Skip to content

Commit f432977

Browse files
ercolessiercolessi
authored
[PWGLF] [WIP] [PWGLF] nuclei-proton correlation: add pp correlation in mc gen (#10591)
1 parent 402cb1c commit f432977

File tree

1 file changed

+126
-19
lines changed

1 file changed

+126
-19
lines changed

‎PWGLF/Tasks/Nuspex/hadronnucleicorrelation.cxx‎

Lines changed: 126 additions & 19 deletions
Original file line numberDiff line numberDiff line change
@@ -162,6 +162,9 @@ struct hadronnucleicorrelation {
162162
std::vector<std::shared_ptr<TH3>> hPIDEtaPhiRec_AntiDeAntiPr_ME;
163163
std::vector<std::shared_ptr<TH3>> hPIDEtaPhiGen_AntiDeAntiPr_ME;
164164

165+
std::vector<std::shared_ptr<TH3>> hEtaPhiGen_AntiPrAntiPr_SE;
166+
std::vector<std::shared_ptr<TH3>> hEtaPhiGen_AntiPrAntiPr_ME;
167+
165168
int nBinspT;
166169
TH2F* hEffpTEta_proton;
167170
TH2F* hEffpTEta_antiproton;
@@ -238,6 +241,13 @@ struct hadronnucleicorrelation {
238241
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),
239242
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}});
240243
hEtaPhiGen_AntiDeAntiPr_ME.push_back(std::move(htempMEGen_AntiDeAntiPr));
244+
245+
auto htempSEGen_AntiPrAntiPr = registry.add<TH3>(Form("hEtaPhiGen_AntiPrAntiPr_SE_pt%02.0f%02.0f", pTBins.value.at(i) * 10, pTBins.value.at(i + 1) * 10),
246+
Form("Gen #Delta#eta#Delta#phi (%.1f<p_{T} #bar{p} <%.1f GeV/c)", pTBins.value.at(i), pTBins.value.at(i + 1)), {HistType::kTH3F, {DeltaEtaAxis, DeltaPhiAxis, ptBinnedAxis}});
247+
hEtaPhiGen_AntiPrAntiPr_SE.push_back(std::move(htempSEGen_AntiPrAntiPr));
248+
auto htempMEGen_AntiPrAntiPr = registry.add<TH3>(Form("hEtaPhiGen_AntiPrAntiPr_ME_pt%02.0f%02.0f", pTBins.value.at(i) * 10, pTBins.value.at(i + 1) * 10),
249+
Form("Gen #Delta#eta#Delta#phi (%.1f<p_{T} #bar{p} <%.1f GeV/c)", pTBins.value.at(i), pTBins.value.at(i + 1)), {HistType::kTH3F, {DeltaEtaAxis, DeltaPhiAxis, ptBinnedAxis}});
250+
hEtaPhiGen_AntiPrAntiPr_ME.push_back(std::move(htempMEGen_AntiPrAntiPr));
241251
}
242252
}
243253

@@ -363,6 +373,18 @@ struct hadronnucleicorrelation {
363373
if (isMCGen) {
364374
registry.add("Generated/hNEventsMC", "hNEventsMC", {HistType::kTH1D, {{1, 0.f, 1.f}}});
365375
registry.get<TH1>(HIST("Generated/hNEventsMC"))->GetXaxis()->SetBinLabel(1, "All");
376+
377+
registry.add("Generated/hQAProtons", "hQAProtons", {HistType::kTH1D, {{5, 0.f, 5.f}}});
378+
registry.get<TH1>(HIST("Generated/hQAProtons"))->GetXaxis()->SetBinLabel(1, "All");
379+
registry.get<TH1>(HIST("Generated/hQAProtons"))->GetXaxis()->SetBinLabel(2, "PhysicalPrimary");
380+
registry.get<TH1>(HIST("Generated/hQAProtons"))->GetXaxis()->SetBinLabel(3, "|#eta|<0.8");
381+
registry.get<TH1>(HIST("Generated/hQAProtons"))->GetXaxis()->SetBinLabel(4, "no daughters");
382+
registry.get<TH1>(HIST("Generated/hQAProtons"))->GetXaxis()->SetBinLabel(5, "d daughter");
383+
384+
registry.add("Generated/hQADeuterons", "hQADeuterons", {HistType::kTH1D, {{3, 0.f, 3.f}}});
385+
registry.get<TH1>(HIST("Generated/hQADeuterons"))->GetXaxis()->SetBinLabel(1, "All");
386+
registry.get<TH1>(HIST("Generated/hQADeuterons"))->GetXaxis()->SetBinLabel(2, "PhysicalPrimary");
387+
registry.get<TH1>(HIST("Generated/hQADeuterons"))->GetXaxis()->SetBinLabel(3, "|#eta|<0.8");
366388
}
367389
}
368390

