Skip to content

[PWGCF] flowMc: add Ev/Tr sel; flowTask/Dihadron: add multCorr config, roll back to wo PID #12245

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
shouqiye merged 4 commits into from
Jul 26, 2025
Merged

[PWGCF] flowMc: add Ev/Tr sel; flowTask/Dihadron: add multCorr config, roll back to wo PID #12245

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
4 commits
Select commit Hold shift + click to select a range
1b15091
add event/track selections
Luzhiyongg Jul 25, 2025
8ea3701
o2 linter modified
Luzhiyongg Jul 25, 2025
f456028
add configurables for multiplicity correlations; dihadron roll back t…
Luzhiyongg Jul 25, 2025
0b3a11e
add centrality weight
Luzhiyongg Jul 25, 2025
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
100 changes: 89 additions & 11 deletions PWGCF/Flow/Tasks/flowMc.cxx
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -54,13 +54,19 @@ struct FlowMc {

Configurable<float> minB{"minB", 0.0f, "min impact parameter"};
Configurable<float> maxB{"maxB", 20.0f, "max impact parameter"};
O2_DEFINE_CONFIGURABLE(cfgCutVertex, float, 10.0f, "Accepted z-vertex range")
O2_DEFINE_CONFIGURABLE(cfgCutPtMin, float, 0.2f, "Minimal pT for tracks")
O2_DEFINE_CONFIGURABLE(cfgCutPtMax, float, 1000.0f, "Maximal pT for tracks")
O2_DEFINE_CONFIGURABLE(cfgCutEta, float, 0.8f, "Eta range for tracks")
O2_DEFINE_CONFIGURABLE(cfgOutputNUAWeights, bool, false, "Fill and output NUA weights")
O2_DEFINE_CONFIGURABLE(cfgCutPtRefMin, float, 0.2f, "Minimal pT for ref tracks")
O2_DEFINE_CONFIGURABLE(cfgCutPtRefMax, float, 3.0f, "Maximal pT for ref tracks")
O2_DEFINE_CONFIGURABLE(cfgCutPtPOIMin, float, 0.2f, "Minimal pT for poi tracks")
O2_DEFINE_CONFIGURABLE(cfgCutPtPOIMax, float, 10.0f, "Maximal pT for poi tracks")
O2_DEFINE_CONFIGURABLE(cfgCutTPCclu, float, 70.0f, "minimum TPC clusters")
O2_DEFINE_CONFIGURABLE(cfgCutITSclu, float, 6.0f, "minimum ITS clusters")
O2_DEFINE_CONFIGURABLE(cfgCutTPCclu, float, 50.0f, "minimum TPC clusters")
O2_DEFINE_CONFIGURABLE(cfgCutITSclu, float, 5.0f, "minimum ITS clusters")
O2_DEFINE_CONFIGURABLE(cfgCutChi2prTPCcls, float, 2.5f, "max chi2 per TPC clusters")
O2_DEFINE_CONFIGURABLE(cfgCutTPCcrossedrows, float, 70.0f, "minimum TPC crossed rows")
O2_DEFINE_CONFIGURABLE(cfgFlowAcceptance, std::string, "", "CCDB path to acceptance object")
O2_DEFINE_CONFIGURABLE(cfgFlowEfficiency, std::string, "", "CCDB path to efficiency object")
O2_DEFINE_CONFIGURABLE(cfgCentVsIPTruth, std::string, "", "CCDB path to centrality vs IP truth")
Expand All @@ -70,6 +76,10 @@ struct FlowMc {
O2_DEFINE_CONFIGURABLE(cfgFlowCumulantNbootstrap, int, 30, "Number of subsamples")
O2_DEFINE_CONFIGURABLE(cfgTrackDensityCorrUse, bool, false, "Use track density efficiency correction")
O2_DEFINE_CONFIGURABLE(cfgTrackDensityCorrSlopeFactor, float, 1.0f, "A factor to scale the track density efficiency slope")
O2_DEFINE_CONFIGURABLE(cfgRecoEvRejectMC, bool, false, "reject both MC and Reco events when reco do not pass")
O2_DEFINE_CONFIGURABLE(cfgRecoEvSel8, bool, false, "require sel8 for reconstruction events")
O2_DEFINE_CONFIGURABLE(cfgRecoEvkIsGoodITSLayersAll, bool, false, "require kIsGoodITSLayersAll for reconstruction events")
O2_DEFINE_CONFIGURABLE(cfgRecoEvkNoSameBunchPileup, bool, false, "require kNoSameBunchPileup for reconstruction events")
Configurable<std::vector<double>> cfgTrackDensityP0{"cfgTrackDensityP0", std::vector<double>{0.6003720411, 0.6152630970, 0.6288860646, 0.6360694031, 0.6409494798, 0.6450540203, 0.6482117301, 0.6512592056, 0.6640008690, 0.6862631416, 0.7005738691, 0.7106567432, 0.7170728333}, "parameter 0 for track density efficiency correction"};
Configurable<std::vector<double>> cfgTrackDensityP1{"cfgTrackDensityP1", std::vector<double>{-1.007592e-05, -8.932635e-06, -9.114538e-06, -1.054818e-05, -1.220212e-05, -1.312304e-05, -1.376433e-05, -1.412813e-05, -1.289562e-05, -1.050065e-05, -8.635725e-06, -7.380821e-06, -6.201250e-06}, "parameter 1 for track density efficiency correction"};
float maxEta = 0.8;
Expand All @@ -81,6 +91,18 @@ struct FlowMc {

ConfigurableAxis axisPt{"axisPt", {VARIABLE_WIDTH, 0.0f, 0.1f, 0.2f, 0.3f, 0.4f, 0.5f, 0.6f, 0.7f, 0.8f, 0.9f, 1.0f, 1.1f, 1.2f, 1.3f, 1.4f, 1.5f, 1.6f, 1.7f, 1.8f, 1.9f, 2.0f, 2.2f, 2.4f, 2.6f, 2.8f, 3.0f, 3.2f, 3.4f, 3.6f, 3.8f, 4.0f, 4.4f, 4.8f, 5.2f, 5.6f, 6.0f, 6.5f, 7.0f, 7.5f, 8.0f, 9.0f, 10.0f, 11.0f, 12.0f}, "pt axis"};

// Filter for MCcollisions
Filter mccollisionFilter = nabs(aod::mccollision::posZ) < cfgCutVertex;
using FilteredMcCollisions = soa::Filtered<aod::McCollisions>;
// Filter for MCParticle
Filter particleFilter = (nabs(aod::mcparticle::eta) < cfgCutEta) && (aod::mcparticle::pt > cfgCutPtMin) && (aod::mcparticle::pt < cfgCutPtMax);
using FilteredMcParticles = soa::Filtered<soa::Join<aod::McParticles, aod::ParticlesToTracks>>;
// Filter for reco tracks
Filter trackFilter = (nabs(aod::track::eta) < cfgCutEta) && (aod::track::pt > cfgCutPtMin) && (aod::track::pt < cfgCutPtMax) && (aod::track::tpcChi2NCl < cfgCutChi2prTPCcls);
using FilteredTracks = soa::Filtered<soa::Join<aod::Tracks, aod::TracksExtra, aod::TrackSelection, aod::TracksDCA, aod::McTrackLabels>>;

// using FilteredTracks = soa::Join<aod::TracksIU, aod::TracksExtra, aod::TrackSelection>;

// Cent vs IP
TH1D* mCentVsIPTruth = nullptr;
bool centVsIPTruthLoaded = false;
Expand All @@ -98,7 +120,6 @@ struct FlowMc {

// Connect to ccdb
Service<ccdb::BasicCCDBManager> ccdb;
Configurable<int64_t> ccdbNoLaterThan{"ccdbNoLaterThan", std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::system_clock::now().time_since_epoch()).count(), "latest acceptable timestamp of creation for the object"};
Configurable<std::string> ccdbUrl{"ccdbUrl", "http://alice-ccdb.cern.ch", "url of the ccdb repository"};

OutputObj<GFWWeights> fWeights{GFWWeights("weights")};
Expand All @@ -110,10 +131,22 @@ struct FlowMc {
std::vector<GFW::CorrConfig> corrconfigsTruth;
std::vector<GFW::CorrConfig> corrconfigsReco;
TRandom3* fRndm = new TRandom3(0);
double epsilon = 1e-6;

void init(InitContext&)
{
ccdb->setURL(ccdbUrl.value);
ccdb->setCaching(true);
auto now = std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::system_clock::now().time_since_epoch()).count();
ccdb->setCreatedNotAfter(now);

const AxisSpec axisVertex{20, -10, 10, "Vtxz (cm)"};
const AxisSpec axisEta{20, -1., 1., "#eta"};
const AxisSpec axisCounter{1, 0, +1, ""};
// QA histograms
histos.add<TH1>("mcEventCounter", "Monte Carlo Truth EventCounter", HistType::kTH1F, {axisCounter});
histos.add<TH1>("numberOfRecoCollisions", "numberOfRecoCollisions", HistType::kTH1F, {{10, -0.5f, 9.5f}});
histos.add<TH1>("RecoEventCounter", "Reconstruction EventCounter", HistType::kTH1F, {axisCounter});
histos.add<TH1>("hnTPCClu", "Number of found TPC clusters", HistType::kTH1D, {{100, 40, 180}});
histos.add<TH1>("hnITSClu", "Number of found ITS clusters", HistType::kTH1D, {{100, 0, 20}});
// pT histograms
Expand Down Expand Up @@ -156,7 +189,9 @@ struct FlowMc {
histos.add<TH2>("hPtNchGeneratedLambda", "Reco production; pT (GeV/c); multiplicity", HistType::kTH2D, {axisPt, axisNch});
histos.add<TH2>("hPtNchGlobalLambda", "Global production; pT (GeV/c); multiplicity", HistType::kTH2D, {axisPt, axisNch});
histos.add<TH1>("hPtMCGen", "Monte Carlo Truth; pT (GeV/c);", {HistType::kTH1D, {axisPt}});
histos.add<TH3>("hEtaPtVtxzMCGen", "Monte Carlo Truth; #eta; p_{T} (GeV/c); V_{z} (cm);", {HistType::kTH3D, {axisEta, axisPt, axisVertex}});
histos.add<TH1>("hPtMCGlobal", "Monte Carlo Global; pT (GeV/c);", {HistType::kTH1D, {axisPt}});
histos.add<TH3>("hEtaPtVtxzMCGlobal", "Monte Carlo Global; #eta; p_{T} (GeV/c); V_{z} (cm);", {HistType::kTH3D, {axisEta, axisPt, axisVertex}});
histos.add<TH1>("hPhiWeightedTrDen", "corrected #phi distribution, considering track density", {HistType::kTH1D, {axisPhi}});

o2::framework::AxisSpec axis = axisPt;
Expand Down Expand Up @@ -318,9 +353,35 @@ struct FlowMc {
}
}

using RecoTracks = soa::Join<aod::TracksIU, aod::TracksExtra, aod::TrackSelection>;
template <typename TCollision>
bool eventSelected(TCollision collision)
{
if (std::fabs(collision.posZ()) > cfgCutVertex) {
return 0;
}
if (cfgRecoEvSel8 && !collision.sel8()) {
return 0;
}
if (cfgRecoEvkNoSameBunchPileup && !collision.selection_bit(o2::aod::evsel::kNoSameBunchPileup)) {
// rejects collisions which are associated with the same "found-by-T0" bunch crossing
// https://indico.cern.ch/event/1396220/#1-event-selection-with-its-rof
return 0;
}
if (cfgRecoEvkIsGoodITSLayersAll && !collision.selection_bit(o2::aod::evsel::kIsGoodITSLayersAll)) {
// from Jan 9 2025 AOT meeting
// cut time intervals with dead ITS staves
return 0;
}
return 1;
}

template <typename TTrack>
bool trackSelected(TTrack track)
{
return ((track.tpcNClsFound() >= cfgCutTPCclu) && (track.tpcNClsCrossedRows() >= cfgCutTPCcrossedrows) && (track.itsNCls() >= cfgCutITSclu));
}

void process(aod::McCollision const& mcCollision, aod::BCsWithTimestamps const&, soa::Join<aod::McParticles, aod::ParticlesToTracks> const& mcParticles, RecoTracks const&)
void process(FilteredMcCollisions::iterator const& mcCollision, aod::BCsWithTimestamps const&, soa::SmallGroups<soa::Join<aod::McCollisionLabels, aod::Collisions, aod::EvSels>> const& collisions, FilteredMcParticles const& mcParticles, FilteredTracks const&)
{

float imp = mcCollision.impactParameter();
Expand All @@ -337,6 +398,21 @@ struct FlowMc {
auto bc = mcCollision.bc_as<aod::BCsWithTimestamps>();
loadCorrections(bc.timestamp());

if (collisions.size() > -1) {
histos.fill(HIST("mcEventCounter"), 0.5);
histos.fill(HIST("numberOfRecoCollisions"), collisions.size()); // number of times coll was reco-ed
if (cfgRecoEvRejectMC) {
if (collisions.size() != 1) { // only pass those have one reconstruction event
return;
}
for (auto const& collision : collisions) {
if (!eventSelected(collision))
return;
}
}
histos.fill(HIST("RecoEventCounter"), 0.5);
}

if (imp > minB && imp < maxB) {
// event within range
histos.fill(HIST("hImpactParameter"), imp);
Expand All @@ -363,9 +439,9 @@ struct FlowMc {
if (std::fabs(mcParticle.eta()) > maxEta) // main acceptance
continue;
if (mcParticle.has_tracks()) {
auto const& tracks = mcParticle.tracks_as<RecoTracks>();
auto const& tracks = mcParticle.tracks_as<FilteredTracks>();
for (auto const& track : tracks) {
if (!((track.tpcNClsFound() >= cfgCutTPCclu) && (track.itsNCls() >= cfgCutITSclu))) {
if (!trackSelected(track)) {
continue;
}
if (cfgIsGlobalTrack && track.isGlobalTrack()) {
Expand Down Expand Up @@ -418,6 +494,7 @@ struct FlowMc {
histos.fill(HIST("hBVsPtVsPhiGenerated"), imp, deltaPhi, mcParticle.pt());
histos.fill(HIST("hPtNchGenerated"), mcParticle.pt(), nChGlobal);
histos.fill(HIST("hPtMCGen"), mcParticle.pt());
histos.fill(HIST("hEtaPtVtxzMCGen"), mcParticle.eta(), mcParticle.pt(), vtxz);
if (pdgCode == PDG_t::kPiPlus)
histos.fill(HIST("hPtNchGeneratedPion"), mcParticle.pt(), nChGlobal);
if (pdgCode == PDG_t::kKPlus)
Expand All @@ -437,9 +514,9 @@ struct FlowMc {
bool validITSTrack = false;
bool validITSABTrack = false;
if (mcParticle.has_tracks()) {
auto const& tracks = mcParticle.tracks_as<RecoTracks>();
auto const& tracks = mcParticle.tracks_as<FilteredTracks>();
for (auto const& track : tracks) {
if (!((track.tpcNClsFound() >= cfgCutTPCclu) && (track.itsNCls() >= cfgCutITSclu))) {
if (!trackSelected(track)) {
continue;
}
histos.fill(HIST("hnTPCClu"), track.tpcNClsFound());
Expand All @@ -456,10 +533,10 @@ struct FlowMc {
if (track.hasTPC()) {
validTPCTrack = true;
}
if (track.hasITS() && track.itsChi2NCl() > -1e-6) {
if (track.hasITS() && track.itsChi2NCl() > -1. * epsilon) {
validITSTrack = true;
}
if (track.hasITS() && track.itsChi2NCl() < -1e-6) {
if (track.hasITS() && track.itsChi2NCl() < -1. * epsilon) {
validITSABTrack = true;
}
}
Expand Down Expand Up @@ -519,6 +596,7 @@ struct FlowMc {
histos.fill(HIST("hBVsPtVsPhiGlobal"), imp, deltaPhi, mcParticle.pt(), wacc * weff);
histos.fill(HIST("hPtNchGlobal"), mcParticle.pt(), nChGlobal);
histos.fill(HIST("hPtMCGlobal"), mcParticle.pt());
histos.fill(HIST("hEtaPtVtxzMCGlobal"), mcParticle.eta(), mcParticle.pt(), vtxz);
if (pdgCode == PDG_t::kPiPlus)
histos.fill(HIST("hPtNchGlobalPion"), mcParticle.pt(), nChGlobal);
if (pdgCode == PDG_t::kKPlus)
Expand Down
Loading
Loading