From 1d86b98785979e09a322ef19155f90217f011ec7 Mon Sep 17 00:00:00 2001 From: Nicolas Strangmann Date: Fri, 28 Feb 2025 17:49:43 +0100 Subject: [PATCH 1/7] Add first draft of heavy neutral meson software trigger --- EventFiltering/CMakeLists.txt | 5 + .../PWGEM/HeavyNeutralMesonFilter.cxx | 265 ++++++++++++++++++ EventFiltering/cefpTask.cxx | 1 + EventFiltering/filterTables.h | 65 +++-- 4 files changed, 310 insertions(+), 26 deletions(-) create mode 100644 EventFiltering/PWGEM/HeavyNeutralMesonFilter.cxx diff --git a/EventFiltering/CMakeLists.txt b/EventFiltering/CMakeLists.txt index 31cd58f40ed..ec9b151bc66 100644 --- a/EventFiltering/CMakeLists.txt +++ b/EventFiltering/CMakeLists.txt @@ -102,6 +102,11 @@ o2physics_add_dpl_workflow(em-photon-filter-qc PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2::DetectorsBase O2Physics::PWGEMPhotonMesonCore COMPONENT_NAME Analysis) +o2physics_add_dpl_workflow(heavy-neutral-meson-filter + SOURCES PWGEM/HeavyNeutralMesonFilter.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2::EMCALBase O2::EMCALCalib O2Physics::AnalysisCore + COMPONENT_NAME Analysis) + o2physics_add_dpl_workflow(lf-f1proton-filter SOURCES PWGLF/filterf1proton.cxx PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2::DetectorsBase diff --git a/EventFiltering/PWGEM/HeavyNeutralMesonFilter.cxx b/EventFiltering/PWGEM/HeavyNeutralMesonFilter.cxx new file mode 100644 index 00000000000..e5d785c40bb --- /dev/null +++ b/EventFiltering/PWGEM/HeavyNeutralMesonFilter.cxx @@ -0,0 +1,265 @@ +// 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 HeavyNeutralMesonFilter.cxx +/// +/// \brief This code loops over collisions to filter events contaning heavy mesons (omega or eta') using EMCal clusters +/// +/// \author Nicolas Strangmann (nicolas.strangmann@cern.ch) - Goethe University Frankfurt +/// + +#include +#include +#include +#include "TVector3.h" + +#include "Framework/runDataProcessing.h" +#include "Framework/AnalysisTask.h" +#include "Framework/HistogramRegistry.h" +#include "Common/DataModel/EventSelection.h" +#include "EMCALBase/Geometry.h" +#include "PWGJE/DataModel/EMCALClusters.h" +#include "Common/Core/trackUtilities.h" +#include "Common/DataModel/TrackSelectionTables.h" +#include "EventFiltering/filterTables.h" + +using namespace o2; +using namespace o2::framework; +using namespace o2::framework::expressions; + +using MyBCs = soa::Join; +using MyCollisions = soa::Join; + +using SelectedTracks = soa::Filtered>; +using SelectedClusters = soa::Filtered; // Clusters from collisions with only one collision in the BC + +struct Photon { // Struct to store photons (unique and ambiguous clusters that passed the cuts) + Photon(float eta, float phi, float energy) : eta(eta), phi(phi), e(energy), theta(2 * std::atan2(std::exp(-eta), 1)), px(e * std::sin(theta) * std::cos(phi)), py(e * std::sin(theta) * std::sin(phi)), pz(e * std::cos(theta)), pt(std::sqrt(px * px + py * py)) + { + photon.SetPxPyPzE(px, py, pz, e); + } + + ROOT::Math::PxPyPzEVector photon; + float eta, phi, e; + float theta; + float px, py, pz, pt; +}; + +namespace HeavyNeutralMeson +{ +enum MesonType { + kPi0 = 1 << 0, // 0001 + kEta = 1 << 1, // 0010 + kOmega = 1 << 2, // 0100 + kEtaPrime = 1 << 3 // 1000 +}; +} + +struct Meson { + Meson(Photon p1, Photon p2) : p1(p1), p2(p2), mesonType(0) + { + pMeson = p1.photon + p2.photon; + } + Meson(Meson gg, float eTracks, float pxTracks, float pyTracks, float pzTracks) : p1(gg.p1), p2(gg.p2), mesonType(0) + { + pMeson.SetPxPyPzE(gg.pMeson.Px() + pxTracks, gg.pMeson.Py() + pyTracks, gg.pMeson.Pz() + pzTracks, gg.pMeson.e() + eTracks); + } + Photon p1, p2; + ROOT::Math::PxPyPzEVector pMeson; + + uint32_t mesonType; + + void addMesonType(uint32_t type) { mesonType |= type; } + bool hasMesonType(uint32_t type) const { return mesonType & type; } + + float m() const { return pMeson.M(); } + float pT() const { return pMeson.Pt(); } +}; + +struct HeavyNeutralMesonFilter { + + Produces tags; + + HistogramRegistry mHistManager{"HeavyNeutralMesonFilterHistograms", {}, OutputObjHandlingPolicy::QAObject}; + + Configurable cfgClusterDefinition{"cfgClusterDefinition", 13, "Clusterizer to be selected, e.g. 13 for kV3MostSplitLowSeed"}; + Configurable cfgMinClusterEnergy{"cfgMinClusterEnergy", 0.5, "Minimum energy of selected clusters (GeV)"}; + Configurable cfgMinM02{"cfgMinM02", 0., "Minimum M02 of selected clusters"}; + Configurable cfgMaxM02{"cfgMaxM02", 1, "Maximum M02 of selected clusters"}; + Configurable cfgMinTime{"cfgMinTime", -25, "Minimum time of selected clusters (ns)"}; + Configurable cfgMaxTime{"cfgMaxTime", 25, "Maximum time of selected clusters (ns)"}; + + Configurable cfgTrackMinPt{"cfgTrackMinPt", 0.1, "Minimum momentum of tracks (GeV/c)"}; + + static constexpr float defaultMassWindows[2][4] = {{0., 0.4, 0.6, 1.}, {0.4, 0.8, 0.8, 1.2}}; + Configurable> massWindowOmega{"massWindowOmega", {defaultMassWindows[0], 4, {"pi0_min", "pi0_max", "omega_min", "omega_max"}}, "Mass window for selected omegas and their decay pi0"}; + Configurable> massWindowEtaPrime{"massWindowEtaPrime", {defaultMassWindows[1], 4, {"eta_min", "eta_max", "etaprime_min", "etaprime_max"}}, "Mass window for selected eta' and their decay eta"}; + + static constexpr float defaultMinPts[4] = {0., 5., 0., 5.}; // LowPtOmega, HighPtOmega, LowPtEtaPrime, HighPtEtaPrime + Configurable> minHNMPt{"minHNMPt", {defaultMinPts, 4, {"omega_lowPt", "omega_highPt", "etaPrime_lowPt", "etaPrime_highPt"}}, "Minimum pT values for the heavy neutral mesons of the low and high pT trigger"}; + + Filter clusterDefinitionFilter = aod::emcalcluster::definition == cfgClusterDefinition; + Filter energyFilter = aod::emcalcluster::energy > cfgMinClusterEnergy; + Filter m02Filter = (aod::emcalcluster::nCells == 1 || (aod::emcalcluster::m02 > cfgMinM02 && aod::emcalcluster::m02 < cfgMaxM02)); + Filter timeFilter = (aod::emcalcluster::time > cfgMinTime && aod::emcalcluster::time < cfgMaxTime); + + Filter trackPtFilter = aod::track::pt > cfgTrackMinPt; + + std::vector mPhotons; + std::vector mGGs; + + bool colContainsLowPtOmega = false; + bool colContainsHighPtOmega = false; + bool colContainsLowPtEtaPrime = false; + bool colContainsHighPtEtaPrime = false; + + emcal::Geometry* emcalGeom; + + void init(InitContext const&) + { + emcalGeom = emcal::Geometry::GetInstanceFromRunNumber(300000); + auto hCollisionCounter = mHistManager.add("Event/hCollisionCounter", "Number of collisions;;#bf{#it{N}_{Coll}}", HistType::kTH1F, {{8, -0.5, 7.5}}); + hCollisionCounter->GetXaxis()->SetBinLabel(1, "all"); + hCollisionCounter->GetXaxis()->SetBinLabel(2, "kTVXinEMC"); + hCollisionCounter->GetXaxis()->SetBinLabel(3, "L0Triggered"); + hCollisionCounter->GetXaxis()->SetBinLabel(4, "2+ tracks&clusters"); + hCollisionCounter->GetXaxis()->SetBinLabel(5, "Low #it{p}_{T} #omega"); + hCollisionCounter->GetXaxis()->SetBinLabel(6, "High #it{p}_{T} #omega"); + hCollisionCounter->GetXaxis()->SetBinLabel(7, "Low #it{p}_{T} #eta'"); + hCollisionCounter->GetXaxis()->SetBinLabel(8, "High #it{p}_{T} #eta'"); + + mHistManager.add("Event/nClusters", "Number of BCs (with (k)TVX);#bf{TVX triggered};#bf{#it{N}_{BCs}}", HistType::kTH1F, {{51, -0.5, 50.5}}); + mHistManager.add("Event/nGGs", "Number of BCs (with (k)TVX);#bf{TVX triggered};#bf{#it{N}_{BCs}}", HistType::kTH1F, {{51, -0.5, 50.5}}); + mHistManager.add("Event/nTracks", "Number of BCs (with (k)TVX);#bf{TVX triggered};#bf{#it{N}_{BCs}}", HistType::kTH1F, {{51, -0.5, 50.5}}); + mHistManager.add("Event/nHeavyNeutralMesons", "Number of BCs (with (k)TVX);#bf{TVX triggered};#bf{#it{N}_{BCs}}", HistType::kTH1F, {{51, -0.5, 50.5}}); + + mHistManager.add("nCollisionsVsClusters", "Number of collisions vs number of clusters;N_{Collisions};N_{Clusters}", HistType::kTH2F, {{4, -0.5, 3.5}, {21, -0.5, 20.5}}); + + mHistManager.add("Track/trackP", "Momentum of tracks;#bf{#it{p} (GeV/#it{c})};#bf{#it{N}_{tracks}}", HistType::kTH1F, {{200, 0, 20}}); + + mHistManager.add("Cluster/clusterE", "Energy of cluster;#bf{#it{E} (GeV)};#bf{#it{N}_{clusters}}", HistType::kTH1F, {{200, 0, 20}}); + mHistManager.add("Cluster/clusterM02", "Shape of cluster;#bf{#it{M}_{02}};#bf{#it{N}_{clusters}}", HistType::kTH1F, {{200, 0, 2}}); + mHistManager.add("Cluster/clusterTime", "Time of cluster;#bf{#it{t} (ns)};#bf{#it{N}_{clusters}}", HistType::kTH1F, {{200, -100, 100}}); + + mHistManager.add("GG/invMassVsPt", "Invariant mass and pT of meson candidates", HistType::kTH2F, {{400, massWindowOmega->get("pi0_min"), massWindowEtaPrime->get("eta_max")}, {250, 0., 25.}}); + + mHistManager.add("HeavyNeutralMeson/invMassVsPt", "Invariant mass and pT of meson candidates", HistType::kTH2F, {{600, massWindowOmega->get("omega_min"), massWindowEtaPrime->get("etaprime_max")}, {250, 0., 25.}}); + } + + void process(MyCollisions::iterator const& collision, MyBCs const&, SelectedClusters const& clusters, SelectedTracks const& tracks) + { + mHistManager.fill(HIST("Event/hCollisionCounter"), 0.); + + colContainsLowPtOmega = false; + colContainsHighPtOmega = false; + colContainsLowPtEtaPrime = false; + colContainsHighPtEtaPrime = false; + + if (isBCSelected(collision.foundBC_as(), clusters.size(), tracks.size() > 2)) { + processClusters(clusters); + + processGGs(); + processHeavyNeutralMesons(tracks); + + mPhotons.clear(); + mGGs.clear(); + } + + if (colContainsLowPtOmega) + mHistManager.fill(HIST("Event/hCollisionCounter"), 4.); + if (colContainsHighPtOmega) + mHistManager.fill(HIST("Event/hCollisionCounter"), 5.); + if (colContainsLowPtEtaPrime) + mHistManager.fill(HIST("Event/hCollisionCounter"), 6.); + if (colContainsHighPtEtaPrime) + mHistManager.fill(HIST("Event/hCollisionCounter"), 7.); + tags(colContainsLowPtOmega, colContainsHighPtOmega, colContainsLowPtEtaPrime, colContainsHighPtEtaPrime); + } + + bool isBCSelected(MyBCs::iterator const& bc, int nClusters, bool enoughTracks) + { + if (bc.alias_bit(kTVXinEMC)) + mHistManager.fill(HIST("Event/hCollisionCounter"), 1.); + else if (bc.alias_bit(kEMC7) || bc.alias_bit(kEG1) || bc.alias_bit(kEG2) || bc.alias_bit(kDG1) || bc.alias_bit(kDG2) || bc.alias_bit(kEJ1) || bc.alias_bit(kEJ2) || bc.alias_bit(kDJ1) || bc.alias_bit(kDJ2)) + mHistManager.fill(HIST("Event/hCollisionCounter"), 2.); + else + return false; + if (nClusters < 2 || !enoughTracks) + return false; + mHistManager.fill(HIST("Event/hCollisionCounter"), 3.); + + return true; + } + + void processClusters(SelectedClusters const& clusters) + { + for (const auto& cluster : clusters) { + mHistManager.fill(HIST("Cluster/clusterE"), cluster.energy()); + mHistManager.fill(HIST("Cluster/clusterM02"), cluster.m02()); + mHistManager.fill(HIST("Cluster/clusterTime"), cluster.time()); + mPhotons.push_back(Photon(cluster.eta(), cluster.phi(), cluster.energy())); + } + mHistManager.fill(HIST("Event/nClusters"), clusters.size()); + } + + /// \brief Process light neutral meson candidates (pi0 or eta), calculate invariant mass and pT and fill histograms + void processGGs() + { + // loop over all photon combinations and build meson candidates + for (unsigned int ig1 = 0; ig1 < mPhotons.size(); ++ig1) { + for (unsigned int ig2 = ig1 + 1; ig2 < mPhotons.size(); ++ig2) { + Meson lightMeson(mPhotons[ig1], mPhotons[ig2]); // build lightMeson from photons + mHistManager.fill(HIST("GG/invMassVsPt"), lightMeson.m(), lightMeson.pT()); + if (lightMeson.m() > massWindowOmega->get("pi0_min") && lightMeson.m() < massWindowOmega->get("pi0_max")) + lightMeson.addMesonType(HeavyNeutralMeson::kPi0); + else if (lightMeson.m() > massWindowEtaPrime->get("eta_min") && lightMeson.m() < massWindowEtaPrime->get("eta_max")) + lightMeson.addMesonType(HeavyNeutralMeson::kEta); + else + continue; + mGGs.push_back(lightMeson); + } + } + mHistManager.fill(HIST("Event/nGGs"), mGGs.size()); + } + + void processHeavyNeutralMesons(SelectedTracks const& tracks) + { + for (const auto& posTrack : tracks) { + if (!posTrack.isGlobalTrack() || posTrack.sign() < 0) + continue; + for (const auto& negTrack : tracks) { + if (!negTrack.isGlobalTrack() || negTrack.sign() > 0) + continue; + for (unsigned long iGG = 0; iGG < mGGs.size(); iGG++) { + Meson heavyNeutralMeson(mGGs.at(iGG), posTrack.energy(constants::physics::MassPiPlus) + negTrack.energy(constants::physics::MassPiPlus), posTrack.px() + negTrack.px(), posTrack.py() + negTrack.py(), posTrack.pz() + negTrack.pz()); + mHistManager.fill(HIST("HeavyNeutralMeson/invMassVsPt"), heavyNeutralMeson.m(), heavyNeutralMeson.pT()); + if (mGGs.at(iGG).hasMesonType(HeavyNeutralMeson::kPi0) && heavyNeutralMeson.m() > massWindowOmega->get("omega_min") && heavyNeutralMeson.m() < massWindowOmega->get("omega_max")) { + if(minHNMPt->get("omega_lowPt") < heavyNeutralMeson.pT()) + colContainsLowPtOmega = true; + if(minHNMPt->get("omega_highPt") < heavyNeutralMeson.pT()) + colContainsHighPtOmega = true; + } + if (mGGs.at(iGG).hasMesonType(HeavyNeutralMeson::kEta) && heavyNeutralMeson.m() > massWindowEtaPrime->get("etaprime_min") && heavyNeutralMeson.m() < massWindowEtaPrime->get("etaprime_max")) { + if(minHNMPt->get("etaPrime_lowPt") < heavyNeutralMeson.pT()) + colContainsLowPtEtaPrime = true; + if(minHNMPt->get("etaPrime_highPt") < heavyNeutralMeson.pT()) + colContainsHighPtEtaPrime = true; + } + } + } + } + } +}; + +WorkflowSpec defineDataProcessing(o2::framework::ConfigContext const& cfgc) +{ + return WorkflowSpec{adaptAnalysisTask(cfgc)}; +} diff --git a/EventFiltering/cefpTask.cxx b/EventFiltering/cefpTask.cxx index b5cfa90449a..56c75c2a047 100644 --- a/EventFiltering/cefpTask.cxx +++ b/EventFiltering/cefpTask.cxx @@ -227,6 +227,7 @@ struct centralEventFilterTask { FILTER_CONFIGURABLE(MultFilters); FILTER_CONFIGURABLE(FullJetFilters); FILTER_CONFIGURABLE(PhotonFilters); + FILTER_CONFIGURABLE(HeavyNeutralMesonFilters); void init(o2::framework::InitContext& initc) { diff --git a/EventFiltering/filterTables.h b/EventFiltering/filterTables.h index 2155e6528eb..32574ed65f9 100644 --- a/EventFiltering/filterTables.h +++ b/EventFiltering/filterTables.h @@ -46,9 +46,9 @@ namespace o2::aod { namespace filtering { -DECLARE_SOA_COLUMN(H2, hasH2, bool); //! deuteron trigger for the helium normalisation (to be downscaled) -DECLARE_SOA_COLUMN(He, hasHe, bool); //! helium -DECLARE_SOA_COLUMN(H3L3Body, hasH3L3Body, bool); //! hypertriton 3body +DECLARE_SOA_COLUMN(H2, hasH2, bool); //! deuteron trigger for the helium normalisation (to be downscaled) +DECLARE_SOA_COLUMN(He, hasHe, bool); //! helium +DECLARE_SOA_COLUMN(H3L3Body, hasH3L3Body, bool); //! hypertriton 3body DECLARE_SOA_COLUMN(ITSextremeIonisation, hasITSextremeIonisation, bool); //! ITS extreme ionisation DECLARE_SOA_COLUMN(ITSmildIonisation, hasITSmildIonisation, bool); //! ITS mild ionisation (normalisation of the extreme ionisation), to be downscaled @@ -71,25 +71,25 @@ DECLARE_SOA_COLUMN(LMeeHMR, hasLMeeHMR, bool); //! dielectron trigger for high m DECLARE_SOA_COLUMN(ElectronMuon, hasElectronMuon, bool); //! dimuon trigger with low pT on muons // heavy flavours -DECLARE_SOA_COLUMN(HfHighPt2P, hasHfHighPt2P, bool); //! high-pT 2-prong charm hadron -DECLARE_SOA_COLUMN(HfHighPt3P, hasHfHighPt3P, bool); //! high-pT 3-prong charm hadron -DECLARE_SOA_COLUMN(HfBeauty3P, hasHfBeauty3P, bool); //! 3-prong beauty hadron -DECLARE_SOA_COLUMN(HfBeauty4P, hasHfBeauty4P, bool); //! 4-prong beauty hadron -DECLARE_SOA_COLUMN(HfFemto2P, hasHfFemto2P, bool); //! 2-prong charm-hadron - N pair -DECLARE_SOA_COLUMN(HfFemto3P, hasHfFemto3P, bool); //! 3-prong charm-hadron - N pair -DECLARE_SOA_COLUMN(HfDoubleCharm2P, hasHfDoubleCharm2P, bool); //! at least two 2-prong charm-hadron candidates -DECLARE_SOA_COLUMN(HfDoubleCharm3P, hasHfDoubleCharm3P, bool); //! at least two 3-prong charm-hadron candidates -DECLARE_SOA_COLUMN(HfDoubleCharmMix, hasHfDoubleCharmMix, bool); //! at least one 2-prong and one 3-prong charm-hadron candidates -DECLARE_SOA_COLUMN(HfV0Charm2P, hasHfV0Charm2P, bool); //! V0 with 2-prong charm hadron -DECLARE_SOA_COLUMN(HfV0Charm3P, hasHfV0Charm3P, bool); //! V0 with 3-prong charm hadron -DECLARE_SOA_COLUMN(HfCharmBarToXiBach, hasHfCharmBarToXiBach, bool); //! Charm baryon to Xi + bachelor +DECLARE_SOA_COLUMN(HfHighPt2P, hasHfHighPt2P, bool); //! high-pT 2-prong charm hadron +DECLARE_SOA_COLUMN(HfHighPt3P, hasHfHighPt3P, bool); //! high-pT 3-prong charm hadron +DECLARE_SOA_COLUMN(HfBeauty3P, hasHfBeauty3P, bool); //! 3-prong beauty hadron +DECLARE_SOA_COLUMN(HfBeauty4P, hasHfBeauty4P, bool); //! 4-prong beauty hadron +DECLARE_SOA_COLUMN(HfFemto2P, hasHfFemto2P, bool); //! 2-prong charm-hadron - N pair +DECLARE_SOA_COLUMN(HfFemto3P, hasHfFemto3P, bool); //! 3-prong charm-hadron - N pair +DECLARE_SOA_COLUMN(HfDoubleCharm2P, hasHfDoubleCharm2P, bool); //! at least two 2-prong charm-hadron candidates +DECLARE_SOA_COLUMN(HfDoubleCharm3P, hasHfDoubleCharm3P, bool); //! at least two 3-prong charm-hadron candidates +DECLARE_SOA_COLUMN(HfDoubleCharmMix, hasHfDoubleCharmMix, bool); //! at least one 2-prong and one 3-prong charm-hadron candidates +DECLARE_SOA_COLUMN(HfV0Charm2P, hasHfV0Charm2P, bool); //! V0 with 2-prong charm hadron +DECLARE_SOA_COLUMN(HfV0Charm3P, hasHfV0Charm3P, bool); //! V0 with 3-prong charm hadron +DECLARE_SOA_COLUMN(HfCharmBarToXiBach, hasHfCharmBarToXiBach, bool); //! Charm baryon to Xi + bachelor DECLARE_SOA_COLUMN(HfCharmBarToXi2Bach, hasHfCharmBarToXi2Bach, bool); //! Charm baryon to Xi + 2 bachelors -DECLARE_SOA_COLUMN(HfSigmaCPPK, hasHfSigmaCPPK, bool); //! SigmaC(2455)++K- and SigmaC(2520)++K- + c.c. -DECLARE_SOA_COLUMN(HfSigmaC0K0, hasHfSigmaC0K0, bool); //! SigmaC(2455)0KS0 and SigmaC(2520)0KS0 -DECLARE_SOA_COLUMN(HfPhotonCharm2P, hasHfPhotonCharm2P, bool); //! photon with 2-prong charm hadron -DECLARE_SOA_COLUMN(HfPhotonCharm3P, hasHfPhotonCharm3P, bool); //! photon with 3-prong charm hadron -DECLARE_SOA_COLUMN(HfSingleCharm2P, hasHfSingleCharm2P, bool); //! 2-prong charm hadron (for efficiency studies) -DECLARE_SOA_COLUMN(HfSingleCharm3P, hasHfSingleCharm3P, bool); //! 3-prong charm hadron (for efficiency studies) +DECLARE_SOA_COLUMN(HfSigmaCPPK, hasHfSigmaCPPK, bool); //! SigmaC(2455)++K- and SigmaC(2520)++K- + c.c. +DECLARE_SOA_COLUMN(HfSigmaC0K0, hasHfSigmaC0K0, bool); //! SigmaC(2455)0KS0 and SigmaC(2520)0KS0 +DECLARE_SOA_COLUMN(HfPhotonCharm2P, hasHfPhotonCharm2P, bool); //! photon with 2-prong charm hadron +DECLARE_SOA_COLUMN(HfPhotonCharm3P, hasHfPhotonCharm3P, bool); //! photon with 3-prong charm hadron +DECLARE_SOA_COLUMN(HfSingleCharm2P, hasHfSingleCharm2P, bool); //! 2-prong charm hadron (for efficiency studies) +DECLARE_SOA_COLUMN(HfSingleCharm3P, hasHfSingleCharm3P, bool); //! 3-prong charm hadron (for efficiency studies) DECLARE_SOA_COLUMN(HfSingleNonPromptCharm2P, hasHfSingleNonPromptCharm2P, bool); //! 2-prong charm hadron (for efficiency studies) DECLARE_SOA_COLUMN(HfSingleNonPromptCharm3P, hasHfSingleNonPromptCharm3P, bool); //! 3-prong charm hadron (for efficiency studies) DECLARE_SOA_COLUMN(HfBtoJPsiKa, hasHfBtoJPsiKa, bool); //! B+ -> JPsi(->mumu)K+ @@ -174,6 +174,12 @@ DECLARE_SOA_COLUMN(PCMHighPtPhoton, hasPCMHighPtPhoton, bool); //! PCM high pT p // DECLARE_SOA_COLUMN(PCMEtaDalitz, hasPCMEtaDalitz, bool); //! PCM eta -> ee gamma // DECLARE_SOA_COLUMN(PCMEtaGG, hasPCMEtaGG, bool); //! PCM eta -> ee gamma DECLARE_SOA_COLUMN(PCMandEE, hasPCMandEE, bool); //! PCM and ee + +// heavy meson filters +DECLARE_SOA_COLUMN(LowPtOmegaMeson, hasLowPtOmegaMeson, bool); //! Omega meson candidate (3pi) in the collision +DECLARE_SOA_COLUMN(HighPtOmegaMeson, hasHighPtOmegaMeson, bool); //! Omega meson candidate (3pi) in the collision +DECLARE_SOA_COLUMN(LowPtEtaPrimeMeson, hasLowPtEtaPrimeMeson, bool); //! Eta prime candidate (pi+pi-eta) in the collision +DECLARE_SOA_COLUMN(HighPtEtaPrimeMeson, hasHighPtEtaPrimeMeson, bool); //! Eta prime candidate (pi+pi-eta) in the collision } // namespace filtering namespace decision @@ -298,6 +304,13 @@ DECLARE_SOA_TABLE(PhotonFilters, "AOD", "PhotonFilters", //! using PhotonFilter = PhotonFilters::iterator; +// heavy mesons +DECLARE_SOA_TABLE(HeavyNeutralMesonFilters, "AOD", "HeavyNeutralMesonFilters", //! + filtering::LowPtOmegaMeson, filtering::HighPtOmegaMeson, + filtering::LowPtEtaPrimeMeson, filtering::HighPtEtaPrimeMeson); + +using HeavyNeutralMesonFilter = HeavyNeutralMesonFilters::iterator; + // cefp decision DECLARE_SOA_TABLE(CefpDecisions, "AOD", "CefpDecision", //! decision::BCId, decision::GlobalBCId, decision::EvSelBC, decision::CollisionTime, decision::CollisionTimeRes, decision::CefpTriggered0, decision::CefpTriggered1, decision::CefpSelected0, decision::CefpSelected1); @@ -309,11 +322,11 @@ DECLARE_SOA_TABLE(BCRanges, "AOD", "BCRanges", //! using BCRange = BCRanges::iterator; /// List of the available filters, the description of their tables and the name of the tasks -constexpr int NumberOfFilters{12}; -constexpr std::array AvailableFilters{"NucleiFilters", "DiffractionFilters", "DqFilters", "HfFilters", "CFFilters", "JetFilters", "JetHFFilters", "FullJetFilters", "StrangenessFilters", "MultFilters", "PhotonFilters", "F1ProtonFilters"}; -constexpr std::array FilterDescriptions{"NucleiFilters", "DiffFilters", "DqFilters", "HfFilters", "CFFilters", "JetFilters", "JetHFFilters", "FullJetFilters", "LFStrgFilters", "MultFilters", "PhotonFilters", "F1ProtonFilters"}; -constexpr std::array FilteringTaskNames{"o2-analysis-nuclei-filter", "o2-analysis-diffraction-filter", "o2-analysis-dq-filter-pp-with-association", "o2-analysis-hf-filter", "o2-analysis-cf-filter", "o2-analysis-je-filter", "o2-analysis-je-hf-filter", "o2-analysis-fje-filter", "o2-analysis-lf-strangeness-filter", "o2-analysis-mult-filter", "o2-analysis-em-photon-filter", "o2-analysis-lf-f1proton-filter"}; -constexpr o2::framework::pack FiltersPack; +constexpr int NumberOfFilters{13}; +constexpr std::array AvailableFilters{"NucleiFilters", "DiffractionFilters", "DqFilters", "HfFilters", "CFFilters", "JetFilters", "JetHFFilters", "FullJetFilters", "StrangenessFilters", "MultFilters", "PhotonFilters", "F1ProtonFilters", "HeavyNeutralMesonFilters"}; +constexpr std::array FilterDescriptions{"NucleiFilters", "DiffFilters", "DqFilters", "HfFilters", "CFFilters", "JetFilters", "JetHFFilters", "FullJetFilters", "LFStrgFilters", "MultFilters", "PhotonFilters", "F1ProtonFilters", "HNMesonFilters"}; +constexpr std::array FilteringTaskNames{"o2-analysis-nuclei-filter", "o2-analysis-diffraction-filter", "o2-analysis-dq-filter-pp-with-association", "o2-analysis-hf-filter", "o2-analysis-cf-filter", "o2-analysis-je-filter", "o2-analysis-je-hf-filter", "o2-analysis-fje-filter", "o2-analysis-lf-strangeness-filter", "o2-analysis-mult-filter", "o2-analysis-em-photon-filter", "o2-analysis-lf-f1proton-filter", "o2-analysis-heavy-meson-filter"}; +constexpr o2::framework::pack FiltersPack; static_assert(o2::framework::pack_size(FiltersPack) == NumberOfFilters); template From c7eb32feaf409b9cd50e98b4f944459bb1af76cc Mon Sep 17 00:00:00 2001 From: Nicolas Strangmann Date: Sun, 9 Mar 2025 20:24:31 +0100 Subject: [PATCH 2/7] Add PCM to omega trigger task --- EventFiltering/CMakeLists.txt | 2 +- .../PWGEM/HeavyNeutralMesonFilter.cxx | 281 ++++++------------ EventFiltering/filterTables.h | 12 +- .../TableProducer/skimmerGammaCalo.cxx | 58 +++- PWGEM/PhotonMeson/Utils/HNMUtilities.h | 175 +++++++++++ 5 files changed, 320 insertions(+), 208 deletions(-) create mode 100644 PWGEM/PhotonMeson/Utils/HNMUtilities.h diff --git a/EventFiltering/CMakeLists.txt b/EventFiltering/CMakeLists.txt index ec9b151bc66..af42d63b16c 100644 --- a/EventFiltering/CMakeLists.txt +++ b/EventFiltering/CMakeLists.txt @@ -104,7 +104,7 @@ o2physics_add_dpl_workflow(em-photon-filter-qc o2physics_add_dpl_workflow(heavy-neutral-meson-filter SOURCES PWGEM/HeavyNeutralMesonFilter.cxx - PUBLIC_LINK_LIBRARIES O2::Framework O2::EMCALBase O2::EMCALCalib O2Physics::AnalysisCore + PUBLIC_LINK_LIBRARIES O2::Framework O2::EMCALBase O2::EMCALCalib O2Physics::AnalysisCore O2Physics::PWGEMPhotonMesonCore COMPONENT_NAME Analysis) o2physics_add_dpl_workflow(lf-f1proton-filter diff --git a/EventFiltering/PWGEM/HeavyNeutralMesonFilter.cxx b/EventFiltering/PWGEM/HeavyNeutralMesonFilter.cxx index e5d785c40bb..b09fc89f014 100644 --- a/EventFiltering/PWGEM/HeavyNeutralMesonFilter.cxx +++ b/EventFiltering/PWGEM/HeavyNeutralMesonFilter.cxx @@ -16,110 +16,35 @@ /// \author Nicolas Strangmann (nicolas.strangmann@cern.ch) - Goethe University Frankfurt /// -#include -#include -#include -#include "TVector3.h" - -#include "Framework/runDataProcessing.h" -#include "Framework/AnalysisTask.h" -#include "Framework/HistogramRegistry.h" -#include "Common/DataModel/EventSelection.h" -#include "EMCALBase/Geometry.h" -#include "PWGJE/DataModel/EMCALClusters.h" -#include "Common/Core/trackUtilities.h" -#include "Common/DataModel/TrackSelectionTables.h" -#include "EventFiltering/filterTables.h" +#include "PWGEM/PhotonMeson/Utils/HNMUtilities.h" using namespace o2; using namespace o2::framework; using namespace o2::framework::expressions; +using namespace o2::aod::pwgem::photonmeson::utils::hnmutilities; using MyBCs = soa::Join; using MyCollisions = soa::Join; using SelectedTracks = soa::Filtered>; -using SelectedClusters = soa::Filtered; // Clusters from collisions with only one collision in the BC - -struct Photon { // Struct to store photons (unique and ambiguous clusters that passed the cuts) - Photon(float eta, float phi, float energy) : eta(eta), phi(phi), e(energy), theta(2 * std::atan2(std::exp(-eta), 1)), px(e * std::sin(theta) * std::cos(phi)), py(e * std::sin(theta) * std::sin(phi)), pz(e * std::cos(theta)), pt(std::sqrt(px * px + py * py)) - { - photon.SetPxPyPzE(px, py, pz, e); - } - - ROOT::Math::PxPyPzEVector photon; - float eta, phi, e; - float theta; - float px, py, pz, pt; -}; - -namespace HeavyNeutralMeson -{ -enum MesonType { - kPi0 = 1 << 0, // 0001 - kEta = 1 << 1, // 0010 - kOmega = 1 << 2, // 0100 - kEtaPrime = 1 << 3 // 1000 -}; -} - -struct Meson { - Meson(Photon p1, Photon p2) : p1(p1), p2(p2), mesonType(0) - { - pMeson = p1.photon + p2.photon; - } - Meson(Meson gg, float eTracks, float pxTracks, float pyTracks, float pzTracks) : p1(gg.p1), p2(gg.p2), mesonType(0) - { - pMeson.SetPxPyPzE(gg.pMeson.Px() + pxTracks, gg.pMeson.Py() + pyTracks, gg.pMeson.Pz() + pzTracks, gg.pMeson.e() + eTracks); - } - Photon p1, p2; - ROOT::Math::PxPyPzEVector pMeson; - - uint32_t mesonType; - - void addMesonType(uint32_t type) { mesonType |= type; } - bool hasMesonType(uint32_t type) const { return mesonType & type; } - - float m() const { return pMeson.M(); } - float pT() const { return pMeson.Pt(); } -}; struct HeavyNeutralMesonFilter { - Produces tags; HistogramRegistry mHistManager{"HeavyNeutralMesonFilterHistograms", {}, OutputObjHandlingPolicy::QAObject}; - Configurable cfgClusterDefinition{"cfgClusterDefinition", 13, "Clusterizer to be selected, e.g. 13 for kV3MostSplitLowSeed"}; - Configurable cfgMinClusterEnergy{"cfgMinClusterEnergy", 0.5, "Minimum energy of selected clusters (GeV)"}; - Configurable cfgMinM02{"cfgMinM02", 0., "Minimum M02 of selected clusters"}; - Configurable cfgMaxM02{"cfgMaxM02", 1, "Maximum M02 of selected clusters"}; - Configurable cfgMinTime{"cfgMinTime", -25, "Minimum time of selected clusters (ns)"}; - Configurable cfgMaxTime{"cfgMaxTime", 25, "Maximum time of selected clusters (ns)"}; - Configurable cfgTrackMinPt{"cfgTrackMinPt", 0.1, "Minimum momentum of tracks (GeV/c)"}; - + Configurable cfgHNMMassCorrection{"cfgHNMMassCorrection", 1, "Use GG PDG mass to correct HNM mass (0 = off, 1 = subDeltaPi0, 2 = subLambda)"}; static constexpr float defaultMassWindows[2][4] = {{0., 0.4, 0.6, 1.}, {0.4, 0.8, 0.8, 1.2}}; Configurable> massWindowOmega{"massWindowOmega", {defaultMassWindows[0], 4, {"pi0_min", "pi0_max", "omega_min", "omega_max"}}, "Mass window for selected omegas and their decay pi0"}; Configurable> massWindowEtaPrime{"massWindowEtaPrime", {defaultMassWindows[1], 4, {"eta_min", "eta_max", "etaprime_min", "etaprime_max"}}, "Mass window for selected eta' and their decay eta"}; - static constexpr float defaultMinPts[4] = {0., 5., 0., 5.}; // LowPtOmega, HighPtOmega, LowPtEtaPrime, HighPtEtaPrime - Configurable> minHNMPt{"minHNMPt", {defaultMinPts, 4, {"omega_lowPt", "omega_highPt", "etaPrime_lowPt", "etaPrime_highPt"}}, "Minimum pT values for the heavy neutral mesons of the low and high pT trigger"}; - - Filter clusterDefinitionFilter = aod::emcalcluster::definition == cfgClusterDefinition; - Filter energyFilter = aod::emcalcluster::energy > cfgMinClusterEnergy; - Filter m02Filter = (aod::emcalcluster::nCells == 1 || (aod::emcalcluster::m02 > cfgMinM02 && aod::emcalcluster::m02 < cfgMaxM02)); - Filter timeFilter = (aod::emcalcluster::time > cfgMinTime && aod::emcalcluster::time < cfgMaxTime); - Filter trackPtFilter = aod::track::pt > cfgTrackMinPt; - std::vector mPhotons; - std::vector mGGs; + std::vector vGGs; + std::vector vHNMs; - bool colContainsLowPtOmega = false; - bool colContainsHighPtOmega = false; - bool colContainsLowPtEtaPrime = false; - bool colContainsHighPtEtaPrime = false; + bool colContainsPCMOmega, colContainsEMCOmega, colContainsPCMEtaPrime, colContainsEMCEtaPrime = false; emcal::Geometry* emcalGeom; @@ -131,131 +56,117 @@ struct HeavyNeutralMesonFilter { hCollisionCounter->GetXaxis()->SetBinLabel(2, "kTVXinEMC"); hCollisionCounter->GetXaxis()->SetBinLabel(3, "L0Triggered"); hCollisionCounter->GetXaxis()->SetBinLabel(4, "2+ tracks&clusters"); - hCollisionCounter->GetXaxis()->SetBinLabel(5, "Low #it{p}_{T} #omega"); - hCollisionCounter->GetXaxis()->SetBinLabel(6, "High #it{p}_{T} #omega"); - hCollisionCounter->GetXaxis()->SetBinLabel(7, "Low #it{p}_{T} #eta'"); - hCollisionCounter->GetXaxis()->SetBinLabel(8, "High #it{p}_{T} #eta'"); - - mHistManager.add("Event/nClusters", "Number of BCs (with (k)TVX);#bf{TVX triggered};#bf{#it{N}_{BCs}}", HistType::kTH1F, {{51, -0.5, 50.5}}); - mHistManager.add("Event/nGGs", "Number of BCs (with (k)TVX);#bf{TVX triggered};#bf{#it{N}_{BCs}}", HistType::kTH1F, {{51, -0.5, 50.5}}); - mHistManager.add("Event/nTracks", "Number of BCs (with (k)TVX);#bf{TVX triggered};#bf{#it{N}_{BCs}}", HistType::kTH1F, {{51, -0.5, 50.5}}); - mHistManager.add("Event/nHeavyNeutralMesons", "Number of BCs (with (k)TVX);#bf{TVX triggered};#bf{#it{N}_{BCs}}", HistType::kTH1F, {{51, -0.5, 50.5}}); - - mHistManager.add("nCollisionsVsClusters", "Number of collisions vs number of clusters;N_{Collisions};N_{Clusters}", HistType::kTH2F, {{4, -0.5, 3.5}, {21, -0.5, 20.5}}); - - mHistManager.add("Track/trackP", "Momentum of tracks;#bf{#it{p} (GeV/#it{c})};#bf{#it{N}_{tracks}}", HistType::kTH1F, {{200, 0, 20}}); - - mHistManager.add("Cluster/clusterE", "Energy of cluster;#bf{#it{E} (GeV)};#bf{#it{N}_{clusters}}", HistType::kTH1F, {{200, 0, 20}}); - mHistManager.add("Cluster/clusterM02", "Shape of cluster;#bf{#it{M}_{02}};#bf{#it{N}_{clusters}}", HistType::kTH1F, {{200, 0, 2}}); - mHistManager.add("Cluster/clusterTime", "Time of cluster;#bf{#it{t} (ns)};#bf{#it{N}_{clusters}}", HistType::kTH1F, {{200, -100, 100}}); - - mHistManager.add("GG/invMassVsPt", "Invariant mass and pT of meson candidates", HistType::kTH2F, {{400, massWindowOmega->get("pi0_min"), massWindowEtaPrime->get("eta_max")}, {250, 0., 25.}}); - - mHistManager.add("HeavyNeutralMeson/invMassVsPt", "Invariant mass and pT of meson candidates", HistType::kTH2F, {{600, massWindowOmega->get("omega_min"), massWindowEtaPrime->get("etaprime_max")}, {250, 0., 25.}}); + hCollisionCounter->GetXaxis()->SetBinLabel(5, "PCM #omega"); + hCollisionCounter->GetXaxis()->SetBinLabel(6, "EMC #omega"); + hCollisionCounter->GetXaxis()->SetBinLabel(7, "PCM #eta'"); + hCollisionCounter->GetXaxis()->SetBinLabel(8, "EMC #eta'"); + + mHistManager.add("Event/nGGs", "Number of (selected) #gamma#gamma paris;#bf{#it{N}_{#gamma#gamma}};#bf{#it{N}_{#gamma#gamma}^{selected}}", HistType::kTH2F, {{51, -0.5, 50.5}, {51, -0.5, 50.5}}); + mHistManager.add("Event/nTracks", "Number of tracks;#bf{N_{tracks}};#bf{#it{N}_{Coll}}", HistType::kTH1F, {{51, -0.5, 50.5}}); + mHistManager.add("Event/nHeavyNeutralMesons", "Number of (selected) HNM candidates;#bf{#it{N}_{HNM}};#bf{#it{N}_{HNM}^{selected}}", HistType::kTH2F, {{51, -0.5, 50.5}, {51, -0.5, 50.5}}); + mHistManager.add("Event/nClustersVsV0s", "Number of clusters and V0s in the collision;#bf{#it{N}_{clusters}};#bf{#it{N}_{V0s}}", HistType::kTH2F, {{26, -0.5, 25.5}, {26, -0.5, 25.5}}); + + mHistManager.add("GG/invMassVsPt_PCM", "Invariant mass and pT of gg candidates", HistType::kTH2F, {{400, 0., 0.8}, {250, 0., 25.}}); + mHistManager.add("GG/invMassVsPt_PCMEMC", "Invariant mass and pT of gg candidates", HistType::kTH2F, {{400, 0., 0.8}, {250, 0., 25.}}); + mHistManager.add("GG/invMassVsPt_EMC", "Invariant mass and pT of gg candidates", HistType::kTH2F, {{400, 0., 0.8}, {250, 0., 25.}}); + + mHistManager.add("HeavyNeutralMeson/invMassVsPt_PCM", "Invariant mass and pT of HNM candidates", HistType::kTH2F, {{600, 0.6, 1.2}, {250, 0., 25.}}); + mHistManager.add("HeavyNeutralMeson/invMassVsPt_PCMEMC", "Invariant mass and pT of HNM candidates", HistType::kTH2F, {{600, 0.6, 1.2}, {250, 0., 25.}}); + mHistManager.add("HeavyNeutralMeson/invMassVsPt_EMC", "Invariant mass and pT of HNM candidates", HistType::kTH2F, {{600, 0.6, 1.2}, {250, 0., 25.}}); } - void process(MyCollisions::iterator const& collision, MyBCs const&, SelectedClusters const& clusters, SelectedTracks const& tracks) + Preslice perCollision_pcm = aod::v0photonkf::collisionId; + Preslice perCollision_emc = aod::skimmedcluster::collisionId; + + void process(MyCollisions::iterator const& collision, MyBCs const&, aod::SkimEMCClusters const& clusters, aod::V0PhotonsKF const& v0s, SelectedTracks const& tracks) { mHistManager.fill(HIST("Event/hCollisionCounter"), 0.); - colContainsLowPtOmega = false; - colContainsHighPtOmega = false; - colContainsLowPtEtaPrime = false; - colContainsHighPtEtaPrime = false; + if (collision.foundBC_as().alias_bit(kTVXinEMC)) + mHistManager.fill(HIST("Event/hCollisionCounter"), 1.); - if (isBCSelected(collision.foundBC_as(), clusters.size(), tracks.size() > 2)) { - processClusters(clusters); + auto v0sInThisCollision = v0s.sliceBy(perCollision_pcm, collision.globalIndex()); + auto clustersInThisCollision = clusters.sliceBy(perCollision_emc, collision.globalIndex()); - processGGs(); - processHeavyNeutralMesons(tracks); + mHistManager.fill(HIST("Event/nClustersVsV0s"), clustersInThisCollision.size(), v0sInThisCollision.size()); + mHistManager.fill(HIST("Event/nTracks"), tracks.size()); - mPhotons.clear(); - mGGs.clear(); - } - - if (colContainsLowPtOmega) - mHistManager.fill(HIST("Event/hCollisionCounter"), 4.); - if (colContainsHighPtOmega) - mHistManager.fill(HIST("Event/hCollisionCounter"), 5.); - if (colContainsLowPtEtaPrime) - mHistManager.fill(HIST("Event/hCollisionCounter"), 6.); - if (colContainsHighPtEtaPrime) - mHistManager.fill(HIST("Event/hCollisionCounter"), 7.); - tags(colContainsLowPtOmega, colContainsHighPtOmega, colContainsLowPtEtaPrime, colContainsHighPtEtaPrime); - } + colContainsPCMOmega = colContainsEMCOmega = colContainsPCMEtaPrime = colContainsEMCEtaPrime = false; - bool isBCSelected(MyBCs::iterator const& bc, int nClusters, bool enoughTracks) - { - if (bc.alias_bit(kTVXinEMC)) - mHistManager.fill(HIST("Event/hCollisionCounter"), 1.); - else if (bc.alias_bit(kEMC7) || bc.alias_bit(kEG1) || bc.alias_bit(kEG2) || bc.alias_bit(kDG1) || bc.alias_bit(kDG2) || bc.alias_bit(kEJ1) || bc.alias_bit(kEJ2) || bc.alias_bit(kDJ1) || bc.alias_bit(kDJ2)) - mHistManager.fill(HIST("Event/hCollisionCounter"), 2.); - else - return false; - if (nClusters < 2 || !enoughTracks) - return false; - mHistManager.fill(HIST("Event/hCollisionCounter"), 3.); + reconstructGGs(clustersInThisCollision, v0sInThisCollision, vGGs); + processGGs(vGGs); + reconstructHeavyNeutralMesons(tracks, vGGs, vHNMs); + processHNMs(vHNMs); - return true; + tags(colContainsPCMOmega, colContainsEMCOmega, colContainsPCMEtaPrime, colContainsEMCEtaPrime); } - void processClusters(SelectedClusters const& clusters) + void processGGs(std::vector& vGGs) { - for (const auto& cluster : clusters) { - mHistManager.fill(HIST("Cluster/clusterE"), cluster.energy()); - mHistManager.fill(HIST("Cluster/clusterM02"), cluster.m02()); - mHistManager.fill(HIST("Cluster/clusterTime"), cluster.time()); - mPhotons.push_back(Photon(cluster.eta(), cluster.phi(), cluster.energy())); - } - mHistManager.fill(HIST("Event/nClusters"), clusters.size()); - } + int nGGsBeforeMassCuts = vGGs.size(); + for (unsigned int iGG = 0; iGG < vGGs.size(); iGG++) { + auto lightMeson = &vGGs.at(iGG); + + if (lightMeson->reconstructionType == ReconstructionType::kPCM) { + mHistManager.fill(HIST("GG/invMassVsPt_PCM"), lightMeson->m(), lightMeson->pT()); + } else if (lightMeson->reconstructionType == ReconstructionType::kEMC) { + mHistManager.fill(HIST("GG/invMassVsPt_EMC"), lightMeson->m(), lightMeson->pT()); + } else { + mHistManager.fill(HIST("GG/invMassVsPt_PCMEMC"), lightMeson->m(), lightMeson->pT()); + } - /// \brief Process light neutral meson candidates (pi0 or eta), calculate invariant mass and pT and fill histograms - void processGGs() - { - // loop over all photon combinations and build meson candidates - for (unsigned int ig1 = 0; ig1 < mPhotons.size(); ++ig1) { - for (unsigned int ig2 = ig1 + 1; ig2 < mPhotons.size(); ++ig2) { - Meson lightMeson(mPhotons[ig1], mPhotons[ig2]); // build lightMeson from photons - mHistManager.fill(HIST("GG/invMassVsPt"), lightMeson.m(), lightMeson.pT()); - if (lightMeson.m() > massWindowOmega->get("pi0_min") && lightMeson.m() < massWindowOmega->get("pi0_max")) - lightMeson.addMesonType(HeavyNeutralMeson::kPi0); - else if (lightMeson.m() > massWindowEtaPrime->get("eta_min") && lightMeson.m() < massWindowEtaPrime->get("eta_max")) - lightMeson.addMesonType(HeavyNeutralMeson::kEta); - else - continue; - mGGs.push_back(lightMeson); + if (lightMeson->m() > massWindowOmega->get("pi0_min") && lightMeson->m() < massWindowOmega->get("pi0_max")) { + lightMeson->isPi0 = true; + } else if (lightMeson->m() > massWindowEtaPrime->get("eta_min") && lightMeson->m() < massWindowEtaPrime->get("eta_max")) { + lightMeson->isEta = true; + } else { + vGGs.erase(vGGs.begin() + iGG); + iGG--; } } - mHistManager.fill(HIST("Event/nGGs"), mGGs.size()); + mHistManager.fill(HIST("Event/nGGs"), nGGsBeforeMassCuts, vGGs.size()); } - void processHeavyNeutralMesons(SelectedTracks const& tracks) + void processHNMs(std::vector& vHNMs) { - for (const auto& posTrack : tracks) { - if (!posTrack.isGlobalTrack() || posTrack.sign() < 0) - continue; - for (const auto& negTrack : tracks) { - if (!negTrack.isGlobalTrack() || negTrack.sign() > 0) - continue; - for (unsigned long iGG = 0; iGG < mGGs.size(); iGG++) { - Meson heavyNeutralMeson(mGGs.at(iGG), posTrack.energy(constants::physics::MassPiPlus) + negTrack.energy(constants::physics::MassPiPlus), posTrack.px() + negTrack.px(), posTrack.py() + negTrack.py(), posTrack.pz() + negTrack.pz()); - mHistManager.fill(HIST("HeavyNeutralMeson/invMassVsPt"), heavyNeutralMeson.m(), heavyNeutralMeson.pT()); - if (mGGs.at(iGG).hasMesonType(HeavyNeutralMeson::kPi0) && heavyNeutralMeson.m() > massWindowOmega->get("omega_min") && heavyNeutralMeson.m() < massWindowOmega->get("omega_max")) { - if(minHNMPt->get("omega_lowPt") < heavyNeutralMeson.pT()) - colContainsLowPtOmega = true; - if(minHNMPt->get("omega_highPt") < heavyNeutralMeson.pT()) - colContainsHighPtOmega = true; - } - if (mGGs.at(iGG).hasMesonType(HeavyNeutralMeson::kEta) && heavyNeutralMeson.m() > massWindowEtaPrime->get("etaprime_min") && heavyNeutralMeson.m() < massWindowEtaPrime->get("etaprime_max")) { - if(minHNMPt->get("etaPrime_lowPt") < heavyNeutralMeson.pT()) - colContainsLowPtEtaPrime = true; - if(minHNMPt->get("etaPrime_highPt") < heavyNeutralMeson.pT()) - colContainsHighPtEtaPrime = true; - } - } + int nHNMsBeforeMassCuts = vHNMs.size(); + for (unsigned int iHNM = 0; iHNM < vHNMs.size(); iHNM++) { + auto heavyNeutralMeson = vHNMs.at(iHNM); + + float massHNM = heavyNeutralMeson.m(cfgHNMMassCorrection); + if (heavyNeutralMeson.gg->reconstructionType == ReconstructionType::kPCM) { + mHistManager.fill(HIST("HeavyNeutralMeson/invMassVsPt_PCM"), massHNM, heavyNeutralMeson.pT()); + } else if (heavyNeutralMeson.gg->reconstructionType == ReconstructionType::kEMC) { + mHistManager.fill(HIST("HeavyNeutralMeson/invMassVsPt_EMC"), massHNM, heavyNeutralMeson.pT()); + } else { + mHistManager.fill(HIST("HeavyNeutralMeson/invMassVsPt_PCMEMC"), massHNM, heavyNeutralMeson.pT()); + } + + if (heavyNeutralMeson.gg->isPi0 && massHNM > massWindowOmega->get("omega_min") && massHNM < massWindowOmega->get("omega_max")) { + if (heavyNeutralMeson.gg->reconstructionType == ReconstructionType::kPCM) + colContainsPCMOmega = true; + else if (heavyNeutralMeson.gg->reconstructionType == ReconstructionType::kEMC) + colContainsEMCOmega = true; + } else if (heavyNeutralMeson.gg->isEta && massHNM > massWindowEtaPrime->get("etaprime_min") && massHNM < massWindowEtaPrime->get("etaprime_max")) { + if (heavyNeutralMeson.gg->reconstructionType == ReconstructionType::kPCM) + colContainsPCMEtaPrime = true; + else if (heavyNeutralMeson.gg->reconstructionType == ReconstructionType::kEMC) + colContainsEMCEtaPrime = true; + } else { + vHNMs.erase(vHNMs.begin() + iHNM); + iHNM--; } } + mHistManager.fill(HIST("Event/nHeavyNeutralMesons"), nHNMsBeforeMassCuts, vHNMs.size()); + + if (colContainsPCMOmega) + mHistManager.fill(HIST("Event/hCollisionCounter"), 4.); + if (colContainsEMCOmega) + mHistManager.fill(HIST("Event/hCollisionCounter"), 5.); + if (colContainsPCMEtaPrime) + mHistManager.fill(HIST("Event/hCollisionCounter"), 6.); + if (colContainsEMCEtaPrime) + mHistManager.fill(HIST("Event/hCollisionCounter"), 7.); } }; diff --git a/EventFiltering/filterTables.h b/EventFiltering/filterTables.h index 32574ed65f9..b3a30d24292 100644 --- a/EventFiltering/filterTables.h +++ b/EventFiltering/filterTables.h @@ -176,10 +176,10 @@ DECLARE_SOA_COLUMN(PCMHighPtPhoton, hasPCMHighPtPhoton, bool); //! PCM high pT p DECLARE_SOA_COLUMN(PCMandEE, hasPCMandEE, bool); //! PCM and ee // heavy meson filters -DECLARE_SOA_COLUMN(LowPtOmegaMeson, hasLowPtOmegaMeson, bool); //! Omega meson candidate (3pi) in the collision -DECLARE_SOA_COLUMN(HighPtOmegaMeson, hasHighPtOmegaMeson, bool); //! Omega meson candidate (3pi) in the collision -DECLARE_SOA_COLUMN(LowPtEtaPrimeMeson, hasLowPtEtaPrimeMeson, bool); //! Eta prime candidate (pi+pi-eta) in the collision -DECLARE_SOA_COLUMN(HighPtEtaPrimeMeson, hasHighPtEtaPrimeMeson, bool); //! Eta prime candidate (pi+pi-eta) in the collision +DECLARE_SOA_COLUMN(PCMOmegaMeson, hasPCMOmegaMeson, bool); //! Omega meson candidate (3pi) in the collision +DECLARE_SOA_COLUMN(EMCOmegaMeson, hasEMCOmegaMeson, bool); //! Omega meson candidate (3pi) in the collision +DECLARE_SOA_COLUMN(PCMEtaPrimeMeson, hasPCMEtaPrimeMeson, bool); //! Eta' meson candidate (3pi) in the collision +DECLARE_SOA_COLUMN(EMCEtaPrimeMeson, hasEMCEtaPrimeMeson, bool); //! Eta' meson candidate (3pi) in the collision } // namespace filtering namespace decision @@ -306,8 +306,8 @@ using PhotonFilter = PhotonFilters::iterator; // heavy mesons DECLARE_SOA_TABLE(HeavyNeutralMesonFilters, "AOD", "HeavyNeutralMesonFilters", //! - filtering::LowPtOmegaMeson, filtering::HighPtOmegaMeson, - filtering::LowPtEtaPrimeMeson, filtering::HighPtEtaPrimeMeson); + filtering::PCMOmegaMeson, filtering::EMCOmegaMeson, + filtering::PCMEtaPrimeMeson, filtering::EMCEtaPrimeMeson); using HeavyNeutralMesonFilter = HeavyNeutralMesonFilters::iterator; diff --git a/PWGEM/PhotonMeson/TableProducer/skimmerGammaCalo.cxx b/PWGEM/PhotonMeson/TableProducer/skimmerGammaCalo.cxx index c7da3fa270f..701aa70e195 100644 --- a/PWGEM/PhotonMeson/TableProducer/skimmerGammaCalo.cxx +++ b/PWGEM/PhotonMeson/TableProducer/skimmerGammaCalo.cxx @@ -50,6 +50,7 @@ struct skimmerGammaCalo { Configurable minM02{"minM02", 0.0, "Minimum M02 for M02 cut"}; Configurable maxM02{"maxM02", 1.0, "Maximum M02 for M02 cut"}; Configurable minE{"minE", 0.5, "Minimum energy for energy cut"}; + Configurable> clusterDefinitions{"clusterDefinitions", {0, 1, 2, 10, 11, 12, 13, 20, 21, 22, 30, 40, 41, 42, 43, 44, 45}, "Cluster definitions to be accepted (e.g. 13 for kV3MostSplitLowSeed)"}; Configurable maxdEta{"maxdEta", 0.1, "Set a maximum difference in eta for tracks and cluster to still count as matched"}; Configurable maxdPhi{"maxdPhi", 0.1, "Set a maximum difference in phi for tracks and cluster to still count as matched"}; Configurable applyEveSel_at_skimming{"applyEveSel_at_skimming", false, "flag to apply minimal event selection at the skimming level"}; @@ -59,15 +60,23 @@ struct skimmerGammaCalo { void init(o2::framework::InitContext& initContext) { - historeg.add("hCaloClusterEIn", "hCaloClusterEIn", gHistoSpec_clusterE); - historeg.add("hCaloClusterEOut", "hCaloClusterEOut", gHistoSpec_clusterE); - historeg.add("hMTEtaPhi", "hMTEtaPhi", gHistoSpec_clusterTM_dEtadPhi); - auto hCaloClusterFilter = historeg.add("hCaloClusterFilter", "hCaloClusterFilter", kTH1I, {{5, 0, 5}}); + historeg.add("DefinitionIn", "Cluster definitions before cuts;#bf{Cluster definition};#bf{#it{N}_{clusters}}", HistType::kTH1F, {{51, -0.5, 50.5}}); + historeg.add("DefinitionOut", "Cluster definitions after cuts;#bf{Cluster definition};#bf{#it{N}_{clusters}}", HistType::kTH1F, {{51, -0.5, 50.5}}); + historeg.add("EIn", "Energy of clusters before cuts", gHistoSpec_clusterE); + historeg.add("EOut", "Energy of clusters after cuts", gHistoSpec_clusterE); + historeg.add("MTEtaPhi", "Eta phi of matched tracks", gHistoSpec_clusterTM_dEtadPhi); + historeg.add("M02In", "Shape of cluster before cuts;#bf{#it{M}_{02}};#bf{#it{N}_{clusters}}", HistType::kTH1F, {{200, 0, 2}}); + historeg.add("M02Out", "Shape of cluster after cuts;#bf{#it{M}_{02}};#bf{#it{N}_{clusters}}", HistType::kTH1F, {{200, 0, 2}}); + historeg.add("TimeIn", "Time of cluster before cuts;#bf{#it{t} (ns)};#bf{#it{N}_{clusters}}", HistType::kTH1F, {{200, -100, 100}}); + historeg.add("TimeOut", "Time of cluster after cuts;#bf{#it{t} (ns)};#bf{#it{N}_{clusters}}", HistType::kTH1F, {{200, -100, 100}}); + + auto hCaloClusterFilter = historeg.add("hCaloClusterFilter", "hCaloClusterFilter", kTH1I, {{6, 0, 6}}); hCaloClusterFilter->GetXaxis()->SetBinLabel(1, "in"); - hCaloClusterFilter->GetXaxis()->SetBinLabel(2, "E cut"); - hCaloClusterFilter->GetXaxis()->SetBinLabel(3, "time cut"); - hCaloClusterFilter->GetXaxis()->SetBinLabel(4, "M02 cut"); - hCaloClusterFilter->GetXaxis()->SetBinLabel(5, "out"); + hCaloClusterFilter->GetXaxis()->SetBinLabel(2, "Definition cut"); + hCaloClusterFilter->GetXaxis()->SetBinLabel(3, "E cut"); + hCaloClusterFilter->GetXaxis()->SetBinLabel(4, "time cut"); + hCaloClusterFilter->GetXaxis()->SetBinLabel(5, "M02 cut"); + hCaloClusterFilter->GetXaxis()->SetBinLabel(6, "out"); if (inherit_from_emevent_photon) { getTaskOptionValue(initContext, "create-emevent-photon", "applyEveSel_at_skimming", applyEveSel_at_skimming.value, true); // for EM users. @@ -84,24 +93,40 @@ struct skimmerGammaCalo { return; } for (const auto& emccluster : emcclusters) { - historeg.fill(HIST("hCaloClusterEIn"), emccluster.energy()); historeg.fill(HIST("hCaloClusterFilter"), 0); + historeg.fill(HIST("DefinitionIn"), emccluster.definition()); + historeg.fill(HIST("EIn"), emccluster.energy()); + historeg.fill(HIST("M02In"), emccluster.m02()); + historeg.fill(HIST("TimeIn"), emccluster.time()); + + // Definition cut + if (!(std::find(clusterDefinitions.value.begin(), clusterDefinitions.value.end(), emccluster.definition()) != clusterDefinitions.value.end())) { + historeg.fill(HIST("hCaloClusterFilter"), 1); + continue; + } + historeg.fill(HIST("EIn"), emccluster.energy()); // Energy cut if (emccluster.energy() < minE) { - historeg.fill(HIST("hCaloClusterFilter"), 1); + historeg.fill(HIST("hCaloClusterFilter"), 2); continue; } // timing cut if (emccluster.time() > maxTime || emccluster.time() < minTime) { - historeg.fill(HIST("hCaloClusterFilter"), 2); + historeg.fill(HIST("hCaloClusterFilter"), 3); continue; } // M02 cut if (emccluster.nCells() > 1 && (emccluster.m02() > maxM02 || emccluster.m02() < minM02)) { - historeg.fill(HIST("hCaloClusterFilter"), 3); + historeg.fill(HIST("hCaloClusterFilter"), 4); continue; } + historeg.fill(HIST("hCaloClusterFilter"), 5); + + historeg.fill(HIST("DefinitionOut"), emccluster.definition()); + historeg.fill(HIST("EOut"), emccluster.energy()); + historeg.fill(HIST("M02Out"), emccluster.m02()); + historeg.fill(HIST("TimeOut"), emccluster.time()); // Skimmed cell table auto groupedCells = emcclustercells.sliceBy(CellperCluster, emccluster.globalIndex()); @@ -125,7 +150,7 @@ struct skimmerGammaCalo { if (std::abs(emccluster.eta() - emcmatchedtrack.track_as().trackEtaEmcal()) >= maxdEta || std::abs(emccluster.phi() - emcmatchedtrack.track_as().trackPhiEmcal()) >= maxdPhi) { continue; } - historeg.fill(HIST("hMTEtaPhi"), emccluster.eta() - emcmatchedtrack.track_as().trackEtaEmcal(), emccluster.phi() - emcmatchedtrack.track_as().trackPhiEmcal()); + historeg.fill(HIST("MTEtaPhi"), emccluster.eta() - emcmatchedtrack.track_as().trackEtaEmcal(), emccluster.phi() - emcmatchedtrack.track_as().trackPhiEmcal()); vTrackIds.emplace_back(emcmatchedtrack.trackId()); vEta.emplace_back(emcmatchedtrack.track_as().trackEtaEmcal()); vPhi.emplace_back(emcmatchedtrack.track_as().trackPhiEmcal()); @@ -135,9 +160,6 @@ struct skimmerGammaCalo { emcmatchedtrack.track_as().p(), emcmatchedtrack.track_as().pt()); } - historeg.fill(HIST("hCaloClusterEOut"), emccluster.energy()); - historeg.fill(HIST("hCaloClusterFilter"), 4); - tableGammaEMCReco(emccluster.collisionId(), emccluster.definition(), emccluster.energy(), emccluster.eta(), emccluster.phi(), emccluster.m02(), emccluster.nCells(), emccluster.time(), emccluster.isExotic(), vEta, vPhi, vP, vPt); } @@ -149,6 +171,10 @@ struct skimmerGammaCalo { } for (const auto& emccluster : emcclusters) { + // Definition cut + if (!(std::find(clusterDefinitions.value.begin(), clusterDefinitions.value.end(), emccluster.definition()) != clusterDefinitions.value.end())) { + continue; + } // Energy cut if (emccluster.energy() < minE) { continue; diff --git a/PWGEM/PhotonMeson/Utils/HNMUtilities.h b/PWGEM/PhotonMeson/Utils/HNMUtilities.h new file mode 100644 index 00000000000..fcc424370a4 --- /dev/null +++ b/PWGEM/PhotonMeson/Utils/HNMUtilities.h @@ -0,0 +1,175 @@ +// 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 HNMUtilities.h +/// +/// \brief This code provides helper functions for the reconstruction of heavy neutral mesons (omega and eta meson) via their three pion decay +/// +/// \author Nicolas Strangmann (nicolas.strangmann@cern.ch) - Goethe University Frankfurt +/// + +#ifndef PWGEM_PHOTONMESON_UTILS_HNMUTILITIES_H_ +#define PWGEM_PHOTONMESON_UTILS_HNMUTILITIES_H_ + +#include +#include +#include +#include "TVector3.h" + +#include "Framework/runDataProcessing.h" +#include "Framework/AnalysisTask.h" +#include "Framework/HistogramRegistry.h" +#include "Common/DataModel/EventSelection.h" +#include "EMCALBase/Geometry.h" +#include "PWGJE/DataModel/EMCALClusters.h" +#include "PWGEM/PhotonMeson/DataModel/gammaTables.h" +#include "Common/Core/trackUtilities.h" +#include "Common/DataModel/TrackSelectionTables.h" +#include "EventFiltering/filterTables.h" + +namespace o2::aod::pwgem::photonmeson::utils::hnmutilities +{ +struct Photon { + Photon(float px, float py, float pz, bool isFromConversion) : px(px), py(py), pz(pz), pt(std::sqrt(px * px + py * py)), isFromConversion(isFromConversion) + { + photon.SetPxPyPzE(px, py, pz, std::sqrt(px * px + py * py + pz * pz)); + } + + static Photon fromPxPyPz(float px, float py, float pz) + { + Photon photon(px, py, pz, true); + return photon; + } + + static Photon fromEtaPhiEnergy(float eta, float phi, float energy) + { + float theta = 2 * std::atan2(std::exp(-eta), 1); + float px = energy * std::sin(theta) * std::cos(phi); + float py = energy * std::sin(theta) * std::sin(phi); + float pz = energy * std::cos(theta); + Photon photon(px, py, pz, false); + return photon; + } + + ROOT::Math::PxPyPzEVector photon; + float px, py, pz, pt; + bool isFromConversion; +}; + +namespace ReconstructionType +{ +enum ReconstructionType { + kPCM = 0, + kPCMEMC = 1, + kEMC = 2 +}; +} // namespace ReconstructionType + +enum MassCorrectionType { + kNoHNMMassCorrection = 0, + kSubDeltaPi0 = 1, + kSubLambda = 2 +}; + +struct GammaGammaPair { + GammaGammaPair(Photon p1, Photon p2) : p1(p1), p2(p2) + { + vGG = p1.photon + p2.photon; + } + Photon p1, p2; + ROOT::Math::PxPyPzEVector vGG; + + bool isPi0 = false; + bool isEta = false; + ushort reconstructionType; + void setReconstructionType(ushort type) { reconstructionType = type; } + + float m() const { return vGG.M(); } + float pT() const { return vGG.Pt(); } +}; + +struct HeavyNeutralMeson { + HeavyNeutralMeson(GammaGammaPair* gg, float eTracks, float pxTracks, float pyTracks, float pzTracks) : gg(gg) + { + vHeavyNeutralMeson.SetPxPyPzE(gg->vGG.Px() + pxTracks, gg->vGG.Py() + pyTracks, gg->vGG.Pz() + pzTracks, gg->vGG.e() + eTracks); + } + + GammaGammaPair* gg = nullptr; + ROOT::Math::PxPyPzEVector vHeavyNeutralMeson; + + float m(int massCorrectionType) const + { + float massHNM = vHeavyNeutralMeson.M(); + switch (massCorrectionType) { + case kSubDeltaPi0: + massHNM -= gg->m(); + massHNM += (gg->isPi0 ? constants::physics::MassPi0 : 0.547862); + break; + case kSubLambda: + LOGF(warning, "SubLambdaMassCorrection not yet implemented!"); + break; + default: + LOGF(fatal, "Unknown mass correction type %d", massCorrectionType); + } + return massHNM; + } + float pT() const { return vHeavyNeutralMeson.Pt(); } +}; + +/// \brief Process light neutral meson candidates (pi0 or eta), calculate invariant mass and pT and fill histograms +template +void reconstructGGs(C clusters, V v0s, std::vector& vGGs) +{ + std::vector vPhotons; + for (const auto& cluster : clusters) + vPhotons.push_back(Photon::fromEtaPhiEnergy(cluster.eta(), cluster.phi(), cluster.e())); + + for (const auto& v0 : v0s) + vPhotons.push_back(Photon::fromPxPyPz(v0.px(), v0.py(), v0.pz())); + + vGGs.clear(); + // loop over all photon combinations and build meson candidates + for (unsigned int ig1 = 0; ig1 < vPhotons.size(); ++ig1) { + for (unsigned int ig2 = ig1 + 1; ig2 < vPhotons.size(); ++ig2) { + GammaGammaPair lightMeson(vPhotons[ig1], vPhotons[ig2]); // build lightMeson from photons + if (vPhotons[ig1].isFromConversion && vPhotons[ig2].isFromConversion) + lightMeson.setReconstructionType(ReconstructionType::kPCM); + else if (!vPhotons[ig1].isFromConversion && !vPhotons[ig2].isFromConversion) + lightMeson.setReconstructionType(ReconstructionType::kEMC); + else + lightMeson.setReconstructionType(ReconstructionType::kPCMEMC); + + vGGs.push_back(lightMeson); + } + } +} + +template +void reconstructHeavyNeutralMesons(Track const& tracks, std::vector& vGGs, std::vector& vHNMs) +{ + vHNMs.clear(); + for (const auto& posTrack : tracks) { + if (!posTrack.isGlobalTrack() || posTrack.sign() < 0) + continue; + for (const auto& negTrack : tracks) { + if (!negTrack.isGlobalTrack() || negTrack.sign() > 0) + continue; + for (unsigned long iGG = 0; iGG < vGGs.size(); iGG++) { + HeavyNeutralMeson heavyNeutralMeson(&vGGs.at(iGG), posTrack.energy(constants::physics::MassPiPlus) + negTrack.energy(constants::physics::MassPiMinus), posTrack.px() + negTrack.px(), posTrack.py() + negTrack.py(), posTrack.pz() + negTrack.pz()); + vHNMs.push_back(heavyNeutralMeson); + } + } + } +} + +} // namespace o2::aod::pwgem::photonmeson::utils::hnmutilities + +#endif // PWGEM_PHOTONMESON_UTILS_HNMUTILITIES_H_ From 9e3b146fe910ec2a1d78419d2063585d2cfdb392 Mon Sep 17 00:00:00 2001 From: ALICE Action Bot Date: Sun, 9 Mar 2025 19:31:30 +0000 Subject: [PATCH 3/7] Please consider the following formatting changes --- EventFiltering/filterTables.h | 40 +++++++++++++++++------------------ 1 file changed, 20 insertions(+), 20 deletions(-) diff --git a/EventFiltering/filterTables.h b/EventFiltering/filterTables.h index b3a30d24292..5c450b02eb2 100644 --- a/EventFiltering/filterTables.h +++ b/EventFiltering/filterTables.h @@ -71,25 +71,25 @@ DECLARE_SOA_COLUMN(LMeeHMR, hasLMeeHMR, bool); //! dielectron trigger for high m DECLARE_SOA_COLUMN(ElectronMuon, hasElectronMuon, bool); //! dimuon trigger with low pT on muons // heavy flavours -DECLARE_SOA_COLUMN(HfHighPt2P, hasHfHighPt2P, bool); //! high-pT 2-prong charm hadron -DECLARE_SOA_COLUMN(HfHighPt3P, hasHfHighPt3P, bool); //! high-pT 3-prong charm hadron -DECLARE_SOA_COLUMN(HfBeauty3P, hasHfBeauty3P, bool); //! 3-prong beauty hadron -DECLARE_SOA_COLUMN(HfBeauty4P, hasHfBeauty4P, bool); //! 4-prong beauty hadron -DECLARE_SOA_COLUMN(HfFemto2P, hasHfFemto2P, bool); //! 2-prong charm-hadron - N pair -DECLARE_SOA_COLUMN(HfFemto3P, hasHfFemto3P, bool); //! 3-prong charm-hadron - N pair -DECLARE_SOA_COLUMN(HfDoubleCharm2P, hasHfDoubleCharm2P, bool); //! at least two 2-prong charm-hadron candidates -DECLARE_SOA_COLUMN(HfDoubleCharm3P, hasHfDoubleCharm3P, bool); //! at least two 3-prong charm-hadron candidates -DECLARE_SOA_COLUMN(HfDoubleCharmMix, hasHfDoubleCharmMix, bool); //! at least one 2-prong and one 3-prong charm-hadron candidates -DECLARE_SOA_COLUMN(HfV0Charm2P, hasHfV0Charm2P, bool); //! V0 with 2-prong charm hadron -DECLARE_SOA_COLUMN(HfV0Charm3P, hasHfV0Charm3P, bool); //! V0 with 3-prong charm hadron -DECLARE_SOA_COLUMN(HfCharmBarToXiBach, hasHfCharmBarToXiBach, bool); //! Charm baryon to Xi + bachelor +DECLARE_SOA_COLUMN(HfHighPt2P, hasHfHighPt2P, bool); //! high-pT 2-prong charm hadron +DECLARE_SOA_COLUMN(HfHighPt3P, hasHfHighPt3P, bool); //! high-pT 3-prong charm hadron +DECLARE_SOA_COLUMN(HfBeauty3P, hasHfBeauty3P, bool); //! 3-prong beauty hadron +DECLARE_SOA_COLUMN(HfBeauty4P, hasHfBeauty4P, bool); //! 4-prong beauty hadron +DECLARE_SOA_COLUMN(HfFemto2P, hasHfFemto2P, bool); //! 2-prong charm-hadron - N pair +DECLARE_SOA_COLUMN(HfFemto3P, hasHfFemto3P, bool); //! 3-prong charm-hadron - N pair +DECLARE_SOA_COLUMN(HfDoubleCharm2P, hasHfDoubleCharm2P, bool); //! at least two 2-prong charm-hadron candidates +DECLARE_SOA_COLUMN(HfDoubleCharm3P, hasHfDoubleCharm3P, bool); //! at least two 3-prong charm-hadron candidates +DECLARE_SOA_COLUMN(HfDoubleCharmMix, hasHfDoubleCharmMix, bool); //! at least one 2-prong and one 3-prong charm-hadron candidates +DECLARE_SOA_COLUMN(HfV0Charm2P, hasHfV0Charm2P, bool); //! V0 with 2-prong charm hadron +DECLARE_SOA_COLUMN(HfV0Charm3P, hasHfV0Charm3P, bool); //! V0 with 3-prong charm hadron +DECLARE_SOA_COLUMN(HfCharmBarToXiBach, hasHfCharmBarToXiBach, bool); //! Charm baryon to Xi + bachelor DECLARE_SOA_COLUMN(HfCharmBarToXi2Bach, hasHfCharmBarToXi2Bach, bool); //! Charm baryon to Xi + 2 bachelors -DECLARE_SOA_COLUMN(HfSigmaCPPK, hasHfSigmaCPPK, bool); //! SigmaC(2455)++K- and SigmaC(2520)++K- + c.c. -DECLARE_SOA_COLUMN(HfSigmaC0K0, hasHfSigmaC0K0, bool); //! SigmaC(2455)0KS0 and SigmaC(2520)0KS0 -DECLARE_SOA_COLUMN(HfPhotonCharm2P, hasHfPhotonCharm2P, bool); //! photon with 2-prong charm hadron -DECLARE_SOA_COLUMN(HfPhotonCharm3P, hasHfPhotonCharm3P, bool); //! photon with 3-prong charm hadron -DECLARE_SOA_COLUMN(HfSingleCharm2P, hasHfSingleCharm2P, bool); //! 2-prong charm hadron (for efficiency studies) -DECLARE_SOA_COLUMN(HfSingleCharm3P, hasHfSingleCharm3P, bool); //! 3-prong charm hadron (for efficiency studies) +DECLARE_SOA_COLUMN(HfSigmaCPPK, hasHfSigmaCPPK, bool); //! SigmaC(2455)++K- and SigmaC(2520)++K- + c.c. +DECLARE_SOA_COLUMN(HfSigmaC0K0, hasHfSigmaC0K0, bool); //! SigmaC(2455)0KS0 and SigmaC(2520)0KS0 +DECLARE_SOA_COLUMN(HfPhotonCharm2P, hasHfPhotonCharm2P, bool); //! photon with 2-prong charm hadron +DECLARE_SOA_COLUMN(HfPhotonCharm3P, hasHfPhotonCharm3P, bool); //! photon with 3-prong charm hadron +DECLARE_SOA_COLUMN(HfSingleCharm2P, hasHfSingleCharm2P, bool); //! 2-prong charm hadron (for efficiency studies) +DECLARE_SOA_COLUMN(HfSingleCharm3P, hasHfSingleCharm3P, bool); //! 3-prong charm hadron (for efficiency studies) DECLARE_SOA_COLUMN(HfSingleNonPromptCharm2P, hasHfSingleNonPromptCharm2P, bool); //! 2-prong charm hadron (for efficiency studies) DECLARE_SOA_COLUMN(HfSingleNonPromptCharm3P, hasHfSingleNonPromptCharm3P, bool); //! 3-prong charm hadron (for efficiency studies) DECLARE_SOA_COLUMN(HfBtoJPsiKa, hasHfBtoJPsiKa, bool); //! B+ -> JPsi(->mumu)K+ @@ -178,8 +178,8 @@ DECLARE_SOA_COLUMN(PCMandEE, hasPCMandEE, bool); //! PCM and ee // heavy meson filters DECLARE_SOA_COLUMN(PCMOmegaMeson, hasPCMOmegaMeson, bool); //! Omega meson candidate (3pi) in the collision DECLARE_SOA_COLUMN(EMCOmegaMeson, hasEMCOmegaMeson, bool); //! Omega meson candidate (3pi) in the collision -DECLARE_SOA_COLUMN(PCMEtaPrimeMeson, hasPCMEtaPrimeMeson, bool); //! Eta' meson candidate (3pi) in the collision -DECLARE_SOA_COLUMN(EMCEtaPrimeMeson, hasEMCEtaPrimeMeson, bool); //! Eta' meson candidate (3pi) in the collision +DECLARE_SOA_COLUMN(PCMEtaPrimeMeson, hasPCMEtaPrimeMeson, bool); //! Eta' meson candidate (3pi) in the collision +DECLARE_SOA_COLUMN(EMCEtaPrimeMeson, hasEMCEtaPrimeMeson, bool); //! Eta' meson candidate (3pi) in the collision } // namespace filtering namespace decision From fc5972341ba5747bf91c2616e811bf3ebba4ebb6 Mon Sep 17 00:00:00 2001 From: Nicolas Strangmann Date: Sun, 9 Mar 2025 20:38:56 +0100 Subject: [PATCH 4/7] Make MegaLinter happy --- EventFiltering/PWGEM/HeavyNeutralMesonFilter.cxx | 2 ++ PWGEM/PhotonMeson/Utils/HNMUtilities.h | 2 +- 2 files changed, 3 insertions(+), 1 deletion(-) diff --git a/EventFiltering/PWGEM/HeavyNeutralMesonFilter.cxx b/EventFiltering/PWGEM/HeavyNeutralMesonFilter.cxx index b09fc89f014..83731276fe2 100644 --- a/EventFiltering/PWGEM/HeavyNeutralMesonFilter.cxx +++ b/EventFiltering/PWGEM/HeavyNeutralMesonFilter.cxx @@ -16,6 +16,8 @@ /// \author Nicolas Strangmann (nicolas.strangmann@cern.ch) - Goethe University Frankfurt /// +#include + #include "PWGEM/PhotonMeson/Utils/HNMUtilities.h" using namespace o2; diff --git a/PWGEM/PhotonMeson/Utils/HNMUtilities.h b/PWGEM/PhotonMeson/Utils/HNMUtilities.h index fcc424370a4..5d58216cfc8 100644 --- a/PWGEM/PhotonMeson/Utils/HNMUtilities.h +++ b/PWGEM/PhotonMeson/Utils/HNMUtilities.h @@ -162,7 +162,7 @@ void reconstructHeavyNeutralMesons(Track const& tracks, std::vector 0) continue; - for (unsigned long iGG = 0; iGG < vGGs.size(); iGG++) { + for (const auto iGG = 0; iGG < vGGs.size(); iGG++) { HeavyNeutralMeson heavyNeutralMeson(&vGGs.at(iGG), posTrack.energy(constants::physics::MassPiPlus) + negTrack.energy(constants::physics::MassPiMinus), posTrack.px() + negTrack.px(), posTrack.py() + negTrack.py(), posTrack.pz() + negTrack.pz()); vHNMs.push_back(heavyNeutralMeson); } From a21d36f69345b4d9c7696608e50b56cc5e0d8dee Mon Sep 17 00:00:00 2001 From: Nicolas Strangmann Date: Sun, 9 Mar 2025 21:44:07 +0100 Subject: [PATCH 5/7] Remove const in for loop --- PWGEM/PhotonMeson/Utils/HNMUtilities.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PWGEM/PhotonMeson/Utils/HNMUtilities.h b/PWGEM/PhotonMeson/Utils/HNMUtilities.h index 5d58216cfc8..6e85e8e91a8 100644 --- a/PWGEM/PhotonMeson/Utils/HNMUtilities.h +++ b/PWGEM/PhotonMeson/Utils/HNMUtilities.h @@ -162,7 +162,7 @@ void reconstructHeavyNeutralMesons(Track const& tracks, std::vector 0) continue; - for (const auto iGG = 0; iGG < vGGs.size(); iGG++) { + for (auto iGG = 0; iGG < vGGs.size(); iGG++) { HeavyNeutralMeson heavyNeutralMeson(&vGGs.at(iGG), posTrack.energy(constants::physics::MassPiPlus) + negTrack.energy(constants::physics::MassPiMinus), posTrack.px() + negTrack.px(), posTrack.py() + negTrack.py(), posTrack.pz() + negTrack.pz()); vHNMs.push_back(heavyNeutralMeson); } From 87fd3b94da23d844389732a99addc24383fb12d1 Mon Sep 17 00:00:00 2001 From: Nicolas Strangmann Date: Sun, 9 Mar 2025 22:01:53 +0100 Subject: [PATCH 6/7] Make MegaLinter happy --- PWGEM/PhotonMeson/Utils/HNMUtilities.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/PWGEM/PhotonMeson/Utils/HNMUtilities.h b/PWGEM/PhotonMeson/Utils/HNMUtilities.h index 6e85e8e91a8..7dea4725f08 100644 --- a/PWGEM/PhotonMeson/Utils/HNMUtilities.h +++ b/PWGEM/PhotonMeson/Utils/HNMUtilities.h @@ -162,7 +162,7 @@ void reconstructHeavyNeutralMesons(Track const& tracks, std::vector 0) continue; - for (auto iGG = 0; iGG < vGGs.size(); iGG++) { + for (size_t iGG = 0; iGG < vGGs.size(); iGG++) { HeavyNeutralMeson heavyNeutralMeson(&vGGs.at(iGG), posTrack.energy(constants::physics::MassPiPlus) + negTrack.energy(constants::physics::MassPiMinus), posTrack.px() + negTrack.px(), posTrack.py() + negTrack.py(), posTrack.pz() + negTrack.pz()); vHNMs.push_back(heavyNeutralMeson); } From d4c65ad789add812d4164b2d7c3cbde624371473 Mon Sep 17 00:00:00 2001 From: Nicolas Strangmann Date: Mon, 10 Mar 2025 11:00:21 +0100 Subject: [PATCH 7/7] Change reconstruction type enum and general cleanup --- .../PWGEM/HeavyNeutralMesonFilter.cxx | 69 ++++++++++--------- PWGEM/PhotonMeson/Utils/HNMUtilities.h | 46 ++++++------- 2 files changed, 59 insertions(+), 56 deletions(-) diff --git a/EventFiltering/PWGEM/HeavyNeutralMesonFilter.cxx b/EventFiltering/PWGEM/HeavyNeutralMesonFilter.cxx index 83731276fe2..1446943c279 100644 --- a/EventFiltering/PWGEM/HeavyNeutralMesonFilter.cxx +++ b/EventFiltering/PWGEM/HeavyNeutralMesonFilter.cxx @@ -11,7 +11,7 @@ /// /// \file HeavyNeutralMesonFilter.cxx /// -/// \brief This code loops over collisions to filter events contaning heavy mesons (omega or eta') using EMCal clusters +/// \brief This code loops over collisions to filter events contaning heavy mesons (omega or eta') using EMCal clusters and V0s (PCM) /// /// \author Nicolas Strangmann (nicolas.strangmann@cern.ch) - Goethe University Frankfurt /// @@ -23,7 +23,7 @@ using namespace o2; using namespace o2::framework; using namespace o2::framework::expressions; -using namespace o2::aod::pwgem::photonmeson::utils::hnmutilities; +using namespace o2::aod::pwgem::photonmeson; using MyBCs = soa::Join; using MyCollisions = soa::Join; @@ -41,10 +41,13 @@ struct HeavyNeutralMesonFilter { Configurable> massWindowOmega{"massWindowOmega", {defaultMassWindows[0], 4, {"pi0_min", "pi0_max", "omega_min", "omega_max"}}, "Mass window for selected omegas and their decay pi0"}; Configurable> massWindowEtaPrime{"massWindowEtaPrime", {defaultMassWindows[1], 4, {"eta_min", "eta_max", "etaprime_min", "etaprime_max"}}, "Mass window for selected eta' and their decay eta"}; + static constexpr float defaultMinPts[4] = {1.8, 1.8, 2.6, 2.6}; + Configurable> minHNMPts{"minHNMPts", {defaultMinPts, 4, {"PCM_omega", "PCM_etaprime", "EMC_omega", "EMC_etaprime"}}, "Minimum pT values for the trigger decisions (GeV/c)"}; + Filter trackPtFilter = aod::track::pt > cfgTrackMinPt; - std::vector vGGs; - std::vector vHNMs; + std::vector vGGs; + std::vector vHNMs; bool colContainsPCMOmega, colContainsEMCOmega, colContainsPCMEtaPrime, colContainsEMCEtaPrime = false; @@ -53,28 +56,26 @@ struct HeavyNeutralMesonFilter { void init(InitContext const&) { emcalGeom = emcal::Geometry::GetInstanceFromRunNumber(300000); - auto hCollisionCounter = mHistManager.add("Event/hCollisionCounter", "Number of collisions;;#bf{#it{N}_{Coll}}", HistType::kTH1F, {{8, -0.5, 7.5}}); + auto hCollisionCounter = mHistManager.add("Event/hCollisionCounter", "Number of collisions;;#bf{#it{N}_{Coll}}", HistType::kTH1F, {{6, -0.5, 5.5}}); hCollisionCounter->GetXaxis()->SetBinLabel(1, "all"); hCollisionCounter->GetXaxis()->SetBinLabel(2, "kTVXinEMC"); - hCollisionCounter->GetXaxis()->SetBinLabel(3, "L0Triggered"); - hCollisionCounter->GetXaxis()->SetBinLabel(4, "2+ tracks&clusters"); - hCollisionCounter->GetXaxis()->SetBinLabel(5, "PCM #omega"); - hCollisionCounter->GetXaxis()->SetBinLabel(6, "EMC #omega"); - hCollisionCounter->GetXaxis()->SetBinLabel(7, "PCM #eta'"); - hCollisionCounter->GetXaxis()->SetBinLabel(8, "EMC #eta'"); + hCollisionCounter->GetXaxis()->SetBinLabel(3, "PCM #omega"); + hCollisionCounter->GetXaxis()->SetBinLabel(4, "EMC #omega"); + hCollisionCounter->GetXaxis()->SetBinLabel(5, "PCM #eta'"); + hCollisionCounter->GetXaxis()->SetBinLabel(6, "EMC #eta'"); mHistManager.add("Event/nGGs", "Number of (selected) #gamma#gamma paris;#bf{#it{N}_{#gamma#gamma}};#bf{#it{N}_{#gamma#gamma}^{selected}}", HistType::kTH2F, {{51, -0.5, 50.5}, {51, -0.5, 50.5}}); mHistManager.add("Event/nTracks", "Number of tracks;#bf{N_{tracks}};#bf{#it{N}_{Coll}}", HistType::kTH1F, {{51, -0.5, 50.5}}); mHistManager.add("Event/nHeavyNeutralMesons", "Number of (selected) HNM candidates;#bf{#it{N}_{HNM}};#bf{#it{N}_{HNM}^{selected}}", HistType::kTH2F, {{51, -0.5, 50.5}, {51, -0.5, 50.5}}); mHistManager.add("Event/nClustersVsV0s", "Number of clusters and V0s in the collision;#bf{#it{N}_{clusters}};#bf{#it{N}_{V0s}}", HistType::kTH2F, {{26, -0.5, 25.5}, {26, -0.5, 25.5}}); - mHistManager.add("GG/invMassVsPt_PCM", "Invariant mass and pT of gg candidates", HistType::kTH2F, {{400, 0., 0.8}, {250, 0., 25.}}); - mHistManager.add("GG/invMassVsPt_PCMEMC", "Invariant mass and pT of gg candidates", HistType::kTH2F, {{400, 0., 0.8}, {250, 0., 25.}}); - mHistManager.add("GG/invMassVsPt_EMC", "Invariant mass and pT of gg candidates", HistType::kTH2F, {{400, 0., 0.8}, {250, 0., 25.}}); + mHistManager.add("GG/invMassVsPt_PCM", "Invariant mass and pT of gg candidates;#bf{#it{M}_{#gamma#gamma}};#bf{#it{N}_{#gamma#gamma}}", HistType::kTH2F, {{400, 0., 0.8}, {250, 0., 25.}}); + mHistManager.add("GG/invMassVsPt_PCMEMC", "Invariant mass and pT of gg candidates;#bf{#it{M}_{#gamma#gamma}};#bf{#it{N}_{#gamma#gamma}}", HistType::kTH2F, {{400, 0., 0.8}, {250, 0., 25.}}); + mHistManager.add("GG/invMassVsPt_EMC", "Invariant mass and pT of gg candidates;#bf{#it{M}_{#gamma#gamma}};#bf{#it{N}_{#gamma#gamma}}", HistType::kTH2F, {{400, 0., 0.8}, {250, 0., 25.}}); - mHistManager.add("HeavyNeutralMeson/invMassVsPt_PCM", "Invariant mass and pT of HNM candidates", HistType::kTH2F, {{600, 0.6, 1.2}, {250, 0., 25.}}); - mHistManager.add("HeavyNeutralMeson/invMassVsPt_PCMEMC", "Invariant mass and pT of HNM candidates", HistType::kTH2F, {{600, 0.6, 1.2}, {250, 0., 25.}}); - mHistManager.add("HeavyNeutralMeson/invMassVsPt_EMC", "Invariant mass and pT of HNM candidates", HistType::kTH2F, {{600, 0.6, 1.2}, {250, 0., 25.}}); + mHistManager.add("HeavyNeutralMeson/invMassVsPt_PCM", "Invariant mass and pT of HNM candidates;#bf{#it{M}_{#pi^{+}#pi^{-}#gamma#gamma}};#bf{#it{N}_{#pi^{+}#pi^{-}#gamma#gamma}}", HistType::kTH2F, {{600, 0.6, 1.2}, {250, 0., 25.}}); + mHistManager.add("HeavyNeutralMeson/invMassVsPt_PCMEMC", "Invariant mass and pT of HNM candidates;#bf{#it{M}_{#pi^{+}#pi^{-}#gamma#gamma}};#bf{#it{N}_{#pi^{+}#pi^{-}#gamma#gamma}}", HistType::kTH2F, {{600, 0.6, 1.2}, {250, 0., 25.}}); + mHistManager.add("HeavyNeutralMeson/invMassVsPt_EMC", "Invariant mass and pT of HNM candidates;#bf{#it{M}_{#pi^{+}#pi^{-}#gamma#gamma}};#bf{#it{N}_{#pi^{+}#pi^{-}#gamma#gamma}}", HistType::kTH2F, {{600, 0.6, 1.2}, {250, 0., 25.}}); } Preslice perCollision_pcm = aod::v0photonkf::collisionId; @@ -95,23 +96,24 @@ struct HeavyNeutralMesonFilter { colContainsPCMOmega = colContainsEMCOmega = colContainsPCMEtaPrime = colContainsEMCEtaPrime = false; - reconstructGGs(clustersInThisCollision, v0sInThisCollision, vGGs); + hnmutilities::reconstructGGs(clustersInThisCollision, v0sInThisCollision, vGGs); processGGs(vGGs); - reconstructHeavyNeutralMesons(tracks, vGGs, vHNMs); + hnmutilities::reconstructHeavyNeutralMesons(tracks, vGGs, vHNMs); processHNMs(vHNMs); tags(colContainsPCMOmega, colContainsEMCOmega, colContainsPCMEtaPrime, colContainsEMCEtaPrime); } - void processGGs(std::vector& vGGs) + /// \brief Loop over the GG candidates, fill the mass/pt histograms and set the isPi0/isEta flags based on the reconstructed mass + void processGGs(std::vector& vGGs) { int nGGsBeforeMassCuts = vGGs.size(); for (unsigned int iGG = 0; iGG < vGGs.size(); iGG++) { auto lightMeson = &vGGs.at(iGG); - if (lightMeson->reconstructionType == ReconstructionType::kPCM) { + if (lightMeson->reconstructionType == photonpair::kPCMPCM) { mHistManager.fill(HIST("GG/invMassVsPt_PCM"), lightMeson->m(), lightMeson->pT()); - } else if (lightMeson->reconstructionType == ReconstructionType::kEMC) { + } else if (lightMeson->reconstructionType == photonpair::kEMCEMC) { mHistManager.fill(HIST("GG/invMassVsPt_EMC"), lightMeson->m(), lightMeson->pT()); } else { mHistManager.fill(HIST("GG/invMassVsPt_PCMEMC"), lightMeson->m(), lightMeson->pT()); @@ -129,30 +131,31 @@ struct HeavyNeutralMesonFilter { mHistManager.fill(HIST("Event/nGGs"), nGGsBeforeMassCuts, vGGs.size()); } - void processHNMs(std::vector& vHNMs) + /// \brief Loop over the heavy neutral meson candidates, fill the mass/pt histograms and set the trigger flags based on the reconstructed mass + void processHNMs(std::vector& vHNMs) { int nHNMsBeforeMassCuts = vHNMs.size(); for (unsigned int iHNM = 0; iHNM < vHNMs.size(); iHNM++) { auto heavyNeutralMeson = vHNMs.at(iHNM); float massHNM = heavyNeutralMeson.m(cfgHNMMassCorrection); - if (heavyNeutralMeson.gg->reconstructionType == ReconstructionType::kPCM) { + if (heavyNeutralMeson.gg->reconstructionType == photonpair::kPCMPCM) { mHistManager.fill(HIST("HeavyNeutralMeson/invMassVsPt_PCM"), massHNM, heavyNeutralMeson.pT()); - } else if (heavyNeutralMeson.gg->reconstructionType == ReconstructionType::kEMC) { + } else if (heavyNeutralMeson.gg->reconstructionType == photonpair::kEMCEMC) { mHistManager.fill(HIST("HeavyNeutralMeson/invMassVsPt_EMC"), massHNM, heavyNeutralMeson.pT()); } else { mHistManager.fill(HIST("HeavyNeutralMeson/invMassVsPt_PCMEMC"), massHNM, heavyNeutralMeson.pT()); } if (heavyNeutralMeson.gg->isPi0 && massHNM > massWindowOmega->get("omega_min") && massHNM < massWindowOmega->get("omega_max")) { - if (heavyNeutralMeson.gg->reconstructionType == ReconstructionType::kPCM) + if (heavyNeutralMeson.gg->reconstructionType == photonpair::kPCMPCM && heavyNeutralMeson.pT() > minHNMPts->get("PCM_omega")) colContainsPCMOmega = true; - else if (heavyNeutralMeson.gg->reconstructionType == ReconstructionType::kEMC) + else if (heavyNeutralMeson.gg->reconstructionType == photonpair::kEMCEMC && heavyNeutralMeson.pT() > minHNMPts->get("EMC_omega")) colContainsEMCOmega = true; } else if (heavyNeutralMeson.gg->isEta && massHNM > massWindowEtaPrime->get("etaprime_min") && massHNM < massWindowEtaPrime->get("etaprime_max")) { - if (heavyNeutralMeson.gg->reconstructionType == ReconstructionType::kPCM) + if (heavyNeutralMeson.gg->reconstructionType == photonpair::kPCMPCM && heavyNeutralMeson.pT() > minHNMPts->get("PCM_etaprime")) colContainsPCMEtaPrime = true; - else if (heavyNeutralMeson.gg->reconstructionType == ReconstructionType::kEMC) + else if (heavyNeutralMeson.gg->reconstructionType == photonpair::kEMCEMC && heavyNeutralMeson.pT() > minHNMPts->get("EMC_etaprime")) colContainsEMCEtaPrime = true; } else { vHNMs.erase(vHNMs.begin() + iHNM); @@ -162,13 +165,13 @@ struct HeavyNeutralMesonFilter { mHistManager.fill(HIST("Event/nHeavyNeutralMesons"), nHNMsBeforeMassCuts, vHNMs.size()); if (colContainsPCMOmega) - mHistManager.fill(HIST("Event/hCollisionCounter"), 4.); + mHistManager.fill(HIST("Event/hCollisionCounter"), 2.); if (colContainsEMCOmega) - mHistManager.fill(HIST("Event/hCollisionCounter"), 5.); + mHistManager.fill(HIST("Event/hCollisionCounter"), 3.); if (colContainsPCMEtaPrime) - mHistManager.fill(HIST("Event/hCollisionCounter"), 6.); + mHistManager.fill(HIST("Event/hCollisionCounter"), 4.); if (colContainsEMCEtaPrime) - mHistManager.fill(HIST("Event/hCollisionCounter"), 7.); + mHistManager.fill(HIST("Event/hCollisionCounter"), 5.); } }; diff --git a/PWGEM/PhotonMeson/Utils/HNMUtilities.h b/PWGEM/PhotonMeson/Utils/HNMUtilities.h index 7dea4725f08..279b7a0b683 100644 --- a/PWGEM/PhotonMeson/Utils/HNMUtilities.h +++ b/PWGEM/PhotonMeson/Utils/HNMUtilities.h @@ -31,11 +31,15 @@ #include "EMCALBase/Geometry.h" #include "PWGJE/DataModel/EMCALClusters.h" #include "PWGEM/PhotonMeson/DataModel/gammaTables.h" +#include "PWGEM/PhotonMeson/Utils/PairUtilities.h" #include "Common/Core/trackUtilities.h" #include "Common/DataModel/TrackSelectionTables.h" #include "EventFiltering/filterTables.h" -namespace o2::aod::pwgem::photonmeson::utils::hnmutilities +using namespace o2::aod::pwgem::photonmeson; + +// -------> Struct to store photons from EMC clusters or V0s +namespace o2::aod::pwgem::photonmeson::hnmutilities { struct Photon { Photon(float px, float py, float pz, bool isFromConversion) : px(px), py(py), pz(pz), pt(std::sqrt(px * px + py * py)), isFromConversion(isFromConversion) @@ -64,21 +68,7 @@ struct Photon { bool isFromConversion; }; -namespace ReconstructionType -{ -enum ReconstructionType { - kPCM = 0, - kPCMEMC = 1, - kEMC = 2 -}; -} // namespace ReconstructionType - -enum MassCorrectionType { - kNoHNMMassCorrection = 0, - kSubDeltaPi0 = 1, - kSubLambda = 2 -}; - +// -------> Struct to store gamma gamma pairs (pi0 or eta meson candidates) struct GammaGammaPair { GammaGammaPair(Photon p1, Photon p2) : p1(p1), p2(p2) { @@ -96,6 +86,13 @@ struct GammaGammaPair { float pT() const { return vGG.Pt(); } }; +// -------> Enum to specify how the heavy neutral meson mass should be corrected based on the PDG mass of its light neutral meson decay daughter +enum MassCorrectionType { + kNoHNMMassCorrection = 0, + kSubDeltaPi0 = 1, + kSubLambda = 2 +}; + struct HeavyNeutralMeson { HeavyNeutralMeson(GammaGammaPair* gg, float eTracks, float pxTracks, float pyTracks, float pzTracks) : gg(gg) { @@ -109,11 +106,13 @@ struct HeavyNeutralMeson { { float massHNM = vHeavyNeutralMeson.M(); switch (massCorrectionType) { - case kSubDeltaPi0: + case kNoHNMMassCorrection: // No mass correction + break; + case kSubDeltaPi0: // Subtract the mass difference of the reconstructed light neutral meson mass to the PDG mass massHNM -= gg->m(); massHNM += (gg->isPi0 ? constants::physics::MassPi0 : 0.547862); break; - case kSubLambda: + case kSubLambda: // Subtract an opening angle dependent mass correction (see https://arxiv.org/abs/2502.19956 for details) LOGF(warning, "SubLambdaMassCorrection not yet implemented!"); break; default: @@ -124,7 +123,7 @@ struct HeavyNeutralMeson { float pT() const { return vHeavyNeutralMeson.Pt(); } }; -/// \brief Process light neutral meson candidates (pi0 or eta), calculate invariant mass and pT and fill histograms +/// \brief Reconstruct light neutral mesons from EMC clusters and V0s and fill them into the vGGs vector template void reconstructGGs(C clusters, V v0s, std::vector& vGGs) { @@ -141,17 +140,18 @@ void reconstructGGs(C clusters, V v0s, std::vector& vGGs) for (unsigned int ig2 = ig1 + 1; ig2 < vPhotons.size(); ++ig2) { GammaGammaPair lightMeson(vPhotons[ig1], vPhotons[ig2]); // build lightMeson from photons if (vPhotons[ig1].isFromConversion && vPhotons[ig2].isFromConversion) - lightMeson.setReconstructionType(ReconstructionType::kPCM); + lightMeson.setReconstructionType(photonpair::kPCMPCM); else if (!vPhotons[ig1].isFromConversion && !vPhotons[ig2].isFromConversion) - lightMeson.setReconstructionType(ReconstructionType::kEMC); + lightMeson.setReconstructionType(photonpair::kEMCEMC); else - lightMeson.setReconstructionType(ReconstructionType::kPCMEMC); + lightMeson.setReconstructionType(photonpair::kPCMEMC); vGGs.push_back(lightMeson); } } } +/// \brief Reconstruct heavy neutral mesons from tracks and GG candidates and fill them into the vHNMs vector template void reconstructHeavyNeutralMesons(Track const& tracks, std::vector& vGGs, std::vector& vHNMs) { @@ -170,6 +170,6 @@ void reconstructHeavyNeutralMesons(Track const& tracks, std::vector