Skip to content

Commit e9f8ed5

Browse files
YazhenLinYazhenLin
authored
Add files via upload
1 parent da70704 commit e9f8ed5

1 file changed

Lines changed: 67 additions & 16 deletions

File tree

‎PWGDQ/Tasks/dqEnergyCorrelator_direct.cxx‎

Lines changed: 67 additions & 16 deletions
Original file line numberDiff line numberDiff line change
@@ -130,8 +130,10 @@ 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"};
137139
Configurable<bool> fConfigContainlepton{"cfgContainlepton", false, "If true, require the hadron to contain the lepton in its decay tree for the energy correlator study"};
@@ -167,6 +169,7 @@ struct AnalysisEnergyCorrelator {
167169
std::vector<MCSignal*> fRecMCTrackSignals; // MC signals for reconstructed tracks
168170
std::vector<MCSignal*> fRecMCSignals; // MC signals for reconstructed pairs
169171
std::vector<MCSignal*> fGenMCSignals;
172+
std::vector<MCSignal*> fGenMCPairSignals;
170173
std::vector<MCSignal*> fRecMCTripleSignals; // MC signals for reconstructed triples
171174

172175
Service<o2::ccdb::BasicCCDBManager> fCCDB;
@@ -183,6 +186,7 @@ struct AnalysisEnergyCorrelator {
183186
TH2F* hEfficiency_dilepton;
184187
TH2F* hEfficiency_hadron;
185188
TH1F* hMasswindow;
189+
TH2F* hReweighthadron;
186190

187191
void init(o2::framework::InitContext& context)
188192
{
@@ -337,6 +341,29 @@ struct AnalysisEnergyCorrelator {
337341
}
338342
}
339343

344+
// Add histogram classes for each specified MCsignal at the generator level
345+
// TODO: create a std::vector of hist classes to be used at Fill time, to avoid using Form in the process function
346+
TString sigpairGenNamesStr = fConfigDileptonHadronOptions.fConfigMCGenPairSignals.value;
347+
std::unique_ptr<TObjArray> objGenPairSigArray(sigpairGenNamesStr.Tokenize(","));
348+
for (int isig = 0; isig < objGenPairSigArray->GetEntries(); isig++) {
349+
MCSignal* sig = o2::aod::dqmcsignals::GetMCSignal(objGenPairSigArray->At(isig)->GetName());
350+
if (sig) {
351+
fGenMCPairSignals.push_back(sig);
352+
}
353+
}
354+
355+
// Add the MCSignals from the JSON config
356+
TString addMCSignalsPairGenStr = fConfigDileptonHadronOptions.fConfigMCGenPairSignalsJSON.value;
357+
if (addMCSignalsPairGenStr != "") {
358+
std::vector<MCSignal*> addMCSignals = dqmcsignals::GetMCSignalsFromJSON(addMCSignalsPairGenStr.Data());
359+
for (auto& mcIt : addMCSignals) {
360+
if (mcIt->GetNProngs() != 2) { // NOTE: only 2 prong signals
361+
continue;
362+
}
363+
fGenMCPairSignals.push_back(mcIt);
364+
}
365+
}
366+
340367
// TODO: create a std::vector of hist classes to be used at Fill time, to avoid using Form in the process function
341368
TString sigRecTrackNamesStr = fConfigTrackOptions.fConfigMCRecTrackSignals.value;
342369
std::unique_ptr<TObjArray> objRecTrackSigArray(sigRecTrackNamesStr.Tokenize(","));
@@ -440,6 +467,14 @@ struct AnalysisEnergyCorrelator {
440467
if (sig->GetNProngs() == 1) {
441468
if (isMCGen_energycorrelators) {
442469
DefineHistograms(fHistMan, Form("MCTruthGenSel_%s", sig->GetName()), "");
470+
}
471+
}
472+
}
473+
474+
for (auto& sig : fGenMCPairSignals) {
475+
LOG(info) << "Defining histograms for pair signal: " << sig->GetNProngs();
476+
if (sig->GetNProngs() == 2) {
477+
if (isMCGen_energycorrelators) {
443478
DefineHistograms(fHistMan, Form("MCTruthEenergyCorrelators_%s", sig->GetName()), "");
444479
DefineHistograms(fHistMan, Form("MCTruthEenergyCorrelators_Pion_%s", sig->GetName()), "");
445480
}
@@ -472,7 +507,8 @@ struct AnalysisEnergyCorrelator {
472507
hAcceptance_rec = static_cast<TH2F*>(listAccs->FindObject("hAcceptance_rec"));
473508
hAcceptance_gen = static_cast<TH2F*>(listAccs->FindObject("hAcceptance_gen"));
474509
hMasswindow = static_cast<TH1F*>(listAccs->FindObject("hMasswindow"));
475-
if (!hAcceptance_rec || !hAcceptance_gen || !hEfficiency_dilepton || !hEfficiency_hadron || !hMasswindow) {
510+
hReweighthadron = static_cast<TH2F*>(listAccs->FindObject("hReweighthadron"));
511+
if (!hAcceptance_rec || !hAcceptance_gen || !hEfficiency_dilepton || !hEfficiency_hadron || !hMasswindow || !hReweighthadron) {
476512
LOG(fatal) << "Problem getting histograms from the TList object with efficiencies!";
477513
}
478514
}
@@ -533,11 +569,22 @@ struct AnalysisEnergyCorrelator {
533569
float Effdilepton = GetSafeInterpolationWeight(hEfficiency_dilepton, dilepton_rap, dilepton_pt);
534570
float Effhadron = GetSafeInterpolationWeight(hEfficiency_hadron, hadron_eta, hadron.pt());
535571
float Masswindow = hMasswindow->Interpolate(dilepton_pt);
536-
Accweight_gen = Accweight_gen * Effdilepton * Effhadron;
572+
float Reweighthadron = 1.0f;
573+
if (abs(hadronMC.pdgCode()) == PDG_t::kPiPlus) {
574+
int bin = hReweighthadron->FindBin(0, hadron.pt());
575+
Reweighthadron = hReweighthadron->GetBinContent(bin);
576+
} else if (abs(hadronMC.pdgCode()) == PDG_t::kProton) {
577+
int bin = hReweighthadron->FindBin(1, hadron.pt());
578+
Reweighthadron = hReweighthadron->GetBinContent(bin);
579+
} else if (abs(hadronMC.pdgCode()) == PDG_t::kKPlus) {
580+
int bin = hReweighthadron->FindBin(2, hadron.pt());
581+
Reweighthadron = hReweighthadron->GetBinContent(bin);
582+
}
583+
Accweight_gen = Accweight_gen * Effdilepton * Effhadron * Reweighthadron;
537584
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
585+
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
539586
} else {
540-
Effweight_rec = Effweight_rec * Effdilepton * Effhadron * Masswindow; // apply acceptance and efficiency correction for the real pairs
587+
Effweight_rec = Effweight_rec * Effdilepton * Effhadron * Masswindow * Reweighthadron; // apply acceptance and efficiency correction for the real pairs
541588
}
542589
}
543590

@@ -936,25 +983,16 @@ struct AnalysisEnergyCorrelator {
936983
continue;
937984
// for the energy correlators
938985
for (auto& t2 : groupedMCTracks2) {
939-
auto t2_raw = groupedMCTracks2.rawIteratorAt(t2.globalIndex());
986+
auto t2_raw = mcTracks.rawIteratorAt(t2.globalIndex());
940987
if (t2.mcCollisionId() != event2.mcCollisionId()) { // check that the mc track belongs to the same mc collision as the reconstructed event
941988
continue;
942989
}
943-
if (!t2_raw.isPhysicalPrimary()) {
944-
continue;
945-
}
946990
if (t2_raw.has_mothers()) {
947991
auto mother_raw = t2_raw.template mothers_first_as<McParticles>();
948992
if (mother_raw.globalIndex() == t1_raw.globalIndex()) {
949993
continue;
950994
}
951995
}
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-
}
958996
if (t2_raw.pt() < fConfigDileptonHadronOptions.fConfigMCGenHadronPtMin.value || std::abs(t2_raw.eta()) > fConfigDileptonHadronOptions.fConfigMCGenHadronEtaAbs.value) {
959997
continue;
960998
}
@@ -966,11 +1004,24 @@ struct AnalysisEnergyCorrelator {
9661004
float hadron_phi = t2_raw.phi();
9671005
float deltaphi = RecoDecay::constrainAngle(dilepton_phi - hadron_phi, -0.5 * o2::constants::math::PI);
9681006
acceptance = hAcceptance_gen->Interpolate(dilepton_eta - hadron_eta, deltaphi);
1007+
float Reweighthadron = 1.0f;
1008+
if (abs(t2_raw.pdgCode()) == PDG_t::kPiPlus) {
1009+
int bin = hReweighthadron->FindBin(0, t2_raw.pt());
1010+
Reweighthadron = hReweighthadron->GetBinContent(bin);
1011+
} else if (abs(t2_raw.pdgCode()) == PDG_t::kProton) {
1012+
int bin = hReweighthadron->FindBin(1, t2_raw.pt());
1013+
Reweighthadron = hReweighthadron->GetBinContent(bin);
1014+
} else if (abs(t2_raw.pdgCode()) == PDG_t::kKPlus) {
1015+
int bin = hReweighthadron->FindBin(2, t2_raw.pt());
1016+
Reweighthadron = hReweighthadron->GetBinContent(bin);
1017+
}
1018+
acceptance = acceptance * Reweighthadron;
9691019
}
9701020
std::vector<float> fTransRange = fConfigDileptonHadronOptions.fConfigTransRange;
9711021
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)) {
1022+
for (auto& sig : fGenMCPairSignals) {
1023+
1024+
if (sig->CheckSignal(true, t1_raw, t2_raw)) {
9741025
if (!MixedEvent && !PionMass) {
9751026
fHistMan->FillHistClass(Form("MCTruthEenergyCorrelators_%s", sig->GetName()), VarManager::fgValues);
9761027
}

0 commit comments

Comments
 (0)