Skip to content

Commit e4e94eb

Browse files
LuzhiyonggLuzhiyongg
authored
[PWGCF] add MC process (#11482)
1 parent 98d6b15 commit e4e94eb

File tree

2 files changed

+289
-76
lines changed

2 files changed

+289
-76
lines changed

‎PWGCF/Flow/Tasks/flowTask.cxx‎

Lines changed: 115 additions & 28 deletions
Original file line numberDiff line numberDiff line change
@@ -113,6 +113,16 @@ struct FlowTask {
113113

114114
Filter collisionFilter = (nabs(aod::collision::posZ) < cfgCutVertex) && (aod::cent::centFT0C > cfgCentFT0CMin) && (aod::cent::centFT0C < cfgCentFT0CMax);
115115
Filter trackFilter = ((requireGlobalTrackInFilter()) || (aod::track::isGlobalTrackSDD == (uint8_t) true)) && (nabs(aod::track::eta) < cfgCutEta) && (aod::track::pt > cfgCutPtMin) && (aod::track::pt < cfgCutPtMax) && (aod::track::tpcChi2NCl < cfgCutChi2prTPCcls) && (nabs(aod::track::dcaZ) < cfgCutDCAz);
116+
using FilteredCollisions = soa::Filtered<soa::Join<aod::Collisions, aod::EvSels, aod::CentFT0Cs, aod::CentFT0CVariant1s, aod::CentFT0Ms, aod::CentFV0As, aod::Mults>>;
117+
using FilteredTracks = soa::Filtered<soa::Join<aod::Tracks, aod::TrackSelection, aod::TracksExtra, aod::TracksDCA>>;
118+
// Filter for MCcollisions
119+
Filter mccollisionFilter = nabs(aod::mccollision::posZ) < cfgCutVertex;
120+
using FilteredMcCollisions = soa::Filtered<aod::McCollisions>;
121+
// Filter for MCParticle
122+
Filter particleFilter = (nabs(aod::mcparticle::eta) < cfgCutEta) && (aod::mcparticle::pt > cfgCutPtMin) && (aod::mcparticle::pt < cfgCutPtMax);
123+
using FilteredMcParticles = soa::Filtered<aod::McParticles>;
124+
125+
using FilteredSmallGroupMcCollisions = soa::SmallGroups<soa::Join<aod::McCollisionLabels, aod::Collisions, aod::EvSel, aod::CentFT0Cs, aod::CentFT0CVariant1s, aod::CentFT0Ms, aod::CentFV0As, aod::Mults>>;
116126

117127
// Corrections
118128
TH1D* mEfficiency = nullptr;
@@ -127,6 +137,7 @@ struct FlowTask {
127137

128138
// Define output
129139
OutputObj<FlowContainer> fFC{FlowContainer("FlowContainer")};
140+
OutputObj<FlowContainer> fFCgen{FlowContainer("FlowContainer_gen")};
130141
OutputObj<GFWWeights> fWeights{GFWWeights("weights")};
131142
HistogramRegistry registry{"registry"};
132143

@@ -143,6 +154,10 @@ struct FlowTask {
143154
// Count the total number of enum
144155
kCount_CentEstimators
145156
};
157+
enum DataType {
158+
kReco,
159+
kGen
160+
};
146161
int mRunNumber{-1};
147162
uint64_t mSOR{0};
148163
double mMinSeconds{-1.};
@@ -157,9 +172,6 @@ struct FlowTask {
157172
TF1* funcV3;
158173
TF1* funcV4;
159174

160-
using AodCollisions = soa::Filtered<soa::Join<aod::Collisions, aod::EvSels, aod::CentFT0Cs, aod::CentFT0CVariant1s, aod::CentFT0Ms, aod::CentFV0As, aod::Mults>>;
161-
using AodTracks = soa::Filtered<soa::Join<aod::Tracks, aod::TrackSelection, aod::TracksExtra, aod::TracksDCA>>;
162-
163175
// Track selection
164176
TrackSelection myTrackSel;
165177
TF1* fPhiCutLow = nullptr;
@@ -222,6 +234,11 @@ struct FlowTask {
222234
registry.add("hMult", "Multiplicity distribution", {HistType::kTH1D, {{3000, 0.5, 3000.5}}});
223235
std::string hCentTitle = "Centrality distribution, Estimator " + std::to_string(cfgCentEstimator);
224236
registry.add("hCent", hCentTitle.c_str(), {HistType::kTH1D, {{90, 0, 90}}});
237+
if (doprocessMCGen) {
238+
registry.add("MCGen/MChVtxZ", "Vexter Z distribution", {HistType::kTH1D, {axisVertex}});
239+
registry.add("MCGen/MChMult", "Multiplicity distribution", {HistType::kTH1D, {{3000, 0.5, 3000.5}}});
240+
registry.add("MCGen/MChCent", hCentTitle.c_str(), {HistType::kTH1D, {{90, 0, 90}}});
241+
}
225242
if (!cfgUseSmallMemory) {
226243
registry.add("BeforeSel8_globalTracks_centT0C", "before sel8;Centrality T0C;mulplicity global tracks", {HistType::kTH2D, {axisCentForQA, axisNch}});
227244
registry.add("BeforeCut_globalTracks_centT0C", "before cut;Centrality T0C;mulplicity global tracks", {HistType::kTH2D, {axisCentForQA, axisNch}});
@@ -256,6 +273,11 @@ struct FlowTask {
256273
registry.add("hDCAz", "DCAz after cuts; DCAz (cm); Pt", {HistType::kTH2D, {{200, -0.5, 0.5}, {200, 0, 5}}});
257274
registry.add("hDCAxy", "DCAxy after cuts; DCAxy (cm); Pt", {HistType::kTH2D, {{200, -0.5, 0.5}, {200, 0, 5}}});
258275
registry.add("hTrackCorrection2d", "Correlation table for number of tracks table; uncorrected track; corrected track", {HistType::kTH2D, {axisNch, axisNch}});
276+
if (doprocessMCGen) {
277+
registry.add("MCGen/MChPhi", "#phi distribution", {HistType::kTH1D, {axisPhi}});
278+
registry.add("MCGen/MChEta", "#eta distribution", {HistType::kTH1D, {axisEta}});
279+
registry.add("MCGen/MChPtRef", "p_{T} distribution after cut", {HistType::kTH1D, {axisPtHist}});
280+
}
259281

260282
o2::framework::AxisSpec axis = axisPt;
261283
int nPtBins = axis.binEdges.size() - 1;
@@ -327,6 +349,11 @@ struct FlowTask {
327349
fFC->SetName("FlowContainer");
328350
fFC->SetXAxis(fPtAxis);
329351
fFC->Initialize(oba, axisIndependent, cfgNbootstrap);
352+
if (doprocessMCGen) {
353+
fFCgen->SetName("FlowContainer_gen");
354+
fFCgen->SetXAxis(fPtAxis);
355+
fFCgen->Initialize(oba, axisIndependent, cfgNbootstrap);
356+
}
330357
delete oba;
331358

332359
// eta region
@@ -471,6 +498,7 @@ struct FlowTask {
471498
return;
472499
}
473500

501+
template <DataType dt>
474502
void fillFC(const GFW::CorrConfig& corrconf, const double& cent, const double& rndm)
475503
{
476504
double dnx, val;
@@ -479,17 +507,19 @@ struct FlowTask {
479507
if (dnx == 0)
480508
return;
481509
val = fGFW->Calculate(corrconf, 0, kFALSE).real() / dnx;
482-
if (std::fabs(val) < 1)
483-
fFC->FillProfile(corrconf.Head.c_str(), cent, val, dnx, rndm);
510+
if (std::fabs(val) < 1) {
511+
(dt == kGen) ? fFCgen->FillProfile(corrconf.Head.c_str(), cent, val, dnx, rndm) : fFC->FillProfile(corrconf.Head.c_str(), cent, val, dnx, rndm);
512+
}
484513
return;
485514
}
486515
for (auto i = 1; i <= fPtAxis->GetNbins(); i++) {
487516
dnx = fGFW->Calculate(corrconf, i - 1, kTRUE).real();
488517
if (dnx == 0)
489518
continue;
490519
val = fGFW->Calculate(corrconf, i - 1, kFALSE).real() / dnx;
491-
if (std::fabs(val) < 1)
492-
fFC->FillProfile(Form("%s_pt_%i", corrconf.Head.c_str(), i), cent, val, dnx, rndm);
520+
if (std::fabs(val) < 1) {
521+
(dt == kGen) ? fFCgen->FillProfile(Form("%s_pt_%i", corrconf.Head.c_str(), i), cent, val, dnx, rndm) : fFC->FillProfile(Form("%s_pt_%i", corrconf.Head.c_str(), i), cent, val, dnx, rndm);
522+
}
493523
}
494524
return;
495525
}
@@ -610,7 +640,8 @@ struct FlowTask {
610640
registry.fill(HIST("hEventCountSpecific"), 8.5);
611641

612642
// V0A T0A 5 sigma cut
613-
if (cfgEvSelV0AT0ACut && (std::fabs(collision.multFV0A() - fT0AV0AMean->Eval(collision.multFT0A())) > 5 * fT0AV0ASigma->Eval(collision.multFT0A())))
643+
float sigma = 5.0;
644+
if (cfgEvSelV0AT0ACut && (std::fabs(collision.multFV0A() - fT0AV0AMean->Eval(collision.multFT0A())) > sigma * fT0AV0ASigma->Eval(collision.multFT0A())))
614645
return 0;
615646
if (cfgEvSelV0AT0ACut)
616647
registry.fill(HIST("hEventCountSpecific"), 9.5);
@@ -641,7 +672,8 @@ struct FlowTask {
641672
registry.fill(HIST("hEventCountTentative"), 7.5);
642673
if (!((multNTracksPV < fMultPVCutLow->Eval(centrality)) || (multNTracksPV > fMultPVCutHigh->Eval(centrality)) || (multTrk < fMultCutLow->Eval(centrality)) || (multTrk > fMultCutHigh->Eval(centrality))))
643674
registry.fill(HIST("hEventCountTentative"), 8.5);
644-
if (!(std::fabs(collision.multFV0A() - fT0AV0AMean->Eval(collision.multFT0A())) > 5 * fT0AV0ASigma->Eval(collision.multFT0A())))
675+
float sigma = 5.0;
676+
if (!(std::fabs(collision.multFV0A() - fT0AV0AMean->Eval(collision.multFT0A())) > sigma * fT0AV0ASigma->Eval(collision.multFT0A())))
645677
registry.fill(HIST("hEventCountTentative"), 9.5);
646678
}
647679

@@ -675,7 +707,30 @@ struct FlowTask {
675707
gCurrentHadronicRate = gHadronicRate[mRunNumber];
676708
}
677709

678-
void process(AodCollisions::iterator const& collision, aod::BCsWithTimestamps const&, AodTracks const& tracks)
710+
template <typename TCollision>
711+
float getCentrality(TCollision const& collision)
712+
{
713+
float cent;
714+
switch (cfgCentEstimator) {
715+
case kCentFT0C:
716+
cent = collision.centFT0C();
717+
break;
718+
case kCentFT0CVariant1:
719+
cent = collision.centFT0CVariant1();
720+
break;
721+
case kCentFT0M:
722+
cent = collision.centFT0M();
723+
break;
724+
case kCentFV0A:
725+
cent = collision.centFV0A();
726+
break;
727+
default:
728+
cent = collision.centFT0C();
729+
}
730+
return cent;
731+
}
732+
733+
void process(FilteredCollisions::iterator const& collision, aod::BCsWithTimestamps const&, FilteredTracks const& tracks)
679734
{
680735
registry.fill(HIST("hEventCount"), 0.5);
681736
if (!cfgUseSmallMemory && tracks.size() >= 1) {
@@ -703,23 +758,8 @@ struct FlowTask {
703758
registry.fill(HIST("BeforeCut_multV0A_multT0A"), collision.multFT0A(), collision.multFV0A());
704759
registry.fill(HIST("BeforeCut_multT0C_centT0C"), collision.centFT0C(), collision.multFT0C());
705760
}
706-
float cent;
707-
switch (cfgCentEstimator) {
708-
case kCentFT0C:
709-
cent = collision.centFT0C();
710-
break;
711-
case kCentFT0CVariant1:
712-
cent = collision.centFT0CVariant1();
713-
break;
714-
case kCentFT0M:
715-
cent = collision.centFT0M();
716-
break;
717-
case kCentFV0A:
718-
cent = collision.centFV0A();
719-
break;
720-
default:
721-
cent = collision.centFT0C();
722-
}
761+
float cent = getCentrality(collision);
762+
723763
if (cfgUseTentativeEventCounter)
724764
eventCounterQA(collision, tracks.size(), cent);
725765
if (cfgUseAdditionalEventCut && !eventSelected(collision, tracks.size(), cent))
@@ -844,9 +884,56 @@ struct FlowTask {
844884

845885
// Filling Flow Container
846886
for (uint l_ind = 0; l_ind < corrconfigs.size(); l_ind++) {
847-
fillFC(corrconfigs.at(l_ind), independent, lRandom);
887+
fillFC<kReco>(corrconfigs.at(l_ind), independent, lRandom);
888+
}
889+
}
890+
891+
void processMCGen(FilteredMcCollisions::iterator const& mcCollision, FilteredSmallGroupMcCollisions const& collisions, FilteredMcParticles const& mcParticles)
892+
{
893+
if (collisions.size() != 1)
894+
return;
895+
896+
float cent = -1.;
897+
for (const auto& collision : collisions) {
898+
cent = getCentrality(collision);
899+
}
900+
901+
float lRandom = fRndm->Rndm();
902+
float vtxz = mcCollision.posZ();
903+
registry.fill(HIST("MCGen/MChVtxZ"), vtxz);
904+
registry.fill(HIST("MCGen/MChMult"), mcParticles.size());
905+
registry.fill(HIST("MCGen/MChCent"), cent);
906+
float independent = cent;
907+
if (cfgUseNch)
908+
independent = static_cast<float>(mcParticles.size());
909+
910+
fGFW->Clear();
911+
912+
for (const auto& mcParticle : mcParticles) {
913+
if (!mcParticle.isPhysicalPrimary())
914+
continue;
915+
bool withinPtPOI = (cfgCutPtPOIMin < mcParticle.pt()) && (mcParticle.pt() < cfgCutPtPOIMax); // within POI pT range
916+
bool withinPtRef = (cfgCutPtRefMin < mcParticle.pt()) && (mcParticle.pt() < cfgCutPtRefMax); // within RF pT range
917+
918+
if (withinPtRef) {
919+
registry.fill(HIST("MCGen/MChPhi"), mcParticle.phi());
920+
registry.fill(HIST("MCGen/MChEta"), mcParticle.eta());
921+
registry.fill(HIST("MCGen/MChPtRef"), mcParticle.pt());
922+
}
923+
if (withinPtRef)
924+
fGFW->Fill(mcParticle.eta(), fPtAxis->FindBin(mcParticle.pt()) - 1, mcParticle.phi(), 1., 1);
925+
if (withinPtPOI)
926+
fGFW->Fill(mcParticle.eta(), fPtAxis->FindBin(mcParticle.pt()) - 1, mcParticle.phi(), 1., 2);
927+
if (withinPtPOI && withinPtRef)
928+
fGFW->Fill(mcParticle.eta(), fPtAxis->FindBin(mcParticle.pt()) - 1, mcParticle.phi(), 1., 4);
929+
}
930+
931+
// Filling Flow Container
932+
for (uint l_ind = 0; l_ind < corrconfigs.size(); l_ind++) {
933+
fillFC<kGen>(corrconfigs.at(l_ind), independent, lRandom);
848934
}
849935
}
936+
PROCESS_SWITCH(FlowTask, processMCGen, "Process analysis for MC generated events", false);
850937
};
851938

852939
WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)

0 commit comments

Comments
 (0)