diff --git a/PWGCF/Flow/Tasks/flowMc.cxx b/PWGCF/Flow/Tasks/flowMc.cxx index ea3223dc153..e370a7a9847 100644 --- a/PWGCF/Flow/Tasks/flowMc.cxx +++ b/PWGCF/Flow/Tasks/flowMc.cxx @@ -78,7 +78,7 @@ struct FlowMc { O2_DEFINE_CONFIGURABLE(cfgFlowEfficiency, std::string, "", "CCDB path to efficiency object") O2_DEFINE_CONFIGURABLE(cfgCentVsIPTruth, std::string, "", "CCDB path to centrality vs IP truth") O2_DEFINE_CONFIGURABLE(cfgIsGlobalTrack, bool, false, "Use global tracks instead of hasTPC&&hasITS") - O2_DEFINE_CONFIGURABLE(cfgK0Lambda0Enabled, bool, false, "Add K0 and Lambda0") + O2_DEFINE_CONFIGURABLE(cfgK0Lambda0Enabled, bool, true, "Add K0 and Lambda0") O2_DEFINE_CONFIGURABLE(cfgFlowCumulantEnabled, bool, false, "switch of calculating flow") O2_DEFINE_CONFIGURABLE(cfgFlowCumulantNbootstrap, int, 30, "Number of subsamples") O2_DEFINE_CONFIGURABLE(cfgTrackDensityCorrUse, bool, false, "Use track density efficiency correction") @@ -96,7 +96,6 @@ struct FlowMc { Configurable> cfgTrackDensityP0{"cfgTrackDensityP0", std::vector{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> cfgTrackDensityP1{"cfgTrackDensityP1", std::vector{-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; ConfigurableAxis axisB{"axisB", {100, 0.0f, 20.0f}, ""}; ConfigurableAxis axisCentrality{"axisCentrality", {VARIABLE_WIDTH, 0, 5, 10, 20, 30, 40, 50, 60, 70, 80, 90}, "X axis for histograms"}; @@ -162,9 +161,10 @@ struct FlowMc { const AxisSpec axisEta{20, -1., 1., "#eta"}; const AxisSpec axisCounter{1, 0, +1, ""}; // QA histograms - histos.add("mcEventCounter", "Monte Carlo Truth EventCounter", HistType::kTH1F, {axisCounter}); - histos.add("numberOfRecoCollisions", "numberOfRecoCollisions", HistType::kTH1F, {{10, -0.5f, 9.5f}}); - histos.add("RecoEventCounter", "Reconstruction EventCounter", HistType::kTH1F, {axisCounter}); + histos.add("mcEventCounter", "Monte Carlo Truth EventCounter", HistType::kTH1F, {{5, 0, 5}}); + histos.add("numberOfRecoCollisions", "numberOfRecoCollisions", HistType::kTH1F, {{10, -1.5f, 8.5f}}); + histos.add("RecoEventCounter", "Reconstruction EventCounter", HistType::kTH1F, {{5, 0, 5}}); + histos.add("RecoProcessEventCounter", "Reconstruction EventCounter", HistType::kTH1F, {{5, 0, 5}}); histos.add("hnTPCClu", "Number of found TPC clusters", HistType::kTH1D, {{100, 40, 180}}); histos.add("hnITSClu", "Number of found ITS clusters", HistType::kTH1D, {{100, 0, 20}}); // pT histograms @@ -209,6 +209,8 @@ struct FlowMc { histos.add("hPtMCGen", "Monte Carlo Truth; pT (GeV/c);", {HistType::kTH1D, {axisPt}}); histos.add("hEtaPtVtxzMCGen", "Monte Carlo Truth; #eta; p_{T} (GeV/c); V_{z} (cm);", {HistType::kTH3D, {axisEta, axisPt, axisVertex}}); histos.add("hPtMCGlobal", "Monte Carlo Global; pT (GeV/c);", {HistType::kTH1D, {axisPt}}); + histos.add("hPtReco", "Monte Carlo Reco Global; pT (GeV/c);", {HistType::kTH1D, {axisPt}}); + histos.add("hPtReco_PhysicalPrimary", "Monte Carlo Reco Global; pT (GeV/c);", {HistType::kTH1D, {axisPt}}); histos.add("hEtaPtVtxzMCGlobal", "Monte Carlo Global; #eta; p_{T} (GeV/c); V_{z} (cm);", {HistType::kTH3D, {axisEta, axisPt, axisVertex}}); histos.add("hPhiWeightedTrDen", "corrected #phi distribution, considering track density", {HistType::kTH1D, {axisPhi}}); @@ -447,9 +449,9 @@ struct FlowMc { return myTrackSel.IsSelected(track); } - void process(FilteredMcCollisions::iterator const& mcCollision, aod::BCsWithTimestamps const&, soa::SmallGroups> const& collisions, FilteredMcParticles const& mcParticles, FilteredTracks const&) + void processMCTrue(FilteredMcCollisions::iterator const& mcCollision, aod::BCsWithTimestamps const&, soa::SmallGroups> const& collisions, FilteredMcParticles const& mcParticles, FilteredTracks const&) { - + histos.fill(HIST("mcEventCounter"), 0.5); float imp = mcCollision.impactParameter(); float evPhi = mcCollision.eventPlaneAngle(); float vtxz = mcCollision.posZ(); @@ -465,21 +467,23 @@ struct FlowMc { 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 + histos.fill(HIST("RecoEventCounter"), 0.5); if (cfgRecoEvRejectMC) { if (collisions.size() != 1) { // only pass those have one reconstruction event return; } + histos.fill(HIST("RecoEventCounter"), 1.5); for (auto const& collision : collisions) { if (!eventSelected(collision)) return; } + histos.fill(HIST("RecoEventCounter"), 2.5); } - histos.fill(HIST("RecoEventCounter"), 0.5); } + histos.fill(HIST("RecoEventCounter"), 3.5); - if (imp > minB && imp < maxB) { + if (imp >= minB && imp <= maxB) { // event within range histos.fill(HIST("hImpactParameter"), imp); histos.fill(HIST("hEventPlaneAngle"), evPhi); @@ -502,7 +506,7 @@ struct FlowMc { continue; if (!mcParticle.isPhysicalPrimary()) continue; - if (std::fabs(mcParticle.eta()) > maxEta) // main acceptance + if (std::fabs(mcParticle.eta()) > cfgCutEta) // main acceptance continue; if (mcParticle.has_tracks()) { auto const& tracks = mcParticle.tracks_as(); @@ -551,7 +555,7 @@ struct FlowMc { if (!mcParticle.isPhysicalPrimary()) continue; - if (std::fabs(mcParticle.eta()) > maxEta) // main acceptance + if (std::fabs(mcParticle.eta()) > cfgCutEta) // main acceptance continue; float deltaPhi = mcParticle.phi() - mcCollision.eventPlaneAngle(); @@ -696,6 +700,37 @@ struct FlowMc { } histos.fill(HIST("hNchVsImpactParameter"), imp, nCh); } + PROCESS_SWITCH(FlowMc, processMCTrue, "process pure simulation information", true); + + void processReco(soa::Join::iterator const& collision, aod::BCsWithTimestamps const&, FilteredTracks const& tracks, aod::McParticles const&, aod::McCollisions const&) + { + histos.fill(HIST("RecoProcessEventCounter"), 0.5); + if (!eventSelected(collision)) + return; + histos.fill(HIST("RecoProcessEventCounter"), 1.5); + if (tracks.size() < 1) + return; + histos.fill(HIST("RecoProcessEventCounter"), 2.5); + for (const auto& track : tracks) { + if (!trackSelected(track)) + continue; + histos.fill(HIST("hPtReco"), track.pt()); + auto mcParticle = track.mcParticle(); + int pdgCode = std::abs(mcParticle.pdgCode()); + bool extraPDGType = true; + if (cfgK0Lambda0Enabled) { + extraPDGType = (pdgCode != PDG_t::kK0Short && pdgCode != PDG_t::kLambda0); + } + if (extraPDGType && pdgCode != PDG_t::kElectron && pdgCode != PDG_t::kMuonMinus && pdgCode != PDG_t::kPiPlus && pdgCode != kKPlus && pdgCode != PDG_t::kProton) + continue; + if (!mcParticle.isPhysicalPrimary()) + continue; + if (std::fabs(mcParticle.eta()) > cfgCutEta) // main acceptance + continue; + histos.fill(HIST("hPtReco_PhysicalPrimary"), track.pt()); + } + } + PROCESS_SWITCH(FlowMc, processReco, "process reconstructed information", true); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) diff --git a/PWGCF/Flow/Tasks/flowTask.cxx b/PWGCF/Flow/Tasks/flowTask.cxx index da88271ee0b..a38fac52843 100644 --- a/PWGCF/Flow/Tasks/flowTask.cxx +++ b/PWGCF/Flow/Tasks/flowTask.cxx @@ -41,6 +41,7 @@ #include "TList.h" #include #include +#include #include #include @@ -300,6 +301,7 @@ struct FlowTask { std::string hCentTitle = "Centrality distribution, Estimator " + std::to_string(cfgCentEstimator); registry.add("hCent", hCentTitle.c_str(), {HistType::kTH1D, {{100, 0, 100}}}); if (doprocessMCGen) { + registry.add("MCGen/hMCEventCount", "Number of Event;; Count", {HistType::kTH1D, {{5, 0, 5}}}); registry.add("MCGen/MChVtxZ", "Vexter Z distribution", {HistType::kTH1D, {axisVertex}}); registry.add("MCGen/MChMult", "Multiplicity distribution", {HistType::kTH1D, {{3000, 0.5, 3000.5}}}); registry.add("MCGen/MChCent", hCentTitle.c_str(), {HistType::kTH1D, {{100, 0, 100}}}); @@ -360,7 +362,7 @@ struct FlowTask { if (doprocessMCGen) { registry.add("MCGen/MChPhi", "#phi distribution", {HistType::kTH1D, {axisPhi}}); registry.add("MCGen/MChEta", "#eta distribution", {HistType::kTH1D, {axisEta}}); - registry.add("MCGen/MChPtRef", "p_{T} distribution after cut", {HistType::kTH1D, {axisPtHist}}); + registry.add("MCGen/MChPt", "p_{T} distribution after cut", {HistType::kTH1D, {axisPtHist}}); registry.add("hMeanPtWithinGap08_MC", "mean p_{T}", {HistType::kTProfile, {axisIndependent}}); for (auto i = 0; i < cfgNbootstrap; i++) { bootstrapArray[i][kMeanPtWithinGap08_MC] = registry.add(Form("BootstrapContainer_%d/hMeanPtWithinGap08_MC", i), "", {HistType::kTProfile, {axisIndependent}}); @@ -1242,14 +1244,26 @@ struct FlowTask { void processMCGen(FilteredMcCollisions::iterator const& mcCollision, FilteredSmallGroupMcCollisions const& collisions, FilteredMcParticles const& mcParticles) { + registry.fill(HIST("MCGen/hMCEventCount"), 0.5); if (collisions.size() != 1) return; + registry.fill(HIST("MCGen/hMCEventCount"), 1.5); float cent = -1.; for (const auto& collision : collisions) { cent = getCentrality(collision); } + if (cfgUseAdditionalEventCut) { + for (auto const& collision : collisions) { + if (!collision.sel8()) + return; + if (!eventSelected(collision, mcParticles.size(), cent)) + return; + } + } + registry.fill(HIST("MCGen/hMCEventCount"), 2.5); + float lRandom = fRndm->Rndm(); float vtxz = mcCollision.posZ(); registry.fill(HIST("MCGen/MChVtxZ"), vtxz); @@ -1262,26 +1276,25 @@ struct FlowTask { fGFW->Clear(); fFCptgen->clearVector(); - if (cfgUseAdditionalEventCut) { - for (auto const& collision : collisions) { - if (!eventSelected(collision, mcParticles.size(), cent)) - return; - } - } - double ptSum_Gap08 = 0.; double count_Gap08 = 0.; for (const auto& mcParticle : mcParticles) { + int pdgCode = std::abs(mcParticle.pdgCode()); + bool extraPDGType = (pdgCode != PDG_t::kK0Short && pdgCode != PDG_t::kLambda0); + if (extraPDGType && pdgCode != PDG_t::kElectron && pdgCode != PDG_t::kMuonMinus && pdgCode != PDG_t::kPiPlus && pdgCode != kKPlus && pdgCode != PDG_t::kProton) + continue; if (!mcParticle.isPhysicalPrimary()) continue; + if (std::fabs(mcParticle.eta()) > cfgCutEta) // main acceptance + continue; bool withinPtPOI = (cfgCutPtPOIMin < mcParticle.pt()) && (mcParticle.pt() < cfgCutPtPOIMax); // within POI pT range bool withinPtRef = (cfgCutPtRefMin < mcParticle.pt()) && (mcParticle.pt() < cfgCutPtRefMax); // within RF pT range bool withinEtaGap08 = (std::abs(mcParticle.eta()) < cfgCutEta); + registry.fill(HIST("MCGen/MChPt"), mcParticle.pt()); if (withinPtRef) { registry.fill(HIST("MCGen/MChPhi"), mcParticle.phi()); registry.fill(HIST("MCGen/MChEta"), mcParticle.eta()); - registry.fill(HIST("MCGen/MChPtRef"), mcParticle.pt()); if (withinEtaGap08) { ptSum_Gap08 += mcParticle.pt(); count_Gap08 += 1.;