diff --git a/PWGLF/Tasks/GlobalEventProperties/CMakeLists.txt b/PWGLF/Tasks/GlobalEventProperties/CMakeLists.txt index ed61c203672..258a43a91d0 100644 --- a/PWGLF/Tasks/GlobalEventProperties/CMakeLists.txt +++ b/PWGLF/Tasks/GlobalEventProperties/CMakeLists.txt @@ -23,3 +23,8 @@ o2physics_add_dpl_workflow(dndeta-mft-pp SOURCES dndeta-mft-pp.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(flattenicty-pikp + SOURCES flattenictyPikp.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::AnalysisCCDB + COMPONENT_NAME Analysis) diff --git a/PWGLF/Tasks/GlobalEventProperties/flattenictyPikp.cxx b/PWGLF/Tasks/GlobalEventProperties/flattenictyPikp.cxx new file mode 100644 index 00000000000..df99f8b7a6f --- /dev/null +++ b/PWGLF/Tasks/GlobalEventProperties/flattenictyPikp.cxx @@ -0,0 +1,1939 @@ +// Copyright 2020-2022 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// 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. + +/// \file flattenictyPikp.cxx +/// \author Gyula Bencedi, gyula.bencedi@cern.ch +/// \brief Task to produce pion, kaon, proton high-pT +/// distributions as a function of charged-particle flattenicity +/// \since 26 June 2025 + +#include "PWGLF/DataModel/LFParticleIdentification.h" +#include "PWGLF/DataModel/LFStrangenessTables.h" +#include "PWGLF/DataModel/mcCentrality.h" +#include "PWGLF/Utils/inelGt.h" + +#include "Common/Core/TrackSelection.h" +#include "Common/Core/TrackSelectionDefaults.h" +#include "Common/DataModel/Centrality.h" +#include "Common/DataModel/EventSelection.h" +#include "Common/DataModel/FT0Corrected.h" +#include "Common/DataModel/McCollisionExtra.h" +#include "Common/DataModel/Multiplicity.h" +#include "Common/DataModel/PIDResponse.h" +#include "Common/DataModel/TrackSelectionTables.h" + +#include "CommonConstants/MathConstants.h" +#include "DataFormatsFIT/Triggers.h" +#include "DetectorsCommonDataFormats/AlignParam.h" +#include "Framework/ASoAHelpers.h" +#include "Framework/AnalysisDataModel.h" +#include "Framework/AnalysisTask.h" +#include "Framework/Configurable.h" +#include "Framework/HistogramRegistry.h" +#include "Framework/O2DatabasePDGPlugin.h" +#include "Framework/StaticFor.h" +#include "Framework/runDataProcessing.h" +#include "ReconstructionDataFormats/PID.h" +#include "ReconstructionDataFormats/Track.h" +#include +#include +#include + +#include "TEfficiency.h" +#include "THashList.h" +#include "TPDGCode.h" +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace o2; +using namespace o2::framework; +using namespace o2::constants::math; + +auto static constexpr kMinCharge = 3.f; +// FV0 specific constants +static constexpr int kMaxRingsFV0 = 5; +static constexpr int kNCellsFV0 = 48; +static constexpr int kInnerFV0 = 32; +static constexpr float kFV0IndexPhi[5] = {0., 8., 16., 24., 32.}; +static constexpr float kEtaMaxFV0 = 5.1; +static constexpr float kEtaMinFV0 = 2.2; +static constexpr float kDEtaFV0 = (kEtaMaxFV0 - kEtaMinFV0) / kMaxRingsFV0; +// PID names +static constexpr int kProcessIdWeak = 4; +static constexpr int Ncharges = 2; +static constexpr o2::track::PID::ID Npart = 5; +static constexpr o2::track::PID::ID NpartChrg = Npart * Ncharges; +static constexpr int PDGs[] = {11, 13, 211, 321, 2212}; +static constexpr int PidSgn[NpartChrg] = {11, 13, 211, 321, 2212, -11, -13, -211, -321, -2212}; +static constexpr const char* Pid[Npart] = {"el", "mu", "pi", "ka", "pr"}; +static constexpr const char* PidChrg[NpartChrg] = {"e^{-}", "#mu^{-}", "#pi^{+}", "K^{+}", "p", "e^{+}", "#mu^{+}", "#pi^{-}", "K^{-}", "#bar{p}"}; +static constexpr std::string_view kSpecies[NpartChrg] = {"Elminus", "Muplus", "PiPlus", "KaPlus", "Pr", "ElPlus", "MuMinus", "PiMinus", "KaMinus", "PrBar"}; +static constexpr std::string_view kSpeciesAll[Npart] = {"El", "Mu", "Pi", "Ka", "Pr"}; +// histogram naming +static constexpr std::string_view kCharge[] = {"all/", "pos/", "neg/"}; +static constexpr std::string_view kPrefix = "Tracks/"; +static constexpr std::string_view kQA = "QA/"; +static constexpr std::string_view kPrefixCleanTof = "Tracks/CleanTof/"; +static constexpr std::string_view kPrefixCleanV0 = "Tracks/CleanV0/"; +static constexpr std::string_view kStatus[] = {"preSel/", "postSel/"}; +static constexpr std::string_view kStatCalib[] = {"preCalib/", "postCalib/"}; +static constexpr std::string_view kPtGenPrimSgn = "/hPtGenPrimSgn"; +static constexpr std::string_view kPtGenPrimSgnF = "Tracks/{}/hPtGenPrimSgn"; +static constexpr std::string_view kPtGenPrimSgnINEL = "/hPtGenPrimSgnINEL"; +static constexpr std::string_view kPtGenPrimSgnINELF = "Tracks/{}/hPtGenPrimSgnINEL"; +static constexpr std::string_view kPtRecCollPrimSgn = "/hPtRecCollPrimSgn"; +static constexpr std::string_view kPtRecCollPrimSgnF = "Tracks/{}/hPtRecCollPrimSgn"; +static constexpr std::string_view kPtRecCollPrimSgnINEL = "/hPtRecCollPrimSgnINEL"; +static constexpr std::string_view kPtRecCollPrimSgnINELF = "Tracks/{}/hPtRecCollPrimSgnINEL"; +static constexpr std::string_view kPtGenRecCollPrimSgn = "/hPtGenRecCollPrimSgn"; +static constexpr std::string_view kPtGenRecCollPrimSgnF = "Tracks/{}/hPtGenRecCollPrimSgn"; +static constexpr std::string_view kPtGenRecCollPrimSgnINEL = "/hPtGenRecCollPrimSgnINEL"; +static constexpr std::string_view kPtGenRecCollPrimSgnINELF = "Tracks/{}/hPtGenRecCollPrimSgnINEL"; +static constexpr std::string_view kPtMCclosurePrim = "/hPtMCclosurePrim"; +static constexpr std::string_view kPtMCclosurePrimF = "Tracks/{}/hPtMCclosurePrim"; + +enum V0sSel { + kNaN = -1, + kGa = 0, + kKz = 1, + kLam = 2, + kaLam = 3 +}; + +enum TrkSelNoFilt { + trkSelEta, + trkSelPt, + trkSelDCA, + trkSelNCls, + trkSelTPCBndr, + nTrkSel +}; + +enum EvtSel { + evtSelAll, + evtSelSel8, + evtSelNoITSROFrameBorder, + evtSelkNoTimeFrameBorder, + evtSelkNoSameBunchPileup, + evtSelkIsGoodZvtxFT0vsPV, + evtSelkIsVertexITSTPC, + evtSelkIsVertexTOFmatched, + evtSelVtxZ, + evtSelINELgt0, + nEvtSel +}; + +enum FillType { + kBefore, + kAfter +}; + +enum ChargeType { + kAll, + kPos, + kNeg +}; + +struct MultE { + static constexpr int kNoMult = 0; + static constexpr int kMultFT0M = 1; + static constexpr int kMultTPC = 2; +}; + +std::array rhoLatticeFV0{0}; +std::array fv0AmplitudeWoCalib{0}; + +std::array, NpartChrg> hPtEffRecPrim{}; +std::array, NpartChrg> hPtEffRecWeak{}; +std::array, NpartChrg> hPtEffRecMat{}; +std::array, NpartChrg> hPtEffRec{}; +std::array, NpartChrg> hPtEffGen{}; +std::array, NpartChrg> hPtGenRecEvt{}; +std::array, NpartChrg> hPtGenPrimRecEvt{}; +std::array, NpartChrg> hPtEffGenPrim{}; +std::array, NpartChrg> hPtEffGenWeak{}; +std::array, NpartChrg> hPtEffGenMat{}; +std::array, NpartChrg> hPtVsDCAxyPrim{}; +std::array, NpartChrg> hPtVsDCAxyWeak{}; +std::array, NpartChrg> hPtVsDCAxyMat{}; +std::array, NpartChrg> hDCAxyBadCollPrim{}; +std::array, NpartChrg> hDCAxyBadCollWeak{}; +std::array, NpartChrg> hDCAxyBadCollMat{}; +std::array, NpartChrg> hPtNsigmaTPC{}; +std::array, NpartChrg> hThPtNsigmaTPC{}; +std::array, NpartChrg> hPtNsigmaTOF{}; +std::array, NpartChrg> hPtNsigmaTPCTOF{}; + +struct FlattenictyPikp { + + HistogramRegistry flatchrg{"flatchrg", {}, OutputObjHandlingPolicy::AnalysisObject, true, false}; + OutputObj listEfficiency{"Efficiency"}; + Service pdg; + + std::vector fv0AmplCorr{}; + TProfile2D* zVtxMap = nullptr; + int runNumber{-1}; + + Configurable multEst{"multEst", 1, "0: without multiplicity; 1: MultFT0M; 2: MultTPC"}; + Configurable applyCalibGain{"applyCalibGain", true, "equalize detector amplitudes"}; + Configurable applyCalibVtx{"applyCalibVtx", false, "equalize Amp vs vtx"}; + Configurable applyCalibDeDx{"applyCalibDeDx", true, "calibration of dedx signal"}; + Configurable cfgFillQAHisto{"cfgFillQAHisto", true, "fill QA histograms"}; + Configurable cfgFillChrgType{"cfgFillChrgType", true, "fill histograms per charge types"}; + Configurable> paramsFuncMIPposEta{"paramsFuncMIPposEta", std::vector{-1.f}, "parameters of pol2"}; + Configurable> paramsFuncMIPnegEta{"paramsFuncMIPnegEta", std::vector{-1.f}, "parameters of pol2"}; + Configurable> paramsFuncPlateuPosEta{"paramsFuncPlateuPosEta", std::vector{-1.f}, "parameters of pol2"}; + Configurable> paramsFuncPlateuNegEta{"paramsFuncPlateuNegEta", std::vector{-1.f}, "parameters of pol2"}; + Configurable cfgGainEqCcdbPath{"cfgGainEqCcdbPath", "Users/g/gbencedi/flattenicity/GainEq", "CCDB path for gain equalization constants"}; + Configurable cfgVtxEqCcdbPath{"cfgVtxEqCcdbPath", "Users/g/gbencedi/flattenicity/ZvtxEq", "CCDB path for z-vertex equalization constants"}; + Configurable cfgUseCcdbForRun{"cfgUseCcdbForRun", true, "Get ccdb object based on run number instead of timestamp"}; + Configurable cfgStoreThnSparse{"cfgStoreThnSparse", false, "Store histograms as THnSparse"}; + + struct : ConfigurableGroup { + Configurable cfgCustomTVX{"cfgCustomTVX", false, "Ask for custom TVX instead of sel8"}; + Configurable cfgRemoveNoTimeFrameBorder{"cfgRemoveNoTimeFrameBorder", true, "Bunch crossing is far from Time Frame borders"}; + Configurable cfgRemoveITSROFrameBorder{"cfgRemoveITSROFrameBorder", true, "Bunch crossing is far from ITS RO Frame border"}; + Configurable cfgCutVtxZ{"cfgCutVtxZ", 10.0f, "Accepted z-vertex range"}; + Configurable cfgINELCut{"cfgINELCut", true, "INEL event selection"}; + Configurable cfgRemoveNoSameBunchPileup{"cfgRemoveNoSameBunchPileup", false, "Reject collisions in case of pileup with another collision in the same foundBC"}; + Configurable cfgRequireIsGoodZvtxFT0vsPV{"cfgRequireIsGoodZvtxFT0vsPV", false, "Small difference between z-vertex from PV and from FT0"}; + Configurable cfgRequireIsVertexITSTPC{"cfgRequireIsVertexITSTPC", false, "At least one ITS-TPC track (reject vertices built from ITS-only tracks)"}; + Configurable cfgRequirekIsVertexTOFmatched{"cfgRequirekIsVertexTOFmatched", false, "Require kIsVertexTOFmatched: at least one of vertex contributors is matched to TOF"}; + } evtSelOpt; + + struct : ConfigurableGroup { + Configurable useFlatData{"useFlatData", true, "use flattenicity from rec collisions"}; + } flatSelOpt; + + struct : ConfigurableGroup { + ConfigurableAxis axisPt{"axisPt", {VARIABLE_WIDTH, 0.0, 0.1, 0.12, 0.14, 0.16, 0.18, 0.2, 0.25, 0.3, 0.35, 0.4, 0.45, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9, 0.95, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0, 2.2, 2.4, 2.6, 2.8, 3.0, 3.2, 3.4, 3.6, 3.8, 4.0, 4.2, 4.4, 4.6, 4.8, 5.0}, "pT binning"}; + ConfigurableAxis axisMultPerc{"axisMultPerc", {100, 0, 100}, "Multiplicity percentiles binning"}; + ConfigurableAxis axisVertexZ{"axisVertexZ", {60, -15., 15.}, "Vertex z binning"}; + ConfigurableAxis axisMult{"axisMult", {301, -0.5, 300.5}, "Multiplicity binning"}; + ConfigurableAxis axisDCAxy{"axisDCAxy", {200, -5, 5}, "DCAxy binning"}; + ConfigurableAxis axisDCAz{"axisDCAz", {200, -5, 5}, "DCAz binning"}; + ConfigurableAxis axisPhi = {"axisPhi", {60, 0, constants::math::TwoPI}, "#varphi binning"}; + ConfigurableAxis axisPhiMod = {"axisPhiMod", {100, 0, constants::math::PI / 9}, "fmod(#varphi,#pi/9)"}; + ConfigurableAxis axisEta = {"axisEta", {8, -0.8, 0.8}, "#eta binning"}; + ConfigurableAxis axisDedx{"axisDedx", {1000, 0, 1000}, "dE/dx binning"}; + ConfigurableAxis axisNsigmaTPC{"axisNsigmaTPC", {200, -10, 10}, "nsigmaTPC binning"}; + ConfigurableAxis axisNsigmaTOF{"axisNsigmaTOF", {200, -10, 10}, "nsigmaTOF binning"}; + } binOpt; + + struct : ConfigurableGroup { + Configurable cfgTrkEtaMax{"cfgTrkEtaMax", 0.8f, "Eta range for tracks"}; + Configurable cfgRapMax{"cfgRapMax", 0.5f, "Maximum range of rapidity for tracks"}; + Configurable cfgTrkPtMin{"cfgTrkPtMin", 0.15f, "Minimum pT of tracks"}; + Configurable cfgPhiCutPtMin{"cfgPhiCutPtMin", 2.0f, "Minimum pT for phi cut"}; + Configurable cfgTOFBetaPion{"cfgTOFBetaPion", 1.0f, "Minimum beta for TOF pions"}; + Configurable cfgUseExtraTrkCut{"cfgUseExtraTrkCut", true, "Use extra track cut"}; + Configurable cfgMomMIPMax{"cfgMomMIPMax", 0.6f, "Maximum momentum of MIP pions"}; + Configurable cfgMomMIPMin{"cfgMomMIPMin", 0.4f, "Minimum momentum of MIP pions"}; + Configurable cfgDeDxMIPMax{"cfgDeDxMIPMax", 60.0f, "Maximum range of MIP dedx"}; + Configurable cfgDeDxMIPMin{"cfgDeDxMIPMin", 40.0f, "Maximum range of MIP dedx"}; + Configurable cfgNsigmaMax{"cfgNsigmaMax", 100.0f, "Maximum range of nsgima for tracks"}; + Configurable cfgMomSelPiTOF{"cfgMomSelPiTOF", 0.7f, "Momentum cut for TOF pions"}; + Configurable cfgNsigSelKaTOF{"cfgNsigSelKaTOF", 3.0f, "Nsigma cut for TOF kaons"}; + Configurable cfgBetaPlateuMax{"cfgBetaPlateuMax", 0.1f, "Beta max for Plateau electrons"}; + } trkSelOpt; + + struct : ConfigurableGroup { + Configurable cfgDCAv0daughter{"cfgDCAv0daughter", 0.3, "DCA of V0 daughter tracks"}; + Configurable cfgV0etamax{"cfgV0etamax", 0.9f, "max eta of V0s"}; + Configurable cfgV0DaughterTpcMomMax{"cfgV0DaughterTpcMomMax", 0.6f, "Maximum momentum of V0 daughter tracks in TPC"}; + Configurable cfgTPCnClsmin{"cfgTPCnClsmin", 50.0f, "cfgTPCnClsmin"}; + Configurable cfgDCAposToPV{"cfgDCAposToPV", 0.05f, "cfgDCAposToPV"}; + Configurable cfgv0cospa{"cfgv0cospa", 0.998, "V0 CosPA"}; + Configurable cfgv0Rmin{"cfgv0Rmin", 0.0, "cfgv0Rmin"}; + Configurable cfgv0Rmax{"cfgv0Rmax", 90.0, "cfgv0Rmax"}; + Configurable cfgdmassG{"cfgdmassG", 0.1, "max mass for photon"}; + Configurable cfgdmassK{"cfgdmassK", 0.1, "max mass for K0s"}; + Configurable cfgdmassL{"cfgdmassL", 0.1, "max mass for Lambda"}; + Configurable cfgdmassAL{"cfgdmassAL", 0.1, "max mass for anti Lambda"}; + Configurable cfgNsigmaElTPC{"cfgNsigmaElTPC", 5.0, "max nsigma of TPC for electorn"}; + Configurable cfgNsigmaPiTPC{"cfgNsigmaPiTPC", 5.0, "max nsigma of TPC for pion"}; + Configurable cfgNsigmaPrTPC{"cfgNsigmaPrTPC", 5.0, "max nsigma of TPC for proton"}; + Configurable cfgNsigmaElTOF{"cfgNsigmaElTOF", 3.0, "max nsigma of TOF for electorn"}; + Configurable cfgNsigmaPiTOF{"cfgNsigmaPiTOF", 3.0, "max nsigma of TOF for pion"}; + Configurable cfgNsigmaPrTOF{"cfgNsigmaPrTOF", 3.0, "max nsigma of TOF for proton"}; + } v0SelOpt; + + Service ccdb; + struct : ConfigurableGroup { + Configurable cfgMagField{"cfgMagField", 99999, "Configurable magnetic field;default CCDB will be queried"}; + Configurable ccdbUrl{"ccdbUrl", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; + Configurable grpmagPath{"grpmagPath", "GLO/Config/GRPMagField", "CCDB path of the GRPMagField object"}; + Configurable grpPath{"grpPath", "GLO/GRP/GRP", "Path of the grp file"}; + } ccdbConf; + + TrackSelection mTrackSelector; + Configurable isCustomTracks{"isCustomTracks", false, "Use custom track cuts"}; + Configurable requirePt{"requirePt", 0.15f, "Set minimum pT of tracks"}; + Configurable requireEta{"requireEta", 0.8f, "Set eta range of tracks"}; + Configurable setITSreq{"setITSreq", 1, "0 = Run3ITSibAny, 1 = Run3ITSallAny, 2 = Run3ITSall7Layers, 3 = Run3ITSibTwo"}; + Configurable requireITS{"requireITS", true, "Additional cut on the ITS requirement"}; + Configurable requireTPC{"requireTPC", true, "Additional cut on the TPC requirement"}; + Configurable requireGoldenChi2{"requireGoldenChi2", true, "Additional cut on the GoldenChi2"}; + Configurable minNCrossedRowsTPC{"minNCrossedRowsTPC", 70.f, "Additional cut on the minimum number of crossed rows in the TPC"}; + Configurable minNCrossedRowsOverFindableClustersTPC{"minNCrossedRowsOverFindableClustersTPC", 0.8f, "Additional cut on the minimum value of the ratio between crossed rows and findable clusters in the TPC"}; + Configurable maxChi2PerClusterTPC{"maxChi2PerClusterTPC", 4.f, "Additional cut on the maximum value of the chi2 per cluster in the TPC"}; + Configurable maxChi2PerClusterITS{"maxChi2PerClusterITS", 36.f, "Additional cut on the maximum value of the chi2 per cluster in the ITS"}; + Configurable minITSnClusters{"minITSnClusters", 5, "minimum number of found ITS clusters"}; + Configurable maxDcaXYFactor{"maxDcaXYFactor", 1.f, "Multiplicative factor on the maximum value of the DCA xy"}; + Configurable maxDcaZ{"maxDcaZ", 2.f, "Additional cut on the maximum value of the DCA z"}; + Configurable minTPCNClsFound{"minTPCNClsFound", 0.f, "Additional cut on the minimum value of the number of found clusters in the TPC"}; + + TF1* fPhiCutLow = nullptr; + TF1* fPhiCutHigh = nullptr; + + std::unique_ptr fDeDxVsEtaPos = nullptr; + std::unique_ptr fDeDxVsEtaNeg = nullptr; + std::unique_ptr fEDeDxVsEtaPos = nullptr; + std::unique_ptr fEDeDxVsEtaNeg = nullptr; + + void init(InitContext&) + { + auto vecParamsMIPposEta = (std::vector)paramsFuncMIPposEta; + auto vecParamsMIPnegEta = (std::vector)paramsFuncMIPnegEta; + auto vecParamsPlateuPosEta = (std::vector)paramsFuncPlateuPosEta; + auto vecParamsPlateuNegEta = (std::vector)paramsFuncPlateuNegEta; + + fDeDxVsEtaPos = setFuncPars(vecParamsMIPposEta); + fDeDxVsEtaNeg = setFuncPars(vecParamsMIPnegEta); + fEDeDxVsEtaPos = setFuncPars(vecParamsPlateuPosEta); + fEDeDxVsEtaNeg = setFuncPars(vecParamsPlateuNegEta); + + ccdb->setURL(ccdbConf.ccdbUrl.value); + ccdb->setCaching(true); + ccdb->setLocalObjectValidityChecking(); + ccdb->setCreatedNotAfter(std::chrono::duration_cast(std::chrono::system_clock::now().time_since_epoch()).count()); + ccdb->setFatalWhenNull(false); + + if (isCustomTracks.value) { + mTrackSelector = getGlobalTrackSelectionRun3ITSMatch(setITSreq.value); + mTrackSelector.SetPtRange(requirePt.value, 1e10f); + mTrackSelector.SetEtaRange(-requireEta.value, requireEta.value); + mTrackSelector.SetRequireITSRefit(requireITS.value); + mTrackSelector.SetRequireTPCRefit(requireTPC.value); + mTrackSelector.SetRequireGoldenChi2(requireGoldenChi2.value); + mTrackSelector.SetMaxChi2PerClusterTPC(maxChi2PerClusterTPC.value); + mTrackSelector.SetMaxChi2PerClusterITS(maxChi2PerClusterITS.value); + mTrackSelector.SetMinNClustersITS(minITSnClusters.value); + mTrackSelector.SetMinNCrossedRowsTPC(minNCrossedRowsTPC.value); + mTrackSelector.SetMinNClustersTPC(minTPCNClsFound.value); + mTrackSelector.SetMinNCrossedRowsOverFindableClustersTPC(minNCrossedRowsOverFindableClustersTPC.value); + // // mTrackSelector.SetMaxDcaXYPtDep([](float pt) { return 0.0105f + 0.0350f / pow(pt, 1.1f); }); + mTrackSelector.SetMaxDcaXYPtDep([](float /*pt*/) { return 10000.f; }); + mTrackSelector.SetMaxDcaZ(maxDcaZ.value); + mTrackSelector.print(); + } + + const AxisSpec chargeAxis{2, -2.f, 2.f, "Charge"}; + const AxisSpec dEdxAxis{binOpt.axisDedx, "TPC dEdx (a.u.)"}; + const AxisSpec vtxzAxis{100, -20, 20, "Z_{vtx} (cm)"}; + const AxisSpec flatAxis{102, -0.01, 1.01, "Flat FV0"}; + const AxisSpec etaAxis{binOpt.axisEta, "#eta"}; + const AxisSpec phiAxis{binOpt.axisPhi, "#varphi"}; + const AxisSpec phiAxisMod{binOpt.axisPhiMod, "fmod(#varphi,#pi/9)"}; + const AxisSpec pAxis{binOpt.axisPt, "#it{p} (GeV/#it{c})"}; + const AxisSpec ptAxis{binOpt.axisPt, "#it{p}_{T} (GeV/#it{c})"}; + const AxisSpec dcaXYAxis{binOpt.axisDCAxy, "DCA_{xy} (cm)"}; + const AxisSpec dcaZAxis{binOpt.axisDCAz, "DCA_{z} (cm)"}; + const AxisSpec shCluserAxis{200, 0, 1, "Fraction of shared TPC clusters"}; + const AxisSpec clTpcAxis{160, 0, 160, "Number of clusters in TPC"}; + const AxisSpec nSigmaTPCAxis{binOpt.axisNsigmaTPC, "n#sigma_{TPC}"}; + const AxisSpec nSigmaTOFAxis{binOpt.axisNsigmaTOF, "n#sigma_{TOF}"}; + const AxisSpec amplitudeFT0 = {5000, 0, 10000, "FT0 amplitude"}; + const AxisSpec channelFT0Axis = {220, 0.0, 220.0, "FT0 channel"}; + + AxisSpec multAxis{binOpt.axisMultPerc, "multiplicity estimator"}; + + switch (multEst) { + case MultE::kNoMult: + break; + case MultE::kMultFT0M: + multAxis.name = "multFT0M"; + break; + case MultE::kMultTPC: + multAxis.name = "multTPC"; + break; + default: + LOG(fatal) << "No valid option for mult estimator " << multEst; + } + + // Event counter + flatchrg.add("Events/hEvtSel", "Number of events; Cut; #Events Passed Cut", {HistType::kTH1D, {{nEvtSel, 0, nEvtSel}}}); + flatchrg.get(HIST("Events/hEvtSel"))->GetXaxis()->SetBinLabel(evtSelAll + 1, "Events read"); + flatchrg.get(HIST("Events/hEvtSel"))->GetXaxis()->SetBinLabel(evtSelSel8 + 1, "Evt. sel8"); + flatchrg.get(HIST("Events/hEvtSel"))->GetXaxis()->SetBinLabel(evtSelNoITSROFrameBorder + 1, "NoITSROFrameBorder"); + flatchrg.get(HIST("Events/hEvtSel"))->GetXaxis()->SetBinLabel(evtSelkNoTimeFrameBorder + 1, "NoTimeFrameBorder"); + flatchrg.get(HIST("Events/hEvtSel"))->GetXaxis()->SetBinLabel(evtSelkNoSameBunchPileup + 1, "NoSameBunchPileup"); + flatchrg.get(HIST("Events/hEvtSel"))->GetXaxis()->SetBinLabel(evtSelkIsGoodZvtxFT0vsPV + 1, "IsGoodZvtxFT0vsPV"); + flatchrg.get(HIST("Events/hEvtSel"))->GetXaxis()->SetBinLabel(evtSelkIsVertexITSTPC + 1, "IsVertexITSTPC"); + flatchrg.get(HIST("Events/hEvtSel"))->GetXaxis()->SetBinLabel(evtSelkIsVertexTOFmatched + 1, "IsVertexTOFmatched"); + flatchrg.get(HIST("Events/hEvtSel"))->GetXaxis()->SetBinLabel(evtSelVtxZ + 1, "Vtx-z pos"); + flatchrg.get(HIST("Events/hEvtSel"))->GetXaxis()->SetBinLabel(evtSelINELgt0 + 1, "INEL>0"); + // Track counter + flatchrg.add("Tracks/hTrkSel", "Number of tracks; Cut; #Tracks Passed Cut", {HistType::kTH1D, {{nTrkSel, 0, nTrkSel}}}); + flatchrg.get(HIST("Tracks/hTrkSel"))->GetXaxis()->SetBinLabel(trkSelEta + 1, "Eta"); + flatchrg.get(HIST("Tracks/hTrkSel"))->GetXaxis()->SetBinLabel(trkSelPt + 1, "Pt"); + flatchrg.get(HIST("Tracks/hTrkSel"))->GetXaxis()->SetBinLabel(trkSelDCA + 1, "DCA"); + flatchrg.get(HIST("Tracks/hTrkSel"))->GetXaxis()->SetBinLabel(trkSelNCls + 1, "NCls"); + flatchrg.get(HIST("Tracks/hTrkSel"))->GetXaxis()->SetBinLabel(trkSelTPCBndr + 1, "TPC Boundary"); + + if (trkSelOpt.cfgUseExtraTrkCut) { + fPhiCutLow = new TF1("fPhiCutLow", "0.06/x+pi/18.0-0.06", 0, 100); + fPhiCutHigh = new TF1("fPhiCutHigh", "0.1/x+pi/18.0+0.06", 0, 100); + } + + if (doprocessFlat) { + flatchrg.add("Events/hVtxZ", "Measured vertex z position", HistType::kTH1D, {vtxzAxis}); + flatchrg.add("Events/hFlatVsMultEst", "hFlatVsMultEst", HistType::kTH2D, {flatAxis, multAxis}); + + if (cfgFillQAHisto) { + flatchrg.add("Tracks/postSel/hPtPhi", "; #it{p}_{T} (GeV/#it{c}); fmod(#varphi,#pi/9)", {HistType::kTH2D, {ptAxis, phiAxisMod}}); + // P vs PT vs ETA + flatchrg.add("Tracks/postSel/hPVsPtEta", "; #it{p} (GeV/#it{c}); #it{p}_{T} (GeV/#it{c}); #eta;", {HistType::kTH3D, {pAxis, ptAxis, etaAxis}}); + flatchrg.addClone("Tracks/postSel/", "Tracks/preSel/"); + // dEdx MIP, Plateau + flatchrg.add("Tracks/postCalib/all/hMIP", "; mult; flat; #eta; #LT dE/dx #GT_{MIP, primary tracks};", {HistType::kTHnSparseD, {multAxis, flatAxis, etaAxis, dEdxAxis}}); + flatchrg.add("Tracks/postCalib/all/hPlateau", "; mult; flat; #eta; #LT dE/dx #GT_{Plateau, primary tracks};", {HistType::kTHnSparseD, {multAxis, flatAxis, etaAxis, dEdxAxis}}); + flatchrg.add("Tracks/postCalib/all/hMIPVsEta", "; #eta; #LT dE/dx #GT_{MIP, primary tracks};", {HistType::kTH2D, {etaAxis, dEdxAxis}}); + flatchrg.add("Tracks/postCalib/all/pMIPVsEta", "; #eta; #LT dE/dx #GT_{MIP, primary tracks};", {HistType::kTProfile, {etaAxis}}); + flatchrg.add("Tracks/postCalib/all/hMIPVsPhi", "; #varphi; #LT dE/dx #GT_{MIP, primary tracks};", {HistType::kTH2D, {phiAxis, dEdxAxis}}); + flatchrg.add("Tracks/postCalib/all/pMIPVsPhi", "; #varphi; #LT dE/dx #GT_{MIP, primary tracks};", {HistType::kTProfile, {phiAxis}}); + flatchrg.add("Tracks/postCalib/all/hMIPVsPhiVsEta", "; #varphi; #LT dE/dx #GT_{MIP, primary tracks}; #eta;", {HistType::kTH3D, {phiAxis, dEdxAxis, etaAxis}}); + flatchrg.add("Tracks/postCalib/all/hPlateauVsEta", "; #eta; #LT dE/dx #GT_{Plateau, primary tracks};", {HistType::kTH2D, {etaAxis, dEdxAxis}}); + flatchrg.add("Tracks/postCalib/all/pPlateauVsEta", "; #eta; #LT dE/dx #GT_{Plateau, primary tracks};", {HistType::kTProfile, {etaAxis}}); + flatchrg.add("Tracks/postCalib/all/hPlateauVsPhi", "; #varphi; #LT dE/dx #GT_{Plateau, primary tracks};", {HistType::kTH2D, {phiAxis, dEdxAxis}}); + flatchrg.add("Tracks/postCalib/all/pPlateauVsPhi", "; #varphi; #LT dE/dx #GT_{Plateau, primary tracks};", {HistType::kTProfile, {phiAxis}}); + flatchrg.add("Tracks/postCalib/all/hPlateauVsPhiVsEta", "; #varphi; #LT dE/dx #GT_{Plateau, primary tracks}; #eta;", {HistType::kTH3D, {phiAxis, dEdxAxis, etaAxis}}); + if (cfgFillChrgType) { + flatchrg.addClone("Tracks/postCalib/all/", "Tracks/postCalib/pos/"); + flatchrg.addClone("Tracks/postCalib/all/", "Tracks/postCalib/neg/"); + flatchrg.addClone("Tracks/postCalib/all/", "Tracks/preCalib/all/"); + flatchrg.addClone("Tracks/preCalib/all/", "Tracks/preCalib/pos/"); + flatchrg.addClone("Tracks/preCalib/all/", "Tracks/preCalib/neg/"); + } + // fillTrackQA + flatchrg.add("Tracks/QA/hPtVsWOcutDCA", "hPtVsWOcutDCA", HistType::kTH2D, {ptAxis, dcaXYAxis}); + flatchrg.add("Tracks/QA/hPt", "", HistType::kTH1D, {ptAxis}); + flatchrg.add("Tracks/QA/hPhi", "", HistType::kTH1D, {phiAxis}); + flatchrg.add("Tracks/QA/hEta", "", HistType::kTH1D, {etaAxis}); + flatchrg.add("Tracks/QA/hDCAXYvsPt", "", HistType::kTH2D, {ptAxis, dcaXYAxis}); + flatchrg.add("Tracks/QA/hDCAZvsPt", "", HistType::kTH2D, {ptAxis, dcaZAxis}); + // tpc + flatchrg.add("Tracks/QA/hShTpcClvsPt", "", {HistType::kTH2D, {ptAxis, shCluserAxis}}); + flatchrg.add("Tracks/QA/hCrossTPCvsPt", "", {HistType::kTH2D, {ptAxis, clTpcAxis}}); + flatchrg.add("Tracks/QA/hTPCCluster", "N_{cluster}", HistType::kTH1D, {{200, -0.5, 199.5}}); + flatchrg.add("Tracks/QA/tpcNClsShared", " ; # shared TPC clusters TPC", HistType::kTH1D, {{165, -0.5, 164.5}}); + flatchrg.add("Tracks/QA/tpcCrossedRows", " ; # crossed TPC rows", HistType::kTH1D, {{165, -0.5, 164.5}}); + flatchrg.add("Tracks/QA/tpcCrossedRowsOverFindableCls", " ; crossed rows / findable TPC clusters", HistType::kTH1D, {{60, 0.7, 1.3}}); + // its + flatchrg.add("Tracks/QA/itsNCls", " ; # ITS clusters", HistType::kTH1D, {{8, -0.5, 7.5}}); + flatchrg.add("Tracks/QA/hChi2ITSTrkSegment", "chi2ITS", HistType::kTH1D, {{100, -0.5, 99.5}}); + // tof + flatchrg.add("Tracks/QA/hTOFPvsBeta", "Beta from TOF; #it{p} (GeV/#it{c}); #beta", {HistType::kTH2D, {pAxis, {120, 0.0, 1.2}}}); + flatchrg.add("Tracks/QA/hTOFpi", "Primary Pions from TOF; #eta; #it{p} (GeV/#it{c}); dEdx", {HistType::kTHnSparseD, {etaAxis, pAxis, dEdxAxis}}); + // nsigma + for (int i = 0; i < NpartChrg; i++) { + const std::string strID = Form("/%s/%s", (i < Npart) ? "pos" : "neg", Pid[i % Npart]); + hPtNsigmaTPC[i] = flatchrg.add("Tracks/hPtNsigmaTPC" + strID, " ; p_{T} (GeV/c)", HistType::kTH2D, {ptAxis, nSigmaTPCAxis}); + if (cfgStoreThnSparse) { + hThPtNsigmaTPC[i] = flatchrg.add("Tracks/hThPtNsigmaTPC" + strID, " ; p_{T} (GeV/c)", HistType::kTHnSparseD, {ptAxis, nSigmaTPCAxis, multAxis, flatAxis}); + } + hPtNsigmaTOF[i] = flatchrg.add("Tracks/hPtNsigmaTOF" + strID, " ; p_{T} (GeV/c)", HistType::kTH2D, {ptAxis, nSigmaTOFAxis}); + hPtNsigmaTPCTOF[i] = flatchrg.add("Tracks/hPtNsigmaTPCTOF" + strID, PidChrg[i], HistType::kTH2D, {nSigmaTPCAxis, nSigmaTOFAxis}); + } + } + // FV0 QA + flatchrg.add("FV0/hFV0AmplWCalib", "", HistType::kTH2D, {{48, -0.5, 47.5, "channel"}, {500, -0.5, +19999.5, "FV0 amplitude"}}); + flatchrg.add("FV0/hFV0AmplvsVtxzWoCalib", "", HistType::kTH2D, {{30, -15.0, +15.0, "z vtx (cm)"}, {1000, -0.5, +39999.5, "FV0 amplitude"}}); + flatchrg.add("FV0/hFV0AmplvsVtxzCalib", "", HistType::kTH2D, {{30, -15.0, +15.0, "z vtx (cm)"}, {1000, -0.5, +39999.5, "FV0 amplitude"}}); + flatchrg.add("FV0/hFV0amp", "", {HistType::kTH2D, {channelFT0Axis, amplitudeFT0}}); + flatchrg.add("FV0/pFV0amp", "", HistType::kTProfile, {channelFT0Axis}); + flatchrg.add("FV0/hFV0ampCorr", "", {HistType::kTH2D, {channelFT0Axis, amplitudeFT0}}); + // V0's QA + flatchrg.add("Tracks/V0qa/hV0Pt", "pT", HistType::kTH1D, {{100, 0.0f, 10}}); + flatchrg.add("Tracks/V0qa/hV0ArmPod", ";#alpha; #it{q}_T", HistType::kTH2D, {{200, -1.0f, +1.0f}, {250, 0.0f, 0.25f}}); + // dEdx PID + flatchrg.add({"Tracks/all/hdEdx", "; #eta; mult; flat; #it{p} (GeV/#it{c}); dEdx", {HistType::kTHnSparseD, {etaAxis, multAxis, flatAxis, pAxis, dEdxAxis}}}); + // Clean samples + flatchrg.add({"Tracks/CleanTof/all/hPiTof", "; #eta; mult; flat; #it{p} (GeV/#it{c}); dEdx", {HistType::kTHnSparseD, {etaAxis, multAxis, flatAxis, pAxis, dEdxAxis}}}); + flatchrg.add({"Tracks/CleanV0/pos/hEV0", "; #eta; mult; flat; #it{p} (GeV/#it{c}); dEdx", {HistType::kTHnSparseD, {etaAxis, multAxis, flatAxis, pAxis, dEdxAxis}}}); + flatchrg.add({"Tracks/CleanV0/pos/hPiV0", "; #eta; mult; flat; #it{p} (GeV/#it{c}); dEdx", {HistType::kTHnSparseD, {etaAxis, multAxis, flatAxis, pAxis, dEdxAxis}}}); + flatchrg.add({"Tracks/CleanV0/pos/hPV0", "; #eta; mult; flat; #it{p} (GeV/#it{c}); dEdx", {HistType::kTHnSparseD, {etaAxis, multAxis, flatAxis, pAxis, dEdxAxis}}}); + flatchrg.addClone("Tracks/CleanV0/pos/", "Tracks/CleanV0/neg/"); + if (cfgFillChrgType) { + flatchrg.addClone("Tracks/all/", "Tracks/pos/"); + flatchrg.addClone("Tracks/all/", "Tracks/neg/"); + flatchrg.addClone("Tracks/CleanTof/all/", "Tracks/CleanTof/pos/"); + flatchrg.addClone("Tracks/CleanTof/all/", "Tracks/CleanTof/neg/"); + } + } + + if (doprocessMC) { + auto h = flatchrg.add("hEvtGenRec", "Generated and Reconstructed MC Collisions", kTH1D, {{3, 0.5, 3.5}}); + h->GetXaxis()->SetBinLabel(1, "Gen coll"); + h->GetXaxis()->SetBinLabel(2, "Rec coll"); + h->GetXaxis()->SetBinLabel(3, "INEL>0"); + + flatchrg.add("hEvtMcGenColls", "Number of events; Cut; #Events Passed Cut", {HistType::kTH1D, {{5, 0.5, 5.5}}}); + flatchrg.get(HIST("hEvtMcGenColls"))->GetXaxis()->SetBinLabel(1, "Gen. coll"); + flatchrg.get(HIST("hEvtMcGenColls"))->GetXaxis()->SetBinLabel(2, "At least 1 reco"); + flatchrg.get(HIST("hEvtMcGenColls"))->GetXaxis()->SetBinLabel(3, "Reco. coll."); + flatchrg.get(HIST("hEvtMcGenColls"))->GetXaxis()->SetBinLabel(4, "Reco. good coll."); + + for (int i = 0; i < NpartChrg; i++) { + const std::string strID = Form("/%s/%s", (i < Npart) ? "pos" : "neg", Pid[i % Npart]); + hPtGenRecEvt[i] = flatchrg.add("Tracks/hPtGenRecEvt" + strID, " ; p_{T} (GeV/c)", HistType::kTH1D, {ptAxis}); + hPtGenPrimRecEvt[i] = flatchrg.add("Tracks/hPtGenPrimRecEvt" + strID, " ; p_{T} (GeV/c)", HistType::kTH1D, {ptAxis}); + hPtEffGenPrim[i] = flatchrg.add("Tracks/hPtEffGenPrim" + strID, " ; p_{T} (GeV/c)", HistType::kTH3D, {multAxis, flatAxis, ptAxis}); + hPtEffGenWeak[i] = flatchrg.add("Tracks/hPtEffGenWeak" + strID, " ; p_{T} (GeV/c)", HistType::kTH3D, {multAxis, flatAxis, ptAxis}); + hPtEffGenMat[i] = flatchrg.add("Tracks/hPtEffGenMat" + strID, " ; p_{T} (GeV/c)", HistType::kTH3D, {multAxis, flatAxis, ptAxis}); + hPtEffRecPrim[i] = flatchrg.add("Tracks/hPtEffRecPrim" + strID, " ; p_{T} (GeV/c)", HistType::kTH3D, {multAxis, flatAxis, ptAxis}); + hPtEffRecWeak[i] = flatchrg.add("Tracks/hPtEffRecWeak" + strID, " ; p_{T} (GeV/c)", HistType::kTH3D, {multAxis, flatAxis, ptAxis}); + hPtEffRecMat[i] = flatchrg.add("Tracks/hPtEffRecMat" + strID, " ; p_{T} (GeV/c)", HistType::kTH3D, {multAxis, flatAxis, ptAxis}); + hDCAxyBadCollPrim[i] = flatchrg.add("Tracks/hDCAxyBadCollPrim" + strID, " ; p_{T} (GeV/c)", HistType::kTH2D, {ptAxis, dcaXYAxis}); + hDCAxyBadCollWeak[i] = flatchrg.add("Tracks/hDCAxyBadCollWeak" + strID, " ; p_{T} (GeV/c)", HistType::kTH2D, {ptAxis, dcaXYAxis}); + hDCAxyBadCollMat[i] = flatchrg.add("Tracks/hDCAxyBadCollMat" + strID, " ; p_{T} (GeV/c)", HistType::kTH2D, {ptAxis, dcaXYAxis}); + hPtVsDCAxyPrim[i] = flatchrg.add("Tracks/hPtVsDCAxyPrim" + strID, " ; p_{T} (GeV/c)", HistType::kTH2D, {ptAxis, dcaXYAxis}); + hPtVsDCAxyWeak[i] = flatchrg.add("Tracks/hPtVsDCAxyWeak" + strID, " ; p_{T} (GeV/c)", HistType::kTH2D, {ptAxis, dcaXYAxis}); + hPtVsDCAxyMat[i] = flatchrg.add("Tracks/hPtVsDCAxyMat" + strID, " ; p_{T} (GeV/c)", HistType::kTH2D, {ptAxis, dcaXYAxis}); + } + + flatchrg.add({"hPtOut", " ; p_{T} (GeV/c)", {HistType::kTH1D, {ptAxis}}}); + flatchrg.add({"hPtOutPrim", " ; p_{T} (GeV/c)", {HistType::kTH1D, {ptAxis}}}); + flatchrg.add({"hPtOutNoEtaCut", " ; p_{T} (GeV/c)", {HistType::kTH1D, {ptAxis}}}); + flatchrg.add({"PtOutFakes", " ; p_{T} (GeV/c)", {HistType::kTH1D, {ptAxis}}}); + flatchrg.add("hPtVsDCAxyPrimAll", "hPtVsDCAxyPrimAll", HistType::kTH2D, {ptAxis, dcaXYAxis}); + flatchrg.add("hPtVsDCAxyWeakAll", "hPtVsDCAxyWeakAll", HistType::kTH2D, {ptAxis, dcaXYAxis}); + flatchrg.add("hPtVsDCAxyMatAll", "hPtVsDCAxyMatAll", HistType::kTH2D, {ptAxis, dcaXYAxis}); + flatchrg.add("hPtVsDCAxyAll", "hPtVsDCAxyAll", HistType::kTH2D, {ptAxis, dcaXYAxis}); + flatchrg.add({"ResponseGen", " ; N_{part}; F_{FV0};", {HistType::kTHnSparseD, {multAxis, flatAxis}}}); + flatchrg.add("h1flatencityFV0MCGen", "", HistType::kTH1D, {{102, -0.01, 1.01, "1-flatencityFV0"}}); + } + + if (doprocessMCclosure) { + for (int i = 0; i < Npart; i++) { + flatchrg.add({fmt::format(kPtMCclosurePrimF.data(), kSpeciesAll[i]).c_str(), " ; p_{T} (GeV/c)", {HistType::kTH3D, {multAxis, flatAxis, ptAxis}}}); + } + } + + if (doprocessSgnLoss) { + flatchrg.add("hFlatMCGenRecColl", "hFlatMCGenRecColl", {HistType::kTH1D, {flatAxis}}); + flatchrg.add("hFlatMCGen", "hFlatMCGen", {HistType::kTH1D, {flatAxis}}); + // Event counter + flatchrg.add("hEvtMcGen", "hEvtMcGen", {HistType::kTH1D, {{4, 0.f, 4.f}}}); + flatchrg.get(HIST("hEvtMcGen"))->GetXaxis()->SetBinLabel(1, "all"); + flatchrg.get(HIST("hEvtMcGen"))->GetXaxis()->SetBinLabel(2, "z-vtx"); + flatchrg.get(HIST("hEvtMcGen"))->GetXaxis()->SetBinLabel(3, "INELgt0"); + flatchrg.add("hEvtMCRec", "hEvtMCRec", {HistType::kTH1D, {{4, 0.f, 4.f}}}); + flatchrg.get(HIST("hEvtMCRec"))->GetXaxis()->SetBinLabel(1, "all"); + flatchrg.get(HIST("hEvtMCRec"))->GetXaxis()->SetBinLabel(2, "evt sel"); + flatchrg.get(HIST("hEvtMCRec"))->GetXaxis()->SetBinLabel(3, "INELgt0"); + flatchrg.add("hEvtMcGenRecColl", "hEvtMcGenRecColl", {HistType::kTH1D, {{2, 0.f, 2.f}}}); + flatchrg.get(HIST("hEvtMcGenRecColl"))->GetXaxis()->SetBinLabel(1, "INEL"); + flatchrg.get(HIST("hEvtMcGenRecColl"))->GetXaxis()->SetBinLabel(2, "INELgt0"); + + flatchrg.add("hFlatGenINELgt0", "hFlatGenINELgt0", {HistType::kTH1D, {flatAxis}}); + + for (int i = 0; i < NpartChrg; ++i) { + flatchrg.add({fmt::format(kPtGenPrimSgnF.data(), kSpecies[i]).c_str(), " ; p_{T} (GeV/c)", {HistType::kTH3D, {multAxis, flatAxis, ptAxis}}}); + flatchrg.add({fmt::format(kPtGenPrimSgnINELF.data(), kSpecies[i]).c_str(), " ; p_{T} (GeV/c)", {HistType::kTH3D, {multAxis, flatAxis, ptAxis}}}); + flatchrg.add({fmt::format(kPtRecCollPrimSgnF.data(), kSpecies[i]).c_str(), " ; p_{T} (GeV/c)", {HistType::kTH3D, {multAxis, flatAxis, ptAxis}}}); + flatchrg.add({fmt::format(kPtRecCollPrimSgnINELF.data(), kSpecies[i]).c_str(), " ; p_{T} (GeV/c)", {HistType::kTH3D, {multAxis, flatAxis, ptAxis}}}); + flatchrg.add({fmt::format(kPtGenRecCollPrimSgnF.data(), kSpecies[i]).c_str(), " ; p_{T} (GeV/c)", {HistType::kTH3D, {multAxis, flatAxis, ptAxis}}}); + flatchrg.add({fmt::format(kPtGenRecCollPrimSgnINELF.data(), kSpecies[i]).c_str(), " ; p_{T} (GeV/c)", {HistType::kTH3D, {multAxis, flatAxis, ptAxis}}}); + } + } + + // Hash list for efficiency + listEfficiency.setObject(new THashList); + static_for<0, 1>([&](auto pidSgn) { + bookMcHist(); + bookMcHist(); + bookMcHist(); + initEfficiency(); + initEfficiency(); + initEfficiency(); + }); + } + + void initCCDB(aod::BCsWithTimestamps::iterator const& bc) + { + fv0AmplCorr.clear(); + fv0AmplCorr = {}; + std::string fullPathCalibGain; + std::string fullPathCalibVtx; + auto timestamp = bc.timestamp(); + auto runnumber = bc.runNumber(); + + if (applyCalibGain) { + fullPathCalibGain = cfgGainEqCcdbPath; + fullPathCalibGain += "/FV0"; + auto objfv0Gain = getForTsOrRun>(fullPathCalibGain, timestamp, runnumber); + if (!objfv0Gain) { + for (auto i{0u}; i < kNCellsFV0; i++) { + fv0AmplCorr.push_back(1.); + LOGF(warning, "Setting FV0 calibration object values to 1"); + } + } else { + fv0AmplCorr = *(objfv0Gain); + } + } + + if (applyCalibVtx) { + fullPathCalibVtx = cfgVtxEqCcdbPath; + fullPathCalibVtx += "/FV0"; + zVtxMap = getForTsOrRun(fullPathCalibVtx, timestamp, runnumber); + } + } + + using MyCollisions = soa::Join; + using Colls = soa::Join; + using CollsGen = soa::Join; + using MCColls = soa::Join; + using CollsGenSgn = soa::SmallGroups>; + using MyPIDTracks = soa::Join; + using MyLabeledTracks = soa::Join; + using MyFiltLabeledTracks = soa::Filtered; + using MyLabeledPIDTracks = soa::Join; + + template + std::unique_ptr setFuncPars(T const& vecPars) + { + std::unique_ptr fCalibFunc(new TF1("fCalibFunc", "pol2", -1., 1.)); + if (vecPars.size() >= 1) { + for (typename T::size_type i = 0; i < vecPars.size(); i++) { + fCalibFunc->SetParameter(i, vecPars[i]); + } + } + return fCalibFunc; + } + + template + bool isPID(const P& mcParticle) + { + static_assert(pidSgn == 0 || pidSgn == 1); + static_assert(id > 0 || id < Npart); + constexpr int kIdx = id + pidSgn * Npart; + return mcParticle.pdgCode() == PidSgn[kIdx]; + } + + template + bool selTPCtrack(T const& track) + { + return (track.hasTPC() && + track.passedTPCCrossedRowsOverNCls() && + track.passedTPCChi2NDF() && + track.passedTPCNCls() && + track.passedTPCCrossedRows() && + track.passedTPCRefit()); + } + + template + bool selTOFPi(T const& track) + { + if (track.p() < trkSelOpt.cfgMomSelPiTOF) { + if (track.hasTOF()) { + if (std::abs(track.tpcNSigmaPi()) < v0SelOpt.cfgNsigmaPiTPC && std::abs(track.tofNSigmaPi()) < v0SelOpt.cfgNsigmaPiTOF) { + return true; + } + } + } + return false; + } + + template + void fillNsigma(T const& tracks, const C& collision) + { + const float mult = getMult(collision); + const float flat = fillFlat(collision); + for (const auto& track : tracks) { + checkNsigma(track, mult, flat); + } + } + + template + void filldEdx(T const& tracks, V const& v0s, const C& collision, aod::BCsWithTimestamps const& /*bcs*/) + { + auto bc = collision.template bc_as(); + auto magField = (ccdbConf.cfgMagField == 0) ? getMagneticField(bc.timestamp()) : ccdbConf.cfgMagField; + + if (applyCalibGain || applyCalibVtx) { + int currentRun = bc.runNumber(); + if (runNumber != currentRun) { + initCCDB(bc); + runNumber = currentRun; + } + } + + const float mult = getMult(collision); + const float flat = fillFlat(collision); + flatchrg.fill(HIST("Events/hFlatVsMultEst"), flat, mult); + + for (const auto& track : tracks) { + float dEdx = track.tpcSignal(); + bool posP = (track.sign() * track.tpcInnerParam() > 0) ? true : false; + if (cfgFillQAHisto) { + fillTrackQA(track); + } + if (!isGoodTrack(track, magField)) { + continue; + } + if (cfgFillQAHisto) { + fillTrackQA(track); + filldEdxQA(track, collision); + if (posP) { + filldEdxQA(track, collision); + } else { + filldEdxQA(track, collision); + } + } + if (applyCalibDeDx) { + dEdx *= (50.0 / getCalibration(true, track.eta())); + } + if (cfgFillQAHisto) { + filldEdxQA(track, collision); + if (posP) { + filldEdxQA(track, collision); + } else { + filldEdxQA(track, collision); + } + } + + // PID TPC dEdx + flatchrg.fill(HIST(kPrefix) + HIST(kCharge[kAll]) + HIST("hdEdx"), track.eta(), mult, flat, track.tpcInnerParam(), dEdx); + if (posP) { + flatchrg.fill(HIST(kPrefix) + HIST(kCharge[kPos]) + HIST("hdEdx"), track.eta(), mult, flat, track.tpcInnerParam(), dEdx); + } else { + flatchrg.fill(HIST(kPrefix) + HIST(kCharge[kNeg]) + HIST("hdEdx"), track.eta(), mult, flat, track.tpcInnerParam(), dEdx); + } + + // TOF pions + if (track.hasTOF() && track.beta() > 1.) { + if (selTOFPi(track)) { + flatchrg.fill(HIST(kPrefixCleanTof) + HIST(kCharge[kAll]) + HIST("hPiTof"), track.eta(), mult, flat, track.tpcInnerParam(), dEdx); + if (posP) { + flatchrg.fill(HIST(kPrefixCleanTof) + HIST(kCharge[kPos]) + HIST("hPiTof"), track.eta(), mult, flat, track.tpcInnerParam(), dEdx); + } else { + flatchrg.fill(HIST(kPrefixCleanTof) + HIST(kCharge[kNeg]) + HIST("hPiTof"), track.eta(), mult, flat, track.tpcInnerParam(), dEdx); + } + } + } + } + + // V0s + for (const auto& v0 : v0s) { + if (!isGoodV0Track(v0, tracks)) { + continue; + } + + const auto& posTrack = v0.template posTrack_as(); + const auto& negTrack = v0.template negTrack_as(); + float dEdxPos = posTrack.tpcSignal(); + float dEdxNeg = negTrack.tpcSignal(); + + if (applyCalibDeDx) { + dEdxPos *= (50.0 / getCalibration(true, posTrack.eta())); + dEdxNeg *= (50.0 / getCalibration(true, negTrack.eta())); + } + + if (selectTypeV0s(v0, posTrack, negTrack) == kGa) { // Gamma selection + flatchrg.fill(HIST(kPrefixCleanV0) + HIST(kCharge[kPos]) + HIST("hEV0"), posTrack.eta(), mult, flat, posTrack.sign() * posTrack.tpcInnerParam(), dEdxPos); + flatchrg.fill(HIST(kPrefixCleanV0) + HIST(kCharge[kNeg]) + HIST("hEV0"), negTrack.eta(), mult, flat, negTrack.sign() * negTrack.tpcInnerParam(), dEdxNeg); + } + if (selectTypeV0s(v0, posTrack, negTrack) == kKz) { // K0S -> pi + pi + flatchrg.fill(HIST(kPrefixCleanV0) + HIST(kCharge[kPos]) + HIST("hPiV0"), posTrack.eta(), mult, flat, posTrack.sign() * posTrack.tpcInnerParam(), dEdxPos); + flatchrg.fill(HIST(kPrefixCleanV0) + HIST(kCharge[kNeg]) + HIST("hPiV0"), negTrack.eta(), mult, flat, negTrack.sign() * negTrack.tpcInnerParam(), dEdxNeg); + } + if (selectTypeV0s(v0, posTrack, negTrack) == kLam) { // L -> p + pi- + flatchrg.fill(HIST(kPrefixCleanV0) + HIST(kCharge[kPos]) + HIST("hPV0"), posTrack.eta(), mult, flat, posTrack.sign() * posTrack.tpcInnerParam(), dEdxPos); + flatchrg.fill(HIST(kPrefixCleanV0) + HIST(kCharge[kNeg]) + HIST("hPV0"), negTrack.eta(), mult, flat, negTrack.sign() * negTrack.tpcInnerParam(), dEdxNeg); + } + if (selectTypeV0s(v0, posTrack, negTrack) == kaLam) { // L -> p + pi- + flatchrg.fill(HIST(kPrefixCleanV0) + HIST(kCharge[kPos]) + HIST("hPiV0"), posTrack.eta(), mult, flat, posTrack.sign() * posTrack.tpcInnerParam(), dEdxPos); + flatchrg.fill(HIST(kPrefixCleanV0) + HIST(kCharge[kNeg]) + HIST("hPiV0"), negTrack.eta(), mult, flat, negTrack.sign() * negTrack.tpcInnerParam(), dEdxNeg); + } + } + } + + float getCalibration(bool isMIP, const float& eta) + { + double valCalib = 999.0; + if (eta < 0.) { + if (isMIP) { + valCalib = fDeDxVsEtaNeg->Eval(eta); + } else { + valCalib = fEDeDxVsEtaNeg->Eval(eta); + } + } else { + if (isMIP) { + valCalib = fDeDxVsEtaPos->Eval(eta); + } else { + valCalib = fEDeDxVsEtaPos->Eval(eta); + } + } + return valCalib; + } + + bool isChrgParticle(int code) + { + auto p = pdg->GetParticle(code); + auto charge = 0.; + if (p != nullptr) { + charge = p->Charge(); + } + return std::abs(charge) >= kMinCharge; + } + + template + bool phiCut(T const& track, float mag, TF1* fphiCutLow, TF1* fphiCutHigh) + { + if (track.pt() < trkSelOpt.cfgPhiCutPtMin) + return true; + // cut to remove tracks at TPC boundaries + double phimodn = track.phi(); + if (mag < 0) // for negative polarity field + phimodn = o2::constants::math::TwoPI - phimodn; + if (track.sign() < 0) // for negative charge + phimodn = o2::constants::math::TwoPI - phimodn; + if (phimodn < 0) + LOGF(warning, "phi < 0: %g", phimodn); + + phimodn += o2::constants::math::PI / 18.0; // to center gap in the middle + phimodn = std::fmod(phimodn, o2::constants::math::PI / 9.0); + + if (cfgFillQAHisto) { + flatchrg.fill(HIST("Tracks/preSel/hPtPhi"), track.pt(), phimodn); + } + if (phimodn < fphiCutHigh->Eval(track.pt()) && phimodn > fphiCutLow->Eval(track.pt())) { + return false; + } + if (cfgFillQAHisto) { + flatchrg.fill(HIST("Tracks/postSel/hPtPhi"), track.pt(), phimodn); + } + return true; + } + + template + bool isGoodTrack(T const& track, int const magfield) + { + if (std::abs(track.eta()) > trkSelOpt.cfgTrkEtaMax) { + return false; + } + flatchrg.fill(HIST("Tracks/hTrkSel"), trkSelEta); + if (track.pt() < trkSelOpt.cfgTrkPtMin) { + return false; + } + flatchrg.fill(HIST("Tracks/hTrkSel"), trkSelPt); + if (!isDCAxyCut(track)) { + return false; + } + flatchrg.fill(HIST("Tracks/hTrkSel"), trkSelDCA); + if (track.tpcNClsCrossedRows() < minNCrossedRowsTPC) { + return false; + } + flatchrg.fill(HIST("Tracks/hTrkSel"), trkSelNCls); + if (trkSelOpt.cfgUseExtraTrkCut && !phiCut(track, magfield, fPhiCutLow, fPhiCutHigh)) { + return false; + } + flatchrg.fill(HIST("Tracks/hTrkSel"), trkSelTPCBndr); + return true; + } + + template + int selectTypeV0s(T1 const& v0, T2 const& postrk, T2 const& negtrk) + { + // Gamma selection + if (v0.mGamma() < v0SelOpt.cfgdmassG) { + if (postrk.tpcInnerParam() < v0SelOpt.cfgV0DaughterTpcMomMax && postrk.hasTPC() && std::abs(postrk.tpcNSigmaEl()) < v0SelOpt.cfgNsigmaElTPC) { + if (postrk.hasTOF() && std::abs(postrk.tofNSigmaEl()) < v0SelOpt.cfgNsigmaElTOF) { + return kGa; + } + } + if (negtrk.tpcInnerParam() < v0SelOpt.cfgV0DaughterTpcMomMax && negtrk.hasTPC() && std::abs(negtrk.tpcNSigmaEl()) < v0SelOpt.cfgNsigmaElTPC) { + if (negtrk.hasTOF() && std::abs(negtrk.tofNSigmaEl()) < v0SelOpt.cfgNsigmaElTOF) { + return kGa; + } + } + } + // K0S selection, K0S -> pi + pi + if (std::abs(v0.mK0Short() - o2::constants::physics::MassK0Short) < v0SelOpt.cfgdmassK) { + if (postrk.tpcInnerParam() < v0SelOpt.cfgV0DaughterTpcMomMax && postrk.hasTPC() && std::abs(postrk.tpcNSigmaPi()) < v0SelOpt.cfgNsigmaPiTPC) { + if (postrk.hasTOF() && std::abs(postrk.tofNSigmaPi()) < v0SelOpt.cfgNsigmaElTOF) { + return kKz; + } + } + if (negtrk.tpcInnerParam() < v0SelOpt.cfgV0DaughterTpcMomMax && negtrk.hasTPC() && std::abs(negtrk.tpcNSigmaPi()) < v0SelOpt.cfgNsigmaPiTPC) { + if (negtrk.hasTOF() && std::abs(negtrk.tofNSigmaPi()) < v0SelOpt.cfgNsigmaPiTOF) { + return kKz; + } + } + } + // Lambda selection, L -> p + pi- + if (std::abs(v0.mLambda() - o2::constants::physics::MassLambda0) < v0SelOpt.cfgdmassL) { + if (postrk.tpcInnerParam() < v0SelOpt.cfgV0DaughterTpcMomMax && postrk.hasTPC() && std::abs(postrk.tpcNSigmaPr()) < v0SelOpt.cfgNsigmaPrTPC && negtrk.hasTPC() && std::abs(negtrk.tpcNSigmaPi()) < v0SelOpt.cfgNsigmaPiTPC) { + if (postrk.hasTOF() && std::abs(postrk.tofNSigmaPr()) < v0SelOpt.cfgNsigmaPrTOF && negtrk.hasTOF() && std::abs(negtrk.tofNSigmaPi()) < v0SelOpt.cfgNsigmaPiTOF) { + return kLam; + } + } + } + // antiLambda -> pbar + pi+ + if (std::abs(v0.mAntiLambda() - o2::constants::physics::MassLambda0) < v0SelOpt.cfgdmassL) { + if (postrk.tpcInnerParam() < v0SelOpt.cfgV0DaughterTpcMomMax && postrk.hasTPC() && std::abs(postrk.tpcNSigmaPi()) < v0SelOpt.cfgNsigmaPiTPC && negtrk.hasTPC() && std::abs(negtrk.tpcNSigmaPr()) < v0SelOpt.cfgNsigmaPrTPC) { + if (postrk.hasTOF() && std::abs(postrk.tofNSigmaPi()) < v0SelOpt.cfgNsigmaPiTOF && negtrk.hasTOF() && std::abs(negtrk.tofNSigmaPr()) < v0SelOpt.cfgNsigmaPrTOF) { + return kaLam; + } + } + } + return kNaN; + } + + template + bool isGoodV0Track(T1 const& v0, T2 const& /*track*/) + { + const auto& posTrack = v0.template posTrack_as(); + const auto& negTrack = v0.template negTrack_as(); + + if (std::abs(posTrack.eta()) > v0SelOpt.cfgV0etamax || std::abs(negTrack.eta()) > v0SelOpt.cfgV0etamax) { + return false; + } + if (posTrack.tpcNClsFound() < v0SelOpt.cfgTPCnClsmin || negTrack.tpcNClsFound() < v0SelOpt.cfgTPCnClsmin) { + return false; + } + if (posTrack.sign() * negTrack.sign() > 0) { // reject same sign pair + return false; + } + if (v0.dcaV0daughters() > v0SelOpt.cfgDCAv0daughter) { + return false; + } + if (v0.v0cosPA() < v0SelOpt.cfgv0cospa) { + return false; + } + if (v0.v0radius() < v0SelOpt.cfgv0Rmin || v0.v0radius() > v0SelOpt.cfgv0Rmax) { + return false; + } + if (std::abs(v0.dcapostopv()) < v0SelOpt.cfgDCAposToPV || std::abs(v0.dcanegtopv()) < v0SelOpt.cfgDCAposToPV) { + return false; + } + if constexpr (fillHist) { + flatchrg.fill(HIST("Tracks/V0qa/hV0Pt"), v0.pt()); + flatchrg.fill(HIST("Tracks/V0qa/hV0ArmPod"), v0.alpha(), v0.qtarm()); + } + return true; + } + + template + void checkNsigma(const T& track, const float mult, const float flat) + { + if (std::abs(track.rapidity(o2::track::PID::getMass(pid))) > trkSelOpt.cfgRapMax) { + return; + } + + float valTPCnsigma = -999, valTOFnsigma = -999; + switch (pid) { + case o2::track::PID::Pion: + valTPCnsigma = track.tpcNSigmaPi(); + valTOFnsigma = track.tofNSigmaPi(); + break; + case o2::track::PID::Kaon: + valTPCnsigma = track.tpcNSigmaKa(); + valTOFnsigma = track.tofNSigmaKa(); + break; + case o2::track::PID::Proton: + valTPCnsigma = track.tpcNSigmaPr(); + valTOFnsigma = track.tofNSigmaPr(); + break; + case o2::track::PID::Electron: + valTPCnsigma = track.tpcNSigmaEl(); + valTOFnsigma = track.tofNSigmaEl(); + break; + case o2::track::PID::Muon: + valTPCnsigma = track.tpcNSigmaMu(); + valTOFnsigma = track.tofNSigmaMu(); + break; + default: + valTPCnsigma = -999, valTOFnsigma = -999; + break; + } + + if (track.sign() > 0) { + hPtNsigmaTPC[pid]->Fill(track.pt(), valTPCnsigma); + if (cfgStoreThnSparse) { + hThPtNsigmaTPC[pid]->Fill(track.pt(), valTPCnsigma, mult, flat); + } + } else { + hPtNsigmaTPC[pid + Npart]->Fill(track.pt(), valTPCnsigma); + if (cfgStoreThnSparse) { + hThPtNsigmaTPC[pid + Npart]->Fill(track.pt(), valTPCnsigma, mult, flat); + } + } + if (!track.hasTOF()) { + return; + } + if (track.sign() > 0) { + hPtNsigmaTOF[pid]->Fill(track.pt(), valTOFnsigma); + hPtNsigmaTPCTOF[pid]->Fill(valTPCnsigma, valTOFnsigma); + } else { + hPtNsigmaTOF[pid + Npart]->Fill(track.pt(), valTOFnsigma); + hPtNsigmaTPCTOF[pid + Npart]->Fill(valTPCnsigma, valTOFnsigma); + } + } + + template + inline void fillTrackQA(T const& track) + { + if constexpr (fillHist) { + flatchrg.fill(HIST(kPrefix) + HIST(kQA) + HIST("hPt"), track.pt()); + flatchrg.fill(HIST(kPrefix) + HIST(kQA) + HIST("hPhi"), track.phi()); + flatchrg.fill(HIST(kPrefix) + HIST(kQA) + HIST("hEta"), track.eta()); + flatchrg.fill(HIST(kPrefix) + HIST(kQA) + HIST("hDCAXYvsPt"), track.pt(), track.dcaXY()); + flatchrg.fill(HIST(kPrefix) + HIST(kQA) + HIST("hDCAZvsPt"), track.pt(), track.dcaZ()); + + if (track.hasTPC() && track.hasITS()) { + int nFindable = track.tpcNClsFindable(); + int nMinusFound = track.tpcNClsFindableMinusFound(); + int nCluster = nFindable - nMinusFound; + flatchrg.fill(HIST(kPrefix) + HIST(kQA) + HIST("hTPCCluster"), nCluster); + flatchrg.fill(HIST(kPrefix) + HIST(kQA) + HIST("hShTpcClvsPt"), track.pt(), track.tpcFractionSharedCls()); + flatchrg.fill(HIST(kPrefix) + HIST(kQA) + HIST("hCrossTPCvsPt"), track.pt(), track.tpcNClsFound()); + flatchrg.fill(HIST(kPrefix) + HIST(kQA) + HIST("tpcNClsShared"), track.tpcNClsShared()); + flatchrg.fill(HIST(kPrefix) + HIST(kQA) + HIST("tpcCrossedRows"), track.tpcNClsCrossedRows()); + flatchrg.fill(HIST(kPrefix) + HIST(kQA) + HIST("tpcCrossedRowsOverFindableCls"), track.tpcCrossedRowsOverFindableCls()); + flatchrg.fill(HIST(kPrefix) + HIST(kQA) + HIST("hChi2ITSTrkSegment"), track.itsChi2NCl()); + flatchrg.fill(HIST(kPrefix) + HIST(kQA) + HIST("itsNCls"), track.itsNCls()); + } + + flatchrg.fill(HIST(kPrefix) + HIST(kQA) + HIST("hTOFPvsBeta"), track.tpcInnerParam(), track.beta()); + if (track.beta() > trkSelOpt.cfgTOFBetaPion && track.beta() < trkSelOpt.cfgTOFBetaPion + 0.05) { // TOF pions + flatchrg.fill(HIST(kPrefix) + HIST(kQA) + HIST("hTOFpi"), track.eta(), track.tpcInnerParam(), track.tpcSignal()); + } + + if (std::abs(track.eta()) < trkSelOpt.cfgTrkEtaMax) { + if (isDCAxyWoCut(track)) { + flatchrg.fill(HIST(kPrefix) + HIST(kQA) + HIST("hPtVsWOcutDCA"), track.pt(), track.dcaXY()); + } + } + } + + flatchrg.fill(HIST(kPrefix) + HIST(kStatus[ft]) + HIST("hPVsPtEta"), track.tpcInnerParam(), track.pt(), track.eta()); + } + + template + inline void filldEdxQA(T const& track, C const& collision) + { + const float mult = getMult(collision); + const float flat = fillFlat(collision); + float dEdx = track.tpcSignal(); + if constexpr (fillHist) { + if (track.tpcInnerParam() >= trkSelOpt.cfgMomMIPMin && track.tpcInnerParam() <= trkSelOpt.cfgMomMIPMax) { + if (dEdx > trkSelOpt.cfgDeDxMIPMin && dEdx < trkSelOpt.cfgDeDxMIPMax) { // MIP pions + flatchrg.fill(HIST(kPrefix) + HIST(kStatCalib[ft]) + HIST(kCharge[chrg]) + HIST("hMIP"), mult, flat, track.eta(), dEdx); + flatchrg.fill(HIST(kPrefix) + HIST(kStatCalib[ft]) + HIST(kCharge[chrg]) + HIST("hMIPVsEta"), track.eta(), dEdx); + flatchrg.fill(HIST(kPrefix) + HIST(kStatCalib[ft]) + HIST(kCharge[chrg]) + HIST("pMIPVsEta"), track.eta(), dEdx); + flatchrg.fill(HIST(kPrefix) + HIST(kStatCalib[ft]) + HIST(kCharge[chrg]) + HIST("hMIPVsPhi"), track.phi(), dEdx); + flatchrg.fill(HIST(kPrefix) + HIST(kStatCalib[ft]) + HIST(kCharge[chrg]) + HIST("hMIPVsPhiVsEta"), track.phi(), dEdx, track.eta()); + flatchrg.fill(HIST(kPrefix) + HIST(kStatCalib[ft]) + HIST(kCharge[chrg]) + HIST("pMIPVsPhi"), track.phi(), dEdx); + } + if (dEdx > trkSelOpt.cfgDeDxMIPMax + 10. && dEdx < trkSelOpt.cfgDeDxMIPMax + 30.) { // Plateau electrons + if (std::abs(track.beta() - 1) < trkSelOpt.cfgBetaPlateuMax) { + flatchrg.fill(HIST(kPrefix) + HIST(kStatCalib[ft]) + HIST(kCharge[chrg]) + HIST("hPlateau"), mult, flat, track.eta(), dEdx); + flatchrg.fill(HIST(kPrefix) + HIST(kStatCalib[ft]) + HIST(kCharge[chrg]) + HIST("hPlateauVsEta"), track.eta(), dEdx); + flatchrg.fill(HIST(kPrefix) + HIST(kStatCalib[ft]) + HIST(kCharge[chrg]) + HIST("pPlateauVsEta"), track.eta(), dEdx); + flatchrg.fill(HIST(kPrefix) + HIST(kStatCalib[ft]) + HIST(kCharge[chrg]) + HIST("hPlateauVsPhi"), track.phi(), dEdx); + flatchrg.fill(HIST(kPrefix) + HIST(kStatCalib[ft]) + HIST(kCharge[chrg]) + HIST("hPlateauVsPhiVsEta"), track.phi(), dEdx, track.eta()); + flatchrg.fill(HIST(kPrefix) + HIST(kStatCalib[ft]) + HIST(kCharge[chrg]) + HIST("pPlateauVsPhi"), track.phi(), dEdx); + } + } + } + } + } + + template + bool isGoodEvent(C const& collision) + { + if constexpr (fillHist) { + flatchrg.fill(HIST("Events/hEvtSel"), evtSelAll); + } + if (evtSelOpt.cfgCustomTVX) { + if (!collision.selection_bit(aod::evsel::kIsTriggerTVX)) { + return false; + } + } else { + if (!collision.sel8()) { + return false; + } + } + if constexpr (fillHist) { + flatchrg.fill(HIST("Events/hEvtSel"), evtSelSel8); + } + if (evtSelOpt.cfgRemoveITSROFrameBorder && !collision.selection_bit(aod::evsel::kNoITSROFrameBorder)) { + return false; + } + if constexpr (fillHist) { + flatchrg.fill(HIST("Events/hEvtSel"), evtSelNoITSROFrameBorder); + } + if (evtSelOpt.cfgRemoveNoTimeFrameBorder && !collision.selection_bit(aod::evsel::kNoTimeFrameBorder)) { + return false; + } + if constexpr (fillHist) { + flatchrg.fill(HIST("Events/hEvtSel"), evtSelkNoTimeFrameBorder); + } + if (evtSelOpt.cfgRemoveNoSameBunchPileup && !collision.selection_bit(aod::evsel::kNoSameBunchPileup)) { + return false; + } + if constexpr (fillHist) { + flatchrg.fill(HIST("Events/hEvtSel"), evtSelkNoSameBunchPileup); + } + if (evtSelOpt.cfgRequireIsGoodZvtxFT0vsPV && !collision.selection_bit(aod::evsel::kIsGoodZvtxFT0vsPV)) { + return false; + } + if constexpr (fillHist) { + flatchrg.fill(HIST("Events/hEvtSel"), evtSelkIsGoodZvtxFT0vsPV); + } + if (evtSelOpt.cfgRequireIsVertexITSTPC && !collision.selection_bit(aod::evsel::kIsVertexITSTPC)) { + return false; + } + if constexpr (fillHist) { + flatchrg.fill(HIST("Events/hEvtSel"), evtSelkIsVertexITSTPC); + } + if (evtSelOpt.cfgRequirekIsVertexTOFmatched && !collision.selection_bit(aod::evsel::kIsVertexTOFmatched)) { + return false; + } + if constexpr (fillHist) { + flatchrg.fill(HIST("Events/hEvtSel"), evtSelkIsVertexTOFmatched); + } + if (std::abs(collision.posZ()) > evtSelOpt.cfgCutVtxZ) { + return false; + } + if constexpr (fillHist) { + flatchrg.fill(HIST("Events/hEvtSel"), evtSelVtxZ); + flatchrg.fill(HIST("Events/hVtxZ"), collision.posZ()); + } + if (evtSelOpt.cfgINELCut && !collision.isInelGt0()) { + return false; + } + if constexpr (fillHist) { + flatchrg.fill(HIST("Events/hEvtSel"), evtSelINELgt0); + } + return true; + } + + int getFV0IndexPhi(int i_ch) + { + int iRing = -1; + + if (i_ch >= kFV0IndexPhi[0] && i_ch < kFV0IndexPhi[1]) { + if (i_ch < kFV0IndexPhi[0] + 4) { + iRing = i_ch; + } else { + if (i_ch == kFV0IndexPhi[1] - 1) { + iRing = i_ch - 3; // 4; + } else if (i_ch == kFV0IndexPhi[1] - 2) { + iRing = i_ch - 1; // 5; + } else if (i_ch == kFV0IndexPhi[1] - 3) { + iRing = i_ch + 1; // 6; + } else if (i_ch == kFV0IndexPhi[1] - 4) { + iRing = i_ch + 3; // 7; + } + } + } else if (i_ch >= kFV0IndexPhi[1] && i_ch < kFV0IndexPhi[2]) { + if (i_ch < kFV0IndexPhi[2] - 4) { + iRing = i_ch; + } else { + if (i_ch == kFV0IndexPhi[2] - 1) { + iRing = i_ch - 3; // 12; + } else if (i_ch == kFV0IndexPhi[2] - 2) { + iRing = i_ch - 1; // 13; + } else if (i_ch == kFV0IndexPhi[2] - 3) { + iRing = i_ch + 1; // 14; + } else if (i_ch == kFV0IndexPhi[2] - 4) { + iRing = i_ch + 3; // 15; + } + } + } else if (i_ch >= kFV0IndexPhi[2] && i_ch < kFV0IndexPhi[3]) { + if (i_ch < kFV0IndexPhi[3] - 4) { + iRing = i_ch; + } else { + if (i_ch == kFV0IndexPhi[3] - 1) { + iRing = i_ch - 3; // 20; + } else if (i_ch == kFV0IndexPhi[3] - 2) { + iRing = i_ch - 1; // 21; + } else if (i_ch == kFV0IndexPhi[3] - 3) { + iRing = i_ch + 1; // 22; + } else if (i_ch == kFV0IndexPhi[3] - 4) { + iRing = i_ch + 3; // 23; + } + } + } else if (i_ch >= kFV0IndexPhi[3] && i_ch < kFV0IndexPhi[4]) { + if (i_ch < kFV0IndexPhi[3] + 4) { + iRing = i_ch; + } else { + if (i_ch == kFV0IndexPhi[4] - 5) { + iRing = i_ch - 3; // 28; + } else if (i_ch == kFV0IndexPhi[4] - 6) { + iRing = i_ch - 1; // 29; + } else if (i_ch == kFV0IndexPhi[4] - 7) { + iRing = i_ch + 1; // 30; + } else if (i_ch == kFV0IndexPhi[4] - 8) { + iRing = i_ch + 3; // 31; + } + } + } else if (i_ch == kFV0IndexPhi[4]) { + iRing = kFV0IndexPhi[4]; + } else if (i_ch == kFV0IndexPhi[4] + 8) { + iRing = i_ch - 7; // 33; + } else if (i_ch == kFV0IndexPhi[4] - 3) { + iRing = i_ch + 1; // 34; + } else if (i_ch == kFV0IndexPhi[4] + 5) { + iRing = i_ch - 6; // 35; + } else if (i_ch == kFV0IndexPhi[4] - 2) { + iRing = i_ch + 2; // 36; + } else if (i_ch == kFV0IndexPhi[4] + 6) { + iRing = i_ch - 5; // 37; + } else if (i_ch == kFV0IndexPhi[4] - 1) { + iRing = i_ch + 3; // 38; + } else if (i_ch == kFV0IndexPhi[4] + 7) { + iRing = i_ch - 4; // 39; + } else if (i_ch == kFV0IndexPhi[4] + 11) { + iRing = i_ch + 7; // 40; + } else if (i_ch == kFV0IndexPhi[4] + 3) { + iRing = i_ch + 2; // 41; + } else if (i_ch == kFV0IndexPhi[4] + 10) { + iRing = i_ch - 4; // 42; + } else if (i_ch == kFV0IndexPhi[4] + 1) { + iRing = i_ch + 5; // 43; + } else if (i_ch == kFV0IndexPhi[4] + 9) { + iRing = i_ch - 1; // 44; + } else if (i_ch == kFV0IndexPhi[4] + 1) { + iRing = i_ch + 8; // 45; + } else if (i_ch == kFV0IndexPhi[4] + 8) { + iRing = i_ch + 2; // 46; + } else if (i_ch == kFV0IndexPhi[4]) { + iRing = i_ch + 11; // 47; + } + return iRing; + } + + int getMagneticField(uint64_t timestamp) + { + o2::parameters::GRPMagField* grpmag = nullptr; + grpmag = ccdb->getForTimeStamp(ccdbConf.grpmagPath, timestamp); + if (!grpmag) { + return 0; + } + return grpmag->getNominalL3Field(); + } + + template + float getMult(C const& collision) + { + float val = -999.0; + switch (multEst) { + case MultE::kNoMult: + return val; + break; + case MultE::kMultFT0M: + return collision.centFT0M(); + break; + case MultE::kMultTPC: + if constexpr (!isMC) { + return collision.multTPC(); + } else { + LOG(fatal) << "No valid multiplicity estimator: " << multEst; + return val; + } + break; + default: + return collision.centFT0M(); + break; + } + } + + float getMultMC(MCColls::iterator const& collision) + { + return getMult(collision); + } + + template + float calcFlatenicity(std::array const& signals) + { + static_assert(S != 0); + + int entries = signals.size(); + float flat{-1}; + float mRho{0}; + for (int iCell = 0; iCell < entries; ++iCell) { + if (signals[iCell] > 0.) { + mRho += 1.0 * signals[iCell]; + } + } + // average activity per cell + mRho /= (1.0 * entries); + if (mRho <= 0) { + return -1; + } + // get sigma + float sRhoTmp{0}; + float sRho{0}; + for (int iCell = 0; iCell < entries; ++iCell) { + if (signals[iCell] > 0.) { + sRhoTmp += std::pow(1.0 * signals[iCell] - mRho, 2); + } + } + sRhoTmp /= (1.0 * entries * entries); + sRho = std::sqrt(sRhoTmp); + if (mRho > 0.) { + flat = sRho / mRho; + } else { + flat = -1; + } + return flat; + } + + template + float fillFlat(C const& collision) + { + rhoLatticeFV0.fill(0); + fv0AmplitudeWoCalib.fill(0); + bool isOkFV0OrA = false; + if (collision.has_foundFV0()) { + auto fv0 = collision.foundFV0(); + std::bitset<8> fV0Triggers = fv0.triggerMask(); + isOkFV0OrA = fV0Triggers[o2::fit::Triggers::bitA]; + if (isOkFV0OrA) { + for (std::size_t ich = 0; ich < fv0.channel().size(); ich++) { + float amplCh = fv0.amplitude()[ich]; + int chv0 = fv0.channel()[ich]; + int chv0phi = getFV0IndexPhi(chv0); + if constexpr (fillHist) { + flatchrg.fill(HIST("FV0/hFV0amp"), chv0, amplCh); + flatchrg.fill(HIST("FV0/pFV0amp"), chv0, amplCh); + if (applyCalibGain) { + flatchrg.fill(HIST("FV0/hFV0ampCorr"), chv0, amplCh / fv0AmplCorr[chv0]); + } + } + if (amplCh > 0.) { + if (applyCalibGain) { // equalize gain channel-by-channel + amplCh /= fv0AmplCorr[chv0]; + } + if (chv0phi > 0) { + fv0AmplitudeWoCalib[chv0phi] = amplCh; + if constexpr (fillHist) { + flatchrg.fill(HIST("FV0/hFV0AmplWCalib"), ich, fv0AmplitudeWoCalib[ich]); + } + if (chv0 < kInnerFV0) { + rhoLatticeFV0[chv0phi] += amplCh; + } else { // two channels per bin + rhoLatticeFV0[chv0phi] += amplCh / 2.; + } + if constexpr (fillHist) { + flatchrg.fill(HIST("FV0/hFV0AmplvsVtxzWoCalib"), collision.posZ(), rhoLatticeFV0[chv0phi]); + } + if (applyCalibVtx) { + rhoLatticeFV0[chv0phi] *= zVtxMap->GetBinContent(zVtxMap->GetXaxis()->FindBin(chv0phi), zVtxMap->GetYaxis()->FindBin(collision.posZ())); + if constexpr (fillHist) { + flatchrg.fill(HIST("FV0/hFV0AmplvsVtxzCalib"), collision.posZ(), rhoLatticeFV0[chv0phi]); + } + } + } + } + } + float flattenicityFV0 = calcFlatenicity(rhoLatticeFV0); + return 1. - flattenicityFV0; + } else { + return 9999; + } + } else { + return 9999; + } + } + + template + bool isDCAxyCut(T const& track) const + { + if (isCustomTracks.value) { + for (int i = 0; i < static_cast(TrackSelection::TrackCuts::kNCuts); i++) { + if (i == static_cast(TrackSelection::TrackCuts::kDCAxy)) { + continue; + } + if (!mTrackSelector.IsSelected(track, static_cast(i))) { + return false; + } + } + return (std::abs(track.dcaXY()) <= (maxDcaXYFactor.value * (0.0105f + 0.0350f / std::pow(track.pt(), 1.1f)))); + } + return track.isGlobalTrack(); + } + + template + bool isDCAxyWoCut(T const& track) const + { + if (isCustomTracks.value) { + for (int i = 0; i < static_cast(TrackSelection::TrackCuts::kNCuts); i++) { + if (i == static_cast(TrackSelection::TrackCuts::kDCAxy)) { + continue; + } + if (i == static_cast(TrackSelection::TrackCuts::kDCAz)) { + continue; + } + if (!mTrackSelector.IsSelected(track, static_cast(i))) { + return false; + } + } + return true; + } + return track.isGlobalTrackWoDCA(); + } + + Preslice perCol = aod::track::collisionId; + Preslice perColV0s = aod::v0::collisionId; + + template + void processData(C const& collisions, + MyPIDTracks const& tracks, + aod::V0Datas const& v0s, + aod::BCsWithTimestamps const& bcs) + { + for (const auto& collision : collisions) { + if (!isGoodEvent(collision)) { + continue; + } + auto tracksPerCollision = tracks.sliceBy(perCol, collision.globalIndex()); + auto v0sPerCollision = v0s.sliceBy(perColV0s, collision.globalIndex()); + v0sPerCollision.bindExternalIndices(&tracks); + filldEdx(tracksPerCollision, v0sPerCollision, collision, bcs); + if (cfgFillQAHisto) { + fillNsigma(tracksPerCollision, collision); + fillNsigma(tracksPerCollision, collision); + fillNsigma(tracksPerCollision, collision); + } + } + } + + void processFlat( + Colls const& collisions, + MyPIDTracks const& tracks, + aod::V0Datas const& v0s, + aod::BCsWithTimestamps const& bcs, + aod::FT0s const&, aod::FV0As const&) + { + processData(collisions, tracks, v0s, bcs); + } + PROCESS_SWITCH(FlattenictyPikp, processFlat, "process Flat data inclusive", true); + + template + float fillFlatMC(McPart const& mcparts) + { + int nFV0sectors = 0; + float minPhi = 0; + float maxPhi = 0; + float dPhi = 0; + + double etaMinFV0bins[kMaxRingsFV0] = {0.0}; + double etaMaxFV0bins[kMaxRingsFV0] = {0.0}; + for (int i = 0; i < kMaxRingsFV0; ++i) { + etaMaxFV0bins[i] = kEtaMaxFV0 - i * kDEtaFV0; + if (i < kMaxRingsFV0 - 1) { + etaMinFV0bins[i] = kEtaMaxFV0 - (i + 1) * kDEtaFV0; + } else { + etaMinFV0bins[i] = kEtaMinFV0; + } + } + + rhoLatticeFV0.fill(0); + std::vector vNch; + float nCharged{0}; + for (const auto& mcPart : mcparts) { + if (!isChrgParticle(mcPart.pdgCode())) { + continue; + } + auto etaMc = mcPart.eta(); + auto phiMc = mcPart.phi(); + + int isegment = 0; + for (int ieta = 0; ieta < kMaxRingsFV0; ieta++) { + + nFV0sectors = kNCellsFV0 / 6.; + if (ieta == kMaxRingsFV0 - 1) { + nFV0sectors = kNCellsFV0 / 3.; + } + + for (int iphi = 0; iphi < nFV0sectors; iphi++) { + + minPhi = iphi * TwoPI / nFV0sectors; + maxPhi = (iphi + 1) * TwoPI / nFV0sectors; + dPhi = std::abs(maxPhi - minPhi); + + if (etaMc >= etaMinFV0bins[ieta] && etaMc < etaMaxFV0bins[ieta] && phiMc >= minPhi && phiMc < maxPhi) { + rhoLatticeFV0[isegment] += 1. / std::abs(kDEtaFV0 * dPhi); + } + isegment++; + } + } + nCharged++; + } + + vNch.push_back(nCharged); + auto flatFV0 = calcFlatenicity(rhoLatticeFV0); + + if constexpr (fillHist) { + flatchrg.fill(HIST("ResponseGen"), vNch[0], 1. - flatFV0); + flatchrg.fill(HIST("h1flatencityFV0MCGen"), 1. - flatFV0); + } + vNch.clear(); + return 1. - flatFV0; + } + + template + void bookMcHist() + { + AxisSpec ptAxis{binOpt.axisPt, "#it{p}_{T} (GeV/#it{c})"}; + constexpr int kHistIdx = id + pidSgn * Npart; + auto kIdx = static_cast(id); + const std::string strID = Form("/%s/%s", (pidSgn == 0 && id < Npart) ? "pos" : "neg", Pid[kIdx]); + hPtEffRec[kHistIdx] = flatchrg.add("Tracks/hPtEffRec" + strID, " ; p_{T} (GeV/c)", HistType::kTH1D, {ptAxis}); + hPtEffGen[kHistIdx] = flatchrg.add("Tracks/hPtEffGen" + strID, " ; p_{T} (GeV/c)", HistType::kTH1D, {ptAxis}); + } + + template + void initEfficiency() + { + static_assert(pidSgn == 0 || pidSgn == 1); + static_assert(id > 0 || id < Npart); + constexpr int kIdx = id + pidSgn * Npart; + const TString partName = PidChrg[kIdx]; + THashList* lhash = new THashList(); + lhash->SetName(partName); + listEfficiency->Add(lhash); + + auto bookEff = [&](const TString eName, auto h) { + const TAxis* axis = h->GetXaxis(); + TString eTitle = h->GetTitle(); + eTitle.ReplaceAll("Numerator", "").Strip(TString::kBoth); + eTitle = Form("%s;%s;Efficiency", eTitle.Data(), axis->GetTitle()); + lhash->Add(new TEfficiency(eName, eTitle, axis->GetNbins(), axis->GetXbins()->GetArray())); + }; + + const int kHistIdx = id + pidSgn * Npart; + bookEff("hEffvsPt", hPtEffRec[kHistIdx]); + } + + template + void fillEfficiency() + { + static_assert(pidSgn == 0 || pidSgn == 1); + constexpr int kHistIdx = id + pidSgn * Npart; + const char* partName = PidChrg[kHistIdx]; + THashList* lhash = static_cast(listEfficiency->FindObject(partName)); + if (!lhash) { + LOG(warning) << "No efficiency object found for particle " << partName; + return; + } + + auto fillEff = [&](const TString eName, auto num, auto den) { + TEfficiency* eff = static_cast(lhash->FindObject(eName)); + if (!eff) { + LOG(warning) << "Cannot find TEfficiency " << eName; + return; + } + eff->SetTotalHistogram(*den, "f"); + eff->SetPassedHistogram(*num, "f"); + }; + fillEff("hEffvsPt", hPtEffRec[kHistIdx], hPtEffGen[kHistIdx]); + } + + template + void fillMCRecTrack(MyLabeledPIDTracks::iterator const& track, const float mult, const float flat) + { + static_assert(pidSgn == 0 || pidSgn == 1); + constexpr int kHistIdx = id + pidSgn * Npart; + const aod::McParticles::iterator& mcParticle = track.mcParticle(); + const CollsGen::iterator& collision = track.collision_as(); + + if (!isPID(mcParticle)) { + return; + } + flatchrg.fill(HIST("hPtOutNoEtaCut"), track.pt()); + if (std::abs(track.eta()) > trkSelOpt.cfgTrkEtaMax) { + return; + } + if (std::abs(mcParticle.y()) > trkSelOpt.cfgRapMax) { + return; + } + + if ((collision.has_mcCollision() && (mcParticle.mcCollisionId() != collision.mcCollisionId())) || !collision.has_mcCollision()) { + if (!mcParticle.isPhysicalPrimary()) { + if (mcParticle.getProcess() == kProcessIdWeak) { + hDCAxyBadCollWeak[kHistIdx]->Fill(track.pt(), track.dcaXY()); + } else { + hDCAxyBadCollMat[kHistIdx]->Fill(track.pt(), track.dcaXY()); + } + } else { + hDCAxyBadCollPrim[kHistIdx]->Fill(track.pt(), track.dcaXY()); + } + } + + if (!isDCAxyCut(track)) { + return; + } + flatchrg.fill(HIST("hPtVsDCAxyAll"), track.pt(), track.dcaXY()); + + if (selTPCtrack(track)) { + hPtEffRec[kHistIdx]->Fill(mcParticle.pt()); + } + + if (!mcParticle.isPhysicalPrimary()) { + if (mcParticle.getProcess() == kProcessIdWeak) { + hPtEffRecWeak[kHistIdx]->Fill(mult, flat, track.pt()); + hPtVsDCAxyWeak[kHistIdx]->Fill(track.pt(), track.dcaXY()); + flatchrg.fill(HIST("hPtVsDCAxyWeakAll"), track.pt(), track.dcaXY()); + } else { + hPtEffRecMat[kHistIdx]->Fill(mult, flat, track.pt()); + hPtVsDCAxyMat[kHistIdx]->Fill(track.pt(), track.dcaXY()); + flatchrg.fill(HIST("hPtVsDCAxyMatAll"), track.pt(), track.dcaXY()); + } + } else { + hPtEffRecPrim[kHistIdx]->Fill(mult, flat, track.pt()); + hPtVsDCAxyPrim[kHistIdx]->Fill(track.pt(), track.dcaXY()); + flatchrg.fill(HIST("hPtVsDCAxyPrimAll"), track.pt(), track.dcaXY()); + } + } + + template + void fillMCGen(aod::McParticles::iterator const& mcParticle, const float mult, const float flat) + { + static_assert(pidSgn == 0 || pidSgn == 1); + constexpr int kHistIdx = id + pidSgn * Npart; + + if (!isPID(mcParticle)) { + return; + } + + if constexpr (recoEvt) { + hPtGenRecEvt[kHistIdx]->Fill(mcParticle.pt()); + if (mcParticle.isPhysicalPrimary()) { + hPtGenPrimRecEvt[kHistIdx]->Fill(mcParticle.pt()); + } + return; + } + + if (!mcParticle.isPhysicalPrimary()) { + if (mcParticle.getProcess() == kProcessIdWeak) { + hPtEffGenWeak[kHistIdx]->Fill(mult, flat, mcParticle.pt()); + } else { + hPtEffGenMat[kHistIdx]->Fill(mult, flat, mcParticle.pt()); + } + } else { + hPtEffGenPrim[kHistIdx]->Fill(mult, flat, mcParticle.pt()); + hPtEffGen[kHistIdx]->Fill(mcParticle.pt()); + } + } + + void processSgnLoss(MCColls::iterator const& mcCollision, + CollsGenSgn const& collisions, + aod::FV0As const& /*fv0s*/, + aod::McParticles const& particles) + { + float flat; + float mult; + if (flatSelOpt.useFlatData) { + float flatRec = 999.0; + float multRec = 999.0; + for (const auto& collision : collisions) { + multRec = getMult(collision); + flatRec = fillFlat(collision); + } + flat = flatRec; + mult = multRec; + flatchrg.fill(HIST("hFlatMCGenRecColl"), flatRec); + } else { + float flatGen = fillFlatMC(particles); + flat = flatGen; + flatchrg.fill(HIST("hFlatMCGen"), flatGen); + float multGen = getMultMC(mcCollision); + mult = multGen; + } + + // Evt loss den + flatchrg.fill(HIST("hEvtMcGen"), 0.5); + if (std::abs(mcCollision.posZ()) > evtSelOpt.cfgCutVtxZ) { + return; + } + flatchrg.fill(HIST("hEvtMcGen"), 1.5); + + bool isINELgt0mc = false; + if (pwglf::isINELgtNmc(particles, 0, pdg)) { + isINELgt0mc = true; + flatchrg.fill(HIST("hEvtMcGen"), 2.5); + flatchrg.fill(HIST("hFlatGenINELgt0"), flat); + } + + // Sgn loss den + for (const auto& particle : particles) { + if (!particle.isPhysicalPrimary()) { + continue; + } + if (std::abs(particle.y()) > trkSelOpt.cfgRapMax) { + continue; + } + static_for<0, 5>([&](auto i) { + constexpr int kIdx = i.value; + if (particle.pdgCode() == PidSgn[kIdx]) { + flatchrg.fill(HIST(kPrefix) + HIST(kSpecies[kIdx]) + HIST(kPtGenPrimSgn), mult, flat, particle.pt()); + if (isINELgt0mc) { + flatchrg.fill(HIST(kPrefix) + HIST(kSpecies[kIdx]) + HIST(kPtGenPrimSgnINEL), mult, flat, particle.pt()); + } + } + }); + } + + int nRecCollINEL = 0; + int nRecCollINELgt0 = 0; + for (const auto& collision : collisions) { + // Evt split num + flatchrg.fill(HIST("hEvtMCRec"), 0.5); + if (!isGoodEvent(collision)) { + continue; + } + flatchrg.fill(HIST("hEvtMCRec"), 1.5); + + nRecCollINEL++; + + if (collision.isInelGt0() && isINELgt0mc) { + flatchrg.fill(HIST("hEvtMCRec"), 2.5); + nRecCollINELgt0++; + } + // Sgn split num + for (const auto& particle : particles) { + if (!particle.isPhysicalPrimary()) { + continue; + } + if (std::abs(particle.y()) > trkSelOpt.cfgRapMax) { + continue; + } + static_for<0, 5>([&](auto i) { + constexpr int kIdx = i.value; + if (particle.pdgCode() == PidSgn[kIdx]) { + flatchrg.fill(HIST(kPrefix) + HIST(kSpecies[kIdx]) + HIST(kPtRecCollPrimSgn), mult, flat, particle.pt()); + if (nRecCollINELgt0) { + flatchrg.fill(HIST(kPrefix) + HIST(kSpecies[kIdx]) + HIST(kPtRecCollPrimSgnINEL), mult, flat, particle.pt()); + } + } + }); + } + } + + if (nRecCollINEL < 1) { + return; + } + // Evt loss num + flatchrg.fill(HIST("hEvtMcGenRecColl"), 0.5); + if (nRecCollINELgt0 > 0) { + flatchrg.fill(HIST("hEvtMcGenRecColl"), 1.5); + } + + // Sgn loss num + for (const auto& particle : particles) { + if (!particle.isPhysicalPrimary()) { + continue; + } + if (std::abs(particle.y()) > trkSelOpt.cfgRapMax) { + continue; + } + static_for<0, 5>([&](auto i) { + constexpr int kIdx = i.value; + if (particle.pdgCode() == PidSgn[kIdx]) { + flatchrg.fill(HIST(kPrefix) + HIST(kSpecies[kIdx]) + HIST(kPtGenRecCollPrimSgn), mult, flat, particle.pt()); + if (nRecCollINELgt0) { + flatchrg.fill(HIST(kPrefix) + HIST(kSpecies[kIdx]) + HIST(kPtGenRecCollPrimSgnINEL), mult, flat, particle.pt()); + } + } + }); + } + } + PROCESS_SWITCH(FlattenictyPikp, processSgnLoss, "process to calcuate signal/event lossses", false); + + // using Particles = soa::Filtered; + // expressions::Filter primaries = (aod::mcparticle::flags & (uint8_t)o2::aod::mcparticle::enums::PhysicalPrimary) == (uint8_t)o2::aod::mcparticle::enums::PhysicalPrimary; + + void processMCclosure(Colls::iterator const& collision, + MyPIDTracks const& tracks, + MyLabeledTracks const& mcTrackLabels, + aod::McParticles const& particles, + // Particles const& particles, + aod::FV0As const& /*fv0s*/, + aod::BCsWithTimestamps const& /*bcs*/) + { + const float multRec = getMult(collision); + const float flatRec = fillFlat(collision); + for (const auto& track : tracks) { + if (!track.has_collision()) { + continue; + } + const auto& coll = track.collision_as(); + auto bc = coll.template bc_as(); + auto magField = (ccdbConf.cfgMagField == 0) ? getMagneticField(bc.timestamp()) : ccdbConf.cfgMagField; + + if (!isGoodEvent(coll)) { + continue; + } + if (!isGoodTrack(track, magField)) { + continue; + } + if (!isDCAxyCut(track)) { + continue; + } + const auto& mcLabel = mcTrackLabels.iteratorAt(track.globalIndex()); + const auto& mcParticle = particles.iteratorAt(mcLabel.mcParticleId()); + + static_for<0, 4>([&](auto i) { + constexpr int kIdx = i.value; + if ((std::abs(o2::aod::pidutils::tpcNSigma(track)) < trkSelOpt.cfgNsigmaMax) && std::abs(track.rapidity(o2::track::PID::getMass(kIdx))) <= trkSelOpt.cfgRapMax) { + if (std::fabs(mcParticle.pdgCode()) == PDGs[kIdx]) { + flatchrg.fill(HIST(kPrefix) + HIST(kSpeciesAll[kIdx]) + HIST(kPtMCclosurePrim), multRec, flatRec, track.pt()); + } + } + }); + } + } + PROCESS_SWITCH(FlattenictyPikp, processMCclosure, "process MC closure test", false); + + Preslice perCollTrk = aod::track::collisionId; + PresliceUnsorted perCollMcLabel = aod::mccollisionlabel::mcCollisionId; + Preslice perCollMcPart = aod::mcparticle::mcCollisionId; + + void processMC(MCColls const& mcCollisions, + CollsGen const& collisions, + MyLabeledPIDTracks const& tracks, + aod::McParticles const& mcparticles) + { + flatchrg.fill(HIST("hEvtGenRec"), 1.f, mcCollisions.size()); + flatchrg.fill(HIST("hEvtGenRec"), 2.f, collisions.size()); + + for (const auto& mcCollision : mcCollisions) { + if (mcCollision.isInelGt0()) { + flatchrg.fill(HIST("hEvtGenRec"), 3.f); + } + flatchrg.fill(HIST("hEvtMcGenColls"), 1); + const auto groupedColls = collisions.sliceBy(perCollMcLabel, mcCollision.globalIndex()); + const auto groupedParts = mcparticles.sliceBy(perCollMcPart, mcCollision.globalIndex()); + const float flatMC = fillFlatMC(groupedParts); + const float multMC = getMultMC(mcCollision); + if (groupedColls.size() < 1) { // if MC events have no rec collisions + continue; + } + flatchrg.fill(HIST("hEvtMcGenColls"), 2); + for (const auto& collision : groupedColls) { + flatchrg.fill(HIST("hEvtMcGenColls"), 3); + if (!isGoodEvent(collision)) { + continue; + } + flatchrg.fill(HIST("hEvtMcGenColls"), 4); + const auto groupedTrks = tracks.sliceBy(perCollTrk, collision.globalIndex()); + for (const auto& track : groupedTrks) { + if (!isDCAxyWoCut(track)) { + continue; + } + if (!track.has_mcParticle()) { + flatchrg.fill(HIST("PtOutFakes"), track.pt()); + continue; + } + const auto& mcParticle = track.mcParticle(); + if (std::abs(mcParticle.y()) > trkSelOpt.cfgRapMax) { + continue; + } + if (!track.has_collision()) { + continue; + } + static_for<0, 1>([&](auto pidSgn) { + fillMCRecTrack(track, multMC, flatMC); + fillMCRecTrack(track, multMC, flatMC); + fillMCRecTrack(track, multMC, flatMC); + }); + } + + if (std::abs(mcCollision.posZ()) > evtSelOpt.cfgCutVtxZ) { + continue; + } + + if (evtSelOpt.cfgINELCut.value) { + if (!o2::pwglf::isINELgt0mc(groupedParts, pdg)) { + continue; + } + } + + for (const auto& particle : groupedParts) { + if (std::abs(particle.y()) > trkSelOpt.cfgRapMax) { + continue; + } + static_for<0, 1>([&](auto pidSgn) { + fillMCGen(particle, multMC, flatMC); + fillMCGen(particle, multMC, flatMC); + fillMCGen(particle, multMC, flatMC); + }); + } + } // reco collisions + + if (evtSelOpt.cfgINELCut.value) { + if (!o2::pwglf::isINELgt0mc(groupedParts, pdg)) { + continue; + } + } + + for (const auto& mcParticle : groupedParts) { + if (std::abs(mcParticle.y()) > trkSelOpt.cfgRapMax) { + continue; + } + static_for<0, 1>([&](auto pidSgn) { + fillMCGen(mcParticle, multMC, flatMC); + fillMCGen(mcParticle, multMC, flatMC); + fillMCGen(mcParticle, multMC, flatMC); + }); + } + + static_for<0, 1>([&](auto pidSgn) { + fillEfficiency(); + fillEfficiency(); + fillEfficiency(); + }); + + } // gen collisions + } + PROCESS_SWITCH(FlattenictyPikp, processMC, "process MC", false); + + template + ObjType* getForTsOrRun(std::string const& fullPath, int64_t timestamp, int runNumber) + { + if (cfgUseCcdbForRun) { + return ccdb->getForRun(fullPath, runNumber); + } else { + return ccdb->getForTimeStamp(fullPath, timestamp); + } + } +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{adaptAnalysisTask(cfgc)}; +}