Skip to content

Commit 71bb524

Browse files
GiorgioAlbertoLuciaGiorgioAlbertoLucia
authored
[PWGLF] refactor nucleiqc logic, fixed previous errors (#15980)
1 parent eb920af commit 71bb524

2 files changed

Lines changed: 57 additions & 47 deletions

File tree

‎PWGLF/TableProducer/QC/nucleiQC.cxx‎

Lines changed: 40 additions & 45 deletions
Original file line numberDiff line numberDiff line change
@@ -78,10 +78,11 @@ enum trackQuality {
7878
kNCrossedRowsCut = 5,
7979
kTpcChi2Cut = 6,
8080
kItsChi2Cut = 7,
81-
kNtrackQuality = 8
81+
kDcaCuts = 8,
82+
kNtrackQuality = 9
8283
};
8384

84-
std::array<std::string, trackQuality::kNtrackQuality> trackQualityLabels{"All", "#eta cut", "#it{p}_{TPC}^{min} cut", "#it{N}_{cls}^{ITS} cut", "#it{N}_{cls}^{TPC} cut", "Crossed rows cut", "#chi^{2}_{TPC} cut", "#chi^{2}_{ITS} cut"};
85+
std::array<std::string, trackQuality::kNtrackQuality> trackQualityLabels{"All", "#eta cut", "#it{p}_{TPC}^{min} cut", "#it{N}_{cls}^{ITS} cut", "#it{N}_{cls}^{TPC} cut", "Crossed rows cut", "#chi^{2}_{TPC} cut", "#chi^{2}_{ITS} cut", "DCA cuts"};
8586

8687
} // namespace
8788

