Skip to content

Commit d801b84

Browse files
GiorgioAlbertoLuciaGiorgioAlbertoLucia
committed
moved to unordered_set to fix an issue with reconstructed collisions and processed mc particles
1 parent 0bba10f commit d801b84

1 file changed

Lines changed: 92 additions & 102 deletions

File tree

‎PWGLF/TableProducer/QC/nucleiQC.cxx‎

Lines changed: 92 additions & 102 deletions
Original file line numberDiff line numberDiff line change
@@ -156,6 +156,7 @@ struct nucleiQC {
156156
std::shared_ptr<TH1> mHistTrackTunedTracks = mHistograms.get<TH1>(HIST("hTrackTunedTracks"));
157157

158158
std::vector<int> mSpeciesToProcess;
159+
std::array<bool, nuclei::Species::kNspecies> mFillSpecies{false};
159160
Produces<aod::NucleiTableRed> mNucleiTableRed;
160161
Produces<aod::NucleiTableExt> mNucleiTableExt;
161162

@@ -169,6 +170,7 @@ struct nucleiQC {
169170
o2::dataformats::VertexBase mVtx;
170171
o2::track::TrackParametrizationWithError<float> mTrackParCov;
171172
std::array<nuclei::PidManager, static_cast<int>(nuclei::Species::kNspecies)> mPidManagers;
173+
bool mAnyTrackTuner = false;
172174

173175
void init(o2::framework::InitContext&)
174176
{
@@ -177,22 +179,24 @@ struct nucleiQC {
177179
mCcdb->setCaching(true);
178180
mCcdb->setLocalObjectValidityChecking();
179181
mCcdb->setFatalWhenNull(false);
180-
nuclei::lut = o2::base::MatLayerCylSet::rectifyPtrFromFile(mCcdb->get<o2::base::MatLayerCylSet>("GLO/Param/MatLUT"));
182+
// nuclei::lut = o2::base::MatLayerCylSet::rectifyPtrFromFile(mCcdb->get<o2::base::MatLayerCylSet>("GLO/Param/MatLUT"));
181183

182184
for (int iSel = 0; iSel < nuclei::evSel::kNevSels; iSel++) {
183185
mHistograms.get<TH1>(HIST("hEventSelections"))->GetXaxis()->SetBinLabel(iSel + 1, nuclei::eventSelectionLabels[iSel].c_str());
184186
}
185187

186188
for (int iSpecies = 0; iSpecies < static_cast<int>(nuclei::Species::kNspecies); iSpecies++) {
187-
if (cfgSpeciesToProcess->get(iSpecies) == 1)
189+
if (cfgSpeciesToProcess->get(iSpecies) == 1) {
188190
mSpeciesToProcess.emplace_back(iSpecies);
191+
mFillSpecies[iSpecies] = true;
192+
}
189193
}
190194

191195
static_for<0, nuclei::kNspecies - 1>([&](auto iSpecies) {
192196
constexpr int kSpeciesCt = decltype(iSpecies)::value;
193197
const int kSpeciesRt = kSpeciesCt;
194198

195-
if (std::find(mSpeciesToProcess.begin(), mSpeciesToProcess.end(), kSpeciesCt) == mSpeciesToProcess.end())
199+
if (!mFillSpecies[kSpeciesCt])
196200
return;
197201

198202
float tpcBetheBlochParams[6];
@@ -214,12 +218,12 @@ struct nucleiQC {
214218
});
215219

216220
/// TrackTuner initialization
217-
bool anyTrackTuner = false;
221+
mAnyTrackTuner = false;
218222
for (int iSpecies = 0; iSpecies < static_cast<int>(nuclei::Species::kNspecies); iSpecies++) {
219-
anyTrackTuner = anyTrackTuner || cfgUseTrackTuner->get(iSpecies);
223+
mAnyTrackTuner = mAnyTrackTuner || cfgUseTrackTuner->get(iSpecies);
220224
}
221225

222-
if (anyTrackTuner) {
226+
if (mAnyTrackTuner) {
223227
std::string outputStringParams = "";
224228
switch (cfgTrackTunerConfigSource) {
225229
case aod::track_tuner::InputString:
@@ -251,12 +255,16 @@ struct nucleiQC {
251255

252256
o2::parameters::GRPMagField* grpmag = mCcdb->getForTimeStamp<o2::parameters::GRPMagField>("GLO/Config/GRPMagField", timestamp);
253257
o2::base::Propagator::initFieldFromGRP(grpmag);
254-
o2::base::Propagator::Instance()->setMatLUT(nuclei::lut);
258+
// o2::base::Propagator::Instance()->setMatLUT(nuclei::lut);
259+
if (!o2::base::Propagator::Instance()->getMatLUT()) {
260+
nuclei::lut = o2::base::MatLayerCylSet::rectifyPtrFromFile(mCcdb->get<o2::base::MatLayerCylSet>("GLO/Param/MatLUT"));
261+
o2::base::Propagator::Instance()->setMatLUT(nuclei::lut);
262+
}
255263
mBz = static_cast<float>(grpmag->getNominalL3Field());
256264
LOGF(info, "Retrieved GRP for timestamp %ull (%i) with magnetic field of %1.2f kZG", timestamp, mRunNumber, mBz);
257265
}
258266

259-
template <int iSpecies, typename Ttrack, typename Tcollision>
267+
template <int iSpecies, bool isMc, typename Ttrack, typename Tcollision>
260268
bool trackSelection(const Ttrack& track, const Tcollision& collision)
261269
{
262270
mHistograms.fill(HIST(nuclei::cNames[iSpecies]) + HIST("/hTrackQuality"), track.sign() * track.pt(), trackQuality::kNoCuts);
@@ -289,13 +297,24 @@ struct nucleiQC {
289297
return false;
290298
mHistograms.fill(HIST(nuclei::cNames[iSpecies]) + HIST("/hTrackQuality"), track.sign() * track.pt(), trackQuality::kItsChi2Cut);
291299

292-
const o2::math_utils::Point3D<float> collisionVertex{collision.posX(), collision.posY(), collision.posZ()};
293300
mDcaInfoCov.set(999, 999, 999, 999, 999);
294301
setTrackParCov(track, mTrackParCov);
295302
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)
303+
304+
if constexpr (isMc) {
305+
if (track.has_mcParticle() && cfgUseTrackTuner->get(iSpecies)) {
306+
const auto& particle = track.mcParticle();
307+
mHistTrackTunedTracks->Fill(1.);
308+
mTrackTuner.tuneTrackParams(particle, mTrackParCov, mMatCorr, &mDcaInfoCov, mHistTrackTunedTracks);
309+
}
310+
} else {
311+
mMatCorr = static_cast<o2::base::Propagator::MatCorrType>(cfgMaterialCorrection.value);
312+
}
313+
314+
o2::base::Propagator::Instance()->propagateToDCABxByBz(mVtx, mTrackParCov, 2.f, mMatCorr, &mDcaInfoCov);
315+
mDcaInfo[0] = mDcaInfoCov.getY();
316+
mDcaInfo[1] = mDcaInfoCov.getZ();
317+
if (std::abs(mDcaInfo[0]) > cfgCutDCAxy || std::abs(mDcaInfo[1]) > cfgCutDCAz)
299318
return false;
300319
mHistograms.fill(HIST(nuclei::cNames[iSpecies]) + HIST("/hTrackQuality"), track.sign() * track.pt(), trackQuality::kDcaCuts);
301320

@@ -365,9 +384,9 @@ struct nucleiQC {
365384
}
366385

367386
template <typename Tparticle>
368-
void fillCollisionFlag(const Tparticle& particle, nuclei::SlimCandidate& candidate, const std::vector<bool>& reconstructedCollision)
387+
void fillCollisionFlag(const Tparticle& particle, nuclei::SlimCandidate& candidate, const std::unordered_set<int>& reconstructedCollision)
369388
{
370-
if (reconstructedCollision[particle.mcCollisionId()]) {
389+
if (reconstructedCollision.count(particle.mcCollisionId()) > 0) {
371390
candidate.flags |= nuclei::QcFlags::kQcHasReconstructedCollision;
372391
}
373392
}
@@ -396,31 +415,6 @@ struct nucleiQC {
396415
candidate.phiGenerated = particle.phi();
397416
}
398417

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-
424418
template <const bool isMc, typename Tcollision, typename Ttrack>
425419
nuclei::SlimCandidate fillCandidate(const int iSpecies, Tcollision const& collision, Ttrack const& track)
426420
{
@@ -434,8 +428,8 @@ struct nucleiQC {
434428
.clusterSizesITS = track.itsClusterSizes(),
435429
.TPCsignal = track.tpcSignal(),
436430
.beta = mPidManagers[iSpecies].getBetaTOF(track),
437-
.DCAxy = 0.f,
438-
.DCAz = 0.f,
431+
.DCAxy = mDcaInfo[0], // set in the track selection function
432+
.DCAz = mDcaInfo[1], // set in the track selection function
439433
.flags = 0,
440434
.pdgCode = 0,
441435
.motherPdgCode = 0,
@@ -450,19 +444,15 @@ struct nucleiQC {
450444

451445
fillNucleusFlagsPdgs(collision, track, candidate);
452446

453-
aod::McParticles::iterator particle;
454-
455447
if constexpr (isMc) {
456448
if (track.has_mcParticle()) {
457449

458-
particle = track.mcParticle();
450+
const auto& particle = track.mcParticle();
459451
fillNucleusFlagsPdgsMc(particle, candidate);
460452
fillNucleusGeneratedVariables(particle, candidate);
461453
}
462454
}
463455

464-
fillDcaInformation<isMc>(iSpecies, collision, track, candidate, particle);
465-
466456
return candidate;
467457
}
468458

@@ -510,12 +500,36 @@ struct nucleiQC {
510500
}
511501
}
512502

503+
void writeCandidate(const nuclei::SlimCandidate& candidate)
504+
{
505+
if (!cfgFillTable)
506+
return;
507+
508+
mNucleiTableRed(
509+
candidate.pt,
510+
candidate.eta,
511+
candidate.phi,
512+
candidate.tpcInnerParam,
513+
candidate.clusterSizesITS,
514+
candidate.TPCsignal,
515+
candidate.beta,
516+
candidate.DCAxy,
517+
candidate.DCAz,
518+
candidate.flags,
519+
candidate.ptGenerated,
520+
candidate.mcProcess,
521+
candidate.pdgCode,
522+
candidate.motherPdgCode);
523+
mNucleiTableExt(
524+
candidate.nsigmaTpc,
525+
candidate.nsigmaTof);
526+
}
527+
513528
void processMc(const Collisions& collisions, const TrackCandidatesMC& tracks, const aod::BCsWithTimestamps&, const aod::McParticles& mcParticles, const aod::McCollisions& mcCollisions)
514529
{
515530
gRandom->SetSeed(67);
516-
mNucleiCandidates.clear();
517-
std::vector<bool> reconstructedMcParticles(mcParticles.size(), false);
518-
std::vector<bool> reconstructedCollisions(mcCollisions.size(), false);
531+
std::unordered_set<int> reconstructedMcParticles;
532+
std::unordered_set<int> reconstructedCollisions;
519533

520534
for (const auto& collision : collisions) {
521535

@@ -525,13 +539,11 @@ struct nucleiQC {
525539
if (!nuclei::eventSelection(collision, mHistograms, cfgEventSelections, cfgCutVertex))
526540
continue;
527541
mHistograms.fill(HIST("hCentrality"), nuclei::getCentrality(collision, cfgCentralityEstimator, mHistFailCentrality));
528-
reconstructedCollisions[collision.mcCollisionId()] = true;
542+
reconstructedCollisions.insert(collision.mcCollisionId());
543+
mVtx.setPos({collision.posX(), collision.posY(), collision.posZ()});
544+
mVtx.setCov(collision.covXX(), collision.covXY(), collision.covYY(), collision.covXZ(), collision.covYZ(), collision.covZZ());
529545

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) {
546+
if (mAnyTrackTuner && mTrackTuner.autoDetectDcaCalib && !mTrackTuner.areGraphsConfigured) {
535547

536548
mTrackTuner.setRunNumber(mRunNumber);
537549

@@ -544,38 +556,41 @@ struct nucleiQC {
544556
auto tracksThisCollision = tracks.sliceBy(mTracksPerCollision, collision.globalIndex());
545557
for (const auto& track : tracksThisCollision) {
546558

547-
static_for<0, nuclei::kNspecies - 1>([&](auto iSpecies) {
548-
constexpr int kSpeciesCt = decltype(iSpecies)::value;
549-
const int kSpeciesRt = kSpeciesCt;
559+
// selections shared among all species
560+
if (!track.has_mcParticle())
561+
continue;
550562

551-
if (std::find(mSpeciesToProcess.begin(), mSpeciesToProcess.end(), kSpeciesRt) == mSpeciesToProcess.end())
552-
return;
563+
if (track.mcParticleId() < -1 || track.mcParticleId() >= mcParticles.size())
564+
continue;
565+
const auto& particle = mcParticles.iteratorAt(track.mcParticleId());
553566

554-
if (!track.has_mcParticle())
555-
return;
567+
if (cfgRapidityToggle && ((particle.y() - cfgRapidityCenterMass) < cfgRapidityMin || (particle.y() - cfgRapidityCenterMass) > cfgRapidityMax))
568+
continue;
569+
570+
if (cfgFillOnlyPhysicalPrimaries && !particle.isPhysicalPrimary())
571+
continue;
556572

557-
if (track.mcParticleId() < -1 || track.mcParticleId() >= mcParticles.size())
573+
// species-specific selections and filling
574+
static_for<0, nuclei::kNspecies - 1>([&](auto iSpecies) {
575+
constexpr int kSpeciesCt = decltype(iSpecies)::value;
576+
const int kSpeciesRt = kSpeciesCt;
577+
if (!mFillSpecies[kSpeciesCt])
558578
return;
559-
const auto& particle = mcParticles.iteratorAt(track.mcParticleId());
560579

561580
if (cfgDoCheckPdgCode) {
562581
if (std::abs(particle.pdgCode()) != nuclei::pdgCodes[kSpeciesRt])
563582
return;
564583
}
565584

585+
mDcaInfo = {-999.f, -999.f};
586+
566587
if (cfgDownscalingFactor->get(kSpeciesRt) < 1.) {
567588
if ((gRandom->Uniform()) > cfgDownscalingFactor->get(kSpeciesRt))
568589
return;
569590
}
570591

571-
if (cfgRapidityToggle && ((particle.y() - cfgRapidityCenterMass) < cfgRapidityMin || (particle.y() - cfgRapidityCenterMass) > cfgRapidityMax))
572-
return;
573-
574-
if (cfgFillOnlyPhysicalPrimaries && !particle.isPhysicalPrimary())
575-
return;
576-
577592
mHistograms.fill(HIST(nuclei::cNames[kSpeciesCt]) + HIST("/hTrackSelections"), nuclei::trackSelection::kNoCuts);
578-
if (!trackSelection<kSpeciesRt>(track, collision))
593+
if (!trackSelection<kSpeciesRt, /*isMc*/ true>(track, collision))
579594
return;
580595
mHistograms.fill(HIST(nuclei::cNames[kSpeciesCt]) + HIST("/hTrackSelections"), nuclei::trackSelection::kTrackCuts);
581596

@@ -587,8 +602,8 @@ struct nucleiQC {
587602
candidate = fillCandidate</*isMc*/ true>(kSpeciesCt, collision, track);
588603
fillCollisionFlag(particle, candidate, reconstructedCollisions);
589604

590-
mNucleiCandidates.emplace_back(candidate);
591-
reconstructedMcParticles[particle.globalIndex()] = true;
605+
writeCandidate(candidate);
606+
reconstructedMcParticles.insert(particle.globalIndex());
592607

593608
dispatchFillHistograms</*isGenerated*/ true>(kSpeciesRt, candidate);
594609
dispatchFillHistograms</*isGenerated*/ false>(kSpeciesRt, candidate);
@@ -600,15 +615,14 @@ struct nucleiQC {
600615
for (const auto& particle : mcParticles) {
601616

602617
mcIndex++;
603-
604618
int iSpecies = nuclei::getSpeciesFromPdg(particle.pdgCode());
605-
if (std::find(mSpeciesToProcess.begin(), mSpeciesToProcess.end(), iSpecies) == mSpeciesToProcess.end())
619+
if (!mFillSpecies[iSpecies])
606620
continue;
607621

608-
if ((particle.y() - cfgRapidityCenterMass) < cfgRapidityMin || (particle.y() - cfgRapidityCenterMass) > cfgRapidityMax)
622+
if (reconstructedMcParticles.count(particle.globalIndex()) > 0)
609623
continue;
610624

611-
if (reconstructedMcParticles[mcIndex])
625+
if ((particle.y() - cfgRapidityCenterMass) < cfgRapidityMin || (particle.y() - cfgRapidityCenterMass) > cfgRapidityMax)
612626
continue;
613627

614628
if (cfgFillOnlyPhysicalPrimaries && !particle.isPhysicalPrimary())
@@ -626,34 +640,10 @@ struct nucleiQC {
626640
fillNucleusFlagsPdgsMc(particle, candidate);
627641
fillNucleusGeneratedVariables(particle, candidate);
628642

629-
mNucleiCandidates.emplace_back(candidate);
643+
writeCandidate(candidate);
630644
mFilledMcParticleIds.emplace_back(particle.globalIndex());
631645
dispatchFillHistograms</*isGenerated*/ true>(iSpecies, candidate);
632646
}
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-
}
657647
}
658648
PROCESS_SWITCH(nucleiQC, processMc, "Mc analysis", false);
659649
};

0 commit comments

Comments
 (0)