Skip to content

Commit 19494c1

Browse files
mapalharesmapalhares
committed
Add process function to compute event loss
1 parent f6e49e8 commit 19494c1

File tree

1 file changed

+192
-15
lines changed

1 file changed

+192
-15
lines changed

‎PWGLF/TableProducer/Nuspex/hyperRecoTask.cxx‎

Lines changed: 192 additions & 15 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)
@@ -345,6 +402,8 @@ struct hyperRecoTask {
345402
return static_cast<float>((candidate.tpcSignal() - expTPCSignal) / resoTPC);
346403
}
347404

405+
406+
348407
template <class Tcoll>
349408
void selectGoodCollisions(const Tcoll& collisions)
350409
{
@@ -366,12 +425,18 @@ struct hyperRecoTask {
366425
}
367426
}
368427

369-
if (!collision.selection_bit(aod::evsel::kIsTriggerTVX) || !collision.selection_bit(aod::evsel::kNoTimeFrameBorder) || std::abs(collision.posZ()) > 10) {
428+
if (!collision.selection_bit(aod::evsel::kIsTriggerTVX) || !collision.selection_bit(aod::evsel::kNoTimeFrameBorder)) {
370429
continue;
371430
}
372431

373432
hEvents->Fill(1.);
374433

434+
if(std::abs(collision.posZ()) > 10){
435+
hEvents->Fill(2.);
436+
continue;
437+
}
438+
439+
375440
if (zorroSelected) {
376441
hEventsZorro->Fill(1.);
377442
}
@@ -380,14 +445,14 @@ struct hyperRecoTask {
380445
if (!collision.selection_bit(aod::evsel::kNoSameBunchPileup)) {
381446
continue;
382447
}
383-
hEvents->Fill(2.);
384-
}
448+
hEvents->Fill(3.);
385449

450+
}
386451
if (cfgEvSelkIsGoodZvtxFT0vsPV) {
387452
if (!collision.selection_bit(aod::evsel::kIsGoodZvtxFT0vsPV)) {
388453
continue;
389454
}
390-
hEvents->Fill(3.);
455+
hEvents->Fill(4.);
391456
}
392457

393458
goodCollision[collision.globalIndex()] = true;
@@ -408,23 +473,27 @@ struct hyperRecoTask {
408473
if (collision.has_mcCollision()) {
409474
recoCollisionIds[collision.mcCollisionId()] = collision.globalIndex();
410475
}
411-
if (!collision.selection_bit(aod::evsel::kIsTriggerTVX) || !collision.selection_bit(aod::evsel::kNoTimeFrameBorder) || std::abs(collision.posZ()) > 10)
476+
if (!collision.selection_bit(aod::evsel::kIsTriggerTVX) || !collision.selection_bit(aod::evsel::kNoTimeFrameBorder))
412477
continue;
413-
478+
414479
hEvents->Fill(1.);
415480

481+
if(std::abs(collision.posZ()) > 10){
482+
hEvents->Fill(2.);
483+
continue;
484+
}
485+
416486
if (cfgEvSelkNoSameBunchPileup) {
417487
if (!collision.selection_bit(aod::evsel::kNoSameBunchPileup)) {
418488
continue;
419489
}
420-
hEvents->Fill(2.);
421-
}
422-
490+
hEvents->Fill(3.);
491+
}
423492
if (cfgEvSelkIsGoodZvtxFT0vsPV) {
424493
if (!collision.selection_bit(aod::evsel::kIsGoodZvtxFT0vsPV)) {
425494
continue;
426495
}
427-
hEvents->Fill(3.);
496+
hEvents->Fill(4.);
428497
}
429498

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

9271104
WorkflowSpec

0 commit comments

Comments
 (0)