From 6a473966c6afed8ed086dc2b35d36c90cf5f61ac Mon Sep 17 00:00:00 2001 From: balwan Date: Fri, 16 Oct 2026 14:27:30 +0530 Subject: [PATCH] merge-safe slot subsampling replaces the sampleSize path --- .../Tasks/FactorialMomentsTask.cxx | 715 ++++++++---------- 1 file changed, 330 insertions(+), 385 deletions(-) diff --git a/PWGCF/EbyEFluctuations/Tasks/FactorialMomentsTask.cxx b/PWGCF/EbyEFluctuations/Tasks/FactorialMomentsTask.cxx index 26331d34127..8bca9cda29e 100644 --- a/PWGCF/EbyEFluctuations/Tasks/FactorialMomentsTask.cxx +++ b/PWGCF/EbyEFluctuations/Tasks/FactorialMomentsTask.cxx @@ -8,8 +8,8 @@ // In applying this license CERN does not waive the privileges and immunities // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. -/// -/// \brief This task is for Normalized Factorial Moments Analysis: PRC 85, 044914 (2012), nucl-ex:1411.6083 +/// \file FactorialMomentsTask.cxx +/// \brief This task is for Normalized Factorial Moments Analysis: Hwa and Yang, Phys. Rev. C 85, 044914 (2012) /// \author Salman Malik /// \author Balwan Singh @@ -29,101 +29,71 @@ #include #include #include +#include #include -#include #include #include #include #include -#include #include #include -#include #include #include #include -#include +#include #include #include using namespace o2; using namespace o2::framework; using namespace o2::framework::expressions; + struct FactorialMomentsTask { Configurable useITS{"useITS", false, "Select tracks with ITS"}; Configurable useTPC{"useTPC", false, "Select tracks with TPC"}; Configurable useGlobal{"useGlobal", true, "Select global tracks"}; Configurable applyCheckPtForRec{"applyCheckPtForRec", false, "Apply checkpT for reconstructed tracks"}; Configurable applyCheckPtForMC{"applyCheckPtForMC", true, "Apply checkpT for MC-generated tracks"}; - Configurable cfgEvSelkNoITSROFrameBorder{"cfgEvSelkNoITSROFrameBorder", true, "ITSROFrame border event selection cut"}; - Configurable cfgEvSelkNoTimeFrameBorder{"cfgEvSelkNoTimeFrameBorder", true, "TimeFrame border event selection cut"}; + Configurable smearPhi{"smearPhi", false, "Randomise track phi with a Gaussian of sigma=2pi before filling the eta-phi lattices"}; + Configurable evSelNoITSROFrameBorder{"evSelNoITSROFrameBorder", true, "ITSROFrame border event selection cut"}; + Configurable evSelNoTimeFrameBorder{"evSelNoTimeFrameBorder", true, "TimeFrame border event selection cut"}; Configurable centralEta{"centralEta", 0.9, "eta limit for tracks"}; Configurable numPt{"numPt", 5, "number of pT bins"}; Configurable ptMin{"ptMin", 0.2f, "lower pT cut"}; Configurable dcaXY{"dcaXY", 0.1f, "DCA xy cut"}; Configurable dcaZ{"dcaZ", 1.0f, "DCA z cut"}; - Configurable cfgCutTpcChi2NCl{"cfgCutTpcChi2NCl", 2.5f, "Maximum TPCchi2NCl"}; - Configurable cfgCutItsChi2NCl{"cfgCutItsChi2NCl", 40.0f, "Maximum ITSchi2NCl"}; - Configurable mintPCCls{"mintPCCls", 70.0f, "minimum number of TPC clusters"}; + Configurable maxTPCChi2NCl{"maxTPCChi2NCl", 2.5f, "Maximum TPCchi2NCl"}; + Configurable maxITSChi2NCl{"maxITSChi2NCl", 40.0f, "Maximum ITSchi2NCl"}; + Configurable minTPCNCls{"minTPCNCls", 70.0f, "minimum number of TPC clusters"}; Configurable> centLimits{"centLimits", {0, 5}, "centrality min and max"}; Configurable> vertexXYZ{"vertexXYZ", {0.3f, 0.4f, 10.0f}, "vertex cuts"}; Configurable> ptCuts{"ptCuts", {0.2f, 2.0f}, "pT cuts"}; - Configurable isApplySameBunchPileup{"isApplySameBunchPileup", true, "Enable SameBunchPileup cut"}; - Configurable isApplyGoodZvtxFT0vsPV{"isApplyGoodZvtxFT0vsPV", true, "Enable GoodZvtxFT0vsPV cut"}; - Configurable cfgUseGoodITSLayerAllCut{"cfgUseGoodITSLayerAllCut", true, "Remove time interval with dead ITS zone"}; - Configurable isApplyVertexITSTPC{"isApplyVertexITSTPC", true, "Enable VertexITSTPC cut"}; - Configurable isApplyVertexTOFmatched{"isApplyVertexTOFmatched", true, "Enable VertexTOFmatched cut"}; - Configurable isApplyVertexTRDmatched{"isApplyVertexTRDmatched", true, "Enable VertexTRDmatched cut"}; - Configurable isApplyExtraCorrCut{"isApplyExtraCorrCut", false, "Enable extra NPVtracks vs FTOC correlation cut"}; - Configurable isApplyExtraPhiCut{"isApplyExtraPhiCut", false, "Enable extra phi cut"}; + Configurable evSelNoSameBunchPileup{"evSelNoSameBunchPileup", true, "Enable SameBunchPileup cut"}; + Configurable evSelGoodZvtxFT0vsPV{"evSelGoodZvtxFT0vsPV", true, "Enable GoodZvtxFT0vsPV cut"}; + Configurable evSelGoodITSLayersAll{"evSelGoodITSLayersAll", true, "Remove time interval with dead ITS zone"}; + Configurable evSelVertexITSTPC{"evSelVertexITSTPC", true, "Enable VertexITSTPC cut"}; + Configurable evSelVertexTOFmatched{"evSelVertexTOFmatched", true, "Enable VertexTOFmatched cut"}; + Configurable evSelVertexTRDmatched{"evSelVertexTRDmatched", true, "Enable VertexTRDmatched cut"}; + Configurable applyExtraCorrCut{"applyExtraCorrCut", false, "Enable extra NPVtracks vs FTOC correlation cut"}; + Configurable applyExtraPhiCut{"applyExtraPhiCut", false, "Enable extra phi cut"}; Configurable includeGlobalTracks{"includeGlobalTracks", false, "Enable Global Tracks"}; Configurable includeTPCTracks{"includeTPCTracks", false, "TPC Tracks"}; Configurable includeITSTracks{"includeITSTracks", false, "ITS Tracks"}; - Configurable samplesize{"samplesize", 100, "Sample size"}; + Configurable nSlots{"nSlots", 1250, "Number of hash slots for merge-safe subsampling (slot = splitmix64(globalBC) % nSlots); must be identical for every job"}; Configurable useMC{"useMC", false, "Use MC information"}; Configurable useGlobalTrack{"useGlobalTrack", true, "Require global track in filter"}; - Configurable cfgITScluster{"cfgITScluster", 6, "Minimum Number of ITS cluster"}; - Configurable cfgTPCcluster{"cfgTPCcluster", 80, "Minimum Number of TPC cluster"}; - Configurable cfgTPCnCrossedRows{"cfgTPCnCrossedRows", 70, "Minimum Number of TPC crossed-rows"}; - Configurable cfgTPCnCrossedRowsOverFindableCls{"cfgTPCnCrossedRowsOverFindableCls", 0.8, "Minimum ratio of crossed rows over findable clusters TPC"}; + Configurable minITSClusters{"minITSClusters", 6, "Minimum Number of ITS cluster"}; + Configurable minTPCClusters{"minTPCClusters", 80, "Minimum Number of TPC cluster"}; + Configurable minTPCCrossedRows{"minTPCCrossedRows", 70, "Minimum Number of TPC crossed-rows"}; + Configurable minTPCCrossedRowsOverFindableCls{"minTPCCrossedRowsOverFindableCls", 0.8, "Minimum ratio of crossed rows over findable clusters TPC"}; Configurable reduceOutput{"reduceOutput", 0, "Suppress info level output (0 = all output, 1 = per collision, 2 = none)"}; Filter filterTracks = (nabs(aod::track::eta) < centralEta) && (aod::track::pt >= ptMin); Filter filterCollisions = (nabs(aod::collision::posZ) < vertexXYZ.value[2]) && (nabs(aod::collision::posX) < vertexXYZ.value[0]) && (nabs(aod::collision::posY) < vertexXYZ.value[1]); Service pdg{}; // Histograms - HistogramRegistry histos1{ - "histos1", - { - {"h3DVtxZetaPhi", "VtxZ vs #eta vs #phi;VtxZ;#eta;#phi", {HistType::kTH3F, {{20, -10, 10}, {16, -0.8, 0.8}, {100, 0., o2::constants::math::TwoPI}}}}, - {"hRecoPtBefore", "Reco pT before cuts;pt (GeV/c);Counts", {HistType::kTH1F, {{1000, 0.0, 20.0}}}}, - {"hGenPtBefore", "Gen pT before cuts;pt (GeV/c);Counts", {HistType::kTH1F, {{1000, 0.0, 20.0}}}}, - {"hRecoPtAfter", "Reco pT after cuts;pt (GeV/c);Counts", {HistType::kTH1F, {{1000, 0.0, 20.0}}}}, - {"hGenPtAfter", "Gen pT after cuts;pt (GeV/c);Counts", {HistType::kTH1F, {{1000, 0.0, 20.0}}}}, - {"hRecoEtaBefore", "Reco #eta before cuts;#eta;Counts", {HistType::kTH1F, {{200, -2.0, 2.0}}}}, - {"mCentMCImpactFull", "Impact Parameter (MC All Events);Impact Parameter b [fm];Counts", {HistType::kTH1F, {{200, 0.0, 50.0}}}}, - {"mCentMCImpact005", "Impact Parameter (MC 0–5%);Impact Parameter b [fm];Counts", {HistType::kTH1F, {{200, 0.0, 20.0}}}}, - {"hGenEtaBefore", "Gen #eta before cuts;#eta;Counts", {HistType::kTH1F, {{200, -2.0, 2.0}}}}, - {"hRecoEtaAfter", "Reco #eta after cuts;#eta;Counts", {HistType::kTH1F, {{200, -2.0, 2.0}}}}, - {"hGenEtaAfter", "Gen #eta after cuts;#eta;Counts", {HistType::kTH1F, {{200, -2.0, 2.0}}}}, - {"hRecoPhiBefore", "Reco #phi before cuts;#phi;Counts", {HistType::kTH1F, {{100, 0, 6.3}}}}, - {"hGenPhiBefore", "Gen #phi before cuts;#phi;Counts", {HistType::kTH1F, {{100, 0, 6.3}}}}, - {"hRecoPhiAfter", "Reco #phi after cuts;#phi;Counts", {HistType::kTH1F, {{100, 0, 6.3}}}}, - {"hGenCharge", "Gen particle charge;Charge;Counts", {HistType::kTH1F, {{13, -6.5, 6.5}}}}, - {"hGenPhiAfter", "Gen #phi after cuts;#phi;Counts", {HistType::kTH1F, {{100, 0, 6.3}}}}, - {"mChargeDist", "MC particle charge;Charge;Counts", {HistType::kTH1F, {{9, -4.5, 4.5}}}}, - {"mDiffPhi", "dphi between selected MC particles;dphi;Counts", {HistType::kTH1F, {{2000, -0.01, 0.01}}}}, - {"mPrimariesPerEvent", "Primary MC particles per event;N_{primaries};Counts", {HistType::kTH1I, {{20000, 0, 20000}}}}, - {"mDiffPt", "dpt between selected MC particles;dpt (GeV/c);Counts", {HistType::kTH1F, {{2000, -0.01, 0.01}}}}, - {"mDiffEta", "deta between selected MC particles;deta;Counts", {HistType::kTH1F, {{2000, -0.01, 0.01}}}}, - {"mADiffPhi", "dphi between selected MC particles;dphi;Counts", {HistType::kTH1F, {{2000, -0.01, 0.01}}}}, - {"mADiffPt", "dpt between selected MC particles;dpt (GeV/c);Counts", {HistType::kTH1F, {{2000, -0.01, 0.01}}}}, - {"mADiffEta", "deta between selected MC particles;deta;Counts", {HistType::kTH1F, {{2000, -0.01, 0.01}}}}, - }, - }; - HistogramRegistry histos{ "histos", { @@ -132,7 +102,7 @@ struct FactorialMomentsTask { {"mChargeBefore", "Charge before MC cuts;charge;entries", {HistType::kTH1F, {{7, -3.5, 3.5}}}}, {"mChargeAfter", "Charge after MC cuts;charge;entries", {HistType::kTH1F, {{7, -3.5, 3.5}}}}, {"mCollID", "collisionID", {HistType::kTH1I, {{1000, -10000, 10000}}}}, - {"mCentFT0M", "centFT0C_0A", {HistType::kTH1F, {{100, 0, 100}}}}, + {"mCentFT0M", "centFT0M / Run2 V0M", {HistType::kTH1F, {{100, 0, 100}}}}, {"mCentFV0A", "centFV0A", {HistType::kTH1F, {{100, 0, 100}}}}, {"mCentFT0A", "centFT0A", {HistType::kTH1F, {{100, 0, 100}}}}, {"mCentFT0C", "centFT0C", {HistType::kTH1F, {{100, 0, 100}}}}, @@ -142,7 +112,6 @@ struct FactorialMomentsTask { {"mEta", "#eta", {HistType::kTH1F, {{1000, -2, 2}}}}, {"mPt", "#pt", {HistType::kTH1F, {{1000, -0.01, 50}}}}, {"mPhi", "#phi", {HistType::kTH1F, {{100, 0, o2::constants::math::TwoPI}}}}, - {"mEvents", "events", {HistType::kTH1D, {{5, -0.5, 4.5}}}}, {"mNFindableClsTPC", "findable TPC clusters;findable clusters", {HistType::kTH1F, {{100, 0, 200}}}}, {"mNClsTPC", "number of clusters TPC; nClusters TPC", {HistType::kTH1F, {{100, 0, 200}}}}, {"mNClsITS", "number of clusters ITS; nClusters ITS", {HistType::kTH1F, {{100, 0, 10}}}}, @@ -150,10 +119,8 @@ struct FactorialMomentsTask { {"mChi2ITS", "chi2 ITS", {HistType::kTH1F, {{100, 0, 10}}}}, {"mChi2TRD", "chi2 TRD", {HistType::kTH1F, {{100, 0, 100}}}}, {"mDCAxy", "DCA xy", {HistType::kTH1F, {{500, -0.8, 0.8}}}}, - {"mDCAx", "DCA z", {HistType::kTH1F, {{500, -2.0, 2.0}}}}, + {"mDCAz", "DCA z", {HistType::kTH1F, {{500, -2.0, 2.0}}}}, {"mDCAxyPt", "DCA xy vs #pt;#pt;DCAxy", {HistType::kTH2F, {{100, 0, 20}, {500, -0.5, 0.5}}}}, - {"mDCAxyPtbcut", "DCA xy vs #pt;#pt;DCAxycut", {HistType::kTH2F, {{100, 0, 20}, {500, -0.5, 0.5}}}}, - {"mDCAzPtbcut", "DCA z vs #pt;#pt;DCAzcut", {HistType::kTH2F, {{100, 0, 20}, {100, -2.0, 2.0}}}}, {"mDCAzPt", "DCA z vs #pt;#pt;DCAz", {HistType::kTH2F, {{100, 0, 20}, {100, -2.0, 2.0}}}}, {"mNSharedClsTPC", "shared clusters in TPC", {HistType::kTH1F, {{100, 0, 10}}}}, {"mCrossedRowsTPC", "crossedrows in TPC", {HistType::kTH1F, {{100, 0, 200}}}}, @@ -165,88 +132,110 @@ struct FactorialMomentsTask { OutputObjHandlingPolicy::AnalysisObject, true}; float collisionZ = 0.f; - double dcaxyMaxTrackPar0 = 0.0105; - double dcaxyMaxTrackPar1 = 0.035; - double dcaxyMaxTrackPar2 = 1.1; + static constexpr double kDcaXY0 = 0.0105; + static constexpr double kDcaXY1 = 0.035; + static constexpr double kDcaXY2 = 1.1; static const int nBins = 52; - double kMinCharge = 1e-6; + static constexpr double kMinCharge = 1e-6; static const int nfqOrder = 6; - int countSamples = 0; - int testc1 = 0, testc2 = 0, testc3 = 0; + static constexpr int kEvAll = 1; + static constexpr int kEvSel8 = 2; + static constexpr int kEvITSROFrame = 3; + static constexpr int kEvTimeFrame = 4; + static constexpr int kEvSameBunch = 5; + static constexpr int kEvGoodITS = 6; + static constexpr int kEvGoodZvtx = 7; + static constexpr int kEvVertexITSTPC = 8; + static constexpr int kEvCent = 9; + static constexpr int kEvAccepted = 10; + static constexpr int kNEventSteps = 10; std::array binningM{}; - std::array countTracks{0, 0, 0, 0, 0}; - std::array, 5>, 6> fqEvent{}; - std::array, 5>, 6> fqEventSampled{}; + std::array countTracks{}; + std::array, 5>, nfqOrder> fqEvent{}; std::array, 5> binConEvent{}; - std::array, 5>, 6> binConEventSampled{}; - std::array, 5>, nfqOrder> nSubsamples{}; - std::array, 5>, nfqOrder> tmpFqErr{}; std::vector> mHistArrReset; - std::vector> mHistArrQA; + std::vector> mHistEtaQA; + std::vector> mHistPtQA; + std::vector> mHistPhiQA; + std::vector> mHistMultQA; std::vector> mHistArrEff; std::vector> mFqBinFinal; std::vector> mBinConFinal; - std::vector> mFqBinFinalSampled; - std::vector> mBinConFinalSampled; - std::vector> mFqError; + // merge-safe slot output: additive per-slot sums, hadd-exact + std::shared_ptr mSlotNev; + std::vector> mSlotFq; + std::vector> mSlotOcc; + uint64_t mEventGlobalBC = 0; void init(o2::framework::InitContext&) { - // NOTE: check to make number of pt and the vector consistent - if (numPt != static_cast(ptCuts.value.size()) / 2) { - for (int i = numPt; i < static_cast(ptCuts.value.size() / 2); i++) { - ptCuts.value[2 * i] = 0; - ptCuts.value[2 * i + 1] = 0; + if (numPt.value < 1 || numPt.value > 5) { + LOG(fatal) << "numPt must be between 1 and 5, got " << numPt.value; + } + if (static_cast(ptCuts.value.size()) < 2 * numPt.value) { + LOG(fatal) << "ptCuts needs " << 2 * numPt.value << " values for numPt=" << numPt.value + << ", got " << ptCuts.value.size(); + } + if (nSlots.value < 1 || nSlots.value > 1000000) { + LOG(fatal) << "nSlots must be in [1, 1000000], got " << nSlots.value; + } + if (centLimits.value.size() != 2 || vertexXYZ.value.size() != 3) { + LOG(fatal) << "centLimits needs 2 and vertexXYZ needs 3 values"; + } + if (applyCheckPtForRec.value == applyCheckPtForMC.value) { + LOG(fatal) << "exactly one of applyCheckPtForRec / applyCheckPtForMC must be enabled"; + } + TString cfg = Form("FMtask centLimits=%d,%d numPt=%d centralEta=%g ptMin=%g dcaXY=%g dcaZ=%g useMC=%d nfqOrder=%d smearPhi=%d nSlots=%d", + centLimits.value[0], centLimits.value[1], numPt.value, centralEta.value, + ptMin.value, dcaXY.value, dcaZ.value, useMC.value ? 1 : 0, nfqOrder, smearPhi.value ? 1 : 0, nSlots.value); + cfg += Form(" vertexXYZ=%g,%g,%g", vertexXYZ.value[0], vertexXYZ.value[1], vertexXYZ.value[2]); + cfg += " ptCuts="; + for (int i = 0; i < static_cast(ptCuts.value.size()); i++) { + cfg += Form("%s%g", i > 0 ? "," : "", ptCuts.value[i]); + } + histos.add("metaConfig", cfg.Data(), HistType::kTH1D, {{1, 0.5, 1.5}}); + auto mEventSelected = std::get>(histos.add("mEventSelected", "eventSelected", HistType::kTH1D, {{kNEventSteps, 0.5, kNEventSteps + 0.5}})); + const std::array evLabels = {"all", "sel8", "kNoITSROFrameBorder", "kNoTimeFrameBorder", "kNoSameBunchPileup", "kIsGoodITSLayersAll", "kIsGoodZvtxFT0vsPV", "kIsVertexITSTPC", "centrality", "accepted"}; + for (int iStep = 0; iStep < kNEventSteps; ++iStep) { + mEventSelected->GetXaxis()->SetBinLabel(iStep + 1, evLabels[iStep]); + } + // merge-safe slot output: additive per-slot sums, hadd-exact (see slotpool.C) + const AxisSpec slotAxis{nSlots.value, -0.5, nSlots.value - 0.5, ""}; + const AxisSpec slotMAxis{nBins, -0.5, nBins - 0.5, ""}; + mSlotNev = std::get>(histos.add("slots/nEv", "events per slot;slot;events", HistType::kTH1F, {slotAxis})); + // bin1 = merged jobs, bin2 = sum of their nSlots: mismatch witness (hadd keeps the first axis silently) + auto slotParams = std::get>(histos.add("slots/params", "slot campaign parameters;jobs;entries", HistType::kTH1D, {{2, 0.5, 2.5}})); + slotParams->SetBinContent(1, 1.0); + slotParams->SetBinContent(2, static_cast(nSlots.value)); + for (int iPt = 0; iPt < numPt; ++iPt) { + mSlotOcc.push_back(std::get>(histos.add(Form("slots/sumOcc_pt%d", iPt + 1), + Form("sum of occupancy per slot for bin %.2f-%.2f;slot;M", ptCuts.value[2 * iPt], ptCuts.value[2 * iPt + 1]), + HistType::kTH2F, {slotAxis, slotMAxis}))); + for (int iq = 0; iq < nfqOrder; ++iq) { + mSlotFq.push_back(std::get>(histos.add(Form("slots/sumFq_pt%d_o%d", iPt + 1, iq + 2), + Form("sum of f_%d per slot for bin %.2f-%.2f;slot;M", iq + 2, ptCuts.value[2 * iPt], ptCuts.value[2 * iPt + 1]), + HistType::kTH2F, {slotAxis, slotMAxis}))); } } - std::array axisPt = {AxisSpec{100, -0.01, 3 * ptCuts.value[1], ""}, AxisSpec{100, -0.01, 3 * ptCuts.value[3], ""}, AxisSpec{100, -0.01, 3 * ptCuts.value[5], ""}, AxisSpec{100, -0.01, 3 * ptCuts.value[7], ""}, AxisSpec{100, -0.01, 3 * ptCuts.value[9], ""}}; // pT axis - auto mEventSelected = std::get>(histos.add("mEventSelected", "eventSelected", HistType::kTH1D, {{7, 0.5, 8.5}})); - mEventSelected->GetXaxis()->SetBinLabel(0, "all"); - mEventSelected->GetXaxis()->SetBinLabel(1, "sel8"); - mEventSelected->GetXaxis()->SetBinLabel(2, "kNoITSROFrameBorder"); - mEventSelected->GetXaxis()->SetBinLabel(3, "kNoTimeFrameBorder"); - mEventSelected->GetXaxis()->SetBinLabel(4, "sameBunchPileup"); - mEventSelected->GetXaxis()->SetBinLabel(5, "kIsGoodITSLayersAll"); - mEventSelected->GetXaxis()->SetBinLabel(6, "kIsGoodZvtxFT0vsPV"); - mEventSelected->GetXaxis()->SetBinLabel(7, "FTOC"); - auto mTrackSelected = std::get>(histos.add( - "mTrackSelected", "Track Selection Steps", HistType::kTH1D, {{5, 0.5, 5.5}})); - mTrackSelected->GetXaxis()->SetBinLabel(1, "all"); - mTrackSelected->GetXaxis()->SetBinLabel(2, "charge"); - mTrackSelected->GetXaxis()->SetBinLabel(3, "PIDs"); - mTrackSelected->GetXaxis()->SetBinLabel(4, "Pids exclude pions"); - mTrackSelected->GetXaxis()->SetBinLabel(5, "Pids exclude pions + kions"); - if (useMC) { - auto mMcTrackSelected = std::get>(histos.add("mMcTrackSelected", "mcTrackSelection", HistType::kTH1D, {{5, 0.5, 5.5}})); - } for (int iM = 0; iM < nBins; ++iM) { binningM[iM] = 2 * (iM + 2); } for (int iPt = 0; iPt < numPt; ++iPt) { + const AxisSpec ptAxis{100, -0.01, 3 * ptCuts.value[2 * iPt + 1], ""}; mHistArrEff.push_back(std::get>(histos.add(Form("bin%i/m3DVtxZetaPhi", iPt + 1), Form("#eta #phi #vtxz for bin %.2f-%.2f;vz;#eta;#phi", ptCuts.value[2 * iPt], ptCuts.value[2 * iPt + 1]), HistType::kTH3F, {{20, -10, 10}, {16, -0.8, +0.8}, {100, 0., o2::constants::math::TwoPI}}))); - mHistArrQA.push_back(std::get>(histos.add(Form("bin%i/mEta", iPt + 1), Form("#eta for bin %.2f-%.2f;#eta", ptCuts.value[2 * iPt], ptCuts.value[2 * iPt + 1]), HistType::kTH1F, {{1000, -2, 2}}))); - mHistArrQA.push_back(std::get>(histos.add(Form("bin%i/mPt", iPt + 1), Form("pT for bin %.2f-%.2f;pT", ptCuts.value[2 * iPt], ptCuts.value[2 * iPt + 1]), HistType::kTH1F, {axisPt[iPt]}))); - mHistArrQA.push_back(std::get>(histos.add(Form("bin%i/mPhi", iPt + 1), Form("#phi for bin %.2f-%.2f;#phi", ptCuts.value[2 * iPt], ptCuts.value[2 * iPt + 1]), HistType::kTH1F, {{1000, 0, o2::constants::math::TwoPI}}))); - mHistArrQA.push_back(std::get>(histos.add(Form("bin%i/mMultiplicity", iPt + 1), Form("Multiplicity for bin %.2f-%.2f;Multiplicity", ptCuts.value[2 * iPt], ptCuts.value[2 * iPt + 1]), HistType::kTH1F, {{1000, 0, 15000}}))); + mHistEtaQA.push_back(std::get>(histos.add(Form("bin%i/mEta", iPt + 1), Form("#eta for bin %.2f-%.2f;#eta", ptCuts.value[2 * iPt], ptCuts.value[2 * iPt + 1]), HistType::kTH1F, {{1000, -2, 2}}))); + mHistPtQA.push_back(std::get>(histos.add(Form("bin%i/mPt", iPt + 1), Form("pT for bin %.2f-%.2f;pT", ptCuts.value[2 * iPt], ptCuts.value[2 * iPt + 1]), HistType::kTH1F, {ptAxis}))); + mHistPhiQA.push_back(std::get>(histos.add(Form("bin%i/mPhi", iPt + 1), Form("#phi for bin %.2f-%.2f;#phi", ptCuts.value[2 * iPt], ptCuts.value[2 * iPt + 1]), HistType::kTH1F, {{1000, 0, o2::constants::math::TwoPI}}))); + mHistMultQA.push_back(std::get>(histos.add(Form("bin%i/mMultiplicity", iPt + 1), Form("Multiplicity for bin %.2f-%.2f;Multiplicity", ptCuts.value[2 * iPt], ptCuts.value[2 * iPt + 1]), HistType::kTH1F, {{1000, 0, 15000}}))); for (int iM = 0; iM < nBins; ++iM) { auto mHistsR = std::get>(histos.add(Form("bin%i/Reset/mEtaPhi%i", iPt + 1, iM), Form("#eta#phi_%i for bin %.2f-%.2f;#eta;#phi", iM, ptCuts.value[2 * iPt], ptCuts.value[2 * iPt + 1]), HistType::kTH2F, {{binningM[iM], -0.8, 0.8}, {binningM[iM], 0, o2::constants::math::TwoPI}})); mHistArrReset.push_back(mHistsR); - for (int iq = 0; iq < nfqOrder; ++iq) { - tmpFqErr[iq][iPt][iM] = new TH1D(Form("tmpFqErr%i%i%i", iq, iPt, iM), Form("tmpFqErr%i%i%i", iq, iPt, iM), 100, 0, 10); - } } for (int i = 0; i < nfqOrder; ++i) { auto mHistFq = std::get>(histos.add(Form("mFinalFq%i_bin%i", i + 2, iPt + 1), Form("Final F_%i for bin %.2f-%.2f;M", i + 2, ptCuts.value[2 * iPt], ptCuts.value[2 * iPt + 1]), HistType::kTH1F, {{nBins, -0.5, nBins - 0.5}})); mFqBinFinal.push_back(mHistFq); auto mHistAv = std::get>(histos.add(Form("mFinalAvBin%i_bin%i", i + 2, iPt + 1), Form("Final AvBin_%i for bin %.2f-%.2f;M", i + 2, ptCuts.value[2 * iPt], ptCuts.value[2 * iPt + 1]), HistType::kTH1F, {{nBins, -0.5, nBins - 0.5}})); mBinConFinal.push_back(mHistAv); - auto mHistFqSampled = std::get>(histos.add(Form("mFinalFq%iSampled_bin%i", i + 2, iPt + 1), Form("Final F_%i for bin %.2f-%.2f;M", i + 2, ptCuts.value[2 * iPt], ptCuts.value[2 * iPt + 1]), HistType::kTH1F, {{nBins, -0.5, nBins - 0.5}})); - mFqBinFinalSampled.push_back(mHistFqSampled); - auto mHistAvSampled = std::get>(histos.add(Form("mFinalAvBin%iSampled_bin%i", i + 2, iPt + 1), Form("Final AvBin_%i for bin %.2f-%.2f;M", i + 2, ptCuts.value[2 * iPt], ptCuts.value[2 * iPt + 1]), HistType::kTH1F, {{nBins, -0.5, nBins - 0.5}})); - mBinConFinalSampled.push_back(mHistAvSampled); - - auto mHistError = std::get>(histos.add(Form("mFqError%i_bin%i", i + 2, iPt + 1), Form("Error for F_%i for bin %.2f-%.2f;M", i + 2, ptCuts.value[2 * iPt], ptCuts.value[2 * iPt + 1]), HistType::kTH1F, {{nBins, -0.5, nBins - 0.5}})); - mFqError.push_back(mHistError); } } } @@ -256,12 +245,14 @@ struct FactorialMomentsTask { for (auto iPt = 0; iPt < numPt; ++iPt) { if (track.pt() > ptCuts.value[2 * iPt] && track.pt() < ptCuts.value[2 * iPt + 1]) { float iphi = track.phi(); - iphi = gRandom->Gaus(iphi, o2::constants::math::TwoPI); - iphi = RecoDecay::constrainAngle(iphi, 0.); + if (smearPhi) { + iphi = gRandom->Gaus(iphi, o2::constants::math::TwoPI); + iphi = RecoDecay::constrainAngle(iphi, 0.); + } mHistArrEff[iPt]->Fill(collisionZ, track.eta(), iphi); - mHistArrQA[iPt * 4]->Fill(track.eta()); - mHistArrQA[iPt * 4 + 1]->Fill(track.pt()); - mHistArrQA[iPt * 4 + 2]->Fill(iphi); + mHistEtaQA[iPt]->Fill(track.eta()); + mHistPtQA[iPt]->Fill(track.pt()); + mHistPhiQA[iPt]->Fill(iphi); countTracks[iPt]++; for (auto iM = 0; iM < nBins; ++iM) { mHistArrReset[iPt * nBins + iM]->Fill(track.eta(), iphi); @@ -269,106 +260,224 @@ struct FactorialMomentsTask { } } } - void calculateMoments(std::vector> hist) + // splitmix64: deterministic stateless hash of the global BC + static uint64_t slotOf(uint64_t x) + { + x += 0x9e3779b97f4a7c15ULL; + x = (x ^ (x >> 30)) * 0xbf58476d1ce4e5b9ULL; + x = (x ^ (x >> 27)) * 0x94d049bb133111ebULL; + return x ^ (x >> 31); + } + // additive fold of the current accepted event into its slot (never Fill with weights: no fSumw2, hadd-exact) + void accumulateSlot() + { + const int slot = static_cast(slotOf(mEventGlobalBC) % static_cast(nSlots.value)); + mSlotNev->SetBinContent(slot + 1, mSlotNev->GetBinContent(slot + 1) + 1.); + for (int iPt = 0; iPt < numPt; ++iPt) { + auto* occ = mSlotOcc[iPt].get(); + for (int iM = 0; iM < nBins; ++iM) { + occ->SetBinContent(slot + 1, iM + 1, occ->GetBinContent(slot + 1, iM + 1) + binConEvent[iPt][iM]); + for (int iq = 0; iq < nfqOrder; ++iq) { + auto* fq = mSlotFq[iPt * nfqOrder + iq].get(); + fq->SetBinContent(slot + 1, iM + 1, fq->GetBinContent(slot + 1, iM + 1) + fqEvent[iq][iPt][iM]); + } + } + } + } + void calculateMoments(const std::vector>& hist) { - double binContent = 0; - countSamples++; - bool compSample = kFALSE; - if (countSamples == samplesize) { - compSample = kTRUE; - countSamples = 0; - } - // Calculate the normalized factorial moments + // normalized factorial moments per accepted event, per pT bin and M for (int iPt = 0; iPt < numPt; ++iPt) { for (int iM = 0; iM < nBins; ++iM) { - binContent = 0; - std::array sumfqBin{0.}; + double binContent = 0; + std::array sumfqBin{0.}; + const double m2 = static_cast(binningM[iM]) * static_cast(binningM[iM]); + const auto& histM = hist[iPt * nBins + iM]; - for (int iEta = 1; iEta <= hist[iPt * nBins + iM]->GetNbinsX(); ++iEta) { - for (int iPhi = 1; iPhi <= hist[iPt * nBins + iM]->GetNbinsY(); ++iPhi) { - double binconVal = 0; - binconVal = hist[iPt * nBins + iM]->GetBinContent(iEta, iPhi); + for (int iEta = 1; iEta <= histM->GetNbinsX(); ++iEta) { + for (int iPhi = 1; iPhi <= histM->GetNbinsY(); ++iPhi) { + const double binconVal = histM->GetBinContent(iEta, iPhi); binContent += binconVal; for (int iq = 0; iq < nfqOrder; ++iq) { double fqBin = 0; if (binconVal >= iq + 2) { - fqBin = TMath::Factorial(static_cast(binconVal)) / (TMath::Factorial(static_cast(binconVal) - (iq + 2))); - } - if (std::isnan(fqBin)) { - break; + const int nOcc = static_cast(binconVal); + fqBin = 1.0; + for (int k = 0; k < iq + 2; ++k) { + fqBin *= (nOcc - k); + } } sumfqBin[iq] += fqBin; } } } - binConEvent[iPt][iM] = binContent / (std::pow(binningM[iM], 2)); + binConEvent[iPt][iM] = binContent / m2; for (int iq = 0; iq < nfqOrder; ++iq) { - // Always accumulate raw fq numerator and binCon unconditionally so - // both sums run over exactly the same set of events. - fqEvent[iq][iPt][iM] = sumfqBin[iq] / (std::pow(binningM[iM], 2)); - fqEventSampled[iq][iPt][iM] += fqEvent[iq][iPt][iM]; - binConEventSampled[iq][iPt][iM] += binConEvent[iPt][iM]; - mFqBinFinal[iPt * 6 + iq]->Fill(iM, fqEvent[iq][iPt][iM]); - mBinConFinal[iPt * 6 + iq]->Fill(iM, binConEvent[iPt][iM]); - if (compSample) { - const double avM = binConEventSampled[iq][iPt][iM] / samplesize; - mBinConFinalSampled[iPt * 6 + iq]->Fill(iM, avM); - const double den = std::pow(avM, iq + 2); - const double fqM = ((fqEventSampled[iq][iPt][iM]) / samplesize); + fqEvent[iq][iPt][iM] = sumfqBin[iq] / m2; + mFqBinFinal[iPt * nfqOrder + iq]->Fill(iM, fqEvent[iq][iPt][iM]); + mBinConFinal[iPt * nfqOrder + iq]->Fill(iM, binConEvent[iPt][iM]); + } + } + } + accumulateSlot(); + } - const double tmp = (den > 0.) ? fqM / den : 0.; - mFqBinFinalSampled[iPt * 6 + iq]->Fill(iM, tmp); + template + bool passDcaCut(const T& track) const + { + return std::fabs(track.dcaXY()) < (kDcaXY0 + kDcaXY1 / std::pow(track.pt(), kDcaXY2)); + } - tmpFqErr[iq][iPt][iM]->Fill(tmp); - ++nSubsamples[iq][iPt][iM]; - const double stdErr = (nSubsamples[iq][iPt][iM] > 1) - ? tmpFqErr[iq][iPt][iM]->GetStdDev() / std::sqrt(static_cast(nSubsamples[iq][iPt][iM])) - : 0.; - mFqError[iPt * 6 + iq]->SetBinContent(iM + 1, stdErr); - fqEventSampled[iq][iPt][iM] = 0; - binConEventSampled[iq][iPt][iM] = 0; - } - } // end of loop over order of F - } // end of loop over M bins - } // end of loop over pT bins + template + void fillTrackQA(const T& track) + { + histos.fill(HIST("mCollID"), track.collisionId()); + histos.fill(HIST("mEta"), track.eta()); + histos.fill(HIST("mPt"), track.pt()); + histos.fill(HIST("mPhi"), track.phi()); + histos.fill(HIST("mNFindableClsTPC"), track.tpcNClsFindable()); + histos.fill(HIST("mNClsTPC"), track.tpcNClsFound()); + histos.fill(HIST("mNClsITS"), track.itsNCls()); + histos.fill(HIST("mChi2TPC"), track.tpcChi2NCl()); + histos.fill(HIST("mChi2ITS"), track.itsChi2NCl()); + histos.fill(HIST("mChi2TRD"), track.trdChi2()); + histos.fill(HIST("mDCAxy"), track.dcaXY()); + histos.fill(HIST("mDCAz"), track.dcaZ()); + histos.fill(HIST("mtpcsignalvspt"), track.pt(), track.tpcSignal()); + histos.fill(HIST("mDCAxyPt"), track.pt(), track.dcaXY()); + histos.fill(HIST("mDCAzPt"), track.pt(), track.dcaZ()); + histos.fill(HIST("mNSharedClsTPC"), track.tpcNClsShared()); + histos.fill(HIST("mCrossedRowsTPC"), track.tpcNClsCrossedRows()); + histos.fill(HIST("mNFinClsminusCRows"), track.tpcNClsFindableMinusCrossedRows()); + histos.fill(HIST("mNFractionShClsTPC"), track.tpcFractionSharedCls()); + histos.fill(HIST("mSharedClsvsPt"), track.pt(), track.tpcNClsShared()); + const double crossedRows = track.tpcNClsCrossedRows(); + histos.fill(HIST("mSharedClsProbvsPt"), track.pt(), crossedRows > 0. ? track.tpcFractionSharedCls() / crossedRows : 0.); } - // write a template that takes coll and tracks from processRun3 and processMCRec - using TracksFMs = soa::Filtered>; - Preslice perMcCollision = aod::mcparticle::mcCollisionId; - void processRun3(soa::Filtered>::iterator const& coll, TracksFMs const& tracks) + template + void fillMcQA(const T& mcParts) + { + for (auto const& mc : mcParts) { + int pdgCode = mc.pdgCode(); + auto pdgInfo = pdg->GetParticle(pdgCode); + if (!pdgInfo) { + continue; + } + double charge = pdgInfo->Charge(); + double physCharge = charge / 3.0; + histos.fill(HIST("mChargeBefore"), physCharge); + if (mc.isPhysicalPrimary() && std::abs(mc.eta()) < centralEta && std::abs(physCharge) >= kMinCharge) { + histos.fill(HIST("mChargeAfter"), physCharge); + histos.fill(HIST("mEta"), mc.eta()); + histos.fill(HIST("mPt"), mc.pt()); + histos.fill(HIST("mPhi"), mc.phi()); + if (applyCheckPtForMC && !applyCheckPtForRec) { + checkpT(mc); + } + } + } + } + + template + bool selectEventRun3(const T& coll) { - // selection of events - histos.fill(HIST("mEventSelected"), 0); + histos.fill(HIST("mEventSelected"), kEvAll); if (!coll.sel8()) { - return; + return false; } - if (cfgEvSelkNoITSROFrameBorder && !(coll.selection_bit(o2::aod::evsel::kNoITSROFrameBorder))) { - return; + histos.fill(HIST("mEventSelected"), kEvSel8); + if (evSelNoITSROFrameBorder && !(coll.selection_bit(o2::aod::evsel::kNoITSROFrameBorder))) { + return false; + } + histos.fill(HIST("mEventSelected"), kEvITSROFrame); + if (evSelNoTimeFrameBorder && !(coll.selection_bit(o2::aod::evsel::kNoTimeFrameBorder))) { + return false; + } + histos.fill(HIST("mEventSelected"), kEvTimeFrame); + if (evSelNoSameBunchPileup && !coll.selection_bit(o2::aod::evsel::kNoSameBunchPileup)) { + return false; + } + histos.fill(HIST("mEventSelected"), kEvSameBunch); + if (evSelGoodITSLayersAll && !(coll.selection_bit(o2::aod::evsel::kIsGoodITSLayersAll))) { + return false; + } + histos.fill(HIST("mEventSelected"), kEvGoodITS); + if (evSelGoodZvtxFT0vsPV && !coll.selection_bit(o2::aod::evsel::kIsGoodZvtxFT0vsPV)) { + return false; + } + histos.fill(HIST("mEventSelected"), kEvGoodZvtx); + if (coll.centFT0C() < centLimits.value[0] || coll.centFT0C() > centLimits.value[1]) { + return false; } + collisionZ = coll.posZ(); + histos.fill(HIST("mEventSelected"), kEvCent); + histos.fill(HIST("mEventSelected"), kEvAccepted); + return true; + } - histos.fill(HIST("mEventSelected"), 2); - if (cfgEvSelkNoTimeFrameBorder && !(coll.selection_bit(o2::aod::evsel::kNoTimeFrameBorder))) { - return; + template + bool selectEventMCRec(const T& coll) + { + histos.fill(HIST("mEventSelected"), kEvAll); + if (!coll.sel8()) { + return false; } - histos.fill(HIST("mEventSelected"), 3); - if (isApplySameBunchPileup && !coll.selection_bit(o2::aod::evsel::kNoSameBunchPileup)) { - return; + histos.fill(HIST("mEventSelected"), kEvSel8); + if (evSelNoSameBunchPileup && !coll.selection_bit(o2::aod::evsel::kNoSameBunchPileup)) { + return false; } - histos.fill(HIST("mEventSelected"), 4); - if (cfgUseGoodITSLayerAllCut && !(coll.selection_bit(o2::aod::evsel::kIsGoodITSLayersAll))) { - return; + histos.fill(HIST("mEventSelected"), kEvSameBunch); + if (evSelGoodZvtxFT0vsPV && !coll.selection_bit(o2::aod::evsel::kIsGoodZvtxFT0vsPV)) { + return false; } - histos.fill(HIST("mEventSelected"), 5); - if (isApplyGoodZvtxFT0vsPV && !coll.selection_bit(o2::aod::evsel::kIsGoodZvtxFT0vsPV)) { - return; + histos.fill(HIST("mEventSelected"), kEvGoodZvtx); + if (evSelVertexITSTPC && !coll.selection_bit(o2::aod::evsel::kIsVertexITSTPC)) { + return false; } - histos.fill(HIST("mEventSelected"), 6); + histos.fill(HIST("mEventSelected"), kEvVertexITSTPC); if (coll.centFT0C() < centLimits.value[0] || coll.centFT0C() > centLimits.value[1]) { - return; + return false; } collisionZ = coll.posZ(); - histos.fill(HIST("mEventSelected"), 7); + histos.fill(HIST("mEventSelected"), kEvCent); + histos.fill(HIST("mEventSelected"), kEvAccepted); + return true; + } + + void beginEvent() + { + for (int iPt = 0; iPt < numPt; ++iPt) { + if (countTracks[iPt] == 0) { + continue; + } + for (int iM = 0; iM < nBins; ++iM) { + mHistArrReset[iPt * nBins + iM]->Reset(); + } + } + countTracks.fill(0); + fqEvent = {}; + binConEvent = {}; + } + + void endEvent() + { + for (int iPt = 0; iPt < numPt; ++iPt) { + if (countTracks[iPt] > 0) { + mHistMultQA[iPt]->Fill(countTracks[iPt]); + } + } + calculateMoments(mHistArrReset); + } + + using TracksFMs = soa::Filtered>; + Preslice perMcCollision = aod::mcparticle::mcCollisionId; + void processRun3(soa::Filtered>::iterator const& coll, TracksFMs const& tracks, aod::BCs const&) + { + mEventGlobalBC = coll.bc_as().globalBC(); + if (!selectEventRun3(coll)) { + return; + } histos.fill(HIST("mVertexX"), coll.posX()); histos.fill(HIST("mVertexY"), coll.posY()); histos.fill(HIST("mVertexZ"), coll.posZ()); @@ -376,12 +485,7 @@ struct FactorialMomentsTask { histos.fill(HIST("mCentFV0A"), coll.centFV0A()); histos.fill(HIST("mCentFT0A"), coll.centFT0A()); histos.fill(HIST("mCentFT0C"), coll.centFT0C()); - for (auto const& h : mHistArrReset) { - h->Reset(); - } - countTracks = {0, 0, 0, 0, 0}; - fqEvent = {{{{{0, 0, 0, 0, 0, 0}}}}}; - binConEvent = {{{0, 0, 0, 0, 0}}}; + beginEvent(); for (auto const& track : tracks) { if (useITS && !track.hasITS()) { continue; @@ -392,79 +496,31 @@ struct FactorialMomentsTask { if (useGlobal && !track.isGlobalTrack()) { continue; } - - if (std::fabs(track.dcaXY()) < (dcaxyMaxTrackPar0 + dcaxyMaxTrackPar1 / std::pow(track.pt(), dcaxyMaxTrackPar2))) { - histos.fill(HIST("mCollID"), track.collisionId()); - histos.fill(HIST("mEta"), track.eta()); - histos.fill(HIST("mPt"), track.pt()); - histos.fill(HIST("mPhi"), track.phi()); - histos.fill(HIST("mNFindableClsTPC"), track.tpcNClsFindable()); - histos.fill(HIST("mNClsTPC"), track.tpcNClsFound()); - histos.fill(HIST("mNClsITS"), track.itsNCls()); - histos.fill(HIST("mChi2TPC"), track.tpcChi2NCl()); - histos.fill(HIST("mChi2ITS"), track.itsChi2NCl()); - histos.fill(HIST("mChi2TRD"), track.trdChi2()); - histos.fill(HIST("mDCAxy"), track.dcaXY()); - histos.fill(HIST("mtpcsignalvspt"), track.pt(), track.tpcSignal()); - histos.fill(HIST("mDCAxyPt"), track.pt(), track.dcaXY()); - histos.fill(HIST("mDCAzPt"), track.pt(), track.dcaZ()); - histos.fill(HIST("mDCAzPt"), track.pt(), track.dcaZ()); - histos.fill(HIST("mNSharedClsTPC"), track.tpcNClsShared()); - histos.fill(HIST("mCrossedRowsTPC"), track.tpcNClsCrossedRows()); - histos.fill(HIST("mNFinClsminusCRows"), track.tpcNClsFindableMinusCrossedRows()); - histos.fill(HIST("mNFractionShClsTPC"), track.tpcFractionSharedCls()); - histos.fill(HIST("mSharedClsvsPt"), track.pt(), track.tpcNClsShared()); - histos.fill(HIST("mSharedClsProbvsPt"), track.pt(), track.tpcFractionSharedCls() / track.tpcNClsCrossedRows()); - checkpT(track); - } - } - for (int iPt = 0; iPt < numPt; ++iPt) { - if (countTracks[iPt] > 0) { - mHistArrQA[iPt * 4 + 3]->Fill(countTracks[iPt]); + if (!passDcaCut(track)) { + continue; } + fillTrackQA(track); + checkpT(track); } - calculateMoments(mHistArrReset); + endEvent(); } PROCESS_SWITCH(FactorialMomentsTask, processRun3, "main process function", false); using CollisionCandidateMCRec = soa::Join; using TracksMc = soa::Filtered>; - void processMCRec(soa::Filtered::iterator const& coll, TracksMc const& colltracks, aod::McParticles const& mcParticles, aod::McCollisions const&) + void processMCRec(soa::Filtered::iterator const& coll, TracksMc const& colltracks, aod::McParticles const& mcParticles, aod::McCollisions const&, aod::BCs const&) { + mEventGlobalBC = coll.bc_as().globalBC(); if (!coll.has_mcCollision()) { return; } - histos.fill(HIST("mEventSelected"), 0); - if (!coll.sel8()) { - return; - } - histos.fill(HIST("mEventSelected"), 1); - if (isApplySameBunchPileup && !coll.selection_bit(o2::aod::evsel::kNoSameBunchPileup)) { - return; - } - histos.fill(HIST("mEventSelected"), 2); - if (isApplyGoodZvtxFT0vsPV && !coll.selection_bit(o2::aod::evsel::kIsGoodZvtxFT0vsPV)) { + if (!selectEventMCRec(coll)) { return; } - histos.fill(HIST("mEventSelected"), 3); - if (isApplyVertexITSTPC && !coll.selection_bit(o2::aod::evsel::kIsVertexITSTPC)) { - return; - } - histos.fill(HIST("mEventSelected"), 4); - if (coll.centFT0C() < centLimits.value[0] || coll.centFT0C() > centLimits.value[1]) { - return; - } - - histos.fill(HIST("mEventSelected"), 5); histos.fill(HIST("mVertexX"), coll.posX()); histos.fill(HIST("mVertexY"), coll.posY()); histos.fill(HIST("mVertexZ"), coll.posZ()); histos.fill(HIST("mCentFT0C"), coll.centFT0C()); - for (auto const& h : mHistArrReset) { - h->Reset(); - } - countTracks = {0, 0, 0, 0, 0}; - fqEvent = {{{{{0, 0, 0, 0, 0, 0}}}}}; - binConEvent = {{{0, 0, 0, 0, 0}}}; + beginEvent(); for (auto const& track : colltracks) { if (useITS && !track.hasITS()) { continue; @@ -475,60 +531,17 @@ struct FactorialMomentsTask { if (useGlobal && !track.isGlobalTrack()) { continue; } - if (std::fabs(track.dcaXY()) < (dcaxyMaxTrackPar0 + dcaxyMaxTrackPar1 / std::pow(track.pt(), dcaxyMaxTrackPar2))) { - histos.fill(HIST("mCollID"), track.collisionId()); - histos.fill(HIST("mEta"), track.eta()); - histos.fill(HIST("mPt"), track.pt()); - histos.fill(HIST("mPhi"), track.phi()); - histos.fill(HIST("mNFindableClsTPC"), track.tpcNClsFindable()); - histos.fill(HIST("mNClsTPC"), track.tpcNClsFound()); - histos.fill(HIST("mNClsITS"), track.itsNCls()); - histos.fill(HIST("mChi2TPC"), track.tpcChi2NCl()); - histos.fill(HIST("mChi2ITS"), track.itsChi2NCl()); - histos.fill(HIST("mChi2TRD"), track.trdChi2()); - histos.fill(HIST("mDCAxy"), track.dcaXY()); - histos.fill(HIST("mDCAx"), track.dcaZ()); - histos.fill(HIST("mtpcsignalvspt"), track.pt(), track.tpcSignal()); - histos.fill(HIST("mDCAxyPt"), track.pt(), track.dcaXY()); - histos.fill(HIST("mDCAzPt"), track.pt(), track.dcaZ()); - histos.fill(HIST("mNSharedClsTPC"), track.tpcNClsShared()); - histos.fill(HIST("mCrossedRowsTPC"), track.tpcNClsCrossedRows()); - histos.fill(HIST("mNFinClsminusCRows"), track.tpcNClsFindableMinusCrossedRows()); - histos.fill(HIST("mNFractionShClsTPC"), track.tpcFractionSharedCls()); - histos.fill(HIST("mSharedClsvsPt"), track.pt(), track.tpcNClsShared()); - histos.fill(HIST("mSharedClsProbvsPt"), track.pt(), track.tpcFractionSharedCls() / track.tpcNClsCrossedRows()); - if (applyCheckPtForRec && !applyCheckPtForMC) { - checkpT(track); - } - } - } - auto mcParts = mcParticles.sliceBy(perMcCollision, coll.mcCollision().globalIndex()); - for (auto const& mc : mcParts) { - int pdgCode = mc.pdgCode(); - auto pdgInfo = pdg->GetParticle(pdgCode); - if (!pdgInfo) { + if (!passDcaCut(track)) { continue; } - double charge = pdgInfo->Charge(); - double physCharge = charge / 3.0; - histos.fill(HIST("mChargeBefore"), physCharge); - if (mc.isPhysicalPrimary() && std::abs(mc.eta()) < centralEta && std::abs(physCharge) >= kMinCharge) { - histos.fill(HIST("mChargeAfter"), physCharge); - histos.fill(HIST("mEta"), mc.eta()); - histos.fill(HIST("mPt"), mc.pt()); - histos.fill(HIST("mPhi"), mc.phi()); - if (applyCheckPtForMC && !applyCheckPtForRec) { - checkpT(mc); - } - } - } - for (auto iPt = 0; iPt < numPt; ++iPt) { - if (countTracks[iPt] > 0) { - mHistArrQA[iPt * 4 + 3]->Fill(countTracks[iPt]); + fillTrackQA(track); + if (applyCheckPtForRec && !applyCheckPtForMC) { + checkpT(track); } } - histos.fill(HIST("mEventSelected"), 6); - calculateMoments(mHistArrReset); + auto mcParts = mcParticles.sliceBy(perMcCollision, coll.mcCollision().globalIndex()); + fillMcQA(mcParts); + endEvent(); } PROCESS_SWITCH(FactorialMomentsTask, processMCRec, "main process function", false); using EventSelectionrun2 = soa::Join; @@ -543,129 +556,61 @@ struct FactorialMomentsTask { BCsWithRun2Info const&) { auto bc = coll.bc_as(); + mEventGlobalBC = bc.globalBC(); if (!(static_cast(bc.eventCuts() & BIT(aod::Run2EventCuts::kAliEventCutsAccepted)))) { return; } if (coll.centRun2V0M() < centLimits.value[0] || coll.centRun2V0M() > centLimits.value[1]) { return; } + collisionZ = coll.posZ(); histos.fill(HIST("mVertexX"), coll.posX()); histos.fill(HIST("mVertexY"), coll.posY()); histos.fill(HIST("mVertexZ"), coll.posZ()); histos.fill(HIST("mCentFT0M"), coll.centRun2V0M()); - for (auto const& h : mHistArrReset) { - h->Reset(); - } - countTracks = {0, 0, 0, 0, 0}; - fqEvent = {{{{{0, 0, 0, 0, 0, 0}}}}}; - binConEvent = {{{0, 0, 0, 0, 0}}}; + beginEvent(); for (auto const& track : tracks) { double recoCharge = (track.sign() != 0) ? track.sign() : 0.; if (std::abs(track.eta()) < centralEta && track.isGlobalTrack() && std::abs(recoCharge) >= kMinCharge) { - histos.fill(HIST("mCollID"), track.collisionId()); - histos.fill(HIST("mNFindableClsTPC"), track.tpcNClsFindable()); - histos.fill(HIST("mNClsTPC"), track.tpcNClsFound()); - histos.fill(HIST("mNClsITS"), track.itsNCls()); - histos.fill(HIST("mChi2TPC"), track.tpcChi2NCl()); - histos.fill(HIST("mChi2ITS"), track.itsChi2NCl()); - histos.fill(HIST("mChi2TRD"), track.trdChi2()); - histos.fill(HIST("mNSharedClsTPC"), track.tpcNClsShared()); - histos.fill(HIST("mCrossedRowsTPC"), track.tpcNClsCrossedRows()); - histos.fill(HIST("mNFinClsminusCRows"), track.tpcNClsFindableMinusCrossedRows()); - histos.fill(HIST("mNFractionShClsTPC"), track.tpcFractionSharedCls()); - histos.fill(HIST("mSharedClsvsPt"), track.pt(), track.tpcNClsShared()); - histos.fill(HIST("mSharedClsProbvsPt"), track.pt(), track.tpcFractionSharedCls() / track.tpcNClsCrossedRows()); + fillTrackQA(track); if (applyCheckPtForRec && !applyCheckPtForMC) { checkpT(track); } } } auto mcParts = mcParticles.sliceBy(perMcCollision, coll.mcCollision().globalIndex()); - for (auto const& mc : mcParts) { - int pdgCode = mc.pdgCode(); - auto pdgInfo = pdg->GetParticle(pdgCode); - if (!pdgInfo) { - continue; - } - double charge = pdgInfo->Charge(); - double physCharge = charge / 3.0; - histos.fill(HIST("mChargeBefore"), physCharge); - if (mc.isPhysicalPrimary() && std::abs(mc.eta()) < centralEta && std::abs(physCharge) >= kMinCharge) { - histos.fill(HIST("mChargeAfter"), physCharge); - histos.fill(HIST("mEta"), mc.eta()); - histos.fill(HIST("mPt"), mc.pt()); - histos.fill(HIST("mPhi"), mc.phi()); - if (applyCheckPtForMC && !applyCheckPtForRec) { - checkpT(mc); - } - } - } - - for (int iPt = 0; iPt < numPt; ++iPt) { - if (countTracks[iPt] > 0) { - mHistArrQA[iPt * 4 + 3]->Fill(countTracks[iPt]); - } - } - - calculateMoments(mHistArrReset); + fillMcQA(mcParts); + endEvent(); } - PROCESS_SWITCH(FactorialMomentsTask, processMcRun2, "process MC Run2", true); - void processRun2(soa::Filtered>::iterator const& coll, TracksFMs const& tracks) + PROCESS_SWITCH(FactorialMomentsTask, processMcRun2, "process MC Run2", false); + void processRun2(soa::Filtered>::iterator const& coll, TracksFMs const& tracks, aod::BCs const&) { + mEventGlobalBC = coll.bc_as().globalBC(); if ((!coll.alias_bit(kINT7)) || (!coll.sel7())) { return; } if (coll.centRun2V0M() < centLimits.value[0] || coll.centRun2V0M() > centLimits.value[1]) { return; } + collisionZ = coll.posZ(); histos.fill(HIST("mVertexX"), coll.posX()); histos.fill(HIST("mVertexY"), coll.posY()); histos.fill(HIST("mVertexZ"), coll.posZ()); histos.fill(HIST("mCentFT0M"), coll.centRun2V0M()); - - for (auto const& h : mHistArrReset) { - h->Reset(); - } - countTracks = {0, 0, 0, 0, 0}; - fqEvent = {{{{{0, 0, 0, 0, 0, 0}}}}}; - binConEvent = {{{0, 0, 0, 0, 0}}}; + beginEvent(); for (auto const& track : tracks) { - if ((track.pt() < ptMin) || (!track.isGlobalTrack()) || (track.tpcNClsFindable() < mintPCCls)) { + if ((track.pt() < ptMin) || (!track.isGlobalTrack()) || (track.tpcNClsFindable() < minTPCNCls)) { continue; } - histos.fill(HIST("mCollID"), track.collisionId()); - histos.fill(HIST("mEta"), track.eta()); - histos.fill(HIST("mPt"), track.pt()); - histos.fill(HIST("mPhi"), track.phi()); - histos.fill(HIST("mNFindableClsTPC"), track.tpcNClsFindable()); - histos.fill(HIST("mNClsTPC"), track.tpcNClsFound()); - histos.fill(HIST("mNClsITS"), track.itsNCls()); - histos.fill(HIST("mChi2TPC"), track.tpcChi2NCl()); - histos.fill(HIST("mChi2ITS"), track.itsChi2NCl()); - histos.fill(HIST("mChi2TRD"), track.trdChi2()); - histos.fill(HIST("mDCAxy"), track.dcaXY()); - histos.fill(HIST("mDCAx"), track.dcaZ()); - histos.fill(HIST("mDCAxyPt"), track.pt(), track.dcaXY()); - histos.fill(HIST("mDCAzPt"), track.pt(), track.dcaZ()); - histos.fill(HIST("mNSharedClsTPC"), track.tpcNClsShared()); - histos.fill(HIST("mCrossedRowsTPC"), track.tpcNClsCrossedRows()); - histos.fill(HIST("mNFinClsminusCRows"), track.tpcNClsFindableMinusCrossedRows()); - histos.fill(HIST("mNFractionShClsTPC"), track.tpcFractionSharedCls()); - histos.fill(HIST("mSharedClsvsPt"), track.pt(), track.tpcNClsShared()); - histos.fill(HIST("mSharedClsProbvsPt"), track.pt(), track.tpcFractionSharedCls() / track.tpcNClsCrossedRows()); + fillTrackQA(track); checkpT(track); } - for (int iPt = 0; iPt < numPt; ++iPt) { - if (countTracks[iPt] > 0) { - mHistArrQA[iPt * 4 + 3]->Fill(countTracks[iPt]); - } - } - // Calculate the normalized factorial moments - calculateMoments(mHistArrReset); + endEvent(); } PROCESS_SWITCH(FactorialMomentsTask, processRun2, "for RUN2", false); }; + WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { return WorkflowSpec{