PID implementation

include/FlavorTag.h

@@ -47,6 +47,8 @@ class FlavorTag : public Algorithm {
   int _nhitsJointProbZ0;
   int _nhitsMostSignificantTrack;
 
+  double _dEdxGausWidth, _dEdxMinMom, _dEdxMaxAng;
+
   ClassDef(FlavorTag,1);
 };
 

include/algoSigProb.h

@@ -28,6 +28,10 @@ extern double jointProbZ0(const Jet* jet, const Vertex* pri, int minhitcut, doub
 extern double jointProb2D0(const Jet* jet, const Vertex* pri, int minhitcut, double maxd0sigcut, bool useVertexTracks, const TH1* jh1, const TH1* jh2);
 extern double jointProb2Z0(const Jet* jet, const Vertex* pri, int minhitcut, double maxz0sigcut, bool useVertexTracks, const TH1* jh1, const TH1* jh2);
 
+extern double dEdxNPartPri(const Vertex* pri, string particle, float GausWidth, float MaxMom, float MaxAngle);
+extern double dEdxNPartSec(const VertexVec& vtxList, string particle, float GausWidth, float MaxMom, float MaxAngle);
+extern double dEdxRatioPri(const Vertex* pri, string P1overP2, float GausWidth, float MaxMom, float MaxAngle);
+extern double dEdxRatioSec(const VertexVec& vtxList, string P1overP2, float GausWidth, float MaxMom, float MaxAngle);
 }
 }
 

include/flavtag.h

@@ -42,6 +42,7 @@ class FTAlgo {
   void setNHitsJointProbZ0(int value);
   void setNHitsMostSignificantTrack(int value);
   float getValue();
+
   const string& getName() const {
     return _name;
   }
@@ -87,6 +88,8 @@ class FTManager {
   void openFile(const char* filename);
   void closeFile();
   void process(const Event* event, const Vertex* privtx, int nhitsJointProbD0, int nhitsJointProbZ0, int nhitsMostSignificantTrack, JetVec& jets);
+  //dEdx-KDS
+  void setDEDXParameters(const double gaus, const double mom, const double ang);
 
   float* getVarAddress(const string& varname);
   void setEval(bool seteval, bool exportAllVars) {
@@ -114,6 +117,18 @@ class FTManager {
   void setAuxiliary(double aux) {
     _aux = aux;
   }
+
+  double getDEDXGausWidth()const {
+    return _dEdxGausWidth;
+  }
+
+  double getDEDXMinMom()const {
+    return _dEdxMinMom;
+  }
+
+  double getDEDXMaxAng()const {
+    return _dEdxMaxAng;
+  }
 
  private:
   void add(FTAlgo* v); // 121214 moved to private
@@ -136,9 +151,11 @@ class FTManager {
   // variables for flavor tagging
   FtIPProbHolder* _holder;
   double _aux;
+  // dEdx-KDS
+  double _dEdxGausWidth, _dEdxMinMom, _dEdxMaxAng;
 };
 
-// historgram holder for d0/z0 probability
+// histogram holder for d0/z0 probability
 class FtD0bProb;
 class FtD0cProb;
 class FtD0qProb;

include/lcfiplus.h

@@ -414,6 +414,13 @@ class Track : public TLorentzVector {//, protected TrackData {//, public EventPo
     _pdg = pdg;
   }
 
+  bool isMultiTrack() const {
+    return _isMultiTrack;
+  }
+  void setMultiTrack(bool multi) {
+    _isMultiTrack = multi;
+  }
+
   double getCharge() const {
     return _charge;
   }
@@ -584,6 +591,8 @@ class Track : public TLorentzVector {//, protected TrackData {//, public EventPo
   //ParticleID posterior probability
   map<string, double> _pidProbability;
   double _correnergy;
+  //If a track is imported from a PFO having multiple tracks or not
+  bool _isMultiTrack;
 
   //BNess
   double _bness, _cness;

src/FlavorTag.cc

@@ -1648,12 +1648,9 @@ class FtBNess0 : public FTAlgo {
 	bness.push_back(tracks[i]->getBNess());
  	//bness.push_back(0.5*tracks[i]->getBNess()+0.5*tracks[i]->getCNess());
       }
-
       std::sort(bness.begin(), bness.end());
-
       _result = bness[tracks.size()-1];
     }
-
   }
 };
 
@@ -1668,13 +1665,10 @@ class FtBNess1 : public FTAlgo {
       for(unsigned int i=0;i<tracks.size();i++){
 	bness.push_back(tracks[i]->getBNess());
  	//bness.push_back(0.5*tracks[i]->getBNess()+0.5*tracks[i]->getCNess());
-      }
-
+      } 
       std::sort(bness.begin(), bness.end());
-
       _result = bness[tracks.size()-2];
     }
-
   }
 };
 
@@ -1690,12 +1684,9 @@ class FtBNess2 : public FTAlgo {
 	bness.push_back(tracks[i]->getBNess());
  	//bness.push_back(0.5*tracks[i]->getBNess()+0.5*tracks[i]->getCNess());
       }
