diff --git a/PWGLF/Tasks/Nuspex/nucleitpcpbpb.cxx b/PWGLF/Tasks/Nuspex/nucleitpcpbpb.cxx deleted file mode 100644 index 52dd21fdeef..00000000000 --- a/PWGLF/Tasks/Nuspex/nucleitpcpbpb.cxx +++ /dev/null @@ -1,1416 +0,0 @@ -// Copyright 2019-2020 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 nucleitpcpbpb.cxx -/// \brief nuclei analysis -/// \note under work -/// -/// \author Jaideep Tanwar , Panjab University - -#include "Common/CCDB/EventSelectionParams.h" -#include "Common/Core/PID/PIDTOF.h" -#include "Common/Core/trackUtilities.h" -#include "Common/DataModel/Centrality.h" -#include "Common/DataModel/EventSelection.h" -#include "Common/DataModel/PIDResponseITS.h" -#include "Common/DataModel/PIDResponseTOF.h" -#include "Common/DataModel/PIDResponseTPC.h" -#include "Common/DataModel/TrackSelectionTables.h" - -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#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; -using CollisionsFull = soa::Join; -using TracksFull = soa::Join; -using CollisionsFullMC = soa::Join; -//--------------------------------------------------------------------------------------------------------------------------------- -namespace -{ -static const int nParticles = 6; -static const std::vector particleNames{"pion", "proton", "deuteron", "triton", "helion", "alpha"}; -static const std::vector correctedparticleNames{"helion", "antihelion", "alpha", "antialpha"}; -static const std::vector particlePdgCodes{211, 2212, o2::constants::physics::kDeuteron, o2::constants::physics::kTriton, o2::constants::physics::kHelium3, o2::constants::physics::kAlpha}; -static const std::vector particleMasses{o2::constants::physics::MassPionCharged, o2::constants::physics::MassProton, o2::constants::physics::MassDeuteron, o2::constants::physics::MassTriton, o2::constants::physics::MassHelium3, o2::constants::physics::MassAlpha}; -static const std::vector particleCharge{1, 1, 1, 1, 2, 2}; -const int nBetheParams = 6; -std::vector hfMothCodes = {511, 521, 531, 541, 5122}; // b-mesons + Lambda_b -static const std::vector betheBlochParNames{"p0", "p1", "p2", "p3", "p4", "resolution"}; -constexpr double kBetheBlochDefault[nParticles][nBetheParams]{ - {13.611469, 3.598765, -0.021138, 2.039562, 0.651040, 0.09}, // pion - {5.393020, 7.859534, 0.004048, 2.323197, 1.609307, 0.09}, // proton - {5.393020, 7.859534, 0.004048, 2.323197, 1.609307, 0.09}, // deuteron - {5.393020, 7.859534, 0.004048, 2.323197, 1.609307, 0.09}, // triton - {-126.557359, -0.858569, 1.111643, 1.210323, 2.656374, 0.09}, // helion - {-126.557359, -0.858569, 1.111643, 1.210323, 2.656374, 0.09}}; // alpha -const int nTrkSettings = 13; -static const std::vector trackPIDsettingsNames{"useBBparams", "minITSnCls", "minITSnClscos", "minTPCnCls", "maxTPCchi2", "minTPCchi2", "maxITSchi2", "maxTPCnSigma", "maxDcaXY", "maxDcaZ", "minITSclsSize", "minTPCnClsCrossedRows", "minReqClusterITSib"}; -constexpr double kTrackPIDSettings[nParticles][nTrkSettings]{ - {0, 0, 4, 60, 4.0, 0.5, 100, 2.5, 2., 2., 0., 70, 1}, - {1, 0, 4, 70, 4.0, 0.5, 100, 3.0, 2., 2., 0., 70, 1}, - {1, 0, 4, 70, 4.0, 0.5, 100, 3.0, 2., 2., 0., 70, 1}, - {1, 0, 4, 70, 4.0, 0.5, 100, 3.0, 2., 2., 0., 70, 1}, - {1, 0, 4, 75, 4.0, 0.5, 100, 5.0, 2., 2., 0., 70, 1}, - {1, 0, 4, 70, 4.0, 0.5, 100, 5.0, 2., 2., 0., 70, 1}}; - -const int nTrkSettings2 = 7; -static const std::vector trackPIDsettingsNames2{"useITSnsigma", "minITSnsigma", "maxITSnsigma", "fillsparsh", "useTPCnsigmaTOF", "maxTPCnsigmaTOF", "maxTOFnsigma"}; -constexpr double kTrackPIDSettings2[nParticles][nTrkSettings2]{ - {1, -5, 4, 0, 1, 2, 2}, - {1, -5, 4, 0, 1, 2, 2}, - {1, -5, 4, 0, 1, 2, 2}, - {1, -5, 4, 1, 1, 2, 2}, - {1, -5, 4, 1, 1, 2, 2}, - {1, -5, 4, 1, 1, 2, 2}}; - -static const int nfittingparticle = 4; -const int nfittingparameters = 4; -static const std::vector trackcorrectionNames{"a", "b", "c", "d"}; -constexpr double ktrackcorrection[nfittingparticle][nfittingparameters]{ - {0.464215, 0.195771, 0.0183111, 0.0}, // He3 - {0.464215, 0.195771, 0.0183111, 0.0}, // anti-He3 - {0.00765, 0.503791, -1.10517, 0.0}, // He4 - {0.00765, 0.503791, -1.10517, 0.0}}; // anti-He4 - -// DCA parameters - simple 2x3 array (DCAxy, DCAz) with p0, p1, p2 -static const int nDCAtypes = 2; -const int nDCAparameters = 3; -static const std::vector dcaNames{"p0", "p1", "p2"}; -static const std::vector dcaTypeNames{"DCAxy", "DCAz"}; -constexpr double kDCAcorrection[nDCAtypes][nDCAparameters]{ - {0.0118, -0.6889, 0.0017}, // DCAxy: p0*exp(p1*pt) + p2 - {0.1014, 1.7512, 0.0024}}; // DCAz: p0*exp(-p1*pt) + p2 - -struct PrimParticles { - TString name; - int pdgCode, charge; - double mass, resolution; - std::vector betheParams; - bool active; - PrimParticles(std::string name_, int pdgCode_, double mass_, int charge_, LabeledArray bethe) : name(name_), pdgCode(pdgCode_), charge(charge_), mass(mass_), active(false) - { - resolution = bethe.get(name, "resolution"); - betheParams.clear(); - constexpr unsigned int kNSpecies = 5; - for (unsigned int i = 0; i < kNSpecies; i++) - betheParams.push_back(bethe.get(name, i)); - } -}; // struct PrimParticles - -//---------------------------------------------------------------------------------------------------------------- -std::vector> hmass; -std::vector> hmassnsigma; -} // namespace -//---------------------------------------------------------------------------------------------------------------- -struct NucleitpcPbPb { - - Preslice tracksPerCollision = aod::track::collisionId; - - HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject}; - HistogramRegistry histomc{"histomc", {}, OutputObjHandlingPolicy::AnalysisObject, true, true}; - - // Event Selection Configurables - Configurable removeNoSameBunchPileup{"removeNoSameBunchPileup", false, "Remove no same bunch pileup"}; - Configurable requireIsGoodZvtxFT0vsPV{"requireIsGoodZvtxFT0vsPV", false, "Require is good Zvtx FT0 vs PV"}; - Configurable requireIsVertexITSTPC{"requireIsVertexITSTPC", false, "Require is vertex ITS TPC"}; - Configurable removeNoTimeFrameBorder{"removeNoTimeFrameBorder", false, "Remove no time frame border"}; - Configurable cfgsel8Require{"cfgsel8Require", true, "sel8 cut require"}; - Configurable cfgZvertex{"cfgZvertex", 10, "Min Z Vertex"}; - - // Track Selection Configurables - Configurable cfgUsePVcontributors{"cfgUsePVcontributors", true, "use tracks that are PV contibutors"}; - Configurable cfgPassedITSRefit{"cfgPassedITSRefit", true, "Require ITS refit"}; - Configurable cfgPassedTPCRefit{"cfgPassedTPCRefit", true, "Require TPC refit"}; - Configurable cfgetaRequire{"cfgetaRequire", true, "eta cut require"}; - Configurable cfgetaRequireMC{"cfgetaRequireMC", true, "eta cut require for generated particles"}; - Configurable cfgRapidityRequire{"cfgRapidityRequire", true, "Require Rapidity cut"}; - Configurable cfgRapidityRequireMC{"cfgRapidityRequireMC", true, "rapidity cut require for generated particles"}; - Configurable cfgCutEta{"cfgCutEta", 0.9f, "Eta range for tracks"}; - Configurable cfgCutRapiditymin{"cfgCutRapiditymin", -0.5f, "min Rapidity range"}; - Configurable cfgCutRapiditymax{"cfgCutRapiditymax", 0.5f, "max Rapidity range"}; - Configurable cfgtpcNClsFindable{"cfgtpcNClsFindable", 0.8f, "tpcNClsFindable over crossedRows"}; - Configurable centcut{"centcut", 80.0f, "centrality cut"}; - Configurable cfgZvertexRequireMC{"cfgZvertexRequireMC", true, "Pos Z cut in MC"}; - - // Track Cut Configurables - Configurable cfgTPCNClsCrossedRowsRequire{"cfgTPCNClsCrossedRowsRequire", true, "Require TPCNClsCrossedRows Cut"}; - Configurable cfgmaxTPCchi2Require{"cfgmaxTPCchi2Require", true, "Require maxTPCchi2 Cut"}; - Configurable cfgminTPCchi2Require{"cfgminTPCchi2Require", true, "Require minTPCchi2 Cut"}; - Configurable cfgminITSnClsRequire{"cfgminITSnClsRequire", false, "Require minITSnCls Cut"}; - Configurable cfgminITSnClscosRequire{"cfgminITSnClscosRequire", true, "Require minITSnCls / cosh(eta) Cut"}; - Configurable cfgminReqClusterITSibRequire{"cfgminReqClusterITSibRequire", true, " Require min number of clusters required in ITS inner barrel"}; - Configurable cfgmaxITSchi2Require{"cfgmaxITSchi2Require", true, "Require maxITSchi2 Cut"}; - Configurable cfgmaxTPCnSigmaRequire{"cfgmaxTPCnSigmaRequire", true, "Require maxTPCnSigma Cut"}; - Configurable cfgminGetMeanItsClsSizeRequire{"cfgminGetMeanItsClsSizeRequire", true, "Require minGetMeanItsClsSize Cut"}; - Configurable cfgmaxGetMeanItsClsSizeRequire{"cfgmaxGetMeanItsClsSizeRequire", true, "Require maxGetMeanItsClsSize Cut"}; - - // He4 Configurables - Configurable cfghe3massrejreq{"cfghe3massrejreq", true, "Require mass cut on He4 particles"}; - Configurable cfgminmassrejection{"cfgminmassrejection", 6.5, "Min side of He3 particle rejection"}; - Configurable cfgmaxmassrejection{"cfgmaxmassrejection", 9.138, "Max side of He3 particle rejection"}; - Configurable deuteronsigmarejection{"deuteronsigmarejection", 2.0f, "Deuteron TPC nsigma rejection for He4"}; - - // MC Configurables - Configurable cfgmccorrectionhe4Require{"cfgmccorrectionhe4Require", true, "MC correction for pp he4 particle"}; - - // DCA Configurables - booleans to enable/disable pT-dependent cuts - Configurable cfgUseDCAxyCorrection{"cfgUseDCAxyCorrection", true, "Use pT-dependent DCAxy cut"}; - Configurable cfgUseDCAzCorrection{"cfgUseDCAzCorrection", true, "Use pT-dependent DCAz cut"}; - - Configurable cfgDebug{"cfgDebug", 1, "debug level"}; - Configurable cfgRigidityCorrection{"cfgRigidityCorrection", false, "apply rigidity correction"}; - Configurable cfgRequirebetaplot{"cfgRequirebetaplot", true, "Require beta plot"}; - Configurable cfgmass2{"cfgmass2", true, "Fill mass square difference"}; - Configurable cfgFillmass{"cfgFillmass", false, "Fill mass histograms"}; - Configurable cfgFillmassnsigma{"cfgFillmassnsigma", true, "Fill mass vs nsigma histograms"}; - - Configurable> cfgBetheBlochParams{"cfgBetheBlochParams", {kBetheBlochDefault[0], nParticles, nBetheParams, particleNames, betheBlochParNames}, "TPC Bethe-Bloch parameterisation for light nuclei"}; - Configurable> cfgTrackPIDsettings{"cfgTrackPIDsettings", {kTrackPIDSettings[0], nParticles, nTrkSettings, particleNames, trackPIDsettingsNames}, "track selection and PID criteria"}; - Configurable> cfgTrackPIDsettings2{"cfgTrackPIDsettings2", {kTrackPIDSettings2[0], nParticles, nTrkSettings2, particleNames, trackPIDsettingsNames2}, "track selection and PID criteria"}; - Configurable> cfgktrackcorrection{"cfgktrackcorrection", {ktrackcorrection[0], nfittingparticle, nfittingparameters, correctedparticleNames, trackcorrectionNames}, "fitting parameters"}; - - // DCA correction parameters - simple 2x3 array - Configurable> cfgDCAcorrection{"cfgDCAcorrection", {kDCAcorrection[0], nDCAtypes, nDCAparameters, dcaTypeNames, dcaNames}, "DCA parameters: DCAxy: p0*exp(p1*pt)+p2, DCAz: p0*exp(-p1*pt)+p2"}; - - o2::track::TrackParametrizationWithError mTrackParCov; - // Binning configuration - ConfigurableAxis axisMagField{"axisMagField", {10, -10., 10.}, "magnetic field"}; - ConfigurableAxis axisNev{"axisNev", {10, 0., 10.}, "Number of events"}; - ConfigurableAxis axisRigidity{"axisRigidity", {4000, -10., 10.}, "#it{p}^{TPC}/#it{z}"}; - ConfigurableAxis axisdEdx{"axisdEdx", {4000, 0, 4000}, "d#it{E}/d#it{x}"}; - ConfigurableAxis axisCent{"axisCent", {100, 0, 100}, "centrality"}; - ConfigurableAxis axisOccupancy{"axisOccupancy", {5000, 0, 50000}, "axis for Occupancy of event"}; - ConfigurableAxis axisVtxZ{"axisVtxZ", {120, -20, 20}, "z"}; - ConfigurableAxis ptAxis{"ptAxis", {200, 0, 10}, "#it{p}_{T} (GeV/#it{c})"}; - ConfigurableAxis axiseta{"axiseta", {100, -1, 1}, "eta"}; - ConfigurableAxis axisrapidity{"axisrapidity", {100, -2, 2}, "rapidity"}; - ConfigurableAxis axismass{"axismass", {1200, 0, 12}, "mass"}; - ConfigurableAxis axismassnsigma{"axismassnsigma", {1200, 0, 12}, "nsigma mass"}; - ConfigurableAxis nsigmaAxis{"nsigmaAxis", {160, -10, 10}, "n#sigma_{#pi^{+}}"}; - ConfigurableAxis axisDCA{"axisDCA", {400, -10., 10.}, "DCA axis"}; - ConfigurableAxis particleAntiAxis{"particleAntiAxis", {2, -0.5, 1.5}, "Particle/Anti-particle"}; // 0 = particle, 1 = anti-particle - ConfigurableAxis decayTypeAxis{"decayTypeAxis", {3, -0.5, 2.5}, "Decay type"}; // 0 = primary, 1 = from decay, 2 = material - - // CCDB - Service ccdb; - Configurable bField{"bField", -999, "bz field, -999 is automatic"}; - Configurable ccdbUrl{"ccdbUrl", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; - Configurable grpPath{"grpPath", "GLO/GRP/GRP", "Path of the grp file"}; - Configurable grpmagPath{"grpmagPath", "GLO/Config/GRPMagField", "CCDB path of the GRPMagField object"}; - Configurable lutPath{"lutPath", "GLO/Param/MatLUT", "Path of the Lut parametrization"}; - Configurable geoPath{"geoPath", "GLO/Config/GeometryAligned", "Path of the geometry file"}; - Configurable pidPath{"pidPath", "", "Path to the PID response object"}; - - std::vector primaryParticles; - std::vector primVtx, cents; - bool collHasCandidate = false; - bool collPassedEvSel = false; - int mRunNumber = 0; - int occupancy = 0; - float dBz = 0.0f; - TRandom3 rand; - float proton = 1; - float de = 2; - float triton = 3; - float he3 = 4; - float he4 = 5; - - //---------------------------------------------------------------------------------------------------------------- - void init(InitContext const&) - { - mRunNumber = 0; - dBz = 0; - rand.SetSeed(0); - ccdb->setURL(ccdbUrl); - ccdb->setCaching(true); - ccdb->setLocalObjectValidityChecking(); - ccdb->setFatalWhenNull(false); - for (int i = 0; i < nParticles; i++) { // create primaryParticles - primaryParticles.push_back(PrimParticles(particleNames.at(i), particlePdgCodes.at(i), particleMasses.at(i), particleCharge.at(i), cfgBetheBlochParams)); - } - // create histograms - if (doprocessData) { - histos.add("histNev", "histNev", kTH1F, {axisNev}); - histos.add("histVtxZ", "histVtxZ", kTH1F, {axisVtxZ}); - histos.add("histCentFT0C", "histCentFT0C", kTH1F, {axisCent}); - histos.add("histCentFT0M", "histCentFT0M", kTH1F, {axisCent}); - histos.add("histCentFTOC_cut", "histCentFTOC_cut", kTH1F, {axisCent}); - histos.add("hSpectra", " ", HistType::kTHnSparseF, {ptAxis, nsigmaAxis, {5, -2.5, 2.5}, axisCent}); - histos.add("hSpectratof", " ", HistType::kTHnSparseF, {ptAxis, nsigmaAxis, {5, -2.5, 2.5}, axisCent}); - histos.add("DCAxy_vs_pT_data", "DCA_{xy} vs p_{T} for He3 (Data);p_{T} (GeV/c);DCA_{xy} (cm)", - {HistType::kTH3F, {ptAxis, axisDCA, axisCent}}); - histos.add("DCAxy_vs_pT_anti_data", "DCA_{xy} vs p_{T} for anti-He3 (Data);p_{T} (GeV/c);DCA_{xy} (cm)", - {HistType::kTH3F, {ptAxis, axisDCA, axisCent}}); - - hmass.resize(2 * nParticles + 2); - hmassnsigma.resize(2 * nParticles + 2); - - for (int i = 0; i < nParticles; i++) { - TString histName = primaryParticles[i].name; - if (cfgFillmass) { - hmass[2 * i] = histos.add(Form("histmass_pt/histmass_%s", histName.Data()), ";p_T{TPC} (GeV/#it{c}); mass^{2}; centrality(%)", HistType::kTH3F, {ptAxis, axismass, axisCent}); - hmass[2 * i + 1] = histos.add(Form("histmass_ptanti/histmass_%s", histName.Data()), ";p_T{TPC} (GeV/#it{c}); mass^{2}; centrality(%)", HistType::kTH3F, {ptAxis, axismass, axisCent}); - } - } - for (int i = 0; i < nParticles; i++) { - TString histName = primaryParticles[i].name; - if (cfgFillmassnsigma) { - hmassnsigma[2 * i] = histos.add(Form("histmass_nsigma/histmass_%s", histName.Data()), ";p_T{TPC} (GeV/#it{c}); mass^{2}/z^{2}; Centrality(%)", HistType::kTH3F, {ptAxis, axismassnsigma, axisCent}); - hmassnsigma[2 * i + 1] = histos.add(Form("histmass_nsigmaanti/histmass_%s", histName.Data()), ";p_T{TPC} (GeV/#it{c}); mass^{2}/z^{2}; centrality(%)", HistType::kTH3F, {ptAxis, axismassnsigma, axisCent}); - } - } - - histos.add("histeta", "histeta", kTH1F, {axiseta}); - - histos.add("histEvents", "histEvents", kTH2F, {axisCent, axisOccupancy}); - } - - histos.add("dcaZ", "dcaZ", kTH2F, {ptAxis, axisDCA}); - histos.add("dcaXY", "dcaXY", kTH2F, {ptAxis, axisDCA}); - histos.add("Tofsignal", "Tofsignal", kTH2F, {axisRigidity, {4000, 0.2, 1.2, "#beta"}}); - histos.add("Tpcsignal", "Tpcsignal", kTH2F, {axisRigidity, axisdEdx}); - - if (doprocessMC) { - histomc.add("hSpectramc", " ", HistType::kTHnSparseF, {{5, -2.5, 2.5}, axisCent, ptAxis, ptAxis}); - - // Efficiency x Acceptance - histomc.add("hDenomEffAcc", "Denominator for Efficiency x Acceptance", - {HistType::kTHnSparseF, {ptAxis, axisCent, particleAntiAxis, decayTypeAxis}}); - histomc.add("hNumerEffAcc", "Numerator for Efficiency x Acceptance", - {HistType::kTHnSparseF, {ptAxis, axisCent, particleAntiAxis, decayTypeAxis}}); - - // The Signal loss correction - histomc.add("hSignalLossDenom", "Signal Loss Denominator", kTH2F, {axisCent, ptAxis}); - histomc.add("hSignalLossNumer", "Signal Loss Numerator", kTH2F, {axisCent, ptAxis}); - - histomc.add("haSignalLossDenom", "anti particle Signal Loss Denominator", kTH2F, {axisCent, ptAxis}); - histomc.add("haSignalLossNumer", "antiparticle Signal Loss Numerator", kTH2F, {axisCent, ptAxis}); - - // The event loss correction - histomc.add("hEventLossDenom", "Event loss denominator", kTH1F, {axisCent}); - histomc.add("hEventLossNumer", "Event loss numerator", kTH1F, {axisCent}); - - histomc.add("histVtxZgen", "histVtxZgen", kTH1F, {axisVtxZ}); - histomc.add("histVtxZReco", "histVtxZReco", kTH1F, {axisVtxZ}); - - histomc.add("histDeltaPtVsPtGenanti", " delta pt vs pt rec", HistType::kTH2F, {{1000, 0, 10}, {1000, -0.5, 0.5, "p_{T}(reco) - p_{T}(gen);p_{T}(reco)"}}); - histomc.add("histDeltaPtVsPtGen", " delta pt vs pt rec", HistType::kTH2F, {{1000, 0, 10}, {1000, -0.5, 0.5, "p_{T}(reco) - p_{T}(gen);p_{T}(reco)"}}); - histomc.add("histPIDtrack", " delta pt vs pt rec", HistType::kTH2F, {{1000, 0, 10, "p_{T}(reco)"}, {9, -0.5, 8.5, "p_{T}(reco) - p_{T}(gen)"}}); - histomc.add("histPIDtrackanti", " delta pt vs pt rec", HistType::kTH2F, {{1000, 0, 10, "p_{T}(reco)"}, {9, -0.5, 8.5, "p_{T}(reco) - p_{T}(gen)"}}); - - // Mass²/z² 3D plots for MC (pT, mass²/z², centrality) - separate for particles and anti-particles - histomc.add("hMassVsPtMC", "mass^{2}/z^{2} vs p_{T} (MC Particles);p_{T} (GeV/c);mass^{2}/z^{2};Centrality", - {HistType::kTH3F, {ptAxis, axismassnsigma, axisCent}}); - histomc.add("hMassVsPtAntiMC", "mass^{2}/z^{2} vs p_{T} (MC Anti-particles);p_{T} (GeV/c);mass^{2}/z^{2};Centrality", - {HistType::kTH3F, {ptAxis, axismassnsigma, axisCent}}); - } - - if (doprocessDCA) { - - // NEW: Add DCAxy vs pT histograms for MC - histomc.add("DCAxy_vs_pT_transport", "DCA_{xy} vs p_{T} (Transport);p_{T} (GeV/c);DCA_{xy} (cm)", - {HistType::kTH3F, {ptAxis, axisDCA, axisCent}}); - histomc.add("DCAxy_vs_pT_primary", "DCA_{xy} vs p_{T} (primary);p_{T} (GeV/c);DCA_{xy} (cm)", - {HistType::kTH3F, {ptAxis, axisDCA, axisCent}}); - histomc.add("DCAxy_vs_pT_weakdecay", "DCA_{xy} vs p_{T} (Weak Decay);p_{T} (GeV/c);DCA_{xy} (cm)", - {HistType::kTH3F, {ptAxis, axisDCA, axisCent}}); - histomc.add("DCAxy_vs_pT_total", "DCA_{xy} vs p_{T} (Total);p_{T} (GeV/c);DCA_{xy} (cm)", - {HistType::kTH3F, {ptAxis, axisDCA, axisCent}}); - histomc.add("DCAxy_vs_pT_anti_total", "DCA_{xy} vs p_{T} for anti (Total);p_{T} (GeV/c);DCA_{xy} (cm)", - {HistType::kTH3F, {ptAxis, axisDCA, axisCent}}); - } - } - //---------------------------------------------------------------------------------------------------------------- - void processData(CollisionsFull const& collisions, - TracksFull const& tracks, - aod::BCsWithTimestamps const&) - { - // Get DCA parameters from configurable - double dcaXY_p0 = cfgDCAcorrection->get(0U, "p0"); - double dcaXY_p1 = cfgDCAcorrection->get(0U, "p1"); - double dcaXY_p2 = cfgDCAcorrection->get(0U, "p2"); - double dcaZ_p0 = cfgDCAcorrection->get(1U, "p0"); - double dcaZ_p1 = cfgDCAcorrection->get(1U, "p1"); - double dcaZ_p2 = cfgDCAcorrection->get(1U, "p2"); - - for (const auto& collision : collisions) { - auto bc = collision.bc_as(); - initCCDB(bc); - initCollision(collision); - - if (!collPassedEvSel) - continue; - if (removeNoSameBunchPileup && !collision.selection_bit(aod::evsel::kNoSameBunchPileup)) - continue; - histos.fill(HIST("histNev"), 2.5); - - if (requireIsGoodZvtxFT0vsPV && !collision.selection_bit(aod::evsel::kIsGoodZvtxFT0vsPV)) - continue; - histos.fill(HIST("histNev"), 3.5); - if (requireIsVertexITSTPC && !collision.selection_bit(aod::evsel::kIsVertexITSTPC)) - continue; - - histos.fill(HIST("histNev"), 4.5); - if (removeNoTimeFrameBorder && !collision.selection_bit(aod::evsel::kNoTimeFrameBorder)) - continue; - histos.fill(HIST("histNev"), 5.5); - histos.fill(HIST("histEvents"), collision.centFT0C(), occupancy); - histos.fill(HIST("histVtxZ"), collision.posZ()); - histos.fill(HIST("histCentFT0C"), collision.centFT0C()); - histos.fill(HIST("histCentFT0M"), collision.centFT0M()); - if (collision.centFT0C() > centcut) - continue; - histos.fill(HIST("histNev"), 6.5); - histos.fill(HIST("histCentFTOC_cut"), collision.centFT0C()); - // new slicing - auto tracksInColl = tracks.sliceBy(tracksPerCollision, collision.globalIndex()); - // loop over sliced tracks - for (const auto& track : tracksInColl) { - if (!track.isPVContributor() && cfgUsePVcontributors) - continue; - if (!track.hasTPC()) - continue; - if (!track.passedITSRefit() && cfgPassedITSRefit) - continue; - if (!track.passedTPCRefit() && cfgPassedTPCRefit) - continue; - if (std::abs(track.eta()) > cfgCutEta && cfgetaRequire) - continue; - for (size_t i = 0; i < primaryParticles.size(); i++) { - if (cfgTrackPIDsettings2->get(i, "fillsparsh") != 1) - continue; - float ptMomn; - setTrackParCov(track, mTrackParCov); - mTrackParCov.setPID(track.pidForTracking()); - ptMomn = (i == he3 || i == he4) ? 2 * mTrackParCov.getPt() : mTrackParCov.getPt(); - float tpcNsigma = getTPCnSigma(track, primaryParticles.at(i)); - if ((std::abs(tpcNsigma) > cfgTrackPIDsettings->get(i, "maxTPCnSigma")) && cfgmaxTPCnSigmaRequire) - continue; - int sign = (track.sign() > 0) ? 1 : ((track.sign() < 0) ? -1 : 0); - float rapidity = getRapidity(track, i); - if ((rapidity < cfgCutRapiditymin || rapidity > cfgCutRapiditymax) && cfgRapidityRequire) - continue; - if (track.tpcNClsFound() < cfgTrackPIDsettings->get(i, "minTPCnCls")) - continue; - if (((track.tpcNClsCrossedRows() < cfgTrackPIDsettings->get(i, "minTPCnClsCrossedRows")) || - track.tpcNClsCrossedRows() < cfgtpcNClsFindable * track.tpcNClsFindable()) && - cfgTPCNClsCrossedRowsRequire) - continue; - if (track.tpcChi2NCl() > cfgTrackPIDsettings->get(i, "maxTPCchi2") && cfgmaxTPCchi2Require) - continue; - if (track.tpcChi2NCl() < cfgTrackPIDsettings->get(i, "minTPCchi2") && cfgminTPCchi2Require) - continue; - if (track.itsNCls() < cfgTrackPIDsettings->get(i, "minITSnCls") && cfgminITSnClsRequire) - continue; - double cosheta = std::cosh(track.eta()); - if ((track.itsNCls() / cosheta) < cfgTrackPIDsettings->get(i, "minITSnClscos") && cfgminITSnClscosRequire) - continue; - if ((track.itsNClsInnerBarrel() < cfgTrackPIDsettings->get(i, "minReqClusterITSib")) && cfgminReqClusterITSibRequire) - continue; - if (track.itsChi2NCl() > cfgTrackPIDsettings->get(i, "maxITSchi2") && cfgmaxITSchi2Require) - continue; - if (getMeanItsClsSize(track) < cfgTrackPIDsettings->get(i, "minITSclsSize") && cfgminGetMeanItsClsSizeRequire) - continue; - - // DCAxy cut - bool insideDCAxy = false; - if (cfgUseDCAxyCorrection) { - // Use pT-dependent DCAxy cut - double sigmaFactor = cfgTrackPIDsettings->get(i, "maxDcaXY"); - double sigma_base = dcaXY_p0 * std::exp(dcaXY_p1 * ptMomn) + dcaXY_p2; - double sigma_new = sigmaFactor * sigma_base; - insideDCAxy = (std::abs(track.dcaXY()) <= sigma_new); - } else { - // Use simple DCAxy cut (no pT dependence) - insideDCAxy = (std::abs(track.dcaXY()) <= cfgTrackPIDsettings->get(i, "maxDcaXY")); - } - - // DCAz cut - bool insideDCAz = false; - if (cfgUseDCAzCorrection) { - // Use pT-dependent DCAz cut - double sigmaFactorZ = cfgTrackPIDsettings->get(i, "maxDcaZ"); - double sigma_base_z = dcaZ_p0 * std::exp(-dcaZ_p1 * ptMomn) + dcaZ_p2; - double sigma_new_z = sigmaFactorZ * sigma_base_z; - insideDCAz = (std::abs(track.dcaZ()) <= sigma_new_z); - } else { - // Use simple DCAz cut (no pT dependence) - insideDCAz = (std::abs(track.dcaZ()) <= cfgTrackPIDsettings->get(i, "maxDcaZ")); - } - - if ((!insideDCAxy || !insideDCAz)) { - continue; - } - - float itsSigma = getITSnSigma(track, primaryParticles.at(i)); - if (itsSigma < cfgTrackPIDsettings2->get(i, "minITSnsigma") && cfgTrackPIDsettings2->get(i, "useITSnsigma") < 1) - continue; - if (itsSigma > cfgTrackPIDsettings2->get(i, "maxITSnsigma") && cfgTrackPIDsettings2->get(i, "useITSnsigma") < 1) - continue; - - histos.fill(HIST("Tpcsignal"), getRigidity(track) * track.sign(), track.tpcSignal()); - histos.fill(HIST("dcaXY"), ptMomn, track.dcaXY()); - histos.fill(HIST("dcaZ"), ptMomn, track.dcaZ()); - float tofnsigma = -999; - if (i == proton) { - tofnsigma = track.tofNSigmaPr(); - } - if (track.sign() > 0) { - histos.fill(HIST("DCAxy_vs_pT_data"), ptMomn, track.dcaXY(), collision.centFT0C()); - } else if (track.sign() < 0) { - histos.fill(HIST("DCAxy_vs_pT_anti_data"), ptMomn, track.dcaXY(), collision.centFT0C()); - } - // Get deuteron TPC nsigma for He4 selection - float tpcNsigmaDe = -999; - tpcNsigmaDe = track.tpcNSigmaDe(); - - if (i != he4) { - histos.fill(HIST("hSpectra"), ptMomn, tpcNsigma, sign, collision.centFT0C()); - if (track.hasTOF() && i == proton) { - histos.fill(HIST("hSpectratof"), ptMomn, tofnsigma, sign, collision.centFT0C()); - } - } else { - - if (!track.hasTOF()) { - if (std::abs(tpcNsigmaDe) > deuteronsigmarejection) { - histos.fill(HIST("hSpectra"), ptMomn, tpcNsigma, sign, collision.centFT0C()); - } - } else { - // Has TOF - apply mass cut - float beta = o2::pid::tof::Beta::GetBeta(track); - const float eps = 1e-6f; - if (beta < eps || beta > 1.0f - eps) - continue; - - float charge = 2.f; // he4 has charge 2 - float p = getRigidity(track); - float massTOF = p * charge * std::sqrt(1.f / (beta * beta) - 1.f); - - // Apply mass cut for he4 (mass^2 around 3.73^2 = 13.9) - if (cfghe3massrejreq && (massTOF * massTOF > cfgminmassrejection && massTOF * massTOF < cfgmaxmassrejection)) { - continue; // Skip if mass cut fails - } - if (std::abs(tpcNsigmaDe) > deuteronsigmarejection) { - histos.fill(HIST("hSpectra"), ptMomn, tpcNsigma, sign, collision.centFT0C()); - } - } - } - - if ((std::abs(tpcNsigma) > cfgTrackPIDsettings2->get(i, "maxTPCnsigmaTOF"))) { - if (i == he4) { - if (std::abs(tpcNsigmaDe) > deuteronsigmarejection) { - fillhmassnsigma(track, i, collision.centFT0C()); - } - } else { - fillhmassnsigma(track, i, collision.centFT0C()); - } - } - - if ((std::abs(tpcNsigma) > cfgTrackPIDsettings2->get(i, "maxTPCnsigmaTOF")) && cfgTrackPIDsettings2->get(i, "useTPCnsigmaTOF") < 1) - continue; - - if (i == he4) { - if (std::abs(tpcNsigmaDe) > deuteronsigmarejection) { - fillhmass(track, i, collision.centFT0C()); - } - } else { - if (i == proton && (std::abs(tofnsigma) < cfgTrackPIDsettings2->get(i, "maxTOFnsigma"))) { - fillhmassnsigma(track, i, collision.centFT0C()); - } else { - fillhmass(track, i, collision.centFT0C()); - } - } - - if (cfgRequirebetaplot) { - histos.fill(HIST("Tofsignal"), getRigidity(track) * track.sign(), o2::pid::tof::Beta::GetBeta(track)); - } - } // loop primaryParticles - - histos.fill(HIST("histeta"), track.eta()); - } // loop sliced tracks - } // collision loop - } - PROCESS_SWITCH(NucleitpcPbPb, processData, "data analysis", false); - - //---------------------------------------------------------------------------------------------------------------- - // MC particles - Efficiency x Acceptance and Signal Loss calculations - //-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=- - struct McCollInfo { - bool passedEvSel = false; - float centrality = -1.0f; - bool passedEvSelVtZ = false; - }; - std::vector mcCollInfos; - - void processMC(CollisionsFullMC const& collisions, - aod::McCollisions const& mcCollisions, - soa::Join const& tracks, - aod::McParticles const& particlesMC, - aod::BCsWithTimestamps const&) - - { - // Get DCA parameters from configurable - double dcaXY_p0 = cfgDCAcorrection->get(0U, "p0"); - double dcaXY_p1 = cfgDCAcorrection->get(0U, "p1"); - double dcaXY_p2 = cfgDCAcorrection->get(0U, "p2"); - double dcaZ_p0 = cfgDCAcorrection->get(1U, "p0"); - double dcaZ_p1 = cfgDCAcorrection->get(1U, "p1"); - double dcaZ_p2 = cfgDCAcorrection->get(1U, "p2"); - - mcCollInfos.clear(); - mcCollInfos.resize(mcCollisions.size()); - - // First pass: Store centrality and apply event selection - for (auto const& collision : collisions) { - int mcCollIdx = collision.mcCollisionId(); - if (mcCollIdx < 0 || mcCollIdx >= static_cast(mcCollisions.size())) { - continue; - } - - // STORE CENTRALITY WITHOUt CUTS - mcCollInfos[mcCollIdx].centrality = collision.centFT0C(); - - if (!collision.sel8() && cfgsel8Require) - continue; - if (collision.centFT0C() > centcut) - continue; - - if (removeNoSameBunchPileup && !collision.selection_bit(aod::evsel::kNoSameBunchPileup)) - continue; - if (requireIsGoodZvtxFT0vsPV && !collision.selection_bit(aod::evsel::kIsGoodZvtxFT0vsPV)) - continue; - if (requireIsVertexITSTPC && !collision.selection_bit(aod::evsel::kIsVertexITSTPC)) - continue; - if (removeNoTimeFrameBorder && !collision.selection_bit(aod::evsel::kNoTimeFrameBorder)) - continue; - - // Mark this MC collision as passing event selection - mcCollInfos[mcCollIdx].passedEvSel = true; - - // Apply event selection cuts - if (std::abs(collision.posZ()) > cfgZvertex && cfgZvertexRequireMC) - continue; - - mcCollInfos[mcCollIdx].passedEvSelVtZ = true; - } - - // FILL EVENT LOSS AND SIGNAL LOSS: Combined loop per MC collision - for (size_t i = 0; i < mcCollInfos.size(); i++) { - if (mcCollInfos[i].centrality >= 0) { // Only if we found a matching collision - // Event loss denominator (always filled regardless of fillsparsh) - histomc.fill(HIST("hEventLossDenom"), mcCollInfos[i].centrality); - - // Event loss numerator (if passed selection) - if (mcCollInfos[i].passedEvSel) { - histomc.fill(HIST("hEventLossNumer"), mcCollInfos[i].centrality); - } - - // Fill signal loss for all primary particles in this MC collision - for (auto const& mcParticle : particlesMC) { - if (mcParticle.mcCollisionId() != static_cast(i)) { - continue; - } - if (!mcParticle.isPhysicalPrimary()) { - continue; - } - - // Check if this particle type has fillsparsh enabled - int particleType = -1; - for (size_t j = 0; j < primaryParticles.size(); j++) { - if (std::abs(mcParticle.pdgCode()) == std::abs(particlePdgCodes.at(j))) { - particleType = j; - break; - } - } - - // Skip if fillsparsh is not enabled for this particle - if (particleType < 0 || cfgTrackPIDsettings2->get(particleType, "fillsparsh") != 1) - continue; - - // Signal loss denominator - if (mcParticle.pdgCode() == particlePdgCodes.at(particleType)) { // particle - histomc.fill(HIST("hSignalLossDenom"), mcCollInfos[i].centrality, mcParticle.pt()); - } else if (mcParticle.pdgCode() == -particlePdgCodes.at(particleType)) { // anti-particle - histomc.fill(HIST("haSignalLossDenom"), mcCollInfos[i].centrality, mcParticle.pt()); - } - // Signal loss numerator (if event passed selection) - if (mcCollInfos[i].passedEvSel) { - if (mcParticle.pdgCode() == particlePdgCodes.at(particleType)) { // particle - histomc.fill(HIST("hSignalLossNumer"), mcCollInfos[i].centrality, mcParticle.pt()); - } else if (mcParticle.pdgCode() == -particlePdgCodes.at(particleType)) { // anti-particle - histomc.fill(HIST("haSignalLossNumer"), mcCollInfos[i].centrality, mcParticle.pt()); - } - } - } - } - } - - // Process MC collisions for efficiency and reconstructed collisions - for (auto const& mcCollision : mcCollisions) { - size_t idx = mcCollision.globalIndex(); - if (idx >= mcCollInfos.size()) - continue; - - float centrality = mcCollInfos[idx].centrality; - bool passedEvSelVtZ = mcCollInfos[idx].passedEvSelVtZ; - - // Process generated particles for efficiency denominators - histomc.fill(HIST("histVtxZgen"), mcCollision.posZ()); - - for (auto const& mcParticle : particlesMC) { - if (mcParticle.mcCollisionId() != mcCollision.globalIndex()) - continue; - - int pdgCode = mcParticle.pdgCode(); - - // Check which particle type this is - int particleType = -1; - for (size_t j = 0; j < primaryParticles.size(); j++) { - if (std::abs(pdgCode) == std::abs(particlePdgCodes.at(j))) { - particleType = j; - break; - } - } - - // Only process if this particle type has fillsparsh enabled - if (particleType < 0 || cfgTrackPIDsettings2->get(particleType, "fillsparsh") != 1) - continue; - - if (std::abs(mcParticle.eta()) > cfgCutEta && cfgetaRequireMC) - continue; - float rapidity = mcParticle.y(); - if ((rapidity < cfgCutRapiditymin || rapidity > cfgCutRapiditymax) && cfgRapidityRequireMC) - continue; - - int decayType = 0; - int particleAnti = (pdgCode > 0) ? 0 : 1; - - if (mcParticle.isPhysicalPrimary()) { - decayType = 0; - if (mcParticle.has_mothers()) { - for (const auto& motherparticle : mcParticle.mothers_as()) { - if (std::find(hfMothCodes.begin(), hfMothCodes.end(), - std::abs(motherparticle.pdgCode())) != hfMothCodes.end()) { - decayType = 1; - break; - } - } - } - } else if (mcParticle.has_mothers()) { - decayType = 1; - } else { - decayType = 2; - continue; - } - - // Efficiency x Acceptance histograms - Denominator - if (passedEvSelVtZ) { - histomc.fill(HIST("hDenomEffAcc"), mcParticle.pt(), centrality, particleAnti, decayType); - } - } - - // Process reconstructed collisions for this MC collision - if (passedEvSelVtZ) { - // Find the corresponding reconstructed collision - for (auto const& collision : collisions) { - if (collision.mcCollisionId() != static_cast(idx)) - continue; - - auto bc = collision.bc_as(); - initCCDB(bc); - - histomc.fill(HIST("histVtxZReco"), collision.posZ()); - - auto tracksInColl = tracks.sliceBy(tracksPerCollision, collision.globalIndex()); - - for (auto const& track : tracksInColl) { - if (!track.has_mcParticle()) - continue; // skip un-matched reco tracks - - auto const& matchedMCParticle = track.mcParticle_as(); - - // Only process particles from this MC collision - if (matchedMCParticle.mcCollisionId() != mcCollision.globalIndex()) - continue; - - int pdg = matchedMCParticle.pdgCode(); - - int decayType = 0; - - if (matchedMCParticle.isPhysicalPrimary()) { - decayType = 0; - if (matchedMCParticle.has_mothers()) { - for (const auto& motherparticle : matchedMCParticle.mothers_as()) { - if (std::find(hfMothCodes.begin(), hfMothCodes.end(), - std::abs(motherparticle.pdgCode())) != hfMothCodes.end()) { - decayType = 1; - break; - } - } - } - } else if (matchedMCParticle.has_mothers()) { - decayType = 1; - } else { - decayType = 2; - } - - if (!track.isPVContributor() && cfgUsePVcontributors) - continue; - - if (!track.hasTPC()) - continue; - if (!track.passedITSRefit() && cfgPassedITSRefit) - continue; - if (!track.passedTPCRefit() && cfgPassedTPCRefit) - continue; - if (std::abs(track.eta()) > cfgCutEta && cfgetaRequire) - continue; - - for (size_t i = 0; i < primaryParticles.size(); i++) { - if (std::abs(pdg) != std::abs(particlePdgCodes.at(i))) - continue; - - // CRITICAL: Only process if fillsparsh is enabled for this particle - if (cfgTrackPIDsettings2->get(i, "fillsparsh") != 1) - continue; - - float ptReco; - setTrackParCov(track, mTrackParCov); - mTrackParCov.setPID(track.pidForTracking()); - - ptReco = (std::abs(pdg) == particlePdgCodes.at(4) || std::abs(pdg) == particlePdgCodes.at(5)) ? 2 * mTrackParCov.getPt() : mTrackParCov.getPt(); - - int particleAnti = (pdg > 0) ? 0 : 1; - - double a = 0, b = 0, c = 0; - - int param = -1; - if (i == he3) { - param = (pdg > 0) ? 0 : 1; - } else if (i == he4) { - param = (pdg > 0) ? 2 : 3; - } - - if (param >= 0) { - a = cfgktrackcorrection->get(param, "a"); - b = cfgktrackcorrection->get(param, "b"); - c = cfgktrackcorrection->get(param, "c"); - } - - if (std::abs(pdg) == particlePdgCodes.at(5) && cfgmccorrectionhe4Require) { - ptReco = ptReco + a + b * std::exp(c * ptReco); - } - - if (std::abs(pdg) == particlePdgCodes.at(4) && cfgmccorrectionhe4Require) { - int pidGuess = track.pidForTracking(); - int antitriton = 6; - if (pidGuess == antitriton) { - ptReco = ptReco - a + b * ptReco - c * ptReco * ptReco; - } - } - - float rapidity = getRapidity(track, i); - if ((rapidity < cfgCutRapiditymin || rapidity > cfgCutRapiditymax) && cfgRapidityRequire) - continue; - - if (track.tpcNClsFound() < cfgTrackPIDsettings->get(i, "minTPCnCls")) - continue; - if (((track.tpcNClsCrossedRows() < cfgTrackPIDsettings->get(i, "minTPCnClsCrossedRows")) || track.tpcNClsCrossedRows() < cfgtpcNClsFindable * track.tpcNClsFindable()) && cfgTPCNClsCrossedRowsRequire) - continue; - if (track.tpcChi2NCl() > cfgTrackPIDsettings->get(i, "maxTPCchi2") && cfgmaxTPCchi2Require) - continue; - if (track.tpcChi2NCl() < cfgTrackPIDsettings->get(i, "minTPCchi2") && cfgminTPCchi2Require) - continue; - if (track.itsNCls() < cfgTrackPIDsettings->get(i, "minITSnCls") && cfgminITSnClsRequire) - continue; - double cosheta = std::cosh(track.eta()); - if ((track.itsNCls() / cosheta) < cfgTrackPIDsettings->get(i, "minITSnClscos") && cfgminITSnClscosRequire) - continue; - if ((track.itsNClsInnerBarrel() < cfgTrackPIDsettings->get(i, "minReqClusterITSib")) && cfgminReqClusterITSibRequire) - continue; - if (track.itsChi2NCl() > cfgTrackPIDsettings->get(i, "maxITSchi2") && cfgmaxITSchi2Require) - continue; - if (getMeanItsClsSize(track) < cfgTrackPIDsettings->get(i, "minITSclsSize") && cfgminGetMeanItsClsSizeRequire) - continue; - - // DCAxy cut - bool insideDCAxy = false; - if (cfgUseDCAxyCorrection) { - // Use pT-dependent DCAxy cut - double sigmaFactor = cfgTrackPIDsettings->get(i, "maxDcaXY"); - double sigma_base = dcaXY_p0 * std::exp(dcaXY_p1 * ptReco) + dcaXY_p2; - double sigma_new = sigmaFactor * sigma_base; - insideDCAxy = (std::abs(track.dcaXY()) <= sigma_new); - } else { - // Use simple DCAxy cut (no pT dependence) - insideDCAxy = (std::abs(track.dcaXY()) <= cfgTrackPIDsettings->get(i, "maxDcaXY")); - } - - // DCAz cut - bool insideDCAz = false; - if (cfgUseDCAzCorrection) { - // Use pT-dependent DCAz cut - double sigmaFactorZ = cfgTrackPIDsettings->get(i, "maxDcaZ"); - double sigma_base_z = dcaZ_p0 * std::exp(-dcaZ_p1 * ptReco) + dcaZ_p2; - double sigma_new_z = sigmaFactorZ * sigma_base_z; - insideDCAz = (std::abs(track.dcaZ()) <= sigma_new_z); - } else { - // Use simple DCAz cut (no pT dependence) - insideDCAz = (std::abs(track.dcaZ()) <= cfgTrackPIDsettings->get(i, "maxDcaZ")); - } - - if ((!insideDCAxy || !insideDCAz)) { - continue; - } - - float tpcNsigma = getTPCnSigma(track, primaryParticles.at(i)); - if ((std::abs(tpcNsigma) > cfgTrackPIDsettings->get(i, "maxTPCnSigma")) && cfgmaxTPCnSigmaRequire) - continue; - - // Efficiency x Acceptance - Numerator - histomc.fill(HIST("hNumerEffAcc"), ptReco, collision.centFT0C(), particleAnti, decayType); - - // Sparse histogram - float ptTOF = -1.0; - if (track.hasTOF()) { - ptTOF = ptReco; - } - histomc.fill(HIST("hSpectramc"), particleAnti, collision.centFT0C(), ptReco, ptTOF); - - // Basic track histograms - if (decayType == 0) { - histos.fill(HIST("dcaXY"), ptReco, track.dcaXY()); - histos.fill(HIST("dcaZ"), ptReco, track.dcaZ()); - histos.fill(HIST("Tpcsignal"), getRigidity(track) * track.sign(), track.tpcSignal()); - } - // Delta Pt histograms - float ptGen = matchedMCParticle.pt(); - float deltaPt = ptReco - ptGen; - - if (pdg == -particlePdgCodes.at(i) && decayType == 0) { // Anti-particle - histomc.fill(HIST("histDeltaPtVsPtGenanti"), ptReco, deltaPt); - histomc.fill(HIST("histPIDtrackanti"), ptReco, track.pidForTracking()); - } - - if (pdg == particlePdgCodes.at(i) && decayType == 0) { // Particle - histomc.fill(HIST("histDeltaPtVsPtGen"), ptReco, deltaPt); - histomc.fill(HIST("histPIDtrack"), ptReco, track.pidForTracking()); - } - - // Fill mass²/z² for MC - separate for particles and anti-particles - if (track.hasTOF()) { - float beta = o2::pid::tof::Beta::GetBeta(track); - const float eps = 1e-6f; - if (beta >= eps && beta <= 1.0f - eps) { - float charge = (i == he3 || i == he4) ? 2.f : 1.f; - float p = getRigidity(track); - float massTOF = p * charge * std::sqrt(1.f / (beta * beta) - 1.f); - float massSquareOverChargeSquare = (massTOF * massTOF) / (charge * charge); - - // Apply He4 rejection if needed (to match data selection) - bool skipHe4 = false; - if (std::abs(pdg) == particlePdgCodes.at(5)) { // He4 - if (cfghe3massrejreq && (massTOF * massTOF > cfgminmassrejection && massTOF * massTOF < cfgmaxmassrejection)) { - skipHe4 = true; - } - } - - if (!skipHe4 && decayType == 0) { - if (pdg > 0) { - histomc.fill(HIST("hMassVsPtMC"), ptReco, massSquareOverChargeSquare, collision.centFT0C()); - } else { - histomc.fill(HIST("hMassVsPtAntiMC"), ptReco, massSquareOverChargeSquare, collision.centFT0C()); - } - } - } - } - } - } - break; // Found the matching collision, break out of collision loop - } - } - } - } - PROCESS_SWITCH(NucleitpcPbPb, processMC, "MC reco+gen analysis with efficiency corrections", false); - //=-=-=-==-=-=-==-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= - void processDCA(CollisionsFullMC const& collisions, - aod::McCollisions const& mcCollisions, - aod::McParticles const& particlesMC, - soa::Join const& tracks, - aod::BCsWithTimestamps const&) - - { - (void)particlesMC; - mcCollInfos.clear(); - mcCollInfos.resize(mcCollisions.size()); - - // First pass: Store centrality and apply event selection - for (auto const& collision : collisions) { - int mcCollIdx = collision.mcCollisionId(); - if (mcCollIdx < 0 || mcCollIdx >= static_cast(mcCollisions.size())) { - continue; - } - - // STORE CENTRALITY WITHOUT CUTS - mcCollInfos[mcCollIdx].centrality = collision.centFT0C(); - - if (!collision.sel8() && cfgsel8Require) - continue; - if (collision.centFT0C() > centcut) - continue; - - if (removeNoSameBunchPileup && !collision.selection_bit(aod::evsel::kNoSameBunchPileup)) - continue; - if (requireIsGoodZvtxFT0vsPV && !collision.selection_bit(aod::evsel::kIsGoodZvtxFT0vsPV)) - continue; - if (requireIsVertexITSTPC && !collision.selection_bit(aod::evsel::kIsVertexITSTPC)) - continue; - if (removeNoTimeFrameBorder && !collision.selection_bit(aod::evsel::kNoTimeFrameBorder)) - continue; - - // Mark this MC collision as passing event selection - mcCollInfos[mcCollIdx].passedEvSel = true; - - // Apply event selection cuts - if (std::abs(collision.posZ()) > cfgZvertex && cfgZvertexRequireMC) - continue; - - mcCollInfos[mcCollIdx].passedEvSelVtZ = true; - } - - // Process MC collisions for efficiency and reconstructed collisions - for (auto const& mcCollision : mcCollisions) { - size_t idx = mcCollision.globalIndex(); - if (idx >= mcCollInfos.size()) - continue; - - bool passedEvSelVtZ = mcCollInfos[idx].passedEvSelVtZ; - - // Process reconstructed collisions for this MC collision - if (passedEvSelVtZ) { - // Find the corresponding reconstructed collision - for (auto const& collision : collisions) { - if (collision.mcCollisionId() != static_cast(idx)) - continue; - - auto bc = collision.bc_as(); - initCCDB(bc); - auto tracksInColl = tracks.sliceBy(tracksPerCollision, collision.globalIndex()); - - for (auto const& track : tracksInColl) { - if (!track.has_mcParticle()) - continue; // skip un-matched reco tracks - - auto const& matchedMCParticle = track.mcParticle_as(); - - // Only process particles from this MC collision - if (matchedMCParticle.mcCollisionId() != mcCollision.globalIndex()) - continue; - - int pdg = matchedMCParticle.pdgCode(); - - if (!track.isPVContributor() && cfgUsePVcontributors) - continue; - - if (!track.hasTPC()) - continue; - if (!track.passedITSRefit() && cfgPassedITSRefit) - continue; - if (!track.passedTPCRefit() && cfgPassedTPCRefit) - continue; - if (std::abs(track.eta()) > cfgCutEta && cfgetaRequire) - continue; - - for (size_t i = 0; i < primaryParticles.size(); i++) { - if (std::abs(pdg) != std::abs(particlePdgCodes.at(i))) - continue; - - if (cfgTrackPIDsettings2->get(i, "fillsparsh") != 1) - continue; - - float ptReco; - setTrackParCov(track, mTrackParCov); - mTrackParCov.setPID(track.pidForTracking()); - - ptReco = (std::abs(pdg) == particlePdgCodes.at(4) || std::abs(pdg) == particlePdgCodes.at(5)) ? 2 * mTrackParCov.getPt() : mTrackParCov.getPt(); - - double a = 0, b = 0, c = 0; - - int param = -1; - if (i == he3) { - param = (pdg > 0) ? 0 : 1; - } else if (i == he4) { - param = (pdg > 0) ? 2 : 3; - } - - if (param >= 0) { - a = cfgktrackcorrection->get(param, "a"); - b = cfgktrackcorrection->get(param, "b"); - c = cfgktrackcorrection->get(param, "c"); - } - - if (std::abs(pdg) == particlePdgCodes.at(5) && cfgmccorrectionhe4Require) { - ptReco = ptReco + a + b * std::exp(c * ptReco); - } - - if (std::abs(pdg) == particlePdgCodes.at(4) && cfgmccorrectionhe4Require) { - int pidGuess = track.pidForTracking(); - int antitriton = 6; - if (pidGuess == antitriton) { - ptReco = ptReco - a + b * ptReco - c * ptReco * ptReco; - } - } - - float rapidity = getRapidity(track, i); - if ((rapidity < cfgCutRapiditymin || rapidity > cfgCutRapiditymax) && cfgRapidityRequire) - continue; - - if (track.tpcNClsFound() < cfgTrackPIDsettings->get(i, "minTPCnCls")) - continue; - if (((track.tpcNClsCrossedRows() < cfgTrackPIDsettings->get(i, "minTPCnClsCrossedRows")) || track.tpcNClsCrossedRows() < cfgtpcNClsFindable * track.tpcNClsFindable()) && cfgTPCNClsCrossedRowsRequire) - continue; - if (track.tpcChi2NCl() > cfgTrackPIDsettings->get(i, "maxTPCchi2") && cfgmaxTPCchi2Require) - continue; - if (track.tpcChi2NCl() < cfgTrackPIDsettings->get(i, "minTPCchi2") && cfgminTPCchi2Require) - continue; - if (track.itsNCls() < cfgTrackPIDsettings->get(i, "minITSnCls") && cfgminITSnClsRequire) - continue; - double cosheta = std::cosh(track.eta()); - if ((track.itsNCls() / cosheta) < cfgTrackPIDsettings->get(i, "minITSnClscos") && cfgminITSnClscosRequire) - continue; - if ((track.itsNClsInnerBarrel() < cfgTrackPIDsettings->get(i, "minReqClusterITSib")) && cfgminReqClusterITSibRequire) - continue; - if (track.itsChi2NCl() > cfgTrackPIDsettings->get(i, "maxITSchi2") && cfgmaxITSchi2Require) - continue; - if (getMeanItsClsSize(track) < cfgTrackPIDsettings->get(i, "minITSclsSize") && cfgminGetMeanItsClsSizeRequire) - continue; - float tpcNsigma = getTPCnSigma(track, primaryParticles.at(i)); - if ((std::abs(tpcNsigma) > cfgTrackPIDsettings->get(i, "maxTPCnSigma")) && cfgmaxTPCnSigmaRequire) - continue; - - int decayType = -999; // 0 = primary, 1 = weak decay, 2 = material - if (matchedMCParticle.isPhysicalPrimary()) { - // ---- Primary particles ---- - decayType = 0; - if (matchedMCParticle.has_mothers()) { - for (auto& motherparticle : matchedMCParticle.mothers_as()) { - if (std::find(hfMothCodes.begin(), hfMothCodes.end(), std::abs(motherparticle.pdgCode())) != hfMothCodes.end()) { - decayType = 1; - break; - } - } - } - } else if (matchedMCParticle.getProcess() == TMCProcess::kPDecay) { - // ---- Secondary from weak decay ---- - if (!matchedMCParticle.has_mothers()) { - continue; // skip secondaries from weak decay without mothers - } - decayType = 1; - - // Check if it's from HF decay - for (auto& motherparticle : matchedMCParticle.mothers_as()) { - if (std::find(hfMothCodes.begin(), hfMothCodes.end(), std::abs(motherparticle.pdgCode())) != hfMothCodes.end()) { - break; - } - } - } else { - // ---- Secondary from material interaction ---- - decayType = 2; - } - - float ptDCA; - - if (i == he3 || i == he4) { - ptDCA = ptReco; - } else { - ptDCA = 2 * ptReco; - } - - if (pdg == particlePdgCodes.at(i)) { // He3 - histomc.fill(HIST("DCAxy_vs_pT_total"), ptDCA, track.dcaXY(), collision.centFT0C()); - if (matchedMCParticle.isPhysicalPrimary()) { - histomc.fill(HIST("DCAxy_vs_pT_primary"), ptDCA, track.dcaXY(), collision.centFT0C()); - } - if (decayType == 2) { // Transport/Material - histomc.fill(HIST("DCAxy_vs_pT_transport"), ptDCA, track.dcaXY(), collision.centFT0C()); - } else if (decayType == 1) { // Weak decay (including HF) - histomc.fill(HIST("DCAxy_vs_pT_weakdecay"), ptDCA, track.dcaXY(), collision.centFT0C()); - } - } else if (pdg == -particlePdgCodes.at(i)) { // anti-He3 - histomc.fill(HIST("DCAxy_vs_pT_anti_total"), ptDCA, track.dcaXY(), collision.centFT0C()); - } - } - } - break; // Found the matching collision, break out of collision loop - } - } - } - } - PROCESS_SWITCH(NucleitpcPbPb, processDCA, "MC DCA analysis For secondary correction", false); - //=-=-=-==-=-=-==-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= - - void initCCDB(aod::BCsWithTimestamps::iterator const& bc) - { - if (mRunNumber == bc.runNumber()) { - return; - } - constexpr float kInvalidBField = -990.f; - auto run3grpTimestamp = bc.timestamp(); - dBz = 0; - o2::parameters::GRPObject* grpo = ccdb->getForTimeStamp(grpPath, run3grpTimestamp); - o2::parameters::GRPMagField* grpmag = 0x0; - if (grpo) { - o2::base::Propagator::initFieldFromGRP(grpo); - if (bField < kInvalidBField) { - // Fetch magnetic field from ccdb for current collision - dBz = grpo->getNominalL3Field(); - LOG(info) << "Retrieved GRP for timestamp " << run3grpTimestamp << " with magnetic field of " << dBz << " kZG"; - } else { - dBz = bField; - } - } else { - grpmag = ccdb->getForTimeStamp(grpmagPath, run3grpTimestamp); - if (!grpmag) { - LOG(fatal) << "Got nullptr from CCDB for path " << grpmagPath << " of object GRPMagField and " << grpPath << " of object GRPObject for timestamp " << run3grpTimestamp; - } - o2::base::Propagator::initFieldFromGRP(grpmag); - if (bField < kInvalidBField) { - // Fetch magnetic field from ccdb for current collision - dBz = std::lround(5.f * grpmag->getL3Current() / 30000.f); - LOG(info) << "Retrieved GRP for timestamp " << run3grpTimestamp << " with magnetic field of " << dBz << " kZG"; - } else { - dBz = bField; - } - } - mRunNumber = bc.runNumber(); - } - //---------------------------------------------------------------------------------------------------------------- - template - void initCollision(const T& collision) - { - collHasCandidate = false; - histos.fill(HIST("histNev"), 0.5); - collPassedEvSel = collision.sel8() && std::abs(collision.posZ()) < cfgZvertex; - occupancy = collision.trackOccupancyInTimeRange(); - if (collPassedEvSel) { - histos.fill(HIST("histNev"), 1.5); - } - primVtx.assign({collision.posX(), collision.posY(), collision.posZ()}); - cents.assign({collision.centFT0A(), collision.centFT0C(), collision.centFT0M()}); - } - //---------------------------------------------------------------------------------------------------------------- - template - void fillhmass(T const& track, int species, float cent) - { - if (!track.hasTOF() || !cfgFillmass) - return; - - float beta{o2::pid::tof::Beta::GetBeta(track)}; - const float eps = 1e-6f; - if (beta < eps || beta > 1.0f - eps) - return; - float charge = (species == he3 || species == he4) ? 2.f : 1.f; - float p = getRigidity(track); // assuming this is the momentum from inner TPC - float massTOF = p * charge * std::sqrt(1.f / (beta * beta) - 1.f); - float massDiff = 0.0; - if (species != he4) { - if (cfgmass2) { - // Compare squared masses - massDiff = (massTOF * massTOF) / (charge * charge); - } else { - // Compare linear masses - massDiff = massTOF / charge; - } - } - if (species == he4) { - if (cfghe3massrejreq && (massTOF * massTOF > cfgminmassrejection && massTOF * massTOF < cfgmaxmassrejection)) - return; - if (cfgmass2) { - // Compare squared masses - massDiff = (massTOF * massTOF) / (charge * charge); - } else { - // Compare linear masses - massDiff = massTOF / charge; - } - } - - float ptMomn; - setTrackParCov(track, mTrackParCov); - mTrackParCov.setPID(track.pidForTracking()); - ptMomn = (species == he3 || species == he4) ? 2 * mTrackParCov.getPt() : mTrackParCov.getPt(); - if (track.sign() > 0) { - hmass[2 * species]->Fill(ptMomn, massDiff, cent); - } else if (track.sign() < 0) { - hmass[2 * species + 1]->Fill(ptMomn, massDiff, cent); - } - } - //---------------------------------------------------------------------------------------------------------------- - template - void fillhmassnsigma(T const& track, int species, float cent) - { - if (!track.hasTOF() || !cfgFillmassnsigma) - return; - - float beta{o2::pid::tof::Beta::GetBeta(track)}; - const float eps = 1e-6f; - if (beta < eps || beta > 1.0f - eps) - return; - - float charge = (species == he3 || species == he4) ? 2.f : 1.f; - float p = getRigidity(track); - float massTOF = p * charge * std::sqrt(1.f / (beta * beta) - 1.f); - - // Calculate mass²/z² - float massSquareOverChargeSquare = (massTOF * massTOF) / (charge * charge); - - // Apply helium-4 rejection if needed - if (species == he4) { - if (cfghe3massrejreq && (massTOF * massTOF > cfgminmassrejection && massTOF * massTOF < cfgmaxmassrejection)) - return; - } - - float ptMomn; - setTrackParCov(track, mTrackParCov); - mTrackParCov.setPID(track.pidForTracking()); - ptMomn = (species == he3 || species == he4) ? 2 * mTrackParCov.getPt() : mTrackParCov.getPt(); - - // Fill histogram: mass²/z² vs pT vs centrality - if (track.sign() > 0) { - hmassnsigma[2 * species]->Fill(ptMomn, massSquareOverChargeSquare, cent); - } else if (track.sign() < 0) { - hmassnsigma[2 * species + 1]->Fill(ptMomn, massSquareOverChargeSquare, cent); - } - } - //---------------------------------------------------------------------------------------------------------------- - template - float getTPCnSigma(T const& track, PrimParticles& particle) - { - const float rigidity = getRigidity(track); - if (!track.hasTPC()) - return -999; - if (particle.name == "pion" && cfgTrackPIDsettings->get("pion", "useBBparams") < 1) - return cfgTrackPIDsettings->get("pion", "useBBparams") == 0 ? track.tpcNSigmaPi() : 0; - if (particle.name == "proton" && cfgTrackPIDsettings->get("proton", "useBBparams") < 1) - return cfgTrackPIDsettings->get("proton", "useBBparams") == 0 ? track.tpcNSigmaPr() : 0; - if (particle.name == "deuteron" && cfgTrackPIDsettings->get("deuteron", "useBBparams") < 1) - return cfgTrackPIDsettings->get("deuteron", "useBBparams") == 0 ? track.tpcNSigmaDe() : 0; - if (particle.name == "triton" && cfgTrackPIDsettings->get("triton", "useBBparams") < 1) - return cfgTrackPIDsettings->get("triton", "useBBparams") == 0 ? track.tpcNSigmaTr() : 0; - if (particle.name == "helion" && cfgTrackPIDsettings->get("helion", "useBBparams") < 1) - return cfgTrackPIDsettings->get("helion", "useBBparams") == 0 ? track.tpcNSigmaHe() : 0; - if (particle.name == "alpha" && cfgTrackPIDsettings->get("alpha", "useBBparams") < 1) - return cfgTrackPIDsettings->get("alpha", "useBBparams") == 0 ? track.tpcNSigmaAl() : 0; - - double expBethe{common::BetheBlochAleph(static_cast(particle.charge * rigidity / particle.mass), particle.betheParams[0], particle.betheParams[1], particle.betheParams[2], particle.betheParams[3], particle.betheParams[4])}; - double expSigma{expBethe * particle.resolution}; - float sigmaTPC = static_cast((track.tpcSignal() - expBethe) / expSigma); - return sigmaTPC; - } - //---------------------------------------------------------------------------------------------------------------- - template - float getITSnSigma(T const& track, PrimParticles& particle) - { - if (!track.hasITS()) - return -999; - o2::aod::ITSResponse itsResponse; - if (particle.name == "pion") - return itsResponse.nSigmaITS(track); - if (particle.name == "proton") - return itsResponse.nSigmaITS(track); - if (particle.name == "deuteron") - return itsResponse.nSigmaITS(track); - if (particle.name == "triton") - return itsResponse.nSigmaITS(track); - if (particle.name == "helion") - return itsResponse.nSigmaITS(track); - if (particle.name == "alpha") - return itsResponse.nSigmaITS(track); - return -999; // fallback if no match - } - //---------------------------------------------------------------------------------------------------------------- - template - float getMeanItsClsSize(T const& track) - { - int sum = 0, n = 0; - constexpr int kNITSLayers = 8; - for (int i = 0; i < kNITSLayers; i++) { - sum += (track.itsClusterSizes() >> (4 * i) & 15); - if (track.itsClusterSizes() >> (4 * i) & 15) - n++; - } - return n > 0 ? static_cast(sum) / n : 0.f; - } - //---------------------------------------------------------------------------------------------------------------- - template - float getRigidity(T const& track) - { - if (!cfgRigidityCorrection) - return track.tpcInnerParam(); - bool hePID = track.pidForTracking() == o2::track::PID::Helium3 || track.pidForTracking() == o2::track::PID::Alpha; - return hePID ? track.tpcInnerParam() / 2 : track.tpcInnerParam(); - } - - //---------------------------------------------------------------------------------------------------------------- - template - float getRapidity(T const& track, int species) - { - using PtEtaPhiMVector = ROOT::Math::LorentzVector>; - double momn; - int speciesHe3 = 4; - int speciesHe4 = 5; - if (species == speciesHe3 || species == speciesHe4) { - momn = 2 * track.pt(); - } else { - momn = track.pt(); - } - PtEtaPhiMVector lorentzVectorParticle(momn, track.eta(), track.phi(), particleMasses[species]); - return lorentzVectorParticle.Rapidity(); - } -}; // end of the task here -//---------------------------------------------------------------------------------------------------------------- -//---------------------------------------------------------------------------------------------------------------- -WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) -{ - return WorkflowSpec{ - adaptAnalysisTask(cfgc)}; -}