diff --git a/PWGDQ/Core/VarManager.h b/PWGDQ/Core/VarManager.h index a3b9d9d367b..99c049374e8 100644 --- a/PWGDQ/Core/VarManager.h +++ b/PWGDQ/Core/VarManager.h @@ -156,9 +156,11 @@ class VarManager : public TObject enum PairCandidateType { // TODO: need to agree on a scheme to incorporate all various hypotheses (e.g. e - mu, jpsi - K+, Jpsi - pipi,...) - kDecayToEE = 0, // e.g. J/psi -> e+ e- - kDecayToMuMu, // e.g. J/psi -> mu+ mu- - kDecayToPiPi, + kDecayToEE = 0, // e.g. J/psi -> e+ e- + kDecayToMuMu, // e.g. J/psi -> mu+ mu- + kDecayToPiPi, // decay channel to pi+ pi- + kDecayToKchKch, // decay channel to K+ K- + kDecayToPrAPr, // decay channel to p+ p- kElectronMuon, // e.g. Electron - muon correlations kBcToThreeMuons, // e.g. Bc -> mu+ mu- mu+ kBtoJpsiEEK, // e.g. B+ -> e+ e- K+ @@ -166,6 +168,7 @@ class VarManager : public TObject kXtoJpsiPiPi, // e.g. X(3872) -> J/psi pi+ pi- kPsi2StoJpsiPiPi, // e.g. Psi(2S) -> J/psi pi+ pi- kChictoJpsiEE, // e.g. Chi_c1 -> J/psi e+ e- + kNonResonantToKpKmPiPi, // e.g. J/psi -> K+ K- pi+ pi+ kDstarToD0KPiPi, // e.g. D*+ -> D0 pi+ -> K- pi+ pi+ kTripleCandidateToEEPhoton, // e.g. chi_c -> e+ e- photon or pi0 -> e+ e- photon kDecayToKPi, // e.g. D0 -> K+ pi- or cc. @@ -3841,6 +3844,16 @@ void VarManager::FillPair(T1 const& t1, T2 const& t2, float* values) m2 = o2::constants::physics::MassPionCharged; } + if constexpr (pairType == kDecayToKchKch) { + m1 = o2::constants::physics::MassKaonCharged; + m2 = o2::constants::physics::MassKaonCharged; + } + + if constexpr (pairType == kDecayToPrAPr) { + m1 = o2::constants::physics::MassProton; + m2 = o2::constants::physics::MassProton; + } + if constexpr (pairType == kDecayToKPi) { m1 = o2::constants::physics::MassKaonCharged; m2 = o2::constants::physics::MassPionCharged; @@ -6910,6 +6923,11 @@ void VarManager::FillDileptonTrackTrack(T1 const& dilepton, T2 const& hadron1, T hadronMass1 = o2::constants::physics::MassElectron; hadronMass2 = o2::constants::physics::MassElectron; } + if (candidateType == kNonResonantToKpKmPiPi) { + defaultDileptonMass = dilepton.mass(); + hadronMass1 = o2::constants::physics::MassPionCharged; + hadronMass2 = o2::constants::physics::MassPionCharged; + } ROOT::Math::PtEtaPhiMVector v1(dilepton.pt(), dilepton.eta(), dilepton.phi(), dilepton.mass()); ROOT::Math::PtEtaPhiMVector v2(hadron1.pt(), hadron1.eta(), hadron1.phi(), hadronMass1); diff --git a/PWGDQ/Tasks/tableReader_withAssoc.h b/PWGDQ/Tasks/tableReader_withAssoc.h index 354e9fe815d..1aff98249a4 100644 --- a/PWGDQ/Tasks/tableReader_withAssoc.h +++ b/PWGDQ/Tasks/tableReader_withAssoc.h @@ -815,12 +815,13 @@ struct AnalysisTrackSelection { fTrackCuts.push_back(o2::aod::dqcuts::GetCompositeCut(objArray->At(icut)->GetName())); } } + // Extra cuts via JSON TString addTrackCutsStr = fConfigCutsJSON.value; if (addTrackCutsStr != "") { std::vector addTrackCuts = o2::aod::dqcuts::GetCutsFromJSON(addTrackCutsStr.Data()); for (auto const& t : addTrackCuts) { - fTrackCuts.push_back(dynamic_cast(t)); + fTrackCuts.push_back(reinterpret_cast(t)); } } @@ -1573,7 +1574,7 @@ struct AnalysisSameEventPairing { void init(o2::framework::InitContext& context) { - fEnableBarrelHistos = context.mOptions.get("processAllSkimmed") || context.mOptions.get("processBarrelOnlySkimmed") || context.mOptions.get("processBarrelOnlyWithCollSkimmed") || context.mOptions.get("processBarrelOnlySkimmedNoCov") || context.mOptions.get("processBarrelOnlySkimmedNoCovWithMultExtra") || context.mOptions.get("processBarrelOnlyWithQvectorCentrSkimmedNoCov") || context.mOptions.get("processBarrelOnlyWithQvectorCentrSkimmed"); + fEnableBarrelHistos = context.mOptions.get("processAllSkimmed") || context.mOptions.get("processBarrelOnlySkimmed") || context.mOptions.get("processBarrelOnlyWithCollSkimmed") || context.mOptions.get("processBarrelOnlySkimmedNoCov") || context.mOptions.get("processBarrelOnlySkimmedNoCovWithMultExtra") || context.mOptions.get("processBarrelOnlyWithQvectorCentrSkimmedNoCov") || context.mOptions.get("processBarrelOnlyWithQvectorCentrSkimmed") || context.mOptions.get("processBarrelOnlyDiKaons") || context.mOptions.get("processBarrelOnlyDiProtons"); fEnableBarrelMixingHistos = context.mOptions.get("processMixingAllSkimmed") || context.mOptions.get("processMixingBarrelSkimmed") || context.mOptions.get("processMixingBarrelSkimmedFlow") || context.mOptions.get("processMixingBarrelWithQvectorCentrSkimmedNoCov"); fEnableBarrelMixingHistos |= fConfigRunMixingAcrossTFs; fEnableMuonHistos = context.mOptions.get("processAllSkimmed") || context.mOptions.get("processMuonOnlySkimmed") || context.mOptions.get("processMuonOnlySkimmedMultExtra") || context.mOptions.get("processMuonOnlySkimmedFlow") || context.mOptions.get("processMuonOnlyVertexingSkimmed"); @@ -2048,7 +2049,7 @@ struct AnalysisSameEventPairing { auto groupedAssocs = assocs.sliceBy(preslice, event.globalIndex()); size_t nGood = 0; for (auto const& t : groupedAssocs) { - if constexpr (TPairType == VarManager::kDecayToEE) { + if constexpr (TPairType == VarManager::kDecayToEE || TPairType == VarManager::kDecayToPiPi || TPairType == VarManager::kDecayToKchKch || TPairType == VarManager::kDecayToPrAPr) { if (t.isBarrelSelected_raw()) { nGood++; } @@ -2126,7 +2127,7 @@ struct AnalysisSameEventPairing { bool isFirst = true; for (auto const& [a1, a2] : o2::soa::combinations(groupedAssocs, groupedAssocs)) { - if constexpr (TPairType == VarManager::kDecayToEE) { + if constexpr (TPairType == VarManager::kDecayToEE || TPairType == VarManager::kDecayToPiPi || TPairType == VarManager::kDecayToKchKch || TPairType == VarManager::kDecayToPrAPr) { 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 @@ -2411,7 +2412,7 @@ struct AnalysisSameEventPairing { isUnambiguous = !(isAmbiInBunch || isAmbiOutOfBunch); isLeg1Ambi = (twoTrackFilter & (static_cast(1) << 28) || (twoTrackFilter & (static_cast(1) << 30))); isLeg2Ambi = (twoTrackFilter & (static_cast(1) << 29) || (twoTrackFilter & (static_cast(1) << 31))); - if constexpr (TPairType == VarManager::kDecayToEE) { + if constexpr (TPairType == VarManager::kDecayToEE || TPairType == VarManager::kDecayToPiPi || TPairType == VarManager::kDecayToKchKch || TPairType == VarManager::kDecayToPrAPr) { if (isLeg1Ambi && isLeg2Ambi) { std::pair iPair(a1.reducedtrackId(), a2.reducedtrackId()); if (fAmbiguousPairs.contains(iPair)) { @@ -2461,7 +2462,7 @@ struct AnalysisSameEventPairing { } } } - if constexpr (TPairType == VarManager::kDecayToEE) { + if constexpr (TPairType == VarManager::kDecayToEE || TPairType == VarManager::kDecayToPiPi || TPairType == VarManager::kDecayToKchKch || TPairType == VarManager::kDecayToPrAPr) { if (fConfigQA) { fHistMan->FillHistClass(histNames[icut][0].Data(), dqtablereader_helpers::varValues()); if (isAmbiExtra) { @@ -2485,7 +2486,7 @@ struct AnalysisSameEventPairing { } } } - if constexpr (TPairType == VarManager::kDecayToEE) { + if constexpr (TPairType == VarManager::kDecayToEE || TPairType == VarManager::kDecayToPiPi || TPairType == VarManager::kDecayToKchKch || TPairType == VarManager::kDecayToPrAPr) { fHistMan->FillHistClass(histNames[icut][1].Data(), dqtablereader_helpers::varValues()); if (isAmbiExtra) { fHistMan->FillHistClass(histNames[icut][4].Data(), dqtablereader_helpers::varValues()); @@ -2506,7 +2507,7 @@ struct AnalysisSameEventPairing { } } } - if constexpr (TPairType == VarManager::kDecayToEE) { + if constexpr (TPairType == VarManager::kDecayToEE || TPairType == VarManager::kDecayToPiPi || TPairType == VarManager::kDecayToKchKch || TPairType == VarManager::kDecayToPrAPr) { fHistMan->FillHistClass(histNames[icut][2].Data(), dqtablereader_helpers::varValues()); if (isAmbiExtra) { fHistMan->FillHistClass(histNames[icut][5].Data(), dqtablereader_helpers::varValues()); @@ -2537,7 +2538,7 @@ struct AnalysisSameEventPairing { // edit // rotation 20 times if (fConfigTRPairs) { - if constexpr (TPairType == VarManager::kDecayToEE) { + if constexpr (TPairType == VarManager::kDecayToEE || TPairType == VarManager::kDecayToPiPi || TPairType == VarManager::kDecayToKchKch || TPairType == VarManager::kDecayToPrAPr) { 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 @@ -2568,40 +2569,35 @@ struct AnalysisSameEventPairing { isUnambiguous = !(isAmbiInBunch || isAmbiOutOfBunch); isLeg1Ambi = (twoTrackFilter & (static_cast(1) << 28) || (twoTrackFilter & (static_cast(1) << 30))); isLeg2Ambi = (twoTrackFilter & (static_cast(1) << 29) || (twoTrackFilter & (static_cast(1) << 31))); - if constexpr (TPairType == VarManager::kDecayToEE) { - if (isLeg1Ambi && isLeg2Ambi) { - std::pair iPair(a1.reducedtrackId(), a2.reducedtrackId()); - if (fAmbiguousPairs.contains(iPair)) { - if (fAmbiguousPairs[iPair] & (static_cast(1) << icut)) { // if this pair is already stored with this cut - isAmbiExtra = true; - } else { - fAmbiguousPairs[iPair] |= static_cast(1) << icut; - } + if (isLeg1Ambi && isLeg2Ambi) { + std::pair iPair(a1.reducedtrackId(), a2.reducedtrackId()); + if (fAmbiguousPairs.contains(iPair)) { + if (fAmbiguousPairs[iPair] & (static_cast(1) << icut)) { // if this pair is already stored with this cut + isAmbiExtra = true; } else { - fAmbiguousPairs[iPair] = static_cast(1) << icut; + fAmbiguousPairs[iPair] |= static_cast(1) << icut; } + } else { + fAmbiguousPairs[iPair] = static_cast(1) << icut; } } + if (sign1 * sign2 < 0) { if (fConfigNRotations.value == 1) { VarManager::FillPairRotation(t1, t2, fConfigNRotations.value); - if constexpr (TPairType == VarManager::kDecayToEE) { - if (fConfigQA) { - fHistMan->FillHistClass(Form("PairsBarrelTRPM_%s", fTrackCuts[icut].Data()), dqtablereader_helpers::varValues()); - if (isAmbiExtra) { - fHistMan->FillHistClass(Form("PairsBarrelTRPM_ambiguousextra_%s", fTrackCuts[icut].Data()), dqtablereader_helpers::varValues()); - } + if (fConfigQA) { + fHistMan->FillHistClass(Form("PairsBarrelTRPM_%s", fTrackCuts[icut].Data()), dqtablereader_helpers::varValues()); + if (isAmbiExtra) { + fHistMan->FillHistClass(Form("PairsBarrelTRPM_ambiguousextra_%s", fTrackCuts[icut].Data()), dqtablereader_helpers::varValues()); } } } else if (fConfigNRotations.value == 3) { for (int irot = 1; irot <= fConfigNRotations.value; irot++) { VarManager::FillPairRotation(t1, t2, irot); - if constexpr (TPairType == VarManager::kDecayToEE) { - if (fConfigQA) { - fHistMan->FillHistClass(Form("PairsBarrelTRPM_%s", fTrackCuts[icut].Data()), dqtablereader_helpers::varValues()); - if (isAmbiExtra) { - fHistMan->FillHistClass(Form("PairsBarrelTRPM_ambiguousextra_%s", fTrackCuts[icut].Data()), dqtablereader_helpers::varValues()); - } + if (fConfigQA) { + fHistMan->FillHistClass(Form("PairsBarrelTRPM_%s", fTrackCuts[icut].Data()), dqtablereader_helpers::varValues()); + if (isAmbiExtra) { + fHistMan->FillHistClass(Form("PairsBarrelTRPM_ambiguousextra_%s", fTrackCuts[icut].Data()), dqtablereader_helpers::varValues()); } } } @@ -2629,7 +2625,7 @@ struct AnalysisSameEventPairing { MixingHandler::MixingEvent mixingEvent; uint32_t trackFilterForMixing = 0; for (auto const& assoc : groupedAssocs) { - if constexpr (TPairType == VarManager::kDecayToEE) { + if constexpr (TPairType == VarManager::kDecayToEE || TPairType == VarManager::kDecayToPiPi || TPairType == VarManager::kDecayToKchKch || TPairType == VarManager::kDecayToPrAPr) { trackFilterForMixing = assoc.isBarrelSelected_raw() & assoc.isBarrelSelectedPrefilter_raw() & fTrackFilterMask; if (!trackFilterForMixing) { continue; @@ -2761,7 +2757,7 @@ struct AnalysisSameEventPairing { auto twoTrackFilter = static_cast(0); for (auto const& a1 : assocs1) { for (auto const& a2 : assocs2) { - if constexpr (TPairType == VarManager::kDecayToEE) { + if constexpr (TPairType == VarManager::kDecayToEE || TPairType == VarManager::kDecayToPiPi || TPairType == VarManager::kDecayToKchKch || TPairType == VarManager::kDecayToPrAPr) { 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 continue; @@ -2876,7 +2872,7 @@ struct AnalysisSameEventPairing { } } } - if constexpr (TPairType == VarManager::kDecayToEE) { + if constexpr (TPairType == VarManager::kDecayToEE || TPairType == VarManager::kDecayToPiPi || TPairType == VarManager::kDecayToKchKch || TPairType == VarManager::kDecayToPrAPr) { fHistMan->FillHistClass(Form("PairsBarrelMEPM_%s", fTrackCuts[icut].Data()), dqtablereader_helpers::varValues()); } } else { @@ -2895,7 +2891,7 @@ struct AnalysisSameEventPairing { } } } - if constexpr (TPairType == VarManager::kDecayToEE) { + if constexpr (TPairType == VarManager::kDecayToEE || TPairType == VarManager::kDecayToPiPi || TPairType == VarManager::kDecayToKchKch || TPairType == VarManager::kDecayToPrAPr) { fHistMan->FillHistClass(Form("PairsBarrelMEPP_%s", fTrackCuts[icut].Data()), dqtablereader_helpers::varValues()); } } else { @@ -2913,7 +2909,7 @@ struct AnalysisSameEventPairing { } } } - if constexpr (TPairType == VarManager::kDecayToEE) { + if constexpr (TPairType == VarManager::kDecayToEE || TPairType == VarManager::kDecayToPiPi || TPairType == VarManager::kDecayToKchKch || TPairType == VarManager::kDecayToPrAPr) { fHistMan->FillHistClass(Form("PairsBarrelMEMM_%s", fTrackCuts[icut].Data()), dqtablereader_helpers::varValues()); } } @@ -3252,6 +3248,20 @@ struct AnalysisSameEventPairing { runSameEventPairing(events, trackAssocsPerCollision, barrelAssocs, barrelTracks); } + void processBarrelOnlyDiKaons(MyEventsVtxCovSelected const& events, + o2::soa::Join const& barrelAssocs, + MyBarrelTracksWithCovWithAmbiguities const& barrelTracks) + { + runSameEventPairing(events, trackAssocsPerCollision, barrelAssocs, barrelTracks); + } + + void processBarrelOnlyDiProtons(MyEventsVtxCovSelected const& events, + o2::soa::Join const& barrelAssocs, + MyBarrelTracksWithCovWithAmbiguities const& barrelTracks) + { + runSameEventPairing(events, trackAssocsPerCollision, barrelAssocs, barrelTracks); + } + void processMuonOnlySkimmed(MyEventsVtxCovSelected const& events, o2::soa::Join const& muonAssocs, MyMuonTracksWithCovWithAmbiguities const& muons) { @@ -3341,6 +3351,8 @@ struct AnalysisSameEventPairing { PROCESS_SWITCH(AnalysisSameEventPairing, processBarrelOnlyWithQvectorCentrSkimmedNoCov, "Run barrel only pairing (no covariances), with skimmed tracks, with Qvector from central framework", false); PROCESS_SWITCH(AnalysisSameEventPairing, processBarrelOnlyWithQvectorCentrSkimmed, "Run barrel only pairing (no covariances), with skimmed tracks, with collision information, with Qvector from central framework", false); PROCESS_SWITCH(AnalysisSameEventPairing, processBarrelOnlySkimmedFlow, "Run barrel only pairing, with skimmed tracks and with flow", false); + PROCESS_SWITCH(AnalysisSameEventPairing, processBarrelOnlyDiKaons, "Run barrel only pairing for di-kaons", false); + PROCESS_SWITCH(AnalysisSameEventPairing, processBarrelOnlyDiProtons, "Run barrel only pairing for di-protons", false); PROCESS_SWITCH(AnalysisSameEventPairing, processMuonOnlySkimmed, "Run muon only pairing, with skimmed tracks", false); PROCESS_SWITCH(AnalysisSameEventPairing, processMuonOnlySkimmedMultExtra, "Run muon only pairing, with skimmed tracks", false); PROCESS_SWITCH(AnalysisSameEventPairing, processMuonOnlySkimmedFlow, "Run muon only pairing, with skimmed tracks and flow", false); @@ -4645,6 +4657,7 @@ struct AnalysisDileptonTrack { o2::framework::Preslice dielectronsPerCollision = o2::aod::reducedpair::reducedeventId; o2::framework::Preslice ditracksPerCollision = o2::aod::reducedpair::reducedeventId; + // Charged B mesons to J/psi+Kch void processBarrelSkimmed(o2::soa::Filtered const& events, o2::soa::Filtered> const& assocs, MyBarrelTracksWithCov const& tracks, o2::soa::Filtered const& dileptons) @@ -4871,8 +4884,7 @@ struct AnalysisDileptonTrack { struct AnalysisDileptonTrackTrack { o2::framework::OutputObj fOutputList{"output"}; - o2::framework::Configurable fConfigTrackCut1{"cfgTrackCut1", "pionPIDCut1", "track1 cut"}; // used for select the tracks from SelectedTracks - o2::framework::Configurable fConfigTrackCut2{"cfgTrackCut2", "pionPIDCut2", "track2 cut"}; // used for select the tracks from SelectedTracks + o2::framework::Configurable fConfigTrackCuts{"cfgTrackCuts", "", "doublets of track cuts in the format cutName1:cutName2,cutName3:cutName4"}; // used for select the tracks from SelectedTracks o2::framework::Configurable fConfigDileptonCut{"cfgDiLeptonCut", "pairJpsi2", "Dilepton cut"}; o2::framework::Configurable fConfigQuadrupletCuts{"cfgQuadrupletCuts", "pairX3872Cut1", "Comma separated list of Dilepton-Track-Track cut"}; o2::framework::Configurable fConfigAddDileptonHistogram{"cfgAddDileptonHistogram", "barrel", "Comma separated list of histograms"}; @@ -4885,14 +4897,17 @@ struct AnalysisDileptonTrackTrack { o2::framework::Produces DileptonTrackTrackTable; int fCurrentRun = -1; // needed to detect if the run changed and trigger update of calibrations etc. - // uint32_t fTrackCutBitMap; // track cut bit mask to be used in the selection of tracks associated with dileptons - // cut name setting - TString fTrackCutName1; - TString fTrackCutName2; - bool fIsSameTrackCut = false; + uint32_t fTrackCutMask1 = 0; // track cut bit mask to be used in the selection of tracks + uint32_t fTrackCutMask2 = 0; // track cut bit mask to be used in the selection of tracks + uint32_t fNDoubletCuts = 0; // number of doublet track cuts + std::vector fDoubletCutNames; + std::vector> fDoubletCutMasks; + AnalysisCompositeCut fDileptonCut{}; + bool fApplyDileptonCut = false; std::vector fQuadrupletCutNames; std::vector fQuadrupletCuts; + uint32_t fNQuadrupletCuts = 0; // number of quadruplet track cuts o2::framework::Service fCCDB{}; @@ -4908,48 +4923,127 @@ struct AnalysisDileptonTrackTrack { void init(o2::framework::InitContext& context) { - // bool isBarrel = context.mOptions.get("processJpsiPiPi"); - if (context.mOptions.get("processDummy")) { return; } + // if the processKKPiPi function is enabled, then the analysis-same-event-pairing must have the processBarrelOnlyDiKaons enabled as well + bool processKKPiPi = context.mOptions.get("processKKPiPi"); + bool processBarrelOnlyDiKaons = false; + o2::common::core::getTaskOptionValue(context, "analysis-same-event-pairing", "processBarrelOnlyDiKaons", processBarrelOnlyDiKaons, false); + bool processDummy = false; + o2::common::core::getTaskOptionValue(context, "analysis-same-event-pairing", "processDummy", processDummy, false); + if (processKKPiPi) { + if (processDummy) { + LOGF(fatal, "processKKPiPi is enabled, but processDummy should not be enabled in analysis-same-event-pairing!"); + } + if (!processBarrelOnlyDiKaons) { + LOGF(fatal, "processKKPiPi is enabled, but processBarrelOnlyDiKaons is not enabled in analysis-same-event-pairing"); + } + } + // the cfgFlatTables option must be enabled in the analysis-same-event-pairing task because the dileptonsAll table has to be filled + bool flatTables = false; + o2::common::core::getTaskOptionValue(context, "analysis-same-event-pairing", "cfgFlatTables", flatTables, false); + if (!flatTables) { + LOGF(fatal, "Flat tables are not enabled in analysis-same-event-pairing! Enable the cfgFlatTables option!"); + } + + // check that there are track cuts defined otherwise throw a FATAL error + TString tempCutsStr = fConfigTrackCuts.value; + if (tempCutsStr.IsNull() || tempCutsStr.IsWhitespace()) { + LOGF(fatal, "No track cuts defined in analysis-track-selection task"); + } + VarManager::SetDefaultVarNames(); + + // get the list of track cuts computed in the AnalysisTrackSelection task + std::string tempCuts; + o2::common::core::getTaskOptionValue(context, "analysis-track-selection", "cfgTrackCuts", tempCuts, false); + tempCutsStr = tempCuts; + // check also the cuts added via JSON and add them to the string of cuts + o2::common::core::getTaskOptionValue(context, "analysis-track-selection", "cfgBarrelTrackCutsJSON", tempCuts, false); + TString addTrackCutsStr = tempCuts; + if (addTrackCutsStr != "") { + std::vector addTrackCuts = o2::aod::dqcuts::GetCutsFromJSON(addTrackCutsStr.Data()); + for (auto const& t : addTrackCuts) { + tempCutsStr += Form(",%s", t->GetName()); + } + } + std::unique_ptr trackSelectionsUpstream(tempCutsStr.Tokenize(",")); + + // interpret the list of track cuts and build the filter mask for each of the two prongs + TString trackCutsStr = fConfigTrackCuts.value; + std::unique_ptr doubletTrackSelections(trackCutsStr.Tokenize(",")); + for (Int_t icut = 0; icut < doubletTrackSelections->GetEntries(); ++icut) { + TString cutName = doubletTrackSelections->At(icut)->GetName(); + std::unique_ptr doubletCuts(cutName.Tokenize(":")); + if (doubletCuts->GetEntries() != 2) { + LOGF(fatal, "Doublet cut %s does not have exactly 2 prongs", cutName.Data()); + } + TObject* obj = trackSelectionsUpstream->FindObject(doubletCuts->At(0)->GetName()); + if (!obj) { + LOGF(fatal, "Track cut %s not found in upstream track selections", doubletCuts->At(0)->GetName()); + } + uint32_t mask1 = (static_cast(1) << trackSelectionsUpstream->IndexOf(obj)); + obj = trackSelectionsUpstream->FindObject(doubletCuts->At(1)->GetName()); + if (!obj) { + LOGF(fatal, "Track cut %s not found in upstream track selections", doubletCuts->At(1)->GetName()); + } + uint32_t mask2 = (static_cast(1) << trackSelectionsUpstream->IndexOf(obj)); + + fTrackCutMask1 |= mask1; + fTrackCutMask2 |= mask2; + fDoubletCutMasks.push_back(std::make_pair(mask1, mask2)); + fDoubletCutNames.push_back(cutName); + fNDoubletCuts++; + } + fCurrentRun = 0; fValuesQuadruplet = new float[VarManager::kNVars]; - VarManager::SetDefaultVarNames(); fHistMan = new HistogramManager("analysisHistos", "aa", VarManager::kNVars); fHistMan->SetUseDefaultVariableNames(kTRUE); fHistMan->SetDefaultVarNames(dqtablereader_helpers::varNames(), dqtablereader_helpers::varUnits()); - // define cuts - fTrackCutName1 = fConfigTrackCut1.value; - fTrackCutName2 = fConfigTrackCut2.value; - if (fTrackCutName1 == fTrackCutName2) { - fIsSameTrackCut = true; - } + // get the dilepton cut TString configDileptonCutNamesStr = fConfigDileptonCut.value; - fDileptonCut = *o2::aod::dqcuts::GetCompositeCut(configDileptonCutNamesStr.Data()); + if (!configDileptonCutNamesStr.IsNull()) { + fApplyDileptonCut = true; + fDileptonCut = *o2::aod::dqcuts::GetCompositeCut(configDileptonCutNamesStr.Data()); + } + TString configQuadruletCutNamesStr = fConfigQuadrupletCuts.value; std::unique_ptr objArray(configQuadruletCutNamesStr.Tokenize(",")); for (Int_t icut = 0; icut < objArray->GetEntries(); ++icut) { TString cutName = objArray->At(icut)->GetName(); fQuadrupletCutNames.push_back(cutName); fQuadrupletCuts.push_back(*o2::aod::dqcuts::GetCompositeCut(cutName.Data())); + fNQuadrupletCuts++; } + // define histogram directories + // One directory for the selected pairs (dileptons) if a cut was requested + // One directory for each combination of track doublet and quadruplet cuts if (!context.mOptions.get("processDummy")) { - DefineHistograms(fHistMan, Form("Pairs_%s", configDileptonCutNamesStr.Data()), fConfigAddDileptonHistogram.value.data()); - if (!configQuadruletCutNamesStr.IsNull()) { - for (std::size_t icut = 0; icut < fQuadrupletCutNames.size(); ++icut) { - if (fIsSameTrackCut) { - DefineHistograms(fHistMan, Form("QuadrupletSEPM_%s", fQuadrupletCutNames[icut].Data()), fConfigAddQuadrupletHistogram.value.data()); - } else { - DefineHistograms(fHistMan, Form("QuadrupletSEPM_%s", fQuadrupletCutNames[icut].Data()), fConfigAddQuadrupletHistogram.value.data()); - DefineHistograms(fHistMan, Form("QuadrupletSEMP_%s", fQuadrupletCutNames[icut].Data()), fConfigAddQuadrupletHistogram.value.data()); + if (fApplyDileptonCut) { + DefineHistograms(fHistMan, Form("Pairs_%s", configDileptonCutNamesStr.Data()), fConfigAddDileptonHistogram.value.data()); + } + for (uint32_t idcut = 0; idcut < fNDoubletCuts; ++idcut) { + for (uint32_t icut = 0; icut < fNQuadrupletCuts; ++icut) { + DefineHistograms(fHistMan, Form("QuadrupletSEPM_%s_%s", fDoubletCutNames[idcut].Data(), fQuadrupletCutNames[icut].Data()), fConfigAddQuadrupletHistogram.value.data()); + // if the doublet cut is asymmetric, add also the opposite ordering of the doublet cut + if (fDoubletCutMasks[idcut].first != fDoubletCutMasks[idcut].second) { + DefineHistograms(fHistMan, Form("QuadrupletSEMP_%s_%s", fDoubletCutNames[idcut].Data(), fQuadrupletCutNames[icut].Data()), fConfigAddQuadrupletHistogram.value.data()); } - DefineHistograms(fHistMan, Form("QuadrupletSEPP_%s", fQuadrupletCutNames[icut].Data()), fConfigAddQuadrupletHistogram.value.data()); - DefineHistograms(fHistMan, Form("QuadrupletSEMM_%s", fQuadrupletCutNames[icut].Data()), fConfigAddQuadrupletHistogram.value.data()); + DefineHistograms(fHistMan, Form("QuadrupletSEPP_%s_%s", fDoubletCutNames[idcut].Data(), fQuadrupletCutNames[icut].Data()), fConfigAddQuadrupletHistogram.value.data()); + DefineHistograms(fHistMan, Form("QuadrupletSEMM_%s_%s", fDoubletCutNames[idcut].Data(), fQuadrupletCutNames[icut].Data()), fConfigAddQuadrupletHistogram.value.data()); + } + if (fNQuadrupletCuts == 0) { + DefineHistograms(fHistMan, Form("QuadrupletSEPM_%s", fDoubletCutNames[idcut].Data()), fConfigAddQuadrupletHistogram.value.data()); + if (fDoubletCutMasks[idcut].first != fDoubletCutMasks[idcut].second) { + DefineHistograms(fHistMan, Form("QuadrupletSEMP_%s", fDoubletCutNames[idcut].Data()), fConfigAddQuadrupletHistogram.value.data()); + } + DefineHistograms(fHistMan, Form("QuadrupletSEPP_%s", fDoubletCutNames[idcut].Data()), fConfigAddQuadrupletHistogram.value.data()); + DefineHistograms(fHistMan, Form("QuadrupletSEMM_%s", fDoubletCutNames[idcut].Data()), fConfigAddQuadrupletHistogram.value.data()); } } } @@ -4983,12 +5077,9 @@ struct AnalysisDileptonTrackTrack { { VarManager::ResetValues(0, VarManager::kNVars, fValuesQuadruplet); VarManager::FillEvent(event, fValuesQuadruplet); - // int indexOffset = -999; - // std::vector trackGlobalIndexes; for (auto const& dilepton : dileptons) { // get full track info of tracks based on the index - int indexLepton1 = dilepton.index0Id(); int indexLepton2 = dilepton.index1Id(); auto lepton1 = tracks.rawIteratorAt(dilepton.index0Id()); @@ -4998,25 +5089,24 @@ struct AnalysisDileptonTrackTrack { continue; } VarManager::FillTrack(dilepton, fValuesQuadruplet); - // LOGP(info, "is dilepton selected: {}", fDileptonCut.IsSelected(fValuesQuadruplet)); // apply the dilepton cut - if (!fDileptonCut.IsSelected(fValuesQuadruplet)) { + if (fApplyDileptonCut && !fDileptonCut.IsSelected(fValuesQuadruplet)) { continue; } - fHistMan->FillHistClass(Form("Pairs_%s", fDileptonCut.GetName()), fValuesQuadruplet); + if (fApplyDileptonCut) { + fHistMan->FillHistClass(Form("Pairs_%s", fDileptonCut.GetName()), fValuesQuadruplet); + } - // loop over hadrons pairs + // loop over track pairs for (auto const& [a1, a2] : o2::soa::combinations(assocs, assocs)) { - uint32_t trackSelection = 0; - if (fIsSameTrackCut) { - trackSelection = ((a1.isBarrelSelected_raw() & (static_cast(1) << 1)) || (a2.isBarrelSelected_raw() & (static_cast(1) << 1))); - } else { - trackSelection = ((a1.isBarrelSelected_raw() & (static_cast(1) << 1)) && (a2.isBarrelSelected_raw() & (static_cast(1) << 2))); - } - // LOGP(info, "trackSelection: {}, a1: {}, a2: {}", trackSelection, a1.isBarrelSelected_raw(), a2.isBarrelSelected_raw()); - if (trackSelection == 0) { + uint32_t trackSelection1 = 0; + uint32_t trackSelection2 = 0; + trackSelection1 = (a1.isBarrelSelected_raw() & fTrackCutMask1) | (a2.isBarrelSelected_raw() & fTrackCutMask1); + trackSelection2 = (a2.isBarrelSelected_raw() & fTrackCutMask2) | (a1.isBarrelSelected_raw() & fTrackCutMask2); + // skip if no tracks were selected + if (trackSelection1 == 0 || trackSelection2 == 0) { continue; } @@ -5031,39 +5121,81 @@ struct AnalysisDileptonTrackTrack { // fill variables VarManager::FillDileptonTrackTrack(dilepton, track1, track2, fValuesQuadruplet); if (fConfigSetupFourProngFitter || fConfigUseKFVertexing) { - // LOGP(info, "Using KF or DCA fitter for secondary vertexing"); VarManager::FillDileptonTrackTrackVertexing(event, lepton1, lepton2, track1, track2, fValuesQuadruplet); } - int iCut = 0; + // first loop over the doublet track cuts and make a decision map + int idCut = 0; + uint32_t idDoubletCutDecision = 0; + for (auto cutmask = fDoubletCutMasks.begin(); cutmask != fDoubletCutMasks.end(); cutmask++, idCut++) { + // check the doublet cuts + if ((trackSelection1 & cutmask->first) && (trackSelection2 & cutmask->second)) { + idDoubletCutDecision |= (static_cast(1) << idCut); + } + // check the reverse double cuts + if ((trackSelection2 & cutmask->first) && (trackSelection1 & cutmask->second)) { + idDoubletCutDecision |= (static_cast(1) << idCut); + } + } + uint32_t CutDecision = 0; - for (auto cutname = fQuadrupletCutNames.begin(); cutname != fQuadrupletCutNames.end(); cutname++, iCut++) { + for (uint32_t iCut = 0; iCut < fNQuadrupletCuts; iCut++) { // apply dilepton-track-track cut if (fQuadrupletCuts[iCut].IsSelected(fValuesQuadruplet)) { - CutDecision |= (1 << iCut); - if (fIsSameTrackCut) { - if (track1.sign() * track2.sign() < 0) { - fHistMan->FillHistClass(Form("QuadrupletSEPM_%s", fQuadrupletCutNames[iCut].Data()), fValuesQuadruplet); + CutDecision |= (static_cast(1) << iCut); + } + } + if (fNQuadrupletCuts > 0 && !CutDecision) { + continue; + } + + // now make a double loop over doublet and quadruplet cuts and fill histograms + for (uint32_t iDoubletCut = 0; iDoubletCut < fNDoubletCuts; iDoubletCut++) { + if (!(idDoubletCutDecision & (static_cast(1) << iDoubletCut))) { + continue; + } + // check if the prong1 and prong2 cuts are the same + bool isSameTrackCut = fDoubletCutMasks[iDoubletCut].first == fDoubletCutMasks[iDoubletCut].second; + for (uint32_t iQuadrupletCut = 0; iQuadrupletCut < fNQuadrupletCuts; iQuadrupletCut++) { + if (CutDecision & (static_cast(1) << iQuadrupletCut)) { + if (isSameTrackCut) { + if (track1.sign() * track2.sign() < 0) { + fHistMan->FillHistClass(Form("QuadrupletSEPM_%s_%s", fDoubletCutNames[iDoubletCut].Data(), fQuadrupletCutNames[iQuadrupletCut].Data()), fValuesQuadruplet); + } + } else { + if ((track1.sign() < 0) && (track2.sign() > 0)) { + fHistMan->FillHistClass(Form("QuadrupletSEMP_%s_%s", fDoubletCutNames[iDoubletCut].Data(), fQuadrupletCutNames[iQuadrupletCut].Data()), fValuesQuadruplet); + } else if ((track1.sign() > 0) && (track2.sign() < 0)) { + fHistMan->FillHistClass(Form("QuadrupletSEPM_%s_%s", fDoubletCutNames[iDoubletCut].Data(), fQuadrupletCutNames[iQuadrupletCut].Data()), fValuesQuadruplet); + } } - } else { - if ((track1.sign() < 0) && (track2.sign() > 0)) { - fHistMan->FillHistClass(Form("QuadrupletSEMP_%s", fQuadrupletCutNames[iCut].Data()), fValuesQuadruplet); - } else if ((track1.sign() > 0) && (track2.sign() < 0)) { - fHistMan->FillHistClass(Form("QuadrupletSEPM_%s", fQuadrupletCutNames[iCut].Data()), fValuesQuadruplet); + if ((track1.sign() > 0) && (track2.sign() > 0)) { + fHistMan->FillHistClass(Form("QuadrupletSEPP_%s_%s", fDoubletCutNames[iDoubletCut].Data(), fQuadrupletCutNames[iQuadrupletCut].Data()), fValuesQuadruplet); + } else if ((track1.sign() < 0) && (track2.sign() < 0)) { + fHistMan->FillHistClass(Form("QuadrupletSEMM_%s_%s", fDoubletCutNames[iDoubletCut].Data(), fQuadrupletCutNames[iQuadrupletCut].Data()), fValuesQuadruplet); } } - if ((track1.sign() > 0) && (track2.sign() > 0)) { - fHistMan->FillHistClass(Form("QuadrupletSEPP_%s", fQuadrupletCutNames[iCut].Data()), fValuesQuadruplet); - } else if ((track1.sign() < 0) && (track2.sign() < 0)) { - fHistMan->FillHistClass(Form("QuadrupletSEMM_%s", fQuadrupletCutNames[iCut].Data()), fValuesQuadruplet); + } // end of loop over quadruplet cuts + // in case there are no quadruplet cuts specified by the user + if (isSameTrackCut) { + if (track1.sign() * track2.sign() < 0) { + fHistMan->FillHistClass(Form("QuadrupletSEPM_%s", fDoubletCutNames[iDoubletCut].Data()), fValuesQuadruplet); + } + } else { + if ((track1.sign() < 0) && (track2.sign() > 0)) { + fHistMan->FillHistClass(Form("QuadrupletSEMP_%s", fDoubletCutNames[iDoubletCut].Data()), fValuesQuadruplet); + } else if ((track1.sign() > 0) && (track2.sign() < 0)) { + fHistMan->FillHistClass(Form("QuadrupletSEPM_%s", fDoubletCutNames[iDoubletCut].Data()), fValuesQuadruplet); } } - } // loop over dilepton-track-track cuts + if ((track1.sign() > 0) && (track2.sign() > 0)) { + fHistMan->FillHistClass(Form("QuadrupletSEPP_%s", fDoubletCutNames[iDoubletCut].Data()), fValuesQuadruplet); + } else if ((track1.sign() < 0) && (track2.sign() < 0)) { + fHistMan->FillHistClass(Form("QuadrupletSEMM_%s", fDoubletCutNames[iDoubletCut].Data()), fValuesQuadruplet); + } + } // end of loop over doublet cuts // fill table - if (!CutDecision) { - continue; - } DileptonTrackTrackTable(fValuesQuadruplet[VarManager::kQuadMass], fValuesQuadruplet[VarManager::kQuadPt], fValuesQuadruplet[VarManager::kQuadEta], fValuesQuadruplet[VarManager::kQuadPhi], fValuesQuadruplet[VarManager::kRap], fValuesQuadruplet[VarManager::kQ], fValuesQuadruplet[VarManager::kDeltaR1], fValuesQuadruplet[VarManager::kDeltaR2], fValuesQuadruplet[VarManager::kDeltaR], dilepton.mass(), dilepton.pt(), dilepton.eta(), dilepton.phi(), dilepton.sign(), @@ -5099,12 +5231,31 @@ struct AnalysisDileptonTrackTrack { } } + void processKKPiPi(o2::soa::Filtered const& events, + o2::soa::Filtered> const& assocs, + MyBarrelTracksWithCov const& tracks, o2::soa::Filtered const& dileptons) + { + if (events.size() == 0) { + return; + } + if (fCurrentRun != events.begin().runNumber()) { // start: runNumber + initParamsFromCCDB(events.begin().timestamp()); + fCurrentRun = events.begin().runNumber(); + } // end: runNumber + for (auto const& event : events) { + auto groupedBarrelAssocs = assocs.sliceBy(trackAssocsPerCollision, event.globalIndex()); + auto groupedDielectrons = dileptons.sliceBy(dielectronsPerCollision, event.globalIndex()); + runDileptonTrackTrack(event, groupedBarrelAssocs, tracks, groupedDielectrons); + } + } + void processDummy(MyEventsBasic const&) { // do nothing } PROCESS_SWITCH(AnalysisDileptonTrackTrack, processJpsiPiPi, "Run barrel pairing of J/psi with pion candidate", false); + PROCESS_SWITCH(AnalysisDileptonTrackTrack, processKKPiPi, "Run barrel pairing of non-resonant K+ K- pi+ pi+ candidate", false); PROCESS_SWITCH(AnalysisDileptonTrackTrack, processDummy, "Dummy function", true); };