Skip to content

[PWGCF] Added technical event cuts and QA histograms #18615

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
abilandz merged 1 commit into from
Nov 6, 2026
Merged

[PWGCF] Added technical event cuts and QA histograms #18615

Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
234 changes: 212 additions & 22 deletions PWGCF/MultiparticleCorrelations/Tasks/multiparticleCorrelationsDm.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -15,6 +15,7 @@

// O2:

#include "Common/CCDB/EventSelectionParams.h"
#include "Common/DataModel/Centrality.h"
#include "Common/DataModel/EventSelection.h"
#include "Common/DataModel/Multiplicity.h"
Expand All @@ -39,10 +40,12 @@
#include <TH2.h>
#include <TList.h>
#include <TObject.h>
#include <TProfile.h>
#include <TString.h>
#include <TSystem.h>

#include <array>
#include <cmath>
#include <cstdint>
#include <string>
#include <vector>
Expand Down Expand Up @@ -129,8 +132,17 @@ struct MultiparticleCorrelationsDm { // this name is used in lower-case format t

Configurable<bool> cfDryRun{"cfDryRun", false, "book all histos and run without filling and calculating anything"};
Configurable<bool> cfFillQA{"cfFillQA", false, "book and fill QA histograms (only for development and validation)"}; // switch for QA histograms

// Define cuts:
Configurable<bool> cfUseVertexZCut{"cfUseVertexZCut", true, "apply the vertex Z event cut"};
Configurable<bool> cfUsePtCut{"cfUsePtCut", true, "apply the pt particle cut"};
Configurable<bool> cfUseTimeRange{"cfUseTimeRange", true, "rejects events with overlapping collisions"};
Configurable<bool> cfUseRof{"cfUseRof", true, "removes events containing multiple collisions in the same readout frame"};
Configurable<bool> cfUsePileUp{"cfUsePileUp", true, "apply pile up cuts"};
Configurable<bool> cfUseITSTPC{"cfUseITSTPC", true, "use reconstructed from the ITS (Inner Tracking System) and TPC (Time Porjection Chamber) tracks"};
Configurable<bool> cfUseITSLayers{"cfUseITSLayers", true, "use the ITS layers cut"};
Configurable<bool> cfUseZFT0vsPV{"cfUseZFT0vsPV", true, "use this cut"};
Configurable<bool> cfUsePrevRof{"cfUsePrevRof", true, " use this cut"};

Configurable<std::string> cfCentralityEstimator{"cfCentralityEstimator", "FT0C", "centrality estimator: FT0C, FT0M, FT0A"}; // centrality estimator with valid values: FT0C (default)
Configurable<std::string> cfReferenceMultiplicityEstimator{"cfReferenceMultiplicityEstimator", "FT0C", "reference multiplicity estimator: TPC, FT0C, FT0M"}; // ref. multiplicity estimator with valid values: TPC (default) FT0C Ft0M
Expand All @@ -147,14 +159,19 @@ struct MultiparticleCorrelationsDm { // this name is used in lower-case format t
Configurable<std::vector<float>> cfVertexXBins{"cfVertexXBins", {200, -1., 1.}, "nvxBins, vxMin, vxMax"};
Configurable<std::vector<float>> cfVertexYBins{"cfVertexYBins", {200, -1, 1.}, "nvyBins, vyMin, vyMax"};
Configurable<std::vector<float>> cfVertexZBins{"cfVertexZBins", {300, -15., 15.}, "nvzBins, vzMin, vzMax"};
Configurable<std::vector<float>> cfNContribBins{"cfNContribBins", {500, 0., 5000.}, "nContribBins, contribMin, contribMax"};
Configurable<std::vector<float>> cfVertexZ{"cfVertexZ", {-10., 10.}, "{min, max} [cm]"};
Configurable<std::vector<float>> cfPt{"cfPt", {0.2, 5.0}, "{min, max} [GeV]"};
Configurable<bool> cfUseEtaCut{"cfUseEtaCut", true, "apply the eta particle cut (pseudorapidity"};
Configurable<std::vector<float>> cfEta{"cfEta", {-0.8, 0.8}, "eta range: {min, max}, with convention: min <= eta < max"};

// for weights (efficiency corrections)
Configurable<bool> cfUseWeights{"cfUseWeights", false, "load weights from the file in cfFileWithWeights"};
Configurable<bool> cfUsePhiWeights{"cfUsePhiWeights", false, "use phi weights from the file in cfFileWithWeights"};
Configurable<bool> cfUsePtWeights{"cfUsePtWeights", false, "use pt weights from the file in cfFileWithWeights"};
Configurable<bool> cfUseEtaWeights{"cfUseEtaWeights", false, "use eta weights from the file in cfFileWithWeights"};
Configurable<std::string> cfFileWithWeights{"cfFileWithWeights", "", "path to external ROOT file with weights; on Hyperloop only CCDB paths (/alice-ccdb.cern.ch/...) work"};

Configurable<bool> cfVerbose{"cfVerbose", false, "print \run number, centrality and vertex for each collision"};
Configurable<bool> cfVerbose{"cfVerbose", false, "print run number, centrality and vertex for each collision"};

// *) Define and initialize all data members to be called in the main process* functions:
// **) Task configuration:
Expand All @@ -164,16 +181,27 @@ struct MultiparticleCorrelationsDm { // this name is used in lower-case format t
bool fFillQA = false; // book and fill QA histos
bool fUseVertexZCut = true;
bool fUsePtCut = true;
bool fUseTimeRange = true;
bool fUseRof = true;
bool fUsePileUp = true;
bool fUseITSTPC = true;
bool fUseITSLayers = true;
bool fUseZFT0vsPV = true;
bool fUsePrevRof = true;

int fCentralityEstimator = ECentFT0C;
int fReferenceMultiplicityEstimator = ERefMultTPC;

float fSigmaInel = 771.; // sigma inel [fm^2]
std::array<float, 2> fVertexZ = {-10., 10.}; // {min, max} [cm]
std::array<float, 2> fPt = {0.2, 5.0}; // {min, max} [GeV]
bool fWeightsLoaded = false; // weights are loaded once, at the first collision
bool fUseWeights = false; // load weights from external file
bool fVerbose = false; // print the values / collision
bool fUseEtaCut = true;
std::array<float, 2> fEta = {-0.8, 0.8}; // {min, max}
int fCurrentRunNumber = -1; // run number of the weights currently in use
bool fUsePhiWeights = false; // load phi weights
bool fUsePtWeights = false; // load pt weights
bool fUseEtaWeights = false; // load eta weights
bool fVerbose = false; // print the values / collision

} tc; // you have to prepend "tc." for all objects name in this group later in the code

Expand All @@ -197,11 +225,20 @@ struct MultiparticleCorrelationsDm { // this name is used in lower-case format t

struct QAHistograms {
TList* fQAHistogramsList = nullptr;
TH2F* fHistCentralityRecVsSim = nullptr; // x = centrality sim (from impact parameter), y = centrality rec
} qa; // you have to prepend "qa." for all objects name in this group later in the code
TH2F* fHistCentralityRecVsSim = nullptr; // x = centrality sim (from impact parameter), y = centrality rec
std::array<TH2F*, EBeforeAfter_N> fHistRefMultVsNContrib{}; // x = reference multiplicity, y = number of vertex contributors [before, after]
} qa; // you have to prepend "qa." for all objects name in this group later in the code

// **) Correlations:
struct Correlations {
TList* fCorrelationsList = nullptr;
TProfile* fProfTwo = nullptr; // <<2>> = v2{2}^2 vs. centrality (rec)
} cc;

// **) Histogram with weights from an external file
TH1F* histWeights = nullptr;
TH1F* histPhiWeights = nullptr; // phi weights for the current run
TH1F* histPtWeights = nullptr; // pt weights for the current run
TH1F* histEtaWeights = nullptr; // eta weights for the current run

// *) Define all member functions to be called in the main process* functions:

Expand Down Expand Up @@ -291,8 +328,7 @@ struct MultiparticleCorrelationsDm { // this name is used in lower-case format t
return nullptr;

} // TObject* getObjectFromList(...)

TH1F* getHistogramWithWeights(const char* filePath, const char* runNumber)
TH1F* getHistogramWithWeights(const char* filePath, const char* runNumber, const char* histName)
{
// Last update: 20261018

Expand Down Expand Up @@ -427,7 +463,7 @@ struct MultiparticleCorrelationsDm { // this name is used in lower-case format t
// listWithRuns->ls();

// get pointer to "hist" from "listWithRuns"
hist = dynamic_cast<TH1F*>(getObjectFromList(listWithRuns, "histWithEfficiencyCorrections_1"));
hist = dynamic_cast<TH1F*>(getObjectFromList(listWithRuns, histName));
if (!hist) {
LOGF(fatal, "\033[1;31m%s at line %d\033[0m", __FUNCTION__, __LINE__);
}
Expand All @@ -443,6 +479,36 @@ struct MultiparticleCorrelationsDm { // this name is used in lower-case format t

} // end of TH1F* getHistogramWithWeights(...)

// load one kind of weights for this run (or reuse it if already loaded), and save it in the output
TH1F* loadWeights(const char* nameInFile, const char* nameInOutput, int runNumber)
{
TString name = Form("%s_%d", nameInOutput, runNumber);
TH1F* hist = dynamic_cast<TH1F*>(fBaseList->FindObject(name.Data()));
if (!hist) {
hist = getHistogramWithWeights(cfFileWithWeights.value.c_str(), Form("%d", runNumber), nameInFile);
hist->SetName(name.Data()); // e.g. histPhiWeights_544116
hist->SetTitle(cfFileWithWeights.value.c_str()); // book-keeping: title = file path
fBaseList->Add(hist); // saved in AnalysisResults.root
}
return hist;
}

// weight of one particle: product of all enabled weights, 1 if none is enabled
float getParticleWeight(float phi, float pt, float eta)
{
float weight = 1.;
if (tc.fUsePhiWeights && histPhiWeights) {
weight *= histPhiWeights->GetBinContent(histPhiWeights->FindBin(phi));
}
if (tc.fUsePtWeights && histPtWeights) {
weight *= histPtWeights->GetBinContent(histPtWeights->FindBin(pt));
}
if (tc.fUseEtaWeights && histEtaWeights) {
weight *= histEtaWeights->GetBinContent(histEtaWeights->FindBin(eta));
}
return weight;
}

//__ EVENT CUTS __
template <ERecSim rs, typename T1>
bool eventCuts(T1 const& collision)
Expand All @@ -454,6 +520,55 @@ struct MultiparticleCorrelationsDm { // this name is used in lower-case format t
return false; // if at least one event cut is not satisfied
}
}
// *) No collision in time range cut
if (tc.fUseTimeRange) {
if (!collision.selection_bit(o2::aod::evsel::kNoCollInTimeRangeStandard)) {
return false; // if at least one event cut is not satisfied
}
}

// *) No collision in the same readout frame cut
if (tc.fUseRof) {
if (!collision.selection_bit(o2::aod::evsel::kNoCollInRofStandard)) {
return false; // if at least one event cut is not satisfied
}
}

// *) No same-bunch pile-up cut
if (tc.fUsePileUp) {
if (!collision.selection_bit(o2::aod::evsel::kNoSameBunchPileup)) {
return false; // if at least one event cut is not satisfied
}
}

// *) Vertex from ITS and TPC tracks cut
if (tc.fUseITSTPC) {
if (!collision.selection_bit(o2::aod::evsel::kIsVertexITSTPC)) {
return false; // if at least one event cut is not satisfied
}
}

// *) All ITS layers good cut
if (tc.fUseITSLayers) {
if (!collision.selection_bit(o2::aod::evsel::kIsGoodITSLayersAll)) {
return false; // if at least one event cut is not satisfied
}
}

// *) z vertex from FT0 and from tracks agree cut
if (tc.fUseZFT0vsPV) {
if (!collision.selection_bit(o2::aod::evsel::kIsGoodZvtxFT0vsPV)) {
return false; // if at least one event cut is not satisfied
}
}

// *) No high-multiplicity collision in previous readout frame cut
if (tc.fUsePrevRof) {
if (!collision.selection_bit(o2::aod::evsel::kNoHighMultCollInPrevRof)) {
return false; // if at least one event cut is not satisfied
}
}

return true; // all event cuts passed
}

Expand All @@ -469,6 +584,13 @@ struct MultiparticleCorrelationsDm { // this name is used in lower-case format t
}
}

// *) eta cut
if (tc.fUseEtaCut) {
if (!(tc.fEta[0] <= track.eta() && track.eta() < tc.fEta[1])) {
return false; // if at least one particle cut is not satisfied
}
}

return true; // all particle cuts passed
}

