Skip to content

Commit 503819b

Browse files
omvazqueomvazque
committed
MC closure
1 parent 88678fd commit 503819b

File tree

1 file changed

+91
-64
lines changed

1 file changed

+91
-64
lines changed

‎PWGLF/Tasks/Nuspex/piKpRAA.cxx‎

Lines changed: 91 additions & 64 deletions
Original file line numberDiff line numberDiff line change
@@ -198,7 +198,7 @@ struct PiKpRAA {
198198
Configurable<bool> isCentSel{"isCentSel", true, "Centrality selection?"};
199199
Configurable<bool> isT0Ccent{"isT0Ccent", true, "Use T0C-based centrality?"};
200200
Configurable<bool> isZvtxPosSel{"isZvtxPosSel", true, "Zvtx position selection?"};
201-
Configurable<bool> isINELgt0{"isINELgt0", true, "Apply INEL > 0?"};
201+
Configurable<bool> selINELgt0{"selINELgt0", true, "Select INEL > 0?"};
202202
Configurable<bool> isApplyFT0CbasedOccupancy{"isApplyFT0CbasedOccupancy", false, "T0C Occu cut"};
203203
Configurable<bool> applyNchSel{"applyNchSel", false, "Use mid-rapidity-based Nch selection"};
204204
Configurable<bool> skipRecoColGTOne{"skipRecoColGTOne", true, "Remove collisions if reconstructed more than once"};
@@ -551,11 +551,19 @@ struct PiKpRAA {
551551
registry.add("PtPiVsNchMC_AllGen", "All Generated Events;;Gen. Nch;", kTH2F, {{axisPt, {nBinsNch, minNch, maxNch}}});
552552
registry.add("PtKaVsNchMC_AllGen", "All Generated Events;;Gen. Nch;", kTH2F, {{axisPt, {nBinsNch, minNch, maxNch}}});
553553
registry.add("PtPrVsNchMC_AllGen", "All Generated Events;;Gen. Nch;", kTH2F, {{axisPt, {nBinsNch, minNch, maxNch}}});
554+
555+
registry.add("MCclosure_PtMCPiVsNchMC", "All Generated Events 4 MC closure;;Gen. Nch;", kTH2F, {{axisPt, {nBinsNch, minNch, maxNch}}});
556+
registry.add("MCclosure_PtMCKaVsNchMC", "All Generated Events 4 MC closure;;Gen. Nch;", kTH2F, {{axisPt, {nBinsNch, minNch, maxNch}}});
557+
registry.add("MCclosure_PtMCPrVsNchMC", "All Generated Events 4 MC closure;;Gen. Nch;", kTH2F, {{axisPt, {nBinsNch, minNch, maxNch}}});
558+
559+
registry.add("MCclosure_PtPiVsCent", "Gen Evts With at least one Rec. Coll. + Sel. criteria 4 MC closure;;Gen. Nch;", kTH2F, {axisPt, axisCent});
560+
registry.add("MCclosure_PtKaVsCent", "Gen Evts With at least one Rec. Coll. + Sel. criteria 4 MC closure;;Gen. Nch;", kTH2F, {axisPt, axisCent});
561+
registry.add("MCclosure_PtPrVsCent", "Gen Evts With at least one Rec. Coll. + Sel. criteria 4 MC closure;;Gen. Nch;", kTH2F, {axisPt, axisCent});
554562
}
555563

556564
LOG(info) << "\tccdbNoLaterThan=" << ccdbNoLaterThan.value;
557565
LOG(info) << "\tapplyNchSel=" << applyNchSel.value;
558-
LOG(info) << "\tisINELgt0=" << isINELgt0.value;
566+
LOG(info) << "\tselINELgt0=" << selINELgt0.value;
559567
LOG(info) << "\tdetector4Calibration=" << detector4Calibration.value;
560568
LOG(info) << "\tv0TypeSelection=" << static_cast<int>(v0Selections.v0TypeSelection);
561569
LOG(info) << "\tselElecFromGammas=" << v0Selections.selElecFromGammas;
@@ -1369,45 +1377,49 @@ struct PiKpRAA {
13691377
void processSim(aod::McCollisions::iterator const& mccollision, soa::SmallGroups<ColEvSelsMC> const& collisions, BCsRun3 const& /*bcs*/, aod::FT0s const& /*ft0s*/, aod::McParticles const& mcParticles, TracksMC const& tracksMC)
13701378
{
13711379

1380+
//---------------------------
13721381
// Only INEL > 0 generated collisions
13731382
// By counting number of primary charged particles in |eta| < 1
1374-
if (isINELgt0) {
1375-
int nChMC{0};
1376-
for (const auto& particle : mcParticles) {
1383+
//---------------------------
1384+
int nChMC{0};
1385+
for (const auto& particle : mcParticles) {
13771386

1378-
if (std::abs(particle.eta()) > kOne)
1379-
continue;
1387+
if (std::abs(particle.eta()) > kOne)
1388+
continue;
13801389

1381-
auto charge{0.};
1382-
// Get the MC particle
1383-
const auto* pdgParticle = pdg->GetParticle(particle.pdgCode());
1384-
if (pdgParticle != nullptr) {
1385-
charge = pdgParticle->Charge();
1386-
} else {
1387-
continue;
1388-
}
1390+
auto charge{0.};
1391+
// Get the MC particle
1392+
const auto* pdgParticle = pdg->GetParticle(particle.pdgCode());
1393+
if (pdgParticle != nullptr) {
1394+
charge = pdgParticle->Charge();
1395+
} else {
1396+
continue;
1397+
}
13891398

1390-
// Is it a charged particle?
1391-
if (std::abs(charge) < kMinCharge)
1392-
continue;
1399+
// Is it a charged particle?
1400+
if (std::abs(charge) < kMinCharge)
1401+
continue;
13931402

1394-
// Is it a primary particle?
1395-
if (!particle.isPhysicalPrimary())
1396-
continue;
1403+
// Is it a primary particle?
1404+
if (!particle.isPhysicalPrimary())
1405+
continue;
13971406

1398-
nChMC++;
1399-
}
1407+
nChMC++;
1408+
}
14001409

1401-
// Only INEL > 0 generated events
1410+
//---------------------------
1411+
// Only INEL > 0 generated events
1412+
//---------------------------
1413+
if (selINELgt0) {
14021414
if (!(nChMC > kZeroInt)) {
14031415
return;
14041416
}
1405-
}
1417+
} // selINELgt0 condition: Rejects all NON-INEL > 0
14061418

14071419
const auto& nRecColls{collisions.size()};
1408-
14091420
registry.fill(HIST("NumberOfRecoCollisions"), nRecColls);
14101421

1422+
// Only Generated evets with at least one reconstrued collision
14111423
if (nRecColls > kZeroInt) {
14121424

14131425
// Finds the collisions with the largest number of contributors
@@ -1427,6 +1439,10 @@ struct PiKpRAA {
14271439
}
14281440
}
14291441

1442+
//---------------------------
1443+
// Loop over the reconstructed collisions
1444+
// Only that one with the largest number of contributors is considered
1445+
//---------------------------
14301446
for (const auto& collision : collisions) {
14311447

14321448
// Choose the collisions with the largest number of contributors
@@ -1452,34 +1468,9 @@ struct PiKpRAA {
14521468
return;
14531469
}
14541470

1471+
//---------------------------
14551472
// Needed to construct the correlation between MC Nch v.s. centrality
1456-
int nChMC{0};
1457-
for (const auto& particle : mcParticles) {
1458-
if (particle.eta() < v0Selections.minEtaDaughter || particle.eta() > v0Selections.maxEtaDaughter)
1459-
continue;
1460-
1461-
if (particle.pt() < v0Selections.minPt || particle.pt() > v0Selections.maxPt)
1462-
continue;
1463-
1464-
auto charge{0.};
1465-
// Get the MC particle
1466-
const auto* pdgParticle = pdg->GetParticle(particle.pdgCode());
1467-
if (pdgParticle != nullptr) {
1468-
charge = pdgParticle->Charge();
1469-
} else {
1470-
continue;
1471-
}
1472-
1473-
// Is it a charged particle?
1474-
if (std::abs(charge) < kMinCharge)
1475-
continue;
1476-
1477-
// Is it a primary particle?
1478-
if (!particle.isPhysicalPrimary())
1479-
continue;
1480-
1481-
nChMC++;
1482-
}
1473+
//---------------------------
14831474

14841475
const float centrality{isT0Ccent ? collision.centFT0C() : collision.centFT0M()};
14851476
registry.fill(HIST("Centrality_WRecoEvt"), centrality);
@@ -1488,6 +1479,7 @@ struct PiKpRAA {
14881479
//---------------------------
14891480
// All Generated events with at least one associated reconstructed collision
14901481
// The Generated events are not subjected to any selection criteria
1482+
// This histograms are used for the denominator of the tracking efficiency
14911483
//---------------------------
14921484
for (const auto& particle : mcParticles) {
14931485
if (particle.eta() < v0Selections.minEtaDaughter || particle.eta() > v0Selections.maxEtaDaughter)
@@ -1547,7 +1539,10 @@ struct PiKpRAA {
15471539

15481540
const auto& groupedTracks{tracksMC.sliceBy(perCollision, collision.globalIndex())};
15491541

1542+
//---------------------------
15501543
// Track selection with Open DCAxy
1544+
// This is needed for the Secondary Particle Correction
1545+
//---------------------------
15511546
for (const auto& track : groupedTracks) {
15521547
// Track Selection
15531548
if (track.eta() < v0Selections.minEtaDaughter || track.eta() > v0Selections.maxEtaDaughter)
@@ -1634,7 +1629,11 @@ struct PiKpRAA {
16341629
}
16351630
}
16361631

1632+
//---------------------------
16371633
// Global track + DCAxy selections
1634+
// This is needed for the number of the Tracking Efficiency
1635+
// and the spectra to be corrected
1636+
//---------------------------
16381637
for (const auto& track : groupedTracks) {
16391638
// Track Selection
16401639
if (track.eta() < v0Selections.minEtaDaughter || track.eta() > v0Selections.maxEtaDaughter)
@@ -1716,25 +1715,33 @@ struct PiKpRAA {
17161715
continue;
17171716
}
17181717

1719-
if (isPi && !isKa && !isPr)
1718+
if (isPi && !isKa && !isPr) {
17201719
registry.fill(HIST("PtPiVsCent_WithRecoEvt"), track.pt(), centrality);
1721-
if (isKa && !isPi && !isPr)
1720+
registry.fill(HIST("MCclosure_PtPiVsCent"), track.pt(), centrality);
1721+
}
1722+
if (isKa && !isPi && !isPr) {
17221723
registry.fill(HIST("PtKaVsCent_WithRecoEvt"), track.pt(), centrality);
1723-
if (isPr && !isPi && !isKa)
1724+
registry.fill(HIST("MCclosure_PtKaVsCent"), track.pt(), centrality);
1725+
}
1726+
if (isPr && !isPi && !isKa) {
17241727
registry.fill(HIST("PtPrVsCent_WithRecoEvt"), track.pt(), centrality);
1725-
1728+
registry.fill(HIST("MCclosure_PtPrVsCent"), track.pt(), centrality);
1729+
}
17261730
registry.fill(HIST("PtResolution"), particle.pt(), (track.pt() - particle.pt()) / particle.pt());
17271731
} // Loop over reconstructed tracks
1728-
} // Loop over Reco. Collisions: These collisions are not required to pass the event selection
1732+
} // Loop over Reco. Collisions: Only the collisions with the largest number of contributors
17291733
} // If condition: Only simulated evets with at least one reconstrued collision
17301734

17311735
//---------------------------
17321736
// All Generated events irrespective of whether there is an associated reconstructed collision
17331737
// Consequently, the centrality being a reconstructed quantity, might not always be available
1734-
// Therefore it is expressed as a function of the generated pT and the generated Nch in ∣eta∣ < 0.5
1738+
// Therefore it is expressed as a function of the generated pT and the generated Nch in ∣eta∣ < 0.8
17351739
//---------------------------
17361740

1737-
int nChMC{0};
1741+
//---------------------------
1742+
// To perform MC closure
1743+
// True Pt vs Generated Nch
1744+
//---------------------------
17381745
for (const auto& particle : mcParticles) {
17391746
if (particle.eta() < v0Selections.minEtaDaughter || particle.eta() > v0Selections.maxEtaDaughter)
17401747
continue;
@@ -1756,12 +1763,28 @@ struct PiKpRAA {
17561763
continue;
17571764

17581765
// Is it a primary particle?
1766+
bool isPrimary{true};
17591767
if (!particle.isPhysicalPrimary())
1760-
continue;
1768+
isPrimary = false;
17611769

1762-
nChMC++;
1763-
}
1770+
if (isPrimary) {
1771+
if (particle.pdgCode() == PDG_t::kPiPlus || particle.pdgCode() == PDG_t::kPiMinus) {
1772+
registry.fill(HIST("MCclosure_PtMCPiVsNchMC"), particle.pt(), nChMC);
1773+
} else if (particle.pdgCode() == PDG_t::kKPlus || particle.pdgCode() == PDG_t::kKMinus) {
1774+
registry.fill(HIST("MCclosure_PtMCKaVsNchMC"), particle.pt(), nChMC);
1775+
} else if (particle.pdgCode() == PDG_t::kProton || particle.pdgCode() == PDG_t::kProtonBar) {
1776+
registry.fill(HIST("MCclosure_PtMCPrVsNchMC"), particle.pt(), nChMC);
1777+
} else {
1778+
continue;
1779+
}
1780+
}
1781+
} // Loop over generated particles per Gen Event 4 MC closure
17641782

1783+
//---------------------------
1784+
// Generated Pt spectra of all INEL > 0 Generated evets
1785+
// irrespective of whether there is a reconstructed collision
1786+
// This is used for the denominator of the signal loss correction
1787+
//---------------------------
17651788
for (const auto& particle : mcParticles) {
17661789
if (particle.eta() < v0Selections.minEtaDaughter || particle.eta() > v0Selections.maxEtaDaughter)
17671790
continue;
@@ -1799,6 +1822,10 @@ struct PiKpRAA {
17991822
}
18001823
}
18011824
} // Loop over Generated Particles
1825+
1826+
//---------------------------
1827+
// This is used for the denominator of the event loss correction
1828+
//---------------------------
18021829
registry.fill(HIST("NchMC_AllGen"), nChMC);
18031830
}
18041831
PROCESS_SWITCH(PiKpRAA, processSim, "Process Sim", false);
@@ -2180,7 +2207,7 @@ struct PiKpRAA {
21802207
registry.fill(HIST("EventCounter"), EvCutLabel::VtxZ);
21812208
}
21822209

2183-
if (isINELgt0) {
2210+
if (selINELgt0) {
21842211
if (!col.isInelGt0()) {
21852212
return false;
21862213
}

0 commit comments

Comments
 (0)