Skip to content

Commit cc910ed

Browse files
YazhenLinCape
YazhenLin
and
Cape
authored
[PWGDQ] Add and change the code for the energy correlator study (#16199)
Co-authored-by: Cape <cape@lindeMacBook-Pro.local>
1 parent c400d8f commit cc910ed

1 file changed

Lines changed: 67 additions & 17 deletions

File tree

‎PWGDQ/Tasks/dqEnergyCorrelator_direct.cxx‎

Lines changed: 67 additions & 17 deletions
Original file line numberDiff line numberDiff line change
@@ -130,11 +130,12 @@ struct AnalysisEnergyCorrelator {
130130
Configurable<std::string> fConfigAddDileptonHadronHistogram{"cfgAddDileptonHadronHistogram", "", "Dilepton-hadron histograms"};
131131
Configurable<std::string> fConfigMCRecSignals{"cfgMCRecDileptonHadronSignals", "", "Comma separated list of MC signals (reconstructed)"};
132132
Configurable<std::string> fConfigMCGenSignals{"cfgMCGenDileptonHadronSignals", "", "Comma separated list of MC signals (generated)"};
133+
Configurable<std::string> fConfigMCGenPairSignals{"cfgMCGenDileptonHadronPairSignals", "", "Comma separated list of MC pair signals (generated)"};
133134
Configurable<std::string> fConfigMCRecSignalsJSON{"cfgMCRecDileptonHadronSignalsJSON", "", "Additional list of MC signals (reconstructed) via JSON"};
134135
Configurable<std::string> fConfigMCGenSignalsJSON{"cfgMCGenDileptonHadronSignalsJSON", "", "Comma separated list of MC signals (generated) via JSON"};
136+
Configurable<std::string> fConfigMCGenPairSignalsJSON{"cfgMCGenDileptonHadronPairSignalsJSON", "", "Comma separated list of MC pair signals (generated) via JSON"};
135137
Configurable<float> fConfigMCGenHadronEtaAbs{"cfgMCGenHadronEtaAbs", 0.9f, "eta abs range for the hadron"};
136138
Configurable<float> fConfigMCGenHadronPtMin{"cfgMCGenHadronPtMin", 0.1f, "minimum pt for the hadron"};
137-
Configurable<bool> fConfigContainlepton{"cfgContainlepton", false, "If true, require the hadron to contain the lepton in its decay tree for the energy correlator study"};
138139
Configurable<bool> fConfigUsePionMass{"cfgUsePionMass", false, "If true, use pion mass for the hadron in the energy correlator study"};
139140
Configurable<bool> fConfigApplyEfficiency{"cfgApplyEfficiency", false, "If true, apply efficiency correction for the energy correlator study"};
140141
Configurable<bool> fConfigApplyEfficiencyME{"cfgApplyEfficiencyME", false, "If true, apply efficiency correction for the energy correlator study"};
@@ -167,6 +168,7 @@ struct AnalysisEnergyCorrelator {
167168
std::vector<MCSignal*> fRecMCTrackSignals; // MC signals for reconstructed tracks
168169
std::vector<MCSignal*> fRecMCSignals; // MC signals for reconstructed pairs
169170
std::vector<MCSignal*> fGenMCSignals;
171+
std::vector<MCSignal*> fGenMCPairSignals;
170172
std::vector<MCSignal*> fRecMCTripleSignals; // MC signals for reconstructed triples
171173

172174
Service<o2::ccdb::BasicCCDBManager> fCCDB;
@@ -183,6 +185,7 @@ struct AnalysisEnergyCorrelator {
183185
TH2F* hEfficiency_dilepton;
184186
TH2F* hEfficiency_hadron;
185187
TH1F* hMasswindow;
188+
TH2F* hReweighthadron;
186189

187190
void init(o2::framework::InitContext& context)
188191
{
@@ -337,6 +340,29 @@ struct AnalysisEnergyCorrelator {
337340
}
338341
}
339342

343+
// Add histogram classes for each specified MCsignal at the generator level
344+
// TODO: create a std::vector of hist classes to be used at Fill time, to avoid using Form in the process function
345+
TString sigpairGenNamesStr = fConfigDileptonHadronOptions.fConfigMCGenPairSignals.value;
346+
std::unique_ptr<TObjArray> objGenPairSigArray(sigpairGenNamesStr.Tokenize(","));
347+
for (int isig = 0; isig < objGenPairSigArray->GetEntries(); isig++) {
348+
MCSignal* sig = o2::aod::dqmcsignals::GetMCSignal(objGenPairSigArray->At(isig)->GetName());
349+
if (sig) {
350+
fGenMCPairSignals.push_back(sig);
351+
}
352+
}
353+
354+
// Add the MCSignals from the JSON config
355+
TString addMCSignalsPairGenStr = fConfigDileptonHadronOptions.fConfigMCGenPairSignalsJSON.value;
356+
if (addMCSignalsPairGenStr != "") {
357+
std::vector<MCSignal*> addMCSignals = dqmcsignals::GetMCSignalsFromJSON(addMCSignalsPairGenStr.Data());
358+
for (auto& mcIt : addMCSignals) {
359+
if (mcIt->GetNProngs() != 2) { // NOTE: only 2 prong signals
360+
continue;
361+
}
362+
fGenMCPairSignals.push_back(mcIt);
363+
}
364+
}
365+
340366
// TODO: create a std::vector of hist classes to be used at Fill time, to avoid using Form in the process function
341367
TString sigRecTrackNamesStr = fConfigTrackOptions.fConfigMCRecTrackSignals.value;
342368
std::unique_ptr<TObjArray> objRecTrackSigArray(sigRecTrackNamesStr.Tokenize(","));
@@ -440,6 +466,14 @@ struct AnalysisEnergyCorrelator {
440466
if (sig->GetNProngs() == 1) {
441467
if (isMCGen_energycorrelators) {
442468
DefineHistograms(fHistMan, Form("MCTruthGenSel_%s", sig->GetName()), "");
469+
}
470+
}
471+
}
472+
473+
for (auto& sig : fGenMCPairSignals) {
474+
LOG(info) << "Defining histograms for pair signal: " << sig->GetNProngs();
475+
if (sig->GetNProngs() == 2) {
476+
if (isMCGen_energycorrelators) {
443477
DefineHistograms(fHistMan, Form("MCTruthEenergyCorrelators_%s", sig->GetName()), "");
444478
DefineHistograms(fHistMan, Form("MCTruthEenergyCorrelators_Pion_%s", sig->GetName()), "");
445479
}
@@ -472,7 +506,8 @@ struct AnalysisEnergyCorrelator {
472506
hAcceptance_rec = static_cast<TH2F*>(listAccs->FindObject("hAcceptance_rec"));
473507
hAcceptance_gen = static_cast<TH2F*>(listAccs->FindObject("hAcceptance_gen"));
474508
hMasswindow = static_cast<TH1F*>(listAccs->FindObject("hMasswindow"));
475-
if (!hAcceptance_rec || !hAcceptance_gen || !hEfficiency_dilepton || !hEfficiency_hadron || !hMasswindow) {
509+
hReweighthadron = static_cast<TH2F*>(listAccs->FindObject("hReweighthadron"));
510+
if (!hAcceptance_rec || !hAcceptance_gen || !hEfficiency_dilepton || !hEfficiency_hadron || !hMasswindow || !hReweighthadron) {
476511
LOG(fatal) << "Problem getting histograms from the TList object with efficiencies!";
477512
}
478513
}
@@ -533,11 +568,22 @@ struct AnalysisEnergyCorrelator {
533568
float Effdilepton = GetSafeInterpolationWeight(hEfficiency_dilepton, dilepton_rap, dilepton_pt);
534569
float Effhadron = GetSafeInterpolationWeight(hEfficiency_hadron, hadron_eta, hadron.pt());
535570
float Masswindow = hMasswindow->Interpolate(dilepton_pt);
536-
Accweight_gen = Accweight_gen * Effdilepton * Effhadron;
571+
float Reweighthadron = 1.0f;
572+
if (std::abs(hadronMC.pdgCode()) == PDG_t::kPiPlus) {
573+
int bin = hReweighthadron->FindBin(0, hadron.pt());
574+
Reweighthadron = hReweighthadron->GetBinContent(bin);
575+
} else if (std::abs(hadronMC.pdgCode()) == PDG_t::kProton) {
576+
int bin = hReweighthadron->FindBin(1, hadron.pt());
577+
Reweighthadron = hReweighthadron->GetBinContent(bin);
578+
} else if (std::abs(hadronMC.pdgCode()) == PDG_t::kKPlus) {
579+
int bin = hReweighthadron->FindBin(2, hadron.pt());
580+
Reweighthadron = hReweighthadron->GetBinContent(bin);
581+
}
582+
Accweight_gen = Accweight_gen * Effdilepton * Effhadron * Reweighthadron;
537583
if (fConfigDileptonHadronOptions.fConfigApplyEfficiencyME) {
538-
Effweight_rec = Effdilepton * Effhadron * Masswindow; // for the moment, apply the efficiency correction also for the mixed event pairs, but this can be changed in case we want to apply it only for the same event pairs
584+
Effweight_rec = Effdilepton * Effhadron * Masswindow * Reweighthadron; // for the moment, apply the efficiency correction also for the mixed event pairs, but this can be changed in case we want to apply it only for the same event pairs
539585
} else {
540-
Effweight_rec = Effweight_rec * Effdilepton * Effhadron * Masswindow; // apply acceptance and efficiency correction for the real pairs
586+
Effweight_rec = Effweight_rec * Effdilepton * Effhadron * Masswindow * Reweighthadron; // apply acceptance and efficiency correction for the real pairs
541587
}
542588
}
543589

@@ -936,25 +982,16 @@ struct AnalysisEnergyCorrelator {
936982
continue;
937983
// for the energy correlators
938984
for (auto& t2 : groupedMCTracks2) {
939-
auto t2_raw = groupedMCTracks2.rawIteratorAt(t2.globalIndex());
985+
auto t2_raw = mcTracks.rawIteratorAt(t2.globalIndex());
940986
if (t2.mcCollisionId() != event2.mcCollisionId()) { // check that the mc track belongs to the same mc collision as the reconstructed event
941987
continue;
942988
}
943-
if (!t2_raw.isPhysicalPrimary()) {
944-
continue;
945-
}
946989
if (t2_raw.has_mothers()) {
947990
auto mother_raw = t2_raw.template mothers_first_as<McParticles>();
948991
if (mother_raw.globalIndex() == t1_raw.globalIndex()) {
949992
continue;
950993
}
951994
}
952-
if (fConfigDileptonHadronOptions.fConfigContainlepton && std::abs(t2_raw.pdgCode()) != PDG_t::kPiPlus && std::abs(t2_raw.pdgCode()) != PDG_t::kKPlus && std::abs(t2_raw.pdgCode()) != PDG_t::kProton && std::abs(t2_raw.pdgCode()) != PDG_t::kElectron && std::abs(t2_raw.pdgCode()) != PDG_t::kMuonMinus) {
953-
continue;
954-
}
955-
if (!fConfigDileptonHadronOptions.fConfigContainlepton && std::abs(t2_raw.pdgCode()) != PDG_t::kPiPlus && std::abs(t2_raw.pdgCode()) != PDG_t::kKPlus && std::abs(t2_raw.pdgCode()) != PDG_t::kProton) {
956-
continue;
957-
}
958995
if (t2_raw.pt() < fConfigDileptonHadronOptions.fConfigMCGenHadronPtMin.value || std::abs(t2_raw.eta()) > fConfigDileptonHadronOptions.fConfigMCGenHadronEtaAbs.value) {
959996
continue;
960997
}
@@ -966,11 +1003,24 @@ struct AnalysisEnergyCorrelator {
9661003
float hadron_phi = t2_raw.phi();
9671004
float deltaphi = RecoDecay::constrainAngle(dilepton_phi - hadron_phi, -0.5 * o2::constants::math::PI);
9681005
acceptance = hAcceptance_gen->Interpolate(dilepton_eta - hadron_eta, deltaphi);
1006+
float Reweighthadron = 1.0f;
1007+
if (std::abs(t2_raw.pdgCode()) == PDG_t::kPiPlus) {
1008+
int bin = hReweighthadron->FindBin(0, t2_raw.pt());
1009+
Reweighthadron = hReweighthadron->GetBinContent(bin);
1010+
} else if (std::abs(t2_raw.pdgCode()) == PDG_t::kProton) {
1011+
int bin = hReweighthadron->FindBin(1, t2_raw.pt());
1012+
Reweighthadron = hReweighthadron->GetBinContent(bin);
1013+
} else if (std::abs(t2_raw.pdgCode()) == PDG_t::kKPlus) {
1014+
int bin = hReweighthadron->FindBin(2, t2_raw.pt());
1015+
Reweighthadron = hReweighthadron->GetBinContent(bin);
1016+
}
1017+
acceptance = acceptance * Reweighthadron;
9691018
}
9701019
std::vector<float> fTransRange = fConfigDileptonHadronOptions.fConfigTransRange;
9711020
VarManager::FillEnergyCorrelatorsMC<THadronMassType>(t1_raw, t2_raw, VarManager::fgValues, fTransRange[0], fTransRange[1], 1. / acceptance);
972-
for (auto& sig : fGenMCSignals) {
973-
if (sig->CheckSignal(true, t1_raw)) {
1021+
for (auto& sig : fGenMCPairSignals) {
1022+
1023+
if (sig->CheckSignal(true, t1_raw, t2_raw)) {
9741024
if (!MixedEvent && !PionMass) {
9751025
fHistMan->FillHistClass(Form("MCTruthEenergyCorrelators_%s", sig->GetName()), VarManager::fgValues);
9761026
}

0 commit comments

Comments
 (0)