Skip to content

Commit dd055b2

Browse files
authored
Update HadNucleiFemto.cxx
1 parent 46423fd commit dd055b2

1 file changed

Lines changed: 291 additions & 4 deletions

File tree

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

Lines changed: 291 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -43,11 +43,16 @@
4343
#include <Framework/AnalysisTask.h>
4444
#include <Framework/Array2D.h>
4545
#include <Framework/BinningPolicy.h>
46+
#include <Framework/CallbackService.h>
4647
#include <Framework/Configurable.h>
48+
#include <Framework/DeviceSpec.h>
49+
#include <Framework/EndOfStreamContext.h>
4750
#include <Framework/HistogramRegistry.h>
4851
#include <Framework/HistogramSpec.h>
4952
#include <Framework/InitContext.h>
5053
#include <Framework/OutputObjHeader.h>
54+
#include <Framework/ProcessingContext.h>
55+
#include <Framework/TimingInfo.h>
5156
#include <Framework/runDataProcessing.h>
5257
#include <MathUtils/BetheBlochAleph.h>
5358
#include <MathUtils/Primitive2D.h>
@@ -62,6 +67,7 @@
6267
#include <TMCProcess.h>
6368
#include <TPDGCode.h>
6469
#include <TString.h>
70+
#include <TUUID.h>
6571

