Skip to content

Commit 4d4093e

Browse files
blacwovieblacwovie
authored
Update HadNucleiFemto.cxx
1 parent 705ec57 commit 4d4093e

1 file changed

Lines changed: 205 additions & 22 deletions

File tree

‎PWGCF/Femto/FemtoNuclei/TableProducer/HadNucleiFemto.cxx‎

Lines changed: 205 additions & 22 deletions
Original file line numberDiff line numberDiff line change
@@ -78,6 +78,7 @@ using CollBracket = o2::math_utils::Bracket<int>;
7878
using CollisionsFull = soa::Join<aod::Collisions, aod::EvSels, aod::CentFT0Cs, aod::FT0Mults>;
7979
using CollisionsFullMC = soa::Join<aod::Collisions, aod::McCollisionLabels, aod::EvSels, aod::CentFT0Cs, aod::FT0Mults>;
8080
using TrackCandidates = soa::Join<aod::TracksIU, aod::TracksExtra, aod::TracksCovIU, aod::TracksDCA, aod::TrackSelection, aod::pidTPCFullDe, aod::pidTPCFullTr, aod::pidTOFFullDe, aod::pidTOFFullTr, aod::pidTOFFullHe, aod::pidTPCFullPr, aod::pidTOFFullPr, aod::pidTPCFullPi, aod::pidTOFFullPi, aod::pidTPCFullKa, aod::pidTOFFullKa, aod::TOFSignal, aod::TOFEvTime>;
81+
using TrackCandidatesMC = soa::Join<aod::TracksIU, aod::TracksExtra, aod::TracksCovIU, aod::TracksDCA, aod::TrackSelection, aod::pidTPCFullDe, aod::pidTPCFullTr, aod::pidTOFFullDe, aod::pidTOFFullTr, aod::pidTOFFullHe, aod::pidTPCFullPr, aod::pidTOFFullPr, aod::pidTPCFullPi, aod::pidTOFFullPi, aod::pidTPCFullKa, aod::pidTOFFullKa, aod::TOFSignal, aod::TOFEvTime, aod::McTrackLabels>;
8182

