Skip to content

Commit 212a6e1

Browse files
wenyaCernalibuild
wenyaCern
and
alibuild
authored
[PWGCF] Add MC gen process (AliceO2Group#16765)
Co-authored-by: ALICE Action Bot <alibuild@cern.ch>
1 parent bbf97ee commit 212a6e1

1 file changed

Lines changed: 142 additions & 33 deletions

File tree

‎PWGCF/Flow/Tasks/flowFlucGfwPp.cxx‎

Lines changed: 142 additions & 33 deletions
Original file line numberDiff line numberDiff line change
@@ -121,7 +121,7 @@ struct FlowFlucGfwPp {
121121

122122
O2_DEFINE_CONFIGURABLE(cfgNbootstrap, int, 10, "Number of subsamples")
123123
O2_DEFINE_CONFIGURABLE(cfgIsMC, bool, false, "Is MC event")
124-
O2_DEFINE_CONFIGURABLE(cfgCentEstimator, int, 0, "0:FT0C; 1:FT0CVariant1; 2:FT0M; 3:FV0A, 4:NTPV, 5:NGlobals, 6:MFT")
124+
O2_DEFINE_CONFIGURABLE(cfgCentEstimator, int, 0, "0:FT0C; 1:FT0CVariant1; 2:FT0M; 3:FV0A, 4:NTPV, 5:NGlobals")
125125
O2_DEFINE_CONFIGURABLE(cfgUseNch, bool, false, "Do correlations as function of Nch")
126126
O2_DEFINE_CONFIGURABLE(cfgQvecQA, bool, false, "Enable filling QA for q-Vec of TPC")
127127
O2_DEFINE_CONFIGURABLE(cfgFillWeights, bool, false, "Fill NUA weights")
@@ -171,6 +171,7 @@ struct FlowFlucGfwPp {
171171
O2_DEFINE_CONFIGURABLE(cfgUseNegativeEtaHalfForq2, bool, true, "If true, use -eta half for qn selection; otherwise use +eta half");
172172
O2_DEFINE_CONFIGURABLE(cfgQnSelectionHarmonic, int, 2, "Harmonic n used to build the reduced q_n vector for event shape selection, use 2 for q2 and 3 for q3");
173173
O2_DEFINE_CONFIGURABLE(cfgQnHistMax, float, 6., "Upper range for q_n calibration histograms");
174+
O2_DEFINE_CONFIGURABLE(cfgQnTrkAbsEtaMax, float, 0.5, "Upper range for abs eta of tracks contributing to q_n");
174175
O2_DEFINE_CONFIGURABLE(cfgBypassQnSelection, bool, false, "Bypass q_n event shape selection and fill one integral q-bin");
175176
O2_DEFINE_CONFIGURABLE(cfgMinPtOnTPC, float, 0.2, "minimum transverse momentum selection for TPC tracks participating in Q-vector reconstruction");
176177
O2_DEFINE_CONFIGURABLE(cfgMaxPtOnTPC, float, 5., "maximum transverse momentum selection for TPC tracks participating in Q-vector reconstruction");
@@ -277,8 +278,7 @@ struct FlowFlucGfwPp {
277278
kCentFT0M,
278279
kCentFV0A,
279280
kCentNTPV,
280-
kCentNGlobal,
281-
kCentMFT
281+
kCentNGlobal
282282
};
283283
enum EventSelFlags {
284284
kFilteredEvent = 1,
@@ -337,11 +337,11 @@ struct FlowFlucGfwPp {
337337
o2::framework::expressions::Filter collisionFilter = nabs(aod::collision::posZ) < cfgVtxZ;
338338
o2::framework::expressions::Filter trackFilter = nabs(aod::track::eta) < cfgEta && aod::track::pt > cfgPtmin&& aod::track::pt < cfgPtmax && (aod::track::itsChi2NCl < cfgChi2PrITSCls) && (aod::track::tpcChi2NCl < cfgChi2PrTPCCls) && nabs(aod::track::dcaZ) < cfgDCAz;
339339

340-
Preslice<aod::Tracks> perCollision = aod::track::collisionId;
341340
o2::framework::expressions::Filter mcCollFilter = nabs(aod::mccollision::posZ) < cfgVtxZ;
342341
o2::framework::expressions::Filter mcParticlesFilter = (aod::mcparticle::eta > o2::analysis::gfwflowflucpp::etalow && aod::mcparticle::eta < o2::analysis::gfwflowflucpp::etaup && aod::mcparticle::pt > o2::analysis::gfwflowflucpp::ptlow && aod::mcparticle::pt < o2::analysis::gfwflowflucpp::ptup);
343342

344343
using GFWTracks = soa::Filtered<soa::Join<aod::Tracks, aod::TracksExtra, aod::TrackSelection, aod::TrackSelectionExtension, aod::TracksDCA>>;
344+
using GFWTracksMC = soa::Filtered<soa::Join<aod::Tracks, aod::TracksExtra, aod::TrackSelection, aod::TrackSelectionExtension, aod::TracksDCA, aod::McTrackLabels>>;
345345

346346
void init(InitContext const&)
347347
{
@@ -415,8 +415,7 @@ struct FlowFlucGfwPp {
415415
{kCentFT0M, "FT0M"},
416416
{kCentFV0A, "FV0A"},
417417
{kCentNTPV, "NTPV"},
418-
{kCentNGlobal, "NGlobals"},
419-
{kCentMFT, "MFT"}};
418+
{kCentNGlobal, "NGlobals"}};
420419
sCentralityEstimator = centEstimatorMap.at(cfgCentEstimator);
421420
sCentralityEstimator += " centrality (%)";
422421
AxisSpec centAxis = {o2::analysis::gfwflowflucpp::centbinning, sCentralityEstimator.c_str()};
@@ -503,7 +502,7 @@ struct FlowFlucGfwPp {
503502
}
504503
}
505504

506-
if (doprocessData) {
505+
if (doprocessData || doprocessMC) {
507506
registry.add("trackQA/before/phi_eta_vtxZ", "", {HistType::kTH3D, {phiAxis, etaAxis, vtxAxis}});
508507
registry.add("trackQA/before/pt_dcaXY_dcaZ", "", {HistType::kTH3D, {ptAxis, dcaXYAXis, dcaZAXis}});
509508
registry.add("trackQA/before/nch_pt", "#it{p}_{T} vs multiplicity; N_{ch}; #it{p}_{T}", {HistType::kTH2D, {nchAxis, ptAxis}});
@@ -530,24 +529,20 @@ struct FlowFlucGfwPp {
530529
registry.add("eventQA/before/globalTracks_multV0A", "", {HistType::kTH2D, {v0aAxis, nchAxis}});
531530
registry.add("eventQA/before/multV0A_multT0A", "", {HistType::kTH2D, {t0aAxis, v0aAxis}});
532531

533-
if (doprocessData) {
534-
registry.add("eventQA/before/centrality", "", {HistType::kTH1D, {centAxis}});
535-
registry.add("eventQA/before/globalTracks_centT0C", "", {HistType::kTH2D, {centAxis, nchAxis}});
536-
registry.add("eventQA/before/PVTracks_centT0C", "", {HistType::kTH2D, {centAxis, multpvAxis}});
537-
registry.add("eventQA/before/multT0C_centT0C", "", {HistType::kTH2D, {centAxis, t0cAxis}});
538-
539-
registry.add("eventQA/before/centT0M_centT0C", "", {HistType::kTH2D, {centAxis, centAxis}});
540-
registry.add("eventQA/before/centV0A_centT0C", "", {HistType::kTH2D, {centAxis, centAxis}});
541-
registry.add("eventQA/before/centGlobal_centT0C", "", {HistType::kTH2D, {centAxis, centAxis}});
542-
registry.add("eventQA/before/centNTPV_centT0C", "", {HistType::kTH2D, {centAxis, centAxis}});
543-
registry.add("eventQA/before/centMFT_centT0C", "", {HistType::kTH2D, {centAxis, centAxis}});
544-
545-
if (cfgIsMC) {
546-
registry.add("MCGen/trackQA/phi_eta_vtxZ", "", {HistType::kTH3D, {phiAxis, etaAxis, vtxAxis}});
547-
registry.add("MCGen/trackQA/nch_pt", "#it{p}_{T} vs multiplicity; N_{ch}; #it{p}_{T}", {HistType::kTH2D, {nchAxis, ptAxis}});
548-
registry.add("MCGen/trackQA/pt_ref", "", {HistType::kTH1D, {{100, o2::analysis::gfwflowflucpp::ptreflow, o2::analysis::gfwflowflucpp::ptrefup}}});
549-
registry.add("MCGen/trackQA/pt_poi", "", {HistType::kTH1D, {{100, o2::analysis::gfwflowflucpp::ptpoilow, o2::analysis::gfwflowflucpp::ptpoiup}}});
550-
}
532+
registry.add("eventQA/before/centrality", "", {HistType::kTH1D, {centAxis}});
533+
registry.add("eventQA/before/globalTracks_centT0C", "", {HistType::kTH2D, {centAxis, nchAxis}});
534+
registry.add("eventQA/before/PVTracks_centT0C", "", {HistType::kTH2D, {centAxis, multpvAxis}});
535+
registry.add("eventQA/before/multT0C_centT0C", "", {HistType::kTH2D, {centAxis, t0cAxis}});
536+
537+
registry.add("eventQA/before/centT0M_centT0C", "", {HistType::kTH2D, {centAxis, centAxis}});
538+
registry.add("eventQA/before/centV0A_centT0C", "", {HistType::kTH2D, {centAxis, centAxis}});
539+
registry.add("eventQA/before/centGlobal_centT0C", "", {HistType::kTH2D, {centAxis, centAxis}});
540+
registry.add("eventQA/before/centNTPV_centT0C", "", {HistType::kTH2D, {centAxis, centAxis}});
541+
if (cfgIsMC || doprocessMC) {
542+
registry.add("MCGen/trackQA/phi_eta_vtxZ", "", {HistType::kTH3D, {phiAxis, etaAxis, vtxAxis}});
543+
registry.add("MCGen/trackQA/nch_pt", "#it{p}_{T} vs multiplicity; N_{ch}; #it{p}_{T}", {HistType::kTH2D, {nchAxis, ptAxis}});
544+
registry.add("MCGen/trackQA/pt_ref", "", {HistType::kTH1D, {{100, o2::analysis::gfwflowflucpp::ptreflow, o2::analysis::gfwflowflucpp::ptrefup}}});
545+
registry.add("MCGen/trackQA/pt_poi", "", {HistType::kTH1D, {{100, o2::analysis::gfwflowflucpp::ptpoilow, o2::analysis::gfwflowflucpp::ptpoiup}}});
551546
}
552547

553548
registry.addClone("eventQA/before/", "eventQA/after/");
@@ -1020,14 +1015,14 @@ struct FlowFlucGfwPp {
10201015
registry.fill(HIST("qvecQA/ChTracks"), trk.pt(), trk.eta(), trk.phi());
10211016
}
10221017

1023-
if (trk.eta() > 0) {
1018+
if (trk.eta() > 0 && std::fabs(trk.eta()) < cfgQnTrkAbsEtaMax) {
10241019
// In qVectorsTable this branch is additionally guarded by useDetector["QvectorTPCposs"] || useDetector["QvectorBPoss"].
10251020
// Here TPCpos is always computed because the downstream ESE selector can require it.
10261021
qvec.qVectTPCPos[0] += trk.pt() * std::cos(trk.phi() * nMode);
10271022
qvec.qVectTPCPos[1] += trk.pt() * std::sin(trk.phi() * nMode);
10281023
qvec.trkTPCPosLabel.push_back(trk.globalIndex());
10291024
qvec.nTrkTPCPos++;
1030-
} else if (trk.eta() < 0) {
1025+
} else if (trk.eta() < 0 && std::fabs(trk.eta()) < cfgQnTrkAbsEtaMax) {
10311026
// In qVectorsTable this branch is additionally guarded by useDetector["QvectorTPCnegs"] || useDetector["QvectorBNegs"].
10321027
// Here TPCneg is always computed because the downstream ESE selector can require it.
10331028
qvec.qVectTPCNeg[0] += trk.pt() * std::cos(trk.phi() * nMode);
@@ -1172,6 +1167,48 @@ struct FlowFlucGfwPp {
11721167
lRandom, qPtmp, run);
11731168
}
11741169

1170+
template <typename TCollision, typename TParticles>
1171+
void processGenCollision(TCollision collision, TParticles particles, const int& mcCollisionId, const XAxis& xaxis, const int& run, const int& qPtmp)
1172+
{
1173+
if (xaxis.multiplicity < cfgFixedMultMin || xaxis.multiplicity > cfgFixedMultMax)
1174+
return;
1175+
1176+
if (cfgFillQA && xaxis.centrality >= 0)
1177+
registry.fill(HIST("eventQA/after/centrality"), xaxis.centrality);
1178+
if (cfgFillQA)
1179+
registry.fill(HIST("eventQA/after/multiplicity"), xaxis.multiplicity);
1180+
1181+
fGFW->Clear();
1182+
float lRandom = fRndm->Rndm();
1183+
float vtxz = collision.posZ();
1184+
1185+
AcceptedTracks acceptedTracks{0, 0, 0, 0};
1186+
for (const auto& particle : particles) {
1187+
if (particle.mcCollisionId() != mcCollisionId)
1188+
continue;
1189+
processTrack(particle, vtxz, xaxis.multiplicity, run, acceptedTracks);
1190+
}
1191+
1192+
if (cfgConsistentEventFlag & kRequireBothEtaSides)
1193+
if (!acceptedTracks.nPos || !acceptedTracks.nNeg)
1194+
return;
1195+
if (cfgConsistentEventFlag & kRequireFullFourParticleTracks)
1196+
if (acceptedTracks.nFull < kMinTracksForFourParticleCorrelation)
1197+
return;
1198+
if (cfgConsistentEventFlag & kRequireTwoTracksInBothEtaSides)
1199+
if (acceptedTracks.nPos < kMinTracksPerEtaSideForGapCorrelation ||
1200+
acceptedTracks.nNeg < kMinTracksPerEtaSideForGapCorrelation)
1201+
return;
1202+
if (cfgConsistentEventFlag & kRequireTwoTracksInThreeEtaRegions)
1203+
if (acceptedTracks.nPos < kMinTracksPerEtaRegionForThreeSubevents ||
1204+
acceptedTracks.nMid < kMinTracksPerEtaRegionForThreeSubevents ||
1205+
acceptedTracks.nNeg < kMinTracksPerEtaRegionForThreeSubevents)
1206+
return;
1207+
1208+
fillOutputContainers<kGen>(cfgUseNch ? static_cast<float>(xaxis.multiplicity) : xaxis.centrality,
1209+
lRandom, qPtmp, run);
1210+
}
1211+
11751212
bool isStable(int pdg)
11761213
{
11771214
if (std::abs(pdg) == PDG_t::kPiPlus)
@@ -1334,8 +1371,6 @@ struct FlowFlucGfwPp {
13341371
return collision.centNTPV();
13351372
case kCentNGlobal:
13361373
return collision.centNGlobal();
1337-
case kCentMFT:
1338-
return collision.centMFT();
13391374
default:
13401375
return collision.centFT0C();
13411376
}
@@ -1352,7 +1387,6 @@ struct FlowFlucGfwPp {
13521387
registry.fill(HIST("eventQA/") + HIST(FillTimeName[ft]) + HIST("centV0A_centT0C"), collision.centFT0C(), collision.centFV0A());
13531388
registry.fill(HIST("eventQA/") + HIST(FillTimeName[ft]) + HIST("centGlobal_centT0C"), collision.centFT0C(), collision.centNGlobal());
13541389
registry.fill(HIST("eventQA/") + HIST(FillTimeName[ft]) + HIST("centNTPV_centT0C"), collision.centFT0C(), collision.centNTPV());
1355-
registry.fill(HIST("eventQA/") + HIST(FillTimeName[ft]) + HIST("centMFT_centT0C"), collision.centFT0C(), collision.centMFT());
13561390
}
13571391
registry.fill(HIST("eventQA/") + HIST(FillTimeName[ft]) + HIST("globalTracks_PVTracks"), collision.multNTracksPV(), xaxis.multiplicity);
13581392
registry.fill(HIST("eventQA/") + HIST(FillTimeName[ft]) + HIST("globalTracks_multT0A"), collision.multFT0A(), xaxis.multiplicity);
@@ -1406,8 +1440,7 @@ struct FlowFlucGfwPp {
14061440

14071441
void processData(soa::Filtered<soa::Join<aod::Collisions, aod::EvSels, aod::Mults,
14081442
aod::CentFT0Cs, aod::CentFT0CVariant1s, aod::CentFT0Ms,
1409-
aod::CentFV0As, aod::CentNTPVs, aod::CentNGlobals,
1410-
aod::CentMFTs>>::iterator const& collision,
1443+
aod::CentFV0As, aod::CentNTPVs, aod::CentNGlobals>>::iterator const& collision,
14111444
aod::BCsWithTimestamps const&, GFWTracks const& tracks)
14121445
{
14131446
auto bc = collision.bc_as<aod::BCsWithTimestamps>();
@@ -1487,7 +1520,7 @@ struct FlowFlucGfwPp {
14871520
}
14881521
PROCESS_SWITCH(FlowFlucGfwPp, processData, "Process analysis for non-derived data", false);
14891522

1490-
void processq2(soa::Filtered<soa::Join<aod::Collisions, aod::EvSels, aod::Mults, aod::CentFT0Cs, aod::CentFT0CVariant1s, aod::CentFT0Ms, aod::CentFV0As, aod::CentNTPVs, aod::CentNGlobals, aod::CentMFTs>>::iterator const& collision, aod::BCsWithTimestamps const&, GFWTracks const& tracks)
1523+
void processq2(soa::Filtered<soa::Join<aod::Collisions, aod::EvSels, aod::Mults, aod::CentFT0Cs, aod::CentFT0CVariant1s, aod::CentFT0Ms, aod::CentFV0As, aod::CentNTPVs, aod::CentNGlobals>>::iterator const& collision, aod::BCsWithTimestamps const&, GFWTracks const& tracks)
14911524
{
14921525
float count{0.5};
14931526
fillQnEventCounter(count++);
@@ -1521,6 +1554,82 @@ struct FlowFlucGfwPp {
15211554
fillQnCalibrationHistograms(centr, multi, qvecPos, qvecNeg);
15221555
}
15231556
PROCESS_SWITCH(FlowFlucGfwPp, processq2, "Process analysis for filling q_n-vector calibration histograms", true);
1557+
1558+
void processMC(soa::Filtered<soa::Join<aod::Collisions, aod::EvSels, aod::Mults,
1559+
aod::CentFT0Cs, aod::CentFT0CVariant1s, aod::CentFT0Ms,
1560+
aod::CentFV0As, aod::CentNTPVs, aod::CentNGlobals,
1561+
aod::McCollisionLabels>>::iterator const& collision,
1562+
aod::BCsWithTimestamps const&, GFWTracksMC const& tracks, aod::McCollisions const&, aod::McParticles const& mcParticles)
1563+
{
1564+
auto bc = collision.bc_as<aod::BCsWithTimestamps>();
1565+
int run = bc.runNumber();
1566+
if (run != lastRun) {
1567+
lastRun = run;
1568+
LOGF(info, "run = %d", run);
1569+
if (cfgRunByRun) {
1570+
if (std::find(runNumbers.begin(), runNumbers.end(), run) == runNumbers.end()) {
1571+
LOGF(info, "Creating histograms for run %d", run);
1572+
createRunByRunHistograms(run);
1573+
runNumbers.push_back(run);
1574+
}
1575+
if (!cfgFillWeights)
1576+
loadCorrections(bc);
1577+
}
1578+
}
1579+
if (!cfgFillWeights && !cfgRunByRun)
1580+
loadCorrections(bc);
1581+
1582+
registry.fill(HIST("eventQA/eventSel"), kFilteredEvent);
1583+
if (cfgRunByRun)
1584+
th1sList[run][hEventSel]->Fill(kFilteredEvent);
1585+
1586+
if (!collision.sel8())
1587+
return;
1588+
1589+
registry.fill(HIST("eventQA/eventSel"), kSel8);
1590+
if (cfgRunByRun)
1591+
th1sList[run][hEventSel]->Fill(kSel8);
1592+
1593+
if (cfgDoOccupancySel) {
1594+
int occupancy = collision.trackOccupancyInTimeRange();
1595+
if (occupancy < 0 || occupancy > cfgOccupancySelection)
1596+
return;
1597+
}
1598+
1599+
registry.fill(HIST("eventQA/eventSel"), kOccupancy);
1600+
if (cfgRunByRun)
1601+
th1sList[run][hEventSel]->Fill(kOccupancy);
1602+
1603+
const XAxis xaxis{
1604+
getCentrality(collision),
1605+
collision.multNTracksPV(),
1606+
(cfgTimeDependent) ? getTimeSinceStartOfFill(bc.timestamp(), *firstRunOfCurrentFill) : -1.0};
1607+
1608+
if (cfgTimeDependent && run == *firstRunOfCurrentFill &&
1609+
firstRunOfCurrentFill != o2::analysis::gfwflowflucpp::firstRunsOfFill.end() - 1)
1610+
++firstRunOfCurrentFill;
1611+
1612+
if (cfgFillQA)
1613+
fillEventQA<kBefore>(collision, xaxis);
1614+
1615+
registry.fill(HIST("eventQA/before/centrality"), xaxis.centrality);
1616+
registry.fill(HIST("eventQA/before/multiplicity"), xaxis.multiplicity);
1617+
1618+
if (!eventSelected(collision, xaxis.multiplicity, xaxis.centrality, run))
1619+
return;
1620+
1621+
if (cfgFillQA)
1622+
fillEventQA<kAfter>(collision, xaxis);
1623+
1624+
processCollision<kReco>(collision, tracks, xaxis, run, 0);
1625+
1626+
if (!collision.has_mcCollision())
1627+
return;
1628+
1629+
const auto genCollision = collision.template mcCollision_as<aod::McCollisions>();
1630+
processGenCollision(genCollision, mcParticles, collision.mcCollisionId(), xaxis, run, 0);
1631+
}
1632+
PROCESS_SWITCH(FlowFlucGfwPp, processMC, "Process analysis for Monte-Carlo data", false);
15241633
};
15251634

15261635
WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)

0 commit comments

Comments
 (0)