diff --git a/ALICE3/Tasks/alice3-dq-efficiency.cxx b/ALICE3/Tasks/alice3-dq-efficiency.cxx index 99dcd281727..ca4c117cabe 100644 --- a/ALICE3/Tasks/alice3-dq-efficiency.cxx +++ b/ALICE3/Tasks/alice3-dq-efficiency.cxx @@ -41,8 +41,12 @@ #include #include +#include #include +#include #include +#include +#include #include #include #include @@ -133,8 +137,15 @@ constexpr static uint32_t gkEventFillMapWithCov = VarManager::ObjTypes::ReducedE constexpr static uint32_t gkTrackFillMapWithCov = VarManager::ObjTypes::ReducedTrack | VarManager::ObjTypes::ReducedTrackBarrel | VarManager::ObjTypes::ReducedTrackBarrelCov | VarManager::ObjTypes::ReducedTrackBarrelPID; constexpr static uint32_t gkTrackFillMap = VarManager::ObjTypes::ReducedTrack | VarManager::ObjTypes::ReducedTrackBarrel | VarManager::ObjTypes::ReducedTrackBarrelPID; +namespace dqefficiency_helpers +{ +inline float* varValues() { return static_cast(VarManager::fgValues); } +inline TString* varNames() { return static_cast(VarManager::fgVariableNames); } +inline TString* varUnits() { return static_cast(VarManager::fgVariableUnits); } +} // namespace dqefficiency_helpers + // Global function used to define needed histogram classes -void DefineHistograms(HistogramManager* histMan, TString histClasses, const char* histGroups); // defines histograms for all tasks +void DefineHistograms(HistogramManager* histMan, const TString& histClasses, const char* histGroups); // defines histograms for all tasks constexpr int TWO_PRONG = 2; constexpr int THREE_PRONG = 3; @@ -155,7 +166,7 @@ struct AnalysisEventSelection { HistogramManager* fHistMan = nullptr; MixingHandler* fMixHandler = nullptr; - AnalysisCompositeCut* fEventCut; + AnalysisCompositeCut* fEventCut = nullptr; std::map fSelMap; // key: reduced event global index, value: event selection decision @@ -188,7 +199,7 @@ struct AnalysisEventSelection { if (fConfigQA) { fHistMan = new HistogramManager("analysisHistos", "", VarManager::kNVars); fHistMan->SetUseDefaultVariableNames(true); - fHistMan->SetDefaultVarNames(VarManager::fgVariableNames, VarManager::fgVariableUnits); + fHistMan->SetDefaultVarNames(dqefficiency_helpers::varNames(), dqefficiency_helpers::varUnits()); DefineHistograms(fHistMan, "Event_BeforeCuts;Event_AfterCuts;", fConfigAddEventHistogram.value.data()); DefineHistograms(fHistMan, "EventsMC", fConfigAddEventMCHistogram.value.data()); dqhistograms::AddHistogramsFromJSON(fHistMan, fConfigAddJSONHistograms.value.c_str()); // aditional histograms via JSON @@ -222,17 +233,17 @@ struct AnalysisEventSelection { bool decision = false; // if QA is requested fill histograms before event selections if (fConfigQA) { - fHistMan->FillHistClass("Event_BeforeCuts", VarManager::fgValues); // automatically fill all the histograms in the class Event + fHistMan->FillHistClass("Event_BeforeCuts", dqefficiency_helpers::varValues()); // automatically fill all the histograms in the class Event } - if (fEventCut->IsSelected(VarManager::fgValues)) { + if (fEventCut->IsSelected(dqefficiency_helpers::varValues())) { if (fConfigQA) { - fHistMan->FillHistClass("Event_AfterCuts", VarManager::fgValues); + fHistMan->FillHistClass("Event_AfterCuts", dqefficiency_helpers::varValues()); } decision = true; } fSelMap[event.globalIndex()] = decision; if (fMixHandler != nullptr) { - int hh = fMixHandler->FindEventCategory(VarManager::fgValues); + int hh = fMixHandler->FindEventCategory(dqefficiency_helpers::varValues()); hash(hh); } } @@ -242,7 +253,7 @@ struct AnalysisEventSelection { VarManager::ResetValues(0, VarManager::kNEventWiseVariables); VarManager::FillEventAlice3(event); if (fConfigQA) { - fHistMan->FillHistClass("EventsMC", VarManager::fgValues); + fHistMan->FillHistClass("EventsMC", dqefficiency_helpers::varValues()); } } } @@ -250,7 +261,7 @@ struct AnalysisEventSelection { void publishSelections(MyEventsVtxCov const& events) { // publish the table - uint32_t evSel = static_cast(0); + auto evSel = static_cast(0); for (const auto& event : events) { evSel = 0; if (fSelMap[event.globalIndex()]) { // event passed the user cuts @@ -292,7 +303,7 @@ struct AnalysisTrackSelection { Configurable fConfigMCSignals{"cfgTrackMCSignals", "", "Comma separated list of MC signals"}; Configurable fConfigMCSignalsJSON{"cfgTrackMCsignalsJSON", "", "Additional list of MC signals via JSON"}; - HistogramManager* fHistMan; + HistogramManager* fHistMan = nullptr; std::vector fTrackCuts; std::vector fMCSignals; // list of signals to be checked std::vector fHistNamesReco; @@ -320,7 +331,7 @@ struct AnalysisTrackSelection { if (addTrackCutsStr != "") { std::vector addTrackCuts = dqcuts::GetCutsFromJSON(addTrackCutsStr.Data()); for (const auto& t : addTrackCuts) { - fTrackCuts.push_back(reinterpret_cast(t)); + fTrackCuts.push_back(dynamic_cast(t)); } } VarManager::SetUseVars(AnalysisCut::fgUsedVars); // provide the list of required variables so that VarManager knows what to fill @@ -352,7 +363,7 @@ struct AnalysisTrackSelection { if (fConfigQA) { fHistMan = new HistogramManager("analysisHistos", "aa", VarManager::kNVars); fHistMan->SetUseDefaultVariableNames(true); - fHistMan->SetDefaultVarNames(VarManager::fgVariableNames, VarManager::fgVariableUnits); + fHistMan->SetDefaultVarNames(dqefficiency_helpers::varNames(), dqefficiency_helpers::varUnits()); // Configure histogram classes for each track cut; // Add histogram classes for each track cut and for each requested MC signal (reconstructed tracks with MC truth) @@ -421,16 +432,16 @@ struct AnalysisTrackSelection { } if (fConfigQA) { - fHistMan->FillHistClass("AssocsBarrel_BeforeCuts", VarManager::fgValues); + fHistMan->FillHistClass("AssocsBarrel_BeforeCuts", dqefficiency_helpers::varValues()); } int iCut = 0; - uint32_t filterMap = static_cast(0); + auto filterMap = static_cast(0); for (auto cut = fTrackCuts.begin(); cut != fTrackCuts.end(); cut++, iCut++) { - if ((*cut)->IsSelected(VarManager::fgValues)) { + if ((*cut)->IsSelected(dqefficiency_helpers::varValues())) { filterMap |= (static_cast(1) << iCut); if (fConfigQA) { - fHistMan->FillHistClass(fHistNamesReco[iCut], VarManager::fgValues); + fHistMan->FillHistClass(fHistNamesReco[iCut], dqefficiency_helpers::varValues()); } } } // end loop over cuts @@ -446,11 +457,11 @@ struct AnalysisTrackSelection { // mcDecision |= (static_cast(1) << isig); // loop over cuts and fill histograms for the cuts that are fulfilled for (unsigned int icut = 0; icut < fTrackCuts.size(); icut++) { - if (filterMap & (static_cast(1) << icut)) { + if ((filterMap & (static_cast(1) << icut)) != 0u) { if (isCorrectAssoc) { - fHistMan->FillHistClass(fHistNamesMCMatched[icut * 2 * fMCSignals.size() + 2 * isig].Data(), VarManager::fgValues); + fHistMan->FillHistClass(fHistNamesMCMatched[icut * 2 * fMCSignals.size() + 2 * isig].Data(), dqefficiency_helpers::varValues()); } else { - fHistMan->FillHistClass(fHistNamesMCMatched[icut * 2 * fMCSignals.size() + 2 * isig + 1].Data(), VarManager::fgValues); + fHistMan->FillHistClass(fHistNamesMCMatched[icut * 2 * fMCSignals.size() + 2 * isig + 1].Data(), dqefficiency_helpers::varValues()); } } } // end loop over cuts @@ -462,7 +473,7 @@ struct AnalysisTrackSelection { if (fConfigPublishAmbiguity && filterMap > 0) { if (event.isEventSelected_bit(1)) { // for this track, count the number of associated collisions with in-bunch pileup and out of bunch associations - if (fNAssocsInBunch.find(track.globalIndex()) == fNAssocsInBunch.end()) { + if (!fNAssocsInBunch.contains(track.globalIndex())) { std::vector evVector = {event.globalIndex()}; fNAssocsInBunch[track.globalIndex()] = evVector; } else { @@ -470,7 +481,7 @@ struct AnalysisTrackSelection { evVector.push_back(event.globalIndex()); } } else { - if (fNAssocsOutOfBunch.find(track.globalIndex()) == fNAssocsOutOfBunch.end()) { + if (!fNAssocsOutOfBunch.contains(track.globalIndex())) { std::vector evVector = {event.globalIndex()}; fNAssocsOutOfBunch[track.globalIndex()] = evVector; } else { @@ -494,7 +505,7 @@ struct AnalysisTrackSelection { VarManager::ResetValues(0, VarManager::kNBarrelTrackVariables); VarManager::FillTrackAlice3(track); VarManager::fgValues[VarManager::kBarrelNAssocsInBunch] = static_cast(evIndices.size()); - fHistMan->FillHistClass("TrackBarrel_AmbiguityInBunch", VarManager::fgValues); + fHistMan->FillHistClass("TrackBarrel_AmbiguityInBunch", dqefficiency_helpers::varValues()); } // end loop over in-bunch ambiguous tracks for (const auto& [trackIdx, evIndices] : fNAssocsOutOfBunch) { @@ -505,18 +516,18 @@ struct AnalysisTrackSelection { VarManager::ResetValues(0, VarManager::kNBarrelTrackVariables); VarManager::FillTrackAlice3(track); VarManager::fgValues[VarManager::kBarrelNAssocsOutOfBunch] = static_cast(evIndices.size()); - fHistMan->FillHistClass("TrackBarrel_AmbiguityOutOfBunch", VarManager::fgValues); + fHistMan->FillHistClass("TrackBarrel_AmbiguityOutOfBunch", dqefficiency_helpers::varValues()); } // end loop over out-of-bunch ambiguous tracks } // publish the ambiguity table for (const auto& track : tracks) { int8_t nInBunch = 0; - if (fNAssocsInBunch.find(track.globalIndex()) != fNAssocsInBunch.end()) { + if (!fNAssocsInBunch.contains(track.globalIndex())) { nInBunch = fNAssocsInBunch[track.globalIndex()].size(); } int8_t nOutOfBunch = 0; - if (fNAssocsOutOfBunch.find(track.globalIndex()) != fNAssocsOutOfBunch.end()) { + if (!fNAssocsOutOfBunch.contains(track.globalIndex())) { nOutOfBunch = fNAssocsOutOfBunch[track.globalIndex()].size(); } trackAmbiguities(nInBunch, nOutOfBunch); @@ -549,9 +560,9 @@ struct AnalysisPrefilterSelection { Configurable fPropTrack{"cfgPropTrack", false, "Propgate tracks to associated collision to recalculate DCA and momentum vector"}; std::map fPrefilterMap; - AnalysisCompositeCut* fPairCut; - uint32_t fPrefilterMask; - int fPrefilterCutBit; + AnalysisCompositeCut* fPairCut = nullptr; + uint32_t fPrefilterMask = 0; + int fPrefilterCutBit = -1; PresliceUnsorted trackAssocsPerCollision = aod::reducedA3track_association::reducedA3eventId; @@ -655,7 +666,7 @@ struct AnalysisPrefilterSelection { bool track1Loose = assoc1.isBarrelSelected_bit(fPrefilterCutBit); bool track2Loose = assoc2.isBarrelSelected_bit(fPrefilterCutBit); - if (!((track1Candidate > 0 && track2Loose) || (track2Candidate > 0 && track1Loose))) { + if ((track1Candidate == 0 || !track2Loose) && (track2Candidate == 0 || !track1Loose)) { continue; } @@ -665,11 +676,11 @@ struct AnalysisPrefilterSelection { VarManager::FillPairCollision(event, track1, track2); } // if the pair fullfils the criteria, add an entry into the prefilter map for the two tracks - if (fPairCut->IsSelected(VarManager::fgValues)) { - if (fPrefilterMap.find(track1.globalIndex()) == fPrefilterMap.end() && track1Candidate > 0) { + if (fPairCut->IsSelected(dqefficiency_helpers::varValues())) { + if (!fPrefilterMap.contains(track1.globalIndex()) && track1Candidate > 0) { fPrefilterMap[track1.globalIndex()] = track1Candidate; } - if (fPrefilterMap.find(track2.globalIndex()) == fPrefilterMap.end() && track2Candidate > 0) { + if (!fPrefilterMap.contains(track2.globalIndex()) && track2Candidate > 0) { fPrefilterMap[track2.globalIndex()] = track2Candidate; } } @@ -699,7 +710,7 @@ struct AnalysisPrefilterSelection { // TODO: just use the index from the assoc (no need to cast the whole track) auto track = assoc.template reducedA3track_as(); mymap = -1; - if (fPrefilterMap.find(track.globalIndex()) != fPrefilterMap.end()) { + if (!fPrefilterMap.contains(track.globalIndex())) { // NOTE: publish the bitwise negated bits (~), so there will be zeroes for cuts that failed the prefiltering and 1 everywhere else mymap = ~fPrefilterMap[track.globalIndex()]; prefilter(mymap); @@ -765,7 +776,7 @@ struct AnalysisSameEventPairing { // Filter filterEventSelected = aod::dqanalysisflags::isEventSelected & uint32_t(1); Filter eventFilter = aod::dqanalysisflags::isEventSelected > static_cast(0); - HistogramManager* fHistMan; + HistogramManager* fHistMan = nullptr; // keep histogram class names in maps, so we don't have to buld their names in the pair loops std::map> fTrackHistNames; @@ -777,12 +788,12 @@ struct AnalysisSameEventPairing { AnalysisCompositeCut fMCGenAccCut; bool fUseMCGenAccCut = false; - uint32_t fTrackFilterMask; // mask for the track cuts required in this task to be applied on the barrel cuts produced upstream - int fNCutsBarrel; - int fNPairCuts; + uint32_t fTrackFilterMask = 0; // mask for the track cuts required in this task to be applied on the barrel cuts produced upstream + int fNCutsBarrel = 0; + int fNPairCuts = 0; bool fHasTwoProngGenMCsignals = false; - bool fEnableBarrelHistos; + bool fEnableBarrelHistos = false; PresliceUnsorted trackAssocsPerCollision = aod::reducedA3track_association::reducedA3eventId; @@ -898,9 +909,9 @@ struct AnalysisSameEventPairing { // if there are pair cuts specified, assign hist directories for each barrel cut - pair cut combination // NOTE: This could possibly lead to large histogram outputs. It is strongly advised to use pair cuts only // if you know what you are doing. - TString cutNamesStr = fConfigCuts.pair.value; - if (!cutNamesStr.IsNull()) { // if pair cuts - std::unique_ptr objArrayPair(cutNamesStr.Tokenize(",")); + TString pairCutNamesStr = fConfigCuts.pair.value; + if (!pairCutNamesStr.IsNull()) { // if pair cuts + std::unique_ptr objArrayPair(pairCutNamesStr.Tokenize(",")); fNPairCuts = objArrayPair->GetEntries(); for (int iPairCut = 0; iPairCut < fNPairCuts; ++iPairCut) { // loop over pair cuts names = { @@ -980,7 +991,7 @@ struct AnalysisSameEventPairing { fHistMan = new HistogramManager("analysisHistos", "aa", VarManager::kNVars); fHistMan->SetUseDefaultVariableNames(true); - fHistMan->SetDefaultVarNames(VarManager::fgVariableNames, VarManager::fgVariableUnits); + fHistMan->SetDefaultVarNames(dqefficiency_helpers::varNames(), dqefficiency_helpers::varUnits()); VarManager::SetCollisionSystem((TString)fConfigOptions.collisionSystem, fConfigOptions.centerMassEnergy); // set collision system and center of mass energy @@ -1003,17 +1014,32 @@ struct AnalysisSameEventPairing { std::map> histNamesMC = fBarrelHistNamesMCmatched; int ncuts = fNCutsBarrel; - uint32_t twoTrackFilter = static_cast(0); + auto twoTrackFilter = static_cast(0); int sign1 = 0; int sign2 = 0; - uint32_t mcDecision = static_cast(0); + auto mcDecision = static_cast(0); bool isCorrectAssoc_leg1 = false; bool isCorrectAssoc_leg2 = false; - dielectronList.reserve(1); - dielectronsExtraList.reserve(1); + + int64_t reserveSize = 0; + for (auto const& event : events) { + if (event.isEventSelected_bit(0)) { + auto groupedAssocs = assocs.sliceBy(preslice, event.globalIndex()); + size_t nGood = 0; + for (auto const& t : groupedAssocs) { + if (t.isBarrelSelected_raw() > 0u) { + nGood++; + } + } + reserveSize += nGood * (nGood - 1) / 2; + } + } + + dielectronList.reserve(reserveSize); + dielectronsExtraList.reserve(reserveSize); if (fConfigOptions.flatTables.value) { - dielectronAllList.reserve(1); + dielectronAllList.reserve(reserveSize); } for (const auto& event : events) { @@ -1022,8 +1048,8 @@ struct AnalysisSameEventPairing { } // Reset the fValues array VarManager::ResetValues(0, VarManager::kNVars); - VarManager::FillEventAlice3(event, VarManager::fgValues); - VarManager::FillEventAlice3(event.reducedA3MCEvent(), VarManager::fgValues); + VarManager::FillEventAlice3(event, dqefficiency_helpers::varValues()); + VarManager::FillEventAlice3(event.reducedA3MCEvent(), dqefficiency_helpers::varValues()); auto groupedAssocs = assocs.sliceBy(preslice, event.globalIndex()); if (groupedAssocs.size() == 0) { @@ -1034,7 +1060,7 @@ struct AnalysisSameEventPairing { twoTrackFilter = a1.isBarrelSelected_raw() & a2.isBarrelSelected_raw() & a1.isBarrelSelectedPrefilter_raw() & a2.isBarrelSelectedPrefilter_raw() & fTrackFilterMask; - if (!twoTrackFilter) { // the tracks must have at least one filter bit in common to continue + if (twoTrackFilter == 0u) { // the tracks must have at least one filter bit in common to continue continue; } @@ -1057,12 +1083,12 @@ struct AnalysisSameEventPairing { } // run MC matching for this pair - int isig = 0; + int iSigMc = 0; mcDecision = 0; - for (auto sig = fRecMCSignals.begin(); sig != fRecMCSignals.end(); sig++, isig++) { + for (auto sig = fRecMCSignals.begin(); sig != fRecMCSignals.end(); sig++, iSigMc++) { if (t1.has_reducedA3MCTrack() && t2.has_reducedA3MCTrack()) { if ((*sig)->CheckSignal(true, t1.reducedA3MCTrack(), t2.reducedA3MCTrack())) { - mcDecision |= (static_cast(1) << isig); + mcDecision |= (static_cast(1) << iSigMc); } } } // end loop over MC signals @@ -1075,11 +1101,9 @@ struct AnalysisSameEventPairing { if (fPropTrack) { VarManager::FillPairCollision(event, t1, t2); } - /* TODO: Reimplement Pair vertexing when secondary vertexing is available - if constexpr (TTwoProngFitter) { - // VarManager::FillPairVertexing(event, t1, t2, fConfigOptions.propToPCA); - }*/ - if (!fConfigMC.skimSignalOnly || (fConfigMC.skimSignalOnly && mcDecision > 0)) { + + VarManager::FillPairVertexingAlice3(event, t1, t2, true); + if (!fConfigMC.skimSignalOnly || mcDecision > 0) { dielectronList(event.globalIndex(), VarManager::fgValues[VarManager::kMass], VarManager::fgValues[VarManager::kPt], VarManager::fgValues[VarManager::kEta], VarManager::fgValues[VarManager::kPhi], t1.sign() + t2.sign(), twoTrackFilter, mcDecision); @@ -1090,91 +1114,92 @@ struct AnalysisSameEventPairing { bool isAmbiOutOfBunch = false; for (int icut = 0; icut < ncuts; icut++) { - if (twoTrackFilter & (static_cast(1) << icut)) { - isAmbiInBunch = (twoTrackFilter & (static_cast(1) << 28)) || (twoTrackFilter & (static_cast(1) << 29)); - isAmbiOutOfBunch = (twoTrackFilter & (static_cast(1) << 30)) || (twoTrackFilter & (static_cast(1) << 31)); - if (sign1 * sign2 < 0) { // +- pairs - fHistMan->FillHistClass(histNames[icut][0].Data(), VarManager::fgValues); // reconstructed, unmatched - for (unsigned int isig = 0; isig < fRecMCSignals.size(); isig++) { // loop over MC signals - if (mcDecision & (static_cast(1) << isig)) { - fHistMan->FillHistClass(histNamesMC[icut * fRecMCSignals.size() + isig][0].Data(), VarManager::fgValues); // matched signal + if ((twoTrackFilter & (static_cast(1) << icut)) != 0u) { + isAmbiInBunch = ((twoTrackFilter & (static_cast(1) << 28)) != 0u) || ((twoTrackFilter & (static_cast(1) << 29)) != 0u); + isAmbiOutOfBunch = ((twoTrackFilter & (static_cast(1) << 30)) != 0u) || ((twoTrackFilter & (static_cast(1) << 31)) != 0u); + if (sign1 * sign2 < 0) { // +- pairs + fHistMan->FillHistClass(histNames[icut][0].Data(), dqefficiency_helpers::varValues()); // reconstructed, unmatched + for (unsigned int isig = 0; isig < fRecMCSignals.size(); isig++) { // loop over MC signals + if ((mcDecision & (static_cast(1) << isig)) != 0u) { + fHistMan->FillHistClass(histNamesMC[icut * fRecMCSignals.size() + isig][0].Data(), dqefficiency_helpers::varValues()); // matched signal if (fConfigQA) { if (isCorrectAssoc_leg1 && isCorrectAssoc_leg2) { // correct track-collision association - fHistMan->FillHistClass(histNamesMC[icut * fRecMCSignals.size() + isig][3].Data(), VarManager::fgValues); + fHistMan->FillHistClass(histNamesMC[icut * fRecMCSignals.size() + isig][3].Data(), dqefficiency_helpers::varValues()); } else { // incorrect track-collision association - fHistMan->FillHistClass(histNamesMC[icut * fRecMCSignals.size() + isig][4].Data(), VarManager::fgValues); + fHistMan->FillHistClass(histNamesMC[icut * fRecMCSignals.size() + isig][4].Data(), dqefficiency_helpers::varValues()); } if (isAmbiInBunch) { // ambiguous in bunch - fHistMan->FillHistClass(histNamesMC[icut * fRecMCSignals.size() + isig][5].Data(), VarManager::fgValues); + fHistMan->FillHistClass(histNamesMC[icut * fRecMCSignals.size() + isig][5].Data(), dqefficiency_helpers::varValues()); if (isCorrectAssoc_leg1 && isCorrectAssoc_leg2) { - fHistMan->FillHistClass(histNamesMC[icut * fRecMCSignals.size() + isig][6].Data(), VarManager::fgValues); + fHistMan->FillHistClass(histNamesMC[icut * fRecMCSignals.size() + isig][6].Data(), dqefficiency_helpers::varValues()); } else { - fHistMan->FillHistClass(histNamesMC[icut * fRecMCSignals.size() + isig][7].Data(), VarManager::fgValues); + fHistMan->FillHistClass(histNamesMC[icut * fRecMCSignals.size() + isig][7].Data(), dqefficiency_helpers::varValues()); } } if (isAmbiOutOfBunch) { // ambiguous out of bunch - fHistMan->FillHistClass(histNamesMC[icut * fRecMCSignals.size() + isig][8].Data(), VarManager::fgValues); + fHistMan->FillHistClass(histNamesMC[icut * fRecMCSignals.size() + isig][8].Data(), dqefficiency_helpers::varValues()); if (isCorrectAssoc_leg1 && isCorrectAssoc_leg2) { - fHistMan->FillHistClass(histNamesMC[icut * fRecMCSignals.size() + isig][9].Data(), VarManager::fgValues); + fHistMan->FillHistClass(histNamesMC[icut * fRecMCSignals.size() + isig][9].Data(), dqefficiency_helpers::varValues()); } else { - fHistMan->FillHistClass(histNamesMC[icut * fRecMCSignals.size() + isig][10].Data(), VarManager::fgValues); + fHistMan->FillHistClass(histNamesMC[icut * fRecMCSignals.size() + isig][10].Data(), dqefficiency_helpers::varValues()); } } } } if (fConfigQA) { if (isAmbiInBunch) { - fHistMan->FillHistClass(histNames[icut][3].Data(), VarManager::fgValues); + fHistMan->FillHistClass(histNames[icut][3].Data(), dqefficiency_helpers::varValues()); } if (isAmbiOutOfBunch) { - fHistMan->FillHistClass(histNames[icut][3 + 3].Data(), VarManager::fgValues); + fHistMan->FillHistClass(histNames[icut][3 + 3].Data(), dqefficiency_helpers::varValues()); } } } } else { if (sign1 > 0) { // ++ pairs - fHistMan->FillHistClass(histNames[icut][1].Data(), VarManager::fgValues); + fHistMan->FillHistClass(histNames[icut][1].Data(), dqefficiency_helpers::varValues()); for (unsigned int isig = 0; isig < fRecMCSignals.size(); isig++) { // loop over MC signals - if (mcDecision & (static_cast(1) << isig)) { - fHistMan->FillHistClass(histNamesMC[icut * fRecMCSignals.size() + isig][1].Data(), VarManager::fgValues); + if ((mcDecision & (static_cast(1) << isig)) != 0u) { + fHistMan->FillHistClass(histNamesMC[icut * fRecMCSignals.size() + isig][1].Data(), dqefficiency_helpers::varValues()); } } if (fConfigQA) { if (isAmbiInBunch) { - fHistMan->FillHistClass(histNames[icut][4].Data(), VarManager::fgValues); + fHistMan->FillHistClass(histNames[icut][4].Data(), dqefficiency_helpers::varValues()); } if (isAmbiOutOfBunch) { - fHistMan->FillHistClass(histNames[icut][4 + 3].Data(), VarManager::fgValues); + fHistMan->FillHistClass(histNames[icut][4 + 3].Data(), dqefficiency_helpers::varValues()); } } } else { // -- pairs - fHistMan->FillHistClass(histNames[icut][2].Data(), VarManager::fgValues); + fHistMan->FillHistClass(histNames[icut][2].Data(), dqefficiency_helpers::varValues()); for (unsigned int isig = 0; isig < fRecMCSignals.size(); isig++) { // loop over MC signals - if (mcDecision & (static_cast(1) << isig)) { - fHistMan->FillHistClass(histNamesMC[icut * fRecMCSignals.size() + isig][2].Data(), VarManager::fgValues); + if ((mcDecision & (static_cast(1) << isig)) != 0u) { + fHistMan->FillHistClass(histNamesMC[icut * fRecMCSignals.size() + isig][2].Data(), dqefficiency_helpers::varValues()); } } if (fConfigQA) { if (isAmbiInBunch) { - fHistMan->FillHistClass(histNames[icut][5].Data(), VarManager::fgValues); + fHistMan->FillHistClass(histNames[icut][5].Data(), dqefficiency_helpers::varValues()); } if (isAmbiOutOfBunch) { - fHistMan->FillHistClass(histNames[icut][5 + 3].Data(), VarManager::fgValues); + fHistMan->FillHistClass(histNames[icut][5 + 3].Data(), dqefficiency_helpers::varValues()); } } } } for (unsigned int iPairCut = 0; iPairCut < fPairCuts.size(); iPairCut++) { AnalysisCompositeCut cut = fPairCuts.at(iPairCut); - if (!(cut.IsSelected(VarManager::fgValues))) // apply pair cuts + if (!(cut.IsSelected(dqefficiency_helpers::varValues()))) { // apply pair cuts continue; + } if (sign1 * sign2 < 0) { - fHistMan->FillHistClass(histNames[ncuts + icut * ncuts + iPairCut][0].Data(), VarManager::fgValues); + fHistMan->FillHistClass(histNames[ncuts + icut * ncuts + iPairCut][0].Data(), dqefficiency_helpers::varValues()); } else { if (sign1 > 0) { - fHistMan->FillHistClass(histNames[ncuts + icut * ncuts + iPairCut][1].Data(), VarManager::fgValues); + fHistMan->FillHistClass(histNames[ncuts + icut * ncuts + iPairCut][1].Data(), dqefficiency_helpers::varValues()); } else { - fHistMan->FillHistClass(histNames[ncuts + icut * ncuts + iPairCut][2].Data(), VarManager::fgValues); + fHistMan->FillHistClass(histNames[ncuts + icut * ncuts + iPairCut][2].Data(), dqefficiency_helpers::varValues()); } } } // end loop (pair cuts) @@ -1188,14 +1213,11 @@ struct AnalysisSameEventPairing { void runMCGenWithGrouping(MyEventsVtxCovSelected const& events, ReducedA3MCEvents const& /*mcEvents*/, ReducedA3MCTracks const& mcTracks) { - [[maybe_unused]] uint32_t mcDecision = 0; - int isig = 0; - for (const auto& mctrack : mcTracks) { VarManager::FillTrackMC(mcTracks, mctrack); // if we have a mc generated acceptance cut, apply it here if (fUseMCGenAccCut) { - if (!fMCGenAccCut.IsSelected(VarManager::fgValues)) { + if (!fMCGenAccCut.IsSelected(dqefficiency_helpers::varValues())) { continue; } } @@ -1204,7 +1226,7 @@ struct AnalysisSameEventPairing { // TODO: Use the mcReducedFlags to select signals for (const auto& sig : fGenMCSignals) { if (sig->CheckSignal(true, mctrack)) { - fHistMan->FillHistClass(Form("MCTruthGen_%s", sig->GetName()), VarManager::fgValues); + fHistMan->FillHistClass(Form("MCTruthGen_%s", sig->GetName()), dqefficiency_helpers::varValues()); } } } @@ -1224,20 +1246,16 @@ struct AnalysisSameEventPairing { VarManager::FillTrackMC(mcTracks, track); // if we have a mc generated acceptance cut, apply it here if (fUseMCGenAccCut) { - if (!fMCGenAccCut.IsSelected(VarManager::fgValues)) { + if (!fMCGenAccCut.IsSelected(dqefficiency_helpers::varValues())) { continue; } } auto track_raw = mcTracks.rawIteratorAt(track.globalIndex()); - mcDecision = 0; - isig = 0; for (const auto& sig : fGenMCSignals) { if (sig->CheckSignal(true, track_raw)) { - mcDecision |= (static_cast(1) << isig); - fHistMan->FillHistClass(Form("MCTruthGenSel_%s", sig->GetName()), VarManager::fgValues); + fHistMan->FillHistClass(Form("MCTruthGenSel_%s", sig->GetName()), dqefficiency_helpers::varValues()); MCTruthTableEffi(VarManager::fgValues[VarManager::kMCPt], VarManager::fgValues[VarManager::kMCEta], VarManager::fgValues[VarManager::kMCY], VarManager::fgValues[VarManager::kMCPhi], VarManager::fgValues[VarManager::kMCVz], VarManager::fgValues[VarManager::kMCVtxZ], VarManager::fgValues[VarManager::kMultFT0A], VarManager::fgValues[VarManager::kMultFT0C], VarManager::fgValues[VarManager::kCentFT0M], VarManager::fgValues[VarManager::kVtxNcontribReal]); } - isig++; } } } // end loop over reconstructed events @@ -1254,11 +1272,11 @@ struct AnalysisSameEventPairing { if (sig->CheckSignal(true, t1_raw, t2_raw)) { VarManager::FillPairMC(t1, t2); if (fUseMCGenAccCut) { - if (!fMCGenAccCut.IsSelected(VarManager::fgValues)) { + if (!fMCGenAccCut.IsSelected(dqefficiency_helpers::varValues())) { continue; } } - fHistMan->FillHistClass(Form("MCTruthGenPair_%s", sig->GetName()), VarManager::fgValues); + fHistMan->FillHistClass(Form("MCTruthGenPair_%s", sig->GetName()), dqefficiency_helpers::varValues()); } } } @@ -1279,23 +1297,19 @@ struct AnalysisSameEventPairing { auto t1_raw = mcTracks.rawIteratorAt(t1.globalIndex()); auto t2_raw = mcTracks.rawIteratorAt(t2.globalIndex()); if (t1_raw.reducedA3MCEventId() == t2_raw.reducedA3MCEventId()) { - mcDecision = 0; - isig = 0; for (const auto& sig : fGenMCSignals) { if (sig->GetNProngs() != TWO_PRONG) { // NOTE: 2-prong signals required here continue; } if (sig->CheckSignal(true, t1_raw, t2_raw)) { - mcDecision |= (static_cast(1) << isig); VarManager::FillPairMC(t1, t2); if (fUseMCGenAccCut) { - if (!fMCGenAccCut.IsSelected(VarManager::fgValues)) { + if (!fMCGenAccCut.IsSelected(dqefficiency_helpers::varValues())) { continue; } } - fHistMan->FillHistClass(Form("MCTruthGenPairSel_%s", sig->GetName()), VarManager::fgValues); + fHistMan->FillHistClass(Form("MCTruthGenPairSel_%s", sig->GetName()), dqefficiency_helpers::varValues()); } - isig++; } } } @@ -1315,10 +1329,6 @@ struct AnalysisSameEventPairing { void processMCGen(soa::Filtered const& events, ReducedA3MCEvents const& /*mcEvents*/, ReducedA3MCTracks const& mcTracks) { - // Fill Generated histograms taking into account all generated tracks - [[maybe_unused]] uint32_t mcDecision = 0; - int isig = 0; - for (const auto& mctrack : mcTracks) { VarManager::FillTrackMC(mcTracks, mctrack); // NOTE: Signals are checked here mostly based on the skimmed MC stack, so depending on the requested signal, the stack could be incomplete. @@ -1326,7 +1336,7 @@ struct AnalysisSameEventPairing { // TODO: Use the mcReducedFlags to select signals for (const auto& sig : fGenMCSignals) { if (sig->CheckSignal(true, mctrack)) { - fHistMan->FillHistClass(Form("MCTruthGen_%s", sig->GetName()), VarManager::fgValues); + fHistMan->FillHistClass(Form("MCTruthGen_%s", sig->GetName()), dqefficiency_helpers::varValues()); } } } @@ -1339,27 +1349,20 @@ struct AnalysisSameEventPairing { if (!event.has_reducedA3MCEvent()) { continue; } - VarManager::FillEventAlice3(event, VarManager::fgValues); - VarManager::FillEventAlice3(event.reducedA3MCEvent(), VarManager::fgValues); - // auto groupedMCTracks = mcTracks.sliceBy(perReducedMcGenEvent, event.reducedA3MCEventId()); - // groupedMCTracks.bindInternalIndicesTo(&mcTracks); - // for (const auto& track : groupedMCTracks) { + VarManager::FillEventAlice3(event, dqefficiency_helpers::varValues()); + VarManager::FillEventAlice3(event.reducedA3MCEvent(), dqefficiency_helpers::varValues()); + for (const auto& track : mcTracks) { if (track.reducedA3MCEventId() != event.reducedA3MCEventId()) { continue; } VarManager::FillTrackMC(mcTracks, track); auto track_raw = mcTracks.rawIteratorAt(track.globalIndex()); - // auto track_raw = groupedMCTracks.rawIteratorAt(track.globalIndex()); - mcDecision = 0; - isig = 0; for (const auto& sig : fGenMCSignals) { if (sig->CheckSignal(true, track_raw)) { - mcDecision |= (static_cast(1) << isig); - fHistMan->FillHistClass(Form("MCTruthGenSel_%s", sig->GetName()), VarManager::fgValues); + fHistMan->FillHistClass(Form("MCTruthGenSel_%s", sig->GetName()), dqefficiency_helpers::varValues()); MCTruthTableEffi(VarManager::fgValues[VarManager::kMCPt], VarManager::fgValues[VarManager::kMCEta], VarManager::fgValues[VarManager::kMCY], VarManager::fgValues[VarManager::kMCPhi], VarManager::fgValues[VarManager::kMCVz], VarManager::fgValues[VarManager::kMCVtxZ], VarManager::fgValues[VarManager::kMultFT0A], VarManager::fgValues[VarManager::kMultFT0C], VarManager::fgValues[VarManager::kCentFT0M], VarManager::fgValues[VarManager::kVtxNcontribReal]); } - isig++; } } } // end loop over reconstructed events @@ -1373,7 +1376,7 @@ struct AnalysisSameEventPairing { continue; } if (sig->CheckSignal(true, t1_raw, t2_raw)) { - fHistMan->FillHistClass(Form("MCTruthGenPair_%s", sig->GetName()), VarManager::fgValues); + fHistMan->FillHistClass(Form("MCTruthGenPair_%s", sig->GetName()), dqefficiency_helpers::varValues()); } } } @@ -1399,17 +1402,13 @@ struct AnalysisSameEventPairing { auto t1_raw = mcTracks.rawIteratorAt(t1.globalIndex()); auto t2_raw = mcTracks.rawIteratorAt(t2.globalIndex()); if (t1_raw.reducedA3MCEventId() == t2_raw.reducedA3MCEventId()) { - mcDecision = 0; - isig = 0; for (const auto& sig : fGenMCSignals) { if (sig->GetNProngs() != TWO_PRONG) { // NOTE: 2-prong signals required here continue; } if (sig->CheckSignal(true, t1_raw, t2_raw)) { - mcDecision |= (static_cast(1) << isig); - fHistMan->FillHistClass(Form("MCTruthGenPairSel_%s", sig->GetName()), VarManager::fgValues); + fHistMan->FillHistClass(Form("MCTruthGenPairSel_%s", sig->GetName()), dqefficiency_helpers::varValues()); } - isig++; } } } @@ -1419,9 +1418,6 @@ struct AnalysisSameEventPairing { void processMCGenWithGrouping(soa::Filtered const& events, ReducedA3MCEvents const& /*mcEvents*/, ReducedA3MCTracks const& mcTracks) { - [[maybe_unused]] uint32_t mcDecision = 0; - int isig = 0; - for (const auto& mctrack : mcTracks) { VarManager::FillTrackMC(mcTracks, mctrack); // NOTE: Signals are checked here mostly based on the skimmed MC stack, so depending on the requested signal, the stack could be incomplete. @@ -1429,7 +1425,7 @@ struct AnalysisSameEventPairing { // TODO: Use the mcReducedFlags to select signals for (const auto& sig : fGenMCSignals) { if (sig->CheckSignal(true, mctrack)) { - fHistMan->FillHistClass(Form("MCTruthGen_%s", sig->GetName()), VarManager::fgValues); + fHistMan->FillHistClass(Form("MCTruthGen_%s", sig->GetName()), dqefficiency_helpers::varValues()); } } } @@ -1448,15 +1444,11 @@ struct AnalysisSameEventPairing { } VarManager::FillTrackMC(mcTracks, track); auto track_raw = mcTracks.rawIteratorAt(track.globalIndex()); - mcDecision = 0; - isig = 0; for (const auto& sig : fGenMCSignals) { if (sig->CheckSignal(true, track_raw)) { - mcDecision |= (static_cast(1) << isig); - fHistMan->FillHistClass(Form("MCTruthGenSel_%s", sig->GetName()), VarManager::fgValues); + fHistMan->FillHistClass(Form("MCTruthGenSel_%s", sig->GetName()), dqefficiency_helpers::varValues()); MCTruthTableEffi(VarManager::fgValues[VarManager::kMCPt], VarManager::fgValues[VarManager::kMCEta], VarManager::fgValues[VarManager::kMCY], VarManager::fgValues[VarManager::kMCPhi], VarManager::fgValues[VarManager::kMCVz], VarManager::fgValues[VarManager::kMCVtxZ], VarManager::fgValues[VarManager::kMultFT0A], VarManager::fgValues[VarManager::kMultFT0C], VarManager::fgValues[VarManager::kCentFT0M], VarManager::fgValues[VarManager::kVtxNcontribReal]); } - isig++; } } } // end loop over reconstructed events @@ -1470,7 +1462,7 @@ struct AnalysisSameEventPairing { continue; } if (sig->CheckSignal(true, t1_raw, t2_raw)) { - fHistMan->FillHistClass(Form("MCTruthGenPair_%s", sig->GetName()), VarManager::fgValues); + fHistMan->FillHistClass(Form("MCTruthGenPair_%s", sig->GetName()), dqefficiency_helpers::varValues()); } } } @@ -1491,17 +1483,13 @@ struct AnalysisSameEventPairing { auto t1_raw = groupedMCTracks.rawIteratorAt(t1.globalIndex()); auto t2_raw = groupedMCTracks.rawIteratorAt(t2.globalIndex()); if (t1_raw.reducedA3MCEventId() == t2_raw.reducedA3MCEventId()) { - mcDecision = 0; - isig = 0; for (const auto& sig : fGenMCSignals) { if (sig->GetNProngs() != TWO_PRONG) { // NOTE: 2-prong signals required here continue; } if (sig->CheckSignal(true, t1_raw, t2_raw)) { - mcDecision |= (static_cast(1) << isig); - fHistMan->FillHistClass(Form("MCTruthGenPairSel_%s", sig->GetName()), VarManager::fgValues); + fHistMan->FillHistClass(Form("MCTruthGenPairSel_%s", sig->GetName()), dqefficiency_helpers::varValues()); } - isig++; } } } @@ -1549,30 +1537,30 @@ struct AnalysisAsymmetricPairing { Configurable fConfigMCRecSignalsJSON{"cfgMCRecSignalsJSON", "", "Additional list of MC signals (reconstructed) via JSON"}; Configurable fConfigMCGenSignalsJSON{"cfgMCGenSignalsJSON", "", "Comma separated list of MC signals (generated) via JSON"}; - HistogramManager* fHistMan; + HistogramManager* fHistMan = nullptr; std::vector fPairCuts; - int fNPairHistPrefixes; + int fNPairHistPrefixes = 0; std::vector fRecMCSignals; std::vector fGenMCSignals; // Filter masks to find legs in BarrelTrackCuts table - uint32_t fLegAFilterMask; - uint32_t fLegBFilterMask; - uint32_t fLegCFilterMask; + uint32_t fLegAFilterMask = 0; + uint32_t fLegBFilterMask = 0; + uint32_t fLegCFilterMask = 0; // Maps tracking which combination of leg cuts the track cuts participate in std::map fConstructedLegAFilterMasksMap; std::map fConstructedLegBFilterMasksMap; std::map fConstructedLegCFilterMasksMap; // Filter map for common track cuts - uint32_t fCommonTrackCutMask; + uint32_t fCommonTrackCutMask = 0; // Map tracking which common track cut the track cuts correspond to std::map fCommonTrackCutFilterMasks; - int fNLegCuts; + int fNLegCuts = 0; int fNPairCuts = 0; - int fNCommonTrackCuts; + int fNCommonTrackCuts = 0; // vectors for cut names and signal names, for easy access when calling FillHistogramList() std::vector fLegCutNames; std::vector fPairCutNames; @@ -1602,7 +1590,7 @@ struct AnalysisAsymmetricPairing { VarManager::SetDefaultVarNames(); fHistMan = new HistogramManager("analysisHistos", "aa", VarManager::kNVars); fHistMan->SetUseDefaultVariableNames(true); - fHistMan->SetDefaultVarNames(VarManager::fgVariableNames, VarManager::fgVariableUnits); + fHistMan->SetDefaultVarNames(dqefficiency_helpers::varNames(), dqefficiency_helpers::varUnits()); // Get the leg cut filter masks fLegAFilterMask = fConfigLegAFilterMask.value; @@ -1622,7 +1610,7 @@ struct AnalysisAsymmetricPairing { if (addPairCutsStr != "") { std::vector addPairCuts = dqcuts::GetCutsFromJSON(addPairCutsStr.Data()); for (const auto& t : addPairCuts) { - fPairCuts.push_back(reinterpret_cast(t)); + fPairCuts.push_back(dynamic_cast(t)); cutNamesStr += Form(",%s", t->GetName()); } } @@ -1674,7 +1662,7 @@ struct AnalysisAsymmetricPairing { } std::unique_ptr objArray(tempCutsStr.Tokenize(",")); // Get the common leg cuts - int commonCutIdx; + int commonCutIdx = -1; TString commonNamesStr = fConfigCommonTrackCuts.value; if (!commonNamesStr.IsNull()) { // if common track cuts std::unique_ptr objArrayCommon(commonNamesStr.Tokenize(",")); @@ -1712,9 +1700,9 @@ struct AnalysisAsymmetricPairing { } fNLegCuts = objArrayLegs->GetEntries(); std::vector isThreeProng; - int legAIdx; - int legBIdx; - int legCIdx; + int legAIdx = -1; + int legBIdx = -1; + int legCIdx = -1; // Loop over leg defining cuts for (int icut = 0; icut < fNLegCuts; ++icut) { TString legsStr = objArrayLegs->At(icut)->GetName(); @@ -1949,7 +1937,7 @@ struct AnalysisAsymmetricPairing { } // Reset the fValues array VarManager::ResetValues(0, VarManager::kNVars); - VarManager::FillEventAlice3(event, VarManager::fgValues); + VarManager::FillEventAlice3(event, dqefficiency_helpers::varValues()); auto groupedLegAAssocs = legACandidateAssocs.sliceBy(preslice, event.globalIndex()); if (groupedLegAAssocs.size() == 0) { @@ -1961,25 +1949,24 @@ struct AnalysisAsymmetricPairing { } for (const auto& [a1, a2] : combinations(soa::CombinationsFullIndexPolicy(groupedLegAAssocs, groupedLegBAssocs))) { - uint32_t twoTrackFilter = 0; uint32_t twoTrackCommonFilter = 0; uint32_t pairFilter = 0; bool isPairIdWrong = false; for (int icut = 0; icut < fNLegCuts; ++icut) { // Find leg pair definitions both candidates participate in - if ((a1.isBarrelSelected_raw() & fConstructedLegAFilterMasksMap[icut]) && (a2.isBarrelSelected_raw() & fConstructedLegBFilterMasksMap[icut])) { + if (((a1.isBarrelSelected_raw() & fConstructedLegAFilterMasksMap[icut]) != 0u) && ((a2.isBarrelSelected_raw() & fConstructedLegBFilterMasksMap[icut]) != 0u)) { twoTrackFilter |= static_cast(1) << icut; // If the supposed pion passes a kaon cut, this is a K+K-. Skip it. if (fConfigSkipAmbiguousIdCombinations.value) { - if (a2.isBarrelSelected_raw() & fLegAFilterMask) { + if ((a2.isBarrelSelected_raw() & fLegAFilterMask) != 0u) { isPairIdWrong = true; } } } } - if (!twoTrackFilter || isPairIdWrong) { + if (twoTrackFilter == 0u || isPairIdWrong) { continue; } @@ -1996,12 +1983,12 @@ struct AnalysisAsymmetricPairing { bool isReflected = false; std::pair trackIds(t1.globalIndex(), t2.globalIndex()); - if (fPairCount.find(trackIds) != fPairCount.end()) { + if (fPairCount.contains(trackIds)) { // Double counting is possible due to track-collision ambiguity. Skip pairs which were counted before fPairCount[trackIds] += 1; continue; } - if (fPairCount.find(std::pair(trackIds.second, trackIds.first)) != fPairCount.end()) { + if (fPairCount.contains(std::pair(trackIds.second, trackIds.first))) { isReflected = true; } fPairCount[trackIds] += 1; @@ -2017,97 +2004,94 @@ struct AnalysisAsymmetricPairing { } // run MC matching for this pair - int isig = 0; + int iSigMc = 0; mcDecision = 0; - for (auto sig = fRecMCSignals.begin(); sig != fRecMCSignals.end(); sig++, isig++) { + for (auto sig = fRecMCSignals.begin(); sig != fRecMCSignals.end(); sig++, iSigMc++) { if (t1.has_reducedA3MCTrack() && t2.has_reducedA3MCTrack()) { VarManager::FillPairMC(t1.reducedA3MCTrack(), t2.reducedA3MCTrack()); if ((*sig)->CheckSignal(true, t1.reducedA3MCTrack(), t2.reducedA3MCTrack())) { - mcDecision |= static_cast(1) << isig; + mcDecision |= static_cast(1) << iSigMc; } } } // end loop over MC signals VarManager::FillPairAlice3(t1, t2); - /*TODO: Reimplement when secondary vertexing is available - if constexpr (TTwoProngFitter) { - VarManager::FillPairVertexing(event, t1, t2, fConfigPropToPCA); - }*/ + VarManager::FillPairVertexingAlice3(event, t1, t2, true); // Fill histograms bool isAmbi = false; for (int icut = 0; icut < fNLegCuts; icut++) { - if (twoTrackFilter & (static_cast(1) << icut)) { - isAmbi = (twoTrackFilter & (static_cast(1) << 30)) || (twoTrackFilter & (static_cast(1) << 31)); - if (sign1 * sign2 < 0) { // +- pairs - fHistMan->FillHistClass(Form("PairsBarrelSEPM_%s", fLegCutNames[icut].Data()), VarManager::fgValues); // reconstructed, unmatched + if ((twoTrackFilter & (static_cast(1) << icut)) != 0u) { + isAmbi = ((twoTrackFilter & (static_cast(1) << 30)) != 0u) || ((twoTrackFilter & (static_cast(1) << 31)) != 0u); + if (sign1 * sign2 < 0) { // +- pairs + fHistMan->FillHistClass(Form("PairsBarrelSEPM_%s", fLegCutNames[icut].Data()), dqefficiency_helpers::varValues()); // reconstructed, unmatched if (isAmbi && fConfigQA) { - fHistMan->FillHistClass(Form("PairsBarrelSEPM_ambiguous_%s", fLegCutNames[icut].Data()), VarManager::fgValues); + fHistMan->FillHistClass(Form("PairsBarrelSEPM_ambiguous_%s", fLegCutNames[icut].Data()), dqefficiency_helpers::varValues()); } if (isReflected && fConfigReflectedHistograms.value) { - fHistMan->FillHistClass(Form("PairsBarrelSEPM_reflected_%s", fLegCutNames[icut].Data()), VarManager::fgValues); + fHistMan->FillHistClass(Form("PairsBarrelSEPM_reflected_%s", fLegCutNames[icut].Data()), dqefficiency_helpers::varValues()); } } else if (fConfigSameSignHistograms.value) { if (sign1 > 0) { // ++ pairs - fHistMan->FillHistClass(Form("PairsBarrelSEPP_%s", fLegCutNames[icut].Data()), VarManager::fgValues); + fHistMan->FillHistClass(Form("PairsBarrelSEPP_%s", fLegCutNames[icut].Data()), dqefficiency_helpers::varValues()); if (isAmbi && fConfigQA) { - fHistMan->FillHistClass(Form("PairsBarrelSEPP_ambiguous_%s", fLegCutNames[icut].Data()), VarManager::fgValues); + fHistMan->FillHistClass(Form("PairsBarrelSEPP_ambiguous_%s", fLegCutNames[icut].Data()), dqefficiency_helpers::varValues()); } if (isReflected && fConfigReflectedHistograms.value) { - fHistMan->FillHistClass(Form("PairsBarrelSEPP_reflected_%s", fLegCutNames[icut].Data()), VarManager::fgValues); + fHistMan->FillHistClass(Form("PairsBarrelSEPP_reflected_%s", fLegCutNames[icut].Data()), dqefficiency_helpers::varValues()); } } else { // -- pairs - fHistMan->FillHistClass(Form("PairsBarrelSEMM_%s", fLegCutNames[icut].Data()), VarManager::fgValues); + fHistMan->FillHistClass(Form("PairsBarrelSEMM_%s", fLegCutNames[icut].Data()), dqefficiency_helpers::varValues()); if (isAmbi && fConfigQA) { - fHistMan->FillHistClass(Form("PairsBarrelSEMM_ambiguous_%s", fLegCutNames[icut].Data()), VarManager::fgValues); + fHistMan->FillHistClass(Form("PairsBarrelSEMM_ambiguous_%s", fLegCutNames[icut].Data()), dqefficiency_helpers::varValues()); } if (isReflected && fConfigReflectedHistograms) { - fHistMan->FillHistClass(Form("PairsBarrelSEMM_reflected_%s", fLegCutNames[icut].Data()), VarManager::fgValues); + fHistMan->FillHistClass(Form("PairsBarrelSEMM_reflected_%s", fLegCutNames[icut].Data()), dqefficiency_helpers::varValues()); } } } for (unsigned int isig = 0; isig < fRecMCSignals.size(); isig++) { // loop over MC signals - if (mcDecision & (static_cast(1) << isig)) { + if ((mcDecision & (static_cast(1) << isig)) != 0u) { if (sign1 * sign2 < 0) { - fHistMan->FillHistClass(Form("PairsBarrelSEPM_%s_%s", fLegCutNames[icut].Data(), fRecMCSignalNames[isig].Data()), VarManager::fgValues); + fHistMan->FillHistClass(Form("PairsBarrelSEPM_%s_%s", fLegCutNames[icut].Data(), fRecMCSignalNames[isig].Data()), dqefficiency_helpers::varValues()); if (isReflected && fConfigReflectedHistograms.value) { - fHistMan->FillHistClass(Form("PairsBarrelSEPM_reflected_%s_%s", fLegCutNames[icut].Data(), fRecMCSignalNames[isig].Data()), VarManager::fgValues); + fHistMan->FillHistClass(Form("PairsBarrelSEPM_reflected_%s_%s", fLegCutNames[icut].Data(), fRecMCSignalNames[isig].Data()), dqefficiency_helpers::varValues()); } } else if (fConfigSameSignHistograms.value) { if (sign1 > 0) { - fHistMan->FillHistClass(Form("PairsBarrelSEPP_%s_%s", fLegCutNames[icut].Data(), fRecMCSignalNames[isig].Data()), VarManager::fgValues); + fHistMan->FillHistClass(Form("PairsBarrelSEPP_%s_%s", fLegCutNames[icut].Data(), fRecMCSignalNames[isig].Data()), dqefficiency_helpers::varValues()); if (isReflected && fConfigReflectedHistograms.value) { - fHistMan->FillHistClass(Form("PairsBarrelSEPP_reflected_%s_%s", fLegCutNames[icut].Data(), fRecMCSignalNames[isig].Data()), VarManager::fgValues); + fHistMan->FillHistClass(Form("PairsBarrelSEPP_reflected_%s_%s", fLegCutNames[icut].Data(), fRecMCSignalNames[isig].Data()), dqefficiency_helpers::varValues()); } } else { - fHistMan->FillHistClass(Form("PairsBarrelSEMM_%s_%s", fLegCutNames[icut].Data(), fRecMCSignalNames[isig].Data()), VarManager::fgValues); + fHistMan->FillHistClass(Form("PairsBarrelSEMM_%s_%s", fLegCutNames[icut].Data(), fRecMCSignalNames[isig].Data()), dqefficiency_helpers::varValues()); if (isReflected && fConfigReflectedHistograms.value) { - fHistMan->FillHistClass(Form("PairsBarrelSEMM_reflected_%s_%s", fLegCutNames[icut].Data(), fRecMCSignalNames[isig].Data()), VarManager::fgValues); + fHistMan->FillHistClass(Form("PairsBarrelSEMM_reflected_%s_%s", fLegCutNames[icut].Data(), fRecMCSignalNames[isig].Data()), dqefficiency_helpers::varValues()); } } } } } for (int iCommonCut = 0; iCommonCut < fNCommonTrackCuts; iCommonCut++) { - if (twoTrackCommonFilter & fCommonTrackCutFilterMasks[iCommonCut]) { + if ((twoTrackCommonFilter & fCommonTrackCutFilterMasks[iCommonCut]) != 0u) { if (sign1 * sign2 < 0) { - fHistMan->FillHistClass(Form("PairsBarrelSEPM_%s_%s", fLegCutNames[icut].Data(), fCommonCutNames[iCommonCut].Data()), VarManager::fgValues); + fHistMan->FillHistClass(Form("PairsBarrelSEPM_%s_%s", fLegCutNames[icut].Data(), fCommonCutNames[iCommonCut].Data()), dqefficiency_helpers::varValues()); } else if (fConfigSameSignHistograms.value) { if (sign1 > 0) { - fHistMan->FillHistClass(Form("PairsBarrelSEPP_%s_%s", fLegCutNames[icut].Data(), fCommonCutNames[iCommonCut].Data()), VarManager::fgValues); + fHistMan->FillHistClass(Form("PairsBarrelSEPP_%s_%s", fLegCutNames[icut].Data(), fCommonCutNames[iCommonCut].Data()), dqefficiency_helpers::varValues()); } else { - fHistMan->FillHistClass(Form("PairsBarrelSEMM_%s_%s", fLegCutNames[icut].Data(), fCommonCutNames[iCommonCut].Data()), VarManager::fgValues); + fHistMan->FillHistClass(Form("PairsBarrelSEMM_%s_%s", fLegCutNames[icut].Data(), fCommonCutNames[iCommonCut].Data()), dqefficiency_helpers::varValues()); } } for (unsigned int isig = 0; isig < fRecMCSignals.size(); isig++) { // loop over MC signals - if (mcDecision & (static_cast(1) << isig)) { + if ((mcDecision & (static_cast(1) << isig)) != 0u) { if (sign1 * sign2 < 0) { - fHistMan->FillHistClass(Form("PairsBarrelSEPM_%s_%s_%s", fLegCutNames[icut].Data(), fCommonCutNames[iCommonCut].Data(), fRecMCSignalNames[isig].Data()), VarManager::fgValues); + fHistMan->FillHistClass(Form("PairsBarrelSEPM_%s_%s_%s", fLegCutNames[icut].Data(), fCommonCutNames[iCommonCut].Data(), fRecMCSignalNames[isig].Data()), dqefficiency_helpers::varValues()); } else if (fConfigSameSignHistograms.value) { if (sign1 > 0) { - fHistMan->FillHistClass(Form("PairsBarrelSEPP_%s_%s_%s", fLegCutNames[icut].Data(), fCommonCutNames[iCommonCut].Data(), fRecMCSignalNames[isig].Data()), VarManager::fgValues); + fHistMan->FillHistClass(Form("PairsBarrelSEPP_%s_%s_%s", fLegCutNames[icut].Data(), fCommonCutNames[iCommonCut].Data(), fRecMCSignalNames[isig].Data()), dqefficiency_helpers::varValues()); } else { - fHistMan->FillHistClass(Form("PairsBarrelSEMM_%s_%s_%s", fLegCutNames[icut].Data(), fCommonCutNames[iCommonCut].Data(), fRecMCSignalNames[isig].Data()), VarManager::fgValues); + fHistMan->FillHistClass(Form("PairsBarrelSEMM_%s_%s_%s", fLegCutNames[icut].Data(), fCommonCutNames[iCommonCut].Data(), fRecMCSignalNames[isig].Data()), dqefficiency_helpers::varValues()); } } } @@ -2116,53 +2100,54 @@ struct AnalysisAsymmetricPairing { } // end loop (common cuts) int iPairCut = 0; for (auto cut = fPairCuts.begin(); cut != fPairCuts.end(); cut++, iPairCut++) { - if (!((*cut)->IsSelected(VarManager::fgValues))) // apply pair cuts + if (!((*cut)->IsSelected(dqefficiency_helpers::varValues()))) { // apply pair cuts continue; + } pairFilter |= (static_cast(1) << iPairCut); // Histograms with pair cuts if (sign1 * sign2 < 0) { - fHistMan->FillHistClass(Form("PairsBarrelSEPM_%s_%s", fLegCutNames[icut].Data(), fPairCutNames[iPairCut].Data()), VarManager::fgValues); + fHistMan->FillHistClass(Form("PairsBarrelSEPM_%s_%s", fLegCutNames[icut].Data(), fPairCutNames[iPairCut].Data()), dqefficiency_helpers::varValues()); } else if (fConfigSameSignHistograms.value) { if (sign1 > 0) { - fHistMan->FillHistClass(Form("PairsBarrelSEPP_%s_%s", fLegCutNames[icut].Data(), fPairCutNames[iPairCut].Data()), VarManager::fgValues); + fHistMan->FillHistClass(Form("PairsBarrelSEPP_%s_%s", fLegCutNames[icut].Data(), fPairCutNames[iPairCut].Data()), dqefficiency_helpers::varValues()); } else { - fHistMan->FillHistClass(Form("PairsBarrelSEMM_%s_%s", fLegCutNames[icut].Data(), fPairCutNames[iPairCut].Data()), VarManager::fgValues); + fHistMan->FillHistClass(Form("PairsBarrelSEMM_%s_%s", fLegCutNames[icut].Data(), fPairCutNames[iPairCut].Data()), dqefficiency_helpers::varValues()); } } for (unsigned int isig = 0; isig < fRecMCSignals.size(); isig++) { // loop over MC signals - if (mcDecision & (static_cast(1) << isig)) { + if ((mcDecision & (static_cast(1) << isig)) != 0u) { if (sign1 * sign2 < 0) { - fHistMan->FillHistClass(Form("PairsBarrelSEPM_%s_%s_%s", fLegCutNames[icut].Data(), fPairCutNames[iPairCut].Data(), fRecMCSignalNames[isig].Data()), VarManager::fgValues); + fHistMan->FillHistClass(Form("PairsBarrelSEPM_%s_%s_%s", fLegCutNames[icut].Data(), fPairCutNames[iPairCut].Data(), fRecMCSignalNames[isig].Data()), dqefficiency_helpers::varValues()); } else if (fConfigSameSignHistograms.value) { if (sign1 > 0) { - fHistMan->FillHistClass(Form("PairsBarrelSEPP_%s_%s_%s", fLegCutNames[icut].Data(), fPairCutNames[iPairCut].Data(), fRecMCSignalNames[isig].Data()), VarManager::fgValues); + fHistMan->FillHistClass(Form("PairsBarrelSEPP_%s_%s_%s", fLegCutNames[icut].Data(), fPairCutNames[iPairCut].Data(), fRecMCSignalNames[isig].Data()), dqefficiency_helpers::varValues()); } else { - fHistMan->FillHistClass(Form("PairsBarrelSEMM_%s_%s_%s", fLegCutNames[icut].Data(), fPairCutNames[iPairCut].Data(), fRecMCSignalNames[isig].Data()), VarManager::fgValues); + fHistMan->FillHistClass(Form("PairsBarrelSEMM_%s_%s_%s", fLegCutNames[icut].Data(), fPairCutNames[iPairCut].Data(), fRecMCSignalNames[isig].Data()), dqefficiency_helpers::varValues()); } } } } // Histograms with pair cuts and common track cuts for (int iCommonCut = 0; iCommonCut < fNCommonTrackCuts; ++iCommonCut) { - if (twoTrackCommonFilter & fCommonTrackCutFilterMasks[iCommonCut]) { + if ((twoTrackCommonFilter & fCommonTrackCutFilterMasks[iCommonCut]) != 0u) { if (sign1 * sign2 < 0) { - fHistMan->FillHistClass(Form("PairsBarrelSEPM_%s_%s_%s", fLegCutNames[icut].Data(), fCommonCutNames[iCommonCut].Data(), fPairCutNames[iPairCut].Data()), VarManager::fgValues); + fHistMan->FillHistClass(Form("PairsBarrelSEPM_%s_%s_%s", fLegCutNames[icut].Data(), fCommonCutNames[iCommonCut].Data(), fPairCutNames[iPairCut].Data()), dqefficiency_helpers::varValues()); } else if (fConfigSameSignHistograms.value) { if (sign1 > 0) { - fHistMan->FillHistClass(Form("PairsBarrelSEPP_%s_%s_%s", fLegCutNames[icut].Data(), fCommonCutNames[iCommonCut].Data(), fPairCutNames[iPairCut].Data()), VarManager::fgValues); + fHistMan->FillHistClass(Form("PairsBarrelSEPP_%s_%s_%s", fLegCutNames[icut].Data(), fCommonCutNames[iCommonCut].Data(), fPairCutNames[iPairCut].Data()), dqefficiency_helpers::varValues()); } else { - fHistMan->FillHistClass(Form("PairsBarrelSEMM_%s_%s_%s", fLegCutNames[icut].Data(), fCommonCutNames[iCommonCut].Data(), fPairCutNames[iPairCut].Data()), VarManager::fgValues); + fHistMan->FillHistClass(Form("PairsBarrelSEMM_%s_%s_%s", fLegCutNames[icut].Data(), fCommonCutNames[iCommonCut].Data(), fPairCutNames[iPairCut].Data()), dqefficiency_helpers::varValues()); } } for (unsigned int isig = 0; isig < fRecMCSignals.size(); isig++) { // loop over MC signals - if (mcDecision & (static_cast(1) << isig)) { + if ((mcDecision & (static_cast(1) << isig)) != 0u) { if (sign1 * sign2 < 0) { - fHistMan->FillHistClass(Form("PairsBarrelSEPM_%s_%s_%s_%s", fLegCutNames[icut].Data(), fCommonCutNames[iCommonCut].Data(), fPairCutNames[iPairCut].Data(), fRecMCSignalNames[isig].Data()), VarManager::fgValues); + fHistMan->FillHistClass(Form("PairsBarrelSEPM_%s_%s_%s_%s", fLegCutNames[icut].Data(), fCommonCutNames[iCommonCut].Data(), fPairCutNames[iPairCut].Data(), fRecMCSignalNames[isig].Data()), dqefficiency_helpers::varValues()); } else if (fConfigSameSignHistograms.value) { if (sign1 > 0) { - fHistMan->FillHistClass(Form("PairsBarrelSEPP_%s_%s_%s_%s", fLegCutNames[icut].Data(), fCommonCutNames[iCommonCut].Data(), fPairCutNames[iPairCut].Data(), fRecMCSignalNames[isig].Data()), VarManager::fgValues); + fHistMan->FillHistClass(Form("PairsBarrelSEPP_%s_%s_%s_%s", fLegCutNames[icut].Data(), fCommonCutNames[iCommonCut].Data(), fPairCutNames[iPairCut].Data(), fRecMCSignalNames[isig].Data()), dqefficiency_helpers::varValues()); } else { - fHistMan->FillHistClass(Form("PairsBarrelSEMM_%s_%s_%s_%s", fLegCutNames[icut].Data(), fCommonCutNames[iCommonCut].Data(), fPairCutNames[iPairCut].Data(), fRecMCSignalNames[isig].Data()), VarManager::fgValues); + fHistMan->FillHistClass(Form("PairsBarrelSEMM_%s_%s_%s_%s", fLegCutNames[icut].Data(), fCommonCutNames[iCommonCut].Data(), fPairCutNames[iPairCut].Data(), fRecMCSignalNames[isig].Data()), dqefficiency_helpers::varValues()); } } } @@ -2180,7 +2165,8 @@ struct AnalysisAsymmetricPairing { } // Function to run same event triplets (e.g. D+->K-pi+pi+) - void runThreeProng(MyEventsVtxCovSelected const& events, PresliceUnsorted& preslice, MyBarrelAssocs const& /*assocs*/, MyBarrelTracksWithCovWithAmbiguities const& tracks, ReducedA3MCEvents const& /*mcEvents*/, ReducedA3MCTracks const& /*mcTracks*/, VarManager::PairCandidateType tripletType) + template + void runThreeProng(TEvents const& events, PresliceUnsorted& preslice, TTrackAssocs const& /*assocs*/, TTracks const& tracks, ReducedA3MCEvents const& /*mcEvents*/, ReducedA3MCTracks const& /*mcTracks*/, VarManager::PairCandidateType tripletType) { for (const auto& event : events) { if (!event.isEventSelected_bit(0)) { @@ -2188,7 +2174,7 @@ struct AnalysisAsymmetricPairing { } // Reset the fValues array VarManager::ResetValues(0, VarManager::kNVars); - VarManager::FillEventAlice3(event, VarManager::fgValues); + VarManager::FillEventAlice3(event, dqefficiency_helpers::varValues()); auto groupedLegAAssocs = legACandidateAssocs.sliceBy(preslice, event.globalIndex()); if (groupedLegAAssocs.size() == 0) { @@ -2206,12 +2192,12 @@ struct AnalysisAsymmetricPairing { // Based on triplet type, make suitable combinations of the partitions if (tripletType == VarManager::kTripleCandidateToPKPi) { for (const auto& [a1, a2, a3] : combinations(soa::CombinationsFullIndexPolicy(groupedLegAAssocs, groupedLegBAssocs, groupedLegCAssocs))) { - readTriplet(a1, a2, a3, tracks, event, tripletType); + readTriplet(a1, a2, a3, tracks, event, tripletType); } } else if (tripletType == VarManager::kTripleCandidateToKPiPi) { for (const auto& a1 : groupedLegAAssocs) { for (const auto& [a2, a3] : combinations(groupedLegBAssocs, groupedLegCAssocs)) { - readTriplet(a1, a2, a3, tracks, event, tripletType); + readTriplet(a1, a2, a3, tracks, event, tripletType); } } } else { @@ -2220,8 +2206,8 @@ struct AnalysisAsymmetricPairing { } // end event loop } - // Helper function to process triplet - void readTriplet(MyBarrelAssocs::iterator const& a1, MyBarrelAssocs::iterator const& a2, MyBarrelAssocs::iterator const& a3, MyBarrelTracksWithCovWithAmbiguities const& /*tracks*/, MyEventsVtxCovSelected::iterator const& /*event*/, VarManager::PairCandidateType tripletType) + template + void readTriplet(TTrackAssoc const& a1, TTrackAssoc const& a2, TTrackAssoc const& a3, TTracks const& /*tracks*/, TEvent const& event, VarManager::PairCandidateType tripletType) { uint32_t mcDecision = 0; @@ -2256,9 +2242,9 @@ struct AnalysisAsymmetricPairing { // Find common track cuts all candidates pass threeTrackCommonFilter |= a1.isBarrelSelected_raw() & a2.isBarrelSelected_raw() & a3.isBarrelSelected_raw() & fCommonTrackCutMask; - auto t1 = a1.template reducedA3track_as(); - auto t2 = a2.template reducedA3track_as(); - auto t3 = a3.template reducedA3track_as(); + auto t1 = a1.template reducedA3track_as(); + auto t2 = a2.template reducedA3track_as(); + auto t3 = a3.template reducedA3track_as(); // Avoid self-pairs if (t1 == t2 || t1 == t3 || t2 == t3) { @@ -2288,64 +2274,64 @@ struct AnalysisAsymmetricPairing { } // run MC matching for this triplet - int isig = 0; + int iSigMc = 0; mcDecision = 0; - for (auto sig = fRecMCSignals.begin(); sig != fRecMCSignals.end(); sig++, isig++) { + for (auto sig = fRecMCSignals.begin(); sig != fRecMCSignals.end(); sig++, iSigMc++) { if (t1.has_reducedA3MCTrack() && t2.has_reducedA3MCTrack() && t3.has_reducedA3MCTrack()) { if ((*sig)->CheckSignal(true, t1.reducedA3MCTrack(), t2.reducedA3MCTrack(), t3.reducedA3MCTrack())) { - mcDecision |= (static_cast(1) << isig); + mcDecision |= (static_cast(1) << iSigMc); } } } // end loop over MC signals - VarManager::FillTriple(t1, t2, t3, VarManager::fgValues, tripletType); - /* TODO: Reimplement when secondary vertexing is available + VarManager::FillTriple(t1, t2, t3, dqefficiency_helpers::varValues(), tripletType); if constexpr (TThreeProngFitter) { - VarManager::FillTripletVertexing(event, t1, t2, t3, tripletType); - }*/ + VarManager::FillTripletVertexingALICE3(event, t1, t2, t3, tripletType); + } // Fill histograms bool isAmbi = false; for (int icut = 0; icut < fNLegCuts; icut++) { isAmbi = (threeTrackFilter & (static_cast(1) << 29)) || (threeTrackFilter & (static_cast(1) << 30)) || (threeTrackFilter & (static_cast(1) << 31)); if (threeTrackFilter & (static_cast(1) << icut)) { - fHistMan->FillHistClass(Form("TripletsBarrelSE_%s", fLegCutNames[icut].Data()), VarManager::fgValues); + fHistMan->FillHistClass(Form("TripletsBarrelSE_%s", fLegCutNames[icut].Data()), dqefficiency_helpers::varValues()); for (unsigned int isig = 0; isig < fRecMCSignals.size(); isig++) { // loop over MC signals if (mcDecision & (static_cast(1) << isig)) { - fHistMan->FillHistClass(Form("TripletsBarrelSE_%s_%s", fLegCutNames[icut].Data(), fRecMCSignalNames[isig].Data()), VarManager::fgValues); // matched signal + fHistMan->FillHistClass(Form("TripletsBarrelSE_%s_%s", fLegCutNames[icut].Data(), fRecMCSignalNames[isig].Data()), dqefficiency_helpers::varValues()); // matched signal } } // end loop (MC signals) if (fConfigQA && isAmbi) { - fHistMan->FillHistClass(Form("TripletsBarrelSE_ambiguous_%s", fLegCutNames[icut].Data()), VarManager::fgValues); + fHistMan->FillHistClass(Form("TripletsBarrelSE_ambiguous_%s", fLegCutNames[icut].Data()), dqefficiency_helpers::varValues()); } for (int iCommonCut = 0; iCommonCut < fNCommonTrackCuts; iCommonCut++) { if (threeTrackCommonFilter & fCommonTrackCutFilterMasks[iCommonCut]) { - fHistMan->FillHistClass(Form("TripletsBarrelSE_%s_%s", fLegCutNames[icut].Data(), fCommonCutNames[iCommonCut].Data()), VarManager::fgValues); + fHistMan->FillHistClass(Form("TripletsBarrelSE_%s_%s", fLegCutNames[icut].Data(), fCommonCutNames[iCommonCut].Data()), dqefficiency_helpers::varValues()); for (unsigned int isig = 0; isig < fRecMCSignals.size(); isig++) { // loop over MC signals if (mcDecision & (static_cast(1) << isig)) { - fHistMan->FillHistClass(Form("TripletsBarrelSE_%s_%s_%s", fLegCutNames[icut].Data(), fCommonCutNames[iCommonCut].Data(), fRecMCSignalNames[isig].Data()), VarManager::fgValues); // matched signal + fHistMan->FillHistClass(Form("TripletsBarrelSE_%s_%s_%s", fLegCutNames[icut].Data(), fCommonCutNames[iCommonCut].Data(), fRecMCSignalNames[isig].Data()), dqefficiency_helpers::varValues()); // matched signal } } // end loop (MC signals) } } // end loop (common cuts) int iPairCut = 0; for (auto cut = fPairCuts.begin(); cut != fPairCuts.end(); cut++, iPairCut++) { - if (!((*cut)->IsSelected(VarManager::fgValues))) // apply pair cuts + if (!((*cut)->IsSelected(dqefficiency_helpers::varValues()))) { // apply pair cuts continue; + } // Histograms with pair cuts - fHistMan->FillHistClass(Form("TripletsBarrelSE_%s_%s", fLegCutNames[icut].Data(), fPairCutNames[iPairCut].Data()), VarManager::fgValues); + fHistMan->FillHistClass(Form("TripletsBarrelSE_%s_%s", fLegCutNames[icut].Data(), fPairCutNames[iPairCut].Data()), dqefficiency_helpers::varValues()); for (unsigned int isig = 0; isig < fRecMCSignals.size(); isig++) { // loop over MC signals if (mcDecision & (static_cast(1) << isig)) { - fHistMan->FillHistClass(Form("TripletsBarrelSE_%s_%s_%s", fLegCutNames[icut].Data(), fPairCutNames[iPairCut].Data(), fRecMCSignalNames[isig].Data()), VarManager::fgValues); // matched signal + fHistMan->FillHistClass(Form("TripletsBarrelSE_%s_%s_%s", fLegCutNames[icut].Data(), fPairCutNames[iPairCut].Data(), fRecMCSignalNames[isig].Data()), dqefficiency_helpers::varValues()); // matched signal } } // end loop (MC signals) // Histograms with pair cuts and common track cuts for (int iCommonCut = 0; iCommonCut < fNCommonTrackCuts; ++iCommonCut) { if (threeTrackCommonFilter & fCommonTrackCutFilterMasks[iCommonCut]) { - fHistMan->FillHistClass(Form("TripletsBarrelSE_%s_%s_%s", fLegCutNames[icut].Data(), fCommonCutNames[iCommonCut].Data(), fPairCutNames[iPairCut].Data()), VarManager::fgValues); + fHistMan->FillHistClass(Form("TripletsBarrelSE_%s_%s_%s", fLegCutNames[icut].Data(), fCommonCutNames[iCommonCut].Data(), fPairCutNames[iPairCut].Data()), dqefficiency_helpers::varValues()); for (unsigned int isig = 0; isig < fRecMCSignals.size(); isig++) { // loop over MC signals if (mcDecision & (static_cast(1) << isig)) { - fHistMan->FillHistClass(Form("TripletsBarrelSE_%s_%s_%s_%s", fLegCutNames[icut].Data(), fCommonCutNames[iCommonCut].Data(), fPairCutNames[iPairCut].Data(), fRecMCSignalNames[isig].Data()), VarManager::fgValues); // matched signal + fHistMan->FillHistClass(Form("TripletsBarrelSE_%s_%s_%s_%s", fLegCutNames[icut].Data(), fCommonCutNames[iCommonCut].Data(), fPairCutNames[iPairCut].Data(), fRecMCSignalNames[isig].Data()), dqefficiency_helpers::varValues()); // matched signal } } // end loop (MC signals) } @@ -2368,7 +2354,15 @@ struct AnalysisAsymmetricPairing { MyBarrelTracksWithCovWithAmbiguities const& barrelTracks, ReducedA3MCEvents const& mcEvents, ReducedA3MCTracks const& mcTracks) { - runThreeProng(events, trackAssocsPerCollision, barrelAssocs, barrelTracks, mcEvents, mcTracks, VarManager::kTripleCandidateToKPiPi); + runThreeProng(events, trackAssocsPerCollision, barrelAssocs, barrelTracks, mcEvents, mcTracks, VarManager::kTripleCandidateToKPiPi); + } + + void processProtonKaonPionSkimmed(MyEventsVtxCovSelected const& events, + MyBarrelAssocs const& barrelAssocs, + MyBarrelTracksWithCovWithAmbiguities const& barrelTracks, + ReducedA3MCEvents const& mcEvents, ReducedA3MCTracks const& mcTracks) + { + runThreeProng(events, trackAssocsPerCollision, barrelAssocs, barrelTracks, mcEvents, mcTracks, VarManager::kTripleCandidateToPKPi); } void processMCGen(ReducedA3MCTracks const& mcTracks) @@ -2384,7 +2378,7 @@ struct AnalysisAsymmetricPairing { // TODO: Use the mcReducedFlags to select signals for (const auto& sig : fGenMCSignals) { if (sig->CheckSignal(true, mctrack)) { - fHistMan->FillHistClass(Form("MCTruthGen_%s", sig->GetName()), VarManager::fgValues); + fHistMan->FillHistClass(Form("MCTruthGen_%s", sig->GetName()), dqefficiency_helpers::varValues()); } } } @@ -2412,7 +2406,7 @@ struct AnalysisAsymmetricPairing { auto track_raw = groupedMCTracks.rawIteratorAt(track.globalIndex()); for (const auto& sig : fGenMCSignals) { if (sig->CheckSignal(true, track_raw)) { - fHistMan->FillHistClass(Form("MCTruthGenSel_%s", sig->GetName()), VarManager::fgValues); + fHistMan->FillHistClass(Form("MCTruthGenSel_%s", sig->GetName()), dqefficiency_helpers::varValues()); } } } @@ -2426,6 +2420,7 @@ struct AnalysisAsymmetricPairing { PROCESS_SWITCH(AnalysisAsymmetricPairing, processKaonPionSkimmed, "Run kaon pion pairing, with skimmed tracks", false); PROCESS_SWITCH(AnalysisAsymmetricPairing, processKaonPionPionSkimmed, "Run kaon pion pion triplets, with skimmed tracks", false); + PROCESS_SWITCH(AnalysisAsymmetricPairing, processProtonKaonPionSkimmed, "Run proton kaon pion triplets, with skimmed tracks", false); PROCESS_SWITCH(AnalysisAsymmetricPairing, processMCGen, "Loop over MC particle stack and fill generator level histograms", false); PROCESS_SWITCH(AnalysisAsymmetricPairing, processMCGenWithEventSelection, "Loop over MC particle stack and fill generator level histograms", false); PROCESS_SWITCH(AnalysisAsymmetricPairing, processDummy, "Dummy function, enabled only if none of the others are enabled", true); @@ -2441,7 +2436,7 @@ WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) adaptAnalysisTask(cfgc)}; } -void DefineHistograms(HistogramManager* histMan, TString histClasses, const char* histGroups) +void DefineHistograms(HistogramManager* histMan, const TString& histClasses, const char* histGroups) { // // Define here the histograms for all the classes required in analysis. diff --git a/PWGDQ/Core/CutsLibrary.cxx b/PWGDQ/Core/CutsLibrary.cxx index 8c069d64920..fae1754b529 100644 --- a/PWGDQ/Core/CutsLibrary.cxx +++ b/PWGDQ/Core/CutsLibrary.cxx @@ -4012,9 +4012,7 @@ AnalysisCompositeCut* o2::aod::dqcuts::GetCompositeCut(const char* cutName) if (nameStr == "alice3DielectronPID") { cut->AddCut(GetAnalysisCut("alice3JpsiKine")); cut->AddCut(GetAnalysisCut("alice3TrackQuality")); - cut->AddCut(GetAnalysisCut("alice3iTOFPIDEl")); - cut->AddCut(GetAnalysisCut("alice3oTOFPIDEl")); - cut->AddCut(GetAnalysisCut("alice3RICHPIDEl")); + cut->AddCut(GetAnalysisCut("alice3CharmoniumPID")); return cut; } @@ -4060,6 +4058,12 @@ AnalysisCompositeCut* o2::aod::dqcuts::GetCompositeCut(const char* cutName) return cut; } + if (nameStr == "alice3LambdaCQualityCuts") { + cut->AddCut(GetAnalysisCut("alice3LambdaCKine")); + cut->AddCut(GetAnalysisCut("alice3TrackQuality")); + return cut; + } + delete cut; LOGF(fatal, Form("Did not find cut %s. Returning nullptr", cutName)); return nullptr; @@ -7264,29 +7268,33 @@ AnalysisCut* o2::aod::dqcuts::GetAnalysisCut(const char* cutName) return cut; } + if (nameStr == "alice3LambdaCKine") { + cut->AddCut(VarManager::kPt, 0.2, 1000.0); + cut->AddCut(VarManager::kEta, -2.0, 2.0); + return cut; + } + if (nameStr == "alice3JpsiKine") { cut->AddCut(VarManager::kPt, 1.0, 1000.0); - cut->AddCut(VarManager::kEta, -2.5, 2.5); // Total tracker acceptance in v3b geomety + cut->AddCut(VarManager::kEta, -2.5, 2.5); return cut; } if (nameStr == "alice3JpsiKineTOFAcceptance") { cut->AddCut(VarManager::kPt, 1.0, 1000.0); - cut->AddCut(VarManager::kEta, -2.0, 2.0); // TOF acceptance in v3b geomety + cut->AddCut(VarManager::kEta, -2.0, 2.0); return cut; } if (nameStr == "alice3JpsiKineRICHAcceptance") { cut->AddCut(VarManager::kPt, 1.0, 1000.0); - cut->AddCut(VarManager::kEta, -0.8, 0.8); // RICH acceptance in v3b geomety + cut->AddCut(VarManager::kEta, -0.8, 0.8); return cut; } if (nameStr == "alice3TrackQuality") { cut->AddCut(VarManager::kIsReconstructed, 0.5, 1.5); - cut->AddCut(VarManager::kNSiliconHits, 6.0, 12.0); - cut->AddCut(VarManager::kTrackDCAxy, -3.0, 3.0); - cut->AddCut(VarManager::kTrackDCAz, -3.0, 3.0); + cut->AddCut(VarManager::kNSiliconHits, 5.0, 12.0); return cut; } @@ -7363,27 +7371,29 @@ AnalysisCut* o2::aod::dqcuts::GetAnalysisCut(const char* cutName) return cut; } + if (nameStr == "alice3CharmoniumPID") { + cut->AddCut(VarManager::kOuterTOFnSigmaEl, -2.0, 3.0, false, VarManager::kP, 0.0, 1.2); + cut->AddCut(VarManager::kRICHnSigmaEl, -2.0, 3.0, false, VarManager::kHasRICHSigEl, 0.5, 1.5, false, VarManager::kP, 0.7, 1000.0); + return cut; + } + if (nameStr == "alice3RICHPIDEl") { - cut->AddCut(VarManager::kRICHnSigmaEl, -3.0, 3.0); - cut->AddCut(VarManager::kHasRICHSigEl, 0.5, 1.5); + cut->AddCut(VarManager::kRICHnSigmaEl, -3.0, 3.0, false, VarManager::kP, 0.7, 1000.0, false, VarManager::kHasRICHSigEl, 0.5, 1.5); return cut; } if (nameStr == "alice3RICHPIDPi") { - cut->AddCut(VarManager::kRICHnSigmaPi, -3.0, 3.0); - cut->AddCut(VarManager::kHasRICHSigPi, 0.5, 1.5); + cut->AddCut(VarManager::kRICHnSigmaPi, -3.0, 3.0, false, VarManager::kP, 0.57, 1000.0, false, VarManager::kHasRICHSigPi, 0.5, 1.5); return cut; } if (nameStr == "alice3RICHPIDKa") { - cut->AddCut(VarManager::kRICHnSigmaKa, -3.0, 3.0); - cut->AddCut(VarManager::kHasRICHSigKa, 0.5, 1.5); + cut->AddCut(VarManager::kRICHnSigmaKa, -3.0, 3.0, false, VarManager::kP, 2.0, 1000.0, false, VarManager::kHasRICHSigKa, 0.5, 1.5); return cut; } if (nameStr == "alice3RICHPIDPr") { - cut->AddCut(VarManager::kRICHnSigmaPr, -3.0, 3.0); - cut->AddCut(VarManager::kHasRICHSigPr, 0.5, 1.5); + cut->AddCut(VarManager::kRICHnSigmaPr, -3.0, 3.0, false, VarManager::kP, 3.8, 1000.0, false, VarManager::kHasRICHSigPr, 0.5, 1.5); return cut; } diff --git a/PWGDQ/Core/MCSignalLibrary.cxx b/PWGDQ/Core/MCSignalLibrary.cxx index 19dcbbc8aa9..bdd86156c7d 100644 --- a/PWGDQ/Core/MCSignalLibrary.cxx +++ b/PWGDQ/Core/MCSignalLibrary.cxx @@ -1717,6 +1717,13 @@ MCSignal* o2::aod::dqmcsignals::GetMCSignal(const char* name) signal = new MCSignal(name, "Lambda_c", {prong}, {-1}); return signal; } + if (nameStr == "PrKPiFromLambdaC") { + MCProng prongProton(2, {2212, Pdg::kLambdaCPlus}, {true, true}, {false, false}, {0, 0}, {0, 0}, {false, false}); + MCProng prongKaon(2, {321, Pdg::kLambdaCPlus}, {true, true}, {false, false}, {0, 0}, {0, 0}, {false, false}); + MCProng prongPion(2, {211, Pdg::kLambdaCPlus}, {true, true}, {false, false}, {0, 0}, {0, 0}, {false, false}); + signal = new MCSignal(name, "Proton kaon pion triplet from Lambda_c+/-", {prongProton, prongKaon, prongPion}, {1, 1, 1}); + return signal; + } //-------------------------------------------------------------------------------- diff --git a/PWGDQ/Core/VarManager.h b/PWGDQ/Core/VarManager.h index b9ac54a5616..45006f3d030 100644 --- a/PWGDQ/Core/VarManager.h +++ b/PWGDQ/Core/VarManager.h @@ -1495,6 +1495,10 @@ class VarManager : public TObject static void FillTrackAlice3(T const& track, float* values = nullptr); template static void FillResolutions(M const& mcTrack, T const& track, float* values = nullptr); + template + static void FillPairVertexingAlice3(C const& collision, T const& t1, T const& t2, bool propToSV = false, float* values = nullptr); + template + static void FillTripletVertexingALICE3(C const& collision, T const& t1, T const& t2, T const& t3, VarManager::PairCandidateType tripletType, float* values = nullptr); static void SetCalibrationObject(CalibObjects calib, TObject* obj) { @@ -1606,10 +1610,10 @@ class VarManager : public TObject static o2::vertexing::FwdDCAFitterN<3> fgFitterThreeProngFwd; static o2::globaltracking::MatchGlobalFwd mMatching; - static std::map fgCalibs; // map of calibration histograms + static std::map fgCalibs; // map of calibration histograms static std::array fgRunTPCPostCalibration; // 0-electron, 1-pion, 2-kaon, 3-proton - static int fgCalibrationType; // 0 - no calibration, 1 - calibration vs (TPCncls,pIN,eta) typically for pp, 2 - calibration vs (eta,nPV,nLong,tLong) typically for PbPb - static bool fgUseInterpolatedCalibration; // use interpolated calibration histograms (default: true) + static int fgCalibrationType; // 0 - no calibration, 1 - calibration vs (TPCncls,pIN,eta) typically for pp, 2 - calibration vs (eta,nPV,nLong,tLong) typically for PbPb + static bool fgUseInterpolatedCalibration; // use interpolated calibration histograms (default: true) static int fgEfficiencyType; // type of efficiency correction to apply static TObject* fgEfficiencyHist; // histogram for efficiency correction @@ -7542,4 +7546,308 @@ void VarManager::FillResolutions(M const& mcTrack, T const& track, float* values values[kEtaResolution] = track.eta() - mcTrack.eta(); } +template +void VarManager::FillPairVertexingAlice3(C const& collision, T const& t1, T const& t2, bool propToSV, float* values) +{ + // check at compile time that the event and cov matrix have the cov matrix + constexpr bool eventHasVtxCov = ((collFillMap & Collision) > 0 || (collFillMap & ReducedEventVtxCov) > 0); + constexpr bool trackHasCov = ((fillMap & TrackCov) > 0 || (fillMap & ReducedTrackBarrelCov) > 0); + constexpr bool muonHasCov = ((fillMap & MuonCov) > 0 || (fillMap & ReducedMuonCov) > 0); + + if (!values) { + values = fgValues; + } + float m1 = o2::constants::physics::MassElectron; + float m2 = o2::constants::physics::MassElectron; + if constexpr (pairType == kDecayToKPi) { + m1 = o2::constants::physics::MassKaonCharged; + m2 = o2::constants::physics::MassPionCharged; + } + if constexpr (pairType == kDecayToMuMu && muonHasCov) { + m1 = o2::constants::physics::MassMuon; + m2 = o2::constants::physics::MassMuon; + } + ROOT::Math::PtEtaPhiMVector v1(t1.pt(), t1.eta(), t1.phi(), m1); + ROOT::Math::PtEtaPhiMVector v2(t2.pt(), t2.eta(), t2.phi(), m2); + ROOT::Math::PtEtaPhiMVector v12 = v1 + v2; + + values[kUsedKF] = static_cast(fgUsedKF); + if (!fgUsedKF) { + int procCode = 0; + + // TODO: use trackUtilities functions to initialize the various matrices to avoid code duplication + // auto pars1 = getTrackParCov(t1); + // auto pars2 = getTrackParCov(t2); + // We need to hide the cov data members from the cases when no cov table is provided + if constexpr ((pairType == kDecayToEE || pairType == kDecayToKPi) && trackHasCov) { + std::array t1pars = {t1.y(), t1.z(), t1.snp(), t1.tgl(), t1.signed1Pt()}; + std::array t1covs = {t1.cYY(), t1.cZY(), t1.cZZ(), t1.cSnpY(), t1.cSnpZ(), + t1.cSnpSnp(), t1.cTglY(), t1.cTglZ(), t1.cTglSnp(), t1.cTglTgl(), + t1.c1PtY(), t1.c1PtZ(), t1.c1PtSnp(), t1.c1PtTgl(), t1.c1Pt21Pt2()}; + o2::track::TrackParCov pars1{t1.x(), t1.alpha(), t1pars, t1covs}; + std::array t2pars = {t2.y(), t2.z(), t2.snp(), t2.tgl(), t2.signed1Pt()}; + std::array t2covs = {t2.cYY(), t2.cZY(), t2.cZZ(), t2.cSnpY(), t2.cSnpZ(), + t2.cSnpSnp(), t2.cTglY(), t2.cTglZ(), t2.cTglSnp(), t2.cTglTgl(), + t2.c1PtY(), t2.c1PtZ(), t2.c1PtSnp(), t2.c1PtTgl(), t2.c1Pt21Pt2()}; + o2::track::TrackParCov pars2{t2.x(), t2.alpha(), t2pars, t2covs}; + procCode = fgFitterTwoProngBarrel.process(pars1, pars2); + } else if constexpr ((pairType == kDecayToMuMu) && muonHasCov) { + // Initialize track parameters for forward + o2::track::TrackParCovFwd pars1 = FwdToTrackPar(t1, t1); + o2::track::TrackParCovFwd pars2 = FwdToTrackPar(t2, t2); + procCode = fgFitterTwoProngFwd.process(pars1, pars2); + } else { + return; + } + + values[kVertexingProcCode] = procCode; + if (procCode == 0) { + // TODO: set the other variables to appropriate values and return + values[kVertexingChi2PCA] = -999.; + values[kVertexingLxy] = -999.; + values[kVertexingLxyz] = -999.; + values[kVertexingLz] = -999.; + values[kVertexingLxyErr] = -999.; + values[kVertexingLxyzErr] = -999.; + values[kVertexingLzErr] = -999.; + + values[kVertexingTauxy] = -999.; + values[kVertexingTauz] = -999.; + values[kVertexingTauxyErr] = -999.; + values[kVertexingTauzErr] = -999.; + values[kVertexingPz] = -999.; + values[kVertexingSV] = -999.; + return; + } + + Vec3D secondaryVertex; + o2::dataformats::VertexBase primaryVertexNew; + + if constexpr (eventHasVtxCov) { + + std::array covMatrixPCA{}; + // get track impact parameters + // This modifies track momenta! + o2::math_utils::Point3D vtxXYZ(collision.posX(), collision.posY(), collision.posZ()); + std::array vtxCov{collision.covXX(), collision.covXY(), collision.covYY(), collision.covXZ(), collision.covYZ(), collision.covZZ()}; + o2::dataformats::VertexBase primaryVertex = {vtxXYZ, vtxCov}; + // auto primaryVertex = getPrimaryVertex(collision); + auto covMatrixPV = primaryVertex.getCov(); + + if constexpr ((pairType == kDecayToEE || pairType == kDecayToKPi) && trackHasCov) { + secondaryVertex = fgFitterTwoProngBarrel.getPCACandidate(); + // printf("secVtx (first) %f %f %f \n",secondaryVertex[0],secondaryVertex[1],secondaryVertex[2]); + covMatrixPCA = fgFitterTwoProngBarrel.calcPCACovMatrixFlat(); + auto chi2PCA = fgFitterTwoProngBarrel.getChi2AtPCACandidate(); + auto trackParVar0 = fgFitterTwoProngBarrel.getTrack(0); + auto trackParVar1 = fgFitterTwoProngBarrel.getTrack(1); + values[kVertexingChi2PCA] = chi2PCA; + v1 = {trackParVar0.getPt(), trackParVar0.getEta(), trackParVar0.getPhi(), m1}; + v2 = {trackParVar1.getPt(), trackParVar1.getEta(), trackParVar1.getPhi(), m2}; + v12 = v1 + v2; + } + double phi = std::atan2(secondaryVertex[1] - collision.posY(), secondaryVertex[0] - collision.posX()); + double theta = std::atan2(secondaryVertex[2] - collision.posZ(), + std::sqrt((secondaryVertex[0] - collision.posX()) * (secondaryVertex[0] - collision.posX()) + + (secondaryVertex[1] - collision.posY()) * (secondaryVertex[1] - collision.posY()))); + + values[kVertexingLxyzErr] = std::sqrt(getRotatedCovMatrixXX(covMatrixPV, phi, theta) + getRotatedCovMatrixXX(covMatrixPCA, phi, theta)); + values[kVertexingLxyErr] = std::sqrt(getRotatedCovMatrixXX(covMatrixPV, phi, 0.) + getRotatedCovMatrixXX(covMatrixPCA, phi, 0.)); + values[kVertexingLzErr] = std::sqrt(getRotatedCovMatrixXX(covMatrixPV, 0, theta) + getRotatedCovMatrixXX(covMatrixPCA, 0, theta)); + + values[kVertexingLxy] = (collision.posX() - secondaryVertex[0]) * (collision.posX() - secondaryVertex[0]) + + (collision.posY() - secondaryVertex[1]) * (collision.posY() - secondaryVertex[1]); + values[kVertexingLz] = (collision.posZ() - secondaryVertex[2]) * (collision.posZ() - secondaryVertex[2]); + values[kVertexingLxyz] = values[kVertexingLxy] + values[kVertexingLz]; + values[kVertexingLxy] = std::sqrt(values[kVertexingLxy]); + values[kVertexingLz] = std::sqrt(values[kVertexingLz]); + values[kVertexingLxyz] = std::sqrt(values[kVertexingLxyz]); + + values[kVertexingTauz] = (collision.posZ() - secondaryVertex[2]) * v12.M() / (TMath::Abs(v12.Pz()) * o2::constants::physics::LightSpeedCm2NS); + values[kVertexingTauxy] = values[kVertexingLxy] * v12.M() / (v12.Pt() * o2::constants::physics::LightSpeedCm2NS); + + values[kVertexingPz] = TMath::Abs(v12.Pz()); + values[kVertexingSV] = secondaryVertex[2]; + + values[kVertexingTauzErr] = values[kVertexingLzErr] * v12.M() / (TMath::Abs(v12.Pz()) * o2::constants::physics::LightSpeedCm2NS); + values[kVertexingTauxyErr] = values[kVertexingLxyErr] * v12.M() / (v12.Pt() * o2::constants::physics::LightSpeedCm2NS); + + values[kCosPointingAngle] = ((secondaryVertex[0] - collision.posX()) * v12.Px() + + (secondaryVertex[1] - collision.posY()) * v12.Py() + + (secondaryVertex[2] - collision.posZ()) * v12.Pz()) / + (v12.P() * values[VarManager::kVertexingLxyz]); + // Decay length defined as in Run 2 + values[kVertexingLzProjected] = ((secondaryVertex[2] - collision.posZ()) * v12.Pz()) / TMath::Sqrt(v12.Pz() * v12.Pz()); + values[kVertexingLxyProjected] = ((secondaryVertex[0] - collision.posX()) * v12.Px()) + ((secondaryVertex[1] - collision.posY()) * v12.Py()); + values[kVertexingLxyProjected] = values[kVertexingLxyProjected] / TMath::Sqrt((v12.Px() * v12.Px()) + (v12.Py() * v12.Py())); + values[kVertexingLxyzProjected] = ((secondaryVertex[0] - collision.posX()) * v12.Px()) + ((secondaryVertex[1] - collision.posY()) * v12.Py()) + ((secondaryVertex[2] - collision.posZ()) * v12.Pz()); + values[kVertexingLxyzProjected] = values[kVertexingLxyzProjected] / TMath::Sqrt((v12.Px() * v12.Px()) + (v12.Py() * v12.Py()) + (v12.Pz() * v12.Pz())); + if (fgPVrecalKF) { + values[kVertexingLxyProjectedRecalculatePV] = (secondaryVertex[0] - primaryVertexNew.getX()) * v12.Px() + (secondaryVertex[1] - primaryVertexNew.getY()) * v12.Py(); + values[kVertexingLxyProjectedRecalculatePV] = values[kVertexingLxyProjectedRecalculatePV] / v12.Pt(); + } + values[kVertexingTauxyProjected] = values[kVertexingLxyProjected] * v12.M() / (v12.Pt()); + values[kVertexingTauxyProjectedPoleJPsiMass] = values[kVertexingLxyProjected] * o2::constants::physics::MassJPsi / (v12.Pt()); + values[kVertexingTauxyProjectedNs] = values[kVertexingTauxyProjected] / o2::constants::physics::LightSpeedCm2NS; + if (fgPVrecalKF) { + values[kVertexingTauxyProjectedPoleJPsiMassRecalculatePV] = values[kVertexingLxyProjectedRecalculatePV] * o2::constants::physics::MassJPsi / (v12.Pt()); + } + values[kVertexingTauzProjected] = values[kVertexingLzProjected] * v12.M() / TMath::Abs(v12.Pz()); + values[kVertexingTauxyzProjected] = values[kVertexingLxyzProjected] * v12.M() / (v12.P()); + } + } + if (propToSV) { + values[kMass] = v12.M(); + values[kPt] = v12.Pt(); + values[kEta] = v12.Eta(); + // values[kPhi] = v12.Phi(); + values[kPhi] = RecoDecay::constrainAngle(v12.Phi()); + } else { + values[kPt1] = t1.pt(); + values[kEta1] = t1.eta(); + values[kPhi1] = t1.phi(); + + values[kPt2] = t2.pt(); + values[kEta2] = t2.eta(); + values[kPhi2] = t2.phi(); + } +} + +template +void VarManager::FillTripletVertexingALICE3(C const& collision, T const& t1, T const& t2, T const& t3, VarManager::PairCandidateType tripletType, float* values) +{ + // TODO: Vertexing error variables + constexpr bool eventHasVtxCov = ((collFillMap & Collision) > 0 || (collFillMap & ReducedEventVtxCov) > 0); + bool trackHasCov = ((fillMap & ReducedTrackBarrelCov) > 0); + + if (!values) { + values = fgValues; + } + + float m1 = o2::constants::physics::MassKaonCharged; + float m2 = o2::constants::physics::MassPionCharged; + float m3 = o2::constants::physics::MassPionCharged; + + if (tripletType == kTripleCandidateToKPiPi) { + m1 = o2::constants::physics::MassKaonCharged; + m2 = o2::constants::physics::MassPionCharged; + m3 = o2::constants::physics::MassPionCharged; + } + if (tripletType == kTripleCandidateToPKPi) { + m1 = o2::constants::physics::MassProton; + m2 = o2::constants::physics::MassKaonCharged; + m3 = o2::constants::physics::MassPionCharged; + } + ROOT::Math::PtEtaPhiMVector v1(t1.pt(), t1.eta(), t1.phi(), m1); + ROOT::Math::PtEtaPhiMVector v2(t2.pt(), t2.eta(), t2.phi(), m2); + ROOT::Math::PtEtaPhiMVector v3(t3.pt(), t3.eta(), t3.phi(), m3); + ROOT::Math::PtEtaPhiMVector v123 = v1 + v2 + v3; + + int procCode = 0; + + if (trackHasCov) { + std::array t1pars = {t1.y(), t1.z(), t1.snp(), t1.tgl(), t1.signed1Pt()}; + std::array t1covs = {t1.cYY(), t1.cZY(), t1.cZZ(), t1.cSnpY(), t1.cSnpZ(), + t1.cSnpSnp(), t1.cTglY(), t1.cTglZ(), t1.cTglSnp(), t1.cTglTgl(), + t1.c1PtY(), t1.c1PtZ(), t1.c1PtSnp(), t1.c1PtTgl(), t1.c1Pt21Pt2()}; + o2::track::TrackParCov pars1{t1.x(), t1.alpha(), t1pars, t1covs}; + std::array t2pars = {t2.y(), t2.z(), t2.snp(), t2.tgl(), t2.signed1Pt()}; + std::array t2covs = {t2.cYY(), t2.cZY(), t2.cZZ(), t2.cSnpY(), t2.cSnpZ(), + t2.cSnpSnp(), t2.cTglY(), t2.cTglZ(), t2.cTglSnp(), t2.cTglTgl(), + t2.c1PtY(), t2.c1PtZ(), t2.c1PtSnp(), t2.c1PtTgl(), t2.c1Pt21Pt2()}; + o2::track::TrackParCov pars2{t2.x(), t2.alpha(), t2pars, t2covs}; + std::array t3pars = {t3.y(), t3.z(), t3.snp(), t3.tgl(), t3.signed1Pt()}; + std::array t3covs = {t3.cYY(), t3.cZY(), t3.cZZ(), t3.cSnpY(), t3.cSnpZ(), + t3.cSnpSnp(), t3.cTglY(), t3.cTglZ(), t3.cTglSnp(), t3.cTglTgl(), + t3.c1PtY(), t3.c1PtZ(), t3.c1PtSnp(), t3.c1PtTgl(), t3.c1Pt21Pt2()}; + o2::track::TrackParCov pars3{t3.x(), t3.alpha(), t3pars, t3covs}; + procCode = VarManager::fgFitterThreeProngBarrel.process(pars1, pars2, pars3); + } else { + return; + } + + values[VarManager::kVertexingProcCode] = procCode; + if (procCode == 0) { + // TODO: set the other variables to appropriate values and return + values[kVertexingChi2PCA] = -999.; + values[kVertexingLxy] = -999.; + values[kVertexingLxyz] = -999.; + values[kVertexingLz] = -999.; + values[kVertexingLxyErr] = -999.; + values[kVertexingLxyzErr] = -999.; + values[kVertexingLzErr] = -999.; + + values[kVertexingTauxy] = -999.; + values[kVertexingTauz] = -999.; + values[kVertexingTauxyErr] = -999.; + values[kVertexingTauzErr] = -999.; + + values[kVertexingLzProjected] = -999.; + values[kVertexingLxyProjected] = -999.; + values[kVertexingLxyzProjected] = -999.; + values[kVertexingTauzProjected] = -999.; + values[kVertexingTauxyProjected] = -999.; + values[kVertexingTauxyzProjected] = -999.; + + return; + } + + Vec3D secondaryVertex; + + if constexpr (eventHasVtxCov) { + secondaryVertex = fgFitterThreeProngBarrel.getPCACandidate(); + + std::array covMatrixPCA = fgFitterThreeProngBarrel.calcPCACovMatrixFlat(); + + o2::math_utils::Point3D vtxXYZ(collision.posX(), collision.posY(), collision.posZ()); + std::array vtxCov{collision.covXX(), collision.covXY(), collision.covYY(), collision.covXZ(), collision.covYZ(), collision.covZZ()}; + o2::dataformats::VertexBase primaryVertex = {vtxXYZ, vtxCov}; + auto covMatrixPV = primaryVertex.getCov(); + + if (fgUsedVars[kVertexingChi2PCA]) { + auto chi2PCA = fgFitterThreeProngBarrel.getChi2AtPCACandidate(); + values[VarManager::kVertexingChi2PCA] = chi2PCA; + } + + double phi = std::atan2(secondaryVertex[1] - collision.posY(), secondaryVertex[0] - collision.posX()); + double theta = std::atan2(secondaryVertex[2] - collision.posZ(), + std::sqrt((secondaryVertex[0] - collision.posX()) * (secondaryVertex[0] - collision.posX()) + + (secondaryVertex[1] - collision.posY()) * (secondaryVertex[1] - collision.posY()))); + + values[kVertexingLxy] = (collision.posX() - secondaryVertex[0]) * (collision.posX() - secondaryVertex[0]) + + (collision.posY() - secondaryVertex[1]) * (collision.posY() - secondaryVertex[1]); + values[kVertexingLz] = (collision.posZ() - secondaryVertex[2]) * (collision.posZ() - secondaryVertex[2]); + values[kVertexingLxyz] = values[kVertexingLxy] + values[kVertexingLz]; + values[kVertexingLxy] = std::sqrt(values[kVertexingLxy]); + values[kVertexingLz] = std::sqrt(values[kVertexingLz]); + values[kVertexingLxyz] = std::sqrt(values[kVertexingLxyz]); + + values[kVertexingLxyzErr] = std::sqrt(getRotatedCovMatrixXX(covMatrixPV, phi, theta) + getRotatedCovMatrixXX(covMatrixPCA, phi, theta)); + values[kVertexingLxyErr] = std::sqrt(getRotatedCovMatrixXX(covMatrixPV, phi, 0.) + getRotatedCovMatrixXX(covMatrixPCA, phi, 0.)); + values[kVertexingLzErr] = std::sqrt(getRotatedCovMatrixXX(covMatrixPV, 0, theta) + getRotatedCovMatrixXX(covMatrixPCA, 0, theta)); + + values[kVertexingTauz] = (collision.posZ() - secondaryVertex[2]) * v123.M() / (TMath::Abs(v123.Pz()) * o2::constants::physics::LightSpeedCm2NS); + values[kVertexingTauxy] = values[kVertexingLxy] * v123.M() / (v123.Pt() * o2::constants::physics::LightSpeedCm2NS); + + values[kVertexingTauzErr] = values[kVertexingLzErr] * v123.M() / (TMath::Abs(v123.Pz()) * o2::constants::physics::LightSpeedCm2NS); + values[kVertexingTauxyErr] = values[kVertexingLxyErr] * v123.M() / (v123.Pt() * o2::constants::physics::LightSpeedCm2NS); + + values[kCosPointingAngle] = ((secondaryVertex[0] - collision.posX()) * v123.Px() + + (secondaryVertex[1] - collision.posY()) * v123.Py() + + (secondaryVertex[2] - collision.posZ()) * v123.Pz()) / + (v123.P() * values[VarManager::kVertexingLxyz]); + // run 2 definitions: Decay length projected onto the momentum vector of the candidate + values[kVertexingLzProjected] = (secondaryVertex[2] - collision.posZ()) * v123.Pz(); + values[kVertexingLzProjected] = values[kVertexingLzProjected] / TMath::Sqrt(v123.Pz() * v123.Pz()); + values[kVertexingLxyProjected] = ((secondaryVertex[0] - collision.posX()) * v123.Px()) + ((secondaryVertex[1] - collision.posY()) * v123.Py()); + values[kVertexingLxyProjected] = values[kVertexingLxyProjected] / TMath::Sqrt((v123.Px() * v123.Px()) + (v123.Py() * v123.Py())); + values[kVertexingLxyzProjected] = ((secondaryVertex[0] - collision.posX()) * v123.Px()) + ((secondaryVertex[1] - collision.posY()) * v123.Py()) + ((secondaryVertex[2] - collision.posZ()) * v123.Pz()); + values[kVertexingLxyzProjected] = values[kVertexingLxyzProjected] / TMath::Sqrt((v123.Px() * v123.Px()) + (v123.Py() * v123.Py()) + (v123.Pz() * v123.Pz())); + + values[kVertexingTauzProjected] = values[kVertexingLzProjected] * v123.M() / TMath::Abs(v123.Pz()); + values[kVertexingTauxyProjected] = values[kVertexingLxyProjected] * v123.M() / (v123.Pt()); + values[kVertexingTauxyzProjected] = values[kVertexingLxyzProjected] * v123.M() / (v123.P()); + } +} + #endif // PWGDQ_CORE_VARMANAGER_H_