Skip to content

Commit 1c6d714

Browse files
YazhenLinCape
YazhenLin
and
Cape
authored
[PWGDQ] Add and fix some code for the energy correlator study (#15912)
Co-authored-by: Cape <cape@lindeMacBook-Pro.local>
1 parent 848b308 commit 1c6d714

4 files changed

Lines changed: 77 additions & 21 deletions

File tree

‎PWGDQ/Core/VarManager.cxx‎

Lines changed: 4 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -2253,6 +2253,7 @@ void VarManager::SetDefaultVarNames()
22532253
fgVarNamesMap["kMCCosChi"] = kMCCosChi;
22542254
fgVarNamesMap["kMCHadronPt"] = kMCHadronPt;
22552255
fgVarNamesMap["kMCWeight_before"] = kMCWeight_before;
2256+
fgVarNamesMap["kMCEWeight_before"] = kMCEWeight_before;
22562257
fgVarNamesMap["kMCdeltaeta"] = kMCdeltaeta;
22572258
fgVarNamesMap["kMCHadronPt"] = kMCHadronPt;
22582259
fgVarNamesMap["kMCHadronEta"] = kMCHadronEta;
@@ -2267,6 +2268,7 @@ void VarManager::SetDefaultVarNames()
22672268
fgVarNamesMap["kMCdeltaphi_randomPhi_toward"] = kMCdeltaphi_randomPhi_toward;
22682269
fgVarNamesMap["kMCdeltaphi_randomPhi_away"] = kMCdeltaphi_randomPhi_away;
22692270
fgVarNamesMap["kMCdeltaphi_randomPhi_trans"] = kMCdeltaphi_randomPhi_trans;
2271+
fgVarNamesMap["kMCHadronpt_randomPhi_trans"] = kMCHadronpt_randomPhi_trans;
22702272
fgVarNamesMap["kMCCosChi_gen"] = kMCCosChi_gen;
22712273
fgVarNamesMap["kMCWeight_gen"] = kMCWeight_gen;
22722274
fgVarNamesMap["kMCdeltaeta_gen"] = kMCdeltaeta_gen;
@@ -2504,6 +2506,7 @@ void VarManager::SetDefaultVarNames()
25042506
fgVarNamesMap["kWeight"] = kWeight;
25052507
fgVarNamesMap["kECWeight"] = kECWeight;
25062508
fgVarNamesMap["kEWeight_before"] = kEWeight_before;
2509+
fgVarNamesMap["kWeight_before"] = kWeight_before;
25072510
fgVarNamesMap["kPtDau"] = kPtDau;
25082511
fgVarNamesMap["kEtaDau"] = kEtaDau;
25092512
fgVarNamesMap["kPhiDau"] = kPhiDau;
@@ -2516,6 +2519,7 @@ void VarManager::SetDefaultVarNames()
25162519
fgVarNamesMap["kdeltaphi_randomPhi_trans"] = kdeltaphi_randomPhi_trans;
25172520
fgVarNamesMap["kdeltaphi_randomPhi_toward"] = kdeltaphi_randomPhi_toward;
25182521
fgVarNamesMap["kdeltaphi_randomPhi_away"] = kdeltaphi_randomPhi_away;
2522+
fgVarNamesMap["kPtDau_randomPhi_trans"] = kPtDau_randomPhi_trans;
25192523
fgVarNamesMap["kdileptonmass"] = kdileptonmass;
25202524
fgVarNamesMap["kNCorrelationVariables"] = kNCorrelationVariables;
25212525
fgVarNamesMap["kQuadMass"] = kQuadMass;

‎PWGDQ/Core/VarManager.h‎

Lines changed: 10 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -692,7 +692,9 @@ class VarManager : public TObject
692692
kMCdeltaphi_randomPhi_toward,
693693
kMCdeltaphi_randomPhi_away,
694694
kMCdeltaphi_randomPhi_trans,
695+
kMCHadronpt_randomPhi_trans,
695696
kMCWeight_before,
697+
kMCEWeight_before,
696698
kMCCosChi_gen,
697699
kMCWeight_gen,
698700
kMCdeltaeta_gen,
@@ -945,6 +947,7 @@ class VarManager : public TObject
945947
kPtDau,
946948
kCosTheta,
947949
kEWeight_before,
950+
kWeight_before,
948951
kCosChi_randomPhi_trans,
949952
kCosChi_randomPhi_toward,
950953
kCosChi_randomPhi_away,
@@ -955,6 +958,7 @@ class VarManager : public TObject
955958
kdeltaphi_randomPhi_toward,
956959
kdeltaphi_randomPhi_away,
957960
kdileptonmass,
961+
kPtDau_randomPhi_trans,
958962

959963
// Dilepton-track-track variables
960964
kQuadMass,
@@ -3370,6 +3374,7 @@ void VarManager::FillEnergyCorrelatorsMC(T const& track, T1 const& t1, float* va
33703374
values[kMCAccweight] = Accweight;
33713375
values[kMCCosChi] = CosChi;
33723376
values[kMCWeight_before] = t1.pt() / o2::constants::physics::MassJPsi * Accweight;
3377+
values[kMCEWeight_before] = t1.e() / o2::constants::physics::MassJPsi * Accweight;
33733378
values[kMCCosTheta] = CosTheta;
33743379
values[kMCdeltaphi] = deltaphi;
33753380
values[kMCdeltaeta] = deltaeta;
@@ -3397,6 +3402,7 @@ void VarManager::FillEnergyCorrelatorsMC(T const& track, T1 const& t1, float* va
33973402
randomPhi_toward = gRandom->Uniform(-o2::constants::math::PIHalf, 3. * o2::constants::math::PIHalf);
33983403
randomPhi_away = gRandom->Uniform(-o2::constants::math::PIHalf, 3. * o2::constants::math::PIHalf);
33993404

3405+
values[kMCHadronpt_randomPhi_trans] = v2.pt();
34003406
ROOT::Math::PtEtaPhiMVector v2_randomPhi_trans(v2.pt(), v2.eta(), randomPhi_trans, MassHadron);
34013407
values[kMCCosChi_randomPhi_trans] = LorentzTransformJpsihadroncosChi("coschi", v1, v2_randomPhi_trans);
34023408
values[kMCWeight_randomPhi_trans] = LorentzTransformJpsihadroncosChi("weight_boost", v1, v2_randomPhi_trans) / v1.M() * Accweight;
@@ -5877,7 +5883,8 @@ void VarManager::FillEnergyCorrelatorTriple(T1 const& lepton1, T2 const& lepton2
58775883
values[kWeight] = weight;
58785884
values[kECWeight] = E_boost / v1.M() * weight;
58795885
values[kCosTheta] = LorentzTransformJpsihadroncosChi("costheta", v1, v2);
5880-
values[kEWeight_before] = v2.Pt() / v1.M() * weight;
5886+
values[kEWeight_before] = v2.E() / v1.M() * weight;
5887+
values[kWeight_before] = v2.Pt() / v1.M() * weight;
58815888
values[kPtDau] = v2.pt();
58825889
values[kEtaDau] = v2.eta();
58835890
values[kPhiDau] = RecoDecay::constrainAngle(v2.phi(), -o2::constants::math::PIHalf);
@@ -5901,7 +5908,7 @@ void VarManager::FillEnergyCorrelatorTriple(T1 const& lepton1, T2 const& lepton2
59015908
randomPhi_trans = gRandom->Uniform(-o2::constants::math::PIHalf, 3. * o2::constants::math::PIHalf);
59025909
randomPhi_toward = gRandom->Uniform(-o2::constants::math::PIHalf, 3. * o2::constants::math::PIHalf);
59035910
randomPhi_away = gRandom->Uniform(-o2::constants::math::PIHalf, 3. * o2::constants::math::PIHalf);
5904-
5911+
values[kPtDau_randomPhi_trans] = v2.pt();
59055912
ROOT::Math::PtEtaPhiMVector v2_randomPhi_trans(v2.pt(), v2.eta(), randomPhi_trans, o2::constants::physics::MassPionCharged);
59065913
values[kCosChi_randomPhi_trans] = LorentzTransformJpsihadroncosChi("coschi", v1, v2_randomPhi_trans);
59075914
values[kWeight_randomPhi_trans] = LorentzTransformJpsihadroncosChi("weight_boost", v1, v2_randomPhi_trans) / v1.M() * weight;
@@ -5969,8 +5976,7 @@ void VarManager::FillEnergyCorrelatorsUnfoldingTriple(T1 const& lepton1, T2 cons
59695976
ROOT::Math::PtEtaPhiMVector v2_randomPhi_trans_rec(v2_rec.pt(), v2_rec.eta(), randomPhi_trans_rec, o2::constants::physics::MassPionCharged);
59705977
values[kMCCosChi_randomPhi_trans_rec] = LorentzTransformJpsihadroncosChi("coschi", v1_rec, v2_randomPhi_trans_rec);
59715978
values[kMCWeight_randomPhi_trans_rec] = LorentzTransformJpsihadroncosChi("weight_boost", v1_rec, v2_randomPhi_trans_rec) / v1_rec.M() * Effweight_rec;
5972-
float randomPhi_trans_gen = gRandom->Uniform(-o2::constants::math::PIHalf, 3. * o2::constants::math::PIHalf);
5973-
ROOT::Math::PtEtaPhiMVector v2_randomPhi_trans_gen(v2_gen.pt(), v2_gen.eta(), randomPhi_trans_gen, MassHadron);
5979+
ROOT::Math::PtEtaPhiMVector v2_randomPhi_trans_gen(v2_gen.pt(), v2_gen.eta(), randomPhi_trans_rec, MassHadron);
59745980
values[kMCCosChi_randomPhi_trans_gen] = LorentzTransformJpsihadroncosChi("coschi", v1_gen, v2_randomPhi_trans_gen);
59755981
values[kMCWeight_randomPhi_trans_gen] = LorentzTransformJpsihadroncosChi("weight_boost", v1_gen, v2_randomPhi_trans_gen) / v1_gen.M() * Accweight_gen;
59765982
}

‎PWGDQ/Tasks/dqEnergyCorrelator_direct.cxx‎

Lines changed: 54 additions & 8 deletions
Original file line numberDiff line numberDiff line change
@@ -104,6 +104,8 @@ struct AnalysisEnergyCorrelator {
104104
Configurable<std::string> fConfigTrackCuts{"cfgTrackCuts", "electronSelection1_ionut", "Comma separated list of barrel track cuts for electrons"};
105105
Configurable<std::string> fConfigTrackCutsJSON{"cfgTrackCutsJSON", "", "Additional track cuts in JSON"};
106106
Configurable<std::string> fConfigAddTrackHistogram{"cfgAddTrackHistogram", "", "Track histograms"};
107+
Configurable<std::string> fConfigMCRecTrackSignals{"cfgMCRecTrackSignals", "", "Comma separated list of MC signals (reconstructed)"};
108+
Configurable<std::string> fConfigMCRecTrackSignalsJSON{"cfgMCRecTrackSignalsJSON", "", "Additional list of MC signals (reconstructed) via JSON"};
107109
Configurable<bool> fConfigTrackQA{"cfgTrackQA", false, "If true, fill Track QA histograms"};
108110
} fConfigTrackOptions;
109111

@@ -155,12 +157,14 @@ struct AnalysisEnergyCorrelator {
155157
std::vector<TString> fTrackCutNames;
156158
std::vector<TString> fHadronCutNames;
157159
std::vector<TString> fHistNamesReco;
160+
std::vector<TString> fHistNamesMCMatched;
158161

159162
std::map<int, std::vector<TString>> fTrackHistNames;
160163
std::map<int, std::vector<TString>> fBarrelHistNamesMCmatched;
161164
std::map<int64_t, bool> fSelMap;
162165

163-
std::vector<MCSignal*> fRecMCSignals; // MC signals for reconstructed pairs
166+
std::vector<MCSignal*> fRecMCTrackSignals; // MC signals for reconstructed tracks
167+
std::vector<MCSignal*> fRecMCSignals; // MC signals for reconstructed pairs
164168
std::vector<MCSignal*> fGenMCSignals;
165169
std::vector<MCSignal*> fRecMCTripleSignals; // MC signals for reconstructed triples
166170

@@ -175,8 +179,8 @@ struct AnalysisEnergyCorrelator {
175179

176180
TH2F* hAcceptance_rec;
177181
TH2F* hAcceptance_gen;
178-
TH1F* hEfficiency_dilepton;
179-
TH1F* hEfficiency_hadron;
182+
TH2F* hEfficiency_dilepton;
183+
TH2F* hEfficiency_hadron;
180184
TH1F* hMasswindow;
181185

182186
void init(o2::framework::InitContext& context)
@@ -332,6 +336,28 @@ struct AnalysisEnergyCorrelator {
332336
}
333337
}
334338

339+
// TODO: create a std::vector of hist classes to be used at Fill time, to avoid using Form in the process function
340+
TString sigRecTrackNamesStr = fConfigTrackOptions.fConfigMCRecTrackSignals.value;
341+
std::unique_ptr<TObjArray> objRecTrackSigArray(sigRecTrackNamesStr.Tokenize(","));
342+
for (int isig = 0; isig < objRecTrackSigArray->GetEntries(); isig++) {
343+
MCSignal* sig_RecTrack = o2::aod::dqmcsignals::GetMCSignal(objRecTrackSigArray->At(isig)->GetName());
344+
if (sig_RecTrack) {
345+
fRecMCTrackSignals.push_back(sig_RecTrack);
346+
}
347+
}
348+
349+
// Add the MCSignals from the JSON config
350+
TString addRecTrackSignalsGenStr = fConfigTrackOptions.fConfigMCRecTrackSignalsJSON.value;
351+
if (addRecTrackSignalsGenStr != "") {
352+
std::vector<MCSignal*> addMCRecTrackSignals = dqmcsignals::GetMCSignalsFromJSON(addRecTrackSignalsGenStr.Data());
353+
for (auto& mcIt : addMCRecTrackSignals) {
354+
if (mcIt->GetNProngs() > 2) { // NOTE: only 2 prong signals
355+
continue;
356+
}
357+
fRecMCTrackSignals.push_back(mcIt);
358+
}
359+
}
360+
335361
VarManager::SetUseVars(AnalysisCut::fgUsedVars);
336362

337363
fHistMan = new HistogramManager("analysisHistos", "", VarManager::kNVars);
@@ -352,6 +378,11 @@ struct AnalysisEnergyCorrelator {
352378
TString nameStr = Form("AssocsBarrel_%s", cut->GetName());
353379
fHistNamesReco.push_back(nameStr);
354380
histClasses += Form("%s;", nameStr.Data());
381+
for (auto& sig : fRecMCTrackSignals) {
382+
TString nameStr2 = Form("AssocsBarrelMatched_%s_%s", cut->GetName(), sig->GetName());
383+
fHistNamesMCMatched.push_back(nameStr2);
384+
histClasses += Form("%s;", nameStr2.Data());
385+
}
355386
}
356387
DefineHistograms(fHistMan, histClasses.Data(), fConfigTrackOptions.fConfigAddTrackHistogram.value.data());
357388
}
@@ -435,8 +466,8 @@ struct AnalysisEnergyCorrelator {
435466
if (!listAccs) {
436467
LOG(fatal) << "Problem getting TList object with efficiencies!";
437468
}
438-
hEfficiency_dilepton = static_cast<TH1F*>(listAccs->FindObject("hEfficiency_dilepton"));
439-
hEfficiency_hadron = static_cast<TH1F*>(listAccs->FindObject("hEfficiency_hadron"));
469+
hEfficiency_dilepton = static_cast<TH2F*>(listAccs->FindObject("hEfficiency_dilepton"));
470+
hEfficiency_hadron = static_cast<TH2F*>(listAccs->FindObject("hEfficiency_hadron"));
440471
hAcceptance_rec = static_cast<TH2F*>(listAccs->FindObject("hAcceptance_rec"));
441472
hAcceptance_gen = static_cast<TH2F*>(listAccs->FindObject("hAcceptance_gen"));
442473
hMasswindow = static_cast<TH1F*>(listAccs->FindObject("hMasswindow"));
@@ -479,9 +510,9 @@ struct AnalysisEnergyCorrelator {
479510
float deltaphi = RecoDecay::constrainAngle(dilepton_phi - hadron_phi, -0.5 * o2::constants::math::PI);
480511
Effweight_rec = hAcceptance_rec->Interpolate(dilepton_eta - hadron_eta, deltaphi);
481512
Accweight_gen = hAcceptance_gen->Interpolate(dilepton_eta - hadron_eta, deltaphi);
482-
float Effdilepton = hEfficiency_dilepton->Interpolate(VarManager::fgValues[VarManager::kPt]);
513+
float Effdilepton = hEfficiency_dilepton->Interpolate(VarManager::fgValues[VarManager::kPt], dilepton_eta);
483514
float Masswindow = hMasswindow->Interpolate(VarManager::fgValues[VarManager::kPt]);
484-
float Effhadron = hEfficiency_hadron->Interpolate(hadron.pt());
515+
float Effhadron = hEfficiency_hadron->Interpolate(hadron.pt(), hadron.eta());
485516
Accweight_gen = Accweight_gen * Effdilepton * Effhadron;
486517
if (fConfigDileptonHadronOptions.fConfigApplyEfficiencyME) {
487518
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
@@ -647,17 +678,32 @@ struct AnalysisEnergyCorrelator {
647678
fHistMan->FillHistClass("AssocsBarrel_BeforeCuts", VarManager::fgValues);
648679
}
649680

681+
uint32_t mcDecision = static_cast<uint32_t>(0);
682+
// run MC matching for this pair
683+
int isig = 0;
684+
mcDecision = 0;
685+
for (auto sig = fRecMCTrackSignals.begin(); sig != fRecMCTrackSignals.end(); sig++, isig++) {
686+
if (t1.has_mcParticle()) {
687+
if ((*sig)->CheckSignal(true, t1.mcParticle())) {
688+
mcDecision |= (static_cast<uint32_t>(1) << isig);
689+
}
690+
}
691+
}
650692
// Apply electron cuts and fill histograms
651693
int iCut1 = 0;
652694
for (auto cut1 = fTrackCuts.begin(); cut1 != fTrackCuts.end(); cut1++, iCut1++) {
653695
if ((*cut1)->IsSelected(VarManager::fgValues)) {
654696
filter1 |= (static_cast<uint32_t>(1) << iCut1);
655697
if (fConfigTrackOptions.fConfigTrackQA) {
656698
fHistMan->FillHistClass(fHistNamesReco[iCut1], VarManager::fgValues);
699+
for (size_t isig = 0; isig < fRecMCTrackSignals.size(); isig++) { // loop over MC signals
700+
if (mcDecision & (static_cast<uint32_t>(1) << isig)) {
701+
fHistMan->FillHistClass(fHistNamesMCMatched[iCut1 * fRecMCTrackSignals.size() + isig], VarManager::fgValues); // matched signal
702+
}
703+
}
657704
}
658705
}
659706
}
660-
661707
// Check opposite charge with t2
662708
for (auto& a2 : groupedAssocs) {
663709
auto t2 = a2.template track_as<MyBarrelTracksWithCov>();

‎PWGDQ/Tasks/tableReader_withAssoc.cxx‎

Lines changed: 9 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -3493,8 +3493,8 @@ struct AnalysisDileptonTrack {
34933493

34943494
TH2F* hAcceptance_rec;
34953495
TH2F* hAcceptance_gen;
3496-
TH1F* hEfficiency_dilepton;
3497-
TH1F* hEfficiency_hadron;
3496+
TH2F* hEfficiency_dilepton;
3497+
TH2F* hEfficiency_hadron;
34983498
TH1F* hMasswindow;
34993499

35003500
void init(o2::framework::InitContext& context)
@@ -3760,8 +3760,8 @@ struct AnalysisDileptonTrack {
37603760
if (!listAccs) {
37613761
LOG(fatal) << "Problem getting TList object with efficiencies!";
37623762
}
3763-
hEfficiency_dilepton = static_cast<TH1F*>(listAccs->FindObject("hEfficiency_dilepton"));
3764-
hEfficiency_hadron = static_cast<TH1F*>(listAccs->FindObject("hEfficiency_hadron"));
3763+
hEfficiency_dilepton = static_cast<TH2F*>(listAccs->FindObject("hEfficiency_dilepton"));
3764+
hEfficiency_hadron = static_cast<TH2F*>(listAccs->FindObject("hEfficiency_hadron"));
37653765
hAcceptance_rec = static_cast<TH2F*>(listAccs->FindObject("hAcceptance_rec"));
37663766
hAcceptance_gen = static_cast<TH2F*>(listAccs->FindObject("hAcceptance_gen"));
37673767
hMasswindow = static_cast<TH1F*>(listAccs->FindObject("hMasswindow"));
@@ -3854,9 +3854,9 @@ struct AnalysisDileptonTrack {
38543854
float hadron_phi = track.phi();
38553855
float deltaphi = RecoDecay::constrainAngle(dilepton_phi - hadron_phi, -0.5 * o2::constants::math::PI);
38563856
Effweight_rec = hAcceptance_rec->Interpolate(dilepton_eta - hadron_eta, deltaphi);
3857-
float Effdilepton = hEfficiency_dilepton->Interpolate(dilepton.pt());
3857+
float Effdilepton = hEfficiency_dilepton->Interpolate(dilepton.pt(), dilepton_eta);
38583858
float Masswindow = hMasswindow->Interpolate(dilepton.pt());
3859-
float Effhadron = hEfficiency_hadron->Interpolate(track.pt());
3859+
float Effhadron = hEfficiency_hadron->Interpolate(track.pt(), hadron_eta);
38603860
Effweight_rec = Effweight_rec * Effdilepton * Effhadron * Masswindow;
38613861
}
38623862
std::vector<float> fTransRange = fConfigTransRange;
@@ -4088,17 +4088,17 @@ struct AnalysisDileptonTrack {
40884088
float hadron_phi = track.phi();
40894089
float deltaphi = RecoDecay::constrainAngle(dilepton_phi - hadron_phi, -0.5 * o2::constants::math::PI);
40904090
Effweight_rec = hAcceptance_rec->Interpolate(dilepton_eta - hadron_eta, deltaphi);
4091-
float Effdilepton = hEfficiency_dilepton->Interpolate(dilepton.pt());
4091+
float Effdilepton = hEfficiency_dilepton->Interpolate(dilepton.pt(), dilepton_eta);
40924092
float Masswindow = hMasswindow->Interpolate(dilepton.pt());
4093-
float Effhadron = hEfficiency_hadron->Interpolate(track.pt());
4093+
float Effhadron = hEfficiency_hadron->Interpolate(track.pt(), hadron_eta);
40944094
if (fConfigApplyEfficiencyME) {
40954095
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
40964096
} else {
40974097
Effweight_rec = Effweight_rec * Effdilepton * Effhadron * Masswindow; // apply acceptance and efficiency correction for the real pairs
40984098
}
40994099
}
41004100
std::vector<float> fTransRange = fConfigTransRange;
4101-
VarManager::FillEnergyCorrelatorTriple(lepton1, lepton2, track, fValuesHadron, fTransRange[0], fTransRange[1], fConfigApplyMassEC, fMassBkg->GetRandom(), 1. / Effweight_rec);
4101+
VarManager::FillEnergyCorrelatorTriple(lepton1, lepton2, track, VarManager::fgValues, fTransRange[0], fTransRange[1], fConfigApplyMassEC, fMassBkg->GetRandom(), 1. / Effweight_rec);
41024102

41034103
// loop over dilepton leg cuts and track cuts and fill histograms separately for each combination
41044104
for (int icut = 0; icut < fNCuts; icut++) {

0 commit comments

Comments
 (0)