Skip to content
69 changes: 39 additions & 30 deletions PWGCF/Femto/Core/collisionHistManager.h
Original file line number Diff line number Diff line change
Expand Up @@ -65,9 +65,17 @@ enum ColHist {
kColHistLast
};

constexpr std::string_view AnalysisDir = "Collisions/Analysis/";
constexpr std::string_view QaDir = "Collisions/QA/";
constexpr std::string_view McDir = "Collisions/MC/";
constexpr std::string_view AnalysisDir = "Analysis/";
constexpr std::string_view QaDir = "QA/";
constexpr std::string_view McDir = "MC/";

// prefixes: dedicated to the collision histograms, independent of whatever pair/triplet-type
// prefix the enclosing builder uses for its other managers -- keeps every collision histogram
// under one shared "Collision/" directory regardless of which pair/triplet type produced it.
// Standalone QA tasks have no SE/ME split, so they use PrefixQa directly.
constexpr char PrefixQa[] = "Collision/";
constexpr char PrefixSe[] = "Collision/SE/";
constexpr char PrefixMe[] = "Collision/ME/";

constexpr std::array<histmanager::HistInfo<ColHist>, kColHistLast> HistTable = {
{
Expand Down Expand Up @@ -192,6 +200,7 @@ struct ConfCollisionQaBinning : o2::framework::ConfigurableGroup {
o2::framework::ConfigurableAxis occupancy{"occupancy", {500, 0.f, 5000.f}, "Occupancy Binning"};
};

template <auto& prefix>
class CollisionHistManager
{
public:
Expand Down Expand Up @@ -270,7 +279,7 @@ class CollisionHistManager
private:
void initAnalysis(std::map<ColHist, std::vector<o2::framework::AxisSpec>> const& Specs)
{
std::string analysisDir = std::string(AnalysisDir);
std::string analysisDir = std::string(prefix) + std::string(AnalysisDir);
mHistogramRegistry->add(analysisDir + getHistNameV2(kPosZ, HistTable), getHistDesc(kPosZ, HistTable), getHistType(kPosZ, HistTable), {Specs.at(kPosZ)});
mHistogramRegistry->add(analysisDir + getHistNameV2(kMult, HistTable), getHistDesc(kMult, HistTable), getHistType(kMult, HistTable), {Specs.at(kMult)});
mHistogramRegistry->add(analysisDir + getHistNameV2(kCent, HistTable), getHistDesc(kCent, HistTable), getHistType(kCent, HistTable), {Specs.at(kCent)});
Expand All @@ -284,7 +293,7 @@ class CollisionHistManager

void initQa(std::map<ColHist, std::vector<o2::framework::AxisSpec>> const& Specs)
{
std::string qaDir = std::string(QaDir);
std::string qaDir = std::string(prefix) + std::string(QaDir);
mHistogramRegistry->add(qaDir + getHistNameV2(kPosX, HistTable), getHistDesc(kPosX, HistTable), getHistType(kPosX, HistTable), {Specs.at(kPosX)});
mHistogramRegistry->add(qaDir + getHistNameV2(kPosY, HistTable), getHistDesc(kPosY, HistTable), getHistType(kPosY, HistTable), {Specs.at(kPosY)});
mHistogramRegistry->add(qaDir + getHistNameV2(kPos, HistTable), getHistDesc(kPos, HistTable), getHistType(kPos, HistTable), {Specs.at(kPos)});
Expand All @@ -303,7 +312,7 @@ class CollisionHistManager
// reco-vs-truth correlation histograms (kReco and kMc both set)
void initMc(std::map<ColHist, std::vector<o2::framework::AxisSpec>> const& Specs)
{
std::string mcDir = std::string(McDir);
std::string mcDir = std::string(prefix) + std::string(McDir);
mHistogramRegistry->add(mcDir + getHistNameV2(kTruePosZVsPosZ, HistTable), getHistDesc(kTruePosZVsPosZ, HistTable), getHistType(kTruePosZVsPosZ, HistTable), {Specs.at(kTruePosZVsPosZ)});
mHistogramRegistry->add(mcDir + getHistNameV2(kTrueCentVsCent, HistTable), getHistDesc(kTrueCentVsCent, HistTable), getHistType(kTrueCentVsCent, HistTable), {Specs.at(kTrueCentVsCent)});
mHistogramRegistry->add(mcDir + getHistNameV2(kTrueMultVsMult, HistTable), getHistDesc(kTrueMultVsMult, HistTable), getHistType(kTrueMultVsMult, HistTable), {Specs.at(kTrueMultVsMult)});
Expand All @@ -312,7 +321,7 @@ class CollisionHistManager
// pure mc-truth collision: 1D only, no reco collision exists (kMc without kReco)
void initMcTruth(std::map<ColHist, std::vector<o2::framework::AxisSpec>> const& Specs)
{
std::string mcDir = std::string(McDir);
std::string mcDir = std::string(prefix) + std::string(McDir);
mHistogramRegistry->add(mcDir + getHistNameV2(kTruePosZ, HistTable), getHistDesc(kTruePosZ, HistTable), getHistType(kTruePosZ, HistTable), {Specs.at(kTruePosZ)});
mHistogramRegistry->add(mcDir + getHistNameV2(kTrueCent, HistTable), getHistDesc(kTrueCent, HistTable), getHistType(kTrueCent, HistTable), {Specs.at(kTrueCent)});
mHistogramRegistry->add(mcDir + getHistNameV2(kTrueMult, HistTable), getHistDesc(kTrueMult, HistTable), getHistType(kTrueMult, HistTable), {Specs.at(kTrueMult)});
Expand All @@ -321,34 +330,34 @@ class CollisionHistManager
template <typename T>
void fillAnalysis(T const& col)
{
mHistogramRegistry->fill(HIST(AnalysisDir) + HIST(getHistName(kPosZ, HistTable)), col.posZ());
mHistogramRegistry->fill(HIST(AnalysisDir) + HIST(getHistName(kMult, HistTable)), col.mult());
mHistogramRegistry->fill(HIST(AnalysisDir) + HIST(getHistName(kCent, HistTable)), col.cent());
mHistogramRegistry->fill(HIST(AnalysisDir) + HIST(getHistName(kMagField, HistTable)), col.magField());
mHistogramRegistry->fill(HIST(prefix) + HIST(AnalysisDir) + HIST(getHistName(kPosZ, HistTable)), col.posZ());
mHistogramRegistry->fill(HIST(prefix) + HIST(AnalysisDir) + HIST(getHistName(kMult, HistTable)), col.mult());
mHistogramRegistry->fill(HIST(prefix) + HIST(AnalysisDir) + HIST(getHistName(kCent, HistTable)), col.cent());
mHistogramRegistry->fill(HIST(prefix) + HIST(AnalysisDir) + HIST(getHistName(kMagField, HistTable)), col.magField());

if (mPlotEventShape) {
if constexpr (utils::HasEventShape<T>) {
mHistogramRegistry->fill(HIST(AnalysisDir) + HIST(getHistName(kQvector, HistTable)), col.qvec());
mHistogramRegistry->fill(HIST(AnalysisDir) + HIST(getHistName(kEventPlaneAngle, HistTable)), col.eventPlaneAngle());
mHistogramRegistry->fill(HIST(prefix) + HIST(AnalysisDir) + HIST(getHistName(kQvector, HistTable)), col.qvec());
mHistogramRegistry->fill(HIST(prefix) + HIST(AnalysisDir) + HIST(getHistName(kEventPlaneAngle, HistTable)), col.eventPlaneAngle());
}
}
}

template <typename T>
void fillQa(T const& col)
{
mHistogramRegistry->fill(HIST(QaDir) + HIST(getHistName(kPosX, HistTable)), col.posX());
mHistogramRegistry->fill(HIST(QaDir) + HIST(getHistName(kPosY, HistTable)), col.posY());
mHistogramRegistry->fill(HIST(QaDir) + HIST(getHistName(kPos, HistTable)), std::hypot(col.posX(), col.posY(), col.posZ()));
mHistogramRegistry->fill(HIST(QaDir) + HIST(getHistName(kSphericity, HistTable)), col.sphericity());
mHistogramRegistry->fill(HIST(QaDir) + HIST(getHistName(kOccupancy, HistTable)), col.trackOccupancyInTimeRange());
mHistogramRegistry->fill(HIST(prefix) + HIST(QaDir) + HIST(getHistName(kPosX, HistTable)), col.posX());
mHistogramRegistry->fill(HIST(prefix) + HIST(QaDir) + HIST(getHistName(kPosY, HistTable)), col.posY());
mHistogramRegistry->fill(HIST(prefix) + HIST(QaDir) + HIST(getHistName(kPos, HistTable)), std::hypot(col.posX(), col.posY(), col.posZ()));
mHistogramRegistry->fill(HIST(prefix) + HIST(QaDir) + HIST(getHistName(kSphericity, HistTable)), col.sphericity());
mHistogramRegistry->fill(HIST(prefix) + HIST(QaDir) + HIST(getHistName(kOccupancy, HistTable)), col.trackOccupancyInTimeRange());
if (mPlot2d) {
mHistogramRegistry->fill(HIST(QaDir) + HIST(getHistName(kPoszVsMult, HistTable)), col.posZ(), col.mult());
mHistogramRegistry->fill(HIST(QaDir) + HIST(getHistName(kPoszVsCent, HistTable)), col.posZ(), col.cent());
mHistogramRegistry->fill(HIST(QaDir) + HIST(getHistName(kCentVsMult, HistTable)), col.cent(), col.mult());
mHistogramRegistry->fill(HIST(QaDir) + HIST(getHistName(kMultVsSphericity, HistTable)), col.mult(), col.sphericity());
mHistogramRegistry->fill(HIST(QaDir) + HIST(getHistName(kCentVsSphericity, HistTable)), col.cent(), col.sphericity());
mHistogramRegistry->fill(HIST(QaDir) + HIST(getHistName(kFT0AvsFT0C, HistTable)), col.centFT0A(), col.centFT0C());
mHistogramRegistry->fill(HIST(prefix) + HIST(QaDir) + HIST(getHistName(kPoszVsMult, HistTable)), col.posZ(), col.mult());
mHistogramRegistry->fill(HIST(prefix) + HIST(QaDir) + HIST(getHistName(kPoszVsCent, HistTable)), col.posZ(), col.cent());
mHistogramRegistry->fill(HIST(prefix) + HIST(QaDir) + HIST(getHistName(kCentVsMult, HistTable)), col.cent(), col.mult());
mHistogramRegistry->fill(HIST(prefix) + HIST(QaDir) + HIST(getHistName(kMultVsSphericity, HistTable)), col.mult(), col.sphericity());
mHistogramRegistry->fill(HIST(prefix) + HIST(QaDir) + HIST(getHistName(kCentVsSphericity, HistTable)), col.cent(), col.sphericity());
mHistogramRegistry->fill(HIST(prefix) + HIST(QaDir) + HIST(getHistName(kFT0AvsFT0C, HistTable)), col.centFT0A(), col.centFT0C());
}
}

Expand All @@ -360,18 +369,18 @@ class CollisionHistManager
return;
}
auto mcCol = col.template fMcCol_as<T2>();
mHistogramRegistry->fill(HIST(McDir) + HIST(getHistName(kTruePosZVsPosZ, HistTable)), mcCol.posZ(), col.posZ());
mHistogramRegistry->fill(HIST(McDir) + HIST(getHistName(kTrueCentVsCent, HistTable)), mcCol.cent(), col.cent());
mHistogramRegistry->fill(HIST(McDir) + HIST(getHistName(kTrueMultVsMult, HistTable)), mcCol.mult(), col.mult());
mHistogramRegistry->fill(HIST(prefix) + HIST(McDir) + HIST(getHistName(kTruePosZVsPosZ, HistTable)), mcCol.posZ(), col.posZ());
mHistogramRegistry->fill(HIST(prefix) + HIST(McDir) + HIST(getHistName(kTrueCentVsCent, HistTable)), mcCol.cent(), col.cent());
mHistogramRegistry->fill(HIST(prefix) + HIST(McDir) + HIST(getHistName(kTrueMultVsMult, HistTable)), mcCol.mult(), col.mult());
}

// pure mc-truth collision: 'col' here IS the truth collision, no reco object exists
template <typename T>
void fillMc(T const& col)
{
mHistogramRegistry->fill(HIST(McDir) + HIST(getHistName(kTruePosZ, HistTable)), col.posZ());
mHistogramRegistry->fill(HIST(McDir) + HIST(getHistName(kTrueCent, HistTable)), col.cent());
mHistogramRegistry->fill(HIST(McDir) + HIST(getHistName(kTrueMult, HistTable)), col.mult());
mHistogramRegistry->fill(HIST(prefix) + HIST(McDir) + HIST(getHistName(kTruePosZ, HistTable)), col.posZ());
mHistogramRegistry->fill(HIST(prefix) + HIST(McDir) + HIST(getHistName(kTrueCent, HistTable)), col.cent());
mHistogramRegistry->fill(HIST(prefix) + HIST(McDir) + HIST(getHistName(kTrueMult, HistTable)), col.mult());
}

o2::framework::HistogramRegistry* mHistogramRegistry = nullptr;
Expand Down
44 changes: 29 additions & 15 deletions PWGCF/Femto/Core/dndetaBuilder.h
Original file line number Diff line number Diff line change
Expand Up @@ -173,9 +173,13 @@ class DndetaBuilder
// data
// ---------------------------------------------------------------------------

/// \param triggerKinematic k*/Q3 of the triggering pair/triplet (minimum over the event), written
/// by the trigger callable as a side effect on itself returning true; -1 (or whatever the
/// caller initializes it to) for minimum bias. Read only after trigger() has run, i.e.
/// only meaningful once this function returns true.
/// \return true if the collision was selected (including the trigger)
template <typename T1, typename T2, typename T3, typename T4, typename T5>
bool processData(T1 const& col, T2& partitionGlobal, T3& partitionItsOnly, T4& cache, T5 const& trigger)
bool processData(T1 const& col, T2& partitionGlobal, T3& partitionItsOnly, T4& cache, T5 const& trigger, float triggerKinematic = -1.f)
{
if (!selectCollision(col, trigger)) {
return false;
Expand All @@ -194,15 +198,15 @@ class DndetaBuilder
}
nch++;
mTrackHistManagerGlobal.template fill<modes::Mode::kReco_Qa>(track, tracksGlobal);
mHistManager.fillDataTrack(posZ, cent, track.eta(), track.phi(), true);
mHistManager.fillDataTrack(posZ, cent, track.eta(), track.phi(), true, triggerKinematic);
}
for (auto const& track : tracksItsOnly) {
if (!isTrackAccepted(track, false)) {
continue;
}
nch++;
mTrackHistManagerItsOnly.template fill<modes::Mode::kReco_Qa>(track, tracksItsOnly);
mHistManager.fillDataTrack(posZ, cent, track.eta(), track.phi(), false);
mHistManager.fillDataTrack(posZ, cent, track.eta(), track.phi(), false, triggerKinematic);
}
if (mDoCorrelation) {
mHistManager.fillCorrelation(static_cast<float>(nch), col.mult(), col.multFT0A(), col.multFT0C());
Expand Down Expand Up @@ -234,7 +238,11 @@ class DndetaBuilder
/// - gen: efficiency histograms for all generated collisions, loss histograms for the generated event class
/// (optional INEL>0, optional generator-level trigger); the loss block's reco requirement is Sel8
/// only (no vtxZ, see isSelectedByMaskOnly), the vtxZ cut only gates the efficiency block above
template <typename T1, typename T2, typename T3, typename T4, typename T5, typename T6, typename T7, typename T8, typename T9, typename T10, typename T11>
/// \param triggerKinematicScratch scratch variable the recoTrigger callable writes the
/// triggering k*/Q3 into as a side effect when it returns true (minimum over the event);
/// read back immediately after each recoTrigger() call, so its value right before the
/// call is irrelevant and it is freely reused across collisions
template <typename T1, typename T2, typename T3, typename T4, typename T5, typename T6, typename T7, typename T8, typename T9, typename T10, typename T11, typename T12>
void processMc(T1 const& cols,
T2 const& mcCols,
T3& partitionGlobal,
Expand All @@ -245,7 +253,8 @@ class DndetaBuilder
T8 const& perMcColParticles,
T9& pdg,
T10 const& recoTrigger,
T11 const& genTrigger)
T11 const& genTrigger,
T12& triggerKinematicScratch)
{
const auto nMcCols = static_cast<std::size_t>(mcCols.size());

Expand All @@ -270,6 +279,7 @@ class DndetaBuilder
// not implemented (every Run3 event with a vertex is already within the stored range), so the loss-block
// reco requirement must not carry the vtxZ cut like the Ctrk-side "selvtx" one does
std::vector<int64_t> selectedReco(nMcCols, -1);
std::vector<float> selectedRecoTriggerKinematic(nMcCols, -1.f);
std::vector<bool> hasSel8Reco(nMcCols, false);
for (auto const& col : cols) {
if (!col.has_fMcCol() || bestReco[static_cast<std::size_t>(col.fMcColId())] != col.globalIndex()) {
Expand All @@ -284,10 +294,11 @@ class DndetaBuilder
continue;
}
selectedReco[static_cast<std::size_t>(col.fMcColId())] = col.globalIndex();
selectedRecoTriggerKinematic[static_cast<std::size_t>(col.fMcColId())] = triggerKinematicScratch;
if (mDoMcEfficiency) {
auto tracksGlobal = partitionGlobal->sliceByCached(o2::aod::femtobase::stored::fColId, col.globalIndex(), cache);
auto tracksItsOnly = partitionItsOnly->sliceByCached(o2::aod::femtobase::stored::fColId, col.globalIndex(), cache);
fillRecoMc(col, tracksGlobal, tracksItsOnly, mcParticles, mcMothers);
fillRecoMc(col, tracksGlobal, tracksItsOnly, mcParticles, mcMothers, triggerKinematicScratch);
}
}

Expand All @@ -300,13 +311,14 @@ class DndetaBuilder
// efficiency: centrality of the selected reconstructed collision (underflow if there is none)
if (mDoMcEfficiency) {
const float recoCent = hasSelectedReco ? cols.rawIteratorAt(selectedReco[mcColIndex]).cent() : -1.f;
const float triggerKinematic = hasSelectedReco ? selectedRecoTriggerKinematic[mcColIndex] : -1.f;
mHistManager.fillGenCollision(mcCol.posZ(), recoCent, hasSelectedReco);
for (auto const& particle : particles) {
if (!isGenParticleSelected(particle, pdg)) {
continue;
}
mHistManager.fillGenParticle(mcCol.posZ(), recoCent, particle.eta(), particle.phi(), genType(particle.pdgCode()), hasSelectedReco,
ptVariationWeight(particle.pt(), mPtVariationUp), ptVariationWeight(particle.pt(), mPtVariationDown));
ptVariationWeight(particle.pt(), mPtVariationUp), ptVariationWeight(particle.pt(), mPtVariationDown), triggerKinematic);
}
}

Expand Down Expand Up @@ -492,20 +504,22 @@ class DndetaBuilder
}

template <typename T1, typename T2, typename T3, typename T4>
void fillRecoMcTrack(T1 const& col, T2 const& track, T3 const& /*mcParticles*/, T4 const& /*mcMothers*/, std::vector<int64_t>& usedLabels)
void fillRecoMcTrack(T1 const& col, T2 const& track, T3 const& /*mcParticles*/, T4 const& /*mcMothers*/, std::vector<int64_t>& usedLabels, float triggerKinematic)
{
const float posZ = col.posZ();
const float cent = col.cent();
if (!track.has_fMcParticle()) {
mHistManager.fillRecoTrack(posZ, cent, track.eta(), track.phi(), dndetahistmanager::kRecoBkg);
mHistManager.fillRecoTrack(posZ, cent, track.eta(), track.phi(), dndetahistmanager::kRecoBkg, triggerKinematic);
return;
}
auto mcParticle = track.template fMcParticle_as<T3>();
// tracks matched to a particle from another mc collision are skipped (as in the legacy task)
// a track whose matched mc particle belongs to a different mc collision is a mis-association
// (e.g. an ambiguous/pileup-adjacent match); counting it here would attribute it to the wrong
// collision's efficiency, so it is skipped
if (mcParticle.fMcColId() != col.fMcColId()) {
return;
}
mHistManager.fillRecoTrack(posZ, cent, track.eta(), track.phi(), dndetahistmanager::kRecoAll);
mHistManager.fillRecoTrack(posZ, cent, track.eta(), track.phi(), dndetahistmanager::kRecoAll, triggerKinematic);

dndetahistmanager::RecoType type = dndetahistmanager::kRecoSecondary;
if (isPhysicalPrimary(mcParticle)) {
Expand All @@ -521,25 +535,25 @@ class DndetaBuilder
type = dndetahistmanager::kRecoFake;
}
usedLabels.push_back(track.fMcParticleId());
mHistManager.fillRecoTrack(posZ, cent, mcParticle.eta(), mcParticle.phi(), type);
mHistManager.fillRecoTrack(posZ, cent, mcParticle.eta(), mcParticle.phi(), type, triggerKinematic);
}

template <typename T1, typename T2, typename T3, typename T4, typename T5>
void fillRecoMc(T1 const& col, T2 const& tracksGlobal, T3 const& tracksItsOnly, T4 const& mcParticles, T5 const& mcMothers)
void fillRecoMc(T1 const& col, T2 const& tracksGlobal, T3 const& tracksItsOnly, T4 const& mcParticles, T5 const& mcMothers, float triggerKinematic)
{
mHistManager.fillRecoCollision(col.posZ(), col.cent());
// labels are tracked over both partitions, so a fake is also found if the other track comes from the other partition
std::vector<int64_t> usedLabels;
for (auto const& track : tracksGlobal) {
if (isTrackAccepted(track, true)) {
mTrackHistManagerGlobal.template fill<modes::Mode::kReco_Qa>(track, tracksGlobal);
fillRecoMcTrack(col, track, mcParticles, mcMothers, usedLabels);
fillRecoMcTrack(col, track, mcParticles, mcMothers, usedLabels, triggerKinematic);
}
}
for (auto const& track : tracksItsOnly) {
if (isTrackAccepted(track, false)) {
mTrackHistManagerItsOnly.template fill<modes::Mode::kReco_Qa>(track, tracksItsOnly);
fillRecoMcTrack(col, track, mcParticles, mcMothers, usedLabels);
fillRecoMcTrack(col, track, mcParticles, mcMothers, usedLabels, triggerKinematic);
}
}
}
Expand Down
Loading
Loading