4141#include < TMath.h>
4242#include < TPDGCode.h>
4343
44+ #include < cstdlib>
4445#include < vector>
4546
4647using namespace o2 ;
@@ -55,13 +56,14 @@ struct Alice3Lepton {
5556 kCe = 0 ,
5657 kBe = 1 ,
5758 kBCe = 2 ,
58- kPPi0 = 4 ,
59+ kPPi0 = 3 ,
60+ kHFE = 4 ,
5961 };
6062
6163 Service<o2::framework::O2DatabasePDG> inspdg;
6264
6365 Configurable<int > pdg{" pdg" , 11 , " pdg code for analysis. dielectron:11, dimuon:13" };
64- Configurable<int > requireHFE{" requireHFE" , -1 , " -1: no selection, 0: charm, 1: direct beauty, 2: beauty->charm->e, 3: HFE , 4: promptPi0 " };
66+ Configurable<int > requireHFE{" requireHFE" , -1 , " -1: no selection, 0: charm, 1: direct beauty, 2: beauty->charm->e, 3: promptPi0 , 4: HFE " };
6567 Configurable<float > ptMin{" ptMin" , 0 .f , " Lower limit in pT" };
6668 Configurable<float > ptMax{" ptMax" , 5 .f , " Upper limit in pT" };
6769 Configurable<float > etaMin{" etaMin" , -5 .f , " Lower limit in eta" };
@@ -92,7 +94,7 @@ struct Alice3Lepton {
9294 const AxisSpec axisTrackLengthOuterTOF{300 , 0 ., 300 ., " Track length (cm)" };
9395 const AxisSpec axisEta{1000 , -5 , 5 , " #it{#eta}" };
9496 const AxisSpec axisDCAxysigma{1000 , -20 , 20 , " DCA_{xy} (#sigma)" };
95- const AxisSpec axisDCAxy{1000 , -1000 , 1000 , " DCA_{xy} (#micro m)" };
97+ const AxisSpec axisDCAxy{1200 , -300 , 300 , " DCA_{xy} (#micro m)" };
9698 const AxisSpec axisPhi{360 , 0 , TMath::TwoPi (), " #it{#varphi} (rad.)" };
9799 const AxisSpec axisProdx{2000 , -100 , 100 , " Prod. Vertex X (cm)" };
98100 const AxisSpec axisPrody{2000 , -100 , 100 , " Prod. Vertex Y (cm)" };
@@ -162,8 +164,9 @@ struct Alice3Lepton {
162164 template <typename TMCParticle1, typename TMCParticles>
163165 int IsHF (TMCParticle1 const & p1, TMCParticles const & mcparticles)
164166 {
165- if (!p1.has_mothers ())
167+ if (!p1.has_mothers ()) {
166168 return HFType::kUndef ;
169+ }
167170
168171 int motherid_p1 = p1.mothersIds ()[0 ];
169172 if (motherid_p1 > -1 ) {
@@ -285,15 +288,12 @@ struct Alice3Lepton {
285288 continue ;
286289 }
287290 }
288- if (std::abs (mcParticle.pdgCode ()) != pdg) {
289- continue ;
290- }
291291 if (requireHFE > -1 ) {
292292 int typehfe = IsHF (mcParticle, mcParticles);
293- if (requireHFE < 3 && typehfe != requireHFE) {
293+ if (requireHFE < HFType:: kHFE && typehfe != requireHFE) {
294294 continue ;
295295 }
296- if (requireHFE == 3 && typehfe == HFType::kUndef ) {
296+ if (requireHFE == HFType:: kHFE && ( typehfe == HFType::kUndef || typehfe == kPPi0 ) ) {
297297 continue ;
298298 }
299299 }
@@ -431,10 +431,10 @@ struct Alice3Lepton {
431431 }
432432 if (requireHFE > -1 ) {
433433 int typehfe = IsHF (mcParticle, mcParticles);
434- if (requireHFE < 3 && typehfe != requireHFE) {
434+ if (requireHFE < HFType:: kHFE && typehfe != requireHFE) {
435435 continue ;
436436 }
437- if (requireHFE == 3 && typehfe == HFType::kUndef ) {
437+ if (requireHFE == HFType:: kHFE && ( typehfe == HFType::kUndef || typehfe == kPPi0 ) ) {
438438 continue ;
439439 }
440440 }
@@ -528,7 +528,8 @@ struct Alice3Dilepton {
528528 SliceCache cache_rec;
529529
530530 Configurable<int > pdg{" pdg" , 11 , " pdg code for analysis. dielectron:11, dimuon:13" };
531- Configurable<bool > requireHFEid{" requireHFEid" , true , " Require HFE identification for both leptons" };
531+ Configurable<bool > requireHFEid{" requireHFEid" , true , " Require HFE identification" };
532+ Configurable<bool > contamination{" contamination" , false , " Fill only pairs with one misidentifixed electrons" };
532533 Configurable<float > ptMin{" ptMin" , 0 .f , " Lower limit in pT" };
533534 Configurable<float > ptMax{" ptMax" , 5 .f , " Upper limit in pT" };
534535 Configurable<float > etaMin{" etaMin" , -5 .f , " Lower limit in eta" };
@@ -629,13 +630,40 @@ struct Alice3Dilepton {
629630 }
630631 }
631632
633+ template <typename TTrack, typename TMCParticles>
634+ bool IsHF (TTrack const & track, TMCParticles const & mcparticles)
635+ {
636+ if (!track.has_mcParticle ()) {
637+ return false ;
638+ }
639+ const auto p1 = track.template mcParticle_as <aod::McParticles>();
640+ if (!p1.has_mothers ()) {
641+ return false ;
642+ }
643+ int motherId = p1.mothersIds ()[0 ];
644+ while (motherId > -1 ) {
645+ auto mp = mcparticles.rawIteratorAt (motherId);
646+ if (((500 < std::abs (mp.pdgCode ()) && std::abs (mp.pdgCode ()) < 599 ) || (5000 < std::abs (mp.pdgCode ()) && std::abs (mp.pdgCode ()) < 5999 )) || ((400 < std::abs (mp.pdgCode ()) && std::abs (mp.pdgCode ()) < 499 ) || (4000 < std::abs (mp.pdgCode ()) && std::abs (mp.pdgCode ()) < 4999 ))) {
647+ return true ;
648+ }
649+ if (mp.has_mothers ()) {
650+ motherId = mp.mothersIds ()[0 ];
651+ } else {
652+ motherId = -999 ;
653+ }
654+ }
655+ return false ;
656+ }
657+
632658 template <typename TMCParticle1, typename TMCParticle2, typename TMCParticles>
633659 int IsSameMother (TMCParticle1 const & p1, TMCParticle2 const & p2, TMCParticles const & mcparticles)
634660 {
635- if (!p1.has_mothers ())
661+ if (!p1.has_mothers ()) {
636662 return -1 ;
637- if (!p2.has_mothers ())
663+ }
664+ if (!p2.has_mothers ()) {
638665 return -1 ;
666+ }
639667
640668 int motherid1 = p1.mothersIds ()[0 ];
641669 auto mother1 = mcparticles.iteratorAt (motherid1);
@@ -645,10 +673,12 @@ struct Alice3Dilepton {
645673 auto mother2 = mcparticles.iteratorAt (motherid2);
646674 int mother2_pdg = mother2.pdgCode ();
647675
648- if (motherid1 != motherid2)
676+ if (motherid1 != motherid2) {
649677 return -1 ;
650- if (mother1_pdg != mother2_pdg)
678+ }
679+ if (mother1_pdg != mother2_pdg) {
651680 return -1 ;
681+ }
652682
653683 if (std::abs (mother1_pdg) != PDG_t::kGamma // photon
654684 && std::abs (mother1_pdg) != PDG_t::kPi0 // pi0
@@ -674,10 +704,12 @@ struct Alice3Dilepton {
674704 // 1. b->e- and bbar->e+ (different b and bbar)
675705 // 2. b->c->e+ and bbar->cbar->e- (different b and bbar)
676706 // 3. b->c->e+ and b->e- (1 same b (or bbar))
677- if (!p1.has_mothers ())
707+ if (!p1.has_mothers ()) {
678708 return HFllType::kUndef ;
679- if (!p2.has_mothers ())
709+ }
710+ if (!p2.has_mothers ()) {
680711 return HFllType::kUndef ;
712+ }
681713
682714 int motherid_p1 = p1.mothersIds ()[0 ];
683715 int motherid_p2 = p2.mothersIds ()[0 ];
@@ -727,10 +759,12 @@ struct Alice3Dilepton {
727759 {
728760 // in total, 1 case for LS pairs
729761 // 4. b->c->e+ and bbar->e+
730- if (!p1.has_mothers ())
762+ if (!p1.has_mothers ()) {
731763 return HFllType::kUndef ;
732- if (!p2.has_mothers ())
764+ }
765+ if (!p2.has_mothers ()) {
733766 return HFllType::kUndef ;
767+ }
734768
735769 int motherid_p1 = p1.mothersIds ()[0 ];
736770 int motherid_p2 = p2.mothersIds ()[0 ];
@@ -788,13 +822,27 @@ struct Alice3Dilepton {
788822 }
789823
790824 template <bool isWithSmearing, PairType pairtype, typename TTracks>
791- void FillPairRecAll (TTracks const & tracks1, TTracks const & tracks2)
825+ void FillPairRecAll (TTracks const & tracks1, TTracks const & tracks2, const aod::McParticles& mcParticles )
792826 {
793827 if constexpr (pairtype == PairType::kULS ) {
794828 for (const auto & [t1, t2] : combinations (soa::CombinationsFullIndexPolicy (tracks1, tracks2))) {
795829 if (!IsInAcceptance<isWithSmearing>(t1) || !IsInAcceptance<isWithSmearing>(t2)) {
796830 continue ;
797831 }
832+ if (contamination) {
833+ if (t1.has_mcParticle () && t2.has_mcParticle ()) {
834+ auto mct1 = t1.template mcParticle_as <aod::McParticles>();
835+ auto mct2 = t2.template mcParticle_as <aod::McParticles>();
836+ if (std::abs (mct1.pdgCode ()) == pdg && std::abs (mct2.pdgCode ()) == pdg) {
837+ continue ;
838+ }
839+ }
840+ }
841+ if (requireHFEid) {
842+ if (!IsHF (t1, mcParticles) && !IsHF (t2, mcParticles)) {
843+ continue ;
844+ }
845+ }
798846 float pair_dca_xy = 999 .f ;
799847 ROOT ::Math::PtEtaPhiMVector v12 = buildPairDCA<isWithSmearing>(t1, t2, pair_dca_xy);
800848
@@ -811,6 +859,20 @@ struct Alice3Dilepton {
811859 if (!IsInAcceptance<isWithSmearing>(t1) || !IsInAcceptance<isWithSmearing>(t2)) {
812860 continue ;
813861 }
862+ if (contamination) {
863+ if (t1.has_mcParticle () && t2.has_mcParticle ()) {
864+ auto mct1 = t1.template mcParticle_as <aod::McParticles>();
865+ auto mct2 = t2.template mcParticle_as <aod::McParticles>();
866+ if (std::abs (mct1.pdgCode ()) == pdg && std::abs (mct2.pdgCode ()) == pdg) {
867+ continue ;
868+ }
869+ }
870+ }
871+ if (requireHFEid) {
872+ if (!IsHF (t1, mcParticles) && !IsHF (t2, mcParticles)) {
873+ continue ;
874+ }
875+ }
814876 float pair_dca_xy = 999 .f ;
815877 ROOT ::Math::PtEtaPhiMVector v12 = buildPairDCA<isWithSmearing>(t1, t2, pair_dca_xy);
816878
@@ -1084,7 +1146,7 @@ struct Alice3Dilepton {
10841146
10851147 void processRecAll (MyFilteredAlice3Collision const & collisions,
10861148 MyFilteredTracksMC const &,
1087- const aod::McParticles&)
1149+ const aod::McParticles& mcParticles )
10881150 {
10891151 for (const auto & collision : collisions) {
10901152 registry.fill (HIST (" Reconstructed/Event/VtxZ" ), collision.posZ ());
@@ -1125,9 +1187,9 @@ struct Alice3Dilepton {
11251187 registry.fill (HIST (" Reconstructed/Track/Pre" ), track.isTrackPrefilter ());
11261188 }
11271189
1128- FillPairRecAll<false , PairType::kULS >(negTracks_coll, posTracks_coll);
1129- FillPairRecAll<false , PairType::kLSpp >(posTracks_coll, posTracks_coll);
1130- FillPairRecAll<false , PairType::kLSnn >(negTracks_coll, negTracks_coll);
1190+ FillPairRecAll<false , PairType::kULS >(negTracks_coll, posTracks_coll, mcParticles );
1191+ FillPairRecAll<false , PairType::kLSpp >(posTracks_coll, posTracks_coll, mcParticles );
1192+ FillPairRecAll<false , PairType::kLSnn >(negTracks_coll, negTracks_coll, mcParticles );
11311193
11321194 } // end of collision loop
11331195 } // end of processRec
@@ -1150,7 +1212,7 @@ struct Alice3Dilepton {
11501212
11511213 void processRecAllWithSmearing (MyFilteredAlice3Collision const & collisions,
11521214 MyFilteredTracksWithSmearing const &,
1153- const aod::McParticles&)
1215+ const aod::McParticles& mcParticles )
11541216 {
11551217 for (const auto & collision : collisions) {
11561218 registry.fill (HIST (" Reconstructed/Event/VtxZ" ), collision.posZ ());
@@ -1191,9 +1253,9 @@ struct Alice3Dilepton {
11911253 registry.fill (HIST (" Reconstructed/Track/Pre" ), track.isTrackPrefilter ());
11921254 }
11931255
1194- FillPairRecAll<true , PairType::kULS >(negTracks_coll, posTracks_coll);
1195- FillPairRecAll<true , PairType::kLSpp >(posTracks_coll, posTracks_coll);
1196- FillPairRecAll<true , PairType::kLSnn >(negTracks_coll, negTracks_coll);
1256+ FillPairRecAll<true , PairType::kULS >(negTracks_coll, posTracks_coll, mcParticles );
1257+ FillPairRecAll<true , PairType::kLSpp >(posTracks_coll, posTracks_coll, mcParticles );
1258+ FillPairRecAll<true , PairType::kLSnn >(negTracks_coll, negTracks_coll, mcParticles );
11971259
11981260 } // end of collision loop
11991261 } // end of processRec
0 commit comments