From 536cd5ad9b7f4a741a839b581c8bc955b086cc30 Mon Sep 17 00:00:00 2001 From: sigurd Date: Mon, 9 Dec 2024 16:37:49 +0100 Subject: [PATCH 1/7] Add open charm MC analyses. Remove KKPi process function. --- PWGDQ/Core/CutsLibrary.cxx | 36 + PWGDQ/Core/HistogramsLibrary.cxx | 2 + PWGDQ/Core/MCSignalLibrary.cxx | 97 ++ PWGDQ/Core/VarManager.cxx | 2 + PWGDQ/Core/VarManager.h | 35 +- PWGDQ/Tasks/dqEfficiency_withAssoc.cxx | 1159 +++++++++++++++++++++++- PWGDQ/Tasks/tableReader_withAssoc.cxx | 10 +- 7 files changed, 1297 insertions(+), 44 deletions(-) diff --git a/PWGDQ/Core/CutsLibrary.cxx b/PWGDQ/Core/CutsLibrary.cxx index 20ae5000bc8..2f6ad8a958e 100644 --- a/PWGDQ/Core/CutsLibrary.cxx +++ b/PWGDQ/Core/CutsLibrary.cxx @@ -3172,6 +3172,27 @@ AnalysisCompositeCut* o2::aod::dqcuts::GetCompositeCut(const char* cutName) return cut; } + if (!nameStr.compare("pairCosPointingPos")) { + cut->AddCut(GetAnalysisCut("pairCosPointingPos")); + return cut; + } + + if (!nameStr.compare("pairCosPointingNeg90")) { + cut->AddCut(GetAnalysisCut("pairCosPointingNeg90")); + return cut; + } + + if (!nameStr.compare("pairCosPointingNeg85")) { + cut->AddCut(GetAnalysisCut("pairCosPointingNeg85")); + return cut; + } + + if (!nameStr.compare("pairTauxyzProjectedCosPointing1")) { + cut->AddCut(GetAnalysisCut("pairCosPointingNeg")); + cut->AddCut(GetAnalysisCut("pairTauxyzProjected1")); + return cut; + } + // ------------------------------------------------------------------------------------------------- // // Below are a list of single electron single muon and in order or optimize the trigger @@ -5962,6 +5983,21 @@ AnalysisCut* o2::aod::dqcuts::GetAnalysisCut(const char* cutName) return cut; } + if (!nameStr.compare("pairCosPointingPos")) { + cut->AddCut(VarManager::kCosPointingAngle, 0.9, 1000.); + return cut; + } + + if (!nameStr.compare("pairCosPointingNeg90")) { + cut->AddCut(VarManager::kCosPointingAngle, -1000., -0.9); + return cut; + } + + if (!nameStr.compare("pairCosPointingNeg85")) { + cut->AddCut(VarManager::kCosPointingAngle, -1000., -0.85); + return cut; + } + // ------------------------------------------------------------------------------------------------- // // Below are a list of single electron single muon and pair selection in order or optimize the trigger diff --git a/PWGDQ/Core/HistogramsLibrary.cxx b/PWGDQ/Core/HistogramsLibrary.cxx index c55f34bf1c9..c8eac97666b 100644 --- a/PWGDQ/Core/HistogramsLibrary.cxx +++ b/PWGDQ/Core/HistogramsLibrary.cxx @@ -941,6 +941,7 @@ void o2::aod::dqhistograms::DefineHistograms(HistogramManager* hm, const char* h hm->AddHistogram(histClass, "LxyProj_Pt", "", false, 10, 0.0, 20.0, VarManager::kPt, 1000, -1.0, 1.0, VarManager::kVertexingLxyProjected); hm->AddHistogram(histClass, "LxyProj_Mass_Pt", "", false, 50, 2.0, 4.0, VarManager::kMass, 10, 0.0, 20.0, VarManager::kPt, 1000, -1.0, 1.0, VarManager::kVertexingLxyProjected); hm->AddHistogram(histClass, "LzProj_Mass_Pt", "", false, 50, 2.0, 4.0, VarManager::kMass, 10, 0.0, 20.0, VarManager::kPt, 1000, -1.0, 1.0, VarManager::kVertexingLzProjected); + hm->AddHistogram(histClass, "CosPointingAngle", "", false, 200, -1.0, 1.0, VarManager::kCosPointingAngle); } if (subGroupStr.Contains("kalman-filter")) { @@ -1378,6 +1379,7 @@ void o2::aod::dqhistograms::DefineHistograms(HistogramManager* hm, const char* h hm->AddHistogram(histClass, "MassD0region_Rapidity", "", false, 140, 1.5, 2.2, VarManager::kMass, 10, -0.8, 0.8, VarManager::kRap); hm->AddHistogram(histClass, "MassD0region_eta", "", false, 140, 1.5, 2.2, VarManager::kMass, 40, -2., 2., VarManager::kEta); hm->AddHistogram(histClass, "MassD0region_TauxyzProj", "", false, 140, 1.5, 2.2, VarManager::kMass, 200, -0.03, 0.03, VarManager::kVertexingTauxyzProjected); + hm->AddHistogram(histClass, "MassD0region_CosPointing", "", false, 140, 1.5, 2.2, VarManager::kMass, 200, -1.0, 1.0, VarManager::kCosPointingAngle); hm->AddHistogram(histClass, "MassD0region_VtxNContribReal", "", false, 140, 1.5, 2.2, VarManager::kMass, 50, 0, 50, VarManager::kVtxNcontribReal); hm->AddHistogram(histClass, "MassD0region_Rapidity_AveragePt", "", true, 140, 1.5, 2.2, VarManager::kMass, 10, -0.8, 0.8, VarManager::kRap, 150, 0.0, 30.0, VarManager::kPt); } diff --git a/PWGDQ/Core/MCSignalLibrary.cxx b/PWGDQ/Core/MCSignalLibrary.cxx index 22f445252f9..ca390993dd9 100644 --- a/PWGDQ/Core/MCSignalLibrary.cxx +++ b/PWGDQ/Core/MCSignalLibrary.cxx @@ -1185,6 +1185,103 @@ MCSignal* o2::aod::dqmcsignals::GetMCSignal(const char* name) return signal; } + //------------------------------------------------------------------------------------ + + if (!nameStr.compare("D0")) { + MCProng prong(1, {Pdg::kD0}, {true}, {false}, {0}, {0}, {false}); + signal = new MCSignal(name, "D0", {prong}, {-1}); + return signal; + } + if (!nameStr.compare("KPiFromD0")) { + MCProng prongKaon(2, {321, Pdg::kD0}, {true, true}, {false, false}, {0, 0}, {0, 0}, {false, false}); + MCProng prongPion(2, {211, Pdg::kD0}, {true, true}, {false, false}, {0, 0}, {0, 0}, {false, false}); + signal = new MCSignal(name, "Kaon and pion pair from D0", {prongKaon, prongPion}, {1, 1}); + return signal; + } + if (!nameStr.compare("Dcharged")) { + MCProng prong(1, {Pdg::kDPlus}, {true}, {false}, {0}, {0}, {false}); + signal = new MCSignal(name, "D+/-", {prong}, {-1}); + return signal; + } + if (!nameStr.compare("Dplus")) { + MCProng prong(1, {Pdg::kDPlus}, {false}, {false}, {0}, {0}, {false}); + signal = new MCSignal(name, "D+", {prong}, {-1}); + return signal; + } + if (!nameStr.compare("Dminus")) { + MCProng prong(1, {-Pdg::kDPlus}, {false}, {false}, {0}, {0}, {false}); + signal = new MCSignal(name, "D+", {prong}, {-1}); + return signal; + } + if (!nameStr.compare("KPiPiFromDcharged")) { + MCProng prongKaon(2, {321, Pdg::kDPlus}, {true, true}, {false, false}, {0, 0}, {0, 0}, {false, false}); + MCProng prongPion(2, {211, Pdg::kDPlus}, {true, true}, {false, false}, {0, 0}, {0, 0}, {false, false}); + signal = new MCSignal(name, "Kaon pion pion triplet from D+/-", {prongKaon, prongPion, prongPion}, {1, 1, 1}); + return signal; + } + if (!nameStr.compare("KPiPiFromDplus")) { + MCProng prongKaon(2, {-321, Pdg::kDPlus}, {false, false}, {false, false}, {0, 0}, {0, 0}, {false, false}); + MCProng prongPion(2, {211, Pdg::kDPlus}, {false, false}, {false, false}, {0, 0}, {0, 0}, {false, false}); + signal = new MCSignal(name, "Kaon pion pion triplet from D+", {prongKaon, prongPion, prongPion}, {1, 1, 1}); + return signal; + } + if (!nameStr.compare("KPiPiFromDminus")) { + MCProng prongKaon(2, {321, -Pdg::kDPlus}, {false, false}, {false, false}, {0, 0}, {0, 0}, {false, false}); + MCProng prongPion(2, {-211, -Pdg::kDPlus}, {false, false}, {false, false}, {0, 0}, {0, 0}, {false, false}); + signal = new MCSignal(name, "Kaon pion pion triplet from D-", {prongKaon, prongPion, prongPion}, {1, 1, 1}); + return signal; + } + if (!nameStr.compare("Dstar")) { + MCProng prong(1, {Pdg::kDStar}, {true}, {false}, {0}, {0}, {false}); + signal = new MCSignal(name, "D*", {prong}, {-1}); + return signal; + } + if (!nameStr.compare("pionFromDstar")) { + MCProng prong(2, {211, Pdg::kDStar}, {true, true}, {false, false}, {0, 0}, {0, 0}, {false, false}); + signal = new MCSignal(name, "Pions from D* decays", {prong}, {1}); + return signal; + } + if (!nameStr.compare("D0FromDstar")) { + MCProng prong(2, {Pdg::kD0, Pdg::kDStar}, {true, true}, {false, false}, {0, 0}, {0, 0}, {false, false}); + signal = new MCSignal(name, "D0 from D* decays", {prong}, {1}); + return signal; + } + if (!nameStr.compare("KFromD0FromDstar")) { + MCProng prong(3, {321, Pdg::kD0, Pdg::kDStar}, {true, true, true}, {false, false, false}, {0, 0, 0}, {0, 0, 0}, {false, false, false}); + signal = new MCSignal(name, "Kaons from D0 from D* decays", {prong}, {1}); + return signal; + } + if (!nameStr.compare("PiFromD0FromDstar")) { + MCProng prong(3, {211, Pdg::kD0, Pdg::kDStar}, {true, true, true}, {false, false, false}, {0, 0, 0}, {0, 0, 0}, {false, false, false}); + signal = new MCSignal(name, "Pions from D0 from D* decays", {prong}, {1}); + return signal; + } + if (!nameStr.compare("KPiFromD0FromDstar")) { + MCProng prongKaon(3, {321, Pdg::kD0, Pdg::kDStar}, {true, true, true}, {false, false, false}, {0, 0, 0}, {0, 0, 0}, {false, false, false}); + MCProng prongPion(3, {321, Pdg::kD0, Pdg::kDStar}, {true, true, true}, {false, false, false}, {0, 0, 0}, {0, 0, 0}, {false, false, false}); + signal = new MCSignal(name, "Kaon and pion pair from D0 from D* decay", {prongKaon, prongPion}, {1, 1}); + return signal; + } + if (!nameStr.compare("KPiPiFromD0FromDstar")) { + MCProng prongKaon(3, {321, Pdg::kD0, Pdg::kDStar}, {true, true, true}, {false, false, false}, {0, 0, 0}, {0, 0, 0}, {false, false, false}); + MCProng prongPionSecondary(3, {211, Pdg::kD0, Pdg::kDStar}, {true, true, true}, {false, false, false}, {0, 0, 0}, {0, 0, 0}, {false, false, false}); + MCProng prongPion(2, {211, Pdg::kDStar}, {true, true}, {false, false}, {0, 0}, {0, 0}, {false, false}); + signal = new MCSignal(name, "Kaon pion pion triplet from D*", {prongKaon, prongPionSecondary, prongPion}, {2, 2, 1}); + return signal; + } + if (!nameStr.compare("PiPiPiFromD0FromDstar")) { + MCProng prongPionSecondary(3, {211, Pdg::kD0, Pdg::kDStar}, {true, true, true}, {false, false, false}, {0, 0, 0}, {0, 0, 0}, {false, false, false}); + MCProng prongPion(2, {211, Pdg::kDStar}, {true, true}, {false, false}, {0, 0}, {0, 0}, {false, false}); + signal = new MCSignal(name, "Kaon pion pion triplet from D*", {prongPionSecondary, prongPionSecondary, prongPion}, {2, 2, 1}); + return signal; + } + if (!nameStr.compare("KFromDplus")) { + MCProng prong(2, {321, Pdg::kDPlus}, {true, true}, {false, false}, {0, 0}, {0, 0}, {false, false}, false, {502}, {true}); + prong.SetSourceBit(0, MCProng::kPhysicalPrimary); + signal = new MCSignal(name, "Kaons from D+/- decays", {prong}, {-1}); + return signal; + } + //-------------------------------------------------------------------------------- if (!nameStr.compare("JpsiFromChic0")) { diff --git a/PWGDQ/Core/VarManager.cxx b/PWGDQ/Core/VarManager.cxx index 77de7838f8e..91315f7069c 100644 --- a/PWGDQ/Core/VarManager.cxx +++ b/PWGDQ/Core/VarManager.cxx @@ -642,6 +642,8 @@ void VarManager::SetDefaultVarNames() fgVariableUnits[kVertexingLzOverErr] = ""; fgVariableNames[kVertexingLxyzOverErr] = "Pair Lxyz/DLxyz"; fgVariableUnits[kVertexingLxyzOverErr] = ""; + fgVariableNames[kCosPointingAngle] = "Cos #theta_{pointing}"; + fgVariableUnits[kCosPointingAngle] = ""; fgVariableNames[kKFTrack0DCAxyz] = "Daughter0 DCAxyz"; fgVariableUnits[kKFTrack0DCAxyz] = "cm"; fgVariableNames[kKFTrack1DCAxyz] = "Daughter1 DCAxyz"; diff --git a/PWGDQ/Core/VarManager.h b/PWGDQ/Core/VarManager.h index bf328b6fa0c..6a2443dc3bb 100644 --- a/PWGDQ/Core/VarManager.h +++ b/PWGDQ/Core/VarManager.h @@ -146,7 +146,6 @@ class VarManager : public TObject kDecayToKPi, // e.g. D0 -> K+ pi- or cc. kTripleCandidateToKPiPi, // e.g. D+ -> K- pi+ pi+ kTripleCandidateToPKPi, // e.g. Lambda_c -> p K- pi+ - kTripleCandidateToKKPi, // e.g. D_s -> K+ K- pi+ kNMaxCandidateTypes }; @@ -2885,21 +2884,6 @@ void VarManager::FillTriple(T1 const& t1, T2 const& t2, T3 const& t3, float* val values[kPhi] = v123.Phi(); values[kRap] = -v123.Rapidity(); } - - if (pairType == kTripleCandidateToKKPi) { - float m1 = o2::constants::physics::MassKaonCharged; - float m2 = o2::constants::physics::MassPionCharged; - - ROOT::Math::PtEtaPhiMVector v1(t1.pt(), t1.eta(), t1.phi(), m1); - ROOT::Math::PtEtaPhiMVector v2(t2.pt(), t2.eta(), t2.phi(), m1); - ROOT::Math::PtEtaPhiMVector v3(t3.pt(), t3.eta(), t3.phi(), m2); - ROOT::Math::PtEtaPhiMVector v123 = v1 + v2 + v3; - values[kMass] = v123.M(); - values[kPt] = v123.Pt(); - values[kEta] = v123.Eta(); - values[kPhi] = v123.Phi(); - values[kRap] = -v123.Rapidity(); - } } template @@ -3063,6 +3047,21 @@ void VarManager::FillTripleMC(T1 const& t1, T2 const& t2, T3 const& t3, float* v values[kPt1] = t1.pt(); values[kPt2] = t2.pt(); } + + if (pairType == kTripleCandidateToKPiPi) { + float m1 = o2::constants::physics::MassKaonCharged; + float m2 = o2::constants::physics::MassPionCharged; + + ROOT::Math::PtEtaPhiMVector v1(t1.pt(), t1.eta(), t1.phi(), m1); + ROOT::Math::PtEtaPhiMVector v2(t2.pt(), t2.eta(), t2.phi(), m2); + ROOT::Math::PtEtaPhiMVector v3(t3.pt(), t3.eta(), t3.phi(), m2); + ROOT::Math::PtEtaPhiMVector v123 = v1 + v2 + v3; + values[kMass] = v123.M(); + values[kPt] = v123.Pt(); + values[kEta] = v123.Eta(); + values[kPhi] = v123.Phi(); + values[kRap] = -v123.Rapidity(); + } } template @@ -3574,6 +3573,10 @@ void VarManager::FillTripletVertexing(C const& collision, T const& t1, T const& values[kVertexingTauzErr] = values[kVertexingLzErr] * v123.M() / (TMath::Abs(v123.Pz()) * o2::constants::physics::LightSpeedCm2NS); values[kVertexingTauxyErr] = values[kVertexingLxyErr] * v123.M() / (v123.Pt() * o2::constants::physics::LightSpeedCm2NS); + values[kCosPointingAngle] = ((collision.posX() - secondaryVertex[0]) * v123.Px() + + (collision.posY() - secondaryVertex[1]) * v123.Py() + + (collision.posZ() - secondaryVertex[2]) * v123.Pz()) / + (v123.P() * values[VarManager::kVertexingLxyz]); // run 2 definitions: Decay length projected onto the momentum vector of the candidate values[kVertexingLzProjected] = (secondaryVertex[2] - collision.posZ()) * v123.Pz(); values[kVertexingLzProjected] = values[kVertexingLzProjected] / TMath::Sqrt(v123.Pz() * v123.Pz()); diff --git a/PWGDQ/Tasks/dqEfficiency_withAssoc.cxx b/PWGDQ/Tasks/dqEfficiency_withAssoc.cxx index 62dfcd46652..6391e4d20fd 100644 --- a/PWGDQ/Tasks/dqEfficiency_withAssoc.cxx +++ b/PWGDQ/Tasks/dqEfficiency_withAssoc.cxx @@ -29,6 +29,7 @@ #include "Framework/AnalysisTask.h" #include "Framework/AnalysisDataModel.h" #include "Framework/ASoAHelpers.h" +#include "Framework/AnalysisHelpers.h" #include "PWGDQ/DataModel/ReducedInfoTables.h" #include "PWGDQ/Core/VarManager.h" #include "PWGDQ/Core/HistogramManager.h" @@ -106,6 +107,7 @@ using MyBarrelTracksWithCov = soa::Join; using MyBarrelTracksWithCovWithAmbiguitiesWithColl = soa::Join; using MyDielectronCandidates = soa::Join; +using MyDitrackCandidates = soa::Join; using MyDimuonCandidates = soa::Join; using MyMuonTracks = soa::Join; using MyMuonTracksWithCov = soa::Join; @@ -1878,6 +1880,975 @@ struct AnalysisSameEventPairing { PROCESS_SWITCH(AnalysisSameEventPairing, processDummy, "Dummy function, enabled only if none of the others are enabled", false); }; +// Run pairing for resonance with legs fulfilling separate cuts (asymmetric decay channel) +struct AnalysisAsymmetricPairing { + + Produces ditrackList; + Produces ditrackExtraList; + + o2::base::MatLayerCylSet* fLUT = nullptr; + int fCurrentRun; // needed to detect if the run changed and trigger update of calibrations etc. + + // Output objects + OutputObj fOutputList{"output"}; + + // Configurables + Configurable fConfigLegCuts{"cfgLegCuts", "", ":[:],[:[:],...]"}; + Configurable fConfigLegAFilterMask{"cfgLegAFilterMask", 0, "Filter mask corresponding to cuts in event-selection"}; + Configurable fConfigLegBFilterMask{"cfgLegBFilterMask", 0, "Filter mask corresponding to cuts in event-selection"}; + Configurable fConfigLegCFilterMask{"cfgLegCFilterMask", 0, "Filter mask corresponding to cuts in event-selection"}; + Configurable fConfigCommonTrackCuts{"cfgCommonTrackCuts", "", "Comma separated list of cuts to be applied to all legs"}; + Configurable fConfigPairCuts{"cfgPairCuts", "", "Comma separated list of pair cuts"}; + Configurable fConfigSkipAmbiguousIdCombinations{"cfgSkipAmbiguousIdCombinations", true, "Choose whether to skip pairs/triples which pass a stricter combination of cuts, e.g. KKPi triplets for D+ -> KPiPi"}; + + Configurable fConfigHistogramSubgroups{"cfgAsymmetricPairingHistogramsSubgroups", "barrel,vertexing", "Comma separated list of asymmetric-pairing histogram subgroups"}; + Configurable fConfigSameSignHistograms{"cfgSameSignHistograms", false, "Include same sign pair histograms for 2-prong decays"}; + // Configurable fConfigAmbiguousHistograms{"cfgAmbiguousHistograms", false, "Include separate histograms for pairs/triplets with ambiguous tracks"}; + Configurable fConfigQA{"cfgQA", false, "If true, fill QA histograms"}; + + Configurable fConfigCcdbUrl{"ccdb-url", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; + Configurable fConfigGRPMagPath{"grpmagPath", "GLO/Config/GRPMagField", "CCDB path of the GRPMagField object"}; + Configurable fConfigUseRemoteField{"cfgUseRemoteField", false, "Choose whether to fetch the magnetic field from ccdb or set it manually"}; + Configurable fConfigMagField{"cfgMagField", 5.0f, "Manually set magnetic field"}; + + Configurable fConfigUseKFVertexing{"cfgUseKFVertexing", false, "Use KF Particle for secondary vertex reconstruction (DCAFitter is used by default)"}; + Configurable fConfigUseAbsDCA{"cfgUseAbsDCA", false, "Use absolute DCA minimization instead of chi^2 minimization in secondary vertexing"}; + Configurable fConfigPropToPCA{"cfgPropToPCA", false, "Propagate tracks to secondary vertex"}; + Configurable fConfigLutPath{"lutPath", "GLO/Param/MatLUT", "Path of the Lut parametrization"}; + + Configurable fConfigRunMCGenPair{"cfgRunMCGenPair", false, "Do pairing of true MC particles"}; + Configurable fConfigMCGenSignals{"cfgBarrelMCGenSignals", "", "Comma separated list of MC signals (generated)"}; + Configurable fConfigMCRecSignals{"cfgBarrelMCRecSignals", "", "Comma separated list of MC signals (reconstructed)"}; + + Service fCCDB; + + HistogramManager* fHistMan; + + std::map> fTrackHistNames; + std::map> fBarrelHistNamesMCmatched; + std::vector fPairCuts; + + std::vector fRecMCSignals; + std::vector fGenMCSignals; + + // Filter masks to find legs in BarrelTrackCuts table + uint32_t fLegAFilterMask; + uint32_t fLegBFilterMask; + uint32_t fLegCFilterMask; + // Map tracking which pair of leg cuts the track cuts participate in + std::map fTrackCutFilterMasks; + // Filter map for common track cuts + uint32_t fCommonTrackCutMask; + // Map tracking which common track cut the track cuts correspond to + std::map fCommonTrackCutFilterMasks; + + int fNLegCuts; + int fNPairCuts = 0; + int fNCommonTrackCuts; + bool fHasTwoProngGenMCsignals = false; + bool fHasThreeProngGenMCsignals = false; + + Preslice> trackAssocsPerCollision = aod::reducedtrack_association::reducedeventId; + + // Partitions for triplets and asymmetric pairs + Partition> legACandidateAssocs = (o2::aod::dqanalysisflags::isBarrelSelected & fConfigLegAFilterMask) > static_cast(0); + Partition> legBCandidateAssocs = (o2::aod::dqanalysisflags::isBarrelSelected & fConfigLegBFilterMask) > static_cast(0); + Partition> legCCandidateAssocs = (o2::aod::dqanalysisflags::isBarrelSelected & fConfigLegCFilterMask) > static_cast(0); + + void init(o2::framework::InitContext& context) + { + if (context.mOptions.get("processDummy")) { + return; + } + + TString histNames = ""; + std::vector names; + + // Get the leg cut filter maps + fLegAFilterMask = fConfigLegAFilterMask.value; + fLegBFilterMask = fConfigLegBFilterMask.value; + fLegCFilterMask = fConfigLegCFilterMask.value; + // Get the pair cuts + TString cutNamesStr = fConfigPairCuts.value; + if (!cutNamesStr.IsNull()) { + std::unique_ptr objArray(cutNamesStr.Tokenize(",")); + for (int icut = 0; icut < objArray->GetEntries(); ++icut) { + fPairCuts.push_back(*dqcuts::GetCompositeCut(objArray->At(icut)->GetName())); + } + } + + // Setting the MC rec signal names + TString sigNamesStr = fConfigMCRecSignals.value; + std::unique_ptr objRecSigArray(sigNamesStr.Tokenize(",")); + + for (int isig = 0; isig < objRecSigArray->GetEntries(); ++isig) { + MCSignal* sig = o2::aod::dqmcsignals::GetMCSignal(objRecSigArray->At(isig)->GetName()); + if (sig) { + fRecMCSignals.push_back(*sig); + } + } + + // Get the barrel track selection cuts + string tempCuts; + getTaskOptionValue(context, "analysis-track-selection", "cfgTrackCuts", tempCuts, false); + TString tempCutsStr = tempCuts; + std::unique_ptr objArray(tempCutsStr.Tokenize(",")); + // Get the common leg cuts + int commonCutIdx; + TString commonNamesStr = fConfigCommonTrackCuts.value; + if (!commonNamesStr.IsNull()) { // if common track cuts + std::unique_ptr objArrayCommon(commonNamesStr.Tokenize(",")); + fNCommonTrackCuts = objArrayCommon->GetEntries(); + for (int icut = 0; icut < fNCommonTrackCuts; ++icut) { + commonCutIdx = objArray->IndexOf(objArrayCommon->At(icut)); + if (commonCutIdx >= 0) { + fCommonTrackCutMask |= static_cast(1) << objArray->IndexOf(objArrayCommon->At(icut)); + fCommonTrackCutFilterMasks[icut] = static_cast(1) << objArray->IndexOf(objArrayCommon->At(icut)); + } else { + LOGF(fatal, "Common track cut %s was not calculated upstream. Check the config!", objArrayCommon->At(icut)->GetName()); + } + } + } + // Check that the leg cut masks make sense + if (static_cast(std::floor(TMath::Log2(fLegAFilterMask))) + 1 > objArray->GetEntries()) { + LOGF(fatal, "fConfigLegAFilterMask has highest bit at position %d, but track-selection only has %d cuts!", static_cast(std::floor(TMath::Log2(fLegAFilterMask))) + 1, objArray->GetEntries()); + } + if (static_cast(std::floor(TMath::Log2(fLegBFilterMask))) + 1 > objArray->GetEntries()) { + LOGF(fatal, "fConfigLegBFilterMask has highest bit at position %d, but track-selection only has %d cuts!", static_cast(std::floor(TMath::Log2(fLegBFilterMask))) + 1, objArray->GetEntries()); + } + if (static_cast(std::floor(TMath::Log2(fLegCFilterMask))) + 1 > objArray->GetEntries()) { + LOGF(fatal, "fConfigLegCFilterMask has highest bit at position %d, but track-selection only has %d cuts!", static_cast(std::floor(TMath::Log2(fLegCFilterMask))) + 1, objArray->GetEntries()); + } + + // Get the cuts defining the legs + uint32_t fConstructedLegAFilterMask = 0; + uint32_t fConstructedLegBFilterMask = 0; + uint32_t fConstructedLegCFilterMask = 0; + TString legCutsStr = fConfigLegCuts.value; + std::unique_ptr objArrayLegs(legCutsStr.Tokenize(",")); + if (objArrayLegs->GetEntries() == 0) { + LOG(fatal) << "No cuts defining legs. Check the config!"; + } + fNLegCuts = objArrayLegs->GetEntries(); + std::vector isThreeProng; + int legAIdx; + int legBIdx; + int legCIdx; + // Loop over leg defining cuts + for (int icut = 0; icut < fNLegCuts; ++icut) { + TString legsStr = objArrayLegs->At(icut)->GetName(); + std::unique_ptr legs(legsStr.Tokenize(":")); + if (legs->GetEntries() == 3) { + isThreeProng.push_back(true); + } else if (legs->GetEntries() == 2) { + isThreeProng.push_back(false); + } else { + LOGF(fatal, "Leg cuts %s has the wrong format and could not be parsed!", legsStr.Data()); + continue; + } + // Find leg cuts in the track selection cuts + legAIdx = objArray->IndexOf(legs->At(0)); + if (legAIdx >= 0) { + fConstructedLegAFilterMask |= static_cast(1) << legAIdx; + fTrackCutFilterMasks[icut] |= static_cast(1) << legAIdx; + } else { + LOGF(fatal, "Leg A cut %s was not calculated upstream. Check the config!", legs->At(0)->GetName()); + continue; + } + legBIdx = objArray->IndexOf(legs->At(1)); + if (legBIdx >= 0) { + fConstructedLegBFilterMask |= static_cast(1) << legBIdx; + fTrackCutFilterMasks[icut] |= static_cast(1) << legBIdx; + } else { + LOGF(fatal, "Leg B cut %s was not calculated upstream. Check the config!", legs->At(1)->GetName()); + continue; + } + if (isThreeProng[icut]) { + legCIdx = objArray->IndexOf(legs->At(2)); + if (legCIdx >= 0) { + fConstructedLegCFilterMask |= static_cast(1) << legCIdx; + fTrackCutFilterMasks[icut] |= static_cast(1) << legCIdx; + } else { + LOGF(fatal, "Leg C cut %s was not calculated upstream. Check the config!", legs->At(2)->GetName()); + continue; + } + } + if (isThreeProng[icut]) { + names = { + Form("TripletsBarrelSE_%s", legsStr.Data())}; + histNames += Form("%s;", names[0].Data()); + if (fConfigQA) { + names.push_back(Form("TripletsBarrelSE_ambiguous_%s", legsStr.Data())); + histNames += Form("%s;", names[1].Data()); + } + fTrackHistNames[icut] = names; + + std::unique_ptr objArrayCommon(commonNamesStr.Tokenize(",")); + for (int iCommonCut = 0; iCommonCut < fNCommonTrackCuts; ++iCommonCut) { + names = {}; + names.push_back(Form("TripletsBarrelSE_%s_%s", legsStr.Data(), objArrayCommon->At(iCommonCut)->GetName())); + histNames += Form("%s;", names[0].Data()); + fTrackHistNames[fNLegCuts + icut * fNCommonTrackCuts + iCommonCut] = names; + } + + TString cutNamesStr = fConfigPairCuts.value; + if (!cutNamesStr.IsNull()) { // if pair cuts + std::unique_ptr objArrayPair(cutNamesStr.Tokenize(",")); + fNPairCuts = objArrayPair->GetEntries(); + for (int iPairCut = 0; iPairCut < fNPairCuts; ++iPairCut) { // loop over pair cuts + names = {}; + names.push_back(Form("TripletsBarrelSE_%s_%s", legsStr.Data(), objArrayPair->At(iPairCut)->GetName())); + histNames += Form("%s;", names[0].Data()); + fTrackHistNames[fNLegCuts * (fNCommonTrackCuts + 1) + icut * fNPairCuts + iPairCut] = names; + for (int iCommonCut = 0; iCommonCut < fNCommonTrackCuts; ++iCommonCut) { + names = {}; + names.push_back(Form("TripletsBarrelSE_%s_%s_%s", legsStr.Data(), objArrayCommon->At(iCommonCut)->GetName(), objArrayPair->At(iPairCut)->GetName())); + histNames += Form("%s;", names[0].Data()); + fTrackHistNames[(fNLegCuts * (fNCommonTrackCuts + 1) + fNLegCuts * fNPairCuts) + icut * (fNPairCuts * fNCommonTrackCuts + 1) + iCommonCut * (1 + fNPairCuts) + iPairCut] = names; + } // end loop (common cuts) + } // end loop (pair cuts) + } // end if (pair cuts) + + // TODO: assign hist directories for the MC matched triplets for each (leg cut combo,MCsignal) combination + if (!sigNamesStr.IsNull()) { + for (unsigned int isig = 0; isig < fRecMCSignals.size(); isig++) { + auto sig = fRecMCSignals.at(isig); + int offset = fNLegCuts * isig * (1 + fNCommonTrackCuts + fNPairCuts + fNCommonTrackCuts * fNPairCuts); + names = { + Form("TripletsBarrelSE_%s_%s", legsStr.Data(), sig.GetName())}; + histNames += Form("%s;", names[0].Data()); + fBarrelHistNamesMCmatched[offset + icut] = names; + + for (int iCommonCut = 0; iCommonCut < fNCommonTrackCuts; ++iCommonCut) { + names = {}; + names.push_back(Form("TripletsBarrelSE_%s_%s_%s", legsStr.Data(), objArrayCommon->At(iCommonCut)->GetName(), sig.GetName())); + histNames += Form("%s;", names[0].Data()); + fBarrelHistNamesMCmatched[offset + fNLegCuts + icut * fNCommonTrackCuts + iCommonCut] = names; + } + + if (!cutNamesStr.IsNull()) { // if pair cuts + std::unique_ptr objArrayPair(cutNamesStr.Tokenize(",")); + for (int iPairCut = 0; iPairCut < fNPairCuts; ++iPairCut) { // loop over pair cuts + names = {}; + names.push_back(Form("TripletsBarrelSE_%s_%s_%s", legsStr.Data(), objArrayPair->At(iPairCut)->GetName(), sig.GetName())); + histNames += Form("%s;", names[0].Data()); + fBarrelHistNamesMCmatched[offset + fNLegCuts * (fNCommonTrackCuts + 1) + icut * fNPairCuts + iPairCut] = names; + for (int iCommonCut = 0; iCommonCut < fNCommonTrackCuts; ++iCommonCut) { + names = {}; + names.push_back(Form("TripletsBarrelSE_%s_%s_%s_%s", legsStr.Data(), objArrayCommon->At(iCommonCut)->GetName(), objArrayPair->At(iPairCut)->GetName(), sig.GetName())); + histNames += Form("%s;", names[0].Data()); + fBarrelHistNamesMCmatched[offset + (fNLegCuts * (fNCommonTrackCuts + 1) + fNLegCuts * fNPairCuts) + icut * (fNPairCuts * fNCommonTrackCuts) + iCommonCut * (1 + fNPairCuts) + iPairCut] = names; + } // end loop (common cuts) + } // end loop (pair cuts) + } // end if (pair cuts) + } // end loop over MC signals + } // end if (MC signals) + } else { + names = {}; + std::vector pairHistPrefixes = {"PairsBarrelSEPM"}; + if (fConfigSameSignHistograms.value) { + pairHistPrefixes.push_back("PairsBarrelSEPP"); + pairHistPrefixes.push_back("PairsBarrelSEMM"); + } + int fNPairHistPrefixes = pairHistPrefixes.size(); + + for (int iPrefix = 0; iPrefix < fNPairHistPrefixes; ++iPrefix) { + names.push_back(Form("%s_%s", pairHistPrefixes[iPrefix].Data(), legsStr.Data())); + histNames += Form("%s;", names[iPrefix].Data()); + } + if (fConfigQA) { + for (int iPrefix = 0; iPrefix < fNPairHistPrefixes; ++iPrefix) { + names.push_back(Form("%s_ambiguous_%s", pairHistPrefixes[iPrefix].Data(), legsStr.Data())); + histNames += Form("%s;", names[fNPairHistPrefixes + iPrefix].Data()); + } + } + fTrackHistNames[icut] = names; + + std::unique_ptr objArrayCommon(commonNamesStr.Tokenize(",")); + for (int iCommonCut = 0; iCommonCut < fNCommonTrackCuts; ++iCommonCut) { + names = {}; + for (int iPrefix = 0; iPrefix < fNPairHistPrefixes; ++iPrefix) { + names.push_back(Form("%s_%s_%s", pairHistPrefixes[iPrefix].Data(), legsStr.Data(), objArrayCommon->At(iCommonCut)->GetName())); + histNames += Form("%s;", names[iPrefix].Data()); + } + fTrackHistNames[fNLegCuts + icut * fNCommonTrackCuts + iCommonCut] = names; + } + + if (!cutNamesStr.IsNull()) { // if pair cuts + std::unique_ptr objArrayPair(cutNamesStr.Tokenize(",")); + fNPairCuts = objArrayPair->GetEntries(); + for (int iPairCut = 0; iPairCut < fNPairCuts; ++iPairCut) { // loop over pair cuts + names = {}; + for (int iPrefix = 0; iPrefix < fNPairHistPrefixes; ++iPrefix) { + names.push_back(Form("%s_%s_%s", pairHistPrefixes[iPrefix].Data(), legsStr.Data(), objArrayPair->At(iPairCut)->GetName())); + histNames += Form("%s;", names[iPrefix].Data()); + } + fTrackHistNames[fNLegCuts * (fNCommonTrackCuts + 1) + icut * fNPairCuts + iPairCut] = names; + for (int iCommonCut = 0; iCommonCut < fNCommonTrackCuts; ++iCommonCut) { + names = {}; + for (int iPrefix = 0; iPrefix < fNPairHistPrefixes; ++iPrefix) { + names.push_back(Form("%s_%s_%s_%s", pairHistPrefixes[iPrefix].Data(), legsStr.Data(), objArrayCommon->At(iCommonCut)->GetName(), objArrayPair->At(iPairCut)->GetName())); + histNames += Form("%s;", names[iPrefix].Data()); + } + fTrackHistNames[(fNLegCuts * (fNCommonTrackCuts + 1) + fNLegCuts * fNPairCuts) + icut * (fNPairCuts * fNCommonTrackCuts) + iCommonCut * (1 + fNPairCuts) + iPairCut] = names; + } // end loop (common cuts) + } // end loop (pair cuts) + } // end if (pair cuts) + + // assign hist directories for the MC matched triplets for each (leg cut combo,MCsignal) combination + if (!sigNamesStr.IsNull()) { + for (unsigned int isig = 0; isig < fRecMCSignals.size(); isig++) { + auto sig = fRecMCSignals.at(isig); + names = {}; + int offset = fNLegCuts * isig * (1 + fNCommonTrackCuts + fNPairCuts + fNCommonTrackCuts * fNPairCuts); + for (int iPrefix = 0; iPrefix < fNPairHistPrefixes; ++iPrefix) { + names.push_back(Form("%s_%s_%s", pairHistPrefixes[iPrefix].Data(), legsStr.Data(), sig.GetName())); + histNames += Form("%s;", names[iPrefix].Data()); + } + fBarrelHistNamesMCmatched[offset + icut] = names; + + for (int iCommonCut = 0; iCommonCut < fNCommonTrackCuts; ++iCommonCut) { + names = {}; + for (int iPrefix = 0; iPrefix < fNPairHistPrefixes; ++iPrefix) { + names.push_back(Form("%s_%s_%s_%s", pairHistPrefixes[iPrefix].Data(), legsStr.Data(), objArrayCommon->At(iCommonCut)->GetName(), sig.GetName())); + histNames += Form("%s;", names[iPrefix].Data()); + } + fBarrelHistNamesMCmatched[offset + fNLegCuts + icut * fNCommonTrackCuts + iCommonCut] = names; + } + + if (!cutNamesStr.IsNull()) { // if pair cuts + std::unique_ptr objArrayPair(cutNamesStr.Tokenize(",")); + for (int iPairCut = 0; iPairCut < fNPairCuts; ++iPairCut) { // loop over pair cuts + names = {}; + for (int iPrefix = 0; iPrefix < fNPairHistPrefixes; ++iPrefix) { + names.push_back(Form("%s_%s_%s_%s", pairHistPrefixes[iPrefix].Data(), legsStr.Data(), objArrayPair->At(iPairCut)->GetName(), sig.GetName())); + histNames += Form("%s;", names[iPrefix].Data()); + } + fBarrelHistNamesMCmatched[offset + fNLegCuts * (fNCommonTrackCuts + 1) + icut * fNPairCuts + iPairCut] = names; + for (int iCommonCut = 0; iCommonCut < fNCommonTrackCuts; ++iCommonCut) { + names = {}; + for (int iPrefix = 0; iPrefix < fNPairHistPrefixes; ++iPrefix) { + names.push_back(Form("%s_%s_%s_%s_%s", pairHistPrefixes[iPrefix].Data(), legsStr.Data(), objArrayCommon->At(iCommonCut)->GetName(), objArrayPair->At(iPairCut)->GetName(), sig.GetName())); + histNames += Form("%s;", names[iPrefix].Data()); + } + fBarrelHistNamesMCmatched[offset + (fNLegCuts * (fNCommonTrackCuts + 1) + fNLegCuts * fNPairCuts) + icut * (fNPairCuts * fNCommonTrackCuts) + iCommonCut * (1 + fNPairCuts) + iPairCut] = names; + } // end loop (common cuts) + } // end loop (pair cuts) + } // end if (pair cuts) + } // end loop over MC signals + } // end if (MC signals) + } + } + + // Add histogram classes for each specified MCsignal at the generator level + // TODO: create a std::vector of hist classes to be used at Fill time, to avoid using Form in the process function + TString sigGenNamesStr = fConfigMCGenSignals.value; + std::unique_ptr objGenSigArray(sigGenNamesStr.Tokenize(",")); + for (int isig = 0; isig < objGenSigArray->GetEntries(); isig++) { + MCSignal* sig = o2::aod::dqmcsignals::GetMCSignal(objGenSigArray->At(isig)->GetName()); + if (sig) { + if (sig->GetNProngs() == 1) { // NOTE: 1-prong signals required + fGenMCSignals.push_back(*sig); + histNames += Form("MCTruthGen_%s;", sig->GetName()); // TODO: Add these names to a std::vector to avoid using Form in the process function + } else if (sig->GetNProngs() == 2) { // NOTE: 2-prong signals required + fGenMCSignals.push_back(*sig); + histNames += Form("MCTruthGenPair_%s;", sig->GetName()); + fHasTwoProngGenMCsignals = true; + } else if (sig->GetNProngs() == 3) { // NOTE: 3-prong signals required + fGenMCSignals.push_back(*sig); + histNames += Form("MCTruthGenTriplet_%s;", sig->GetName()); + fHasThreeProngGenMCsignals = true; + } + } + } + + // Make sure the leg cuts are covered by the configured filter masks + if (fLegAFilterMask != fConstructedLegAFilterMask) { + LOGF(fatal, "cfgLegAFilterMask (%d) is not equal to the mask constructed by the cuts specified in cfgLegCuts (%d)!", fLegAFilterMask, fConstructedLegAFilterMask); + } + if (fLegBFilterMask != fConstructedLegBFilterMask) { + LOGF(fatal, "cfgLegBFilterMask (%d) is not equal to the mask constructed by the cuts specified in cfgLegCuts (%d)!", fLegBFilterMask, fConstructedLegBFilterMask); + } + if (fLegCFilterMask != fConstructedLegCFilterMask) { + LOGF(fatal, "cfgLegCFilterMask (%d) is not equal to the mask constructed by the cuts specified in cfgLegCuts (%d)!", fLegCFilterMask, fConstructedLegCFilterMask); + } + // Make sure only pairs or only triplets of leg cuts were given + int tripletCheckSum = std::count(isThreeProng.begin(), isThreeProng.end(), true); + if (tripletCheckSum != 0 && tripletCheckSum != fNLegCuts) { + LOGF(fatal, "A mix of pairs and triplets was given as leg cuts. Check your config!"); + } + + fCurrentRun = 0; + + fCCDB->setURL(fConfigCcdbUrl.value); + fCCDB->setCaching(true); + fCCDB->setLocalObjectValidityChecking(); + + fLUT = o2::base::MatLayerCylSet::rectifyPtrFromFile(fCCDB->get(fConfigLutPath)); + VarManager::SetupMatLUTFwdDCAFitter(fLUT); + + VarManager::SetDefaultVarNames(); + fHistMan = new HistogramManager("analysisHistos", "aa", VarManager::kNVars); + fHistMan->SetUseDefaultVariableNames(kTRUE); + fHistMan->SetDefaultVarNames(VarManager::fgVariableNames, VarManager::fgVariableUnits); + + DefineHistograms(fHistMan, histNames.Data(), fConfigHistogramSubgroups.value.data()); // define all histograms + VarManager::SetUseVars(fHistMan->GetUsedVars()); // provide the list of required variables so that VarManager knows what to fill + fOutputList.setObject(fHistMan->GetMainHistogramList()); + } + + void initParamsFromCCDB(uint64_t timestamp, bool isTriplets) + { + if (fConfigUseRemoteField.value) { + o2::parameters::GRPMagField* grpmag = fCCDB->getForTimeStamp(fConfigGRPMagPath, timestamp); + float magField = 0.0; + if (grpmag != nullptr) { + magField = grpmag->getNominalL3Field(); + } else { + LOGF(fatal, "GRP object is not available in CCDB at timestamp=%llu", timestamp); + } + if (isTriplets) { + if (fConfigUseKFVertexing.value) { + VarManager::SetupThreeProngKFParticle(magField); + } else { + VarManager::SetupThreeProngDCAFitter(magField, true, 200.0f, 4.0f, 1.0e-3f, 0.9f, fConfigUseAbsDCA.value); + } + } else { if (fConfigUseKFVertexing.value) { + VarManager::SetupTwoProngKFParticle(magField); + } else { + VarManager::SetupTwoProngDCAFitter(magField, true, 200.0f, 4.0f, 1.0e-3f, 0.9f, fConfigUseAbsDCA.value); // TODO: get these parameters from Configurables + } + } + } else { + if (isTriplets) { + if (fConfigUseKFVertexing.value) { + VarManager::SetupThreeProngKFParticle(fConfigMagField.value); + } else { + VarManager::SetupThreeProngDCAFitter(fConfigMagField.value, true, 200.0f, 4.0f, 1.0e-3f, 0.9f, fConfigUseAbsDCA.value); + } + } else { + if (fConfigUseKFVertexing.value) { + VarManager::SetupTwoProngKFParticle(fConfigMagField.value); + } else { + VarManager::SetupTwoProngDCAFitter(fConfigMagField.value, true, 200.0f, 4.0f, 1.0e-3f, 0.9f, fConfigUseAbsDCA.value); // TODO: get these parameters from Configurables + } + } + } + } + + // Template function to run same event pairing with asymmetric pairs (e.g. kaon-pion) + template + void runAsymmetricPairing(TEvents const& events, Preslice& preslice, TTrackAssocs const& /*assocs*/, TTracks const& /*tracks*/, ReducedMCEvents const& /*mcEvents*/, ReducedMCTracks const& /*mcTracks*/) + { + if (events.size() > 0) { // Additional protection to avoid crashing of events.begin().runNumber() + if (fCurrentRun != events.begin().runNumber()) { + initParamsFromCCDB(events.begin().timestamp(), false); + fCurrentRun = events.begin().runNumber(); + } + } + + std::map> histNamesMC = fBarrelHistNamesMCmatched; + std::map> histNames = fTrackHistNames; + + int sign1 = 0; + int sign2 = 0; + uint32_t mcDecision = 0; + ditrackList.reserve(1); + ditrackExtraList.reserve(1); + + constexpr bool trackHasCov = ((TTrackFillMap & VarManager::ObjTypes::TrackCov) > 0 || (TTrackFillMap & VarManager::ObjTypes::ReducedTrackBarrelCov) > 0); + + for (auto& event : events) { + if (!event.isEventSelected_bit(0)) { + continue; + } + // Reset the fValues array + VarManager::ResetValues(0, VarManager::kNVars); + VarManager::FillEvent(event, VarManager::fgValues); + + auto groupedLegAAssocs = legACandidateAssocs.sliceBy(preslice, event.globalIndex()); + if (groupedLegAAssocs.size() == 0) { + continue; + } + auto groupedLegBAssocs = legBCandidateAssocs.sliceBy(preslice, event.globalIndex()); + if (groupedLegBAssocs.size() == 0) { + continue; + } + + // TODO: Think about double counting + std::set> globIdxPairs; + for (auto& [a1, a2] : combinations(soa::CombinationsFullIndexPolicy(groupedLegAAssocs, groupedLegBAssocs))) { + + uint32_t twoTrackFilter = 0; + uint32_t twoTrackCommonFilter = 0; + uint32_t pairFilter = 0; + bool isPairIdWrong = false; + for (int icut = 0; icut < fNLegCuts; ++icut) { + // Find leg pair definitions both candidates participate in + if ((((a1.isBarrelSelected_raw() & fLegAFilterMask) | (a2.isBarrelSelected_raw() & fLegBFilterMask)) & fTrackCutFilterMasks[icut]) == fTrackCutFilterMasks[icut]) { + twoTrackFilter |= static_cast(1) << icut; + // If the supposed pion passes a kaon cut, this is a K+K-. Skip it. + if (TPairType == VarManager::kDecayToKPi && fConfigSkipAmbiguousIdCombinations.value) { + if (a2.isBarrelSelected_raw() & fLegAFilterMask) { + isPairIdWrong = true; + } + } + } + } + + if (!twoTrackFilter || isPairIdWrong) { + continue; + } + + // Find common track cuts both candidates pass + twoTrackCommonFilter |= a1.isBarrelSelected_raw() & a2.isBarrelSelected_raw() & fCommonTrackCutMask; + + auto t1 = a1.template reducedtrack_as(); + auto t2 = a2.template reducedtrack_as(); + + // Avoid self-pairs + if (t1.globalIndex() == t2.globalIndex()) { + continue; + } + + sign1 = t1.sign(); + sign2 = t2.sign(); + // store the ambiguity number of the two dilepton legs in the last 4 digits of the two-track filter + if (t1.barrelAmbiguityInBunch() > 1 || t1.barrelAmbiguityOutOfBunch() > 1) { + twoTrackFilter |= static_cast(1) << 30; + } + if (t2.barrelAmbiguityInBunch() > 1 || t2.barrelAmbiguityOutOfBunch() > 1) { + twoTrackFilter |= static_cast(1) << 31; + } + + // run MC matching for this pair + int isig = 0; + mcDecision = 0; + for (auto sig = fRecMCSignals.begin(); sig != fRecMCSignals.end(); sig++, isig++) { + if (t1.has_reducedMCTrack() && t2.has_reducedMCTrack()) { + if ((*sig).CheckSignal(true, t1.reducedMCTrack(), t2.reducedMCTrack())) { + mcDecision |= static_cast(1) << isig; + } + } + } // end loop over MC signals + + VarManager::FillPair(t1, t2); + if constexpr (TTwoProngFitter) { + VarManager::FillPairVertexing(event, t1, t2, fConfigPropToPCA); + } + + // Fill histograms + bool isAmbi = false; + for (int icut = 0; icut < fNLegCuts; icut++) { + if (twoTrackFilter & (static_cast(1) << icut)) { + isAmbi = (twoTrackFilter & (static_cast(1) << 30)) || (twoTrackFilter & (static_cast(1) << 31)); + if (sign1 * sign2 < 0) { // +- pairs + fHistMan->FillHistClass(histNames[icut][0].Data(), VarManager::fgValues); // reconstructed, unmatched + if (isAmbi && fConfigQA) { + fHistMan->FillHistClass(histNames[icut][3].Data(), VarManager::fgValues); + } + } else if (fConfigSameSignHistograms.value) { + if (sign1 > 0) { // ++ pairs + fHistMan->FillHistClass(histNames[icut][1].Data(), VarManager::fgValues); + if (isAmbi && fConfigQA) { + fHistMan->FillHistClass(histNames[icut][4].Data(), VarManager::fgValues); + } + } else { // -- pairs + fHistMan->FillHistClass(histNames[icut][2].Data(), VarManager::fgValues); + if (isAmbi && fConfigQA) { + fHistMan->FillHistClass(histNames[icut][5].Data(), VarManager::fgValues); + } + } + } + for (unsigned int isig = 0; isig < fRecMCSignals.size(); isig++) { // loop over MC signals + int offset = fNLegCuts * isig * (1 + fNCommonTrackCuts + fNPairCuts + fNCommonTrackCuts * fNPairCuts); + if (mcDecision & (static_cast(1) << isig)) { + if (sign1 * sign2 < 0) { + fHistMan->FillHistClass(histNamesMC[offset + icut][0].Data(), VarManager::fgValues); + } else if (fConfigSameSignHistograms.value) { + if (sign1 > 0) { + fHistMan->FillHistClass(histNamesMC[offset + icut][1].Data(), VarManager::fgValues); + } else { + fHistMan->FillHistClass(histNamesMC[offset + icut][2].Data(), VarManager::fgValues); + } + } + } + } + for (int iCommonCut = 0; iCommonCut < fNCommonTrackCuts; iCommonCut++) { + if (twoTrackCommonFilter & fCommonTrackCutFilterMasks[iCommonCut]) { + if (sign1 * sign2 < 0) { + fHistMan->FillHistClass(histNames[fNLegCuts + icut * fNCommonTrackCuts + iCommonCut][0].Data(), VarManager::fgValues); + } else if (fConfigSameSignHistograms.value) { + if (sign1 > 0) { + fHistMan->FillHistClass(histNames[fNLegCuts + icut * fNCommonTrackCuts + iCommonCut][1].Data(), VarManager::fgValues); + } else { + fHistMan->FillHistClass(histNames[fNLegCuts + icut * fNCommonTrackCuts + iCommonCut][2].Data(), VarManager::fgValues); + } + } + for (unsigned int isig = 0; isig < fRecMCSignals.size(); isig++) { // loop over MC signals + int offset = fNLegCuts * isig * (1 + fNCommonTrackCuts + fNPairCuts + fNCommonTrackCuts * fNPairCuts); + if (mcDecision & (static_cast(1) << isig)) { + if (sign1 * sign2 < 0) { + fHistMan->FillHistClass(histNamesMC[offset + fNLegCuts + icut * fNCommonTrackCuts + iCommonCut][0].Data(), VarManager::fgValues); + } else if (fConfigSameSignHistograms.value) { + if (sign1 > 0) { + fHistMan->FillHistClass(histNamesMC[offset + fNLegCuts + icut * fNCommonTrackCuts + iCommonCut][1].Data(), VarManager::fgValues); + } else { + fHistMan->FillHistClass(histNamesMC[offset + fNLegCuts + icut * fNCommonTrackCuts + iCommonCut][2].Data(), VarManager::fgValues); + } + } + } + } + } + } // end loop (common cuts) + for (unsigned int iPairCut = 0; iPairCut < fPairCuts.size(); iPairCut++) { + AnalysisCompositeCut cut = fPairCuts.at(iPairCut); + if (!(cut.IsSelected(VarManager::fgValues))) // apply pair cuts + continue; + pairFilter |= (static_cast(1) << iPairCut); + // Histograms with pair cuts + if (sign1 * sign2 < 0) { + fHistMan->FillHistClass(histNames[fNLegCuts * (fNCommonTrackCuts + 1) + icut * fNPairCuts + iPairCut][0].Data(), VarManager::fgValues); + } else if (fConfigSameSignHistograms.value) { + if (sign1 > 0) { + fHistMan->FillHistClass(histNames[fNLegCuts * (fNCommonTrackCuts + 1) + icut * fNPairCuts + iPairCut][1].Data(), VarManager::fgValues); + } else { + fHistMan->FillHistClass(histNames[fNLegCuts * (fNCommonTrackCuts + 1) + icut * fNPairCuts + iPairCut][2].Data(), VarManager::fgValues); + } + } + for (unsigned int isig = 0; isig < fRecMCSignals.size(); isig++) { // loop over MC signals + int offset = fNLegCuts * isig * (1 + fNCommonTrackCuts + fNPairCuts + fNCommonTrackCuts * fNPairCuts); + if (mcDecision & (static_cast(1) << isig)) { + if (sign1 * sign2 < 0) { + fHistMan->FillHistClass(histNamesMC[offset + fNLegCuts * (fNCommonTrackCuts + 1) + icut * fNPairCuts + iPairCut][0].Data(), VarManager::fgValues); + } else if (fConfigSameSignHistograms.value) { + if (sign1 > 0) { + fHistMan->FillHistClass(histNamesMC[offset + fNLegCuts * (fNCommonTrackCuts + 1) + icut * fNPairCuts + iPairCut][1].Data(), VarManager::fgValues); + } else { + fHistMan->FillHistClass(histNamesMC[offset + fNLegCuts * (fNCommonTrackCuts + 1) + icut * fNPairCuts + iPairCut][2].Data(), VarManager::fgValues); + } + } + } + } + // Histograms with pair cuts and common track cuts + for (int iCommonCut = 0; iCommonCut < fNCommonTrackCuts; ++iCommonCut) { + if (twoTrackCommonFilter & fCommonTrackCutFilterMasks[iCommonCut]) { + if (sign1 * sign2 < 0) { + fHistMan->FillHistClass(histNames[(fNLegCuts * (fNCommonTrackCuts + 1) + fNLegCuts * fNPairCuts) + icut * (fNPairCuts * fNCommonTrackCuts) + iCommonCut * (1 + fNPairCuts) + iPairCut][0].Data(), VarManager::fgValues); + } else if (fConfigSameSignHistograms.value) { + if (sign1 > 0) { + fHistMan->FillHistClass(histNames[(fNLegCuts * (fNCommonTrackCuts + 1) + fNLegCuts * fNPairCuts) + icut * (fNPairCuts * fNCommonTrackCuts) + iCommonCut * (1 + fNPairCuts) + iPairCut][1].Data(), VarManager::fgValues); + } else { + fHistMan->FillHistClass(histNames[(fNLegCuts * (fNCommonTrackCuts + 1) + fNLegCuts * fNPairCuts) + icut * (fNPairCuts * fNCommonTrackCuts) + iCommonCut * (1 + fNPairCuts) + iPairCut][2].Data(), VarManager::fgValues); + } + } + for (unsigned int isig = 0; isig < fRecMCSignals.size(); isig++) { // loop over MC signals + int offset = fNLegCuts * isig * (1 + fNCommonTrackCuts + fNPairCuts + fNCommonTrackCuts * fNPairCuts); + if (mcDecision & (static_cast(1) << isig)) { + if (sign1 * sign2 < 0) { + fHistMan->FillHistClass(histNamesMC[offset + (fNLegCuts * (fNCommonTrackCuts + 1) + fNLegCuts * fNPairCuts) + icut * (fNPairCuts * fNCommonTrackCuts) + iCommonCut * (1 + fNPairCuts) + iPairCut][0].Data(), VarManager::fgValues); + } else if (fConfigSameSignHistograms.value) { + if (sign1 > 0) { + fHistMan->FillHistClass(histNamesMC[offset + (fNLegCuts * (fNCommonTrackCuts + 1) + fNLegCuts * fNPairCuts) + icut * (fNPairCuts * fNCommonTrackCuts) + iCommonCut * (1 + fNPairCuts) + iPairCut][1].Data(), VarManager::fgValues); + } else { + fHistMan->FillHistClass(histNamesMC[offset + (fNLegCuts * (fNCommonTrackCuts + 1) + fNLegCuts * fNPairCuts) + icut * (fNPairCuts * fNCommonTrackCuts) + iCommonCut * (1 + fNPairCuts) + iPairCut][2].Data(), VarManager::fgValues); + } + } + } + } + } + } + } // end loop (pair cuts) + } + } // end loop (cuts) + ditrackList(event.globalIndex(), VarManager::fgValues[VarManager::kMass], + VarManager::fgValues[VarManager::kPt], VarManager::fgValues[VarManager::kEta], VarManager::fgValues[VarManager::kPhi], + t1.sign() + t2.sign(), twoTrackFilter, pairFilter, twoTrackCommonFilter); + if constexpr (trackHasCov && TTwoProngFitter) { + ditrackExtraList(t1.globalIndex(), t2.globalIndex(), VarManager::fgValues[VarManager::kVertexingTauzProjected], VarManager::fgValues[VarManager::kVertexingLzProjected], VarManager::fgValues[VarManager::kVertexingLxyProjected]); + } + } // end inner assoc loop (leg A) + } // end event loop + } + + // Template function to run same event triplets (e.g. D+->K-pi+pi+) + template + void runThreeProng(TEvents const& events, Preslice& preslice, TTrackAssocs const& /*assocs*/, TTracks const& tracks, ReducedMCEvents const& /*mcEvents*/, ReducedMCTracks const& /*mcTracks*/, VarManager::PairCandidateType tripletType) + { + if (events.size() > 0) { // Additional protection to avoid crashing of events.begin().runNumber() + if (fCurrentRun != events.begin().runNumber()) { + initParamsFromCCDB(events.begin().timestamp(), true); + fCurrentRun = events.begin().runNumber(); + } + } + + std::map> histNames = fTrackHistNames; + std::map> histNamesMC = fBarrelHistNamesMCmatched; + + for (auto& event : events) { + if (!event.isEventSelected_bit(0)) { + continue; + } + // Reset the fValues array + VarManager::ResetValues(0, VarManager::kNVars); + VarManager::FillEvent(event, VarManager::fgValues); + + auto groupedLegAAssocs = legACandidateAssocs.sliceBy(preslice, event.globalIndex()); + if (groupedLegAAssocs.size() == 0) { + continue; + } + auto groupedLegBAssocs = legBCandidateAssocs.sliceBy(preslice, event.globalIndex()); + if (groupedLegBAssocs.size() == 0) { + continue; + } + auto groupedLegCAssocs = legCCandidateAssocs.sliceBy(preslice, event.globalIndex()); + if (groupedLegCAssocs.size() == 0) { + continue; + } + + std::set> globIdxTriplets; + // Based on triplet type, make suitable combinations of the partitions + if (tripletType == VarManager::kTripleCandidateToPKPi) { + for (auto& [a1, a2, a3] : combinations(soa::CombinationsFullIndexPolicy(groupedLegAAssocs, groupedLegBAssocs, groupedLegCAssocs))) { + readTriplet(a1, a2, a3, tracks, event, tripletType, histNames, histNamesMC); + } + } else if (tripletType == VarManager::kTripleCandidateToKPiPi) { + for (auto& a1 : groupedLegAAssocs) { + for (auto& [a2, a3] : combinations(groupedLegBAssocs, groupedLegCAssocs)) { + readTriplet(a1, a2, a3, tracks, event, tripletType, histNames, histNamesMC); + } + } + } else { + LOG(fatal) << "Given tripletType not recognized. Don't know how to make combinations!" << endl; + } + } // end event loop + } + + // Helper function to process triplet + template + void readTriplet(TTrackAssoc const& a1, TTrackAssoc const& a2, TTrackAssoc const& a3, TTracks const& /*tracks*/, TEvent const& event, VarManager::PairCandidateType tripletType, std::map> histNames, std::map> histNamesMC) + { + uint32_t mcDecision = 0; + + uint32_t threeTrackFilter = 0; + uint32_t threeTrackCommonFilter = 0; + for (int icut = 0; icut < fNLegCuts; ++icut) { + // Find out which leg cut combination the triplet passes + if ((((a1.isBarrelSelected_raw() & fLegAFilterMask) | (a2.isBarrelSelected_raw() & fLegBFilterMask) | (a3.isBarrelSelected_raw() & fLegCFilterMask)) & fTrackCutFilterMasks[icut]) == fTrackCutFilterMasks[icut]) { + threeTrackFilter |= (static_cast(1) << icut); + if (tripletType == VarManager::kTripleCandidateToPKPi && fConfigSkipAmbiguousIdCombinations.value) { + // Check if the supposed pion passes as a proton or kaon, if so, skip this triplet. It is pKp or pKK. + if ((a3.isBarrelSelected_raw() & fLegAFilterMask) || (a3.isBarrelSelected_raw() & fLegBFilterMask)) { + return; + } + // Check if the supposed kaon passes as a proton, if so, skip this triplet. It is ppPi. + if (a2.isBarrelSelected_raw() & fLegAFilterMask) { + return; + } + } + if (tripletType == VarManager::kTripleCandidateToKPiPi && fConfigSkipAmbiguousIdCombinations.value) { + // Check if one of the supposed pions pass as a kaon, if so, skip this triplet. It is KKPi. + if ((a2.isBarrelSelected_raw() & fLegAFilterMask) || (a3.isBarrelSelected_raw() & fLegAFilterMask)) { + return; + } + } + } + } + if (!threeTrackFilter) { + return; + } + + // Find common track cuts all candidates pass + threeTrackCommonFilter |= a1.isBarrelSelected_raw() & a2.isBarrelSelected_raw() & a3.isBarrelSelected_raw() & fCommonTrackCutMask; + + auto t1 = a1.template reducedtrack_as(); + auto t2 = a2.template reducedtrack_as(); + auto t3 = a3.template reducedtrack_as(); + + // Avoid self-pairs + if (t1 == t2 || t1 == t3 || t2 == t3) { + return; + } + + // store the ambiguity of the three legs in the last 3 digits of the two-track filter + if (t1.barrelAmbiguityInBunch() > 1 || t1.barrelAmbiguityOutOfBunch() > 1) { + threeTrackFilter |= (static_cast(1) << 29); + } + if (t2.barrelAmbiguityInBunch() > 1 || t2.barrelAmbiguityOutOfBunch() > 1) { + threeTrackFilter |= (static_cast(1) << 30); + } + if (t3.barrelAmbiguityInBunch() > 1 || t3.barrelAmbiguityOutOfBunch() > 1) { + threeTrackFilter |= (static_cast(1) << 31); + } + + // run MC matching for this triplet + int isig = 0; + mcDecision = 0; + for (auto sig = fRecMCSignals.begin(); sig != fRecMCSignals.end(); sig++, isig++) { + if (t1.has_reducedMCTrack() && t2.has_reducedMCTrack() && t3.has_reducedMCTrack()) { + if ((*sig).CheckSignal(true, t1.reducedMCTrack(), t2.reducedMCTrack(), t3.reducedMCTrack())) { + mcDecision |= (static_cast(1) << isig); + } + } + } // end loop over MC signals + + VarManager::FillTriple(t1, t2, t3, VarManager::fgValues, tripletType); + if constexpr (TThreeProngFitter) { + VarManager::FillTripletVertexing(event, t1, t2, t3, tripletType); + } + + // Fill histograms + bool isAmbi = false; + for (int icut = 0; icut < fNLegCuts; icut++) { + isAmbi = (threeTrackFilter & (static_cast(1) << 29)) || (threeTrackFilter & (static_cast(1) << 30)) || (threeTrackFilter & (static_cast(1) << 31)); + if (threeTrackFilter & (static_cast(1) << icut)) { + fHistMan->FillHistClass(histNames[icut][0].Data(), VarManager::fgValues); + // TODO: loop over MC signals + for (unsigned int isig = 0; isig < fRecMCSignals.size(); isig++) { // loop over MC signals + int offset = fNLegCuts * isig * (1 + fNCommonTrackCuts + fNPairCuts + fNCommonTrackCuts * fNPairCuts); + if (mcDecision & (static_cast(1) << isig)) { + fHistMan->FillHistClass(histNamesMC[offset + icut][0].Data(), VarManager::fgValues); // matched signal + } + } // end loop (MC signals) + if (fConfigQA && isAmbi) { + fHistMan->FillHistClass(histNames[icut][1].Data(), VarManager::fgValues); + } + for (int iCommonCut = 0; iCommonCut < fNCommonTrackCuts; iCommonCut++) { + if (threeTrackCommonFilter & fCommonTrackCutFilterMasks[iCommonCut]) { + fHistMan->FillHistClass(histNames[fNLegCuts + icut * fNCommonTrackCuts + iCommonCut][0].Data(), VarManager::fgValues); + for (unsigned int isig = 0; isig < fRecMCSignals.size(); isig++) { // loop over MC signals + int offset = fNLegCuts * isig * (1 + fNCommonTrackCuts + fNPairCuts + fNCommonTrackCuts * fNPairCuts); + if (mcDecision & (static_cast(1) << isig)) { + fHistMan->FillHistClass(histNamesMC[offset + fNLegCuts + icut * fNCommonTrackCuts + iCommonCut][0].Data(), VarManager::fgValues); // matched signal + } + } // end loop (MC signals) + } + } // end loop (common cuts) + for (unsigned int iPairCut = 0; iPairCut < fPairCuts.size(); iPairCut++) { + AnalysisCompositeCut cut = fPairCuts.at(iPairCut); + if (!(cut.IsSelected(VarManager::fgValues))) { // apply pair cuts + continue; + } + // Histograms with pair cuts + fHistMan->FillHistClass(histNames[fNLegCuts * (fNCommonTrackCuts + 1) + icut * fNPairCuts + iPairCut][0].Data(), VarManager::fgValues); + for (unsigned int isig = 0; isig < fRecMCSignals.size(); isig++) { // loop over MC signals + int offset = fNLegCuts * isig * (1 + fNCommonTrackCuts + fNPairCuts + fNCommonTrackCuts * fNPairCuts); + if (mcDecision & (static_cast(1) << isig)) { + fHistMan->FillHistClass(histNamesMC[offset + fNLegCuts * (fNCommonTrackCuts + 1) + icut * fNPairCuts + iPairCut][0].Data(), VarManager::fgValues); // matched signal + } + } // end loop (MC signals) + // Histograms with pair cuts and common track cuts + for (int iCommonCut = 0; iCommonCut < fNCommonTrackCuts; ++iCommonCut) { + if (threeTrackCommonFilter & fCommonTrackCutFilterMasks[iCommonCut]) { + fHistMan->FillHistClass(histNames[(fNLegCuts * (fNCommonTrackCuts + 1) + fNLegCuts * fNPairCuts) + icut * (fNPairCuts * fNCommonTrackCuts + 1) + iCommonCut * (1 + fNPairCuts) + iPairCut][0].Data(), VarManager::fgValues); + for (unsigned int isig = 0; isig < fRecMCSignals.size(); isig++) { // loop over MC signals + int offset = fNLegCuts * isig * (1 + fNCommonTrackCuts + fNPairCuts + fNCommonTrackCuts * fNPairCuts); + if (mcDecision & (static_cast(1) << isig)) { + fHistMan->FillHistClass(histNamesMC[offset + (fNLegCuts * (fNCommonTrackCuts + 1) + fNLegCuts * fNPairCuts) + icut * (fNPairCuts * fNCommonTrackCuts) + iCommonCut * (1 + fNPairCuts) + iPairCut][0].Data(), VarManager::fgValues); // matched signal + } + } // end loop (MC signals) + } + } + } // end loop (pair cuts) + } + } // end loop (cuts) + } + + PresliceUnsorted perReducedMcEvent = aod::reducedtrackMC::reducedMCeventId; + + void runMCGen(ReducedMCEvents const& mcEvents, ReducedMCTracks const& mcTracks, VarManager::PairCandidateType pairType) + { + // loop over mc stack and fill histograms for pure MC truth signals + // group all the MC tracks which belong to the MC event corresponding to the current reconstructed event + for (auto& mctrack : mcTracks) { + VarManager::FillTrackMC(mcTracks, mctrack); + // NOTE: Signals are checked here mostly based on the skimmed MC stack, so depending on the requested signal, the stack could be incomplete. + // NOTE: However, the working model is that the decisions on MC signals are precomputed during skimming and are stored in the mcReducedFlags member. + // TODO: Use the mcReducedFlags to select signals + for (auto& sig : fGenMCSignals) { + if (sig.GetNProngs() != 1) { // NOTE: 1-prong signals required here + continue; + } + bool checked = false; + checked = sig.CheckSignal(true, mctrack); + if (checked) { + fHistMan->FillHistClass(Form("MCTruthGen_%s", sig.GetName()), VarManager::fgValues); + } + } + } + + if (fHasTwoProngGenMCsignals) { + for (auto& event : mcEvents) { + auto groupedMCTracks = mcTracks.sliceBy(perReducedMcEvent, event.globalIndex()); + groupedMCTracks.bindInternalIndicesTo(&mcTracks); + for (auto& [t1, t2] : combinations(groupedMCTracks, groupedMCTracks)) { + auto t1_raw = groupedMCTracks.rawIteratorAt(t1.globalIndex()); + auto t2_raw = groupedMCTracks.rawIteratorAt(t2.globalIndex()); + for (auto& sig : fGenMCSignals) { + if (sig.GetNProngs() != 2) { // NOTE: 2-prong signals required here + continue; + } + if (sig.CheckSignal(true, t1_raw, t2_raw)) { + VarManager::FillPairMC(t1, t2, VarManager::fgValues, pairType); + fHistMan->FillHistClass(Form("MCTruthGenPair_%s", sig.GetName()), VarManager::fgValues); + } + } // end loop over MC signals + } // end loop over pairs + } // end loop over events + } + + if (fHasThreeProngGenMCsignals) { + for (auto& event : mcEvents) { + auto groupedMCTracks = mcTracks.sliceBy(perReducedMcEvent, event.globalIndex()); + groupedMCTracks.bindInternalIndicesTo(&mcTracks); + for (auto& [t1, t2, t3] : combinations(groupedMCTracks, groupedMCTracks, groupedMCTracks)) { + auto t1_raw = groupedMCTracks.rawIteratorAt(t1.globalIndex()); + auto t2_raw = groupedMCTracks.rawIteratorAt(t2.globalIndex()); + auto t3_raw = groupedMCTracks.rawIteratorAt(t3.globalIndex()); + for (auto& sig : fGenMCSignals) { + if (sig.GetNProngs() != 3) { // NOTE: 3-prong signals required here + continue; + } + if (sig.CheckSignal(true, t1_raw, t2_raw, t3_raw)) { + VarManager::FillTripleMC(t1, t2, t3, VarManager::fgValues, pairType); + fHistMan->FillHistClass(Form("MCTruthGenTriplet_%s", sig.GetName()), VarManager::fgValues); + } + } // end loop over MC signals + } // end loop over pairs + } // end loop over events + } + } // end runMCGen + + void processKaonPionSkimmed(MyEventsVtxCovSelected const& events, + soa::Join const& barrelAssocs, + MyBarrelTracksWithCovWithAmbiguities const& barrelTracks, + ReducedMCEvents const& mcEvents, ReducedMCTracks const& mcTracks) + { + runAsymmetricPairing(events, trackAssocsPerCollision, barrelAssocs, barrelTracks, mcEvents, mcTracks); + if (fConfigRunMCGenPair) + runMCGen(mcEvents, mcTracks, VarManager::kDecayToKPi); + } + + void processKaonPionPionSkimmed(MyEventsVtxCovSelected const& events, + soa::Join const& barrelAssocs, + MyBarrelTracksWithCovWithAmbiguities const& barrelTracks, + ReducedMCEvents const& mcEvents, ReducedMCTracks const& mcTracks) + { + runThreeProng(events, trackAssocsPerCollision, barrelAssocs, barrelTracks, mcEvents, mcTracks, VarManager::kTripleCandidateToKPiPi); + if (fConfigRunMCGenPair) + runMCGen(mcEvents, mcTracks, VarManager::kTripleCandidateToKPiPi); + } + + void processDummy(MyEvents&) + { + // do nothing + } + + PROCESS_SWITCH(AnalysisAsymmetricPairing, processKaonPionSkimmed, "Run kaon pion pairing, with skimmed tracks", false); + PROCESS_SWITCH(AnalysisAsymmetricPairing, processKaonPionPionSkimmed, "Run kaon pion pion triplets, with skimmed tracks", false); + PROCESS_SWITCH(AnalysisAsymmetricPairing, processDummy, "Dummy function, enabled only if none of the others are enabled", true); +}; + // Combines dileptons with barrel or muon tracks for either resonance or correlation analyses // Dileptons produced with all the selection cuts specified in the same-event pairing task are combined with the // tracks passing the fConfigTrackCut cut. The dileptons cuts from the same-event pairing task are auto-detected @@ -1903,8 +2874,12 @@ struct AnalysisDileptonTrack { int fCurrentRun; // needed to detect if the run changed and trigger update of calibrations etc. int fNCuts; + int fNPairCuts; + int fNCommonTrackCuts; + std::map fCommonTrackCutMap; int fTrackCutBit; std::map fHistNamesDileptonTrack; + // std::map fHistNamesDileptonTrackMCmatched; std::map> fHistNamesDileptonTrackMCmatched; std::vector fHistNamesMCgen; std::map fHistNamesDileptons; @@ -1934,8 +2909,9 @@ struct AnalysisDileptonTrack { } bool isBarrel = context.mOptions.get("processBarrelSkimmed"); + bool isBarrelAsymmetric = context.mOptions.get("processDstarToD0Pi"); bool isMuon = context.mOptions.get("processMuonSkimmed"); - bool isAnyProcessEnabled = isBarrel || isMuon; + bool isAnyProcessEnabled = isBarrel || isMuon || isBarrelAsymmetric; bool isDummy = context.mOptions.get("processDummy"); if (isDummy) { if (isAnyProcessEnabled) { @@ -1972,10 +2948,10 @@ struct AnalysisDileptonTrack { // For each track/muon selection used to produce dileptons, create a separate histogram directory using the // name of the track/muon cut. // Also, create a map which will hold the name of the histogram directories so they can be accessed directly in the pairing loop - if (isBarrel || isMuon) { + if (isBarrel || isMuon || isBarrelAsymmetric) { // get the list of single track and muon cuts computed in the dedicated tasks upstream string tempCutsSingle; - if (isBarrel) { + if (isBarrel || isBarrelAsymmetric) { getTaskOptionValue(context, "analysis-track-selection", "cfgTrackCuts", tempCutsSingle, false); } else { getTaskOptionValue(context, "analysis-muon-selection", "cfgMuonCuts", tempCutsSingle, false); @@ -1988,6 +2964,7 @@ struct AnalysisDileptonTrack { if (objArraySingleCuts->FindObject(fConfigTrackCut.value.data()) == nullptr) { LOG(fatal) << " Track cut chosen for the correlation task was not computed in the single-track task! Check it out!"; } + // Loop over single-track/muon task cuts and find the cuts used for the track to be combined with dileptons for (int icut = 0; icut < objArraySingleCuts->GetEntries(); ++icut) { TString tempStr = objArraySingleCuts->At(icut)->GetName(); if (tempStr.CompareTo(fConfigTrackCut.value.data()) == 0) { @@ -1996,31 +2973,99 @@ struct AnalysisDileptonTrack { } // get the cuts employed for same-event pairing string tempCutsPair; + string tempCutsAsymPair; + string tempCutsAsymCommon; if (isBarrel) { getTaskOptionValue(context, "analysis-same-event-pairing", "cfgTrackCuts", tempCutsPair, false); - } else { + } else if (isMuon) { getTaskOptionValue(context, "analysis-same-event-pairing", "cfgMuonCuts", tempCutsPair, false); + } else if (isBarrelAsymmetric) { + getTaskOptionValue(context, "analysis-asymmetric-pairing", "cfgLegCuts", tempCutsPair, false); + getTaskOptionValue(context, "analysis-asymmetric-pairing", "cfgPairCuts", tempCutsAsymPair, false); + getTaskOptionValue(context, "analysis-asymmetric-pairing", "cfgCommonTrackCuts", tempCutsAsymCommon, false); + } + + // If asymmetric pair is used, it may have common track cuts + TString tempCutsAsymCommonStr = tempCutsAsymCommon; + if (!tempCutsAsymCommonStr.IsNull()) { // if common track cuts + std::unique_ptr objArrayCommon(tempCutsAsymCommonStr.Tokenize(",")); + fNCommonTrackCuts = objArrayCommon->GetEntries(); + for (int icut = 0; icut < fNCommonTrackCuts; ++icut) { + for (int iicut = 0; iicut < objArraySingleCuts->GetEntries(); ++iicut) { + if (std::strcmp(objArrayCommon->At(icut)->GetName(), objArraySingleCuts->At(iicut)->GetName()) == 0) { + fCommonTrackCutMap[icut] = iicut; + } + } + } } + TString tempCutsPairStr = tempCutsPair; if (!tempCutsSingleStr.IsNull() && !tempCutsPairStr.IsNull()) { std::unique_ptr objArray(tempCutsPairStr.Tokenize(",")); fNCuts = objArray->GetEntries(); - for (int icut = 0; icut < objArray->GetEntries(); ++icut) { + for (int icut = 0; icut < fNCuts; ++icut) { TString tempStr = objArray->At(icut)->GetName(); - if (objArray->FindObject(tempStr.Data()) != nullptr) { - fHistNamesDileptonTrack[icut] = Form("DileptonTrack_%s_%s", tempStr.Data(), fConfigTrackCut.value.data()); - fHistNamesDileptons[icut] = Form("DileptonsSelected_%s", tempStr.Data()); - DefineHistograms(fHistMan, fHistNamesDileptonTrack[icut], fConfigHistogramSubgroups.value.data()); // define dilepton-track histograms - DefineHistograms(fHistMan, fHistNamesDileptons[icut], "barrel,vertexing"); // define dilepton histograms - std::vector mcHistNames; - - for (auto& sig : fRecMCSignals) { - mcHistNames.push_back(Form("DileptonTrackMCMatched_%s_%s_%s", tempStr.Data(), fConfigTrackCut.value.data(), sig.GetName())); - DefineHistograms(fHistMan, mcHistNames[mcHistNames.size() - 1], fConfigHistogramSubgroups.value.data()); + fHistNamesDileptonTrack[icut] = Form("DileptonTrack_%s_%s", tempStr.Data(), fConfigTrackCut.value.data()); + fHistNamesDileptons[icut] = Form("DileptonsSelected_%s", tempStr.Data()); + DefineHistograms(fHistMan, fHistNamesDileptonTrack[icut], fConfigHistogramSubgroups.value.data()); // define dilepton-track histograms + DefineHistograms(fHistMan, fHistNamesDileptons[icut], "barrel,vertexing"); // define dilepton histograms + if (!tempCutsAsymCommonStr.IsNull()) { + std::unique_ptr objArrayCommon(tempCutsAsymCommonStr.Tokenize(",")); + for (int iCommonCut = 0; iCommonCut < fNCommonTrackCuts; ++iCommonCut) { + fHistNamesDileptonTrack[fNCuts + icut * fNCommonTrackCuts + iCommonCut] = Form("DileptonTrack_%s_%s_%s", tempStr.Data(), objArrayCommon->At(iCommonCut)->GetName(), fConfigTrackCut.value.data()); + fHistNamesDileptons[fNCuts + icut * fNCommonTrackCuts + iCommonCut] = Form("DileptonsSelected_%s_%s", tempStr.Data(), objArrayCommon->At(iCommonCut)->GetName()); + DefineHistograms(fHistMan, fHistNamesDileptonTrack[fNCuts + icut * fNCommonTrackCuts + iCommonCut], fConfigHistogramSubgroups.value.data()); // define dilepton-track histograms + DefineHistograms(fHistMan, fHistNamesDileptons[fNCuts + icut * fNCommonTrackCuts + iCommonCut], "barrel,vertexing"); // define dilepton histograms } - fHistNamesDileptonTrackMCmatched[icut] = mcHistNames; } - } + TString tempCutsAsymPairStr = tempCutsAsymPair; + if (!tempCutsAsymPairStr.IsNull()) { + std::unique_ptr objArrayPairCuts(tempCutsAsymPairStr.Tokenize(",")); + fNPairCuts = objArrayPairCuts->GetEntries(); + for (int iPairCut = 0; iPairCut < fNPairCuts; ++iPairCut) { + fHistNamesDileptonTrack[fNCuts * (fNCommonTrackCuts + 1) + icut * fNPairCuts + iPairCut] = Form("DileptonTrack_%s_%s_%s", tempStr.Data(), objArrayPairCuts->At(iPairCut)->GetName(), fConfigTrackCut.value.data()); + fHistNamesDileptons[fNCuts * (fNCommonTrackCuts + 1) + icut * fNPairCuts + iPairCut] = Form("DileptonsSelected_%s_%s", tempStr.Data(), objArrayPairCuts->At(iPairCut)->GetName()); + DefineHistograms(fHistMan, fHistNamesDileptonTrack[fNCuts * (fNCommonTrackCuts + 1) + icut * fNPairCuts + iPairCut], fConfigHistogramSubgroups.value.data()); // define dilepton-track histograms + DefineHistograms(fHistMan, fHistNamesDileptons[fNCuts * (fNCommonTrackCuts + 1) + icut * fNPairCuts + iPairCut], "barrel,vertexing"); // define dilepton histograms + if (!tempCutsAsymCommonStr.IsNull()) { + std::unique_ptr objArrayCommon(tempCutsAsymCommonStr.Tokenize(",")); + for (int iCommonCut = 0; iCommonCut < fNCommonTrackCuts; ++iCommonCut) { + fHistNamesDileptonTrack[(fNCuts * (fNCommonTrackCuts + 1) + fNCuts * fNPairCuts) + icut * (fNPairCuts * fNCommonTrackCuts) + iCommonCut * fNPairCuts + iPairCut] = Form("DileptonTrack_%s_%s_%s_%s", tempStr.Data(), objArrayCommon->At(iCommonCut)->GetName(), objArrayPairCuts->At(iPairCut)->GetName(), fConfigTrackCut.value.data()); + fHistNamesDileptons[(fNCuts * (fNCommonTrackCuts + 1) + fNCuts * fNPairCuts) + icut * (fNPairCuts * fNCommonTrackCuts) + iCommonCut * fNPairCuts + iPairCut] = Form("DileptonsSelected_%s_%s_%s", tempStr.Data(), objArrayCommon->At(iCommonCut)->GetName(), objArrayPairCuts->At(iPairCut)->GetName()); + DefineHistograms(fHistMan, fHistNamesDileptonTrack[(fNCuts * (fNCommonTrackCuts + 1) + fNCuts * fNPairCuts) + icut * (fNPairCuts * fNCommonTrackCuts) + iCommonCut * fNPairCuts + iPairCut], fConfigHistogramSubgroups.value.data()); // define dilepton-track histograms + DefineHistograms(fHistMan, fHistNamesDileptons[(fNCuts * (fNCommonTrackCuts + 1) + fNCuts * fNPairCuts) + icut * (fNPairCuts * fNCommonTrackCuts) + iCommonCut * fNPairCuts + iPairCut], "barrel,vertexing"); // define dilepton histograms + } + } + } // end loop (pair cuts) + } + for (unsigned int isig = 0; isig < fRecMCSignals.size(); isig++) { + auto sig = fRecMCSignals.at(isig); + fHistNamesDileptonTrackMCmatched[icut].push_back(Form("DileptonTrackMCMatched_%s_%s_%s", tempStr.Data(), fConfigTrackCut.value.data(), sig.GetName())); + DefineHistograms(fHistMan, fHistNamesDileptonTrackMCmatched[icut].back(), fConfigHistogramSubgroups.value.data()); + if (!tempCutsAsymCommonStr.IsNull()) { + std::unique_ptr objArrayCommon(tempCutsAsymCommonStr.Tokenize(",")); + for (int iCommonCut = 0; iCommonCut < fNCommonTrackCuts; ++iCommonCut) { + fHistNamesDileptonTrackMCmatched[fNCuts + icut * fNCommonTrackCuts + iCommonCut].push_back(Form("DileptonTrackMCMatched_%s_%s_%s_%s", tempStr.Data(), objArrayCommon->At(iCommonCut)->GetName(), fConfigTrackCut.value.data(), sig.GetName())); + DefineHistograms(fHistMan, fHistNamesDileptonTrackMCmatched[fNCuts + icut * fNCommonTrackCuts + iCommonCut].back(), fConfigHistogramSubgroups.value.data()); // define dilepton-track histograms + } + } + if (!tempCutsAsymPairStr.IsNull()) { + std::unique_ptr objArrayPairCuts(tempCutsAsymPairStr.Tokenize(",")); + fNPairCuts = objArrayPairCuts->GetEntries(); + for (int iPairCut = 0; iPairCut < fNPairCuts; ++iPairCut) { + fHistNamesDileptonTrackMCmatched[fNCuts * (fNCommonTrackCuts + 1) + icut * fNPairCuts + iPairCut].push_back(Form("DileptonTrackMCMatched_%s_%s_%s_%s", tempStr.Data(), objArrayPairCuts->At(iPairCut)->GetName(), fConfigTrackCut.value.data(), sig.GetName())); + DefineHistograms(fHistMan, fHistNamesDileptonTrackMCmatched[fNCuts * (fNCommonTrackCuts + 1) + icut * fNPairCuts + iPairCut].back(), fConfigHistogramSubgroups.value.data()); // define dilepton-track histograms + if (!tempCutsAsymCommonStr.IsNull()) { + std::unique_ptr objArrayCommon(tempCutsAsymCommonStr.Tokenize(",")); + for (int iCommonCut = 0; iCommonCut < fNCommonTrackCuts; ++iCommonCut) { + fHistNamesDileptonTrackMCmatched[(fNCuts * (fNCommonTrackCuts + 1) + fNCuts * fNPairCuts) + icut * (fNPairCuts * fNCommonTrackCuts) + iCommonCut * fNPairCuts + iPairCut].push_back(Form("DileptonTrackMCMatched_%s_%s_%s_%s_%s", tempStr.Data(), objArrayCommon->At(iCommonCut)->GetName(), objArrayPairCuts->At(iPairCut)->GetName(), fConfigTrackCut.value.data(), sig.GetName())); + DefineHistograms(fHistMan, fHistNamesDileptonTrackMCmatched[(fNCuts * (fNCommonTrackCuts + 1) + fNCuts * fNPairCuts) + icut * (fNPairCuts * fNCommonTrackCuts) + iCommonCut * fNPairCuts + iPairCut].back(), fConfigHistogramSubgroups.value.data()); // define dilepton-track histograms + } // end loop (common cuts) + } + } // end loop (pair cuts) + } + } // end loop (MC signals) + } // end loop (leg defining cuts) } // Add histogram classes for each specified MCsignal at the generator level // TODO: create a std::vector of hist classes to be used at Fill time, to avoid using Form in the process function @@ -2099,12 +3144,29 @@ struct AnalysisDileptonTrack { for (int icut = 0; icut < fNCuts; icut++) { if (dilepton.filterMap_bit(icut)) { fHistMan->FillHistClass(fHistNamesDileptons[icut].Data(), fValuesDilepton); + if constexpr (TCandidateType == VarManager::kDstarToD0KPiPi) { // Dielectrons and Dimuons don't have the PairFilterMap column + for (int iCommonCut = 0; iCommonCut < fNCommonTrackCuts; iCommonCut++) { + if (dilepton.commonFilterMap_bit(fCommonTrackCutMap[iCommonCut])) { + fHistMan->FillHistClass(fHistNamesDileptons[fNCuts + icut * fNCommonTrackCuts + iCommonCut].Data(), fValuesDilepton); + } + } + for (int iPairCut = 0; iPairCut < fNPairCuts; iPairCut++) { + if (dilepton.pairFilterMap_bit(iPairCut)) { + fHistMan->FillHistClass(fHistNamesDileptons[fNCuts * (fNCommonTrackCuts + 1) + icut * fNPairCuts + iPairCut].Data(), fValuesDilepton); + for (int iCommonCut = 0; iCommonCut < fNCommonTrackCuts; iCommonCut++) { + if (dilepton.commonFilterMap_bit(fCommonTrackCutMap[iCommonCut])) { + fHistMan->FillHistClass(fHistNamesDileptons[(fNCuts * (fNCommonTrackCuts + 1) + fNCuts * fNPairCuts) + icut * (fNPairCuts * fNCommonTrackCuts) + iCommonCut * fNPairCuts + iPairCut].Data(), fValuesDilepton); + } + } + } + } + } } } // loop over hadrons for (auto& assoc : assocs) { - if constexpr (TCandidateType == VarManager::kBtoJpsiEEK) { + if constexpr (TCandidateType == VarManager::kBtoJpsiEEK || TCandidateType == VarManager::kDstarToD0KPiPi) { if (!assoc.isBarrelSelected_bit(fTrackCutBit)) { continue; } @@ -2155,7 +3217,39 @@ struct AnalysisDileptonTrack { fHistMan->FillHistClass(fHistNamesDileptonTrackMCmatched[icut][isig], fValuesHadron); } } - } + if constexpr (TCandidateType == VarManager::kDstarToD0KPiPi) { // Dielectrons and Dimuons don't have the PairFilterMap column + for (int iCommonCut = 0; iCommonCut < fNCommonTrackCuts; iCommonCut++) { + if (dilepton.commonFilterMap_bit(fCommonTrackCutMap[iCommonCut])) { + fHistMan->FillHistClass(fHistNamesDileptonTrack[fNCuts + icut * fNCommonTrackCuts + iCommonCut].Data(), fValuesHadron); + for (isig = 0; isig < fRecMCSignals.size(); isig++) { + if (mcDecision & (static_cast(1) << isig)) { + fHistMan->FillHistClass(fHistNamesDileptonTrackMCmatched[fNCuts + icut * fNCommonTrackCuts + iCommonCut][isig], fValuesHadron); + } + } // end loop (MC signals) + } + } // end loop (common track cuts) + for (int iPairCut = 0; iPairCut < fNPairCuts; iPairCut++) { + if (dilepton.pairFilterMap_bit(iPairCut)) { + fHistMan->FillHistClass(fHistNamesDileptonTrack[fNCuts * (fNCommonTrackCuts + 1) + icut * fNPairCuts + iPairCut].Data(), fValuesHadron); + for (isig = 0; isig < fRecMCSignals.size(); isig++) { + if (mcDecision & (static_cast(1) << isig)) { + fHistMan->FillHistClass(fHistNamesDileptonTrackMCmatched[fNCuts * (fNCommonTrackCuts + 1) + icut * fNPairCuts + iPairCut][isig], fValuesHadron); + } + } + for (int iCommonCut = 0; iCommonCut < fNCommonTrackCuts; iCommonCut++) { + if (dilepton.commonFilterMap_bit(fCommonTrackCutMap[iCommonCut])) { + fHistMan->FillHistClass(fHistNamesDileptonTrack[(fNCuts * (fNCommonTrackCuts + 1) + fNCuts * fNPairCuts) + icut * (fNPairCuts * fNCommonTrackCuts) + iCommonCut * fNPairCuts + iPairCut].Data(), fValuesHadron); + for (isig = 0; isig < fRecMCSignals.size(); isig++) { + if (mcDecision & (static_cast(1) << isig)) { + fHistMan->FillHistClass(fHistNamesDileptonTrackMCmatched[(fNCuts * (fNCommonTrackCuts + 1) + fNCuts * fNPairCuts) + icut * (fNPairCuts * fNCommonTrackCuts) + iCommonCut * fNPairCuts + iPairCut][isig], fValuesHadron); + } + } // end loop (MC signals) + } + } // end loop (common track cuts) + } + } // end loop (pair cuts) + } + } // end loop (cuts) } // table to be written out for ML analysis BmesonsTable(fValuesHadron[VarManager::kPairMass], fValuesHadron[VarManager::kPairPt], fValuesHadron[VarManager::kVertexingLxy], fValuesHadron[VarManager::kVertexingLxyz], fValuesHadron[VarManager::kVertexingLz], fValuesHadron[VarManager::kVertexingTauxy], fValuesHadron[VarManager::kVertexingTauz], fValuesHadron[VarManager::kKFDCAxyzBetweenProngs], fValuesHadron[VarManager::kCosPointingAngle], fValuesHadron[VarManager::kVertexingChi2PCA], mcDecision); @@ -2165,6 +3259,7 @@ struct AnalysisDileptonTrack { Preslice trackAssocsPerCollision = aod::reducedtrack_association::reducedeventId; Preslice dielectronsPerCollision = aod::reducedpair::reducedeventId; + Preslice ditracksPerCollision = aod::reducedpair::reducedeventId; void processBarrelSkimmed(soa::Filtered const& events, soa::Filtered> const& assocs, @@ -2189,6 +3284,26 @@ struct AnalysisDileptonTrack { } } + void processDstarToD0Pi(soa::Filtered const& events, + soa::Filtered> const& assocs, + MyBarrelTracksWithCov const& tracks, soa::Filtered const& ditracks, + ReducedMCEvents const& mcEvents, ReducedMCTracks const& mcTracks) + { + // set up KF or DCAfitter + if (events.size() == 0) { + return; + } + if (fCurrentRun != events.begin().runNumber()) { // start: runNumber + initParamsFromCCDB(events.begin().timestamp()); + fCurrentRun = events.begin().runNumber(); + } // end: runNumber + for (auto& event : events) { + auto groupedBarrelAssocs = assocs.sliceBy(trackAssocsPerCollision, event.globalIndex()); + auto groupedDitracks = ditracks.sliceBy(ditracksPerCollision, event.globalIndex()); + runDileptonHadron(event, groupedBarrelAssocs, tracks, groupedDitracks, mcEvents, mcTracks); + } + } + Preslice muonAssocsPerCollision = aod::reducedtrack_association::reducedeventId; Preslice dimuonsPerCollision = aod::reducedpair::reducedeventId; @@ -2241,6 +3356,7 @@ struct AnalysisDileptonTrack { } PROCESS_SWITCH(AnalysisDileptonTrack, processBarrelSkimmed, "Run barrel dilepton-track pairing, using skimmed data", false); + PROCESS_SWITCH(AnalysisDileptonTrack, processDstarToD0Pi, "Run barrel pairing of D0 daughters with pion candidate, using skimmed data", false); PROCESS_SWITCH(AnalysisDileptonTrack, processMuonSkimmed, "Run muon dilepton-track pairing, using skimmed data", false); PROCESS_SWITCH(AnalysisDileptonTrack, processMCGen, "Loop over MC particle stack and fill generator level histograms", false); PROCESS_SWITCH(AnalysisDileptonTrack, processDummy, "Dummy function", false); @@ -2254,6 +3370,7 @@ WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) adaptAnalysisTask(cfgc), adaptAnalysisTask(cfgc), adaptAnalysisTask(cfgc), + adaptAnalysisTask(cfgc), adaptAnalysisTask(cfgc)}; } @@ -2308,6 +3425,10 @@ void DefineHistograms(HistogramManager* histMan, TString histClasses, const char dqhistograms::DefineHistograms(histMan, objArray->At(iclass)->GetName(), "pair", histName); } + if (classStr.Contains("Triplets")) { + dqhistograms::DefineHistograms(histMan, objArray->At(iclass)->GetName(), "pair", histName); + } + if (classStr.Contains("MCTruthGenPair")) { dqhistograms::DefineHistograms(histMan, objArray->At(iclass)->GetName(), "mctruth_pair"); } diff --git a/PWGDQ/Tasks/tableReader_withAssoc.cxx b/PWGDQ/Tasks/tableReader_withAssoc.cxx index 739a9ff1957..f3cb779db40 100644 --- a/PWGDQ/Tasks/tableReader_withAssoc.cxx +++ b/PWGDQ/Tasks/tableReader_withAssoc.cxx @@ -2397,7 +2397,7 @@ struct AnalysisAsymmetricPairing { std::set> globIdxTriplets; // Based on triplet type, make suitable combinations of the partitions - if (tripletType == VarManager::kTripleCandidateToPKPi || tripletType == VarManager::kTripleCandidateToKKPi) { + if (tripletType == VarManager::kTripleCandidateToPKPi) { for (auto& [a1, a2, a3] : combinations(soa::CombinationsFullIndexPolicy(groupedLegAAssocs, groupedLegBAssocs, groupedLegCAssocs))) { readTriplet(a1, a2, a3, tracks, event, tripletType, histNames); } @@ -2517,13 +2517,6 @@ struct AnalysisAsymmetricPairing { runThreeProng(events, trackAssocsPerCollision, barrelAssocs, barrelTracks, VarManager::kTripleCandidateToKPiPi); } - void processKaonKaonPionSkimmed(MyEventsVtxCovZdcSelected const& events, - soa::Join const& barrelAssocs, - MyBarrelTracksWithCovWithAmbiguities const& barrelTracks) - { - runThreeProng(events, trackAssocsPerCollision, barrelAssocs, barrelTracks, VarManager::kTripleCandidateToKKPi); - } - void processProtonKaonPionSkimmed(MyEventsVtxCovZdcSelected const& events, soa::Join const& barrelAssocs, MyBarrelTracksWithCovWithAmbiguities const& barrelTracks) @@ -2538,7 +2531,6 @@ struct AnalysisAsymmetricPairing { PROCESS_SWITCH(AnalysisAsymmetricPairing, processKaonPionSkimmed, "Run kaon pion pairing, with skimmed tracks", false); PROCESS_SWITCH(AnalysisAsymmetricPairing, processKaonPionPionSkimmed, "Run kaon pion pion triplets, with skimmed tracks", false); - PROCESS_SWITCH(AnalysisAsymmetricPairing, processKaonKaonPionSkimmed, "Run kaon kaon pion triplets, with skimmed tracks", false); PROCESS_SWITCH(AnalysisAsymmetricPairing, processProtonKaonPionSkimmed, "Run proton kaon pion triplets, with skimmed tracks", false); PROCESS_SWITCH(AnalysisAsymmetricPairing, processDummy, "Dummy function, enabled only if none of the others are enabled", true); }; From 839b2a922fcbcc956a088954076ef2c230b0aa56 Mon Sep 17 00:00:00 2001 From: sigurd Date: Tue, 10 Dec 2024 11:19:33 +0100 Subject: [PATCH 2/7] Remove unnecessary code for 2- and 3-prong signals --- PWGDQ/Tasks/dqEfficiency_withAssoc.cxx | 51 -------------------------- 1 file changed, 51 deletions(-) diff --git a/PWGDQ/Tasks/dqEfficiency_withAssoc.cxx b/PWGDQ/Tasks/dqEfficiency_withAssoc.cxx index 6391e4d20fd..4ca955d1d28 100644 --- a/PWGDQ/Tasks/dqEfficiency_withAssoc.cxx +++ b/PWGDQ/Tasks/dqEfficiency_withAssoc.cxx @@ -1945,8 +1945,6 @@ struct AnalysisAsymmetricPairing { int fNLegCuts; int fNPairCuts = 0; int fNCommonTrackCuts; - bool fHasTwoProngGenMCsignals = false; - bool fHasThreeProngGenMCsignals = false; Preslice> trackAssocsPerCollision = aod::reducedtrack_association::reducedeventId; @@ -2250,14 +2248,6 @@ struct AnalysisAsymmetricPairing { if (sig->GetNProngs() == 1) { // NOTE: 1-prong signals required fGenMCSignals.push_back(*sig); histNames += Form("MCTruthGen_%s;", sig->GetName()); // TODO: Add these names to a std::vector to avoid using Form in the process function - } else if (sig->GetNProngs() == 2) { // NOTE: 2-prong signals required - fGenMCSignals.push_back(*sig); - histNames += Form("MCTruthGenPair_%s;", sig->GetName()); - fHasTwoProngGenMCsignals = true; - } else if (sig->GetNProngs() == 3) { // NOTE: 3-prong signals required - fGenMCSignals.push_back(*sig); - histNames += Form("MCTruthGenTriplet_%s;", sig->GetName()); - fHasThreeProngGenMCsignals = true; } } } @@ -2776,47 +2766,6 @@ struct AnalysisAsymmetricPairing { } } } - - if (fHasTwoProngGenMCsignals) { - for (auto& event : mcEvents) { - auto groupedMCTracks = mcTracks.sliceBy(perReducedMcEvent, event.globalIndex()); - groupedMCTracks.bindInternalIndicesTo(&mcTracks); - for (auto& [t1, t2] : combinations(groupedMCTracks, groupedMCTracks)) { - auto t1_raw = groupedMCTracks.rawIteratorAt(t1.globalIndex()); - auto t2_raw = groupedMCTracks.rawIteratorAt(t2.globalIndex()); - for (auto& sig : fGenMCSignals) { - if (sig.GetNProngs() != 2) { // NOTE: 2-prong signals required here - continue; - } - if (sig.CheckSignal(true, t1_raw, t2_raw)) { - VarManager::FillPairMC(t1, t2, VarManager::fgValues, pairType); - fHistMan->FillHistClass(Form("MCTruthGenPair_%s", sig.GetName()), VarManager::fgValues); - } - } // end loop over MC signals - } // end loop over pairs - } // end loop over events - } - - if (fHasThreeProngGenMCsignals) { - for (auto& event : mcEvents) { - auto groupedMCTracks = mcTracks.sliceBy(perReducedMcEvent, event.globalIndex()); - groupedMCTracks.bindInternalIndicesTo(&mcTracks); - for (auto& [t1, t2, t3] : combinations(groupedMCTracks, groupedMCTracks, groupedMCTracks)) { - auto t1_raw = groupedMCTracks.rawIteratorAt(t1.globalIndex()); - auto t2_raw = groupedMCTracks.rawIteratorAt(t2.globalIndex()); - auto t3_raw = groupedMCTracks.rawIteratorAt(t3.globalIndex()); - for (auto& sig : fGenMCSignals) { - if (sig.GetNProngs() != 3) { // NOTE: 3-prong signals required here - continue; - } - if (sig.CheckSignal(true, t1_raw, t2_raw, t3_raw)) { - VarManager::FillTripleMC(t1, t2, t3, VarManager::fgValues, pairType); - fHistMan->FillHistClass(Form("MCTruthGenTriplet_%s", sig.GetName()), VarManager::fgValues); - } - } // end loop over MC signals - } // end loop over pairs - } // end loop over events - } } // end runMCGen void processKaonPionSkimmed(MyEventsVtxCovSelected const& events, From f059d990661839235121a54a9636f16158a04ccd Mon Sep 17 00:00:00 2001 From: ALICE Builder Date: Tue, 10 Dec 2024 11:23:05 +0100 Subject: [PATCH 3/7] Please consider the following formatting changes (#9) --- PWGDQ/Tasks/dqEfficiency_withAssoc.cxx | 19 ++++++++++--------- 1 file changed, 10 insertions(+), 9 deletions(-) diff --git a/PWGDQ/Tasks/dqEfficiency_withAssoc.cxx b/PWGDQ/Tasks/dqEfficiency_withAssoc.cxx index 4ca955d1d28..f1e08248226 100644 --- a/PWGDQ/Tasks/dqEfficiency_withAssoc.cxx +++ b/PWGDQ/Tasks/dqEfficiency_withAssoc.cxx @@ -2105,7 +2105,7 @@ struct AnalysisAsymmetricPairing { fTrackHistNames[(fNLegCuts * (fNCommonTrackCuts + 1) + fNLegCuts * fNPairCuts) + icut * (fNPairCuts * fNCommonTrackCuts + 1) + iCommonCut * (1 + fNPairCuts) + iPairCut] = names; } // end loop (common cuts) } // end loop (pair cuts) - } // end if (pair cuts) + } // end if (pair cuts) // TODO: assign hist directories for the MC matched triplets for each (leg cut combo,MCsignal) combination if (!sigNamesStr.IsNull()) { @@ -2191,7 +2191,7 @@ struct AnalysisAsymmetricPairing { fTrackHistNames[(fNLegCuts * (fNCommonTrackCuts + 1) + fNLegCuts * fNPairCuts) + icut * (fNPairCuts * fNCommonTrackCuts) + iCommonCut * (1 + fNPairCuts) + iPairCut] = names; } // end loop (common cuts) } // end loop (pair cuts) - } // end if (pair cuts) + } // end if (pair cuts) // assign hist directories for the MC matched triplets for each (leg cut combo,MCsignal) combination if (!sigNamesStr.IsNull()) { @@ -2303,7 +2303,8 @@ struct AnalysisAsymmetricPairing { } else { VarManager::SetupThreeProngDCAFitter(magField, true, 200.0f, 4.0f, 1.0e-3f, 0.9f, fConfigUseAbsDCA.value); } - } else { if (fConfigUseKFVertexing.value) { + } else { + if (fConfigUseKFVertexing.value) { VarManager::SetupTwoProngKFParticle(magField); } else { VarManager::SetupTwoProngDCAFitter(magField, true, 200.0f, 4.0f, 1.0e-3f, 0.9f, fConfigUseAbsDCA.value); // TODO: get these parameters from Configurables @@ -2557,8 +2558,8 @@ struct AnalysisAsymmetricPairing { if constexpr (trackHasCov && TTwoProngFitter) { ditrackExtraList(t1.globalIndex(), t2.globalIndex(), VarManager::fgValues[VarManager::kVertexingTauzProjected], VarManager::fgValues[VarManager::kVertexingLzProjected], VarManager::fgValues[VarManager::kVertexingLxyProjected]); } - } // end inner assoc loop (leg A) - } // end event loop + } // end inner assoc loop (leg A) + } // end event loop } // Template function to run same event triplets (e.g. D+->K-pi+pi+) @@ -2694,7 +2695,7 @@ struct AnalysisAsymmetricPairing { if (threeTrackFilter & (static_cast(1) << icut)) { fHistMan->FillHistClass(histNames[icut][0].Data(), VarManager::fgValues); // TODO: loop over MC signals - for (unsigned int isig = 0; isig < fRecMCSignals.size(); isig++) { // loop over MC signals + for (unsigned int isig = 0; isig < fRecMCSignals.size(); isig++) { // loop over MC signals int offset = fNLegCuts * isig * (1 + fNCommonTrackCuts + fNPairCuts + fNCommonTrackCuts * fNPairCuts); if (mcDecision & (static_cast(1) << isig)) { fHistMan->FillHistClass(histNamesMC[offset + icut][0].Data(), VarManager::fgValues); // matched signal @@ -2706,7 +2707,7 @@ struct AnalysisAsymmetricPairing { for (int iCommonCut = 0; iCommonCut < fNCommonTrackCuts; iCommonCut++) { if (threeTrackCommonFilter & fCommonTrackCutFilterMasks[iCommonCut]) { fHistMan->FillHistClass(histNames[fNLegCuts + icut * fNCommonTrackCuts + iCommonCut][0].Data(), VarManager::fgValues); - for (unsigned int isig = 0; isig < fRecMCSignals.size(); isig++) { // loop over MC signals + for (unsigned int isig = 0; isig < fRecMCSignals.size(); isig++) { // loop over MC signals int offset = fNLegCuts * isig * (1 + fNCommonTrackCuts + fNPairCuts + fNCommonTrackCuts * fNPairCuts); if (mcDecision & (static_cast(1) << isig)) { fHistMan->FillHistClass(histNamesMC[offset + fNLegCuts + icut * fNCommonTrackCuts + iCommonCut][0].Data(), VarManager::fgValues); // matched signal @@ -2721,7 +2722,7 @@ struct AnalysisAsymmetricPairing { } // Histograms with pair cuts fHistMan->FillHistClass(histNames[fNLegCuts * (fNCommonTrackCuts + 1) + icut * fNPairCuts + iPairCut][0].Data(), VarManager::fgValues); - for (unsigned int isig = 0; isig < fRecMCSignals.size(); isig++) { // loop over MC signals + for (unsigned int isig = 0; isig < fRecMCSignals.size(); isig++) { // loop over MC signals int offset = fNLegCuts * isig * (1 + fNCommonTrackCuts + fNPairCuts + fNCommonTrackCuts * fNPairCuts); if (mcDecision & (static_cast(1) << isig)) { fHistMan->FillHistClass(histNamesMC[offset + fNLegCuts * (fNCommonTrackCuts + 1) + icut * fNPairCuts + iPairCut][0].Data(), VarManager::fgValues); // matched signal @@ -2731,7 +2732,7 @@ struct AnalysisAsymmetricPairing { for (int iCommonCut = 0; iCommonCut < fNCommonTrackCuts; ++iCommonCut) { if (threeTrackCommonFilter & fCommonTrackCutFilterMasks[iCommonCut]) { fHistMan->FillHistClass(histNames[(fNLegCuts * (fNCommonTrackCuts + 1) + fNLegCuts * fNPairCuts) + icut * (fNPairCuts * fNCommonTrackCuts + 1) + iCommonCut * (1 + fNPairCuts) + iPairCut][0].Data(), VarManager::fgValues); - for (unsigned int isig = 0; isig < fRecMCSignals.size(); isig++) { // loop over MC signals + for (unsigned int isig = 0; isig < fRecMCSignals.size(); isig++) { // loop over MC signals int offset = fNLegCuts * isig * (1 + fNCommonTrackCuts + fNPairCuts + fNCommonTrackCuts * fNPairCuts); if (mcDecision & (static_cast(1) << isig)) { fHistMan->FillHistClass(histNamesMC[offset + (fNLegCuts * (fNCommonTrackCuts + 1) + fNLegCuts * fNPairCuts) + icut * (fNPairCuts * fNCommonTrackCuts) + iCommonCut * (1 + fNPairCuts) + iPairCut][0].Data(), VarManager::fgValues); // matched signal From 2ddde3e8cb9917bcbeefbe361fdc22c7bcefe804 Mon Sep 17 00:00:00 2001 From: sigurd Date: Tue, 10 Dec 2024 11:37:45 +0100 Subject: [PATCH 4/7] MegaLinter: Remove unused variables --- PWGDQ/Tasks/dqEfficiency_withAssoc.cxx | 3 --- PWGDQ/Tasks/tableReader_withAssoc.cxx | 3 --- 2 files changed, 6 deletions(-) diff --git a/PWGDQ/Tasks/dqEfficiency_withAssoc.cxx b/PWGDQ/Tasks/dqEfficiency_withAssoc.cxx index f1e08248226..c8bc424a2fc 100644 --- a/PWGDQ/Tasks/dqEfficiency_withAssoc.cxx +++ b/PWGDQ/Tasks/dqEfficiency_withAssoc.cxx @@ -2366,8 +2366,6 @@ struct AnalysisAsymmetricPairing { continue; } - // TODO: Think about double counting - std::set> globIdxPairs; for (auto& [a1, a2] : combinations(soa::CombinationsFullIndexPolicy(groupedLegAAssocs, groupedLegBAssocs))) { uint32_t twoTrackFilter = 0; @@ -2597,7 +2595,6 @@ struct AnalysisAsymmetricPairing { continue; } - std::set> globIdxTriplets; // Based on triplet type, make suitable combinations of the partitions if (tripletType == VarManager::kTripleCandidateToPKPi) { for (auto& [a1, a2, a3] : combinations(soa::CombinationsFullIndexPolicy(groupedLegAAssocs, groupedLegBAssocs, groupedLegCAssocs))) { diff --git a/PWGDQ/Tasks/tableReader_withAssoc.cxx b/PWGDQ/Tasks/tableReader_withAssoc.cxx index f3cb779db40..7d4c90fcefc 100644 --- a/PWGDQ/Tasks/tableReader_withAssoc.cxx +++ b/PWGDQ/Tasks/tableReader_withAssoc.cxx @@ -2232,8 +2232,6 @@ struct AnalysisAsymmetricPairing { continue; } - // TODO: Think about double counting - std::set> globIdxPairs; for (auto& [a1, a2] : combinations(soa::CombinationsFullIndexPolicy(groupedLegAAssocs, groupedLegBAssocs))) { uint32_t twoTrackFilter = static_cast(0); @@ -2395,7 +2393,6 @@ struct AnalysisAsymmetricPairing { continue; } - std::set> globIdxTriplets; // Based on triplet type, make suitable combinations of the partitions if (tripletType == VarManager::kTripleCandidateToPKPi) { for (auto& [a1, a2, a3] : combinations(soa::CombinationsFullIndexPolicy(groupedLegAAssocs, groupedLegBAssocs, groupedLegCAssocs))) { From b48c7cb78270cbf7abebb6aec6d6994ed2c4155a Mon Sep 17 00:00:00 2001 From: sigurd Date: Tue, 10 Dec 2024 12:04:11 +0100 Subject: [PATCH 5/7] Remove MC signal used for debugging --- PWGDQ/Core/MCSignalLibrary.cxx | 6 ------ 1 file changed, 6 deletions(-) diff --git a/PWGDQ/Core/MCSignalLibrary.cxx b/PWGDQ/Core/MCSignalLibrary.cxx index ca390993dd9..acbbbf418d4 100644 --- a/PWGDQ/Core/MCSignalLibrary.cxx +++ b/PWGDQ/Core/MCSignalLibrary.cxx @@ -1269,12 +1269,6 @@ MCSignal* o2::aod::dqmcsignals::GetMCSignal(const char* name) signal = new MCSignal(name, "Kaon pion pion triplet from D*", {prongKaon, prongPionSecondary, prongPion}, {2, 2, 1}); return signal; } - if (!nameStr.compare("PiPiPiFromD0FromDstar")) { - MCProng prongPionSecondary(3, {211, Pdg::kD0, Pdg::kDStar}, {true, true, true}, {false, false, false}, {0, 0, 0}, {0, 0, 0}, {false, false, false}); - MCProng prongPion(2, {211, Pdg::kDStar}, {true, true}, {false, false}, {0, 0}, {0, 0}, {false, false}); - signal = new MCSignal(name, "Kaon pion pion triplet from D*", {prongPionSecondary, prongPionSecondary, prongPion}, {2, 2, 1}); - return signal; - } if (!nameStr.compare("KFromDplus")) { MCProng prong(2, {321, Pdg::kDPlus}, {true, true}, {false, false}, {0, 0}, {0, 0}, {false, false}, false, {502}, {true}); prong.SetSourceBit(0, MCProng::kPhysicalPrimary); From c71c909fde971e663c13505001b09033566c0ebe Mon Sep 17 00:00:00 2001 From: sigurd Date: Wed, 11 Dec 2024 15:45:30 +0100 Subject: [PATCH 6/7] Add charge-specific D* signals and LambdaC signal --- PWGDQ/Core/MCSignalLibrary.cxx | 29 +++++++++++++++++++++++++++++ 1 file changed, 29 insertions(+) diff --git a/PWGDQ/Core/MCSignalLibrary.cxx b/PWGDQ/Core/MCSignalLibrary.cxx index acbbbf418d4..6cef5603b9a 100644 --- a/PWGDQ/Core/MCSignalLibrary.cxx +++ b/PWGDQ/Core/MCSignalLibrary.cxx @@ -1236,6 +1236,16 @@ MCSignal* o2::aod::dqmcsignals::GetMCSignal(const char* name) signal = new MCSignal(name, "D*", {prong}, {-1}); return signal; } + if (!nameStr.compare("DstarPlus")) { + MCProng prong(1, {Pdg::kDStar}, {false}, {false}, {0}, {0}, {false}); + signal = new MCSignal(name, "D*+", {prong}, {-1}); + return signal; + } + if (!nameStr.compare("DstarMinus")) { + MCProng prong(1, {-Pdg::kDStar}, {false}, {false}, {0}, {0}, {false}); + signal = new MCSignal(name, "D*-", {prong}, {-1}); + return signal; + } if (!nameStr.compare("pionFromDstar")) { MCProng prong(2, {211, Pdg::kDStar}, {true, true}, {false, false}, {0, 0}, {0, 0}, {false, false}); signal = new MCSignal(name, "Pions from D* decays", {prong}, {1}); @@ -1269,12 +1279,31 @@ MCSignal* o2::aod::dqmcsignals::GetMCSignal(const char* name) signal = new MCSignal(name, "Kaon pion pion triplet from D*", {prongKaon, prongPionSecondary, prongPion}, {2, 2, 1}); return signal; } + if (!nameStr.compare("KPiPiFromD0FromDstarPlus")) { + MCProng prongKaon(3, {-321, Pdg::kD0, Pdg::kDStar}, {false, false, false}, {false, false, false}, {0, 0, 0}, {0, 0, 0}, {false, false, false}); + MCProng prongPionSecondary(3, {211, Pdg::kD0, Pdg::kDStar}, {false, false, false}, {false, false, false}, {0, 0, 0}, {0, 0, 0}, {false, false, false}); + MCProng prongPion(2, {211, Pdg::kDStar}, {false, false}, {false, false}, {0, 0}, {0, 0}, {false, false}); + signal = new MCSignal(name, "Kaon pion pion triplet from D*+", {prongKaon, prongPionSecondary, prongPion}, {2, 2, 1}); + return signal; + } + if (!nameStr.compare("KPiPiFromD0FromDstarMinus")) { + MCProng prongKaon(3, {321, Pdg::kD0, Pdg::kDStar}, {false, false, false}, {false, false, false}, {0, 0, 0}, {0, 0, 0}, {false, false, false}); + MCProng prongPionSecondary(3, {-211, Pdg::kD0, Pdg::kDStar}, {false, false, false}, {false, false, false}, {0, 0, 0}, {0, 0, 0}, {false, false, false}); + MCProng prongPion(2, {-211, Pdg::kDStar}, {false, false}, {false, false}, {0, 0}, {0, 0}, {false, false}); + signal = new MCSignal(name, "Kaon pion pion triplet from D*-", {prongKaon, prongPionSecondary, prongPion}, {2, 2, 1}); + return signal; + } if (!nameStr.compare("KFromDplus")) { MCProng prong(2, {321, Pdg::kDPlus}, {true, true}, {false, false}, {0, 0}, {0, 0}, {false, false}, false, {502}, {true}); prong.SetSourceBit(0, MCProng::kPhysicalPrimary); signal = new MCSignal(name, "Kaons from D+/- decays", {prong}, {-1}); return signal; } + if (!nameStr.compare("LambdaC")) { + MCProng prong(1, {Pdg::kLambdaCPlus}, {true}, {false}, {0}, {0}, {false}); + signal = new MCSignal(name, "Lambda_c", {prong}, {-1}); + return signal; + } //-------------------------------------------------------------------------------- From 54dee398ef3fc46b35a1632864fbfceffb6ec492 Mon Sep 17 00:00:00 2001 From: sigurd Date: Fri, 13 Dec 2024 11:27:07 +0100 Subject: [PATCH 7/7] Remove unused parameters in function call --- PWGDQ/Tasks/dqEfficiency_withAssoc.cxx | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/PWGDQ/Tasks/dqEfficiency_withAssoc.cxx b/PWGDQ/Tasks/dqEfficiency_withAssoc.cxx index c8bc424a2fc..7cdb426d55e 100644 --- a/PWGDQ/Tasks/dqEfficiency_withAssoc.cxx +++ b/PWGDQ/Tasks/dqEfficiency_withAssoc.cxx @@ -2744,7 +2744,7 @@ struct AnalysisAsymmetricPairing { PresliceUnsorted perReducedMcEvent = aod::reducedtrackMC::reducedMCeventId; - void runMCGen(ReducedMCEvents const& mcEvents, ReducedMCTracks const& mcTracks, VarManager::PairCandidateType pairType) + void runMCGen(ReducedMCTracks const& mcTracks) { // loop over mc stack and fill histograms for pure MC truth signals // group all the MC tracks which belong to the MC event corresponding to the current reconstructed event @@ -2773,7 +2773,7 @@ struct AnalysisAsymmetricPairing { { runAsymmetricPairing(events, trackAssocsPerCollision, barrelAssocs, barrelTracks, mcEvents, mcTracks); if (fConfigRunMCGenPair) - runMCGen(mcEvents, mcTracks, VarManager::kDecayToKPi); + runMCGen(mcTracks); } void processKaonPionPionSkimmed(MyEventsVtxCovSelected const& events, @@ -2783,7 +2783,7 @@ struct AnalysisAsymmetricPairing { { runThreeProng(events, trackAssocsPerCollision, barrelAssocs, barrelTracks, mcEvents, mcTracks, VarManager::kTripleCandidateToKPiPi); if (fConfigRunMCGenPair) - runMCGen(mcEvents, mcTracks, VarManager::kTripleCandidateToKPiPi); + runMCGen(mcTracks); } void processDummy(MyEvents&)