@@ -54,13 +54,19 @@ struct FlowMc {
5454
5555 Configurable<float> minB{"minB", 0.0f, "min impact parameter"};
5656 Configurable<float> maxB{"maxB", 20.0f, "max impact parameter"};
57+ O2_DEFINE_CONFIGURABLE(cfgCutVertex, float, 10.0f, "Accepted z-vertex range")
58+ O2_DEFINE_CONFIGURABLE(cfgCutPtMin, float, 0.2f, "Minimal pT for tracks")
59+ O2_DEFINE_CONFIGURABLE(cfgCutPtMax, float, 1000.0f, "Maximal pT for tracks")
60+ O2_DEFINE_CONFIGURABLE(cfgCutEta, float, 0.8f, "Eta range for tracks")
5761 O2_DEFINE_CONFIGURABLE(cfgOutputNUAWeights, bool, false, "Fill and output NUA weights")
5862 O2_DEFINE_CONFIGURABLE(cfgCutPtRefMin, float, 0.2f, "Minimal pT for ref tracks")
5963 O2_DEFINE_CONFIGURABLE(cfgCutPtRefMax, float, 3.0f, "Maximal pT for ref tracks")
6064 O2_DEFINE_CONFIGURABLE(cfgCutPtPOIMin, float, 0.2f, "Minimal pT for poi tracks")
6165 O2_DEFINE_CONFIGURABLE(cfgCutPtPOIMax, float, 10.0f, "Maximal pT for poi tracks")
62- O2_DEFINE_CONFIGURABLE(cfgCutTPCclu, float, 70.0f, "minimum TPC clusters")
63- O2_DEFINE_CONFIGURABLE(cfgCutITSclu, float, 6.0f, "minimum ITS clusters")
66+ O2_DEFINE_CONFIGURABLE(cfgCutTPCclu, float, 50.0f, "minimum TPC clusters")
67+ O2_DEFINE_CONFIGURABLE(cfgCutITSclu, float, 5.0f, "minimum ITS clusters")
68+ O2_DEFINE_CONFIGURABLE(cfgCutChi2prTPCcls, float, 2.5f, "max chi2 per TPC clusters")
69+ O2_DEFINE_CONFIGURABLE(cfgCutTPCcrossedrows, float, 70.0f, "minimum TPC crossed rows")
6470 O2_DEFINE_CONFIGURABLE(cfgFlowAcceptance, std::string, "", "CCDB path to acceptance object")
6571 O2_DEFINE_CONFIGURABLE(cfgFlowEfficiency, std::string, "", "CCDB path to efficiency object")
6672 O2_DEFINE_CONFIGURABLE(cfgCentVsIPTruth, std::string, "", "CCDB path to centrality vs IP truth")
@@ -70,6 +76,10 @@ struct FlowMc {
7076 O2_DEFINE_CONFIGURABLE(cfgFlowCumulantNbootstrap, int, 30, "Number of subsamples")
7177 O2_DEFINE_CONFIGURABLE(cfgTrackDensityCorrUse, bool, false, "Use track density efficiency correction")
7278 O2_DEFINE_CONFIGURABLE(cfgTrackDensityCorrSlopeFactor, float, 1.0f, "A factor to scale the track density efficiency slope")
79+ O2_DEFINE_CONFIGURABLE(cfgRecoEvRejectMC, bool, false, "reject both MC and Reco events when reco do not pass")
80+ O2_DEFINE_CONFIGURABLE(cfgRecoEvSel8, bool, false, "require sel8 for reconstruction events")
81+ O2_DEFINE_CONFIGURABLE(cfgRecoEvkIsGoodITSLayersAll, bool, false, "require kIsGoodITSLayersAll for reconstruction events")
82+ O2_DEFINE_CONFIGURABLE(cfgRecoEvkNoSameBunchPileup, bool, false, "require kNoSameBunchPileup for reconstruction events")
7383 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"};
7484 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"};
7585 float maxEta = 0.8;
@@ -81,6 +91,18 @@ struct FlowMc {
8191
8292 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"};
8393
94+ // Filter for MCcollisions
95+ Filter mccollisionFilter = nabs(aod::mccollision::posZ) < cfgCutVertex;
96+ using FilteredMcCollisions = soa::Filtered<aod::McCollisions>;
97+ // Filter for MCParticle
98+ Filter particleFilter = (nabs(aod::mcparticle::eta) < cfgCutEta) && (aod::mcparticle::pt > cfgCutPtMin) && (aod::mcparticle::pt < cfgCutPtMax);
99+ using FilteredMcParticles = soa::Filtered<soa::Join<aod::McParticles, aod::ParticlesToTracks>>;
100+ // Filter for reco tracks
101+ Filter trackFilter = (nabs(aod::track::eta) < cfgCutEta) && (aod::track::pt > cfgCutPtMin) && (aod::track::pt < cfgCutPtMax) && (aod::track::tpcChi2NCl < cfgCutChi2prTPCcls);
102+ using FilteredTracks = soa::Filtered<soa::Join<aod::Tracks, aod::TracksExtra, aod::TrackSelection, aod::TracksDCA, aod::McTrackLabels>>;
103+
104+ // using FilteredTracks = soa::Join<aod::TracksIU, aod::TracksExtra, aod::TrackSelection>;
105+
84106 // Cent vs IP
85107 TH1D* mCentVsIPTruth = nullptr;
86108 bool centVsIPTruthLoaded = false;
@@ -98,7 +120,6 @@ struct FlowMc {
98120
99121 // Connect to ccdb
100122 Service<ccdb::BasicCCDBManager> ccdb;
101- 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"};
102123 Configurable<std::string> ccdbUrl{"ccdbUrl", "http://alice-ccdb.cern.ch", "url of the ccdb repository"};
103124
104125 OutputObj<GFWWeights> fWeights{GFWWeights("weights")};
@@ -110,10 +131,22 @@ struct FlowMc {
110131 std::vector<GFW::CorrConfig> corrconfigsTruth;
111132 std::vector<GFW::CorrConfig> corrconfigsReco;
112133 TRandom3* fRndm = new TRandom3(0);
134+ double epsilon = 1e-6;
113135
114136 void init(InitContext&)
115137 {
138+ ccdb->setURL(ccdbUrl.value);
139+ ccdb->setCaching(true);
140+ auto now = std::chrono::duration_cast<std::chrono::milliseconds>(std::chrono::system_clock::now().time_since_epoch()).count();
141+ ccdb->setCreatedNotAfter(now);
142+
143+ const AxisSpec axisVertex{20, -10, 10, "Vtxz (cm)"};
144+ const AxisSpec axisEta{20, -1., 1., "#eta"};
145+ const AxisSpec axisCounter{1, 0, +1, ""};
116146 // QA histograms
147+ histos.add<TH1>("mcEventCounter", "Monte Carlo Truth EventCounter", HistType::kTH1F, {axisCounter});
148+ histos.add<TH1>("numberOfRecoCollisions", "numberOfRecoCollisions", HistType::kTH1F, {{10, -0.5f, 9.5f}});
149+ histos.add<TH1>("RecoEventCounter", "Reconstruction EventCounter", HistType::kTH1F, {axisCounter});
117150 histos.add<TH1>("hnTPCClu", "Number of found TPC clusters", HistType::kTH1D, {{100, 40, 180}});
118151 histos.add<TH1>("hnITSClu", "Number of found ITS clusters", HistType::kTH1D, {{100, 0, 20}});
119152 // pT histograms
@@ -156,7 +189,9 @@ struct FlowMc {
156189 histos.add<TH2>("hPtNchGeneratedLambda", "Reco production; pT (GeV/c); multiplicity", HistType::kTH2D, {axisPt, axisNch});
157190 histos.add<TH2>("hPtNchGlobalLambda", "Global production; pT (GeV/c); multiplicity", HistType::kTH2D, {axisPt, axisNch});
158191 histos.add<TH1>("hPtMCGen", "Monte Carlo Truth; pT (GeV/c);", {HistType::kTH1D, {axisPt}});
192+ histos.add<TH3>("hEtaPtVtxzMCGen", "Monte Carlo Truth; #eta; p_{T} (GeV/c); V_{z} (cm);", {HistType::kTH3D, {axisEta, axisPt, axisVertex}});
159193 histos.add<TH1>("hPtMCGlobal", "Monte Carlo Global; pT (GeV/c);", {HistType::kTH1D, {axisPt}});
194+ histos.add<TH3>("hEtaPtVtxzMCGlobal", "Monte Carlo Global; #eta; p_{T} (GeV/c); V_{z} (cm);", {HistType::kTH3D, {axisEta, axisPt, axisVertex}});
160195 histos.add<TH1>("hPhiWeightedTrDen", "corrected #phi distribution, considering track density", {HistType::kTH1D, {axisPhi}});
161196
162197 o2::framework::AxisSpec axis = axisPt;
@@ -318,9 +353,35 @@ struct FlowMc {
318353 }
319354 }
320355
321- using RecoTracks = soa::Join<aod::TracksIU, aod::TracksExtra, aod::TrackSelection>;
356+ template <typename TCollision>
357+ bool eventSelected(TCollision collision)
358+ {
359+ if (std::fabs(collision.posZ()) > cfgCutVertex) {
360+ return 0;
361+ }
362+ if (cfgRecoEvSel8 && !collision.sel8()) {
363+ return 0;
364+ }
365+ if (cfgRecoEvkNoSameBunchPileup && !collision.selection_bit(o2::aod::evsel::kNoSameBunchPileup)) {
366+ // rejects collisions which are associated with the same "found-by-T0" bunch crossing
367+ // https://indico.cern.ch/event/1396220/#1-event-selection-with-its-rof
368+ return 0;
369+ }
370+ if (cfgRecoEvkIsGoodITSLayersAll && !collision.selection_bit(o2::aod::evsel::kIsGoodITSLayersAll)) {
371+ // from Jan 9 2025 AOT meeting
372+ // cut time intervals with dead ITS staves
373+ return 0;
374+ }
375+ return 1;
376+ }
377+
378+ template <typename TTrack>
379+ bool trackSelected(TTrack track)
380+ {
381+ return ((track.tpcNClsFound() >= cfgCutTPCclu) && (track.tpcNClsCrossedRows() >= cfgCutTPCcrossedrows) && (track.itsNCls() >= cfgCutITSclu));
382+ }
322383
323- void process(
aod::McCollision const& mcCollision, aod::BCsWithTimestamps const&, soa::Join<aod::
McParticles , aod::
ParticlesToTracks> const& mcParticles,
RecoTracks const&)
384+ 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&)
324385 {
325386
326387 float imp = mcCollision.impactParameter();
@@ -337,6 +398,21 @@ struct FlowMc {
337398 auto bc = mcCollision.bc_as<aod::BCsWithTimestamps>();
338399 loadCorrections(bc.timestamp());
339400
401+ if (collisions.size() > -1) {
402+ histos.fill(HIST("mcEventCounter"), 0.5);
403+ histos.fill(HIST("numberOfRecoCollisions"), collisions.size()); // number of times coll was reco-ed
404+ if (cfgRecoEvRejectMC) {
405+ if (collisions.size() != 1) { // only pass those have one reconstruction event
406+ return;
407+ }
408+ for (auto const& collision : collisions) {
409+ if (!eventSelected(collision))
410+ return;
411+ }
412+ }
413+ histos.fill(HIST("RecoEventCounter"), 0.5);
414+ }
415+
340416 if (imp > minB && imp < maxB) {
341417 // event within range
342418 histos.fill(HIST("hImpactParameter"), imp);
@@ -363,9 +439,9 @@ struct FlowMc {
363439 if (std::fabs(mcParticle.eta()) > maxEta) // main acceptance
364440 continue;
365441 if (mcParticle.has_tracks()) {
366- auto const& tracks = mcParticle.tracks_as<RecoTracks >();
442+ auto const& tracks = mcParticle.tracks_as<FilteredTracks >();
367443 for (auto const& track : tracks) {
368- if (!(( track.tpcNClsFound() >= cfgCutTPCclu) && (track.itsNCls() >= cfgCutITSclu) )) {
444+ if (!trackSelected( track)) {
369445 continue;
370446 }
371447 if (cfgIsGlobalTrack && track.isGlobalTrack()) {
@@ -418,6 +494,7 @@ struct FlowMc {
418494 histos.fill(HIST("hBVsPtVsPhiGenerated"), imp, deltaPhi, mcParticle.pt());
419495 histos.fill(HIST("hPtNchGenerated"), mcParticle.pt(), nChGlobal);
420496 histos.fill(HIST("hPtMCGen"), mcParticle.pt());
497+ histos.fill(HIST("hEtaPtVtxzMCGen"), mcParticle.eta(), mcParticle.pt(), vtxz);
421498 if (pdgCode == PDG_t::kPiPlus)
422499 histos.fill(HIST("hPtNchGeneratedPion"), mcParticle.pt(), nChGlobal);
423500 if (pdgCode == PDG_t::kKPlus)
@@ -437,9 +514,9 @@ struct FlowMc {
437514 bool validITSTrack = false;
438515 bool validITSABTrack = false;
439516 if (mcParticle.has_tracks()) {
440- auto const& tracks = mcParticle.tracks_as<RecoTracks >();
517+ auto const& tracks = mcParticle.tracks_as<FilteredTracks >();
441518 for (auto const& track : tracks) {
442- if (!(( track.tpcNClsFound() >= cfgCutTPCclu) && (track.itsNCls() >= cfgCutITSclu) )) {
519+ if (!trackSelected( track)) {
443520 continue;
444521 }
445522 histos.fill(HIST("hnTPCClu"), track.tpcNClsFound());
@@ -456,10 +533,10 @@ struct FlowMc {
456533 if (track.hasTPC()) {
457534 validTPCTrack = true;
458535 }
459- if (track.hasITS() && track.itsChi2NCl() > -1e-6 ) {
536+ if (track.hasITS() && track.itsChi2NCl() > -1. * epsilon ) {
460537 validITSTrack = true;
461538 }
462- if (track.hasITS() && track.itsChi2NCl() < -1e-6 ) {
539+ if (track.hasITS() && track.itsChi2NCl() < -1. * epsilon ) {
463540 validITSABTrack = true;
464541 }
465542 }
@@ -519,6 +596,7 @@ struct FlowMc {
519596 histos.fill(HIST("hBVsPtVsPhiGlobal"), imp, deltaPhi, mcParticle.pt(), wacc * weff);
520597 histos.fill(HIST("hPtNchGlobal"), mcParticle.pt(), nChGlobal);
521598 histos.fill(HIST("hPtMCGlobal"), mcParticle.pt());
599+ histos.fill(HIST("hEtaPtVtxzMCGlobal"), mcParticle.eta(), mcParticle.pt(), vtxz);
522600 if (pdgCode == PDG_t::kPiPlus)
523601 histos.fill(HIST("hPtNchGlobalPion"), mcParticle.pt(), nChGlobal);
524602 if (pdgCode == PDG_t::kKPlus)
0 commit comments