diff --git a/PWGCF/FemtoUniverse/Core/FemtoUniverseParticleHisto.h b/PWGCF/FemtoUniverse/Core/FemtoUniverseParticleHisto.h index 508cb8994c4..4b860faae66 100644 --- a/PWGCF/FemtoUniverse/Core/FemtoUniverseParticleHisto.h +++ b/PWGCF/FemtoUniverse/Core/FemtoUniverseParticleHisto.h @@ -70,7 +70,7 @@ class FemtoUniverseParticleHisto // comment template - void init_debug(std::string folderName, T& tempFitVarMomAxis) // o2-linter: disable=name/function-variable + void init_debug(std::string folderName, T& tempFitVarMomAxis, bool isFillITSNsigma) // o2-linter: disable=name/function-variable { std::string folderSuffix = static_cast(o2::aod::femtouniverse_mc_particle::MCTypeName[mc]).c_str(); if constexpr (mParticleType == o2::aod::femtouniverseparticle::ParticleType::kTrack || mParticleType == o2::aod::femtouniverseparticle::ParticleType::kV0Child || mParticleType == o2::aod::femtouniverseparticle::ParticleType::kCascadeBachelor || mParticleType == o2::aod::femtouniverseparticle::ParticleType::kMCTruthTrack) { @@ -102,6 +102,18 @@ class FemtoUniverseParticleHisto mHistogramRegistry->add((folderName + folderSuffix + "/nSigmaComb_K").c_str(), "; #it{p} (GeV/#it{c}); n#sigma_{comb}^{K}", kTH2F, {{tempFitVarMomAxis}, {100, 0, 5}}); mHistogramRegistry->add((folderName + folderSuffix + "/nSigmaComb_p").c_str(), "; #it{p} (GeV/#it{c}); n#sigma_{comb}^{p}", kTH2F, {{tempFitVarMomAxis}, {100, 0, 5}}); mHistogramRegistry->add((folderName + folderSuffix + "/nSigmaComb_d").c_str(), "; #it{p} (GeV/#it{c}); n#sigma_{comb}^{d}", kTH2F, {{tempFitVarMomAxis}, {100, 0, 5}}); + if (isFillITSNsigma) { + mHistogramRegistry->add((folderName + folderSuffix + "/nSigmaITS_el").c_str(), "; #it{p} (GeV/#it{c}); n#sigma_{ITS}^{e}", kTH2F, {{tempFitVarMomAxis}, {100, 0, 5}}); + mHistogramRegistry->add((folderName + folderSuffix + "/nSigmaITS_pi").c_str(), "; #it{p} (GeV/#it{c}); n#sigma_{ITS}^{#pi}", kTH2F, {{tempFitVarMomAxis}, {100, 0, 5}}); + mHistogramRegistry->add((folderName + folderSuffix + "/nSigmaITS_K").c_str(), "; #it{p} (GeV/#it{c}); n#sigma_{ITS}^{K}", kTH2F, {{tempFitVarMomAxis}, {100, 0, 5}}); + mHistogramRegistry->add((folderName + folderSuffix + "/nSigmaITS_p").c_str(), "; #it{p} (GeV/#it{c}); n#sigma_{ITS}^{p}", kTH2F, {{tempFitVarMomAxis}, {100, 0, 5}}); + mHistogramRegistry->add((folderName + folderSuffix + "/nSigmaITS_d").c_str(), "; #it{p} (GeV/#it{c}); n#sigma_{ITS}^{d}", kTH2F, {{tempFitVarMomAxis}, {100, 0, 5}}); + mHistogramRegistry->add((folderName + folderSuffix + "/nSigmaCombITSTPC_el").c_str(), "; #it{p} (GeV/#it{c}); n#sigma_{ITSTPC}^{e}", kTH2F, {{tempFitVarMomAxis}, {100, 0, 5}}); + mHistogramRegistry->add((folderName + folderSuffix + "/nSigmaCombITSTPC_pi").c_str(), "; #it{p} (GeV/#it{c}); n#sigma_{ITSTPC}^{#pi}", kTH2F, {{tempFitVarMomAxis}, {100, 0, 5}}); + mHistogramRegistry->add((folderName + folderSuffix + "/nSigmaCombITSTPC_K").c_str(), "; #it{p} (GeV/#it{c}); n#sigma_{ITSTPC}^{K}", kTH2F, {{tempFitVarMomAxis}, {100, 0, 5}}); + mHistogramRegistry->add((folderName + folderSuffix + "/nSigmaCombITSTPC_p").c_str(), "; #it{p} (GeV/#it{c}); n#sigma_{ITSTPC}^{p}", kTH2F, {{tempFitVarMomAxis}, {100, 0, 5}}); + mHistogramRegistry->add((folderName + folderSuffix + "/nSigmaCombITSTPC_d").c_str(), "; #it{p} (GeV/#it{c}); n#sigma_{ITSTPC}^{d}", kTH2F, {{tempFitVarMomAxis}, {100, 0, 5}}); + } } else if constexpr (mParticleType == o2::aod::femtouniverseparticle::ParticleType::kV0) { mHistogramRegistry->add((folderName + folderSuffix + "/hDaughDCA").c_str(), "; DCA^{daugh} (cm); Entries", kTH1F, {{1000, 0, 10}}); mHistogramRegistry->add((folderName + folderSuffix + "/hTransRadius").c_str(), "; #it{r}_{xy} (cm); Entries", kTH1F, {{1500, 0, 150}}); @@ -208,7 +220,7 @@ class FemtoUniverseParticleHisto /// \param tempFitVarBins binning of the tempFitVar (DCA_xy in case of tracks, CPA in case of V0s, etc.) /// \param isMC add Monte Carlo truth histograms to the output file template - void init(HistogramRegistry* registry, T& tempFitVarpTBins, T& tempFitVarBins, bool isMC, int pdgCode, bool isDebug = false, std::optional flexibleFolder = std::nullopt) + void init(HistogramRegistry* registry, T& tempFitVarpTBins, T& tempFitVarBins, bool isMC, int pdgCode, bool isDebug = false, std::optional flexibleFolder = std::nullopt, bool isFillITSNsigma = false) { mPDG = pdgCode; if (registry) { @@ -247,7 +259,7 @@ class FemtoUniverseParticleHisto // Fill here the actual histogramms by calling init_base and init_MC init_base(folderName, tempFitVarAxisTitle, tempFitVarpTAxis, tempFitVarAxis); if (isDebug) { - init_debug(folderName, tempFitVarMomAxis); + init_debug(folderName, tempFitVarMomAxis, isFillITSNsigma); } if (isMC) { init_base(folderName, tempFitVarAxisTitle, tempFitVarpTAxis, tempFitVarAxis); @@ -512,6 +524,28 @@ class FemtoUniverseParticleHisto } } + template + void fillQAITSPID(T const& part) + { + fillQABaseITSPID(part, HIST(o2::aod::femtouniverseparticle::ParticleTypeName[mParticleType]) + HIST(mFolderSuffix[mFolderSuffixType])); + } + + template + void fillQABaseITSPID(T const& part, H const& histFolder) + { + // std::string tempFitVarName; + mHistogramRegistry->fill(histFolder + HIST(o2::aod::femtouniverse_mc_particle::MCTypeName[o2::aod::femtouniverse_mc_particle::MCType::kRecon]) + HIST("/nSigmaITS_el"), part.p(), part.itsNSigmaEl()); + mHistogramRegistry->fill(histFolder + HIST(o2::aod::femtouniverse_mc_particle::MCTypeName[o2::aod::femtouniverse_mc_particle::MCType::kRecon]) + HIST("/nSigmaITS_pi"), part.p(), part.itsNSigmaPi()); + mHistogramRegistry->fill(histFolder + HIST(o2::aod::femtouniverse_mc_particle::MCTypeName[o2::aod::femtouniverse_mc_particle::MCType::kRecon]) + HIST("/nSigmaITS_K"), part.p(), part.itsNSigmaKa()); + mHistogramRegistry->fill(histFolder + HIST(o2::aod::femtouniverse_mc_particle::MCTypeName[o2::aod::femtouniverse_mc_particle::MCType::kRecon]) + HIST("/nSigmaITS_p"), part.p(), part.itsNSigmaPr()); + mHistogramRegistry->fill(histFolder + HIST(o2::aod::femtouniverse_mc_particle::MCTypeName[o2::aod::femtouniverse_mc_particle::MCType::kRecon]) + HIST("/nSigmaITS_d"), part.p(), part.itsNSigmaDe()); + mHistogramRegistry->fill(histFolder + HIST(o2::aod::femtouniverse_mc_particle::MCTypeName[o2::aod::femtouniverse_mc_particle::MCType::kRecon]) + HIST("/nSigmaCombITSTPC_el"), part.p(), std::sqrt(part.tpcNSigmaEl() * part.tpcNSigmaEl() + part.itsNSigmaEl() * part.itsNSigmaEl())); + mHistogramRegistry->fill(histFolder + HIST(o2::aod::femtouniverse_mc_particle::MCTypeName[o2::aod::femtouniverse_mc_particle::MCType::kRecon]) + HIST("/nSigmaCombITSTPC_pi"), part.p(), std::sqrt(part.tpcNSigmaPi() * part.tpcNSigmaPi() + part.itsNSigmaPi() * part.itsNSigmaPi())); + mHistogramRegistry->fill(histFolder + HIST(o2::aod::femtouniverse_mc_particle::MCTypeName[o2::aod::femtouniverse_mc_particle::MCType::kRecon]) + HIST("/nSigmaCombITSTPC_K"), part.p(), std::sqrt(part.tpcNSigmaKa() * part.tpcNSigmaKa() + part.itsNSigmaKa() * part.itsNSigmaKa())); + mHistogramRegistry->fill(histFolder + HIST(o2::aod::femtouniverse_mc_particle::MCTypeName[o2::aod::femtouniverse_mc_particle::MCType::kRecon]) + HIST("/nSigmaCombITSTPC_p"), part.p(), std::sqrt(part.tpcNSigmaPr() * part.tpcNSigmaPr() + part.itsNSigmaPr() * part.itsNSigmaPr())); + mHistogramRegistry->fill(histFolder + HIST(o2::aod::femtouniverse_mc_particle::MCTypeName[o2::aod::femtouniverse_mc_particle::MCType::kRecon]) + HIST("/nSigmaCombITSTPC_d"), part.p(), std::sqrt(part.tpcNSigmaDe() * part.tpcNSigmaDe() + part.itsNSigmaDe() * part.itsNSigmaDe())); + } + /// Templated function to fill particle histograms for data/ Monte Carlo reconstructed and Monte Carlo truth /// Always calls fillQA_base fill histogramms with data/ Monte Carlo reconstructed /// In case of Monte Carlo, calls fillQA_base with Monte Carlo truth info and specialized function fillQA_MC for additional histogramms diff --git a/PWGCF/FemtoUniverse/Core/FemtoUniverseTrackSelection.h b/PWGCF/FemtoUniverse/Core/FemtoUniverseTrackSelection.h index 3aee759b29b..a7145af2c71 100644 --- a/PWGCF/FemtoUniverse/Core/FemtoUniverseTrackSelection.h +++ b/PWGCF/FemtoUniverse/Core/FemtoUniverseTrackSelection.h @@ -23,6 +23,7 @@ #include "Common/Core/TrackSelection.h" #include "Common/Core/TrackSelectionDefaults.h" +#include "Common/DataModel/PIDResponseITS.h" #include "Common/DataModel/PIDResponseTOF.h" #include "Common/DataModel/PIDResponseTPC.h" #include "Common/DataModel/TrackSelectionTables.h" @@ -135,6 +136,14 @@ class FemtoUniverseTrackSelection : public FemtoUniverseObjectSelection auto getNsigmaTOF(T const& track, o2::track::PID pid); + /// Computes the n_sigma for a track and a particle-type hypothesis in the TOF + /// \tparam T Data type of the track + /// \param track Track for which PID is evaluated + /// \param pid Particle species for which PID is evaluated + /// \return Value of n_{sigma, TOF} + template + auto getNsigmaITS(T const& track, o2::track::PID pid); + /// Checks whether the most open combination of all selection criteria is fulfilled /// \tparam T Data type of the track /// \param track Track @@ -151,6 +160,18 @@ class FemtoUniverseTrackSelection : public FemtoUniverseObjectSelection std::array getCutContainer(T const& track); + /// Obtain the bit-wise container for the selections + /// \todo For the moment, PID is separated from the other selections, hence + /// instead of a single value an std::array of size two is returned + /// \tparam CutContainerType Data type of the bit-wise container for the + /// selections + /// \tparam T Data type of the track + /// \param track Track + /// \return The bit-wise container for the selections, separately with all + /// selection criteria, and the PID + template + std::array getCutContainerWithITS(T const& track); + /// Some basic QA histograms /// \tparam part Type of the particle for proper naming of the folders for QA /// \tparam tracktype Type of track (track, positive child, negative child) for proper naming of the folders for QA @@ -329,6 +350,11 @@ void FemtoUniverseTrackSelection::init(HistogramRegistry* registry) mHistogramRegistry->add((folderName + "/hDCAz").c_str(), "; #it{p}_{T} (GeV/#it{c}); DCA_{z} (cm)", kTH2F, {{100, 0, 10}, {500, -5, 5}}); mHistogramRegistry->add((folderName + "/hDCA").c_str(), "; #it{p}_{T} (GeV/#it{c}); DCA (cm)", kTH2F, {{100, 0, 10}, {301, 0., 1.5}}); mHistogramRegistry->add((folderName + "/hTPCdEdX").c_str(), "; #it{p} (GeV/#it{c}); TPC Signal", kTH2F, {{100, 0, 10}, {1000, 0, 1000}}); + mHistogramRegistry->add((folderName + "/nSigmaITS_el").c_str(), "; #it{p} (GeV/#it{c}); n#sigma_{ITS}^{e}", kTH2F, {{100, 0, 10}, {200, -4.975, 5.025}}); + mHistogramRegistry->add((folderName + "/nSigmaITS_pi").c_str(), "; #it{p} (GeV/#it{c}); n#sigma_{ITS}^{#pi}", kTH2F, {{100, 0, 10}, {200, -4.975, 5.025}}); + mHistogramRegistry->add((folderName + "/nSigmaITS_K").c_str(), "; #it{p} (GeV/#it{c}); n#sigma_{ITS}^{K}", kTH2F, {{100, 0, 10}, {200, -4.975, 5.025}}); + mHistogramRegistry->add((folderName + "/nSigmaITS_p").c_str(), "; #it{p} (GeV/#it{c}); n#sigma_{ITS}^{p}", kTH2F, {{100, 0, 10}, {200, -4.975, 5.025}}); + mHistogramRegistry->add((folderName + "/nSigmaITS_d").c_str(), "; #it{p} (GeV/#it{c}); n#sigma_{ITS}^{d}", kTH2F, {{100, 0, 10}, {200, -4.975, 5.025}}); mHistogramRegistry->add((folderName + "/nSigmaTPC_el").c_str(), "; #it{p} (GeV/#it{c}); n#sigma_{TPC}^{e}", kTH2F, {{100, 0, 10}, {200, -4.975, 5.025}}); mHistogramRegistry->add((folderName + "/nSigmaTPC_pi").c_str(), "; #it{p} (GeV/#it{c}); n#sigma_{TPC}^{#pi}", kTH2F, {{100, 0, 10}, {200, -4.975, 5.025}}); mHistogramRegistry->add((folderName + "/nSigmaTPC_K").c_str(), "; #it{p} (GeV/#it{c}); n#sigma_{TPC}^{K}", kTH2F, {{100, 0, 10}, {200, -4.975, 5.025}}); @@ -394,6 +420,24 @@ auto FemtoUniverseTrackSelection::getNsigmaTOF(T const& track, o2::track::PID pi return o2::aod::pidutils::tofNSigma(pid, track); } +template +auto FemtoUniverseTrackSelection::getNsigmaITS(T const& track, o2::track::PID pid) +{ + if (pid == o2::track::PID::Electron) { + return track.itsNSigmaEl(); + } else if (pid == o2::track::PID::Pion) { + return track.itsNSigmaPi(); + } else if (pid == o2::track::PID::Kaon) { + return track.itsNSigmaKa(); + } else if (pid == o2::track::PID::Proton) { + return track.itsNSigmaPr(); + } else if (pid == o2::track::PID::Deuteron) { + return track.itsNSigmaDe(); + } + // if nothing matched, return default value + return -999.f; +} + template bool FemtoUniverseTrackSelection::isSelectedMinimal(T const& track) { @@ -565,6 +609,101 @@ std::array FemtoUniverseTrackSelection::getCutContainer(T c return {output, outputPID}; } +template +std::array + FemtoUniverseTrackSelection::getCutContainerWithITS(T const& track) +{ + CutContainerType output = 0; + size_t counter = 0; + CutContainerType outputPID = 0; + const auto sign = track.sign(); + const auto pT = track.pt(); + const auto eta = track.eta(); + const auto tpcNClsF = track.tpcNClsFound(); + const auto tpcRClsC = track.tpcCrossedRowsOverFindableCls(); + const auto tpcNClsC = track.tpcNClsCrossedRows(); + const auto tpcNClsS = track.tpcNClsShared(); + const auto tpcNClsFracS = track.tpcFractionSharedCls(); + const auto itsNCls = track.itsNCls(); + const auto itsNClsIB = track.itsNClsInnerBarrel(); + const auto dcaXY = track.dcaXY(); + const auto dcaZ = track.dcaZ(); + const auto dca = std::sqrt(std::pow(dcaXY, 2.) + std::pow(dcaZ, 2.)); + + std::vector pidTPC, pidTOF, pidITS; + for (auto it : kPIDspecies) { + pidTPC.push_back(getNsigmaTPC(track, it)); + pidTOF.push_back(getNsigmaTOF(track, it)); + pidITS.push_back(getNsigmaITS(track, it)); + } + + float observable = 0.; + for (auto& sel : mSelections) { + const auto selVariable = sel.getSelectionVariable(); + if (selVariable == femto_universe_track_selection::kPIDnSigmaMax) { + /// PID needs to be handled a bit differently since we may need more than + /// one species + for (size_t i = 0; i < kPIDspecies.size(); ++i) { + auto pidITSVal = pidITS.at(i); + auto pidTPCVal = pidTPC.at(i) - nSigmaPIDOffsetTPC; + auto pidTOFVal = pidTOF.at(i) - nSigmaPIDOffsetTOF; + auto pidComb = std::sqrt(pidTPCVal * pidTPCVal + pidTOFVal * pidTOFVal); + sel.checkSelectionSetBitPID(pidTPCVal, outputPID); + sel.checkSelectionSetBitPID(pidComb, outputPID); + sel.checkSelectionSetBitPID(pidITSVal, outputPID); + } + } else { + /// for the rest it's all the same + switch (selVariable) { + case (femto_universe_track_selection::kSign): + observable = sign; + break; + case (femto_universe_track_selection::kpTMin): + case (femto_universe_track_selection::kpTMax): + observable = pT; + break; + case (femto_universe_track_selection::kEtaMax): + observable = eta; + break; + case (femto_universe_track_selection::kTPCnClsMin): + observable = tpcNClsF; + break; + case (femto_universe_track_selection::kTPCfClsMin): + observable = tpcRClsC; + break; + case (femto_universe_track_selection::kTPCcRowsMin): + observable = tpcNClsC; + break; + case (femto_universe_track_selection::kTPCsClsMax): + observable = tpcNClsS; + break; + case (femto_universe_track_selection::kTPCfracsClsMax): + observable = tpcNClsFracS; + break; + case (femto_universe_track_selection::kITSnClsMin): + observable = itsNCls; + break; + case (femto_universe_track_selection::kITSnClsIbMin): + observable = itsNClsIB; + break; + case (femto_universe_track_selection::kDCAxyMax): + observable = dcaXY; + break; + case (femto_universe_track_selection::kDCAzMax): + observable = dcaZ; + break; + case (femto_universe_track_selection::kDCAMin): + observable = dca; + break; + case (femto_universe_track_selection::kPIDnSigmaMax): + break; + } + sel.checkSelectionSetBit(observable, output, counter); + } + } + return {output, outputPID}; +} + template void FemtoUniverseTrackSelection::fillQA(T const& track) { diff --git a/PWGCF/FemtoUniverse/DataModel/FemtoDerived.h b/PWGCF/FemtoUniverse/DataModel/FemtoDerived.h index c9b158ce1b5..8c5614b12ad 100644 --- a/PWGCF/FemtoUniverse/DataModel/FemtoDerived.h +++ b/PWGCF/FemtoUniverse/DataModel/FemtoDerived.h @@ -123,6 +123,15 @@ DECLARE_SOA_DYNAMIC_COLUMN(P, p, //! Compute the overall momentum in GeV/c [](float pt, float eta) -> float { return pt * std::cosh(eta); }); + +DECLARE_SOA_COLUMN(ITSNSigmaEl, itsNSigmaEl, float); //! Nsigma separation with the Its detector for electron +DECLARE_SOA_COLUMN(ITSNSigmaPi, itsNSigmaPi, float); //! Nsigma separation with the Its detector for pion +DECLARE_SOA_COLUMN(ITSNSigmaKa, itsNSigmaKa, float); //! Nsigma separation with the Its detector for kaon +DECLARE_SOA_COLUMN(ITSNSigmaPr, itsNSigmaPr, float); //! Nsigma separation with the Its detector for proton +DECLARE_SOA_COLUMN(ITSNSigmaDe, itsNSigmaDe, float); //! Nsigma separation with the Its detector for deuteron +DECLARE_SOA_COLUMN(ITSNSigmaTr, itsNSigmaTr, float); //! Nsigma separation with the Its detector for triton +DECLARE_SOA_COLUMN(ITSNSigmaHe, itsNSigmaHe, float); //! Nsigma separation with the Its detector for helium3 + // debug variables DECLARE_SOA_COLUMN(Sign, sign, int8_t); //! Sign of the track charge DECLARE_SOA_COLUMN(TpcNClsFound, tpcNClsFound, uint8_t); //! Number of TPC clusters @@ -222,6 +231,14 @@ DECLARE_SOA_TABLE(FDExtParticles, "AOD", "FDEXTPARTICLE", pidtof_tiny::TOFNSigmaDe); using FDFullParticle = FDExtParticles::iterator; +DECLARE_SOA_TABLE(FDItsParticles, "AOD", "FDITSPARTICLE", + femtouniverseparticle::ITSNSigmaEl, + femtouniverseparticle::ITSNSigmaPi, + femtouniverseparticle::ITSNSigmaKa, + femtouniverseparticle::ITSNSigmaPr, + femtouniverseparticle::ITSNSigmaDe); +using FDItsParticle = FDItsParticles::iterator; + DECLARE_SOA_TABLE(FDCascParticles, "AOD", "FDCASCPARTICLE", o2::soa::Index<>, femtouniverseparticle::FdCollisionId, diff --git a/PWGCF/FemtoUniverse/TableProducer/femtoUniverseProducerTask.cxx b/PWGCF/FemtoUniverse/TableProducer/femtoUniverseProducerTask.cxx index d73fc6794e5..a49e9c0cdf7 100644 --- a/PWGCF/FemtoUniverse/TableProducer/femtoUniverseProducerTask.cxx +++ b/PWGCF/FemtoUniverse/TableProducer/femtoUniverseProducerTask.cxx @@ -39,6 +39,7 @@ #include "Common/DataModel/Centrality.h" #include "Common/DataModel/EventSelection.h" #include "Common/DataModel/Multiplicity.h" +#include "Common/DataModel/PIDResponseITS.h" #include "Common/DataModel/PIDResponseTOF.h" #include "Common/DataModel/PIDResponseTPC.h" #include "Common/DataModel/TrackSelectionTables.h" @@ -117,11 +118,13 @@ struct FemtoUniverseProducerTask { Produces outputParts; Produces outputPartsMC; Produces outputDebugParts; + Produces outputDebugITSParts; Produces outputPartsMCLabels; Produces outputDebugPartsMC; Produces outputCascParts; Configurable confIsDebug{"confIsDebug", true, "Enable Debug tables"}; + Configurable confFillITSPid{"confFillITSPid", false, "Fill ITSPid information"}; Configurable confIsUseCutculator{"confIsUseCutculator", true, "Enable cutculator for track cuts"}; // Choose if filtering or skimming version is run // Configurable confIsTrigger{"confIsTrigger", false, "Store all collisions"}; //Commented: not used configurable @@ -201,6 +204,12 @@ struct FemtoUniverseProducerTask { Configurable> confTrkPIDspecies{"confTrkPIDspecies", std::vector{o2::track::PID::Pion, o2::track::PID::Kaon, o2::track::PID::Proton, o2::track::PID::Deuteron}, "Trk sel: Particles species for PID (Pion=2, Kaon=3, Proton=4, Deuteron=5)"}; Configurable confIsOnlyMCTrack{"confIsOnlyMCTrack", false, "Enable filling of only MC Tracks"}; // Numbers from ~/alice/O2/DataFormats/Reconstruction/include/ReconstructionDataFormats/PID.h //static constexpr ID Pion = 2; static constexpr ID Kaon = 3; static constexpr ID Proton = 4; static constexpr ID Deuteron = 5; + + Configurable confTrkMinChi2PerClusterTPC{"confTrkMinChi2PerClusterTPC", 0.0, "Lower limit of the Chi2 per TPC Cluster"}; + Configurable confTrkMaxChi2PerClusterTPC{"confTrkMaxChi2PerClusterTPC", 100.0, "Upper limit of the Chi2 per TPC Cluster"}; + Configurable confTrkMaxChi2PerClusterITS{"confTrkMaxChi2PerClusterITS", 100.0, "Upper limit of the Chi2 per ITS Cluster"}; + Configurable confTrkTPCRefit{"confTrkTPCRefit", false, "Enable TPC refit"}; + Configurable confTrkITSRefit{"confTrkITSRefit", false, "Enable ITS refit"}; } ConfTrkSelection; Configurable confTrkPIDnSigmaOffsetTPC{"confTrkPIDnSigmaOffsetTPC", 0., "Offset for TPC nSigma because of bad calibration"}; @@ -1220,6 +1229,13 @@ struct FemtoUniverseProducerTask { continue; } + if (track.tpcChi2NCl() < ConfTrkSelection.confTrkMinChi2PerClusterTPC || track.tpcChi2NCl() > ConfTrkSelection.confTrkMaxChi2PerClusterTPC) { + continue; + } + if (ConfTrkSelection.confTrkTPCRefit && !track.hasTPC()) { + continue; + } + if (track.pt() > confTOFpTmin) { if (!track.hasTOF()) { continue; @@ -1254,6 +1270,61 @@ struct FemtoUniverseProducerTask { } } + template + void fillTracksITSTPCTOF(TrackType const& tracks) + { + std::vector childIDs = {0, 0}; // these IDs are necessary to keep track of the children + std::vector tmpIDtrack; // this vector keeps track of the matching of the primary track table row <-> aod::track table global index + + for (const auto& track : tracks) { + /// if the most open selection criteria are not fulfilled there is no + /// point looking further at the track + + if (!trackCuts.isSelectedMinimal(track)) { + continue; + } + + if (track.tpcChi2NCl() < ConfTrkSelection.confTrkMinChi2PerClusterTPC || track.tpcChi2NCl() > ConfTrkSelection.confTrkMaxChi2PerClusterTPC) { + continue; + } + if (track.itsChi2NCl() > ConfTrkSelection.confTrkMaxChi2PerClusterITS) { + continue; + } + if ((ConfTrkSelection.confTrkTPCRefit && !track.hasTPC()) || (ConfTrkSelection.confTrkITSRefit && !track.hasITS())) { + continue; + } + + if (track.pt() > confTOFpTmin) { + if (!track.hasTOF()) { + continue; + } + } + + trackCuts.fillQA(track); + // the bit-wise container of the systematic variations is obtained + auto cutContainerITS = trackCuts.getCutContainerWithITS(track); + + // now the table is filled + outputParts(outputCollision.lastIndex(), track.pt(), track.eta(), track.phi(), + aod::femtouniverseparticle::ParticleType::kTrack, + cutContainerITS.at( + femto_universe_track_selection::TrackContainerPosition::kCuts), + cutContainerITS.at( + femto_universe_track_selection::TrackContainerPosition::kPID), + track.dcaXY(), childIDs, 0, + track.sign()); // sign getter is mAntiLambda() + + tmpIDtrack.push_back(track.globalIndex()); + fillDebugParticle(track); + outputDebugITSParts(track.itsNSigmaEl(), track.itsNSigmaPi(), + track.itsNSigmaKa(), track.itsNSigmaPr(), + track.itsNSigmaDe()); + if constexpr (isMC) { + fillMCParticle(track, o2::aod::femtouniverseparticle::ParticleType::kTrack); + } + } + } + template void fillV0(CollisionType const& col, V0Type const& fullV0s, TrackType const&) { @@ -2838,11 +2909,17 @@ struct FemtoUniverseProducerTask { getMagneticFieldTesla(bc); const auto ir = mRateFetcher.fetch(ccdb.service, bc.timestamp(), mRunNumber, "ZNC hadronic") * 1.e-3; // fetch IR + auto tracksWithItsPid = soa::Attach, aod::pidits::ITSNSigmaEl, aod::pidits::ITSNSigmaPi, aod::pidits::ITSNSigmaKa, aod::pidits::ITSNSigmaPr, aod::pidits::ITSNSigmaDe, aod::pidits::ITSNSigmaTr, aod::pidits::ITSNSigmaHe>(tracks); + // fill the tables const auto colcheck = fillCollisionsCentRun3(col); if (colcheck) { fillCollisionsCentRun3ColExtra(col, ir); - fillTracks(tracks); + if (!confFillITSPid) { + fillTracks(tracks); + } else { + fillTracksITSTPCTOF(tracksWithItsPid); + } } } PROCESS_SWITCH(FemtoUniverseProducerTask, processTrackCentRun3Data, "Provide experimental data for Run 3 with centrality for track track", false); diff --git a/PWGCF/FemtoUniverse/Tasks/femtoUniversePairTaskTrackTrackSpherHarMultKtExtended.cxx b/PWGCF/FemtoUniverse/Tasks/femtoUniversePairTaskTrackTrackSpherHarMultKtExtended.cxx index a1891f20ecf..f03ad811289 100644 --- a/PWGCF/FemtoUniverse/Tasks/femtoUniversePairTaskTrackTrackSpherHarMultKtExtended.cxx +++ b/PWGCF/FemtoUniverse/Tasks/femtoUniversePairTaskTrackTrackSpherHarMultKtExtended.cxx @@ -72,7 +72,9 @@ struct femtoUniversePairTaskTrackTrackSpherHarMultKtExtended { struct : o2::framework::ConfigurableGroup { Configurable confNsigmaCombined{"confNsigmaCombined", 3.0f, "TPC and TOF Pion Sigma (combined) for momentum > confTOFPtMin"}; Configurable confNsigmaTPC{"confNsigmaTPC", 3.0f, "TPC Pion Sigma for momentum < confTOFPtMin"}; + Configurable confNsigmaITSTPCCombined{"confNsigmaITSTPCCombined", 3.0f, "ITS and TPC Pion Sigma (combined) for momentum < confTOFPtMin"}; Configurable confIsElReject{"confIsElReject", false, "Is electron rejection activated"}; + Configurable confIsAddITSNsigma{"confIsAddITSNsigma", false, "Is ITS Pions nsigma added"}; Configurable confNsigmaTPCElRejectMin{"confNsigmaTPCElRejectMin", 2.0f, "TPC Electron SigmaMin for momentum < confTOFPtMin"}; Configurable confNsigmaTPCElRejectMax{"confNsigmaTPCElRejectMax", 2.0f, "TPC Electron SigmaMax for momentum < confTOFPtMin"}; Configurable confTOFPtMin{"confTOFPtMin", 0.5f, "Min. Pt for which TOF is required for PID."}; @@ -133,7 +135,7 @@ struct femtoUniversePairTaskTrackTrackSpherHarMultKtExtended { ConfigurableAxis confDeltaPhiStarAxis{"confDeltaPhiStarAxis", {100, -0.15, 0.15}, "DeltaPhiStar"}; } twotracksconfigs; - using FemtoFullParticles = soa::Join; + using FemtoFullParticles = soa::Join; // Filters for selecting particles (both p1 and p2) Filter trackAdditionalfilter = ((nabs(aod::femtouniverseparticle::eta) < twotracksconfigs.confEtaMax) && (aod::track::dcaXY <= twotracksconfigs.confTrkDCAxyMax) && @@ -145,13 +147,13 @@ struct femtoUniversePairTaskTrackTrackSpherHarMultKtExtended { using FilteredFemtoFullParticles = soa::Filtered; // using FilteredFemtoFullParticles = FemtoFullParticles; //if no filtering is applied uncomment this optionconfIsCPRkT - using FemtoRecoParticles = soa::Join; + using FemtoRecoParticles = soa::Join; using FilteredFemtoRecoParticles = soa::Filtered; SliceCache cache; Preslice perCol = aod::femtouniverseparticle::fdCollisionId; - using FemtoTruthParticles = soa::Join; + using FemtoTruthParticles = soa::Join; Preslice perColMCTruth = aod::femtouniverseparticle::fdCollisionId; /// Particle 1 @@ -361,7 +363,7 @@ struct femtoUniversePairTaskTrackTrackSpherHarMultKtExtended { } } - bool IsPionNSigma(float mom, float nsigmaTPCPi, float nsigmaTOFPi, float nsigmaTPCElReject) + bool IsPionNSigma(float mom, float nsigmaITSPi, float nsigmaTPCPi, float nsigmaTOFPi, float nsigmaTPCElReject) { //|nsigma_TPC| < 3 for p < 0.5 GeV/c //|nsigma_combined| < 3 for p > 0.5 @@ -378,6 +380,12 @@ struct femtoUniversePairTaskTrackTrackSpherHarMultKtExtended { } else { return false; } + } else if (twotracksconfigs.confIsAddITSNsigma) { + if (std::hypot(nsigmaITSPi, nsigmaTPCPi) < twotracksconfigs.confNsigmaITSTPCCombined) { + return true; + } else { + return false; + } } else { if ((std::abs(nsigmaTPCPi) < twotracksconfigs.confNsigmaTPC)) { return true; @@ -396,7 +404,7 @@ struct femtoUniversePairTaskTrackTrackSpherHarMultKtExtended { return false; } - bool IsParticleNSigma(int8_t particle_number, float mom, float nsigmaTPCPr, float nsigmaTOFPr, float nsigmaTPCPi, float nsigmaTOFPi, float nsigmaTPCK, float nsigmaTOFK, float nsigmaTPCElReject) + bool IsParticleNSigma(int8_t particle_number, float mom, float nsigmaITSPi, float nsigmaTPCPr, float nsigmaTOFPr, float nsigmaTPCPi, float nsigmaTOFPi, float nsigmaTPCK, float nsigmaTOFK, float nsigmaTPCElReject) { if (particle_number == 1) { switch (trackonefilter.confPDGCodePartOne) { @@ -406,7 +414,7 @@ struct femtoUniversePairTaskTrackTrackSpherHarMultKtExtended { break; case 211: // Pion+ case -211: // Pion- - return IsPionNSigma(mom, nsigmaTPCPi, nsigmaTOFPi, nsigmaTPCElReject); + return IsPionNSigma(mom, nsigmaITSPi, nsigmaTPCPi, nsigmaTOFPi, nsigmaTPCElReject); break; case 321: // Kaon+ case -321: // Kaon- @@ -424,7 +432,7 @@ struct femtoUniversePairTaskTrackTrackSpherHarMultKtExtended { break; case 211: // Pion+ case -211: // Pion- - return IsPionNSigma(mom, nsigmaTPCPi, nsigmaTOFPi, nsigmaTPCElReject); + return IsPionNSigma(mom, nsigmaITSPi, nsigmaTPCPi, nsigmaTOFPi, nsigmaTPCElReject); break; case 321: // Kaon+ case -321: // Kaon- @@ -443,9 +451,9 @@ struct femtoUniversePairTaskTrackTrackSpherHarMultKtExtended { void init(InitContext&) { eventHisto.init(&qaRegistry); - trackHistoPartOne.init(&qaRegistry, confTempFitVarpTBins, confTempFitVarBins, twotracksconfigs.confIsMC, trackonefilter.confPDGCodePartOne, true); + trackHistoPartOne.init(&qaRegistry, confTempFitVarpTBins, confTempFitVarBins, twotracksconfigs.confIsMC, trackonefilter.confPDGCodePartOne, true, std::nullopt, twotracksconfigs.confIsAddITSNsigma); - trackHistoPartTwo.init(&qaRegistry, confTempFitVarpTBins, confTempFitVarBins, twotracksconfigs.confIsMC, tracktwofilter.confPDGCodePartTwo, true); + trackHistoPartTwo.init(&qaRegistry, confTempFitVarpTBins, confTempFitVarBins, twotracksconfigs.confIsMC, tracktwofilter.confPDGCodePartTwo, true, std::nullopt, twotracksconfigs.confIsAddITSNsigma); MixQaRegistry.add("MixingQA/hSECollisionBins", ";bin;Entries", kTH1F, {{120, -0.5, 119.5}}); MixQaRegistry.add("MixingQA/hMECollisionBins", ";bin;Entries", kTH1F, {{120, -0.5, 119.5}}); @@ -535,20 +543,26 @@ struct femtoUniversePairTaskTrackTrackSpherHarMultKtExtended { /// Histogramming same event if ((contType == PairType::PlusMinus || contType == PairType::PlusPlus) && fillQA) { for (const auto& part : groupPartsOne) { - if (!IsParticleNSigma((int8_t)1, part.p(), trackCuts.getNsigmaTPC(part, o2::track::PID::Proton), trackCuts.getNsigmaTOF(part, o2::track::PID::Proton), trackCuts.getNsigmaTPC(part, o2::track::PID::Pion), trackCuts.getNsigmaTOF(part, o2::track::PID::Pion), trackCuts.getNsigmaTPC(part, o2::track::PID::Kaon), trackCuts.getNsigmaTOF(part, o2::track::PID::Kaon), trackCuts.getNsigmaTPC(part, o2::track::PID::Electron))) { + if (!IsParticleNSigma((int8_t)1, part.p(), trackCuts.getNsigmaITS(part, o2::track::PID::Pion), trackCuts.getNsigmaTPC(part, o2::track::PID::Proton), trackCuts.getNsigmaTOF(part, o2::track::PID::Proton), trackCuts.getNsigmaTPC(part, o2::track::PID::Pion), trackCuts.getNsigmaTOF(part, o2::track::PID::Pion), trackCuts.getNsigmaTPC(part, o2::track::PID::Kaon), trackCuts.getNsigmaTOF(part, o2::track::PID::Kaon), trackCuts.getNsigmaTPC(part, o2::track::PID::Electron))) { continue; } trackHistoPartOne.fillQA(part); + if (twotracksconfigs.confIsAddITSNsigma) { + trackHistoPartOne.fillQAITSPID<>(part); + } trackHistoPartOne.fillQAMisIden(part, trackonefilter.confPDGCodePartOne); } } if ((contType == PairType::PlusMinus || contType == PairType::MinusMinus) && fillQA) { for (const auto& part : groupPartsTwo) { - if (!IsParticleNSigma((int8_t)2, part.p(), trackCuts.getNsigmaTPC(part, o2::track::PID::Proton), trackCuts.getNsigmaTOF(part, o2::track::PID::Proton), trackCuts.getNsigmaTPC(part, o2::track::PID::Pion), trackCuts.getNsigmaTOF(part, o2::track::PID::Pion), trackCuts.getNsigmaTPC(part, o2::track::PID::Kaon), trackCuts.getNsigmaTOF(part, o2::track::PID::Kaon), trackCuts.getNsigmaTPC(part, o2::track::PID::Electron))) { + if (!IsParticleNSigma((int8_t)2, part.p(), trackCuts.getNsigmaITS(part, o2::track::PID::Pion), trackCuts.getNsigmaTPC(part, o2::track::PID::Proton), trackCuts.getNsigmaTOF(part, o2::track::PID::Proton), trackCuts.getNsigmaTPC(part, o2::track::PID::Pion), trackCuts.getNsigmaTOF(part, o2::track::PID::Pion), trackCuts.getNsigmaTPC(part, o2::track::PID::Kaon), trackCuts.getNsigmaTOF(part, o2::track::PID::Kaon), trackCuts.getNsigmaTPC(part, o2::track::PID::Electron))) { continue; } trackHistoPartTwo.fillQA(part); + if (twotracksconfigs.confIsAddITSNsigma) { + trackHistoPartTwo.fillQAITSPID<>(part); + } } } @@ -557,11 +571,11 @@ struct femtoUniversePairTaskTrackTrackSpherHarMultKtExtended { /// Now build the combinations for non-identical particle pairs for (const auto& [p1, p2] : combinations(CombinationsFullIndexPolicy(groupPartsOne, groupPartsTwo))) { - if (!IsParticleNSigma((int8_t)1, p1.p(), trackCuts.getNsigmaTPC(p1, o2::track::PID::Proton), trackCuts.getNsigmaTOF(p1, o2::track::PID::Proton), trackCuts.getNsigmaTPC(p1, o2::track::PID::Pion), trackCuts.getNsigmaTOF(p1, o2::track::PID::Pion), trackCuts.getNsigmaTPC(p1, o2::track::PID::Kaon), trackCuts.getNsigmaTOF(p1, o2::track::PID::Kaon), trackCuts.getNsigmaTPC(p1, o2::track::PID::Electron))) { + if (!IsParticleNSigma((int8_t)1, p1.p(), trackCuts.getNsigmaITS(p1, o2::track::PID::Pion), trackCuts.getNsigmaTPC(p1, o2::track::PID::Proton), trackCuts.getNsigmaTOF(p1, o2::track::PID::Proton), trackCuts.getNsigmaTPC(p1, o2::track::PID::Pion), trackCuts.getNsigmaTOF(p1, o2::track::PID::Pion), trackCuts.getNsigmaTPC(p1, o2::track::PID::Kaon), trackCuts.getNsigmaTOF(p1, o2::track::PID::Kaon), trackCuts.getNsigmaTPC(p1, o2::track::PID::Electron))) { continue; } - if (!IsParticleNSigma((int8_t)2, p2.p(), trackCuts.getNsigmaTPC(p2, o2::track::PID::Proton), trackCuts.getNsigmaTOF(p2, o2::track::PID::Proton), trackCuts.getNsigmaTPC(p2, o2::track::PID::Pion), trackCuts.getNsigmaTOF(p2, o2::track::PID::Pion), trackCuts.getNsigmaTPC(p2, o2::track::PID::Kaon), trackCuts.getNsigmaTOF(p2, o2::track::PID::Kaon), trackCuts.getNsigmaTPC(p2, o2::track::PID::Electron))) { + if (!IsParticleNSigma((int8_t)2, p2.p(), trackCuts.getNsigmaITS(p2, o2::track::PID::Pion), trackCuts.getNsigmaTPC(p2, o2::track::PID::Proton), trackCuts.getNsigmaTOF(p2, o2::track::PID::Proton), trackCuts.getNsigmaTPC(p2, o2::track::PID::Pion), trackCuts.getNsigmaTOF(p2, o2::track::PID::Pion), trackCuts.getNsigmaTPC(p2, o2::track::PID::Kaon), trackCuts.getNsigmaTOF(p2, o2::track::PID::Kaon), trackCuts.getNsigmaTPC(p2, o2::track::PID::Electron))) { continue; } @@ -603,11 +617,11 @@ struct femtoUniversePairTaskTrackTrackSpherHarMultKtExtended { } else { for (const auto& [p1, p2] : combinations(CombinationsStrictlyUpperIndexPolicy(groupPartsOne, groupPartsOne))) { - if (!IsParticleNSigma((int8_t)2, p1.p(), trackCuts.getNsigmaTPC(p1, o2::track::PID::Proton), trackCuts.getNsigmaTOF(p1, o2::track::PID::Proton), trackCuts.getNsigmaTPC(p1, o2::track::PID::Pion), trackCuts.getNsigmaTOF(p1, o2::track::PID::Pion), trackCuts.getNsigmaTPC(p1, o2::track::PID::Kaon), trackCuts.getNsigmaTOF(p1, o2::track::PID::Kaon), trackCuts.getNsigmaTPC(p1, o2::track::PID::Electron))) { + if (!IsParticleNSigma((int8_t)2, p1.p(), trackCuts.getNsigmaITS(p1, o2::track::PID::Pion), trackCuts.getNsigmaTPC(p1, o2::track::PID::Proton), trackCuts.getNsigmaTOF(p1, o2::track::PID::Proton), trackCuts.getNsigmaTPC(p1, o2::track::PID::Pion), trackCuts.getNsigmaTOF(p1, o2::track::PID::Pion), trackCuts.getNsigmaTPC(p1, o2::track::PID::Kaon), trackCuts.getNsigmaTOF(p1, o2::track::PID::Kaon), trackCuts.getNsigmaTPC(p1, o2::track::PID::Electron))) { continue; } - if (!IsParticleNSigma((int8_t)2, p2.p(), trackCuts.getNsigmaTPC(p2, o2::track::PID::Proton), trackCuts.getNsigmaTOF(p2, o2::track::PID::Proton), trackCuts.getNsigmaTPC(p2, o2::track::PID::Pion), trackCuts.getNsigmaTOF(p2, o2::track::PID::Pion), trackCuts.getNsigmaTPC(p2, o2::track::PID::Kaon), trackCuts.getNsigmaTOF(p2, o2::track::PID::Kaon), trackCuts.getNsigmaTPC(p2, o2::track::PID::Electron))) { + if (!IsParticleNSigma((int8_t)2, p2.p(), trackCuts.getNsigmaITS(p2, o2::track::PID::Pion), trackCuts.getNsigmaTPC(p2, o2::track::PID::Proton), trackCuts.getNsigmaTOF(p2, o2::track::PID::Proton), trackCuts.getNsigmaTPC(p2, o2::track::PID::Pion), trackCuts.getNsigmaTOF(p2, o2::track::PID::Pion), trackCuts.getNsigmaTPC(p2, o2::track::PID::Kaon), trackCuts.getNsigmaTOF(p2, o2::track::PID::Kaon), trackCuts.getNsigmaTPC(p2, o2::track::PID::Electron))) { continue; } @@ -972,11 +986,11 @@ struct femtoUniversePairTaskTrackTrackSpherHarMultKtExtended { for (const auto& [p1, p2] : combinations(CombinationsFullIndexPolicy(groupPartsOne, groupPartsTwo))) { - if (!IsParticleNSigma((int8_t)2, p1.p(), trackCuts.getNsigmaTPC(p1, o2::track::PID::Proton), trackCuts.getNsigmaTOF(p1, o2::track::PID::Proton), trackCuts.getNsigmaTPC(p1, o2::track::PID::Pion), trackCuts.getNsigmaTOF(p1, o2::track::PID::Pion), trackCuts.getNsigmaTPC(p1, o2::track::PID::Kaon), trackCuts.getNsigmaTOF(p1, o2::track::PID::Kaon), trackCuts.getNsigmaTPC(p1, o2::track::PID::Electron))) { + if (!IsParticleNSigma((int8_t)2, p1.p(), trackCuts.getNsigmaITS(p1, o2::track::PID::Pion), trackCuts.getNsigmaTPC(p1, o2::track::PID::Proton), trackCuts.getNsigmaTOF(p1, o2::track::PID::Proton), trackCuts.getNsigmaTPC(p1, o2::track::PID::Pion), trackCuts.getNsigmaTOF(p1, o2::track::PID::Pion), trackCuts.getNsigmaTPC(p1, o2::track::PID::Kaon), trackCuts.getNsigmaTOF(p1, o2::track::PID::Kaon), trackCuts.getNsigmaTPC(p1, o2::track::PID::Electron))) { continue; } - if (!IsParticleNSigma((int8_t)2, p2.p(), trackCuts.getNsigmaTPC(p2, o2::track::PID::Proton), trackCuts.getNsigmaTOF(p2, o2::track::PID::Proton), trackCuts.getNsigmaTPC(p2, o2::track::PID::Pion), trackCuts.getNsigmaTOF(p2, o2::track::PID::Pion), trackCuts.getNsigmaTPC(p2, o2::track::PID::Kaon), trackCuts.getNsigmaTOF(p2, o2::track::PID::Kaon), trackCuts.getNsigmaTPC(p2, o2::track::PID::Electron))) { + if (!IsParticleNSigma((int8_t)2, p2.p(), trackCuts.getNsigmaITS(p2, o2::track::PID::Pion), trackCuts.getNsigmaTPC(p2, o2::track::PID::Proton), trackCuts.getNsigmaTOF(p2, o2::track::PID::Proton), trackCuts.getNsigmaTPC(p2, o2::track::PID::Pion), trackCuts.getNsigmaTOF(p2, o2::track::PID::Pion), trackCuts.getNsigmaTPC(p2, o2::track::PID::Kaon), trackCuts.getNsigmaTOF(p2, o2::track::PID::Kaon), trackCuts.getNsigmaTPC(p2, o2::track::PID::Electron))) { continue; }