diff --git a/PWGLF/Tasks/Nuspex/nucleitpcpbpb.cxx b/PWGLF/Tasks/Nuspex/nucleitpcpbpb.cxx index 110afcc0a19..f56eee2ebc9 100644 --- a/PWGLF/Tasks/Nuspex/nucleitpcpbpb.cxx +++ b/PWGLF/Tasks/Nuspex/nucleitpcpbpb.cxx @@ -8,11 +8,13 @@ // 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 -/// \author Jaideep Tanwar +/// \note under work /// +/// \author Jaideep Tanwar , Panjab University + #include "Common/Core/PID/PIDTOF.h" #include "Common/Core/PID/TPCPIDResponse.h" #include "Common/Core/RecoDecay.h" @@ -47,6 +49,7 @@ using namespace o2::framework; using namespace o2::framework::expressions; using CollisionsFull = soa::Join; using TracksFull = soa::Join; +using CollisionsFullMC = soa::Join; //--------------------------------------------------------------------------------------------------------------------------------- namespace { @@ -57,22 +60,22 @@ static const std::vector particleMasses{o2::constants::physics::MassPion static const std::vector particleCharge{1, 1, 1, 1, 2, 2}; const int nBetheParams = 6; static const std::vector betheBlochParNames{"p0", "p1", "p2", "p3", "p4", "resolution"}; -constexpr double betheBlochDefault[nParticles][nBetheParams]{ +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 = 18; -static const std::vector trackPIDsettingsNames{"useBBparams", "minITSnCls", "minITSnClscos", "minTPCnCls", "maxTPCchi2", "minTPCchi2", "maxITSchi2", "minRigidity", "maxRigidity", "maxTPCnSigma", "TOFrequiredabove", "minTOFmass", "maxTOFmass", "maxDcaXY", "maxDcaZ", "minITSclsSize", "maxITSclsSize", "minTPCnClsCrossedRows"}; -constexpr double trackPIDsettings[nParticles][nTrkSettings]{ - {0, 0, 4, 60, 4.0, 0.5, 100, 0.15, 1.2, 2.5, -1, 0, 100, 2., 2., 0., 1000, 70}, - {1, 0, 4, 70, 4.0, 0.5, 100, 0.20, 4.0, 3.0, -1, 0, 100, 2., 2., 0., 1000, 70}, - {1, 0, 4, 70, 4.0, 0.5, 100, 0.50, 5.0, 3.0, -1, 0, 100, 2., 2., 0., 1000, 70}, - {1, 0, 4, 70, 4.0, 0.5, 100, 0.50, 5.0, 3.0, -1, 0, 100, 2., 2., 0., 1000, 70}, - {1, 0, 4, 75, 4.0, 0.5, 100, 0.50, 5.0, 5.0, -1, 0, 100, 2., 2., 0., 1000, 70}, - {1, 0, 4, 70, 4.0, 0.5, 100, 0.50, 5.0, 5.0, -1, 0, 100, 2., 2., 0., 1000, 70}}; +const int nTrkSettings = 19; +static const std::vector trackPIDsettingsNames{"useBBparams", "minITSnCls", "minITSnClscos", "minTPCnCls", "maxTPCchi2", "minTPCchi2", "maxITSchi2", "minRigidity", "maxRigidity", "maxTPCnSigma", "TOFrequiredabove", "minTOFmass", "maxTOFmass", "maxDcaXY", "maxDcaZ", "minITSclsSize", "maxITSclsSize", "minTPCnClsCrossedRows", "minReqClusterITSib"}; +constexpr double kTrackPIDSettings[nParticles][nTrkSettings]{ + {0, 0, 4, 60, 4.0, 0.5, 100, 0.15, 1.2, 2.5, -1, 0, 100, 2., 2., 0., 1000, 70, 1}, + {1, 0, 4, 70, 4.0, 0.5, 100, 0.20, 4.0, 3.0, -1, 0, 100, 2., 2., 0., 1000, 70, 1}, + {1, 0, 4, 70, 4.0, 0.5, 100, 0.50, 5.0, 3.0, -1, 0, 100, 2., 2., 0., 1000, 70, 1}, + {1, 0, 4, 70, 4.0, 0.5, 100, 0.50, 5.0, 3.0, -1, 0, 100, 2., 2., 0., 1000, 70, 1}, + {1, 0, 4, 75, 4.0, 0.5, 100, 0.50, 5.0, 5.0, -1, 0, 100, 2., 2., 0., 1000, 70, 1}, + {1, 0, 4, 70, 4.0, 0.5, 100, 0.50, 5.0, 5.0, -1, 0, 100, 2., 2., 0., 1000, 70, 1}}; struct PrimParticles { TString name; int pdgCode, charge; @@ -97,6 +100,7 @@ std::vector> hmass; struct NucleitpcPbPb { Preslice perCollision = aod::track_association::collisionId; HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject}; + HistogramRegistry histomc{"histomc", {}, OutputObjHandlingPolicy::AnalysisObject, true, true}; Configurable cfgDebug{"cfgDebug", 1, "debug level"}; // event Selections cuts Configurable removeITSROFrameBorder{"removeITSROFrameBorder", false, "Remove TF border"}; @@ -118,7 +122,9 @@ struct NucleitpcPbPb { Configurable cfgmaxTPCchi2Require{"cfgmaxTPCchi2Require", true, "Require maxTPCchi2 Cut"}; Configurable cfgminTPCchi2Require{"cfgminTPCchi2Require", true, "Require minTPCchi2 Cut"}; Configurable cfgminITSnClsRequire{"cfgminITSnClsRequire", false, "Require minITSnCls Cut"}; + Configurable cfgmccorrectionhe4Require{"cfgmccorrectionhe4Require", true, "MC correction for pp he4 particle"}; 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 cfgmaxITSnSigmaRequire{"cfgmaxITSnSigmaRequire", true, "Require maxITSnSigma Cut for helium"}; @@ -129,8 +135,11 @@ struct NucleitpcPbPb { Configurable cfgDCAwithptRequire{"cfgDCAwithptRequire", true, "Require DCA cuts with pt dependance"}; Configurable cfgDCAnopt{"cfgDCAnopt", true, "Require DCA cuts without pt dependance"}; Configurable cfgTwicemass{"cfgTwicemass", true, "multiply mass by its charge"}; - Configurable> cfgBetheBlochParams{"cfgBetheBlochParams", {betheBlochDefault[0], nParticles, nBetheParams, particleNames, betheBlochParNames}, "TPC Bethe-Bloch parameterisation for light nuclei"}; - Configurable> cfgTrackPIDsettings{"cfgTrackPIDsettings", {trackPIDsettings[0], nParticles, nTrkSettings, particleNames, trackPIDsettingsNames}, "track selection and PID criteria"}; + Configurable cfgRequirebetaplot{"cfgRequirebetaplot", true, "Require beta plot"}; + Configurable cfgRequireMCposZ{"cfgRequireMCposZ", true, "Require beta plot"}; + + 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 cfgFillDeDxWithCut{"cfgFillDeDxWithCut", true, "Fill with cut beth bloch"}; Configurable cfgFillnsigma{"cfgFillnsigma", false, "Fill n-sigma histograms"}; Configurable cfgFillmass{"cfgFillmass", false, "Fill mass histograms"}; @@ -141,6 +150,19 @@ struct NucleitpcPbPb { Configurable cfgtpcNClsFound{"cfgtpcNClsFound", 100.0f, "min. no. of tpcNClsFound"}; Configurable cfgitsNCls{"cfgitsNCls", 2.0f, "min. no. of itsNCls"}; o2::track::TrackParametrizationWithError mTrackParCov; + // Binning configuration + ConfigurableAxis axisMagField{"axisMagField", {10, -10., 10.}, "magnetic field"}; + ConfigurableAxis axisNev{"axisNev", {3, 0., 3.}, "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 axisVtxZ{"axisVtxZ", {120, -20, 20}, "z"}; + ConfigurableAxis axisImpt{"axisImpt", {100, 0, 20}, "impact parameter"}; + 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", {100, -10, 10}, "mass^{2}"}; + ConfigurableAxis nsigmaAxis{"nsigmaAxis", {160, -20, 20}, "n#sigma_{#pi^{+}}"}; // CCDB Service ccdb; Configurable bField{"bField", -999, "bz field, -999 is automatic"}; @@ -153,12 +175,12 @@ struct NucleitpcPbPb { std::vector primaryParticles; std::vector primVtx, cents; - bool collHasCandidate, collPassedEvSel; + bool collHasCandidate, collPassedEvSel, collPassedEvSelMc; int mRunNumber, occupancy; float dBz, momn; TRandom3 rand; - float He3 = 4; - float He4 = 5; + float he3 = 4; + float he4 = 5; //---------------------------------------------------------------------------------------------------------------- void init(InitContext const&) { @@ -172,31 +194,21 @@ struct NucleitpcPbPb { 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)); } - // define histogram axes - const AxisSpec axisMagField{10, -10., 10., "magnetic field"}; - const AxisSpec axisNev{3, 0., 3., "Number of events"}; - const AxisSpec axisRigidity{4000, -10., 10., "#it{p}^{TPC}/#it{z}"}; - const AxisSpec axisdEdx{4000, 0, 4000, "d#it{E}/d#it{x}"}; - const AxisSpec axisCent{100, 0, 100, "centrality"}; - const AxisSpec axisVtxZ{100, -20, 20, "z"}; - const AxisSpec ptAxis{200, 0, 20, "#it{p}_{T} (GeV/#it{c})"}; - const AxisSpec axiseta{100, -1, 1, "eta"}; - const AxisSpec axisrapidity{100, -2, 2, "rapidity"}; - const AxisSpec axismass{100, -10, 10, "mass^{2}"}; - AxisSpec nsigmaAxis = {160, -20, 20, "n#sigma_{#pi^{+}}"}; // create histograms + if (doprocessData) { + histos.add("histMagField", "histMagField", kTH1F, {axisMagField}); + 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("histeta", "histeta", kTH1F, {axiseta}); - histos.add("histCentFTOC_cut", "histCentFTOC_cut", kTH1F, {axisCent}); - histos.add("Tof_signal", "Tof_signal", kTH2F, {axisRigidity, {4000, 0.2, 1.2, "#beta"}}); + histos.add("Tofsignal", "Tofsignal", kTH2F, {axisRigidity, {4000, 0.2, 1.2, "#beta"}}); histos.add("histDcaZVsPtData_particle", "dcaZ vs Pt (particle)", HistType::kTH2F, {{1000, 0, 20}, {1000, -2.5, 2.5, "dca"}}); histos.add("histDcaXYVsPtData_particle", "dcaXY vs Pt (particle)", HistType::kTH2F, {{1000, 0, 20}, {1000, -2.0, 2.0, "dca"}}); histos.add("histDcaZVsPtData_antiparticle", "dcaZ vs Pt (antiparticle)", HistType::kTH2F, {{1000, 0, 20}, {1000, -2.5, 2.5, "dca"}}); histos.add("histDcaXYVsPtData_antiparticle", "dcaXY vs Pt (antiparticle)", HistType::kTH2F, {{1000, 0, 20}, {1000, -2.0, 2.0, "dca"}}); - histos.add("histMagField", "histMagField", kTH1F, {axisMagField}); - histos.add("histNev", "histNev", kTH1F, {axisNev}); - histos.add("histVtxZ", "histVtxZ", kTH1F, {axisVtxZ}); - histos.add("histCentFT0C", "histCentFT0C", kTH1F, {axisCent}); - histos.add("histCentFT0M", "histCentFT0M", kTH1F, {axisCent}); hDeDx.resize(2 * nParticles + 2); hNsigmaPt.resize(2 * nParticles + 2); hmass.resize(2 * nParticles + 2); @@ -217,6 +229,27 @@ struct NucleitpcPbPb { hmass[2 * i + 1] = histos.add(Form("histmass_ptanti/histmass_%s", histName.Data()), ";p_T{TPC} (GeV/#it{c}); mass^{2}", HistType::kTH2F, {ptAxis, axismass}); } } + + if (doprocessMC) { + histomc.add("histVtxZgen", "histVtxZgen", kTH1F, {axisVtxZ}); + histomc.add("histEtagen", "histEtagen", kTH1F, {axiseta}); + histomc.add("histPtgenHe3", "histPtgenHe3", kTH1F, {ptAxis}); + histomc.add("histPtgenAntiHe3", "histPtgenAntiHe3", kTH1F, {ptAxis}); + histomc.add("histPtgenHe4", "histPtgenHe4", kTH1F, {ptAxis}); + histomc.add("histPtgenAntiHe4", "histPtgenAntiHe4", kTH1F, {ptAxis}); + // Reconstrcuted eta + histomc.add("histNevReco", "histNevReco", kTH1F, {axisNev}); + histomc.add("histVtxZReco", "histVtxZReco", kTH1F, {axisVtxZ}); + histomc.add("histCentFT0CReco", "histCentFT0CReco", kTH1F, {axisCent}); + histomc.add("histCentFT0MReco", "histCentFT0MReco", kTH1F, {axisCent}); + histomc.add("histPtRecoHe3", "histPtgenHe3", kTH1F, {ptAxis}); + histomc.add("histPtRecoAntiHe3", "histPtgenAntiHe3", kTH1F, {ptAxis}); + histomc.add("histPtRecoHe4", "histPtgenHe4", kTH1F, {ptAxis}); + histomc.add("histPtRecoAntiHe4", "histPtgenAntiHe4", kTH1F, {ptAxis}); + 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("histDeltaPtVsPtGenHe4", " delta pt vs pt rec", HistType::kTH2F, {{1000, 0, 10}, {1000, -0.5, 0.5, "p_{T}(reco) - p_{T}(gen);p_{T}(reco)"}}); + } } //---------------------------------------------------------------------------------------------------------------- void findprimaryParticles(aod::TrackAssoc const& tracksByColl, TracksFull const& tracks) @@ -251,6 +284,8 @@ struct NucleitpcPbPb { 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) @@ -262,9 +297,9 @@ struct NucleitpcPbPb { float ptMomn; setTrackParCov(track, mTrackParCov); mTrackParCov.setPID(track.pidForTracking()); - ptMomn = (i == He3 || i == He4) ? 2 * mTrackParCov.getPt() : mTrackParCov.getPt(); + ptMomn = (i == he3 || i == he4) ? 2 * mTrackParCov.getPt() : mTrackParCov.getPt(); bool insideDCAxy = (std::abs(track.dcaXY()) <= (cfgTrackPIDsettings->get(i, "maxDcaXY") * (0.0105f + 0.0350f / std::pow(ptMomn, 1.1f)))); // o2-linter: disable=magic-number (To be checked) - if ((!(insideDCAxy) || std::abs(track.dcaZ()) > DCAzSigma(ptMomn, cfgTrackPIDsettings->get(i, "maxDcaZ"))) && cfgDCAwithptRequire) + if ((!(insideDCAxy) || std::abs(track.dcaZ()) > dcazSigma(ptMomn, cfgTrackPIDsettings->get(i, "maxDcaZ"))) && cfgDCAwithptRequire) continue; if (track.sign() > 0) { histos.fill(HIST("histDcaZVsPtData_particle"), ptMomn, track.dcaZ()); @@ -286,18 +321,20 @@ struct NucleitpcPbPb { if (!track.hasTOF() && cfgTrackPIDsettings->get(i, "TOFrequiredabove") < 1) continue; float beta{o2::pid::tof::Beta::GetBeta(track)}; - float charge{1.f + static_cast(i == He3 || i == He4)}; + float charge{1.f + static_cast(i == he3 || i == he4)}; float tofMasses = getRigidity(track) * charge * std::sqrt(1.f / (beta * beta) - 1.f); if ((getRigidity(track) > cfgTrackPIDsettings->get(i, "TOFrequiredabove") && (tofMasses < cfgTrackPIDsettings->get(i, "minTOFmass") || tofMasses > cfgTrackPIDsettings->get(i, "maxTOFmass"))) && cfgmassRequire) continue; fillhmass(track, i); } histos.fill(HIST("histeta"), track.eta()); + if (cfgRequirebetaplot) { + histos.fill(HIST("Tofsignal"), getRigidity(track), o2::pid::tof::Beta::GetBeta(track)); + } filldedx(track, nParticles); } // track loop } //---------------------------------------------------------------------------------------------------------------- - //---------------------------------------------------------------------------------------------------------------- void processData(CollisionsFull const& collisions, TracksFull const& tracks, aod::BCsWithTimestamps const&, aod::TrackAssoc const& tracksColl) { for (const auto& collision : collisions) { @@ -326,8 +363,209 @@ struct NucleitpcPbPb { continue; } } - PROCESS_SWITCH(NucleitpcPbPb, processData, "data analysis", true); + PROCESS_SWITCH(NucleitpcPbPb, processData, "data analysis", false); //---------------------------------------------------------------------------------------------------------------- + // MC particles + //-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= + struct McCollInfo { + bool passedEvSel = false; + }; + std::vector mcCollInfos; + void processMC(CollisionsFullMC const& collisions, aod::McCollisions const& mcCollisions, TracksFull const& tracks, aod::BCsWithTimestamps const& bcs, aod::McParticles const& particlesMC, aod::McTrackLabels const& trackLabelsMC, aod::McCollisionLabels const& collLabels, aod::TrackAssoc const& tracksColl, CollisionsFull const& colls) + { + (void)colls; + (void)collLabels; + (void)bcs; + mcCollInfos.clear(); + mcCollInfos.resize(mcCollisions.size()); + // ----------------------------- Generated particles loop ----------------------------- + for (auto const& mcColl : mcCollisions) { + if (std::abs(mcColl.posZ()) > cfgZvertex) + continue; + histomc.fill(HIST("histVtxZgen"), mcColl.posZ()); + for (auto const& mcParticle : particlesMC) { + if (mcParticle.mcCollisionId() != mcColl.globalIndex()) + continue; + if (!mcParticle.isPhysicalPrimary()) + continue; + int pdgCode = mcParticle.pdgCode(); + if (std::abs(mcParticle.y()) > cfgCutRapidity) + continue; + if (std::abs(mcParticle.eta()) > cfgCutEta) + continue; + histomc.fill(HIST("histEtagen"), mcParticle.eta()); + float ptScaled = mcParticle.pt(); + if (pdgCode == particlePdgCodes.at(4)) { + histomc.fill(HIST("histPtgenHe3"), ptScaled); + } else if (pdgCode == -particlePdgCodes.at(4)) { + histomc.fill(HIST("histPtgenAntiHe3"), ptScaled); + } else if (pdgCode == particlePdgCodes.at(5)) { + histomc.fill(HIST("histPtgenHe4"), ptScaled); + } else if (pdgCode == -particlePdgCodes.at(5)) { + histomc.fill(HIST("histPtgenAntiHe4"), ptScaled); + } + } // mc track loop generated + } // mc collision loop generated + // ----------------------------- Reconstructed track loop ----------------------------- + for (auto const& collision : collisions) { + auto mcCollIdx = collision.mcCollisionId(); + if (mcCollIdx < 0 || mcCollIdx >= mcCollisions.size()) + continue; + auto bc = collision.bc_as(); + initCCDB(bc); + collHasCandidate = false; + histomc.fill(HIST("histNevReco"), 0.5); + collPassedEvSel = collision.sel8() && std::abs(collision.posZ()) < cfgZvertex; + occupancy = collision.trackOccupancyInTimeRange(); + if (!collPassedEvSel) + continue; + histomc.fill(HIST("histNevReco"), 1.5); + histomc.fill(HIST("histVtxZReco"), collision.posZ()); + histomc.fill(HIST("histCentFT0CReco"), collision.centFT0C()); + histomc.fill(HIST("histCentFT0MReco"), collision.centFT0M()); + if (collision.centFT0C() > centcut) + continue; + if (removeITSROFrameBorder && !collision.selection_bit(aod::evsel::kNoITSROFrameBorder)) + 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; + // Loop through TrackAssoc and check if this track belongs to the current collision + for (const auto& assoc : tracksColl) { + if (assoc.collisionId() != collision.globalIndex()) + continue; + const auto& track = tracks.rawIteratorAt(assoc.trackId()); + auto labelRow = trackLabelsMC.iteratorAt(track.globalIndex()); + int label = labelRow.mcParticleId(); + if (label < 0 || label >= particlesMC.size()) + continue; + auto const& matchedMCParticle = particlesMC.iteratorAt(label); + int pdg = matchedMCParticle.pdgCode(); + if (!track.isPVContributor() && cfgUsePVcontributors) + continue; + if (!track.hasITS() && cfgITSrequire) + continue; + if (!track.hasTPC() && cfgTPCrequire) + continue; + if (!track.passedITSRefit() && cfgPassedITSRefit) + continue; + if (!track.passedTPCRefit() && cfgPassedTPCRefit) + continue; + if (std::abs(track.eta()) > cfgCutEta) + continue; + if (!matchedMCParticle.isPhysicalPrimary()) + continue; + for (size_t i = 0; i < primaryParticles.size(); i++) { + if (std::abs(pdg) != std::abs(particlePdgCodes.at(i))) + continue; + if (std::abs(getRapidity(track, i)) > cfgCutRapidity && cfgRapidityRequire) + continue; + if (track.tpcNClsFound() < cfgTrackPIDsettings->get(i, "minTPCnCls") && cfgTPCNClsfoundRequire) + continue; + if (((track.tpcNClsCrossedRows() < cfgTrackPIDsettings->get(i, "minTPCnClsCrossedRows")) || track.tpcNClsCrossedRows() < 0.8 * track.tpcNClsFindable()) && cfgTPCNClsCrossedRowsRequire) // o2-linter: disable=magic-number (To be checked) + 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; + if (getMeanItsClsSize(track) > cfgTrackPIDsettings->get(i, "maxITSclsSize") && cfgmaxGetMeanItsClsSizeRequire) + continue; + if ((getRigidity(track) < cfgTrackPIDsettings->get(i, "minRigidity") || getRigidity(track) > cfgTrackPIDsettings->get(i, "maxRigidity")) && cfgRigidityCutRequire) + continue; + float ptMomn; + setTrackParCov(track, mTrackParCov); + mTrackParCov.setPID(track.pidForTracking()); + ptMomn = (i == he3 || i == he4) ? 2 * mTrackParCov.getPt() : mTrackParCov.getPt(); + bool insideDCAxy = (std::abs(track.dcaXY()) <= (cfgTrackPIDsettings->get(i, "maxDcaXY") * (0.0105f + 0.0350f / std::pow(ptMomn, 1.1f)))); // o2-linter: disable=magic-number (To be checked) + if ((!(insideDCAxy) || std::abs(track.dcaZ()) > dcazSigma(ptMomn, cfgTrackPIDsettings->get(i, "maxDcaZ"))) && cfgDCAwithptRequire) + continue; + if (track.sign() > 0) { + histos.fill(HIST("histDcaZVsPtData_particle"), ptMomn, track.dcaZ()); + histos.fill(HIST("histDcaXYVsPtData_particle"), ptMomn, track.dcaXY()); + } + if (track.sign() < 0) { + histos.fill(HIST("histDcaZVsPtData_antiparticle"), ptMomn, track.dcaZ()); + histos.fill(HIST("histDcaXYVsPtData_antiparticle"), ptMomn, track.dcaXY()); + } + float tpcNsigma = getTPCnSigma(track, primaryParticles.at(i)); + if ((std::abs(tpcNsigma) > cfgTrackPIDsettings->get(i, "maxTPCnSigma")) && cfgmaxTPCnSigmaRequire) + continue; + float itsSigma = getITSnSigma(track, primaryParticles.at(i)); + if ((std::abs(itsSigma) > cfgITSnsigma) && cfgmaxITSnSigmaRequire) + continue; + fillnsigma(track, i); + filldedx(track, i); + if (!track.hasTOF() && cfgTrackPIDsettings->get(i, "TOFrequiredabove") < 1) + continue; + float beta{o2::pid::tof::Beta::GetBeta(track)}; + float charge{1.f + static_cast(i == he3 || i == he4)}; + float tofMasses = getRigidity(track) * charge * std::sqrt(1.f / (beta * beta) - 1.f); + if ((getRigidity(track) > cfgTrackPIDsettings->get(i, "TOFrequiredabove") && (tofMasses < cfgTrackPIDsettings->get(i, "minTOFmass") || tofMasses > cfgTrackPIDsettings->get(i, "maxTOFmass"))) && cfgmassRequire) + continue; + fillhmass(track, i); + } + histos.fill(HIST("histeta"), track.eta()); + if (cfgRequirebetaplot) { + histos.fill(HIST("Tofsignal"), getRigidity(track), o2::pid::tof::Beta::GetBeta(track)); + } + filldedx(track, nParticles); + /*----------------------------------------------------------------------------------------------------------------*/ + 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(); + if (pdg == -particlePdgCodes.at(5) && cfgmccorrectionhe4Require) { + ptReco = ptReco + 0.00765 + 0.503791 * std::exp(-1.10517 * ptReco); + } + + if (pdg == -particlePdgCodes.at(4) && cfgmccorrectionhe4Require) { + int pidGuess = track.pidForTracking(); + int antitriton = 6; + if (pidGuess == antitriton) { + ptReco = ptReco - 0.464215 + 0.195771 * ptReco - 0.0183111 * ptReco * ptReco; + // LOG(info) << "we have he3" << pidGuess; + } + } + float ptGen = matchedMCParticle.pt(); + float deltaPt = ptReco - ptGen; + + if (pdg == -particlePdgCodes.at(4)) { + histomc.fill(HIST("histDeltaPtVsPtGen"), ptReco, deltaPt); + histomc.fill(HIST("histPIDtrack"), ptReco, track.pidForTracking()); + } + if (pdg == -particlePdgCodes.at(5)) { + histomc.fill(HIST("histDeltaPtVsPtGenHe4"), ptReco, deltaPt); + } + if (pdg == particlePdgCodes.at(4)) { + histomc.fill(HIST("histPtRecoHe3"), ptReco); + } else if (pdg == -particlePdgCodes.at(4)) { + histomc.fill(HIST("histPtRecoAntiHe3"), ptReco); + } else if (pdg == particlePdgCodes.at(5)) { + histomc.fill(HIST("histPtRecoHe4"), ptReco); + } else if (pdg == -particlePdgCodes.at(5)) { + histomc.fill(HIST("histPtRecoAntiHe4"), ptReco); + } + } // Track loop + } // Collision loop + } + PROCESS_SWITCH(NucleitpcPbPb, processMC, "MC reco+gen analysis", true); + //-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-=-= void initCCDB(aod::BCsWithTimestamps::iterator const& bc) { if (mRunNumber == bc.runNumber()) { @@ -404,12 +642,12 @@ struct NucleitpcPbPb { float ptMomn; setTrackParCov(track, mTrackParCov); mTrackParCov.setPID(track.pidForTracking()); - ptMomn = (i == He3 || i == He4) ? 2 * mTrackParCov.getPt() : mTrackParCov.getPt(); + ptMomn = (i == he3 || i == he4) ? 2 * mTrackParCov.getPt() : mTrackParCov.getPt(); if (track.sign() > 0) { hNsigmaPt[2 * species]->Fill(ptMomn, tpcNsigma); } if (track.sign() < 0) { - hNsigmaPt[2 * species]->Fill(ptMomn, tpcNsigma); + hNsigmaPt[2 * species + 1]->Fill(ptMomn, tpcNsigma); } } } @@ -424,7 +662,7 @@ struct NucleitpcPbPb { float beta{o2::pid::tof::Beta::GetBeta(track)}; if (beta <= 0.f || beta >= 1.f) return; - float charge = (species == He3 || species == He4) ? 2.f : 1.f; + 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); // get PDG mass @@ -433,14 +671,13 @@ struct NucleitpcPbPb { float ptMomn; setTrackParCov(track, mTrackParCov); mTrackParCov.setPID(track.pidForTracking()); - ptMomn = (species == He3 || species == He4) ? 2 * mTrackParCov.getPt() : mTrackParCov.getPt(); + ptMomn = (species == he3 || species == he4) ? 2 * mTrackParCov.getPt() : mTrackParCov.getPt(); if (track.sign() > 0) { hmass[2 * species]->Fill(ptMomn, massDiff * massDiff); } else if (track.sign() < 0) { hmass[2 * species + 1]->Fill(ptMomn, massDiff * massDiff); } } - //---------------------------------------------------------------------------------------------------------------- template float getTPCnSigma(T const& track, PrimParticles& particle) @@ -487,7 +724,6 @@ struct NucleitpcPbPb { return itsResponse.nSigmaITS(track); return -999; // fallback if no match } - //---------------------------------------------------------------------------------------------------------------- template float getMeanItsClsSize(T const& track) @@ -510,7 +746,8 @@ struct NucleitpcPbPb { bool hePID = track.pidForTracking() == o2::track::PID::Helium3 || track.pidForTracking() == o2::track::PID::Alpha; return hePID ? track.tpcInnerParam() / 2 : track.tpcInnerParam(); } - float DCAzSigma(double pt, float dcasigma) + //---------------------------------------------------------------------------------------------------------------- + float dcazSigma(double pt, float dcasigma) { float invPt = 1.f / pt; return (5.00000e-04 + 8.73690e-03 * invPt + 9.62329e-04 * invPt * invPt) * dcasigma; // o2-linter: disable=magic-number (To be checked)