Skip to content

Commit eba2097

Browse files
ariedel-cernariedel-cern
committed
Feat: update pair/triplet efficiency task
1 parent b629bc8 commit eba2097

File tree

1 file changed

+232
-25
lines changed

1 file changed

+232
-25
lines changed

‎PWGCF/FemtoDream/Tasks/femtoDreamPairEfficiency.cxx‎

Lines changed: 232 additions & 25 deletions
Original file line numberDiff line numberDiff line change
@@ -83,13 +83,20 @@ struct FemtoDreamPairEfficiency {
8383
SliceCache cache;
8484
Preslice<McParticles> perColMc = mcparticle::mcCollisionId;
8585

86+
struct : ConfigurableGroup {
87+
std::string prefix = std::string("Mode");
88+
Configurable<bool> countPairs{"countPairs", true, "Count Pairs"};
89+
Configurable<bool> countTriplets{"countTriplets", false, "Count Triplets"};
90+
} mode;
91+
8692
// Event Selections
8793
struct : ConfigurableGroup {
8894
std::string prefix = std::string("EventSelection");
8995
Configurable<float> zvtxAbsMax{"zvtxAbsMax", 10.f, "|z-Vertex| max"};
9096
Configurable<bool> offlineCheck{"offlineCheck", true, "Check for Sel8"};
9197
Configurable<float> etaAbsMax{"etaAbsMax", 0.8f, "Common eta cut for particles in dNdEta distribution"};
9298
Configurable<float> kstarMax{"kstarMax", 999.f, "Cut on kstar"};
99+
Configurable<float> q3Max{"q3Max", 999.f, "Cut on Q3"};
93100
} eventSelection;
94101

95102
struct : ConfigurableGroup {
@@ -146,13 +153,39 @@ struct FemtoDreamPairEfficiency {
146153
Configurable<float> tpctofNsigmaMax{"tpctofNsigmaMax", 3.f, "TPCTOC nsigma max"};
147154
} trackCuts2;
148155

156+
struct : ConfigurableGroup {
157+
std::string prefix = std::string("SelectionTrack3");
158+
Configurable<int> sign{"sign", 1, "Sign of charge"};
159+
Configurable<float> ptMin{"ptMin", 0.0f, "pt min"};
160+
Configurable<float> ptMax{"ptMax", 3.0f, "pt max"};
161+
Configurable<float> etaAbsMax{"etaAbsMax", 0.8f, "|eta| max"};
162+
Configurable<float> dcazAbsMax{"dcazAbsMax", 0.1f, "|dca_z| max"};
163+
Configurable<bool> useDcaxyPtDepCut{"useDcaxyPtDepCut", true, "|dca_z| max"};
164+
Configurable<float> tpcClusterMin{"tpcClusterMin", 80.f, "TPC clusters min"};
165+
Configurable<float> tpcCrossedOverClusterMin{"tpcCrossedOverClusterMin", 0.83f, "TPC clusters/TPC crossed rows min"};
166+
Configurable<float> tpcCrossedMin{"tpcCrossedMin", 70.f, "TPC crossed rows min"};
167+
Configurable<float> tpcSharedMax{"tpcSharedMax", 160.f, "TPC shared clusters max"};
168+
Configurable<float> itsClusterMin{"itsClusterMin", 0.f, "ITS clusters min"};
169+
Configurable<float> itsIbClusterMin{"itsIbClusterMin", 0.f, "ITS inner barrle min"};
170+
Configurable<int> pdgCode{"pdgCode", 2212, "PDG code"};
171+
Configurable<float> pidThreshold{"pidThreshold", 0.75f, "Momentum threshold for PID"};
172+
Configurable<float> itsNsigmaMax{"itsNsigmaMax", 99.f, "its nsigma max"};
173+
Configurable<float> tpcNsigmaMax{"tpcNsigmaMax", 3.f, "TPC nsigma max"};
174+
Configurable<float> tpctofNsigmaMax{"tpctofNsigmaMax", 3.f, "TPCTOC nsigma max"};
175+
} trackCuts3;
176+
149177
HistogramRegistry dataEventHist{"dataEventHist", {}, OutputObjHandlingPolicy::AnalysisObject, false, true};
150178
HistogramRegistry mcEventHist{"mcEventHist", {}, OutputObjHandlingPolicy::AnalysisObject, false, true};
151179

152180
HistogramRegistry dataTrackHist{"dataTrackHist", {}, OutputObjHandlingPolicy::AnalysisObject, false, true};
153181

154182
void init(InitContext&)
155183
{
184+
185+
if (mode.countPairs.value && mode.countTriplets.value) {
186+
LOG(fatal) << "We cannot count pairs and triplets at the same time. Turn one of the off.";
187+
}
188+
156189
dataEventHist.add("hDataEventSelection", "hDataEventSelection", kTH1F, {{6, -0.5f, 5.5f}});
157190
dataEventHist.get<TH1>(HIST("hDataEventSelection"))->GetXaxis()->SetBinLabel(1, "All collisions");
158191
dataEventHist.get<TH1>(HIST("hDataEventSelection"))->GetXaxis()->SetBinLabel(2, "sel8 cut");
@@ -193,21 +226,35 @@ struct FemtoDreamPairEfficiency {
193226
dataTrackHist.add("Track2/tpcNsigma", "Track 2 nsigma tpc", kTH2F, {{600, 0, 6}, {1000, -5, 5}});
194227
dataTrackHist.add("Track2/tpctofNsigma", "Track 2 nsigma tpctof", kTH2F, {{600, 0, 6}, {500, 0, 5}});
195228

229+
if (mode.countTriplets) {
230+
dataTrackHist.add("Track3/pt", "Track 3 pt", kTH1F, {{600, 0, 6}});
231+
dataTrackHist.add("Track3/eta", "Track 3 eta", kTH1F, {{200, -1, 1}});
232+
dataTrackHist.add("Track3/phi", "Track 3 phi", kTH1F, {{720, 0, o2::constants::math::TwoPI}});
233+
dataTrackHist.add("Track3/tpcCluster", "Track 3 cluster", kTH1F, {{160, 0, 160}});
234+
dataTrackHist.add("Track3/tpcCrossed", "Track 3 crossed rows", kTH1F, {{160, 0, 160}});
235+
dataTrackHist.add("Track3/tpcShared", "Track 3 cluster shared", kTH1F, {{160, 0, 160}});
236+
dataTrackHist.add("Track3/itsCluster", "Track 3 its cluster", kTH1F, {{0, 0, 7}});
237+
dataTrackHist.add("Track3/itsIbCluster", "Track 3 its cluster inner barrel", kTH1F, {{3, 0, 3}});
238+
dataTrackHist.add("Track3/itsNsigma", "Track 3 nsigma its", kTH2F, {{600, 0, 6}, {1000, -5, 5}});
239+
dataTrackHist.add("Track3/tpcNsigma", "Track 3 nsigma tpc", kTH2F, {{600, 0, 6}, {1000, -5, 5}});
240+
dataTrackHist.add("Track3/tpctofNsigma", "Track 3 nsigma tpctof", kTH2F, {{600, 0, 6}, {500, 0, 5}});
241+
}
242+
196243
mcEventHist.add("hRecoEventSelection", "hRecoEventSelection", kTH1F, {{7, -0.5f, 6.5f}});
197244
mcEventHist.get<TH1>(HIST("hRecoEventSelection"))->GetXaxis()->SetBinLabel(1, "All collisions");
198245
mcEventHist.get<TH1>(HIST("hRecoEventSelection"))->GetXaxis()->SetBinLabel(2, "Sel8 cut");
199246
mcEventHist.get<TH1>(HIST("hRecoEventSelection"))->GetXaxis()->SetBinLabel(3, "posZ cut");
200247
mcEventHist.get<TH1>(HIST("hRecoEventSelection"))->GetXaxis()->SetBinLabel(4, "INEL>0 cut");
201248
mcEventHist.get<TH1>(HIST("hRecoEventSelection"))->GetXaxis()->SetBinLabel(5, "With at least a gen coll");
202-
mcEventHist.get<TH1>(HIST("hRecoEventSelection"))->GetXaxis()->SetBinLabel(6, "With at least a pair");
203-
mcEventHist.get<TH1>(HIST("hRecoEventSelection"))->GetXaxis()->SetBinLabel(7, "With at least a lowkstar pair");
249+
mcEventHist.get<TH1>(HIST("hRecoEventSelection"))->GetXaxis()->SetBinLabel(6, "With at least a pair/triplet");
250+
mcEventHist.get<TH1>(HIST("hRecoEventSelection"))->GetXaxis()->SetBinLabel(7, "With at least a lowkstar pair/lowq3 triplet");
204251

205252
mcEventHist.add("hGenEventSelection", "hGenEventSelection", kTH1F, {{6, -0.5f, 5.5f}});
206253
mcEventHist.get<TH1>(HIST("hGenEventSelection"))->GetXaxis()->SetBinLabel(1, "All collisions");
207254
mcEventHist.get<TH1>(HIST("hGenEventSelection"))->GetXaxis()->SetBinLabel(2, "posZ cut");
208255
mcEventHist.get<TH1>(HIST("hGenEventSelection"))->GetXaxis()->SetBinLabel(3, "INEL>0 cut");
209-
mcEventHist.get<TH1>(HIST("hGenEventSelection"))->GetXaxis()->SetBinLabel(4, "With at least a pair");
210-
mcEventHist.get<TH1>(HIST("hGenEventSelection"))->GetXaxis()->SetBinLabel(5, "With at least a lowkstar pair");
256+
mcEventHist.get<TH1>(HIST("hGenEventSelection"))->GetXaxis()->SetBinLabel(4, "With at least a pair/triplet");
257+
mcEventHist.get<TH1>(HIST("hGenEventSelection"))->GetXaxis()->SetBinLabel(5, "With at least a lowkstar pair/lowq3 triplet");
211258
mcEventHist.get<TH1>(HIST("hGenEventSelection"))->GetXaxis()->SetBinLabel(6, "With at least a reco coll");
212259

213260
// MC Event information for Rec and Gen
@@ -236,6 +283,9 @@ struct FemtoDreamPairEfficiency {
236283
int pairsFound = 0;
237284
std::vector<float> vecKstar{};
238285

286+
int tripletsFound = 0;
287+
std::vector<float> vecq3{};
288+
239289
template <typename T1, typename T2>
240290
bool checkRecoTrackSelections(const T1& track, const T2& sel)
241291
{
@@ -417,6 +467,72 @@ struct FemtoDreamPairEfficiency {
417467
return pairs;
418468
}
419469

470+
template <typename T>
471+
int countRecoTriplets(const T& tracks)
472+
{
473+
int triplets = 0;
474+
for (auto track1 = tracks.begin(); track1 != tracks.end(); track1++) {
475+
if (!checkRecoTrackSelections(track1, trackCuts1) || !checkRecoTrackPidSelections(track1, trackCuts1)) {
476+
continue;
477+
}
478+
for (auto track2 = track1 + 1; track2 != tracks.end(); track2++) {
479+
if (!checkRecoTrackSelections(track2, trackCuts2) || !checkRecoTrackPidSelections(track2, trackCuts2)) {
480+
continue;
481+
}
482+
483+
for (auto track3 = track2 + 1; track3 != tracks.end(); track3++) {
484+
if (!checkRecoTrackSelections(track3, trackCuts3) || !checkRecoTrackPidSelections(track3, trackCuts3)) {
485+
continue;
486+
}
487+
488+
auto nsigma1 = getNSigmaValues(track1, trackCuts1.pdgCode.value);
489+
dataTrackHist.fill(HIST("Track1/pt"), track1.pt());
490+
dataTrackHist.fill(HIST("Track1/eta"), track1.eta());
491+
dataTrackHist.fill(HIST("Track1/phi"), track1.phi());
492+
dataTrackHist.fill(HIST("Track1/tpcCluster"), track1.tpcNClsFound());
493+
dataTrackHist.fill(HIST("Track1/tpcCrossed"), track1.tpcNClsCrossedRows());
494+
dataTrackHist.fill(HIST("Track1/tpcShared"), track1.tpcNClsShared());
495+
dataTrackHist.fill(HIST("Track1/itsCluster"), track1.itsNCls());
496+
dataTrackHist.fill(HIST("Track1/itsIbCluster"), track1.itsNClsInnerBarrel());
497+
dataTrackHist.fill(HIST("Track1/itsNsigma"), track1.p(), nsigma1[0]);
498+
dataTrackHist.fill(HIST("Track1/tpcNsigma"), track1.p(), nsigma1[1]);
499+
dataTrackHist.fill(HIST("Track1/tpctofNsigma"), track1.p(), std::hypot(nsigma1[1], nsigma1[2]));
500+
501+
auto nsigma2 = getNSigmaValues(track2, trackCuts2.pdgCode.value);
502+
dataTrackHist.fill(HIST("Track2/pt"), track2.pt());
503+
dataTrackHist.fill(HIST("Track2/eta"), track2.eta());
504+
dataTrackHist.fill(HIST("Track2/phi"), track2.phi());
505+
dataTrackHist.fill(HIST("Track2/tpcCluster"), track2.tpcNClsFound());
506+
dataTrackHist.fill(HIST("Track2/tpcCrossed"), track2.tpcNClsCrossedRows());
507+
dataTrackHist.fill(HIST("Track2/tpcShared"), track2.tpcNClsShared());
508+
dataTrackHist.fill(HIST("Track2/itsCluster"), track2.itsNCls());
509+
dataTrackHist.fill(HIST("Track2/itsIbCluster"), track2.itsNClsInnerBarrel());
510+
dataTrackHist.fill(HIST("Track2/itsNsigma"), track2.p(), nsigma2[0]);
511+
dataTrackHist.fill(HIST("Track2/tpcNsigma"), track2.p(), nsigma2[1]);
512+
dataTrackHist.fill(HIST("Track2/tpctofNsigma"), track2.p(), std::hypot(nsigma2[1], nsigma2[2]));
513+
514+
auto nsigma3 = getNSigmaValues(track3, trackCuts3.pdgCode.value);
515+
dataTrackHist.fill(HIST("Track3/pt"), track3.pt());
516+
dataTrackHist.fill(HIST("Track3/eta"), track3.eta());
517+
dataTrackHist.fill(HIST("Track3/phi"), track3.phi());
518+
dataTrackHist.fill(HIST("Track3/tpcCluster"), track3.tpcNClsFound());
519+
dataTrackHist.fill(HIST("Track3/tpcCrossed"), track3.tpcNClsCrossedRows());
520+
dataTrackHist.fill(HIST("Track3/tpcShared"), track3.tpcNClsShared());
521+
dataTrackHist.fill(HIST("Track3/itsCluster"), track3.itsNCls());
522+
dataTrackHist.fill(HIST("Track3/itsIbCluster"), track3.itsNClsInnerBarrel());
523+
dataTrackHist.fill(HIST("Track3/itsNsigma"), track3.p(), nsigma2[0]);
524+
dataTrackHist.fill(HIST("Track3/tpcNsigma"), track3.p(), nsigma2[1]);
525+
dataTrackHist.fill(HIST("Track3/tpctofNsigma"), track3.p(), std::hypot(nsigma3[1], nsigma3[2]));
526+
527+
triplets++;
528+
529+
vecq3.push_back(femtoDream::FemtoDreamMath::getQ3(track1, femtoDream::getMass(trackCuts1.pdgCode.value), track2, femtoDream::getMass(trackCuts2.pdgCode.value), track3, femtoDream::getMass(trackCuts3.pdgCode.value)));
530+
}
531+
}
532+
}
533+
return triplets;
534+
}
535+
420536
template <typename T>
421537
int countGenPairs(const T& tracks)
422538
{
@@ -444,6 +560,42 @@ struct FemtoDreamPairEfficiency {
444560
return pairs;
445561
}
446562

563+
template <typename T>
564+
int countGenTriplets(const T& tracks)
565+
{
566+
int triplets = 0;
567+
for (auto track1 = tracks.begin(); track1 != tracks.end(); track1++) {
568+
if (!track1.isPhysicalPrimary() ||
569+
track1.pdgCode() != trackCuts1.pdgCode.value ||
570+
track1.pt() < trackCuts1.ptMin.value ||
571+
track1.pt() > trackCuts1.ptMax.value ||
572+
std::fabs(track1.eta()) > trackCuts1.etaAbsMax.value) {
573+
continue;
574+
}
575+
for (auto track2 = track1 + 1; track2 != tracks.end(); track2++) {
576+
if (!track2.isPhysicalPrimary() ||
577+
track2.pdgCode() != trackCuts2.pdgCode.value ||
578+
track2.pt() < trackCuts2.ptMin.value ||
579+
track2.pt() > trackCuts2.ptMax.value ||
580+
std::fabs(track2.eta()) > trackCuts2.etaAbsMax.value) {
581+
continue;
582+
}
583+
for (auto track3 = track2 + 1; track3 != tracks.end(); track3++) {
584+
if (!track3.isPhysicalPrimary() ||
585+
track3.pdgCode() != trackCuts3.pdgCode.value ||
586+
track3.pt() < trackCuts3.ptMin.value ||
587+
track3.pt() > trackCuts3.ptMax.value ||
588+
std::fabs(track3.eta()) > trackCuts3.etaAbsMax.value) {
589+
continue;
590+
}
591+
vecq3.push_back(femtoDream::FemtoDreamMath::getQ3(track1, femtoDream::getMass(trackCuts1.pdgCode.value), track2, femtoDream::getMass(trackCuts2.pdgCode.value), track3, femtoDream::getMass(trackCuts3.pdgCode.value)));
592+
triplets++;
593+
}
594+
}
595+
}
596+
return triplets;
597+
}
598+
447599
void processdNdetaData(Collisions::iterator const& collision, FilteredFullTracks const& tracks)
448600
{
449601
if (!checkEventSelections<false>(collision, true))
@@ -452,16 +604,33 @@ struct FemtoDreamPairEfficiency {
452604
auto tracksWithItsPid = soa::Attach<FilteredFullTracks, aod::pidits::ITSNSigmaEl, aod::pidits::ITSNSigmaPi, aod::pidits::ITSNSigmaKa,
453605
aod::pidits::ITSNSigmaPr, aod::pidits::ITSNSigmaDe, aod::pidits::ITSNSigmaTr, aod::pidits::ITSNSigmaHe>(tracks);
454606

455-
vecKstar.clear();
456-
pairsFound = countRecoPairs(tracksWithItsPid);
457-
458-
if (pairsFound == 0) {
459-
return;
607+
if (mode.countPairs.value) {
608+
vecKstar.clear();
609+
pairsFound = countRecoPairs(tracksWithItsPid);
610+
if (pairsFound == 0) {
611+
return;
612+
}
613+
}
614+
if (mode.countTriplets.value) {
615+
vecq3.clear();
616+
tripletsFound = countRecoTriplets(tracksWithItsPid);
617+
if (tripletsFound == 0) {
618+
return;
619+
}
460620
}
621+
461622
dataEventHist.fill(HIST("hDataEventSelection"), 4); // found a pair
462623

463-
if (*std::min_element(vecKstar.begin(), vecKstar.end()) > eventSelection.kstarMax) {
464-
return;
624+
if (mode.countPairs.value) {
625+
if (*std::min_element(vecKstar.begin(), vecKstar.end()) > eventSelection.kstarMax) {
626+
return;
627+
}
628+
}
629+
630+
if (mode.countTriplets.value) {
631+
if (*std::min_element(vecq3.begin(), vecq3.end()) > eventSelection.q3Max) {
632+
return;
633+
}
465634
}
466635

467636
dataEventHist.fill(HIST("hDataEventSelection"), 5); // found a lowkstar pair
@@ -491,17 +660,35 @@ struct FemtoDreamPairEfficiency {
491660
const auto& mcCollision = collision.mcCollision_as<GenCollisions>();
492661
auto mcParticlesThisColl = mcParticles.sliceByCached(mcparticle::mcCollisionId, mcCollision.globalIndex(), cache);
493662

494-
vecKstar.clear();
495-
pairsFound = countGenPairs(mcParticlesThisColl);
496-
497-
if (pairsFound == 0) {
498-
return;
663+
if (mode.countPairs.value) {
664+
vecKstar.clear();
665+
pairsFound = countGenPairs(mcParticlesThisColl);
666+
if (pairsFound == 0) {
667+
return;
668+
}
499669
}
670+
if (mode.countTriplets.value) {
671+
vecq3.clear();
672+
tripletsFound = countGenTriplets(mcParticlesThisColl);
673+
if (tripletsFound == 0) {
674+
return;
675+
}
676+
}
677+
500678
mcEventHist.fill(HIST("hRecoEventSelection"), 5);
501679

502-
if (*std::min_element(vecKstar.begin(), vecKstar.end()) > eventSelection.kstarMax) {
503-
return;
680+
if (mode.countPairs.value) {
681+
if (*std::min_element(vecKstar.begin(), vecKstar.end()) > eventSelection.kstarMax) {
682+
return;
683+
}
684+
}
685+
686+
if (mode.countTriplets.value) {
687+
if (*std::min_element(vecq3.begin(), vecq3.end()) > eventSelection.q3Max) {
688+
return;
689+
}
504690
}
691+
505692
mcEventHist.fill(HIST("hRecoEventSelection"), 6);
506693

507694
float genCentrality = mcCollision.centFT0M();
@@ -551,17 +738,37 @@ struct FemtoDreamPairEfficiency {
551738
return;
552739
}
553740
mcEventHist.fill(HIST("hGenEventSelection"), 2);
554-
vecKstar.clear();
555-
pairsFound = countGenPairs(mcParticles);
556-
if (pairsFound == 0) {
557-
return;
741+
742+
if (mode.countPairs.value) {
743+
vecKstar.clear();
744+
pairsFound = countGenPairs(mcParticles);
745+
if (pairsFound == 0) {
746+
return;
747+
}
748+
}
749+
if (mode.countTriplets.value) {
750+
vecq3.clear();
751+
tripletsFound = countGenTriplets(mcParticles);
752+
if (tripletsFound == 0) {
753+
return;
754+
}
558755
}
756+
559757
mcEventHist.fill(HIST("hGenEventSelection"), 3);
560-
if (*std::min_element(vecKstar.begin(), vecKstar.end()) > eventSelection.kstarMax) {
561-
return;
758+
759+
if (mode.countPairs.value) {
760+
if (*std::min_element(vecKstar.begin(), vecKstar.end()) > eventSelection.kstarMax) {
761+
return;
762+
}
763+
}
764+
765+
if (mode.countTriplets.value) {
766+
if (*std::min_element(vecq3.begin(), vecq3.end()) > eventSelection.q3Max) {
767+
return;
768+
}
562769
}
563-
mcEventHist.fill(HIST("hGenEventSelection"), 4);
564770

771+
mcEventHist.fill(HIST("hGenEventSelection"), 4);
565772
mcEventHist.fill(HIST("hGenMcVertexZ"), mcCollision.posZ());
566773
mcEventHist.fill(HIST("hGenMcMultiplicity"), mcCollision.multMCNParticlesEta08());
567774
mcEventHist.fill(HIST("hGenMcCentrality"), mcCollision.centFT0M());

0 commit comments

Comments
 (0)