diff --git a/PWGLF/TableProducer/Nuspex/nucleiFlowTree.cxx b/PWGLF/TableProducer/Nuspex/nucleiFlowTree.cxx index c336bd8a74c..beb4da99806 100644 --- a/PWGLF/TableProducer/Nuspex/nucleiFlowTree.cxx +++ b/PWGLF/TableProducer/Nuspex/nucleiFlowTree.cxx @@ -19,54 +19,48 @@ // o2-analysis-pid-tof-base, o2-analysis-multiplicity-table, o2-analysis-event-selection // (to add flow: o2-analysis-qvector-table, o2-analysis-centrality-table) -#include -#include -#include -#include -#include - -#include "Math/Vector4D.h" - -#include "CCDB/BasicCCDBManager.h" +#include "PWGLF/DataModel/EPCalibrationTables.h" +#include "PWGLF/DataModel/LFSlimNucleiTables.h" +#include "PWGLF/Utils/nucleiUtils.h" +#include "Common/Core/EventPlaneHelper.h" +#include "Common/Core/PID/PIDTOF.h" +#include "Common/Core/RecoDecay.h" #include "Common/Core/TrackSelection.h" #include "Common/Core/trackUtilities.h" #include "Common/DataModel/Centrality.h" -#include "Common/DataModel/Multiplicity.h" #include "Common/DataModel/EventSelection.h" +#include "Common/DataModel/Multiplicity.h" #include "Common/DataModel/PIDResponse.h" #include "Common/DataModel/PIDResponseITS.h" +#include "Common/DataModel/Qvectors.h" #include "Common/DataModel/TrackSelectionTables.h" -#include "Common/Core/PID/PIDTOF.h" #include "Common/TableProducer/PID/pidTOFBase.h" -#include "Common/Core/EventPlaneHelper.h" -#include "Common/DataModel/Qvectors.h" #include "Common/Tools/TrackTuner.h" -#include "Common/Core/RecoDecay.h" +#include "EventFiltering/Zorro.h" +#include "EventFiltering/ZorroSummary.h" +#include "CCDB/BasicCCDBManager.h" #include "DataFormatsParameters/GRPMagField.h" #include "DataFormatsParameters/GRPObject.h" #include "DataFormatsTPC/BetheBlochAleph.h" #include "DetectorsBase/GeometryManager.h" #include "DetectorsBase/Propagator.h" - -#include "EventFiltering/Zorro.h" -#include "EventFiltering/ZorroSummary.h" - +#include "Framework/ASoAHelpers.h" #include "Framework/AnalysisDataModel.h" #include "Framework/AnalysisTask.h" -#include "Framework/ASoAHelpers.h" #include "Framework/HistogramRegistry.h" #include "Framework/runDataProcessing.h" - #include "ReconstructionDataFormats/Track.h" -#include "PWGLF/DataModel/EPCalibrationTables.h" -#include "PWGLF/DataModel/LFSlimNucleiTables.h" - +#include "Math/Vector4D.h" #include "TRandom3.h" -#include "PWGLF/Utils/nucleiUtils.h" +#include +#include +#include +#include +#include using namespace o2; using namespace o2::framework; diff --git a/PWGLF/TableProducer/QC/nucleiQC.cxx b/PWGLF/TableProducer/QC/nucleiQC.cxx index afd71da2204..2672ef1029d 100644 --- a/PWGLF/TableProducer/QC/nucleiQC.cxx +++ b/PWGLF/TableProducer/QC/nucleiQC.cxx @@ -19,61 +19,54 @@ // o2-analysis-pid-tof-base, o2-analysis-multiplicity-table, o2-analysis-event-selection // (to add flow: o2-analysis-qvector-table, o2-analysis-centrality-table) -#include -#include -#include -#include -#include -#include - -#include "Math/Vector4D.h" - -#include "CCDB/BasicCCDBManager.h" +#include "PWGLF/DataModel/EPCalibrationTables.h" +#include "PWGLF/DataModel/LFSlimNucleiTables.h" +#include "PWGLF/Utils/nucleiUtils.h" +#include "Common/Core/EventPlaneHelper.h" +#include "Common/Core/PID/PIDTOF.h" +#include "Common/Core/RecoDecay.h" #include "Common/Core/TrackSelection.h" #include "Common/Core/trackUtilities.h" #include "Common/DataModel/Centrality.h" -#include "Common/DataModel/Multiplicity.h" #include "Common/DataModel/EventSelection.h" +#include "Common/DataModel/Multiplicity.h" #include "Common/DataModel/PIDResponse.h" #include "Common/DataModel/PIDResponseITS.h" +#include "Common/DataModel/Qvectors.h" #include "Common/DataModel/TrackSelectionTables.h" -#include "Common/Core/PID/PIDTOF.h" #include "Common/TableProducer/PID/pidTOFBase.h" -#include "Common/Core/EventPlaneHelper.h" -#include "Common/DataModel/Qvectors.h" #include "Common/Tools/TrackTuner.h" -#include "Common/Core/RecoDecay.h" +#include "EventFiltering/Zorro.h" +#include "EventFiltering/ZorroSummary.h" +#include "CCDB/BasicCCDBManager.h" #include "DataFormatsParameters/GRPMagField.h" #include "DataFormatsParameters/GRPObject.h" #include "DataFormatsTPC/BetheBlochAleph.h" #include "DetectorsBase/GeometryManager.h" #include "DetectorsBase/Propagator.h" - -#include "EventFiltering/Zorro.h" -#include "EventFiltering/ZorroSummary.h" - +#include "Framework/ASoAHelpers.h" #include "Framework/AnalysisDataModel.h" #include "Framework/AnalysisTask.h" -#include "Framework/ASoAHelpers.h" #include "Framework/HistogramRegistry.h" -#include "Framework/runDataProcessing.h" #include "Framework/StaticFor.h" - +#include "Framework/runDataProcessing.h" #include "ReconstructionDataFormats/Track.h" -#include "PWGLF/DataModel/EPCalibrationTables.h" -#include "PWGLF/DataModel/LFSlimNucleiTables.h" - +#include "Math/Vector4D.h" #include "TRandom3.h" -#include "PWGLF/Utils/nucleiUtils.h" +#include +#include +#include +#include +#include +#include using namespace o2; using namespace o2::framework; - struct nucleiQC { using TrackCandidates = soa::Join; @@ -109,7 +102,7 @@ struct nucleiQC { int mRunNumber = 0; float mBz = 0.f; - HistogramRegistry mHistograms { + HistogramRegistry mHistograms{ "histos", { {"hEventSelections", "Event selections; Selection step; Counts", {HistType::kTH1D, {{nuclei::evSel::kNevSels + 1, -0.5f, static_cast(nuclei::evSel::kNevSels) + 0.5f}}}}, @@ -118,19 +111,19 @@ struct nucleiQC { }, OutputObjHandlingPolicy::AnalysisObject, false, - true - }; + true}; std::vector mSpeciesToProcess; Produces mNucleiTableRed; - + std::vector mNucleiCandidates; std::vector mFilledMcParticleIds; - + o2::dataformats::DCA mDcaInfoCov; o2::track::TrackParametrizationWithError mTrackParCov; std::array(nuclei::Species::kNspecies)> mPidManagers; - void init(o2::framework::InitContext&) { + void init(o2::framework::InitContext&) + { mCcdb->setURL(cfgCCDBurl); mCcdb->setCaching(true); @@ -148,8 +141,7 @@ struct nucleiQC { } } - static_for<0, nuclei::kNspecies-1> ([&] (auto iSpecies) { - + static_for<0, nuclei::kNspecies - 1>([&](auto iSpecies) { constexpr int iSpeciesCt = decltype(iSpecies)::value; const int iSpeciesRt = iSpeciesCt; @@ -163,7 +155,7 @@ struct nucleiQC { } nuclei::createHistogramRegistryNucleus(mHistograms); - + if (cfgUseCentralTpcCalibration->get(uint32_t(iSpeciesRt), uint32_t(0)) == 0) { mPidManagers[iSpeciesRt] = nuclei::PidManager(iSpeciesRt, tpcBetheBlochParams); } else { @@ -188,8 +180,9 @@ struct nucleiQC { } template - bool trackSelection(const Ttrack& track) { - if (std::abs(track.eta()) > cfgCutEta || + bool trackSelection(const Ttrack& track) + { + if (std::abs(track.eta()) > cfgCutEta || track.tpcInnerParam() < cfgCutTpcMom || track.itsNCls() < cfgCutNclusITS || track.tpcNClsFound() < cfgCutNclusTPC || @@ -201,10 +194,10 @@ struct nucleiQC { } return true; } - + template - bool pidSelection(const Ttrack& track, const Tcollision& collision) - { + bool pidSelection(const Ttrack& track, const Tcollision& collision) + { constexpr int index = iSpecies; if (!nuclei::checkSpeciesValidity(index)) { std::runtime_error("species contains invalid nucleus index"); @@ -213,7 +206,8 @@ struct nucleiQC { const float centrality = nuclei::getCentrality(collision, cfgCentralityEstimator); const float nsigmaTPC = mPidManagers[index].getNSigmaTPC(track); mHistograms.fill(HIST(nuclei::cNames[index]) + HIST("/h3NsigmaTPC_preselectionVsCentrality"), track.pt(), nsigmaTPC, centrality); - if (std::abs(nsigmaTPC) > cfgNsigmaTPC->get(index, 1)) return false; + if (std::abs(nsigmaTPC) > cfgNsigmaTPC->get(index, 1)) + return false; mHistograms.fill(HIST(nuclei::cNames[index]) + HIST("/h3NsigmaTPCVsCentrality"), track.pt(), nsigmaTPC, centrality); const float nsigmaITS = mPidManagers[index].getNSigmaITS(track); @@ -223,54 +217,55 @@ struct nucleiQC { const float nsigmaTOF = mPidManagers[index].getNSigmaTOF(track); mHistograms.fill(HIST(nuclei::cNames[index]) + HIST("/h3NsigmaTOF_preselectionVsCentrality"), track.pt(), nsigmaTOF, centrality); - if (std::abs(nsigmaTOF) > cfgNsigmaTOF->get(index, 1)) return false; + if (std::abs(nsigmaTOF) > cfgNsigmaTOF->get(index, 1)) + return false; mHistograms.fill(HIST(nuclei::cNames[index]) + HIST("/h3NsigmaTOFVsCentrality"), track.pt(), nsigmaTOF, centrality); return true; } template - void fillNucleusFlagsPdgsMc(const Tparticle& particle, nuclei::SlimCandidate& candidate) + void fillNucleusFlagsPdgsMc(const Tparticle& particle, nuclei::SlimCandidate& candidate) { candidate.pdgCode = particle.pdgCode(); if (particle.isPhysicalPrimary()) { - candidate.flags |= nuclei::Flags::kIsPhysicalPrimary; - - // heavy flavour mother - //if (particle.has_mothers()) { - // for (auto& motherparticle : particle.mothers_as()) { - // if (std::find(nuclei::hfMothCodes.begin(), nuclei::hfMothCodes.end(), std::abs(motherparticle.pdgCode())) != nuclei::hfMothCodes.end()) { - // flags |= kIsSecondaryFromWeakDecay; - // motherPdgCode = motherparticle.pdgCode(); - // break; - // } - // } - //} - - } else if (particle.has_mothers()) { - candidate.flags |= nuclei::Flags::kIsSecondaryFromWeakDecay; - for (auto& motherparticle : particle.template mothers_as()) { - candidate.motherPdgCode = motherparticle.pdgCode(); - } - - } else { - candidate.flags |= nuclei::Flags::kIsSecondaryFromMaterial; + candidate.flags |= nuclei::Flags::kIsPhysicalPrimary; + + // heavy flavour mother + // if (particle.has_mothers()) { + // for (auto& motherparticle : particle.mothers_as()) { + // if (std::find(nuclei::hfMothCodes.begin(), nuclei::hfMothCodes.end(), std::abs(motherparticle.pdgCode())) != nuclei::hfMothCodes.end()) { + // flags |= kIsSecondaryFromWeakDecay; + // motherPdgCode = motherparticle.pdgCode(); + // break; + // } + // } + //} + + } else if (particle.has_mothers()) { + candidate.flags |= nuclei::Flags::kIsSecondaryFromWeakDecay; + for (auto& motherparticle : particle.template mothers_as()) { + candidate.motherPdgCode = motherparticle.pdgCode(); } - mFilledMcParticleIds.emplace_back(particle.globalIndex()); + } else { + candidate.flags |= nuclei::Flags::kIsSecondaryFromMaterial; + } + + mFilledMcParticleIds.emplace_back(particle.globalIndex()); } template - void fillNucleusFlagsPdgs(const int iSpecies, const Tcollision& collision, const Ttrack& track, nuclei::SlimCandidate& candidate) + void fillNucleusFlagsPdgs(const int iSpecies, const Tcollision& collision, const Ttrack& track, nuclei::SlimCandidate& candidate) { candidate.flags = static_cast((track.pidForTracking() & 0xF) << 12); - candidate.flags |= iSpecies == nuclei::Species::kPr ? nuclei::Flags::kProton : - iSpecies == nuclei::Species::kDe ? nuclei::Flags::kDeuteron : - iSpecies == nuclei::Species::kTr ? nuclei::Flags::kTriton : - iSpecies == nuclei::Species::kHe ? nuclei::Flags::kHe3 : - iSpecies == nuclei::Species::kAl ? nuclei::Flags::kHe4 : 0; - + candidate.flags |= iSpecies == nuclei::Species::kPr ? nuclei::Flags::kProton : iSpecies == nuclei::Species::kDe ? nuclei::Flags::kDeuteron + : iSpecies == nuclei::Species::kTr ? nuclei::Flags::kTriton + : iSpecies == nuclei::Species::kHe ? nuclei::Flags::kHe3 + : iSpecies == nuclei::Species::kAl ? nuclei::Flags::kHe4 + : 0; + if (track.hasTOF()) { candidate.flags |= nuclei::Flags::kHasTOF; } @@ -291,7 +286,8 @@ struct nucleiQC { } template - void fillDcaInformation(const Tcollision& collision, const Ttrack& track, nuclei::SlimCandidate& candidate) { + void fillDcaInformation(const Tcollision& collision, const Ttrack& track, nuclei::SlimCandidate& candidate) + { const o2::math_utils::Point3D collisionVertex{collision.posX(), collision.posY(), collision.posZ()}; @@ -303,7 +299,6 @@ struct nucleiQC { candidate.DCAxy = dcaInfo[0]; candidate.DCAz = dcaInfo[1]; - } template @@ -327,35 +322,40 @@ struct nucleiQC { .ptGenerated = 0.f, // to be filled for mc .etaGenerated = 0.f, .phiGenerated = 0.f, - .centrality = nuclei::getCentrality(collision, cfgCentralityEstimator) - }; - + .centrality = nuclei::getCentrality(collision, cfgCentralityEstimator)}; + fillDcaInformation(collision, track, candidate); fillNucleusFlagsPdgs(iSpecies, collision, track, candidate); - + if constexpr (isMc) { if (track.has_mcParticle()) { - + const auto& particle = track.mcParticle(); fillNucleusFlagsPdgsMc(particle, candidate); fillNucleusGeneratedVariables(particle, candidate); } } - + mNucleiCandidates.emplace_back(candidate); return candidate; - } - template - void dispatchFillHistograms(const int iSpecies, const nuclei::SlimCandidate& candidate) { + template + void dispatchFillHistograms(const int iSpecies, const nuclei::SlimCandidate& candidate) + { switch (iSpecies) { - case nuclei::Species::kPr: return fillHistograms(candidate); - case nuclei::Species::kDe: return fillHistograms(candidate); - case nuclei::Species::kTr: return fillHistograms(candidate); - case nuclei::Species::kHe: return fillHistograms(candidate); - case nuclei::Species::kAl: return fillHistograms(candidate); - default: return; + case nuclei::Species::kPr: + return fillHistograms(candidate); + case nuclei::Species::kDe: + return fillHistograms(candidate); + case nuclei::Species::kTr: + return fillHistograms(candidate); + case nuclei::Species::kHe: + return fillHistograms(candidate); + case nuclei::Species::kAl: + return fillHistograms(candidate); + default: + return; } } @@ -366,7 +366,7 @@ struct nucleiQC { if (!nuclei::checkSpeciesValidity(index)) { std::runtime_error("species contains invalid nucleus index"); } - + if (isGenerated) { const float ptGenerated = (index == nuclei::Species::kPr || index == nuclei::Species::kDe || index == nuclei::Species::kTr) ? candidate.ptGenerated : candidate.ptGenerated / 2.f; mHistograms.fill(HIST(nuclei::cNames[index]) + HIST("/hPtGenerated"), ptGenerated); @@ -388,7 +388,7 @@ struct nucleiQC { { mNucleiCandidates.clear(); mFilledMcParticleIds.clear(); - + auto bc = collision.template bc_as(); initCCDB(bc); @@ -396,11 +396,9 @@ struct nucleiQC { return; } - for (const auto& track: tracks) - { - - static_for<0, nuclei::kNspecies-1> ([&] (auto iSpecies) { + for (const auto& track : tracks) { + static_for<0, nuclei::kNspecies - 1>([&](auto iSpecies) { constexpr int iSpeciesCt = decltype(iSpecies)::value; const int iSpeciesRt = iSpeciesCt; @@ -411,37 +409,39 @@ struct nucleiQC { if (track.has_mcParticle()) { const auto& particle = track.mcParticle(); if (cfgDoCheckPdgCode) { - if (particle.pdgCode() != nuclei::pdgCodes[iSpeciesRt]) return; + if (particle.pdgCode() != nuclei::pdgCodes[iSpeciesRt]) + return; } if (particle.y() - cfgRapidityCenterMass < cfgRapidityMin || particle.y() - cfgRapidityCenterMass > cfgRapidityMax) { return; } - } + } mHistograms.fill(HIST(nuclei::cNames[iSpeciesCt]) + HIST("/hTrackSelections"), nuclei::trackSelection::kNoCuts); - if (!trackSelection(track)) return; + if (!trackSelection(track)) + return; mHistograms.fill(HIST(nuclei::cNames[iSpeciesCt]) + HIST("/hTrackSelections"), nuclei::trackSelection::kTrackCuts); - if (!pidSelection(track, collision)) return; + if (!pidSelection(track, collision)) + return; mHistograms.fill(HIST(nuclei::cNames[iSpeciesCt]) + HIST("/hTrackSelections"), nuclei::trackSelection::kPidCuts); nuclei::SlimCandidate candidate; if (track.has_mcParticle()) { - candidate = fillCandidate(iSpeciesCt, collision, track); - dispatchFillHistograms(iSpeciesRt, candidate); + candidate = fillCandidate(iSpeciesCt, collision, track); + dispatchFillHistograms(iSpeciesRt, candidate); } else { - candidate = fillCandidate(iSpeciesCt, collision, track); + candidate = fillCandidate(iSpeciesCt, collision, track); } - dispatchFillHistograms(iSpeciesRt, candidate); + dispatchFillHistograms(iSpeciesRt, candidate); }); } - + const int mcCollisionId = collision.mcCollisionId(); auto mcParticlesThisCollision = mcParticles.sliceBy(mMcParticlesPerCollision, mcCollisionId); mcParticlesThisCollision.bindExternalIndices(&mcParticles); - for (const auto& particle: mcParticlesThisCollision) - { + for (const auto& particle : mcParticlesThisCollision) { if (std::find(mFilledMcParticleIds.begin(), mFilledMcParticleIds.end(), particle.globalIndex()) != mFilledMcParticleIds.end()) { continue; } @@ -455,13 +455,13 @@ struct nucleiQC { fillNucleusGeneratedVariables(particle, candidate); mNucleiCandidates.emplace_back(candidate); - dispatchFillHistograms(iSpecies, candidate); + dispatchFillHistograms(iSpecies, candidate); } - if (!cfgFillTable) return; + if (!cfgFillTable) + return; - for (const auto& candidate: mNucleiCandidates) - { + for (const auto& candidate : mNucleiCandidates) { mNucleiTableRed( candidate.pt, candidate.eta, @@ -474,8 +474,7 @@ struct nucleiQC { candidate.DCAz, candidate.flags, candidate.pdgCode, - candidate.motherPdgCode - ); + candidate.motherPdgCode); } } PROCESS_SWITCH(nucleiQC, processMc, "Mc analysis", false); diff --git a/PWGLF/Utils/nucleiUtils.h b/PWGLF/Utils/nucleiUtils.h index eeb70aa51b2..e954ef76c1e 100644 --- a/PWGLF/Utils/nucleiUtils.h +++ b/PWGLF/Utils/nucleiUtils.h @@ -12,12 +12,12 @@ #ifndef PWGLF_UTILS_NUCLEIUTILS_H_ #define PWGLF_UTILS_NUCLEIUTILS_H_ -#include -#include - #include "Framework/HistogramRegistry.h" #include "Framework/HistogramSpec.h" +#include +#include + using namespace o2; using namespace o2::framework; using namespace o2::framework::expressions; @@ -119,18 +119,26 @@ enum Flags { kIsSecondaryFromWeakDecay = BIT(11) /// the last 4 bits are reserved for the PID in tracking }; -constexpr int getSpeciesFromPdg(int pdg) { +constexpr int getSpeciesFromPdg(int pdg) +{ switch (std::abs(pdg)) { - case PDG_t::kProton: return Species::kPr; - case o2::constants::physics::Pdg::kDeuteron: return Species::kDe; - case o2::constants::physics::Pdg::kTriton: return Species::kTr; - case o2::constants::physics::Pdg::kHelium3: return Species::kHe; - case o2::constants::physics::Pdg::kAlpha: return Species::kAl; - default: return -1; + case PDG_t::kProton: + return Species::kPr; + case o2::constants::physics::Pdg::kDeuteron: + return Species::kDe; + case o2::constants::physics::Pdg::kTriton: + return Species::kTr; + case o2::constants::physics::Pdg::kHelium3: + return Species::kHe; + case o2::constants::physics::Pdg::kAlpha: + return Species::kAl; + default: + return -1; } } -bool checkSpeciesValidity(const int species) { +bool checkSpeciesValidity(const int species) +{ if (species < 0 || species > Species::kNspecies) { return false; } @@ -356,8 +364,9 @@ enum trackSelection { }; static const std::array(trackSelection::kNtrackSelections)> trackSelectionLabels{"All", "Track cuts", "PID cuts"}; -template -void createHistogramRegistryNucleus(HistogramRegistry& registry) { +template +void createHistogramRegistryNucleus(HistogramRegistry& registry) +{ constexpr int index = iSpecies; if (!checkSpeciesValidity(index)) { @@ -369,7 +378,7 @@ void createHistogramRegistryNucleus(HistogramRegistry& registry) { registry.add(fmt::format("{}/h3PtVsEtaVsCentralityReconstructed", cNames[index]).c_str(), (fmt::format("{} - reconstructed variables;", cNames[index]) + std::string("#it{p}_{T} / |#it{Z}| (GeV/#it{c}); #eta; CentralityFT0C (%)")).c_str(), HistType::kTH3F, {{240, -6.0f, 6.0f}, {40, -1.0f, 1.f}, {20, 0.0f, 100.0f}}); registry.add(fmt::format("{}/h3PhiVsEtaVsCentralityReconstructed", cNames[index]).c_str(), (fmt::format("{} - reconstructed variables;", cNames[index]) + std::string("#phi (radians); #eta; CentralityFT0C (%)")).c_str(), HistType::kTH3F, {{240, -6.0f, 6.0f}, {40, 0.0f, o2::constants::math::TwoPI}, {20, 0.0f, 100.0f}}); registry.add(fmt::format("{}/h3DCAxyVsPtVsCentrality", cNames[index]).c_str(), (fmt::format(";", cNames[index]) + std::string("#it{p}_{T} / |#it{Z}| (GeV/#it{c}); DCA_{xy} (cm); CentralityFT0C (%)")).c_str(), HistType::kTH3F, {{240, -6.0f, 6.0f}, {200, -5.0f, 5.0f}, {20, 0.0f, 100.0f}}); - registry.add(fmt::format("{}/h3DCAzVsPtVsCentrality", cNames[index]).c_str(), (fmt::format("{};", cNames[index]) + std::string("#it{p}_{T} / |#it{Z}| (GeV/#it{c}); DCA_{z} (cm); CentralityFT0C (%)")).c_str() , HistType::kTH3F, {{240, -6.0f, 6.0f}, {200, -5.0f, 5.0f}, {20, 0.0f, 100.0f}}); + registry.add(fmt::format("{}/h3DCAzVsPtVsCentrality", cNames[index]).c_str(), (fmt::format("{};", cNames[index]) + std::string("#it{p}_{T} / |#it{Z}| (GeV/#it{c}); DCA_{z} (cm); CentralityFT0C (%)")).c_str(), HistType::kTH3F, {{240, -6.0f, 6.0f}, {200, -5.0f, 5.0f}, {20, 0.0f, 100.0f}}); registry.add(fmt::format("{}/h3NsigmaTPC_preselectionVsCentrality", cNames[index]).c_str(), (fmt::format("Nsigma{} TPC distribution;", cNames[index]) + std::string("#it{p}_{T} / |#it{Z}| (GeV/#it{c});") + fmt::format("n#sigma_{{TPC}}({}); CentralityFT0C (%)", cNames[index])).c_str(), HistType::kTH3F, {{100, -5.0f, 5.0f}, {400, -10.0f, 10.0f}, {20, 0.0f, 100.0f}}); registry.add(fmt::format("{}/h3NsigmaTPCVsCentrality", cNames[index]).c_str(), (fmt::format("Nsigma{} TPC distribution;", cNames[index]) + std::string("#it{p}_{T} / |#it{Z}| (GeV/#it{c});") + fmt::format("n#sigma_{{TPC}}({}); Centrality FT0C (%)", cNames[index])).c_str(), HistType::kTH3F, {{20, -5.0f, 5.0f}, {200, -5.0f, 5.0f}, {20, 0.0f, 100.0f}}); registry.add(fmt::format("{}/h3NsigmaITS_preselectionVsCentrality", cNames[index]).c_str(), (fmt::format("Nsigma{} ITS distribution;", cNames[index]) + std::string("signed #it{p}_{T} / |#it{Z}| (GeV/#it{c});") + fmt::format("n#sigma_{{ITS}}({}); Centrality FT0C (%)", cNames[index])).c_str(), HistType::kTH3F, {{50, -5.0f, 5.0f}, {120, -3.0f, 3.0f}, {20, 0.0f, 100.0f}}); @@ -390,141 +399,159 @@ void createHistogramRegistryNucleus(HistogramRegistry& registry) { // PID manager class -class PidManager { - - public: - PidManager(const int species, const float* tpcBetheBlochParams = nullptr) +class PidManager +{ + + public: + PidManager(const int species, const float* tpcBetheBlochParams = nullptr) : mSpecies(species) - { - if (!checkSpeciesValidity(species)) { - std::runtime_error("species contains invalid nucleus index"); - } - - if (!tpcBetheBlochParams) { - mUseTpcCentralCalibration = true; - return; - } - - for (int i = 0; i < 6; i++) { - mTpcBetheBlochParams[i] = tpcBetheBlochParams[i]; - } + { + if (!checkSpeciesValidity(species)) { + std::runtime_error("species contains invalid nucleus index"); } - PidManager() = default; - ~PidManager() = default; - - // TOF - template - float getBetaTOF(const Ttrack& track) - { - if (!track.hasTOF()) { - return -999.f; - } - float beta = o2::pid::tof::Beta::GetBeta(track); - return std::min(1.f - 1.e-6f, std::max(1.e-4f, beta)); /// sometimes beta > 1 or < 0, to be checked + if (!tpcBetheBlochParams) { + mUseTpcCentralCalibration = true; + return; } - template - float getMassTOF(const Ttrack& track) - { - if (!track.hasTOF()) { - return -999.f; - } - const float charge{1.f + static_cast(mSpecies == Species::kHe || mSpecies == Species::kAl)}; - const float beta = getBetaTOF(track); - return track.tpcInnerParam() * charge * std::sqrt(1.f / (beta * beta) - 1.f); + for (int i = 0; i < 6; i++) { + mTpcBetheBlochParams[i] = tpcBetheBlochParams[i]; } - - template - float getNSigmaTOF(const Ttrack& track) - { - if (!track.hasTOF()) { - return -999.f; - } - - switch (mSpecies) { - case Species::kPr: return track.tofNSigmaPr(); - case Species::kDe: return track.tofNSigmaDe(); - case Species::kTr: return track.tofNSigmaTr(); - case Species::kHe: return track.tofNSigmaHe(); - case Species::kAl: return track.tofNSigmaAl(); - default: return -999.f; - } + } + PidManager() = default; + ~PidManager() = default; + + // TOF + template + float getBetaTOF(const Ttrack& track) + { + if (!track.hasTOF()) { + return -999.f; } + float beta = o2::pid::tof::Beta::GetBeta(track); + return std::min(1.f - 1.e-6f, std::max(1.e-4f, beta)); /// sometimes beta > 1 or < 0, to be checked + } - template - float getClusterSizeCosLambdaITS(const Ttrack& track) - { - return mResponseITS.averageClusterSize(track.itsClusterSizes()) / std::cosh(track.eta()); + template + float getMassTOF(const Ttrack& track) + { + if (!track.hasTOF()) { + return -999.f; } + const float charge{1.f + static_cast(mSpecies == Species::kHe || mSpecies == Species::kAl)}; + const float beta = getBetaTOF(track); + return track.tpcInnerParam() * charge * std::sqrt(1.f / (beta * beta) - 1.f); + } - float getClusterSizeCosLambdaITS(const u_int32_t clusterSizesITS, const float eta) - { - return mResponseITS.averageClusterSize(clusterSizesITS) / std::cosh(eta); + template + float getNSigmaTOF(const Ttrack& track) + { + if (!track.hasTOF()) { + return -999.f; } - template - float getNSigmaITS(const Ttrack& track) - { - switch (mSpecies) { - case Species::kPr: return mResponseITS.nSigmaITS(track.itsClusterSizes(), track.p(), track.eta()); - case Species::kDe: return mResponseITS.nSigmaITS(track.itsClusterSizes(), track.p(), track.eta()); - case Species::kTr: return mResponseITS.nSigmaITS(track.itsClusterSizes(), track.p(), track.eta()); - case Species::kHe: return mResponseITS.nSigmaITS(track.itsClusterSizes(), 2. * track.p(), track.eta()); - case Species::kAl: return mResponseITS.nSigmaITS(track.itsClusterSizes(), 2. * track.p(), track.eta()); - default: return -999.f; - } + switch (mSpecies) { + case Species::kPr: + return track.tofNSigmaPr(); + case Species::kDe: + return track.tofNSigmaDe(); + case Species::kTr: + return track.tofNSigmaTr(); + case Species::kHe: + return track.tofNSigmaHe(); + case Species::kAl: + return track.tofNSigmaAl(); + default: + return -999.f; } + } - // TPC - float getExpectedTPCsignal(const float p) - { - if (!mUseTpcCentralCalibration) { + template + float getClusterSizeCosLambdaITS(const Ttrack& track) + { + return mResponseITS.averageClusterSize(track.itsClusterSizes()) / std::cosh(track.eta()); + } + + float getClusterSizeCosLambdaITS(const u_int32_t clusterSizesITS, const float eta) + { + return mResponseITS.averageClusterSize(clusterSizesITS) / std::cosh(eta); + } + + template + float getNSigmaITS(const Ttrack& track) + { + switch (mSpecies) { + case Species::kPr: + return mResponseITS.nSigmaITS(track.itsClusterSizes(), track.p(), track.eta()); + case Species::kDe: + return mResponseITS.nSigmaITS(track.itsClusterSizes(), track.p(), track.eta()); + case Species::kTr: + return mResponseITS.nSigmaITS(track.itsClusterSizes(), track.p(), track.eta()); + case Species::kHe: + return mResponseITS.nSigmaITS(track.itsClusterSizes(), 2. * track.p(), track.eta()); + case Species::kAl: + return mResponseITS.nSigmaITS(track.itsClusterSizes(), 2. * track.p(), track.eta()); + default: return -999.f; - } - float pScaled = p * mMomScaling[0] + mMomScaling[1]; - float betaGamma = pScaled / masses[mSpecies]; - return tpc::BetheBlochAleph(betaGamma, - mTpcBetheBlochParams[0], - mTpcBetheBlochParams[1], - mTpcBetheBlochParams[2], - mTpcBetheBlochParams[3], - mTpcBetheBlochParams[4]); } - - template - float getNSigmaTPC(const Ttrack& track) - { - if (!mUseTpcCentralCalibration) { - return getNSigmaTPCcentral(track); - } - float expectedSignal = getExpectedTPCsignal(track.tpcInnerParam()); - float resolution = mTpcBetheBlochParams[5]; - return (track.tpcSignal() - expectedSignal) / (expectedSignal * resolution); + } + + // TPC + float getExpectedTPCsignal(const float p) + { + if (!mUseTpcCentralCalibration) { + return -999.f; } + float pScaled = p * mMomScaling[0] + mMomScaling[1]; + float betaGamma = pScaled / masses[mSpecies]; + return tpc::BetheBlochAleph(betaGamma, + mTpcBetheBlochParams[0], + mTpcBetheBlochParams[1], + mTpcBetheBlochParams[2], + mTpcBetheBlochParams[3], + mTpcBetheBlochParams[4]); + } - protected: - // TPC - template - float getNSigmaTPCcentral(const Ttrack& track) { - switch (mSpecies) { - case Species::kPr: return track.tpcNSigmaPr(); - case Species::kDe: return track.tpcNSigmaDe(); - case Species::kTr: return track.tpcNSigmaTr(); - case Species::kHe: return track.tpcNSigmaHe(); - case Species::kAl: return track.tpcNSigmaAl(); - default: return -999.f; - } + template + float getNSigmaTPC(const Ttrack& track) + { + if (!mUseTpcCentralCalibration) { + return getNSigmaTPCcentral(track); } + float expectedSignal = getExpectedTPCsignal(track.tpcInnerParam()); + float resolution = mTpcBetheBlochParams[5]; + return (track.tpcSignal() - expectedSignal) / (expectedSignal * resolution); + } - private: - float mTpcBetheBlochParams[6]; - bool mUseTpcCentralCalibration = true; // this just becomes a check for the null pointer in the parameters - o2::aod::ITSResponse mResponseITS; - float mMomScaling[2]{1., 0.}; - int mSpecies; -}; + protected: + // TPC + template + float getNSigmaTPCcentral(const Ttrack& track) + { + switch (mSpecies) { + case Species::kPr: + return track.tpcNSigmaPr(); + case Species::kDe: + return track.tpcNSigmaDe(); + case Species::kTr: + return track.tpcNSigmaTr(); + case Species::kHe: + return track.tpcNSigmaHe(); + case Species::kAl: + return track.tpcNSigmaAl(); + default: + return -999.f; + } + } + private: + float mTpcBetheBlochParams[6]; + bool mUseTpcCentralCalibration = true; // this just becomes a check for the null pointer in the parameters + o2::aod::ITSResponse mResponseITS; + float mMomScaling[2]{1., 0.}; + int mSpecies; +}; } // namespace nuclei