Expand All @@ -488,6 +610,11 @@ struct MultiparticleCorrelationsDm { // this name is used in lower-case format t
ec.fEventHistograms[EVertexY][ERec][ba]->Fill(collision.posY());
ec.fEventHistograms[EVertexZ][ERec][ba]->Fill(collision.posZ());

// QA: reference multiplicity vs. number of contributors
if (tc.fFillQA) {
qa.fHistRefMultVsNContrib[ba]->Fill(getReferenceMultiplicity(collision), collision.numContrib());
}

if constexpr (rs == ERecAndSim) {
if (!collision.has_mcCollision()) {
LOGF(debug, " No MC collision for this collision, skip simulated part...");
Expand Down Expand Up @@ -529,14 +656,19 @@ struct MultiparticleCorrelationsDm { // this name is used in lower-case format t
return;
}

// Load the weights once, at the first collision, because the run number is only known here
if (tc.fUseWeights && !tc.fWeightsLoaded) {
histWeights = getHistogramWithWeights(cfFileWithWeights.value.c_str(), Form("%d", collision.bc().runNumber()));
histWeights->SetTitle(cfFileWithWeights.value.c_str()); // book-keeping: title = file path
// value takes the text as std::string and c_str makes it into a constant char*

fBaseList->Add(histWeights); // so it is saved in AnalysisResults.root
tc.fWeightsLoaded = true;
// weights are run-specific: load them again when the run number changes
int runNumber = collision.bc().runNumber();
if (runNumber != tc.fCurrentRunNumber) {
if (tc.fUsePhiWeights) {
histPhiWeights = loadWeights("histWithEfficiencyCorrections_1", "histPhiWeights", runNumber);
}
if (tc.fUsePtWeights) {
histPtWeights = loadWeights("histWithEfficiencyCorrections_2", "histPtWeights", runNumber);
}
if (tc.fUseEtaWeights) {
histEtaWeights = loadWeights("histWithEfficiencyCorrections_3", "histEtaWeights", runNumber);
}
tc.fCurrentRunNumber = runNumber;
}

// Print info for each collision, only if switched on in the config:
Expand All @@ -555,6 +687,9 @@ struct MultiparticleCorrelationsDm { // this name is used in lower-case format t
}

// Main loop over particles:
// Q-vector of harmonic 2 and sums of weights, for <2> of this event
double q2Re = 0., q2Im = 0., sumW = 0., sumW2 = 0.;

auto track = tracks.iteratorAt(0); // set the type and scope from one instance
for (int64_t i = 0; i < tracks.size(); i++) {

Expand Down Expand Up @@ -594,6 +729,13 @@ struct MultiparticleCorrelationsDm { // this name is used in lower-case format t
pc.fParticleHistograms[EPt][ERec][EAfter]->Fill(track.pt());
pc.fParticleHistograms[EPhi][ERec][EAfter]->Fill(track.phi());

// add this particle to the Q-vector
const double w = getParticleWeight(track.phi(), track.pt(), track.eta());
q2Re += w * std::cos(2. * track.phi());
q2Im += w * std::sin(2. * track.phi());
sumW += w;
sumW2 += w * w;

if constexpr (rs == ERecAndSim) {
if (!track.has_mcParticle()) {
continue;
Expand All @@ -610,6 +752,15 @@ struct MultiparticleCorrelationsDm { // this name is used in lower-case format t
} // end of for (int64_t i = 0; i < tracks.size(); i++) {

// Fill event histograms after event and particle cuts
// <2> of this event, filled vs. centrality with event weight sumW^2 - sumW2 (= M(M-1) for unit weights)
if constexpr (rs == ERec || rs == ERecAndSim) {
const double denominator = sumW * sumW - sumW2;
if (denominator > 0.) {
const double two = (q2Re * q2Re + q2Im * q2Im - sumW2) / denominator;
cc.fProfTwo->Fill(getCentrality(collision), two, denominator);
}
}

fillEventHistograms<rs>(collision, tracks, EAfter);

} // end of template <ERecSim rs, typename T1, typename T2> void steer(T1 const& collision, T2 const& tracks) {
Expand All @@ -631,9 +782,11 @@ struct MultiparticleCorrelationsDm { // this name is used in lower-case format t
tc.fVerbose = cfVerbose;
tc.fSigmaInel = cfSigmaInel;

tc.fUseWeights = cfUseWeights;
if (tc.fUseWeights && cfFileWithWeights.value.empty()) {
LOGF(fatal, "cfUseWeights is true, but cfFileWithWeights is empty");
tc.fUsePhiWeights = cfUsePhiWeights;
tc.fUsePtWeights = cfUsePtWeights;
tc.fUseEtaWeights = cfUseEtaWeights;
if ((tc.fUsePhiWeights || tc.fUsePtWeights || tc.fUseEtaWeights) && cfFileWithWeights.value.empty()) {
LOGF(fatal, "weights are switched on, but cfFileWithWeights is empty");
}

// vertex Z cut copy the swtich and range
Expand All @@ -650,6 +803,22 @@ struct MultiparticleCorrelationsDm { // this name is used in lower-case format t
tc.fPt[0] = lPtCut[0];
tc.fPt[1] = lPtCut[1];

// eta cut copy switch and range
tc.fUseEtaCut = cfUseEtaCut;
std::vector<float> lEtaCut = cfEta.value;

tc.fEta[0] = lEtaCut[0];
tc.fEta[1] = lEtaCut[1];

// technical event cuts: copy the switches
tc.fUseTimeRange = cfUseTimeRange;
tc.fUseRof = cfUseRof;
tc.fUsePileUp = cfUsePileUp;
tc.fUseITSTPC = cfUseITSTPC;
tc.fUseITSLayers = cfUseITSLayers;
tc.fUseZFT0vsPV = cfUseZFT0vsPV;
tc.fUsePrevRof = cfUsePrevRof;

// translate the centrality estimator string into the enum
if (cfCentralityEstimator.value == "FT0C") {
tc.fCentralityEstimator = ECentFT0C;
Expand Down Expand Up @@ -688,6 +857,18 @@ struct MultiparticleCorrelationsDm { // this name is used in lower-case format t
ec.fEventHistogramsList->SetOwner(true);
fBaseList->Add(ec.fEventHistogramsList); // any nested TList in the base TList appears as a subdir in the output ROOT file

// correlations: <<2>> vs. centrality, in the centrality classes of the reference paper
cc.fCorrelationsList = new TList();
cc.fCorrelationsList->SetName("Correlations");
cc.fCorrelationsList->SetOwner(true);
fBaseList->Add(cc.fCorrelationsList);
const std::array<double, 10> centBoundaries = {0., 5., 10., 20., 30., 40., 50., 60., 70., 80.};
cc.fProfTwo = new TProfile("fProfTwo", "<<2>> = v_{2}{2}^{2} vs. centrality", centBoundaries.size() - 1, centBoundaries.data());
cc.fProfTwo->GetXaxis()->SetTitle("centrality (%)");
cc.fProfTwo->GetYaxis()->SetTitle("#LT#LT2#GT#GT");
cc.fProfTwo->Sumw2();
cc.fCorrelationsList->Add(cc.fProfTwo);

if (tc.fFillQA) { // only if QA switched on
qa.fQAHistogramsList = new TList();
qa.fQAHistogramsList->SetName("QA");
Expand Down Expand Up @@ -824,6 +1005,15 @@ struct MultiparticleCorrelationsDm { // this name is used in lower-case format t
qa.fHistCentralityRecVsSim->GetXaxis()->SetTitle("centrality simulated (%)");
qa.fHistCentralityRecVsSim->GetYaxis()->SetTitle("centrality reconstructed (%)");
qa.fQAHistogramsList->Add(qa.fHistCentralityRecVsSim);

// reference multiplicity vs. number of contributors, before and after cuts
std::vector<float> lContribBins = cfNContribBins.value;
for (int ba = 0; ba < EBeforeAfter_N; ba++) {
qa.fHistRefMultVsNContrib[ba] = new TH2F(Form("fHistRefMultVsNContrib[%s]", baNames[ba]), Form("reference multiplicity vs. number of contributors (%s cuts)", baNames[ba]), nRefMultBins, refMultMin, refMultMax, static_cast<int>(lContribBins[0]), lContribBins[1], lContribBins[2]);
qa.fHistRefMultVsNContrib[ba]->GetXaxis()->SetTitle("reference multiplicity");
qa.fHistRefMultVsNContrib[ba]->GetYaxis()->SetTitle("number of contributors");
qa.fQAHistogramsList->Add(qa.fHistRefMultVsNContrib[ba]);
}
}

} // end of void init(InitContext&) {
Expand Down
Loading