diff --git a/EventFiltering/CMakeLists.txt b/EventFiltering/CMakeLists.txt index 31cd58f40ed..af42d63b16c 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 O2Physics::PWGEMPhotonMesonCore + 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..1446943c279 --- /dev/null +++ b/EventFiltering/PWGEM/HeavyNeutralMesonFilter.cxx @@ -0,0 +1,181 @@ +// 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 and V0s (PCM) +/// +/// \author Nicolas Strangmann (nicolas.strangmann@cern.ch) - Goethe University Frankfurt +/// + +#include + +#include "PWGEM/PhotonMeson/Utils/HNMUtilities.h" + +using namespace o2; +using namespace o2::framework; +using namespace o2::framework::expressions; +using namespace o2::aod::pwgem::photonmeson; + +using MyBCs = soa::Join; +using MyCollisions = soa::Join; + +using SelectedTracks = soa::Filtered>; + +struct HeavyNeutralMesonFilter { + Produces tags; + + HistogramRegistry mHistManager{"HeavyNeutralMesonFilterHistograms", {}, OutputObjHandlingPolicy::QAObject}; + + 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] = {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; + + bool colContainsPCMOmega, colContainsEMCOmega, colContainsPCMEtaPrime, colContainsEMCEtaPrime = 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, {{6, -0.5, 5.5}}); + hCollisionCounter->GetXaxis()->SetBinLabel(1, "all"); + hCollisionCounter->GetXaxis()->SetBinLabel(2, "kTVXinEMC"); + 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;#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;#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; + 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.); + + if (collision.foundBC_as().alias_bit(kTVXinEMC)) + mHistManager.fill(HIST("Event/hCollisionCounter"), 1.); + + auto v0sInThisCollision = v0s.sliceBy(perCollision_pcm, collision.globalIndex()); + auto clustersInThisCollision = clusters.sliceBy(perCollision_emc, collision.globalIndex()); + + mHistManager.fill(HIST("Event/nClustersVsV0s"), clustersInThisCollision.size(), v0sInThisCollision.size()); + mHistManager.fill(HIST("Event/nTracks"), tracks.size()); + + colContainsPCMOmega = colContainsEMCOmega = colContainsPCMEtaPrime = colContainsEMCEtaPrime = false; + + hnmutilities::reconstructGGs(clustersInThisCollision, v0sInThisCollision, vGGs); + processGGs(vGGs); + hnmutilities::reconstructHeavyNeutralMesons(tracks, vGGs, vHNMs); + processHNMs(vHNMs); + + tags(colContainsPCMOmega, colContainsEMCOmega, colContainsPCMEtaPrime, colContainsEMCEtaPrime); + } + + /// \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 == photonpair::kPCMPCM) { + mHistManager.fill(HIST("GG/invMassVsPt_PCM"), lightMeson->m(), lightMeson->pT()); + } 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()); + } + + 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"), nGGsBeforeMassCuts, vGGs.size()); + } + + /// \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 == photonpair::kPCMPCM) { + mHistManager.fill(HIST("HeavyNeutralMeson/invMassVsPt_PCM"), massHNM, heavyNeutralMeson.pT()); + } 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 == photonpair::kPCMPCM && heavyNeutralMeson.pT() > minHNMPts->get("PCM_omega")) + colContainsPCMOmega = true; + 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 == photonpair::kPCMPCM && heavyNeutralMeson.pT() > minHNMPts->get("PCM_etaprime")) + colContainsPCMEtaPrime = true; + else if (heavyNeutralMeson.gg->reconstructionType == photonpair::kEMCEMC && heavyNeutralMeson.pT() > minHNMPts->get("EMC_etaprime")) + colContainsEMCEtaPrime = true; + } else { + vHNMs.erase(vHNMs.begin() + iHNM); + iHNM--; + } + } + mHistManager.fill(HIST("Event/nHeavyNeutralMesons"), nHNMsBeforeMassCuts, vHNMs.size()); + + if (colContainsPCMOmega) + mHistManager.fill(HIST("Event/hCollisionCounter"), 2.); + if (colContainsEMCOmega) + mHistManager.fill(HIST("Event/hCollisionCounter"), 3.); + if (colContainsPCMEtaPrime) + mHistManager.fill(HIST("Event/hCollisionCounter"), 4.); + if (colContainsEMCEtaPrime) + mHistManager.fill(HIST("Event/hCollisionCounter"), 5.); + } +}; + +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..5c450b02eb2 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 @@ -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(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 @@ -298,6 +304,13 @@ DECLARE_SOA_TABLE(PhotonFilters, "AOD", "PhotonFilters", //! using PhotonFilter = PhotonFilters::iterator; +// heavy mesons +DECLARE_SOA_TABLE(HeavyNeutralMesonFilters, "AOD", "HeavyNeutralMesonFilters", //! + filtering::PCMOmegaMeson, filtering::EMCOmegaMeson, + filtering::PCMEtaPrimeMeson, filtering::EMCEtaPrimeMeson); + +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 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..279b7a0b683 --- /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 "PWGEM/PhotonMeson/Utils/PairUtilities.h" +#include "Common/Core/trackUtilities.h" +#include "Common/DataModel/TrackSelectionTables.h" +#include "EventFiltering/filterTables.h" + +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) + { + 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; +}; + +// -------> Struct to store gamma gamma pairs (pi0 or eta meson candidates) +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(); } +}; + +// -------> 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) + { + 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 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: // Subtract an opening angle dependent mass correction (see https://arxiv.org/abs/2502.19956 for details) + 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 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) +{ + 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(photonpair::kPCMPCM); + else if (!vPhotons[ig1].isFromConversion && !vPhotons[ig2].isFromConversion) + lightMeson.setReconstructionType(photonpair::kEMCEMC); + else + 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) +{ + 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 (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); + } + } + } +} + +} // namespace o2::aod::pwgem::photonmeson::hnmutilities + +#endif // PWGEM_PHOTONMESON_UTILS_HNMUTILITIES_H_