6161#include < memory>
6262#include < stdexcept>
6363#include < string>
64+ #include < unordered_set>
6465#include < vector>
6566
6667using namespace o2 ;
@@ -156,6 +157,7 @@ struct nucleiQC {
156157 std::shared_ptr<TH1> mHistTrackTunedTracks = mHistograms .get<TH1>(HIST(" hTrackTunedTracks" ));
157158
158159 std::vector<int > mSpeciesToProcess ;
160+ std::array<bool , nuclei::Species::kNspecies > mFillSpecies {false };
159161 Produces<aod::NucleiTableRed> mNucleiTableRed ;
160162 Produces<aod::NucleiTableExt> mNucleiTableExt ;
161163
@@ -169,6 +171,7 @@ struct nucleiQC {
169171 o2::dataformats::VertexBase mVtx ;
170172 o2::track::TrackParametrizationWithError<float > mTrackParCov ;
171173 std::array<nuclei::PidManager, static_cast <int >(nuclei::Species::kNspecies )> mPidManagers ;
174+ bool mAnyTrackTuner = false ;
172175
173176 void init (o2::framework::InitContext&)
174177 {
@@ -177,22 +180,24 @@ struct nucleiQC {
177180 mCcdb ->setCaching (true );
178181 mCcdb ->setLocalObjectValidityChecking ();
179182 mCcdb ->setFatalWhenNull (false );
180- nuclei::lut = o2::base::MatLayerCylSet::rectifyPtrFromFile (mCcdb ->get <o2::base::MatLayerCylSet>(" GLO/Param/MatLUT" ));
183+ // nuclei::lut = o2::base::MatLayerCylSet::rectifyPtrFromFile(mCcdb->get<o2::base::MatLayerCylSet>("GLO/Param/MatLUT"));
181184
182185 for (int iSel = 0 ; iSel < nuclei::evSel::kNevSels ; iSel++) {
183186 mHistograms .get <TH1>(HIST (" hEventSelections" ))->GetXaxis ()->SetBinLabel (iSel + 1 , nuclei::eventSelectionLabels[iSel].c_str ());
184187 }
185188
186189 for (int iSpecies = 0 ; iSpecies < static_cast <int >(nuclei::Species::kNspecies ); iSpecies++) {
187- if (cfgSpeciesToProcess->get (iSpecies) == 1 )
190+ if (cfgSpeciesToProcess->get (iSpecies) == 1 ) {
188191 mSpeciesToProcess .emplace_back (iSpecies);
192+ mFillSpecies [iSpecies] = true ;
193+ }
189194 }
190195
191196 static_for<0 , nuclei::kNspecies - 1 >([&](auto iSpecies) {
192197 constexpr int kSpeciesCt = decltype (iSpecies)::value;
193198 const int kSpeciesRt = kSpeciesCt ;
194199
195- if (std::find ( mSpeciesToProcess . begin (), mSpeciesToProcess . end (), kSpeciesCt ) == mSpeciesToProcess . end () )
200+ if (! mFillSpecies [ kSpeciesCt ] )
196201 return ;
197202
198203 float tpcBetheBlochParams[6 ];
@@ -214,12 +219,12 @@ struct nucleiQC {
214219 });
215220
216221 // / TrackTuner initialization
217- bool anyTrackTuner = false ;
222+ mAnyTrackTuner = false ;
218223 for (int iSpecies = 0 ; iSpecies < static_cast <int >(nuclei::Species::kNspecies ); iSpecies++) {
219- anyTrackTuner = anyTrackTuner || cfgUseTrackTuner->get (iSpecies);
224+ mAnyTrackTuner = mAnyTrackTuner || cfgUseTrackTuner->get (iSpecies);
220225 }
221226
222- if (anyTrackTuner ) {
227+ if (mAnyTrackTuner ) {
223228 std::string outputStringParams = " " ;
224229 switch (cfgTrackTunerConfigSource) {
225230 case aod::track_tuner::InputString:
@@ -251,13 +256,17 @@ struct nucleiQC {
251256
252257 o2::parameters::GRPMagField* grpmag = mCcdb ->getForTimeStamp <o2::parameters::GRPMagField>(" GLO/Config/GRPMagField" , timestamp);
253258 o2::base::Propagator::initFieldFromGRP (grpmag);
254- o2::base::Propagator::Instance ()->setMatLUT (nuclei::lut);
259+ // o2::base::Propagator::Instance()->setMatLUT(nuclei::lut);
260+ if (!o2::base::Propagator::Instance ()->getMatLUT ()) {
261+ nuclei::lut = o2::base::MatLayerCylSet::rectifyPtrFromFile (mCcdb ->get <o2::base::MatLayerCylSet>(" GLO/Param/MatLUT" ));
262+ o2::base::Propagator::Instance ()->setMatLUT (nuclei::lut);
263+ }
255264 mBz = static_cast <float >(grpmag->getNominalL3Field ());
256265 LOGF (info, " Retrieved GRP for timestamp %ull (%i) with magnetic field of %1.2f kZG" , timestamp, mRunNumber , mBz );
257266 }
258267
259- template <int iSpecies, typename Ttrack , typename Tcollision >
260- bool trackSelection (const Ttrack& track, const Tcollision& collision )
268+ template <int iSpecies, bool isMc , typename Ttrack >
269+ bool trackSelection (const Ttrack& track)
261270 {
262271 mHistograms .fill (HIST (nuclei::cNames[iSpecies]) + HIST (" /hTrackQuality" ), track.sign () * track.pt (), trackQuality::kNoCuts );
263272
@@ -289,13 +298,24 @@ struct nucleiQC {
289298 return false ;
290299 mHistograms .fill (HIST (nuclei::cNames[iSpecies]) + HIST (" /hTrackQuality" ), track.sign () * track.pt (), trackQuality::kItsChi2Cut );
291300
292- const o2::math_utils::Point3D<float > collisionVertex{collision.posX (), collision.posY (), collision.posZ ()};
293301 mDcaInfoCov .set (999 , 999 , 999 , 999 , 999 );
294302 setTrackParCov (track, mTrackParCov );
295303 mTrackParCov .setPID (track.pidForTracking ());
296- std::array<float , 2 > dcaInfo;
297- o2::base::Propagator::Instance ()->propagateToDCA (collisionVertex, mTrackParCov , mBz , 2 .f , static_cast <o2::base::Propagator::MatCorrType>(cfgMaterialCorrection.value ), &dcaInfo);
298- if (std::abs (dcaInfo[0 ]) > cfgCutDCAxy || std::abs (dcaInfo[1 ]) > cfgCutDCAz)
304+
305+ if constexpr (isMc) {
306+ if (track.has_mcParticle () && cfgUseTrackTuner->get (iSpecies)) {
307+ const auto & particle = track.mcParticle ();
308+ mHistTrackTunedTracks ->Fill (1 .);
309+ mTrackTuner .tuneTrackParams (particle, mTrackParCov , mMatCorr , &mDcaInfoCov , mHistTrackTunedTracks );
310+ }
311+ } else {
312+ mMatCorr = static_cast <o2::base::Propagator::MatCorrType>(cfgMaterialCorrection.value );
313+ }
314+
315+ o2::base::Propagator::Instance ()->propagateToDCABxByBz (mVtx , mTrackParCov , 2 .f , mMatCorr , &mDcaInfoCov );
316+ mDcaInfo [0 ] = mDcaInfoCov .getY ();
317+ mDcaInfo [1 ] = mDcaInfoCov .getZ ();
318+ if (std::abs (mDcaInfo [0 ]) > cfgCutDCAxy || std::abs (mDcaInfo [1 ]) > cfgCutDCAz)
299319 return false ;
300320 mHistograms .fill (HIST (nuclei::cNames[iSpecies]) + HIST (" /hTrackQuality" ), track.sign () * track.pt (), trackQuality::kDcaCuts );
301321
@@ -365,9 +385,9 @@ struct nucleiQC {
365385 }
366386
367387 template <typename Tparticle>
368- void fillCollisionFlag (const Tparticle& particle, nuclei::SlimCandidate& candidate, const std::vector< bool >& reconstructedCollision)
388+ void fillCollisionFlag (const Tparticle& particle, nuclei::SlimCandidate& candidate, const std::unordered_set< int >& reconstructedCollision)
369389 {
370- if (reconstructedCollision[ particle.mcCollisionId ()] ) {
390+ if (reconstructedCollision. count ( particle.mcCollisionId ()) > 0 ) {
371391 candidate.flags |= nuclei::QcFlags::kQcHasReconstructedCollision ;
372392 }
373393 }
@@ -396,31 +416,6 @@ struct nucleiQC {
396416 candidate.phiGenerated = particle.phi ();
397417 }
398418
399- template <const bool isMc, typename Tcollision, typename Ttrack>
400- void fillDcaInformation (const int iSpecies, const Tcollision& collision, const Ttrack& track, nuclei::SlimCandidate& candidate, const aod::McParticles::iterator& particle)
401- {
402-
403- mDcaInfoCov .set (999 , 999 , 999 , 999 , 999 );
404- setTrackParCov (track, mTrackParCov );
405- mTrackParCov .setPID (track.pidForTracking ());
406-
407- if constexpr (isMc) {
408- if (track.has_mcParticle () && cfgUseTrackTuner->get (iSpecies)) {
409- mHistTrackTunedTracks ->Fill (1 .);
410- mTrackTuner .tuneTrackParams (particle, mTrackParCov , mMatCorr , &mDcaInfoCov , mHistTrackTunedTracks );
411- }
412- } else {
413- mMatCorr = static_cast <o2::base::Propagator::MatCorrType>(cfgMaterialCorrection.value );
414- }
415-
416- mVtx .setPos ({collision.posX (), collision.posY (), collision.posZ ()});
417- mVtx .setCov (collision.covXX (), collision.covXY (), collision.covYY (), collision.covXZ (), collision.covYZ (), collision.covZZ ());
418- o2::base::Propagator::Instance ()->propagateToDCABxByBz (mVtx , mTrackParCov , 2 .f , mMatCorr , &mDcaInfoCov );
419-
420- candidate.DCAxy = mDcaInfoCov .getY ();
421- candidate.DCAz = mDcaInfoCov .getZ ();
422- }
423-
424419 template <const bool isMc, typename Tcollision, typename Ttrack>
425420 nuclei::SlimCandidate fillCandidate (const int iSpecies, Tcollision const & collision, Ttrack const & track)
426421 {
@@ -434,8 +429,8 @@ struct nucleiQC {
434429 .clusterSizesITS = track.itsClusterSizes (),
435430 .TPCsignal = track.tpcSignal (),
436431 .beta = mPidManagers [iSpecies].getBetaTOF (track),
437- .DCAxy = 0 . f ,
438- .DCAz = 0 . f ,
432+ .DCAxy = mDcaInfo [ 0 ], // set in the track selection function
433+ .DCAz = mDcaInfo [ 1 ], // set in the track selection function
439434 .flags = 0 ,
440435 .pdgCode = 0 ,
441436 .motherPdgCode = 0 ,
@@ -450,19 +445,15 @@ struct nucleiQC {
450445
451446 fillNucleusFlagsPdgs (collision, track, candidate);
452447
453- aod::McParticles::iterator particle;
454-
455448 if constexpr (isMc) {
456449 if (track.has_mcParticle ()) {
457450
458- particle = track.mcParticle ();
451+ const auto & particle = track.mcParticle ();
459452 fillNucleusFlagsPdgsMc (particle, candidate);
460453 fillNucleusGeneratedVariables (particle, candidate);
461454 }
462455 }
463456
464- fillDcaInformation<isMc>(iSpecies, collision, track, candidate, particle);
465-
466457 return candidate;
467458 }
468459
@@ -510,12 +501,36 @@ struct nucleiQC {
510501 }
511502 }
512503
513- void processMc (const Collisions& collisions, const TrackCandidatesMC& tracks, const aod::BCsWithTimestamps&, const aod::McParticles& mcParticles, const aod::McCollisions& mcCollisions)
504+ void writeCandidate (const nuclei::SlimCandidate& candidate)
505+ {
506+ if (!cfgFillTable)
507+ return ;
508+
509+ mNucleiTableRed (
510+ candidate.pt ,
511+ candidate.eta ,
512+ candidate.phi ,
513+ candidate.tpcInnerParam ,
514+ candidate.clusterSizesITS ,
515+ candidate.TPCsignal ,
516+ candidate.beta ,
517+ candidate.DCAxy ,
518+ candidate.DCAz ,
519+ candidate.flags ,
520+ candidate.ptGenerated ,
521+ candidate.mcProcess ,
522+ candidate.pdgCode ,
523+ candidate.motherPdgCode );
524+ mNucleiTableExt (
525+ candidate.nsigmaTpc ,
526+ candidate.nsigmaTof );
527+ }
528+
529+ void processMc (const Collisions& collisions, const TrackCandidatesMC& tracks, const aod::BCsWithTimestamps&, const aod::McParticles& mcParticles, const aod::McCollisions& /* mcCollisions*/ )
514530 {
515531 gRandom ->SetSeed (67 );
516- mNucleiCandidates .clear ();
517- std::vector<bool > reconstructedMcParticles (mcParticles.size (), false );
518- std::vector<bool > reconstructedCollisions (mcCollisions.size (), false );
532+ std::unordered_set<int > reconstructedMcParticles;
533+ std::unordered_set<int > reconstructedCollisions;
519534
520535 for (const auto & collision : collisions) {
521536
@@ -525,13 +540,11 @@ struct nucleiQC {
525540 if (!nuclei::eventSelection (collision, mHistograms , cfgEventSelections, cfgCutVertex))
526541 continue ;
527542 mHistograms .fill (HIST (" hCentrality" ), nuclei::getCentrality (collision, cfgCentralityEstimator, mHistFailCentrality ));
528- reconstructedCollisions[collision.mcCollisionId ()] = true ;
543+ reconstructedCollisions.insert (collision.mcCollisionId ());
544+ mVtx .setPos ({collision.posX (), collision.posY (), collision.posZ ()});
545+ mVtx .setCov (collision.covXX (), collision.covXY (), collision.covYY (), collision.covXZ (), collision.covYZ (), collision.covZZ ());
529546
530- bool anyTrackTuner = false ;
531- for (int iSpecies = 0 ; iSpecies < static_cast <int >(nuclei::Species::kNspecies ); iSpecies++) {
532- anyTrackTuner = anyTrackTuner || cfgUseTrackTuner->get (iSpecies);
533- }
534- if (anyTrackTuner && mTrackTuner .autoDetectDcaCalib && !mTrackTuner .areGraphsConfigured ) {
547+ if (mAnyTrackTuner && mTrackTuner .autoDetectDcaCalib && !mTrackTuner .areGraphsConfigured ) {
535548
536549 mTrackTuner .setRunNumber (mRunNumber );
537550
@@ -544,38 +557,41 @@ struct nucleiQC {
544557 auto tracksThisCollision = tracks.sliceBy (mTracksPerCollision , collision.globalIndex ());
545558 for (const auto & track : tracksThisCollision) {
546559
547- static_for< 0 , nuclei:: kNspecies - 1 >([&]( auto iSpecies) {
548- constexpr int kSpeciesCt = decltype (iSpecies)::value;
549- const int kSpeciesRt = kSpeciesCt ;
560+ // selections shared among all species
561+ if (!track. has_mcParticle ())
562+ continue ;
550563
551- if (std::find (mSpeciesToProcess .begin (), mSpeciesToProcess .end (), kSpeciesRt ) == mSpeciesToProcess .end ())
552- return ;
564+ if (track.mcParticleId () < -1 || track.mcParticleId () >= mcParticles.size ())
565+ continue ;
566+ const auto & particle = mcParticles.iteratorAt (track.mcParticleId ());
553567
554- if (!track.has_mcParticle ())
555- return ;
568+ if (cfgRapidityToggle && ((particle.y () - cfgRapidityCenterMass) < cfgRapidityMin || (particle.y () - cfgRapidityCenterMass) > cfgRapidityMax))
569+ continue ;
570+
571+ if (cfgFillOnlyPhysicalPrimaries && !particle.isPhysicalPrimary ())
572+ continue ;
556573
557- if (track.mcParticleId () < -1 || track.mcParticleId () >= mcParticles.size ())
574+ // species-specific selections and filling
575+ static_for<0 , nuclei::kNspecies - 1 >([&](auto iSpecies) {
576+ constexpr int kSpeciesCt = decltype (iSpecies)::value;
577+ const int kSpeciesRt = kSpeciesCt ;
578+ if (!mFillSpecies [kSpeciesCt ])
558579 return ;
559- const auto & particle = mcParticles.iteratorAt (track.mcParticleId ());
560580
561581 if (cfgDoCheckPdgCode) {
562582 if (std::abs (particle.pdgCode ()) != nuclei::pdgCodes[kSpeciesRt ])
563583 return ;
564584 }
565585
586+ mDcaInfo = {-999 .f , -999 .f };
587+
566588 if (cfgDownscalingFactor->get (kSpeciesRt ) < 1 .) {
567589 if ((gRandom ->Uniform ()) > cfgDownscalingFactor->get (kSpeciesRt ))
568590 return ;
569591 }
570592
571- if (cfgRapidityToggle && ((particle.y () - cfgRapidityCenterMass) < cfgRapidityMin || (particle.y () - cfgRapidityCenterMass) > cfgRapidityMax))
572- return ;
573-
574- if (cfgFillOnlyPhysicalPrimaries && !particle.isPhysicalPrimary ())
575- return ;
576-
577593 mHistograms .fill (HIST (nuclei::cNames[kSpeciesCt ]) + HIST (" /hTrackSelections" ), nuclei::trackSelection::kNoCuts );
578- if (!trackSelection<kSpeciesRt >(track, collision ))
594+ if (!trackSelection<kSpeciesRt , /* isMc */ true >(track))
579595 return ;
580596 mHistograms .fill (HIST (nuclei::cNames[kSpeciesCt ]) + HIST (" /hTrackSelections" ), nuclei::trackSelection::kTrackCuts );
581597
@@ -587,8 +603,8 @@ struct nucleiQC {
587603 candidate = fillCandidate</* isMc*/ true >(kSpeciesCt , collision, track);
588604 fillCollisionFlag (particle, candidate, reconstructedCollisions);
589605
590- mNucleiCandidates . emplace_back (candidate);
591- reconstructedMcParticles[ particle.globalIndex ()] = true ;
606+ writeCandidate (candidate);
607+ reconstructedMcParticles. insert ( particle.globalIndex ()) ;
592608
593609 dispatchFillHistograms</* isGenerated*/ true >(kSpeciesRt , candidate);
594610 dispatchFillHistograms</* isGenerated*/ false >(kSpeciesRt , candidate);
@@ -600,15 +616,14 @@ struct nucleiQC {
600616 for (const auto & particle : mcParticles) {
601617
602618 mcIndex++;
603-
604619 int iSpecies = nuclei::getSpeciesFromPdg (particle.pdgCode ());
605- if (std::find ( mSpeciesToProcess . begin (), mSpeciesToProcess . end (), iSpecies) == mSpeciesToProcess . end () )
620+ if (! mFillSpecies [ iSpecies] )
606621 continue ;
607622
608- if ((particle. y () - cfgRapidityCenterMass) < cfgRapidityMin || (particle. y () - cfgRapidityCenterMass) > cfgRapidityMax )
623+ if (reconstructedMcParticles. count (mcIndex) > 0 )
609624 continue ;
610625
611- if (reconstructedMcParticles[mcIndex] )
626+ if ((particle. y () - cfgRapidityCenterMass) < cfgRapidityMin || (particle. y () - cfgRapidityCenterMass) > cfgRapidityMax )
612627 continue ;
613628
614629 if (cfgFillOnlyPhysicalPrimaries && !particle.isPhysicalPrimary ())
@@ -626,34 +641,10 @@ struct nucleiQC {
626641 fillNucleusFlagsPdgsMc (particle, candidate);
627642 fillNucleusGeneratedVariables (particle, candidate);
628643
629- mNucleiCandidates . emplace_back (candidate);
644+ writeCandidate (candidate);
630645 mFilledMcParticleIds .emplace_back (particle.globalIndex ());
631646 dispatchFillHistograms</* isGenerated*/ true >(iSpecies, candidate);
632647 }
633-
634- if (!cfgFillTable)
635- return ;
636-
637- for (const auto & candidate : mNucleiCandidates ) {
638- mNucleiTableRed (
639- candidate.pt ,
640- candidate.eta ,
641- candidate.phi ,
642- candidate.tpcInnerParam ,
643- candidate.clusterSizesITS ,
644- candidate.TPCsignal ,
645- candidate.beta ,
646- candidate.DCAxy ,
647- candidate.DCAz ,
648- candidate.flags ,
649- candidate.ptGenerated ,
650- candidate.mcProcess ,
651- candidate.pdgCode ,
652- candidate.motherPdgCode );
653- mNucleiTableExt (
654- candidate.nsigmaTpc ,
655- candidate.nsigmaTof );
656- }
657648 }
658649 PROCESS_SWITCH (nucleiQC, processMc, " Mc analysis" , false );
659650};
0 commit comments