diff --git a/PWGDQ/Core/HistogramsLibrary.cxx b/PWGDQ/Core/HistogramsLibrary.cxx index df3bdc0470b..33b129186d6 100644 --- a/PWGDQ/Core/HistogramsLibrary.cxx +++ b/PWGDQ/Core/HistogramsLibrary.cxx @@ -980,25 +980,24 @@ void o2::aod::dqhistograms::DefineHistograms(HistogramManager* hm, const char* h } if (!groupStr.CompareTo("energy-correlator-gen")) { - hm->AddHistogram(histClass, "MCCostheta", "Cos#theta", false, 40, -1.0, 1.0, VarManager::kMCCosTheta, 0, 0, 0, -1, 0, 0, 0, -1, "", "", "", -1, VarManager::kMCWeight_before); - hm->AddHistogram(histClass, "MCHadronPdgCode", "MCHadronPdgCode", false, 6000, -3000, 3000, VarManager::kMCHadronPdgCode); - hm->AddHistogram(histClass, "MCMotherPdgCode", "MCMotherPdgCode", false, 6000, -3000, 3000, VarManager::kMCMotherPdgCode); - hm->AddHistogram(histClass, "MCPdgCode", "MCPdgCode", false, 1000, -1000, 1000, VarManager::kMCPdgCode); - hm->AddHistogram(histClass, "Coschi", "", false, 40, -1.0, 1.0, VarManager::kMCCosChi, 0, 0, 0, -1, 0, 0, 0, -1, "", "", "", -1, VarManager::kMCWeight); - hm->AddHistogram(histClass, "Pt_Hadron", "", false, 120, 0.0, 30.0, VarManager::kMCHadronPt); - hm->AddHistogram(histClass, "Eta_Hadron", "", false, 120, -2.0, 2.0, VarManager::kMCHadronEta); - hm->AddHistogram(histClass, "Phi_Hadron", "", false, 120, -2.0, 2.0, VarManager::kMCHadronPhi); - hm->AddHistogram(histClass, "DeltaEta", "", false, 20, -2.0, 2.0, VarManager::kMCdeltaeta); - hm->AddHistogram(histClass, "DeltaPhi", "", false, 50, -8.0, 8.0, VarManager::kMCdeltaphi); - hm->AddHistogram(histClass, "DeltaEta_DeltaPhi", "", false, 20, -2.0, 2.0, VarManager::kMCdeltaeta, 50, -8.0, 8.0, VarManager::kMCdeltaphi); - // for bkg - hm->AddHistogram(histClass, "DeltaPhi_randomPhi_trans", "", false, 50, -8.0, 8.0, VarManager::kMCdeltaphi_randomPhi_trans); - hm->AddHistogram(histClass, "DeltaPhi_randomPhi_toward", "", false, 50, -8.0, 8.0, VarManager::kMCdeltaphi_randomPhi_toward); - hm->AddHistogram(histClass, "DeltaPhi_randomPhi_away", "", false, 50, -8.0, 8.0, VarManager::kMCdeltaphi_randomPhi_away); + double coschiBins[26]; + for (int i = 0; i < 26; i++) { + coschiBins[i] = -1.0 + 2.0 * TMath::Power(0.04 * i, 2.0); + } - hm->AddHistogram(histClass, "Coschi_randomPhi_trans", "", false, 40, -1.0, 1.0, VarManager::kMCCosChi_randomPhi_trans, 0, 0, 0, -1, 0, 0, 0, -1, "", "", "", -1, VarManager::kMCWeight_randomPhi_trans); - hm->AddHistogram(histClass, "Coschi_randomPhi_toward", "", false, 40, -1.0, 1.0, VarManager::kMCCosChi_randomPhi_toward, 0, 0, 0, -1, 0, 0, 0, -1, "", "", "", -1, VarManager::kMCWeight_randomPhi_toward); - hm->AddHistogram(histClass, "Coschi_randomPhi_away", "", false, 40, -1.0, 1.0, VarManager::kMCCosChi_randomPhi_away, 0, 0, 0, -1, 0, 0, 0, -1, "", "", "", -1, VarManager::kMCWeight_randomPhi_away); + hm->AddHistogram(histClass, "Coschi", "", false, 25, coschiBins, VarManager::kMCCosChi, 0, nullptr, -1, 0, nullptr, -1, "", "", "", -1, VarManager::kMCWeight); + + // for bkg + hm->AddHistogram(histClass, "Coschi_randomPhi_trans", "", false, 25, coschiBins, VarManager::kMCCosChi_randomPhi_trans, 0, nullptr, -1, 0, nullptr, -1, "", "", "", -1, VarManager::kMCWeight_randomPhi_trans); + hm->AddHistogram(histClass, "Coschi_randomPhi_toward", "", false, 25, coschiBins, VarManager::kMCCosChi_randomPhi_toward, 0, nullptr, -1, 0, nullptr, -1, "", "", "", -1, VarManager::kMCWeight_randomPhi_toward); + hm->AddHistogram(histClass, "Coschi_randomPhi_away", "", false, 25, coschiBins, VarManager::kMCCosChi_randomPhi_away, 0, nullptr, -1, 0, nullptr, -1, "", "", "", -1, VarManager::kMCWeight_randomPhi_away); + } + if (!groupStr.CompareTo("energy-correlator-unfolding")) { + double coschiBins[26]; + for (int i = 0; i < 26; i++) { + coschiBins[i] = -1.0 + 0.08 * i; + } + hm->AddHistogram(histClass, "Coschi_unfolding", "", false, 25, coschiBins, VarManager::kMCCosChi_rec, 25, coschiBins, VarManager::kMCCosChi_gen); } if (!groupStr.CompareTo("polarization-pseudoproper-gen")) { int varspTHE[3] = {VarManager::kMCPt, VarManager::kMCCosThetaHE, VarManager::kMCVertexingTauxyProjected}; @@ -1952,25 +1951,8 @@ void o2::aod::dqhistograms::DefineHistograms(HistogramManager* hm, const char* h for (int i = 0; i < 21; i++) { deltaetaBins[i] = -2.0 + 0.2 * i; } - - hm->AddHistogram(histClass, "dileptonmass", "", false, 125, 2.5, 3.5, VarManager::kdileptonmass); hm->AddHistogram(histClass, "Coschi", "", false, 25, coschiBins, VarManager::kCosChi, 0, nullptr, -1, 0, nullptr, -1, "", "", "", -1, VarManager::kECWeight); - hm->AddHistogram(histClass, "Pt_Hadron", ";P_{T}", false, 120, 0.0, 30.0, VarManager::kPtDau); - hm->AddHistogram(histClass, "Coschi_wo", "", false, 25, coschiBins, VarManager::kCosChi); - hm->AddHistogram(histClass, "Eta_Hadron", ";#eta", false, 120, -2.0, 2.0, VarManager::kEtaDau); - hm->AddHistogram(histClass, "Phi_Hadron", ";#phi", false, 120, -8, 8, VarManager::kPhiDau); - hm->AddHistogram(histClass, "DeltaEta_DeltaPhi_weight", "", false, 20, -2.0, 2.0, VarManager::kDeltaEta, 50, -2.0, 6.0, VarManager::kDeltaPhi, 0, 0, 0, -1, "", "", "", -1, VarManager::kEWeight_before); - hm->AddHistogram(histClass, "DeltaEta_DeltaPhi", "", false, 20, -2.0, 2.0, VarManager::kDeltaEta, 50, -2.0, 6.0, VarManager::kDeltaPhi); - hm->AddHistogram(histClass, "Coschi_DeltaEta", "", false, 25, coschiBins, VarManager::kCosChi, 20, deltaetaBins, VarManager::kDeltaEta, 0, nullptr, -1, "", "", "", -1, VarManager::kECWeight); - hm->AddHistogram(histClass, "Coschi_wo_DeltaEta", "", false, 25, coschiBins, VarManager::kCosChi, 20, deltaetaBins, VarManager::kDeltaEta); - // for bkg - hm->AddHistogram(histClass, "Coschi_randomPhi_trans", "", false, 25, coschiBins, VarManager::kCosChi_randomPhi_trans, 0, nullptr, -1, 0, nullptr, -1, "", "", "", -1, VarManager::kWeight_randomPhi_trans); - hm->AddHistogram(histClass, "Coschi_randomPhi_toward", "", false, 25, coschiBins, VarManager::kCosChi_randomPhi_toward, 0, nullptr, -1, 0, nullptr, -1, "", "", "", -1, VarManager::kWeight_randomPhi_toward); - hm->AddHistogram(histClass, "Coschi_randomPhi_away", "", false, 25, coschiBins, VarManager::kCosChi_randomPhi_away, 0, nullptr, -1, 0, nullptr, -1, "", "", "", -1, VarManager::kWeight_randomPhi_away); - - hm->AddHistogram(histClass, "Coschi_wo_randomPhi_trans", "", false, 25, coschiBins, VarManager::kCosChi_randomPhi_trans); - hm->AddHistogram(histClass, "Coschi_wo_randomPhi_toward", "", false, 25, coschiBins, VarManager::kCosChi_randomPhi_toward); - hm->AddHistogram(histClass, "Coschi_wo_randomPhi_away", "", false, 25, coschiBins, VarManager::kCosChi_randomPhi_away); + hm->AddHistogram(histClass, "DeltaEta_DeltaPhi_weight", "", false, 20, -2.0, 2.0, VarManager::kDeltaEta, 50, -2.0, 6.0, VarManager::kDeltaPhi, 0, 0, 0, -1, "", "", "", -1, VarManager::kPtDau); } } diff --git a/PWGDQ/Core/VarManager.cxx b/PWGDQ/Core/VarManager.cxx index ecbec0dc442..e69d0794ece 100644 --- a/PWGDQ/Core/VarManager.cxx +++ b/PWGDQ/Core/VarManager.cxx @@ -1824,7 +1824,29 @@ void VarManager::SetDefaultVarNames() fgVarNamesMap["kMCCosChi"] = kMCCosChi; fgVarNamesMap["kMCHadronPt"] = kMCHadronPt; fgVarNamesMap["kMCWeight_before"] = kMCWeight_before; + fgVarNamesMap["kMCdeltaeta"] = kMCdeltaeta; + fgVarNamesMap["kMCHadronPt"] = kMCHadronPt; + fgVarNamesMap["kMCHadronEta"] = kMCHadronEta; + fgVarNamesMap["kMCHadronPhi"] = kMCHadronPhi; + fgVarNamesMap["kMCWeight"] = kMCWeight; + fgVarNamesMap["kMCCosChi_randomPhi_toward"] = kMCCosChi_randomPhi_toward; + fgVarNamesMap["kMCWeight_randomPhi_toward"] = kMCWeight_randomPhi_toward; + fgVarNamesMap["kMCCosChi_randomPhi_away"] = kMCCosChi_randomPhi_away; + fgVarNamesMap["kMCWeight_randomPhi_away"] = kMCWeight_randomPhi_away; + fgVarNamesMap["kMCCosChi_randomPhi_trans"] = kMCCosChi_randomPhi_trans; + fgVarNamesMap["kMCWeight_randomPhi_trans"] = kMCWeight_randomPhi_trans; + fgVarNamesMap["kMCdeltaphi_randomPhi_toward"] = kMCdeltaphi_randomPhi_toward; + fgVarNamesMap["kMCdeltaphi_randomPhi_away"] = kMCdeltaphi_randomPhi_away; + fgVarNamesMap["kMCdeltaphi_randomPhi_trans"] = kMCdeltaphi_randomPhi_trans; + fgVarNamesMap["kMCCosChi_gen"] = kMCCosChi_gen; + fgVarNamesMap["kMCWeight_gen"] = kMCWeight_gen; + fgVarNamesMap["kMCdeltaeta_gen"] = kMCdeltaeta_gen; + fgVarNamesMap["kMCCosChi_rec"] = kMCCosChi_rec; + fgVarNamesMap["kMCWeight_rec"] = kMCWeight_rec; + fgVarNamesMap["kMCdeltaeta_rec"] = kMCdeltaeta_rec; fgVarNamesMap["kMCParticleWeight"] = kMCParticleWeight; + fgVarNamesMap["kMCCosTheta"] = kMCCosTheta; + fgVarNamesMap["kMCJpsiPt"] = kMCJpsiPt; fgVarNamesMap["kMCPx"] = kMCPx; fgVarNamesMap["kMCPy"] = kMCPy; fgVarNamesMap["kMCPz"] = kMCPz; @@ -2045,6 +2067,16 @@ void VarManager::SetDefaultVarNames() fgVarNamesMap["kPtDau"] = kPtDau; fgVarNamesMap["kEtaDau"] = kEtaDau; fgVarNamesMap["kPhiDau"] = kPhiDau; + fgVarNamesMap["kCosChi_randomPhi_trans"] = kCosChi_randomPhi_trans; + fgVarNamesMap["kCosChi_randomPhi_toward"] = kCosChi_randomPhi_toward; + fgVarNamesMap["kCosChi_randomPhi_away"] = kCosChi_randomPhi_away; + fgVarNamesMap["kWeight_randomPhi_trans"] = kWeight_randomPhi_trans; + fgVarNamesMap["kWeight_randomPhi_toward"] = kWeight_randomPhi_toward; + fgVarNamesMap["kWeight_randomPhi_away"] = kWeight_randomPhi_away; + fgVarNamesMap["kdeltaphi_randomPhi_trans"] = kdeltaphi_randomPhi_trans; + fgVarNamesMap["kdeltaphi_randomPhi_toward"] = kdeltaphi_randomPhi_toward; + fgVarNamesMap["kdeltaphi_randomPhi_away"] = kdeltaphi_randomPhi_away; + fgVarNamesMap["kdileptonmass"] = kdileptonmass; fgVarNamesMap["kNCorrelationVariables"] = kNCorrelationVariables; fgVarNamesMap["kQuadMass"] = kQuadMass; fgVarNamesMap["kQuadDefaultDileptonMass"] = kQuadDefaultDileptonMass; diff --git a/PWGDQ/Core/VarManager.h b/PWGDQ/Core/VarManager.h index 6097dcf3ef8..341e24d6829 100644 --- a/PWGDQ/Core/VarManager.h +++ b/PWGDQ/Core/VarManager.h @@ -655,6 +655,12 @@ class VarManager : public TObject kMCdeltaphi_randomPhi_away, kMCdeltaphi_randomPhi_trans, kMCWeight_before, + kMCCosChi_gen, + kMCWeight_gen, + kMCdeltaeta_gen, + kMCCosChi_rec, + kMCWeight_rec, + kMCdeltaeta_rec, // MC mother particle variables kMCMotherPdgCode, @@ -1208,6 +1214,8 @@ class VarManager : public TObject static void FillTrackMC(const U& mcStack, T const& track, float* values = nullptr); template static void FillEnergyCorrelatorsMC(T const& track, T1 const& t1, float* values = nullptr); + template + static void FillEnergyCorrelatorsMCUnfolding(T1 const& dilepton, T2 const& hadron, T const& track, T3 const& t1, float* values = nullptr); template static void FillPairPropagateMuon(T1 const& muon1, T2 const& muon2, const C& collision, float* values = nullptr); template @@ -3004,6 +3012,33 @@ void VarManager::FillEnergyCorrelatorsMC(T const& track, T1 const& t1, float* va } } +template +void VarManager::FillEnergyCorrelatorsMCUnfolding(T1 const& dilepton, T2 const& hadron, T const& track, T3 const& t1, float* values) +{ + // energy correlators + float MassHadron; + if constexpr (pairType == kJpsiHadronMass) { + MassHadron = TMath::Sqrt(t1.e() * t1.e() - t1.p() * t1.p()); + } + if constexpr (pairType == kJpsiPionMass) { + MassHadron = o2::constants::physics::MassPionCharged; + } + ROOT::Math::PtEtaPhiMVector v1_gen(track.pt(), track.eta(), track.phi(), o2::constants::physics::MassJPsi); + ROOT::Math::PtEtaPhiMVector v2_gen(t1.pt(), t1.eta(), t1.phi(), MassHadron); + float E_boost_gen = LorentzTransformJpsihadroncosChi("weight_boost", v1_gen, v2_gen); + float CosChi_gen = LorentzTransformJpsihadroncosChi("coschi", v1_gen, v2_gen); + values[kMCCosChi_gen] = CosChi_gen; + values[kMCWeight_gen] = E_boost_gen / o2::constants::physics::MassJPsi; + values[kMCdeltaeta_gen] = track.eta() - t1.eta(); + + ROOT::Math::PtEtaPhiMVector v1_rec(dilepton.pt(), dilepton.eta(), dilepton.phi(), dilepton.mass()); + ROOT::Math::PtEtaPhiMVector v2_rec(hadron.pt(), hadron.eta(), hadron.phi(), o2::constants::physics::MassPionCharged); + values[kMCCosChi_rec] = LorentzTransformJpsihadroncosChi("coschi", v1_rec, v2_rec); + float E_boost_rec = LorentzTransformJpsihadroncosChi("weight_boost", v1_rec, v2_rec); + values[kMCWeight_rec] = E_boost_rec / v1_rec.M(); + values[kMCdeltaeta_rec] = dilepton.eta() - hadron.eta(); +} + template void VarManager::FillPairPropagateMuon(T1 const& muon1, T2 const& muon2, const C& collision, float* values) { diff --git a/PWGDQ/Tasks/dqEfficiency_withAssoc_direct.cxx b/PWGDQ/Tasks/dqEfficiency_withAssoc_direct.cxx index 44656cd7476..d9f1aa16061 100644 --- a/PWGDQ/Tasks/dqEfficiency_withAssoc_direct.cxx +++ b/PWGDQ/Tasks/dqEfficiency_withAssoc_direct.cxx @@ -2065,6 +2065,7 @@ struct AnalysisDileptonTrack { Configurable fConfigAddJSONHistograms{"cfgAddJSONHistograms", "", "Histograms in JSON format"}; Configurable fConfigMixingDepth{"cfgMixingDepth", 5, "Event mixing pool depth"}; Configurable fConfigPublishTripletTable{"cfgPublishTripletTable", false, "Publish the triplet tables, BmesonCandidates"}; + Configurable fConfigApplyMassEC{"cfgApplyMassEC", false, "Apply fit mass for sideband for the energy correlator study"}; } fConfigOptions; struct : ConfigurableGroup { @@ -2137,6 +2138,7 @@ struct AnalysisDileptonTrack { bool isDummy = context.mOptions.get("processDummy"); bool isMCGen_energycorrelators = context.mOptions.get("processMCGenEnergyCorrelators") || context.mOptions.get("processMCGenEnergyCorrelatorsPion"); bool isMCGen_energycorrelatorsME = context.mOptions.get("processMCGenEnergyCorrelatorsME") || context.mOptions.get("processMCGenEnergyCorrelatorsPionME"); + bool isMC_energycorrelatorsUnfolding = context.mOptions.get("processMCEnergyCorrelatorsUnfolding"); if (isDummy) { if (isBarrel || isMCGen /*|| isBarrelAsymmetric*/ /*|| isMuon*/) { @@ -2393,6 +2395,9 @@ struct AnalysisDileptonTrack { DefineHistograms(fHistMan, Form("DileptonTrack_%s_%s", pairLegCutName.Data(), fTrackCutNames[iCutTrack].Data()), fConfigOptions.fConfigHistogramSubgroups.value.data()); for (auto& sig : fRecMCSignals) { DefineHistograms(fHistMan, Form("DileptonTrackMCMatched_%s_%s_%s", pairLegCutName.Data(), fTrackCutNames[iCutTrack].Data(), sig->GetName()), fConfigOptions.fConfigHistogramSubgroups.value.data()); + if (isMC_energycorrelatorsUnfolding) { + DefineHistograms(fHistMan, Form("DileptonTrackMCUnfold_%s_%s_%s", pairLegCutName.Data(), fTrackCutNames[iCutTrack].Data(), sig->GetName()), ""); + } } if (!cfgPairing_strCommonTrackCuts.IsNull()) { @@ -2580,7 +2585,8 @@ struct AnalysisDileptonTrack { // compute needed quantities VarManager::FillDileptonHadron(dilepton, track, fValuesHadron); VarManager::FillDileptonTrackVertexing(event, lepton1, lepton2, track, fValuesHadron); - + // for the energy correlator analysis + VarManager::FillEnergyCorrelator(dilepton, track, fValuesHadron, fConfigOptions.fConfigApplyMassEC); if (!track.has_mcParticle()) { continue; } @@ -2938,29 +2944,28 @@ struct AnalysisDileptonTrack { { auto groupedMCTracks = mcTracks.sliceBy(perReducedMcEvent, event.mcCollisionId()); groupedMCTracks.bindInternalIndicesTo(&mcTracks); - for (auto& [t1, t2] : combinations(groupedMCTracks, groupedMCTracks)) { + for (auto& t1 : groupedMCTracks) { auto t1_raw = mcTracks.rawIteratorAt(t1.globalIndex()); // apply kinematic cuts for signal - if ((t1_raw.pt() < fConfigOptions.fConfigDileptonLowpTCut || t1_raw.pt() > fConfigOptions.fConfigDileptonHighpTCut)) { - continue; - } - if (abs(t1_raw.y()) > fConfigOptions.fConfigDileptonRapCutAbs) { + if ((t1_raw.pt() < fConfigOptions.fConfigDileptonLowpTCut || t1_raw.pt() > fConfigOptions.fConfigDileptonHighpTCut)) continue; - } - auto t2_raw = mcTracks.rawIteratorAt(t2.globalIndex()); - if (TMath::Abs(t2_raw.pdgCode()) == 443 || TMath::Abs(t2_raw.pdgCode()) == 11 || TMath::Abs(t2_raw.pdgCode()) == 22) { - continue; - } - if (t2_raw.pt() < fConfigMCOptions.fConfigMCGenHadronPtMin.value || std::abs(t2_raw.eta()) > fConfigMCOptions.fConfigMCGenHadronEtaAbs.value) { - continue; - } - if (t2_raw.getGenStatusCode() <= 0) { + if (abs(t1_raw.y()) > fConfigOptions.fConfigDileptonRapCutAbs) continue; - } - for (auto& sig : fGenMCSignals) { - if (sig->CheckSignal(true, t1_raw)) { - VarManager::FillEnergyCorrelatorsMC(t1_raw, t2_raw, VarManager::fgValues); - fHistMan->FillHistClass(Form("MCTruthEenergyCorrelators_%s", sig->GetName()), VarManager::fgValues); + // for the energy correlators + for (auto& t2 : groupedMCTracks) { + auto t2_raw = groupedMCTracks.rawIteratorAt(t2.globalIndex()); + if (TMath::Abs(t2_raw.pdgCode()) == 443 || TMath::Abs(t2_raw.pdgCode()) == 11 || TMath::Abs(t2_raw.pdgCode()) == 22) + continue; + if (t2_raw.pt() < fConfigMCOptions.fConfigMCGenHadronPtMin.value || std::abs(t2_raw.eta()) > fConfigMCOptions.fConfigMCGenHadronEtaAbs.value) { + continue; + } + if (t2_raw.getGenStatusCode() <= 0) + continue; + VarManager::FillEnergyCorrelatorsMC(t1_raw, t2_raw, VarManager::fgValues); + for (auto& sig : fGenMCSignals) { + if (sig->CheckSignal(true, t1_raw)) { + fHistMan->FillHistClass(Form("MCTruthEenergyCorrelators_%s", sig->GetName()), VarManager::fgValues); + } } } } @@ -3078,6 +3083,99 @@ struct AnalysisDileptonTrack { } } + void processMCEnergyCorrelatorsUnfolding(soa::Filtered const& events, BCsWithTimestamps const& bcs, + soa::Join const& assocs, + MyBarrelTracksWithCov const& tracks, soa::Filtered const& dileptons, + McCollisions const& /* mcEvents*/, McParticles const& mcTracks) + { + // set up KF or DCAfitter + if (events.size() == 0) { + return; + } + if (fCurrentRun != bcs.begin().runNumber()) { // start: runNumber + initParamsFromCCDB(bcs.begin().timestamp()); + fCurrentRun = bcs.begin().runNumber(); + } // end: runNumber + for (auto& event : events) { + if (!event.isEventSelected_bit(0)) { + continue; + } + auto groupedBarrelAssocs = assocs.sliceBy(trackAssocsPerCollision, event.globalIndex()); + auto groupedDielectrons = dileptons.sliceBy(dielectronsPerCollision, event.globalIndex()); + + uint32_t mcDecision = static_cast(0); + size_t isig = 0; + + for (auto dilepton : groupedDielectrons) { + // get full track info of tracks based on the index + auto lepton1 = tracks.rawIteratorAt(dilepton.index0Id()); + auto lepton2 = tracks.rawIteratorAt(dilepton.index1Id()); + if (!lepton1.has_mcParticle() || !lepton2.has_mcParticle()) { + continue; + } + auto lepton1MC = lepton1.mcParticle(); + auto lepton2MC = lepton2.mcParticle(); + if (!lepton1MC.has_mothers()) + continue; + const auto& motherParticle = lepton1MC.template mothers_first_as(); + // Check that the dilepton has zero charge + if (dilepton.sign() != 0) { + continue; + } + // loop over track associations + for (auto& assoc : groupedBarrelAssocs) { + uint32_t trackSelection = 0; + // check the cuts fulfilled by this candidate track; if none just continue + trackSelection = (assoc.isBarrelSelected_raw() & fTrackCutBitMap); + if (!trackSelection) { + continue; + } + // get the track from this association + auto track = tracks.rawIteratorAt(assoc.trackId()); + // check that this track is not included in the current dilepton + if (track.globalIndex() == dilepton.index0Id() || track.globalIndex() == dilepton.index1Id()) { + continue; + } + + if (!track.has_mcParticle()) { + continue; + } + auto trackMC = track.mcParticle(); + VarManager::FillEnergyCorrelatorsMCUnfolding(dilepton, track, motherParticle, trackMC, VarManager::fgValues); + + mcDecision = 0; + isig = 0; + for (auto sig = fRecMCSignals.begin(); sig != fRecMCSignals.end(); sig++, isig++) { + if ((*sig)->CheckSignal(true, lepton1MC, lepton2MC, trackMC)) { + mcDecision |= (static_cast(1) << isig); + } + } + // Fill histograms for the triplets + // loop over dilepton / ditrack cuts and MC signals + for (int icut = 0; icut < fNCuts; icut++) { + + if (!dilepton.filterMap_bit(icut)) { + continue; + } + + // loop over specified track cuts (the tracks to be combined with the dileptons) + for (int iTrackCut = 0; iTrackCut < fNCuts; iTrackCut++) { + + if (!(trackSelection & (static_cast(1) << iTrackCut))) { + continue; + } + for (uint32_t isig = 0; isig < fRecMCSignals.size(); isig++) { + if (mcDecision & (static_cast(1) << isig)) { + fHistMan->FillHistClass(Form("DileptonTrackMCUnfold_%s_%s_%s", fLegCutNames[icut].Data(), fTrackCutNames[iTrackCut].Data(), fRecMCSignals[isig]->GetName()), VarManager::fgValues); + } + } + } + } // end loop over track cuts + } // end loop over dilepton cuts + } // end loop over dileptons + } // end loop over events + } + void processDummy(MyEvents&) { // do nothing @@ -3091,6 +3189,7 @@ struct AnalysisDileptonTrack { PROCESS_SWITCH(AnalysisDileptonTrack, processMCGenEnergyCorrelatorsPion, "Loop over MC particle stack and fill generator level histograms(energy correlators)", false); PROCESS_SWITCH(AnalysisDileptonTrack, processMCGenEnergyCorrelatorsME, "Loop over MC particle stack and fill generator level histograms(energy correlators)", false); PROCESS_SWITCH(AnalysisDileptonTrack, processMCGenEnergyCorrelatorsPionME, "Loop over MC particle stack and fill generator level histograms(energy correlators)", false); + PROCESS_SWITCH(AnalysisDileptonTrack, processMCEnergyCorrelatorsUnfolding, "Loop over MC particle stack and fill response matrix(energy correlators)", false); PROCESS_SWITCH(AnalysisDileptonTrack, processDummy, "Dummy function", true); }; @@ -3205,6 +3304,8 @@ void DefineHistograms(HistogramManager* histMan, TString histClasses, const char if (classStr.Contains("MCTruthEenergyCorrelators")) { dqhistograms::DefineHistograms(histMan, objArray->At(iclass)->GetName(), "energy-correlator-gen"); } - + if (classStr.Contains("DileptonTrackMCUnfold")) { + dqhistograms::DefineHistograms(histMan, objArray->At(iclass)->GetName(), "energy-correlator-unfolding"); + } } // end loop over histogram classes }