Skip to content

Commit 237c89d

Browse files
mapalharesalibuild
mapalhares
and
alibuild
authored
[PWGLF/NuSpEx] Added process function to compute event loss MC (#15373)
Co-authored-by: ALICE Action Bot <alibuild@cern.ch>
1 parent 44df812 commit 237c89d

File tree

1 file changed

+182
-13
lines changed

1 file changed

+182
-13
lines changed

‎PWGLF/TableProducer/Nuspex/hyperRecoTask.cxx‎

Lines changed: 182 additions & 13 deletions
Original file line numberDiff line numberDiff line change
@@ -40,6 +40,8 @@
4040
#include "MathUtils/BetheBlochAleph.h"
4141
#include "ReconstructionDataFormats/Track.h"
4242

43+
#include "Math/Vector4D.h"
44+
4345
#include <algorithm>
4446
#include <array>
4547
#include <memory>
@@ -49,13 +51,18 @@
4951
using namespace o2;
5052
using namespace o2::framework;
5153
using namespace o2::framework::expressions;
54+
using std::array;
55+
5256
using CollBracket = o2::math_utils::Bracket<int>;
5357
using TracksFull = soa::Join<aod::TracksIU, aod::TracksExtra, aod::TracksCovIU, aod::TOFSignal, aod::TOFEvTime>;
5458
using CollisionsFull = soa::Join<aod::Collisions, aod::EvSels, aod::CentFT0As, aod::CentFT0Cs, aod::CentFT0Ms>;
5559
using CollisionsFullMC = soa::Join<aod::Collisions, aod::McCollisionLabels, aod::EvSels, aod::CentFT0As, aod::CentFT0Cs, aod::CentFT0Ms>;
5660

5761
using CollisionsFullWithFlow = soa::Join<aod::Collisions, aod::EvSels, aod::CentFT0As, aod::CentFT0Cs, aod::CentFT0Ms, aod::FT0Mults, aod::FV0Mults, aod::TPCMults, aod::EPCalibrationTables>;
5862

63+
using McCollisionMults = soa::Join<aod::McCollisions, aod::MultMCExtras>;
64+
using EventCandidatesMC = soa::Join<aod::Collisions, aod::EvSels, aod::McCollisionLabels, aod::CentFT0Ms, aod::CentFT0As, aod::CentFT0Cs, aod::CentFV0As, aod::Mults>;
65+
5966
namespace
6067
{
6168
constexpr double betheBlochDefault[1][6]{{-1.e32, -1.e32, -1.e32, -1.e32, -1.e32, -1.e32}};
@@ -77,6 +84,20 @@ std::shared_ptr<TH1> hH4LMassTracked;
7784
std::shared_ptr<TH1> hDecayChannel;
7885
std::shared_ptr<TH1> hIsMatterGen;
7986
std::shared_ptr<TH1> hIsMatterGenTwoBody;
87+
std::shared_ptr<TH1> hEvtMC;
88+
std::shared_ptr<TH1> hImpactParamGen;
89+
std::shared_ptr<TH1> hImpactParamReco;
90+
std::shared_ptr<TH1> hGen3HLBeforeEvtSel;
91+
std::shared_ptr<TH1> hGen3HLAfterSel;
92+
std::shared_ptr<TH2> hRecoCentralityColvsMultiplicityRecoEta05;
93+
std::shared_ptr<TH2> hRecoCentralityColvsImpactParamReco;
94+
std::shared_ptr<TH2> hGen3HLvsImpactParameterBeforeEvtSel;
95+
std::shared_ptr<TH2> hGen3HLvsImpactParameterAfterSel;
96+
std::shared_ptr<TH2> hGen3HLvsMultiplicityGenEta05BeforeEvtSel;
97+
std::shared_ptr<TH2> hGen3HLvsMultiplicityGenEta05AfterSel;
98+
std::shared_ptr<TH2> hGen3HLvsMultiplicityFT0CBeforeEvtSel;
99+
std::shared_ptr<TH2> hGen3HLvsMultiplicityFT0CAfterSel;
100+
80101
} // namespace
81102

82103
struct hyperCandidate {
@@ -147,6 +168,7 @@ struct hyperRecoTask {
147168
// PDG codes
148169
Configurable<int> hyperPdg{"hyperPDG", 1010010030, "PDG code of the hyper-mother (could be 3LamH or 4LamH)"};
149170
Configurable<int> heDauPdg{"heDauPDG", 1000020030, "PDG code of the helium (could be 3He or 4He)"};
171+
Configurable<int> piDauPdg{"piDauPdg", 211, "PDG code of pion"};
150172

151173
// Selection criteria
152174
Configurable<double> v0cospacut{"hypcospa", 0.95, "V0 CosPA"};
@@ -199,6 +221,13 @@ struct hyperRecoTask {
199221
ConfigurableAxis zVtxBins{"zVtxBins", {100, -20.f, 20.f}, "Binning for n sigma"};
200222
ConfigurableAxis centBins{"centBins", {100, 0.f, 100.f}, "Binning for centrality"};
201223

224+
// histogram axes for EvtLossMC
225+
ConfigurableAxis binsImpactPar{"binsImpactPar", {80, 0, 16}, "Binning of the impact parameter axis"};
226+
ConfigurableAxis binsCent{"binsCent", {10, 0.0, 100.0}, "Binning of the centrality axis"};
227+
ConfigurableAxis binsPt{"binsPt", {20, 0, 10}, "Binning of the pt"};
228+
ConfigurableAxis binsFT0CMult{"binsFT0CMult", {500, 0.0f, +500.0f}, "Binning of the FT0C multiplicity"};
229+
ConfigurableAxis binsMult{"binsMult", {500, 0.0f, +500.0f}, ""};
230+
202231
// std vector of candidates
203232
std::vector<hyperCandidate> hyperCandidates;
204233
// vector to keep track of MC mothers already filled
@@ -219,7 +248,6 @@ struct hyperRecoTask {
219248

220249
void init(InitContext const&)
221250
{
222-
223251
zorroSummary.setObject(zorro.getZorroSummary());
224252

225253
mRunNumber = 0;
@@ -250,6 +278,11 @@ struct hyperRecoTask {
250278
const AxisSpec nSigma3HeAxis{nSigmaBins, "n_{#sigma}({}^{3}He)"};
251279
const AxisSpec zVtxAxis{zVtxBins, "z_{vtx} (cm)"};
252280
const AxisSpec centAxis{centBins, "Centrality"};
281+
const AxisSpec impactParamAxis{binsImpactPar, "Impact Parameter (b)"};
282+
const AxisSpec centFT0CAxis{binsCent, "Centrality (FT0C %)"};
283+
const AxisSpec binsFT0CMultAxis{binsFT0CMult, "FT0C multiplicity"};
284+
const AxisSpec ptAxis{binsPt, "#it{p}_{T} (GeV/#it{c})"};
285+
const AxisSpec multAxis = {binsMult, "Multiplicity #eta <0.5"};
253286

254287
hNsigma3HeSel = qaRegistry.add<TH2>("hNsigma3HeSel", "; p_{TPC}/z (GeV/#it{c}); n_{#sigma} ({}^{3}He)", HistType::kTH2F, {rigidityAxis, nSigma3HeAxis});
255288
hDeDx3HeSel = qaRegistry.add<TH2>("hDeDx3HeSel", ";p_{TPC}/z (GeV/#it{c}); dE/dx", HistType::kTH2F, {rigidityAxis, dedxAxis});
@@ -259,11 +292,12 @@ struct hyperRecoTask {
259292
hH4LMassBefSel = qaRegistry.add<TH1>("hH4LMassBefSel", ";M (GeV/#it{c}^{2}); ", HistType::kTH1D, {{60, 3.76, 3.84}});
260293
hH4LMassTracked = qaRegistry.add<TH1>("hH4LMassTracked", ";M (GeV/#it{c}^{2}); ", HistType::kTH1D, {{60, 3.76, 3.84}});
261294

262-
hEvents = qaRegistry.add<TH1>("hEvents", ";Events; ", HistType::kTH1D, {{4, -0.5, 3.5}});
295+
hEvents = qaRegistry.add<TH1>("hEvents", ";Events; ", HistType::kTH1D, {{5, -0.5, 4.5}});
263296
hEvents->GetXaxis()->SetBinLabel(1, "All");
264297
hEvents->GetXaxis()->SetBinLabel(2, "sel8");
265-
hEvents->GetXaxis()->SetBinLabel(3, "kNoSameBunchPileup");
266-
hEvents->GetXaxis()->SetBinLabel(4, "kIsGoodZvtxFT0vsPV");
298+
hEvents->GetXaxis()->SetBinLabel(3, "z_{vtx}");
299+
hEvents->GetXaxis()->SetBinLabel(4, "kNoSameBunchPileup");
300+
hEvents->GetXaxis()->SetBinLabel(5, "kIsGoodZvtxFT0vsPV");
267301

268302
hEventsZorro = qaRegistry.add<TH1>("hEventsZorro", ";Events; ", HistType::kTH1D, {{2, -0.5, 1.5}});
269303
hEventsZorro->GetXaxis()->SetBinLabel(1, "Zorro before evsel");
@@ -284,6 +318,29 @@ struct hyperRecoTask {
284318
hCentFT0A = qaRegistry.add<TH1>("hCentFT0A", ";Centrality; ", HistType::kTH1D, {{100, 0, 100}});
285319
hCentFT0C = qaRegistry.add<TH1>("hCentFT0C", ";Centrality; ", HistType::kTH1D, {{100, 0, 100}});
286320
hCentFT0M = qaRegistry.add<TH1>("hCentFT0M", ";Centrality; ", HistType::kTH1D, {{100, 0, 100}});
321+
322+
if (doprocessEventLossMC) {
323+
hEvtMC = qaRegistry.add<TH1>("QAEvent/hEvtMC", ";; ", HistType::kTH1D, {{3, -0.5, 2.5}});
324+
hEvtMC->GetXaxis()->SetBinLabel(1, "All gen evts");
325+
hEvtMC->GetXaxis()->SetBinLabel(2, "Gen evts with al least one reconstructed");
326+
hEvtMC->GetXaxis()->SetBinLabel(3, "Gen evts with no reconstructed collisions");
327+
// Infomation for all generated collisions collisions
328+
hImpactParamGen = qaRegistry.add<TH1>("QAEvent/McColAll/hImpactParamGen", "Impact parameter of generated MC events; Impact Parameter (b); Counts", HistType::kTH1D, {impactParamAxis});
329+
// Infomation for generated collisions collisions with at least one rec. event
330+
hImpactParamReco = qaRegistry.add<TH1>("QAEvent/McColAll/hImpactParamReco", "Impact parameter of generated MC events with at least one rec. evt; Impact Parameter (b); Counts", HistType::kTH1D, {impactParamAxis});
331+
hRecoCentralityColvsMultiplicityRecoEta05 = qaRegistry.add<TH2>("QAEvent/McColAll/hRecoCentralityColvsMultiplicityRecoEta05", "Correlation between FT0C centrality and charged particle multiplicity in generated MC events with at least one rec. evt; Multiplicity #eta <0.5; Counts", HistType::kTH2D, {centFT0CAxis, multAxis});
332+
hRecoCentralityColvsImpactParamReco = qaRegistry.add<TH2>("QAEvent/McColAll/hRecoCentralityColvsImpactParamReco", "Correlation between FT0C centrality and impact parameter in generated MC events with at least one rec. evt; Impact Parameter (b); Counts", HistType::kTH2D, {centFT0CAxis, impactParamAxis});
333+
// Information of generated 3HL in generated events
334+
hGen3HLBeforeEvtSel = qaRegistry.add<TH1>("QAEvent/McCol3HL/hGen3HLBeforeEvtSel", "3HL generated #it{p}_{T} distribution in all gen evt;#it{p}_{T} (GeV/#it{c}); Counts", HistType::kTH1D, {ptAxis});
335+
hGen3HLvsImpactParameterBeforeEvtSel = qaRegistry.add<TH2>("QAEvent/McCol3HL/hGen3HLvsImpactParameterBeforeEvtSel", "Correlation 3HL generated #it{p}_{T} and impact parameter in all gen evt;#it{p}_{T} (GeV/#it{c}); Impact parameter (b)", HistType::kTH2D, {ptAxis, impactParamAxis});
336+
hGen3HLvsMultiplicityGenEta05BeforeEvtSel = qaRegistry.add<TH2>("QAEvent/McCol3HL/hGen3HLvsMultiplicityGenEta05BeforeEvtSel", "Correlation 3HL generated #it{p}_{T} and charged particle multiplicity in all gen evt;#it{p}_{T} (GeV/#it{c}); Multiplicity #eta <0.5", HistType::kTH2D, {ptAxis, multAxis});
337+
hGen3HLvsMultiplicityFT0CBeforeEvtSel = qaRegistry.add<TH2>("QAEvent/McCol3HL/hGen3HLvsMultiplicityFT0CBeforeEvtSel", "Correlation 3HL generated #it{p}_{T} and FT0C multiplicity in all gen evt;#it{p}_{T} (GeV/#it{c}); FT0C Multiplicity", HistType::kTH2D, {ptAxis, binsFT0CMultAxis});
338+
// Information of generated 3HL in generated events with at least one rec. event
339+
hGen3HLAfterSel = qaRegistry.add<TH1>("QAEvent/McCol3HL/hGen3HLAfterSel", "3HL generated #it{p}_{T} distribution in gen. evts with at least one rec. evt; #it{p}_{T} (GeV/#it{c}); Counts", HistType::kTH1D, {ptAxis});
340+
hGen3HLvsImpactParameterAfterSel = qaRegistry.add<TH2>("QAEvent/McCol3HL/hGen3HLvsImpactParameterAfterSel", "Correlation 3HL generated #it{p}_{T} and impact parameter in gen. evts with at least one rec. evt;#it{p}_{T} (GeV/#it{c}); Impact parameter (b)", HistType::kTH2D, {ptAxis, impactParamAxis});
341+
hGen3HLvsMultiplicityGenEta05AfterSel = qaRegistry.add<TH2>("QAEvent/McCol3HL/hGen3HLvsMultiplicityGenEta05AfterSel", "Correlation 3HL generated #it{p}_{T} and charged particle multiplicity in gen. evts with at least one rec. evt;#it{p}_{T} (GeV/#it{c}); Multiplicity #eta <0.5", HistType::kTH2D, {ptAxis, multAxis});
342+
hGen3HLvsMultiplicityFT0CAfterSel = qaRegistry.add<TH2>("QAEvent/McCol3HL/hGen3HLvsMultiplicityFT0CAfterSel", "Correlation 3HL generated #it{p}_{T} and FT0C multiplicity in gen. evts with at least one rec;#it{p}_{T} (GeV/#it{c}); FT0C Multiplicity", HistType::kTH2D, {ptAxis, binsFT0CMultAxis});
343+
}
287344
}
288345

289346
void initCCDB(aod::BCsWithTimestamps::iterator const& bc)
@@ -366,12 +423,17 @@ struct hyperRecoTask {
366423
}
367424
}
368425

369-
if (!collision.selection_bit(aod::evsel::kIsTriggerTVX) || !collision.selection_bit(aod::evsel::kNoTimeFrameBorder) || std::abs(collision.posZ()) > 10) {
426+
if (!collision.selection_bit(aod::evsel::kIsTriggerTVX) || !collision.selection_bit(aod::evsel::kNoTimeFrameBorder)) {
370427
continue;
371428
}
372429

373430
hEvents->Fill(1.);
374431

432+
if (std::abs(collision.posZ()) > 10) {
433+
hEvents->Fill(2.);
434+
continue;
435+
}
436+
375437
if (zorroSelected) {
376438
hEventsZorro->Fill(1.);
377439
}
@@ -380,14 +442,13 @@ struct hyperRecoTask {
380442
if (!collision.selection_bit(aod::evsel::kNoSameBunchPileup)) {
381443
continue;
382444
}
383-
hEvents->Fill(2.);
445+
hEvents->Fill(3.);
384446
}
385-
386447
if (cfgEvSelkIsGoodZvtxFT0vsPV) {
387448
if (!collision.selection_bit(aod::evsel::kIsGoodZvtxFT0vsPV)) {
388449
continue;
389450
}
390-
hEvents->Fill(3.);
451+
hEvents->Fill(4.);
391452
}
392453

393454
goodCollision[collision.globalIndex()] = true;
@@ -408,23 +469,27 @@ struct hyperRecoTask {
408469
if (collision.has_mcCollision()) {
409470
recoCollisionIds[collision.mcCollisionId()] = collision.globalIndex();
410471
}
411-
if (!collision.selection_bit(aod::evsel::kIsTriggerTVX) || !collision.selection_bit(aod::evsel::kNoTimeFrameBorder) || std::abs(collision.posZ()) > 10)
472+
if (!collision.selection_bit(aod::evsel::kIsTriggerTVX) || !collision.selection_bit(aod::evsel::kNoTimeFrameBorder))
412473
continue;
413474

414475
hEvents->Fill(1.);
415476

477+
if (std::abs(collision.posZ()) > 10) {
478+
hEvents->Fill(2.);
479+
continue;
480+
}
481+
416482
if (cfgEvSelkNoSameBunchPileup) {
417483
if (!collision.selection_bit(aod::evsel::kNoSameBunchPileup)) {
418484
continue;
419485
}
420-
hEvents->Fill(2.);
486+
hEvents->Fill(3.);
421487
}
422-
423488
if (cfgEvSelkIsGoodZvtxFT0vsPV) {
424489
if (!collision.selection_bit(aod::evsel::kIsGoodZvtxFT0vsPV)) {
425490
continue;
426491
}
427-
hEvents->Fill(3.);
492+
hEvents->Fill(4.);
428493
}
429494

430495
if (collision.has_mcCollision()) {
@@ -437,7 +502,6 @@ struct hyperRecoTask {
437502
hCentFT0M->Fill(collision.centFT0M());
438503
}
439504
}
440-
441505
template <class Ttrack, class Tcolls>
442506
void fillHyperCand(Ttrack& heTrack, Ttrack& piTrack, CollBracket collBracket, const Tcolls& collisions, hyperCandidate& hypCand)
443507
{
@@ -922,6 +986,111 @@ struct hyperRecoTask {
922986
processMC(collisions, mcCollisions, V0s, tracks, ambiTracks, bcs, trackLabelsMC, particlesMC);
923987
}
924988
PROCESS_SWITCH(hyperRecoTask, processMCTracked, "MC analysis with tracked V0s", false);
989+
990+
template <typename CollType>
991+
bool passEvtSel(const CollType& collision)
992+
{
993+
if (!collision.sel8())
994+
return false;
995+
996+
if ((std::abs(collision.posZ())) > 10)
997+
return false;
998+
999+
if (cfgEvSelkNoSameBunchPileup && !collision.selection_bit(aod::evsel::kNoSameBunchPileup))
1000+
return false;
1001+
1002+
if (cfgEvSelkIsGoodZvtxFT0vsPV && !collision.selection_bit(o2::aod::evsel::kIsGoodZvtxFT0vsPV))
1003+
return false;
1004+
1005+
return true;
1006+
}
1007+
1008+
void processEventLossMC(McCollisionMults::iterator const& mcCollision, soa::SmallGroups<EventCandidatesMC> const& collisions, aod::McParticles const& GenParticles)
1009+
{
1010+
if (std::abs(mcCollision.posZ()) > 10) {
1011+
return;
1012+
}
1013+
1014+
//////////// Event loss estimation via impact parameter and multiplicity by MCFT0C
1015+
1016+
// Fill all generated events
1017+
hEvtMC->Fill(0);
1018+
hImpactParamGen->Fill(mcCollision.impactParameter());
1019+
1020+
// Fill generated events with no reconstructed collisions
1021+
if (collisions.size() == 0) {
1022+
hEvtMC->Fill(1);
1023+
}
1024+
1025+
// Define the generated events with at least one reconstructed event
1026+
bool atLeastOneRecoEvt = false;
1027+
auto centralityFT0C = -999.;
1028+
1029+
for (auto const& col : collisions) {
1030+
if (!passEvtSel(col)) {
1031+
continue;
1032+
}
1033+
centralityFT0C = col.centFT0C();
1034+
atLeastOneRecoEvt = true;
1035+
}
1036+
1037+
if (atLeastOneRecoEvt) {
1038+
hEvtMC->Fill(2);
1039+
hImpactParamReco->Fill(mcCollision.impactParameter());
1040+
hRecoCentralityColvsMultiplicityRecoEta05->Fill(centralityFT0C, mcCollision.multMCNParticlesEta05());
1041+
hRecoCentralityColvsImpactParamReco->Fill(centralityFT0C, mcCollision.impactParameter());
1042+
}
1043+
// Construct the H3L 4-vector based on the generated daugthers identification by PDG
1044+
ROOT::Math::PxPyPzMVector daugh1, daugh2, mother;
1045+
1046+
for (const auto& genParticle : GenParticles) {
1047+
if (std::abs(genParticle.y()) > 1)
1048+
continue;
1049+
if (std::abs(genParticle.pdgCode()) != hyperPdg)
1050+
continue;
1051+
1052+
auto daughters = genParticle.daughters_as<aod::McParticles>();
1053+
1054+
bool dauHe3 = false, dauPiMinus = false, dauAntiHe3 = false, dauPiPos = true;
1055+
1056+
for (const auto& daughter : daughters) {
1057+
if (daughter.pdgCode() == heDauPdg) {
1058+
dauHe3 = true;
1059+
daugh1 = ROOT::Math::PxPyPzMVector(daughter.px(), daughter.py(), daughter.pz(), he3Mass);
1060+
} else if (daughter.pdgCode() == -piDauPdg) {
1061+
dauPiMinus = true;
1062+
daugh2 = ROOT::Math::PxPyPzMVector(daughter.px(), daughter.py(), daughter.pz(), piMass);
1063+
}
1064+
if (daughter.pdgCode() == -heDauPdg) {
1065+
dauAntiHe3 = true;
1066+
daugh1 = ROOT::Math::PxPyPzMVector(daughter.px(), daughter.py(), daughter.pz(), he3Mass);
1067+
} else if (daughter.pdgCode() == piDauPdg) {
1068+
dauPiPos = true;
1069+
daugh2 = ROOT::Math::PxPyPzMVector(daughter.px(), daughter.py(), daughter.pz(), piMass);
1070+
}
1071+
}
1072+
// Check pairs to avoid wrong charge associations
1073+
if (!((dauHe3 && dauPiMinus) || !(dauAntiHe3 && dauPiPos)))
1074+
continue;
1075+
1076+
mother = daugh1 + daugh2;
1077+
1078+
// Fill informations for generated 3HL in all generated events
1079+
hGen3HLBeforeEvtSel->Fill(mother.pt());
1080+
hGen3HLvsImpactParameterBeforeEvtSel->Fill(mother.pt(), mcCollision.impactParameter());
1081+
hGen3HLvsMultiplicityGenEta05BeforeEvtSel->Fill(mother.pt(), mcCollision.multMCNParticlesEta05());
1082+
hGen3HLvsMultiplicityFT0CBeforeEvtSel->Fill(mother.pt(), mcCollision.multMCFT0C());
1083+
1084+
// Fill informations for generated 3HL in generated events with at least one reconstructed event
1085+
if (atLeastOneRecoEvt) {
1086+
hGen3HLAfterSel->Fill(mother.pt());
1087+
hGen3HLvsImpactParameterAfterSel->Fill(mother.pt(), mcCollision.impactParameter());
1088+
hGen3HLvsMultiplicityGenEta05AfterSel->Fill(mother.pt(), mcCollision.multMCNParticlesEta05());
1089+
hGen3HLvsMultiplicityFT0CAfterSel->Fill(mother.pt(), mcCollision.multMCFT0C());
1090+
}
1091+
}
1092+
}
1093+
PROCESS_SWITCH(hyperRecoTask, processEventLossMC, "Event loss analysis", false);
9251094
};
9261095

9271096
WorkflowSpec

0 commit comments

Comments
 (0)