@@ -125,6 +126,8 @@ struct nucleiQC {
125126
Configurable<float> cfgCutNclusCrossedRowsTPC{"cfgCutNclusCrossedRowsTPC", 70, "Minimum number of TPC clusters crossed rows"};
126127
Configurable<float> cfgCutChi2PerClusterTPC{"cfgCutChi2PerClusterTPC", 4.f, "Maximum chi2 per TPC cluster"};
127128
Configurable<float> cfgCutChi2PerClusterITS{"cfgCutChi2PerClusterITS", 36.f, "Maximum chi2 per ITS cluster"};
129+
Configurable<float> cfgCutDCAxy{"cfgCutDCAxy", 10.f, "Maximum DCA in the transverse plane"};
130+
Configurable<float> cfgCutDCAz{"cfgCutDCAz", 10.f, "Maximum DCA in the longitudinal direction"};
128131

129132
Configurable<LabeledArray<double>> cfgNsigmaTPC{"cfgNsigmaTPC", {nuclei::nSigmaTPCdefault[0], nuclei::Species::kNspecies, 2, nuclei::names, nuclei::nSigmaConfigName}, "TPC nsigma selection for light nuclei"};
130133
Configurable<LabeledArray<double>> cfgNsigmaTOF{"cfgNsigmaTOF", {nuclei::nSigmaTOFdefault[0], nuclei::Species::kNspecies, 2, nuclei::names, nuclei::nSigmaConfigName}, "TPC nsigma selection for light nuclei"};
@@ -142,6 +145,7 @@ struct nucleiQC {
142145
{"hEventSelections", "Event selections; Selection step; Counts", {HistType::kTH1D, {{nuclei::evSel::kNevSels + 1, -0.5f, static_cast<float>(nuclei::evSel::kNevSels) + 0.5f}}}},
143146
{"hVtxZBefore", "Vertex distribution in Z before selections;Z (cm)", {HistType::kTH1F, {{400, -20.0, 20.0}}}},
144147
{"hVtxZ", "Vertex distribution in Z;Z (cm)", {HistType::kTH1F, {{400, -20.0, 20.0}}}},
148+
{"hCentrality", "Centrality distribution;Centrality (%)", {HistType::kTH1F, {{100, 0.0, 100.0}}}},
145149
{"hFailCentrality", "0: all the times the centrality filling function is called - 1: each time it fails ; Bool", {HistType::kTH1F, {{2, -0.5, 1.50}}}},
146150
{"hTrackTunedTracks", "", {HistType::kTH1F, {{1, 0.5, 1.5}}}},
147151
},
@@ -252,8 +256,8 @@ struct nucleiQC {
252256
LOGF(info, "Retrieved GRP for timestamp %ull (%i) with magnetic field of %1.2f kZG", timestamp, mRunNumber, mBz);
253257
}
254258

255-
template <int iSpecies, typename Ttrack>
256-
bool trackSelection(const Ttrack& track)
259+
template <int iSpecies, typename Ttrack, typename Tcollision>
260+
bool trackSelection(const Ttrack& track, const Tcollision& collision)
257261
{
258262
mHistograms.fill(HIST(nuclei::cNames[iSpecies]) + HIST("/hTrackQuality"), track.sign() * track.pt(), trackQuality::kNoCuts);
259263

@@ -285,6 +289,16 @@ struct nucleiQC {
285289
return false;
286290
mHistograms.fill(HIST(nuclei::cNames[iSpecies]) + HIST("/hTrackQuality"), track.sign() * track.pt(), trackQuality::kItsChi2Cut);
287291

292+
const o2::math_utils::Point3D<float> collisionVertex{collision.posX(), collision.posY(), collision.posZ()};
293+
mDcaInfoCov.set(999, 999, 999, 999, 999);
294+
setTrackParCov(track, mTrackParCov);
295+
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)
299+
return false;
300+
mHistograms.fill(HIST(nuclei::cNames[iSpecies]) + HIST("/hTrackQuality"), track.sign() * track.pt(), trackQuality::kDcaCuts);
301+
288302
return true;
289303
}
290304

@@ -293,7 +307,7 @@ struct nucleiQC {
293307
{
294308
constexpr int kIndex = iSpecies;
295309
if (!nuclei::checkSpeciesValidity(kIndex))
296-
std::runtime_error("species contains invalid nucleus kIndex");
310+
throw std::runtime_error("species contains invalid nucleus kIndex");
297311

298312
float centrality = nuclei::getCentrality(collision, cfgCentralityEstimator);
299313
float nsigmaTPC = mPidManagers[kIndex].getNSigmaTPC(track);
@@ -329,7 +343,7 @@ struct nucleiQC {
329343
}
330344

331345
if (particle.isPhysicalPrimary()) {
332-
candidate.flags |= nuclei::Flags::kIsPhysicalPrimary;
346+
candidate.flags |= nuclei::QcFlags::kQcIsPhysicalPrimary;
333347

334348
///< heavy flavour mother
335349
/*if (particle.has_mothers()) {
@@ -344,52 +358,33 @@ struct nucleiQC {
344358

345359
} else if (particle.getProcess() == TMCProcess::kPDecay) {
346360
///< assuming that strong decays are included in the previous step
347-
candidate.flags |= nuclei::Flags::kIsSecondaryFromWeakDecay;
361+
candidate.flags |= nuclei::QcFlags::kQcIsSecondaryFromWeakDecay;
348362
} else {
349-
candidate.flags |= nuclei::Flags::kIsSecondaryFromMaterial;
363+
candidate.flags |= nuclei::QcFlags::kQcIsSecondaryFromMaterial;
350364
}
351365
}
352366

353-
void fillSpeciesFlags(const int iSpecies, nuclei::SlimCandidate& candidate)
367+
template <typename Tparticle>
368+
void fillCollisionFlag(const Tparticle& particle, nuclei::SlimCandidate& candidate, const std::vector<bool>& reconstructedCollision)
354369
{
355-
356-
switch (iSpecies) {
357-
case nuclei::Species::kPr:
358-
candidate.flags |= nuclei::Flags::kProton;
359-
break;
360-
case nuclei::Species::kDe:
361-
candidate.flags |= nuclei::Flags::kDeuteron;
362-
break;
363-
case nuclei::Species::kTr:
364-
candidate.flags |= nuclei::Flags::kTriton;
365-
break;
366-
case nuclei::Species::kHe:
367-
candidate.flags |= nuclei::Flags::kHe3;
368-
break;
369-
case nuclei::Species::kAl:
370-
candidate.flags |= nuclei::Flags::kHe4;
371-
break;
372-
default:
373-
candidate.flags |= 0;
374-
break;
370+
if (reconstructedCollision[particle.mcCollisionId()]) {
371+
candidate.flags |= nuclei::QcFlags::kQcHasReconstructedCollision;
375372
}
376373
}
377374

378375
template <typename Tcollision, typename Ttrack>
379-
void fillNucleusFlagsPdgs(const int iSpecies, const Tcollision& collision, const Ttrack& track, nuclei::SlimCandidate& candidate)
376+
void fillNucleusFlagsPdgs(const Tcollision& collision, const Ttrack& track, nuclei::SlimCandidate& candidate)
380377
{
381378
candidate.flags = static_cast<uint16_t>((track.pidForTracking() & 0xF) << 12);
382379

383-
fillSpeciesFlags(iSpecies, candidate);
384-
385380
if (track.hasTOF())
386-
candidate.flags |= nuclei::Flags::kHasTOF;
381+
candidate.flags |= nuclei::QcFlags::kQcHasTOF;
387382

388383
if (track.hasTRD())
389-
candidate.flags |= nuclei::Flags::kHasTRD;
384+
candidate.flags |= nuclei::QcFlags::kQcHasTRD;
390385

391386
if (!collision.selection_bit(o2::aod::evsel::kNoITSROFrameBorder))
392-
candidate.flags |= nuclei::Flags::kITSrof;
387+
candidate.flags |= nuclei::QcFlags::kQcITSrof;
393388
}
394389

395390
template <typename Tparticle>
@@ -405,8 +400,6 @@ struct nucleiQC {
405400
void fillDcaInformation(const int iSpecies, const Tcollision& collision, const Ttrack& track, nuclei::SlimCandidate& candidate, const aod::McParticles::iterator& particle)
406401
{
407402

408-
const o2::math_utils::Point3D<float> collisionVertex{collision.posX(), collision.posY(), collision.posZ()};
409-
410403
mDcaInfoCov.set(999, 999, 999, 999, 999);
411404
setTrackParCov(track, mTrackParCov);
412405
mTrackParCov.setPID(track.pidForTracking());
@@ -432,7 +425,7 @@ struct nucleiQC {
432425
nuclei::SlimCandidate fillCandidate(const int iSpecies, Tcollision const& collision, Ttrack const& track)
433426
{
434427
if (!nuclei::checkSpeciesValidity(iSpecies))
435-
std::runtime_error("species contains invalid nucleus index");
428+
throw std::runtime_error("species contains invalid nucleus index");
436429

437430
nuclei::SlimCandidate candidate = {.pt = track.pt() * track.sign(),
438431
.eta = track.eta(),
@@ -455,7 +448,7 @@ struct nucleiQC {
455448
.nsigmaTpc = mPidManagers[iSpecies].getNSigmaTPC(track),
456449
.nsigmaTof = mPidManagers[iSpecies].getNSigmaTOF(track)};
457450

458-
fillNucleusFlagsPdgs(iSpecies, collision, track, candidate);
451+
fillNucleusFlagsPdgs(collision, track, candidate);
459452

460453
aod::McParticles::iterator particle;
461454

@@ -497,7 +490,7 @@ struct nucleiQC {
497490
{
498491
constexpr int kIndex = iSpecies;
499492
if (!nuclei::checkSpeciesValidity(kIndex))
500-
std::runtime_error("species contains invalid nucleus kIndex");
493+
throw std::runtime_error("species contains invalid nucleus kIndex");
501494

502495
if (isGenerated) {
503496
const float ptGenerated = (kIndex == nuclei::Species::kPr || kIndex == nuclei::Species::kDe || kIndex == nuclei::Species::kTr) ? candidate.ptGenerated : candidate.ptGenerated / 2.f;
@@ -517,11 +510,12 @@ struct nucleiQC {
517510
}
518511
}
519512

520-
void processMc(const Collisions& collisions, const TrackCandidatesMC& tracks, const aod::BCsWithTimestamps&, const aod::McParticles& mcParticles)
513+
void processMc(const Collisions& collisions, const TrackCandidatesMC& tracks, const aod::BCsWithTimestamps&, const aod::McParticles& mcParticles, const aod::McCollisions& mcCollisions)
521514
{
522515
gRandom->SetSeed(67);
523516
mNucleiCandidates.clear();
524517
std::vector<bool> reconstructedMcParticles(mcParticles.size(), false);
518+
std::vector<bool> reconstructedCollisions(mcCollisions.size(), false);
525519

526520
for (const auto& collision : collisions) {
527521

@@ -530,6 +524,8 @@ struct nucleiQC {
530524

531525
if (!nuclei::eventSelection(collision, mHistograms, cfgEventSelections, cfgCutVertex))
532526
continue;
527+
mHistograms.fill(HIST("hCentrality"), nuclei::getCentrality(collision, cfgCentralityEstimator, mHistFailCentrality));
528+
reconstructedCollisions[collision.mcCollisionId()] = true;
533529

534530
bool anyTrackTuner = false;
535531
for (int iSpecies = 0; iSpecies < static_cast<int>(nuclei::Species::kNspecies); iSpecies++) {
@@ -578,10 +574,8 @@ struct nucleiQC {
578574
if (cfgFillOnlyPhysicalPrimaries && !particle.isPhysicalPrimary())
579575
return;
580576

581-
LOG(info) << "track passed physical primary cut";
582-
583577
mHistograms.fill(HIST(nuclei::cNames[kSpeciesCt]) + HIST("/hTrackSelections"), nuclei::trackSelection::kNoCuts);
584-
if (!trackSelection<kSpeciesRt>(track))
578+
if (!trackSelection<kSpeciesRt>(track, collision))
585579
return;
586580
mHistograms.fill(HIST(nuclei::cNames[kSpeciesCt]) + HIST("/hTrackSelections"), nuclei::trackSelection::kTrackCuts);
587581

@@ -591,6 +585,7 @@ struct nucleiQC {
591585

592586
nuclei::SlimCandidate candidate;
593587
candidate = fillCandidate</*isMc*/ true>(kSpeciesCt, collision, track);
588+
fillCollisionFlag(particle, candidate, reconstructedCollisions);
594589

595590
mNucleiCandidates.emplace_back(candidate);
596591
reconstructedMcParticles[particle.globalIndex()] = true;
@@ -627,7 +622,7 @@ struct nucleiQC {
627622
nuclei::SlimCandidate candidate;
628623
// candidate.centrality = nuclei::getCentrality(collision, cfgCentralityEstimator, mHistFailCentrality);
629624
candidate.centrality = -1.f; // centrality is not well defined for non-reconstructed particles, set to -1 for now
630-
fillSpeciesFlags(iSpecies, candidate);
625+
fillCollisionFlag(particle, candidate, reconstructedCollisions);
631626
fillNucleusFlagsPdgsMc(particle, candidate);
632627
fillNucleusGeneratedVariables(particle, candidate);
633628

‎PWGLF/Utils/nucleiUtils.h‎

Lines changed: 17 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -148,6 +148,21 @@ enum Flags {
148148
kIsSecondaryFromWeakDecay = BIT(11) /// the last 4 bits are reserved for the PID in tracking
149149
};
150150

151+
enum QcFlags {
152+
kQcHasReconstructedCollision = BIT(0),
153+
kQcPlaceholder0 = BIT(1), /// placeholdedrs
154+
kQcPlaceholder1 = BIT(2), /// placeholdedrs
155+
kQcPlaceholder2 = BIT(3), /// placeholdedrs
156+
kQcPlaceholder3 = BIT(4), /// placeholdedrs
157+
kQcHasTOF = BIT(5),
158+
kQcHasTRD = BIT(6),
159+
kQcIsAmbiguous = BIT(7), /// just a placeholder now
160+
kQcITSrof = BIT(8),
161+
kQcIsPhysicalPrimary = BIT(9), /// MC flags starting from the second half of the short
162+
kQcIsSecondaryFromMaterial = BIT(10),
163+
kQcIsSecondaryFromWeakDecay = BIT(11) /// the last 4 bits are reserved for the PID in tracking
164+
};
165+
151166
constexpr int getSpeciesFromPdg(int pdg)
152167
{
153168
switch (std::abs(pdg)) {
@@ -419,7 +434,7 @@ void createHistogramRegistryNucleus(o2::framework::HistogramRegistry& registry)
419434

420435
constexpr int index = iSpecies;
421436
if (!checkSpeciesValidity(index)) {
422-
std::runtime_error("species contains invalid nucleus index");
437+
throw std::runtime_error("species contains invalid nucleus index");
423438
}
424439

425440
registry.add(fmt::format("{}/hTrackSelections", cNames[index]).c_str(), (fmt::format("{} track selections;", cNames[index]) + std::string("Selection step; Counts")).c_str(), o2::framework::HistType::kTH1D, {{trackSelection::kNtrackSelections, -0.5f, static_cast<float>(trackSelection::kNtrackSelections) - 0.5f}});
@@ -457,7 +472,7 @@ class PidManager
457472
: mSpecies(species)
458473
{
459474
if (!checkSpeciesValidity(species)) {
460-
std::runtime_error("species contains invalid nucleus index");
475+
throw std::runtime_error("species contains invalid nucleus index");
461476
}
462477

463478
if (!tpcBetheBlochParams) {

0 commit comments

Comments
 (0)