88// In applying this license CERN does not waive the privileges and immunities
99// granted to it by virtue of its status as an Intergovernmental Organization
1010// or submit itself to any jurisdiction.
11- // / \author Andrea Rossi <andrea.rossi@cern.ch>
1211
12+ // / \file filterTracks.cxx
1313// / \brief Simple task to filter tracks and save infos to trees for DCA-related studies (alignment, HF-related issues, ...)
14+ // /
15+ // / \author Andrea Rossi <andrea.rossi@cern.ch>
1416
1517#include " Common/Core/RecoDecay.h"
1618#include " Common/Core/TrackSelection.h"
2527#include " Framework/AnalysisTask.h"
2628#include " Framework/runDataProcessing.h"
2729
30+ #include < TPDGCode.h>
31+
32+ #include < vector>
33+
2834using namespace o2 ;
2935using namespace o2 ::framework;
3036using namespace o2 ::aod::track;
@@ -33,7 +39,7 @@ using namespace o2::framework::expressions;
3339
3440namespace o2 ::aod
3541{
36- namespace filterTracks
42+ namespace filtertracks
3743{
3844
3945DECLARE_SOA_INDEX_COLUMN (Collision, collision); // ! Collision index
@@ -58,7 +64,7 @@ DECLARE_SOA_COLUMN(NsigmaTOFpr, nsigmaTOFpr, float); //! TOF nsigma w.r.t. proto
5864DECLARE_SOA_COLUMN (TpcNCluster, tpcNCluster, int ); // ! TOF nsigma w.r.t. proton mass hypothesis
5965
6066// /// MC INFO
61- DECLARE_SOA_COLUMN (MainHfMotherPdgCode, mainMotherPdgCode , int ); // ! mother pdg code for particles coming from HF, skipping intermediate resonance states. Not trustable when mother is not HF. Not suited for Sc->Lc decays, since Sc are never pointed to
67+ DECLARE_SOA_COLUMN (MainHfMotherPdgCode, mainHfMotherPdgCode , int ); // ! mother pdg code for particles coming from HF, skipping intermediate resonance states. Not trustable when mother is not HF. Not suited for Sc->Lc decays, since Sc are never pointed to
6268DECLARE_SOA_COLUMN (IsPhysicalPrimary, isPhysicalPrimary, bool ); // ! is phyiscal primary according to ALICE definition
6369DECLARE_SOA_COLUMN (MainBeautyAncestorPdgCode, mainBeautyAncestorPdgCode, int ); // ! pdgcode of beauty particle when this is an ancestor, otherwise -1
6470DECLARE_SOA_COLUMN (MainMotherOrigIndex, mainMotherOrigIndex, int ); // ! original index in MCParticle tree of main mother: needed when checking if particles come from same mother
@@ -69,7 +75,7 @@ DECLARE_SOA_COLUMN(MainMotherY, mainMotherY, float); //! origi
6975DECLARE_SOA_COLUMN (MainBeautyAncestorPt, mainBeautyAncestorPt, float ); // ! original index in MCParticle tree of main mother: needed when chekcing if particles come from same mother
7076DECLARE_SOA_COLUMN (MainBeautyAncestorY, mainBeautyAncestorY, float ); // ! original index in MCParticle tree of main mother: needed when chekcing if particles come from same mother
7177DECLARE_SOA_COLUMN (MaxEtaDaughter, maxEtaDaughter, float ); // ! max (abs) eta of daughter particles, needed to reproduce acceptance cut
72- } // namespace filterTracks
78+ } // namespace filtertracks
7379DECLARE_SOA_TABLE (FilterColl, " AOD" , " FILTERCOLL" ,
7480 o2::aod::collision::BCId,
7581 o2::aod::collision::PosX,
@@ -88,7 +94,7 @@ DECLARE_SOA_TABLE(FilterColl, "AOD", "FILTERCOLL",
8894 o2::aod::collision::CollisionTimeRes);
8995DECLARE_SOA_TABLE (FilterTrack, " AOD" , " FILTERTRACK" ,
9096 o2::aod::track::CollisionId,
91- aod::filterTracks ::IsInsideBeamPipe,
97+ aod::filtertracks ::IsInsideBeamPipe,
9298 o2::aod::track::TrackType,
9399 o2::aod::track::X,
94100 o2::aod::track::Alpha,
@@ -98,43 +104,43 @@ DECLARE_SOA_TABLE(FilterTrack, "AOD", "FILTERTRACK",
98104 o2::aod::track::Tgl,
99105 o2::aod::track::Signed1Pt);
100106DECLARE_SOA_TABLE (FilterTrackExtr, " AOD" , " FILTERTRACKEXTR" ,
101- // aod::filterTracks ::Px,aod::filterTracks ::Py, aod::filterTracks ::Pz,
102- aod::filterTracks ::Pt, o2::aod::track::Eta,
103- o2::aod::filterTracks ::Charge,
107+ // aod::filtertracks ::Px,aod::filtertracks ::Py, aod::filtertracks ::Pz,
108+ aod::filtertracks ::Pt, o2::aod::track::Eta,
109+ o2::aod::filtertracks ::Charge,
104110 o2::aod::track::DcaXY,
105111 o2::aod::track::DcaZ,
106112 o2::aod::track::SigmaDcaXY2,
107113 o2::aod::track::SigmaDcaZ2,
108- aod::filterTracks ::NsigmaTPCpi, aod::filterTracks ::NsigmaTPCka, aod::filterTracks ::NsigmaTPCpr,
109- aod::filterTracks ::NsigmaTOFpi, aod::filterTracks ::NsigmaTOFka, aod::filterTracks ::NsigmaTOFpr);
114+ aod::filtertracks ::NsigmaTPCpi, aod::filtertracks ::NsigmaTPCka, aod::filtertracks ::NsigmaTPCpr,
115+ aod::filtertracks ::NsigmaTOFpi, aod::filtertracks ::NsigmaTOFka, aod::filtertracks ::NsigmaTOFpr);
110116DECLARE_SOA_TABLE (FiltTracExtDet, " AOD" , " FILTTRACEXTDET" ,
111117 o2::aod::track::ITSClusterSizes,
112118 o2::aod::track::ITSChi2NCl,
113119 o2::aod::track::TPCChi2NCl,
114- aod::filterTracks ::TpcNCluster,
120+ aod::filtertracks ::TpcNCluster,
115121 o2::aod::track::TrackTime);
116122DECLARE_SOA_TABLE (FilterTrackMC, " AOD" , " FILTERTRACKMC" ,
117- // aod::filterTracks ::Px,aod::filterTracks ::Py, aod::filterTracks ::Pz,
123+ // aod::filtertracks ::Px,aod::filtertracks ::Py, aod::filtertracks ::Pz,
118124 o2::aod::mcparticle::PdgCode,
119- o2::aod::filterTracks ::IsPhysicalPrimary,
120- o2::aod::filterTracks ::MainHfMotherPdgCode,
121- o2::aod::filterTracks ::MainBeautyAncestorPdgCode,
122- o2::aod::filterTracks ::MainMotherOrigIndex,
123- o2::aod::filterTracks ::MainMotherNfinalStateDaught,
124- o2::aod::filterTracks ::MainMotherPt,
125- o2::aod::filterTracks ::MainMotherY,
126- o2::aod::filterTracks ::MainBeautyAncestorPt,
127- o2::aod::filterTracks ::MainBeautyAncestorY);
125+ o2::aod::filtertracks ::IsPhysicalPrimary,
126+ o2::aod::filtertracks ::MainHfMotherPdgCode,
127+ o2::aod::filtertracks ::MainBeautyAncestorPdgCode,
128+ o2::aod::filtertracks ::MainMotherOrigIndex,
129+ o2::aod::filtertracks ::MainMotherNfinalStateDaught,
130+ o2::aod::filtertracks ::MainMotherPt,
131+ o2::aod::filtertracks ::MainMotherY,
132+ o2::aod::filtertracks ::MainBeautyAncestorPt,
133+ o2::aod::filtertracks ::MainBeautyAncestorY);
128134DECLARE_SOA_TABLE (GenParticles, " AOD" , " GENPARTICLES" ,
129- // aod::filterTracks ::Px,aod::filterTracks ::Py, aod::filterTracks ::Pz,
135+ // aod::filtertracks ::Px,aod::filtertracks ::Py, aod::filtertracks ::Pz,
130136 o2::aod::mcparticle::PdgCode,
131137 o2::aod::mcparticle::McCollisionId,
132- o2::aod::filterTracks ::MainBeautyAncestorPdgCode,
133- o2::aod::filterTracks ::MainMotherPt,
134- o2::aod::filterTracks ::MainMotherY,
135- o2::aod::filterTracks ::MaxEtaDaughter,
136- o2::aod::filterTracks ::MainBeautyAncestorPt,
137- o2::aod::filterTracks ::MainBeautyAncestorY);
138+ o2::aod::filtertracks ::MainBeautyAncestorPdgCode,
139+ o2::aod::filtertracks ::MainMotherPt,
140+ o2::aod::filtertracks ::MainMotherY,
141+ o2::aod::filtertracks ::MaxEtaDaughter,
142+ o2::aod::filtertracks ::MainBeautyAncestorPt,
143+ o2::aod::filtertracks ::MainBeautyAncestorY);
138144} // namespace o2::aod
139145
140146struct FilterTracks {
@@ -162,18 +168,22 @@ struct FilterTracks {
162168 using TracksWithSelAndDcaMc = soa::Join<TracksWithSelAndDca, aod::McTrackLabels>;
163169 using FilterCollisionsWithEvSel = soa::Filtered<CollisionsWithEvSel>;
164170
165- Partition<soa::Filtered<TracksWithSelAndDca>> lowPtTracks = aod::track::pt < 2 .f && (nabs(aod::track::pt * 10000 .f - nround(aod::track::pt * 10000 .f)) < trackPtWeightLowPt.node() * 2 .f);
166- Partition<soa::Filtered<TracksWithSelAndDca>> midPtTracks = aod::track::pt > 2 .f && aod::track::pt < 5 .f && (nabs(aod::track::pt * 10000 .f - nround(aod::track::pt * 10000 .f)) < trackPtWeightMidPt.node() * 2 .f);
167- Partition<soa::Filtered<TracksWithSelAndDca>> highPtTracks = aod::track::pt > 5 .f;
171+ float lowPtThreshold = 2 .;
172+ float midPtThreshold = 5 .;
173+ float nDigitScaleFactor = 10000 .;
174+ Partition<soa::Filtered<TracksWithSelAndDca>> lowPtTracks = aod::track::pt < lowPtThreshold && (nabs(aod::track::pt * nDigitScaleFactor - nround(aod::track::pt * nDigitScaleFactor)) < trackPtWeightLowPt.node() * lowPtThreshold);
175+ Partition<soa::Filtered<TracksWithSelAndDca>> midPtTracks = aod::track::pt > lowPtThreshold&& aod::track::pt < midPtThreshold && (nabs(aod::track::pt * nDigitScaleFactor - nround(aod::track::pt * nDigitScaleFactor)) < trackPtWeightMidPt.node() * lowPtThreshold);
176+ Partition<soa::Filtered<TracksWithSelAndDca>> highPtTracks = aod::track::pt > midPtThreshold;
168177
169- Partition<soa::Filtered<TracksWithSelAndDcaMc>> lowPtTracksMC = aod::track::pt < 2 .f && (nabs(aod::track::pt * 10000 .f - nround(aod::track::pt * 10000 .f )) < trackPtWeightLowPt.node() * 2 .f );
170- Partition<soa::Filtered<TracksWithSelAndDcaMc>> midPtTracksMC = aod::track::pt > 2 .f && aod::track::pt < 5 .f && (nabs(aod::track::pt * 10000 .f - nround(aod::track::pt * 10000 .f )) < trackPtWeightMidPt.node() * 2 .f );
171- Partition<soa::Filtered<TracksWithSelAndDcaMc>> highPtTracksMC = aod::track::pt > 5 .f ;
178+ Partition<soa::Filtered<TracksWithSelAndDcaMc>> lowPtTracksMC = aod::track::pt < lowPtThreshold && (nabs(aod::track::pt * nDigitScaleFactor - nround(aod::track::pt * nDigitScaleFactor )) < trackPtWeightLowPt.node() * lowPtThreshold );
179+ Partition<soa::Filtered<TracksWithSelAndDcaMc>> midPtTracksMC = aod::track::pt > lowPtThreshold && aod::track::pt < midPtThreshold && (nabs(aod::track::pt * nDigitScaleFactor - nround(aod::track::pt * nDigitScaleFactor )) < trackPtWeightMidPt.node() * lowPtThreshold );
180+ Partition<soa::Filtered<TracksWithSelAndDcaMc>> highPtTracksMC = aod::track::pt > midPtThreshold ;
172181
173- std::array<int , 3 > pdgSignalParticleArray = {310 , 421 , 4122 }; // K0s, D0 and Lc
174- std::array<int , 3 > pdgDecayLc = {2212 , -321 , 211 };
175- std::array<int , 2 > pdgDecayDzero = {-321 , 211 };
176- std::array<int , 2 > pdgDecayKzero = {-211 , 211 };
182+ const int nStudiedParticlesMc = 3 ;
183+ std::array<int , 3 > pdgSignalParticleArray = {kK0Short , o2::constants::physics::Pdg::kD0 , o2::constants::physics::Pdg::kLambdaCPlus }; // K0s, D0 and Lc
184+ std::array<int , 3 > pdgDecayLc = {kProton , kKMinus , kPiPlus };
185+ std::array<int , 2 > pdgDecayDzero = {kKMinus , kPiPlus };
186+ std::array<int , 2 > pdgDecayKzero = {kPiMinus , kPiPlus };
177187
178188 void init (InitContext&)
179189 {
@@ -192,23 +202,23 @@ struct FilterTracks {
192202 {
193203
194204 fillTableData (track);
195- bool has_MCparticle = track.has_mcParticle ();
196- if (has_MCparticle ) {
205+ bool has_McParticle = track.has_mcParticle ();
206+ if (has_McParticle ) {
197207 // / the track is not fake
198208
199209 // check whether the particle comes from a charm or beauty hadron and store its index
200210
201211 auto mcparticle = track.mcParticle ();
202212 int pdgParticleMother = 0 ;
203- for (int iSignPart = 0 ; iSignPart < 3 ; iSignPart++) {
213+ for (int iSignPart = 0 ; iSignPart < nStudiedParticlesMc ; iSignPart++) {
204214 pdgParticleMother = pdgSignalParticleArray[iSignPart];
205215 auto motherIndex = RecoDecay::getMother (mcParticles, mcparticle, pdgParticleMother, true ); // check whether mcparticle derives from a particle with pdg = pdgparticlemother, accepting also antiparticle (<- the true parameter)
206216 if (motherIndex != -1 ) {
207217 auto particleMother = mcParticles.rawIteratorAt (motherIndex);
208218 // just for internal check
209219 // double mass=particleMother.e()*particleMother.e()-particleMother.pt()*particleMother.pt()-particleMother.pz()*particleMother.pz();
210220 // filteredTracksMC(mcparticle.pdgCode(),mcparticle.isPhysicalPrimary(),particleMother.pdgCode(),0,motherIndex,0,particleMother.pt(),particleMother.y(),std::sqrt(mass),0);
211- if (pdgParticleMother == 310 ) {
221+ if (pdgParticleMother == kK0Short ) {
212222 auto daughtersSlice = mcparticle.template daughters_as <aod::McParticles>();
213223 int ndaught = daughtersSlice.size (); // might not be accurate in case K0s interact with material before decaying
214224 if (ndaught != 2 )
@@ -220,24 +230,26 @@ struct FilterTracks {
220230
221231 int ndaught = 0 ;
222232 std::vector<int > indxDaughers;
223- if (pdgParticleMother == 421 ) {
233+ if (pdgParticleMother == o2::constants::physics::Pdg:: kD0 ) {
224234 if (RecoDecay::isMatchedMCGen<true , false >(mcParticles, particleMother, pdgParticleMother, pdgDecayDzero, true , nullptr , 3 , &indxDaughers)) {
225235 ndaught = 2 ;
226236 // std::cout<<"######## FOUND D0, MATCHED! pdg: " <<particleMother.pdgCode()<<"################ size array "<<indxDaughers.size()<<std::endl;
227- } else
237+ } else {
228238 ndaught = -indxDaughers.size ();
229- } else if (pdgParticleMother == 4122 ) {
239+ }
240+ } else if (pdgParticleMother == o2::constants::physics::Pdg::kLambdaCPlus ) {
230241 if (RecoDecay::isMatchedMCGen<true , false >(mcParticles, particleMother, pdgParticleMother, pdgDecayLc, true , nullptr , 3 , &indxDaughers)) {
231242 ndaught = 3 ;
232- } else
243+ } else {
233244 ndaught = -indxDaughers.size ();
245+ }
234246 }
235247 // now check whether the charm hadron is prompt or comes from beauty decay
236248 std::vector<int > idxBhadMothers;
237249 if (RecoDecay::getCharmHadronOrigin (mcParticles, particleMother, false , &idxBhadMothers) == RecoDecay::OriginType::NonPrompt) {
238250 if (idxBhadMothers.size () > 1 ) {
239251 LOG (info) << " more than 1 B mother hadron found, should not be: " ;
240- for (unsigned long iBhM = 0 ; iBhM < idxBhadMothers.size (); iBhM++) {
252+ for (int64_t iBhM = 0 ; iBhM < idxBhadMothers.size (); iBhM++) {
241253 auto particleBhadr = mcParticles.rawIteratorAt (idxBhadMothers[iBhM]);
242254 LOG (info) << particleBhadr.pdgCode ();
243255 }
@@ -262,20 +274,20 @@ struct FilterTracks {
262274 void processData (FilterCollisionsWithEvSel::iterator const & collision, soa::Filtered<TracksWithSelAndDca> const & tracks)
263275 {
264276 if (trackPtSampling == 0 ) {
265- for (auto & track : tracks) {
277+ for (auto const & track : tracks) {
266278 fillTableData (track);
267279 }
268280 } else {
269281 auto lowPtTracksThisColl = lowPtTracks->sliceByCached (aod::track::collisionId, collision.globalIndex (), cache);
270- for (auto & track : lowPtTracksThisColl) {
282+ for (auto const & track : lowPtTracksThisColl) {
271283 fillTableData (track);
272284 }
273285 auto midPtTracksThisColl = midPtTracks->sliceByCached (aod::track::collisionId, collision.globalIndex (), cache);
274- for (auto & track : midPtTracksThisColl) {
286+ for (auto const & track : midPtTracksThisColl) {
275287 fillTableData (track);
276288 }
277289 auto highPtTracksThisColl = highPtTracks->sliceByCached (aod::track::collisionId, collision.globalIndex (), cache);
278- for (auto & track : highPtTracksThisColl) {
290+ for (auto const & track : highPtTracksThisColl) {
279291 fillTableData (track);
280292 }
281293 }
@@ -290,38 +302,38 @@ struct FilterTracks {
290302 void processMC (FilterCollisionsWithEvSel::iterator const & collision, soa::Filtered<TracksWithSelAndDcaMc> const & tracks, aod::McParticles const & mcParticles)
291303 {
292304 if (trackPtSampling == 0 ) {
293- for (auto & track : tracks) {
305+ for (auto const & track : tracks) {
294306 fillTableDataMC (track, mcParticles);
295307 }
296308 } else {
297309 auto lowPtTracksMCThisColl = lowPtTracksMC->sliceByCached (aod::track::collisionId, collision.globalIndex (), cache);
298- for (auto & track : lowPtTracksMCThisColl) {
310+ for (auto const & track : lowPtTracksMCThisColl) {
299311 fillTableDataMC (track, mcParticles);
300312 }
301313 auto midPtTracksMCThisColl = midPtTracksMC->sliceByCached (aod::track::collisionId, collision.globalIndex (), cache);
302- for (auto & track : midPtTracksMCThisColl) {
314+ for (auto const & track : midPtTracksMCThisColl) {
303315 fillTableDataMC (track, mcParticles);
304316 }
305317 auto highPtTracksMCThisColl = highPtTracksMC->sliceByCached (aod::track::collisionId, collision.globalIndex (), cache);
306- for (auto & track : highPtTracksMCThisColl) {
318+ for (auto const & track : highPtTracksMCThisColl) {
307319 fillTableDataMC (track, mcParticles);
308320 }
309321 }
310322
311- for (auto & mcpart : mcParticles) { // NOTE THAT OF COURSE IN CASE OF SAMPLING THE GEN TABLE WON'T MATCH THE RECO EVEN CONSIDERING EFFICIENCY
323+ for (auto const & mcpart : mcParticles) { // NOTE THAT OF COURSE IN CASE OF SAMPLING THE GEN TABLE WON'T MATCH THE RECO EVEN CONSIDERING EFFICIENCY
312324 int pdgCode = mcpart.pdgCode ();
313325 // for(int iSignPart=0;iSignPart<3;iSignPart++){
314326
315327 std::vector<int > indxDaughers;
316328 float etamax = 0 ;
317329 bool isMatchedToSignal = false ;
318- if (std::abs (pdgCode) == 310 ) {
319- isMatchedToSignal = RecoDecay::isMatchedMCGen<true , false >(mcParticles, mcpart, 310 , pdgDecayKzero, true , nullptr , 1 , &indxDaughers);
330+ if (std::abs (pdgCode) == kK0Short ) {
331+ isMatchedToSignal = RecoDecay::isMatchedMCGen<true , false >(mcParticles, mcpart, kK0Short , pdgDecayKzero, true , nullptr , 1 , &indxDaughers);
320332 }
321- if (std::abs (pdgCode) == 421 ) {
322- isMatchedToSignal = RecoDecay::isMatchedMCGen<true , false >(mcParticles, mcpart, 421 , pdgDecayDzero, true , nullptr , 3 , &indxDaughers);
323- } else if (std::abs (pdgCode) == 4122 ) {
324- isMatchedToSignal = RecoDecay::isMatchedMCGen<true , false >(mcParticles, mcpart, 4122 , pdgDecayLc, true , nullptr , 3 , &indxDaughers);
333+ if (std::abs (pdgCode) == o2::constants::physics::Pdg:: kD0 ) {
334+ isMatchedToSignal = RecoDecay::isMatchedMCGen<true , false >(mcParticles, mcpart, o2::constants::physics::Pdg:: kD0 , pdgDecayDzero, true , nullptr , 3 , &indxDaughers);
335+ } else if (std::abs (pdgCode) == o2::constants::physics::Pdg:: kLambdaCPlus ) {
336+ isMatchedToSignal = RecoDecay::isMatchedMCGen<true , false >(mcParticles, mcpart, o2::constants::physics::Pdg:: kLambdaCPlus , pdgDecayLc, true , nullptr , 3 , &indxDaughers);
325337 // std::cout<<"Lc found, matched to MC? "<<isMatchedToSignal<<std::endl;
326338 // if(!isMatchedToSignal){
327339 // auto daughtersLxSlice = mcpart.daughters_as<aod::McParticles>();
@@ -339,24 +351,25 @@ struct FilterTracks {
339351 etamax = eta;
340352 }
341353 }
342- if (pdgCode == 310 ) {
354+ if (pdgCode == kK0Short ) {
343355 selectedGenParticles (mcpart.pdgCode (), mcpart.mcCollisionId (), 0 , mcpart.pt (), mcpart.y (), etamax, 0 , 0 );
344356 continue ;
345357 }
346358 std::vector<int > idxBhadMothers;
347359 if (RecoDecay::getCharmHadronOrigin (mcParticles, mcpart, false , &idxBhadMothers) == RecoDecay::OriginType::NonPrompt) {
348360 if (idxBhadMothers.size () > 1 ) {
349361 LOG (info) << " loop on gen particles: more than 1 B mother hadron found, should not be: " ;
350- for (unsigned long iBhM = 0 ; iBhM < idxBhadMothers.size (); iBhM++) {
362+ for (int64_t iBhM = 0 ; iBhM < idxBhadMothers.size (); iBhM++) {
351363 auto particleBhadr = mcParticles.rawIteratorAt (idxBhadMothers[iBhM]);
352364 LOG (info) << particleBhadr.pdgCode ();
353365 }
354366 }
355367 auto particleBhadr = mcParticles.rawIteratorAt (idxBhadMothers[0 ]);
356368 // int pdgBhad=particleBhadr.pdgCode();
357369 selectedGenParticles (mcpart.pdgCode (), mcpart.mcCollisionId (), particleBhadr.pdgCode (), mcpart.pt (), mcpart.y (), etamax, particleBhadr.pt (), particleBhadr.y ());
358- } else
370+ } else {
359371 selectedGenParticles (mcpart.pdgCode (), mcpart.mcCollisionId (), 0 , mcpart.pt (), mcpart.y (), etamax, 0 , 0 );
372+ }
360373 }
361374 //
362375 }
0 commit comments