5656#include < cstdint>
5757#include < cstdlib>
5858#include < string>
59+ #include < utility>
5960#include < vector>
6061
6162using namespace o2 ;
@@ -113,9 +114,10 @@ struct HfTaskElectronWeakBoson {
113114 Configurable<bool > isTHnElectron{" isTHnElectron" , true , " Enables THn for electrons" };
114115 Configurable<float > ptTHnThresh{" ptTHnThresh" , 5.0 , " Threshold for THn make" };
115116
116- // Skimmed dataset processing configurations
117+ // Skimmed (trigger) dataset processing configurations
117118 Configurable<bool > cfgSkimmedProcessing{" cfgSkimmedProcessing" , true , " Enables processing of skimmed datasets" };
118119 Configurable<std::string> cfgTriggerName{" cfgTriggerName" , " fGammaHighPtEMCAL" , " Trigger of interest (comma separated for multiple)" };
120+ Configurable<bool > applySel8{" applySel8" , true , " Apply sel8 filter or not" };
119121
120122 // CCDB service configurations
121123 Configurable<std::string> cfgCCDBPath{" cfgCCDBPath" , " Users/m/mpuccio/EventFiltering/OTS/" , " Path to CCDB for trigger data" };
@@ -130,19 +132,20 @@ struct HfTaskElectronWeakBoson {
130132 Configurable<bool > enableCentralityAnalysis{" enableCentralityAnalysis" , true , " Enable centrality-dependent analysis" };
131133 Configurable<float > centralityMin{" centralityMin" , -1 , " minimum cut on centrality selection" };
132134 Configurable<float > centralityMax{" centralityMax" , 101 , " maximum cut on centrality selection" };
135+ Configurable<std::vector<double >> centralityBins{" centralityBins" , {0 , 20 , 60 , 100 }, " centrality bins" };
133136
137+ // QA for Z->ee
138+ Configurable<bool > enableZeeRecoQA{" enableZeeRecoQA" , false , " Enable QA for Z->ee reconstruction" };
134139 // CCDB service object
135140 Service<o2::ccdb::BasicCCDBManager> ccdb;
136141
137142 struct HfElectronCandidate {
138- float pt, eta, phi, energy;
139- int charge;
140- HfElectronCandidate (float ptr, float e, float ph, float en, int ch)
141- : pt(ptr), eta(e), phi(ph), energy(en), charge(ch) {}
142-
143- int sign () const { return charge; }
143+ float pt, eta, phi, eop, energyIso, momIso, ntrackIso;
144+ HfElectronCandidate (float ptr, float e, float ph, float ep, float eiso, float piso, int ntrkiso)
145+ : pt(ptr), eta(e), phi(ph), eop(ep), energyIso(eiso), momIso(piso), ntrackIso(ntrkiso) {}
144146 };
145147 std::vector<HfElectronCandidate> selectedElectronsIso;
148+ std::vector<HfElectronCandidate> selectedPositronsIso;
146149 std::vector<HfElectronCandidate> selectedElectronsAss;
147150
148151 struct HfZeeCandidate {
@@ -161,8 +164,7 @@ struct HfTaskElectronWeakBoson {
161164 // pp
162165 // using TrackEle = o2::soa::Filtered<o2::soa::Join<o2::aod::Tracks, o2::aod::FullTracks, o2::aod::TracksDCA, o2::aod::TrackSelection, o2::aod::pidTPCEl, o2::aod::pidTOFEl>>;
163166
164- // Filter
165- Filter eventFilter = (o2::aod::evsel::sel8 == true );
167+ Filter eventFilter = (applySel8 ? (o2::aod::evsel::sel8 == true ) : (o2::aod::evsel::sel8 == o2::aod::evsel::sel8));
166168 Filter posZFilter = (nabs(o2::aod::collision::posZ) < vtxZ);
167169
168170 Filter etafilter = (aod::track::eta < etaTrMax) && (aod::track::eta > etaTrMin);
@@ -209,6 +211,9 @@ struct HfTaskElectronWeakBoson {
209211 const AxisSpec axisCounter{1 , 0 , 1 , " events" };
210212 const AxisSpec axisEta{20 , -1.0 , 1.0 , " #eta" };
211213 const AxisSpec axisPt{nBinsPt, 0 , binPtmax, " p_{T}" };
214+ const AxisSpec axisPtZee{60 , 20 , 80 , " p_{T}" };
215+ const AxisSpec axisPtZele{60 , 20 , 80 , " p_{T,ele} (GeV/c)" };
216+ const AxisSpec axisPtZpos{60 , 20 , 80 , " p_{T,pos} (GeV/c)" };
212217 const AxisSpec axisNsigma{100 , -5 , 5 , " N#sigma" };
213218 const AxisSpec axisDedx{150 , 0 , 150 , " dEdx" };
214219 const AxisSpec axisE{nBinsE, 0 , binEmax, " Energy" };
@@ -218,20 +223,30 @@ struct HfTaskElectronWeakBoson {
218223 const AxisSpec axisdR{20 , 0.0 , 0.2 , " dR" };
219224 const AxisSpec axisNcell{50 , 0.0 , 50.0 , " Ncell" };
220225 const AxisSpec axisPhi{350 , 0 , 7 , " Phi" };
221- const AxisSpec axisEop{200 , 0 , 2 , " Eop" };
226+ const AxisSpec axisEop{200 , 0 , 2 , " E/p" };
227+ const AxisSpec axisEopZele{200 , 0 , 2 , " E/p electon" };
228+ const AxisSpec axisEopZpos{200 , 0 , 2 , " E/p positron" };
222229 const AxisSpec axisChi2{250 , 0.0 , 25.0 , " #chi^{2}" };
223230 const AxisSpec axisCluster{100 , 0.0 , 200.0 , " counts" };
224231 const AxisSpec axisITSNCls{10 , 0.0 , 10 , " counts" };
225232 const AxisSpec axisEMCtime{100 , -50.0 , 50 , " EMC time" };
226- const AxisSpec axisIsoEnergy{100 , 0 , 1.0 , " Isolation energy(GeV/C)" };
227- const AxisSpec axisIsoTrack{15 , -0.5 , 14.5 , " Isolation Track" };
228- const AxisSpec axisInvMassZ{150 , 0 , 150 , " M_{ee} (GeV/c^{2})" };
233+ const AxisSpec axisIsoEnergy{100 , 0 , 1.0 , " E_{iso}" };
234+ const AxisSpec axisIsoEnergyZele{100 , 0 , 1.0 , " E_{iso,ele}" };
235+ const AxisSpec axisIsoEnergyZpos{100 , 0 , 1.0 , " E_{iso,pos}" };
236+ const AxisSpec axisIsoMomentum{100 , 0 , 10.0 , " Isolation momentum(GeV/C)" };
237+ const AxisSpec axisIsoMomentumZele{100 , 0 , 10.0 , " p_{iso,ele}" };
238+ const AxisSpec axisIsoMomentumZpos{100 , 0 , 10.0 , " p_{iso,pos}" };
239+ const AxisSpec axisIsoTrack{25 , -0.5 , 24.5 , " Isolation Track" };
240+ const AxisSpec axisIsoTrackZele{25 , -0.5 , 24.5 , " N_{isotrk,ele}" };
241+ const AxisSpec axisIsoTrackZpos{25 , -0.5 , 24.5 , " N_{isotrk,pos}" };
242+ const AxisSpec axisInvMassZgamma{150 , 0 , 150 , " M_{ee} (GeV/c^{2})" };
243+ const AxisSpec axisInvMassZ{130 , 20 , 150 , " M_{ee} (GeV/c^{2})" };
229244 const AxisSpec axisTrigger{3 , -0.5 , 2.5 , " Trigger status of zorro" };
230245 const AxisSpec axisDPhiZh{64 , -o2::constants::math::PIHalf, 3 * o2::constants::math::PIHalf, " #Delta#phi(Z-h)" };
231246 const AxisSpec axisPtHadron{50 , 0 , 50 , " p_{T,hadron} (GeV/c)" };
232247 const AxisSpec axisPtZ{150 , 0 , 150 , " p_{T,Z} (GeV/c)" };
233248 const AxisSpec axisSign{2 , -2 , 2 , " charge sign" };
234- const AxisSpec axisCentrality{10 , 0 , 100 , " Centrality (%) " };
249+ const AxisSpec axisCentrality{centralityBins };
235250 const AxisSpec axisPtRatio{200 , 0 , 2.0 , " pt ratio for h and Z" };
236251
237252 // create registrygrams
@@ -259,14 +274,14 @@ struct HfTaskElectronWeakBoson {
259274 registry.add (" hEopNsigTPC" , " Eop vs. Nsigma" , kTH2F , {{axisNsigma}, {axisEop}});
260275 registry.add (" hEMCtime" , " EMC timing" , kTH1F , {axisEMCtime});
261276 registry.add (" hIsolationEnergy" , " Isolation Energy" , kTH2F , {{axisE}, {axisIsoEnergy}});
262- registry.add (" hIsolationTrack " , " Isolation Track " , kTH2F , {{axisE}, {axisIsoTrack} });
263- registry.add (" hInvMassZee " , " invariant mass for Z ULS pair" , HistType::kTHnSparseF , {axisSign, axisPt, axisInvMassZ });
264- registry.add (" hKfInvMassZee " , " invariant mass for Z ULS pair KFp " , HistType::kTHnSparseF , {axisSign, axisPt, axisInvMassZ });
277+ registry.add (" hInvMassZee " , " invariant mass for Z ULS pair " , HistType:: kTHnSparseF , {axisCentrality, axisSign, axisPt, axisInvMassZgamma });
278+ registry.add (" hKfInvMassZee " , " invariant mass for Z ULS pair KFp " , HistType::kTHnSparseF , {axisCentrality, axisSign, axisPt, axisInvMassZgamma });
279+ registry.add (" hInvMassZeeQA " , " QA for invariant mass for Z" , HistType::kTHnSparseF , {axisInvMassZ, axisPtZele, axisPtZpos, axisEopZele, axisEopZpos, axisIsoEnergyZele, axisIsoEnergyZpos, axisIsoMomentumZele, axisIsoMomentumZpos, axisIsoTrackZele, axisIsoTrackZpos });
265280 registry.add (" hTHnElectrons" , " electron info" , HistType::kTHnSparseF , {axisPt, axisNsigma, axisM02, axisEop, axisIsoEnergy, axisIsoTrack, axisEta, axisDedx});
266281 registry.add (" hTHnTrMatch" , " Track EMC Match" , HistType::kTHnSparseF , {axisPt, axisdPhi, axisdEta});
267282
268283 // Z-hadron correlation histograms
269- registry.add (" hZHadronDphi" , " Z-hadron #Delta#phi correlation" , HistType::kTHnSparseF , {axisSign, axisPtZ, axisDPhiZh, axisPtRatio, axisPtHadron});
284+ registry.add (" hZHadronDphi" , " Z-hadron #Delta#phi correlation" , HistType::kTHnSparseF , {axisCentrality, axisSign, axisPtZ, axisDPhiZh, axisPtRatio, axisPtHadron});
270285 registry.add (" hZptSpectrum" , " Z boson p_{T} spectrum" , kTH2F , {{axisSign}, {axisPtZ}});
271286
272287 // hisotgram for EMCal trigger
@@ -306,12 +321,15 @@ struct HfTaskElectronWeakBoson {
306321
307322 return (isoEnergy);
308323 }
309- int getIsolatedTrack (double etaEle,
310- double phiEle,
311- float ptEle ,
312- TrackEle const & tracks)
324+ std::pair< int , double > getIsolatedTrack (double etaEle,
325+ double phiEle,
326+ float pEle ,
327+ TrackEle const & tracks)
313328 {
314329 int trackCount = 0 ;
330+ double isoMomentum = 10 ;
331+ double pSum = 0.0 ;
332+ // LOG(info) << "track p = " << pEle;
315333
316334 for (const auto & track : tracks) {
317335
@@ -323,16 +341,22 @@ struct HfTaskElectronWeakBoson {
323341
324342 if (deltaR < rIsolation) {
325343 trackCount++;
344+ pSum += track.p ();
326345 }
327346 }
328347
329- registry.fill (HIST (" hIsolationTrack" ), ptEle, trackCount);
348+ // LOG(info) << "momSun = " << pSum;
349+ if (pSum > 0 ) {
350+ isoMomentum = pSum / pEle - 1.0 ;
351+ }
330352
331- return (trackCount);
353+ // LOG(info) << "isop = " << isoMomentum;
354+ return std::make_pair (trackCount - 1 , isoMomentum);
332355 }
333356
334357 void recoMassZee (KFParticle kfpIsoEle,
335358 int charge,
359+ float centrality,
336360 TrackEle const & tracks)
337361 {
338362 // LOG(info) << "Invarimass cal by KF particle ";
@@ -362,7 +386,7 @@ struct HfTaskElectronWeakBoson {
362386 auto child2 = RecoDecayPtEtaPhi::pVector (kfpAssEle.GetPt () * correctionPtElectron, kfpAssEle.GetEta (), kfpAssEle.GetPhi ());
363387 double invMassEE = RecoDecay::m (std::array{child1, child2}, std::array{o2::constants::physics::MassElectron, o2::constants::physics::MassElectron});
364388
365- registry.fill (HIST (" hInvMassZee" ), track.sign () * charge, kfpIsoEle.GetPt (), invMassEE);
389+ registry.fill (HIST (" hInvMassZee" ), centrality, track.sign () * charge, kfpIsoEle.GetPt (), invMassEE);
366390
367391 // reco by KFparticle
368392 const KFParticle* electronPairs[2 ] = {&kfpIsoEle, &kfpAssEle};
@@ -382,7 +406,7 @@ struct HfTaskElectronWeakBoson {
382406 }
383407 float massZee, massZeeErr;
384408 zeeKF.GetMass (massZee, massZeeErr);
385- registry.fill (HIST (" hKfInvMassZee" ), track.sign () * charge, kfpIsoEle.GetPt (), massZee);
409+ registry.fill (HIST (" hKfInvMassZee" ), centrality, track.sign () * charge, kfpIsoEle.GetPt (), massZee);
386410 // LOG(info) << "Invarimass cal by KF particle mass = " << massZee;
387411 // LOG(info) << "Invarimass cal by RecoDecay = " << invMassEE;
388412 reconstructedZ.emplace_back (
@@ -448,6 +472,7 @@ struct HfTaskElectronWeakBoson {
448472 }
449473 // initialze for inclusive-electron
450474 selectedElectronsIso.clear ();
475+ selectedPositronsIso.clear ();
451476 selectedElectronsAss.clear ();
452477 reconstructedZ.clear ();
453478
@@ -460,8 +485,9 @@ struct HfTaskElectronWeakBoson {
460485 registry.fill (HIST (" hZvtx" ), collision.posZ ());
461486
462487 // Calculate centrality
488+ float centrality = 1.0 ;
463489 if (enableCentralityAnalysis) {
464- float centrality = o2::hf_centrality::getCentralityColl (collision, centralityEstimator);
490+ centrality = o2::hf_centrality::getCentralityColl (collision, centralityEstimator);
465491 // LOG(info) << centrality;
466492 if (centrality < centralityMin || centrality > centralityMax) {
467493 return ;
@@ -502,38 +528,46 @@ struct HfTaskElectronWeakBoson {
502528 registry.fill (HIST (" hPt" ), track.pt ());
503529 registry.fill (HIST (" hTPCNsigma" ), track.p (), track.tpcNSigmaEl ());
504530
505- float energyTrk = 0.0 ;
531+ float eop = 0.0 ;
532+ float isoEnergy = 1.0 ;
533+ // track isolation
534+ auto [trackCount, isoMomentum] = getIsolatedTrack (track.eta (), track.phi (), track.p (), tracks);
535+ // LOG(info) << "isoMomentum = " << isoMomentum;
506536
507537 if (track.pt () > ptAssMin) {
508538 selectedElectronsAss.emplace_back (
509539 track.pt (),
510540 track.eta (),
511541 track.phi (),
512- energyTrk,
513- track.sign ());
542+ eop,
543+ isoEnergy,
544+ isoMomentum,
545+ trackCount);
514546 }
515547
516548 if (track.pt () < ptMin) {
517549 continue ;
518550 }
519- // track - match
520551
521- // continue;
522- if (track.phi () < phiEmcMin || track.phi () > phiEmcMax)
523- continue ;
524- if (std::abs (track.eta ()) > etaEmcMax)
525- continue ;
552+ // LOG(info) << "tr phi, eta = " << track.phi() << " ; " << track.eta();
553+ // EMC acc
554+ bool isEMCacceptance = true ;
555+ if (track.phi () < phiEmcMin || track.phi () > phiEmcMax) {
556+ isEMCacceptance = false ;
557+ }
558+ if (std::abs (track.eta ()) > etaEmcMax) {
559+ isEMCacceptance = false ;
560+ }
561+ // LOG(info) << "EMC acc = " << isEMCacceptance;
526562 auto tracksofcluster = matchedtracks.sliceBy (perClusterMatchedTracks, track.globalIndex ());
527563
528- // LOGF(info, "Number of matched track: %d", tracksofcluster.size());
529-
530564 double rMin = 999.9 ;
531565 double dPhiMin = 999.9 ;
532566 double dEtaMin = 999.9 ;
533567 bool isIsolated = false ;
534568 bool isIsolatedTr = false ;
535569
536- if (tracksofcluster.size ()) {
570+ if (tracksofcluster.size () && isEMCacceptance ) {
537571 int nMatch = 0 ;
538572 for (const auto & match : tracksofcluster) {
539573 if (match.emcalcluster_as <SelectedClusters>().time () < timeEmcMin || match.emcalcluster_as <SelectedClusters>().time () > timeEmcMax)
@@ -557,6 +591,7 @@ struct HfTaskElectronWeakBoson {
557591 registry.fill (HIST (" hMatchEta" ), etaEmc, match.track_as <TrackEle>().trackEtaEmcal ());
558592
559593 double r = RecoDecay::sqrtSumOfSquares (dPhi, dEta);
594+ // LOG(info) << "r match = " << r;
560595 if (r < rMin) {
561596 rMin = r;
562597 dPhiMin = dPhi;
@@ -574,11 +609,10 @@ struct HfTaskElectronWeakBoson {
574609
575610 const auto & cluster = match.emcalcluster_as <SelectedClusters>();
576611
577- double eop = energyEmc / match.track_as <TrackEle>().p ();
612+ eop = energyEmc / match.track_as <TrackEle>().p ();
613+ // LOG(info) << "eop = " << eop;
578614
579- double isoEnergy = getIsolatedCluster (cluster, emcClusters);
580-
581- int trackCount = getIsolatedTrack (track.eta (), track.phi (), track.pt (), tracks) - 1 ;
615+ isoEnergy = getIsolatedCluster (cluster, emcClusters);
582616
583617 if (match.track_as <TrackEle>().pt () > ptTHnThresh && isTHnElectron) {
584618 registry.fill (HIST (" hTHnElectrons" ), match.track_as <TrackEle>().pt (), match.track_as <TrackEle>().tpcNSigmaEl (), m02Emc, eop, isoEnergy, trackCount, track.eta (), track.tpcSignal ());
@@ -606,29 +640,44 @@ struct HfTaskElectronWeakBoson {
606640 }
607641 KFPTrack kfpTrackIsoEle = createKFPTrackFromTrack (match.track_as <TrackEle>());
608642 KFParticle kfpIsoEle (kfpTrackIsoEle, pdgIso);
609- recoMassZee (kfpIsoEle, match.track_as <TrackEle>().sign (), tracks);
610-
611- selectedElectronsIso.emplace_back (
612- match.track_as <TrackEle>().pt (),
613- match.track_as <TrackEle>().eta (),
614- match.track_as <TrackEle>().phi (),
615- energyEmc,
616- match.track_as <TrackEle>().sign ());
617- }
618- }
643+ recoMassZee (kfpIsoEle, match.track_as <TrackEle>().sign (), centrality, tracks);
644+
645+ } // end of pt cut for e from Z
646+ } // end if isolation cut
619647 if (isIsolatedTr) {
620648 registry.fill (HIST (" hEopIsolationTr" ), match.track_as <TrackEle>().pt (), eop);
621649 }
622- }
623- }
650+ } // end of PID cut
651+ } // end of nmatch == 0
624652 nMatch++;
625- }
626- }
653+ } // end of cluster match
654+ } // end of cluster
627655
628656 if (rMin < rMatchMax) {
629657 // LOG(info) << "R mim = " << rMin;
630658 registry.fill (HIST (" hTrMatch_mim" ), dPhiMin, dEtaMin);
631659 }
660+ if (enableZeeRecoQA && track.pt () > ptZeeMin) {
661+ if (track.sign () < 0 ) {
662+ selectedElectronsIso.emplace_back (
663+ track.pt (),
664+ track.eta (),
665+ track.phi (),
666+ eop,
667+ isoEnergy,
668+ isoMomentum,
669+ trackCount);
670+ } else {
671+ selectedPositronsIso.emplace_back (
672+ track.pt (),
673+ track.eta (),
674+ track.phi (),
675+ eop,
676+ isoEnergy,
677+ isoMomentum,
678+ trackCount);
679+ }
680+ }
632681
633682 } // end of track loop
634683 // Z-hadron
@@ -649,10 +698,23 @@ struct HfTaskElectronWeakBoson {
649698 // calculate Z-h correlation
650699 double deltaPhi = RecoDecay::constrainAngle (trackAss.phi - zBoson.phi , -o2::constants::math::PIHalf);
651700 double ptRatio = trackAss.pt / zBoson.pt ;
652- registry.fill (HIST (" hZHadronDphi" ), zBoson.charge , zBoson.pt , deltaPhi, ptRatio, trackAss.pt );
701+ registry.fill (HIST (" hZHadronDphi" ), centrality, zBoson.charge , zBoson.pt , deltaPhi, ptRatio, trackAss.pt );
653702 }
654703 }
655704 } // end of Z-hadron correlation
705+ // Z->ee QA
706+ if (enableZeeRecoQA) {
707+ if (selectedElectronsIso.size () > 0 && selectedPositronsIso.size () > 0 ) {
708+ for (const auto & trackEle : selectedElectronsIso) {
709+ for (const auto & trackPos : selectedPositronsIso) {
710+ auto child1 = RecoDecayPtEtaPhi::pVector (trackEle.pt , trackEle.eta , trackEle.phi );
711+ auto child2 = RecoDecayPtEtaPhi::pVector (trackPos.pt , trackPos.eta , trackPos.phi );
712+ double invMass = RecoDecay::m (std::array{child1, child2}, std::array{o2::constants::physics::MassElectron, o2::constants::physics::MassElectron});
713+ registry.fill (HIST (" hInvMassZeeQA" ), invMass, trackEle.pt , trackPos.pt , trackEle.eop , trackPos.eop , trackEle.energyIso , trackPos.energyIso , trackEle.momIso , trackPos.momIso , trackEle.ntrackIso , trackPos.ntrackIso );
714+ }
715+ }
716+ }
717+ } // end of Z->ee QA
656718 }
657719};
658720
0 commit comments