-
-      std::sort(bness.begin(), bness.end());
-
+      std::sort(bness.begin(), bness.end()); 
       _result = bness[tracks.size()-3];
-    }
-
+    }	
   }
 };
 
@@ -1711,12 +1702,9 @@ class FtBNess3 : public FTAlgo {
 	bness.push_back(tracks[i]->getBNess());
  	//bness.push_back(0.5*tracks[i]->getBNess()+0.5*tracks[i]->getCNess());
       }
-
-      std::sort(bness.begin(), bness.end());
-
+      std::sort(bness.begin(), bness.end()); 
       _result = bness[tracks.size()-4];
-    }
-
+    }	
   }
 };
 
@@ -1730,7 +1718,6 @@ class FtBNessMass : public FTAlgo {
     for(unsigned int i=0;i<tracks.size();i++){
       if(tracks[i]->getBNess()>=0.5) bnessvec += (*tracks[i]);
     }
-
     _result = bnessvec.M();
   }
 };
@@ -1744,7 +1731,98 @@ class FtNBNess : public FTAlgo {
     for(unsigned int i=0;i<tracks.size();i++){
       if(tracks[i]->getBNess()>=0.5) _result ++;
     }
-
+  }
+};
+
+  //dEdx
+class dEdxNKaonSec : public FTAlgo {
+ public:
+  dEdxNKaonSec() : FTAlgo("dEdxNKaonSec") {}
+  void process() {
+    _result = 0;
+    const VertexVec& vtxList = _jet->getVertices();
+    if(vtxList.size()>0){
+      double gaus = FTManager::getInstance().getDEDXGausWidth();
+      double mom = FTManager::getInstance().getDEDXMinMom();
+      double ang = FTManager::getInstance().getDEDXMaxAng();
+      _result=dEdxNPartSec( vtxList, string("kaon"), gaus, mom, ang);
+    }
+  }
+};
+
+class dEdxNPionSec : public FTAlgo {
+public:
+  dEdxNPionSec() : FTAlgo("dEdxNPionSec") {}
+  void process() {
+    _result = 0;
+    const VertexVec& vtxList = _jet->getVertices();
+    if(vtxList.size()>0){
+      double gaus = FTManager::getInstance().getDEDXGausWidth();
+      double mom = FTManager::getInstance().getDEDXMinMom();
+      double ang = FTManager::getInstance().getDEDXMaxAng();
+      _result=dEdxNPartSec( vtxList, string("pion"), gaus, mom, ang);
+    }
+  }
+};
+
+class dEdxNProtonSec : public FTAlgo {
+public:
+  dEdxNProtonSec() : FTAlgo("dEdxNProtonSec") {}
+  void process() {
+    _result = 0;
+    const VertexVec& vtxList = _jet->getVertices();
+    if(vtxList.size()>0){
+      double gaus = FTManager::getInstance().getDEDXGausWidth();
+      double mom = FTManager::getInstance().getDEDXMinMom();
+      double ang = FTManager::getInstance().getDEDXMaxAng();
+      _result=dEdxNPartSec( vtxList, string("proton"), gaus, mom, ang);
+    }
+  }
+};
+
+  //dEdx Ratios
+class dEdxRatioPionOverKaonSec : public FTAlgo {
+ public:
+  dEdxRatioPionOverKaonSec() : FTAlgo("dEdxRatioPionOverKaonSec") {}
+  void process() {
+    _result = -2; //Entry in -2 means no secondary vtx
+    const VertexVec& vtxList = _jet->getVertices(); //We also count pseudo-vtx
+    if(vtxList.size()>0){
+      double gaus = FTManager::getInstance().getDEDXGausWidth();
+      double mom = FTManager::getInstance().getDEDXMinMom();
+      double ang = FTManager::getInstance().getDEDXMaxAng();
+      _result=dEdxRatioSec( vtxList, string("PionOverKaon"), gaus, mom, ang);
+    }
+  }
+};
+
+class dEdxRatioKaonOverProtonSec : public FTAlgo {
+ public:
+  dEdxRatioKaonOverProtonSec() : FTAlgo("dEdxRatioKaonOverProtonSec") {}
+  void process() {
+    _result = -2;
+    const VertexVec& vtxList = _jet->getVertices();
+      if(vtxList.size()>0){
+	double gaus = FTManager::getInstance().getDEDXGausWidth();
+	double mom = FTManager::getInstance().getDEDXMinMom();
+	double ang = FTManager::getInstance().getDEDXMaxAng();
+        _result=dEdxRatioSec( vtxList, string("KaonOverProton"), gaus, mom, ang);
+      }
+  }
+};
+
+class dEdxRatioPionOverProtonSec : public FTAlgo {
+ public:
+  dEdxRatioPionOverProtonSec() : FTAlgo("dEdxRatioPionOverProtonSec") {}
+  void process() {
+    _result = -2;
+    const VertexVec& vtxList = _jet->getVertices();
+    if(vtxList.size()>0){
+      double gaus = FTManager::getInstance().getDEDXGausWidth();
+      double mom = FTManager::getInstance().getDEDXMinMom();
+      double ang = FTManager::getInstance().getDEDXMaxAng();
+      _result=dEdxRatioSec( vtxList, string("PionOverProton"), gaus, mom, ang);
+    }
   }
 };
 