6672
#include <algorithm>
6773
#include <array>
@@ -322,6 +328,16 @@ struct HadNucleiFemto {
322328
Produces<aod::HadronHyperTable> mOutputHadHyperDataTable;
323329
Produces<aod::HadronHyperTableMC> mOutputHadHyperMCTable;
324330
Produces<aod::HadronNucleiMult> mOutputMultiplicityTable;
331+
Produces<aod::HadNucleiMEDiag> mOutputMEDiagnostic;
332+
Produces<aod::HadNucleiMEDiagLinks> mOutputMEDiagnosticLinks;
333+
334+
struct : o2::framework::ConfigurableGroup {
335+
std::string prefix{"meDiagnostic"};
336+
Configurable<bool> saveRejectedEvents{"saveRejectedEvents", false, "Save rejected collision rows too (potentially very large); rejection counters are always recorded"};
337+
Configurable<bool> savePartnerLinks{"savePartnerLinks", true, "Save one diagnostic row per undirected mixed-event combination"};
338+
Configurable<std::string> productionTag{"productionTag", "unspecified", "Production label recorded with diagnostic instance UUID"};
339+
Configurable<std::string> codeTag{"codeTag", "unspecified", "Exact source/build revision label for the diagnostic run"};
340+
} meDiagnostic;
325341

326342
struct : o2::framework::ConfigurableGroup {
327343
// cppcheck-suppress unusedStructMember
@@ -647,6 +663,29 @@ struct HadNucleiFemto {
647663
int mMixingRunNumber{-1};
648664
int mHyperMixingRunNumber{-1};
649665
int mHyperMCMixingRunNumber{-1};
666+
// A separate state machine: the ordinary ME pools and counters are untouched.
667+
std::unordered_map<int, std::deque<BufferedCollision>> mDiagnosticPools;
668+
TimingInfo mDiagnosticTiming;
669+
std::string mDiagnosticUUID;
670+
uint64_t mDiagnosticInstanceHi{0}, mDiagnosticInstanceLo{0};
671+
uint64_t mDiagnosticCall{0}, mDiagnosticEpoch{0}, mDiagnosticBuffered{0};
672+
uint64_t mDiagnosticExpected{0}, mDiagnosticAccepted{0}, mDiagnosticCombinations{0};
673+
int64_t mNextDiagnosticEventId{0};
674+
int mDiagnosticRun{-1};
675+
static constexpr uint8_t DiagnosticEventKind = 3;
676+
struct MEDiagnosticRow {
677+
uint8_t kind{DiagnosticEventKind}, status{0};
678+
int run{-1}, previousRun{-1};
679+
int64_t eventId{-1}, inputCollisionId{-1};
680+
uint64_t globalBC{0}, timestamp{0};
681+
float z{0.f}, centrality{0.f}, multiplicity{0.f};
682+
uint32_t numContrib{0}, nNuclei{0}, nHadrons{0};
683+
int poolBin{-1};
684+
uint32_t before{0}, after{0};
685+
uint64_t expected{0}, accepted{0}, combinations{0};
686+
uint64_t inputCollisions{0}, passedEventSelection{0}, passedZorro{0};
687+
};
688+
HistogramRegistry mDiagnosticRegistry{"MEDiagnostic", {}, OutputObjHandlingPolicy::AnalysisObject, false, true};
650689
o2::vertexing::DCAFitterN<2> mFitter;
651690

652691
int mRunNumber{0};
@@ -769,10 +808,10 @@ struct HadNucleiFemto {
769808
{"purity/h2NsigmaNuTPC_preselecComp", "NsigmaNu TPC distribution; #it{p}_{T} (GeV/#it{c}); n#sigma_{TPC}(Nu)", {HistType::kTH2F, {{280, -7.0f, 7.0f}, {400, -10.0f, 10.0f}}}},
770809
{"purity/h2NSigmaNuITS_preselection", "NsigmaNu ITS distribution; signed #it{p}_{T} (GeV/#it{c}); n#sigma_{ITS} Nu", {HistType::kTH2F, {{280, -7.0f, 7.0f}, {120, -5.0f, 5.0f}}}},
771810
{"purity/h2NsigmaNuTOF_preselection", "NsigmaNu TOF distribution; #it{p}_{T} (GeV/#it{c}); n#sigma_{TOF}(Nu)", {HistType::kTH2F, {{280, -7.0f, 7.0f}, {400, -10.0f, 10.0f}}}},
772-
{"purity/h2NsigmaNuComb_preselection", "NsigmaNu TPCTOF comb distribution; #it{p}_{T} (GeV/#it{c}); n#sigma_{comb}(Nu)", {HistType::kTH2F, {{280, -7.0f, 7.0f}, {100, 0.0f, 7.0f}}}},
811+
{"purity/h2NsigmaNuComb_preselection", "NsigmaNu TPCTOF comb distribution; #it{p}_{T} (GeV/#it{c}); n#sigma_{comb}(Nu)", {HistType::kTH2F, {{280, -7.0f, 7.0f}, {100, 0.0f, 10.0f}}}},
773812
{"purity/h2NsigmaHadTPC_preselection", "NsigmaNu TPC distribution; #it{p}_{T} (GeV/#it{c}); n#sigma_{TPC}(Nu)", {HistType::kTH2F, {{280, -7.0f, 7.0f}, {400, -10.0f, 10.0f}}}},
774813
{"purity/h2NsigmaHadTOF_preselection", "NsigmaHad TOF distribution; #iit{p}_{T} (GeV/#it{c}); n#sigma_{TOF}(p)", {HistType::kTH2F, {{280, -7.0f, 7.0f}, {400, -10.0f, 10.0f}}}},
775-
{"purity/h2NsigmaHadComb_preselection", "NsigmaHad TPCTOF comb distribution; #it{p}_{T} (GeV/#it{c}); n#sigma_{comb}(had)", {HistType::kTH2F, {{280, -7.0f, 7.0f}, {100, 0.0f, 7.0f}}}},
814+
{"purity/h2NsigmaHadComb_preselection", "NsigmaHad TPCTOF comb distribution; #it{p}_{T} (GeV/#it{c}); n#sigma_{comb}(had)", {HistType::kTH2F, {{280, -7.0f, 7.0f}, {100, 0.0f, 10.0f}}}},
776815

777816
// Hypertriton
778817
{"hHe3TPCnsigma", "NsigmaHe3 TPC distribution; #it{p}_{T} (GeV/#it{c}); n#sigma_{TPC}(He3)", {HistType::kTH2F, {{280, -7.0f, 7.0f}, {200, -5.0f, 5.0f}}}},
@@ -800,8 +839,44 @@ struct HadNucleiFemto {
800839
false,
801840
true};
802841

803-
void init(o2::framework::InitContext&)
804-
{
842+
void init(o2::framework::InitContext& initContext)
843+
{
844+
if (doprocessMixedEventDiagnostic) {
845+
// AnalysisTask does not automatically invoke a task endOfStream method.
846+
// CallbackService appends this callback without replacing output flushing.
847+
initContext.services().get<CallbackService>().set<CallbackService::Id::EndOfStream>(
848+
[this](EndOfStreamContext& context) { endOfStream(context); });
849+
if (doprocessMixedEvent || doprocessSameEvent || doprocessMC || doprocessHyper || doprocessMCHyper || doprocessPurity || doprocessDcaFractionData || doprocessDcaFractionPurityMC) {
850+
LOG(fatal) << "processMixedEventDiagnostic must run alone; disable the other HadNucleiFemto processes to avoid shared selection-QA contamination";
851+
}
852+
if (eventMixing.settingNoMixedEvents.value <= 0) {
853+
LOG(fatal) << "ME diagnostic requires a positive mixing depth";
854+
}
855+
mDiagnosticUUID = TUUID().AsString();
856+
std::string uuidWords = mDiagnosticUUID;
857+
uuidWords.erase(std::remove(uuidWords.begin(), uuidWords.end(), '-'), uuidWords.end());
858+
mDiagnosticInstanceHi = std::stoull(uuidWords.substr(0, 16), nullptr, 16);
859+
mDiagnosticInstanceLo = std::stoull(uuidWords.substr(16, 16), nullptr, 16);
860+
const auto& device = initContext.services().get<DeviceSpec const>();
861+
const char* alienJob = std::getenv("ALIEN_PROC_ID");
862+
LOG(info) << "ME_DIAG INIT instance=" << mDiagnosticUUID << " device=" << device.id
863+
<< " rank=" << device.rank << " inputTimesliceId=" << device.inputTimesliceId
864+
<< " job=" << (alienJob ? alienJob : "unknown")
865+
<< " production=" << meDiagnostic.productionTag.value << " code=" << meDiagnostic.codeTag.value
866+
<< " depth=" << eventMixing.settingNoMixedEvents.value
867+
<< " requireBoth=" << eventMixing.settingRequireBothSpeciesForMixing.value
868+
<< " multiplicityBinning=" << eventMixing.settingUseMultiplicityBinning.value
869+
<< " applyCPR=" << CPR.settingApplyClosePairRejection.value;
870+
for (const auto& name : {"hLifecycle", "hSelection"}) {
871+
mDiagnosticRegistry.add(name, name, HistType::kTH1D, {{6, -0.5, 5.5}});
872+
}
873+
mDiagnosticRegistry.add("hPoolOccupancy", "Prior partners;partners;events", HistType::kTH1D, {{101, -0.5, 100.5}});
874+
mDiagnosticRegistry.add("hPairs", "0=expected,1=accepted,2=CPR rejected;counter;pairs", HistType::kTH1D, {{3, -0.5, 2.5}});
875+
// The global sink keys by task + histogram name, not registry directory.
876+
for (const auto& name : {"hMEDiagKStar_LS_M", "hMEDiagKStar_LS_A", "hMEDiagKStar_US_M", "hMEDiagKStar_US_A"}) {
877+
mDiagnosticRegistry.add(name, name, HistType::kTH1D, {{300, 0., 3.}});
878+
}
879+
}
805880
const bool processHyperPairs = doprocessHyper || doprocessMCHyper;
806881
if (processHyperPairs && hadHyper.maxOutputKstar.value == 0.f) {
807882
LOG(fatal) << "Hadron-hypertriton mode requires a nonzero output k* range";
@@ -3332,6 +3407,218 @@ struct HadNucleiFemto {
33323407
}
33333408
PROCESS_SWITCH(HadNucleiFemto, processMixedEvent, "Process Mixed event", false);
33343409

3410+
void run(ProcessingContext& context)
3411+
{
3412+
if (doprocessMixedEventDiagnostic) {
3413+
// TimingInfo is stream-scoped; access it through the active context.
3414+
mDiagnosticTiming = context.services().get<TimingInfo>();
3415+
}
3416+
}
3417+
3418+
void writeMEDiagnostic(const MEDiagnosticRow& row)
3419+
{
3420+
mOutputMEDiagnostic(mDiagnosticInstanceHi, mDiagnosticInstanceLo, mDiagnosticCall,
3421+
static_cast<uint64_t>(mDiagnosticTiming.timeslice), mDiagnosticTiming.tfCounter, mDiagnosticTiming.runNumber,
3422+
row.kind, row.status, mDiagnosticEpoch, row.run, row.previousRun, row.eventId,
3423+
row.inputCollisionId, row.globalBC, row.timestamp, row.z, row.centrality, row.multiplicity,
3424+
row.numContrib, row.nNuclei, row.nHadrons, row.poolBin, row.before, row.after,
3425+
mDiagnosticBuffered, static_cast<uint32_t>(mDiagnosticPools.size()),
3426+
row.expected, row.accepted, row.expected - row.accepted, row.combinations,
3427+
row.inputCollisions, row.passedEventSelection, row.passedZorro);
3428+
if (row.kind != DiagnosticEventKind) {
3429+
mDiagnosticRegistry.fill(HIST("hLifecycle"), row.kind);
3430+
}
3431+
}
3432+
3433+
bool fillDiagnosticMixedPair(const BufferedTrack& nucleus, const BufferedTrack& hadron)
3434+
{
3435+
// Same CPR predicate and kinematics as fillMixedPair; separate final QA,
3436+
// and never write the ordinary pair/multiplicity tables from this process.
3437+
if (isClosePair(nucleus, hadron, /*fillQA*/ true)) {
3438+
return false;
3439+
}
3440+
HadNucandidate candidate;
3441+
fillBufferedCandidateInfo(nucleus, hadron, candidate);
3442+
if (candidate.recoPtNu() > 0) {
3443+
if (candidate.isBkgUS) {
3444+
mDiagnosticRegistry.fill(HIST("hMEDiagKStar_US_M"), candidate.kstar);
3445+
} else {
3446+
mDiagnosticRegistry.fill(HIST("hMEDiagKStar_LS_M"), candidate.kstar);
3447+
}
3448+
} else if (candidate.isBkgUS) {
3449+
mDiagnosticRegistry.fill(HIST("hMEDiagKStar_US_A"), candidate.kstar);
3450+
} else {
3451+
mDiagnosticRegistry.fill(HIST("hMEDiagKStar_LS_A"), candidate.kstar);
3452+
}
3453+
return true;
3454+
}
3455+
3456+
void processMixedEventDiagnostic(const CollisionsFullWithPVMult& collisions, const TrackCandidates& tracks, const aod::BCsWithTimestamps&)
3457+
{
3458+
const CentralityBinningType centralityBinningPolicy{{axisVertex, axisCentrality}, true};
3459+
const MultiplicityBinningType multiplicityBinningPolicy{{axisVertex, axisMultiplicity}, true};
3460+
++mDiagnosticCall;
3461+
MEDiagnosticRow call;
3462+
call.run = mDiagnosticRun;
3463+
call.inputCollisions = collisions.size();
3464+
if (mDiagnosticCall == 1) {
3465+
call.kind = 0; // Instance initialization, emitted in a valid table context.
3466+
writeMEDiagnostic(call);
3467+
}
3468+
call.kind = 1;
3469+
writeMEDiagnostic(call);
3470+
LOG(info) << "ME_DIAG BEGIN instance=" << mDiagnosticUUID << " call=" << mDiagnosticCall
3471+
<< " timeslice=" << mDiagnosticTiming.timeslice << " tfCounter=" << mDiagnosticTiming.tfCounter
3472+
<< " run=" << mDiagnosticRun << " buffered=" << mDiagnosticBuffered << " bins=" << mDiagnosticPools.size();
3473+
3474+
for (const auto& collision : collisions) {
3475+
const auto bc = collision.template bc_as<aod::BCsWithTimestamps>();
3476+
MEDiagnosticRow row;
3477+
row.run = bc.runNumber();
3478+
row.inputCollisionId = collision.globalIndex();
3479+
row.globalBC = bc.globalBC();
3480+
row.timestamp = bc.timestamp();
3481+
row.z = collision.posZ();
3482+
row.centrality = collision.centFT0C();
3483+
row.multiplicity = collision.multNTracksPV();
3484+
row.numContrib = collision.numContrib();
3485+
mDiagnosticRegistry.fill(HIST("hSelection"), 0);
3486+
if (!passesEventSelection</*isMC*/ false>(collision)) {
3487+
row.status = 1;
3488+
if (meDiagnostic.saveRejectedEvents) {
3489+
writeMEDiagnostic(row);
3490+
}
3491+
continue;
3492+
}
3493+
++call.passedEventSelection;
3494+
mDiagnosticRegistry.fill(HIST("hSelection"), 1);
3495+
if (!passesZorroSelection(collision)) {
3496+
row.status = 2;
3497+
if (meDiagnostic.saveRejectedEvents) {
3498+
writeMEDiagnostic(row);
3499+
}
3500+
continue;
3501+
}
3502+
++call.passedZorro;
3503+
mDiagnosticRegistry.fill(HIST("hSelection"), 2);
3504+
// Match the ordinary ME ordering: observe the run before species/bin cuts.
3505+
if (mDiagnosticRun != bc.runNumber()) {
3506+
MEDiagnosticRow reset = row;
3507+
reset.kind = 2;
3508+
reset.status = mDiagnosticRun < 0 ? 0 : 1;
3509+
reset.previousRun = mDiagnosticRun;
3510+
++mDiagnosticEpoch;
3511+
writeMEDiagnostic(reset); // Contains the pre-clear buffered totals.
3512+
LOG(info) << "ME_DIAG CLEAR instance=" << mDiagnosticUUID << " call=" << mDiagnosticCall
3513+
<< " reason=" << (mDiagnosticRun < 0 ? "first-run" : "run-change")
3514+
<< " oldRun=" << mDiagnosticRun << " newRun=" << bc.runNumber()
3515+
<< " epoch=" << mDiagnosticEpoch << " buffered=" << mDiagnosticBuffered
3516+
<< " bins=" << mDiagnosticPools.size();
3517+
mDiagnosticPools.clear();
3518+
mDiagnosticBuffered = 0;
3519+
mDiagnosticRun = bc.runNumber();
3520+
}
3521+
BufferedCollision current;
3522+
current.eventId = mNextDiagnosticEventId++;
3523+
current.posZ = collision.posZ();
3524+
current.numContrib = collision.numContrib();
3525+
current.centFT0C = collision.centFT0C();
3526+
current.multFT0C = collision.multFT0C();
3527+
row.eventId = current.eventId;
3528+
auto tracksThisCollision = tracks.sliceBy(mPerCol, collision.globalIndex());
3529+
tracksThisCollision.bindExternalIndices(&tracks);
3530+
for (const auto& track : tracksThisCollision) {
3531+
if (selectTrackNu(track) && selectionPIDNu(track)) {
3532+
current.nuclei.push_back(makeBufferedTrack(track, true, current.eventId));
3533+
}
3534+
if (selectTrackHadron(track) && selectionPIDHadron(track)) {
3535+
current.hadrons.push_back(makeBufferedTrack(track, false, current.eventId));
3536+
}
3537+
}
3538+
row.nNuclei = current.nuclei.size();
3539+
row.nHadrons = current.hadrons.size();
3540+
if (eventMixing.settingRequireBothSpeciesForMixing.value && (current.nuclei.empty() || current.hadrons.empty())) {
3541+
row.status = 3;
3542+
if (meDiagnostic.saveRejectedEvents) {
3543+
writeMEDiagnostic(row);
3544+
}
3545+
continue;
3546+
}
3547+
mDiagnosticRegistry.fill(HIST("hSelection"), 3);
3548+
row.poolBin = eventMixing.settingUseMultiplicityBinning.value
3549+
? multiplicityBinningPolicy.getBin(std::make_tuple(collision.posZ(), collision.multNTracksPV()))
3550+
: centralityBinningPolicy.getBin(std::make_tuple(collision.posZ(), collision.centFT0C()));
3551+
if (row.poolBin < 0) {
3552+
row.status = 4;
3553+
if (meDiagnostic.saveRejectedEvents) {
3554+
writeMEDiagnostic(row);
3555+
}
3556+
continue;
3557+
}
3558+
mDiagnosticRegistry.fill(HIST("hSelection"), 4);
3559+
auto& pool = mDiagnosticPools[row.poolBin];
3560+
row.before = pool.size();
3561+
row.combinations = pool.size();
3562+
mDiagnosticRegistry.fill(HIST("hPoolOccupancy"), row.before);
3563+
for (const auto& partner : pool) {
3564+
const uint64_t expected = static_cast<uint64_t>(current.nuclei.size()) * partner.hadrons.size() +
3565+
static_cast<uint64_t>(partner.nuclei.size()) * current.hadrons.size();
3566+
uint64_t accepted = 0;
3567+
for (const auto& nucleus : current.nuclei) {
3568+
for (const auto& hadron : partner.hadrons) {
3569+
accepted += fillDiagnosticMixedPair(nucleus, hadron);
3570+
}
3571+
}
3572+
for (const auto& nucleus : partner.nuclei) {
3573+
for (const auto& hadron : current.hadrons) {
3574+
accepted += fillDiagnosticMixedPair(nucleus, hadron);
3575+
}
3576+
}
3577+
row.expected += expected;
3578+
row.accepted += accepted;
3579+
if (meDiagnostic.savePartnerLinks) {
3580+
mOutputMEDiagnosticLinks(mDiagnosticInstanceHi, mDiagnosticInstanceLo, mDiagnosticCall,
3581+
mDiagnosticEpoch, mDiagnosticRun, current.eventId, partner.eventId, row.poolBin, expected, accepted);
3582+
}
3583+
}
3584+
if (pool.size() >= static_cast<size_t>(eventMixing.settingNoMixedEvents.value)) {
3585+
pool.pop_front();
3586+
} else {
3587+
++mDiagnosticBuffered;
3588+
}
3589+
pool.push_back(std::move(current));
3590+
row.after = pool.size();
3591+
writeMEDiagnostic(row);
3592+
call.expected += row.expected;
3593+
call.accepted += row.accepted;
3594+
call.combinations += row.combinations;
3595+
mDiagnosticRegistry.fill(HIST("hPairs"), 0., static_cast<double>(row.expected));
3596+
mDiagnosticRegistry.fill(HIST("hPairs"), 1., static_cast<double>(row.accepted));
3597+
mDiagnosticRegistry.fill(HIST("hPairs"), 2., static_cast<double>(row.expected - row.accepted));
3598+
}
3599+
call.kind = 4;
3600+
call.run = mDiagnosticRun;
3601+
writeMEDiagnostic(call);
3602+
mDiagnosticExpected += call.expected;
3603+
mDiagnosticAccepted += call.accepted;
3604+
mDiagnosticCombinations += call.combinations;
3605+
LOG(info) << "ME_DIAG END instance=" << mDiagnosticUUID << " call=" << mDiagnosticCall
3606+
<< " buffered=" << mDiagnosticBuffered << " bins=" << mDiagnosticPools.size()
3607+
<< " combinations=" << call.combinations << " expected=" << call.expected << " accepted=" << call.accepted;
3608+
}
3609+
PROCESS_SWITCH(HadNucleiFemto, processMixedEventDiagnostic, "Independent ME pool lifecycle diagnostic (run alone)", false);
3610+
3611+
void endOfStream(EndOfStreamContext&)
3612+
{
3613+
if (doprocessMixedEventDiagnostic) {
3614+
// EOS table writes are not guaranteed to be flushed; final snapshot in log.
3615+
LOG(info) << "ME_DIAG EOS instance=" << mDiagnosticUUID << " calls=" << mDiagnosticCall
3616+
<< " epochs=" << mDiagnosticEpoch << " buffered=" << mDiagnosticBuffered
3617+
<< " bins=" << mDiagnosticPools.size() << " combinations=" << mDiagnosticCombinations
3618+
<< " expected=" << mDiagnosticExpected << " accepted=" << mDiagnosticAccepted;
3619+
}
3620+
}
3621+
33353622
// Produce the data distributions fitted offline with the MC templates. All
33363623
// nominal quality and PID selections are applied; only DCA is relaxed.
33373624
void processDcaFractionData(const CollisionsFull& collisions, const TrackCandidates& tracks, const aod::BCsWithTimestamps& bcs)

0 commit comments

Comments
 (0)