diff --git a/PWGDQ/Core/HistogramsLibrary.cxx b/PWGDQ/Core/HistogramsLibrary.cxx index ffb02b8ad3f..4cc9f07a0bf 100644 --- a/PWGDQ/Core/HistogramsLibrary.cxx +++ b/PWGDQ/Core/HistogramsLibrary.cxx @@ -969,6 +969,50 @@ void o2::aod::dqhistograms::DefineHistograms(HistogramManager* hm, const char* h hm->AddHistogram(histClass, "MCImpPar_CentFT0CMC", "MC impact param vs MC Cent. FT0C", false, 20, 0.0, 20.0, VarManager::kMCEventImpParam, 100, 0.0, 100.0, VarManager::kMCEventCentrFT0C); } + 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, "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); + hm->AddHistogram(histClass, "Coschi_hadron_kMCWeight_hadron_DeltaEta", "", false, 40, -1.0, 1.0, VarManager::kMCCosChi, 50, -8.0, 8.0, VarManager::kMCdeltaeta, 120, 0.0, 50, VarManager::kMCWeight); + hm->AddHistogram(histClass, "Costheta_hadron_kMCWeight_before_DeltaEta", "", false, 40, -1.0, 1.0, VarManager::kMCCosTheta, 50, -8.0, 8.0, VarManager::kMCdeltaeta, 120, 0.0, 50, VarManager::kMCWeight_before); + // for bkg + hm->AddHistogram(histClass, "DeltaPhi_minus", "", false, 50, -8.0, 8.0, VarManager::kMCdeltaphi_minus); + hm->AddHistogram(histClass, "DeltaPhi_toward_minus", "", false, 50, -8.0, 8.0, VarManager::kMCdeltaphi_toward_minus); + hm->AddHistogram(histClass, "DeltaPhi_away_minus", "", false, 50, -8.0, 8.0, VarManager::kMCdeltaphi_away_minus); + hm->AddHistogram(histClass, "DeltaPhi_trans_minus", "", false, 50, -8.0, 8.0, VarManager::kMCdeltaphi_trans_minus); + + hm->AddHistogram(histClass, "Coschi_minus", "", false, 40, -1.0, 1.0, VarManager::kMCCosChi_minus, 0, 0, 0, -1, 0, 0, 0, -1, "", "", "", -1, VarManager::kMCWeight_minus); + hm->AddHistogram(histClass, "Coschi_toward_minus", "", false, 40, -1.0, 1.0, VarManager::kMCCosChi_toward_minus, 0, 0, 0, -1, 0, 0, 0, -1, "", "", "", -1, VarManager::kMCWeight_toward_minus); + hm->AddHistogram(histClass, "Coschi_away_minus", "", false, 40, -1.0, 1.0, VarManager::kMCCosChi_away_minus, 0, 0, 0, -1, 0, 0, 0, -1, "", "", "", -1, VarManager::kMCWeight_away_minus); + hm->AddHistogram(histClass, "Coschi_trans_minus", "", false, 40, -1.0, 1.0, VarManager::kMCCosChi_trans_minus, 0, 0, 0, -1, 0, 0, 0, -1, "", "", "", -1, VarManager::kMCWeight_trans_minus); + + hm->AddHistogram(histClass, "Costheta_minus", "", false, 40, -1.0, 1.0, VarManager::kMCCosTheta_minus, 0, 0, 0, -1, 0, 0, 0, -1, "", "", "", -1, VarManager::kMCWeight_before); + hm->AddHistogram(histClass, "Costheta_toward_minus", "", false, 40, -1.0, 1.0, VarManager::kMCCosTheta_toward_minus, 0, 0, 0, -1, 0, 0, 0, -1, "", "", "", -1, VarManager::kMCWeight_before); + hm->AddHistogram(histClass, "Costheta_away_minus", "", false, 40, -1.0, 1.0, VarManager::kMCCosTheta_away_minus, 0, 0, 0, -1, 0, 0, 0, -1, "", "", "", -1, VarManager::kMCWeight_before); + hm->AddHistogram(histClass, "Costheta_trans_minus", "", false, 40, -1.0, 1.0, VarManager::kMCCosTheta_trans_minus, 0, 0, 0, -1, 0, 0, 0, -1, "", "", "", -1, VarManager::kMCWeight_before); + + hm->AddHistogram(histClass, "DeltaPhi_plus", "", false, 50, -8.0, 8.0, VarManager::kMCdeltaphi_plus); + hm->AddHistogram(histClass, "DeltaPhi_toward_plus", "", false, 50, -8.0, 8.0, VarManager::kMCdeltaphi_toward_plus); + hm->AddHistogram(histClass, "DeltaPhi_away_plus", "", false, 50, -8.0, 8.0, VarManager::kMCdeltaphi_away_plus); + hm->AddHistogram(histClass, "DeltaPhi_trans_plus", "", false, 50, -8.0, 8.0, VarManager::kMCdeltaphi_trans_plus); + + hm->AddHistogram(histClass, "Coschi_plus", "", false, 40, -1.0, 1.0, VarManager::kMCCosChi_plus, 0, 0, 0, -1, 0, 0, 0, -1, "", "", "", -1, VarManager::kMCWeight_plus); + hm->AddHistogram(histClass, "Coschi_toward_plus", "", false, 40, -1.0, 1.0, VarManager::kMCCosChi_toward_plus, 0, 0, 0, -1, 0, 0, 0, -1, "", "", "", -1, VarManager::kMCWeight_toward_plus); + hm->AddHistogram(histClass, "Coschi_away_plus", "", false, 40, -1.0, 1.0, VarManager::kMCCosChi_away_plus, 0, 0, 0, -1, 0, 0, 0, -1, "", "", "", -1, VarManager::kMCWeight_away_plus); + hm->AddHistogram(histClass, "Coschi_trans_plus", "", false, 40, -1.0, 1.0, VarManager::kMCCosChi_trans_plus, 0, 0, 0, -1, 0, 0, 0, -1, "", "", "", -1, VarManager::kMCWeight_trans_plus); + + hm->AddHistogram(histClass, "Costheta_plus", "", false, 40, -1.0, 1.0, VarManager::kMCCosTheta_plus, 0, 0, 0, -1, 0, 0, 0, -1, "", "", "", -1, VarManager::kMCWeight_before); + hm->AddHistogram(histClass, "Costheta_toward_plus", "", false, 40, -1.0, 1.0, VarManager::kMCCosTheta_toward_plus, 0, 0, 0, -1, 0, 0, 0, -1, "", "", "", -1, VarManager::kMCWeight_before); + hm->AddHistogram(histClass, "Costheta_away_plus", "", false, 40, -1.0, 1.0, VarManager::kMCCosTheta_away_plus, 0, 0, 0, -1, 0, 0, 0, -1, "", "", "", -1, VarManager::kMCWeight_before); + hm->AddHistogram(histClass, "Costheta_trans_plus", "", false, 40, -1.0, 1.0, VarManager::kMCCosTheta_trans_plus, 0, 0, 0, -1, 0, 0, 0, -1, "", "", "", -1, VarManager::kMCWeight_before); + } if (!groupStr.CompareTo("pair")) { if (subGroupStr.Contains("cepf")) { hm->AddHistogram(histClass, "Mass", "", false, 300, 0.0, 12.0, VarManager::kMass); @@ -1008,6 +1052,10 @@ void o2::aod::dqhistograms::DefineHistograms(HistogramManager* hm, const char* h hm->AddHistogram(histClass, "Eta_Pt", "", false, 40, -2.0, 2.0, VarManager::kEta, 40, 0.0, 20.0, VarManager::kPt); hm->AddHistogram(histClass, "Y_Pt", "", false, 40, -2.0, 2.0, VarManager::kRap, 40, 0.0, 20.0, VarManager::kPt); hm->AddHistogram(histClass, "Mass_VtxZ", "", true, 30, -15.0, 15.0, VarManager::kVtxZ, 500, 0.0, 5.0, VarManager::kMass); + if (subGroupStr.Contains("energy-correlator")) { + hm->AddHistogram(histClass, "Mass_Y_Pt", "", false, 125, 0.0, 5.0, VarManager::kMass, 40, -2.0, 2.0, VarManager::kRap, 40, 0.0, 20.0, VarManager::kPt); + hm->AddHistogram(histClass, "Y", ";y", false, 40, -2.0, 2.0, VarManager::kRap); + } if (subGroupStr.Contains("pbpb")) { hm->AddHistogram(histClass, "Mass_CentFT0C", "", false, 125, 0.0, 5.0, VarManager::kMass, 20, 0.0, 100.0, VarManager::kCentFT0C); hm->AddHistogram(histClass, "Pt_CentFT0C", "", false, 100, 0.0, 10.0, VarManager::kPt, 20, 0.0, 100.0, VarManager::kCentFT0C); @@ -1897,6 +1945,18 @@ void o2::aod::dqhistograms::DefineHistograms(HistogramManager* hm, const char* h if (subGroupStr.Contains("opencharm")) { hm->AddHistogram(histClass, "Delta_Mass_DstarD0region", "", false, 50, 0.14, 0.16, VarManager::kDeltaMass); } + if (subGroupStr.Contains("energy-correlator")) { + hm->AddHistogram(histClass, "Coschi", "", false, 40, -1.0, 1.0, VarManager::kCosChi, 0, 0, 0, -1, 0, 0, 0, -1, "", "", "", -1, VarManager::kECWeight); + hm->AddHistogram(histClass, "CosTheta", "", false, 40, -1.0, 1.0, VarManager::kCosTheta, 0, 0, 0, -1, 0, 0, 0, -1, "", "", "", -1, VarManager::kEWeight_before); + hm->AddHistogram(histClass, "Pt_Hadron", ";P_{T}", false, 120, 0.0, 30.0, VarManager::kPtDau); + 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", "", false, 20, -2.0, 2.0, VarManager::kDeltaEta, 50, -8.0, 8.0, VarManager::kDeltaPhi); + hm->AddHistogram(histClass, "DeltaEta", "", false, 20, -2.0, 2.0, VarManager::kDeltaEta); + hm->AddHistogram(histClass, "DeltaPhi", "", false, 50, -8.0, 8.0, VarManager::kDeltaPhi); + hm->AddHistogram(histClass, "Coschi_hadron_kMCWeight_hadron_DeltaEta", "", false, 40, -1.0, 1.0, VarManager::kCosChi, 50, -8.0, 8.0, VarManager::kDeltaEta, 120, 0.0, 50, VarManager::kECWeight); + hm->AddHistogram(histClass, "Costheta_hadron_kMCWeight_before_DeltaEta", "", false, 40, -1.0, 1.0, VarManager::kCosTheta, 50, -8.0, 8.0, VarManager::kDeltaEta, 120, 0.0, 50, VarManager::kEWeight_before); + } } if (!groupStr.CompareTo("dilepton-charmhadron")) { diff --git a/PWGDQ/Core/MCSignalLibrary.cxx b/PWGDQ/Core/MCSignalLibrary.cxx index abf7b18e82d..e335a09fbae 100644 --- a/PWGDQ/Core/MCSignalLibrary.cxx +++ b/PWGDQ/Core/MCSignalLibrary.cxx @@ -1909,6 +1909,33 @@ MCSignal* o2::aod::dqmcsignals::GetMCSignal(const char* name) signal = new MCSignal(name, "Electron pair and pion pair from Psi2S", {pronge, pronge, prongPi, prongPi}, {2, 2, 1, 1}); return signal; } + + if (!nameStr.compare("eeFromPromptJpsiAnyPrimary")) { + MCProng pronge(2, {11, 443}, {true, true}, {false, false}, {0, 0}, {0, 0}, {false, false}, false, {503}, {true}); + pronge.SetSourceBit(0, MCProng::kPhysicalPrimary); + MCProng prongPrimary(1); + prongPrimary.SetSourceBit(0, MCProng::kPhysicalPrimary); + signal = new MCSignal(name, "anyprimary and electron pair from prompt jpsi", {pronge, pronge, prongPrimary}, {1, 1, -1}); + return signal; + } + + if (!nameStr.compare("eeFromJpsiAnyPrimary")) { + MCProng pronge(2, {11, 443}, {true, true}, {false, false}, {0, 0}, {0, 0}, {false, false}); + pronge.SetSourceBit(0, MCProng::kPhysicalPrimary); + MCProng prongPrimary(1); + prongPrimary.SetSourceBit(0, MCProng::kPhysicalPrimary); + signal = new MCSignal(name, "anyprimary and electron pair from prompt jpsi", {pronge, pronge, prongPrimary}, {1, 1, -1}); + return signal; + } + + if (!nameStr.compare("eeFromNonPromptJpsiAnyPrimary")) { + MCProng pronge(3, {11, 443, 503}, {true, true, true}, {false, false, false}, {0, 0, 0}, {0, 0, 0}, {false, false, false}); + pronge.SetSourceBit(0, MCProng::kPhysicalPrimary); + MCProng prongPrimary(1); + prongPrimary.SetSourceBit(0, MCProng::kPhysicalPrimary); + signal = new MCSignal(name, "anyprimary and electron pair from non-prompt jpsi", {pronge, pronge, prongPrimary}, {1, 1, -1}); + return signal; + } return nullptr; } diff --git a/PWGDQ/Core/VarManager.cxx b/PWGDQ/Core/VarManager.cxx index 30f68fdd03e..04445304d56 100644 --- a/PWGDQ/Core/VarManager.cxx +++ b/PWGDQ/Core/VarManager.cxx @@ -8,12 +8,14 @@ // In applying this license CERN does not waive the privileges and immunities // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. +#include "PWGDQ/Core/VarManager.h" + +#include "Tools/KFparticle/KFUtilities.h" + #include #include -#include #include -#include "PWGDQ/Core/VarManager.h" -#include "Tools/KFparticle/KFUtilities.h" +#include using std::cout; using std::endl; @@ -29,7 +31,7 @@ bool VarManager::fgUsedKF = false; float VarManager::fgMagField = 0.5; float VarManager::fgzMatching = -77.5; float VarManager::fgValues[VarManager::kNVars] = {0.0f}; -float VarManager::fgTPCInterSectorBoundary = 1.0; // cm +float VarManager::fgTPCInterSectorBoundary = 1.0; // cm int VarManager::fgITSROFbias = 0; int VarManager::fgITSROFlength = 100; int VarManager::fgITSROFBorderMarginLow = 0; @@ -734,6 +736,24 @@ void VarManager::SetDefaultVarNames() fgVariableUnits[kMuonTimeRes] = "ns"; fgVariableNames[kMCPdgCode] = "MC PDG code"; fgVariableUnits[kMCPdgCode] = ""; + fgVariableNames[kMCCosTheta] = "Cos#theta"; + fgVariableUnits[kMCCosTheta] = ""; + fgVariableNames[kMCHadronPdgCode] = "HadronPdgCode"; + fgVariableUnits[kMCHadronPdgCode] = ""; + fgVariableNames[kMCCosChi] = "Cos#chi"; + fgVariableUnits[kMCCosChi] = ""; + fgVariableNames[kMCJpsiPt] = "Jpsi p_{T}"; + fgVariableUnits[kMCJpsiPt] = "GeV/c"; + fgVariableNames[kMCHadronPt] = "Hadron p_{T}"; + fgVariableUnits[kMCHadronPt] = "GeV/c"; + fgVariableNames[kMCHadronEta] = "Hadron #eta"; + fgVariableUnits[kMCHadronEta] = ""; + fgVariableNames[kMCdeltaphi] = "#Delta#phi"; + fgVariableUnits[kMCdeltaphi] = ""; + fgVariableNames[kMCdeltaeta] = "#Delta#eta"; + fgVariableUnits[kMCdeltaeta] = ""; + fgVariableNames[kNhadron] = "N_{hadron}"; + fgVariableUnits[kNhadron] = ""; fgVariableNames[kMCParticleWeight] = "MC particle weight"; fgVariableUnits[kMCParticleWeight] = ""; fgVariableNames[kMCPx] = "MC px"; @@ -1162,6 +1182,16 @@ void VarManager::SetDefaultVarNames() fgVariableUnits[kDeltaPhi] = "rad."; fgVariableNames[kDeltaPhiSym] = "#Delta#phi"; fgVariableUnits[kDeltaPhiSym] = "rad."; + fgVariableNames[kCosChi] = "Cos#chi"; + fgVariableUnits[kCosChi] = ""; + fgVariableNames[kCosTheta] = "Cos#theta"; + fgVariableUnits[kCosTheta] = ""; + fgVariableNames[kPtDau] = "hadron P_{T}"; + fgVariableUnits[kPtDau] = "GeV/c"; + fgVariableNames[kEtaDau] = "hadron #eta"; + fgVariableUnits[kEtaDau] = ""; + fgVariableNames[kPhiDau] = "hadron #phi"; + fgVariableUnits[kPhiDau] = ""; fgVariableNames[kCosThetaHE] = "cos#it{#theta}"; fgVariableUnits[kCosThetaHE] = ""; fgVariableNames[kPhiHE] = "#varphi_{HE}"; @@ -1691,6 +1721,12 @@ void VarManager::SetDefaultVarNames() fgVarNamesMap["kMuonNAssocsOutOfBunch"] = kMuonNAssocsOutOfBunch; fgVarNamesMap["kNMuonTrackVariables"] = kNMuonTrackVariables; fgVarNamesMap["kMCPdgCode"] = kMCPdgCode; + fgVarNamesMap["kMCCosTheta"] = kMCCosTheta; + fgVarNamesMap["kMCHadronPdgCode"] = kMCHadronPdgCode; + fgVarNamesMap["kMCCosChi"] = kMCCosChi; + fgVarNamesMap["kMCHadronPt"] = kMCHadronPt; + fgVarNamesMap["kMCCosChi_minus"] = kMCCosChi_minus; + fgVarNamesMap["kMCWeight_before"] = kMCWeight_before; fgVarNamesMap["kMCParticleWeight"] = kMCParticleWeight; fgVarNamesMap["kMCPx"] = kMCPx; fgVarNamesMap["kMCPy"] = kMCPy; @@ -1896,6 +1932,13 @@ void VarManager::SetDefaultVarNames() fgVarNamesMap["kDeltaEta"] = kDeltaEta; fgVarNamesMap["kDeltaPhi"] = kDeltaPhi; fgVarNamesMap["kDeltaPhiSym"] = kDeltaPhiSym; + fgVarNamesMap["kCosTheta"] = kCosTheta; + fgVarNamesMap["kCosChi"] = kCosChi; + fgVarNamesMap["kECWeight"] = kECWeight; + fgVarNamesMap["kEWeight_before"] = kEWeight_before; + fgVarNamesMap["kPtDau"] = kPtDau; + fgVarNamesMap["kEtaDau"] = kEtaDau; + fgVarNamesMap["kPhiDau"] = kPhiDau; fgVarNamesMap["kNCorrelationVariables"] = kNCorrelationVariables; fgVarNamesMap["kQuadMass"] = kQuadMass; fgVarNamesMap["kQuadDefaultDileptonMass"] = kQuadDefaultDileptonMass; diff --git a/PWGDQ/Core/VarManager.h b/PWGDQ/Core/VarManager.h index 3bbde1dd281..675f2971554 100644 --- a/PWGDQ/Core/VarManager.h +++ b/PWGDQ/Core/VarManager.h @@ -156,6 +156,8 @@ class VarManager : public TObject kDecayToKPi, // e.g. D0 -> K+ pi- or cc. kTripleCandidateToKPiPi, // e.g. D+ -> K- pi+ pi+ kTripleCandidateToPKPi, // e.g. Lambda_c -> p K- pi+ + kJpsiHadronMass, // using the real hadron mass + kJpsiPionMass, // treat the hadron as pion kNMaxCandidateTypes }; @@ -617,6 +619,49 @@ class VarManager : public TObject kMCY, kMCParticleGeneratorId, kNMCParticleVariables, + kMCHadronPdgCode, + kMCCosTheta, + kMCJpsiPt, + kMCCosChi, + kMCdeltaphi, + kMCdeltaeta, + kMCHadronPt, + kMCHadronEta, + kMCWeight, + kNhadron, + kMCCosChi_minus, + kMCCosTheta_minus, + kMCWeight_minus, + kMCCosChi_toward_minus, + kMCCosTheta_toward_minus, + kMCWeight_toward_minus, + kMCCosChi_away_minus, + kMCCosTheta_away_minus, + kMCWeight_away_minus, + kMCCosChi_trans_minus, + kMCCosTheta_trans_minus, + kMCWeight_trans_minus, + kMCdeltaphi_minus, + kMCdeltaphi_toward_minus, + kMCdeltaphi_away_minus, + kMCdeltaphi_trans_minus, + kMCCosChi_plus, + kMCCosTheta_plus, + kMCWeight_plus, + kMCCosChi_toward_plus, + kMCCosTheta_toward_plus, + kMCWeight_toward_plus, + kMCCosChi_away_plus, + kMCCosTheta_away_plus, + kMCWeight_away_plus, + kMCCosChi_trans_plus, + kMCCosTheta_trans_plus, + kMCWeight_trans_plus, + kMCdeltaphi_plus, + kMCdeltaphi_toward_plus, + kMCdeltaphi_away_plus, + kMCdeltaphi_trans_plus, + kMCWeight_before, // MC mother particle variables kMCMotherPdgCode, @@ -825,6 +870,13 @@ class VarManager : public TObject kDeltaPhiSym, kNCorrelationVariables, kDileptonHadronKstar, + kCosChi, + kEtaDau, + kPhiDau, + kECWeight, + kPtDau, + kCosTheta, + kEWeight_before, // Dilepton-track-track variables kQuadMass, @@ -1100,6 +1152,8 @@ class VarManager : public TObject static void FillTrackCollisionMatCorr(T const& track, C const& collision, M const& materialCorr, P const& propagator, float* values = nullptr); template 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 FillPairPropagateMuon(T1 const& muon1, T2 const& muon2, const C& collision, float* values = nullptr); template @@ -1230,15 +1284,15 @@ class VarManager : public TObject static float fgMagField; static float fgzMatching; - static float fgCenterOfMassEnergy; // collision energy - static float fgMassofCollidingParticle; // mass of the colliding particle - static float fgTPCInterSectorBoundary; // TPC inter-sector border size at the TPC outer radius, in cm - static int fgITSROFbias; // ITS ROF bias (from ALPIDE parameters) - static int fgITSROFlength; // ITS ROF length (from ALPIDE parameters) - static int fgITSROFBorderMarginLow; // ITS ROF border low margin - static int fgITSROFBorderMarginHigh; // ITS ROF border high margin - static uint64_t fgSOR; // Timestamp for start of run - static uint64_t fgEOR; // Timestamp for end of run + static float fgCenterOfMassEnergy; // collision energy + static float fgMassofCollidingParticle; // mass of the colliding particle + static float fgTPCInterSectorBoundary; // TPC inter-sector border size at the TPC outer radius, in cm + static int fgITSROFbias; // ITS ROF bias (from ALPIDE parameters) + static int fgITSROFlength; // ITS ROF length (from ALPIDE parameters) + static int fgITSROFBorderMarginLow; // ITS ROF border low margin + static int fgITSROFBorderMarginHigh; // ITS ROF border high margin + static uint64_t fgSOR; // Timestamp for start of run + static uint64_t fgEOR; // Timestamp for end of run static ROOT::Math::PxPyPzEVector fgBeamA; // beam from A-side 4-momentum vector static ROOT::Math::PxPyPzEVector fgBeamC; // beam from C-side 4-momentum vector @@ -1255,6 +1309,8 @@ class VarManager : public TObject static float calculateCosPA(KFParticle kfp, KFParticle PV); template static float calculatePhiV(const T1& t1, const T2& t2); + template + static float LorentzTransformJpsihadroncosChi(TString Option, const T1& v1, const T2& v2); static o2::vertexing::DCAFitterN<2> fgFitterTwoProngBarrel; static o2::vertexing::DCAFitterN<3> fgFitterThreeProngBarrel; @@ -2757,6 +2813,114 @@ void VarManager::FillTrackMC(const U& mcStack, T const& track, float* values) FillTrackDerived(values); } +template +void VarManager::FillEnergyCorrelatorsMC(T const& track, T1 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(track.pt(), track.eta(), track.phi(), o2::constants::physics::MassJPsi); + float deltaphi = RecoDecay::constrainAngle(track.phi() - t1.phi(), -o2::constants::math::PIHalf); + float deltaeta = t1.eta() - track.eta(); + ROOT::Math::PtEtaPhiMVector v2(t1.pt(), t1.eta(), t1.phi(), MassHadron); + float E_boost = LorentzTransformJpsihadroncosChi("weight_boost", v1, v2); + float CosChi = LorentzTransformJpsihadroncosChi("coschi", v1, v2); + float CosTheta = LorentzTransformJpsihadroncosChi("costheta", v1, v2); + values[kMCCosChi] = CosChi; + values[kMCWeight_before] = t1.pt() / o2::constants::physics::MassJPsi; + values[kMCCosTheta] = CosTheta; + values[kMCdeltaphi] = deltaphi; + values[kMCdeltaeta] = deltaeta; + values[kMCHadronPt] = t1.pt(); + values[kMCHadronEta] = t1.eta(); + values[kMCHadronPdgCode] = t1.pdgCode(); + values[kMCWeight] = E_boost / o2::constants::physics::MassJPsi; + values[kMCCosChi_minus] = -999.9f; + values[kMCCosTheta_minus] = -999.9f; + values[kMCCosChi_toward_minus] = -999.9f; + values[kMCCosTheta_toward_minus] = -999.9f; + values[kMCCosChi_away_minus] = -999.9f; + values[kMCCosTheta_away_minus] = -999.9f; + values[kMCCosChi_trans_minus] = -999.9f; + values[kMCCosTheta_trans_minus] = -999.9f; + values[kMCdeltaphi_minus] = -999.9f; + values[kMCdeltaphi_toward_minus] = -999.9f; + values[kMCdeltaphi_away_minus] = -999.9f; + values[kMCdeltaphi_trans_minus] = -999.9f; + + if ((deltaphi > -0.5 * TMath::Pi() && deltaphi < -0.25 * TMath::Pi()) || (deltaphi > 1.25 * TMath::Pi() && deltaphi < 1.5 * TMath::Pi())) { + + values[kMCCosChi_minus] = CosChi; + values[kMCCosTheta_minus] = CosTheta; + values[kMCWeight_minus] = E_boost / o2::constants::physics::MassJPsi; + + ROOT::Math::PtEtaPhiMVector v2_toward_minus(t1.pt(), t1.eta(), t1.phi() - 0.5 * TMath::Pi(), MassHadron); + values[kMCCosChi_toward_minus] = LorentzTransformJpsihadroncosChi("coschi", v1, v2_toward_minus); + values[kMCCosTheta_toward_minus] = LorentzTransformJpsihadroncosChi("costheta", v1, v2_toward_minus); + values[kMCWeight_toward_minus] = LorentzTransformJpsihadroncosChi("weight_boost", v1, v2_toward_minus) / o2::constants::physics::MassJPsi; + + ROOT::Math::PtEtaPhiMVector v2_away_minus(t1.pt(), t1.eta(), t1.phi() + 0.5 * TMath::Pi(), MassHadron); + values[kMCCosChi_away_minus] = LorentzTransformJpsihadroncosChi("coschi", v1, v2_away_minus); + values[kMCCosTheta_away_minus] = LorentzTransformJpsihadroncosChi("costheta", v1, v2_away_minus); + values[kMCWeight_away_minus] = LorentzTransformJpsihadroncosChi("weight_boost", v1, v2_away_minus) / o2::constants::physics::MassJPsi; + + ROOT::Math::PtEtaPhiMVector v2_trans_minus(t1.pt(), t1.eta(), t1.phi() + TMath::Pi(), MassHadron); + values[kMCCosChi_trans_minus] = LorentzTransformJpsihadroncosChi("coschi", v1, v2_trans_minus); + values[kMCCosTheta_trans_minus] = LorentzTransformJpsihadroncosChi("costheta", v1, v2_trans_minus); + values[kMCWeight_trans_minus] = LorentzTransformJpsihadroncosChi("weight_boost", v1, v2_trans_minus) / o2::constants::physics::MassJPsi; + + values[kMCdeltaphi_minus] = deltaphi; + values[kMCdeltaphi_toward_minus] = RecoDecay::constrainAngle(track.phi() - (t1.phi() - 0.5 * TMath::Pi()), -o2::constants::math::PIHalf); + values[kMCdeltaphi_away_minus] = RecoDecay::constrainAngle(track.phi() - (t1.phi() + 0.5 * TMath::Pi()), -o2::constants::math::PIHalf); + values[kMCdeltaphi_trans_minus] = RecoDecay::constrainAngle(track.phi() - t1.phi() + TMath::Pi(), -o2::constants::math::PIHalf); + } + + values[kMCCosChi_plus] = -999.9f; + values[kMCCosTheta_plus] = -999.9f; + values[kMCCosChi_toward_plus] = -999.9f; + values[kMCCosTheta_toward_plus] = -999.9f; + values[kMCCosChi_away_plus] = -999.9f; + values[kMCCosTheta_away_plus] = -999.9f; + values[kMCCosChi_trans_plus] = -999.9f; + values[kMCCosTheta_trans_plus] = -999.9f; + values[kMCdeltaphi_plus] = -999.9f; + values[kMCdeltaphi_toward_plus] = -999.9f; + values[kMCdeltaphi_away_plus] = -999.9f; + values[kMCdeltaphi_trans_plus] = -999.9f; + + if (deltaphi > 0.25 * TMath::Pi() && deltaphi < 0.75 * TMath::Pi()) { + values[kMCCosChi_plus] = CosChi; + values[kMCCosTheta_plus] = CosTheta; + values[kMCWeight_plus] = E_boost / o2::constants::physics::MassJPsi; + + ROOT::Math::PtEtaPhiMVector v2_toward_plus(t1.pt(), t1.eta(), t1.phi() + 0.5 * TMath::Pi(), MassHadron); + values[kMCCosChi_toward_plus] = LorentzTransformJpsihadroncosChi("coschi", v1, v2_toward_plus); + values[kMCCosTheta_toward_plus] = LorentzTransformJpsihadroncosChi("costheta", v1, v2_toward_plus); + values[kMCWeight_toward_plus] = LorentzTransformJpsihadroncosChi("weight_boost", v1, v2_toward_plus) / o2::constants::physics::MassJPsi; + + ROOT::Math::PtEtaPhiMVector v2_away_plus(t1.pt(), t1.eta(), t1.phi() - 0.5 * TMath::Pi(), MassHadron); + values[kMCCosChi_away_plus] = LorentzTransformJpsihadroncosChi("coschi", v1, v2_away_plus); + values[kMCCosTheta_away_plus] = LorentzTransformJpsihadroncosChi("costheta", v1, v2_away_plus); + values[kMCWeight_away_plus] = LorentzTransformJpsihadroncosChi("weight_boost", v1, v2_away_plus) / o2::constants::physics::MassJPsi; + + ROOT::Math::PtEtaPhiMVector v2_trans_plus(t1.pt(), t1.eta(), t1.phi() + TMath::Pi(), MassHadron); + values[kMCCosChi_trans_plus] = LorentzTransformJpsihadroncosChi("coschi", v1, v2_trans_plus); + values[kMCCosTheta_trans_plus] = LorentzTransformJpsihadroncosChi("costheta", v1, v2_trans_plus); + values[kMCWeight_trans_plus] = LorentzTransformJpsihadroncosChi("weight_boost", v1, v2_trans_plus) / o2::constants::physics::MassJPsi; + + values[kMCdeltaphi_plus] = deltaphi; + values[kMCdeltaphi_toward_plus] = RecoDecay::constrainAngle(track.phi() - (t1.phi() + 0.5 * TMath::Pi()), -o2::constants::math::PIHalf); + values[kMCdeltaphi_away_plus] = RecoDecay::constrainAngle(track.phi() - (t1.phi() - 0.5 * TMath::Pi()), -o2::constants::math::PIHalf); + values[kMCdeltaphi_trans_plus] = RecoDecay::constrainAngle(track.phi() - (t1.phi() - 0.5 * TMath::Pi()), -o2::constants::math::PIHalf); + } +} + template void VarManager::FillPairPropagateMuon(T1 const& muon1, T2 const& muon2, const C& collision, float* values) { @@ -5085,6 +5249,18 @@ void VarManager::FillDileptonHadron(T1 const& dilepton, T2 const& hadron, float* double Q1 = (dilepton.mass() * dilepton.mass() - hadronMass * hadronMass) / Pinv; values[kDileptonHadronKstar] = sqrt(Q1 * Q1 - v12_Qvect.M2()) / 2.0; } + if (fgUsedVars[kCosChi] || fgUsedVars[kECWeight] || fgUsedVars[kCosTheta] || fgUsedVars[kEWeight_before] || fgUsedVars[kPtDau] || fgUsedVars[kEtaDau] || fgUsedVars[kPhiDau]) { + ROOT::Math::PtEtaPhiMVector v1(dilepton.pt(), dilepton.eta(), dilepton.phi(), dilepton.mass()); + ROOT::Math::PtEtaPhiMVector v2(hadron.pt(), hadron.eta(), hadron.phi(), o2::constants::physics::MassPionCharged); + values[kCosChi] = LorentzTransformJpsihadroncosChi("coschi", v1, v2); + float E_boost = LorentzTransformJpsihadroncosChi("weight_boost", v1, v2); + values[kECWeight] = E_boost / v1.M(); + values[kCosTheta] = LorentzTransformJpsihadroncosChi("costheta", v1, v2); + values[kEWeight_before] = v2.Pt() / v1.M(); + values[kPtDau] = v2.pt(); + values[kEtaDau] = v2.eta(); + values[kPhiDau] = v2.phi(); + } if (fgUsedVars[kDeltaPhi]) { double delta = dilepton.phi() - hadron.phi(); if (delta > 3.0 / 2.0 * M_PI) { @@ -5668,5 +5844,28 @@ void VarManager::FillBdtScore(T1 const& bdtScore, float* values) LOG(warning) << "Unexpected number of BDT outputs: " << bdtScore.size(); } } +//__________________________________________________________________ +template +float VarManager::LorentzTransformJpsihadroncosChi(TString Option, T1 const& v1, T2 const& v2) +{ + float value = -999.0f; + auto beta_v1 = v1.BoostToCM(); + ROOT::Math::Boost boostv1{beta_v1}; + auto v2_boost = boostv1(v2); + auto p_v1_lab = v1.Vect(); + float p1_lab = p_v1_lab.R(); + if (Option == "coschi") { + auto p_v2_boost = v2_boost.Vect(); + float p_boost = p_v2_boost.R(); + value = p_v2_boost.Dot(p_v1_lab) / (p1_lab * p_boost); + } else if (Option == "costheta") { + auto p_v2_lab = v2.Vect(); + float p2_lab = p_v2_lab.R(); + value = p_v1_lab.Dot(p_v2_lab) / (p1_lab * p2_lab); + } else if (Option == "weight_boost") { + value = v2_boost.E(); + } + return value; +} #endif // PWGDQ_CORE_VARMANAGER_H_ diff --git a/PWGDQ/Tasks/dqEfficiency_withAssoc.cxx b/PWGDQ/Tasks/dqEfficiency_withAssoc.cxx index 5e01a01cda0..470d32513bd 100644 --- a/PWGDQ/Tasks/dqEfficiency_withAssoc.cxx +++ b/PWGDQ/Tasks/dqEfficiency_withAssoc.cxx @@ -68,6 +68,7 @@ namespace o2::aod { namespace dqanalysisflags { +DECLARE_SOA_COLUMN(MixingHash, mixingHash, int); //! Hash used in event mixing //need to understand DECLARE_SOA_BITMAP_COLUMN(IsEventSelected, isEventSelected, 32); //! Event decision DECLARE_SOA_BITMAP_COLUMN(IsBarrelSelected, isBarrelSelected, 32); //! Barrel track decisions DECLARE_SOA_COLUMN(BarrelAmbiguityInBunch, barrelAmbiguityInBunch, int8_t); //! Barrel track in-bunch ambiguity @@ -163,6 +164,7 @@ DECLARE_SOA_COLUMN(OniaVtxZ, oniaVtxZ, float); } // namespace dqanalysisflags DECLARE_SOA_TABLE(EventCuts, "AOD", "DQANAEVCUTS", dqanalysisflags::IsEventSelected); //! joinable to ReducedEvents +DECLARE_SOA_TABLE(MixingHashes, "AOD", "DQANAMIXHASHA", dqanalysisflags::MixingHash); //! joinable to ReducedEvents DECLARE_SOA_TABLE(BarrelTrackCuts, "AOD", "DQANATRKCUTS", dqanalysisflags::IsBarrelSelected); //! joinable to ReducedTracksAssoc DECLARE_SOA_TABLE(BarrelAmbiguities, "AOD", "DQBARRELAMB", dqanalysisflags::BarrelAmbiguityInBunch, dqanalysisflags::BarrelAmbiguityOutOfBunch); //! joinable to ReducedBarrelTracks DECLARE_SOA_TABLE(MuonTrackCuts, "AOD", "DQANAMUONCUTS", dqanalysisflags::IsMuonSelected); //! joinable to ReducedMuonsAssoc @@ -201,6 +203,7 @@ using MyEventsVtxCovSelected = soa::Join; using MyEventsVtxCovSelectedQvector = soa::Join; using MyEventsQvector = soa::Join; +using MyEventsVtxCovHashSelected = soa::Join; using MyBarrelTracks = soa::Join; using MyBarrelTracksWithAmbiguities = soa::Join; @@ -244,8 +247,9 @@ void PrintBitMap(TMap map, int nbits) // Analysis task that produces event decisions and the Hash table used in event mixing struct AnalysisEventSelection { Produces eventSel; + Produces hash; OutputObj fOutputList{"output"}; - + Configurable fConfigMixingVariables{"cfgMixingVars", "", "Mixing configs separated by a comma, default no mixing"}; Configurable fConfigEventCuts{"cfgEventCuts", "eventStandard", "Event selection"}; Configurable fConfigEventCutsJSON{"cfgEventCutsJSON", "", "Additional event cuts specified in JSON format"}; Configurable fConfigQA{"cfgQA", false, "If true, fill QA histograms"}; @@ -261,6 +265,8 @@ struct AnalysisEventSelection { Configurable fConfigNoLaterThan{"ccdb-no-later-than", std::chrono::duration_cast(std::chrono::system_clock::now().time_since_epoch()).count(), "latest acceptable timestamp of creation for the object"}; HistogramManager* fHistMan = nullptr; + MixingHandler* fMixHandler = nullptr; + AnalysisCompositeCut* fEventCut; Service fCCDB; @@ -311,6 +317,16 @@ struct AnalysisEventSelection { fOutputList.setObject(fHistMan->GetMainHistogramList()); } + TString mixVarsString = fConfigMixingVariables.value; + std::unique_ptr objArray(mixVarsString.Tokenize(",")); + if (objArray->GetEntries() > 0) { + fMixHandler = new MixingHandler("mixingHandler", "mixing handler"); + fMixHandler->Init(); + for (int iVar = 0; iVar < objArray->GetEntries(); ++iVar) { + dqmixing::SetUpMixing(fMixHandler, objArray->At(iVar)->GetName()); + } + } + fCurrentRun = -1; fCCDB->setURL(fConfigCcdbUrl.value); fCCDB->setCaching(true); @@ -359,6 +375,10 @@ struct AnalysisEventSelection { auto& evIndices = fBCCollMap[event.globalBC()]; evIndices.push_back(event.globalIndex()); } + if (fMixHandler != nullptr) { + int hh = fMixHandler->FindEventCategory(VarManager::fgValues); + hash(hh); + } } for (auto& event : mcEvents) { @@ -1733,10 +1753,10 @@ struct AnalysisSameEventPairing { template void runSameEventPairing(TEvents const& events, Preslice& preslice, TTrackAssocs const& assocs, TTracks const& /*tracks*/, ReducedMCEvents const& /*mcEvents*/, ReducedMCTracks const& /*mcTracks*/) { - // if (events.size() == 0) { - // LOG(warning) << "No events in this TF, going to the next one ..."; - // return; - // } + if (events.size() == 0) { + LOG(warning) << "No events in this TF, going to the next one ..."; + return; + } if (fCurrentRun != events.begin().runNumber()) { initParamsFromCCDB(events.begin().timestamp(), TTwoProngFitter); fCurrentRun = events.begin().runNumber(); @@ -3529,7 +3549,9 @@ struct AnalysisDileptonTrack { Configurable fConfigDileptonpTCut{"cfgDileptonpTCut", 0.0, "pT cut for dileptons used in the triplet vertexing"}; Configurable fConfigDileptonLxyCut{"cfgDileptonLxyCut", 0.0, "Lxy cut for dileptons used in the triplet vertexing"}; Configurable fConfigUseKFVertexing{"cfgUseKFVertexing", false, "Use KF Particle for secondary vertex reconstruction (DCAFitter is used by default)"}; - + Configurable fConfigDileptonLowpTCut{"cfgDileptonLowpTCut", 0.0, "Low pT cut for dileptons used in the triplet vertexing"}; + Configurable fConfigDileptonHighpTCut{"cfgDileptonHighpTCut", 1E5, "High pT cut for dileptons used in the triplet vertexing"}; + Configurable fConfigDileptonRapCutAbs{"cfgDileptonRapCutAbs", 1.0, "Rap cut for dileptons used in the triplet vertexing"}; Configurable fConfigHistogramSubgroups{"cfgDileptonTrackHistogramsSubgroups", "invmass,vertexing", "Comma separated list of dilepton-track histogram subgroups"}; Configurable fConfigAddJSONHistograms{"cfgAddJSONHistograms", "", "Histograms in JSON format"}; @@ -3552,6 +3574,8 @@ struct AnalysisDileptonTrack { Configurable fConfigMCGenHadronPtMin{"cfgMCGenHadronPtMin", 1.0f, "minimum pt for the hadron"}; Configurable fConfigMCGenDileptonLegEtaAbs{"cfgMCGenDileptonLegEtaAbs", 0.9f, "eta abs range for the dilepton leg"}; Configurable fConfigMCGenHadronEtaAbs{"cfgMCGenHadronEtaAbs", 0.9f, "eta abs range for the hadron"}; + Configurable fConfigUseMCRapcut{"cfgUseMCRapcut", false, "Use Rap cut for dileptons used in the triplet vertexing(reconstructed)"}; + Configurable fConfigMixingDepth{"cfgMixingDepth", 5, "Event mixing pool depth"}; int fCurrentRun; // needed to detect if the run changed and trigger update of calibrations etc. int fNCuts; @@ -3571,7 +3595,7 @@ struct AnalysisDileptonTrack { // TODO: The filter expressions seem to always use the default value of configurables, not the values from the actual configuration file Filter eventFilter = aod::dqanalysisflags::isEventSelected > static_cast(0); - Filter dileptonFilter = aod::reducedpair::pt > fConfigDileptonpTCut&& aod::reducedpair::mass > fConfigDileptonLowMass&& aod::reducedpair::mass fConfigDileptonLxyCut; + Filter dileptonFilter = aod::reducedpair::pt > fConfigDileptonLowpTCut&& aod::reducedpair::pt fConfigDileptonLowMass&& aod::reducedpair::mass fConfigDileptonLxyCut; Filter filterBarrel = aod::dqanalysisflags::isBarrelSelected > static_cast(0); Filter filterMuon = aod::dqanalysisflags::isMuonSelected > static_cast(0); @@ -3585,6 +3609,8 @@ struct AnalysisDileptonTrack { std::vector fRecMCSignals; std::vector fGenMCSignals; + NoBinningPolicy fHashBin; + void init(o2::framework::InitContext& context) { bool isBarrel = context.mOptions.get("processBarrelSkimmed"); @@ -3592,6 +3618,8 @@ struct AnalysisDileptonTrack { bool isMuon = context.mOptions.get("processMuonSkimmed"); bool isMCGen = context.mOptions.get("processMCGen") || context.mOptions.get("processMCGenWithEventSelection"); 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"); if (isDummy) { if (isBarrel || isMuon || isBarrelAsymmetric || isMCGen) { @@ -3879,6 +3907,18 @@ struct AnalysisDileptonTrack { } } + if (isMCGen_energycorrelators) { + for (auto& sig : fGenMCSignals) { + DefineHistograms(fHistMan, Form("MCTruthEenergyCorrelators_%s", sig->GetName()), ""); + } + } + + if (isMCGen_energycorrelatorsME) { + for (auto& sig : fGenMCSignals) { + DefineHistograms(fHistMan, Form("MCTruthEenergyCorrelatorsME_%s", sig->GetName()), ""); + } + } + TString addHistsStr = fConfigAddJSONHistograms.value; if (addHistsStr != "") { dqhistograms::AddHistogramsFromJSON(fHistMan, addHistsStr.Data()); @@ -3939,6 +3979,10 @@ struct AnalysisDileptonTrack { if (dilepton.sign() != 0) { continue; } + // dilepton rap cut + float rap = dilepton.rap(); + if (fConfigUseMCRapcut && abs(rap) > fConfigDileptonRapCutAbs) + continue; VarManager::FillTrack(dilepton, fValuesDilepton); @@ -4328,6 +4372,143 @@ struct AnalysisDileptonTrack { } // end loop over reconstructed events } + template + void runEnergyCorrelators(TEvent const& event, TMCTracks const& mcTracks) + { + auto groupedMCTracks = mcTracks.sliceBy(perReducedMcEvent, event.reducedMCeventId()); + groupedMCTracks.bindInternalIndicesTo(&mcTracks); + for (auto& t1 : groupedMCTracks) { + auto t1_raw = mcTracks.rawIteratorAt(t1.globalIndex()); + // apply kinematic cuts for signal + if ((t1_raw.pt() < fConfigDileptonLowpTCut || t1_raw.pt() > fConfigDileptonHighpTCut)) + continue; + if (abs(t1_raw.y()) > fConfigDileptonRapCutAbs) + continue; + // 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() < fConfigMCGenHadronPtMin.value || std::abs(t2_raw.eta()) > 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); + } + } + } + } + } + void processMCGenEnergyCorrelators(soa::Filtered const& events, + ReducedMCEvents const& /*mcEvents*/, ReducedMCTracks const& mcTracks) + { + if (events.size() == 0) { + LOG(warning) << "No events in this TF, going to the next one ..."; + return; + } + for (auto& event : events) { + if (!event.isEventSelected_bit(0)) { + continue; + } + if (!event.has_reducedMCevent()) { + continue; + } + runEnergyCorrelators(event, mcTracks); + } + } + + void processMCGenEnergyCorrelatorsPion(soa::Filtered const& events, + ReducedMCEvents const& /*mcEvents*/, ReducedMCTracks const& mcTracks) + { + if (events.size() == 0) { + LOG(warning) << "No events in this TF, going to the next one ..."; + return; + } + for (auto& event : events) { + if (!event.isEventSelected_bit(0)) { + continue; + } + if (!event.has_reducedMCevent()) { + continue; + } + runEnergyCorrelators(event, mcTracks); + } + } + + template + void runEnergyCorrelatorsMixedEvent(TEvent const& event1, TEvent const& event2, TMCTracks const& mcTracks) + { + auto groupedMCTracks1 = mcTracks.sliceBy(perReducedMcEvent, event1.reducedMCeventId()); + auto groupedMCTracks2 = mcTracks.sliceBy(perReducedMcEvent, event2.reducedMCeventId()); + groupedMCTracks1.bindInternalIndicesTo(&mcTracks); + groupedMCTracks2.bindInternalIndicesTo(&mcTracks); + for (auto& t1 : groupedMCTracks1) { + auto t1_raw = mcTracks.rawIteratorAt(t1.globalIndex()); + // apply kinematic cuts for signal + if ((t1_raw.pt() < fConfigDileptonLowpTCut || t1_raw.pt() > fConfigDileptonHighpTCut)) + continue; + if (abs(t1_raw.y()) > fConfigDileptonRapCutAbs) + continue; + // for the energy correlators + for (auto& t2 : groupedMCTracks2) { + auto t2_raw = groupedMCTracks2.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() < fConfigMCGenHadronPtMin.value || std::abs(t2_raw.eta()) > 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("MCTruthEenergyCorrelatorsME_%s", sig->GetName()), VarManager::fgValues); + } + } + } + } + } + + void processMCGenEnergyCorrelatorsME(soa::Filtered const& events, + ReducedMCEvents const& /*mcEvents*/, ReducedMCTracks const& mcTracks) + { + if (events.size() == 0) { + LOG(warning) << "No events in this TF, going to the next one ..."; + return; + } + // loop over two event comibnations + for (auto& [event1, event2] : selfCombinations(fHashBin, fConfigMixingDepth.value, -1, events, events)) { + if (!event1.isEventSelected_bit(0) || !event2.isEventSelected_bit(0)) { + continue; + } + if (!event1.has_reducedMCevent() || !event2.has_reducedMCevent()) { + continue; + } + runEnergyCorrelatorsMixedEvent(event1, event2, mcTracks); + } + } + + void processMCGenEnergyCorrelatorsPionME(soa::Filtered const& events, + ReducedMCEvents const& /*mcEvents*/, ReducedMCTracks const& mcTracks) + { + if (events.size() == 0) { + LOG(warning) << "No events in this TF, going to the next one ..."; + return; + } + // loop over two event comibnations + for (auto& [event1, event2] : selfCombinations(fHashBin, fConfigMixingDepth.value, -1, events, events)) { + if (!event1.isEventSelected_bit(0) || !event2.isEventSelected_bit(0)) { + continue; + } + if (!event1.has_reducedMCevent() || !event2.has_reducedMCevent()) { + continue; + } + runEnergyCorrelatorsMixedEvent(event1, event2, mcTracks); + } + } + void processDummy(MyEvents&) { // do nothing @@ -4338,6 +4519,10 @@ struct AnalysisDileptonTrack { PROCESS_SWITCH(AnalysisDileptonTrack, processMuonSkimmed, "Run muon dilepton-track pairing, using skimmed data", false); PROCESS_SWITCH(AnalysisDileptonTrack, processMCGen, "Loop over MC particle stack and fill generator level histograms", false); PROCESS_SWITCH(AnalysisDileptonTrack, processMCGenWithEventSelection, "Loop over MC particle stack and fill generator level histograms", false); + PROCESS_SWITCH(AnalysisDileptonTrack, processMCGenEnergyCorrelators, "Loop over MC particle stack and fill generator level histograms(energy correlators)", false); + 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, processDummy, "Dummy function", true); }; @@ -4445,5 +4630,10 @@ void DefineHistograms(HistogramManager* histMan, TString histClasses, const char if (classStr.Contains("DileptonHadronCorrelation")) { dqhistograms::DefineHistograms(histMan, objArray->At(iclass)->GetName(), "dilepton-hadron-correlation"); } + + if (classStr.Contains("MCTruthEenergyCorrelators")) { + dqhistograms::DefineHistograms(histMan, objArray->At(iclass)->GetName(), "energy-correlator-gen"); + } + } // end loop over histogram classes }