diff --git a/PWGLF/TableProducer/Nuspex/nucleiSpectra.cxx b/PWGLF/TableProducer/Nuspex/nucleiSpectra.cxx index 1835015f629..90dfcee5aed 100644 --- a/PWGLF/TableProducer/Nuspex/nucleiSpectra.cxx +++ b/PWGLF/TableProducer/Nuspex/nucleiSpectra.cxx @@ -51,9 +51,10 @@ #include "Framework/runDataProcessing.h" #include "ReconstructionDataFormats/Track.h" -#include "Math/Vector4D.h" -#include "TRandom3.h" +#include #include +#include // for PDG codes +#include #include #include @@ -118,6 +119,7 @@ struct NucleusCandidateFlow { namespace nuclei { +constexpr int nITSlayers = 7; constexpr double bbMomScalingDefault[5][2]{ {1., 1.}, {1., 1.}, @@ -272,16 +274,26 @@ struct nucleiSpectra { Configurable cfgCompensatePIDinTracking{"cfgCompensatePIDinTracking", false, "If true, divide tpcInnerParam by the electric charge"}; + struct : o2::framework::ConfigurableGroup { + std::string prefix{"cfgTrackCut"}; + Configurable> DcaMax{"DcaMax", {nuclei::DCAcutDefault[0], 5, 2, nuclei::names, nuclei::nDCAConfigName}, "Max DCAxy and DCAz for light nuclei"}; + Configurable EtaMax{"EtaMax", 0.8f, "Max Eta for tracks"}; + Configurable ITSnClusMin{"ITSnClsMin", 5, "Minimum number of ITS clusters"}; + Configurable ITSchi2ClusMax{"ITSchi2ClusMax", 36.f, "Max ITS Chi2 per cluster"}; + Configurable RapidityMax{"RapidityMax", 1.f, "Maximum rapidity for tracks"}; + Configurable RapidityMin{"RapidityMin", -1.f, "Minimum rapidity for tracks"}; + Configurable RapidityToggle{"RapidityToggle", false, "If true, use rapidity cuts"}; + Configurable TPCchi2ClusMax{"TPCchi2ClusMax", 4.f, "Max TPC Chi2 per cluster"}; + Configurable TPCnCrossedRowsMin{"TPCnCrossedRowsMin", 70, "Minimum number of TPC crossed rows"}; + Configurable TPCnCrossedRowsOverFindableMin{"TPCnCrossedRowsOverFindableMin", 0.8f, "Minimum ratio of crossed rows over findable clusters"}; + Configurable TPCnClsMin{"TPCnClsMin", 80, "Minimum number of TPC clusters"}; + Configurable TPCrigidityMin{"TPCrigidityMin", 0.5f, "Minimum TPC rigidity for tracks"}; + Configurable> TPCnSigmaMax{"TPCnSigmaMax", {nuclei::nSigmaTPCdefault[0], 5, 2, nuclei::names, nuclei::nSigmaConfigName}, "TPC nsigma selection for light nuclei"}; + + } cfgTrackCut; Configurable cfgCentralityEstimator{"cfgCentralityEstimator", 0, "Centrality estimator (FV0A: 0, FT0M: 1, FT0A: 2, FT0C: 3)"}; Configurable cfgCMrapidity{"cfgCMrapidity", 0.f, "Rapidity of the center of mass (only for p-Pb)"}; Configurable cfgCutVertex{"cfgCutVertex", 10.0f, "Accepted z-vertex range"}; - Configurable cfgCutEta{"cfgCutEta", 0.8f, "Eta range for tracks"}; - Configurable cfgCutTpcMom{"cfgCutTpcMom", 0.2f, "Minimum TPC momentum for tracks"}; - Configurable cfgCutRapidityMin{"cfgCutRapidityMin", -1., "Minimum rapidity for tracks"}; - Configurable cfgCutRapidityMax{"cfgCutRapidityMax", 1., "Maximum rapidity for tracks"}; - Configurable cfgCutOnReconstructedRapidity{"cfgCutOnReconstructedRapidity", false, "Cut on reconstructed rapidity"}; - Configurable cfgCutNclusITS{"cfgCutNclusITS", 5, "Minimum number of ITS clusters"}; - Configurable cfgCutNclusTPC{"cfgCutNclusTPC", 70, "Minimum number of TPC clusters"}; Configurable cfgCutPtMinTree{"cfgCutPtMinTree", 0.2f, "Minimum track transverse momentum for tree saving"}; Configurable cfgCutPtMaxTree{"cfgCutPtMaxTree", 15.0f, "Maximum track transverse momentum for tree saving"}; @@ -289,8 +301,6 @@ struct nucleiSpectra { Configurable> cfgMomentumScalingBetheBloch{"cfgMomentumScalingBetheBloch", {nuclei::bbMomScalingDefault[0], 5, 2, nuclei::names, nuclei::chargeLabelNames}, "TPC Bethe-Bloch momentum scaling for light nuclei"}; Configurable> cfgBetheBlochParams{"cfgBetheBlochParams", {nuclei::betheBlochDefault[0], 5, 6, nuclei::names, nuclei::betheBlochParNames}, "TPC Bethe-Bloch parameterisation for light nuclei"}; - Configurable> cfgNsigmaTPC{"cfgNsigmaTPC", {nuclei::nSigmaTPCdefault[0], 5, 2, nuclei::names, nuclei::nSigmaConfigName}, "TPC nsigma selection for light nuclei"}; - Configurable> cfgDCAcut{"cfgDCAcut", {nuclei::DCAcutDefault[0], 5, 2, nuclei::names, nuclei::nDCAConfigName}, "Max DCAxy and DCAz for light nuclei"}; Configurable> cfgDownscaling{"cfgDownscaling", {nuclei::DownscalingDefault[0], 5, 1, nuclei::names, nuclei::DownscalingConfigName}, "Fraction of kept candidates for light nuclei"}; Configurable> cfgTreeConfig{"cfgTreeConfig", {nuclei::TreeConfigDefault[0], 5, 2, nuclei::names, nuclei::treeConfigNames}, "Filtered trees configuration"}; Configurable cfgFillPairTree{"cfgFillPairTree", true, "Fill trees for pairs of light nuclei"}; @@ -331,7 +341,7 @@ struct nucleiSpectra { // CCDB options Configurable cfgMaterialCorrection{"cfgMaterialCorrection", static_cast(o2::base::Propagator::MatCorrType::USEMatCorrLUT), "Type of material correction"}; Configurable cfgCCDBurl{"ccdb-url", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; - Configurable cfgZorroCCDBpath{"cfgZorroCCDBpath", "/Users/m/mpuccio/EventFiltering/OTS/", "path to the zorro ccdb objects"}; + Configurable cfgZorroCCDBpath{"cfgZorroCCDBpath", "EventFiltering/Zorro/", "path to the zorro ccdb objects"}; int mRunNumber = 0; float mBz = 0.f; @@ -356,7 +366,7 @@ struct nucleiSpectra { float dauVtx[3]{0.f, 0.f, 0.f}; auto daughters = particle.daughters_as(); for (const auto& dau : daughters) { - if (abs(dau.pdgCode()) != 22 && abs(dau.pdgCode()) != 11) { + if (std::abs(dau.pdgCode()) != PDG_t::kGamma && std::abs(dau.pdgCode()) != PDG_t::kElectron) { dauVtx[0] = dau.vx(); dauVtx[1] = dau.vy(); dauVtx[2] = dau.vz(); @@ -510,14 +520,14 @@ struct nucleiSpectra { spectra.add("hTpcSignalDataSelected", "Specific energy loss for selected particles", HistType::kTH2F, {{600, -6., 6., "#it{p} (GeV/#it{c})"}, {1400, 0, 1400, "d#it{E} / d#it{X} (a. u.)"}}); spectra.add("hTofSignalData", "TOF beta", HistType::kTH2F, {{500, 0., 5., "#it{p} (GeV/#it{c})"}, {750, 0, 1.5, "TOF #beta"}}); - for (int iC{0}; iC < 2; ++iC) { + for (unsigned int iC{0}; iC < nuclei::matter.size(); ++iC) { nuclei::hGloTOFtracks[iC] = spectra.add(fmt::format("hTPCTOFtracks{}", nuclei::matter[iC]).data(), fmt::format("Global vs TOF matched {} tracks in a collision", nuclei::chargeLabelNames[iC]).data(), HistType::kTH2D, {{300, -0.5, 300.5, "Number of global tracks"}, {300, -0.5, 300.5, "Number of TOF matched tracks"}}); for (int iS{0}; iS < nuclei::species; ++iS) { nuclei::hNsigma[0][iS][iC] = spectra.add(fmt::format("h{}nsigma{}_{}", nuclei::pidName[0], nuclei::matter[iC], nuclei::names[iS]).data(), fmt::format("n#sigma_{{}} {} {}", nuclei::pidName[0], nuclei::matter[iC], nuclei::names[iS]).data(), HistType::kTH3D, {centAxis, ptAxes[iS], nSigmaAxes[0]}); nuclei::hNsigmaEta[0][iS][iC] = spectra.add(fmt::format("h{}nsigmaEta{}_{}", nuclei::pidName[0], nuclei::matter[iC], nuclei::names[iS]).data(), fmt::format("n#sigma_{{}} {} {} vs #eta", nuclei::pidName[0], nuclei::matter[iC], nuclei::names[iS]).data(), HistType::kTH3D, {etaAxis, ptAxes[iS], nSigmaAxes[0]}); - for (int iPID{0}; iPID < 2; ++iPID) { + for (unsigned int iPID{0}; iPID < nuclei::matter.size(); ++iPID) { nuclei::hDCAxy[iPID][iS][iC] = spectra.add(fmt::format("hDCAxy{}_{}_{}", nuclei::pidName[iPID], nuclei::matter[iC], nuclei::names[iS]).data(), fmt::format("DCAxy {} {} {}", nuclei::pidName[iPID], nuclei::matter[iC], nuclei::names[iS]).data(), HistType::kTH3D, {centAxis, ptAxes[iS], dcaxyAxes[iS]}); nuclei::hDCAz[iPID][iS][iC] = spectra.add(fmt::format("hDCAz{}_{}_{}", nuclei::pidName[iPID], nuclei::matter[iC], nuclei::names[iS]).data(), fmt::format("DCAz {} {} {}", nuclei::pidName[iPID], nuclei::matter[iC], nuclei::names[iS]).data(), HistType::kTH3D, {centAxis, ptAxes[iS], dcazAxes[iS]}); } @@ -540,15 +550,15 @@ struct nucleiSpectra { } for (int iS{0}; iS < nuclei::species; ++iS) { - for (int iMax{0}; iMax < 2; ++iMax) { - nuclei::pidCuts[0][iS][iMax] = cfgNsigmaTPC->get(iS, iMax); + for (unsigned int iMax{0}; iMax < nuclei::pidName.size(); ++iMax) { + nuclei::pidCuts[0][iS][iMax] = cfgTrackCut.TPCnSigmaMax->get(iS, iMax); } } if (doprocessMatching) { std::vector occBins{-0.5, 499.5, 999.5, 1999.5, 2999.5, 3999.5, 4999.5, 10000., 50000.}; AxisSpec occAxis{occBins, "Occupancy"}; - for (int iC{0}; iC < 2; ++iC) { + for (unsigned int iC{0}; iC < nuclei::matter.size(); ++iC) { nuclei::hMatchingStudy[iC] = spectra.add(fmt::format("hMatchingStudy{}", nuclei::matter[iC]).data(), ";#it{p}_{T};#phi;#eta;n#sigma_{ITS};n#sigma{TPC};n#sigma_{TOF};Centrality", HistType::kTHnSparseF, {{20, 1., 9.}, {10, 0., o2::constants::math::TwoPI}, {10, -1., 1.}, {50, -5., 5.}, {50, -5., 5.}, {50, 0., 1.}, {8, 0., 80.}}); nuclei::hMatchingStudyHadrons[iC] = spectra.add(fmt::format("hMatchingStudyHadrons{}", nuclei::matter[iC]).data(), ";#it{p}_{T};#phi;#eta;Centrality;Track type; Occupancy", HistType::kTHnF, {{23, 0.4, 5.}, {20, 0., o2::constants::math::TwoPI}, {10, -1., 1.}, {8, 0., 80.}, {2, -0.5, 1.5}, occAxis}); } @@ -605,15 +615,15 @@ struct nucleiSpectra { {nuclei::charges[4] * cfgMomentumScalingBetheBloch->get(3u, 0u) / nuclei::masses[4], nuclei::charges[4] * cfgMomentumScalingBetheBloch->get(3u, 1u) / nuclei::masses[4]}}; int nGloTracks[2]{0, 0}, nTOFTracks[2]{0, 0}; - for (auto& track : tracks) { // start loop over tracks - if (std::abs(track.eta()) > cfgCutEta || - track.tpcInnerParam() < cfgCutTpcMom || - track.itsNCls() < cfgCutNclusITS || - track.tpcNClsFound() < cfgCutNclusTPC || - track.tpcNClsCrossedRows() < 70 || - track.tpcNClsCrossedRows() < 0.8 * track.tpcNClsFindable() || - track.tpcChi2NCl() > 4.f || - track.itsChi2NCl() > 36.f) { + for (const auto& track : tracks) { // start loop over tracks + if (std::abs(track.eta()) > cfgTrackCut.EtaMax || + track.tpcInnerParam() < cfgTrackCut.TPCrigidityMin || + track.itsNCls() < cfgTrackCut.ITSnClusMin || + track.tpcNClsFound() < cfgTrackCut.TPCnClsMin || + track.tpcNClsCrossedRows() < cfgTrackCut.TPCnCrossedRowsMin || + track.tpcNClsCrossedRows() < cfgTrackCut.TPCnCrossedRowsOverFindableMin * track.tpcNClsFindable() || + track.tpcChi2NCl() > cfgTrackCut.TPCchi2ClusMax || + track.itsChi2NCl() > cfgTrackCut.ITSchi2ClusMax) { continue; } // temporary fix: tpcInnerParam() returns the momentum in all the software tags before @@ -678,12 +688,12 @@ struct nucleiSpectra { } for (int iS{0}; iS < nuclei::species; ++iS) { bool selectedTOF{false}; - if (std::abs(dcaInfo[1]) > cfgDCAcut->get(iS, 1)) { + if (std::abs(dcaInfo[1]) > cfgTrackCut.DcaMax->get(iS, 1)) { continue; } ROOT::Math::LorentzVector> fvector{mTrackParCov.getPt() * nuclei::charges[iS], mTrackParCov.getEta(), mTrackParCov.getPhi(), nuclei::masses[iS]}; float y{fvector.Rapidity() + cfgCMrapidity}; - for (int iPID{0}; iPID < 2; ++iPID) { /// 0 TPC, 1 TOF + for (unsigned int iPID{0}; iPID < nuclei::pidName.size(); ++iPID) { /// 0 TPC, 1 TOF if (selectedTPC[iS]) { if (iPID && !track.hasTOF()) { continue; @@ -692,12 +702,12 @@ struct nucleiSpectra { float charge{1.f + static_cast(iS == 3 || iS == 4)}; tofMasses[iS] = correctedTpcInnerParam * charge * std::sqrt(1.f / (beta * beta) - 1.f) - nuclei::masses[iS]; } - if (!cfgCutOnReconstructedRapidity || (y > cfgCutRapidityMin && y < cfgCutRapidityMax)) { + if (!cfgTrackCut.RapidityToggle || (y > cfgTrackCut.RapidityMin && y < cfgTrackCut.RapidityMax)) { if (std::abs(nSigmaTPC[iS]) < cfgNsigmaTPCcutDCAhists && (!iPID || std::abs(tofMasses[iS]) < cfgDeltaTOFmassCutDCAhists)) { nuclei::hDCAxy[iPID][iS][iC]->Fill(centrality, fvector.pt(), dcaInfo[0]); nuclei::hDCAz[iPID][iS][iC]->Fill(centrality, fvector.pt(), dcaInfo[1]); } - if (std::abs(dcaInfo[0]) < cfgDCAcut->get(iS, 0u)) { + if (std::abs(dcaInfo[0]) < cfgTrackCut.DcaMax->get(iS, 0u)) { if (!iPID) { /// temporary exclusion of the TOF nsigma PID for the He3 and Alpha nuclei::hNsigma[iPID][iS][iC]->Fill(centrality, fvector.pt(), nSigma[iPID][iS]); nuclei::hNsigmaEta[iPID][iS][iC]->Fill(fvector.eta(), fvector.pt(), nSigma[iPID][iS]); @@ -857,7 +867,7 @@ struct nucleiSpectra { } } } - for (auto& c : nuclei::candidates_flow) { + for (const auto& c : nuclei::candidates_flow) { nucleiTableFlow(c.centFV0A, c.centFT0M, c.centFT0A, c.centFT0C, c.psiFT0A, c.psiFT0C, c.psiTPC, c.psiTPCl, c.psiTPCr, c.qFT0A, c.qFT0C, c.qTPC, c.qTPCl, c.qTPCr); } } @@ -874,7 +884,7 @@ struct nucleiSpectra { return; } fillDataInfo(collision, tracks); - for (auto& c : nuclei::candidates) { + for (const auto& c : nuclei::candidates) { if (c.fillTree) { nucleiTable(c.pt, c.eta, c.phi, c.tpcInnerParam, c.beta, c.zVertex, c.nContrib, c.DCAxy, c.DCAz, c.TPCsignal, c.ITSchi2, c.TPCchi2, c.TOFchi2, c.flags, c.TPCfindableCls, c.TPCcrossedRows, c.ITSclsMap, c.TPCnCls, c.TPCnClsShared, c.clusterSizesITS); } @@ -886,7 +896,7 @@ struct nucleiSpectra { } } } - for (auto& c : nuclei::candidates_flow) { + for (const auto& c : nuclei::candidates_flow) { nucleiTableFlow(c.centFV0A, c.centFT0M, c.centFT0A, c.centFT0C, c.psiFT0A, c.psiFT0C, c.psiTPC, c.psiTPCl, c.psiTPCr, c.qFT0A, c.qFT0C, c.qTPC, c.qTPCl, c.qTPCr); } } @@ -896,11 +906,11 @@ struct nucleiSpectra { void processMC(soa::Join const& collisions, aod::McCollisions const& mcCollisions, soa::Join const& tracks, aod::McParticles const& particlesMC, aod::BCsWithTimestamps const&) { nuclei::candidates.clear(); - for (auto& c : mcCollisions) { + for (const auto& c : mcCollisions) { spectra.fill(HIST("hGenVtxZ"), c.posZ()); } std::vector goodCollisions(mcCollisions.size(), false); - for (auto& collision : collisions) { + for (const auto& collision : collisions) { if (!eventSelection(collision)) { continue; } @@ -938,7 +948,7 @@ struct nucleiSpectra { if (!storeIt) { continue; } - if (particle.y() < cfgCutRapidityMin || particle.y() > cfgCutRapidityMax) { + if (particle.y() < cfgTrackCut.RapidityMin || particle.y() > cfgTrackCut.RapidityMax) { continue; } @@ -948,7 +958,7 @@ struct nucleiSpectra { if (particle.isPhysicalPrimary()) { c.flags |= kIsPhysicalPrimary; if (particle.has_mothers()) { - for (auto& motherparticle : particle.mothers_as()) { + for (const auto& motherparticle : particle.mothers_as()) { if (std::find(nuclei::hfMothCodes.begin(), nuclei::hfMothCodes.end(), std::abs(motherparticle.pdgCode())) != nuclei::hfMothCodes.end()) { c.flags |= kIsSecondaryFromWeakDecay; motherPdgCode = motherparticle.pdgCode(); @@ -959,7 +969,7 @@ struct nucleiSpectra { } } else if (particle.has_mothers()) { c.flags |= kIsSecondaryFromWeakDecay; - for (auto& motherparticle : particle.mothers_as()) { + for (const auto& motherparticle : particle.mothers_as()) { motherPdgCode = motherparticle.pdgCode(); motherDecRadius = std::hypot(particle.vx() - motherparticle.vx(), particle.vy() - motherparticle.vy()); } @@ -973,13 +983,13 @@ struct nucleiSpectra { } int index{0}; - for (auto& particle : particlesMC) { + for (const auto& particle : particlesMC) { int pdg{std::abs(particle.pdgCode())}; for (int iS{0}; iS < nuclei::species; ++iS) { if (pdg != nuclei::codes[iS]) { continue; } - if (particle.y() < cfgCutRapidityMin || particle.y() > cfgCutRapidityMax) { + if (particle.y() < cfgTrackCut.RapidityMin || particle.y() > cfgTrackCut.RapidityMax) { continue; } @@ -1005,7 +1015,7 @@ struct nucleiSpectra { continue; // skip secondaries from weak decay without mothers } flags |= kIsSecondaryFromWeakDecay; - for (auto& motherparticle : particle.mothers_as()) { + for (const auto& motherparticle : particle.mothers_as()) { motherPdgCode = motherparticle.pdgCode(); motherDecRadius = std::hypot(particle.vx() - motherparticle.vx(), particle.vy() - motherparticle.vy()); } @@ -1038,9 +1048,9 @@ struct nucleiSpectra { const float centrality = collision.centFT0C(); o2::aod::ITSResponse itsResponse; for (const auto& track : tracks) { - if (std::abs(track.eta()) > cfgCutEta || - track.itsNCls() < 7 || - track.itsChi2NCl() > 36.f) { + if (std::abs(track.eta()) > cfgTrackCut.EtaMax || + track.itsNCls() < nuclei::nITSlayers || + track.itsChi2NCl() > cfgTrackCut.ITSchi2ClusMax) { continue; } double expBethe{tpc::BetheBlochAleph(static_cast(track.tpcInnerParam() * 2. / o2::constants::physics::MassHelium3), cfgBetheBlochParams->get(4, 0u), cfgBetheBlochParams->get(4, 1u), cfgBetheBlochParams->get(4, 2u), cfgBetheBlochParams->get(4, 3u), cfgBetheBlochParams->get(4, 4u))}; @@ -1062,7 +1072,7 @@ struct nucleiSpectra { { nuclei::candidates.clear(); std::vector goodCollisions(mcCollisions.size(), false); - for (auto& collision : collisions) { + for (const auto& collision : collisions) { if (!eventSelection(collision)) { continue; } @@ -1085,7 +1095,7 @@ struct nucleiSpectra { if (particle.isPhysicalPrimary()) { c.flags |= kIsPhysicalPrimary; if (particle.has_mothers()) { - for (auto& motherparticle : particle.mothers_as()) { + for (const auto& motherparticle : particle.mothers_as()) { if (std::find(nuclei::hfMothCodes.begin(), nuclei::hfMothCodes.end(), std::abs(motherparticle.pdgCode())) != nuclei::hfMothCodes.end()) { c.flags |= kIsSecondaryFromWeakDecay; motherPdgCode = motherparticle.pdgCode(); @@ -1096,7 +1106,7 @@ struct nucleiSpectra { } } else if (particle.has_mothers()) { c.flags |= kIsSecondaryFromWeakDecay; - for (auto& motherparticle : particle.mothers_as()) { + for (const auto& motherparticle : particle.mothers_as()) { motherPdgCode = motherparticle.pdgCode(); motherDecRadius = std::hypot(particle.vx() - motherparticle.vx(), particle.vy() - motherparticle.vy()); }