8283
namespace
8384
{
@@ -179,6 +180,7 @@ struct HadNucandidate {
179180
struct HadNucleiFemto {
180181

181182
Produces<aod::HadronNucleiTable> mOutputDataTable;
183+
Produces<aod::HadronNucleiTableMC> mOutputMCTable;
182184
Produces<aod::HadronHyperTable> mOutputHyperDataTable;
183185
Produces<aod::HadronNucleiMult> mOutputMultiplicityTable;
184186

@@ -242,7 +244,6 @@ struct HadNucleiFemto {
242244
std::string prefix{"deuteronPid"};
243245
// Deuteron purity and PID cuts
244246
Configurable<float> settingCutPinMinDe{"settingCutPinMinDe", 0.0f, "Minimum Pin for De"};
245-
Configurable<float> settingCutClSizeItsDe{"settingCutClSizeItsDe", 4.0f, "Minimum ITS cluster size for De"};
246247
Configurable<float> settingCutDeptMin{"settingCutDeptMin", 0.6f, "Minimum PT cut on De"};
247248
Configurable<float> settingCutDeptMax{"settingCutDeptMax", 1.6f, "Maximum PT cut on De"};
248249
Configurable<float> settingCutPinMinTOFITSDe{"settingCutPinMinTOFITSDe", 1.2f, "Minimum p to apply the TOF ITS cut on De"};
@@ -289,6 +290,16 @@ struct HadNucleiFemto {
289290
Configurable<float> settingClosePairSpecificRadius{"settingClosePairSpecificRadius", 85.f, "TPC radius in cm used when close pair rejection mode is 2"};
290291
} CPR;
291292

293+
struct : o2::framework::ConfigurableGroup {
294+
// cppcheck-suppress unusedStructMember
295+
std::string prefix{"mc"};
296+
Configurable<bool> settingRequireSel8{"settingRequireSel8", true, "Apply the same sel8 event selection to reconstructed MC as to data"};
297+
Configurable<bool> settingRequireTruthSpecies{"settingRequireTruthSpecies", true, "Store only truth-matched pion-nucleus pairs"};
298+
Configurable<bool> settingRequireSameMCCollision{"settingRequireSameMCCollision", true, "Require the two truth particles to come from the same MC collision"};
299+
Configurable<bool> settingRequireRecoMCCollisionMatch{"settingRequireRecoMCCollisionMatch", true, "Require both truth particles to match the reconstructed collision MC label"};
300+
Configurable<bool> settingRequirePhysicalPrimaries{"settingRequirePhysicalPrimaries", false, "Store only pairs in which both truth particles are physical primaries"};
301+
} mc;
302+
292303
struct : o2::framework::ConfigurableGroup {
293304
// cppcheck-suppress unusedStructMember
294305
std::string prefix{"hypertriton"};
@@ -339,6 +350,7 @@ struct HadNucleiFemto {
339350
} pidCalibration;
340351

341352
Preslice<TrackCandidates> mPerCol = aod::track::collisionId;
353+
Preslice<TrackCandidatesMC> mPerColMC = aod::track::collisionId;
342354
PresliceUnsorted<o2::aod::DataHypCandsWColl> hypPerCol = o2::aod::hyperrec::collisionId;
343355

344356
// binning for EM background
@@ -398,6 +410,12 @@ struct HadNucleiFemto {
398410
{"h2CPRBefore", "Close pair rejection before cut; #Delta#eta; #Delta#phi^{*}", {HistType::kTH2F, {{300, -0.15f, 0.15f}, {400, -0.2f, 0.2f}}}},
399411
{"h2CPRAfter", "Close pair rejection after cut; #Delta#eta; #Delta#phi^{*}", {HistType::kTH2F, {{300, -0.15f, 0.15f}, {400, -0.2f, 0.2f}}}},
400412

413+
// Reconstructed MC pair QA
414+
{"MC/hPairFlow", "MC pair flow;step;counts", {HistType::kTH1F, {{7, -0.5f, 6.5f}}}},
415+
{"MC/hKstarRecVsGen", "Reconstructed versus generated k*;generated k* (GeV/c);reconstructed k* (GeV/c)", {HistType::kTH2F, {{300, 0.f, 3.f}, {300, 0.f, 3.f}}}},
416+
{"MC/hPtNuRecVsGen", "Reconstructed versus generated signed nucleus pT;generated pT (GeV/c);reconstructed pT (GeV/c)", {HistType::kTH2F, {{280, -7.f, 7.f}, {280, -7.f, 7.f}}}},
417+
{"MC/hPtHadRecVsGen", "Reconstructed versus generated signed pion pT;generated pT (GeV/c);reconstructed pT (GeV/c)", {HistType::kTH2F, {{280, -7.f, 7.f}, {280, -7.f, 7.f}}}},
418+
401419
// dE/dx
402420
{"h2dEdxNucandidates", "dEdx distribution; #it{p} (GeV/#it{c}); dE/dx (a.u.)", {HistType::kTH2F, {{200, -5.0f, 5.0f}, {100, 0.0f, 2000.0f}}}},
403421
{"h2dEdxHadcandidates", "dEdx distribution; #it{p} (GeV/#it{c}); dE/dx (a.u.)", {HistType::kTH2F, {{200, -5.0f, 5.0f}, {100, 0.0f, 2000.0f}}}},
@@ -466,6 +484,20 @@ struct HadNucleiFemto {
466484

467485
void init(o2::framework::InitContext&)
468486
{
487+
constexpr int closePairRadiusModePv = 0;
488+
constexpr int closePairRadiusModeSpecificTpc = 2;
489+
if (CPR.settingEnableClosePairRejection.value) {
490+
if (CPR.settingClosePairDeltaEtaMax.value <= 0.f || CPR.settingClosePairDeltaPhiMax.value <= 0.f) {
491+
LOG(fatal) << "Close-pair rejection requires positive delta-eta and delta-phi-star limits";
492+
}
493+
if (CPR.settingClosePairRadiusMode.value < closePairRadiusModePv || CPR.settingClosePairRadiusMode.value > closePairRadiusModeSpecificTpc) {
494+
LOG(fatal) << "Invalid close-pair radius mode " << CPR.settingClosePairRadiusMode.value << "; expected 0, 1, or 2";
495+
}
496+
if (CPR.settingClosePairRadiusMode.value == closePairRadiusModeSpecificTpc && CPR.settingClosePairSpecificRadius.value <= 0.f) {
497+
LOG(fatal) << "Close-pair rejection at a specific TPC radius requires a positive radius";
498+
}
499+
}
500+
469501
mZorroSummary.setObject(mZorro.getZorroSummary());
470502
mRunNumber = 0;
471503

@@ -505,17 +537,31 @@ struct HadNucleiFemto {
505537
}
506538
}
507539

540+
template <bool isMC>
508541
void initCCDB(const aod::BCsWithTimestamps::iterator& bc)
509542
{
510543
if (mRunNumber == bc.runNumber()) {
511544
return;
512545
}
513-
if (zorro.settingSkimmedProcessing) {
514-
mZorro.initCCDB(mCcdb.service, bc.runNumber(), bc.timestamp(), zorroTriggerMask());
515-
mZorro.populateHistRegistry(mQaRegistry, bc.runNumber());
546+
if constexpr (!isMC) {
547+
if (zorro.settingSkimmedProcessing) {
548+
mZorro.initCCDB(mCcdb.service, bc.runNumber(), bc.timestamp(), zorroTriggerMask());
549+
mZorro.populateHistRegistry(mQaRegistry, bc.runNumber());
550+
}
516551
}
517552
mRunNumber = bc.runNumber();
518553
const float defaultBzValue = -999.0f;
554+
555+
// A fixed field is sufficient for CPR and DCAFitter when material
556+
// corrections are disabled. This also makes local MC tests independent
557+
// of an AliEn token when the CCDB payload is stored on Grid.
558+
if (ccdb.settingDbz > defaultBzValue && pidCalibration.settingMaterialCorrection == 0) {
559+
mDbz = ccdb.settingDbz;
560+
mFitter.setBz(mDbz);
561+
LOG(info) << "Using configured magnetic field of " << mDbz << " kZG";
562+
return;
563+
}
564+
519565
auto run3GrpTimestamp = bc.timestamp();
520566
auto* grpo = mCcdb->getForTimeStamp<o2::parameters::GRPObject>(ccdb.settingGrpPath, run3GrpTimestamp);
521567
o2::parameters::GRPMagField* grpmag = nullptr;
@@ -542,31 +588,47 @@ struct HadNucleiFemto {
542588
mDbz = ccdb.settingDbz;
543589
}
544590
}
591+
mFitter.setBz(mDbz);
545592
}
546593

547594
// ==================================================================================================================
548595

549596
template <bool isMC, typename Tcollision>
550-
bool selectCollision(const Tcollision& collision, const aod::BCsWithTimestamps&)
597+
bool passesCollisionSelection(const Tcollision& collision)
551598
{
552-
mQaRegistry.fill(HIST("hEvents"), 0);
599+
// CPR uses phi* and therefore needs the magnetic field for MC as well.
600+
auto bc = collision.template bc_as<aod::BCsWithTimestamps>();
601+
initCCDB<isMC>(bc);
553602

554603
if constexpr (isMC) {
555-
if (/*!collision.sel8() ||*/ std::abs(collision.posZ()) > eventMixing.settingCutVertex) {
604+
if ((mc.settingRequireSel8.value && !collision.sel8()) || std::abs(collision.posZ()) > eventMixing.settingCutVertex) {
556605
return false;
557606
}
558607
} else {
559-
auto bc = collision.template bc_as<aod::BCsWithTimestamps>();
560-
initCCDB(bc);
561-
562608
if (!collision.sel8() || std::abs(collision.posZ()) > eventMixing.settingCutVertex) {
563609
return false;
564610
}
565611
if (zorro.settingSkimmedProcessing) {
566-
const bool zorroSelected = mZorro.isSelected(collision.template bc_as<aod::BCsWithTimestamps>().globalBC());
567-
if (!zorroSelected) {
612+
if (!mZorro.isSelected(bc.globalBC())) {
568613
return false;
569614
}
615+
}
616+
}
617+
618+
return true;
619+
}
620+
621+
template <bool isMC, typename Tcollision>
622+
bool selectCollision(const Tcollision& collision, const aod::BCsWithTimestamps&)
623+
{
624+
mQaRegistry.fill(HIST("hEvents"), 0);
625+
626+
if (!passesCollisionSelection<isMC>(collision)) {
627+
return false;
628+
}
629+
630+
if constexpr (!isMC) {
631+
if (zorro.settingSkimmedProcessing) {
570632
mQaRegistry.fill(HIST("hEvents"), 2);
571633
}
572634
}
@@ -894,7 +956,7 @@ struct HadNucleiFemto {
894956
}
895957

896958
template <typename Ttrack1, typename Ttrack2>
897-
bool isClosePair(const Ttrack1& firstTrack, const Ttrack2& secondTrack)
959+
bool isClosePair(const Ttrack1& firstTrack, const Ttrack2& secondTrack, bool fillQA)
898960
{
899961
if (!CPR.settingEnableClosePairRejection.value) {
900962
return false;
@@ -909,11 +971,13 @@ struct HadNucleiFemto {
909971
return false;
910972
}
911973

912-
mQaRegistry.fill(HIST("h2CPRBefore"), deltaEta, deltaPhi);
974+
if (fillQA) {
975+
mQaRegistry.fill(HIST("h2CPRBefore"), deltaEta, deltaPhi);
976+
}
913977
const bool isRejected = std::pow(deltaPhi, 2.f) / std::pow(CPR.settingClosePairDeltaPhiMax.value, 2.f) +
914978
std::pow(deltaEta, 2.f) / std::pow(CPR.settingClosePairDeltaEtaMax.value, 2.f) <
915979
1.f;
916-
if (!isRejected) {
980+
if (fillQA && !isRejected) {
917981
mQaRegistry.fill(HIST("h2CPRAfter"), deltaEta, deltaPhi);
918982
}
919983
return isRejected;
@@ -1673,7 +1737,7 @@ struct HadNucleiFemto {
16731737
continue;
16741738
}
16751739
hasHadronSelected = true;
1676-
if (isClosePair(track0, track1)) {
1740+
if (isClosePair(track0, track1, /*fillQA*/ true)) {
16771741
continue;
16781742
}
16791743

@@ -1739,7 +1803,7 @@ struct HadNucleiFemto {
17391803
if (!selectTrackHadron(hadCand) || !selectionPIDHadron(hadCand)) {
17401804
continue;
17411805
}
1742-
if (isClosePair(DeCand, hadCand)) {
1806+
if (isClosePair(DeCand, hadCand, /*fillQA*/ true)) {
17431807
continue;
17441808
}
17451809

@@ -1811,6 +1875,35 @@ struct HadNucleiFemto {
18111875
}
18121876
}
18131877

1878+
template <typename TparticleNu, typename TparticleHad>
1879+
void fillMCTable(const HadNucandidate& hadNucand, const TparticleNu& particleNu, const TparticleHad& particleHad, bool sameMCCollision, bool matchesRecoMCCollision)
1880+
{
1881+
const float signedPtNuMC = particleNu.pdgCode() >= 0 ? particleNu.pt() : -particleNu.pt();
1882+
const float signedPtHadMC = particleHad.pdgCode() >= 0 ? particleHad.pt() : -particleHad.pt();
1883+
const std::array<float, 3> momentumNuMC{particleNu.px(), particleNu.py(), particleNu.pz()};
1884+
const std::array<float, 3> momentumHadMC{particleHad.px(), particleHad.py(), particleHad.pz()};
1885+
const float kstarMC = computePairKstar(momentumHadMC, mMassHad, momentumNuMC, nucleusMass());
1886+
1887+
mOutputMCTable(
1888+
signedPtNuMC,
1889+
particleNu.eta(),
1890+
particleNu.phi(),
1891+
signedPtHadMC,
1892+
particleHad.eta(),
1893+
particleHad.phi(),
1894+
kstarMC,
1895+
particleNu.pdgCode(),
1896+
particleHad.pdgCode(),
1897+
particleNu.isPhysicalPrimary(),
1898+
particleHad.isPhysicalPrimary(),
1899+
sameMCCollision,
1900+
matchesRecoMCCollision);
1901+
1902+
mQaRegistry.fill(HIST("MC/hKstarRecVsGen"), kstarMC, hadNucand.kstar);
1903+
mQaRegistry.fill(HIST("MC/hPtNuRecVsGen"), signedPtNuMC, hadNucand.recoPtNu());
1904+
mQaRegistry.fill(HIST("MC/hPtHadRecVsGen"), signedPtHadMC, hadNucand.recoPtHad());
1905+
}
1906+
18141907
template <typename Tcoll>
18151908
void fillTableHyper(const HadNucandidate& hadNucand, const Tcoll& collision)
18161909
{
@@ -1954,6 +2047,89 @@ struct HadNucleiFemto {
19542047

19552048
// ==================================================================================================================
19562049

2050+
void processMC(const CollisionsFullMC& collisions, const TrackCandidatesMC& tracks, const aod::McParticles&, const aod::BCsWithTimestamps& bcs)
2051+
{
2052+
mGoodCollisions.clear();
2053+
mGoodCollisions.resize(collisions.size(), false);
2054+
2055+
for (const auto& collision : collisions) {
2056+
mTrackPairs.clear();
2057+
2058+
if (!selectCollision</*isMC*/ true>(collision, bcs)) {
2059+
continue;
2060+
}
2061+
if (mc.settingRequireRecoMCCollisionMatch.value && !collision.has_mcCollision()) {
2062+
continue;
2063+
}
2064+
2065+
const uint64_t collIdx = collision.globalIndex();
2066+
mGoodCollisions[collIdx] = true;
2067+
auto tracksThisCollision = tracks.sliceBy(mPerColMC, collIdx);
2068+
tracksThisCollision.bindExternalIndices(&tracks);
2069+
2070+
// This is deliberately the same pair builder as for data. It applies
2071+
// nucleus/pion selections and CPR before truth matching.
2072+
pairTracksSameEvent(tracksThisCollision, collision.centFT0C());
2073+
2074+
for (const auto& trackPair : mTrackPairs) {
2075+
mQaRegistry.fill(HIST("MC/hPairFlow"), 0);
2076+
auto trackNu = tracks.rawIteratorAt(trackPair.tr0Idx);
2077+
auto trackHad = tracks.rawIteratorAt(trackPair.tr1Idx);
2078+
2079+
if (!trackNu.has_mcParticle() || !trackHad.has_mcParticle()) {
2080+
continue;
2081+
}
2082+
mQaRegistry.fill(HIST("MC/hPairFlow"), 1);
2083+
2084+
const auto particleNu = trackNu.template mcParticle_as<aod::McParticles>();
2085+
const auto particleHad = trackHad.template mcParticle_as<aod::McParticles>();
2086+
const bool truthSpeciesMatch = std::abs(particleNu.pdgCode()) == std::abs(species.settingNuPDGCode.value) &&
2087+
std::abs(particleHad.pdgCode()) == std::abs(species.settingHadPDGCode.value);
2088+
if (mc.settingRequireTruthSpecies.value && !truthSpeciesMatch) {
2089+
continue;
2090+
}
2091+
mQaRegistry.fill(HIST("MC/hPairFlow"), 2);
2092+
2093+
const bool sameMCCollision = particleNu.mcCollisionId() == particleHad.mcCollisionId();
2094+
if (mc.settingRequireSameMCCollision.value && !sameMCCollision) {
2095+
continue;
2096+
}
2097+
mQaRegistry.fill(HIST("MC/hPairFlow"), 3);
2098+
2099+
const bool matchesRecoMCCollision = collision.has_mcCollision() &&
2100+
particleNu.mcCollisionId() == collision.mcCollisionId() &&
2101+
particleHad.mcCollisionId() == collision.mcCollisionId();
2102+
if (mc.settingRequireRecoMCCollisionMatch.value && !matchesRecoMCCollision) {
2103+
continue;
2104+
}
2105+
mQaRegistry.fill(HIST("MC/hPairFlow"), 4);
2106+
2107+
if (mc.settingRequirePhysicalPrimaries.value &&
2108+
(!particleNu.isPhysicalPrimary() || !particleHad.isPhysicalPrimary())) {
2109+
continue;
2110+
}
2111+
mQaRegistry.fill(HIST("MC/hPairFlow"), 5);
2112+
2113+
HadNucandidate hadNucand;
2114+
if (!fillCandidateInfo(trackNu, trackHad, trackPair.collBracket, collisions, hadNucand, tracks, /*isMixedEvent*/ false)) {
2115+
continue;
2116+
}
2117+
mQaRegistry.fill(HIST("MC/hPairFlow"), 6);
2118+
auto selectedCollision = collisions.rawIteratorAt(hadNucand.collisionID);
2119+
fillKstar(hadNucand, selectedCollision);
2120+
fillHistograms(hadNucand);
2121+
2122+
if (output.settingFillTable && shouldFillOutputTable(hadNucand)) {
2123+
fillTable(hadNucand, selectedCollision);
2124+
fillMCTable(hadNucand, particleNu, particleHad, sameMCCollision, matchesRecoMCCollision);
2125+
}
2126+
}
2127+
}
2128+
}
2129+
PROCESS_SWITCH(HadNucleiFemto, processMC, "Process reconstructed MC same-event pairs", false);
2130+
2131+
// ==================================================================================================================
2132+
19572133
void processSameEvent(const CollisionsFull& collisions, const TrackCandidates& tracks, const aod::BCsWithTimestamps& bcs)
19582134
{
19592135
mGoodCollisions.clear();
@@ -2015,24 +2191,31 @@ struct HadNucleiFemto {
20152191
}
20162192
PROCESS_SWITCH(HadNucleiFemto, processSameEventHyper, "Process Same event", false);
20172193

2018-
void processMixedEvent(const CollisionsFull& collisions, const TrackCandidates& tracks)
2194+
void processMixedEvent(const CollisionsFull& collisions, const TrackCandidates& tracks, const aod::BCsWithTimestamps&)
20192195
{
20202196
LOG(debug) << "Processing mixed event";
2021-
mTrackPairs.clear();
20222197

20232198
for (const auto& [c1, tracks1, c2, tracks2] : mPair) {
2024-
if (!c1.sel8() || !c2.sel8()) {
2199+
if (!passesCollisionSelection</*isMC*/ false>(c1) || !passesCollisionSelection</*isMC*/ false>(c2)) {
20252200
continue;
20262201
}
20272202

20282203
mQaRegistry.fill(HIST("hNcontributor"), c1.numContrib());
20292204
mQaRegistry.fill(HIST("hVtxZ"), c1.posZ());
20302205

2206+
auto bc1 = c1.template bc_as<aod::BCsWithTimestamps>();
2207+
auto bc2 = c2.template bc_as<aod::BCsWithTimestamps>();
2208+
initCCDB</*isMC*/ false>(bc1);
2209+
mTrackPairs.clear();
20312210
pairTracksEventMixing(tracks1, tracks2);
2211+
fillPairs(collisions, tracks, /*isMixedEvent*/ true);
2212+
2213+
initCCDB</*isMC*/ false>(bc2);
2214+
mTrackPairs.clear();
20322215
pairTracksEventMixing(tracks2, tracks1);
2216+
fillPairs(collisions, tracks, /*isMixedEvent*/ true);
20332217
}
2034-
2035-
fillPairs(collisions, tracks, /*isMixedEvent*/ true);
2218+
mTrackPairs.clear();
20362219
}
20372220
PROCESS_SWITCH(HadNucleiFemto, processMixedEvent, "Process Mixed event", false);
20382221

0 commit comments

Comments
 (0)