diff --git a/PWGCF/Femto/Core/cascadeBuilder.h b/PWGCF/Femto/Core/cascadeBuilder.h index 2bdfffd3a0b..dd022f91501 100644 --- a/PWGCF/Femto/Core/cascadeBuilder.h +++ b/PWGCF/Femto/Core/cascadeBuilder.h @@ -90,17 +90,17 @@ struct ConfOmegaBits : o2::framework::ConfigurableGroup { #undef CASCADE_DEFAULT_BITS -#define CASCADE_DEFAULT_SELECTION(defaultMassMin, defaultMassMax, defaultPdgCode) \ - o2::framework::Configurable pdgCode{"pdgCode", defaultPdgCode, "Track PDG code"}; \ - o2::framework::Configurable sign{"sign", 1, "Sign of the charge of the Cascade "}; \ - o2::framework::Configurable ptMin{"ptMin", 0.f, "Minimum pT"}; \ - o2::framework::Configurable ptMax{"ptMax", 999.f, "Maximum pT"}; \ - o2::framework::Configurable etaMin{"etaMin", -10.f, "Minimum eta"}; \ - o2::framework::Configurable etaMax{"etaMax", 10.f, "Maximum eta"}; \ - o2::framework::Configurable phiMin{"phiMin", 0.f, "Minimum eta"}; \ - o2::framework::Configurable phiMax{"phiMax", 1.f * o2::constants::math::TwoPI, "Maximum phi"}; \ - o2::framework::Configurable massMin{"massMin", defaultMassMin, "Minimum invariant mass for Cascade"}; \ - o2::framework::Configurable massMax{"massMax", defaultMassMax, "Maximum invariant mass for Cascade"}; \ +#define CASCADE_DEFAULT_SELECTION(defaultMassMin, defaultMassMax, defaultPdgCode) \ + o2::framework::Configurable pdgCodeAbs{"pdgCodeAbs", defaultPdgCode, "Cascade PDG code. Set sign to +1 to select antiparticle"}; \ + o2::framework::Configurable sign{"sign", -1, "Sign of the charge of the Cascade"}; \ + o2::framework::Configurable ptMin{"ptMin", 0.f, "Minimum pT"}; \ + o2::framework::Configurable ptMax{"ptMax", 999.f, "Maximum pT"}; \ + o2::framework::Configurable etaMin{"etaMin", -10.f, "Minimum eta"}; \ + o2::framework::Configurable etaMax{"etaMax", 10.f, "Maximum eta"}; \ + o2::framework::Configurable phiMin{"phiMin", 0.f, "Minimum eta"}; \ + o2::framework::Configurable phiMax{"phiMax", 1.f * o2::constants::math::TwoPI, "Maximum phi"}; \ + o2::framework::Configurable massMin{"massMin", defaultMassMin, "Minimum invariant mass for Cascade"}; \ + o2::framework::Configurable massMax{"massMax", defaultMassMax, "Maximum invariant mass for Cascade"}; \ o2::framework::Configurable mask{"mask", 0x0, "Bitmask for cascade selection"}; struct ConfXiSelection : o2::framework::ConfigurableGroup { diff --git a/PWGCF/Femto/Core/cascadeHistManager.h b/PWGCF/Femto/Core/cascadeHistManager.h index 05de039359f..e4bd045f7d5 100644 --- a/PWGCF/Femto/Core/cascadeHistManager.h +++ b/PWGCF/Femto/Core/cascadeHistManager.h @@ -173,41 +173,122 @@ class CascadeHistManager CascadeHistManager() = default; ~CascadeHistManager() = default; - template + template void init(o2::framework::HistogramRegistry* registry, std::map> const& cascadeSpecs, + T const& ConfCascadeSelection, std::map> const& BachelorSpecs, std::map> const& PosDauSpecs, std::map> const& NegDauSpecs) { mHistogramRegistry = registry; - mBachelorManager.template init(registry, BachelorSpecs, AbsChargeDaughters); - mPosDauManager.template init(registry, PosDauSpecs, AbsChargeDaughters); - mNegDauManager.template init(registry, NegDauSpecs, AbsChargeDaughters); + mPdgCode = std::abs(ConfCascadeSelection.pdgCodeAbs.value); + + int bachelorPdgCodeAbs = 0; + int posDauPdgCodeAbs = 0; + int negDauPdgCodeAbs = 0; + const int absCharge = 1; + int signBachelor = 0; + const int signPlus = 1; + const int signMinus = -1; + + if (mPdgCode == PDG_t::kXiMinus) { + if (ConfCascadeSelection.sign.value < 0) { + bachelorPdgCodeAbs = std::abs(PDG_t::kPiMinus); + signBachelor = -1; + posDauPdgCodeAbs = std::abs(PDG_t::kProton); + negDauPdgCodeAbs = std::abs(PDG_t::kPiMinus); + } else { + mPdgCode = -1 * mPdgCode; // Xi+ has negative pdg code + bachelorPdgCodeAbs = std::abs(PDG_t::kPiPlus); + signBachelor = 1; + posDauPdgCodeAbs = std::abs(PDG_t::kPiPlus); + negDauPdgCodeAbs = std::abs(PDG_t::kProtonBar); + } + } else if (mPdgCode == PDG_t::kOmegaMinus) { + if (ConfCascadeSelection.sign.value < 0) { + bachelorPdgCodeAbs = std::abs(PDG_t::kKMinus); + signBachelor = -1; + posDauPdgCodeAbs = std::abs(PDG_t::kProton); + negDauPdgCodeAbs = std::abs(PDG_t::kPiMinus); + } else { + mPdgCode = -1 * mPdgCode; // Omega+ has negative pdg code + bachelorPdgCodeAbs = std::abs(PDG_t::kKPlus); + signBachelor = 1; + posDauPdgCodeAbs = std::abs(PDG_t::kPiPlus); + negDauPdgCodeAbs = std::abs(PDG_t::kProtonBar); + } + } else { + LOG(fatal) << "PDG code for Cascade has to be either Xi or Omega"; + } + + mBachelorManager.template init(registry, BachelorSpecs, absCharge, signBachelor, bachelorPdgCodeAbs); + mPosDauManager.template init(registry, PosDauSpecs, absCharge, signPlus, posDauPdgCodeAbs); + mNegDauManager.template init(registry, NegDauSpecs, absCharge, signMinus, negDauPdgCodeAbs); + if constexpr (modes::isFlagSet(mode, modes::Mode::kAnalysis)) { initAnalysis(cascadeSpecs); } - if constexpr (modes::isFlagSet(mode, modes::Mode::kQa)) { - initQa(cascadeSpecs); - } } - template + template void init(o2::framework::HistogramRegistry* registry, std::map> const& cascadeSpecs, - T1 const& CascadeConfBinningQa, + T1 const& ConfCascadeSelection, + T2 const& ConfCascadeBinningQa, std::map> const& BachelorSpecs, - T2 const& BachelorConfBinningQa, + T3 const& ConfBachelorQaBinning, std::map> const& PosDauSpecs, - T3 const& PosDauConfBinningQa, + T4& ConfPosDauQaBinning, std::map> const& NegDauSpecs, - T4 const& NegDauConfBinningQa) + T5& ConfNegDauQaBinning) { mHistogramRegistry = registry; - mBachelorManager.template init(registry, BachelorSpecs, BachelorConfBinningQa, AbsChargeDaughters); - mPosDauManager.template init(registry, PosDauSpecs, PosDauConfBinningQa, AbsChargeDaughters); - mNegDauManager.template init(registry, NegDauSpecs, NegDauConfBinningQa, AbsChargeDaughters); - this->enableOptionalHistograms(CascadeConfBinningQa); + mPdgCode = std::abs(ConfCascadeSelection.pdgCodeAbs.value); + this->enableOptionalHistograms(ConfCascadeBinningQa); + + int bachelorPdgCodeAbs = 0; + int posDauPdgCodeAbs = 0; + int negDauPdgCodeAbs = 0; + const int absCharge = 1; + int signBachelor = 0; + const int signPlus = 1; + const int signMinus = -1; + + if (mPdgCode == PDG_t::kXiMinus) { + if (ConfCascadeSelection.sign.value < 0) { + bachelorPdgCodeAbs = std::abs(PDG_t::kPiMinus); + signBachelor = -1; + posDauPdgCodeAbs = std::abs(PDG_t::kProton); + negDauPdgCodeAbs = std::abs(PDG_t::kPiMinus); + } else { + mPdgCode = -1 * mPdgCode; // Xi+ has negative pdg code + bachelorPdgCodeAbs = std::abs(PDG_t::kPiPlus); + signBachelor = 1; + posDauPdgCodeAbs = std::abs(PDG_t::kPiPlus); + negDauPdgCodeAbs = std::abs(PDG_t::kProtonBar); + } + } else if (mPdgCode == PDG_t::kOmegaMinus) { + if (ConfCascadeSelection.sign.value < 0) { + bachelorPdgCodeAbs = std::abs(PDG_t::kKMinus); + signBachelor = -1; + posDauPdgCodeAbs = std::abs(PDG_t::kProton); + negDauPdgCodeAbs = std::abs(PDG_t::kPiMinus); + } else { + mPdgCode = -1 * mPdgCode; // Omega+ has negative pdg code + bachelorPdgCodeAbs = std::abs(PDG_t::kKPlus); + signBachelor = 1; + posDauPdgCodeAbs = std::abs(PDG_t::kPiPlus); + negDauPdgCodeAbs = std::abs(PDG_t::kProtonBar); + } + } else { + LOG(fatal) << "PDG code for Cascade has to be either Xi or Omega"; + } + + mBachelorManager.template init(registry, BachelorSpecs, absCharge, signBachelor, bachelorPdgCodeAbs, ConfBachelorQaBinning); + mPosDauManager.template init(registry, PosDauSpecs, absCharge, signPlus, posDauPdgCodeAbs, ConfPosDauQaBinning); + mNegDauManager.template init(registry, NegDauSpecs, absCharge, signMinus, negDauPdgCodeAbs, ConfNegDauQaBinning); + if constexpr (modes::isFlagSet(mode, modes::Mode::kAnalysis)) { initAnalysis(cascadeSpecs); } @@ -321,6 +402,7 @@ class CascadeHistManager trackhistmanager::TrackHistManager mBachelorManager; trackhistmanager::TrackHistManager mPosDauManager; trackhistmanager::TrackHistManager mNegDauManager; + int mPdgCode = 0; }; }; // namespace cascadehistmanager }; // namespace o2::analysis::femto diff --git a/PWGCF/Femto/Core/collisionBuilder.h b/PWGCF/Femto/Core/collisionBuilder.h index aef859c5163..0d094bc5abf 100644 --- a/PWGCF/Femto/Core/collisionBuilder.h +++ b/PWGCF/Femto/Core/collisionBuilder.h @@ -89,6 +89,7 @@ struct ConfCcdb : o2::framework::ConfigurableGroup { o2::framework::Configurable ccdbUrl{"ccdbUrl", "http://alice-ccdb.cern.ch", "URL to ccdb"}; o2::framework::Configurable grpPath{"grpPath", "GLO/Config/GRPMagField", "Path to GRP object (Run3 -> GLO/Config/GRPMagField, Run2 -> GLO/GRP/GRP"}; o2::framework::Configurable triggerPath{"triggerPath", "EventFiltering/Zorro/", "CCDB path for trigger information"}; + o2::framework::Configurable magFieldForced{"magFieldForced", 0, "Force value for magnetic field (kG). This will skip calls to the ccdb. Deactivate by setting value to 0"}; }; struct ConfCollisionRctFlags : o2::framework::ConfigurableGroup { @@ -390,6 +391,7 @@ class CollisionBuilder mUseRctFlags = true; mRctFlagsChecker.init(confRct.label.value, confRct.useZdc.value, confRct.treatLimitedAcceptanceAsBad.value); } + mMagFieldForced = confCcdb.magFieldForced.value; mGrpPath = confCcdb.grpPath.value; mSubGeneratorId = confFilter.subGeneratorId.value; @@ -403,12 +405,16 @@ class CollisionBuilder { if (mRunNumber != bc.runNumber()) { mRunNumber = bc.runNumber(); - static o2::parameters::GRPMagField* grpo = nullptr; - grpo = ccdb->template getForRun(mGrpPath, mRunNumber); - if (grpo == nullptr) { - LOG(fatal) << "GRP object not found for Run " << mRunNumber; + if (mMagFieldForced == 0) { + static o2::parameters::GRPMagField* grpo = nullptr; + grpo = ccdb->template getForRun(mGrpPath, mRunNumber); + if (grpo == nullptr) { + LOG(fatal) << "GRP object not found for Run " << mRunNumber; + } + mMagField = static_cast(grpo->getNominalL3Field()); // get magnetic field in kG + } else { + mMagField = mMagFieldForced; } - mMagField = static_cast(grpo->getNominalL3Field()); // get magnetic field in kG if (mUseTrigger) { mZorro.initCCDB(ccdb.service, mRunNumber, bc.timestamp(), mTriggerNames); @@ -534,6 +540,7 @@ class CollisionBuilder bool mUseTrigger = false; int mRunNumber = -1; std::string mGrpPath = std::string(""); + int mMagFieldForced = 0; int mMagField = 0; aod::rctsel::RCTFlagsChecker mRctFlagsChecker; bool mUseRctFlags = false; diff --git a/PWGCF/Femto/Core/collisionHistManager.h b/PWGCF/Femto/Core/collisionHistManager.h index 76f3efb9de4..108921ba5be 100644 --- a/PWGCF/Femto/Core/collisionHistManager.h +++ b/PWGCF/Femto/Core/collisionHistManager.h @@ -52,8 +52,8 @@ enum ColHist { kCentVsSphericity, kMultVsSphericity, // mc - kMcCent, - kMcMult, + kTrueCentVsCent, + kTrueMultVsMult, kColHistLast }; @@ -80,8 +80,8 @@ constexpr std::array, kColHistLast> HistTable = { {kMultVsSphericity, o2::framework::kTH2F, "hMultVsSphericity", "Multiplicity vs Sphericity; Multiplicity; Sphericity"}, {kCentVsSphericity, o2::framework::kTH2F, "hCentVsSphericity", "Centrality vs Sphericity; Centrality (%); Sphericity"}, // mc - {kMcCent, o2::framework::kTH1F, "hMcCent", "Monte Carlo Centrality; Centrality (%); Entries"}, - {kMcMult, o2::framework::kTH1F, "hMcMult", "Monte Carlo Multiplicity; Multiplicity; Entries"}, + {kTrueCentVsCent, o2::framework::kTH2F, "hTrueCentVsCent", "True centrality vs centrality; Centrality_{True} (%); Centrality (%)"}, + {kTrueMultVsMult, o2::framework::kTH2F, "hTrueMultVsMult", "True multiplicity vs multiplicity; Multiplicity_{True}; Multiplicity"}, }}; #define COL_HIST_ANALYSIS_MAP(conf) \ @@ -102,9 +102,9 @@ constexpr std::array, kColHistLast> HistTable = { {kMultVsSphericity, {confAnalysis.mult, confQa.sphericity}}, \ {kCentVsSphericity, {confBinningAnalysis.cent, confQa.sphericity}}, -#define COL_HIST_MC_QA_MAP(conf) \ - {kMcMult, {conf.mult}}, \ - {kMcCent, {conf.cent}}, +#define COL_HIST_MC_MAP(conf) \ + {kTrueMultVsMult, {conf.mult, conf.mult}}, \ + {kTrueCentVsCent, {conf.cent, conf.cent}}, template auto makeColHistSpecMap(const T& confBinningAnalysis) @@ -113,6 +113,14 @@ auto makeColHistSpecMap(const T& confBinningAnalysis) COL_HIST_ANALYSIS_MAP(confBinningAnalysis)}; } +template +auto makeColMcHistSpecMap(const T& confBinningAnalysis) +{ + return std::map>{ + COL_HIST_ANALYSIS_MAP(confBinningAnalysis) + COL_HIST_MC_MAP(confBinningAnalysis)}; +} + template auto makeColQaHistSpecMap(const T1& confBinningAnalysis, const T2& confBinningQa) { @@ -127,12 +135,12 @@ auto makeColMcQaHistSpecMap(const T1& confBinningAnalysis, const T2& confBinning return std::map>{ COL_HIST_ANALYSIS_MAP(confBinningAnalysis) COL_HIST_QA_MAP(confBinningAnalysis, confBinningQa) - COL_HIST_MC_QA_MAP(confBinningAnalysis)}; + COL_HIST_MC_MAP(confBinningAnalysis)}; } #undef COL_HIST_ANALYSIS_MAP #undef COL_HIST_QA_MAP -#undef COL_HIST_MC_QA_MAP +#undef COL_HIST_MC_MAP struct ConfCollisionBinning : o2::framework::ConfigurableGroup { std::string prefix = std::string("CollisionBinning"); @@ -240,8 +248,8 @@ class CollisionHistManager void initMc(std::map> const& Specs) { std::string mcDir = std::string(ColMcDir); - mHistogramRegistry->add(mcDir + getHistNameV2(kMcMult, HistTable), getHistDesc(kMcMult, HistTable), getHistType(kMcMult, HistTable), {Specs.at(kMcMult)}); - mHistogramRegistry->add(mcDir + getHistNameV2(kMcCent, HistTable), getHistDesc(kMcCent, HistTable), getHistType(kMcCent, HistTable), {Specs.at(kMcCent)}); + mHistogramRegistry->add(mcDir + getHistNameV2(kTrueMultVsMult, HistTable), getHistDesc(kTrueMultVsMult, HistTable), getHistType(kTrueMultVsMult, HistTable), {Specs.at(kTrueMultVsMult)}); + mHistogramRegistry->add(mcDir + getHistNameV2(kTrueCentVsCent, HistTable), getHistDesc(kTrueCentVsCent, HistTable), getHistType(kTrueCentVsCent, HistTable), {Specs.at(kTrueCentVsCent)}); } template @@ -273,9 +281,12 @@ class CollisionHistManager template void fillMc(T1 const& col, T2 const& /*mcCols*/) { - auto genCol = col.template fMcCol_as(); - mHistogramRegistry->fill(HIST(ColMcDir) + HIST(getHistName(kMcMult, HistTable)), genCol.multMc()); - mHistogramRegistry->fill(HIST(ColMcDir) + HIST(getHistName(kMcCent, HistTable)), genCol.centMc()); + if (!col.has_fMcCol()) { + return; + } + auto mcCol = col.template fMcCol_as(); + mHistogramRegistry->fill(HIST(ColMcDir) + HIST(getHistName(kTrueMultVsMult, HistTable)), mcCol.mult(), col.mult()); + mHistogramRegistry->fill(HIST(ColMcDir) + HIST(getHistName(kTrueCentVsCent, HistTable)), mcCol.cent(), col.cent()); } o2::framework::HistogramRegistry* mHistogramRegistry = nullptr; diff --git a/PWGCF/Femto/Core/kinkBuilder.h b/PWGCF/Femto/Core/kinkBuilder.h index e4897d88dcd..2e427f79287 100644 --- a/PWGCF/Femto/Core/kinkBuilder.h +++ b/PWGCF/Femto/Core/kinkBuilder.h @@ -94,16 +94,16 @@ struct ConfSigmaPlusBits : o2::framework::ConfigurableGroup { #undef KINK_DEFAULT_BITS // base selection for analysis task for kinks -#define KINK_DEFAULT_SELECTIONS(defaultMassMin, defaultMassMax, defaultPdgCode) \ - o2::framework::Configurable pdgCode{"pdgCode", defaultPdgCode, "Kink PDG code"}; \ - o2::framework::Configurable ptMin{"ptMin", 0.f, "Minimum pT"}; \ - o2::framework::Configurable ptMax{"ptMax", 999.f, "Maximum pT"}; \ - o2::framework::Configurable etaMin{"etaMin", -10.f, "Minimum eta"}; \ - o2::framework::Configurable etaMax{"etaMax", 10.f, "Maximum eta"}; \ - o2::framework::Configurable phiMin{"phiMin", 0.f, "Minimum phi"}; \ - o2::framework::Configurable phiMax{"phiMax", 1.f * o2::constants::math::TwoPI, "Maximum phi"}; \ - o2::framework::Configurable massMin{"massMin", defaultMassMin, "Minimum invariant mass for Sigma"}; \ - o2::framework::Configurable massMax{"massMax", defaultMassMax, "Maximum invariant mass for Sigma"}; \ +#define KINK_DEFAULT_SELECTIONS(defaultMassMin, defaultMassMax, defaultPdgCode) \ + o2::framework::Configurable pdgCodeAbs{"pdgCodeAbs", defaultPdgCode, "PDG code. Select antipartilce via sign"}; \ + o2::framework::Configurable ptMin{"ptMin", 0.f, "Minimum pT"}; \ + o2::framework::Configurable ptMax{"ptMax", 999.f, "Maximum pT"}; \ + o2::framework::Configurable etaMin{"etaMin", -10.f, "Minimum eta"}; \ + o2::framework::Configurable etaMax{"etaMax", 10.f, "Maximum eta"}; \ + o2::framework::Configurable phiMin{"phiMin", 0.f, "Minimum phi"}; \ + o2::framework::Configurable phiMax{"phiMax", 1.f * o2::constants::math::TwoPI, "Maximum phi"}; \ + o2::framework::Configurable massMin{"massMin", defaultMassMin, "Minimum invariant mass for Sigma"}; \ + o2::framework::Configurable massMax{"massMax", defaultMassMax, "Maximum invariant mass for Sigma"}; \ o2::framework::Configurable mask{"mask", 0x0, "Bitmask for kink selection"}; // base selection for analysis task for sigmas @@ -498,11 +498,11 @@ class KinkBuilder if constexpr (modes::isEqual(kinkType, modes::Kink::kSigma)) { fillSigma(collisionProducts, kinkProducts, kink, daughterIndex); - mcBuilder.template fillMcSigmaWithLabel(col, mcCols, kink, daughter, mcParticles, mcProducts); + mcBuilder.template fillMcSigmaWithLabel(col, mcCols, daughter, mcParticles, mcProducts); } if constexpr (modes::isEqual(kinkType, modes::Kink::kSigmaPlus)) { fillSigmaPlus(collisionProducts, kinkProducts, kink, daughterIndex); - mcBuilder.template fillMcSigmaPlusWithLabel(col, mcCols, kink, daughter, mcParticles, mcProducts); + mcBuilder.template fillMcSigmaPlusWithLabel(col, mcCols, daughter, mcParticles, mcProducts); } } } diff --git a/PWGCF/Femto/Core/kinkHistManager.h b/PWGCF/Femto/Core/kinkHistManager.h index e6621f01491..6a3c00c189a 100644 --- a/PWGCF/Femto/Core/kinkHistManager.h +++ b/PWGCF/Femto/Core/kinkHistManager.h @@ -72,6 +72,7 @@ enum KinkHist { kPrimary, kFromWrongCollision, kFromMaterial, + kMissidentified, kSecondary1, kSecondary2, kSecondary3, @@ -87,7 +88,8 @@ constexpr std::size_t MaxSecondary = 3; o2::framework::ConfigurableAxis eta{"eta", {{300, -1.5, 1.5}}, "Eta"}; \ o2::framework::ConfigurableAxis phi{"phi", {{720, 0, 1.f * o2::constants::math::TwoPI}}, "Phi"}; \ o2::framework::ConfigurableAxis mass{"mass", {{200, defaultMassMin, defaultMassMax}}, "Mass"}; \ - o2::framework::ConfigurableAxis sign{"sign", {{3, -1.5, 1.5}}, "Sign"}; + o2::framework::ConfigurableAxis sign{"sign", {{3, -1.5, 1.5}}, "Sign"}; \ + o2::framework::ConfigurableAxis pdgCodes{"pdgCodes", {{8001, -4000.5, 4000.5}}, "PDG codes of selected V0s"}; template struct ConfSigmaBinning : o2::framework::ConfigurableGroup { @@ -116,6 +118,8 @@ struct ConfKinkQaBinning : o2::framework::ConfigurableGroup { o2::framework::ConfigurableAxis dcaDaugToPV{"dcaDaugToPV", {{1000, 0, 100}}, "Daughter DCA to PV (cm)"}; o2::framework::ConfigurableAxis decayVertex{"decayVertex", {{100, 0, 100}}, "Decay vertex position (cm)"}; o2::framework::ConfigurableAxis transRadius{"transRadius", {{100, 0, 100}}, "Transverse radius (cm)"}; + o2::framework::Configurable plotOrigins{"plotOrigins", true, "MC ONLY: Plot pt vs cosPa for different particle origins"}; + o2::framework::Configurable> pdgCodesForMothersOfSecondary{"pdgCodesForMothersOfSecondary", {3312, 3334}, "MC ONLY: PDG codes of mothers of secondaries (Max 3 will be considered)"}; }; constexpr const char PrefixSigmaQaBinning1[] = "SigmaQaBinning1"; @@ -124,23 +128,6 @@ using ConfSigmaQaBinning1 = ConfKinkQaBinning; constexpr const char PrefixSigmaPlusQaBinning1[] = "SigmaPlusQaBinning1"; using ConfSigmaPlusQaBinning1 = ConfKinkQaBinning; -template -struct ConfKinkMcBinning : o2::framework::ConfigurableGroup { - std::string prefix = std::string(Prefix); - o2::framework::ConfigurableAxis pdgCodes{"pdgCodes", {{8001, -4000.5, 4000.5}}, "PDG codes of selected V0s"}; - o2::framework::ConfigurableAxis statusCode{"statusCode", {{21, -0.5, 20.5}}, "Status codes (i.e. Origin)"}; - o2::framework::Configurable plotOrigins{"plotOrigins", true, "Plot pt vs cosPa for different particle origins"}; - o2::framework::Configurable> pdgCodesForMothersOfSecondary{"pdgCodesForMothersOfSecondary", {3312, 3334}, "PDG codes of mothers of secondaries (Max 3 will be considered)"}; - o2::framework::ConfigurableAxis pt{"pt", {{150, 0, 3}}, "Pt"}; - o2::framework::ConfigurableAxis kinkAngle{"kinkAngle", {{315, 0, 3.15}}, "Kink angle"}; -}; - -constexpr const char PrefixSigmaMcBinning1[] = "SigmaMcBinning1"; -using ConfSigmaMcBinning = ConfKinkMcBinning; - -constexpr const char PrefixSigmaPlusMcBinning1[] = "SigmaPlusMcBinning1"; -using ConfSigmaPlusMcBinning = ConfKinkMcBinning; - // must be in sync with enum KinkHist // the enum gives the correct index in the array constexpr std::array, kKinkHistLast> HistTable = { @@ -170,13 +157,14 @@ constexpr std::array, kKinkHistLast> HistTable = {kTrueEta, o2::framework::kTH1F, "hTrueEta", "True pseudorapdity; #eta; Entries"}, {kTruePhi, o2::framework::kTH1F, "hTruePhi", "True azimuthal angle; #varphi; Entries"}, {kNoMcParticle, o2::framework::kTH2F, "hNoMcParticle", "Wrongly reconstructed particles; p_{T} (GeV/#it{c}); cos(#alpha)"}, - {kPrimary, o2::framework::kTH2F, "hPrimary", "Primary particles; p_{T} (GeV/#it{c}); cos(#alpha)"}, - {kFromWrongCollision, o2::framework::kTH2F, "hFromWrongCollision", "Particles associated to wrong collision; p_{T} (GeV/#it{c}); cos(#alpha)"}, - {kFromMaterial, o2::framework::kTH2F, "hFromMaterial", "Particles from material; p_{T} (GeV/#it{c}); cos(#alpha)"}, - {kSecondary1, o2::framework::kTH2F, "hFromSecondary1", "Particles from secondary decay; p_{T} (GeV/#it{c}); cos(#alpha)"}, - {kSecondary2, o2::framework::kTH2F, "hFromSecondary2", "Particles from seconary decay; p_{T} (GeV/#it{c}); cos(#alpha)"}, - {kSecondary3, o2::framework::kTH2F, "hFromSecondary3", "Particles from seconary decay; p_{T} (GeV/#it{c}); cos(#alpha)"}, - {kSecondaryOther, o2::framework::kTH2F, "hFromSecondaryOther", "Particles from every other seconary decay; p_{T} (GeV/#it{c}); cos(#alpha)"}}}; + {kPrimary, o2::framework::kTH2F, "hPrimary", "Primary particles; p_{T} (GeV/#it{c}); kink angle"}, + {kFromWrongCollision, o2::framework::kTH2F, "hFromWrongCollision", "Particles associated to wrong collision; p_{T} (GeV/#it{c}); kink angle"}, + {kFromMaterial, o2::framework::kTH2F, "hFromMaterial", "Particles from material; p_{T} (GeV/#it{c}); kink angle"}, + {kMissidentified, o2::framework::kTH2F, "hMissidentified", "Missidentified particles (fake/wrong PDG code); p_{T} (GeV/#it{c}); kink angle"}, + {kSecondary1, o2::framework::kTH2F, "hFromSecondary1", "Particles from secondary decay; p_{T} (GeV/#it{c}); kink angle"}, + {kSecondary2, o2::framework::kTH2F, "hFromSecondary2", "Particles from seconary decay; p_{T} (GeV/#it{c}); kink angle"}, + {kSecondary3, o2::framework::kTH2F, "hFromSecondary3", "Particles from seconary decay; p_{T} (GeV/#it{c}); kink angle"}, + {kSecondaryOther, o2::framework::kTH2F, "hFromSecondaryOther", "Particles from every other seconary decay; p_{T} (GeV/#it{c}); kink angle"}}}; #define KINK_HIST_ANALYSIS_MAP(conf) \ {kPt, {conf.pt}}, \ @@ -185,6 +173,14 @@ constexpr std::array, kKinkHistLast> HistTable = {kMass, {conf.mass}}, \ {kSign, {conf.sign}}, +#define KINK_HIST_MC_MAP(conf) \ + {kTruePt, {conf.pt}}, \ + {kTrueEta, {conf.eta}}, \ + {kTruePhi, {conf.phi}}, \ + {kPdg, {conf.pdgCodes}}, \ + {kPdgMother, {conf.pdgCodes}}, \ + {kPdgPartonicMother, {conf.pdgCodes}}, + #define KINK_HIST_QA_MAP(confAnalysis, confQa) \ {kKinkAngle, {confQa.kinkAngle}}, \ {kDcaMothToPV, {confQa.dcaMothToPV}}, \ @@ -200,22 +196,16 @@ constexpr std::array, kKinkHistLast> HistTable = {kPtVsKinkAngle, {confAnalysis.pt, confQa.kinkAngle}}, \ {kPtVsDecayRadius, {confAnalysis.pt, confQa.transRadius}}, -#define KINK_HIST_MC_MAP(confAnalysis, confMc) \ - {kTruePt, {confAnalysis.pt}}, \ - {kTrueEta, {confAnalysis.eta}}, \ - {kTruePhi, {confAnalysis.phi}}, \ - {kOrigin, {confMc.statusCode}}, \ - {kPdg, {confMc.pdgCodes}}, \ - {kPdgMother, {confMc.pdgCodes}}, \ - {kPdgPartonicMother, {confMc.pdgCodes}}, \ - {kNoMcParticle, {confMc.pt, confMc.kinkAngle}}, \ - {kPrimary, {confMc.pt, confMc.kinkAngle}}, \ - {kFromWrongCollision, {confMc.pt, confMc.kinkAngle}}, \ - {kFromMaterial, {confMc.pt, confMc.kinkAngle}}, \ - {kSecondary1, {confMc.pt, confMc.kinkAngle}}, \ - {kSecondary2, {confMc.pt, confMc.kinkAngle}}, \ - {kSecondary3, {confMc.pt, confMc.kinkAngle}}, \ - {kSecondaryOther, {confMc.pt, confMc.kinkAngle}}, +#define KINK_HIST_MC_QA_MAP(confAnalysis, confQa) \ + {kNoMcParticle, {confAnalysis.pt, confQa.kinkAngle}}, \ + {kPrimary, {confAnalysis.pt, confQa.kinkAngle}}, \ + {kFromWrongCollision, {confAnalysis.pt, confQa.kinkAngle}}, \ + {kFromMaterial, {confAnalysis.pt, confQa.kinkAngle}}, \ + {kMissidentified, {confAnalysis.pt, confQa.kinkAngle}}, \ + {kSecondary1, {confAnalysis.pt, confQa.kinkAngle}}, \ + {kSecondary2, {confAnalysis.pt, confQa.kinkAngle}}, \ + {kSecondary3, {confAnalysis.pt, confQa.kinkAngle}}, \ + {kSecondaryOther, {confAnalysis.pt, confQa.kinkAngle}}, template auto makeKinkHistSpecMap(const T& confBinningAnalysis) @@ -224,6 +214,14 @@ auto makeKinkHistSpecMap(const T& confBinningAnalysis) KINK_HIST_ANALYSIS_MAP(confBinningAnalysis)}; } +template +auto makeKinkMcHistSpecMap(const T& confBinningAnalysis) +{ + return std::map>{ + KINK_HIST_ANALYSIS_MAP(confBinningAnalysis) + KINK_HIST_MC_MAP(confBinningAnalysis)}; +} + template std::map> makeKinkQaHistSpecMap(T1 const& confBinningAnalysis, T2 const& confBinningQa) { @@ -232,18 +230,20 @@ std::map> makeKinkQaHistSpecMap(T1 co KINK_HIST_QA_MAP(confBinningAnalysis, confBinningQa)}; } -template -std::map> makeKinkMcQaHistSpecMap(T1 const& confBinningAnalysis, T2 const& confBinningQa, T3 const& confBinningMc) +template +std::map> makeKinkMcQaHistSpecMap(T1 const& confBinningAnalysis, T2 const& confBinningQa) { return std::map>{ KINK_HIST_ANALYSIS_MAP(confBinningAnalysis) KINK_HIST_QA_MAP(confBinningAnalysis, confBinningQa) - KINK_HIST_MC_MAP(confBinningAnalysis, confBinningMc)}; + KINK_HIST_MC_MAP(confBinningAnalysis) + KINK_HIST_MC_QA_MAP(confBinningAnalysis, confBinningQa)}; } #undef KINK_HIST_ANALYSIS_MAP #undef KINK_HIST_QA_MAP #undef KINK_HIST_MC_MAP +#undef KINK_HIST_MC_QA_MAP constexpr char PrefixSigmaQa[] = "SigmaQA/"; constexpr char PrefixSigma1[] = "Sigma1/"; @@ -267,68 +267,98 @@ class KinkHistManager KinkHistManager() = default; ~KinkHistManager() = default; - // init for analysis - template + // init for analysis and mc + template void init(o2::framework::HistogramRegistry* registry, std::map> const& KinkSpecs, + T const& ConfKinkSelection, std::map> const& ChaDauSpecs) { mHistogramRegistry = registry; - mChaDauManager.template init(registry, ChaDauSpecs, AbsChargeDaughters); - if constexpr (isFlagSet(mode, modes::Mode::kAnalysis)) { - initAnalysis(KinkSpecs); - } - if constexpr (isFlagSet(mode, modes::Mode::kQa)) { - initQa(KinkSpecs); - } - if constexpr (isFlagSet(mode, modes::Mode::kMc)) { - initMc(KinkSpecs); + mPdgCode = std::abs(ConfKinkSelection.pdgCodeAbs.value) * ConfKinkSelection.sign.value; + + int chaDauPdgCodeAbs = 0; + int chaDauCharge = 0; + const int absCharge = 1; + + if (std::abs(mPdgCode) == PDG_t::kSigmaMinus) { + if (ConfKinkSelection.sign.value < 0) { + chaDauPdgCodeAbs = std::abs(PDG_t::kPiMinus); + chaDauCharge = -1; + } else { + chaDauPdgCodeAbs = std::abs(PDG_t::kPiPlus); + chaDauCharge = 1; + } + } else if (std::abs(mPdgCode) == PDG_t::kSigmaPlus) { + if (ConfKinkSelection.sign.value > 0) { + chaDauPdgCodeAbs = std::abs(PDG_t::kProton); + chaDauCharge = 1; + } else { + chaDauPdgCodeAbs = std::abs(PDG_t::kProtonBar); + chaDauCharge = -1; + } + } else { + LOG(fatal) << "PDG code for Kink has to be either SigmaMinus or SigmaPlus"; } - } - // init for qa - template - void init(o2::framework::HistogramRegistry* registry, - std::map> const& KinkSpecs, - T1 const& KinkConfBinningQa, - std::map> const& ChaDauSpecs, - T2 const& ChaDauConfBinningQa) - { - mHistogramRegistry = registry; - mChaDauManager.template init(registry, ChaDauSpecs, ChaDauConfBinningQa, AbsChargeDaughters); - this->enableOptionalHistograms(KinkConfBinningQa); + mChaDauManager.template init(registry, ChaDauSpecs, absCharge, chaDauCharge, chaDauPdgCodeAbs); if constexpr (isFlagSet(mode, modes::Mode::kAnalysis)) { - initAnalysis(KinkSpecs); + this->initAnalysis(KinkSpecs); } if constexpr (isFlagSet(mode, modes::Mode::kQa)) { - initQa(KinkSpecs); + this->initQa(KinkSpecs); } if constexpr (isFlagSet(mode, modes::Mode::kMc)) { - initMc(KinkSpecs); + this->initMc(KinkSpecs); } } - // init for qa + mc - template + // init for analysis and qa and and mc + template void init(o2::framework::HistogramRegistry* registry, std::map> const& KinkSpecs, - T1 const& KinkConfBinningQa, - T2 const& KinkConfBinningMc, + T1 const& ConfKinkSelection, + T2 const& ConfKinkBinningQa, std::map> const& ChaDauSpecs, - T3 const& ChaDauConfBinningQa, - T4 const& ChaDauConfBinningMc) + T3 const& ConfChaDauBinningQa) { mHistogramRegistry = registry; - mChaDauManager.template init(registry, ChaDauSpecs, ChaDauConfBinningQa, ChaDauConfBinningMc, AbsChargeDaughters); - this->enableOptionalHistograms(KinkConfBinningQa, KinkConfBinningMc); + mPdgCode = std::abs(ConfKinkSelection.pdgCodeAbs.value) * ConfKinkSelection.sign.value; + this->enableOptionalHistograms(ConfKinkBinningQa); + + int chaDauPdgCodeAbs = 0; + int chaDauCharge = 0; + const int absCharge = 1; + + if (std::abs(mPdgCode) == PDG_t::kSigmaMinus) { + if (ConfKinkSelection.sign.value < 0) { + chaDauPdgCodeAbs = std::abs(PDG_t::kPiMinus); + chaDauCharge = -1; + } else { + chaDauPdgCodeAbs = std::abs(PDG_t::kPiPlus); + chaDauCharge = 1; + } + } else if (std::abs(mPdgCode) == PDG_t::kSigmaPlus) { + if (ConfKinkSelection.sign.value > 0) { + chaDauPdgCodeAbs = std::abs(PDG_t::kProton); + chaDauCharge = 1; + } else { + chaDauPdgCodeAbs = std::abs(PDG_t::kProtonBar); + chaDauCharge = -1; + } + } else { + LOG(fatal) << "PDG code for Kink has to be either SigmaMinus or SigmaPlus"; + } + + mChaDauManager.template init(registry, ChaDauSpecs, absCharge, chaDauCharge, chaDauPdgCodeAbs, ConfChaDauBinningQa); if constexpr (isFlagSet(mode, modes::Mode::kAnalysis)) { - initAnalysis(KinkSpecs); + this->initAnalysis(KinkSpecs); } if constexpr (isFlagSet(mode, modes::Mode::kQa)) { - initQa(KinkSpecs); + this->initQa(KinkSpecs); } if constexpr (isFlagSet(mode, modes::Mode::kMc)) { - initMc(KinkSpecs); + this->initMc(KinkSpecs); } } @@ -340,10 +370,10 @@ class KinkHistManager auto chaDaughter = tracks.rawIteratorAt(kinkCandidate.chaDauId() - tracks.offset()); mChaDauManager.template fill(chaDaughter, tracks); if constexpr (isFlagSet(mode, modes::Mode::kAnalysis)) { - fillAnalysis(kinkCandidate); + this->fillAnalysis(kinkCandidate); } if constexpr (isFlagSet(mode, modes::Mode::kQa)) { - fillQa(kinkCandidate); + this->fillQa(kinkCandidate); } } @@ -353,35 +383,28 @@ class KinkHistManager auto chaDaughter = tracks.rawIteratorAt(kinkCandidate.chaDauId() - tracks.offset()); mChaDauManager.template fill(chaDaughter, tracks, mcParticles, mcMothers, mcPartonicMothers); if constexpr (modes::isFlagSet(mode, modes::Mode::kAnalysis)) { - fillAnalysis(kinkCandidate); + this->fillAnalysis(kinkCandidate); } if constexpr (modes::isFlagSet(mode, modes::Mode::kQa)) { - fillQa(kinkCandidate); + this->fillQa(kinkCandidate); } if constexpr (modes::isFlagSet(mode, modes::Mode::kMc)) { - fillMc(kinkCandidate, mcParticles, mcMothers, mcPartonicMothers); + this->template fillMc(kinkCandidate, mcParticles, mcMothers, mcPartonicMothers); } } private: template - void enableOptionalHistograms(T1 const& KinkConfBinningQa) + void enableOptionalHistograms(T1 const& ConfKinkBinningQa) { - mPlot2d = KinkConfBinningQa.plot2d.value; - } - - // for qa and mc - template - void enableOptionalHistograms(T1 const& V0ConfBinningQa, T2 const& V0ConfBinningMc) - { - this->enableOptionalHistograms(V0ConfBinningQa); + mPlot2d = ConfKinkBinningQa.plot2d.value; - mPlotOrigins = V0ConfBinningMc.plotOrigins.value; - mPlotNSecondaries = V0ConfBinningMc.pdgCodesForMothersOfSecondary.value.size(); + mPlotOrigins = ConfKinkBinningQa.plotOrigins.value; + mPlotNSecondaries = ConfKinkBinningQa.pdgCodesForMothersOfSecondary.value.size(); for (std::size_t i = 0; i < MaxSecondary; i++) { - if (i < V0ConfBinningMc.pdgCodesForMothersOfSecondary.value.size()) { - mPdgCodesSecondaryMother.at(i) = std::abs(V0ConfBinningMc.pdgCodesForMothersOfSecondary.value.at(i)); + if (i < ConfKinkBinningQa.pdgCodesForMothersOfSecondary.value.size()) { + mPdgCodesSecondaryMother.at(i) = std::abs(ConfKinkBinningQa.pdgCodesForMothersOfSecondary.value.at(i)); } else { mPdgCodesSecondaryMother.at(i) = 0; } @@ -425,7 +448,17 @@ class KinkHistManager mHistogramRegistry->add(mcDir + getHistNameV2(kTruePt, HistTable), getHistDesc(kTruePt, HistTable), getHistType(kTruePt, HistTable), {KinkSpecs.at(kTruePt)}); mHistogramRegistry->add(mcDir + getHistNameV2(kTrueEta, HistTable), getHistDesc(kTrueEta, HistTable), getHistType(kTrueEta, HistTable), {KinkSpecs.at(kTrueEta)}); mHistogramRegistry->add(mcDir + getHistNameV2(kTruePhi, HistTable), getHistDesc(kTruePhi, HistTable), getHistType(kTruePhi, HistTable), {KinkSpecs.at(kTruePhi)}); - mHistogramRegistry->add(mcDir + getHistNameV2(kOrigin, HistTable), getHistDesc(kOrigin, HistTable), getHistType(kOrigin, HistTable), {KinkSpecs.at(kOrigin)}); + + // mc origin can be configured here + const framework::AxisSpec axisOrigin = {static_cast(modes::McOrigin::kMcOriginLast), -0.5, static_cast(modes::McOrigin::kMcOriginLast) - 0.5}; + mHistogramRegistry->add(mcDir + getHistNameV2(kOrigin, HistTable), getHistDesc(kOrigin, HistTable), getHistType(kOrigin, HistTable), {axisOrigin}); + mHistogramRegistry->get(HIST(kinkPrefix) + HIST(McDir) + HIST(histmanager::getHistName(kOrigin, HistTable)))->GetXaxis()->SetBinLabel(1 + static_cast(modes::McOrigin::kNoMcParticle), modes::mcOriginToString(modes::McOrigin::kNoMcParticle)); + mHistogramRegistry->get(HIST(kinkPrefix) + HIST(McDir) + HIST(histmanager::getHistName(kOrigin, HistTable)))->GetXaxis()->SetBinLabel(1 + static_cast(modes::McOrigin::kFromWrongCollision), modes::mcOriginToString(modes::McOrigin::kFromWrongCollision)); + mHistogramRegistry->get(HIST(kinkPrefix) + HIST(McDir) + HIST(histmanager::getHistName(kOrigin, HistTable)))->GetXaxis()->SetBinLabel(1 + static_cast(modes::McOrigin::kPhysicalPrimary), modes::mcOriginToString(modes::McOrigin::kPhysicalPrimary)); + mHistogramRegistry->get(HIST(kinkPrefix) + HIST(McDir) + HIST(histmanager::getHistName(kOrigin, HistTable)))->GetXaxis()->SetBinLabel(1 + static_cast(modes::McOrigin::kFromSecondaryDecay), modes::mcOriginToString(modes::McOrigin::kFromSecondaryDecay)); + mHistogramRegistry->get(HIST(kinkPrefix) + HIST(McDir) + HIST(histmanager::getHistName(kOrigin, HistTable)))->GetXaxis()->SetBinLabel(1 + static_cast(modes::McOrigin::kFromMaterial), modes::mcOriginToString(modes::McOrigin::kFromMaterial)); + mHistogramRegistry->get(HIST(kinkPrefix) + HIST(McDir) + HIST(histmanager::getHistName(kOrigin, HistTable)))->GetXaxis()->SetBinLabel(1 + static_cast(modes::McOrigin::kMissidentified), modes::mcOriginToString(modes::McOrigin::kMissidentified)); + mHistogramRegistry->add(mcDir + getHistNameV2(kPdg, HistTable), getHistDesc(kPdg, HistTable), getHistType(kPdg, HistTable), {KinkSpecs.at(kPdg)}); mHistogramRegistry->add(mcDir + getHistNameV2(kPdgMother, HistTable), getHistDesc(kPdgMother, HistTable), getHistType(kPdgMother, HistTable), {KinkSpecs.at(kPdgMother)}); mHistogramRegistry->add(mcDir + getHistNameV2(kPdgPartonicMother, HistTable), getHistDesc(kPdgPartonicMother, HistTable), getHistType(kPdgPartonicMother, HistTable), {KinkSpecs.at(kPdgPartonicMother)}); @@ -435,6 +468,7 @@ class KinkHistManager mHistogramRegistry->add(mcDir + getHistNameV2(kPrimary, HistTable), getHistDesc(kPrimary, HistTable), getHistType(kPrimary, HistTable), {KinkSpecs.at(kPrimary)}); mHistogramRegistry->add(mcDir + getHistNameV2(kFromWrongCollision, HistTable), getHistDesc(kFromWrongCollision, HistTable), getHistType(kFromWrongCollision, HistTable), {KinkSpecs.at(kFromWrongCollision)}); mHistogramRegistry->add(mcDir + getHistNameV2(kFromMaterial, HistTable), getHistDesc(kFromMaterial, HistTable), getHistType(kFromMaterial, HistTable), {KinkSpecs.at(kFromMaterial)}); + mHistogramRegistry->add(mcDir + getHistNameV2(kMissidentified, HistTable), getHistDesc(kMissidentified, HistTable), getHistType(kMissidentified, HistTable), {KinkSpecs.at(kMissidentified)}); if (mPlotNSecondaries >= histmanager::kSecondaryPlotLevel1) { mHistogramRegistry->add(mcDir + getHistNameV2(kSecondary1, HistTable), getHistDesc(kSecondary1, HistTable), getHistType(kSecondary1, HistTable), {KinkSpecs.at(kSecondary1)}); @@ -487,28 +521,39 @@ class KinkHistManager } } - template + template void fillMc(T1 const& kinkCandidate, T2 const& /*mcParticles*/, T3 const& /*mcMothers*/, T4 const& /*mcPartonicMothers*/) { // No MC Particle if (!kinkCandidate.has_fMcParticle()) { mHistogramRegistry->fill(HIST(kinkPrefix) + HIST(McDir) + HIST(getHistName(kPdg, HistTable)), 0); mHistogramRegistry->fill(HIST(kinkPrefix) + HIST(McDir) + HIST(getHistName(kOrigin, HistTable)), static_cast(modes::McOrigin::kNoMcParticle)); - if (mPlotOrigins) { - mHistogramRegistry->fill(HIST(kinkPrefix) + HIST(McDir) + HIST(getHistName(kNoMcParticle, HistTable)), kinkCandidate.pt(), kinkCandidate.kinkAngle()); + if constexpr (modes::isFlagSet(mode, modes::Mode::kQa)) { + if (mPlotOrigins) { + mHistogramRegistry->fill(HIST(kinkPrefix) + HIST(McDir) + HIST(getHistName(kNoMcParticle, HistTable)), kinkCandidate.pt(), kinkCandidate.kinkAngle()); + } } return; } // Retrieve MC particle auto mcParticle = kinkCandidate.template fMcParticle_as(); + + // missidentifed particles are special case + // whether a particle is missidentfied or not cannot be known by the producer so we check it here + bool isMissidentified = mcParticle.pdgCode() != mPdgCode; + mHistogramRegistry->fill(HIST(kinkPrefix) + HIST(McDir) + HIST(getHistName(kTruePt, HistTable)), mcParticle.pt()); mHistogramRegistry->fill(HIST(kinkPrefix) + HIST(McDir) + HIST(getHistName(kTrueEta, HistTable)), mcParticle.eta()); mHistogramRegistry->fill(HIST(kinkPrefix) + HIST(McDir) + HIST(getHistName(kTruePhi, HistTable)), mcParticle.phi()); - mHistogramRegistry->fill(HIST(kinkPrefix) + HIST(McDir) + HIST(getHistName(kOrigin, HistTable)), mcParticle.origin()); + if (isMissidentified) { + mHistogramRegistry->fill(HIST(kinkPrefix) + HIST(McDir) + HIST(getHistName(kOrigin, HistTable)), static_cast(modes::McOrigin::kMissidentified)); + } else { + mHistogramRegistry->fill(HIST(kinkPrefix) + HIST(McDir) + HIST(getHistName(kOrigin, HistTable)), mcParticle.origin()); + } mHistogramRegistry->fill(HIST(kinkPrefix) + HIST(McDir) + HIST(getHistName(kPdg, HistTable)), mcParticle.pdgCode()); - // Mother PDG + // get mother if (kinkCandidate.has_fMcMother()) { auto mother = kinkCandidate.template fMcMother_as(); mHistogramRegistry->fill(HIST(kinkPrefix) + HIST(McDir) + HIST(getHistName(kPdgMother, HistTable)), mother.pdgCode()); @@ -516,7 +561,7 @@ class KinkHistManager mHistogramRegistry->fill(HIST(kinkPrefix) + HIST(McDir) + HIST(getHistName(kPdgMother, HistTable)), 0); } - // Partonic Mother PDG + // get partonic mother if (kinkCandidate.has_fMcPartMoth()) { auto partonicMother = kinkCandidate.template fMcPartMoth_as(); mHistogramRegistry->fill(HIST(kinkPrefix) + HIST(McDir) + HIST(getHistName(kPdgPartonicMother, HistTable)), partonicMother.pdgCode()); @@ -524,43 +569,52 @@ class KinkHistManager mHistogramRegistry->fill(HIST(kinkPrefix) + HIST(McDir) + HIST(getHistName(kPdgPartonicMother, HistTable)), 0); } - // Plot origins - if (mPlotOrigins) { - switch (static_cast(mcParticle.origin())) { - case modes::McOrigin::kPhysicalPrimary: - mHistogramRegistry->fill(HIST(kinkPrefix) + HIST(McDir) + HIST(getHistName(kPrimary, HistTable)), kinkCandidate.pt(), kinkCandidate.kinkAngle()); - break; - case modes::McOrigin::kFromWrongCollision: - mHistogramRegistry->fill(HIST(kinkPrefix) + HIST(McDir) + HIST(getHistName(kFromWrongCollision, HistTable)), kinkCandidate.pt(), kinkCandidate.kinkAngle()); - break; - case modes::McOrigin::kFromMaterial: - mHistogramRegistry->fill(HIST(kinkPrefix) + HIST(McDir) + HIST(getHistName(kFromMaterial, HistTable)), kinkCandidate.pt(), kinkCandidate.kinkAngle()); - break; - case modes::McOrigin::kFromSecondaryDecay: - if (kinkCandidate.has_fMcMother()) { - auto mother = kinkCandidate.template fMcMother_as(); - int motherPdgCode = std::abs(mother.pdgCode()); - // Switch on PDG of the mother - if (mPlotNSecondaries >= histmanager::kSecondaryPlotLevel1 && motherPdgCode == mPdgCodesSecondaryMother[0]) { - mHistogramRegistry->fill(HIST(kinkPrefix) + HIST(McDir) + HIST(getHistName(kSecondary1, HistTable)), kinkCandidate.pt(), kinkCandidate.kinkAngle()); - } else if (mPlotNSecondaries >= histmanager::kSecondaryPlotLevel2 && motherPdgCode == mPdgCodesSecondaryMother[1]) { - mHistogramRegistry->fill(HIST(kinkPrefix) + HIST(McDir) + HIST(getHistName(kSecondary2, HistTable)), kinkCandidate.pt(), kinkCandidate.kinkAngle()); - } else if (mPlotNSecondaries >= histmanager::kSecondaryPlotLevel3 && motherPdgCode == mPdgCodesSecondaryMother[2]) { - mHistogramRegistry->fill(HIST(kinkPrefix) + HIST(McDir) + HIST(getHistName(kSecondary3, HistTable)), kinkCandidate.pt(), kinkCandidate.kinkAngle()); - } else { - mHistogramRegistry->fill(HIST(kinkPrefix) + HIST(McDir) + HIST(getHistName(kSecondaryOther, HistTable)), kinkCandidate.pt(), kinkCandidate.kinkAngle()); - } + if constexpr (modes::isFlagSet(mode, modes::Mode::kQa)) { + if (mPlotOrigins) { + // check first if particle is missidentified + if (isMissidentified) { + // if it is, we fill it as such + mHistogramRegistry->fill(HIST(kinkPrefix) + HIST(McDir) + HIST(getHistName(kMissidentified, HistTable)), kinkCandidate.pt(), kinkCandidate.kinkAngle()); + } else { + // if not, we fill it acccoridng to its origin + switch (static_cast(mcParticle.origin())) { + case modes::McOrigin::kPhysicalPrimary: + mHistogramRegistry->fill(HIST(kinkPrefix) + HIST(McDir) + HIST(getHistName(kPrimary, HistTable)), kinkCandidate.pt(), kinkCandidate.kinkAngle()); + break; + case modes::McOrigin::kFromWrongCollision: + mHistogramRegistry->fill(HIST(kinkPrefix) + HIST(McDir) + HIST(getHistName(kFromWrongCollision, HistTable)), kinkCandidate.pt(), kinkCandidate.kinkAngle()); + break; + case modes::McOrigin::kFromMaterial: + mHistogramRegistry->fill(HIST(kinkPrefix) + HIST(McDir) + HIST(getHistName(kFromMaterial, HistTable)), kinkCandidate.pt(), kinkCandidate.kinkAngle()); + break; + case modes::McOrigin::kFromSecondaryDecay: + if (kinkCandidate.has_fMcMother()) { + auto mother = kinkCandidate.template fMcMother_as(); + int motherPdgCode = std::abs(mother.pdgCode()); + // Switch on PDG of the mother + if (mPlotNSecondaries >= histmanager::kSecondaryPlotLevel1 && motherPdgCode == mPdgCodesSecondaryMother[0]) { + mHistogramRegistry->fill(HIST(kinkPrefix) + HIST(McDir) + HIST(getHistName(kSecondary1, HistTable)), kinkCandidate.pt(), kinkCandidate.kinkAngle()); + } else if (mPlotNSecondaries >= histmanager::kSecondaryPlotLevel2 && motherPdgCode == mPdgCodesSecondaryMother[1]) { + mHistogramRegistry->fill(HIST(kinkPrefix) + HIST(McDir) + HIST(getHistName(kSecondary2, HistTable)), kinkCandidate.pt(), kinkCandidate.kinkAngle()); + } else if (mPlotNSecondaries >= histmanager::kSecondaryPlotLevel3 && motherPdgCode == mPdgCodesSecondaryMother[2]) { + mHistogramRegistry->fill(HIST(kinkPrefix) + HIST(McDir) + HIST(getHistName(kSecondary3, HistTable)), kinkCandidate.pt(), kinkCandidate.kinkAngle()); + } else { + mHistogramRegistry->fill(HIST(kinkPrefix) + HIST(McDir) + HIST(getHistName(kSecondaryOther, HistTable)), kinkCandidate.pt(), kinkCandidate.kinkAngle()); + } + } + break; + default: + LOG(warn) << "Encounted partilce with unknown origin!"; + break; } - break; - default: - // Unknown origin → safely ignore - break; + } } } } o2::framework::HistogramRegistry* mHistogramRegistry = nullptr; trackhistmanager::TrackHistManager mChaDauManager; + int mPdgCode = 0; bool mPlot2d = true; bool mPlotOrigins = false; int mPlotNSecondaries = 0; diff --git a/PWGCF/Femto/Core/mcBuilder.h b/PWGCF/Femto/Core/mcBuilder.h index 1a59e35f65b..bf4f023f78d 100644 --- a/PWGCF/Femto/Core/mcBuilder.h +++ b/PWGCF/Femto/Core/mcBuilder.h @@ -105,15 +105,6 @@ class McBuilder LOG(info) << "Initialization done..."; } - template - void fillMcCollision(T1 const& mcCol, T2& mcProducts) - { - mcProducts.producedMcCollisions( - mcCol.multMCNParticlesEta08(), - mcCol.centFT0M()); - mCollisionMap.emplace(mcCol.globalIndex(), mcProducts.producedMcCollisions.lastIndex()); - } - template void fillMcCollisionWithLabel(T1& mcProducts, T2 const& col, T3 const& /*mcCols*/) { @@ -139,424 +130,145 @@ class McBuilder } } - template - modes::McOrigin getOrigin(T1 const& col, T2 const& /*mcCols*/, T3 const& mcParticle) - { - // constants - const int producedByDecay = 4; - // check if reconstructed collision has a generated collision - if (col.has_mcCollision()) { - // now check collision ids, if they do not match, then the track belongs to another collision - if (col.mcCollisionId() != mcParticle.mcCollisionId()) { - return modes::McOrigin::kFromWrongCollision; - } - if (mcParticle.isPhysicalPrimary()) { - return modes::McOrigin::kPhysicalPrimary; - } else if (mcParticle.has_mothers() && mcParticle.getProcess() == producedByDecay) { - return modes::McOrigin::kFromSecondaryDecay; - } else { - // not a primary and not from a decay, we label as material - return modes::McOrigin::kFromMaterial; - } - } - return modes::McOrigin::kFromWrongCollision; - } - - template - int64_t getMcColId(T1 const& mcCol, T2& mcProducts) - { - auto gid = mcCol.globalIndex(); - // Find or create - auto it = mCollisionMap.find(gid); - if (it == mCollisionMap.end()) { - this->fillMcCollision(mcCol, mcProducts); - it = mCollisionMap.find(gid); - return it->second; - } - return it->second; - } - template void fillMcTrackWithLabel(T1 const& col, T2 const& mcCols, T3 const& track, T4 const& mcParticles, T5& mcProducts) { if (!mProduceTrackLabels) { - return; - } - - if (!track.has_mcParticle()) { mcProducts.producedTrackLabels(-1, -1, -1); return; } - - auto mcParticle = track.template mcParticle_as(); - auto mcCol = mcParticle.template mcCollision_as(); - - int64_t particleIndex = mcParticle.globalIndex(); - int64_t trackIndex = track.globalIndex(); - - int64_t mcParticleRow = -1; - - // MC PARTICLE - auto itP = mTrackToMcParticleMap.find(particleIndex); - if (itP != mTrackToMcParticleMap.end()) { - mcParticleRow = itP->second; - } else { - auto origin = this->getOrigin(col, mcCols, mcParticle); - int64_t mcColId = this->getMcColId(mcCol, mcProducts); - - mcProducts.producedMcParticles( - mcColId, - static_cast(origin), - mcParticle.pdgCode(), - mcParticle.pt() * utils::signum(mcParticle.pdgCode()), - mcParticle.eta(), - mcParticle.phi()); - - mcParticleRow = mcProducts.producedMcParticles.lastIndex(); - mTrackToMcParticleMap[particleIndex] = mcParticleRow; - } - - // mothers (fill only if exists) - int64_t mothersRow = -1; - auto itM = mTrackToMcMotherMap.find(trackIndex); - - if (itM != mTrackToMcMotherMap.end()) { - mothersRow = itM->second; - } else { - - auto mothers = mcParticle.template mothers_as(); - bool motherExists = mcParticle.has_mothers() && !mothers.empty(); - - if (motherExists) { - int motherPdg = mothers.front().pdgCode(); // PDG code is ALWAYS valid if the mother exists - - mcProducts.producedMothers(motherPdg); - mothersRow = mcProducts.producedMothers.lastIndex(); - mTrackToMcMotherMap[trackIndex] = mothersRow; - } - } - - // partonic mother (fill only if exists) - int64_t partonicRow = -1; - auto itPM = mTrackToMcPartonicMap.find(trackIndex); - - if (itPM != mTrackToMcPartonicMap.end()) { - partonicRow = itPM->second; - } else { - int partIdx = -1; - if (mcParticle.has_mothers()) { - partIdx = this->findFirstPartonicMother(mcParticle, mcParticles); - } - - bool partonicExists = (partIdx >= 0); - - if (partonicExists) { - int partonicPdg = mcParticles.iteratorAt(partIdx).pdgCode(); - - mcProducts.producedPartonicMothers(partonicPdg); - partonicRow = mcProducts.producedPartonicMothers.lastIndex(); - mTrackToMcPartonicMap[trackIndex] = partonicRow; - } - } - mcProducts.producedTrackLabels(mcParticleRow, mothersRow, partonicRow); + fillMcLabelGeneric(col, mcCols, track, mcParticles, mcProducts, [](auto& prod, int64_t p, int64_t m, int64_t pm) { prod.producedTrackLabels(p, m, pm); }); } template void fillMcLambdaWithLabel(T1 const& col, T2 const& mcCols, T3 const& lambda, T4 const& mcParticles, T5& mcProducts) { if (!mProduceLambdaLabels) { - return; - } - - if (!lambda.has_mcParticle()) { mcProducts.producedLambdaLabels(-1, -1, -1); return; } - - auto mcParticle = lambda.template mcParticle_as(); - auto mcCol = mcParticle.template mcCollision_as(); - - int64_t particleIndex = mcParticle.globalIndex(); - int64_t lambdaIndex = lambda.globalIndex(); - - int64_t mcParticleRow = -1; - - // MC PARTICLE - auto itP = mLambdaToMcParticleMap.find(particleIndex); - if (itP != mLambdaToMcParticleMap.end()) { - mcParticleRow = itP->second; - } else { - auto origin = this->getOrigin(col, mcCols, mcParticle); - int64_t mcColId = this->getMcColId(mcCol, mcProducts); - - mcProducts.producedMcParticles( - mcColId, - static_cast(origin), - mcParticle.pdgCode(), - mcParticle.pt() * utils::signum(mcParticle.pdgCode()), - mcParticle.eta(), - mcParticle.phi()); - - mcParticleRow = mcProducts.producedMcParticles.lastIndex(); - mLambdaToMcParticleMap[particleIndex] = mcParticleRow; - } - - // mothers (fill only if exists) - int64_t mothersRow = -1; - auto itM = mLambdaToMcMotherMap.find(lambdaIndex); - - if (itM != mLambdaToMcMotherMap.end()) { - mothersRow = itM->second; - } else { - - auto mothers = mcParticle.template mothers_as(); - bool motherExists = mcParticle.has_mothers() && !mothers.empty(); - - if (motherExists) { - int motherPdg = mothers.front().pdgCode(); // PDG code is ALWAYS valid if the mother exists - - mcProducts.producedMothers(motherPdg); - mothersRow = mcProducts.producedMothers.lastIndex(); - mLambdaToMcMotherMap[lambdaIndex] = mothersRow; - } - } - - // partonic mother (fill only if exists) - int64_t partonicRow = -1; - auto itPM = mLambdaToMcPartonicMap.find(lambdaIndex); - - if (itPM != mLambdaToMcPartonicMap.end()) { - partonicRow = itPM->second; - } else { - int partIdx = -1; - if (mcParticle.has_mothers()) { - partIdx = this->findFirstPartonicMother(mcParticle, mcParticles); - } - - bool partonicExists = (partIdx >= 0); - - if (partonicExists) { - int partonicPdg = mcParticles.iteratorAt(partIdx).pdgCode(); - - mcProducts.producedPartonicMothers(partonicPdg); - partonicRow = mcProducts.producedPartonicMothers.lastIndex(); - mLambdaToMcPartonicMap[lambdaIndex] = partonicRow; - } - } - mcProducts.producedLambdaLabels(mcParticleRow, mothersRow, partonicRow); + fillMcLabelGeneric(col, mcCols, lambda, mcParticles, mcProducts, [](auto& prod, int64_t p, int64_t m, int64_t pm) { prod.producedLambdaLabels(p, m, pm); }); } template void fillMcK0shortWithLabel(T1 const& col, T2 const& mcCols, T3 const& k0short, T4 const& mcParticles, T5& mcProducts) { if (!mProduceK0shortLabels) { - return; - } - - if (!k0short.has_mcParticle()) { mcProducts.producedK0shortLabels(-1, -1, -1); return; } - - auto mcParticle = k0short.template mcParticle_as(); - auto mcCol = mcParticle.template mcCollision_as(); - - int64_t particleIndex = mcParticle.globalIndex(); - int64_t k0shortIndex = k0short.globalIndex(); - - int64_t mcParticleRow = -1; - - // MC PARTICLE - auto itP = mK0shortToMcParticleMap.find(particleIndex); - if (itP != mK0shortToMcParticleMap.end()) { - mcParticleRow = itP->second; - } else { - auto origin = this->getOrigin(col, mcCols, mcParticle); - int64_t mcColId = this->getMcColId(mcCol, mcProducts); - - mcProducts.producedMcParticles( - mcColId, - static_cast(origin), - mcParticle.pdgCode(), - mcParticle.pt(), - mcParticle.eta(), - mcParticle.phi()); - - mcParticleRow = mcProducts.producedMcParticles.lastIndex(); - mK0shortToMcParticleMap[particleIndex] = mcParticleRow; - } - - // mothers (fill only if exists) - int64_t mothersRow = -1; - auto itM = mK0shortToMcMotherMap.find(k0shortIndex); - - if (itM != mK0shortToMcMotherMap.end()) { - mothersRow = itM->second; - } else { - - auto mothers = mcParticle.template mothers_as(); - bool motherExists = mcParticle.has_mothers() && !mothers.empty(); - - if (motherExists) { - int motherPdg = mothers.front().pdgCode(); // PDG code is ALWAYS valid if the mother exists - - mcProducts.producedMothers(motherPdg); - mothersRow = mcProducts.producedMothers.lastIndex(); - mK0shortToMcMotherMap[k0shortIndex] = mothersRow; - } - } - - // partonic mother (fill only if exists) - int64_t partonicRow = -1; - auto itPM = mK0shortToMcPartonicMap.find(k0shortIndex); - - if (itPM != mK0shortToMcPartonicMap.end()) { - partonicRow = itPM->second; - } else { - int partIdx = -1; - if (mcParticle.has_mothers()) { - partIdx = this->findFirstPartonicMother(mcParticle, mcParticles); - } - - bool partonicExists = (partIdx >= 0); - - if (partonicExists) { - int partonicPdg = mcParticles.iteratorAt(partIdx).pdgCode(); - - mcProducts.producedPartonicMothers(partonicPdg); - partonicRow = mcProducts.producedPartonicMothers.lastIndex(); - mK0shortToMcPartonicMap[k0shortIndex] = partonicRow; - } - } - mcProducts.producedK0shortLabels(mcParticleRow, mothersRow, partonicRow); + fillMcLabelGeneric(col, mcCols, k0short, mcParticles, mcProducts, [](auto& prod, int64_t p, int64_t m, int64_t pm) { prod.producedK0shortLabels(p, m, pm); }); } - template - void fillMcSigmaWithLabel(T1 const& col, T2 const& mcCols, T3 const& sigma, T4 const& sigmaDaughter, T5 const& mcParticles, T6& mcProducts) + template + void fillMcSigmaWithLabel(T1 const& col, T2 const& mcCols, T3 const& sigmaDaughter, T4 const& mcParticles, T5& mcProducts) { if (!mProduceSigmaLabels) { return; } + fillMcLabelGeneric(col, mcCols, sigmaDaughter, mcParticles, mcProducts, [](auto& prod, int64_t p, int64_t m, int64_t pm) { prod.producedSigmaLabels(p, m, pm); }, true); + } - if (!sigmaDaughter.has_mcParticle()) { - mcProducts.producedSigmaLabels(-1, -1, -1); - return; - } - - auto mcKinkDaughter = sigmaDaughter.template mcParticle_as(); - auto mcKinkDaughterMothers = mcKinkDaughter.template mothers_as(); - - if (!mcKinkDaughter.has_mothers() && !mcKinkDaughterMothers.empty()) { - mcProducts.producedSigmaLabels(-1, -1, -1); + template + void fillMcSigmaPlusWithLabel(T1 const& col, T2 const& mcCols, T3 const& sigmaPlusDaughter, T4 const& mcParticles, T5& mcProducts) + { + if (!mProduceSigmaPlusLabels) { return; } + fillMcLabelGeneric(col, mcCols, sigmaPlusDaughter, mcParticles, mcProducts, [](auto& prod, int64_t p, int64_t m, int64_t pm) { prod.producedSigmaPlusLabels(p, m, pm); }, true); + } - // we get the mc kink partilce via the mother of the kink daughter - auto mcParticle = mcKinkDaughterMothers.front(); - auto mcCol = mcParticle.template mcCollision_as(); - - int64_t particleIndex = mcParticle.globalIndex(); - int64_t sigmaIndex = sigma.globalIndex(); - - int64_t mcParticleRow = -1; - - // MC PARTICLE - auto itP = mSigmaToMcParticleMap.find(particleIndex); - if (itP != mSigmaToMcParticleMap.end()) { - mcParticleRow = itP->second; - } else { - auto origin = this->getOrigin(col, mcCols, mcParticle); - int64_t mcColId = this->getMcColId(mcCol, mcProducts); - - mcProducts.producedMcParticles( - mcColId, - static_cast(origin), - mcParticle.pdgCode(), - mcParticle.pt() * utils::signum(mcParticle.pdgCode()), - mcParticle.eta(), - mcParticle.phi()); - - mcParticleRow = mcProducts.producedMcParticles.lastIndex(); - mSigmaToMcParticleMap[particleIndex] = mcParticleRow; - } - - // mothers (fill only if exists) - int64_t mothersRow = -1; - auto itM = mSigmaToMcMotherMap.find(sigmaIndex); - - if (itM != mSigmaToMcMotherMap.end()) { - mothersRow = itM->second; - } else { - - auto mothers = mcParticle.template mothers_as(); - bool motherExists = mcParticle.has_mothers() && !mothers.empty(); - - if (motherExists) { - int motherPdg = mothers.front().pdgCode(); // PDG code is ALWAYS valid if the mother exists + bool fillAnyTable() const { return mFillAnyTable; } - mcProducts.producedMothers(motherPdg); - mothersRow = mcProducts.producedMothers.lastIndex(); - mSigmaToMcMotherMap[sigmaIndex] = mothersRow; - } - } + void reset() + { + mCollisionMap.clear(); + mMcParticleMap.clear(); + mMcMotherMap.clear(); + mMcPartonicMotherMap.clear(); + } - // partonic mother (fill only if exists) - int64_t partonicRow = -1; - auto itPM = mSigmaToMcPartonicMap.find(sigmaIndex); + private: + template + void fillMcCollision(T1 const& mcCol, T2& mcProducts) + { + mcProducts.producedMcCollisions( + mcCol.multMCNParticlesEta08(), + mcCol.centFT0M()); + mCollisionMap.emplace(mcCol.globalIndex(), mcProducts.producedMcCollisions.lastIndex()); + } - if (itPM != mSigmaToMcPartonicMap.end()) { - partonicRow = itPM->second; - } else { - int partIdx = -1; - if (mcParticle.has_mothers()) { - partIdx = this->findFirstPartonicMother(mcParticle, mcParticles); + template + modes::McOrigin getOrigin(T1 const& col, T2 const& /*mcCols*/, T3 const& mcParticle) + { + // constants + const int producedByDecay = 4; + // check if reconstructed collision has a generated collision + if (col.has_mcCollision()) { + // now check collision ids, if they do not match, then the track belongs to another collision + if (col.mcCollisionId() != mcParticle.mcCollisionId()) { + return modes::McOrigin::kFromWrongCollision; } - - bool partonicExists = (partIdx >= 0); - - if (partonicExists) { - int partonicPdg = mcParticles.iteratorAt(partIdx).pdgCode(); - - mcProducts.producedPartonicMothers(partonicPdg); - partonicRow = mcProducts.producedPartonicMothers.lastIndex(); - mSigmaToMcPartonicMap[sigmaIndex] = partonicRow; + if (mcParticle.isPhysicalPrimary()) { + return modes::McOrigin::kPhysicalPrimary; + } else if (mcParticle.has_mothers() && mcParticle.getProcess() == producedByDecay) { + return modes::McOrigin::kFromSecondaryDecay; + } else { + // not a primary and not from a decay, we label as material + return modes::McOrigin::kFromMaterial; } } - mcProducts.producedSigmaLabels(mcParticleRow, mothersRow, partonicRow); + return modes::McOrigin::kFromWrongCollision; } - template - void fillMcSigmaPlusWithLabel(T1 const& col, T2 const& mcCols, T3 const& sigmaPlus, T4 const& sigmaPlusDaughter, T5 const& mcParticles, T6& mcProducts) + template + int64_t getMcColId(T1 const& mcCol, T2& mcProducts) { - if (!mProduceSigmaPlusLabels) { - return; - } - - if (!sigmaPlusDaughter.has_mcParticle()) { - mcProducts.producedSigmaPlusLabels(-1, -1, -1); - return; + auto gid = mcCol.globalIndex(); + // Find or create + auto it = mCollisionMap.find(gid); + if (it == mCollisionMap.end()) { + this->fillMcCollision(mcCol, mcProducts); + it = mCollisionMap.find(gid); + return it->second; } + return it->second; + } - auto mcKinkDaughter = sigmaPlusDaughter.template mcParticle_as(); - auto mcKinkDaughterMothers = mcKinkDaughter.template mothers_as(); - - if (!mcKinkDaughter.has_mothers() && !mcKinkDaughterMothers.empty()) { - mcProducts.producedSigmaPlusLabels(-1, -1, -1); + template + void fillMcLabelGeneric(T1 const& col, + T2 const& mcCols, + T3 const& particle, + T4 const& mcParticles, + T5& mcProducts, + T6 writeLabels, + bool startFromMotherParticle = false) + { + if (!particle.has_mcParticle()) { + writeLabels(mcProducts, -1, -1, -1); return; } - // we get the mc kink partilce via the mother of the kink daughter - auto mcParticle = mcKinkDaughterMothers.front(); + auto mcParticle = particle.template mcParticle_as(); auto mcCol = mcParticle.template mcCollision_as(); - int64_t particleIndex = mcParticle.globalIndex(); - int64_t sigmaIndex = sigmaPlus.globalIndex(); + if (startFromMotherParticle) { + // in case of e.g. sigmas we do not reconstruct the mother but the daughter, so here we want to start from the mother particle + auto mcDaughterParticle = particle.template mcParticle_as(); + if (!mcDaughterParticle.has_mothers()) { + writeLabels(mcProducts, -1, -1, -1); + return; + } + auto mothersOfDaughter = mcDaughterParticle.template mothers_as(); + mcParticle = mothersOfDaughter.front(); + } + int64_t mcParticleIndex = mcParticle.globalIndex(); int64_t mcParticleRow = -1; - // MC PARTICLE - auto itP = mSigmaPlusToMcParticleMap.find(particleIndex); - if (itP != mSigmaPlusToMcParticleMap.end()) { + // MC particle + auto itP = mMcParticleMap.find(mcParticleIndex); + if (itP != mMcParticleMap.end()) { mcParticleRow = itP->second; } else { auto origin = this->getOrigin(col, mcCols, mcParticle); @@ -571,62 +283,51 @@ class McBuilder mcParticle.phi()); mcParticleRow = mcProducts.producedMcParticles.lastIndex(); - mSigmaPlusToMcParticleMap[particleIndex] = mcParticleRow; + mMcParticleMap[mcParticleIndex] = mcParticleRow; } - // mothers (fill only if exists) - int64_t mothersRow = -1; - auto itM = mSigmaPlusToMcMotherMap.find(sigmaIndex); - - if (itM != mSigmaPlusToMcMotherMap.end()) { - mothersRow = itM->second; - } else { - - auto mothers = mcParticle.template mothers_as(); - bool motherExists = mcParticle.has_mothers() && !mothers.empty(); - - if (motherExists) { - int motherPdg = mothers.front().pdgCode(); // PDG code is ALWAYS valid if the mother exists + // MC mother + int64_t mcMotherRow = -1; + if (mcParticle.has_mothers()) { + auto mothers = mcParticle.template mothers_as(); + auto mcMotherIndex = mothers.front().globalIndex(); - mcProducts.producedMothers(motherPdg); - mothersRow = mcProducts.producedMothers.lastIndex(); - mSigmaPlusToMcMotherMap[sigmaIndex] = mothersRow; + auto itM = mMcMotherMap.find(mcMotherIndex); + if (itM != mMcMotherMap.end()) { + mcMotherRow = itM->second; + } else { + mcProducts.producedMothers(mothers.front().pdgCode()); + mcMotherRow = mcProducts.producedMothers.lastIndex(); + mMcMotherMap[mcMotherIndex] = mcMotherRow; } } - // partonic mother (fill only if exists) - int64_t partonicRow = -1; - auto itPM = mSigmaPlusToMcPartonicMap.find(sigmaIndex); - - if (itPM != mSigmaPlusToMcPartonicMap.end()) { - partonicRow = itPM->second; - } else { - int partIdx = -1; - if (mcParticle.has_mothers()) { - partIdx = this->findFirstPartonicMother(mcParticle, mcParticles); - } - - bool partonicExists = (partIdx >= 0); - - if (partonicExists) { - int partonicPdg = mcParticles.iteratorAt(partIdx).pdgCode(); - - mcProducts.producedPartonicMothers(partonicPdg); - partonicRow = mcProducts.producedPartonicMothers.lastIndex(); - mSigmaPlusToMcPartonicMap[sigmaIndex] = partonicRow; + // Partonic mother + int64_t mcPartonicMotherRow = -1; + auto mcPartonicMotherIndex = this->findFirstPartonicMother(mcParticle, mcParticles); + if (mcPartonicMotherIndex >= 0) { + auto itPM = mMcPartonicMotherMap.find(mcPartonicMotherIndex); + if (itPM != mMcPartonicMotherMap.end()) { + mcPartonicMotherRow = itPM->second; + } else { + mcProducts.producedPartonicMothers( + mcParticles.iteratorAt(mcPartonicMotherIndex).pdgCode()); + mcPartonicMotherRow = mcProducts.producedPartonicMothers.lastIndex(); + mMcPartonicMotherMap[mcPartonicMotherIndex] = mcPartonicMotherRow; } } - mcProducts.producedSigmaPlusLabels(mcParticleRow, mothersRow, partonicRow); + + writeLabels(mcProducts, mcParticleRow, mcMotherRow, mcPartonicMotherRow); } - template - int findFirstPartonicMother(const TParticle& p, const TContainer& mcParticles) + template + int findFirstPartonicMother(const T1& mcParticle, const T2& mcParticles) { - if (!p.has_mothers()) { + if (!mcParticle.has_mothers()) { return -1; } - auto motherIds = p.mothersIds(); + auto motherIds = mcParticle.mothersIds(); // adapted these checks from MCUtils in PWGEM const int defaultMotherSize = 2; @@ -650,10 +351,10 @@ class McBuilder continue; const auto& mother = mcParticles.iteratorAt(i); - int apdg = std::abs(mother.pdgCode()); + int pdgAbs = std::abs(mother.pdgCode()); // Is it a parton? (quark or gluon) - if (apdg <= PDG_t::kTop || apdg == PDG_t::kGluon) { + if (pdgAbs <= PDG_t::kTop || pdgAbs == PDG_t::kGluon) { return i; // Found a parton → return index } @@ -667,34 +368,6 @@ class McBuilder return -1; } - void reset() - { - mCollisionMap.clear(); - - mTrackToMcParticleMap.clear(); - mTrackToMcMotherMap.clear(); - mTrackToMcPartonicMap.clear(); - - mLambdaToMcParticleMap.clear(); - mLambdaToMcMotherMap.clear(); - mLambdaToMcPartonicMap.clear(); - - mK0shortToMcParticleMap.clear(); - mK0shortToMcMotherMap.clear(); - mK0shortToMcPartonicMap.clear(); - - mSigmaToMcParticleMap.clear(); - mSigmaToMcMotherMap.clear(); - mSigmaToMcPartonicMap.clear(); - - mSigmaPlusToMcParticleMap.clear(); - mSigmaPlusToMcMotherMap.clear(); - mSigmaPlusToMcPartonicMap.clear(); - } - - bool fillAnyTable() const { return mFillAnyTable; } - - private: bool mPassThrough = false; bool mFillAnyTable = false; bool mProduceMcCollisions = false; @@ -711,25 +384,9 @@ class McBuilder std::unordered_map mCollisionMap; - std::unordered_map mTrackToMcParticleMap; - std::unordered_map mTrackToMcMotherMap; - std::unordered_map mTrackToMcPartonicMap; - - std::unordered_map mLambdaToMcParticleMap; - std::unordered_map mLambdaToMcMotherMap; - std::unordered_map mLambdaToMcPartonicMap; - - std::unordered_map mK0shortToMcParticleMap; - std::unordered_map mK0shortToMcMotherMap; - std::unordered_map mK0shortToMcPartonicMap; - - std::unordered_map mSigmaToMcParticleMap; - std::unordered_map mSigmaToMcMotherMap; - std::unordered_map mSigmaToMcPartonicMap; - - std::unordered_map mSigmaPlusToMcParticleMap; - std::unordered_map mSigmaPlusToMcMotherMap; - std::unordered_map mSigmaPlusToMcPartonicMap; + std::unordered_map mMcParticleMap; + std::unordered_map mMcMotherMap; + std::unordered_map mMcPartonicMotherMap; }; } // namespace mcbuilder diff --git a/PWGCF/Femto/Core/modes.h b/PWGCF/Femto/Core/modes.h index 2fc3dd5e0d9..bf337fab232 100644 --- a/PWGCF/Femto/Core/modes.h +++ b/PWGCF/Femto/Core/modes.h @@ -87,16 +87,37 @@ enum class Particle : o2::aod::femtodatatypes::ParticleType { }; enum class McOrigin : o2::aod::femtodatatypes::McOriginType { - kNoMcParticle, - kFromWrongCollision, - kPhysicalPrimary, - kFromSecondaryDecay, - kFromMaterial, - kMcOriginLast + kNoMcParticle = 0, // no associated mc particle normally indicated a wrongly reconstruced partilce + kFromWrongCollision = 1, // partilce originates from the wrong collision or a collision which was wrongly reconstructed (like a split vertex) + kPhysicalPrimary = 2, // primary particle + kFromSecondaryDecay = 3, // particle from secondary decay + kFromMaterial = 4, // partilce orginates from material + kMissidentified = 5, // partilce was kMissidentified (also know as fake) + kMcOriginLast = 6 // kFromFakeRecoCollision, // kFromUnkown }; +constexpr const char* mcOriginToString(McOrigin origin) +{ + switch (origin) { + case McOrigin::kNoMcParticle: + return "NoMcParticle"; + case McOrigin::kFromWrongCollision: + return "FromWrongCollision"; + case McOrigin::kPhysicalPrimary: + return "PhysicalPrimary"; + case McOrigin::kFromSecondaryDecay: + return "FromSecondaryDecay"; + case McOrigin::kFromMaterial: + return "FromMaterial"; + case McOrigin::kMissidentified: + return "Missidentified"; + default: + return "UnknownMcOrigin"; + } +} + constexpr bool hasMass(Particle p) { diff --git a/PWGCF/Femto/Core/pairBuilder.h b/PWGCF/Femto/Core/pairBuilder.h index 9fef5e04b37..41d8b8d894b 100644 --- a/PWGCF/Femto/Core/pairBuilder.h +++ b/PWGCF/Femto/Core/pairBuilder.h @@ -89,26 +89,27 @@ class PairTrackTrackBuilder mColHistManager.template init(registry, colHistSpec); mPairHistManagerSe.template init(registry, pairHistSpec, confPairBinning, confPairCuts); mPairHistManagerMe.template init(registry, pairHistSpec, confPairBinning, confPairCuts); + mPc.template init(confPairCuts); if (mSameSpecies) { - mTrackHistManager1.template init(registry, trackHistSpec1, confTrackSelection1.chargeAbs.value); + mTrackHistManager1.template init(registry, trackHistSpec1, confTrackSelection1); - mPairHistManagerSe.setMass(confTrackSelection1.pdgCode.value, confTrackSelection1.pdgCode.value); + mPairHistManagerSe.setMass(confTrackSelection1.pdgCodeAbs.value, confTrackSelection1.pdgCodeAbs.value); mPairHistManagerSe.setCharge(confTrackSelection1.chargeAbs.value, confTrackSelection1.chargeAbs.value); mCprSe.init(registry, cprHistSpec, confCpr, confTrackSelection1.chargeAbs.value, confTrackSelection1.chargeAbs.value); - mPairHistManagerMe.setMass(confTrackSelection1.pdgCode.value, confTrackSelection1.pdgCode.value); + mPairHistManagerMe.setMass(confTrackSelection1.pdgCodeAbs.value, confTrackSelection1.pdgCodeAbs.value); mPairHistManagerMe.setCharge(confTrackSelection1.chargeAbs.value, confTrackSelection1.chargeAbs.value); mCprMe.init(registry, cprHistSpec, confCpr, confTrackSelection1.chargeAbs.value, confTrackSelection1.chargeAbs.value); } else { - mTrackHistManager1.template init(registry, trackHistSpec1, confTrackSelection1.chargeAbs.value); - mTrackHistManager2.template init(registry, trackHistSpec2, confTrackSelection2.chargeAbs.value); + mTrackHistManager1.template init(registry, trackHistSpec1, confTrackSelection1); + mTrackHistManager2.template init(registry, trackHistSpec2, confTrackSelection2); - mPairHistManagerSe.setMass(confTrackSelection1.pdgCode.value, confTrackSelection2.pdgCode.value); + mPairHistManagerSe.setMass(confTrackSelection1.pdgCodeAbs.value, confTrackSelection2.pdgCodeAbs.value); mPairHistManagerSe.setCharge(confTrackSelection1.chargeAbs.value, confTrackSelection2.chargeAbs.value); mCprSe.init(registry, cprHistSpec, confCpr, confTrackSelection1.chargeAbs.value, confTrackSelection2.chargeAbs.value); - mPairHistManagerMe.setMass(confTrackSelection1.pdgCode.value, confTrackSelection2.pdgCode.value); + mPairHistManagerMe.setMass(confTrackSelection1.pdgCodeAbs.value, confTrackSelection2.pdgCodeAbs.value); mPairHistManagerMe.setCharge(confTrackSelection1.chargeAbs.value, confTrackSelection2.chargeAbs.value); mCprMe.init(registry, cprHistSpec, confCpr, confTrackSelection1.chargeAbs.value, confTrackSelection2.chargeAbs.value); } @@ -130,6 +131,7 @@ class PairTrackTrackBuilder } } + // data template void processSameEvent(T1 const& col, T2& trackTable, T3& partition1, T4& partition2, T5& cache) { @@ -153,6 +155,30 @@ class PairTrackTrackBuilder } } + // mc + template + void processSameEvent(T1 const& col, T2 const& mcCols, T3& trackTable, T4& partition1, T5& partition2, T6 const& mcParticles, T7 const& mcMothers, T8 const& mcPartonicMothers, T9& cache) + { + if (mSameSpecies) { + auto trackSlice1 = partition1->sliceByCached(o2::aod::femtobase::stored::fColId, col.globalIndex(), cache); + if (trackSlice1.size() == 0) { + return; + } + mColHistManager.template fill(col, mcCols); + mCprSe.setMagField(col.magField()); + pairprocesshelpers::processSameEvent(trackSlice1, trackTable, mcParticles, mcMothers, mcPartonicMothers, col, mcCols, mTrackHistManager1, mPairHistManagerSe, mCprSe, mPc, mRng, mMixIdenticalParticles); + } else { + auto trackSlice1 = partition1->sliceByCached(o2::aod::femtobase::stored::fColId, col.globalIndex(), cache); + auto trackSlice2 = partition2->sliceByCached(o2::aod::femtobase::stored::fColId, col.globalIndex(), cache); + if (trackSlice1.size() == 0 || trackSlice2.size() == 0) { + return; + } + mColHistManager.template fill(col, mcCols); + mCprSe.setMagField(col.magField()); + pairprocesshelpers::processSameEvent(trackSlice1, trackSlice2, trackTable, mcParticles, mcMothers, mcPartonicMothers, col, mcCols, mTrackHistManager1, mTrackHistManager2, mPairHistManagerSe, mCprSe, mPc); + } + } + template void processMixedEvent(T1 const& cols, T2& trackTable, T3& partition1, T4& partition2, T5& cache, T6& binsVtxMult, T7& binsVtxCent, T8& binsVtxMultCent) { @@ -188,6 +214,40 @@ class PairTrackTrackBuilder } } + template + void processMixedEvent(T1 const& cols, T2 const& mcCols, T3& trackTable, T4& partition1, T5& partition2, T6 const& mcParticles, T7& cache, T8& binsVtxMult, T9& binsVtxCent, T10& binsVtxMultCent) + { + if (mSameSpecies) { + switch (mMixingPolicy) { + case static_cast(pairhistmanager::kVtxMult): + pairprocesshelpers::processMixedEvent(cols, mcCols, partition1, trackTable, mcParticles, cache, binsVtxMult, mMixingDepth, mPairHistManagerMe, mCprMe, mPc); + break; + case static_cast(pairhistmanager::kVtxCent): + pairprocesshelpers::processMixedEvent(cols, mcCols, partition1, trackTable, mcParticles, cache, binsVtxCent, mMixingDepth, mPairHistManagerMe, mCprMe, mPc); + break; + case static_cast(pairhistmanager::kVtxMultCent): + pairprocesshelpers::processMixedEvent(cols, mcCols, partition1, trackTable, mcParticles, cache, binsVtxMultCent, mMixingDepth, mPairHistManagerMe, mCprMe, mPc); + break; + default: + LOG(fatal) << "Invalid binning policiy specifed. Breaking..."; + } + } else { + switch (mMixingPolicy) { + case static_cast(pairhistmanager::kVtxMult): + pairprocesshelpers::processMixedEvent(cols, mcCols, partition1, partition2, trackTable, mcParticles, cache, binsVtxMult, mMixingDepth, mPairHistManagerMe, mCprMe, mPc); + break; + case static_cast(pairhistmanager::kVtxCent): + pairprocesshelpers::processMixedEvent(cols, mcCols, partition1, partition2, trackTable, mcParticles, cache, binsVtxCent, mMixingDepth, mPairHistManagerMe, mCprMe, mPc); + break; + case static_cast(pairhistmanager::kVtxMultCent): + pairprocesshelpers::processMixedEvent(cols, mcCols, partition1, partition2, trackTable, mcParticles, cache, binsVtxMultCent, mMixingDepth, mPairHistManagerMe, mCprMe, mPc); + break; + default: + LOG(fatal) << "Invalid binning policiy specifed. Breaking..."; + } + } + } + private: colhistmanager::CollisionHistManager mColHistManager; trackhistmanager::TrackHistManager mTrackHistManager1; @@ -266,24 +326,24 @@ class PairV0V0Builder mPairHistManagerMe.template init(registry, pairHistSpec, confPairBinning, confPairCuts); if (mSameSpecies) { - mV0HistManager1.template init(registry, V0HistSpec1, PosDauHistSpec, NegDauHistSpec); + mV0HistManager1.template init(registry, V0HistSpec1, confV0Selection1, PosDauHistSpec, NegDauHistSpec); - mPairHistManagerSe.setMass(confV0Selection1.pdgCode.value, confV0Selection1.pdgCode.value); + mPairHistManagerSe.setMass(confV0Selection1.pdgCodeAbs.value, confV0Selection1.pdgCodeAbs.value); mPairHistManagerSe.setCharge(1, 1); mCprSe.init(registry, cprHistSpecPos, cprHistSpecNeg, confCprPos, confCprPos); - mPairHistManagerMe.setMass(confV0Selection1.pdgCode.value, confV0Selection1.pdgCode.value); + mPairHistManagerMe.setMass(confV0Selection1.pdgCodeAbs.value, confV0Selection1.pdgCodeAbs.value); mPairHistManagerMe.setCharge(1, 1); mCprMe.init(registry, cprHistSpecPos, cprHistSpecNeg, confCprPos, confCprNeg); } else { - mV0HistManager1.template init(registry, V0HistSpec1, PosDauHistSpec, NegDauHistSpec); - mV0HistManager2.template init(registry, V0HistSpec2, PosDauHistSpec, NegDauHistSpec); + mV0HistManager1.template init(registry, V0HistSpec1, confV0Selection1, PosDauHistSpec, NegDauHistSpec); + mV0HistManager2.template init(registry, V0HistSpec2, confV0Selection2, PosDauHistSpec, NegDauHistSpec); - mPairHistManagerSe.setMass(confV0Selection1.pdgCode.value, confV0Selection2.pdgCode.value); + mPairHistManagerSe.setMass(confV0Selection1.pdgCodeAbs.value, confV0Selection2.pdgCodeAbs.value); mPairHistManagerSe.setCharge(1, 1); mCprSe.init(registry, cprHistSpecPos, cprHistSpecNeg, confCprPos, confCprNeg); - mPairHistManagerMe.setMass(confV0Selection1.pdgCode.value, confV0Selection2.pdgCode.value); + mPairHistManagerMe.setMass(confV0Selection1.pdgCodeAbs.value, confV0Selection2.pdgCodeAbs.value); mPairHistManagerMe.setCharge(1, 1); mCprMe.init(registry, cprHistSpecPos, cprHistSpecNeg, confCprPos, confCprNeg); } @@ -425,18 +485,19 @@ class PairTrackV0Builder { mColHistManager.template init(registry, colHistSpec); - mTrackHistManager.template init(registry, trackHistSpec, confTrackSelection.chargeAbs.value); - mV0HistManager.template init(registry, v0HistSpec, posDauHistSpec, negDauHistSpec); + mTrackHistManager.template init(registry, trackHistSpec, confTrackSelection); + mV0HistManager.template init(registry, v0HistSpec, confV0Selection, posDauHistSpec, negDauHistSpec); mPairHistManagerSe.template init(registry, pairHistSpec, confPairBinning, confPairCuts); - mPairHistManagerSe.setMass(confTrackSelection.pdgCode.value, confV0Selection.pdgCode.value); + mPairHistManagerSe.setMass(confTrackSelection.pdgCodeAbs.value, confV0Selection.pdgCodeAbs.value); mPairHistManagerSe.setCharge(confTrackSelection.chargeAbs.value, 1); mCprSe.init(registry, cprHistSpec, confCpr, confTrackSelection.chargeAbs.value); mPairHistManagerMe.template init(registry, pairHistSpec, confPairBinning, confPairCuts); - mPairHistManagerMe.setMass(confTrackSelection.pdgCode.value, confV0Selection.pdgCode.value); + mPairHistManagerMe.setMass(confTrackSelection.pdgCodeAbs.value, confV0Selection.pdgCodeAbs.value); mPairHistManagerMe.setCharge(confTrackSelection.chargeAbs.value, 1); mCprMe.init(registry, cprHistSpec, confCpr, confTrackSelection.chargeAbs.value); + mPc.template init(confPairCuts); // setup mixing mMixingPolicy = static_cast(confMixing.policy.value); @@ -456,6 +517,19 @@ class PairTrackV0Builder pairprocesshelpers::processSameEvent(trackSlice, v0Slice, trackTable, col, mTrackHistManager, mV0HistManager, mPairHistManagerSe, mCprSe, mPc); } + template + void processSameEvent(T1 const& col, T2 const& mcCols, T3 const& trackTable, T4& trackPartition, T5 const& /*v0table*/, T6& v0Partition, T7 const& mcParticles, T8 const& mcMothers, T9 const& mcPartonicMothers, T10& cache) + { + auto trackSlice = trackPartition->sliceByCached(o2::aod::femtobase::stored::fColId, col.globalIndex(), cache); + auto v0Slice = v0Partition->sliceByCached(o2::aod::femtobase::stored::fColId, col.globalIndex(), cache); + if (trackSlice.size() == 0 || v0Slice.size() == 0) { + return; + } + mColHistManager.template fill(col, mcCols); + mCprSe.setMagField(col.magField()); + pairprocesshelpers::processSameEvent(trackSlice, v0Slice, trackTable, mcParticles, mcMothers, mcPartonicMothers, col, mcCols, mTrackHistManager, mV0HistManager, mPairHistManagerSe, mCprSe, mPc); + } + template void processMixedEvent(T1 const& cols, T2& trackTable, T3& trackPartition, T4& v0Partition, T5& cache, T6& binsVtxMult, T7& binsVtxCent, T8& binsVtxMultCent) { @@ -474,6 +548,24 @@ class PairTrackV0Builder } } + template + void processMixedEvent(T1 const& cols, T2 const& mcCols, T3& trackTable, T4& trackPartition, T5& v0Partition, T6 const& mcParticles, T7& cache, T8& binsVtxMult, T9& binsVtxCent, T10& binsVtxMultCent) + { + switch (mMixingPolicy) { + case static_cast(pairhistmanager::kVtxMult): + pairprocesshelpers::processMixedEvent(cols, mcCols, trackPartition, v0Partition, trackTable, mcParticles, cache, binsVtxMult, mMixingDepth, mPairHistManagerMe, mCprMe, mPc); + break; + case static_cast(pairhistmanager::kVtxCent): + pairprocesshelpers::processMixedEvent(cols, mcCols, trackPartition, v0Partition, trackTable, mcParticles, cache, binsVtxCent, mMixingDepth, mPairHistManagerMe, mCprMe, mPc); + break; + case static_cast(pairhistmanager::kVtxMultCent): + pairprocesshelpers::processMixedEvent(cols, mcCols, trackPartition, v0Partition, trackTable, mcParticles, cache, binsVtxMultCent, mMixingDepth, mPairHistManagerMe, mCprMe, mPc); + break; + default: + LOG(fatal) << "Invalid binning policiy specifed. Breaking..."; + } + } + private: colhistmanager::CollisionHistManager mColHistManager; trackhistmanager::TrackHistManager mTrackHistManager; @@ -533,16 +625,16 @@ class PairTrackTwoTrackResonanceBuilder { mColHistManager.template init(registry, colHistSpec); - mTrackHistManager.template init(registry, trackHistSpec, confTrackSelection.chargeAbs.value); - mResonanceHistManager.template init(registry, resonanceHistSpec, posDauHistSpec, negDauHistSpec); + mTrackHistManager.template init(registry, trackHistSpec, confTrackSelection); + mResonanceHistManager.template init(registry, resonanceHistSpec, confResonanceSelection, posDauHistSpec, negDauHistSpec); mPairHistManagerSe.template init(registry, pairHistSpec, confPairBinning, confPairCuts); - mPairHistManagerSe.setMass(confTrackSelection.pdgCode.value, confResonanceSelection.pdgCode.value); + mPairHistManagerSe.setMass(confTrackSelection.pdgCodeAbs.value, confResonanceSelection.pdgCodeAbs.value); mPairHistManagerSe.setCharge(confTrackSelection.chargeAbs.value, 1); mCprSe.init(registry, cprHistSpec, confCpr, confTrackSelection.chargeAbs.value); mPairHistManagerMe.template init(registry, pairHistSpec, confPairBinning, confPairCuts); - mPairHistManagerMe.setMass(confTrackSelection.pdgCode.value, confResonanceSelection.pdgCode.value); + mPairHistManagerMe.setMass(confTrackSelection.pdgCodeAbs.value, confResonanceSelection.pdgCodeAbs.value); mPairHistManagerMe.setCharge(confTrackSelection.chargeAbs.value, 1); mCprMe.init(registry, cprHistSpec, confCpr, confTrackSelection.chargeAbs.value); @@ -638,18 +730,19 @@ class PairTrackKinkBuilder { mColHistManager.template init(registry, colHistSpec); - mTrackHistManager.template init(registry, trackHistSpec, confTrackSelection.chargeAbs.value); - mKinkHistManager.template init(registry, kinkHistSpec, chaDauHistSpec); + mTrackHistManager.template init(registry, trackHistSpec, confTrackSelection); + mKinkHistManager.template init(registry, kinkHistSpec, confKinkSelection, chaDauHistSpec); mPairHistManagerSe.template init(registry, pairHistSpec, confPairBinning, confPairCuts); - mPairHistManagerSe.setMass(confTrackSelection.pdgCode.value, confKinkSelection.pdgCode.value); + mPairHistManagerSe.setMass(confTrackSelection.pdgCodeAbs.value, confKinkSelection.pdgCodeAbs.value); mPairHistManagerSe.setCharge(confTrackSelection.chargeAbs.value, 1); // abs charge of kink daughter is always 1 mCprSe.init(registry, cprHistSpec, confCpr, confTrackSelection.chargeAbs.value); mPairHistManagerMe.template init(registry, pairHistSpec, confPairBinning, confPairCuts); - mPairHistManagerMe.setMass(confTrackSelection.pdgCode.value, confKinkSelection.pdgCode.value); + mPairHistManagerMe.setMass(confTrackSelection.pdgCodeAbs.value, confKinkSelection.pdgCodeAbs.value); mPairHistManagerMe.setCharge(confTrackSelection.chargeAbs.value, 1); // abs charge of kink daughter is always 1 mCprMe.init(registry, cprHistSpec, confCpr, confTrackSelection.chargeAbs.value); + mPc.template init(confPairCuts); // setup mixing mMixingPolicy = static_cast(confMixing.policy.value); @@ -669,6 +762,19 @@ class PairTrackKinkBuilder pairprocesshelpers::processSameEvent(trackSlice, kinkSlice, trackTable, col, mTrackHistManager, mKinkHistManager, mPairHistManagerSe, mCprSe, mPc); } + template + void processSameEvent(T1 const& col, T2 const& mcCols, T3& trackTable, T4& trackPartition, T5& /*kinktable*/, T6& kinkPartition, T7 const& mcParticles, T8 const& mcMothers, T9 const& mcPartonicMothers, T10& cache) + { + auto trackSlice = trackPartition->sliceByCached(o2::aod::femtobase::stored::fColId, col.globalIndex(), cache); + auto kinkSlice = kinkPartition->sliceByCached(o2::aod::femtobase::stored::fColId, col.globalIndex(), cache); + if (trackSlice.size() == 0 || kinkSlice.size() == 0) { + return; + } + mColHistManager.fill(col, mcCols); + mCprSe.setMagField(col.magField()); + pairprocesshelpers::processSameEvent(trackSlice, kinkSlice, trackTable, mcParticles, mcMothers, mcPartonicMothers, col, mcCols, mTrackHistManager, mKinkHistManager, mPairHistManagerSe, mCprSe, mPc); + } + template void processMixedEvent(T1 const& cols, T2& trackTable, T3& trackPartition, T4& kinkPartition, T5& cache, T6& binsVtxMult, T7& binsVtxCent, T8& binsVtxMultCent) { @@ -687,6 +793,24 @@ class PairTrackKinkBuilder } } + template + void processMixedEvent(T1 const& cols, T2 const& mcCols, T3& trackTable, T4& trackPartition, T5& kinkPartition, T6 const& mcParticles, T7& cache, T8& binsVtxMult, T9& binsVtxCent, T10& binsVtxMultCent) + { + switch (mMixingPolicy) { + case static_cast(pairhistmanager::kVtxMult): + pairprocesshelpers::processMixedEvent(cols, mcCols, trackPartition, kinkPartition, trackTable, mcParticles, cache, binsVtxMult, mMixingDepth, mPairHistManagerMe, mCprMe, mPc); + break; + case static_cast(pairhistmanager::kVtxCent): + pairprocesshelpers::processMixedEvent(cols, mcCols, trackPartition, kinkPartition, trackTable, mcParticles, cache, binsVtxCent, mMixingDepth, mPairHistManagerMe, mCprMe, mPc); + break; + case static_cast(pairhistmanager::kVtxMultCent): + pairprocesshelpers::processMixedEvent(cols, mcCols, trackPartition, kinkPartition, trackTable, mcParticles, cache, binsVtxMultCent, mMixingDepth, mPairHistManagerMe, mCprMe, mPc); + break; + default: + LOG(fatal) << "Invalid binning policiy specifed. Breaking..."; + } + } + private: colhistmanager::CollisionHistManager mColHistManager; trackhistmanager::TrackHistManager mTrackHistManager; @@ -755,16 +879,16 @@ class PairTrackCascadeBuilder { mColHistManager.template init(registry, colHistSpec); - mTrackHistManager.template init(registry, trackHistSpec, confTrackSelection.chargeAbs.value); - mCascadeHistManager.template init(registry, cascadeHistSpec, bachelorHistSpec, posDauHistSpec, negDauHistSpec); + mTrackHistManager.template init(registry, trackHistSpec, confTrackSelection); + mCascadeHistManager.template init(registry, cascadeHistSpec, confCascadeSelection, bachelorHistSpec, posDauHistSpec, negDauHistSpec); mPairHistManagerSe.template init(registry, pairHistSpec, confPairBinning, confPairCuts); - mPairHistManagerSe.setMass(confTrackSelection.pdgCode.value, confCascadeSelection.pdgCode.value); + mPairHistManagerSe.setMass(confTrackSelection.pdgCodeAbs.value, confCascadeSelection.pdgCodeAbs.value); mPairHistManagerSe.setCharge(confTrackSelection.chargeAbs.value, 1); mCprSe.init(registry, cprHistSpecBachelor, cprHistSpecV0Daughter, confCprBachelor, confCprV0Daughter, confTrackSelection.chargeAbs.value); mPairHistManagerMe.template init(registry, pairHistSpec, confPairBinning, confPairCuts); - mPairHistManagerMe.setMass(confTrackSelection.pdgCode.value, confCascadeSelection.pdgCode.value); + mPairHistManagerMe.setMass(confTrackSelection.pdgCodeAbs.value, confCascadeSelection.pdgCodeAbs.value); mPairHistManagerMe.setCharge(confTrackSelection.chargeAbs.value, 1); mCprMe.init(registry, cprHistSpecBachelor, cprHistSpecV0Daughter, confCprBachelor, confCprV0Daughter, confTrackSelection.chargeAbs.value); diff --git a/PWGCF/Femto/Core/pairCleaner.h b/PWGCF/Femto/Core/pairCleaner.h index ac2d3e9e07b..b1076efa0b0 100644 --- a/PWGCF/Femto/Core/pairCleaner.h +++ b/PWGCF/Femto/Core/pairCleaner.h @@ -16,6 +16,10 @@ #ifndef PWGCF_FEMTO_CORE_PAIRCLEANER_H_ #define PWGCF_FEMTO_CORE_PAIRCLEANER_H_ +#include "PWGCF/Femto/Core/modes.h" + +#include "fairlogger/Logger.h" + namespace o2::analysis::femto { namespace paircleaner @@ -27,23 +31,84 @@ class BasePairCleaner BasePairCleaner() = default; virtual ~BasePairCleaner() = default; + template + void init(T const& PairCuts) + { + if constexpr (modes::isFlagSet(mode, modes::Mode::kMc)) { + mMixPairsWithCommonAncestor = PairCuts.mixOnlyCommonAncestor.value; + mMixPairsWithNonCommonAncestor = PairCuts.mixOnlyNonCommonAncestor.value; + if (mMixPairsWithCommonAncestor && mMixPairsWithNonCommonAncestor) { + LOG(fatal) << "Both mixing with common and non-common ancestor is activated. Breaking..."; + } + } + } + protected: template bool isCleanTrackPair(T1 const& track1, T2 const& track2) const { return track1.globalIndex() != track2.globalIndex(); }; + + template + bool pairHasCommonAncestor(T1 const& particle1, T2 const& particle2, T3 const& /*partonicMothers*/) const + { + // if one of the two particles has no associated partonic mother, we cannot now if they have a common anchestor, so we break out with false + if (!particle1.has_fMcPartMoth() || !particle2.has_fMcPartMoth()) { + return false; + } + // get mc particles + auto partonicMother1 = particle1.template fMcPartMoth_as(); + auto partonicMother2 = particle2.template fMcPartMoth_as(); + // get partonic mothers + return partonicMother1.globalIndex() == partonicMother2.globalIndex(); + }; + + template + bool pairHasNonCommonAncestor(T1 const& particle1, T2 const& particle2, T3 const& /*partonicMothers*/) const + { + // if one of the two particles has no associated partonic mother, we cannot now if they have a non-common anchestor, so we break out with false + if (!particle1.has_fMcPartMoth() || !particle2.has_fMcPartMoth()) { + return false; + } + // get mc particles + auto partonicMother1 = particle1.template fMcPartMoth_as(); + auto partonicMother2 = particle2.template fMcPartMoth_as(); + // get partonic mothers + return partonicMother1.globalIndex() != partonicMother2.globalIndex(); + }; + + bool mMixPairsWithCommonAncestor = false; + bool mMixPairsWithNonCommonAncestor = false; }; class TrackTrackPairCleaner : public BasePairCleaner { public: TrackTrackPairCleaner() = default; + template bool isCleanPair(T1 const& track1, T2 const& track2, T3 const& /*trackTable*/) const { return this->isCleanTrackPair(track1, track2); } + + template + bool isCleanPair(T1 const& track1, T2 const& track2, T3 const& /*trackTable*/, T4 const& partonicMothers) const + { + if (!this->isCleanTrackPair(track1, track2)) { + return false; + } + // pair is clean + // no check if we require common or non-common ancestry + if (mMixPairsWithCommonAncestor) { + return this->pairHasCommonAncestor(track1, track2, partonicMothers); + } + if (mMixPairsWithNonCommonAncestor) { + return this->pairHasNonCommonAncestor(track1, track2, partonicMothers); + } + return true; + } }; class V0V0PairCleaner : public BasePairCleaner @@ -72,6 +137,23 @@ class TrackV0PairCleaner : public BasePairCleaner // also works for particles de auto negDaughter = trackTable.rawIteratorAt(v0.negDauId() - trackTable.offset()); return (this->isCleanTrackPair(posDaughter, track) && this->isCleanTrackPair(negDaughter, track)); } + + template + bool isCleanPair(T1 const& track1, T2 const& v0, T3 const& trackTable, T4 const& partonicMothers) const + { + if (!this->isCleanPair(track1, v0, trackTable)) { + return false; + } + // pair is clean + // now check if we require common or non-common ancestry + if (mMixPairsWithCommonAncestor) { + return this->pairHasCommonAncestor(track1, v0, partonicMothers); + } + if (mMixPairsWithNonCommonAncestor) { + return this->pairHasNonCommonAncestor(track1, v0, partonicMothers); + } + return true; + } }; class TrackKinkPairCleaner : public BasePairCleaner @@ -84,6 +166,23 @@ class TrackKinkPairCleaner : public BasePairCleaner auto chaDaughter = trackTable.rawIteratorAt(kink.chaDauId() - trackTable.offset()); return this->isCleanTrackPair(chaDaughter, track); } + + template + bool isCleanPair(T1 const& track1, T2 const& kink, T3 const& trackTable, T4 const& partonicMothers) const + { + if (!this->isCleanPair(track1, kink, trackTable)) { + return false; + } + // pair is clean + // now check if we require common or non-common ancestry + if (mMixPairsWithCommonAncestor) { + return this->pairHasCommonAncestor(track1, kink, partonicMothers); + } + if (mMixPairsWithNonCommonAncestor) { + return this->pairHasNonCommonAncestor(track1, kink, partonicMothers); + } + return true; + } }; class TrackCascadePairCleaner : public BasePairCleaner diff --git a/PWGCF/Femto/Core/pairHistManager.h b/PWGCF/Femto/Core/pairHistManager.h index 370b09bd3bd..21bed35be22 100644 --- a/PWGCF/Femto/Core/pairHistManager.h +++ b/PWGCF/Femto/Core/pairHistManager.h @@ -70,6 +70,12 @@ enum PairHist { kKstarVsMass1VsMult, kKstarVsMass2VsMult, kKstarVsMass1VsMass2VsMult, + // mc + kTrueKstarVsKstar, + kTrueKtVsKt, + kTrueMtVsMt, + kTrueMultVsMult, + kTrueCentVsCent, kPairHistogramLast }; @@ -125,6 +131,8 @@ struct ConfPairCuts : o2::framework::ConfigurableGroup { o2::framework::Configurable ktMin{"ktMin", -1, "Minimal kt (set to -1 to deactivate)"}; o2::framework::Configurable mtMax{"mtMax", -1, "Maximal mt (set to -1 to deactivate)"}; o2::framework::Configurable mtMin{"mtMin", -1, "Minimal mt (set to -1 to deactivate)"}; + o2::framework::Configurable mixOnlyCommonAncestor{"mixOnlyCommonAncestor", false, "Require pair to have common anchestor (in the same event)"}; + o2::framework::Configurable mixOnlyNonCommonAncestor{"mixOnlyNonCommonAncestor", false, "Require pair to have non-common anchestor (in the same event)"}; }; // the enum gives the correct index in the array @@ -161,45 +169,65 @@ constexpr std::array, kPairHistogramLast> {kKstarVsMass1VsMult, o2::framework::kTHnSparseF, "hKstarVsMass1VsMult", "k* vs m_{1} vs multiplicity; k* (GeV/#it{c}); m_{1} (GeV/#it{c}^{2}); Multiplicity"}, {kKstarVsMass2VsMult, o2::framework::kTHnSparseF, "hKstarVsMass2VsMult", "k* vs m_{2} vs multiplicity; k* (GeV/#it{c}); m_{2} (GeV/#it{c}^{2}); Multiplicity"}, {kKstarVsMass1VsMass2VsMult, o2::framework::kTHnSparseF, "hKstarVsMass1VsMass2VsMult", "k* vs m_{1} vs m_{2} vs multiplicity; k* (GeV/#it{c}); m_{1} (GeV/#it{c}^{2}); m_{2} (GeV/#it{c}^{2}); Multiplicity"}, - + {kTrueKstarVsKstar, o2::framework::kTH2F, "hTrueKstarVsKstar", "k*_{True} vs k*; k*_{True} (GeV/#it{c}); k* (GeV/#it{c})"}, + {kTrueKtVsKt, o2::framework::kTH2F, "hTrueKtVsKt", "k_{T,True} vs k_{T}; k_{T,True} (GeV/#it{c}); k_{T} (GeV/#it{c})"}, + {kTrueMtVsMt, o2::framework::kTH2F, "hTrueMtVsMt", "m_{T,True} vs m_{T}; m_{T,True} (GeV/#it{c}^{2}); m_{T,True} (GeV/#it{c}^{2})"}, + {kTrueMultVsMult, o2::framework::kTH2F, "hTrueMultVsMult", "Multiplicity_{True} vs Multiplicity; Multiplicity_{True} ; Multiplicity"}, + {kTrueCentVsCent, o2::framework::kTH2F, "hTrueCentVsCent", "Centrality_{True} vs Centrality; Centrality_{True} (%); Centrality (%)"}, }}; +#define PAIR_HIST_ANALYSIS_MAP(conf) \ + {kKstar, {conf.kstar}}, \ + {kKt, {conf.kt}}, \ + {kMt, {conf.mt}}, \ + {kPt1VsPt2, {conf.pt1, conf.pt2}}, \ + {kPt1VsKstar, {conf.pt1, conf.kstar}}, \ + {kPt2VsKstar, {conf.pt2, conf.kstar}}, \ + {kPt1VsKt, {conf.pt1, conf.kt}}, \ + {kPt2VsKt, {conf.pt2, conf.kt}}, \ + {kPt1VsMt, {conf.pt1, conf.mt}}, \ + {kPt2VsMt, {conf.pt2, conf.mt}}, \ + {kKstarVsKt, {conf.kstar, conf.kt}}, \ + {kKstarVsMt, {conf.kstar, conf.mt}}, \ + {kKstarVsMult, {conf.kstar, conf.multiplicity}}, \ + {kKstarVsCent, {conf.kstar, conf.centrality}}, \ + {kKstarVsMass1, {conf.kstar, conf.mass1}}, \ + {kKstarVsMass2, {conf.kstar, conf.mass2}}, \ + {kMass1VsMass2, {conf.mass1, conf.mass2}}, \ + {kKstarVsMtVsMult, {conf.kstar, conf.mt, conf.multiplicity}}, \ + {kKstarVsMtVsMultVsCent, {conf.kstar, conf.mt, conf.multiplicity, conf.centrality}}, \ + {kKstarVsMtVsPt1VsPt2VsMult, {conf.kstar, conf.mt, conf.pt1, conf.pt2, conf.multiplicity}}, \ + {kKstarVsMtVsPt1VsPt2VsMultVsCent, {conf.kstar, conf.mt, conf.pt1, conf.pt2, conf.multiplicity, conf.centrality}}, \ + {kKstarVsMass1VsMass2, {conf.kstar, conf.mass1, conf.mass2}}, \ + {kKstarVsMass1VsMult, {conf.kstar, conf.mass1, conf.multiplicity}}, \ + {kKstarVsMass2VsMult, {conf.kstar, conf.mass2, conf.multiplicity}}, \ + {kKstarVsMass1VsMass2VsMult, {conf.kstar, conf.mass1, conf.mass2, conf.multiplicity}}, + +#define PAIR_HIST_MC_MAP(conf) \ + {kTrueKstarVsKstar, {conf.kstar, conf.kstar}}, \ + {kTrueKtVsKt, {conf.kt, conf.kt}}, \ + {kTrueMtVsMt, {conf.mt, conf.mt}}, \ + {kTrueMultVsMult, {conf.multiplicity, conf.multiplicity}}, \ + {kTrueCentVsCent, {conf.centrality, conf.centrality}}, + template auto makePairHistSpecMap(const T& confPairBinning) { return std::map>{ - {kKstar, {confPairBinning.kstar}}, - {kKt, {confPairBinning.kt}}, - {kMt, {confPairBinning.mt}}, - {kPt1VsPt2, {confPairBinning.pt1, confPairBinning.pt2}}, - {kPt1VsKstar, {confPairBinning.pt1, confPairBinning.kstar}}, - {kPt2VsKstar, {confPairBinning.pt2, confPairBinning.kstar}}, - {kPt1VsKt, {confPairBinning.pt1, confPairBinning.kt}}, - {kPt2VsKt, {confPairBinning.pt2, confPairBinning.kt}}, - {kPt1VsMt, {confPairBinning.pt1, confPairBinning.mt}}, - {kPt2VsMt, {confPairBinning.pt2, confPairBinning.mt}}, - {kKstarVsKt, {confPairBinning.kstar, confPairBinning.kt}}, - {kKstarVsMt, {confPairBinning.kstar, confPairBinning.mt}}, - {kKstarVsMult, {confPairBinning.kstar, confPairBinning.multiplicity}}, - {kKstarVsCent, {confPairBinning.kstar, confPairBinning.centrality}}, - - {kKstarVsMass1, {confPairBinning.kstar, confPairBinning.mass1}}, - {kKstarVsMass2, {confPairBinning.kstar, confPairBinning.mass2}}, - {kMass1VsMass2, {confPairBinning.mass1, confPairBinning.mass2}}, - - {kKstarVsMtVsMult, {confPairBinning.kstar, confPairBinning.mt, confPairBinning.multiplicity}}, - {kKstarVsMtVsMultVsCent, {confPairBinning.kstar, confPairBinning.mt, confPairBinning.multiplicity, confPairBinning.centrality}}, - {kKstarVsMtVsPt1VsPt2VsMult, {confPairBinning.kstar, confPairBinning.mt, confPairBinning.pt1, confPairBinning.pt2, confPairBinning.multiplicity}}, - {kKstarVsMtVsPt1VsPt2VsMultVsCent, {confPairBinning.kstar, confPairBinning.mt, confPairBinning.pt1, confPairBinning.pt2, confPairBinning.multiplicity, confPairBinning.centrality}}, - - {kKstarVsMass1VsMass2, {confPairBinning.kstar, confPairBinning.mass1, confPairBinning.mass2}}, - {kKstarVsMass1VsMult, {confPairBinning.kstar, confPairBinning.mass1, confPairBinning.multiplicity}}, - {kKstarVsMass2VsMult, {confPairBinning.kstar, confPairBinning.mass2, confPairBinning.multiplicity}}, - {kKstarVsMass1VsMass2VsMult, {confPairBinning.kstar, confPairBinning.mass1, confPairBinning.mass2, confPairBinning.multiplicity}}, - - }; + PAIR_HIST_ANALYSIS_MAP(confPairBinning)}; +}; + +template +auto makePairMcHistSpecMap(const T& confPairBinning) +{ + return std::map>{ + PAIR_HIST_ANALYSIS_MAP(confPairBinning) + PAIR_HIST_MC_MAP(confPairBinning)}; }; +#undef PAIR_HIST_ANALYSIS_MAP +#undef PAIR_HIST_MC_MAP + constexpr char PrefixTrackTrackSe[] = "TrackTrack/SE/"; constexpr char PrefixTrackTrackMe[] = "TrackTrack/ME/"; @@ -220,10 +248,8 @@ constexpr char PrefixTrackKinkMe[] = "TrackKink/ME/"; constexpr std::string_view AnalysisDir = "Analysis/"; constexpr std::string_view QaDir = "QA/"; +constexpr std::string_view McDir = "MC/"; -/// \class FemtoDreamEventHisto -/// \brief Class for histogramming event properties -// template template @@ -269,9 +295,9 @@ class PairHistManager initAnalysis(Specs); } - // if constexpr (isFlagSet(mode, modes::Mode::kQA)) { - // std::string qaDir = std::string(prefix) + std::string(QaDir); - // } + if constexpr (isFlagSet(mode, modes::Mode::kMc)) { + initMc(Specs); + } } void setMass(int PdgParticle1, int PdgParticle2) @@ -284,7 +310,6 @@ class PairHistManager void setCharge(int chargeAbsParticle1, int chargeAbsParticle2) { // the pt stored is actually as pt/z for tracks, so in case of particles with z > 1, we have to rescale the pt (this is so far only for He3 the case) - // similarly, for neutral particles, no reason to rescale so we just set absolute charge to 1 mAbsCharge1 = std::abs(chargeAbsParticle1); mAbsCharge2 = std::abs(chargeAbsParticle2); } @@ -299,10 +324,10 @@ class PairHistManager auto partSum = mParticle1 + mParticle2; // set kT - mKt = partSum.Pt() / 2.f; + mKt = 0.5f * partSum.Pt(); // set mT - computeMt(partSum); + mMt = computeMt(partSum); // Boost particle to the pair rest frame (Prf) and calculate k* (would be equivalent using particle 2) // make a copy of particle 1 @@ -324,17 +349,79 @@ class PairHistManager template void setPair(const T1& particle1, const T2& particle2, const T3& col) { + setPair(particle1, particle2); mMult = col.mult(); mCent = col.cent(); - setPair(particle1, particle2); } template void setPair(const T1& particle1, const T2& particle2, const T3& col1, const T4& col2) { + setPair(particle1, particle2); mMult = 0.5f * (col1.mult() + col2.mult()); // if mixing with multiplicity, should be in the same mixing bin mCent = 0.5f * (col1.cent() + col2.cent()); // if mixing with centrality, should be in the same mixing bin - setPair(particle1, particle2); + } + + template + void setPairMc(const T1& particle1, const T2& particle2, const T3& /*mcParticles*/) + { + if (!particle1.has_fMcParticle() || !particle2.has_fMcParticle()) { + mHasMcPair = false; + return; + } + mHasMcPair = true; + + auto mcParticle1 = particle1.template fMcParticle_as(); + auto mcParticle2 = particle2.template fMcParticle_as(); + + mParticle1 = ROOT::Math::PtEtaPhiMVector{mAbsCharge1 * mcParticle1.pt(), mcParticle1.eta(), mcParticle1.phi(), mPdgMass1}; + mParticle2 = ROOT::Math::PtEtaPhiMVector{mAbsCharge2 * mcParticle2.pt(), mcParticle2.eta(), mcParticle2.phi(), mPdgMass2}; + auto partSum = mParticle1 + mParticle2; + + // set kT + mTrueKt = partSum.Pt() / 2.f; + + // set mT + mTrueMt = computeMt(partSum); + + // Boost particle to the pair rest frame (Prf) and calculate k* (would be equivalent using particle 2) + // make a copy of particle 1 + auto particle1Prf = ROOT::Math::PtEtaPhiMVector(mParticle1); + // get lorentz boost into pair rest frame + ROOT::Math::Boost boostPrf(partSum.BoostToCM()); + // boost particle 1 into pair rest frame and calculate its momentum, which has the same value as k* + mTrueKstar = boostPrf(particle1Prf).P(); + } + + template + void setPairMc(const T1& particle1, const T2& particle2, const T3& mcParticles, const T4& col, const T5& /*mcCols*/) + { + setPair(particle1, particle2, col); + setPairMc(particle1, particle2, mcParticles); + if (!col.has_fMcCol()) { + mHasMcCol = false; + return; + } + mHasMcCol = true; + auto mcCol = col.template fMcCol_as(); + mTrueMult = mcCol.mult(); + mTrueCent = mcCol.cent(); + } + + template + void setPairMc(const T1& particle1, const T2& particle2, const T3& mcParticles, const T4& col1, const T5& col2, const T6& /*mcCols*/) + { + setPair(particle1, particle2, col1, col2); + setPairMc(particle1, particle2, mcParticles); + if (!col1.has_fMcCol() || !col2.has_fMcCol()) { + mHasMcCol = false; + return; + } + mHasMcCol = true; + auto mcCol1 = col1.template fMcCol_as(); + auto mcCol2 = col2.template fMcCol_as(); + mTrueMult = 0.5f * (mcCol1.mult() + mcCol2.mult()); + mTrueCent = 0.5f * (mcCol1.cent() + mcCol2.cent()); } bool checkPairCuts() const @@ -353,6 +440,9 @@ class PairHistManager if constexpr (isFlagSet(mode, modes::Mode::kAnalysis)) { fillAnalysis(); } + if constexpr (isFlagSet(mode, modes::Mode::kMc)) { + fillMc(); + } } float getKstar() const { return mKstar; } @@ -426,6 +516,16 @@ class PairHistManager } } + void initMc(std::map> const& Specs) + { + std::string mcDir = std::string(prefix) + std::string(McDir); + mHistogramRegistry->add(mcDir + getHistNameV2(kTrueKstarVsKstar, HistTable), getHistDesc(kTrueKstarVsKstar, HistTable), getHistType(kTrueKstarVsKstar, HistTable), {Specs.at(kTrueKstarVsKstar)}); + mHistogramRegistry->add(mcDir + getHistNameV2(kTrueKtVsKt, HistTable), getHistDesc(kTrueKtVsKt, HistTable), getHistType(kTrueKtVsKt, HistTable), {Specs.at(kTrueKtVsKt)}); + mHistogramRegistry->add(mcDir + getHistNameV2(kTrueMtVsMt, HistTable), getHistDesc(kTrueMtVsMt, HistTable), getHistType(kTrueMtVsMt, HistTable), {Specs.at(kTrueMtVsMt)}); + mHistogramRegistry->add(mcDir + getHistNameV2(kTrueMultVsMult, HistTable), getHistDesc(kTrueMultVsMult, HistTable), getHistType(kTrueMultVsMult, HistTable), {Specs.at(kTrueMultVsMult)}); + mHistogramRegistry->add(mcDir + getHistNameV2(kTrueCentVsCent, HistTable), getHistDesc(kTrueCentVsCent, HistTable), getHistType(kTrueCentVsCent, HistTable), {Specs.at(kTrueCentVsCent)}); + } + void fillAnalysis() { if (mPlot1d) { @@ -493,21 +593,37 @@ class PairHistManager } } - void computeMt(ROOT::Math::PtEtaPhiMVector const& PairMomentum) + void fillMc() { + if (mHasMcPair) { + mHistogramRegistry->fill(HIST(prefix) + HIST(McDir) + HIST(getHistName(kTrueKstarVsKstar, HistTable)), mTrueKstar, mKstar); + mHistogramRegistry->fill(HIST(prefix) + HIST(McDir) + HIST(getHistName(kTrueKtVsKt, HistTable)), mTrueKt, mKt); + mHistogramRegistry->fill(HIST(prefix) + HIST(McDir) + HIST(getHistName(kTrueMtVsMt, HistTable)), mTrueMt, mMt); + } + if (mHasMcCol) { + mHistogramRegistry->fill(HIST(prefix) + HIST(McDir) + HIST(getHistName(kTrueMultVsMult, HistTable)), mTrueMult, mMult); + mHistogramRegistry->fill(HIST(prefix) + HIST(McDir) + HIST(getHistName(kTrueCentVsCent, HistTable)), mTrueCent, mCent); + } + } + + float computeMt(ROOT::Math::PtEtaPhiMVector const& PairMomentum) + { + float mt = 0; switch (mMtType) { case modes::TransverseMassType::kAveragePdgMass: - mMt = std::hypot(PairMomentum.Pt() / 2.f, mAverageMass); + mt = std::hypot(0.5 * PairMomentum.Pt(), mAverageMass); break; case modes::TransverseMassType::kReducedPdgMass: - mMt = std::hypot(PairMomentum.Pt() / 2.f, mReducedMass); + mt = std::hypot(0.5 * PairMomentum.Pt(), mReducedMass); break; case modes::TransverseMassType::kMt4Vector: - mMt = PairMomentum.Mt() / 2.f; + mt = PairMomentum.Mt() / 2.f; break; default: - mMt = std::hypot(mKt, mAverageMass); + LOG(warn) << "Invalid transverse mass type, falling back to default..."; + mt = std::hypot(0.5 * PairMomentum.Pt(), mAverageMass); } + return mt; } o2::framework::HistogramRegistry* mHistogramRegistry = nullptr; @@ -530,7 +646,16 @@ class PairHistManager float mMult = 0.f; float mCent = 0.f; + // mc + float mTrueKstar = 0.f; + float mTrueKt = 0.f; + float mTrueMt = 0.f; + float mTrueMult = 0.f; + float mTrueCent = 0.f; + // cuts + bool mHasMcPair = false; + bool mHasMcCol = false; float mKstarMin = -1.f; float mKstarMax = -1.f; float mKtMin = -1.f; diff --git a/PWGCF/Femto/Core/pairProcessHelpers.h b/PWGCF/Femto/Core/pairProcessHelpers.h index 74f47920246..6380685bdc2 100644 --- a/PWGCF/Femto/Core/pairProcessHelpers.h +++ b/PWGCF/Femto/Core/pairProcessHelpers.h @@ -80,6 +80,65 @@ void processSameEvent(T1 const& SliceParticle, } } +// process same event for identical particles with mc information +template +void processSameEvent(T1 const& SliceParticle, + T2 const& TrackTable, + T3 const& mcParticles, + T4 const& mcMothers, + T5 const& mcPartonicMothers, + T6 const& Collision, + T7 const& mcCollisions, + T8& ParticleHistManager, + T9& PairHistManager, + T10& CprManager, + T11& PcManager, + T12& rng, + bool randomize) +{ + for (auto const& part : SliceParticle) { + ParticleHistManager.template fill(part, TrackTable, mcParticles, mcMothers, mcPartonicMothers); + } + std::uniform_real_distribution dist(0.f, 1.f); + for (auto const& [p1, p2] : o2::soa::combinations(o2::soa::CombinationsStrictlyUpperIndexPolicy(SliceParticle, SliceParticle))) { + // check if pair is clean + if (!PcManager.isCleanPair(p1, p2, TrackTable, mcPartonicMothers)) { + continue; + } + // check if pair is close + CprManager.setPair(p1, p2, TrackTable); + if (CprManager.isClosePair()) { + continue; + } + // Randomize pair order if enabled + float threshold = 0.5f; + bool swapPair = randomize ? (dist(rng) > threshold) : false; + if (swapPair) { + PairHistManager.setPairMc(p2, p1, mcParticles, Collision, mcCollisions); + } else { + PairHistManager.setPairMc(p1, p2, mcParticles, Collision, mcCollisions); + } + // fill deta-dphi histograms with kstar cutoff + CprManager.fill(PairHistManager.getKstar()); + // if pair cuts are configured check them before filling + if (PairHistManager.checkPairCuts()) { + PairHistManager.template fill(); + } + } +} + // process same event for non-identical particles template +void processSameEvent(T1 const& SliceParticle1, + T2 const& SliceParticle2, + T3 const& TrackTable, + T4 const& mcParticles, + T5 const& mcMothers, + T6 const& mcPartonicMothers, + T7 const& Collision, + T8 const& mcCollisions, + T9& ParticleHistManager1, + T10& ParticleHistManager2, + T11& PairHistManager, + T12& CprManager, + T13& PcManager) +{ + // Fill single particle histograms + for (auto const& part : SliceParticle1) { + ParticleHistManager1.template fill(part, TrackTable, mcParticles, mcMothers, mcPartonicMothers); + } + for (auto const& part : SliceParticle2) { + ParticleHistManager2.template fill(part, TrackTable, mcParticles, mcMothers, mcPartonicMothers); + } + for (auto const& [p1, p2] : o2::soa::combinations(o2::soa::CombinationsFullIndexPolicy(SliceParticle1, SliceParticle2))) { + // pair cleaning + if (!PcManager.isCleanPair(p1, p2, TrackTable, mcPartonicMothers)) { + continue; + } + // Close pair rejection + CprManager.setPair(p1, p2, TrackTable); + if (CprManager.isClosePair()) { + continue; + } + PairHistManager.setPairMc(p1, p2, mcParticles, Collision, mcCollisions); + CprManager.fill(PairHistManager.getKstar()); + if (PairHistManager.checkPairCuts()) { + PairHistManager.template fill(); + } + } +} + // process mixed event for identical particles template -void processMixedEvent(T1& Collisions, +void processMixedEvent(T1 const& Collisions, T2& Partition, - T3& TrackTable, + T3 const& TrackTable, T4& cache, - T5& policy, - T6& depth, + T5 const& policy, + T6 const& depth, T7& PairHistManager, T8& CprManager, T9& PcManager) @@ -176,6 +289,60 @@ void processMixedEvent(T1& Collisions, } } +// process mixed event for identical particles with mc information +template +void processMixedEvent(T1 const& Collisions, + T2 const& mcCollisions, + T3& Partition, + T4 const& TrackTable, + T5 const& mcParticles, + T6& cache, + T7 const& policy, + T8 const& depth, + T9& PairHistManager, + T10& CprManager, + T11& PcManager) +{ + for (auto const& [collision1, collision2] : o2::soa::selfCombinations(policy, depth, -1, Collisions, Collisions)) { + if (!(std::fabs(collision1.magField() - collision2.magField()) < o2::constants::math::Epsilon)) { + continue; + } + CprManager.setMagField(collision1.magField()); + auto sliceParticle1 = Partition->sliceByCached(o2::aod::femtobase::stored::fColId, collision1.globalIndex(), cache); + auto sliceParticle2 = Partition->sliceByCached(o2::aod::femtobase::stored::fColId, collision2.globalIndex(), cache); + if (sliceParticle1.size() == 0 || sliceParticle2.size() == 0) { + continue; + } + for (auto const& [p1, p2] : o2::soa::combinations(o2::soa::CombinationsFullIndexPolicy(sliceParticle1, sliceParticle2))) { + // pair cleaning + if (!PcManager.isCleanPair(p1, p2, TrackTable)) { + continue; + } + // Close pair rejection + CprManager.setPair(p1, p2, TrackTable); + if (CprManager.isClosePair()) { + continue; + } + PairHistManager.setPairMc(p1, p2, mcParticles, collision1, collision2, mcCollisions); + CprManager.fill(PairHistManager.getKstar()); + if (PairHistManager.checkPairCuts()) { + PairHistManager.template fill(); + } + } + } +} + // process mixed event different particles template -void processMixedEvent(T1& Collisions, +void processMixedEvent(T1 const& Collisions, T2& Partition1, T3& Partition2, - T4& TrackTable, + T4 const& TrackTable, T5& cache, - T6& policy, - T7& depth, + T6 const& policy, + T7 const& depth, T8& PairHistManager, T9& CprManager, T10& PcManager) @@ -227,6 +394,63 @@ void processMixedEvent(T1& Collisions, } } } + +// process mixed event different particles with mc information +template +void processMixedEvent(T1 const& Collisions, + T2 const& mcCollisions, + T3& Partition1, + T4& Partition2, + T5 const& TrackTable, + T6 const& mcParticles, + T7& cache, + T8& policy, + T9& depth, + T10& PairHistManager, + T11& CprManager, + T12& PcManager) +{ + for (auto const& [collision1, collision2] : o2::soa::selfCombinations(policy, depth, -1, Collisions, Collisions)) { + if (!(std::fabs(collision1.magField() - collision2.magField()) < o2::constants::math::Epsilon)) { + continue; + } + CprManager.setMagField(collision1.magField()); + auto sliceParticle1 = Partition1->sliceByCached(o2::aod::femtobase::stored::fColId, collision1.globalIndex(), cache); + auto sliceParticle2 = Partition2->sliceByCached(o2::aod::femtobase::stored::fColId, collision2.globalIndex(), cache); + if (sliceParticle1.size() == 0 || sliceParticle2.size() == 0) { + continue; + } + for (auto const& [p1, p2] : o2::soa::combinations(o2::soa::CombinationsFullIndexPolicy(sliceParticle1, sliceParticle2))) { + // pair cleaning + if (!PcManager.isCleanPair(p1, p2, TrackTable)) { + continue; + } + // Close pair rejection + CprManager.setPair(p1, p2, TrackTable); + if (CprManager.isClosePair()) { + continue; + } + PairHistManager.setPairMc(p1, p2, mcParticles, collision1, collision2, mcCollisions); + CprManager.fill(PairHistManager.getKstar()); + if (PairHistManager.checkPairCuts()) { + PairHistManager.template fill(); + } + } + } +} + } // namespace pairprocesshelpers } // namespace o2::analysis::femto diff --git a/PWGCF/Femto/Core/trackBuilder.h b/PWGCF/Femto/Core/trackBuilder.h index 4316a993bcd..7e1c90317bc 100644 --- a/PWGCF/Femto/Core/trackBuilder.h +++ b/PWGCF/Femto/Core/trackBuilder.h @@ -65,7 +65,7 @@ struct ConfTrackBits : o2::framework::ConfigurableGroup { o2::framework::Configurable> dcazMax{"dcazMax", {"0.004 + 0.013*pow(x, -1)"}, "Maximum |dca_z| as a function of pT. Has to be a valid TForumal, where x=pt"}; // Electron PID cuts - o2::framework::Configurable pidIsOptionalElectron{"pidIsOptionalElectron", false, "Make election PID optional"}; + o2::framework::Configurable requirePidElectron{"requirePidElectron", false, "Make election PID optional"}; o2::framework::Configurable minMomTofElectron{"minMomTofElectron", 0.3, "Minimum momentum to required TOF PID for Electron"}; o2::framework::Configurable> itsElectron{"itsElectron", {}, "Maximum |nsigma| for Electron PID"}; o2::framework::Configurable> tpcElectron{"tpcElectron", {}, "Maximum |nsigma| for Electron PID"}; @@ -74,7 +74,7 @@ struct ConfTrackBits : o2::framework::ConfigurableGroup { o2::framework::Configurable> tpctofElectron{"tpctofElectron", {}, "Maximum |nsigma| for Electron PID"}; // Pion PID cuts - o2::framework::Configurable pidIsOptionalPion{"pidIsOptionalPion", false, "Make election PID optional"}; + o2::framework::Configurable requirePidPion{"requirePidPion", false, "Make election PID optional"}; o2::framework::Configurable minMomTofPion{"minMomTofPion", 0.5, "Minimum momentum to required TOF PID for Pion"}; o2::framework::Configurable> itsPion{"itsPion", {}, "Maximum |nsigma| for Pion PID"}; o2::framework::Configurable> tpcPion{"tpcPion", {}, "Maximum |nsigma| for Pion PID"}; @@ -83,7 +83,7 @@ struct ConfTrackBits : o2::framework::ConfigurableGroup { o2::framework::Configurable> tpctofPion{"tpctofPion", {}, "Maximum |nsigma| for Pion PID"}; // Kaon PID cuts - o2::framework::Configurable pidIsOptionalKaon{"pidIsOptionalKaon", false, "Make election PID optional"}; + o2::framework::Configurable requirePidKaon{"requirePidKaon", false, "Make election PID optional"}; o2::framework::Configurable minMomTofKaon{"minMomTofKaon", 0.4, "Minimum momentum to required TOF PID for Kaon"}; o2::framework::Configurable> itsKaon{"itsKaon", {}, "Maximum |nsigma| for Kaon PID"}; o2::framework::Configurable> tpcKaon{"tpcKaon", {}, "Maximum |nsigma| for Kaon PID"}; @@ -92,7 +92,7 @@ struct ConfTrackBits : o2::framework::ConfigurableGroup { o2::framework::Configurable> tpctofKaon{"tpctofKaon", {}, "Maximum |nsigma| for Kaon PID"}; // Proton PID cuts - o2::framework::Configurable pidIsOptionalProton{"pidIsOptionalProton", true, "Make election PID optional"}; + o2::framework::Configurable requirePidProton{"requirePidProton", true, "Make election PID optional"}; o2::framework::Configurable minMomTofProton{"minMomTofProton", 0.75, "Minimum momentum to required TOF PID for Proton"}; o2::framework::Configurable> itsProton{"itsProton", {}, "Maximum |nsigma| for Proton PID"}; o2::framework::Configurable> tpcProton{"tpcProton", {}, "Maximum |nsigma| for Proton PID"}; @@ -101,7 +101,7 @@ struct ConfTrackBits : o2::framework::ConfigurableGroup { o2::framework::Configurable> tpctofProton{"tpctofProton", {3.f}, "Maximum |nsigma| for Proton PID"}; // Deuteron PID cuts - o2::framework::Configurable pidIsOptionalDeuteron{"pidIsOptionalDeuteron", false, "Make election PID optional"}; + o2::framework::Configurable requirePidDeuteron{"requirePidDeuteron", false, "Make election PID optional"}; o2::framework::Configurable minMomTofDeuteron{"minMomTofDeuteron", 1.2, "Minimum momentum to required TOF PID for Deuteron"}; o2::framework::Configurable> itsDeuteron{"itsDeuteron", {}, "Maximum |nsigma| for Deuteron PID"}; o2::framework::Configurable> tpcDeuteron{"tpcDeuteron", {}, "Maximum |nsigma| for Deuteron PID"}; @@ -110,7 +110,7 @@ struct ConfTrackBits : o2::framework::ConfigurableGroup { o2::framework::Configurable> tpctofDeuteron{"tpctofDeuteron", {}, "Maximum |nsigma| for Deuteron PID"}; // Triton PID cuts - o2::framework::Configurable pidIsOptionalTriton{"pidIsOptionalTriton", false, "Make election PID optional"}; + o2::framework::Configurable requirePidTriton{"requirePidTriton", false, "Make election PID optional"}; o2::framework::Configurable minMomTofTriton{"minMomTofTriton", 1.4, "Minimum momentum to required TOF PID for Triton"}; o2::framework::Configurable> itsTriton{"itsTriton", {}, "Maximum |nsigma| for Triton PID"}; o2::framework::Configurable> tpcTriton{"tpcTriton", {}, "Maximum |nsigma| for Triton PID"}; @@ -119,7 +119,7 @@ struct ConfTrackBits : o2::framework::ConfigurableGroup { o2::framework::Configurable> tpctofTriton{"tpctofTriton", {}, "Maximum |nsigma| for Triton PID"}; // Helium PID cuts - o2::framework::Configurable pidIsOptionalHelium{"pidIsOptionalHelium", false, "Make election PID optional"}; + o2::framework::Configurable requirePidHelium{"requirePidHelium", false, "Make election PID optional"}; o2::framework::Configurable minMomTofHelium{"minMomTofHelium", 1.6, "Minimum momentum to required TOF PID for Helium"}; o2::framework::Configurable> itsHelium{"itsHelium", {}, "Maximum |nsigma| for Helium PID"}; o2::framework::Configurable> tpcHelium{"tpcHelium", {}, "Maximum |nsigma| for Helium PID"}; @@ -133,8 +133,8 @@ template struct ConfTrackSelection : public o2::framework::ConfigurableGroup { std::string prefix = Prefix; // Unique prefix based on the template argument // configuration parameters - o2::framework::Configurable pdgCode{"pdgCode", 2212, "Track PDG code"}; - o2::framework::Configurable chargeAbs{"chargeAbs", 1, "Absolute value of charge (e.g. 1 for most tracks, 2 for He3)"}; + o2::framework::Configurable pdgCodeAbs{"pdgCodeAbs", 2212, "Absolute value of PDG code. Set sign of charge to -1 for antiparticle."}; + o2::framework::Configurable chargeAbs{"chargeAbs", 1, "Absolute value of charge (e.g. 1 for most tracks, 2 for He3). Set sign of charge to -1 for antiparticle"}; o2::framework::Configurable chargeSign{"chargeSign", 1, "Track charge sign: +1 for positive, -1 for negative, 0 for both"}; // filters for kinematics o2::framework::Configurable ptMin{"ptMin", 0.2f, "Minimum pT (GeV/c)"}; @@ -302,59 +302,59 @@ class TrackSelection : public BaseSelectionaddSelection(kDCAzMax, trackSelectionNames.at(kDCAzMax), filter.ptMin.value, filter.ptMax.value, config.dcazMax.value, limits::kAbsUpperFunctionLimit, true, true, false); // add selections for Electron pid - this->addSelection(kItsElectron, trackSelectionNames.at(kItsElectron), config.itsElectron.value, limits::kAbsUpperLimit, false, false, config.pidIsOptionalElectron); - this->addSelection(kTpcElectron, trackSelectionNames.at(kTpcElectron), config.tpcElectron.value, limits::kAbsUpperLimit, false, false, config.pidIsOptionalElectron); - this->addSelection(kTofElectron, trackSelectionNames.at(kTofElectron), config.tofElectron.value, limits::kAbsUpperLimit, false, false, config.pidIsOptionalElectron); - this->addSelection(kTpcitsElectron, trackSelectionNames.at(kTpcitsElectron), config.tpcitsElectron.value, limits::kUpperLimit, false, false, config.pidIsOptionalElectron); - this->addSelection(kTpctofElectron, trackSelectionNames.at(kTpctofElectron), config.tpctofElectron.value, limits::kUpperLimit, false, false, config.pidIsOptionalElectron); + this->addSelection(kItsElectron, trackSelectionNames.at(kItsElectron), config.itsElectron.value, limits::kAbsUpperLimit, false, false, config.requirePidElectron); + this->addSelection(kTpcElectron, trackSelectionNames.at(kTpcElectron), config.tpcElectron.value, limits::kAbsUpperLimit, false, false, config.requirePidElectron); + this->addSelection(kTofElectron, trackSelectionNames.at(kTofElectron), config.tofElectron.value, limits::kAbsUpperLimit, false, false, config.requirePidElectron); + this->addSelection(kTpcitsElectron, trackSelectionNames.at(kTpcitsElectron), config.tpcitsElectron.value, limits::kUpperLimit, false, false, config.requirePidElectron); + this->addSelection(kTpctofElectron, trackSelectionNames.at(kTpctofElectron), config.tpctofElectron.value, limits::kUpperLimit, false, false, config.requirePidElectron); mElectronTofThres = config.minMomTofElectron.value; // add selections for Pion pid - this->addSelection(kItsPion, trackSelectionNames.at(kItsPion), config.itsPion.value, limits::kAbsUpperLimit, false, false, config.pidIsOptionalPion); - this->addSelection(kTpcPion, trackSelectionNames.at(kTpcPion), config.tpcPion.value, limits::kAbsUpperLimit, false, false, config.pidIsOptionalPion); - this->addSelection(kTofPion, trackSelectionNames.at(kTofPion), config.tofPion.value, limits::kAbsUpperLimit, false, false, config.pidIsOptionalPion); - this->addSelection(kTpcitsPion, trackSelectionNames.at(kTpcitsPion), config.tpcitsPion.value, limits::kUpperLimit, false, false, config.pidIsOptionalPion); - this->addSelection(kTpctofPion, trackSelectionNames.at(kTpctofPion), config.tpctofPion.value, limits::kUpperLimit, false, false, config.pidIsOptionalPion); + this->addSelection(kItsPion, trackSelectionNames.at(kItsPion), config.itsPion.value, limits::kAbsUpperLimit, false, false, config.requirePidPion); + this->addSelection(kTpcPion, trackSelectionNames.at(kTpcPion), config.tpcPion.value, limits::kAbsUpperLimit, false, false, config.requirePidPion); + this->addSelection(kTofPion, trackSelectionNames.at(kTofPion), config.tofPion.value, limits::kAbsUpperLimit, false, false, config.requirePidPion); + this->addSelection(kTpcitsPion, trackSelectionNames.at(kTpcitsPion), config.tpcitsPion.value, limits::kUpperLimit, false, false, config.requirePidPion); + this->addSelection(kTpctofPion, trackSelectionNames.at(kTpctofPion), config.tpctofPion.value, limits::kUpperLimit, false, false, config.requirePidPion); mPionTofThres = config.minMomTofPion.value; // add selections for Kaon pid - this->addSelection(kItsKaon, trackSelectionNames.at(kItsKaon), config.itsKaon.value, limits::kAbsUpperLimit, false, false, config.pidIsOptionalKaon); - this->addSelection(kTpcKaon, trackSelectionNames.at(kTpcKaon), config.tpcKaon.value, limits::kAbsUpperLimit, false, false, config.pidIsOptionalKaon); - this->addSelection(kTofKaon, trackSelectionNames.at(kTofKaon), config.tofKaon.value, limits::kAbsUpperLimit, false, false, config.pidIsOptionalKaon); - this->addSelection(kTpcitsKaon, trackSelectionNames.at(kTpcitsKaon), config.tpcitsKaon.value, limits::kUpperLimit, false, false, config.pidIsOptionalKaon); - this->addSelection(kTpctofKaon, trackSelectionNames.at(kTpctofKaon), config.tpctofKaon.value, limits::kUpperLimit, false, false, config.pidIsOptionalKaon); + this->addSelection(kItsKaon, trackSelectionNames.at(kItsKaon), config.itsKaon.value, limits::kAbsUpperLimit, false, false, config.requirePidKaon); + this->addSelection(kTpcKaon, trackSelectionNames.at(kTpcKaon), config.tpcKaon.value, limits::kAbsUpperLimit, false, false, config.requirePidKaon); + this->addSelection(kTofKaon, trackSelectionNames.at(kTofKaon), config.tofKaon.value, limits::kAbsUpperLimit, false, false, config.requirePidKaon); + this->addSelection(kTpcitsKaon, trackSelectionNames.at(kTpcitsKaon), config.tpcitsKaon.value, limits::kUpperLimit, false, false, config.requirePidKaon); + this->addSelection(kTpctofKaon, trackSelectionNames.at(kTpctofKaon), config.tpctofKaon.value, limits::kUpperLimit, false, false, config.requirePidKaon); mKaonTofThres = config.minMomTofKaon.value; // add selections for Proton pid - this->addSelection(kItsProton, trackSelectionNames.at(kItsProton), config.itsProton.value, limits::kAbsUpperLimit, false, false, config.pidIsOptionalProton); - this->addSelection(kTpcProton, trackSelectionNames.at(kTpcProton), config.tpcProton.value, limits::kAbsUpperLimit, false, false, config.pidIsOptionalProton); - this->addSelection(kTofProton, trackSelectionNames.at(kTofProton), config.tofProton.value, limits::kAbsUpperLimit, false, false, config.pidIsOptionalProton); - this->addSelection(kTpcitsProton, trackSelectionNames.at(kTpcitsProton), config.tpcitsProton.value, limits::kUpperLimit, false, false, config.pidIsOptionalProton); - this->addSelection(kTpctofProton, trackSelectionNames.at(kTpctofProton), config.tpctofProton.value, limits::kUpperLimit, false, false, config.pidIsOptionalProton); + this->addSelection(kItsProton, trackSelectionNames.at(kItsProton), config.itsProton.value, limits::kAbsUpperLimit, false, false, config.requirePidProton); + this->addSelection(kTpcProton, trackSelectionNames.at(kTpcProton), config.tpcProton.value, limits::kAbsUpperLimit, false, false, config.requirePidProton); + this->addSelection(kTofProton, trackSelectionNames.at(kTofProton), config.tofProton.value, limits::kAbsUpperLimit, false, false, config.requirePidProton); + this->addSelection(kTpcitsProton, trackSelectionNames.at(kTpcitsProton), config.tpcitsProton.value, limits::kUpperLimit, false, false, config.requirePidProton); + this->addSelection(kTpctofProton, trackSelectionNames.at(kTpctofProton), config.tpctofProton.value, limits::kUpperLimit, false, false, config.requirePidProton); mProtonTofThres = config.minMomTofProton.value; // add selections for Deuteron pid - this->addSelection(kItsDeuteron, trackSelectionNames.at(kItsDeuteron), config.itsDeuteron.value, limits::kAbsUpperLimit, false, false, config.pidIsOptionalDeuteron); - this->addSelection(kTpcDeuteron, trackSelectionNames.at(kTpcDeuteron), config.tpcDeuteron.value, limits::kAbsUpperLimit, false, false, config.pidIsOptionalDeuteron); - this->addSelection(kTofDeuteron, trackSelectionNames.at(kTofDeuteron), config.tofDeuteron.value, limits::kAbsUpperLimit, false, false, config.pidIsOptionalDeuteron); - this->addSelection(kTpcitsDeuteron, trackSelectionNames.at(kTpcitsDeuteron), config.tpcitsDeuteron.value, limits::kUpperLimit, false, false, config.pidIsOptionalDeuteron); - this->addSelection(kTpctofDeuteron, trackSelectionNames.at(kTpctofDeuteron), config.tpctofDeuteron.value, limits::kUpperLimit, false, false, config.pidIsOptionalDeuteron); + this->addSelection(kItsDeuteron, trackSelectionNames.at(kItsDeuteron), config.itsDeuteron.value, limits::kAbsUpperLimit, false, false, config.requirePidDeuteron); + this->addSelection(kTpcDeuteron, trackSelectionNames.at(kTpcDeuteron), config.tpcDeuteron.value, limits::kAbsUpperLimit, false, false, config.requirePidDeuteron); + this->addSelection(kTofDeuteron, trackSelectionNames.at(kTofDeuteron), config.tofDeuteron.value, limits::kAbsUpperLimit, false, false, config.requirePidDeuteron); + this->addSelection(kTpcitsDeuteron, trackSelectionNames.at(kTpcitsDeuteron), config.tpcitsDeuteron.value, limits::kUpperLimit, false, false, config.requirePidDeuteron); + this->addSelection(kTpctofDeuteron, trackSelectionNames.at(kTpctofDeuteron), config.tpctofDeuteron.value, limits::kUpperLimit, false, false, config.requirePidDeuteron); mDeuteronTofThres = config.minMomTofDeuteron.value; // add selections for Triton pid - this->addSelection(kItsTriton, trackSelectionNames.at(kItsTriton), config.itsTriton.value, limits::kAbsUpperLimit, false, false, config.pidIsOptionalTriton); - this->addSelection(kTpcTriton, trackSelectionNames.at(kTpcTriton), config.tpcTriton.value, limits::kAbsUpperLimit, false, false, config.pidIsOptionalTriton); - this->addSelection(kTofTriton, trackSelectionNames.at(kTofTriton), config.tofTriton.value, limits::kAbsUpperLimit, false, false, config.pidIsOptionalTriton); - this->addSelection(kTpcitsTriton, trackSelectionNames.at(kTpcitsTriton), config.tpcitsTriton.value, limits::kUpperLimit, false, false, config.pidIsOptionalTriton); - this->addSelection(kTpctofTriton, trackSelectionNames.at(kTpctofTriton), config.tpctofTriton.value, limits::kUpperLimit, false, false, config.pidIsOptionalTriton); + this->addSelection(kItsTriton, trackSelectionNames.at(kItsTriton), config.itsTriton.value, limits::kAbsUpperLimit, false, false, config.requirePidTriton); + this->addSelection(kTpcTriton, trackSelectionNames.at(kTpcTriton), config.tpcTriton.value, limits::kAbsUpperLimit, false, false, config.requirePidTriton); + this->addSelection(kTofTriton, trackSelectionNames.at(kTofTriton), config.tofTriton.value, limits::kAbsUpperLimit, false, false, config.requirePidTriton); + this->addSelection(kTpcitsTriton, trackSelectionNames.at(kTpcitsTriton), config.tpcitsTriton.value, limits::kUpperLimit, false, false, config.requirePidTriton); + this->addSelection(kTpctofTriton, trackSelectionNames.at(kTpctofTriton), config.tpctofTriton.value, limits::kUpperLimit, false, false, config.requirePidTriton); mTritonTofThres = config.minMomTofTriton.value; // add selections for Helium pid - this->addSelection(kItsHelium, trackSelectionNames.at(kItsHelium), config.itsHelium.value, limits::kAbsUpperLimit, false, false, config.pidIsOptionalHelium); - this->addSelection(kTpcHelium, trackSelectionNames.at(kTpcHelium), config.tpcHelium.value, limits::kAbsUpperLimit, false, false, config.pidIsOptionalHelium); - this->addSelection(kTofHelium, trackSelectionNames.at(kTofHelium), config.tofHelium.value, limits::kAbsUpperLimit, false, false, config.pidIsOptionalHelium); - this->addSelection(kTpcitsHelium, trackSelectionNames.at(kTpcitsHelium), config.tpcitsHelium.value, limits::kUpperLimit, false, false, config.pidIsOptionalHelium); - this->addSelection(kTpctofHelium, trackSelectionNames.at(kTpctofHelium), config.tpctofHelium.value, limits::kUpperLimit, false, false, config.pidIsOptionalHelium); + this->addSelection(kItsHelium, trackSelectionNames.at(kItsHelium), config.itsHelium.value, limits::kAbsUpperLimit, false, false, config.requirePidHelium); + this->addSelection(kTpcHelium, trackSelectionNames.at(kTpcHelium), config.tpcHelium.value, limits::kAbsUpperLimit, false, false, config.requirePidHelium); + this->addSelection(kTofHelium, trackSelectionNames.at(kTofHelium), config.tofHelium.value, limits::kAbsUpperLimit, false, false, config.requirePidHelium); + this->addSelection(kTpcitsHelium, trackSelectionNames.at(kTpcitsHelium), config.tpcitsHelium.value, limits::kUpperLimit, false, false, config.requirePidHelium); + this->addSelection(kTpctofHelium, trackSelectionNames.at(kTpctofHelium), config.tpctofHelium.value, limits::kUpperLimit, false, false, config.requirePidHelium); mHeliumTofThres = config.minMomTofHelium.value; this->setupContainers(registry); diff --git a/PWGCF/Femto/Core/trackHistManager.h b/PWGCF/Femto/Core/trackHistManager.h index d8d68abac2a..4977aa4bd30 100644 --- a/PWGCF/Femto/Core/trackHistManager.h +++ b/PWGCF/Femto/Core/trackHistManager.h @@ -119,14 +119,15 @@ enum TrackHist { kPdg, kPdgMother, kPdgPartonicMother, - kTruePt, - kTrueEta, - kTruePhi, + kTruePtVsPt, + kTrueEtaVsEta, + kTruePhiVsPhi, // histograms for fraction estimation of tracks kNoMcParticle, kPrimary, kFromWrongCollision, kFromMaterial, + kMissidentified, kSecondary1, kSecondary2, kSecondary3, @@ -144,6 +145,7 @@ struct ConfTrackBinning : o2::framework::ConfigurableGroup { o2::framework::ConfigurableAxis eta{"eta", {{300, -1.5, 1.5}}, "Eta"}; o2::framework::ConfigurableAxis phi{"phi", {{720, 0, 1.f * o2::constants::math::TwoPI}}, "Phi"}; o2::framework::ConfigurableAxis sign{"sign", {{3, -1.5, 1.5}}, "Sign"}; + o2::framework::ConfigurableAxis pdgCodes{"pdgCodes", {{8001, -4000.5, 4000.5}}, "MC ONLY: PDG codes of selected tracks"}; }; constexpr const char PrefixTrackBinning1[] = "TrackBinning1"; @@ -180,6 +182,8 @@ struct ConfTrackQaBinning : o2::framework::ConfigurableGroup { o2::framework::Configurable plotDeuteronPid{"plotDeuteronPid", true, "Generate plots for Deuteron PID"}; o2::framework::Configurable plotTritonPid{"plotTritonPid", true, "Generate plots for Triton PID"}; o2::framework::Configurable plotHeliumPid{"plotHeliumPid", true, "Generate plots for Helium PID"}; + o2::framework::Configurable plotOrigins{"plotOrigins", true, "MC ONLY: Plot pt vs DCAxy vs DCAz for different particle origins"}; + o2::framework::Configurable> pdgCodesForMothersOfSecondary{"pdgCodesForMothersOfSecondary", {3122}, "MC ONLY: PDG codes of mothers of secondaries (Max 3 will be considered)"}; o2::framework::ConfigurableAxis itsCluster{"itsCluster", {{8, -0.5, 7.5}}, "ITS cluster"}; o2::framework::ConfigurableAxis itsClusterIb{"itsClusterIb", {{4, -0.5, 3.5}}, "ITS cluster in inner barrel"}; o2::framework::ConfigurableAxis tpcCrossedRows{"tpcCrossedRows", {{161, -0.5, 160.5}}, "TPC cluster"}; @@ -254,28 +258,6 @@ using ConfCascadeNegDauQaBinning = ConfTrackQaBinning; using ConfKinkChaDauQaBinning = ConfTrackQaBinning; -template -struct ConfTrackMcBinning : o2::framework::ConfigurableGroup { - std::string prefix = std::string(Prefix); - o2::framework::ConfigurableAxis pdgCodes{"pdgCodes", {{8001, -4000.5, 4000.5}}, "PDG codes of selected tracks"}; - o2::framework::ConfigurableAxis statusCode{"statusCode", {{21, -0.5, 20.5}}, "Status codes (i.e. Origin)"}; - o2::framework::Configurable plotOrigins{"plotOrigins", true, "Plot pt vs DCAxy vs DCAz for different particle origins"}; - o2::framework::Configurable> pdgCodesForMothersOfSecondary{"pdgCodesForMothersOfSecondary", {3122}, "PDG codes of mothers of secondaries (Max 3 will be considered)"}; - o2::framework::ConfigurableAxis pt{"pt", {{150, 0, 3}}, "Pt"}; - o2::framework::ConfigurableAxis dcaXy{"dcaXy", {{50, -0.1, 0.1}}, "DCA_xy"}; - o2::framework::ConfigurableAxis dcaZ{"dcaZ", {{50, -0.1, 0.1}}, "DCA_z"}; -}; - -constexpr const char PrefixTrackMcBinning1[] = "TrackMcBinning1"; -constexpr const char PrefixV0PosDauMcBinning[] = "V0PosDauMcBinning"; -constexpr const char PrefixV0NegDauMcBinning[] = "V0NegDauMcBinning"; -constexpr const char PrefixKinkChaDauMcBinning[] = "KinkChaDauMcBinning"; - -using ConfTrackMcBinning1 = ConfTrackMcBinning; -using ConfV0PosDauMcBinning = ConfTrackMcBinning; -using ConfV0NegDauMcBinning = ConfTrackMcBinning; -using ConfKinkChaDauMcBinning = ConfTrackMcBinning; - // must be in sync with enum TrackVariables // the enum gives the correct index in the array constexpr std::array, kTrackHistLast> @@ -351,13 +333,14 @@ constexpr std::array, kTrackHistLast> {kPdg, o2::framework::kTH1F, "hPdg", "PDG Codes of selected tracks; PDG Code; Entries"}, {kPdgMother, o2::framework::kTH1F, "hPdgMother", "PDG Codes of mother of selected tracks; PDG Code; Entries"}, {kPdgPartonicMother, o2::framework::kTH1F, "hPdgPartonicMother", "PDG Codes of partonic mother selected tracks; PDG Code; Entries"}, - {kTruePt, o2::framework::kTH1F, "hTruePt", "True transverse momentum; p_{T} (GeV/#it{c}); Entries"}, - {kTrueEta, o2::framework::kTH1F, "hTrueEta", "True pseudorapdity; #eta; Entries"}, - {kTruePhi, o2::framework::kTH1F, "hTruePhi", "True azimuthal angle; #varphi; Entries"}, + {kTruePtVsPt, o2::framework::kTH2F, "hTruePtVsPt", "True transverse momentum vs transverse momentum; p_{T,True} (GeV/#it{c}); p_{T,True} (GeV/#it{c})"}, + {kTrueEtaVsEta, o2::framework::kTH2F, "hTrueEtaVsEta", "True pseudorapdity vs pseudorapdity; #eta_{True}; #eta"}, + {kTruePhiVsPhi, o2::framework::kTH2F, "hTruePhiVsPhi", "True azimuthal angle vs azimuthal angle; #varphi_{True}; #varphi"}, {kNoMcParticle, o2::framework::kTHnSparseF, "hNoMcParticle", "Wrongly reconstructed particles; p_{T} (GeV/#it{c}); DCA_{xy} (cm); DCA_{z} (cm)"}, {kPrimary, o2::framework::kTHnSparseF, "hPrimary", "Primary particles; p_{T} (GeV/#it{c}); DCA_{xy} (cm); DCA_{z} (cm)"}, {kFromWrongCollision, o2::framework::kTHnSparseF, "hFromWrongCollision", "Particles associated to wrong collision; p_{T} (GeV/#it{c}); DCA_{xy} (cm); DCA_{z} (cm)"}, {kFromMaterial, o2::framework::kTHnSparseF, "hFromMaterial", "Particles from material; p_{T} (GeV/#it{c}); DCA_{xy} (cm); DCA_{z} (cm)"}, + {kMissidentified, o2::framework::kTHnSparseF, "hMissidentified", "Missidentified particles (fake/wrong PDG code); p_{T} (GeV/#it{c}); DCA_{xy} (cm); DCA_{z} (cm)"}, {kSecondary1, o2::framework::kTHnSparseF, "hFromSecondary1", "Particles from secondary decay; p_{T} (GeV/#it{c}); DCA_{xy} (cm); DCA_{z} (cm)"}, {kSecondary2, o2::framework::kTHnSparseF, "hFromSecondary2", "Particles from seconary decay; p_{T} (GeV/#it{c}); DCA_{xy} (cm); DCA_{z} (cm)"}, {kSecondary3, o2::framework::kTHnSparseF, "hFromSecondary3", "Particles from seconary decay; p_{T} (GeV/#it{c}); DCA_{xy} (cm); DCA_{z} (cm)"}, @@ -434,22 +417,24 @@ constexpr std::array, kTrackHistLast> {kTpctofTriton, {confQa.p, confQa.tpctofTriton}}, \ {kTpctofHelium, {confQa.p, confQa.tpctofHelium}}, -#define TRACK_HIST_MC_MAP(confAnalysis, confMc) \ - {kTruePt, {confAnalysis.pt}}, \ - {kTrueEta, {confAnalysis.eta}}, \ - {kTruePhi, {confAnalysis.phi}}, \ - {kOrigin, {confMc.statusCode}}, \ - {kPdg, {confMc.pdgCodes}}, \ - {kPdgMother, {confMc.pdgCodes}}, \ - {kPdgPartonicMother, {confMc.pdgCodes}}, \ - {kNoMcParticle, {confMc.pt, confMc.dcaXy, confMc.dcaZ}}, \ - {kPrimary, {confMc.pt, confMc.dcaXy, confMc.dcaZ}}, \ - {kFromWrongCollision, {confMc.pt, confMc.dcaXy, confMc.dcaZ}}, \ - {kFromMaterial, {confMc.pt, confMc.dcaXy, confMc.dcaZ}}, \ - {kSecondary1, {confMc.pt, confMc.dcaXy, confMc.dcaZ}}, \ - {kSecondary2, {confMc.pt, confMc.dcaXy, confMc.dcaZ}}, \ - {kSecondary3, {confMc.pt, confMc.dcaXy, confMc.dcaZ}}, \ - {kSecondaryOther, {confMc.pt, confMc.dcaXy, confMc.dcaZ}}, +#define TRACK_HIST_MC_MAP(conf) \ + {kTruePtVsPt, {conf.pt, conf.pt}}, \ + {kTrueEtaVsEta, {conf.eta, conf.eta}}, \ + {kTruePhiVsPhi, {conf.phi, conf.phi}}, \ + {kPdg, {conf.pdgCodes}}, \ + {kPdgMother, {conf.pdgCodes}}, \ + {kPdgPartonicMother, {conf.pdgCodes}}, + +#define TRACK_HIST_MC_QA_MAP(confAnalysis, confQa) \ + {kNoMcParticle, {confAnalysis.pt, confQa.dcaXy, confQa.dcaZ}}, \ + {kPrimary, {confAnalysis.pt, confQa.dcaXy, confQa.dcaZ}}, \ + {kFromWrongCollision, {confAnalysis.pt, confQa.dcaXy, confQa.dcaZ}}, \ + {kFromMaterial, {confAnalysis.pt, confQa.dcaXy, confQa.dcaZ}}, \ + {kMissidentified, {confAnalysis.pt, confQa.dcaXy, confQa.dcaZ}}, \ + {kSecondary1, {confAnalysis.pt, confQa.dcaXy, confQa.dcaZ}}, \ + {kSecondary2, {confAnalysis.pt, confQa.dcaXy, confQa.dcaZ}}, \ + {kSecondary3, {confAnalysis.pt, confQa.dcaXy, confQa.dcaZ}}, \ + {kSecondaryOther, {confAnalysis.pt, confQa.dcaXy, confQa.dcaZ}}, template auto makeTrackHistSpecMap(const T& confBinningAnalysis) @@ -458,6 +443,14 @@ auto makeTrackHistSpecMap(const T& confBinningAnalysis) TRACK_HIST_ANALYSIS_MAP(confBinningAnalysis)}; } +template +auto makeTrackMcHistSpecMap(T const& confBinningAnalysis) +{ + return std::map>{ + TRACK_HIST_ANALYSIS_MAP(confBinningAnalysis) + TRACK_HIST_MC_MAP(confBinningAnalysis)}; +}; + template auto makeTrackQaHistSpecMap(T1 const& confBinningAnalysis, T2 const& confBinningQa) { @@ -466,18 +459,20 @@ auto makeTrackQaHistSpecMap(T1 const& confBinningAnalysis, T2 const& confBinning TRACK_HIST_QA_MAP(confBinningAnalysis, confBinningQa)}; } -template -auto makeTrackMcQaHistSpecMap(T1 const& confBinningAnalysis, T2 const& confBinningQa, T3 const& confBinningMc) +template +auto makeTrackMcQaHistSpecMap(T1 const& confBinningAnalysis, T2 const& confBinningQa) { return std::map>{ TRACK_HIST_ANALYSIS_MAP(confBinningAnalysis) TRACK_HIST_QA_MAP(confBinningAnalysis, confBinningQa) - TRACK_HIST_MC_MAP(confBinningAnalysis, confBinningMc)}; + TRACK_HIST_MC_MAP(confBinningAnalysis) + TRACK_HIST_MC_QA_MAP(confBinningAnalysis, confBinningQa)}; }; #undef TRACK_HIST_ANALYSIS_MAP #undef TRACK_HIST_QA_MAP #undef TRACK_HIST_MC_MAP +#undef TRACK_HIST_MC_QA_MAP constexpr char PrefixTrackQa[] = "TrackQA/"; constexpr char PrefixTrack1[] = "Track1/"; @@ -518,12 +513,15 @@ class TrackHistManager TrackHistManager() = default; ~TrackHistManager() = default; - // init for analysis - template - void init(o2::framework::HistogramRegistry* registry, std::map> const& Specs, int AbsCharge) + // init for analysis and mc + template + void init(o2::framework::HistogramRegistry* registry, + std::map> const& Specs, + T const& ConfTrackSelection) { mHistogramRegistry = registry; - mAbsCharge = std::abs(AbsCharge); + mAbsCharge = std::abs(ConfTrackSelection.chargeAbs.value); + mPdgCode = std::abs(ConfTrackSelection.pdgCodeAbs.value) * ConfTrackSelection.chargeSign.value; if constexpr (isFlagSet(mode, modes::Mode::kAnalysis)) { this->initAnalysis(Specs); } @@ -535,20 +533,48 @@ class TrackHistManager } } - // init for analysis and qa - template - void init(o2::framework::HistogramRegistry* registry, std::map> const& Specs, T const& ConfBinningQa, int AbsCharge) + template + void init(o2::framework::HistogramRegistry* registry, + std::map> const& Specs, + int ChargeAbs, + int ChargeSign, + int PdgCodeAbs) { - this->enableOptionalHistograms(ConfBinningQa); - this->template init(registry, Specs, AbsCharge); + mHistogramRegistry = registry; + mAbsCharge = std::abs(ChargeAbs); + mPdgCode = std::abs(PdgCodeAbs) * ChargeSign; + if constexpr (isFlagSet(mode, modes::Mode::kAnalysis)) { + this->initAnalysis(Specs); + } + if constexpr (isFlagSet(mode, modes::Mode::kQa)) { + this->initQa(Specs); + } + if constexpr (isFlagSet(mode, modes::Mode::kMc)) { + this->initMc(Specs); + } } - // init for analysis and mc and qa + // init for analysis and qa and mc template - void init(o2::framework::HistogramRegistry* registry, std::map> const& Specs, T1 const& ConfBinningQa, T2 const& ConfBinningMc, int AbsCharge) + void init(o2::framework::HistogramRegistry* registry, + std::map> const& Specs, + T1 const& ConfTrackSelection, + T2 const& ConfBinningQa) { - this->enableOptionalHistograms(ConfBinningQa, ConfBinningMc); - this->template init(registry, Specs, AbsCharge); + this->template enableOptionalHistograms(ConfBinningQa); + this->template init(registry, Specs, ConfTrackSelection); + } + + template + void init(o2::framework::HistogramRegistry* registry, + std::map> const& Specs, + int ChargeAbs, + int ChargeSign, + int PdgCodeAbs, + T const& ConfBinningQa) + { + this->template enableOptionalHistograms(ConfBinningQa); + this->template init(registry, Specs, ChargeAbs, ChargeSign, PdgCodeAbs); } template @@ -572,12 +598,11 @@ class TrackHistManager this->fillQa(track); } if constexpr (isFlagSet(mode, modes::Mode::kMc)) { - this->fillMc(track, mcParticles, mcMothers, mcPartonicMothers); + this->template fillMc(track, mcParticles, mcMothers, mcPartonicMothers); } } private: - // enable additional qa histograms template void enableOptionalHistograms(T const& ConfBinningQa) { @@ -590,20 +615,13 @@ class TrackHistManager mPlotTritonPid = ConfBinningQa.plotTritonPid.value; mPlotHeliumPid = ConfBinningQa.plotHeliumPid.value; mMomentumType = static_cast(ConfBinningQa.momentumType.value); - } - // enable additional qa and mc histograms - template - void enableOptionalHistograms(T1 const& ConfBinningQa, T2 const& ConfBinningMc) - { - this->template enableOptionalHistograms(ConfBinningQa); - - mPlotOrigins = ConfBinningMc.plotOrigins.value; - mPlotNSecondaries = ConfBinningMc.pdgCodesForMothersOfSecondary.value.size(); + mPlotOrigins = ConfBinningQa.plotOrigins.value; + mPlotNSecondaries = ConfBinningQa.pdgCodesForMothersOfSecondary.value.size(); for (std::size_t i = 0; i < MaxSecondary; i++) { - if (i < ConfBinningMc.pdgCodesForMothersOfSecondary.value.size()) { - mPdgCodesSecondaryMother.at(i) = std::abs(ConfBinningMc.pdgCodesForMothersOfSecondary.value.at(i)); + if (i < ConfBinningQa.pdgCodesForMothersOfSecondary.value.size()) { + mPdgCodesSecondaryMother.at(i) = std::abs(ConfBinningQa.pdgCodesForMothersOfSecondary.value.at(i)); } else { mPdgCodesSecondaryMother.at(i) = 0; } @@ -719,10 +737,20 @@ class TrackHistManager void initMc(std::map> const& Specs) { std::string mcDir = std::string(prefix) + std::string(McDir); - mHistogramRegistry->add(mcDir + getHistNameV2(kTruePt, HistTable), getHistDesc(kTruePt, HistTable), getHistType(kTruePt, HistTable), {Specs.at(kTruePt)}); - mHistogramRegistry->add(mcDir + getHistNameV2(kTrueEta, HistTable), getHistDesc(kTrueEta, HistTable), getHistType(kTrueEta, HistTable), {Specs.at(kTrueEta)}); - mHistogramRegistry->add(mcDir + getHistNameV2(kTruePhi, HistTable), getHistDesc(kTruePhi, HistTable), getHistType(kTruePhi, HistTable), {Specs.at(kTruePhi)}); - mHistogramRegistry->add(mcDir + getHistNameV2(kOrigin, HistTable), getHistDesc(kOrigin, HistTable), getHistType(kOrigin, HistTable), {Specs.at(kOrigin)}); + mHistogramRegistry->add(mcDir + getHistNameV2(kTruePtVsPt, HistTable), getHistDesc(kTruePtVsPt, HistTable), getHistType(kTruePtVsPt, HistTable), {Specs.at(kTruePtVsPt)}); + mHistogramRegistry->add(mcDir + getHistNameV2(kTrueEtaVsEta, HistTable), getHistDesc(kTrueEtaVsEta, HistTable), getHistType(kTrueEtaVsEta, HistTable), {Specs.at(kTrueEtaVsEta)}); + mHistogramRegistry->add(mcDir + getHistNameV2(kTruePhiVsPhi, HistTable), getHistDesc(kTruePhiVsPhi, HistTable), getHistType(kTruePhiVsPhi, HistTable), {Specs.at(kTruePhiVsPhi)}); + + // mc origin can be configured here + const framework::AxisSpec axisOrigin = {static_cast(modes::McOrigin::kMcOriginLast), -0.5, static_cast(modes::McOrigin::kMcOriginLast) - 0.5}; + mHistogramRegistry->add(mcDir + getHistNameV2(kOrigin, HistTable), getHistDesc(kOrigin, HistTable), getHistType(kOrigin, HistTable), {axisOrigin}); + mHistogramRegistry->get(HIST(prefix) + HIST(McDir) + HIST(histmanager::getHistName(kOrigin, HistTable)))->GetXaxis()->SetBinLabel(1 + static_cast(modes::McOrigin::kNoMcParticle), modes::mcOriginToString(modes::McOrigin::kNoMcParticle)); + mHistogramRegistry->get(HIST(prefix) + HIST(McDir) + HIST(histmanager::getHistName(kOrigin, HistTable)))->GetXaxis()->SetBinLabel(1 + static_cast(modes::McOrigin::kFromWrongCollision), modes::mcOriginToString(modes::McOrigin::kFromWrongCollision)); + mHistogramRegistry->get(HIST(prefix) + HIST(McDir) + HIST(histmanager::getHistName(kOrigin, HistTable)))->GetXaxis()->SetBinLabel(1 + static_cast(modes::McOrigin::kPhysicalPrimary), modes::mcOriginToString(modes::McOrigin::kPhysicalPrimary)); + mHistogramRegistry->get(HIST(prefix) + HIST(McDir) + HIST(histmanager::getHistName(kOrigin, HistTable)))->GetXaxis()->SetBinLabel(1 + static_cast(modes::McOrigin::kFromSecondaryDecay), modes::mcOriginToString(modes::McOrigin::kFromSecondaryDecay)); + mHistogramRegistry->get(HIST(prefix) + HIST(McDir) + HIST(histmanager::getHistName(kOrigin, HistTable)))->GetXaxis()->SetBinLabel(1 + static_cast(modes::McOrigin::kFromMaterial), modes::mcOriginToString(modes::McOrigin::kFromMaterial)); + mHistogramRegistry->get(HIST(prefix) + HIST(McDir) + HIST(histmanager::getHistName(kOrigin, HistTable)))->GetXaxis()->SetBinLabel(1 + static_cast(modes::McOrigin::kMissidentified), modes::mcOriginToString(modes::McOrigin::kMissidentified)); + mHistogramRegistry->add(mcDir + getHistNameV2(kPdg, HistTable), getHistDesc(kPdg, HistTable), getHistType(kPdg, HistTable), {Specs.at(kPdg)}); mHistogramRegistry->add(mcDir + getHistNameV2(kPdgMother, HistTable), getHistDesc(kPdgMother, HistTable), getHistType(kPdgMother, HistTable), {Specs.at(kPdgMother)}); mHistogramRegistry->add(mcDir + getHistNameV2(kPdgPartonicMother, HistTable), getHistDesc(kPdgPartonicMother, HistTable), getHistType(kPdgPartonicMother, HistTable), {Specs.at(kPdgPartonicMother)}); @@ -732,6 +760,7 @@ class TrackHistManager mHistogramRegistry->add(mcDir + getHistNameV2(kPrimary, HistTable), getHistDesc(kPrimary, HistTable), getHistType(kPrimary, HistTable), {Specs.at(kPrimary)}); mHistogramRegistry->add(mcDir + getHistNameV2(kFromWrongCollision, HistTable), getHistDesc(kFromWrongCollision, HistTable), getHistType(kFromWrongCollision, HistTable), {Specs.at(kFromWrongCollision)}); mHistogramRegistry->add(mcDir + getHistNameV2(kFromMaterial, HistTable), getHistDesc(kFromMaterial, HistTable), getHistType(kFromMaterial, HistTable), {Specs.at(kFromMaterial)}); + mHistogramRegistry->add(mcDir + getHistNameV2(kMissidentified, HistTable), getHistDesc(kMissidentified, HistTable), getHistType(kMissidentified, HistTable), {Specs.at(kMissidentified)}); if (mPlotNSecondaries >= histmanager::kSecondaryPlotLevel1) { mHistogramRegistry->add(mcDir + getHistNameV2(kSecondary1, HistTable), getHistDesc(kSecondary1, HistTable), getHistType(kSecondary1, HistTable), {Specs.at(kSecondary1)}); @@ -858,28 +887,40 @@ class TrackHistManager mHistogramRegistry->fill(HIST(prefix) + HIST(PidDir) + HIST(getHistName(kTpctofHelium, HistTable)), momentum, track.tpctofNSigmaHe()); } } - template + + template void fillMc(T1 const& track, T2 const& /*mcParticles*/, T3 const& /*mcMothers*/, T4 const& /*mcPartonicMothers*/) { // No MC Particle if (!track.has_fMcParticle()) { mHistogramRegistry->fill(HIST(prefix) + HIST(McDir) + HIST(getHistName(kPdg, HistTable)), 0); mHistogramRegistry->fill(HIST(prefix) + HIST(McDir) + HIST(getHistName(kOrigin, HistTable)), static_cast(modes::McOrigin::kNoMcParticle)); - if (mPlotOrigins) { - mHistogramRegistry->fill(HIST(prefix) + HIST(McDir) + HIST(getHistName(kNoMcParticle, HistTable)), track.pt(), track.dcaXY(), track.dcaZ()); + if constexpr (modes::isFlagSet(mode, modes::Mode::kQa)) { + if (mPlotOrigins) { + mHistogramRegistry->fill(HIST(prefix) + HIST(McDir) + HIST(getHistName(kNoMcParticle, HistTable)), track.pt(), track.dcaXY(), track.dcaZ()); + } } return; } // Retrieve MC particle auto mcParticle = track.template fMcParticle_as(); - mHistogramRegistry->fill(HIST(prefix) + HIST(McDir) + HIST(getHistName(kTruePt, HistTable)), mcParticle.pt()); - mHistogramRegistry->fill(HIST(prefix) + HIST(McDir) + HIST(getHistName(kTrueEta, HistTable)), mcParticle.eta()); - mHistogramRegistry->fill(HIST(prefix) + HIST(McDir) + HIST(getHistName(kTruePhi, HistTable)), mcParticle.phi()); - mHistogramRegistry->fill(HIST(prefix) + HIST(McDir) + HIST(getHistName(kOrigin, HistTable)), mcParticle.origin()); + + // missidentifed particles are special case + // whether a particle is missidentfied or not cannot be known by the producer so we check it here + bool isMissidentified = mcParticle.pdgCode() != mPdgCode; + + mHistogramRegistry->fill(HIST(prefix) + HIST(McDir) + HIST(getHistName(kTruePtVsPt, HistTable)), mcParticle.pt(), track.pt()); + mHistogramRegistry->fill(HIST(prefix) + HIST(McDir) + HIST(getHistName(kTrueEtaVsEta, HistTable)), mcParticle.eta(), track.eta()); + mHistogramRegistry->fill(HIST(prefix) + HIST(McDir) + HIST(getHistName(kTruePhiVsPhi, HistTable)), mcParticle.phi(), track.phi()); + if (isMissidentified) { + mHistogramRegistry->fill(HIST(prefix) + HIST(McDir) + HIST(getHistName(kOrigin, HistTable)), static_cast(modes::McOrigin::kMissidentified)); + } else { + mHistogramRegistry->fill(HIST(prefix) + HIST(McDir) + HIST(getHistName(kOrigin, HistTable)), mcParticle.origin()); + } mHistogramRegistry->fill(HIST(prefix) + HIST(McDir) + HIST(getHistName(kPdg, HistTable)), mcParticle.pdgCode()); - // Mother PDG + // get mother if (track.has_fMcMother()) { auto mother = track.template fMcMother_as(); mHistogramRegistry->fill(HIST(prefix) + HIST(McDir) + HIST(getHistName(kPdgMother, HistTable)), mother.pdgCode()); @@ -887,7 +928,7 @@ class TrackHistManager mHistogramRegistry->fill(HIST(prefix) + HIST(McDir) + HIST(getHistName(kPdgMother, HistTable)), 0); } - // Partonic Mother PDG + // get partonic mother if (track.has_fMcPartMoth()) { auto partonicMother = track.template fMcPartMoth_as(); mHistogramRegistry->fill(HIST(prefix) + HIST(McDir) + HIST(getHistName(kPdgPartonicMother, HistTable)), partonicMother.pdgCode()); @@ -895,43 +936,52 @@ class TrackHistManager mHistogramRegistry->fill(HIST(prefix) + HIST(McDir) + HIST(getHistName(kPdgPartonicMother, HistTable)), 0); } - // Plot origins - if (mPlotOrigins) { - switch (static_cast(mcParticle.origin())) { - case modes::McOrigin::kPhysicalPrimary: - mHistogramRegistry->fill(HIST(prefix) + HIST(McDir) + HIST(getHistName(kPrimary, HistTable)), track.pt(), track.dcaXY(), track.dcaZ()); - break; - case modes::McOrigin::kFromWrongCollision: - mHistogramRegistry->fill(HIST(prefix) + HIST(McDir) + HIST(getHistName(kFromWrongCollision, HistTable)), track.pt(), track.dcaXY(), track.dcaZ()); - break; - case modes::McOrigin::kFromMaterial: - mHistogramRegistry->fill(HIST(prefix) + HIST(McDir) + HIST(getHistName(kFromMaterial, HistTable)), track.pt(), track.dcaXY(), track.dcaZ()); - break; - case modes::McOrigin::kFromSecondaryDecay: - if (track.has_fMcMother()) { - auto mother = track.template fMcMother_as(); - int motherPdgCode = std::abs(mother.pdgCode()); - // Switch on PDG of the mother - if (mPlotNSecondaries >= histmanager::kSecondaryPlotLevel1 && motherPdgCode == mPdgCodesSecondaryMother[0]) { - mHistogramRegistry->fill(HIST(prefix) + HIST(McDir) + HIST(getHistName(kSecondary1, HistTable)), track.pt(), track.dcaXY(), track.dcaZ()); - } else if (mPlotNSecondaries >= histmanager::kSecondaryPlotLevel2 && motherPdgCode == mPdgCodesSecondaryMother[1]) { - mHistogramRegistry->fill(HIST(prefix) + HIST(McDir) + HIST(getHistName(kSecondary2, HistTable)), track.pt(), track.dcaXY(), track.dcaZ()); - } else if (mPlotNSecondaries >= histmanager::kSecondaryPlotLevel3 && motherPdgCode == mPdgCodesSecondaryMother[2]) { - mHistogramRegistry->fill(HIST(prefix) + HIST(McDir) + HIST(getHistName(kSecondary3, HistTable)), track.pt(), track.dcaXY(), track.dcaZ()); - } else { - mHistogramRegistry->fill(HIST(prefix) + HIST(McDir) + HIST(getHistName(kSecondaryOther, HistTable)), track.pt(), track.dcaXY(), track.dcaZ()); - } + if constexpr (modes::isFlagSet(mode, modes::Mode::kQa)) { + if (mPlotOrigins) { + // check first if particle is missidentified + if (isMissidentified) { + // if it is, we fill it as such + mHistogramRegistry->fill(HIST(prefix) + HIST(McDir) + HIST(getHistName(kMissidentified, HistTable)), track.pt(), track.dcaXY(), track.dcaZ()); + } else { + // if not, we fill it acccoridng to its origin + switch (static_cast(mcParticle.origin())) { + case modes::McOrigin::kPhysicalPrimary: + mHistogramRegistry->fill(HIST(prefix) + HIST(McDir) + HIST(getHistName(kPrimary, HistTable)), track.pt(), track.dcaXY(), track.dcaZ()); + break; + case modes::McOrigin::kFromWrongCollision: + mHistogramRegistry->fill(HIST(prefix) + HIST(McDir) + HIST(getHistName(kFromWrongCollision, HistTable)), track.pt(), track.dcaXY(), track.dcaZ()); + break; + case modes::McOrigin::kFromMaterial: + mHistogramRegistry->fill(HIST(prefix) + HIST(McDir) + HIST(getHistName(kFromMaterial, HistTable)), track.pt(), track.dcaXY(), track.dcaZ()); + break; + case modes::McOrigin::kFromSecondaryDecay: + if (track.has_fMcMother()) { + auto mother = track.template fMcMother_as(); + int motherPdgCode = std::abs(mother.pdgCode()); + // Switch on PDG of the mother + if (mPlotNSecondaries >= histmanager::kSecondaryPlotLevel1 && motherPdgCode == mPdgCodesSecondaryMother[0]) { + mHistogramRegistry->fill(HIST(prefix) + HIST(McDir) + HIST(getHistName(kSecondary1, HistTable)), track.pt(), track.dcaXY(), track.dcaZ()); + } else if (mPlotNSecondaries >= histmanager::kSecondaryPlotLevel2 && motherPdgCode == mPdgCodesSecondaryMother[1]) { + mHistogramRegistry->fill(HIST(prefix) + HIST(McDir) + HIST(getHistName(kSecondary2, HistTable)), track.pt(), track.dcaXY(), track.dcaZ()); + } else if (mPlotNSecondaries >= histmanager::kSecondaryPlotLevel3 && motherPdgCode == mPdgCodesSecondaryMother[2]) { + mHistogramRegistry->fill(HIST(prefix) + HIST(McDir) + HIST(getHistName(kSecondary3, HistTable)), track.pt(), track.dcaXY(), track.dcaZ()); + } else { + mHistogramRegistry->fill(HIST(prefix) + HIST(McDir) + HIST(getHistName(kSecondaryOther, HistTable)), track.pt(), track.dcaXY(), track.dcaZ()); + } + } + break; + default: + LOG(warn) << "Encounted partilce with unknown origin!"; + break; } - break; - default: - // Unknown origin → safely ignore - break; + } } } } o2::framework::HistogramRegistry* mHistogramRegistry = nullptr; int mAbsCharge = 1; + int mPdgCode = 0; bool mPlot2d = false; bool mPlotElectronPid = false; bool mPlotPionPid = false; diff --git a/PWGCF/Femto/Core/twoTrackResonanceBuilder.h b/PWGCF/Femto/Core/twoTrackResonanceBuilder.h index 14cb1b2c5a5..9183527b890 100644 --- a/PWGCF/Femto/Core/twoTrackResonanceBuilder.h +++ b/PWGCF/Femto/Core/twoTrackResonanceBuilder.h @@ -129,7 +129,7 @@ struct ConfKstar0Bits : o2::framework::ConfigurableGroup { #undef TWOTRACKRESONANCE_PIONPID_BITS #define TWOTRACKRESONANCE_DEFAULT_SELECTION(defaultPdgCode, defaultMassMin, defaultMassMax) \ - o2::framework::Configurable pdgCode{"pdgCode", defaultPdgCode, "Resonance PDG code"}; \ + o2::framework::Configurable pdgCodeAbs{"pdgCodeAbs", defaultPdgCode, "Resonance PDG code. Set sign to minus 1 for antiparticle"}; \ o2::framework::Configurable ptMin{"ptMin", 0.f, "Minimum pT"}; \ o2::framework::Configurable ptMax{"ptMax", 6.f, "Maximum pT"}; \ o2::framework::Configurable etaMin{"etaMin", -0.9f, "Minimum eta"}; \ @@ -148,11 +148,13 @@ struct ConfKstar0Bits : o2::framework::ConfigurableGroup { struct ConfPhiSelection : o2::framework::ConfigurableGroup { std::string prefix = std::string("PhiSelection"); TWOTRACKRESONANCE_DEFAULT_SELECTION(333, 0.95f, 1.05f) + o2::framework::Configurable sign{"sign", 1, "Dummy value for compatability"}; }; struct ConfRho0Selection : o2::framework::ConfigurableGroup { std::string prefix = std::string("Rho0Selection"); TWOTRACKRESONANCE_DEFAULT_SELECTION(113, 0.7f, 0.84f) + o2::framework::Configurable sign{"sign", 1, "Dummy value for compatability"}; }; struct ConfKstar0Selection : o2::framework::ConfigurableGroup { diff --git a/PWGCF/Femto/Core/twoTrackResonanceHistManager.h b/PWGCF/Femto/Core/twoTrackResonanceHistManager.h index e9275f15adf..848dd69f521 100644 --- a/PWGCF/Femto/Core/twoTrackResonanceHistManager.h +++ b/PWGCF/Femto/Core/twoTrackResonanceHistManager.h @@ -118,8 +118,6 @@ constexpr char PrefixKstar[] = "Kstar0/"; constexpr std::string_view AnalysisDir = "Kinematics/"; constexpr std::string_view QaDir = "QA/"; -constexpr int AbsChargeDaughers = 1; - template + // init analysis and mc + template void init(o2::framework::HistogramRegistry* registry, std::map> const& ResoSpecs, + T const& ConfResoSelection, std::map> const& PosDauSpecs, std::map> const& NegDauSpecs) { mHistogramRegistry = registry; - mPosDauManager.template init(registry, PosDauSpecs, AbsChargeDaughers); - mNegDauManager.template init(registry, NegDauSpecs, AbsChargeDaughers); + mPdgCode = std::abs(ConfResoSelection.pdgCodeAbs.value); + + int posDauPdgCodeAbs = 0; + int negDauPdgCodeAbs = 0; + const int absCharge = 1; + const int signPlus = 1; + const int signMinus = -1; + + if (mPdgCode == PDG_t::kRho770_0) { + posDauPdgCodeAbs = std::abs(PDG_t::kPiPlus); + negDauPdgCodeAbs = std::abs(PDG_t::kPiMinus); + } else if (mPdgCode == o2::constants::physics::Pdg::kPhi) { + posDauPdgCodeAbs = std::abs(PDG_t::kKPlus); + negDauPdgCodeAbs = std::abs(PDG_t::kKMinus); + } else if (mPdgCode == o2::constants::physics::Pdg::kK0Star892) { + if (ConfResoSelection.sign.value > 0) { + posDauPdgCodeAbs = std::abs(PDG_t::kKPlus); + negDauPdgCodeAbs = std::abs(PDG_t::kPiMinus); + } else { + mPdgCode = -1 * mPdgCode; + posDauPdgCodeAbs = std::abs(PDG_t::kKMinus); + negDauPdgCodeAbs = std::abs(PDG_t::kPiPlus); + } + } else { + LOG(fatal) << "PDG code for V0 has to be either Lambda or K0short"; + } + + mPosDauManager.template init(registry, PosDauSpecs, absCharge, signPlus, posDauPdgCodeAbs); + mNegDauManager.template init(registry, NegDauSpecs, absCharge, signMinus, negDauPdgCodeAbs); if constexpr (modes::isFlagSet(mode, modes::Mode::kAnalysis)) { initAnalysis(ResoSpecs); } @@ -147,17 +174,46 @@ class TwoTrackResonanceHistManager } } - template + // init analysis and qa and mc + template void init(o2::framework::HistogramRegistry* registry, - std::map> ResoSpecs, - std::map> PosDauSpecs, - T1 const& PosDauConfBinningQa, - std::map> NegDauSpecs, - T2 const& NegDauConfBinningQa) + std::map> const& ResoSpecs, + T1 const& ConfResoSelection, + std::map> const& PosDauSpecs, + T2 const& ConfPosDauBinningQa, + std::map> const& NegDauSpecs, + T3 const& ConfNegDauBinningQa) { mHistogramRegistry = registry; - mPosDauManager.template init(registry, PosDauSpecs, PosDauConfBinningQa, AbsChargeDaughers); - mNegDauManager.template init(registry, NegDauSpecs, NegDauConfBinningQa, AbsChargeDaughers); + mPdgCode = std::abs(ConfResoSelection.pdgCodeAbs.value); + + int posDauPdgCodeAbs = 0; + int negDauPdgCodeAbs = 0; + const int absCharge = 1; + const int signPlus = 1; + const int signMinus = -1; + + if (mPdgCode == PDG_t::kRho770_0) { + posDauPdgCodeAbs = std::abs(PDG_t::kPiPlus); + negDauPdgCodeAbs = std::abs(PDG_t::kPiMinus); + } else if (mPdgCode == o2::constants::physics::Pdg::kPhi) { + posDauPdgCodeAbs = std::abs(PDG_t::kKPlus); + negDauPdgCodeAbs = std::abs(PDG_t::kKMinus); + } else if (mPdgCode == o2::constants::physics::Pdg::kK0Star892) { + if (ConfResoSelection.sign.value > 0) { + posDauPdgCodeAbs = std::abs(PDG_t::kKPlus); + negDauPdgCodeAbs = std::abs(PDG_t::kPiMinus); + } else { + mPdgCode = -1 * mPdgCode; + posDauPdgCodeAbs = std::abs(PDG_t::kKMinus); + negDauPdgCodeAbs = std::abs(PDG_t::kPiPlus); + } + } else { + LOG(fatal) << "PDG code for TwoTrackResonance has to be either Rho, Phi or K^0*(892)"; + } + + mPosDauManager.template init(registry, PosDauSpecs, absCharge, signPlus, posDauPdgCodeAbs, ConfPosDauBinningQa); + mNegDauManager.template init(registry, NegDauSpecs, absCharge, signMinus, negDauPdgCodeAbs, ConfNegDauBinningQa); if constexpr (modes::isFlagSet(mode, modes::Mode::kAnalysis)) { initAnalysis(ResoSpecs); } @@ -230,6 +286,7 @@ class TwoTrackResonanceHistManager o2::framework::HistogramRegistry* mHistogramRegistry = nullptr; trackhistmanager::TrackHistManager mPosDauManager; trackhistmanager::TrackHistManager mNegDauManager; + int mPdgCode = 0; }; }; // namespace twotrackresonancehistmanager }; // namespace o2::analysis::femto diff --git a/PWGCF/Femto/Core/v0Builder.h b/PWGCF/Femto/Core/v0Builder.h index b70182cac04..e7907f56f41 100644 --- a/PWGCF/Femto/Core/v0Builder.h +++ b/PWGCF/Femto/Core/v0Builder.h @@ -93,16 +93,16 @@ struct ConfK0shortBits : o2::framework::ConfigurableGroup { #undef V0_DEFAULT_BITS // base selection for analysis task for v0s -#define V0_DEFAULT_SELECTIONS(defaultMassMin, defaultMassMax, defaultPdgCode) \ - o2::framework::Configurable pdgCode{"pdgCode", defaultPdgCode, "V0 PDG code"}; \ - o2::framework::Configurable ptMin{"ptMin", 0.f, "Minimum pT"}; \ - o2::framework::Configurable ptMax{"ptMax", 999.f, "Maximum pT"}; \ - o2::framework::Configurable etaMin{"etaMin", -10.f, "Minimum eta"}; \ - o2::framework::Configurable etaMax{"etaMax", 10.f, "Maximum eta"}; \ - o2::framework::Configurable phiMin{"phiMin", 0.f, "Minimum eta"}; \ - o2::framework::Configurable phiMax{"phiMax", 1.f * o2::constants::math::TwoPI, "Maximum phi"}; \ - o2::framework::Configurable massMin{"massMin", defaultMassMin, "Minimum invariant mass for Lambda"}; \ - o2::framework::Configurable massMax{"massMax", defaultMassMax, "Maximum invariant mass for Lambda"}; \ +#define V0_DEFAULT_SELECTIONS(defaultMassMin, defaultMassMax, defaultPdgCode) \ + o2::framework::Configurable pdgCodeAbs{"pdgCodeAbs", defaultPdgCode, "PDG code. Set sign to -1 for antiparticle"}; \ + o2::framework::Configurable ptMin{"ptMin", 0.f, "Minimum pT"}; \ + o2::framework::Configurable ptMax{"ptMax", 999.f, "Maximum pT"}; \ + o2::framework::Configurable etaMin{"etaMin", -10.f, "Minimum eta"}; \ + o2::framework::Configurable etaMax{"etaMax", 10.f, "Maximum eta"}; \ + o2::framework::Configurable phiMin{"phiMin", 0.f, "Minimum eta"}; \ + o2::framework::Configurable phiMax{"phiMax", 1.f * o2::constants::math::TwoPI, "Maximum phi"}; \ + o2::framework::Configurable massMin{"massMin", defaultMassMin, "Minimum invariant mass for Lambda"}; \ + o2::framework::Configurable massMax{"massMax", defaultMassMax, "Maximum invariant mass for Lambda"}; \ o2::framework::Configurable mask{"mask", 0, "Bitmask for v0 selection"}; // base selection for analysis task for lambdas @@ -118,6 +118,7 @@ template struct ConfK0shortSelection : o2::framework::ConfigurableGroup { std::string prefix = Prefix; V0_DEFAULT_SELECTIONS(0.47, 0.51, 310) + o2::framework::Configurable sign{"sign", 0, "Dummy value. For compatability with Lambda selection"}; }; #undef V0_DEFAULT_SELECTIONS @@ -142,8 +143,8 @@ enum V0Sels { // selection for daughter kDauAbsEtaMax, ///< Max. absolute pseudo rapidity - kDauDcaMin, ///< Min. DCA of the positive daughters at primary vertex - kDauTpcClsMin, ///< Min. number of TPC clusters of positive daughter + kDauDcaMin, ///< Min. DCA of the daughters at primary vertex + kDauTpcClsMin, ///< Min. number of TPC clusters of daughter // pid selection for daughters kPosDaughTpcPion, ///< TPC Pion PID for positive daughter @@ -165,9 +166,9 @@ const std::unordered_map v0SelectionNames = { {kTransRadMin, "Min. transverse radius"}, {kTransRadMax, "Max. transverse radius"}, - {kDauAbsEtaMax, "Max. absolute pseudo rapidity"}, - {kDauDcaMin, "Min. DCA of the positive daughters at primary vertex"}, - {kDauTpcClsMin, "Min. number of TPC clusters of positive daughter"}, + {kDauAbsEtaMax, "Max. absolute pseudo rapidity of daughters"}, + {kDauDcaMin, "Min. DCA of the daughters at primary vertex"}, + {kDauTpcClsMin, "Min. number of TPC clusters of daughters"}, {kPosDaughTpcPion, "TPC Pion PID for positive daughter"}, {kPosDaughTpcProton, "TPC Proton PID for positive daughter"}, diff --git a/PWGCF/Femto/Core/v0HistManager.h b/PWGCF/Femto/Core/v0HistManager.h index 6e2512da15a..b79972437ef 100644 --- a/PWGCF/Femto/Core/v0HistManager.h +++ b/PWGCF/Femto/Core/v0HistManager.h @@ -70,14 +70,15 @@ enum V0Hist { kPdg, kPdgMother, kPdgPartonicMother, - kTruePt, - kTrueEta, - kTruePhi, + kTruePtVsPt, + kTrueEtaVsEta, + kTruePhiVsPhi, // histograms for fraction estimation of v0s kNoMcParticle, kPrimary, kFromWrongCollision, kFromMaterial, + kMissidentified, kSecondary1, kSecondary2, kSecondary3, @@ -93,7 +94,8 @@ constexpr std::size_t MaxSecondary = 3; o2::framework::ConfigurableAxis eta{"eta", {{300, -1.5, 1.5}}, "Eta"}; \ o2::framework::ConfigurableAxis phi{"phi", {{720, 0, 1.f * o2::constants::math::TwoPI}}, "Phi"}; \ o2::framework::ConfigurableAxis mass{"mass", {{200, defaultMassMin, defaultMassMax}}, "Mass"}; \ - o2::framework::ConfigurableAxis sign{"sign", {{3, -1.5, 1.5}}, "Sign"}; + o2::framework::ConfigurableAxis sign{"sign", {{3, -1.5, 1.5}}, "Sign"}; \ + o2::framework::ConfigurableAxis pdgCodes{"pdgCodes", {{8001, -4000.5, 4000.5}}, "MC ONLY: PDG codes of reconstructed V0s"}; template struct ConfLambdaBinning : o2::framework::ConfigurableGroup { @@ -116,6 +118,8 @@ template struct ConfV0QaBinning : o2::framework::ConfigurableGroup { std::string prefix = Prefix; o2::framework::Configurable plot2d{"plot2d", true, "Generate various 2D QA plots"}; + o2::framework::Configurable plotOrigins{"plotOrigins", true, "MC ONLY: Plot pt vs cosPa for different particle origins"}; + o2::framework::Configurable> pdgCodesForMothersOfSecondary{"pdgCodesForMothersOfSecondary", {3312, 3334}, "MC ONLY: PDG codes of mothers of secondaries (Max 3 will be considered)"}; o2::framework::ConfigurableAxis cosPa{"cosPa", {{100, 0.9, 1}}, "Cosine of poiting angle"}; o2::framework::ConfigurableAxis dauDcaAtDecay{"dauDcaAtDecay", {{150, 0, 1.5}}, "Daughter DCA at decay vertex"}; o2::framework::ConfigurableAxis decayVertex{"decayVertex", {{100, 0, 100}}, "Decay vertex"}; @@ -131,23 +135,6 @@ using ConfLambdaQaBinning1 = ConfV0QaBinning; constexpr const char PrefixK0shortQaBinning1[] = "K0shortQaBinning1"; using ConfK0shortQaBinning1 = ConfV0QaBinning; -template -struct ConfV0McBinning : o2::framework::ConfigurableGroup { - std::string prefix = std::string(Prefix); - o2::framework::ConfigurableAxis pdgCodes{"pdgCodes", {{8001, -4000.5, 4000.5}}, "PDG codes of selected V0s"}; - o2::framework::ConfigurableAxis statusCode{"statusCode", {{21, -0.5, 20.5}}, "Status codes (i.e. Origin)"}; - o2::framework::Configurable plotOrigins{"plotOrigins", true, "Plot pt vs cosPa for different particle origins"}; - o2::framework::Configurable> pdgCodesForMothersOfSecondary{"pdgCodesForMothersOfSecondary", {3312, 3334}, "PDG codes of mothers of secondaries (Max 3 will be considered)"}; - o2::framework::ConfigurableAxis pt{"pt", {{150, 0, 3}}, "Pt"}; - o2::framework::ConfigurableAxis cosPa{"cosPa", {{100, 0.8, 1}}, "Cosine pointing angle"}; -}; - -constexpr const char PrefixLambdaMcBinning1[] = "LambdaMcBinning1"; -using ConfLambdaMcBinning1 = ConfV0McBinning; - -constexpr const char PrefixK0shortMcBinning1[] = "K0shortMcBinning1"; -using ConfK0shortMcBinning1 = ConfV0McBinning; - // must be in sync with enum TrackVariables // the enum gives the correct index in the array constexpr std::array, kV0HistLast> HistTable = { @@ -180,13 +167,14 @@ constexpr std::array, kV0HistLast> HistTable = { {kPdg, o2::framework::kTH1F, "hPdg", "PDG Codes of reconstructed v0; PDG Code; Entries"}, {kPdgMother, o2::framework::kTH1F, "hPdgMother", "PDG Codes of mother of reconstructed v0; PDG Code; Entries"}, {kPdgPartonicMother, o2::framework::kTH1F, "hPdgPartonicMother", "PDG Codes of partonic mother of reconstructed v0; PDG Code; Entries"}, - {kTruePt, o2::framework::kTH1F, "hTruePt", "True transverse momentum; p_{T} (GeV/#it{c}); Entries"}, - {kTrueEta, o2::framework::kTH1F, "hTrueEta", "True pseudorapdity; #eta; Entries"}, - {kTruePhi, o2::framework::kTH1F, "hTruePhi", "True azimuthal angle; #varphi; Entries"}, + {kTruePtVsPt, o2::framework::kTH2F, "hTruePtVsPt", "True transverse momentum vs transverse momentum; p_{T,True} (GeV/#it{c}); p_{T,True} (GeV/#it{c})"}, + {kTrueEtaVsEta, o2::framework::kTH2F, "hTrueEtaVsEta", "True pseudorapdity vs pseudorapdity; #eta_{True}; #eta"}, + {kTruePhiVsPhi, o2::framework::kTH2F, "hTruePhiVsPhi", "True azimuthal angle vs azimuthal angle; #varphi_{True}; #varphi"}, {kNoMcParticle, o2::framework::kTH2F, "hNoMcParticle", "Wrongly reconstructed particles; p_{T} (GeV/#it{c}); cos(#alpha)"}, {kPrimary, o2::framework::kTH2F, "hPrimary", "Primary particles; p_{T} (GeV/#it{c}); cos(#alpha)"}, {kFromWrongCollision, o2::framework::kTH2F, "hFromWrongCollision", "Particles associated to wrong collision; p_{T} (GeV/#it{c}); cos(#alpha)"}, {kFromMaterial, o2::framework::kTH2F, "hFromMaterial", "Particles from material; p_{T} (GeV/#it{c}); cos(#alpha)"}, + {kMissidentified, o2::framework::kTH2F, "hMissidentified", "Missidentified particles (fake/wrong PDG code); p_{T} (GeV/#it{c}); cos(#alpha)"}, {kSecondary1, o2::framework::kTH2F, "hFromSecondary1", "Particles from secondary decay; p_{T} (GeV/#it{c}); cos(#alpha)"}, {kSecondary2, o2::framework::kTH2F, "hFromSecondary2", "Particles from seconary decay; p_{T} (GeV/#it{c}); cos(#alpha)"}, {kSecondary3, o2::framework::kTH2F, "hFromSecondary3", "Particles from seconary decay; p_{T} (GeV/#it{c}); cos(#alpha)"}, @@ -200,6 +188,14 @@ constexpr std::array, kV0HistLast> HistTable = { {kMass, {conf.mass}}, \ {kSign, {conf.sign}}, +#define V0_HIST_MC_MAP(conf) \ + {kTruePtVsPt, {conf.pt, conf.pt}}, \ + {kTrueEtaVsEta, {conf.eta, conf.eta}}, \ + {kTruePhiVsPhi, {conf.phi, conf.phi}}, \ + {kPdg, {conf.pdgCodes}}, \ + {kPdgMother, {conf.pdgCodes}}, \ + {kPdgPartonicMother, {conf.pdgCodes}}, + #define V0_HIST_QA_MAP(confAnalysis, confQa) \ {kCosPa, {confQa.cosPa}}, \ {kDecayDauDca, {confQa.dauDcaAtDecay}}, \ @@ -222,22 +218,16 @@ constexpr std::array, kV0HistLast> HistTable = { {kK0shortMassVsLambdaMass, {confQa.massK0short, confQa.massLambda}}, \ {kK0shortMassVsAntiLambdaMass, {confQa.massK0short, confQa.massAntiLambda}}, -#define V0_HIST_MC_MAP(confAnalysis, confMc) \ - {kTruePt, {confAnalysis.pt}}, \ - {kTrueEta, {confAnalysis.eta}}, \ - {kTruePhi, {confAnalysis.phi}}, \ - {kOrigin, {confMc.statusCode}}, \ - {kPdg, {confMc.pdgCodes}}, \ - {kPdgMother, {confMc.pdgCodes}}, \ - {kPdgPartonicMother, {confMc.pdgCodes}}, \ - {kNoMcParticle, {confMc.pt, confMc.cosPa}}, \ - {kPrimary, {confMc.pt, confMc.cosPa}}, \ - {kFromWrongCollision, {confMc.pt, confMc.cosPa}}, \ - {kFromMaterial, {confMc.pt, confMc.cosPa}}, \ - {kSecondary1, {confMc.pt, confMc.cosPa}}, \ - {kSecondary2, {confMc.pt, confMc.cosPa}}, \ - {kSecondary3, {confMc.pt, confMc.cosPa}}, \ - {kSecondaryOther, {confMc.pt, confMc.cosPa}}, +#define V0_HIST_MC_QA_MAP(confAnalysis, confQa) \ + {kNoMcParticle, {confAnalysis.pt, confQa.cosPa}}, \ + {kPrimary, {confAnalysis.pt, confQa.cosPa}}, \ + {kFromWrongCollision, {confAnalysis.pt, confQa.cosPa}}, \ + {kFromMaterial, {confAnalysis.pt, confQa.cosPa}}, \ + {kMissidentified, {confAnalysis.pt, confQa.cosPa}}, \ + {kSecondary1, {confAnalysis.pt, confQa.cosPa}}, \ + {kSecondary2, {confAnalysis.pt, confQa.cosPa}}, \ + {kSecondary3, {confAnalysis.pt, confQa.cosPa}}, \ + {kSecondaryOther, {confAnalysis.pt, confQa.cosPa}}, template auto makeV0HistSpecMap(const T& confBinningAnalysis) @@ -246,6 +236,14 @@ auto makeV0HistSpecMap(const T& confBinningAnalysis) V0_HIST_ANALYSIS_MAP(confBinningAnalysis)}; } +template +auto makeV0McHistSpecMap(const T& confBinningAnalysis) +{ + return std::map>{ + V0_HIST_ANALYSIS_MAP(confBinningAnalysis) + V0_HIST_MC_MAP(confBinningAnalysis)}; +} + template std::map> makeV0QaHistSpecMap(T1 const& confBinningAnalysis, T2 const& confBinningQa) { @@ -254,18 +252,20 @@ std::map> makeV0QaHistSpecMap(T1 const& V0_HIST_QA_MAP(confBinningAnalysis, confBinningQa)}; } -template -std::map> makeV0McQaHistSpecMap(T1 const& confBinningAnalysis, T2 const& confBinningQa, T3 const& confBinningMc) +template +std::map> makeV0McQaHistSpecMap(T1 const& confBinningAnalysis, T2 const& confBinningQa) { return std::map>{ V0_HIST_ANALYSIS_MAP(confBinningAnalysis) V0_HIST_QA_MAP(confBinningAnalysis, confBinningQa) - V0_HIST_MC_MAP(confBinningAnalysis, confBinningMc)}; + V0_HIST_MC_MAP(confBinningAnalysis) + V0_HIST_MC_QA_MAP(confBinningAnalysis, confBinningQa)}; } #undef V0_HIST_ANALYSIS_MAP #undef V0_HIST_QA_MAP #undef V0_HIST_MC_MAP +#undef V0_HIST_MC_QA_MAP constexpr char PrefixLambdaQa[] = "LambdaQA/"; constexpr char PrefixLambda1[] = "Lambda1/"; @@ -280,8 +280,6 @@ constexpr std::string_view AnalysisDir = "Kinematics/"; constexpr std::string_view QaDir = "QA/"; constexpr std::string_view McDir = "MC/"; -constexpr int AbsChargeDaughters = 1; - /// \class FemtoDreamEventHisto /// \brief Class for histogramming event properties // template @@ -295,80 +293,98 @@ class V0HistManager V0HistManager() = default; ~V0HistManager() = default; - // init for analysis - template + // init for analysis and mc + template void init(o2::framework::HistogramRegistry* registry, std::map> const& V0Specs, + T const& ConfV0Selection, std::map> const& PosDauSpecs, std::map> const& NegDauSpecs) { mHistogramRegistry = registry; - mPosDauManager.template init(registry, PosDauSpecs, AbsChargeDaughters); - mNegDauManager.template init(registry, NegDauSpecs, AbsChargeDaughters); - - if constexpr (modes::isFlagSet(mode, modes::Mode::kAnalysis)) { - initAnalysis(V0Specs); - } - if constexpr (modes::isFlagSet(mode, modes::Mode::kQa)) { - initQa(V0Specs); - } - if constexpr (modes::isFlagSet(mode, modes::Mode::kMc)) { - initMc(V0Specs); + mPdgCode = std::abs(ConfV0Selection.pdgCodeAbs.value); + + int posDauPdgCodeAbs = 0; + int negDauPdgCodeAbs = 0; + const int absCharge = 1; + const int signPlus = 1; + const int signMinus = -1; + + if (mPdgCode == PDG_t::kK0Short) { + posDauPdgCodeAbs = std::abs(PDG_t::kPiPlus); + negDauPdgCodeAbs = std::abs(PDG_t::kPiMinus); + } else if (mPdgCode == PDG_t::kLambda0) { + if (ConfV0Selection.sign.value > 0) { + posDauPdgCodeAbs = std::abs(PDG_t::kProton); + negDauPdgCodeAbs = std::abs(PDG_t::kPiMinus); + } else { + mPdgCode = -1 * mPdgCode; // switch sign for antilambda + posDauPdgCodeAbs = std::abs(PDG_t::kPiPlus); + negDauPdgCodeAbs = std::abs(PDG_t::kProtonBar); + } + } else { + LOG(fatal) << "PDG code for V0 has to be either Lambda or K0short"; } - } - // init for qa - template - void init(o2::framework::HistogramRegistry* registry, - std::map> const& V0Specs, - T1 const& V0ConfBinningQa, - std::map> const& PosDauSpecs, - T2 const& PosDauConfBinningQa, - std::map> const& NegDauSpecs, - T3 const& NegDauConfBinningQa) - { - mHistogramRegistry = registry; - mPosDauManager.template init(registry, PosDauSpecs, PosDauConfBinningQa, AbsChargeDaughters); - mNegDauManager.template init(registry, NegDauSpecs, NegDauConfBinningQa, AbsChargeDaughters); - this->enableOptionalHistograms(V0ConfBinningQa); + mPosDauManager.template init(registry, PosDauSpecs, absCharge, signPlus, posDauPdgCodeAbs); + mNegDauManager.template init(registry, NegDauSpecs, absCharge, signMinus, negDauPdgCodeAbs); if constexpr (modes::isFlagSet(mode, modes::Mode::kAnalysis)) { - initAnalysis(V0Specs); - } - if constexpr (modes::isFlagSet(mode, modes::Mode::kQa)) { - initQa(V0Specs); + this->initAnalysis(V0Specs); } if constexpr (modes::isFlagSet(mode, modes::Mode::kMc)) { - initMc(V0Specs); + this->initMc(V0Specs); } } - // init for qa + mc - template + // init for analysis and qa and and mc + template void init(o2::framework::HistogramRegistry* registry, std::map> const& V0Specs, - T1 const& V0ConfBinningQa, - T2 const& V0ConfBinningMc, + T1 const& ConfV0Selection, + T2 const& ConfV0BinningQa, std::map> const& PosDauSpecs, - T3 const& PosDauConfBinningQa, - T4 const& PosDauConfBinningMc, + T3 const& ConfPosDauBinningQa, std::map> const& NegDauSpecs, - T5 const& NegDauConfBinningQa, - T6 const& NegDauConfBinningMc) + T4 const& ConfNegDauBinningQa) { mHistogramRegistry = registry; - mPosDauManager.template init(registry, PosDauSpecs, PosDauConfBinningQa, PosDauConfBinningMc, AbsChargeDaughters); - mNegDauManager.template init(registry, NegDauSpecs, NegDauConfBinningQa, NegDauConfBinningMc, AbsChargeDaughters); - this->enableOptionalHistograms(V0ConfBinningQa, V0ConfBinningMc); + mPdgCode = std::abs(ConfV0Selection.pdgCodeAbs.value); + this->enableOptionalHistograms(ConfV0BinningQa); + + int posDauPdgCode = 0; + int negDauPdgCode = 0; + const int absCharge = 1; + const int signPlus = 1; + const int signMinus = -1; + + if (mPdgCode == PDG_t::kK0Short) { + posDauPdgCode = std::abs(PDG_t::kPiPlus); + negDauPdgCode = std::abs(PDG_t::kPiMinus); + } else if (mPdgCode == PDG_t::kLambda0) { + if (ConfV0Selection.sign.value > 0) { + posDauPdgCode = std::abs(PDG_t::kProton); + negDauPdgCode = std::abs(PDG_t::kPiMinus); + } else { + mPdgCode = -1 * mPdgCode; // set pdg code for antilambda + posDauPdgCode = std::abs(PDG_t::kPiPlus); + negDauPdgCode = std::abs(PDG_t::kProtonBar); + } + } else { + LOG(fatal) << "PDG code for V0 has to be either Lambda or K0short"; + } + + mPosDauManager.template init(registry, PosDauSpecs, absCharge, signPlus, posDauPdgCode, ConfPosDauBinningQa); + mNegDauManager.template init(registry, NegDauSpecs, absCharge, signMinus, negDauPdgCode, ConfNegDauBinningQa); if constexpr (modes::isFlagSet(mode, modes::Mode::kAnalysis)) { - initAnalysis(V0Specs); + this->initAnalysis(V0Specs); } if constexpr (modes::isFlagSet(mode, modes::Mode::kQa)) { - initQa(V0Specs); + this->initQa(V0Specs); } if constexpr (modes::isFlagSet(mode, modes::Mode::kMc)) { - initMc(V0Specs); + this->initMc(V0Specs); } } @@ -381,10 +397,10 @@ class V0HistManager mNegDauManager.template fill(negDaughter, tracks); if constexpr (modes::isFlagSet(mode, modes::Mode::kAnalysis)) { - fillAnalysis(v0candidate); + this->fillAnalysis(v0candidate); } if constexpr (modes::isFlagSet(mode, modes::Mode::kQa)) { - fillQa(v0candidate); + this->fillQa(v0candidate); } } @@ -396,36 +412,26 @@ class V0HistManager auto negDaughter = tracks.rawIteratorAt(v0candidate.negDauId() - tracks.offset()); mNegDauManager.template fill(negDaughter, tracks, mcParticles, mcMothers, mcPartonicMothers); if constexpr (modes::isFlagSet(mode, modes::Mode::kAnalysis)) { - fillAnalysis(v0candidate); + this->fillAnalysis(v0candidate); } if constexpr (modes::isFlagSet(mode, modes::Mode::kQa)) { - fillQa(v0candidate); + this->fillQa(v0candidate); } if constexpr (modes::isFlagSet(mode, modes::Mode::kMc)) { - fillMc(v0candidate, mcParticles, mcMothers, mcPartonicMothers); + this->template fillMc(v0candidate, mcParticles, mcMothers, mcPartonicMothers); } } private: - // for qa template void enableOptionalHistograms(T1 const& V0ConfBinningQa) { mPlot2d = V0ConfBinningQa.plot2d.value; - } - - // for qa and mc - template - void enableOptionalHistograms(T1 const& V0ConfBinningQa, T2 const& V0ConfBinningMc) - { - this->enableOptionalHistograms(V0ConfBinningQa); - - mPlotOrigins = V0ConfBinningMc.plotOrigins.value; - mPlotNSecondaries = V0ConfBinningMc.pdgCodesForMothersOfSecondary.value.size(); - + mPlotOrigins = V0ConfBinningQa.plotOrigins.value; + mPlotNSecondaries = V0ConfBinningQa.pdgCodesForMothersOfSecondary.value.size(); for (std::size_t i = 0; i < MaxSecondary; i++) { - if (i < V0ConfBinningMc.pdgCodesForMothersOfSecondary.value.size()) { - mPdgCodesSecondaryMother.at(i) = std::abs(V0ConfBinningMc.pdgCodesForMothersOfSecondary.value.at(i)); + if (i < V0ConfBinningQa.pdgCodesForMothersOfSecondary.value.size()) { + mPdgCodesSecondaryMother.at(i) = std::abs(V0ConfBinningQa.pdgCodesForMothersOfSecondary.value.at(i)); } else { mPdgCodesSecondaryMother.at(i) = 0; } @@ -475,10 +481,20 @@ class V0HistManager void initMc(std::map> const& V0Specs) { std::string mcDir = std::string(v0Prefix) + std::string(McDir); - mHistogramRegistry->add(mcDir + getHistNameV2(kTruePt, HistTable), getHistDesc(kTruePt, HistTable), getHistType(kTruePt, HistTable), {V0Specs.at(kTruePt)}); - mHistogramRegistry->add(mcDir + getHistNameV2(kTrueEta, HistTable), getHistDesc(kTrueEta, HistTable), getHistType(kTrueEta, HistTable), {V0Specs.at(kTrueEta)}); - mHistogramRegistry->add(mcDir + getHistNameV2(kTruePhi, HistTable), getHistDesc(kTruePhi, HistTable), getHistType(kTruePhi, HistTable), {V0Specs.at(kTruePhi)}); - mHistogramRegistry->add(mcDir + getHistNameV2(kOrigin, HistTable), getHistDesc(kOrigin, HistTable), getHistType(kOrigin, HistTable), {V0Specs.at(kOrigin)}); + mHistogramRegistry->add(mcDir + getHistNameV2(kTruePtVsPt, HistTable), getHistDesc(kTruePtVsPt, HistTable), getHistType(kTruePtVsPt, HistTable), {V0Specs.at(kTruePtVsPt)}); + mHistogramRegistry->add(mcDir + getHistNameV2(kTrueEtaVsEta, HistTable), getHistDesc(kTrueEtaVsEta, HistTable), getHistType(kTrueEtaVsEta, HistTable), {V0Specs.at(kTrueEtaVsEta)}); + mHistogramRegistry->add(mcDir + getHistNameV2(kTruePhiVsPhi, HistTable), getHistDesc(kTruePhiVsPhi, HistTable), getHistType(kTruePhiVsPhi, HistTable), {V0Specs.at(kTruePhiVsPhi)}); + + // mc origin can be configured here + const framework::AxisSpec axisOrigin = {static_cast(modes::McOrigin::kMcOriginLast), -0.5, static_cast(modes::McOrigin::kMcOriginLast) - 0.5}; + mHistogramRegistry->add(mcDir + getHistNameV2(kOrigin, HistTable), getHistDesc(kOrigin, HistTable), getHistType(kOrigin, HistTable), {axisOrigin}); + mHistogramRegistry->get(HIST(v0Prefix) + HIST(McDir) + HIST(histmanager::getHistName(kOrigin, HistTable)))->GetXaxis()->SetBinLabel(1 + static_cast(modes::McOrigin::kNoMcParticle), modes::mcOriginToString(modes::McOrigin::kNoMcParticle)); + mHistogramRegistry->get(HIST(v0Prefix) + HIST(McDir) + HIST(histmanager::getHistName(kOrigin, HistTable)))->GetXaxis()->SetBinLabel(1 + static_cast(modes::McOrigin::kFromWrongCollision), modes::mcOriginToString(modes::McOrigin::kFromWrongCollision)); + mHistogramRegistry->get(HIST(v0Prefix) + HIST(McDir) + HIST(histmanager::getHistName(kOrigin, HistTable)))->GetXaxis()->SetBinLabel(1 + static_cast(modes::McOrigin::kPhysicalPrimary), modes::mcOriginToString(modes::McOrigin::kPhysicalPrimary)); + mHistogramRegistry->get(HIST(v0Prefix) + HIST(McDir) + HIST(histmanager::getHistName(kOrigin, HistTable)))->GetXaxis()->SetBinLabel(1 + static_cast(modes::McOrigin::kFromSecondaryDecay), modes::mcOriginToString(modes::McOrigin::kFromSecondaryDecay)); + mHistogramRegistry->get(HIST(v0Prefix) + HIST(McDir) + HIST(histmanager::getHistName(kOrigin, HistTable)))->GetXaxis()->SetBinLabel(1 + static_cast(modes::McOrigin::kFromMaterial), modes::mcOriginToString(modes::McOrigin::kFromMaterial)); + mHistogramRegistry->get(HIST(v0Prefix) + HIST(McDir) + HIST(histmanager::getHistName(kOrigin, HistTable)))->GetXaxis()->SetBinLabel(1 + static_cast(modes::McOrigin::kMissidentified), modes::mcOriginToString(modes::McOrigin::kMissidentified)); + mHistogramRegistry->add(mcDir + getHistNameV2(kPdg, HistTable), getHistDesc(kPdg, HistTable), getHistType(kPdg, HistTable), {V0Specs.at(kPdg)}); mHistogramRegistry->add(mcDir + getHistNameV2(kPdgMother, HistTable), getHistDesc(kPdgMother, HistTable), getHistType(kPdgMother, HistTable), {V0Specs.at(kPdgMother)}); mHistogramRegistry->add(mcDir + getHistNameV2(kPdgPartonicMother, HistTable), getHistDesc(kPdgPartonicMother, HistTable), getHistType(kPdgPartonicMother, HistTable), {V0Specs.at(kPdgPartonicMother)}); @@ -488,6 +504,7 @@ class V0HistManager mHistogramRegistry->add(mcDir + getHistNameV2(kPrimary, HistTable), getHistDesc(kPrimary, HistTable), getHistType(kPrimary, HistTable), {V0Specs.at(kPrimary)}); mHistogramRegistry->add(mcDir + getHistNameV2(kFromWrongCollision, HistTable), getHistDesc(kFromWrongCollision, HistTable), getHistType(kFromWrongCollision, HistTable), {V0Specs.at(kFromWrongCollision)}); mHistogramRegistry->add(mcDir + getHistNameV2(kFromMaterial, HistTable), getHistDesc(kFromMaterial, HistTable), getHistType(kFromMaterial, HistTable), {V0Specs.at(kFromMaterial)}); + mHistogramRegistry->add(mcDir + getHistNameV2(kMissidentified, HistTable), getHistDesc(kMissidentified, HistTable), getHistType(kMissidentified, HistTable), {V0Specs.at(kMissidentified)}); if (mPlotNSecondaries >= histmanager::kSecondaryPlotLevel1) { mHistogramRegistry->add(mcDir + getHistNameV2(kSecondary1, HistTable), getHistDesc(kSecondary1, HistTable), getHistType(kSecondary1, HistTable), {V0Specs.at(kSecondary1)}); @@ -569,79 +586,99 @@ class V0HistManager } } - template - void fillMc(T1 const& v0candidate, T2 const& /*mcParticles*/, T3 const& /*mcMothers*/, T4 const& /*mcPartonicMothers*/) + template + void fillMc(T1 const& v0Candidate, T2 const& /*mcParticles*/, T3 const& /*mcMothers*/, T4 const& /*mcPartonicMothers*/) { // No MC Particle - if (!v0candidate.has_fMcParticle()) { + if (!v0Candidate.has_fMcParticle()) { mHistogramRegistry->fill(HIST(v0Prefix) + HIST(McDir) + HIST(getHistName(kPdg, HistTable)), 0); mHistogramRegistry->fill(HIST(v0Prefix) + HIST(McDir) + HIST(getHistName(kOrigin, HistTable)), static_cast(modes::McOrigin::kNoMcParticle)); - if (mPlotOrigins) { - mHistogramRegistry->fill(HIST(v0Prefix) + HIST(McDir) + HIST(getHistName(kNoMcParticle, HistTable)), v0candidate.pt(), v0candidate.cosPa()); + if constexpr (modes::isFlagSet(mode, modes::Mode::kQa)) { + if (mPlotOrigins) { + mHistogramRegistry->fill(HIST(v0Prefix) + HIST(McDir) + HIST(getHistName(kNoMcParticle, HistTable)), v0Candidate.pt(), v0Candidate.cosPa()); + } } return; } // Retrieve MC particle - auto mcParticle = v0candidate.template fMcParticle_as(); - mHistogramRegistry->fill(HIST(v0Prefix) + HIST(McDir) + HIST(getHistName(kTruePt, HistTable)), mcParticle.pt()); - mHistogramRegistry->fill(HIST(v0Prefix) + HIST(McDir) + HIST(getHistName(kTrueEta, HistTable)), mcParticle.eta()); - mHistogramRegistry->fill(HIST(v0Prefix) + HIST(McDir) + HIST(getHistName(kTruePhi, HistTable)), mcParticle.phi()); - mHistogramRegistry->fill(HIST(v0Prefix) + HIST(McDir) + HIST(getHistName(kOrigin, HistTable)), mcParticle.origin()); + auto mcParticle = v0Candidate.template fMcParticle_as(); + + // missidentifed particles are special case + // whether a particle is missidentfied or not cannot be known by the producer so we check it here + bool isMissidentified = mcParticle.pdgCode() != mPdgCode; + + mHistogramRegistry->fill(HIST(v0Prefix) + HIST(McDir) + HIST(getHistName(kTruePtVsPt, HistTable)), mcParticle.pt(), v0Candidate.pt()); + mHistogramRegistry->fill(HIST(v0Prefix) + HIST(McDir) + HIST(getHistName(kTrueEtaVsEta, HistTable)), mcParticle.eta(), v0Candidate.eta()); + mHistogramRegistry->fill(HIST(v0Prefix) + HIST(McDir) + HIST(getHistName(kTruePhiVsPhi, HistTable)), mcParticle.phi(), v0Candidate.phi()); + if (isMissidentified) { + mHistogramRegistry->fill(HIST(v0Prefix) + HIST(McDir) + HIST(getHistName(kOrigin, HistTable)), static_cast(modes::McOrigin::kMissidentified)); + } else { + mHistogramRegistry->fill(HIST(v0Prefix) + HIST(McDir) + HIST(getHistName(kOrigin, HistTable)), mcParticle.origin()); + } mHistogramRegistry->fill(HIST(v0Prefix) + HIST(McDir) + HIST(getHistName(kPdg, HistTable)), mcParticle.pdgCode()); - // Mother PDG - if (v0candidate.has_fMcMother()) { - auto mother = v0candidate.template fMcMother_as(); + // get mother + if (v0Candidate.has_fMcMother()) { + auto mother = v0Candidate.template fMcMother_as(); mHistogramRegistry->fill(HIST(v0Prefix) + HIST(McDir) + HIST(getHistName(kPdgMother, HistTable)), mother.pdgCode()); } else { mHistogramRegistry->fill(HIST(v0Prefix) + HIST(McDir) + HIST(getHistName(kPdgMother, HistTable)), 0); } - // Partonic Mother PDG - if (v0candidate.has_fMcPartMoth()) { - auto partonicMother = v0candidate.template fMcPartMoth_as(); + // get partonic mother + if (v0Candidate.has_fMcPartMoth()) { + auto partonicMother = v0Candidate.template fMcPartMoth_as(); mHistogramRegistry->fill(HIST(v0Prefix) + HIST(McDir) + HIST(getHistName(kPdgPartonicMother, HistTable)), partonicMother.pdgCode()); } else { mHistogramRegistry->fill(HIST(v0Prefix) + HIST(McDir) + HIST(getHistName(kPdgPartonicMother, HistTable)), 0); } - // Plot origins - if (mPlotOrigins) { - switch (static_cast(mcParticle.origin())) { - case modes::McOrigin::kPhysicalPrimary: - mHistogramRegistry->fill(HIST(v0Prefix) + HIST(McDir) + HIST(getHistName(kPrimary, HistTable)), v0candidate.pt(), v0candidate.cosPa()); - break; - case modes::McOrigin::kFromWrongCollision: - mHistogramRegistry->fill(HIST(v0Prefix) + HIST(McDir) + HIST(getHistName(kFromWrongCollision, HistTable)), v0candidate.pt(), v0candidate.cosPa()); - break; - case modes::McOrigin::kFromMaterial: - mHistogramRegistry->fill(HIST(v0Prefix) + HIST(McDir) + HIST(getHistName(kFromMaterial, HistTable)), v0candidate.pt(), v0candidate.cosPa()); - break; - case modes::McOrigin::kFromSecondaryDecay: - if (v0candidate.has_fMcMother()) { - auto mother = v0candidate.template fMcMother_as(); - int motherPdgCode = std::abs(mother.pdgCode()); - // Switch on PDG of the mother - if (mPlotNSecondaries >= histmanager::kSecondaryPlotLevel1 && motherPdgCode == mPdgCodesSecondaryMother[0]) { - mHistogramRegistry->fill(HIST(v0Prefix) + HIST(McDir) + HIST(getHistName(kSecondary1, HistTable)), v0candidate.pt(), v0candidate.cosPa()); - } else if (mPlotNSecondaries >= histmanager::kSecondaryPlotLevel2 && motherPdgCode == mPdgCodesSecondaryMother[1]) { - mHistogramRegistry->fill(HIST(v0Prefix) + HIST(McDir) + HIST(getHistName(kSecondary2, HistTable)), v0candidate.pt(), v0candidate.cosPa()); - } else if (mPlotNSecondaries >= histmanager::kSecondaryPlotLevel3 && motherPdgCode == mPdgCodesSecondaryMother[2]) { - mHistogramRegistry->fill(HIST(v0Prefix) + HIST(McDir) + HIST(getHistName(kSecondary3, HistTable)), v0candidate.pt(), v0candidate.cosPa()); - } else { - mHistogramRegistry->fill(HIST(v0Prefix) + HIST(McDir) + HIST(getHistName(kSecondaryOther, HistTable)), v0candidate.pt(), v0candidate.cosPa()); - } + if constexpr (modes::isFlagSet(mode, modes::Mode::kQa)) { + if (mPlotOrigins) { + // check first if particle is missidentified + if (isMissidentified) { + // if it is, we fill it as such + mHistogramRegistry->fill(HIST(v0Prefix) + HIST(McDir) + HIST(getHistName(kMissidentified, HistTable)), v0Candidate.pt(), v0Candidate.cosPa()); + } else { + // if not, we fill it acccoridng to its origin + switch (static_cast(mcParticle.origin())) { + case modes::McOrigin::kPhysicalPrimary: + mHistogramRegistry->fill(HIST(v0Prefix) + HIST(McDir) + HIST(getHistName(kPrimary, HistTable)), v0Candidate.pt(), v0Candidate.cosPa()); + break; + case modes::McOrigin::kFromWrongCollision: + mHistogramRegistry->fill(HIST(v0Prefix) + HIST(McDir) + HIST(getHistName(kFromWrongCollision, HistTable)), v0Candidate.pt(), v0Candidate.cosPa()); + break; + case modes::McOrigin::kFromMaterial: + mHistogramRegistry->fill(HIST(v0Prefix) + HIST(McDir) + HIST(getHistName(kFromMaterial, HistTable)), v0Candidate.pt(), v0Candidate.cosPa()); + break; + case modes::McOrigin::kFromSecondaryDecay: + if (v0Candidate.has_fMcMother()) { + auto mother = v0Candidate.template fMcMother_as(); + int motherPdgCode = std::abs(mother.pdgCode()); + // Switch on PDG of the mother + if (mPlotNSecondaries >= histmanager::kSecondaryPlotLevel1 && motherPdgCode == mPdgCodesSecondaryMother[0]) { + mHistogramRegistry->fill(HIST(v0Prefix) + HIST(McDir) + HIST(getHistName(kSecondary1, HistTable)), v0Candidate.pt(), v0Candidate.cosPa()); + } else if (mPlotNSecondaries >= histmanager::kSecondaryPlotLevel2 && motherPdgCode == mPdgCodesSecondaryMother[1]) { + mHistogramRegistry->fill(HIST(v0Prefix) + HIST(McDir) + HIST(getHistName(kSecondary2, HistTable)), v0Candidate.pt(), v0Candidate.cosPa()); + } else if (mPlotNSecondaries >= histmanager::kSecondaryPlotLevel3 && motherPdgCode == mPdgCodesSecondaryMother[2]) { + mHistogramRegistry->fill(HIST(v0Prefix) + HIST(McDir) + HIST(getHistName(kSecondary3, HistTable)), v0Candidate.pt(), v0Candidate.cosPa()); + } else { + mHistogramRegistry->fill(HIST(v0Prefix) + HIST(McDir) + HIST(getHistName(kSecondaryOther, HistTable)), v0Candidate.pt(), v0Candidate.cosPa()); + } + } + break; + default: + LOG(warn) << "Encounted partilce with unknown origin!"; + break; } - break; - default: - // Unknown origin → safely ignore - break; + } } } } o2::framework::HistogramRegistry* mHistogramRegistry = nullptr; + int mPdgCode = 0; bool mPlot2d = false; bool mPlotOrigins = false; int mPlotNSecondaries = 0; diff --git a/PWGCF/Femto/DataModel/FemtoTables.h b/PWGCF/Femto/DataModel/FemtoTables.h index 81d2e4c7d24..9f742edf96c 100644 --- a/PWGCF/Femto/DataModel/FemtoTables.h +++ b/PWGCF/Femto/DataModel/FemtoTables.h @@ -673,15 +673,15 @@ using FOmegaExtras = FOmegaExtras_001; namespace femtomccollisions { -DECLARE_SOA_COLUMN(MultMc, multMc, int); //! Multiplicity of the event as given by the generator in |eta|<0.8 -DECLARE_SOA_COLUMN(CentMc, centMc, float); //! Multiplicity of the event as given by the generator in |eta|<0.8 - // +DECLARE_SOA_COLUMN(Mult, mult, int); //! Multiplicity of the event as given by the generator in |eta|<0.8 +DECLARE_SOA_COLUMN(Cent, cent, float); //! Multiplicity of the event as given by the generator in |eta|<0.8 + // } // namespace femtomccollisions DECLARE_SOA_TABLE_STAGED_VERSIONED(FMcCols_001, "FMCCOL", 1, //! femto mc collisions o2::soa::Index<>, - femtomccollisions::MultMc, - femtomccollisions::CentMc); + femtomccollisions::Mult, + femtomccollisions::Cent); using FMcCols = FMcCols_001; using FMcCol = FMcCols_001::iterator; @@ -717,7 +717,7 @@ DECLARE_SOA_TABLE_STAGED_VERSIONED(FMcMothers_001, "FMCMOTHER", 1, //! first dir using FMcMothers = FMcMothers_001; using FMcMother = FMcMothers::iterator; -DECLARE_SOA_TABLE_STAGED_VERSIONED(FMcPartMoths_001, "FMCPARTMOTH", 1, //! last partonic mother of the monte carlo particle +DECLARE_SOA_TABLE_STAGED_VERSIONED(FMcPartMoths_001, "FMCPARTMOTH", 1, //! first partonic mother of the monte carlo particle after hadronization o2::soa::Index<>, femtomcparticle::PdgCode); using FMcPartMoths = FMcPartMoths_001; diff --git a/PWGCF/Femto/Macros/cutculator.py b/PWGCF/Femto/Macros/cutculator.py index c1c6a4bf9e3..9e1d771b6da 100755 --- a/PWGCF/Femto/Macros/cutculator.py +++ b/PWGCF/Femto/Macros/cutculator.py @@ -139,7 +139,14 @@ def main(rootfile_path, tdir_path="femto-producer"): print("\nHistograms found in directory:") for i, hname in enumerate(histograms): print(f" [{i}] {hname}") - hidx = int(input("\nSelect histogram index: ")) + while True: + try: + hidx = int(input("\nSelect histogram index: ")) + if 0 <= hidx < len(histograms): + break + except ValueError: + pass + print("Invalid index.") hname = histograms[hidx] hist = d.Get(hname) nbins = hist.GetNbinsX() diff --git a/PWGCF/Femto/TableProducer/femtoProducer.cxx b/PWGCF/Femto/TableProducer/femtoProducer.cxx index f72870f2d28..36897c02afb 100644 --- a/PWGCF/Femto/TableProducer/femtoProducer.cxx +++ b/PWGCF/Femto/TableProducer/femtoProducer.cxx @@ -446,7 +446,7 @@ struct FemtoProducer { } PROCESS_SWITCH(FemtoProducer, processTracksV0sRun3ppMc, "Provide reconstructed and generated tracks and v0s", false); - // process monte carlo tracks and v0s + // process monte carlo tracks and kinks void processTracksKinksRun3ppMc(consumeddata::Run3PpMcRecoCollisions::iterator const& col, consumeddata::Run3PpMcGenCollisions const& mcCols, BCsWithTimestamps const& bcs, @@ -463,6 +463,26 @@ struct FemtoProducer { processMcKinks(col, mcCols, tracks, tracksWithItsPid, kinks, mcParticles); } PROCESS_SWITCH(FemtoProducer, processTracksKinksRun3ppMc, "Provide reconstructed and generated tracks and kinks", false); + + // process monte carlo tracks and v0s and kinks (adding cascades later here) + void processTracksV0sKinksRun3ppMc(consumeddata::Run3PpMcRecoCollisions::iterator const& col, + consumeddata::Run3PpMcGenCollisions const& mcCols, + BCsWithTimestamps const& bcs, + consumeddata::Run3McRecoTracks const& tracks, + consumeddata::Run3RecoVzeros const& v0s, + consumeddata::Run3Kinks const& kinks, + consumeddata::Run3McGenParticles const& mcParticles) + { + if (!processMcCollisions(col, mcCols, bcs, tracks)) { + return; + } + auto tracksWithItsPid = o2::soa::Attach(tracks); + processMcTracks(col, mcCols, tracks, tracksWithItsPid, mcParticles); + processMcV0s(col, mcCols, tracks, tracksWithItsPid, v0s, mcParticles); + processMcKinks(col, mcCols, tracks, tracksWithItsPid, kinks, mcParticles); + } + PROCESS_SWITCH(FemtoProducer, processTracksV0sKinksRun3ppMc, "Provide reconstructed and generated tracks and v0s and kinks", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) diff --git a/PWGCF/Femto/Tasks/CMakeLists.txt b/PWGCF/Femto/Tasks/CMakeLists.txt index 1cee38e9eee..b124fbb05d0 100644 --- a/PWGCF/Femto/Tasks/CMakeLists.txt +++ b/PWGCF/Femto/Tasks/CMakeLists.txt @@ -53,7 +53,7 @@ o2physics_add_dpl_workflow(femto-pair-track-cascade SOURCES femtoPairTrackCascade.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore COMPONENT_NAME Analysis) - +# o2physics_add_dpl_workflow(femto-pair-track-kink SOURCES femtoPairTrackKink.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore diff --git a/PWGCF/Femto/Tasks/femtoCascadeQa.cxx b/PWGCF/Femto/Tasks/femtoCascadeQa.cxx index 1a745e2d41b..2e9fd0c1b8c 100644 --- a/PWGCF/Femto/Tasks/femtoCascadeQa.cxx +++ b/PWGCF/Femto/Tasks/femtoCascadeQa.cxx @@ -105,27 +105,26 @@ struct FemtoCascadeQa { void init(InitContext&) { - // create a map for histogram specs + if ((doprocessXis + doprocessOmegas) > 1) { + LOG(fatal) << "Only one process can be activated"; + } + auto colHistSpec = colhistmanager::makeColQaHistSpecMap(confCollisionBinning, confCollisionQaBinning); colHistManager.init(&hRegistry, colHistSpec, confCollisionQaBinning); - auto bachelorHistSpec = trackhistmanager::makeTrackQaHistSpecMap(confBachelorBinning, confBachelorQaBinning); auto posDaughterHistSpec = trackhistmanager::makeTrackQaHistSpecMap(confPosDaughterBinning, confPosDaughterQaBinning); auto negDaughterHistSpec = trackhistmanager::makeTrackQaHistSpecMap(confNegDaughterBinning, confNegDaughterQaBinning); - if ((doprocessXis + doprocessOmegas) > 1) { - LOG(fatal) << "Only one process can be activated"; - } - if (doprocessXis) { auto xiHistSpec = cascadehistmanager::makeCascadeQaHistSpecMap(confXiBinning, confXiQaBinning); - xiHistManager.init(&hRegistry, xiHistSpec, confXiQaBinning, bachelorHistSpec, confBachelorQaBinning, posDaughterHistSpec, confPosDaughterQaBinning, negDaughterHistSpec, confNegDaughterQaBinning); + xiHistManager.init(&hRegistry, xiHistSpec, confXiSelection, confXiQaBinning, bachelorHistSpec, confBachelorQaBinning, posDaughterHistSpec, confPosDaughterQaBinning, negDaughterHistSpec, confNegDaughterQaBinning); } if (doprocessOmegas) { auto omegaHistSpec = cascadehistmanager::makeCascadeQaHistSpecMap(confOmegaBinning, confOmegaQaBinning); - omegaHistManager.init(&hRegistry, omegaHistSpec, confOmegaQaBinning, bachelorHistSpec, confBachelorQaBinning, posDaughterHistSpec, confPosDaughterQaBinning, negDaughterHistSpec, confNegDaughterQaBinning); + omegaHistManager.init(&hRegistry, omegaHistSpec, confOmegaSelection, confOmegaQaBinning, bachelorHistSpec, confBachelorQaBinning, posDaughterHistSpec, confPosDaughterQaBinning, negDaughterHistSpec, confNegDaughterQaBinning); } + hRegistry.print(); }; void processXis(FilteredFemtoCollision const& col, FemtoXis const& /*xis*/, FemtoTracks const& tracks) diff --git a/PWGCF/Femto/Tasks/femtoKinkQa.cxx b/PWGCF/Femto/Tasks/femtoKinkQa.cxx index 2dd445adf5a..7887b563463 100644 --- a/PWGCF/Femto/Tasks/femtoKinkQa.cxx +++ b/PWGCF/Femto/Tasks/femtoKinkQa.cxx @@ -34,7 +34,6 @@ #include "Framework/OutputObjHeader.h" #include "Framework/runDataProcessing.h" -#include #include #include @@ -84,7 +83,6 @@ struct FemtoKinkQa { kinkhistmanager::ConfSigmaBinning1 confSigmaBinning; kinkhistmanager::ConfSigmaQaBinning1 confSigmaQaBinning; - kinkhistmanager::ConfSigmaMcBinning confSigmaMcBinning; kinkhistmanager::KinkHistManager< kinkhistmanager::PrefixSigmaQa, trackhistmanager::PrefixKinkChaDaughterQa, @@ -102,7 +100,6 @@ struct FemtoKinkQa { kinkhistmanager::ConfSigmaPlusBinning1 confSigmaPlusBinning; kinkhistmanager::ConfSigmaPlusQaBinning1 confSigmaPlusQaBinning; - kinkhistmanager::ConfSigmaPlusMcBinning confSigmaPlusMcBinning; kinkhistmanager::KinkHistManager< kinkhistmanager::PrefixSigmaPlusQa, trackhistmanager::PrefixKinkChaDaughterQa, @@ -112,7 +109,6 @@ struct FemtoKinkQa { // setup for daughters trackhistmanager::ConfKinkChaDauBinning confKinkChaDaughterBinning; trackhistmanager::ConfKinkChaDauQaBinning confKinkChaDaughterQaBinning; - trackhistmanager::ConfKinkChaDauMcBinning confKinkChaDaughterMcBinning; HistogramRegistry hRegistry{"FemtoKinkQa", {}, OutputObjHandlingPolicy::AnalysisObject}; @@ -130,23 +126,23 @@ struct FemtoKinkQa { auto chaDauHistSpec = trackhistmanager::makeTrackQaHistSpecMap(confKinkChaDaughterBinning, confKinkChaDaughterQaBinning); if (doprocessSigma) { auto sigmaHistSpec = kinkhistmanager::makeKinkQaHistSpecMap(confSigmaBinning, confSigmaQaBinning); - sigmaHistManager.init(&hRegistry, sigmaHistSpec, confSigmaQaBinning, chaDauHistSpec, confKinkChaDaughterQaBinning); + sigmaHistManager.init(&hRegistry, sigmaHistSpec, confSigmaSelection, confSigmaQaBinning, chaDauHistSpec, confKinkChaDaughterQaBinning); } if (doprocessSigmaPlus) { auto sigmaPlusHistSpec = kinkhistmanager::makeKinkQaHistSpecMap(confSigmaPlusBinning, confSigmaPlusQaBinning); - sigmaPlusHistManager.init(&hRegistry, sigmaPlusHistSpec, confSigmaPlusQaBinning, chaDauHistSpec, confKinkChaDaughterQaBinning); + sigmaPlusHistManager.init(&hRegistry, sigmaPlusHistSpec, confSigmaPlusSelection, confSigmaPlusQaBinning, chaDauHistSpec, confKinkChaDaughterQaBinning); } } else { auto colHistSpec = colhistmanager::makeColMcQaHistSpecMap(confCollisionBinning, confCollisionQaBinning); colHistManager.init(&hRegistry, colHistSpec, confCollisionQaBinning); - auto chaDauHistSpec = trackhistmanager::makeTrackMcQaHistSpecMap(confKinkChaDaughterBinning, confKinkChaDaughterQaBinning, confKinkChaDaughterMcBinning); + auto chaDauHistSpec = trackhistmanager::makeTrackMcQaHistSpecMap(confKinkChaDaughterBinning, confKinkChaDaughterQaBinning); if (doprocessSigmaMc) { - auto sigmaHistSpec = kinkhistmanager::makeKinkMcQaHistSpecMap(confSigmaBinning, confSigmaQaBinning, confSigmaMcBinning); - sigmaHistManager.init(&hRegistry, sigmaHistSpec, confSigmaQaBinning, confSigmaMcBinning, chaDauHistSpec, confKinkChaDaughterQaBinning, confKinkChaDaughterMcBinning); + auto sigmaHistSpec = kinkhistmanager::makeKinkMcQaHistSpecMap(confSigmaBinning, confSigmaQaBinning); + sigmaHistManager.init(&hRegistry, sigmaHistSpec, confSigmaSelection, confSigmaQaBinning, chaDauHistSpec, confKinkChaDaughterQaBinning); } if (doprocessSigmaPlusMc) { - auto sigmaPlusHistSpec = kinkhistmanager::makeKinkMcQaHistSpecMap(confSigmaPlusBinning, confSigmaPlusQaBinning, confSigmaPlusMcBinning); - sigmaPlusHistManager.init(&hRegistry, sigmaPlusHistSpec, confSigmaPlusQaBinning, confSigmaPlusMcBinning, chaDauHistSpec, confKinkChaDaughterQaBinning, confKinkChaDaughterMcBinning); + auto sigmaPlusHistSpec = kinkhistmanager::makeKinkMcQaHistSpecMap(confSigmaPlusBinning, confSigmaPlusQaBinning); + sigmaPlusHistManager.init(&hRegistry, sigmaPlusHistSpec, confSigmaPlusSelection, confSigmaPlusQaBinning, chaDauHistSpec, confKinkChaDaughterQaBinning); } } hRegistry.print(); diff --git a/PWGCF/Femto/Tasks/femtoPairTrackKink.cxx b/PWGCF/Femto/Tasks/femtoPairTrackKink.cxx index 546cd51df97..438c91b474c 100644 --- a/PWGCF/Femto/Tasks/femtoPairTrackKink.cxx +++ b/PWGCF/Femto/Tasks/femtoPairTrackKink.cxx @@ -51,15 +51,21 @@ using namespace o2::analysis::femto; struct FemtoPairTrackKink { // setup tables - using Collisions = Join; - using Collision = Collisions::iterator; + using FemtoCollisions = Join; + using FilteredFemtoCollisions = o2::soa::Filtered; + using FilteredFemtoCollision = FilteredFemtoCollisions::iterator; - using FilteredCollisions = o2::soa::Filtered; - using FilteredCollision = FilteredCollisions::iterator; + using FemtoCollisionsWithLabel = o2::soa::Join; + using FilteredFemtoCollisionsWithLabel = o2::soa::Filtered; + using FilteredFemtoCollisionWithLabel = FilteredFemtoCollisionsWithLabel::iterator; - using Tracks = o2::soa::Join; - using Sigmas = o2::soa::Join; - using SigmaPlus = o2::soa::Join; + using FemtoTracks = o2::soa::Join; + using FemtoSigmas = o2::soa::Join; + using FemtoSigmaPlus = o2::soa::Join; + + using FemtoTracksWithLabel = o2::soa::Join; + using FemtoSigmasWithLabel = o2::soa::Join; + using FemtoSigmaPlusWithLabel = o2::soa::Join; SliceCache cache; @@ -69,25 +75,37 @@ struct FemtoPairTrackKink { colhistmanager::ConfCollisionBinning confCollisionBinning; // setup tracks - trackbuilder::ConfTrackSelection1 trackSelection; + trackbuilder::ConfTrackSelection1 confTrackSelection; trackhistmanager::ConfTrackBinning1 confTrackBinning; - Partition trackPartition = MAKE_TRACK_PARTITION(trackSelection); - Preslice perColTracks = aod::femtobase::stored::fColId; + + Partition trackPartition = MAKE_TRACK_PARTITION(confTrackSelection); + Preslice perColTracks = aod::femtobase::stored::fColId; + + Partition trackWithLabelPartition = MAKE_TRACK_PARTITION(confTrackSelection); + Preslice perColtracksWithLabel = aod::femtobase::stored::fColId; // setup for daughters trackhistmanager::ConfKinkChaDauBinning confChaDauBinning; // setup sigmas - kinkbuilder::ConfSigmaSelection1 sigmaSelection; + kinkbuilder::ConfSigmaSelection1 confSigmaSelection; kinkhistmanager::ConfSigmaBinning1 confSigmaBinning; - Partition sigmaPartition = MAKE_SIGMA_PARTITION(sigmaSelection); - Preslice perColSigmas = aod::femtobase::stored::fColId; + + Partition sigmaPartition = MAKE_SIGMA_PARTITION(confSigmaSelection); + Preslice perColSigmas = aod::femtobase::stored::fColId; + + Partition sigmaWithLabelPartition = MAKE_SIGMA_PARTITION(confSigmaSelection); + Preslice perColSigmasWithLabel = aod::femtobase::stored::fColId; // setup for sigma plus - kinkbuilder::ConfSigmaPlusSelection1 sigmaPlusSelection; + kinkbuilder::ConfSigmaPlusSelection1 confSigmaPlusSelection; kinkhistmanager::ConfSigmaPlusBinning1 confSigmaPlusBinning; - Partition sigmaPlusPartition = MAKE_SIGMAPLUS_PARTITION(sigmaPlusSelection); - Preslice perColSigmaPlus = aod::femtobase::stored::fColId; + + Partition sigmaPlusPartition = MAKE_SIGMAPLUS_PARTITION(confSigmaPlusSelection); + Preslice perColSigmaPlus = aod::femtobase::stored::fColId; + + Partition sigmaPlusWithLabelPartition = MAKE_SIGMAPLUS_PARTITION(confSigmaPlusSelection); + Preslice perColSigmaPlusWithLabel = aod::femtobase::stored::fColId; // setup pairs pairhistmanager::ConfPairBinning confPairBinning; @@ -131,6 +149,19 @@ struct FemtoPairTrackKink { void init(InitContext&) { + bool processData = doprocessSigmaSameEvent || doprocessSigmaMixedEvent || doprocessSigmaPlusSameEvent || doprocessSigmaPlusMixedEvent; + bool processMc = doprocessSigmaSameEventMc || doprocessSigmaMixedEventMc || doprocessSigmaPlusSameEventMc || doprocessSigmaPlusMixedEventMc; + + if (processData && processMc) { + LOG(fatal) << "Both data and mc processing is enabled. Breaking..."; + } + + bool processSigma = doprocessSigmaSameEvent || doprocessSigmaSameEventMc || doprocessSigmaMixedEvent || doprocessSigmaMixedEventMc; + bool processSigmaPlus = doprocessSigmaPlusSameEvent || doprocessSigmaPlusSameEventMc || doprocessSigmaPlusMixedEvent || doprocessSigmaPlusMixedEventMc; + + if (processSigma && processSigmaPlus) { + LOG(fatal) << "Both Sigma-track and SigmaPlus-track processing is enabled. Breaking..."; + } // setup columnpolicy for binning // default values are used during instantiation, so we need to explicity update them here @@ -138,55 +169,83 @@ struct FemtoPairTrackKink { mixBinsVtxCent = {{confMixing.vtxBins.value, confMixing.centBins.value}, true}; mixBinsVtxMultCent = {{confMixing.vtxBins.value, confMixing.multBins.value, confMixing.centBins.value}, true}; - // setup histograms - auto colHistSpec = colhistmanager::makeColHistSpecMap(confCollisionBinning); - auto trackHistSpec = trackhistmanager::makeTrackHistSpecMap(confTrackBinning); - auto chaDauSpec = trackhistmanager::makeTrackHistSpecMap(confChaDauBinning); - auto pairHistSpec = pairhistmanager::makePairHistSpecMap(confPairBinning); auto cprHistSpec = closepairrejection::makeCprHistSpecMap(confCpr); - // setup for sigma - if (doprocessSigmaSameEvent || doprocessSigmaMixedEvent) { - auto sigmaHistSpec = kinkhistmanager::makeKinkHistSpecMap(confSigmaBinning); - auto pairTrackSigmaHistSpec = pairhistmanager::makePairHistSpecMap(confPairBinning); - pairTrackSigmaBuilder.init(&hRegistry, trackSelection, sigmaSelection, confCpr, confMixing, confPairBinning, confPairCuts, colHistSpec, trackHistSpec, sigmaHistSpec, chaDauSpec, pairTrackSigmaHistSpec, cprHistSpec); - } - - // setup for sigma plus - if (doprocessSigmaPlusSameEvent || doprocessSigmaPlusMixedEvent) { - auto sigmaplusHistSpec = kinkhistmanager::makeKinkHistSpecMap(confSigmaPlusBinning); - auto pairTrackSigmaPlusHistSpec = pairhistmanager::makePairHistSpecMap(confPairBinning); - pairTrackSigmaPlusBuilder.init(&hRegistry, trackSelection, sigmaPlusSelection, confCpr, confMixing, confPairBinning, confPairCuts, colHistSpec, trackHistSpec, sigmaplusHistSpec, chaDauSpec, pairTrackSigmaPlusHistSpec, cprHistSpec); - } - - if (((doprocessSigmaSameEvent || doprocessSigmaMixedEvent) + (doprocessSigmaPlusSameEvent || doprocessSigmaPlusMixedEvent)) > 1) { - LOG(fatal) << "Can only process sigma-tracks Or sigmaplus-tracks"; + if (processData) { + auto colHistSpec = colhistmanager::makeColHistSpecMap(confCollisionBinning); + auto trackHistSpec = trackhistmanager::makeTrackHistSpecMap(confTrackBinning); + auto chaDauSpec = trackhistmanager::makeTrackHistSpecMap(confChaDauBinning); + auto pairTrackKinkHistSpec = pairhistmanager::makePairHistSpecMap(confPairBinning); + if (processSigma) { + auto sigmaHistSpec = kinkhistmanager::makeKinkHistSpecMap(confSigmaBinning); + pairTrackSigmaBuilder.init(&hRegistry, confTrackSelection, confSigmaSelection, confCpr, confMixing, confPairBinning, confPairCuts, colHistSpec, trackHistSpec, sigmaHistSpec, chaDauSpec, pairTrackKinkHistSpec, cprHistSpec); + } else { + auto sigmaplusHistSpec = kinkhistmanager::makeKinkHistSpecMap(confSigmaPlusBinning); + pairTrackSigmaPlusBuilder.init(&hRegistry, confTrackSelection, confSigmaPlusSelection, confCpr, confMixing, confPairBinning, confPairCuts, colHistSpec, trackHistSpec, sigmaplusHistSpec, chaDauSpec, pairTrackKinkHistSpec, cprHistSpec); + } + } else { + auto colHistSpec = colhistmanager::makeColMcHistSpecMap(confCollisionBinning); + auto trackHistSpec = trackhistmanager::makeTrackMcHistSpecMap(confTrackBinning); + auto chaDauSpec = trackhistmanager::makeTrackMcHistSpecMap(confChaDauBinning); + auto pairTrackKinkHistSpec = pairhistmanager::makePairMcHistSpecMap(confPairBinning); + if (processSigma) { + auto sigmaHistSpec = kinkhistmanager::makeKinkMcHistSpecMap(confSigmaBinning); + pairTrackSigmaBuilder.init(&hRegistry, confTrackSelection, confSigmaSelection, confCpr, confMixing, confPairBinning, confPairCuts, colHistSpec, trackHistSpec, sigmaHistSpec, chaDauSpec, pairTrackKinkHistSpec, cprHistSpec); + } else { + auto sigmaplusHistSpec = kinkhistmanager::makeKinkMcHistSpecMap(confSigmaPlusBinning); + pairTrackSigmaPlusBuilder.init(&hRegistry, confTrackSelection, confSigmaPlusSelection, confCpr, confMixing, confPairBinning, confPairCuts, colHistSpec, trackHistSpec, sigmaplusHistSpec, chaDauSpec, pairTrackKinkHistSpec, cprHistSpec); + } } + hRegistry.print(); }; - void processSigmaSameEvent(FilteredCollision const& col, Tracks const& tracks, Sigmas const& sigmas) + void processSigmaSameEvent(FilteredFemtoCollision const& col, FemtoTracks const& tracks, FemtoSigmas const& sigmas) { pairTrackSigmaBuilder.processSameEvent(col, tracks, trackPartition, sigmas, sigmaPartition, cache); } PROCESS_SWITCH(FemtoPairTrackKink, processSigmaSameEvent, "Enable processing same event processing for tracks and sigmas", true); - void processSigmaMixedEvent(FilteredCollisions const& cols, Tracks const& tracks, Sigmas const& /*sigmas*/) + void processSigmaSameEventMc(FilteredFemtoCollisionWithLabel const& col, FMcCols const& mcCols, FemtoTracksWithLabel const& tracks, FemtoSigmasWithLabel const& sigmas, FMcParticles const& mcParticles, FMcMothers const& mcMothers, FMcPartMoths const& mcPartonicMothers) + { + pairTrackSigmaBuilder.processSameEvent(col, mcCols, tracks, trackWithLabelPartition, sigmas, sigmaWithLabelPartition, mcParticles, mcMothers, mcPartonicMothers, cache); + } + PROCESS_SWITCH(FemtoPairTrackKink, processSigmaSameEventMc, "Enable processing same event processing for tracks and sigmas with MC information", false); + + void processSigmaMixedEvent(FilteredFemtoCollisions const& cols, FemtoTracks const& tracks, FemtoSigmas const& /*sigmas*/) { pairTrackSigmaBuilder.processMixedEvent(cols, tracks, trackPartition, sigmaPartition, cache, mixBinsVtxMult, mixBinsVtxCent, mixBinsVtxMultCent); } PROCESS_SWITCH(FemtoPairTrackKink, processSigmaMixedEvent, "Enable processing mixed event processing for tracks and sigmas", true); - // - void processSigmaPlusSameEvent(FilteredCollision const& col, Tracks const& tracks, SigmaPlus const& sigmaplus) + + void processSigmaMixedEventMc(FilteredFemtoCollisionsWithLabel const& cols, FMcCols const& mcCols, FemtoTracksWithLabel const& tracks, FemtoSigmasWithLabel const& /*sigmas*/, FMcParticles const& mcParticles) + { + pairTrackSigmaBuilder.processMixedEvent(cols, mcCols, tracks, trackWithLabelPartition, sigmaWithLabelPartition, mcParticles, cache, mixBinsVtxMult, mixBinsVtxCent, mixBinsVtxMultCent); + } + PROCESS_SWITCH(FemtoPairTrackKink, processSigmaMixedEventMc, "Enable processing mixed event processing for tracks and sigmas with MC information", false); + + void processSigmaPlusSameEvent(FilteredFemtoCollision const& col, FemtoTracks const& tracks, FemtoSigmaPlus const& sigmaplus) { pairTrackSigmaPlusBuilder.processSameEvent(col, tracks, trackPartition, sigmaplus, sigmaPlusPartition, cache); } PROCESS_SWITCH(FemtoPairTrackKink, processSigmaPlusSameEvent, "Enable processing same event processing for tracks and sigma plus", false); - void processSigmaPlusMixedEvent(FilteredCollisions const& cols, Tracks const& tracks, SigmaPlus const& /*sigmaplus*/) + void processSigmaPlusSameEventMc(FilteredFemtoCollisionWithLabel const& col, FMcCols const& mcCols, FemtoTracksWithLabel const& tracks, FemtoSigmaPlusWithLabel const& sigmaplus, FMcParticles const& mcParticles, FMcMothers const& mcMothers, FMcPartMoths const& mcPartonicMothers) + { + pairTrackSigmaPlusBuilder.processSameEvent(col, mcCols, tracks, trackWithLabelPartition, sigmaplus, sigmaPlusWithLabelPartition, mcParticles, mcMothers, mcPartonicMothers, cache); + } + PROCESS_SWITCH(FemtoPairTrackKink, processSigmaPlusSameEventMc, "Enable processing same event processing for tracks and sigma plus with MC information", false); + + void processSigmaPlusMixedEvent(FilteredFemtoCollisions const& cols, FemtoTracks const& tracks, FemtoSigmaPlus const& /*sigmaplus*/) { pairTrackSigmaPlusBuilder.processMixedEvent(cols, tracks, trackPartition, sigmaPlusPartition, cache, mixBinsVtxMult, mixBinsVtxCent, mixBinsVtxMultCent); } PROCESS_SWITCH(FemtoPairTrackKink, processSigmaPlusMixedEvent, "Enable processing mixed event processing for tracks and sigma plus", false); + + void processSigmaPlusMixedEventMc(FilteredFemtoCollisionsWithLabel const& cols, FMcCols const& mcCols, FemtoTracksWithLabel const& tracks, FemtoSigmaPlusWithLabel const& /*sigmaplus*/, FMcParticles const& mcParticles) + { + pairTrackSigmaPlusBuilder.processMixedEvent(cols, mcCols, tracks, trackWithLabelPartition, sigmaPlusWithLabelPartition, mcParticles, cache, mixBinsVtxMult, mixBinsVtxCent, mixBinsVtxMultCent); + } + PROCESS_SWITCH(FemtoPairTrackKink, processSigmaPlusMixedEventMc, "Enable processing mixed event processing for tracks and sigma plus with MC information", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) diff --git a/PWGCF/Femto/Tasks/femtoPairTrackTrack.cxx b/PWGCF/Femto/Tasks/femtoPairTrackTrack.cxx index a14c5b38590..f8a647dfcfe 100644 --- a/PWGCF/Femto/Tasks/femtoPairTrackTrack.cxx +++ b/PWGCF/Femto/Tasks/femtoPairTrackTrack.cxx @@ -48,13 +48,17 @@ using namespace o2::analysis::femto; struct FemtoPairTrackTrack { // setup tables - using Collisions = Join; - using Collision = Collisions::iterator; + using FemtoCollisions = Join; + using FilteredFemtoCollisions = o2::soa::Filtered; + using FilteredFemtoCollision = FilteredFemtoCollisions::iterator; - using FilteredCollisions = o2::soa::Filtered; - using FilteredCollision = FilteredCollisions::iterator; + using FemtoCollisionsWithLabel = o2::soa::Join; + using FilteredFemtoCollisionsWithLabel = o2::soa::Filtered; + using FilteredFemtoCollisionWithLabel = FilteredFemtoCollisionsWithLabel::iterator; - using Tracks = o2::soa::Join; + using FemtoTracks = o2::soa::Join; + + using FemtoTracksWithLabel = o2::soa::Join; SliceCache cache; @@ -64,15 +68,18 @@ struct FemtoPairTrackTrack { colhistmanager::ConfCollisionBinning confCollisionBinning; // setup tracks - trackbuilder::ConfTrackSelection1 trackSelections1; + trackbuilder::ConfTrackSelection1 confTrackSelections1; trackhistmanager::ConfTrackBinning1 confTrackBinning1; - trackbuilder::ConfTrackSelection2 trackSelections2; + trackbuilder::ConfTrackSelection2 confTrackSelections2; trackhistmanager::ConfTrackBinning2 confTrackBinning2; - Partition trackPartition1 = MAKE_TRACK_PARTITION(trackSelections1); - Partition trackPartition2 = MAKE_TRACK_PARTITION(trackSelections2); + Partition trackPartition1 = MAKE_TRACK_PARTITION(confTrackSelections1); + Partition trackPartition2 = MAKE_TRACK_PARTITION(confTrackSelections2); + Preslice perColtracks = aod::femtobase::stored::fColId; - Preslice perColReco = aod::femtobase::stored::fColId; + Partition trackWithLabelPartition1 = MAKE_TRACK_PARTITION(confTrackSelections1); + Partition trackWithLabelPartition2 = MAKE_TRACK_PARTITION(confTrackSelections2); + Preslice perColtracksWithLabel = aod::femtobase::stored::fColId; // setup pairs pairhistmanager::ConfPairBinning confPairBinning; @@ -102,6 +109,15 @@ struct FemtoPairTrackTrack { void init(InitContext&) { + if ((doprocessSameEvent + doprocessSameEventMc) > 1 || (doprocessMixedEvent + doprocessMixedEventMc) > 1) { + LOG(fatal) << "More than 1 same or mixed event process function is activated. Breaking..."; + } + bool processData = doprocessSameEvent || doprocessMixedEvent; + bool processMc = doprocessSameEventMc || doprocessMixedEventMc; + if (processData && processMc) { + LOG(fatal) << "Both data and mc processing is activated. Breaking..."; + } + // setup columnpolicy for binning // default values are used during instantiation, so we need to explicity update them here mixBinsVtxMult = {{confMixing.vtxBins, confMixing.multBins.value}, true}; @@ -109,26 +125,47 @@ struct FemtoPairTrackTrack { mixBinsVtxMultCent = {{confMixing.vtxBins.value, confMixing.multBins.value, confMixing.centBins.value}, true}; // setup histogram specs - auto colHistSpec = colhistmanager::makeColHistSpecMap(confCollisionBinning); - auto trackHistSpec1 = trackhistmanager::makeTrackHistSpecMap(confTrackBinning1); - auto trackHistSpec2 = trackhistmanager::makeTrackHistSpecMap(confTrackBinning2); - auto pairHistSpec = pairhistmanager::makePairHistSpecMap(confPairBinning); auto cprHistSpec = closepairrejection::makeCprHistSpecMap(confCpr); - pairTrackTrackBuilder.init(&hRegistry, trackSelections1, trackSelections2, confCpr, confMixing, confPairBinning, confPairCuts, colHistSpec, trackHistSpec1, trackHistSpec2, pairHistSpec, cprHistSpec); + if (processData) { + auto colHistSpec = colhistmanager::makeColHistSpecMap(confCollisionBinning); + auto trackHistSpec1 = trackhistmanager::makeTrackHistSpecMap(confTrackBinning1); + auto trackHistSpec2 = trackhistmanager::makeTrackHistSpecMap(confTrackBinning2); + auto pairHistSpec = pairhistmanager::makePairHistSpecMap(confPairBinning); + pairTrackTrackBuilder.init(&hRegistry, confTrackSelections1, confTrackSelections2, confCpr, confMixing, confPairBinning, confPairCuts, colHistSpec, trackHistSpec1, trackHistSpec2, pairHistSpec, cprHistSpec); + } else { + auto colHistSpec = colhistmanager::makeColMcHistSpecMap(confCollisionBinning); + auto trackHistSpec1 = trackhistmanager::makeTrackMcHistSpecMap(confTrackBinning1); + auto trackHistSpec2 = trackhistmanager::makeTrackMcHistSpecMap(confTrackBinning2); + auto pairHistSpec = pairhistmanager::makePairMcHistSpecMap(confPairBinning); + pairTrackTrackBuilder.init(&hRegistry, confTrackSelections1, confTrackSelections2, confCpr, confMixing, confPairBinning, confPairCuts, colHistSpec, trackHistSpec1, trackHistSpec2, pairHistSpec, cprHistSpec); + } + hRegistry.print(); }; - void processSameEvent(FilteredCollision const& col, Tracks const& tracks) + void processSameEvent(FilteredFemtoCollision const& col, FemtoTracks const& tracks) { pairTrackTrackBuilder.processSameEvent(col, tracks, trackPartition1, trackPartition2, cache); } PROCESS_SWITCH(FemtoPairTrackTrack, processSameEvent, "Enable processing same event processing", true); - void processMixedEvent(FilteredCollisions const& cols, Tracks const& tracks) + void processSameEventMc(FilteredFemtoCollisionWithLabel const& col, FMcCols const& mcCols, FemtoTracksWithLabel const& tracks, FMcParticles const& mcParticles, FMcMothers const& mcMothers, FMcPartMoths const& mcPartonicMothers) + { + pairTrackTrackBuilder.processSameEvent(col, mcCols, tracks, trackWithLabelPartition1, trackWithLabelPartition2, mcParticles, mcMothers, mcPartonicMothers, cache); + } + PROCESS_SWITCH(FemtoPairTrackTrack, processSameEventMc, "Enable processing same event processing", false); + + void processMixedEvent(FilteredFemtoCollisions const& cols, FemtoTracks const& tracks) { pairTrackTrackBuilder.processMixedEvent(cols, tracks, trackPartition1, trackPartition2, cache, mixBinsVtxMult, mixBinsVtxCent, mixBinsVtxMultCent); } PROCESS_SWITCH(FemtoPairTrackTrack, processMixedEvent, "Enable processing mixed event processing", true); + + void processMixedEventMc(FilteredFemtoCollisionsWithLabel const& cols, FMcCols const& mcCols, FemtoTracksWithLabel const& tracks, FMcParticles const& mcParticles) + { + pairTrackTrackBuilder.processMixedEvent(cols, mcCols, tracks, trackWithLabelPartition1, trackWithLabelPartition2, mcParticles, cache, mixBinsVtxMult, mixBinsVtxCent, mixBinsVtxMultCent); + } + PROCESS_SWITCH(FemtoPairTrackTrack, processMixedEventMc, "Enable processing mixed event processing", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) diff --git a/PWGCF/Femto/Tasks/femtoPairTrackTwoTrackResonance.cxx b/PWGCF/Femto/Tasks/femtoPairTrackTwoTrackResonance.cxx index 28afbed7b9e..e3a43c5b67e 100644 --- a/PWGCF/Femto/Tasks/femtoPairTrackTwoTrackResonance.cxx +++ b/PWGCF/Femto/Tasks/femtoPairTrackTwoTrackResonance.cxx @@ -157,7 +157,7 @@ struct FemtoPairTrackTwoTrackResonance { { if (((doprocessPhiSameEvent || doprocessPhiMixedEvent) + (doprocessKstar0SameEvent || doprocessKstar0MixedEvent)) + (doprocessRho0SameEvent || doprocessRho0MixedEvent) > 1) { - LOG(fatal) << "Can only process lambda-tracks Or k0short-tracks"; + LOG(fatal) << "Can only process phi-tracks, rho-tracks or k0*-tracks"; } // setup columnpolicy for binning diff --git a/PWGCF/Femto/Tasks/femtoPairTrackV0.cxx b/PWGCF/Femto/Tasks/femtoPairTrackV0.cxx index 7ea0cf3846c..550d9778c6b 100644 --- a/PWGCF/Femto/Tasks/femtoPairTrackV0.cxx +++ b/PWGCF/Femto/Tasks/femtoPairTrackV0.cxx @@ -50,15 +50,21 @@ using namespace o2::analysis::femto; struct FemtoPairTrackV0 { // setup tables - using Collisions = Join; - using Collision = Collisions::iterator; + using FemtoCollisions = Join; + using FilteredFemtoCollisions = o2::soa::Filtered; + using FilteredFemtoCollision = FilteredFemtoCollisions::iterator; - using FilteredCollisions = o2::soa::Filtered; - using FilteredCollision = FilteredCollisions::iterator; + using FemtoCollisionsWithLabel = o2::soa::Join; + using FilteredFemtoCollisionsWithLabel = o2::soa::Filtered; + using FilteredFemtoCollisionWithLabel = FilteredFemtoCollisionsWithLabel::iterator; - using Tracks = o2::soa::Join; - using Lambdas = o2::soa::Join; - using K0shorts = o2::soa::Join; + using FemtoTracks = o2::soa::Join; + using FemtoLambdas = o2::soa::Join; + using FemtoK0shorts = o2::soa::Join; + + using FemtoTracksWithLabel = o2::soa::Join; + using FemtoLambdasWithLabel = o2::soa::Join; + using FemtoK0shortsWithLabel = o2::soa::Join; SliceCache cache; @@ -68,10 +74,14 @@ struct FemtoPairTrackV0 { colhistmanager::ConfCollisionBinning confCollisionBinning; // setup tracks - trackbuilder::ConfTrackSelection1 trackSelection; + trackbuilder::ConfTrackSelection1 confTrackSelection; trackhistmanager::ConfTrackBinning1 confTrackBinning; - Partition trackPartition = MAKE_TRACK_PARTITION(trackSelection); - Preslice perColTracks = aod::femtobase::stored::fColId; + + Partition trackPartition = MAKE_TRACK_PARTITION(confTrackSelection); + Preslice perColTracks = aod::femtobase::stored::fColId; + + Partition trackWithLabelPartition = MAKE_TRACK_PARTITION(confTrackSelection); + Preslice perColtracksWithLabel = aod::femtobase::stored::fColId; // setup for daughters trackhistmanager::ConfV0PosDauBinning confPosDauBinning; @@ -80,14 +90,22 @@ struct FemtoPairTrackV0 { // setup lambdas v0builder::ConfLambdaSelection1 lambdaSelection; v0histmanager::ConfLambdaBinning1 confLambdaBinning; - Partition lambdaPartition = MAKE_LAMBDA_PARTITION(lambdaSelection); - Preslice perColLambdas = aod::femtobase::stored::fColId; + + Partition lambdaPartition = MAKE_LAMBDA_PARTITION(lambdaSelection); + Preslice perColLambdas = aod::femtobase::stored::fColId; + + Partition lambdaWithLabelPartition = MAKE_LAMBDA_PARTITION(lambdaSelection); + Preslice perColLambdasWithLabel = aod::femtobase::stored::fColId; // setup k0shorts v0builder::ConfK0shortSelection1 k0shortSelection; v0histmanager::ConfK0shortBinning1 confK0shortBinning; - Partition k0shortPartition = MAKE_K0SHORT_PARTITION(k0shortSelection); - Preslice perColk0shorts = aod::femtobase::stored::fColId; + + Partition k0shortPartition = MAKE_K0SHORT_PARTITION(k0shortSelection); + Preslice perColk0shorts = aod::femtobase::stored::fColId; + + Partition k0shortWithLabelPartition = MAKE_K0SHORT_PARTITION(k0shortSelection); + Preslice perColk0shortsWithLabel = aod::femtobase::stored::fColId; // setup pairs pairhistmanager::ConfPairBinning confPairBinning; @@ -133,6 +151,19 @@ struct FemtoPairTrackV0 { void init(InitContext&) { + bool processData = doprocessLambdaSameEvent || doprocessLambdaMixedEvent || doprocessK0shortSameEvent || doprocessK0shortMixedEvent; + bool processMc = doprocessLambdaSameEventMc || doprocessLambdaMixedEventMc || doprocessK0shortSameEventMc || doprocessK0shortMixedEventMc; + + if (processData && processMc) { + LOG(fatal) << "Both data and mc processing is enabled. Breaking..."; + } + + bool processLambda = doprocessLambdaSameEvent || doprocessLambdaSameEventMc || doprocessLambdaMixedEvent || doprocessLambdaMixedEventMc; + bool processK0short = doprocessK0shortSameEvent || doprocessK0shortSameEventMc || doprocessK0shortMixedEvent || doprocessK0shortMixedEventMc; + + if (processLambda && processK0short) { + LOG(fatal) << "Both lambda-track and k0short-track processing is enabled. Breaking..."; + } // setup columnpolicy for binning // default values are used during instantiation, so we need to explicity update them here @@ -140,55 +171,85 @@ struct FemtoPairTrackV0 { mixBinsVtxCent = {{confMixing.vtxBins.value, confMixing.centBins.value}, true}; mixBinsVtxMultCent = {{confMixing.vtxBins.value, confMixing.multBins.value, confMixing.centBins.value}, true}; - // setup histograms - auto colHistSpec = colhistmanager::makeColHistSpecMap(confCollisionBinning); - auto trackHistSpec = trackhistmanager::makeTrackHistSpecMap(confTrackBinning); - auto posDauSpec = trackhistmanager::makeTrackHistSpecMap(confPosDauBinning); - auto negDauSpec = trackhistmanager::makeTrackHistSpecMap(confNegDauBinning); auto cprHistSpec = closepairrejection::makeCprHistSpecMap(confCpr); - // setup for lambda - if (doprocessLambdaSameEvent || doprocessLambdaMixedEvent) { - auto lambdaHistSpec = v0histmanager::makeV0HistSpecMap(confLambdaBinning); - auto pairTrackLambdaHistSpec = pairhistmanager::makePairHistSpecMap(confPairBinning); - pairTrackLambdaBuilder.init(&hRegistry, trackSelection, lambdaSelection, confCpr, confMixing, confPairBinning, confPairCuts, colHistSpec, trackHistSpec, lambdaHistSpec, posDauSpec, negDauSpec, pairTrackLambdaHistSpec, cprHistSpec); - } - - // setup for k0short - if (doprocessK0shortSameEvent || doprocessK0shortMixedEvent) { - auto k0shortHistSpec = v0histmanager::makeV0HistSpecMap(confK0shortBinning); - auto pairTrackK0shortHistSpec = pairhistmanager::makePairHistSpecMap(confPairBinning); - pairTrackK0shortBuilder.init(&hRegistry, trackSelection, lambdaSelection, confCpr, confMixing, confPairBinning, confPairCuts, colHistSpec, trackHistSpec, k0shortHistSpec, posDauSpec, negDauSpec, pairTrackK0shortHistSpec, cprHistSpec); - } - - if (((doprocessLambdaSameEvent || doprocessLambdaMixedEvent) + (doprocessK0shortSameEvent || doprocessK0shortMixedEvent)) > 1) { - LOG(fatal) << "Can only process lambda-tracks Or k0short-tracks"; + if (processData) { + auto colHistSpec = colhistmanager::makeColHistSpecMap(confCollisionBinning); + auto trackHistSpec = trackhistmanager::makeTrackHistSpecMap(confTrackBinning); + auto posDauSpec = trackhistmanager::makeTrackHistSpecMap(confPosDauBinning); + auto negDauSpec = trackhistmanager::makeTrackHistSpecMap(confNegDauBinning); + auto pairTrackV0HistSpec = pairhistmanager::makePairHistSpecMap(confPairBinning); + if (processLambda) { + auto lambdaHistSpec = v0histmanager::makeV0HistSpecMap(confLambdaBinning); + pairTrackLambdaBuilder.init(&hRegistry, confTrackSelection, lambdaSelection, confCpr, confMixing, confPairBinning, confPairCuts, colHistSpec, trackHistSpec, lambdaHistSpec, posDauSpec, negDauSpec, pairTrackV0HistSpec, cprHistSpec); + } else { + auto k0shortHistSpec = v0histmanager::makeV0HistSpecMap(confK0shortBinning); + pairTrackK0shortBuilder.init(&hRegistry, confTrackSelection, lambdaSelection, confCpr, confMixing, confPairBinning, confPairCuts, colHistSpec, trackHistSpec, k0shortHistSpec, posDauSpec, negDauSpec, pairTrackV0HistSpec, cprHistSpec); + } + } else { + auto colHistSpec = colhistmanager::makeColMcHistSpecMap(confCollisionBinning); + auto trackHistSpec = trackhistmanager::makeTrackMcHistSpecMap(confTrackBinning); + auto posDauSpec = trackhistmanager::makeTrackMcHistSpecMap(confPosDauBinning); + auto negDauSpec = trackhistmanager::makeTrackMcHistSpecMap(confNegDauBinning); + auto pairTrackV0HistSpec = pairhistmanager::makePairMcHistSpecMap(confPairBinning); + if (processLambda) { + auto lambdaHistSpec = v0histmanager::makeV0McHistSpecMap(confLambdaBinning); + pairTrackLambdaBuilder.init(&hRegistry, confTrackSelection, lambdaSelection, confCpr, confMixing, confPairBinning, confPairCuts, colHistSpec, trackHistSpec, lambdaHistSpec, posDauSpec, negDauSpec, pairTrackV0HistSpec, cprHistSpec); + } else { + auto k0shortHistSpec = v0histmanager::makeV0McHistSpecMap(confK0shortBinning); + pairTrackK0shortBuilder.init(&hRegistry, confTrackSelection, lambdaSelection, confCpr, confMixing, confPairBinning, confPairCuts, colHistSpec, trackHistSpec, k0shortHistSpec, posDauSpec, negDauSpec, pairTrackV0HistSpec, cprHistSpec); + } } + hRegistry.print(); }; - void processLambdaSameEvent(FilteredCollision const& col, Tracks const& tracks, Lambdas const& lambdas) + void processLambdaSameEvent(FilteredFemtoCollision const& col, FemtoTracks const& tracks, FemtoLambdas const& lambdas) { pairTrackLambdaBuilder.processSameEvent(col, tracks, trackPartition, lambdas, lambdaPartition, cache); } PROCESS_SWITCH(FemtoPairTrackV0, processLambdaSameEvent, "Enable processing same event processing for tracks and lambdas", true); - void processLambdaMixedEvent(FilteredCollisions const& cols, Tracks const& tracks, Lambdas const& /*lambas*/) + void processLambdaSameEventMc(FilteredFemtoCollisionWithLabel const& col, FMcCols const& mcCols, FemtoTracksWithLabel const& tracks, FemtoLambdasWithLabel const& lambdas, FMcParticles const& mcParticles, FMcMothers const& mcMothers, FMcPartMoths const& mcPartonicMothers) + { + pairTrackLambdaBuilder.processSameEvent(col, mcCols, tracks, trackWithLabelPartition, lambdas, lambdaWithLabelPartition, mcParticles, mcMothers, mcPartonicMothers, cache); + } + PROCESS_SWITCH(FemtoPairTrackV0, processLambdaSameEventMc, "Enable processing same event processing for tracks and lambdas with MC information", false); + + void processLambdaMixedEvent(FilteredFemtoCollisions const& cols, FemtoTracks const& tracks, FemtoLambdas const& /*lambas*/) { pairTrackLambdaBuilder.processMixedEvent(cols, tracks, trackPartition, lambdaPartition, cache, mixBinsVtxMult, mixBinsVtxCent, mixBinsVtxMultCent); } PROCESS_SWITCH(FemtoPairTrackV0, processLambdaMixedEvent, "Enable processing mixed event processing for tracks and lambdas", true); - // - void processK0shortSameEvent(FilteredCollision const& col, Tracks const& tracks, K0shorts const& k0shorts) + + void processLambdaMixedEventMc(FilteredFemtoCollisionsWithLabel const& cols, FMcCols const& mcCols, FemtoTracksWithLabel const& tracks, FemtoLambdasWithLabel const& /*lambas*/, FMcParticles const& mcParticles) + { + pairTrackLambdaBuilder.processMixedEvent(cols, mcCols, tracks, trackWithLabelPartition, lambdaWithLabelPartition, mcParticles, cache, mixBinsVtxMult, mixBinsVtxCent, mixBinsVtxMultCent); + } + PROCESS_SWITCH(FemtoPairTrackV0, processLambdaMixedEventMc, "Enable processing mixed event processing for tracks and lambdas with MC information", false); + + void processK0shortSameEvent(FilteredFemtoCollision const& col, FemtoTracks const& tracks, FemtoK0shorts const& k0shorts) { pairTrackK0shortBuilder.processSameEvent(col, tracks, trackPartition, k0shorts, k0shortPartition, cache); } PROCESS_SWITCH(FemtoPairTrackV0, processK0shortSameEvent, "Enable processing same event processing for tracks and k0shorts", false); - void processK0shortMixedEvent(FilteredCollisions const& cols, Tracks const& tracks, K0shorts const& /*k0shorts*/) + void processK0shortSameEventMc(FilteredFemtoCollisionWithLabel const& col, FMcCols const& mcCols, FemtoTracksWithLabel const& tracks, FemtoK0shortsWithLabel const& k0shorts, FMcParticles const& mcParticles, FMcMothers const& mcMothers, FMcPartMoths const& mcPartonicMothers) + { + pairTrackK0shortBuilder.processSameEvent(col, mcCols, tracks, trackWithLabelPartition, k0shorts, k0shortWithLabelPartition, mcParticles, mcMothers, mcPartonicMothers, cache); + } + PROCESS_SWITCH(FemtoPairTrackV0, processK0shortSameEventMc, "Enable processing same event processing for tracks and k0shorts with MC information", false); + + void processK0shortMixedEvent(FilteredFemtoCollisions const& cols, FemtoTracks const& tracks, FemtoK0shorts const& /*k0shorts*/) { pairTrackK0shortBuilder.processMixedEvent(cols, tracks, trackPartition, k0shortPartition, cache, mixBinsVtxMult, mixBinsVtxCent, mixBinsVtxMultCent); } PROCESS_SWITCH(FemtoPairTrackV0, processK0shortMixedEvent, "Enable processing mixed event processing for tracks and k0shorts", false); + + void processK0shortMixedEventMc(FilteredFemtoCollisionsWithLabel const& cols, FMcCols const& mcCols, FemtoTracksWithLabel const& tracks, FemtoK0shortsWithLabel const& /*k0shorts*/, FMcParticles const& mcParticles) + { + pairTrackK0shortBuilder.processMixedEvent(cols, mcCols, tracks, trackWithLabelPartition, k0shortWithLabelPartition, mcParticles, cache, mixBinsVtxMult, mixBinsVtxCent, mixBinsVtxMultCent); + } + PROCESS_SWITCH(FemtoPairTrackV0, processK0shortMixedEventMc, "Enable processing mixed event processing for tracks and k0shorts with mc information", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) diff --git a/PWGCF/Femto/Tasks/femtoPairV0V0.cxx b/PWGCF/Femto/Tasks/femtoPairV0V0.cxx index c3280579349..b3cf25d3967 100644 --- a/PWGCF/Femto/Tasks/femtoPairV0V0.cxx +++ b/PWGCF/Femto/Tasks/femtoPairV0V0.cxx @@ -139,6 +139,10 @@ struct FemtoPairV0V0 { void init(InitContext&) { + if (((doprocessLambdaLambdaSameEvent || doprocessLambdaLambdaMixedEvent) + (doprocessK0shortK0shortSameEvent || doprocessK0shortK0shortMixedEvent)) > 1) { + LOG(fatal) << "Can only process lambda-tracks Or k0short-tracks"; + } + // setup columnpolicy for binning // default values are used during instantiation, so we need to explicity update them here mixBinsVtxMult = {{confMixing.vtxBins, confMixing.multBins.value}, true}; @@ -165,10 +169,6 @@ struct FemtoPairV0V0 { auto pairK0shortK0shortHistSpec = pairhistmanager::makePairHistSpecMap(confPairBinning); pairK0shortK0shortBuilder.init(&hRegistry, k0shortSelection, k0shortSelection, confCprPos, confCprNeg, confMixing, confPairBinning, confPairCuts, colHistSpec, k0shortHistSpec, k0shortHistSpec, posDauSpec, negDauSpec, pairK0shortK0shortHistSpec, cprHistSpecPos, cprHistSpecNeg); } - - if (((doprocessLambdaLambdaSameEvent || doprocessLambdaLambdaMixedEvent) + (doprocessK0shortK0shortSameEvent || doprocessK0shortK0shortMixedEvent)) > 1) { - LOG(fatal) << "Can only process lambda-tracks Or k0short-tracks"; - } }; void processLambdaLambdaSameEvent(FilteredCollision const& col, Tracks const& tracks, Lambdas const& lambdas) diff --git a/PWGCF/Femto/Tasks/femtoTrackQa.cxx b/PWGCF/Femto/Tasks/femtoTrackQa.cxx index 349296b7cd2..52d082ff7fa 100644 --- a/PWGCF/Femto/Tasks/femtoTrackQa.cxx +++ b/PWGCF/Femto/Tasks/femtoTrackQa.cxx @@ -30,7 +30,6 @@ #include "Framework/OutputObjHeader.h" #include "Framework/runDataProcessing.h" -#include #include using namespace o2::aod; @@ -63,41 +62,54 @@ struct FemtoTrackQa { colhistmanager::CollisionHistManager colHistManager; // setup tracks - trackbuilder::ConfTrackSelection1 trackSelections; + trackbuilder::ConfTrackSelection1 confTrackSelection; trackhistmanager::ConfTrackBinning1 confTrackBinning; trackhistmanager::ConfTrackQaBinning1 confTrackQaBinning; - trackhistmanager::ConfTrackMcBinning1 confTrackMcBinning; trackhistmanager::TrackHistManager trackHistManager; - Partition trackPartition = MAKE_TRACK_PARTITION(trackSelections); + Partition trackPartition = MAKE_TRACK_PARTITION(confTrackSelection); Preslice perColReco = femtobase::stored::fColId; - Partition trackWithLabelPartition = MAKE_TRACK_PARTITION(trackSelections); + Partition trackWithLabelPartition = MAKE_TRACK_PARTITION(confTrackSelection); Preslice perColRecoWithLabel = femtobase::stored::fColId; HistogramRegistry hRegistry{"FemtoTrackQA", {}, OutputObjHandlingPolicy::AnalysisObject}; - void init(InitContext&) + template + void initMode(MakeColSpec&& makeColSpec, + MakeTrackSpec&& makeTrackSpec) { + auto colHistSpec = makeColSpec(confCollisionBinning, confCollisionQaBinning); + colHistManager.init(&hRegistry, colHistSpec, confCollisionQaBinning); + auto trackHistSpec = makeTrackSpec(confTrackBinning, confTrackQaBinning); + + trackHistManager.init( + &hRegistry, + trackHistSpec, + confTrackSelection.chargeAbs.value, + confTrackSelection.chargeSign.value, + confTrackSelection.pdgCodeAbs.value, + confTrackQaBinning); + } + void init(InitContext&) + { if ((doprocessData + doprocessMc) > 1) { LOG(fatal) << "More than 1 process function is activated. Breaking..."; } - bool processData = doprocessData; - if (processData) { - // create a map for histogram specs auto colHistSpec = colhistmanager::makeColQaHistSpecMap(confCollisionBinning, confCollisionQaBinning); colHistManager.init(&hRegistry, colHistSpec, confCollisionQaBinning); auto trackHistSpec = trackhistmanager::makeTrackQaHistSpecMap(confTrackBinning, confTrackQaBinning); - trackHistManager.init(&hRegistry, trackHistSpec, confTrackQaBinning, trackSelections.chargeAbs.value); + trackHistManager.init(&hRegistry, trackHistSpec, confTrackSelection, confTrackQaBinning); } else { - // create a map for histogram specs auto colHistSpec = colhistmanager::makeColMcQaHistSpecMap(confCollisionBinning, confCollisionQaBinning); colHistManager.init(&hRegistry, colHistSpec, confCollisionQaBinning); - auto trackHistSpec = trackhistmanager::makeTrackMcQaHistSpecMap(confTrackBinning, confTrackQaBinning, confTrackMcBinning); - trackHistManager.init(&hRegistry, trackHistSpec, confTrackQaBinning, confTrackMcBinning, trackSelections.chargeAbs.value); + auto trackHistSpec = trackhistmanager::makeTrackMcQaHistSpecMap(confTrackBinning, confTrackQaBinning); + trackHistManager.init(&hRegistry, trackHistSpec, confTrackSelection, confTrackQaBinning); } hRegistry.print(); }; @@ -112,7 +124,6 @@ struct FemtoTrackQa { }; PROCESS_SWITCH(FemtoTrackQa, processData, "Track QA in Data", true); - // void processMc(FilteredFemtoCollisionWithLabel const& col, FemtoTracksWithLabel const& tracks, FMcCols const& mcCols,FMcMothers const& mcMothers, FMcPartMoths const& mcPartonicMothers) void processMc(FilteredFemtoCollisionWithLabel const& col, FMcCols const& mcCols, FemtoTracksWithLabel const& tracks, FMcParticles const& mcParticles, FMcMothers const& mcMothers, FMcPartMoths const& mcPartonicMothers) { colHistManager.fill(col, mcCols); diff --git a/PWGCF/Femto/Tasks/femtoTwotrackresonanceQa.cxx b/PWGCF/Femto/Tasks/femtoTwotrackresonanceQa.cxx index ec285f53c83..537286f1482 100644 --- a/PWGCF/Femto/Tasks/femtoTwotrackresonanceQa.cxx +++ b/PWGCF/Femto/Tasks/femtoTwotrackresonanceQa.cxx @@ -126,16 +126,16 @@ struct FemtoTwotrackresonanceQa { if (doprocessPhis) { auto phiHistSpec = twotrackresonancehistmanager::makeTwoTrackResonanceQaHistSpecMap(confPhiBinning); - phiHistManager.init(&hRegistry, phiHistSpec, posDaughterHistSpec, confPosDaughterQaBinning, negDaughterHistSpec, confNegDaughterQaBinning); + phiHistManager.init(&hRegistry, phiHistSpec, confPhiSelection, posDaughterHistSpec, confPosDaughterQaBinning, negDaughterHistSpec, confNegDaughterQaBinning); } if (doprocessRho0s) { auto rho0HistSpec = twotrackresonancehistmanager::makeTwoTrackResonanceQaHistSpecMap(confRho0Binning); - rho0HistManager.init(&hRegistry, rho0HistSpec, posDaughterHistSpec, confPosDaughterQaBinning, negDaughterHistSpec, confNegDaughterQaBinning); + rho0HistManager.init(&hRegistry, rho0HistSpec, confRho0Selection, posDaughterHistSpec, confPosDaughterQaBinning, negDaughterHistSpec, confNegDaughterQaBinning); } if (doprocessKstar0s) { auto kstar0HistSpec = twotrackresonancehistmanager::makeTwoTrackResonanceQaHistSpecMap(confKstar0Binning); - kstar0HistManager.init(&hRegistry, kstar0HistSpec, posDaughterHistSpec, confPosDaughterQaBinning, negDaughterHistSpec, confNegDaughterQaBinning); + kstar0HistManager.init(&hRegistry, kstar0HistSpec, confKstar0Selection, posDaughterHistSpec, confPosDaughterQaBinning, negDaughterHistSpec, confNegDaughterQaBinning); } }; diff --git a/PWGCF/Femto/Tasks/femtoV0Qa.cxx b/PWGCF/Femto/Tasks/femtoV0Qa.cxx index daeb05964de..1d1e0ec45f9 100644 --- a/PWGCF/Femto/Tasks/femtoV0Qa.cxx +++ b/PWGCF/Femto/Tasks/femtoV0Qa.cxx @@ -79,7 +79,6 @@ struct FemtoV0Qa { v0histmanager::ConfLambdaBinning1 confLambdaBinning; v0histmanager::ConfLambdaQaBinning1 confLambdaQaBinning; - v0histmanager::ConfLambdaMcBinning1 confLambdaMcBinning; v0histmanager::V0HistManager< v0histmanager::PrefixLambdaQa, trackhistmanager::PrefixV0PosDaughterQa, @@ -98,7 +97,6 @@ struct FemtoV0Qa { v0histmanager::ConfK0shortBinning1 confK0shortBinning; v0histmanager::ConfK0shortQaBinning1 confK0shortQaBinning; - v0histmanager::ConfK0shortMcBinning1 confK0shortMcBinning; v0histmanager::V0HistManager< v0histmanager::PrefixK0shortQa, trackhistmanager::PrefixV0PosDaughterQa, @@ -109,11 +107,9 @@ struct FemtoV0Qa { // setup for daughters trackhistmanager::ConfV0PosDauBinning confV0PosDaughterBinning; trackhistmanager::ConfV0PosDauQaBinning confV0PosDaughterQaBinning; - trackhistmanager::ConfV0PosDauMcBinning confV0PosDaughterMcBinning; trackhistmanager::ConfV0NegDauBinning confV0NegDaughterBinning; trackhistmanager::ConfV0NegDauQaBinning confV0NegDaughterQaBinning; - trackhistmanager::ConfV0NegDauMcBinning confV0NegDaughterMcBinning; HistogramRegistry hRegistry{"FemtoV0Qa", {}, OutputObjHandlingPolicy::AnalysisObject}; @@ -128,35 +124,28 @@ struct FemtoV0Qa { if (processData) { auto colHistSpec = colhistmanager::makeColQaHistSpecMap(confCollisionBinning, confCollisionQaBinning); colHistManager.init(&hRegistry, colHistSpec, confCollisionQaBinning); - auto posDaughterHistSpec = trackhistmanager::makeTrackQaHistSpecMap(confV0PosDaughterBinning, confV0PosDaughterQaBinning); auto negDaughterHistSpec = trackhistmanager::makeTrackQaHistSpecMap(confV0NegDaughterBinning, confV0NegDaughterQaBinning); - if (doprocessLambda) { auto lambdaHistSpec = v0histmanager::makeV0QaHistSpecMap(confLambdaBinning, confLambdaQaBinning); - lambdaHistManager.init(&hRegistry, lambdaHistSpec, confLambdaQaBinning, posDaughterHistSpec, confV0PosDaughterQaBinning, negDaughterHistSpec, confV0NegDaughterQaBinning); + lambdaHistManager.init(&hRegistry, lambdaHistSpec, confLambdaSelection, confLambdaQaBinning, posDaughterHistSpec, confV0PosDaughterQaBinning, negDaughterHistSpec, confV0NegDaughterQaBinning); } - if (doprocessK0short) { auto k0shortHistSpec = v0histmanager::makeV0QaHistSpecMap(confK0shortBinning, confK0shortQaBinning); - k0shortHistManager.init(&hRegistry, k0shortHistSpec, confK0shortQaBinning, posDaughterHistSpec, confV0PosDaughterQaBinning, negDaughterHistSpec, confV0NegDaughterQaBinning); + k0shortHistManager.init(&hRegistry, k0shortHistSpec, confK0shortSelection, confK0shortQaBinning, posDaughterHistSpec, confV0PosDaughterQaBinning, negDaughterHistSpec, confV0NegDaughterQaBinning); } } else { - // process mc auto colHistSpec = colhistmanager::makeColMcQaHistSpecMap(confCollisionBinning, confCollisionQaBinning); colHistManager.init(&hRegistry, colHistSpec, confCollisionQaBinning); - - auto posDaughterHistSpec = trackhistmanager::makeTrackMcQaHistSpecMap(confV0PosDaughterBinning, confV0PosDaughterQaBinning, confV0PosDaughterMcBinning); - auto negDaughterHistSpec = trackhistmanager::makeTrackMcQaHistSpecMap(confV0NegDaughterBinning, confV0NegDaughterQaBinning, confV0NegDaughterMcBinning); - + auto posDaughterHistSpec = trackhistmanager::makeTrackMcQaHistSpecMap(confV0PosDaughterBinning, confV0PosDaughterQaBinning); + auto negDaughterHistSpec = trackhistmanager::makeTrackMcQaHistSpecMap(confV0NegDaughterBinning, confV0NegDaughterQaBinning); if (doprocessLambdaMc) { - auto lambdaHistSpec = v0histmanager::makeV0McQaHistSpecMap(confLambdaBinning, confLambdaQaBinning, confLambdaMcBinning); - lambdaHistManager.init(&hRegistry, lambdaHistSpec, confLambdaQaBinning, confLambdaMcBinning, posDaughterHistSpec, confV0PosDaughterQaBinning, confV0PosDaughterMcBinning, negDaughterHistSpec, confV0NegDaughterQaBinning, confV0NegDaughterMcBinning); + auto lambdaHistSpec = v0histmanager::makeV0McQaHistSpecMap(confLambdaBinning, confLambdaQaBinning); + lambdaHistManager.init(&hRegistry, lambdaHistSpec, confLambdaSelection, confLambdaQaBinning, posDaughterHistSpec, confV0PosDaughterQaBinning, negDaughterHistSpec, confV0NegDaughterQaBinning); } - if (doprocessK0shortMc) { - auto k0shortHistSpec = v0histmanager::makeV0McQaHistSpecMap(confK0shortBinning, confK0shortQaBinning, confK0shortMcBinning); - k0shortHistManager.init(&hRegistry, k0shortHistSpec, confK0shortQaBinning, confK0shortMcBinning, posDaughterHistSpec, confV0PosDaughterQaBinning, confV0PosDaughterMcBinning, negDaughterHistSpec, confV0NegDaughterQaBinning, confV0NegDaughterMcBinning); + auto k0shortHistSpec = v0histmanager::makeV0McQaHistSpecMap(confK0shortBinning, confK0shortQaBinning); + k0shortHistManager.init(&hRegistry, k0shortHistSpec, confK0shortSelection, confK0shortQaBinning, posDaughterHistSpec, confV0PosDaughterQaBinning, negDaughterHistSpec, confV0NegDaughterQaBinning); } } hRegistry.print();