@@ -1859,7 +1937,6 @@ void FTManager::initVars() {
   add( new FtZ0qProb(false));
   add( new FtZ0bProbIP());
   add( new FtZ0cProbIP());
-
   add( new FtD0bProb2());
   add( new FtD0cProb2());
   add( new FtD0qProb2());
@@ -1890,6 +1967,15 @@ void FTManager::initVars() {
   add( new FtBNess3() ); 
   add( new FtBNessMass() ); 
   add( new FtNBNess() ); 
+
+  //dEdx Variables
+  add( new dEdxNKaonSec() );
+  add( new dEdxNPionSec() );
+  add( new dEdxNProtonSec() );
+  add( new dEdxRatioPionOverKaonSec() );
+  add( new dEdxRatioPionOverKaonSec() );
+  add( new dEdxRatioKaonOverProtonSec() );
+  add( new dEdxRatioPionOverProtonSec() );
 }
 
 void FlavorTag::init(Parameters* param) {
@@ -1909,7 +1995,11 @@ void FlavorTag::init(Parameters* param) {
   _nhitsJointProbD0 = param->get("FlavorTag.NVTXhitsJointProbD0", int(4));
   _nhitsJointProbZ0 = param->get("FlavorTag.NVTXhitsJointProbZ0", int(4));
   _nhitsMostSignificantTrack = param->get("FlavorTag.NhitsMostSignificantTrack", int(4));
-
+  //dEdx
+  _dEdxGausWidth = param->get("FlavorTag.DEDXGaussianWidth",double(4));
+  _dEdxMinMom = param->get("FlavorTag.DEDXMinMomentum",double(3));
+  _dEdxMaxAng = param->get("FlavorTag.DEDXMaxAngle",double(0.95));
+
   //string outputFilename = param->get("TrainNtupleFile",string("lcfiplus.root"));
   //_nJet = (int)param->get("TrainNJet",float(2));
 
@@ -1943,6 +2033,7 @@ void FlavorTag::process() {
   mgr.setIPProbHolder(_holder);
 
   mgr.process(event, privtx, _nhitsJointProbD0, _nhitsJointProbZ0, _nhitsMostSignificantTrack, jets);
+  mgr.setDEDXParameters(_dEdxGausWidth, _dEdxMinMom, _dEdxMaxAng);
 }
 
 void FlavorTag::end() {

src/LCIOStorer.cc

@@ -438,6 +438,10 @@ void LCIOStorer::SetEvent(lcio::LCEvent* evt) {
         track->SetE(pfo->getEnergy());
         TVector3 pTrack(pfo->getMomentum());
         track->SetVect(pTrack);
+
+	// Check if the PFO has multiple tracks or not
+	if(pfo->getTracks().size() > 1) track->setMultiTrack(true);
+	else track->setMultiTrack(false);
 
 	//PIDs
 	try{

src/LcfiplusProcessor.cc

@@ -53,14 +53,15 @@ LcfiplusProcessor::LcfiplusProcessor() : Processor("LcfiplusProcessor") {
   registerProcessorParameter("IgnoreLackOfVertexRP", "Keep running even if vertex RP collection is not present",
                              _ignoreLackOfVertexRP, int(0));
   registerProcessorParameter("PrintEventNumber", "Event number printing period in std output: 0 = no printing", _printPeriod, int(0));
+  // PID
   registerProcessorParameter("PIDAlgorithmName", "ParticleID Algorithm Name", _pidAlgoName, string("LikelihoodPID"));
 
+  // Beam
   registerOptionalParameter("MagneticField", "Manually set magnetic field, overriding the value from DD4hep [T]", _magneticField, float(0.0));
   registerOptionalParameter("BeamSizeX", "Bunch size in the X direction [mm]", _beamSizeX, float(639e-6));
   registerOptionalParameter("BeamSizeY", "Bunch size in the Y direction [mm]", _beamSizeY, float(5.7e-6));
   registerOptionalParameter("BeamSizeZ", "Bunch size in the Z direction [mm]", _beamSizeZ, float(9.13e-2));
 
-
   // ROOT object
   /*	int argc = 0;
   	if(gROOT->GetApplication() == 0){
@@ -139,6 +140,7 @@ void LcfiplusProcessor::init() {
       _param->add(keys[i].c_str(), vals);
     }
 
+
     // initialize LCIOStorer
     if (!_lcio) {
       _lcio = new LCIOStorer(0,0,true,false,0); // no file
@@ -147,7 +149,6 @@ void LcfiplusProcessor::init() {
       _lcio->setUpdateVertexRPDaughters(_updateVertexRPDaughters);
       _lcio->setIgnoreLackOfVertexRP(_ignoreLackOfVertexRP);
       _lcio->setParticleIDAlgorithmName(_pidAlgoName.c_str());
-
       _lcioowner = true;
     } else {
       _lcioowner = false;