Skip to content

Commit 542ed5b

Browse files
ercolessiercolessi
authored
[PWGLF] nuclei-proton correlations: Add process Gen for simulation studies (#9622)
1 parent 1480d7a commit 542ed5b

File tree

1 file changed

+194
-13
lines changed

1 file changed

+194
-13
lines changed

‎PWGLF/Tasks/Nuspex/hadronnucleicorrelation.cxx‎

Lines changed: 194 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -59,6 +59,8 @@ struct hadronnucleicorrelation {
5959
Configurable<bool> doQA{"doQA", true, "save QA histograms"};
6060
Configurable<bool> doMCQA{"doMCQA", false, "save MC QA histograms"};
6161
Configurable<bool> isMC{"isMC", false, "is MC"};
62+
Configurable<bool> isMCGen{"isMCGen", false, "is isMCGen"};
63+
Configurable<bool> isPrim{"isPrim", true, "is isPrim"};
6264
Configurable<bool> mcCorrelation{"mcCorrelation", false, "true: build the correlation function only for SE"};
6365
Configurable<bool> docorrection{"docorrection", false, "do efficiency correction"};
6466

@@ -101,25 +103,32 @@ struct hadronnucleicorrelation {
101103
ConfigurableAxis AxisNSigma{"AxisNSigma", {35, -7.f, 7.f}, "n#sigma"};
102104

103105
using FilteredCollisions = soa::Filtered<aod::SingleCollSels>;
104-
using FilteredTracks = soa::Filtered<soa::Join<aod::SingleTrackSels_v2, aod::SingleTrkExtras, aod::SinglePIDEls_v0>>; // old tables (v2)
105-
using FilteredTracksMC = soa::Filtered<soa::Join<aod::SingleTrackSels_v2, aod::SingleTrkMCs, aod::SingleTrkExtras, aod::SinglePIDEls_v0>>; // old tables (v2)
106-
107-
// using FilteredTracks = soa::Filtered<soa::Join<aod::SingleTrackSels, aod::SingleTrkExtras, aod::SinglePIDEls, aod::SinglePIDPis, aod::SinglePIDKas, aod::SinglePIDPrs, aod::SinglePIDDes, aod::SinglePIDHes>>; // new tables (v3)
108-
// using FilteredTracksMC = soa::Filtered<soa::Join<aod::SingleTrackSels, aod::SingleTrkMCs, aod::SingleTrkExtras, aod::SinglePIDEls, aod::SinglePIDPis, aod::SinglePIDKas, aod::SinglePIDPrs, aod::SinglePIDDes, aod::SinglePIDHes>>; // new tables (v3)
106+
using SimCollisions = aod::McCollisions;
107+
using SimParticles = aod::McParticles;
108+
using FilteredTracks = soa::Filtered<soa::Join<aod::SingleTrackSels, aod::SingleTrkExtras, aod::SinglePIDEls, aod::SinglePIDPrs, aod::SinglePIDDes, aod::SinglePIDHes>>; // new tables (v3)
109+
using FilteredTracksMC = soa::Filtered<soa::Join<aod::SingleTrackSels, aod::SingleTrkMCs, aod::SingleTrkExtras, aod::SinglePIDEls, aod::SinglePIDPrs, aod::SinglePIDDes, aod::SinglePIDHes>>; // new tables (v3)
109110

110111
HistogramRegistry registry{"registry"};
111112
HistogramRegistry QA{"QA"};
112113

113114
typedef std::shared_ptr<FilteredTracks::iterator> trkType;
114115
typedef std::shared_ptr<FilteredTracksMC::iterator> trkTypeMC;
116+
typedef std::shared_ptr<SimParticles::iterator> partTypeMC;
115117
typedef std::shared_ptr<FilteredCollisions::iterator> colType;
118+
typedef std::shared_ptr<SimCollisions::iterator> MCcolType;
116119

117120
// key: int64_t - value: vector of trkType objects
118121
std::map<int64_t, std::vector<trkType>> selectedtracks_p;
119122
std::map<int64_t, std::vector<trkType>> selectedtracks_d;
120123
std::map<int64_t, std::vector<trkType>> selectedtracks_antid;
121124
std::map<int64_t, std::vector<trkType>> selectedtracks_antip;
122125

126+
// key: int64_t - value: vector of trkType objects
127+
std::map<int64_t, std::vector<partTypeMC>> selectedparticlesMC_d;
128+
std::map<int64_t, std::vector<partTypeMC>> selectedparticlesMC_p;
129+
std::map<int64_t, std::vector<partTypeMC>> selectedparticlesMC_antid;
130+
std::map<int64_t, std::vector<partTypeMC>> selectedparticlesMC_antip;
131+
123132
// key: int64_t - value: vector of trkType objects
124133
std::map<int64_t, std::vector<trkTypeMC>> selectedtracksMC_d;
125134
std::map<int64_t, std::vector<trkTypeMC>> selectedtracksMC_p;
@@ -134,6 +143,8 @@ struct hadronnucleicorrelation {
134143
// for each key I have a vector of collisions
135144
std::map<std::pair<int, float>, std::vector<colType>> mixbins_antidantip;
136145
std::map<std::pair<int, float>, std::vector<colType>> mixbinsPID_antidantip;
146+
std::map<float, std::vector<MCcolType>> mixbinsMC_antidantip;
147+
std::map<float, std::vector<MCcolType>> mixbinsMC_dp;
137148

138149
std::vector<std::shared_ptr<TH3>> hEtaPhi_AntiDeAntiPr_SE;
139150
std::vector<std::shared_ptr<TH3>> hEtaPhi_AntiDeAntiPr_ME;
@@ -197,16 +208,10 @@ struct hadronnucleicorrelation {
197208
for (int i = 0; i < nBinspT; i++) {
198209
auto htempSERec_AntiDeAntiPr = registry.add<TH3>(Form("hEtaPhiRec_AntiDeAntiPr_SE_pt%02.0f%02.0f", pTBins.value.at(i) * 10, pTBins.value.at(i + 1) * 10),
199210
Form("Rec #Delta#eta#Delta#phi (%.1f<p_{T} #bar{d} <%.1f GeV/c)", pTBins.value.at(i), pTBins.value.at(i + 1)), {HistType::kTH3F, {DeltaEtaAxis, DeltaPhiAxis, ptBinnedAxis}});
200-
auto htempSEGen_AntiDeAntiPr = registry.add<TH3>(Form("hEtaPhiGen_AntiDeAntiPr_SE_pt%02.0f%02.0f", pTBins.value.at(i) * 10, pTBins.value.at(i + 1) * 10),
201-
Form("Gen #Delta#eta#Delta#phi (%.1f<p_{T} #bar{d} <%.1f GeV/c)", pTBins.value.at(i), pTBins.value.at(i + 1)), {HistType::kTH3F, {DeltaEtaAxis, DeltaPhiAxis, ptBinnedAxis}});
202211
hEtaPhiRec_AntiDeAntiPr_SE.push_back(std::move(htempSERec_AntiDeAntiPr));
203-
hEtaPhiGen_AntiDeAntiPr_SE.push_back(std::move(htempSEGen_AntiDeAntiPr));
204212
auto htempMERec_AntiDeAntiPr = registry.add<TH3>(Form("hEtaPhiRec_AntiDeAntiPr_ME_pt%02.0f%02.0f", pTBins.value.at(i) * 10, pTBins.value.at(i + 1) * 10),
205213
Form("Rec #Delta#eta#Delta#phi (%.1f<p_{T} #bar{d} <%.1f GeV/c)", pTBins.value.at(i), pTBins.value.at(i + 1)), {HistType::kTH3F, {DeltaEtaAxis, DeltaPhiAxis, ptBinnedAxis}});
206-
auto htempMEGen_AntiDeAntiPr = registry.add<TH3>(Form("hEtaPhiGen_AntiDeAntiPr_ME_pt%02.0f%02.0f", pTBins.value.at(i) * 10, pTBins.value.at(i + 1) * 10),
207-
Form("Gen #Delta#eta#Delta#phi (%.1f<p_{T} #bar{d} <%.1f GeV/c)", pTBins.value.at(i), pTBins.value.at(i + 1)), {HistType::kTH3F, {DeltaEtaAxis, DeltaPhiAxis, ptBinnedAxis}});
208214
hEtaPhiRec_AntiDeAntiPr_ME.push_back(std::move(htempMERec_AntiDeAntiPr));
209-
hEtaPhiGen_AntiDeAntiPr_ME.push_back(std::move(htempMEGen_AntiDeAntiPr));
210215

211216
auto htempPIDSERec_AntiDeAntiPr = registry.add<TH3>(Form("hPIDEtaPhiRec_AntiDeAntiPr_SE_pt%02.0f%02.0f", pTBins.value.at(i) * 10, pTBins.value.at(i + 1) * 10),
212217
Form("Rec #Delta#eta#Delta#phi (%.1f<p_{T} #bar{d} <%.1f GeV/c)", pTBins.value.at(i), pTBins.value.at(i + 1)), {HistType::kTH3F, {DeltaEtaAxis, DeltaPhiAxis, ptBinnedAxis}});
@@ -223,6 +228,17 @@ struct hadronnucleicorrelation {
223228
}
224229
}
225230

231+
if (isMCGen || mcCorrelation) {
232+
for (int i = 0; i < nBinspT; i++) {
233+
auto htempSEGen_AntiDeAntiPr = registry.add<TH3>(Form("hEtaPhiGen_AntiDeAntiPr_SE_pt%02.0f%02.0f", pTBins.value.at(i) * 10, pTBins.value.at(i + 1) * 10),
234+
Form("Gen #Delta#eta#Delta#phi (%.1f<p_{T} #bar{d} <%.1f GeV/c)", pTBins.value.at(i), pTBins.value.at(i + 1)), {HistType::kTH3F, {DeltaEtaAxis, DeltaPhiAxis, ptBinnedAxis}});
235+
hEtaPhiGen_AntiDeAntiPr_SE.push_back(std::move(htempSEGen_AntiDeAntiPr));
236+
auto htempMEGen_AntiDeAntiPr = registry.add<TH3>(Form("hEtaPhiGen_AntiDeAntiPr_ME_pt%02.0f%02.0f", pTBins.value.at(i) * 10, pTBins.value.at(i + 1) * 10),
237+
Form("Gen #Delta#eta#Delta#phi (%.1f<p_{T} #bar{d} <%.1f GeV/c)", pTBins.value.at(i), pTBins.value.at(i + 1)), {HistType::kTH3F, {DeltaEtaAxis, DeltaPhiAxis, ptBinnedAxis}});
238+
hEtaPhiGen_AntiDeAntiPr_ME.push_back(std::move(htempMEGen_AntiDeAntiPr));
239+
}
240+
}
241+
226242
if (!isMC) {
227243
for (int i = 0; i < nBinspT; i++) {
228244
auto htempSE_AntiDeAntiPr = registry.add<TH3>(Form("hEtaPhi_AntiDeAntiPr_SE_pt%02.0f%02.0f", pTBins.value.at(i) * 10, pTBins.value.at(i + 1) * 10), Form("Raw #Delta#eta#Delta#phi (%.1f<p_{T} #bar{d} <%.1f GeV/c)", pTBins.value.at(i), pTBins.value.at(i + 1)), {HistType::kTH3F, {DeltaEtaAxis, DeltaPhiAxis, ptBinnedAxis}});
@@ -341,6 +357,12 @@ struct hadronnucleicorrelation {
341357
registry.add("hReco_Pt_Deuteron_TPCEl_or_TOF", "Reco (anti)protons in reco collisions", {HistType::kTH1F, {pTAxis_small}});
342358
}
343359
}
360+
361+
if (isMCGen) {
362+
registry.add("Generated/hNEventsMC", "hNEventsMC", {HistType::kTH1D, {{2, 0.f, 2.f}}});
363+
registry.get<TH1>(HIST("Generated/hNEventsMC"))->GetXaxis()->SetBinLabel(1, "All");
364+
registry.get<TH1>(HIST("Generated/hNEventsMC"))->GetXaxis()->SetBinLabel(2, "|z_{vtx}| < 10 cm");
365+
}
344366
}
345367

346368
// Filters
@@ -349,8 +371,8 @@ struct hadronnucleicorrelation {
349371
o2::aod::singletrackselector::unPack<singletrackselector::binning::chi2>(o2::aod::singletrackselector::storedTpcChi2NCl) <= max_chi2_TPC &&
350372
o2::aod::singletrackselector::unPack<singletrackselector::binning::rowsOverFindable>(o2::aod::singletrackselector::storedTpcCrossedRowsOverFindableCls) >= min_TPC_nCrossedRowsOverFindableCls &&
351373
o2::aod::singletrackselector::unPack<singletrackselector::binning::chi2>(o2::aod::singletrackselector::storedItsChi2NCl) <= max_chi2_ITS &&
352-
nabs(o2::aod::singletrackselector::unPack<singletrackselector::binning::dca>(o2::aod::singletrackselector::storedDcaXY)) <= max_dcaxy && // For now no filtering on the DCAxy or DCAz (casting not supported)
353-
nabs(o2::aod::singletrackselector::unPack<singletrackselector::binning::dca>(o2::aod::singletrackselector::storedDcaXY)) <= max_dcaz && // For now no filtering on the DCAxy or DCAz (casting not supported)
374+
nabs(o2::aod::singletrackselector::unPack<singletrackselector::binning::dca>(o2::aod::singletrackselector::storedDcaXY)) <= max_dcaxy &&
375+
nabs(o2::aod::singletrackselector::unPack<singletrackselector::binning::dca>(o2::aod::singletrackselector::storedDcaXY)) <= max_dcaz &&
354376
nabs(o2::aod::singletrackselector::eta) <= etacut;
355377

356378
template <typename Type>
@@ -536,6 +558,32 @@ struct hadronnucleicorrelation {
536558
} // tracks 1
537559
}
538560

561+
template <int ME, typename Type>
562+
void mixMCParticles(Type const& particles1, Type const& particles2)
563+
{ // last value: 0 -- SE; 1 -- ME
564+
for (auto it1 : particles1) {
565+
for (auto it2 : particles2) {
566+
567+
// Calculate Delta-eta Delta-phi (gen)
568+
float deltaEtaGen = it2->eta() - it1->eta();
569+
float deltaPhiGen = it2->phi() - it1->phi();
570+
deltaPhiGen = getDeltaPhi(deltaPhiGen);
571+
572+
for (int k = 0; k < nBinspT; k++) {
573+
574+
if (it1->pt() >= pTBins.value.at(k) && it1->pt() < pTBins.value.at(k + 1)) {
575+
576+
if (ME) {
577+
hEtaPhiGen_AntiDeAntiPr_ME[k]->Fill(deltaEtaGen, deltaPhiGen, it2->pt());
578+
} else {
579+
hEtaPhiGen_AntiDeAntiPr_SE[k]->Fill(deltaEtaGen, deltaPhiGen, it2->pt());
580+
} // SE
581+
} // pT condition
582+
} // nBinspT loop
583+
} // particles 2
584+
} // particles 1
585+
}
586+
539587
float getDeltaPhi(float deltaPhi)
540588
{
541589
if (deltaPhi < -TMath::Pi() / 2) {
@@ -1236,6 +1284,139 @@ struct hadronnucleicorrelation {
12361284
mixbins_antidantip.clear(); // clear the map
12371285
}
12381286
PROCESS_SWITCH(hadronnucleicorrelation, processMC, "processMC", false);
1287+
1288+
void processGen(SimCollisions const& mcCollisions,
1289+
SimParticles const& mcParticles)
1290+
{
1291+
for (auto particle : mcParticles) {
1292+
1293+
if (isPrim && !particle.isPhysicalPrimary()) {
1294+
continue;
1295+
}
1296+
1297+
if (particle.pdgCode() == pdgDeuteron) {
1298+
selectedparticlesMC_d[particle.mcCollisionId()].push_back(std::make_shared<decltype(particle)>(particle));
1299+
}
1300+
if (particle.pdgCode() == -pdgDeuteron) {
1301+
selectedparticlesMC_antid[particle.mcCollisionId()].push_back(std::make_shared<decltype(particle)>(particle));
1302+
}
1303+
if (particle.pdgCode() == pdgProton) {
1304+
selectedparticlesMC_p[particle.mcCollisionId()].push_back(std::make_shared<decltype(particle)>(particle));
1305+
}
1306+
if (particle.pdgCode() == -pdgProton) {
1307+
selectedparticlesMC_antip[particle.mcCollisionId()].push_back(std::make_shared<decltype(particle)>(particle));
1308+
}
1309+
}
1310+
1311+
for (auto collision : mcCollisions) {
1312+
1313+
registry.fill(HIST("Generated/hNEventsMC"), 0.5);
1314+
if (std::abs(collision.posZ()) < 10.f) {
1315+
// return;
1316+
registry.fill(HIST("Generated/hNEventsMC"), 1.5);
1317+
}
1318+
1319+
int vertexBinToMix = std::floor((collision.posZ() + cutzvertex) / (2 * cutzvertex / _vertexNbinsToMix));
1320+
// int centBinToMix = std::floor(collision.centFT0M() / (100.0 / _multNsubBins));
1321+
1322+
if (selectedparticlesMC_antid.find(collision.globalIndex()) != selectedparticlesMC_antid.end()) {
1323+
mixbinsMC_antidantip[vertexBinToMix].push_back(std::make_shared<decltype(collision)>(collision));
1324+
}
1325+
if (selectedparticlesMC_d.find(collision.globalIndex()) != selectedparticlesMC_d.end()) {
1326+
mixbinsMC_dp[vertexBinToMix].push_back(std::make_shared<decltype(collision)>(collision));
1327+
}
1328+
} // coll
1329+
1330+
if (!mixbinsMC_antidantip.empty()) {
1331+
1332+
for (auto i = mixbinsMC_antidantip.begin(); i != mixbinsMC_antidantip.end(); i++) { // iterating over all vertex&mult bins
1333+
1334+
std::vector<MCcolType> value = i->second;
1335+
int EvPerBin = value.size(); // number of collisions in each vertex&mult bin
1336+
1337+
for (int indx1 = 0; indx1 < EvPerBin; indx1++) { // loop over all the events in each vertex&mult bin
1338+
1339+
auto col1 = value[indx1];
1340+
1341+
if (selectedparticlesMC_antip.find(col1->index()) != selectedparticlesMC_antip.end()) {
1342+
mixMCParticles<0>(selectedparticlesMC_antid[col1->index()], selectedparticlesMC_antip[col1->index()]); // mixing SE
1343+
}
1344+
1345+
for (int indx2 = indx1 + 1; indx2 < EvPerBin; indx2++) { // nested loop for all the combinations of collisions in a chosen mult/vertex bin
1346+
1347+
auto col2 = (i->second)[indx2];
1348+
1349+
if (col1 == col2) {
1350+
continue;
1351+
}
1352+
1353+
if (selectedparticlesMC_antip.find(col2->index()) != selectedparticlesMC_antip.end()) {
1354+
mixMCParticles<1>(selectedparticlesMC_antid[col1->index()], selectedparticlesMC_antip[col2->index()]); // mixing ME
1355+
}
1356+
}
1357+
}
1358+
}
1359+
}
1360+
1361+
if (!mixbinsMC_dp.empty()) {
1362+
1363+
for (auto i = mixbinsMC_dp.begin(); i != mixbinsMC_dp.end(); i++) { // iterating over all vertex&mult bins
1364+
1365+
std::vector<MCcolType> value = i->second;
1366+
int EvPerBin = value.size(); // number of collisions in each vertex&mult bin
1367+
1368+
for (int indx1 = 0; indx1 < EvPerBin; indx1++) { // loop over all the events in each vertex&mult bin
1369+
1370+
auto col1 = value[indx1];
1371+
1372+
if (selectedparticlesMC_p.find(col1->index()) != selectedparticlesMC_p.end()) {
1373+
mixMCParticles<0>(selectedparticlesMC_d[col1->index()], selectedparticlesMC_p[col1->index()]); // mixing SE
1374+
}
1375+
1376+
for (int indx2 = indx1 + 1; indx2 < EvPerBin; indx2++) { // nested loop for all the combinations of collisions in a chosen mult/vertex bin
1377+
1378+
auto col2 = (i->second)[indx2];
1379+
1380+
if (col1 == col2) {
1381+
continue;
1382+
}
1383+
1384+
if (selectedparticlesMC_p.find(col2->index()) != selectedparticlesMC_p.end()) {
1385+
mixMCParticles<1>(selectedparticlesMC_d[col1->index()], selectedparticlesMC_p[col2->index()]); // mixing ME
1386+
}
1387+
}
1388+
}
1389+
}
1390+
}
1391+
1392+
// clearing up
1393+
for (auto i = selectedparticlesMC_antid.begin(); i != selectedparticlesMC_antid.end(); i++)
1394+
(i->second).clear();
1395+
selectedparticlesMC_antid.clear();
1396+
1397+
for (auto i = selectedparticlesMC_d.begin(); i != selectedparticlesMC_d.end(); i++)
1398+
(i->second).clear();
1399+
selectedparticlesMC_d.clear();
1400+
1401+
for (auto i = selectedparticlesMC_antip.begin(); i != selectedparticlesMC_antip.end(); i++)
1402+
(i->second).clear();
1403+
selectedparticlesMC_antip.clear();
1404+
1405+
for (auto i = selectedparticlesMC_p.begin(); i != selectedparticlesMC_p.end(); i++)
1406+
(i->second).clear();
1407+
selectedparticlesMC_p.clear();
1408+
1409+
for (auto& pair : mixbinsMC_antidantip) {
1410+
pair.second.clear(); // clear the vector associated with the key
1411+
}
1412+
mixbinsMC_antidantip.clear(); // clear the map
1413+
1414+
for (auto& pair : mixbinsMC_dp) {
1415+
pair.second.clear(); // clear the vector associated with the key
1416+
}
1417+
mixbinsMC_dp.clear(); // clear the map
1418+
}
1419+
PROCESS_SWITCH(hadronnucleicorrelation, processGen, "processGen", false);
12391420
};
12401421

12411422
WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)

0 commit comments

Comments
 (0)