diff --git a/PWGLF/Tasks/QC/CMakeLists.txt b/PWGLF/Tasks/QC/CMakeLists.txt index 58248ec3b6a..c640ee82586 100644 --- a/PWGLF/Tasks/QC/CMakeLists.txt +++ b/PWGLF/Tasks/QC/CMakeLists.txt @@ -108,3 +108,8 @@ o2physics_add_dpl_workflow(tracked-cascade-properties SOURCES tracked_cascade_properties.cxx PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(mc-particle-predictions + SOURCES mcParticlePrediction.cxx + PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore + COMPONENT_NAME Analysis) diff --git a/PWGLF/Tasks/QC/mcParticlePrediction.cxx b/PWGLF/Tasks/QC/mcParticlePrediction.cxx new file mode 100644 index 00000000000..0df9dc44c9b --- /dev/null +++ b/PWGLF/Tasks/QC/mcParticlePrediction.cxx @@ -0,0 +1,192 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// +/// \file mcParticlePrediction.cxx +/// \author Nicolò Jacazio nicolo.jacazio@cern.ch +/// \author Francesca Ercolessi francesca.ercolessi@cern.ch +/// \brief Task to build the predictions from the models based on the generated particles +/// + +#include "Framework/runDataProcessing.h" +#include "Framework/AnalysisTask.h" +#include "Framework/HistogramRegistry.h" +#include "Framework/StaticFor.h" +#include "Framework/O2DatabasePDGPlugin.h" +#include "PWGLF/Utils/mcParticle.h" +#include "PWGLF/Utils/inelGt.h" + +#include "TPDGCode.h" + +using namespace o2; +using namespace o2::framework; + +static const std::vector parameterNames{"Enable"}; +static constexpr int nParameters = 1; +static const int defaultParameters[o2::pwglf::PIDExtended::NIDsTot][nParameters]{{0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}, {0}}; +bool enabledArray[o2::pwglf::PIDExtended::NIDsTot]; + +// Histograms +const int FT0A = 0; +const int FT0C = 1; +const int FT0AC = 2; +const int FV0A = 3; +const int FDDA = 4; +const int FDDC = 5; +const int FDDAC = 6; +const int ZNA = 7; +const int ZNC = 8; +// const int ZEM1 = 9; +// const int ZEM2 = 10; +// const int ZPA = 11; +// const int ZPC = 12; +const int nEstimators = 13; + +const char* estimatorNames[nEstimators] = {"FT0A", + "FT0C", + "FT0AC", + "FV0A", + "FDDA", + "FDDC", + "FDDAC", + "ZNA", + "ZNC", + "ZEM1", + "ZEM2", + "ZPA", + "ZPC"}; + +std::array, o2::pwglf::PIDExtended::NIDsTot>, nEstimators> hpt; +std::array, o2::pwglf::PIDExtended::NIDsTot>, nEstimators> hyield; + +struct mcParticlePrediction { + + // Histograms + HistogramRegistry histos{"Histos", {}, OutputObjHandlingPolicy::AnalysisObject}; + ConfigurableAxis binsPt{"binsPt", {100, 0, 10}, "Binning of the Pt axis"}; + ConfigurableAxis binsMultiplicity{"binsMultiplicity", {100, 0, 1000}, "Binning of the Multiplicity axis"}; + Configurable> enabledSpecies{"enabledSpecies", + {defaultParameters[0], o2::pwglf::PIDExtended::NIDsTot, nParameters, o2::pwglf::PIDExtended::arrayNames(), parameterNames}, + "Bethe Bloch parameters"}; + Configurable selectInelGt0{"selectInelGt0", true, "Select only inelastic events"}; + Service pdgDB; + + void init(o2::framework::InitContext&) + { + histos.add("collisions", "collisions", kTH1D, {{10, 0, 10}}); + auto h = histos.add("particles", "particles", kTH1D, {{o2::pwglf::PIDExtended::NIDsTot, -0.5, -0.5 + o2::pwglf::PIDExtended::NIDsTot}}); + for (int i = 0; i < o2::pwglf::PIDExtended::NIDsTot; i++) { + h->GetXaxis()->SetBinLabel(i + 1, o2::pwglf::PIDExtended::getName(i)); + } + h = histos.add("particlesPrim", "particlesPrim", kTH1D, {{o2::pwglf::PIDExtended::NIDsTot, -0.5, -0.5 + o2::pwglf::PIDExtended::NIDsTot}}); + for (int i = 0; i < o2::pwglf::PIDExtended::NIDsTot; i++) { + h->GetXaxis()->SetBinLabel(i + 1, o2::pwglf::PIDExtended::getName(i)); + } + + const AxisSpec axisPt{binsPt, "#it{p}_{T} (GeV/#it{c})"}; + const AxisSpec axisMultiplicity{binsMultiplicity, "Mult"}; + // FT0 + histos.add("multiplicity/FT0A", "FT0A", kTH1D, {axisMultiplicity}); + histos.add("multiplicity/FT0C", "FT0C", kTH1D, {axisMultiplicity}); + histos.add("multiplicity/FT0AC", "FT0AC", kTH1D, {axisMultiplicity}); + + // FV0 + histos.add("multiplicity/FV0A", "FV0A", kTH1D, {axisMultiplicity}); + + // FDD + histos.add("multiplicity/FDDA", "FDDA", kTH1D, {axisMultiplicity}); + histos.add("multiplicity/FDDC", "FDDC", kTH1D, {axisMultiplicity}); + histos.add("multiplicity/FDDAC", "FDDAC", kTH1D, {axisMultiplicity}); + + // ZDC + histos.add("multiplicity/ZNA", "ZNA", kTH1D, {axisMultiplicity}); + histos.add("multiplicity/ZNC", "ZNC", kTH1D, {axisMultiplicity}); + histos.add("multiplicity/ZEM1", "ZEM1", kTH1D, {axisMultiplicity}); + histos.add("multiplicity/ZEM2", "ZEM2", kTH1D, {axisMultiplicity}); + histos.add("multiplicity/ZPA", "ZPA", kTH1D, {axisMultiplicity}); + histos.add("multiplicity/ZPC", "ZPC", kTH1D, {axisMultiplicity}); + + // ITS + histos.add("multiplicity/ITS", "ITS", kTH1D, {axisMultiplicity}); + + for (int i = 0; i < o2::pwglf::PIDExtended::NIDsTot; i++) { + if (enabledSpecies->get(o2::pwglf::PIDExtended::getName(i), "Enable") != 1) { + enabledArray[i] = false; + continue; + } + enabledArray[i] = true; + for (int j = 0; j < nEstimators; j++) { + hpt[j][i] = histos.add(Form("pt/%s/%s", estimatorNames[j], o2::pwglf::PIDExtended::getName(i)), o2::pwglf::PIDExtended::getName(i), kTH2D, {binsPt, axisMultiplicity}); + hyield[j][i] = histos.add(Form("yield/%s/%s", estimatorNames[j], o2::pwglf::PIDExtended::getName(i)), o2::pwglf::PIDExtended::getName(i), kTH1D, {axisMultiplicity}); + } + } + } + + int countInAcceptance(const aod::McParticles& mcParticles, const float etamin, const float etamax) + { + // static_assert(etamin < etamax, "etamin must be smaller than etamax"); + int counter = 0; + for (const auto& particle : mcParticles) { + if (particle.eta() > etamin && particle.eta() < etamax) { + counter++; + } + } + return counter; + } + + int countFT0A(const aod::McParticles& mcParticles) { return countInAcceptance(mcParticles, 3.5f, 4.9f); } + int countFT0C(const aod::McParticles& mcParticles) { return countInAcceptance(mcParticles, -3.3f, -2.1f); } + int countFV0A(const aod::McParticles& mcParticles) { return countInAcceptance(mcParticles, 2.2f, 5.1f); } + int countFDDA(const aod::McParticles& mcParticles) { return countInAcceptance(mcParticles, 4.9f, 6.3f); } + int countFDDC(const aod::McParticles& mcParticles) { return countInAcceptance(mcParticles, -7.f, -4.9f); } + int countZNA(const aod::McParticles& mcParticles) { return countInAcceptance(mcParticles, 8.8f, 100.f); } + int countZNC(const aod::McParticles& mcParticles) { return countInAcceptance(mcParticles, -100.f, -8.8f); } + + void process(aod::McCollision const& mcCollision, + aod::McParticles const& mcParticles) + { + if (selectInelGt0.value && !o2::pwglf::isINELgt0mc(mcParticles, pdgDB)) { + return; + } + + histos.fill(HIST("collisions"), 0.5); + const int nFT0A = countFT0A(mcParticles); + const int nFT0C = countFT0C(mcParticles); + const int nFV0A = countFV0A(mcParticles); + const int nFDDA = countFDDA(mcParticles); + const int nFDDC = countFDDC(mcParticles); + const int nZNA = countZNA(mcParticles); + const int nZNC = countZNC(mcParticles); + + for (const auto& particle : mcParticles) { + particle.pdgCode(); + const auto id = o2::pwglf::PIDExtended::pdgToId(particle); + if (id < 0) { + continue; + } + if (!enabledArray[id]) { + continue; + } + hpt[FT0A][id]->Fill(particle.pt(), nFT0A); + hpt[FT0C][id]->Fill(particle.pt(), nFT0C); + hpt[FT0AC][id]->Fill(particle.pt(), nFT0A + nFT0C); + + hpt[FV0A][id]->Fill(particle.pt(), nFV0A); + hpt[FDDA][id]->Fill(particle.pt(), nFDDA); + hpt[FDDC][id]->Fill(particle.pt(), nFDDC); + hpt[FDDAC][id]->Fill(particle.pt(), nFDDA + nFDDC); + hpt[ZNA][id]->Fill(particle.pt(), nZNA); + hpt[ZNC][id]->Fill(particle.pt(), nZNC); + } + } +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { return WorkflowSpec{adaptAnalysisTask(cfgc)}; } diff --git a/PWGLF/Utils/mcParticle.h b/PWGLF/Utils/mcParticle.h new file mode 100644 index 00000000000..822496f9535 --- /dev/null +++ b/PWGLF/Utils/mcParticle.h @@ -0,0 +1,320 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. + +/// +/// \file mcParticle.h +/// \author Nicolò Jacazio nicolo.jacazio@cern.ch +/// \author Francesca Ercolessi francesca.ercolessi@cern.ch +/// \since 31/05/2024 +/// \brief Utilities to handle the MC information +/// + +#ifndef PWGLF_UTILS_MCPARTICLE_H_ +#define PWGLF_UTILS_MCPARTICLE_H_ + +#include +#include + +#include "ReconstructionDataFormats/PID.h" + +namespace o2 +{ +namespace pwglf +{ + +/// @brief Extended the indices of the PID class with additional particles +class PIDExtended +{ + public: + typedef o2::track::PID::ID ID; + + static constexpr ID Electron = 0; + static constexpr ID Muon = 1; + static constexpr ID Pion = 2; + static constexpr ID Kaon = 3; + static constexpr ID Proton = 4; + static constexpr ID Deuteron = 5; + static constexpr ID Triton = 6; + static constexpr ID Helium3 = 7; + static constexpr ID Alpha = 8; + static constexpr ID PI0 = 9; + static constexpr ID Photon = 10; + static constexpr ID K0 = 11; + static constexpr ID Lambda = 12; + static constexpr ID HyperTriton = 13; + static constexpr ID Hyperhydrog4 = 14; + static constexpr ID XiMinus = 15; + static constexpr ID OmegaMinus = 16; + + static_assert(Electron == o2::track::PID::Electron, "PID::Electron mismatch"); + static_assert(Muon == o2::track::PID::Muon, "PID::Muon mismatch"); + static_assert(Pion == o2::track::PID::Pion, "PID::Pion mismatch"); + static_assert(Kaon == o2::track::PID::Kaon, "PID::Kaon mismatch"); + static_assert(Proton == o2::track::PID::Proton, "PID::Proton mismatch"); + static_assert(Deuteron == o2::track::PID::Deuteron, "PID::Deuteron mismatch"); + static_assert(Triton == o2::track::PID::Triton, "PID::Triton mismatch"); + static_assert(Helium3 == o2::track::PID::Helium3, "PID::Helium3 mismatch"); + static_assert(Alpha == o2::track::PID::Alpha, "PID::Alpha mismatch"); + static_assert(PI0 == o2::track::PID::PI0, "PID::PI0 mismatch"); + static_assert(Photon == o2::track::PID::Photon, "PID::Photon mismatch"); + static_assert(K0 == o2::track::PID::K0, "PID::K0 mismatch"); + static_assert(Lambda == o2::track::PID::Lambda, "PID::Lambda mismatch"); + static_assert(HyperTriton == o2::track::PID::HyperTriton, "PID::HyperTriton mismatch"); + static_assert(Hyperhydrog4 == o2::track::PID::Hyperhydrog4, "PID::Hyperhydrog4 mismatch"); + static_assert(XiMinus == o2::track::PID::XiMinus, "PID::XiMinus mismatch"); + static_assert(OmegaMinus == o2::track::PID::OmegaMinus, "PID::OmegaMinus mismatch"); + + static constexpr ID PIDCounts = 17; // Number of indices defined in PID.h + static_assert(PIDCounts == o2::track::PID::NIDsTot, "PID::NIDsTot mismatch"); + + // Define the antiparticles + static constexpr ID Positron = 17; + static constexpr ID MuonPlus = 18; + static constexpr ID PionMinus = 19; + static constexpr ID KaonMinus = 20; + static constexpr ID AntiProton = 21; + static constexpr ID AntiDeuteron = 22; + static constexpr ID AntiTriton = 23; + static constexpr ID AntiHelium3 = 24; + static constexpr ID AntiAlpha = 25; + static constexpr ID AntiLambda = 26; + static constexpr ID AntiHyperTriton = 27; + static constexpr ID AntiHyperhydrog4 = 28; + static constexpr ID XiPlus = 29; + static constexpr ID OmegaPlus = 30; + + static constexpr ID Neutron = 31; + static constexpr ID AntiNeutron = 32; + static constexpr ID HyperHelium4 = 33; + static constexpr ID AntiHyperHelium4 = 34; + static constexpr ID Phi = 35; + + static constexpr ID BZero = 36; + static constexpr ID BPlus = 37; + static constexpr ID BS = 38; + static constexpr ID D0 = 39; + static constexpr ID DPlus = 40; + static constexpr ID DS = 41; + static constexpr ID DStar = 42; + static constexpr ID ChiC1 = 43; + static constexpr ID JPsi = 44; + static constexpr ID LambdaB0 = 45; + static constexpr ID LambdaCPlus = 46; + static constexpr ID OmegaC0 = 47; + static constexpr ID SigmaC0 = 48; + static constexpr ID SigmaCPlusPlus = 49; + static constexpr ID X3872 = 50; + static constexpr ID Xi0 = 51; + static constexpr ID XiB0 = 52; + static constexpr ID XiCCPlusPlus = 53; + static constexpr ID XiCPlus = 54; + static constexpr ID XiC0 = 55; + static constexpr ID NIDsTot = 56; + + static constexpr const char* sNames[NIDsTot + 1] = { + o2::track::pid_constants::sNames[Electron], // Electron + o2::track::pid_constants::sNames[Muon], // Muon + o2::track::pid_constants::sNames[Pion], // Pion + o2::track::pid_constants::sNames[Kaon], // Kaon + o2::track::pid_constants::sNames[Proton], // Proton + o2::track::pid_constants::sNames[Deuteron], // Deuteron + o2::track::pid_constants::sNames[Triton], // Triton + o2::track::pid_constants::sNames[Helium3], // Helium3 + o2::track::pid_constants::sNames[Alpha], // Alpha + o2::track::pid_constants::sNames[PI0], // PI0 + o2::track::pid_constants::sNames[Photon], // Photon + o2::track::pid_constants::sNames[K0], // K0 + o2::track::pid_constants::sNames[Lambda], // Lambda + o2::track::pid_constants::sNames[HyperTriton], // HyperTriton + o2::track::pid_constants::sNames[Hyperhydrog4], // Hyperhydrog4 + o2::track::pid_constants::sNames[XiMinus], // XiMinus + o2::track::pid_constants::sNames[OmegaMinus], // OmegaMinus + "Positron", // Positron + "MuonPlus", // MuonPlus + "PionMinus", // PionMinus + "KaonMinus", // KaonMinus + "AntiProton", // AntiProton + "AntiDeuteron", // AntiDeuteron + "AntiTriton", // AntiTriton + "AntiHelium3", // AntiHelium3 + "AntiAlpha", // AntiAlpha + "AntiLambda", // AntiLambda + "AntiHyperTriton", // AntiHyperTriton + "AntiHyperhydrog4", // AntiHyperhydrog4 + "XiPlus", // XiPlus + "OmegaPlus", // OmegaPlus + "Neutron", // Neutron + "AntiNeutron", // AntiNeutron + "HyperHelium4", // HyperHelium4 + "AntiHyperHelium4", // AntiHyperHelium4 + "Phi", // Phi + "BZero", // BZero + "BPlus", // BPlus + "BS", // BS + "D0", // D0 + "DPlus", // DPlus + "DS", // DS + "DStar", // DStar + "ChiC1", // ChiC1 + "JPsi", // JPsi + "LambdaB0", // LambdaB0 + "LambdaCPlus", // LambdaCPlus + "OmegaC0", // OmegaC0 + "SigmaC0", // SigmaC0 + "SigmaCPlusPlus", // SigmaCPlusPlus + "X3872", // X3872 + "Xi0", // Xi0 + "XiB0", // XiB0 + "XiCCPlusPlus", // XiCCPlusPlus + "XiCPlus", // XiCPlus + "XiC0", // XiC0 + nullptr}; + + static std::vector arrayNames() + { + std::vector names; + for (int i = 0; i < NIDsTot; i++) { + names.push_back(sNames[i]); + } + return names; + } + + static const char* getName(ID id) { return sNames[id]; } + + /// \brief Convert PDG code to PID index + template + static o2::track::pid_constants::ID pdgToId(const TrackType& particle) + { + switch (particle.pdgCode()) { + case 11: + return Electron; + case -11: + return Positron; + case 13: + return Muon; + case -13: + return MuonPlus; + case 211: + return Pion; + case -211: + return PionMinus; + case 321: + return Kaon; + case -321: + return KaonMinus; + case 2212: + return Proton; + case -2212: + return AntiProton; + case 2112: + return Neutron; + case -2112: + return AntiNeutron; + case o2::constants::physics::Pdg::kDeuteron: + return Deuteron; + case -o2::constants::physics::Pdg::kDeuteron: + return AntiDeuteron; + case o2::constants::physics::Pdg::kTriton: + return Triton; + case -o2::constants::physics::Pdg::kTriton: + return AntiTriton; + case o2::constants::physics::Pdg::kHelium3: + return Helium3; + case -o2::constants::physics::Pdg::kHelium3: + return AntiHelium3; + case o2::constants::physics::Pdg::kAlpha: + return Alpha; + case -o2::constants::physics::Pdg::kAlpha: + return AntiAlpha; + case o2::constants::physics::Pdg::kHyperHelium4: + return HyperHelium4; + case -o2::constants::physics::Pdg::kHyperHelium4: + return AntiHyperHelium4; + case 111: + return PI0; + case 22: + return Photon; + case 130: + return K0; + case 3122: + return Lambda; + case -3122: + return AntiLambda; + case o2::constants::physics::Pdg::kHyperTriton: + return HyperTriton; + case -o2::constants::physics::Pdg::kHyperTriton: + return AntiHyperTriton; + case o2::constants::physics::Pdg::kHyperHydrogen4: + return Hyperhydrog4; + case -o2::constants::physics::Pdg::kHyperHydrogen4: + return AntiHyperhydrog4; + case 3312: + return XiMinus; + case -3312: + return XiPlus; + case 3334: + return OmegaMinus; + case -3334: + return OmegaPlus; + case o2::constants::physics::Pdg::kPhi: + return Phi; + case o2::constants::physics::Pdg::kB0: + return BZero; + case o2::constants::physics::Pdg::kBPlus: + return BPlus; + case o2::constants::physics::Pdg::kBS: + return BS; + case o2::constants::physics::Pdg::kD0: + return D0; + case o2::constants::physics::Pdg::kDPlus: + return DPlus; + case o2::constants::physics::Pdg::kDS: + return DS; + case o2::constants::physics::Pdg::kDStar: + return DStar; + case o2::constants::physics::Pdg::kChiC1: + return ChiC1; + case o2::constants::physics::Pdg::kJPsi: + return JPsi; + case o2::constants::physics::Pdg::kLambdaB0: + return LambdaB0; + case o2::constants::physics::Pdg::kLambdaCPlus: + return LambdaCPlus; + case o2::constants::physics::Pdg::kOmegaC0: + return OmegaC0; + case o2::constants::physics::Pdg::kSigmaC0: + return SigmaC0; + case o2::constants::physics::Pdg::kSigmaCPlusPlus: + return SigmaCPlusPlus; + case o2::constants::physics::Pdg::kX3872: + return X3872; + case o2::constants::physics::Pdg::kXi0: + return Xi0; + case o2::constants::physics::Pdg::kXiB0: + return XiB0; + case o2::constants::physics::Pdg::kXiCCPlusPlus: + return XiCCPlusPlus; + case o2::constants::physics::Pdg::kXiCPlus: + return XiCPlus; + case o2::constants::physics::Pdg::kXiC0: + return XiC0; + default: + LOG(info) << "Cannot identify particle with PDG code " << particle.pdgCode(); + break; + } + return -1; + } +}; + +} // namespace pwglf +} // namespace o2 + +#endif // PWGLF_UTILS_MCPARTICLE_H_