From 55c4e912a38de451b1008b2bc93299630045b0d4 Mon Sep 17 00:00:00 2001 From: Jesus Pedro Marquez Hernandez Date: Thu, 2 Feb 2023 18:21:46 +0100 Subject: [PATCH 01/13] beta version implementing dEdx --- include/algoSigProb.h | 3 ++ include/lcfiplus.h | 14 +++++++ src/FlavorTag.cc | 75 +++++++++++++++++++++++++++++++++ src/LCIOStorer.cc | 5 +++ src/LcfiplusProcessor.cc | 7 +++- src/algoSigProb.cc | 89 +++++++++++++++++++++++++++++++++++++++- 6 files changed, 190 insertions(+), 3 deletions(-) diff --git a/include/algoSigProb.h b/include/algoSigProb.h index 913ce3ca..2c4455ee 100644 --- a/include/algoSigProb.h +++ b/include/algoSigProb.h @@ -28,6 +28,9 @@ 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); +//dEdx + extern double dEdxKDSRatioPri(const Vertex* pri, string P1overP2); + extern double dEdxKDSRatioSec(const VertexVec& vtxList, string P1overP2); } } diff --git a/include/lcfiplus.h b/include/lcfiplus.h index 38deb740..dafc4732 100644 --- a/include/lcfiplus.h +++ b/include/lcfiplus.h @@ -414,6 +414,15 @@ class Track : public TLorentzVector {//, protected TrackData {//, public EventPo _pdg = pdg; } + // + double getCharge() const { return _charge; } @@ -588,6 +597,11 @@ class Track : public TLorentzVector {//, protected TrackData {//, public EventPo //BNess double _bness, _cness; + // + ClassDef(lcfiplus::Track, 2); }; diff --git a/src/FlavorTag.cc b/src/FlavorTag.cc index 5caf3435..dbc5fa14 100644 --- a/src/FlavorTag.cc +++ b/src/FlavorTag.cc @@ -1748,6 +1748,71 @@ class FtNBNess : public FTAlgo { } }; + //getVertices(); //getVerticesForFT()? + // I believe getVertices gets you all and getVerticesForFT don't count pseudovtx + if(vtxList.size()>0){ + _result=dEdxKDSRatioSec( vtxList, string("PionOverKaon")); + } + } + }; + + class dEdxRatioKaonOverProtonSec : public FTAlgo { + public: + dEdxRatioKaonOverProtonSec() : FTAlgo("dEdxRatioKaonOverProtonSec") {} + void process() { + _result = -2; + const VertexVec& vtxList = _jet->getVertices(); + if(vtxList.size()>0){ + _result=dEdxKDSRatioSec( vtxList, string("KaonOverProton")); + } + } + }; + + class dEdxRatioPionOverProtonSec : public FTAlgo { + public: + dEdxRatioPionOverProtonSec() : FTAlgo("dEdxRatioPionOverProtonSec") {} + void process() { + _result = -2; + const VertexVec& vtxList = _jet->getVertices(); + if(vtxList.size()>0){ + _result=dEdxKDSRatioSec( vtxList, string("PionOverProton")); + } + } + }; + + //Added by JP> + void FTManager::initVars() { if (_initialized)return; _initialized = true; @@ -1890,8 +1955,18 @@ void FTManager::initVars() { add( new FtBNess3() ); add( new FtBNessMass() ); add( new FtNBNess() ); + + //dEdx Variables + add( new dEdxRatioPionOverKaonSec() ); + add( new dEdxRatioKaonOverProtonSec() ); + add( new dEdxRatioPionOverProtonSec() ); + add( new dEdxRatioPionOverKaonPri() ); + add( new dEdxRatioKaonOverProtonPri() ); + add( new dEdxRatioPionOverProtonPri() ); } + + void FlavorTag::init(Parameters* param) { Algorithm::init(param); diff --git a/src/LCIOStorer.cc b/src/LCIOStorer.cc index 0515af54..ab273385 100644 --- a/src/LCIOStorer.cc +++ b/src/LCIOStorer.cc @@ -438,6 +438,11 @@ void LCIOStorer::SetEvent(lcio::LCEvent* evt) { track->SetE(pfo->getEnergy()); TVector3 pTrack(pfo->getMomentum()); track->SetVect(pTrack); + + //getTracks().size() > 1) track->setIsUnique(false); + else track->setIsUnique(true); + //JP> //PIDs try{ diff --git a/src/LcfiplusProcessor.cc b/src/LcfiplusProcessor.cc index 37278fa6..fb7bba92 100644 --- a/src/LcfiplusProcessor.cc +++ b/src/LcfiplusProcessor.cc @@ -53,14 +53,17 @@ 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)); + + //registerProcessorParameter("PIDAlgorithmName", "ParticleID Algorithm Name", _pidAlgoName, string("")); registerProcessorParameter("PIDAlgorithmName", "ParticleID Algorithm Name", _pidAlgoName, string("LikelihoodPID")); - + //registerProcessorParameter("PIDAlgorithmName", "ParticleID Algorithm Name", _pidAlgoName, string("dEdxPIDv2")); + 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){ diff --git a/src/algoSigProb.cc b/src/algoSigProb.cc index bb5e87b9..4c987da9 100644 --- a/src/algoSigProb.cc +++ b/src/algoSigProb.cc @@ -496,6 +496,93 @@ void findMostSignificantTrack(const Jet* jet, const Vertex* pri, int minhitcut, } } +//getTracks(); // We load the primary vtx + try{ + 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 isunique=vtxTracks.at(i)->getIsUnique(); + if(KDS == 0) continue; // Initialization issue + if(mom < 3) continue; // Energy cut + if(costheta > 0.95) continue; // Angle cut + if(isunique == false) continue; + + if(KDS > -5 && KDS < -2) neg_counter++; + if(KDS > -2 && KDS < 2) null_counter++; + if(KDS > 2 && KDS < 5) pos_counter++; + } + // Now we fill the ratio: + 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; + + } + + double dEdxKDSRatioSec(const VertexVec& vtxList, string P1overP2) { + + 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 isunique=vtxTracks.at(i)->getIsUnique(); + if(KDS == 0) continue; // Initialization issue + if(mom < 3) continue; // Energy cut + if(costheta > 0.95) continue; // Angle cut + if(isunique == false) continue; + + if(KDS > -5 && KDS < -2) neg_counter++; + if(KDS > -2 && KDS < 2) null_counter++; + if(KDS > 2 && KDS < 5) pos_counter++; + } + } + // Now we fill the ratio: + 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; + } + } } - From f95f7b1b24d34abaf4a36fa19f19b26712729640 Mon Sep 17 00:00:00 2001 From: Jesus Pedro Marquez Hernandez Date: Thu, 9 Feb 2023 00:17:28 +0100 Subject: [PATCH 02/13] Changing IsUnique to isMultiTrack, fixing abs value in angular cut for KDS and turning features for dEdx-KDS into parameters for the FlavourTag processor --- include/FlavorTag.h | 2 ++ include/LcfiplusProcessor.h | 2 ++ include/algoSigProb.h | 5 ++--- include/flavtag.h | 5 +++++ include/lcfiplus.h | 17 ++++++----------- src/FlavorTag.cc | 27 ++++++++++++--------------- src/LCIOStorer.cc | 9 +++++---- src/LcfiplusProcessor.cc | 10 ++++------ src/algoSigProb.cc | 37 ++++++++++++++++++------------------- 9 files changed, 56 insertions(+), 58 deletions(-) diff --git a/include/FlavorTag.h b/include/FlavorTag.h index 5878394d..c1122af7 100644 --- a/include/FlavorTag.h +++ b/include/FlavorTag.h @@ -47,6 +47,8 @@ class FlavorTag : public Algorithm { int _nhitsJointProbZ0; int _nhitsMostSignificantTrack; + double _kdsGausWidth, _kdsMinMom, _kdsMaxAng; + ClassDef(FlavorTag,1); }; diff --git a/include/LcfiplusProcessor.h b/include/LcfiplusProcessor.h index 17b615b4..fa6367f5 100644 --- a/include/LcfiplusProcessor.h +++ b/include/LcfiplusProcessor.h @@ -95,6 +95,8 @@ class LcfiplusProcessor : public Processor, public lcfiplus::EventStoreObserver //Particle ID Algorithm Name string _pidAlgoName; + //dEdx parameters + } ; #endif diff --git a/include/algoSigProb.h b/include/algoSigProb.h index 2c4455ee..4f9dab65 100644 --- a/include/algoSigProb.h +++ b/include/algoSigProb.h @@ -28,9 +28,8 @@ 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); -//dEdx - extern double dEdxKDSRatioPri(const Vertex* pri, string P1overP2); - extern double dEdxKDSRatioSec(const VertexVec& vtxList, string P1overP2); +extern double dEdxKDSRatioPri(const Vertex* pri, string P1overP2, float GausWidth, float MaxMom, float MaxAngle); +extern double dEdxKDSRatioSec(const VertexVec& vtxList, string P1overP2, float GausWidth, float MaxMom, float MaxAngle); } } diff --git a/include/flavtag.h b/include/flavtag.h index 4665a580..3d89fade 100644 --- a/include/flavtag.h +++ b/include/flavtag.h @@ -42,6 +42,10 @@ class FTAlgo { void setNHitsJointProbZ0(int value); void setNHitsMostSignificantTrack(int value); float getValue(); + // + const string& getName() const { return _name; } @@ -56,6 +60,7 @@ class FTAlgo { int _nhitsJointProbD0; int _nhitsJointProbZ0; int _nhitsMostSignificantTrack; + double _kdsGausWidth, _kdsMinMom, _kdsMaxAng; float _result; string _name; diff --git a/include/lcfiplus.h b/include/lcfiplus.h index dafc4732..81c5366e 100644 --- a/include/lcfiplus.h +++ b/include/lcfiplus.h @@ -414,14 +414,12 @@ class Track : public TLorentzVector {//, protected TrackData {//, public EventPo _pdg = pdg; } - // double getCharge() const { return _charge; @@ -593,15 +591,12 @@ class Track : public TLorentzVector {//, protected TrackData {//, public EventPo //ParticleID posterior probability map _pidProbability; double _correnergy; + // For avoid double-counting in dEdx + bool _isMultiTrack; //BNess double _bness, _cness; - // - ClassDef(lcfiplus::Track, 2); }; diff --git a/src/FlavorTag.cc b/src/FlavorTag.cc index dbc5fa14..f1d1bf2c 100644 --- a/src/FlavorTag.cc +++ b/src/FlavorTag.cc @@ -1748,12 +1748,11 @@ class FtNBNess : public FTAlgo { } }; - //getVertices(); //getVerticesForFT()? - // I believe getVertices gets you all and getVerticesForFT don't count pseudovtx + const VertexVec& vtxList = _jet->getVertices(); //We also count pseudo-vtx if(vtxList.size()>0){ - _result=dEdxKDSRatioSec( vtxList, string("PionOverKaon")); + _result=dEdxKDSRatioSec( vtxList, string("PionOverKaon"), _kdsGausWidth, _kdsMinMom, _kdsMaxAng); } } }; @@ -1794,7 +1791,7 @@ class FtNBNess : public FTAlgo { _result = -2; const VertexVec& vtxList = _jet->getVertices(); if(vtxList.size()>0){ - _result=dEdxKDSRatioSec( vtxList, string("KaonOverProton")); + _result=dEdxKDSRatioSec( vtxList, string("KaonOverProton"), _kdsGausWidth, _kdsMinMom, _kdsMaxAng); } } }; @@ -1806,13 +1803,11 @@ class FtNBNess : public FTAlgo { _result = -2; const VertexVec& vtxList = _jet->getVertices(); if(vtxList.size()>0){ - _result=dEdxKDSRatioSec( vtxList, string("PionOverProton")); + _result=dEdxKDSRatioSec( vtxList, string("PionOverProton"), _kdsGausWidth, _kdsMinMom, _kdsMaxAng); } } }; - //Added by JP> - void FTManager::initVars() { if (_initialized)return; _initialized = true; @@ -1965,8 +1960,6 @@ void FTManager::initVars() { add( new dEdxRatioPionOverProtonPri() ); } - - void FlavorTag::init(Parameters* param) { Algorithm::init(param); @@ -1984,7 +1977,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-KDS + _kdsGausWidth = param->get("FlavorTag.KDSGaussianWidth",double(4)); + _kdsMinMom = param->get("FlavorTag.KDSMinMomentum",double(3)); + _kdsMaxAng = param->get("FlavorTag.KDSMaxAngle",double(0.95)); + //string outputFilename = param->get("TrainNtupleFile",string("lcfiplus.root")); //_nJet = (int)param->get("TrainNJet",float(2)); diff --git a/src/LCIOStorer.cc b/src/LCIOStorer.cc index ab273385..39cbcbd0 100644 --- a/src/LCIOStorer.cc +++ b/src/LCIOStorer.cc @@ -439,13 +439,13 @@ void LCIOStorer::SetEvent(lcio::LCEvent* evt) { TVector3 pTrack(pfo->getMomentum()); track->SetVect(pTrack); - //getTracks().size() > 1) track->setIsUnique(false); - else track->setIsUnique(true); - //JP> + // To avoid double-counting in dEdx + if(pfo->getTracks().size() > 1) track->setMultiTrack(true); + else track->setMultiTrack(false); //PIDs try{ + int pidAlgoID = PID.getAlgorithmID(_pidAlgoName); //pdg value track->setPDG(PID.getParticleID(pfo,pidAlgoID).getPDG()); @@ -458,6 +458,7 @@ void LCIOStorer::SetEvent(lcio::LCEvent* evt) { //cal. corrected mass track->setCorrEnergy(pmass[PID.getParticleID(pfo,pidAlgoID).getPDG()]); //track->swapEnergy(); //really temporal need flag... + }catch(UTIL::UnknownAlgorithm e){ } diff --git a/src/LcfiplusProcessor.cc b/src/LcfiplusProcessor.cc index fb7bba92..ba0e102f 100644 --- a/src/LcfiplusProcessor.cc +++ b/src/LcfiplusProcessor.cc @@ -53,16 +53,14 @@ 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)); - - //registerProcessorParameter("PIDAlgorithmName", "ParticleID Algorithm Name", _pidAlgoName, string("")); + // PID registerProcessorParameter("PIDAlgorithmName", "ParticleID Algorithm Name", _pidAlgoName, string("LikelihoodPID")); - //registerProcessorParameter("PIDAlgorithmName", "ParticleID Algorithm Name", _pidAlgoName, string("dEdxPIDv2")); - + + // 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; @@ -142,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 @@ -150,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 4c987da9..6074bfc2 100644 --- a/src/algoSigProb.cc +++ b/src/algoSigProb.cc @@ -496,10 +496,7 @@ void findMostSignificantTrack(const Jet* jet, const Vertex* pri, int minhitcut, } } -//getParticleIDProbability("kaon_dEdxdistance"); // Following my steps, KDS:Kaon Distance Significance double mom=vtxTracks.at(i)->P(); double costheta=vtxTracks.at(i)->CosTheta(); - bool isunique=vtxTracks.at(i)->getIsUnique(); + bool isMultiTrack=vtxTracks.at(i)->isMultiTrack(); if(KDS == 0) continue; // Initialization issue - if(mom < 3) continue; // Energy cut - if(costheta > 0.95) continue; // Angle cut - if(isunique == false) continue; + if(mom < MaxMom) continue; // Energy cut + if(abs(costheta) > MaxAngle) continue; // Angle cut + if(isMultiTrack == true) continue; - if(KDS > -5 && KDS < -2) neg_counter++; - if(KDS > -2 && KDS < 2) null_counter++; - if(KDS > 2 && KDS < 5) pos_counter++; + if(KDS > -1.5*GausWidth && KDS < -0.5*GausWidth) neg_counter++; + if(KDS > -0.5*GausWidth && KDS < 0.5*GausWidth) null_counter++; + if(KDS > 0.5*GausWidth && KDS < 1.5*GausWidth) pos_counter++; } // Now we fill the ratio: if(P1overP2 == "PionOverKaon"){ @@ -534,6 +531,7 @@ void findMostSignificantTrack(const Jet* jet, const Vertex* pri, int minhitcut, if(neg_counter==0)ratio=-1; // Initialize to not divide by 0 else ratio=pos_counter/neg_counter; } + throw(Exception("PID Collection not found. Don't use dEdx observables.")); } catch (lcfiplus::Exception& e) { } @@ -541,7 +539,7 @@ void findMostSignificantTrack(const Jet* jet, const Vertex* pri, int minhitcut, } - double dEdxKDSRatioSec(const VertexVec& vtxList, string P1overP2) { + double dEdxKDSRatioSec(const VertexVec& vtxList, string P1overP2, float GausWidth, float MaxMom, float MaxAngle) { double ratio=-3; double neg_counter=0; //estimated protons @@ -554,15 +552,15 @@ void findMostSignificantTrack(const Jet* jet, const Vertex* pri, int minhitcut, double KDS=vtxTracks.at(i)->getParticleIDProbability("kaon_dEdxdistance"); //Following my steps, KDS:Kaon Distance Significance double mom=vtxTracks.at(i)->P(); double costheta=vtxTracks.at(i)->CosTheta(); - bool isunique=vtxTracks.at(i)->getIsUnique(); + bool isMultiTrack=vtxTracks.at(i)->isMultiTrack(); if(KDS == 0) continue; // Initialization issue - if(mom < 3) continue; // Energy cut - if(costheta > 0.95) continue; // Angle cut - if(isunique == false) continue; + if(mom < MaxMom) continue; // Energy cut + if(abs(costheta) > MaxAngle) continue; // Angle cut + if(isMultiTrack == true) continue; - if(KDS > -5 && KDS < -2) neg_counter++; - if(KDS > -2 && KDS < 2) null_counter++; - if(KDS > 2 && KDS < 5) pos_counter++; + if(KDS > -1.5*GausWidth && KDS < -0.5*GausWidth) neg_counter++; + if(KDS > -0.5*GausWidth && KDS < 0.5*GausWidth) null_counter++; + if(KDS > 0.5*GausWidth && KDS < 1.5*GausWidth) pos_counter++; } } // Now we fill the ratio: @@ -578,6 +576,7 @@ void findMostSignificantTrack(const Jet* jet, const Vertex* pri, int minhitcut, if(neg_counter==0)ratio=-1; // Initialize to not divide by 0 else ratio=pos_counter/neg_counter; } + throw(Exception("PID Collection not found. Don't use dEdx observables.")); } catch (lcfiplus::Exception& e) { } From 09a2c71c8f7865d900f1d9f5fbc0741925e2ed7a Mon Sep 17 00:00:00 2001 From: Jesus Pedro Marquez Hernandez Date: Thu, 9 Feb 2023 00:36:36 +0100 Subject: [PATCH 03/13] Cosmetic changes --- include/LcfiplusProcessor.h | 2 -- include/flavtag.h | 3 --- include/lcfiplus.h | 2 +- src/LCIOStorer.cc | 3 +-- 4 files changed, 2 insertions(+), 8 deletions(-) diff --git a/include/LcfiplusProcessor.h b/include/LcfiplusProcessor.h index fa6367f5..17b615b4 100644 --- a/include/LcfiplusProcessor.h +++ b/include/LcfiplusProcessor.h @@ -95,8 +95,6 @@ class LcfiplusProcessor : public Processor, public lcfiplus::EventStoreObserver //Particle ID Algorithm Name string _pidAlgoName; - //dEdx parameters - } ; #endif diff --git a/include/flavtag.h b/include/flavtag.h index 3d89fade..d2ff309a 100644 --- a/include/flavtag.h +++ b/include/flavtag.h @@ -42,9 +42,6 @@ class FTAlgo { void setNHitsJointProbZ0(int value); void setNHitsMostSignificantTrack(int value); float getValue(); - // const string& getName() const { return _name; diff --git a/include/lcfiplus.h b/include/lcfiplus.h index 81c5366e..2e7182e2 100644 --- a/include/lcfiplus.h +++ b/include/lcfiplus.h @@ -591,7 +591,7 @@ class Track : public TLorentzVector {//, protected TrackData {//, public EventPo //ParticleID posterior probability map _pidProbability; double _correnergy; - // For avoid double-counting in dEdx + //To avoid double-counting in dEdx bool _isMultiTrack; //BNess diff --git a/src/LCIOStorer.cc b/src/LCIOStorer.cc index 39cbcbd0..e297289f 100644 --- a/src/LCIOStorer.cc +++ b/src/LCIOStorer.cc @@ -445,7 +445,6 @@ void LCIOStorer::SetEvent(lcio::LCEvent* evt) { //PIDs try{ - int pidAlgoID = PID.getAlgorithmID(_pidAlgoName); //pdg value track->setPDG(PID.getParticleID(pfo,pidAlgoID).getPDG()); @@ -458,7 +457,7 @@ void LCIOStorer::SetEvent(lcio::LCEvent* evt) { //cal. corrected mass track->setCorrEnergy(pmass[PID.getParticleID(pfo,pidAlgoID).getPDG()]); //track->swapEnergy(); //really temporal need flag... - + throw(Exception("PID Collection not found.")); }catch(UTIL::UnknownAlgorithm e){ } From d7ac95aabdd804905a99157a659404dcda5bc4b3 Mon Sep 17 00:00:00 2001 From: Jesus Pedro Marquez Hernandez Date: Thu, 9 Feb 2023 01:00:02 +0100 Subject: [PATCH 04/13] Minor bug. --- src/LCIOStorer.cc | 1 - 1 file changed, 1 deletion(-) diff --git a/src/LCIOStorer.cc b/src/LCIOStorer.cc index e297289f..d8c137ad 100644 --- a/src/LCIOStorer.cc +++ b/src/LCIOStorer.cc @@ -457,7 +457,6 @@ void LCIOStorer::SetEvent(lcio::LCEvent* evt) { //cal. corrected mass track->setCorrEnergy(pmass[PID.getParticleID(pfo,pidAlgoID).getPDG()]); //track->swapEnergy(); //really temporal need flag... - throw(Exception("PID Collection not found.")); }catch(UTIL::UnknownAlgorithm e){ } From da83398e22b698371e8e1a5658af6cc39fccf195 Mon Sep 17 00:00:00 2001 From: Jesus Pedro Marquez Hernandez Date: Thu, 9 Feb 2023 01:10:05 +0100 Subject: [PATCH 05/13] Minor bug handling exception --- src/algoSigProb.cc | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/algoSigProb.cc b/src/algoSigProb.cc index 6074bfc2..f17034fd 100644 --- a/src/algoSigProb.cc +++ b/src/algoSigProb.cc @@ -531,7 +531,6 @@ void findMostSignificantTrack(const Jet* jet, const Vertex* pri, int minhitcut, if(neg_counter==0)ratio=-1; // Initialize to not divide by 0 else ratio=pos_counter/neg_counter; } - throw(Exception("PID Collection not found. Don't use dEdx observables.")); } catch (lcfiplus::Exception& e) { } @@ -576,7 +575,6 @@ void findMostSignificantTrack(const Jet* jet, const Vertex* pri, int minhitcut, if(neg_counter==0)ratio=-1; // Initialize to not divide by 0 else ratio=pos_counter/neg_counter; } - throw(Exception("PID Collection not found. Don't use dEdx observables.")); } catch (lcfiplus::Exception& e) { } From 2c4598617b6aba5aeada666954fc6eee21dd568f Mon Sep 17 00:00:00 2001 From: Jesus Pedro Marquez Hernandez Date: Thu, 9 Feb 2023 16:00:29 +0100 Subject: [PATCH 06/13] Fixing many bugs. Now the inputs for dEdx-KDS work. --- include/flavtag.h | 21 +++++++++++++++--- src/FlavorTag.cc | 55 +++++++++++++++++++++++++--------------------- src/algoSigProb.cc | 27 ++++++++++++----------- src/flavtag.cc | 7 ++++++ 4 files changed, 69 insertions(+), 41 deletions(-) diff --git a/include/flavtag.h b/include/flavtag.h index d2ff309a..12d77995 100644 --- a/include/flavtag.h +++ b/include/flavtag.h @@ -42,7 +42,7 @@ class FTAlgo { void setNHitsJointProbZ0(int value); void setNHitsMostSignificantTrack(int value); float getValue(); - + const string& getName() const { return _name; } @@ -57,7 +57,6 @@ class FTAlgo { int _nhitsJointProbD0; int _nhitsJointProbZ0; int _nhitsMostSignificantTrack; - double _kdsGausWidth, _kdsMinMom, _kdsMaxAng; float _result; string _name; @@ -89,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 setKDSParameters(const double gaus, const double mom, const double ang); float* getVarAddress(const string& varname); void setEval(bool seteval, bool exportAllVars) { @@ -116,6 +117,18 @@ class FTManager { void setAuxiliary(double aux) { _aux = aux; } + + double getKDSGausWidth()const { + return _kdsGausWidth; + } + + double getKDSMinMom()const { + return _kdsMinMom; + } + + double getKDSMaxAng()const { + return _kdsMaxAng; + } private: void add(FTAlgo* v); // 121214 moved to private @@ -138,9 +151,11 @@ class FTManager { // variables for flavor tagging FtIPProbHolder* _holder; double _aux; + // dEdx-KDS + double _kdsGausWidth, _kdsMinMom, _kdsMaxAng; }; -// historgram holder for d0/z0 probability +// histogram holder for d0/z0 probability class FtD0bProb; class FtD0cProb; class FtD0qProb; diff --git a/src/FlavorTag.cc b/src/FlavorTag.cc index f1d1bf2c..3083b1ef 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,6 @@ class FtNBNess : public FTAlgo { for(unsigned int i=0;igetBNess()>=0.5) _result ++; } - } }; @@ -1752,7 +1738,10 @@ class FtNBNess : public FTAlgo { public: dEdxRatioPionOverKaonPri() : FTAlgo("dEdxRatioPionOverKaonPri") {} void process() { - _result=dEdxKDSRatioPri( _privtx, string("PionOverKaon"), _kdsGausWidth, _kdsMinMom, _kdsMaxAng); + double gaus = FTManager::getInstance().getKDSGausWidth(); + double mom = FTManager::getInstance().getKDSMinMom(); + double ang = FTManager::getInstance().getKDSMaxAng(); + _result=dEdxKDSRatioPri( _privtx, string("PionOverKaon"), gaus, mom, ang); } }; @@ -1760,7 +1749,10 @@ class FtNBNess : public FTAlgo { public: dEdxRatioKaonOverProtonPri() : FTAlgo("dEdxRatioKaonOverProtonPri") {} void process() { - _result=dEdxKDSRatioPri( _privtx, string("KaonOverProton"), _kdsGausWidth, _kdsMinMom, _kdsMaxAng); + double gaus = FTManager::getInstance().getKDSGausWidth(); + double mom = FTManager::getInstance().getKDSMinMom(); + double ang = FTManager::getInstance().getKDSMaxAng(); + _result=dEdxKDSRatioPri( _privtx, string("KaonOverProton"), gaus, mom, ang); } }; @@ -1768,7 +1760,10 @@ class FtNBNess : public FTAlgo { public: dEdxRatioPionOverProtonPri() : FTAlgo("dEdxRatioPionOverProtonPri") {} void process() { - _result=dEdxKDSRatioPri( _privtx, string("PionOverProton"), _kdsGausWidth, _kdsMinMom, _kdsMaxAng); + double gaus = FTManager::getInstance().getKDSGausWidth(); + double mom = FTManager::getInstance().getKDSMinMom(); + double ang = FTManager::getInstance().getKDSMaxAng(); + _result=dEdxKDSRatioPri( _privtx, string("PionOverProton"), gaus, mom, ang); } }; @@ -1779,7 +1774,10 @@ class FtNBNess : public FTAlgo { _result = -2; //Entry in -2 means no secondary vtx const VertexVec& vtxList = _jet->getVertices(); //We also count pseudo-vtx if(vtxList.size()>0){ - _result=dEdxKDSRatioSec( vtxList, string("PionOverKaon"), _kdsGausWidth, _kdsMinMom, _kdsMaxAng); + double gaus = FTManager::getInstance().getKDSGausWidth(); + double mom = FTManager::getInstance().getKDSMinMom(); + double ang = FTManager::getInstance().getKDSMaxAng(); + _result=dEdxKDSRatioSec( vtxList, string("PionOverKaon"), gaus, mom, ang); } } }; @@ -1791,7 +1789,10 @@ class FtNBNess : public FTAlgo { _result = -2; const VertexVec& vtxList = _jet->getVertices(); if(vtxList.size()>0){ - _result=dEdxKDSRatioSec( vtxList, string("KaonOverProton"), _kdsGausWidth, _kdsMinMom, _kdsMaxAng); + double gaus = FTManager::getInstance().getKDSGausWidth(); + double mom = FTManager::getInstance().getKDSMinMom(); + double ang = FTManager::getInstance().getKDSMaxAng(); + _result=dEdxKDSRatioSec( vtxList, string("KaonOverProton"), gaus, mom, ang); } } }; @@ -1803,7 +1804,10 @@ class FtNBNess : public FTAlgo { _result = -2; const VertexVec& vtxList = _jet->getVertices(); if(vtxList.size()>0){ - _result=dEdxKDSRatioSec( vtxList, string("PionOverProton"), _kdsGausWidth, _kdsMinMom, _kdsMaxAng); + double gaus = FTManager::getInstance().getKDSGausWidth(); + double mom = FTManager::getInstance().getKDSMinMom(); + double ang = FTManager::getInstance().getKDSMaxAng(); + _result=dEdxKDSRatioSec( vtxList, string("PionOverProton"), gaus, mom, ang); } } }; @@ -2015,6 +2019,7 @@ void FlavorTag::process() { mgr.setIPProbHolder(_holder); mgr.process(event, privtx, _nhitsJointProbD0, _nhitsJointProbZ0, _nhitsMostSignificantTrack, jets); + mgr.setKDSParameters(_kdsGausWidth, _kdsMinMom, _kdsMaxAng); } void FlavorTag::end() { diff --git a/src/algoSigProb.cc b/src/algoSigProb.cc index f17034fd..be8f9edd 100644 --- a/src/algoSigProb.cc +++ b/src/algoSigProb.cc @@ -495,30 +495,31 @@ void findMostSignificantTrack(const Jet* jet, const Vertex* pri, int minhitcut, sigVec[5] = 0; } } - - double dEdxKDSRatioPri(const Vertex* pri, string P1overP2, float GausWidth, float MaxMom, float MaxAngle) { + // Functions to fill observables using the Kaon dEdx istance in PID processors + // Jesus Marquez Hernandez (IFIC-CSIC/UV) + double dEdxKDSRatioPri(const Vertex* pri, string P1overP2, float GausWidth, float MinMom, float MaxAngle) { double ratio=-3; // If a -3 appear, something went wrong - double neg_counter=0; // Estimated protons + double neg_counter=0; // Estimated protons double null_counter=0; // Estimated kaons double pos_counter=0; // Estimated pion TrackVec& vtxTracks = pri->getTracks(); // We load the primary vtx try{ for(unsigned int i=0;igetParticleIDProbability("kaon_dEdxdistance"); // Following my steps, KDS:Kaon Distance Significance + double KDS=vtxTracks.at(i)->getParticleIDProbability("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(mom < MaxMom) continue; // Energy cut + if(mom < MinMom) continue; // Energy cut if(abs(costheta) > MaxAngle) continue; // Angle cut - if(isMultiTrack == true) continue; - + if(isMultiTrack == true) continue; // Remove possible multitracks + // Fill the particle counters if(KDS > -1.5*GausWidth && KDS < -0.5*GausWidth) neg_counter++; if(KDS > -0.5*GausWidth && KDS < 0.5*GausWidth) null_counter++; if(KDS > 0.5*GausWidth && KDS < 1.5*GausWidth) pos_counter++; } - // Now we fill the ratio: + // Fill the ratios: if(P1overP2 == "PionOverKaon"){ if(null_counter==0)ratio=-1; // Initialize to not divide by 0 else ratio=pos_counter/null_counter; @@ -538,7 +539,7 @@ void findMostSignificantTrack(const Jet* jet, const Vertex* pri, int minhitcut, } - double dEdxKDSRatioSec(const VertexVec& vtxList, string P1overP2, float GausWidth, float MaxMom, float MaxAngle) { + double dEdxKDSRatioSec(const VertexVec& vtxList, string P1overP2, float GausWidth, float MinMom, float MaxAngle) { double ratio=-3; double neg_counter=0; //estimated protons @@ -553,16 +554,16 @@ void findMostSignificantTrack(const Jet* jet, const Vertex* pri, int minhitcut, double costheta=vtxTracks.at(i)->CosTheta(); bool isMultiTrack=vtxTracks.at(i)->isMultiTrack(); if(KDS == 0) continue; // Initialization issue - if(mom < MaxMom) continue; // Energy cut + if(mom < MinMom) continue; // Energy cut if(abs(costheta) > MaxAngle) continue; // Angle cut - if(isMultiTrack == true) continue; - + if(isMultiTrack == true) continue; // Remove possible multitracks + // Fill the particle counters if(KDS > -1.5*GausWidth && KDS < -0.5*GausWidth) neg_counter++; if(KDS > -0.5*GausWidth && KDS < 0.5*GausWidth) null_counter++; if(KDS > 0.5*GausWidth && KDS < 1.5*GausWidth) pos_counter++; } } - // Now we fill the ratio: + // Fill the ratios: if(P1overP2 == "PionOverKaon"){ if(null_counter==0)ratio=-1; // Initialize to not divide by 0 else ratio=pos_counter/null_counter; diff --git a/src/flavtag.cc b/src/flavtag.cc index b1209a58..f3c081cf 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::setKDSParameters(const double gaus, const double mom, const double ang) { + _kdsGausWidth = gaus; + _kdsMinMom = mom; + _kdsMaxAng = ang; +} + float* FTManager::getVarAddress(const string& varname) { std::vector::iterator iter; for (iter = _algoList.begin(); iter != _algoList.end(); ++iter) { From 4f88edd7c7cd488700c8ba39e296fcc085a2270b Mon Sep 17 00:00:00 2001 From: Jesus Pedro Marquez Hernandez Date: Thu, 9 Feb 2023 17:14:36 +0100 Subject: [PATCH 07/13] Changing KDS sigma definition/use. --- src/algoSigProb.cc | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/src/algoSigProb.cc b/src/algoSigProb.cc index be8f9edd..6296d2ca 100644 --- a/src/algoSigProb.cc +++ b/src/algoSigProb.cc @@ -515,9 +515,9 @@ void findMostSignificantTrack(const Jet* jet, const Vertex* pri, int minhitcut, if(abs(costheta) > MaxAngle) continue; // Angle cut if(isMultiTrack == true) continue; // Remove possible multitracks // Fill the particle counters - if(KDS > -1.5*GausWidth && KDS < -0.5*GausWidth) neg_counter++; - if(KDS > -0.5*GausWidth && KDS < 0.5*GausWidth) null_counter++; - if(KDS > 0.5*GausWidth && KDS < 1.5*GausWidth) pos_counter++; + if(KDS > -3*GausWidth && KDS < -GausWidth) neg_counter++; + if(KDS > GausWidth && KDS < GausWidth) null_counter++; + if(KDS > GausWidth && KDS < 3*GausWidth) pos_counter++; } // Fill the ratios: if(P1overP2 == "PionOverKaon"){ @@ -558,9 +558,9 @@ void findMostSignificantTrack(const Jet* jet, const Vertex* pri, int minhitcut, if(abs(costheta) > MaxAngle) continue; // Angle cut if(isMultiTrack == true) continue; // Remove possible multitracks // Fill the particle counters - if(KDS > -1.5*GausWidth && KDS < -0.5*GausWidth) neg_counter++; - if(KDS > -0.5*GausWidth && KDS < 0.5*GausWidth) null_counter++; - if(KDS > 0.5*GausWidth && KDS < 1.5*GausWidth) pos_counter++; + if(KDS > -3*GausWidth && KDS < -GausWidth) neg_counter++; + if(KDS > GausWidth && KDS < GausWidth) null_counter++; + if(KDS > GausWidth && KDS < 3*GausWidth) pos_counter++; } } // Fill the ratios: From 3cf9b2e634da1c2f4a0eaabb7f8207c50174bd21 Mon Sep 17 00:00:00 2001 From: Jesus Pedro Marquez Hernandez Date: Thu, 9 Feb 2023 20:28:37 +0100 Subject: [PATCH 08/13] Fixing a typo. --- src/algoSigProb.cc | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/algoSigProb.cc b/src/algoSigProb.cc index 6296d2ca..a6b3cc9e 100644 --- a/src/algoSigProb.cc +++ b/src/algoSigProb.cc @@ -516,7 +516,7 @@ void findMostSignificantTrack(const Jet* jet, const Vertex* pri, int minhitcut, if(isMultiTrack == true) continue; // Remove possible multitracks // Fill the particle counters if(KDS > -3*GausWidth && KDS < -GausWidth) neg_counter++; - if(KDS > GausWidth && KDS < GausWidth) null_counter++; + if(KDS > -GausWidth && KDS < GausWidth) null_counter++; if(KDS > GausWidth && KDS < 3*GausWidth) pos_counter++; } // Fill the ratios: @@ -559,7 +559,7 @@ void findMostSignificantTrack(const Jet* jet, const Vertex* pri, int minhitcut, if(isMultiTrack == true) continue; // Remove possible multitracks // Fill the particle counters if(KDS > -3*GausWidth && KDS < -GausWidth) neg_counter++; - if(KDS > GausWidth && KDS < GausWidth) null_counter++; + if(KDS > -GausWidth && KDS < GausWidth) null_counter++; if(KDS > GausWidth && KDS < 3*GausWidth) pos_counter++; } } From 6f8607d9338fa06896140c6280f32ca9a56e0863 Mon Sep 17 00:00:00 2001 From: Jesus Pedro Marquez Hernandez Date: Mon, 20 Feb 2023 17:55:21 +0100 Subject: [PATCH 09/13] Change comments for general purposes. Changed limits in KDE counts for new observables. --- include/lcfiplus.h | 2 +- src/LCIOStorer.cc | 2 +- src/algoSigProb.cc | 8 ++++---- 3 files changed, 6 insertions(+), 6 deletions(-) diff --git a/include/lcfiplus.h b/include/lcfiplus.h index 2e7182e2..6015c876 100644 --- a/include/lcfiplus.h +++ b/include/lcfiplus.h @@ -591,7 +591,7 @@ class Track : public TLorentzVector {//, protected TrackData {//, public EventPo //ParticleID posterior probability map _pidProbability; double _correnergy; - //To avoid double-counting in dEdx + //If a track is imported from a PFO having multiple tracks or not bool _isMultiTrack; //BNess diff --git a/src/LCIOStorer.cc b/src/LCIOStorer.cc index d8c137ad..78190846 100644 --- a/src/LCIOStorer.cc +++ b/src/LCIOStorer.cc @@ -439,7 +439,7 @@ void LCIOStorer::SetEvent(lcio::LCEvent* evt) { TVector3 pTrack(pfo->getMomentum()); track->SetVect(pTrack); - // To avoid double-counting in dEdx + // Check if the PFO has multiple tracks or not if(pfo->getTracks().size() > 1) track->setMultiTrack(true); else track->setMultiTrack(false); diff --git a/src/algoSigProb.cc b/src/algoSigProb.cc index a6b3cc9e..568a4333 100644 --- a/src/algoSigProb.cc +++ b/src/algoSigProb.cc @@ -515,9 +515,9 @@ void findMostSignificantTrack(const Jet* jet, const Vertex* pri, int minhitcut, if(abs(costheta) > MaxAngle) continue; // Angle cut if(isMultiTrack == true) continue; // Remove possible multitracks // Fill the particle counters - if(KDS > -3*GausWidth && KDS < -GausWidth) neg_counter++; + if(KDS > -10 && KDS < -GausWidth) neg_counter++; if(KDS > -GausWidth && KDS < GausWidth) null_counter++; - if(KDS > GausWidth && KDS < 3*GausWidth) pos_counter++; + if(KDS > GausWidth && KDS < 10) pos_counter++; } // Fill the ratios: if(P1overP2 == "PionOverKaon"){ @@ -558,9 +558,9 @@ void findMostSignificantTrack(const Jet* jet, const Vertex* pri, int minhitcut, if(abs(costheta) > MaxAngle) continue; // Angle cut if(isMultiTrack == true) continue; // Remove possible multitracks // Fill the particle counters - if(KDS > -3*GausWidth && KDS < -GausWidth) neg_counter++; + if(KDS > -10 && KDS < -GausWidth) neg_counter++; if(KDS > -GausWidth && KDS < GausWidth) null_counter++; - if(KDS > GausWidth && KDS < 3*GausWidth) pos_counter++; + if(KDS > GausWidth && KDS < 10) pos_counter++; } } // Fill the ratios: From 3c4197f9761efd5f4dc67290106fdb885c7eec95 Mon Sep 17 00:00:00 2001 From: Jesus Pedro Marquez Hernandez Date: Thu, 23 Mar 2023 17:09:23 +0100 Subject: [PATCH 10/13] Using Number of Particles from dEdx --- include/algoSigProb.h | 2 + src/FlavorTag.cc | 205 ++++++++++++++++++++++++++++++------------ src/algoSigProb.cc | 49 ++++++++++ 3 files changed, 200 insertions(+), 56 deletions(-) diff --git a/include/algoSigProb.h b/include/algoSigProb.h index 4f9dab65..038d88f7 100644 --- a/include/algoSigProb.h +++ b/include/algoSigProb.h @@ -28,6 +28,8 @@ 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 dEdxKDSRatioPri(const Vertex* pri, string P1overP2, float GausWidth, float MaxMom, float MaxAngle); extern double dEdxKDSRatioSec(const VertexVec& vtxList, string P1overP2, float GausWidth, float MaxMom, float MaxAngle); } diff --git a/src/FlavorTag.cc b/src/FlavorTag.cc index 3083b1ef..39f1571a 100644 --- a/src/FlavorTag.cc +++ b/src/FlavorTag.cc @@ -1734,83 +1734,169 @@ class FtNBNess : public FTAlgo { } }; - class dEdxRatioPionOverKaonPri : public FTAlgo { - public: - dEdxRatioPionOverKaonPri() : FTAlgo("dEdxRatioPionOverKaonPri") {} - void process() { + //dEdx +class dEdxNKaonPri : public FTAlgo { + public: + dEdxNKaonPri() : FTAlgo("dEdxNKaonPri") {} + void process() { + _result = -2; + const VertexVec& vtxList = _jet->getVertices(); + double gaus = FTManager::getInstance().getKDSGausWidth(); + double mom = FTManager::getInstance().getKDSMinMom(); + double ang = FTManager::getInstance().getKDSMaxAng(); + _result=dEdxNPartPri( _privtx, string("Kaon"), gaus, mom, ang); + } +}; + +class dEdxNPionPri : public FTAlgo { + public: + dEdxNPionPri() : FTAlgo("dEdxNPionPri") {} + void process() { + _result = -2; + const VertexVec& vtxList = _jet->getVertices(); + double gaus = FTManager::getInstance().getKDSGausWidth(); + double mom = FTManager::getInstance().getKDSMinMom(); + double ang = FTManager::getInstance().getKDSMaxAng(); + _result=dEdxNPartPri( _privtx, string("Pion"), gaus, mom, ang); + } +}; + +class dEdxNProtonPri : public FTAlgo { + public: + dEdxNProtonPri() : FTAlgo("dEdxNProtonPri") {} + void process() { + _result = -2; + const VertexVec& vtxList = _jet->getVertices(); + double gaus = FTManager::getInstance().getKDSGausWidth(); + double mom = FTManager::getInstance().getKDSMinMom(); + double ang = FTManager::getInstance().getKDSMaxAng(); + _result=dEdxNPartPri( _privtx, string("Proton"), gaus, mom, ang); + } +}; + +class dEdxNKaonSec : public FTAlgo { + public: + dEdxNKaonSec() : FTAlgo("dEdxNKaonSec") {} + void process() { + _result = -2; + const VertexVec& vtxList = _jet->getVertices(); + if(vtxList.size()>0){ double gaus = FTManager::getInstance().getKDSGausWidth(); double mom = FTManager::getInstance().getKDSMinMom(); double ang = FTManager::getInstance().getKDSMaxAng(); - _result=dEdxKDSRatioPri( _privtx, string("PionOverKaon"), gaus, mom, ang); - } - }; - - class dEdxRatioKaonOverProtonPri : public FTAlgo { - public: - dEdxRatioKaonOverProtonPri() : FTAlgo("dEdxRatioKaonOverProtonPri") {} - void process() { + _result=dEdxNPartSec( vtxList, string("Kaon"), gaus, mom, ang); + } + } +}; + +class dEdxNPionSec : public FTAlgo { +public: + dEdxNPionSec() : FTAlgo("dEdxNPionSec") {} + void process() { + _result = -2; + const VertexVec& vtxList = _jet->getVertices(); + if(vtxList.size()>0){ double gaus = FTManager::getInstance().getKDSGausWidth(); double mom = FTManager::getInstance().getKDSMinMom(); double ang = FTManager::getInstance().getKDSMaxAng(); - _result=dEdxKDSRatioPri( _privtx, string("KaonOverProton"), gaus, mom, ang); + _result=dEdxNPartSec( vtxList, string("Pion"), gaus, mom, ang); } - }; - - class dEdxRatioPionOverProtonPri : public FTAlgo { - public: - dEdxRatioPionOverProtonPri() : FTAlgo("dEdxRatioPionOverProtonPri") {} - void process() { + } +}; + +class dEdxNProtonSec : public FTAlgo { +public: + dEdxNProtonSec() : FTAlgo("dEdxNProtonSec") {} + void process() { + _result = -2; + const VertexVec& vtxList = _jet->getVertices(); + if(vtxList.size()>0){ double gaus = FTManager::getInstance().getKDSGausWidth(); double mom = FTManager::getInstance().getKDSMinMom(); double ang = FTManager::getInstance().getKDSMaxAng(); - _result=dEdxKDSRatioPri( _privtx, string("PionOverProton"), gaus, mom, ang); + _result=dEdxNPartSec( vtxList, string("Proton"), gaus, mom, ang); } - }; - - 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().getKDSGausWidth(); - double mom = FTManager::getInstance().getKDSMinMom(); - double ang = FTManager::getInstance().getKDSMaxAng(); - _result=dEdxKDSRatioSec( vtxList, string("PionOverKaon"), gaus, mom, ang); - } + } +}; + + //dEdx Ratios +class dEdxRatioPionOverKaonPri : public FTAlgo { + public: + dEdxRatioPionOverKaonPri() : FTAlgo("dEdxRatioPionOverKaonPri") {} + void process() { + double gaus = FTManager::getInstance().getKDSGausWidth(); + double mom = FTManager::getInstance().getKDSMinMom(); + double ang = FTManager::getInstance().getKDSMaxAng(); + _result=dEdxKDSRatioPri( _privtx, string("PionOverKaon"), gaus, mom, ang); + } +}; + +class dEdxRatioKaonOverProtonPri : public FTAlgo { + public: + dEdxRatioKaonOverProtonPri() : FTAlgo("dEdxRatioKaonOverProtonPri") {} + void process() { + double gaus = FTManager::getInstance().getKDSGausWidth(); + double mom = FTManager::getInstance().getKDSMinMom(); + double ang = FTManager::getInstance().getKDSMaxAng(); + _result=dEdxKDSRatioPri( _privtx, string("KaonOverProton"), gaus, mom, ang); + } +}; + +class dEdxRatioPionOverProtonPri : public FTAlgo { + public: + dEdxRatioPionOverProtonPri() : FTAlgo("dEdxRatioPionOverProtonPri") {} + void process() { + double gaus = FTManager::getInstance().getKDSGausWidth(); + double mom = FTManager::getInstance().getKDSMinMom(); + double ang = FTManager::getInstance().getKDSMaxAng(); + _result=dEdxKDSRatioPri( _privtx, string("PionOverProton"), gaus, mom, ang); } - }; +}; - class dEdxRatioKaonOverProtonSec : public FTAlgo { - public: - dEdxRatioKaonOverProtonSec() : FTAlgo("dEdxRatioKaonOverProtonSec") {} - void process() { - _result = -2; - const VertexVec& vtxList = _jet->getVertices(); +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().getKDSGausWidth(); + double mom = FTManager::getInstance().getKDSMinMom(); + double ang = FTManager::getInstance().getKDSMaxAng(); + _result=dEdxKDSRatioSec( 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().getKDSGausWidth(); double mom = FTManager::getInstance().getKDSMinMom(); double ang = FTManager::getInstance().getKDSMaxAng(); _result=dEdxKDSRatioSec( 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().getKDSGausWidth(); - double mom = FTManager::getInstance().getKDSMinMom(); - double ang = FTManager::getInstance().getKDSMaxAng(); - _result=dEdxKDSRatioSec( vtxList, string("PionOverProton"), 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().getKDSGausWidth(); + double mom = FTManager::getInstance().getKDSMinMom(); + double ang = FTManager::getInstance().getKDSMaxAng(); + _result=dEdxKDSRatioSec( vtxList, string("PionOverProton"), gaus, mom, ang); } - }; + } +}; void FTManager::initVars() { if (_initialized)return; @@ -1956,6 +2042,13 @@ void FTManager::initVars() { add( new FtNBNess() ); //dEdx Variables + add( new dEdxNKaonPri() ); + add( new dEdxNPionPri() ); + add( new dEdxNProtonPri() ); + add( new dEdxNKaonSec() ); + add( new dEdxNPionSec() ); + add( new dEdxNProtonSec() ); + add( new dEdxRatioPionOverKaonSec() ); add( new dEdxRatioPionOverKaonSec() ); add( new dEdxRatioKaonOverProtonSec() ); add( new dEdxRatioPionOverProtonSec() ); diff --git a/src/algoSigProb.cc b/src/algoSigProb.cc index 568a4333..5aefeb60 100644 --- a/src/algoSigProb.cc +++ b/src/algoSigProb.cc @@ -497,6 +497,55 @@ void findMostSignificantTrack(const Jet* jet, const Vertex* pri, int minhitcut, } // Functions to fill observables using the Kaon dEdx istance in PID processors // Jesus Marquez Hernandez (IFIC-CSIC/UV) + double dEdxNPartPri(const Vertex* pri, string particle, float GausWidth, float MinMom, float MaxAngle) { + double counter=0; + TrackVec& vtxTracks = pri->getTracks(); // We load the primary vtx + try{ + for(unsigned int i=0;igetParticleIDProbability(PID_name); + 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(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; + } + + 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(PID_name); + 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(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; + } + + //Ratios double dEdxKDSRatioPri(const Vertex* pri, string P1overP2, float GausWidth, float MinMom, float MaxAngle) { double ratio=-3; // If a -3 appear, something went wrong From 56899192e42faee63318155b84855cd8f7b76f2d Mon Sep 17 00:00:00 2001 From: Jesus Pedro Marquez Hernandez Date: Mon, 3 Apr 2023 13:05:06 +0200 Subject: [PATCH 11/13] small changes for initialization issues --- src/FlavorTag.cc | 27 ++++++++++++--------------- src/algoSigProb.cc | 20 ++++++++++++++------ 2 files changed, 26 insertions(+), 21 deletions(-) diff --git a/src/FlavorTag.cc b/src/FlavorTag.cc index 39f1571a..383c8a4d 100644 --- a/src/FlavorTag.cc +++ b/src/FlavorTag.cc @@ -1739,12 +1739,11 @@ class dEdxNKaonPri : public FTAlgo { public: dEdxNKaonPri() : FTAlgo("dEdxNKaonPri") {} void process() { - _result = -2; - const VertexVec& vtxList = _jet->getVertices(); + _result = 0; double gaus = FTManager::getInstance().getKDSGausWidth(); double mom = FTManager::getInstance().getKDSMinMom(); double ang = FTManager::getInstance().getKDSMaxAng(); - _result=dEdxNPartPri( _privtx, string("Kaon"), gaus, mom, ang); + _result=dEdxNPartPri( _privtx, string("kaon"), gaus, mom, ang); } }; @@ -1752,12 +1751,11 @@ class dEdxNPionPri : public FTAlgo { public: dEdxNPionPri() : FTAlgo("dEdxNPionPri") {} void process() { - _result = -2; - const VertexVec& vtxList = _jet->getVertices(); + _result = 0; double gaus = FTManager::getInstance().getKDSGausWidth(); double mom = FTManager::getInstance().getKDSMinMom(); double ang = FTManager::getInstance().getKDSMaxAng(); - _result=dEdxNPartPri( _privtx, string("Pion"), gaus, mom, ang); + _result=dEdxNPartPri( _privtx, string("pion"), gaus, mom, ang); } }; @@ -1765,12 +1763,11 @@ class dEdxNProtonPri : public FTAlgo { public: dEdxNProtonPri() : FTAlgo("dEdxNProtonPri") {} void process() { - _result = -2; - const VertexVec& vtxList = _jet->getVertices(); + _result = 0; double gaus = FTManager::getInstance().getKDSGausWidth(); double mom = FTManager::getInstance().getKDSMinMom(); double ang = FTManager::getInstance().getKDSMaxAng(); - _result=dEdxNPartPri( _privtx, string("Proton"), gaus, mom, ang); + _result=dEdxNPartPri( _privtx, string("proton"), gaus, mom, ang); } }; @@ -1778,13 +1775,13 @@ class dEdxNKaonSec : public FTAlgo { public: dEdxNKaonSec() : FTAlgo("dEdxNKaonSec") {} void process() { - _result = -2; + _result = 0; const VertexVec& vtxList = _jet->getVertices(); if(vtxList.size()>0){ double gaus = FTManager::getInstance().getKDSGausWidth(); double mom = FTManager::getInstance().getKDSMinMom(); double ang = FTManager::getInstance().getKDSMaxAng(); - _result=dEdxNPartSec( vtxList, string("Kaon"), gaus, mom, ang); + _result=dEdxNPartSec( vtxList, string("kaon"), gaus, mom, ang); } } }; @@ -1793,13 +1790,13 @@ class dEdxNPionSec : public FTAlgo { public: dEdxNPionSec() : FTAlgo("dEdxNPionSec") {} void process() { - _result = -2; + _result = 0; const VertexVec& vtxList = _jet->getVertices(); if(vtxList.size()>0){ double gaus = FTManager::getInstance().getKDSGausWidth(); double mom = FTManager::getInstance().getKDSMinMom(); double ang = FTManager::getInstance().getKDSMaxAng(); - _result=dEdxNPartSec( vtxList, string("Pion"), gaus, mom, ang); + _result=dEdxNPartSec( vtxList, string("pion"), gaus, mom, ang); } } }; @@ -1808,13 +1805,13 @@ class dEdxNProtonSec : public FTAlgo { public: dEdxNProtonSec() : FTAlgo("dEdxNProtonSec") {} void process() { - _result = -2; + _result = 0; const VertexVec& vtxList = _jet->getVertices(); if(vtxList.size()>0){ double gaus = FTManager::getInstance().getKDSGausWidth(); double mom = FTManager::getInstance().getKDSMinMom(); double ang = FTManager::getInstance().getKDSMaxAng(); - _result=dEdxNPartSec( vtxList, string("Proton"), gaus, mom, ang); + _result=dEdxNPartSec( vtxList, string("proton"), gaus, mom, ang); } } }; diff --git a/src/algoSigProb.cc b/src/algoSigProb.cc index 5aefeb60..8fdd1522 100644 --- a/src/algoSigProb.cc +++ b/src/algoSigProb.cc @@ -502,12 +502,15 @@ void findMostSignificantTrack(const Jet* jet, const Vertex* pri, int minhitcut, TrackVec& vtxTracks = pri->getTracks(); // We load the primary vtx try{ for(unsigned int i=0;igetParticleIDProbability(PID_name); + double distance=0; + if(particle=="kaon")distance=vtxTracks.at(i)->getParticleIDProbability("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(KDS == 0) continue; // Initialization issue + if(distance == 0) continue; // Initialization issue + if(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 @@ -526,12 +529,15 @@ void findMostSignificantTrack(const Jet* jet, const Vertex* pri, int minhitcut, for (unsigned int j=0; jgetTracks(); for(unsigned int i=0;igetParticleIDProbability(PID_name); + double distance=0; + if(particle=="kaon")distance=vtxTracks.at(i)->getParticleIDProbability("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(KDS == 0) continue; // Initialization issue + if(distance == 0) continue; // Initialization issue + if(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 @@ -560,6 +566,7 @@ void findMostSignificantTrack(const Jet* jet, const Vertex* pri, int minhitcut, double costheta=vtxTracks.at(i)->CosTheta(); bool isMultiTrack=vtxTracks.at(i)->isMultiTrack(); if(KDS == 0) continue; // Initialization issue + if(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 @@ -603,6 +610,7 @@ void findMostSignificantTrack(const Jet* jet, const Vertex* pri, int minhitcut, double costheta=vtxTracks.at(i)->CosTheta(); bool isMultiTrack=vtxTracks.at(i)->isMultiTrack(); if(KDS == 0) continue; // Initialization issue + if(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 From 813ec455a403feb5d5df82d33ba088d811d508c5 Mon Sep 17 00:00:00 2001 From: Jesus Pedro Marquez Hernandez Date: Thu, 6 Apr 2023 18:52:18 +0200 Subject: [PATCH 12/13] Minor changes --- include/FlavorTag.h | 2 +- include/algoSigProb.h | 4 +- include/flavtag.h | 16 ++++---- src/FlavorTag.cc | 94 +++++++++++++++++++++---------------------- src/algoSigProb.cc | 12 +++--- src/flavtag.cc | 8 ++-- 6 files changed, 68 insertions(+), 68 deletions(-) diff --git a/include/FlavorTag.h b/include/FlavorTag.h index c1122af7..4be54b35 100644 --- a/include/FlavorTag.h +++ b/include/FlavorTag.h @@ -47,7 +47,7 @@ class FlavorTag : public Algorithm { int _nhitsJointProbZ0; int _nhitsMostSignificantTrack; - double _kdsGausWidth, _kdsMinMom, _kdsMaxAng; + double _dEdxGausWidth, _dEdxMinMom, _dEdxMaxAng; ClassDef(FlavorTag,1); }; diff --git a/include/algoSigProb.h b/include/algoSigProb.h index 038d88f7..6a30d85b 100644 --- a/include/algoSigProb.h +++ b/include/algoSigProb.h @@ -30,8 +30,8 @@ extern double jointProb2Z0(const Jet* jet, const Vertex* pri, int minhitcut, dou 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 dEdxKDSRatioPri(const Vertex* pri, string P1overP2, float GausWidth, float MaxMom, float MaxAngle); -extern double dEdxKDSRatioSec(const VertexVec& vtxList, string P1overP2, 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 12d77995..98e02a2c 100644 --- a/include/flavtag.h +++ b/include/flavtag.h @@ -89,7 +89,7 @@ class FTManager { void closeFile(); void process(const Event* event, const Vertex* privtx, int nhitsJointProbD0, int nhitsJointProbZ0, int nhitsMostSignificantTrack, JetVec& jets); //dEdx-KDS - void setKDSParameters(const double gaus, const double mom, const double ang); + void setDEDXParameters(const double gaus, const double mom, const double ang); float* getVarAddress(const string& varname); void setEval(bool seteval, bool exportAllVars) { @@ -118,16 +118,16 @@ class FTManager { _aux = aux; } - double getKDSGausWidth()const { - return _kdsGausWidth; + double getDEDXGausWidth()const { + return _dEdxGausWidth; } - double getKDSMinMom()const { - return _kdsMinMom; + double getDEDXMinMom()const { + return _dEdxMinMom; } - double getKDSMaxAng()const { - return _kdsMaxAng; + double getDEDXMaxAng()const { + return _dEdxMaxAng; } private: @@ -152,7 +152,7 @@ class FTManager { FtIPProbHolder* _holder; double _aux; // dEdx-KDS - double _kdsGausWidth, _kdsMinMom, _kdsMaxAng; + double _dEdxGausWidth, _dEdxMinMom, _dEdxMaxAng; }; // histogram holder for d0/z0 probability diff --git a/src/FlavorTag.cc b/src/FlavorTag.cc index 383c8a4d..effb1506 100644 --- a/src/FlavorTag.cc +++ b/src/FlavorTag.cc @@ -1740,9 +1740,9 @@ class dEdxNKaonPri : public FTAlgo { dEdxNKaonPri() : FTAlgo("dEdxNKaonPri") {} void process() { _result = 0; - double gaus = FTManager::getInstance().getKDSGausWidth(); - double mom = FTManager::getInstance().getKDSMinMom(); - double ang = FTManager::getInstance().getKDSMaxAng(); + double gaus = FTManager::getInstance().getDEDXGausWidth(); + double mom = FTManager::getInstance().getDEDXMinMom(); + double ang = FTManager::getInstance().getDEDXMaxAng(); _result=dEdxNPartPri( _privtx, string("kaon"), gaus, mom, ang); } }; @@ -1752,9 +1752,9 @@ class dEdxNPionPri : public FTAlgo { dEdxNPionPri() : FTAlgo("dEdxNPionPri") {} void process() { _result = 0; - double gaus = FTManager::getInstance().getKDSGausWidth(); - double mom = FTManager::getInstance().getKDSMinMom(); - double ang = FTManager::getInstance().getKDSMaxAng(); + double gaus = FTManager::getInstance().getDEDXGausWidth(); + double mom = FTManager::getInstance().getDEDXMinMom(); + double ang = FTManager::getInstance().getDEDXMaxAng(); _result=dEdxNPartPri( _privtx, string("pion"), gaus, mom, ang); } }; @@ -1764,9 +1764,9 @@ class dEdxNProtonPri : public FTAlgo { dEdxNProtonPri() : FTAlgo("dEdxNProtonPri") {} void process() { _result = 0; - double gaus = FTManager::getInstance().getKDSGausWidth(); - double mom = FTManager::getInstance().getKDSMinMom(); - double ang = FTManager::getInstance().getKDSMaxAng(); + double gaus = FTManager::getInstance().getDEDXGausWidth(); + double mom = FTManager::getInstance().getDEDXMinMom(); + double ang = FTManager::getInstance().getDEDXMaxAng(); _result=dEdxNPartPri( _privtx, string("proton"), gaus, mom, ang); } }; @@ -1778,9 +1778,9 @@ class dEdxNKaonSec : public FTAlgo { _result = 0; const VertexVec& vtxList = _jet->getVertices(); if(vtxList.size()>0){ - double gaus = FTManager::getInstance().getKDSGausWidth(); - double mom = FTManager::getInstance().getKDSMinMom(); - double ang = FTManager::getInstance().getKDSMaxAng(); + double gaus = FTManager::getInstance().getDEDXGausWidth(); + double mom = FTManager::getInstance().getDEDXMinMom(); + double ang = FTManager::getInstance().getDEDXMaxAng(); _result=dEdxNPartSec( vtxList, string("kaon"), gaus, mom, ang); } } @@ -1793,9 +1793,9 @@ class dEdxNPionSec : public FTAlgo { _result = 0; const VertexVec& vtxList = _jet->getVertices(); if(vtxList.size()>0){ - double gaus = FTManager::getInstance().getKDSGausWidth(); - double mom = FTManager::getInstance().getKDSMinMom(); - double ang = FTManager::getInstance().getKDSMaxAng(); + double gaus = FTManager::getInstance().getDEDXGausWidth(); + double mom = FTManager::getInstance().getDEDXMinMom(); + double ang = FTManager::getInstance().getDEDXMaxAng(); _result=dEdxNPartSec( vtxList, string("pion"), gaus, mom, ang); } } @@ -1808,9 +1808,9 @@ class dEdxNProtonSec : public FTAlgo { _result = 0; const VertexVec& vtxList = _jet->getVertices(); if(vtxList.size()>0){ - double gaus = FTManager::getInstance().getKDSGausWidth(); - double mom = FTManager::getInstance().getKDSMinMom(); - double ang = FTManager::getInstance().getKDSMaxAng(); + double gaus = FTManager::getInstance().getDEDXGausWidth(); + double mom = FTManager::getInstance().getDEDXMinMom(); + double ang = FTManager::getInstance().getDEDXMaxAng(); _result=dEdxNPartSec( vtxList, string("proton"), gaus, mom, ang); } } @@ -1821,10 +1821,10 @@ class dEdxRatioPionOverKaonPri : public FTAlgo { public: dEdxRatioPionOverKaonPri() : FTAlgo("dEdxRatioPionOverKaonPri") {} void process() { - double gaus = FTManager::getInstance().getKDSGausWidth(); - double mom = FTManager::getInstance().getKDSMinMom(); - double ang = FTManager::getInstance().getKDSMaxAng(); - _result=dEdxKDSRatioPri( _privtx, string("PionOverKaon"), gaus, mom, ang); + double gaus = FTManager::getInstance().getDEDXGausWidth(); + double mom = FTManager::getInstance().getDEDXMinMom(); + double ang = FTManager::getInstance().getDEDXMaxAng(); + _result=dEdxRatioPri( _privtx, string("PionOverKaon"), gaus, mom, ang); } }; @@ -1832,10 +1832,10 @@ class dEdxRatioKaonOverProtonPri : public FTAlgo { public: dEdxRatioKaonOverProtonPri() : FTAlgo("dEdxRatioKaonOverProtonPri") {} void process() { - double gaus = FTManager::getInstance().getKDSGausWidth(); - double mom = FTManager::getInstance().getKDSMinMom(); - double ang = FTManager::getInstance().getKDSMaxAng(); - _result=dEdxKDSRatioPri( _privtx, string("KaonOverProton"), gaus, mom, ang); + double gaus = FTManager::getInstance().getDEDXGausWidth(); + double mom = FTManager::getInstance().getDEDXMinMom(); + double ang = FTManager::getInstance().getDEDXMaxAng(); + _result=dEdxRatioPri( _privtx, string("KaonOverProton"), gaus, mom, ang); } }; @@ -1843,10 +1843,10 @@ class dEdxRatioPionOverProtonPri : public FTAlgo { public: dEdxRatioPionOverProtonPri() : FTAlgo("dEdxRatioPionOverProtonPri") {} void process() { - double gaus = FTManager::getInstance().getKDSGausWidth(); - double mom = FTManager::getInstance().getKDSMinMom(); - double ang = FTManager::getInstance().getKDSMaxAng(); - _result=dEdxKDSRatioPri( _privtx, string("PionOverProton"), gaus, mom, ang); + double gaus = FTManager::getInstance().getDEDXGausWidth(); + double mom = FTManager::getInstance().getDEDXMinMom(); + double ang = FTManager::getInstance().getDEDXMaxAng(); + _result=dEdxRatioPri( _privtx, string("PionOverProton"), gaus, mom, ang); } }; @@ -1857,10 +1857,10 @@ class dEdxRatioPionOverKaonSec : public FTAlgo { _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().getKDSGausWidth(); - double mom = FTManager::getInstance().getKDSMinMom(); - double ang = FTManager::getInstance().getKDSMaxAng(); - _result=dEdxKDSRatioSec( vtxList, string("PionOverKaon"), gaus, mom, ang); + double gaus = FTManager::getInstance().getDEDXGausWidth(); + double mom = FTManager::getInstance().getDEDXMinMom(); + double ang = FTManager::getInstance().getDEDXMaxAng(); + _result=dEdxRatioSec( vtxList, string("PionOverKaon"), gaus, mom, ang); } } }; @@ -1872,10 +1872,10 @@ class dEdxRatioKaonOverProtonSec : public FTAlgo { _result = -2; const VertexVec& vtxList = _jet->getVertices(); if(vtxList.size()>0){ - double gaus = FTManager::getInstance().getKDSGausWidth(); - double mom = FTManager::getInstance().getKDSMinMom(); - double ang = FTManager::getInstance().getKDSMaxAng(); - _result=dEdxKDSRatioSec( vtxList, string("KaonOverProton"), gaus, mom, ang); + double gaus = FTManager::getInstance().getDEDXGausWidth(); + double mom = FTManager::getInstance().getDEDXMinMom(); + double ang = FTManager::getInstance().getDEDXMaxAng(); + _result=dEdxRatioSec( vtxList, string("KaonOverProton"), gaus, mom, ang); } } }; @@ -1887,10 +1887,10 @@ class dEdxRatioPionOverProtonSec : public FTAlgo { _result = -2; const VertexVec& vtxList = _jet->getVertices(); if(vtxList.size()>0){ - double gaus = FTManager::getInstance().getKDSGausWidth(); - double mom = FTManager::getInstance().getKDSMinMom(); - double ang = FTManager::getInstance().getKDSMaxAng(); - _result=dEdxKDSRatioSec( vtxList, string("PionOverProton"), gaus, mom, ang); + double gaus = FTManager::getInstance().getDEDXGausWidth(); + double mom = FTManager::getInstance().getDEDXMinMom(); + double ang = FTManager::getInstance().getDEDXMaxAng(); + _result=dEdxRatioSec( vtxList, string("PionOverProton"), gaus, mom, ang); } } }; @@ -2071,10 +2071,10 @@ 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-KDS - _kdsGausWidth = param->get("FlavorTag.KDSGaussianWidth",double(4)); - _kdsMinMom = param->get("FlavorTag.KDSMinMomentum",double(3)); - _kdsMaxAng = param->get("FlavorTag.KDSMaxAngle",double(0.95)); + //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)); @@ -2109,7 +2109,7 @@ void FlavorTag::process() { mgr.setIPProbHolder(_holder); mgr.process(event, privtx, _nhitsJointProbD0, _nhitsJointProbZ0, _nhitsMostSignificantTrack, jets); - mgr.setKDSParameters(_kdsGausWidth, _kdsMinMom, _kdsMaxAng); + mgr.setDEDXParameters(_dEdxGausWidth, _dEdxMinMom, _dEdxMaxAng); } void FlavorTag::end() { diff --git a/src/algoSigProb.cc b/src/algoSigProb.cc index 8fdd1522..98b42843 100644 --- a/src/algoSigProb.cc +++ b/src/algoSigProb.cc @@ -510,7 +510,7 @@ void findMostSignificantTrack(const Jet* jet, const Vertex* pri, int minhitcut, double costheta=vtxTracks.at(i)->CosTheta(); bool isMultiTrack=vtxTracks.at(i)->isMultiTrack(); if(distance == 0) continue; // Initialization issue - if(distance > 100) continue; // Initialization issue if rewritten + 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 @@ -537,7 +537,7 @@ void findMostSignificantTrack(const Jet* jet, const Vertex* pri, int minhitcut, double costheta=vtxTracks.at(i)->CosTheta(); bool isMultiTrack=vtxTracks.at(i)->isMultiTrack(); if(distance == 0) continue; // Initialization issue - if(distance > 100) continue; // Initialization issue if rewritten + 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 @@ -552,7 +552,7 @@ void findMostSignificantTrack(const Jet* jet, const Vertex* pri, int minhitcut, } //Ratios - double dEdxKDSRatioPri(const Vertex* pri, string P1overP2, float GausWidth, float MinMom, float MaxAngle) { + double dEdxRatioPri(const Vertex* pri, string P1overP2, float GausWidth, float MinMom, float MaxAngle) { double ratio=-3; // If a -3 appear, something went wrong double neg_counter=0; // Estimated protons @@ -566,7 +566,7 @@ void findMostSignificantTrack(const Jet* jet, const Vertex* pri, int minhitcut, double costheta=vtxTracks.at(i)->CosTheta(); bool isMultiTrack=vtxTracks.at(i)->isMultiTrack(); if(KDS == 0) continue; // Initialization issue - if(KDS > 100) continue; // Initialization issue if rewritten + 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 @@ -595,7 +595,7 @@ void findMostSignificantTrack(const Jet* jet, const Vertex* pri, int minhitcut, } - double dEdxKDSRatioSec(const VertexVec& vtxList, string P1overP2, float GausWidth, float MinMom, float MaxAngle) { + double dEdxRatioSec(const VertexVec& vtxList, string P1overP2, float GausWidth, float MinMom, float MaxAngle) { double ratio=-3; double neg_counter=0; //estimated protons @@ -610,7 +610,7 @@ void findMostSignificantTrack(const Jet* jet, const Vertex* pri, int minhitcut, double costheta=vtxTracks.at(i)->CosTheta(); bool isMultiTrack=vtxTracks.at(i)->isMultiTrack(); if(KDS == 0) continue; // Initialization issue - if(KDS > 100) continue; // Initialization issue if rewritten + 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 diff --git a/src/flavtag.cc b/src/flavtag.cc index f3c081cf..bd8ad53a 100644 --- a/src/flavtag.cc +++ b/src/flavtag.cc @@ -166,10 +166,10 @@ void FTManager::process(const Event* event, const Vertex* privtx, int nhitsJoint } //KDS Parameters -void FTManager::setKDSParameters(const double gaus, const double mom, const double ang) { - _kdsGausWidth = gaus; - _kdsMinMom = mom; - _kdsMaxAng = ang; +void FTManager::setDEDXParameters(const double gaus, const double mom, const double ang) { + _dEdxGausWidth = gaus; + _dEdxMinMom = mom; + _dEdxMaxAng = ang; } float* FTManager::getVarAddress(const string& varname) { From 294931e48f7f8c3fed48d146c33c0148209a49ea Mon Sep 17 00:00:00 2001 From: Jesus Pedro Marquez Hernandez Date: Tue, 6 Jun 2023 17:24:16 +0200 Subject: [PATCH 13/13] Reducing changes to the minimum for an official implementation. --- src/FlavorTag.cc | 76 ---------------------------------------------- src/algoSigProb.cc | 74 +++----------------------------------------- 2 files changed, 4 insertions(+), 146 deletions(-) diff --git a/src/FlavorTag.cc b/src/FlavorTag.cc index effb1506..fc3beb78 100644 --- a/src/FlavorTag.cc +++ b/src/FlavorTag.cc @@ -1735,42 +1735,6 @@ class FtNBNess : public FTAlgo { }; //dEdx -class dEdxNKaonPri : public FTAlgo { - public: - dEdxNKaonPri() : FTAlgo("dEdxNKaonPri") {} - void process() { - _result = 0; - double gaus = FTManager::getInstance().getDEDXGausWidth(); - double mom = FTManager::getInstance().getDEDXMinMom(); - double ang = FTManager::getInstance().getDEDXMaxAng(); - _result=dEdxNPartPri( _privtx, string("kaon"), gaus, mom, ang); - } -}; - -class dEdxNPionPri : public FTAlgo { - public: - dEdxNPionPri() : FTAlgo("dEdxNPionPri") {} - void process() { - _result = 0; - double gaus = FTManager::getInstance().getDEDXGausWidth(); - double mom = FTManager::getInstance().getDEDXMinMom(); - double ang = FTManager::getInstance().getDEDXMaxAng(); - _result=dEdxNPartPri( _privtx, string("pion"), gaus, mom, ang); - } -}; - -class dEdxNProtonPri : public FTAlgo { - public: - dEdxNProtonPri() : FTAlgo("dEdxNProtonPri") {} - void process() { - _result = 0; - double gaus = FTManager::getInstance().getDEDXGausWidth(); - double mom = FTManager::getInstance().getDEDXMinMom(); - double ang = FTManager::getInstance().getDEDXMaxAng(); - _result=dEdxNPartPri( _privtx, string("proton"), gaus, mom, ang); - } -}; - class dEdxNKaonSec : public FTAlgo { public: dEdxNKaonSec() : FTAlgo("dEdxNKaonSec") {} @@ -1817,39 +1781,6 @@ class dEdxNProtonSec : public FTAlgo { }; //dEdx Ratios -class dEdxRatioPionOverKaonPri : public FTAlgo { - public: - dEdxRatioPionOverKaonPri() : FTAlgo("dEdxRatioPionOverKaonPri") {} - void process() { - double gaus = FTManager::getInstance().getDEDXGausWidth(); - double mom = FTManager::getInstance().getDEDXMinMom(); - double ang = FTManager::getInstance().getDEDXMaxAng(); - _result=dEdxRatioPri( _privtx, string("PionOverKaon"), gaus, mom, ang); - } -}; - -class dEdxRatioKaonOverProtonPri : public FTAlgo { - public: - dEdxRatioKaonOverProtonPri() : FTAlgo("dEdxRatioKaonOverProtonPri") {} - void process() { - double gaus = FTManager::getInstance().getDEDXGausWidth(); - double mom = FTManager::getInstance().getDEDXMinMom(); - double ang = FTManager::getInstance().getDEDXMaxAng(); - _result=dEdxRatioPri( _privtx, string("KaonOverProton"), gaus, mom, ang); - } -}; - -class dEdxRatioPionOverProtonPri : public FTAlgo { - public: - dEdxRatioPionOverProtonPri() : FTAlgo("dEdxRatioPionOverProtonPri") {} - void process() { - double gaus = FTManager::getInstance().getDEDXGausWidth(); - double mom = FTManager::getInstance().getDEDXMinMom(); - double ang = FTManager::getInstance().getDEDXMaxAng(); - _result=dEdxRatioPri( _privtx, string("PionOverProton"), gaus, mom, ang); - } -}; - class dEdxRatioPionOverKaonSec : public FTAlgo { public: dEdxRatioPionOverKaonSec() : FTAlgo("dEdxRatioPionOverKaonSec") {} @@ -2006,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()); @@ -2039,9 +1969,6 @@ void FTManager::initVars() { add( new FtNBNess() ); //dEdx Variables - add( new dEdxNKaonPri() ); - add( new dEdxNPionPri() ); - add( new dEdxNProtonPri() ); add( new dEdxNKaonSec() ); add( new dEdxNPionSec() ); add( new dEdxNProtonSec() ); @@ -2049,9 +1976,6 @@ void FTManager::initVars() { add( new dEdxRatioPionOverKaonSec() ); add( new dEdxRatioKaonOverProtonSec() ); add( new dEdxRatioPionOverProtonSec() ); - add( new dEdxRatioPionOverKaonPri() ); - add( new dEdxRatioKaonOverProtonPri() ); - add( new dEdxRatioPionOverProtonPri() ); } void FlavorTag::init(Parameters* param) { diff --git a/src/algoSigProb.cc b/src/algoSigProb.cc index 98b42843..604a9710 100644 --- a/src/algoSigProb.cc +++ b/src/algoSigProb.cc @@ -495,33 +495,10 @@ void findMostSignificantTrack(const Jet* jet, const Vertex* pri, int minhitcut, sigVec[5] = 0; } } - // Functions to fill observables using the Kaon dEdx istance in PID processors + // 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 dEdxNPartPri(const Vertex* pri, string particle, float GausWidth, float MinMom, float MaxAngle) { - double counter=0; - TrackVec& vtxTracks = pri->getTracks(); // We load the primary vtx - try{ - 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; - } double dEdxNPartSec(const VertexVec& vtxList, string particle, float GausWidth, float MinMom, float MaxAngle) { double counter=0; @@ -551,52 +528,9 @@ void findMostSignificantTrack(const Jet* jet, const Vertex* pri, int minhitcut, return counter; } - //Ratios - double dEdxRatioPri(const Vertex* pri, string P1overP2, float GausWidth, float MinMom, float MaxAngle) { - - double ratio=-3; // If a -3 appear, something went wrong - double neg_counter=0; // Estimated protons - double null_counter=0; // Estimated kaons - double pos_counter=0; // Estimated pion - TrackVec& vtxTracks = pri->getTracks(); // We load the primary vtx - try{ - 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; - - } + // 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