@@ -583,6 +605,34 @@ struct hadronnucleicorrelation {
583605
}
584606
}
585607

608+
template <int ME, typename Type>
609+
void mixMCParticlesIdentical(Type const& particles1, Type const& particles2)
610+
{
611+
for (auto it1 : particles1) {
612+
for (auto it2 : particles2) {
613+
// Calculate Delta-eta Delta-phi (gen)
614+
float deltaEtaGen = it2->eta() - it1->eta();
615+
float deltaPhiGen = getDeltaPhi(it2->phi() - it1->phi());
616+
617+
if (!ME && std::abs(deltaPhiGen) < 0.001 && std::abs(deltaEtaGen) < 0.001) {
618+
continue;
619+
}
620+
621+
// Loop over pT bins
622+
for (int k = 0; k < nBinspT; k++) {
623+
if (it1->pt() >= pTBins.value.at(k) && it1->pt() < pTBins.value.at(k + 1)) {
624+
// Use correct histogram based on ME flag
625+
if constexpr (ME) {
626+
hEtaPhiGen_AntiPrAntiPr_ME[k]->Fill(deltaEtaGen, deltaPhiGen, it2->pt());
627+
} else {
628+
hEtaPhiGen_AntiPrAntiPr_SE[k]->Fill(deltaEtaGen, deltaPhiGen, it2->pt());
629+
}
630+
}
631+
}
632+
}
633+
}
634+
}
635+
586636
float getDeltaPhi(float deltaPhi)
587637
{
588638
if (deltaPhi < -TMath::Pi() / 2) {
@@ -1289,13 +1339,32 @@ struct hadronnucleicorrelation {
12891339
{
12901340
for (auto particle : mcParticles) {
12911341

1342+
if (particle.pdgCode() == pdgProton) {
1343+
registry.fill(HIST("Generated/hQAProtons"), 0.5);
1344+
}
1345+
if (particle.pdgCode() == pdgDeuteron) {
1346+
registry.fill(HIST("Generated/hQADeuterons"), 0.5);
1347+
}
1348+
12921349
if (isPrim && !particle.isPhysicalPrimary()) {
12931350
continue;
12941351
}
1352+
if (particle.pdgCode() == pdgProton) {
1353+
registry.fill(HIST("Generated/hQAProtons"), 1.5);
1354+
}
1355+
if (particle.pdgCode() == pdgDeuteron) {
1356+
registry.fill(HIST("Generated/hQADeuterons"), 1.5);
1357+
}
12951358

12961359
if (std::abs(particle.eta()) > etacut) {
12971360
continue;
12981361
}
1362+
if (particle.pdgCode() == pdgProton) {
1363+
registry.fill(HIST("Generated/hQAProtons"), 2.5);
1364+
}
1365+
if (particle.pdgCode() == pdgDeuteron) {
1366+
registry.fill(HIST("Generated/hQADeuterons"), 2.5);
1367+
}
12991368

13001369
if (particle.pdgCode() == pdgDeuteron) {
13011370
selectedparticlesMC_d[particle.mcCollisionId()].push_back(std::make_shared<decltype(particle)>(particle));
@@ -1306,35 +1375,22 @@ struct hadronnucleicorrelation {
13061375
if (particle.pdgCode() == pdgProton) {
13071376
if (!particle.has_daughters()) {
13081377
selectedparticlesMC_p[particle.mcCollisionId()].push_back(std::make_shared<decltype(particle)>(particle));
1378+
registry.fill(HIST("Generated/hQAProtons"), 3.5);
13091379
} else {
1310-
bool isp = false;
1311-
1380+
bool isd = false;
13121381
for (auto& dau : particle.daughters_as<aod::McParticles>()) {
1313-
if (dau.pdgCode() == pdgProton) {
1314-
isp = true;
1382+
if (dau.pdgCode() == pdgDeuteron) {
1383+
isd = true;
13151384
}
13161385
}
1317-
1318-
if (isp) {
1319-
selectedparticlesMC_p[particle.mcCollisionId()].push_back(std::make_shared<decltype(particle)>(particle));
1386+
if (isd) {
1387+
registry.fill(HIST("Generated/hQAProtons"), 4.5);
13201388
}
13211389
}
13221390
}
13231391
if (particle.pdgCode() == -pdgProton) {
13241392
if (!particle.has_daughters()) {
13251393
selectedparticlesMC_antip[particle.mcCollisionId()].push_back(std::make_shared<decltype(particle)>(particle));
1326-
} else {
1327-
bool isantip = false;
1328-
1329-
for (auto& dau : particle.daughters_as<aod::McParticles>()) {
1330-
if (dau.pdgCode() == -pdgProton) {
1331-
isantip = true;
1332-
}
1333-
}
1334-
1335-
if (isantip) {
1336-
selectedparticlesMC_antip[particle.mcCollisionId()].push_back(std::make_shared<decltype(particle)>(particle));
1337-
}
13381394
}
13391395
}
13401396
}
@@ -1396,6 +1452,57 @@ struct hadronnucleicorrelation {
13961452
}
13971453
}
13981454
}
1455+
1456+
// anti-p - anti-p correlation
1457+
if (selectedparticlesMC_antip.find(collision1.globalIndex()) != selectedparticlesMC_antip.end()) {
1458+
1459+
mixMCParticlesIdentical<0>(selectedparticlesMC_antip[collision1.globalIndex()], selectedparticlesMC_antip[collision1.globalIndex()]); // mixing SE
1460+
1461+
int stop3 = 0;
1462+
1463+
for (auto collision2 : mcCollisions) { // nested loop on collisions
1464+
1465+
if (collision1.globalIndex() == collision2.globalIndex()) {
1466+
continue;
1467+
}
1468+
1469+
if (stop3 > maxmixcollsGen) {
1470+
break;
1471+
}
1472+
1473+
if (selectedparticlesMC_antip.find(collision2.globalIndex()) != selectedparticlesMC_antip.end()) {
1474+
mixMCParticlesIdentical<1>(selectedparticlesMC_antip[collision1.globalIndex()], selectedparticlesMC_antip[collision2.globalIndex()]); // mixing ME
1475+
}
1476+
1477+
stop3++;
1478+
}
1479+
}
1480+
// p - p correlation
1481+
if (domatterGen) {
1482+
if (selectedparticlesMC_p.find(collision1.globalIndex()) != selectedparticlesMC_p.end()) {
1483+
1484+
mixMCParticlesIdentical<0>(selectedparticlesMC_p[collision1.globalIndex()], selectedparticlesMC_p[collision1.globalIndex()]); // mixing SE
1485+
1486+
int stop4 = 0;
1487+
1488+
for (auto collision2 : mcCollisions) { // nested loop on collisions
1489+
1490+
if (collision1.globalIndex() == collision2.globalIndex()) {
1491+
continue;
1492+
}
1493+
1494+
if (stop4 > maxmixcollsGen) {
1495+
break;
1496+
}
1497+
1498+
if (selectedparticlesMC_p.find(collision2.globalIndex()) != selectedparticlesMC_p.end()) {
1499+
mixMCParticlesIdentical<1>(selectedparticlesMC_p[collision1.globalIndex()], selectedparticlesMC_p[collision2.globalIndex()]); // mixing ME
1500+
}
1501+
1502+
stop4++;
1503+
}
1504+
}
1505+
}
13991506
}
14001507

14011508
// clearing up

0 commit comments

Comments
 (0)