Skip to content

[PWGDQ] Fixing mu calculation and moving to table reader #18732

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Merged
iarsene merged 2 commits into from
Nov 13, 2026
Merged

[PWGDQ] Fixing mu calculation and moving to table reader #18732

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
2 commits
Select commit Hold shift + click to select a range
5e0d759
Fixing mu calculation and moving to table reader
lucamicheletti Nov 13, 2026
7a5fedb
Fix clang
lucamicheletti Nov 13, 2026
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions PWGDQ/Core/HistogramsLibrary.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -59,6 +59,7 @@ void o2::aod::dqhistograms::DefineHistograms(HistogramManager* hm, const char* h
hm->AddHistogram(histClass, "VtxZ", "Vtx Z", false, 60, -15.0, 15.0, VarManager::kVtxZ);
hm->AddHistogram(histClass, "VtxZ_Run", "Vtx Z", true, 1, -0.5, 0.5, VarManager::kRunNo, 60, -15.0, 15.0, VarManager::kVtxZ, 1, 0, 1, VarManager::kNothing, "", "", "", VarManager::kNothing, VarManager::kNothing, false, true);
hm->AddHistogram(histClass, "BC", "Event per BC", false, 3564, 0.0, 3564.0, VarManager::kBCOrbit);
hm->AddHistogram(histClass, "mu", "#mu distribution", false, 5000, 0.0, 0.5, VarManager::kMuPerBc);
}
if (subGroupStr.Contains("trigger")) {
hm->AddHistogram(histClass, "IsINT7", "Is INT7", false, 2, -0.5, 1.5, VarManager::kIsINT7);
Expand Down
3 changes: 3 additions & 0 deletions PWGDQ/Core/VarManager.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1000,6 +1000,8 @@ void VarManager::SetDefaultVarNames()
fgVariableUnits[kDCAzNPeaksTrimmed3] = "";
fgVariableNames[kInteractionRate] = "Interaction rate";
fgVariableUnits[kInteractionRate] = "kHz";
fgVariableNames[kMuPerBc] = "#mu";
fgVariableUnits[kMuPerBc] = "";
fgVariableNames[kPt] = "p_{T}";
fgVariableUnits[kPt] = "GeV/c";
fgVariableNames[kPt1] = "p_{T1}";
Expand Down Expand Up @@ -2386,6 +2388,7 @@ void VarManager::SetDefaultVarNames()
fgVarNamesMap["kTwoR2EP2"] = kTwoR2EP2;
fgVarNamesMap["kNPairsPerEvent"] = kNPairsPerEvent;
fgVarNamesMap["kInteractionRate"] = kInteractionRate;
fgVarNamesMap["kMuPerBc"] = kMuPerBc;
fgVarNamesMap["kNEventWiseVariables"] = kNEventWiseVariables;
fgVarNamesMap["kX"] = kX;
fgVarNamesMap["kY"] = kY;
Expand Down
1 change: 1 addition & 0 deletions PWGDQ/Core/VarManager.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -525,6 +525,7 @@ class VarManager : public TObject
kA2EPME_FT0A, // A2 of VM to FT0A event plane for ME technique
kNPairsPerEvent, // number of pairs per event in same-event or mixed-event pairing
kInteractionRate,
kMuPerBc,

// Variables for event mixing with cumulant
kV22m,
Expand Down
68 changes: 0 additions & 68 deletions PWGDQ/TableProducer/tableMaker_withAssoc.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -175,7 +175,6 @@ constexpr static uint32_t gkMFTCovFillMap = VarManager::ObjTypes::TrackMFT | Var
// Enum containing the ordering of statistics histograms to be written in the QA file
enum SkimStatsHists {
kStatsEvent = 0,
kStatsBcs,
kStatsTracks,
kStatsMuons,
kStatsOrphanTracks,
Expand Down Expand Up @@ -258,14 +257,12 @@ struct TableMaker {
// Steer QA output
struct : ConfigurableGroup {
Configurable<bool> fConfigQA{"cfgQA", false, "If true, fill QA histograms"};
Configurable<bool> fConfigFillBcStat{"cfgFillBcStat", false, "If true, fill QA histograms for normalization studies (for OO and Pb-Pb)"};
Configurable<bool> fConfigDetailedQA{"cfgDetailedQA", false, "If true, include more QA histograms (BeforeCuts classes)"};
Configurable<std::string> fConfigAddEventHistogram{"cfgAddEventHistogram", "", "Comma separated list of histograms"};
Configurable<std::string> fConfigAddTrackHistogram{"cfgAddTrackHistogram", "", "Comma separated list of histograms"};
Configurable<std::string> fConfigAddMuonHistogram{"cfgAddMuonHistogram", "", "Comma separated list of histograms"};
Configurable<std::string> fConfigAddEMCalHistogram{"cfgAddEMCalHistogram", "emcal", "Comma separated list of histogram subgroups for the EMCal cluster classes"};
Configurable<std::string> fConfigAddJSONHistograms{"cfgAddJSONHistograms", "", "Histograms in JSON format"};
Configurable<std::string> fConfigIrEstimator{"cfgIrEstimator", "", "Estimator of the interaction rate (pp,OO --> T0VTX, Pb-Pb --> ZNC hadronic), to be used with cfgFillBcStat"};
} fConfigHistOutput;

Configurable<bool> fIsRun2{"cfgIsRun2", false, "Whether we analyze Run-2 or Run-3 data"};
Expand Down Expand Up @@ -728,14 +725,6 @@ struct TableMaker {
histEvents->GetYaxis()->SetBinLabel(o2::aod::evsel::kNsel + 1, "Total");
fStatsList->AddAt(histEvents, kStatsEvent);

std::vector<TString> bcLabels{"all", "tvx", "sel8", "sel8 & Cent", "sel8 & Scent", "sel8 & (Cent | Scent)", "sel8 & (ZNA & ZNC)"};
TH2D* histBcs = new TH2D("BcStats", "Bc statistics;;#mu", bcLabels.size(), -0.5, bcLabels.size() - 0.5, 3000, 0, 0.3);
ib = 1;
for (auto label = bcLabels.begin(); label != bcLabels.end(); label++, ib++) {
histBcs->GetXaxis()->SetBinLabel(ib, (*label).Data());
}
fStatsList->AddAt(histBcs, kStatsBcs);

// Track statistics: one bin for each track selection and 5 bins for V0 tags (gamma, K0s, Lambda, anti-Lambda, Omega)
TH1D* histTracks = new TH1D("TrackStats", "Track statistics", fTrackCuts.size() + 5, -0.5, fTrackCuts.size() - 0.5 + 5.0);
ib = 1;
Expand Down Expand Up @@ -951,22 +940,6 @@ struct TableMaker {
} // end loop over collisions
}

// Function to compute the mu for pileup estimation, taken from EM code
double calculateMu(const auto& bc)
{
uint64_t timeStamp = bc.timestamp();
auto bfilling = mLHCIFdata->getBunchFilling();
double nbc = bfilling.getFilledBCs().size();

const double tvxRate = fConfigHistOutput.fConfigIrEstimator.value.empty()
? mRateFetcher.fetch(fCCDB.service, timeStamp, bc.runNumber(), "T0VTX")
: mRateFetcher.fetch(fCCDB.service, timeStamp, bc.runNumber(), fConfigHistOutput.fConfigIrEstimator.value);

double nTriggersPerFilledBC = tvxRate / nbc / o2::constants::lhc::LHCRevFreq;
double mu = -std::log(1 - nTriggersPerFilledBC);
return mu;
}

template <typename TEvents, typename TTracks>
void computeCollMergingTag(TEvents const& collisions, TTracks const& tracks, Preslice<TTracks>& presliceTracks)
{
Expand Down Expand Up @@ -1059,44 +1032,6 @@ struct TableMaker {
float centFT0A = -1.0;
float centFT0M = -1.0;

if (fConfigHistOutput.fConfigFillBcStat) {
for (const auto& bc : bcs) {
double muTVX = calculateMu(bc);

if (bc.has_ft0()) {
std::bitset<8> fT0Triggers = bc.ft0().triggerMask();
bool isTvx = fT0Triggers[o2::ft0::Triggers::bitVertex];
bool isSemiCentral = fT0Triggers[o2::ft0::Triggers::bitSCen];
bool isCentral = fT0Triggers[o2::ft0::Triggers::bitCen];

bool noBorder = bc.selection_bit(aod::evsel::kNoTimeFrameBorder) && bc.selection_bit(aod::evsel::kNoITSROFrameBorder);
// bool isTriggerTVX = bc.selection_bit(aod::evsel::kIsTriggerTVX); // difference w.r.t fT0Triggers[o2::ft0::Triggers::bitVertex] ?
bool isTriggerZNA = bc.selection_bit(aod::evsel::kIsBBZNA);
bool isTriggerZNC = bc.selection_bit(aod::evsel::kIsBBZNC);

(dynamic_cast<TH2D*>(fStatsList->At(kStatsBcs)))->Fill(0.0, muTVX);
if (isTvx) {
(dynamic_cast<TH2D*>(fStatsList->At(kStatsBcs)))->Fill(1.0, muTVX);
if (noBorder) {
(dynamic_cast<TH2D*>(fStatsList->At(kStatsBcs)))->Fill(2.0, muTVX);
if (isCentral) {
(dynamic_cast<TH2D*>(fStatsList->At(kStatsBcs)))->Fill(3.0, muTVX);
}
if (isSemiCentral) {
(dynamic_cast<TH2D*>(fStatsList->At(kStatsBcs)))->Fill(4.0, muTVX);
}
if (isCentral || isSemiCentral) {
(dynamic_cast<TH2D*>(fStatsList->At(kStatsBcs)))->Fill(5.0, muTVX);
}
if (isTriggerZNA && isTriggerZNC) {
(dynamic_cast<TH2D*>(fStatsList->At(kStatsBcs)))->Fill(6.0, muTVX);
}
}
}
}
}
}

for (const auto& collision : collisions) {

for (int i = 0; i < o2::aod::evsel::kNsel; i++) {
Expand Down Expand Up @@ -1894,9 +1829,6 @@ struct TableMaker {
LOG(fatal) << "Unexpected number of shift values from CCDB: " << fFwdShift->size() << ", expected 1 (z), 3 (x, y, z) or 10 (top/bottom x,y,z + slopes)";
}
}
if (fConfigHistOutput.fConfigFillBcStat) {
mLHCIFdata = fCCDB->getSpecific<o2::parameters::GRPLHCIFData>("GLO/Config/GRPLHCIF", bcs.begin().timestamp());
}
if (fConfigVariousOptions.fPropMuon) {
VarManager::SetupMuonMagField();
}
Expand Down
20 changes: 20 additions & 0 deletions PWGDQ/Tasks/tableReader_withAssoc.h
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -342,6 +342,20 @@ void PrintBitMap(TMap map, int nbits)
LOG(info) << bits;
}

// Function to compute the mu for pileup estimation, taken from EM code
double calculateMu(const double interactionRate, const double nbc)
{
// uint64_t timeStamp = bc.timestamp();
// auto grpmag = fCCDB->getForTimeStamp<o2::parameters::GRPMagField>(grpmagPath, events.begin().timestamp());
// auto grpo = fCCDB->getForTimeStamp<o2::parameters::GRPLHCIFData>(fConfigCCDB.GrpLhcIfPath, timestamp);
// auto bfilling = grpo->getBunchFilling();
// double nbc = bfilling.getFilledBCs().size();

double nTriggersPerFilledBC = (interactionRate * 1000) / nbc / o2::constants::lhc::LHCRevFreq;
double mu = -std::log(1 - nTriggersPerFilledBC);
return mu;
}

// Analysis task that produces event decisions (analysis cut, in bunch pileup and split collision check) and the Hash table used in event mixing
struct AnalysisEventSelection {

Expand Down Expand Up @@ -551,6 +565,12 @@ struct AnalysisEventSelection {
// Get the instantaneous IR from the CCDB
if (fConfigFetchInteractionRate.value) {
VarManager::fgValues[VarManager::kInteractionRate] = rateFetcher.fetch(fCCDB.service, event.timestamp(), fCurrentRun, fConfigIRSource.value, true) / 1000.; // kHz
// Compute Mu
auto grpo = fCCDB->getForTimeStamp<o2::parameters::GRPLHCIFData>("GLO/Config/GRPLHCIF", event.timestamp());
auto bfilling = grpo->getBunchFilling();
double nbc = bfilling.getFilledBCs().size();
// IR is now in kHz, but in calcualteMu it is restored to Hz
VarManager::fgValues[VarManager::kMuPerBc] = calculateMu(VarManager::fgValues[VarManager::kInteractionRate], nbc);
}

bool decision = false;
Expand Down
Loading