diff --git a/PWGCF/DataModel/CorrelationsDerived.h b/PWGCF/DataModel/CorrelationsDerived.h index c6b96b77466..5f79a78e974 100644 --- a/PWGCF/DataModel/CorrelationsDerived.h +++ b/PWGCF/DataModel/CorrelationsDerived.h @@ -124,6 +124,10 @@ enum ParticleDecay { JPsiToEE, JPsiToMuMu, Generic2Prong, + PhiToKK, + K0stoPiPi, + LambdatoPPi, + AntiLambdatoPiP }; } // namespace cf2prongtrack DECLARE_SOA_TABLE(CF2ProngTracks, "AOD", "CF2PRONGTRACK", //! Reduced track table diff --git a/PWGCF/TableProducer/filter2Prong.cxx b/PWGCF/TableProducer/filter2Prong.cxx index 2ee72e9b923..105d4d25871 100644 --- a/PWGCF/TableProducer/filter2Prong.cxx +++ b/PWGCF/TableProducer/filter2Prong.cxx @@ -16,6 +16,8 @@ #include "PWGHF/DataModel/CandidateReconstructionTables.h" #include "PWGHF/DataModel/CandidateSelectionTables.h" +#include "Common/DataModel/PIDResponseITS.h" + #include "Framework/ASoAHelpers.h" #include "Framework/AnalysisDataModel.h" #include "Framework/AnalysisTask.h" @@ -34,6 +36,10 @@ using namespace o2::framework; using namespace o2::framework::expressions; using namespace o2::math_utils::detail; +enum LambdaPid { kLambda = 0, + kAntiLambda +}; + // #define FLOAT_PRECISION 0xFFFFFFF0 #define O2_DEFINE_CONFIGURABLE(NAME, TYPE, DEFAULT, HELP) Configurable NAME{#NAME, DEFAULT, HELP}; @@ -44,10 +50,56 @@ struct Filter2Prong { O2_DEFINE_CONFIGURABLE(cfgImPart2Mass, float, o2::constants::physics::MassKMinus, "Daughter particle 2 mass in GeV") O2_DEFINE_CONFIGURABLE(cfgImPart1PID, int, o2::track::PID::Kaon, "PID of daughter particle 1 (O2 PID ID)") O2_DEFINE_CONFIGURABLE(cfgImPart2PID, int, o2::track::PID::Kaon, "PID of daughter particle 2 (O2 PID ID)") + O2_DEFINE_CONFIGURABLE(cfgMomDepPID, bool, 1, "Use mommentum dependent PID for Phi meson") O2_DEFINE_CONFIGURABLE(cfgImCutPt, float, 0.2f, "Minimal pT for candidates") - O2_DEFINE_CONFIGURABLE(cfgImMinInvMass, float, 0.95f, "Minimum invariant mass (GeV)") - O2_DEFINE_CONFIGURABLE(cfgImMaxInvMass, float, 1.07f, "Maximum invariant mass (GeV)") - O2_DEFINE_CONFIGURABLE(cfgImSigmaFormula, std::string, "(z < 0.5 && x < 3.0) || (z >= 0.5 && x < 2.5 && y < 3.0)", "pT dependent daughter track sigma pass condition (x = TPC sigma, y = TOF sigma, z = pT)") + O2_DEFINE_CONFIGURABLE(cfgImMinInvMass, float, 0.95f, "Minimum invariant mass for generic 2 prong") + O2_DEFINE_CONFIGURABLE(cfgImMaxInvMass, float, 1.07f, "Maximum invariant mass for generic 2 prong") + O2_DEFINE_CONFIGURABLE(cfgImSigmaFormula, std::string, "(z < 0.5 && abs(x) < 3.0) || (z >= 0.5 && abs(x) < 2.5 && abs(y) < 3.0)", "pT dependent daughter track sigma pass condition (x = TPC sigma, y = TOF sigma, z = pT)") + + struct : ConfigurableGroup{ + O2_DEFINE_CONFIGURABLE(tpcNClsCrossedRowsTrackMin, float, 70, "Minimum number of crossed rows in TPC") + O2_DEFINE_CONFIGURABLE(etaTrackMax, float, 0.8, "Maximum pseudorapidity") + O2_DEFINE_CONFIGURABLE(ptTrackMin, float, 0.1, "Minimum transverse momentum") + O2_DEFINE_CONFIGURABLE(minV0DCAPr, float, 0.1, "Min V0 proton DCA") + O2_DEFINE_CONFIGURABLE(minV0DCAPiLambda, float, 0.1, "Min V0 pion DCA for lambda") + O2_DEFINE_CONFIGURABLE(minV0DCAPiK0s, float, 0.1, "Min V0 pion DCA for K0s") + O2_DEFINE_CONFIGURABLE(daughPIDCuts, float, 4.0, "PID nsigma for V0s") + O2_DEFINE_CONFIGURABLE(massK0Min, float, 0.4, "Minimum mass for K0") + O2_DEFINE_CONFIGURABLE(massK0Max, float, 0.6, "Maximum mass for K0") + O2_DEFINE_CONFIGURABLE(massLambdaMin, float, 1.0, "Minimum mass for lambda") + O2_DEFINE_CONFIGURABLE(massLambdaMax, float, 1.3, "Maximum mass for lambda") + O2_DEFINE_CONFIGURABLE(radiusMaxLambda, float, 2.3, "Maximum decay radius (cm) for lambda") + O2_DEFINE_CONFIGURABLE(radiusMinLambda, float, 0.0, "Minimum decay radius (cm) for lambda") + O2_DEFINE_CONFIGURABLE(radiusMaxK0s, float, 2.3, "Maximum decay radius (cm) for K0s") + O2_DEFINE_CONFIGURABLE(radiusMinK0s, float, 0.0, "Minimum decay radius (cm) for K0s") + O2_DEFINE_CONFIGURABLE(cosPaMinLambda, float, 0.98, "Minimum cosine of pointing angle for lambda") + O2_DEFINE_CONFIGURABLE(cosPaMinK0s, float, 0.98, "Minimum cosine of pointing angle for K0s") + O2_DEFINE_CONFIGURABLE(dcaV0DaughtersMaxLambda, float, 0.2, "Maximum DCA among the V0 daughters (cm) for lambda") + O2_DEFINE_CONFIGURABLE(dcaV0DaughtersMaxK0s, float, 0.2, "Maximum DCA among the V0 daughters (cm) for K0s") + O2_DEFINE_CONFIGURABLE(qtArmenterosMinForK0s, float, 0.12, "Minimum Armenteros' qt for K0s") + O2_DEFINE_CONFIGURABLE(maxLambdaLifeTime, float, 30, "Maximum lambda lifetime (in cm)") + O2_DEFINE_CONFIGURABLE(maxK0sLifeTime, float, 30, "Maximum K0s lifetime (in cm)")} grpV0; + + struct : ConfigurableGroup{ + O2_DEFINE_CONFIGURABLE(ImMinInvMassPhiMeson, float, 0.98f, "Minimum invariant mass Phi meson (GeV)") + O2_DEFINE_CONFIGURABLE(ImMaxInvMassPhiMeson, float, 1.07f, "Maximum invariant mass Phi meson (GeV)") + O2_DEFINE_CONFIGURABLE(ITSPIDSelection, bool, true, "PID ITS") + O2_DEFINE_CONFIGURABLE(ITSPIDPthreshold, float, 1.0, "Momentum threshold for ITS PID (GeV/c) (only used if ITSPIDSelection is true)") + O2_DEFINE_CONFIGURABLE(lowITSPIDNsigma, float, 3.0, "lower cut on PID nsigma for ITS") + O2_DEFINE_CONFIGURABLE(highITSPIDNsigma, float, 3.0, "higher cut on PID nsigma for ITS") + O2_DEFINE_CONFIGURABLE(ITSclusterPhiMeson, int, 5, "Minimum number of ITS cluster for phi meson track") + O2_DEFINE_CONFIGURABLE(TPCCrossedRowsPhiMeson, int, 80, "Minimum number of TPC Crossed Rows for phi meson track") + O2_DEFINE_CONFIGURABLE(cutDCAxyPhiMeson, float, 0.1, "Maximum DCAxy for phi meson track") + O2_DEFINE_CONFIGURABLE(cutDCAzPhiMeson, float, 0.1, "Maximum DCAz for phi meson track") + O2_DEFINE_CONFIGURABLE(cutEtaPhiMeson, float, 0.8, "Maximum eta for phi meson track") + O2_DEFINE_CONFIGURABLE(cutPTPhiMeson, float, 0.15, "Maximum pt for phi meson track") + O2_DEFINE_CONFIGURABLE(isDeepAngle, bool, true, "Flag for applying deep angle") + O2_DEFINE_CONFIGURABLE(deepAngle, float, 0.04, "Deep angle cut") + O2_DEFINE_CONFIGURABLE(nsigmaCutTPC, float, 2.5, "nsigma tpc") + O2_DEFINE_CONFIGURABLE(nsigmaCutTOF, float, 2.5, "nsigma tof") + O2_DEFINE_CONFIGURABLE(cutTOFBeta, float, 0.5, "TOF beta") + O2_DEFINE_CONFIGURABLE(confFakeKaonCut, float, 0.15, "Cut based on track from momentum difference") + O2_DEFINE_CONFIGURABLE(removefaketrack, bool, true, "Flag to remove fake kaon")} grpPhi; HfHelper hfHelper; Produces output2ProngTracks; @@ -64,6 +116,9 @@ struct Filter2Prong { template using HasMLProb = decltype(std::declval().mlProbD0()); + using PIDTrack = soa::Join; + using ResoV0s = aod::V0Datas; + std::unique_ptr sigmaFormula; void init(InitContext&) @@ -166,9 +221,166 @@ struct Filter2Prong { } PROCESS_SWITCH(Filter2Prong, processMC, "Process MC 2-prong daughters", false); + template + bool selectionTrack(const T& candidate) + { + if (candidate.isGlobalTrack() && candidate.isPVContributor() && candidate.itsNCls() >= grpPhi.ITSclusterPhiMeson && candidate.tpcNClsCrossedRows() > grpPhi.TPCCrossedRowsPhiMeson && std::abs(candidate.dcaXY()) <= grpPhi.cutDCAxyPhiMeson && std::abs(candidate.dcaZ()) <= grpPhi.cutDCAzPhiMeson && std::abs(candidate.eta()) <= grpPhi.cutEtaPhiMeson && candidate.pt() >= grpPhi.cutPTPhiMeson) { + return true; + } + return false; + } + + template + bool selectionPair(const T1& candidate1, const T2& candidate2) + { + double pt1, pt2, pz1, pz2, p1, p2, angle; + pt1 = candidate1.pt(); + pt2 = candidate2.pt(); + pz1 = candidate1.pz(); + pz2 = candidate2.pz(); + p1 = candidate1.p(); + p2 = candidate2.p(); + angle = TMath::ACos((pt1 * pt2 + pz1 * pz2) / (p1 * p2)); + if (grpPhi.isDeepAngle && angle < grpPhi.deepAngle) { + return false; + } + return true; + } + + template + bool isFakeTrack(T const& track) + { + const auto pglobal = track.p(); + const auto ptpc = track.tpcInnerParam(); + if (TMath::Abs(pglobal - ptpc) > grpPhi.confFakeKaonCut) { + return true; + } + return false; + } + + template + bool isSelectedV0AsK0s(Collision const& collision, const V0Cand& v0) + { + const auto& posTrack = v0.template posTrack_as(); + const auto& negTrack = v0.template negTrack_as(); + + float CtauK0s = v0.distovertotmom(collision.posX(), collision.posY(), collision.posZ()) * o2::constants::physics::MassK0; + + if (v0.mK0Short() < grpV0.massK0Min || v0.mK0Short() > grpV0.massK0Max) { + return false; + } + if ((v0.qtarm() / std::abs(v0.alpha())) < grpV0.qtArmenterosMinForK0s) { + return false; + } + if (v0.v0radius() > grpV0.radiusMaxK0s || v0.v0radius() < grpV0.radiusMinK0s) { + return false; + } + if (v0.v0cosPA() < grpV0.cosPaMinK0s) { + return false; + } + if (v0.dcaV0daughters() > grpV0.dcaV0DaughtersMaxK0s) { + return false; + } + if (std::abs(CtauK0s) > grpV0.maxK0sLifeTime) { + return false; + } + if (((std::abs(posTrack.tpcNSigmaPi()) > grpV0.daughPIDCuts) || (std::abs(negTrack.tpcNSigmaPi()) > grpV0.daughPIDCuts))) { + return false; + } + if ((TMath::Abs(v0.dcapostopv()) < grpV0.minV0DCAPiK0s || TMath::Abs(v0.dcanegtopv()) < grpV0.minV0DCAPiK0s)) { + return false; + } + return true; + } + + template + bool isSelectedV0AsLambda(Collision const& collision, const V0Cand& v0) + { + const auto& posTrack = v0.template posTrack_as(); + const auto& negTrack = v0.template negTrack_as(); + + float CtauLambda = v0.distovertotmom(collision.posX(), collision.posY(), collision.posZ()) * o2::constants::physics::MassLambda; + + if ((v0.mLambda() < grpV0.massLambdaMin || v0.mLambda() > grpV0.massLambdaMax) && + (v0.mAntiLambda() < grpV0.massLambdaMin || v0.mAntiLambda() > grpV0.massLambdaMax)) { + return false; + } + if (v0.v0radius() > grpV0.radiusMaxLambda || v0.v0radius() < grpV0.radiusMinLambda) { + return false; + } + if (v0.v0cosPA() < grpV0.cosPaMinLambda) { + return false; + } + if (v0.dcaV0daughters() > grpV0.dcaV0DaughtersMaxLambda) { + return false; + } + if (pid == LambdaPid::kLambda && (TMath::Abs(v0.dcapostopv()) < grpV0.minV0DCAPr || TMath::Abs(v0.dcanegtopv()) < grpV0.minV0DCAPiLambda)) { + return false; + } + if (pid == LambdaPid::kAntiLambda && (TMath::Abs(v0.dcapostopv()) < grpV0.minV0DCAPiLambda || TMath::Abs(v0.dcanegtopv()) < grpV0.minV0DCAPr)) { + return false; + } + if (pid == LambdaPid::kLambda && ((std::abs(posTrack.tpcNSigmaPr()) > grpV0.daughPIDCuts) || (std::abs(negTrack.tpcNSigmaPi()) > grpV0.daughPIDCuts))) { + return false; + } + if (pid == LambdaPid::kAntiLambda && ((std::abs(posTrack.tpcNSigmaPi()) > grpV0.daughPIDCuts) || (std::abs(negTrack.tpcNSigmaPr()) > grpV0.daughPIDCuts))) { + return false; + } + if (std::abs(CtauLambda) > grpV0.maxLambdaLifeTime) { + return false; + } + return true; + } + + template + bool isV0TrackSelected(const T1& v0) + { + const auto& posTrack = v0.template posTrack_as(); + const auto& negTrack = v0.template negTrack_as(); + + if (!posTrack.hasTPC() || !negTrack.hasTPC()) { + return false; + } + if (posTrack.tpcNClsCrossedRows() < grpV0.tpcNClsCrossedRowsTrackMin || negTrack.tpcNClsCrossedRows() < grpV0.tpcNClsCrossedRowsTrackMin) { + return false; + } + if (posTrack.tpcCrossedRowsOverFindableCls() < 0.8 || negTrack.tpcCrossedRowsOverFindableCls() < 0.8) { + return false; + } + if (std::abs(v0.positiveeta()) > grpV0.etaTrackMax || std::abs(v0.negativeeta()) > grpV0.etaTrackMax) { + return false; + } + if (v0.positivept() < grpV0.ptTrackMin || v0.negativept() < grpV0.ptTrackMin) { + return false; + } + return true; + } + + template + bool selectionPID(const T& candidate) + { + if (cfgMomDepPID) { + if (candidate.p() < 0.5) { + if (!candidate.hasTOF() && TMath::Abs(candidate.tpcNSigmaKa()) < grpPhi.nsigmaCutTPC) { + return true; + } else if (candidate.hasTOF() && TMath::Sqrt(candidate.tpcNSigmaKa() * candidate.tpcNSigmaKa() + candidate.tofNSigmaKa() * candidate.tofNSigmaKa()) < grpPhi.nsigmaCutTOF && candidate.beta() > grpPhi.cutTOFBeta) { + return true; + } + } else if (candidate.hasTOF() && TMath::Sqrt(candidate.tpcNSigmaKa() * candidate.tpcNSigmaKa() + candidate.tofNSigmaKa() * candidate.tofNSigmaKa()) < grpPhi.nsigmaCutTOF && candidate.beta() > grpPhi.cutTOFBeta) { + return true; + } + } else { + if (!candidate.hasTOF() && TMath::Abs(candidate.tpcNSigmaKa()) < grpPhi.nsigmaCutTPC) { + return true; + } else if (candidate.hasTOF() && TMath::Sqrt(candidate.tpcNSigmaKa() * candidate.tpcNSigmaKa() + candidate.tofNSigmaKa() * candidate.tofNSigmaKa()) < grpPhi.nsigmaCutTOF && candidate.beta() > grpPhi.cutTOFBeta) { + return true; + } + } + return false; + } + // Generic 2-prong invariant mass method candidate finder. Only works for non-identical daughters of opposite charge for now. - using PIDTrack = soa::Join; - void processDataInvMass(aod::Collisions::iterator const&, aod::BCsWithTimestamps const&, aod::CFCollRefs const& cfcollisions, aod::CFTrackRefs const& cftracks, PIDTrack const& tracks) + void processDataInvMass(aod::Collisions::iterator const&, aod::BCsWithTimestamps const&, aod::CFCollRefs const& cfcollisions, aod::CFTrackRefs const& cftracks, Filter2Prong::PIDTrack const& tracks) { if (cfcollisions.size() <= 0 || cftracks.size() <= 0) return; // rejected collision @@ -199,6 +411,112 @@ struct Filter2Prong { } } PROCESS_SWITCH(Filter2Prong, processDataInvMass, "Process data generic 2-prong candidates with invariant mass method", false); + + // Phi and V0s invariant mass method candidate finder. Only works for non-identical daughters of opposite charge for now. + void processDataV0(aod::Collisions::iterator const& collision, aod::BCsWithTimestamps const&, aod::CFCollRefs const& cfcollisions, aod::CFTrackRefs const& cftracks, Filter2Prong::PIDTrack const&, aod::V0Datas const& V0s) + { + if (cfcollisions.size() <= 0 || cftracks.size() <= 0) + return; // rejected collision + + for (const auto& v0 : V0s) { // Loop over V0 candidates + if (!isV0TrackSelected(v0)) { // Quality selection for V0 prongs + continue; + } + + const auto& posTrack = v0.template posTrack_as(); + const auto& negTrack = v0.template negTrack_as(); + double massV0 = 0.0; + + // K0s + if (isSelectedV0AsK0s(collision, v0)) { // candidate is K0s + output2ProngTracks(cfcollisions.begin().globalIndex(), + posTrack.globalIndex(), negTrack.globalIndex(), + v0.pt(), v0.eta(), v0.phi(), v0.mK0Short(), aod::cf2prongtrack::K0stoPiPi); + } + + // Lambda and Anti-Lambda + bool LambdaTag = isSelectedV0AsLambda(collision, v0); + bool aLambdaTag = isSelectedV0AsLambda(collision, v0); + + // Note: candidate compatible with Lambda and Anti-Lambda hypothesis are counted twice (once for each hypothesis) + if (LambdaTag) { // candidate is Lambda + massV0 = v0.mLambda(); + output2ProngTracks(cfcollisions.begin().globalIndex(), posTrack.globalIndex(), negTrack.globalIndex(), + v0.pt(), v0.eta(), v0.phi(), massV0, aod::cf2prongtrack::LambdatoPPi); + } + if (aLambdaTag) { // candidate is Anti-lambda + massV0 = v0.mAntiLambda(); + output2ProngTracks(cfcollisions.begin().globalIndex(), posTrack.globalIndex(), negTrack.globalIndex(), + v0.pt(), v0.eta(), v0.phi(), massV0, aod::cf2prongtrack::AntiLambdatoPiP); + } // end of Lambda and Anti-Lambda processing + } // end of loop over V0 candidates + } + PROCESS_SWITCH(Filter2Prong, processDataV0, "Process data V0 candidates with invariant mass method", false); + + // Phi and V0s invariant mass method candidate finder. Only works for non-identical daughters of opposite charge for now. + void processDataPhi(aod::Collisions::iterator const&, aod::BCsWithTimestamps const&, aod::CFCollRefs const& cfcollisions, aod::CFTrackRefs const& cftracks, Filter2Prong::PIDTrack const& tracks) + { + if (cfcollisions.size() <= 0 || cftracks.size() <= 0) + return; // rejected collision + + o2::aod::ITSResponse itsResponse; + + for (const auto& cftrack1 : cftracks) { // Loop over first track + const auto& p1 = tracks.iteratorAt(cftrack1.trackId() - tracks.begin().globalIndex()); + if (p1.sign() != 1) { + continue; + } + if (!selectionTrack(p1)) { + continue; + } + if (grpPhi.ITSPIDSelection && p1.p() < grpPhi.ITSPIDPthreshold.value && !(itsResponse.nSigmaITS(p1) > grpPhi.lowITSPIDNsigma.value && itsResponse.nSigmaITS(p1) < grpPhi.highITSPIDNsigma.value)) { // Check ITS PID condition + continue; + } + if (!selectionPID(p1)) { + continue; + } + if (grpPhi.removefaketrack && isFakeTrack(p1)) { // Check if the track is a fake kaon + continue; + } + + for (const auto& cftrack2 : cftracks) { // Loop over second track + if (cftrack2.globalIndex() == cftrack1.globalIndex()) // Skip if it's the same track as the first one + continue; + + const auto& p2 = tracks.iteratorAt(cftrack2.trackId() - tracks.begin().globalIndex()); + if (p2.sign() != -1) { + continue; + } + if (!selectionTrack(p2)) { + continue; + } + if (!selectionPID(p2)) { + continue; + } + if (grpPhi.ITSPIDSelection && p2.p() < grpPhi.ITSPIDPthreshold.value && !(itsResponse.nSigmaITS(p2) > grpPhi.lowITSPIDNsigma.value && itsResponse.nSigmaITS(p2) < grpPhi.highITSPIDNsigma.value)) { // Check ITS PID condition + continue; + } + if (grpPhi.removefaketrack && isFakeTrack(p2)) { // Check if the track is a fake kaon + continue; + } + if (!selectionPair(p1, p2)) { + continue; + } + + ROOT::Math::PtEtaPhiMVector vec1(p1.pt(), p1.eta(), p1.phi(), cfgImPart1Mass); + ROOT::Math::PtEtaPhiMVector vec2(p2.pt(), p2.eta(), p2.phi(), cfgImPart2Mass); + ROOT::Math::PtEtaPhiMVector s = vec1 + vec2; + if (s.M() < grpPhi.ImMinInvMassPhiMeson || s.M() > grpPhi.ImMaxInvMassPhiMeson) { + continue; + } + float phi = RecoDecay::constrainAngle(s.Phi(), 0.0f); + output2ProngTracks(cfcollisions.begin().globalIndex(), + cftrack1.globalIndex(), cftrack2.globalIndex(), s.pt(), s.eta(), phi, s.M(), aod::cf2prongtrack::PhiToKK); + } // end of loop over second track + } // end of loop over first track + } + PROCESS_SWITCH(Filter2Prong, processDataPhi, "Process data Phi candidates with invariant mass method", false); + }; // struct WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) diff --git a/PWGCF/TableProducer/filterCorrelations.cxx b/PWGCF/TableProducer/filterCorrelations.cxx index 7d5e0b17c38..572e4c6c886 100644 --- a/PWGCF/TableProducer/filterCorrelations.cxx +++ b/PWGCF/TableProducer/filterCorrelations.cxx @@ -8,23 +8,27 @@ // In applying this license CERN does not waive the privileges and immunities // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. -#include +#include "PWGCF/DataModel/CorrelationsDerived.h" +#include "PWGHF/DataModel/CandidateReconstructionTables.h" +#include "PWGHF/DataModel/CandidateSelectionTables.h" + +#include "Common/DataModel/Centrality.h" +#include "Common/DataModel/EventSelection.h" +#include "Common/DataModel/PIDResponseITS.h" +#include "Common/DataModel/TrackSelectionTables.h" -#include "Framework/runDataProcessing.h" -#include "Framework/AnalysisTask.h" -#include "Framework/AnalysisDataModel.h" #include "Framework/ASoAHelpers.h" +#include "Framework/AnalysisDataModel.h" +#include "Framework/AnalysisTask.h" #include "Framework/O2DatabasePDGPlugin.h" - +#include "Framework/runDataProcessing.h" #include "MathUtils/detail/TypeTruncation.h" -#include "PWGCF/DataModel/CorrelationsDerived.h" -#include "Common/DataModel/EventSelection.h" -#include "Common/DataModel/TrackSelectionTables.h" -#include "Common/DataModel/Centrality.h" - #include +#include // required for is_detected +#include + using namespace o2; using namespace o2::framework; using namespace o2::framework::expressions; @@ -47,10 +51,14 @@ using CFMultiplicity = CFMultiplicities::iterator; struct FilterCF { Service pdg; - enum TrackSelectionCuts : uint8_t { + enum TrackSelectionCuts1 : uint8_t { kTrackSelected = BIT(0), kITS5Clusters = BIT(1), - kTPC90CrossedRows = BIT(2) + kTPCCrossedRows = BIT(2), + }; + + enum TrackSelectionCuts2 : uint8_t { + kPIDProton = BIT(1) }; // Configuration @@ -65,9 +73,17 @@ struct FilterCF { O2_DEFINE_CONFIGURABLE(cfgMaxOcc, int, 3000, "maximum occupancy selection") O2_DEFINE_CONFIGURABLE(cfgCollisionFlags, uint16_t, aod::collision::CollisionFlagsRun2::Run2VertexerTracks, "Request collision flags if non-zero (0 = off, 1 = Run2VertexerTracks)") O2_DEFINE_CONFIGURABLE(cfgTransientTables, bool, false, "Output transient tables for collision and track IDs to enable successive filtering tasks") - O2_DEFINE_CONFIGURABLE(cfgTrackSelection, int, 0, "Type of track selection (0 = Run 2/3 without systematics | 1 = Run 3 with systematics)") + O2_DEFINE_CONFIGURABLE(cfgTrackSelection, int, 0, "Type of track selection (0 = Run 2/3 without systematics | 1 = Run 3 with systematics | 2 = Run 3 with proton pid selection)") O2_DEFINE_CONFIGURABLE(cfgMinMultiplicity, float, -1, "Minimum multiplicity considered for filtering (if value positive)") O2_DEFINE_CONFIGURABLE(cfgMcSpecialPDGs, std::vector, {}, "Special MC PDG codes to include in the MC primary particle output (additional to charged particles). Empty = charged particles only.") // needed for some neutral particles + O2_DEFINE_CONFIGURABLE(nsigmaCutTPCProton, float, 3, "proton nsigma TPC") + O2_DEFINE_CONFIGURABLE(nsigmaCutTOFProton, float, 3, "proton nsigma TOF") + O2_DEFINE_CONFIGURABLE(ITSProtonselection, bool, false, "flag for ITS proton nsigma selection") + O2_DEFINE_CONFIGURABLE(nsigmaCutITSProton, float, 3, "proton nsigma ITS") + O2_DEFINE_CONFIGURABLE(dcaxymax, float, 999.f, "maximum dcaxy of tracks") + O2_DEFINE_CONFIGURABLE(dcazmax, float, 999.f, "maximum dcaz of tracks") + O2_DEFINE_CONFIGURABLE(itsnclusters, int, 5, "minimum number of ITS clusters for tracks") + O2_DEFINE_CONFIGURABLE(tpcncrossedrows, int, 80, "minimum number of TPC crossed rows for tracks") // Filters and input definitions Filter collisionZVtxFilter = nabs(aod::collision::posZ) < cfgCutVertex; @@ -128,8 +144,49 @@ struct FilterCF { return false; } - template - uint8_t getTrackType(TTrack& track) + using TrackType = soa::Filtered>; + + template + bool selectionPIDProton(const T& candidate) + { + o2::aod::ITSResponse itsResponse; + + if (ITSProtonselection && candidate.pt() <= 0.6 && !(itsResponse.nSigmaITS(candidate) > nsigmaCutITSProton)) { + return false; + } + if (ITSProtonselection && candidate.pt() > 0.6 && candidate.pt() <= 0.8 && !(itsResponse.nSigmaITS(candidate) > nsigmaCutITSProton)) { + return false; + } + + if (candidate.hasTOF()) { + if (candidate.pt() < 0.7 && std::abs(candidate.tpcNSigmaPr()) < nsigmaCutTPCProton) { + return true; + } + if (candidate.p() >= 0.7 && std::abs(candidate.tpcNSigmaPr()) < nsigmaCutTPCProton && std::abs(candidate.tofNSigmaPr()) < nsigmaCutTOFProton) { + return true; + } + } else { + if (std::abs(candidate.tpcNSigmaPr()) < nsigmaCutTPCProton) { + return true; + } + } + return false; + } + + template + struct HasProtonPID : std::false_type { + }; + + template + struct HasProtonPID().tpcNSigmaPr()), + decltype(std::declval().tofNSigmaPr()), + decltype(std::declval().hasTOF())>> + : std::true_type { + }; + + template + uint8_t getTrackType(const TTrack& track) { if (cfgTrackSelection == 0) { if (track.isGlobalTrack()) { @@ -142,20 +199,33 @@ struct FilterCF { uint8_t trackType = 0; if (track.isGlobalTrack()) { trackType |= kTrackSelected; - if (track.itsNCls() >= 5) { + if (track.itsNCls() >= itsnclusters) { trackType |= kITS5Clusters; } - if (track.tpcNClsCrossedRows() >= 90) { - trackType |= kTPC90CrossedRows; + if (track.tpcNClsCrossedRows() >= tpcncrossedrows) { + trackType |= kTPCCrossedRows; + } + } + return trackType; + } else if (cfgTrackSelection == 2) { + uint8_t trackType = 0; + if constexpr (HasProtonPID::value) { + if (track.isGlobalTrack() && (track.itsNCls() >= itsnclusters) && (track.tpcNClsCrossedRows() >= tpcncrossedrows) && selectionPIDProton(track)) { + trackType |= kPIDProton; } } return trackType; } + LOGF(fatal, "Invalid setting for cfgTrackSelection: %d", cfgTrackSelection.value); return 0; } - void processData(soa::Filtered>::iterator const& collision, aod::BCsWithTimestamps const&, soa::Filtered> const& tracks) + /// \brief Templetized process data for a given collision and its associated tracks + /// \param collision The collision object containing information about the collision + /// \param tracks The collection of tracks associated with the collision + template + void processDataT(const C1& collision, const T1& tracks) { if (cfgVerbosity > 0) { LOGF(info, "processData: Tracks for collision: %d | Vertex: %.1f (%d) | INT7: %d | Multiplicity: %.1f", tracks.size(), collision.posZ(), collision.flags(), collision.sel7(), collision.multiplicity()); @@ -165,13 +235,16 @@ struct FilterCF { return; } - auto bc = collision.bc_as(); + auto bc = collision.template bc_as(); outputCollisions(bc.runNumber(), collision.posZ(), collision.multiplicity(), bc.timestamp()); if (cfgTransientTables) outputCollRefs(collision.globalIndex()); - for (auto& track : tracks) { + if ((std::abs(track.dcaXY()) > dcaxymax) || (std::abs(track.dcaZ()) > dcazmax)) { + continue; + } + outputTracks(outputCollisions.lastIndex(), track.pt(), track.eta(), track.phi(), track.sign(), getTrackType(track)); if (cfgTransientTables) outputTrackRefs(collision.globalIndex(), track.globalIndex()); @@ -180,15 +253,31 @@ struct FilterCF { etaphi->Fill(collision.multiplicity(), track.eta(), track.phi()); } } + + void processData(soa::Filtered>::iterator const& collision, aod::BCsWithTimestamps const&, soa::Filtered> const& tracks) + { + processDataT(collision, tracks); + } PROCESS_SWITCH(FilterCF, processData, "Process data", true); - // NOTE not filtering collisions here because in that case there can be tracks referring to MC particles which are not part of the selected MC collisions - Preslice perMcCollision = aod::mcparticle::mcCollisionId; - Preslice perCollision = aod::track::collisionId; - void processMC(aod::McCollisions const& mcCollisions, aod::McParticles const& allParticles, - soa::Join const& allCollisions, - soa::Filtered> const& tracks, - aod::BCsWithTimestamps const&) + void processDataPid(soa::Filtered>::iterator const& collision, aod::BCsWithTimestamps const&, soa::Filtered> const& tracks) + { + processDataT(collision, tracks); + } + PROCESS_SWITCH(FilterCF, processDataPid, "Process data with PID", false); + + /// \brief Process MC data for a given set of MC collisions and associated particles and tracks + /// \param mcCollisions The collection of MC collisions + /// \param allParticles The collection of all MC particles + /// \param allCollisions The collection of all collisions, joined with MC collision labels and + /// event selections + /// \param tracks The collection of tracks, filtered by selection criteria + /// \param bcs The collection of bunch crossings with timestamps + template + void processMCT(aod::McCollisions const& mcCollisions, aod::McParticles const& allParticles, + soa::Join const& allCollisions, + T1 const& tracks, + aod::BCsWithTimestamps const&) { mcReconstructedCache.reserve(allParticles.size()); mcParticleLabelsCache.reserve(allParticles.size()); @@ -286,8 +375,7 @@ struct FilterCF { LOGP(fatal, "processMC: Track {} is referring to a MC particle which we do not store {} {} (reco flag {})", track.index(), track.mcParticleId(), mcParticleId, static_cast(mcReconstructedCache[track.mcParticleId()])); } } - outputTracks(outputCollisions.lastIndex(), - truncateFloatFraction(track.pt()), truncateFloatFraction(track.eta()), truncateFloatFraction(track.phi()), track.sign(), getTrackType(track)); + outputTracks(outputCollisions.lastIndex(), truncateFloatFraction(track.pt()), truncateFloatFraction(track.eta()), truncateFloatFraction(track.phi()), track.sign(), getTrackType(track)); outputTrackLabels(mcParticleId); if (cfgTransientTables) outputTrackRefs(collision.globalIndex(), track.globalIndex()); @@ -297,8 +385,29 @@ struct FilterCF { } } } + + // NOTE not filtering collisions here because in that case there can be tracks referring to MC particles which are not part of the selected MC collisions + Preslice perMcCollision = aod::mcparticle::mcCollisionId; + Preslice perCollision = aod::track::collisionId; + void processMC(aod::McCollisions const& mcCollisions, aod::McParticles const& allParticles, + soa::Join const& allCollisions, + soa::Filtered> const& tracks, + aod::BCsWithTimestamps const& bcs) + { + processMCT(mcCollisions, allParticles, allCollisions, tracks, bcs); + } PROCESS_SWITCH(FilterCF, processMC, "Process MC", false); + // NOTE not filtering collisions here because in that case there can be tracks referring to MC particles which are not part of the selected MC collisions + void processMCPid(aod::McCollisions const& mcCollisions, aod::McParticles const& allParticles, + soa::Join const& allCollisions, + soa::Filtered> const& tracks, + aod::BCsWithTimestamps const& bcs) + { + processMCT(mcCollisions, allParticles, allCollisions, tracks, bcs); + } + PROCESS_SWITCH(FilterCF, processMCPid, "Process MC with PID", false); + void processMCGen(aod::McCollisions::iterator const& mcCollision, aod::McParticles const& particles) { float multiplicity = 0.0f;