diff --git a/include/FlavorTag.h b/include/FlavorTag.h index 5878394..4be54b3 100644 --- a/include/FlavorTag.h +++ b/include/FlavorTag.h @@ -47,6 +47,8 @@ class FlavorTag : public Algorithm { int _nhitsJointProbZ0; int _nhitsMostSignificantTrack; + double _dEdxGausWidth, _dEdxMinMom, _dEdxMaxAng; + ClassDef(FlavorTag,1); }; diff --git a/include/algoSigProb.h b/include/algoSigProb.h index 913ce3c..6a30d85 100644 --- a/include/algoSigProb.h +++ b/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); } } diff --git a/include/flavtag.h b/include/flavtag.h index 4665a58..98e02a2 100644 --- a/include/flavtag.h +++ b/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; diff --git a/include/lcfiplus.h b/include/lcfiplus.h index 38deb74..6015c87 100644 --- a/include/lcfiplus.h +++ b/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 _pidProbability; double _correnergy; + //If a track is imported from a PFO having multiple tracks or not + bool _isMultiTrack; //BNess double _bness, _cness; diff --git a/src/FlavorTag.cc b/src/FlavorTag.cc index 5caf343..fc3beb7 100644 --- a/src/FlavorTag.cc +++ b/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;igetBNess()); //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;igetBNess()>=0.5) bnessvec += (*tracks[i]); } - _result = bnessvec.M(); } }; @@ -1744,7 +1731,98 @@ class FtNBNess : public FTAlgo { for(unsigned int i=0;igetBNess()>=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() { diff --git a/src/LCIOStorer.cc b/src/LCIOStorer.cc index 0515af5..7819084 100644 --- a/src/LCIOStorer.cc +++ b/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{ diff --git a/src/LcfiplusProcessor.cc b/src/LcfiplusProcessor.cc index 37278fa..ba0e102 100644 --- a/src/LcfiplusProcessor.cc +++ b/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; diff --git a/src/algoSigProb.cc b/src/algoSigProb.cc index bb5e87b..604a971 100644 --- a/src/algoSigProb.cc +++ b/src/algoSigProb.cc @@ -495,7 +495,83 @@ void findMostSignificantTrack(const Jet* jet, const Vertex* pri, int minhitcut, sigVec[5] = 0; } } + // Function to fill observables using the significance of the dEdx/dNdx distance wrt Bethe-Bloch estimation + // Reading particleIDProbability from the lcio files + // Two different options: Counting Number of particles or building ratios. + // Jesus Marquez Hernandez (IFIC-CSIC/UV) + + double dEdxNPartSec(const VertexVec& vtxList, string particle, float GausWidth, float MinMom, float MaxAngle) { + double counter=0; + try{ + for (unsigned int j=0; jgetTracks(); + for(unsigned int i=0;igetParticleIDProbability("kaon_dEdxdistance"); + else if(particle=="pion")distance=vtxTracks.at(i)->getParticleIDProbability("pion_dEdxdistance"); + else if(particle=="proton")distance=vtxTracks.at(i)->getParticleIDProbability("proton_dEdxdistance"); + double mom=vtxTracks.at(i)->P(); + double costheta=vtxTracks.at(i)->CosTheta(); + bool isMultiTrack=vtxTracks.at(i)->isMultiTrack(); + if(distance == 0) continue; // Initialization issue + if(abs(distance) > 100) continue; // Initialization issue if rewritten + if(mom < MinMom) continue; // Energy cut + if(abs(costheta) > MaxAngle) continue; // Angle cut + if(isMultiTrack == true) continue; // Remove possible multitracks + // Fill the particle counter + if(distance > -GausWidth && distance < GausWidth) counter++; + } + } + } + catch (lcfiplus::Exception& e) { + } + return counter; + } + // To use ratios of particles + + double dEdxRatioSec(const VertexVec& vtxList, string P1overP2, float GausWidth, float MinMom, float MaxAngle) { + double ratio=-3; + double neg_counter=0; //estimated protons + double null_counter=0; //estimated kaons + double pos_counter=0; //estimated pion + try{ + for (unsigned int j=0; jgetTracks(); + for(unsigned int i=0;igetParticleIDProbability("kaon_dEdxdistance"); //Following my steps, KDS:Kaon Distance Significance + double mom=vtxTracks.at(i)->P(); + double costheta=vtxTracks.at(i)->CosTheta(); + bool isMultiTrack=vtxTracks.at(i)->isMultiTrack(); + if(KDS == 0) continue; // Initialization issue + if(abs(KDS) > 100) continue; // Initialization issue if rewritten + if(mom < MinMom) continue; // Energy cut + if(abs(costheta) > MaxAngle) continue; // Angle cut + if(isMultiTrack == true) continue; // Remove possible multitracks + // Fill the particle counters + if(KDS > -10 && KDS < -GausWidth) neg_counter++; + if(KDS > -GausWidth && KDS < GausWidth) null_counter++; + if(KDS > GausWidth && KDS < 10) pos_counter++; + } + } + // Fill the ratios: + if(P1overP2 == "PionOverKaon"){ + if(null_counter==0)ratio=-1; // Initialize to not divide by 0 + else ratio=pos_counter/null_counter; + } + else if(P1overP2 == "KaonOverProton"){ + if(neg_counter==0)ratio=-1; // Initialize to not divide by 0 + else ratio=null_counter/neg_counter; + } + else if(P1overP2 == "PionOverProton"){ + if(neg_counter==0)ratio=-1; // Initialize to not divide by 0 + else ratio=pos_counter/neg_counter; + } + } + catch (lcfiplus::Exception& e) { + } + return ratio; + } + } } - diff --git a/src/flavtag.cc b/src/flavtag.cc index b1209a5..bd8ad53 100644 --- a/src/flavtag.cc +++ b/src/flavtag.cc @@ -165,6 +165,13 @@ void FTManager::process(const Event* event, const Vertex* privtx, int nhitsJoint } } +//KDS Parameters +void FTManager::setDEDXParameters(const double gaus, const double mom, const double ang) { + _dEdxGausWidth = gaus; + _dEdxMinMom = mom; + _dEdxMaxAng = ang; +} + float* FTManager::getVarAddress(const string& varname) { std::vector::iterator iter; for (iter = _algoList.begin(); iter != _algoList.end(); ++iter) {