diff --git a/PWGEM/PhotonMeson/Core/Pi0EtaToGammaGamma.h b/PWGEM/PhotonMeson/Core/Pi0EtaToGammaGamma.h index 91818c3f3d2..85c74cc14bb 100644 --- a/PWGEM/PhotonMeson/Core/Pi0EtaToGammaGamma.h +++ b/PWGEM/PhotonMeson/Core/Pi0EtaToGammaGamma.h @@ -428,9 +428,9 @@ struct Pi0EtaToGammaGamma { std::vector epBinsForHs(ep_bin_edges.begin(), ep_bin_edges.end()); std::vector vtxBinsForHs(zvtx_bin_edges.begin(), zvtx_bin_edges.end()); if constexpr (pairtype == o2::aod::pwgem::photonmeson::photonpair::PairType::kPCMDalitzEE) { - o2::aod::pwgem::photonmeson::utils::nmhistogram::addNMHistograms(&fRegistry, false, cfgSparseFullAxes.value, occBinsForHs, epBinsForHs, vtxBinsForHs, "ee#gamma"); + o2::aod::pwgem::photonmeson::utils::nmhistogram::addNMHistograms(&fRegistry, false, "ee#gamma", cfgSparseFullAxes.value, occBinsForHs, epBinsForHs, vtxBinsForHs); } else { - o2::aod::pwgem::photonmeson::utils::nmhistogram::addNMHistograms(&fRegistry, false, cfgSparseFullAxes.value, occBinsForHs, epBinsForHs, vtxBinsForHs, "#gamma#gamma"); + o2::aod::pwgem::photonmeson::utils::nmhistogram::addNMHistograms(&fRegistry, false, "#gamma#gamma", cfgSparseFullAxes.value, occBinsForHs, epBinsForHs, vtxBinsForHs); } DefineEMEventCut(); DefinePCMCut(); diff --git a/PWGEM/PhotonMeson/Core/Pi0EtaToGammaGammaMC.h b/PWGEM/PhotonMeson/Core/Pi0EtaToGammaGammaMC.h index dc20e1b3883..351bb86db35 100644 --- a/PWGEM/PhotonMeson/Core/Pi0EtaToGammaGammaMC.h +++ b/PWGEM/PhotonMeson/Core/Pi0EtaToGammaGammaMC.h @@ -637,7 +637,7 @@ struct Pi0EtaToGammaGammaMC { fRegistry.fill(HIST("Event/after/hCollisionCounter"), 12.0, weight); // accepted int photonid1 = -1, photonid2 = -1, pi0id = -1, etaid = -1; - if constexpr (pairtype == o2::aod::pwgem::photonmeson::photonpair::PairType::kPCMPCM || pairtype == o2::aod::pwgem::photonmeson::photonpair::PairType::kPHOSPHOS || pairtype == o2::aod::pwgem::photonmeson::photonpair::PairType::kEMCEMC) { // same kinds pairing + if constexpr (pairtype == o2::aod::pwgem::photonmeson::photonpair::PairType::kPCMPCM || pairtype == o2::aod::pwgem::photonmeson::photonpair::PairType::kPHOSPHOS || pairtype == o2::aod::pwgem::photonmeson::photonpair::PairType::kEMCEMC || pairtype == o2::aod::pwgem::photonmeson::photonpair::PairType::kPCMEMC) { // same kinds pairing and PCM-EMC auto photons1_per_collision = photons1.sliceByCached(TDetectorTag1::perCollision(), collision.globalIndex(), cache); auto photons2_per_collision = photons2.sliceByCached(TDetectorTag2::perCollision(), collision.globalIndex(), cache); @@ -652,7 +652,7 @@ struct Pi0EtaToGammaGammaMC { continue; } } else { - if (!TDetectorTag1::applyCut(*this, g1) || !TDetectorTag2::applyCut(*this, g2)) { + if (!TDetectorTag1::applyCut(*this, g1)) { // g2 is checked below with its own tag (and matched tracks for EMC) continue; } } @@ -696,6 +696,20 @@ struct Pi0EtaToGammaGammaMC { photonid1 = o2::aod::pwgem::photonmeson::utils::mcutil::FindMotherInChain(cluster1mcparticle, mcparticles, std::vector{111, 221}); photonid2 = o2::aod::pwgem::photonmeson::utils::mcutil::FindMotherInChain(cluster2mcparticle, mcparticles, std::vector{111, 221}); + } else if constexpr (pairtype == o2::aod::pwgem::photonmeson::photonpair::PairType::kPCMEMC) { + // PCM: true conversion as for PCM-PCM + auto pos1 = g1.template posTrack_as(); + auto ele1 = g1.template negTrack_as(); + auto pos1mc = pos1.template emmcparticle_as(); + auto ele1mc = ele1.template emmcparticle_as(); + photonid1 = o2::aod::pwgem::dilepton::utils::mcutil::FindCommonMotherFrom2Prongs(pos1mc, ele1mc, -11, 11, 22, mcparticles); + + // EMC: photon below the pi0/eta of the leading contributor as for EMC-EMC + if (g2.emmcparticleIds().empty()) { + continue; // noise cluster without MC label + } + auto cluster2mcparticle = mcparticles.iteratorAt(g2.emmcparticleIds()[0]); + photonid2 = o2::aod::pwgem::photonmeson::utils::mcutil::FindMotherInChain(cluster2mcparticle, mcparticles, std::vector{111, 221}); } else { photonid1 = -1; photonid2 = -1; @@ -711,6 +725,10 @@ struct Pi0EtaToGammaGammaMC { if (!o2::aod::pwgem::photonmeson::utils::mcutil::IsConversionPointInAcceptance(g1mc, maxRgen, maxY_gen, margin_z_mc, mcparticles) || !o2::aod::pwgem::photonmeson::utils::mcutil::IsConversionPointInAcceptance(g2mc, maxRgen, maxY_gen, margin_z_mc, mcparticles)) { continue; } + } else if constexpr (pairtype == o2::aod::pwgem::photonmeson::photonpair::PairType::kPCMEMC) { + if (!o2::aod::pwgem::photonmeson::utils::mcutil::IsConversionPointInAcceptance(g1mc, maxRgen, maxY_gen, margin_z_mc, mcparticles)) { + continue; + } } pi0id = o2::aod::pwgem::dilepton::utils::mcutil::FindCommonMotherFrom2Prongs(g1mc, g2mc, 22, 22, 111, mcparticles); @@ -1024,19 +1042,17 @@ struct Pi0EtaToGammaGammaMC { auto&& [emcclusters, emcMatchedTracks, emcMatchedSecondaries] = std::forward_as_tuple(args...); runTruePairing(collisions, emcclusters, emcclusters, mccollisions, mcparticles, nullptr, emcMatchedTracks, emcMatchedSecondaries); runGenInfo(collisions, mccollisions, mcparticles); + } else if constexpr (pairtype == o2::aod::pwgem::photonmeson::photonpair::PairType::kPCMEMC) { + auto&& [v0photons, emcclusters, v0legs, emcMatchedTracks, emcMatchedSecondaries] = std::forward_as_tuple(args...); + runTruePairing(collisions, v0photons, emcclusters, mccollisions, mcparticles, v0legs, emcMatchedTracks, emcMatchedSecondaries); + runGenInfo(collisions, mccollisions, mcparticles); } // else if constexpr (pairtype == o2::aod::pwgem::photonmeson::photonpair::PairType::kPHOSPHOS) { // auto phosclusters = std::get<0>(std::tie(args...)); // runPairing(collisions, phosclusters, phosclusters, nullptr, nullptr, perCollision_phos, perCollision_phos, fPHOSCut, fPHOSCut, nullptr, nullptr); // } - // else if constexpr (pairtype == o2::aod::pwgem::photonmeson::photonpair::PairType::kPCMEMC) { - // auto v0photons = std::get<0>(std::tie(args...)); - // auto v0legs = std::get<1>(std::tie(args...)); - // auto emcclusters = std::get<2>(std::tie(args...)); - // auto emcmatchedtracks = std::get<3>(std::tie(args...)); - // runPairing(collisions, v0photons, emcclusters, v0legs, nullptr, perCollision_pcm, perCollision_emc, fV0PhotonCut, fEMCCut, emcmatchedtracks, nullptr); - // } else if constexpr (pairtype == o2::aod::pwgem::photonmeson::photonpair::PairType::kPCMPHOS) { + // else if constexpr (pairtype == o2::aod::pwgem::photonmeson::photonpair::PairType::kPCMPHOS) { // auto v0photons = std::get<0>(std::tie(args...)); // auto v0legs = std::get<1>(std::tie(args...)); // auto phosclusters = std::get<2>(std::tie(args...)); @@ -1070,19 +1086,17 @@ struct Pi0EtaToGammaGammaMC { // runTruePairing(collisions, emcclusters, emcclusters, nullptr, nullptr, perCollision_emc, perCollision_emc, fEMCCut, fEMCCut, mccollisions, mcparticles); runTruePairing(collisions, emcclusters, emcclusters, mccollisions, mcparticles, nullptr, emcMatchedTracks, emcMatchedSecondaries); runGenInfo(collisions, mccollisions, mcparticles); + } else if constexpr (pairtype == o2::aod::pwgem::photonmeson::photonpair::PairType::kPCMEMC) { + auto&& [v0photons, emcclusters, v0legs, emcMatchedTracks, emcMatchedSecondaries] = std::forward_as_tuple(args...); + runTruePairing(collisions, v0photons, emcclusters, mccollisions, mcparticles, v0legs, emcMatchedTracks, emcMatchedSecondaries); + runGenInfo(collisions, mccollisions, mcparticles); } // else if constexpr (pairtype == o2::aod::pwgem::photonmeson::photonpair::PairType::kPHOSPHOS) { // auto phosclusters = std::get<0>(std::tie(args...)); // runPairing(collisions, phosclusters, phosclusters, nullptr, nullptr, perCollision_phos, perCollision_phos, fPHOSCut, fPHOSCut, nullptr, nullptr); // } - // else if constexpr (pairtype == o2::aod::pwgem::photonmeson::photonpair::PairType::kPCMEMC) { - // auto v0photons = std::get<0>(std::tie(args...)); - // auto v0legs = std::get<1>(std::tie(args...)); - // auto emcclusters = std::get<2>(std::tie(args...)); - // auto emcmatchedtracks = std::get<3>(std::tie(args...)); - // runPairing(collisions, v0photons, emcclusters, v0legs, nullptr, perCollision_pcm, perCollision_emc, fV0PhotonCut, fEMCCut, emcmatchedtracks, nullptr); - // } else if constexpr (pairtype == o2::aod::pwgem::photonmeson::photonpair::PairType::kPCMPHOS) { + // else if constexpr (pairtype == o2::aod::pwgem::photonmeson::photonpair::PairType::kPCMPHOS) { // auto v0photons = std::get<0>(std::tie(args...)); // auto v0legs = std::get<1>(std::tie(args...)); // auto phosclusters = std::get<2>(std::tie(args...)); diff --git a/PWGEM/PhotonMeson/Tasks/CMakeLists.txt b/PWGEM/PhotonMeson/Tasks/CMakeLists.txt index 729346bab9e..56c32ab6a3b 100644 --- a/PWGEM/PhotonMeson/Tasks/CMakeLists.txt +++ b/PWGEM/PhotonMeson/Tasks/CMakeLists.txt @@ -116,6 +116,11 @@ o2physics_add_dpl_workflow(pi0eta-to-gammagamma-mc-emcemc PUBLIC_LINK_LIBRARIES O2::Framework O2::EMCALBase O2::EMCALCalib O2Physics::AnalysisCore O2Physics::PWGEMPhotonMesonCore COMPONENT_NAME Analysis) +o2physics_add_dpl_workflow(pi0eta-to-gammagamma-mc-pcmemc + SOURCES Pi0EtaToGammaGammaMCPCMEMC.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2::EMCALBase O2::EMCALCalib O2Physics::AnalysisCore O2Physics::PWGEMPhotonMesonCore + COMPONENT_NAME Analysis) + o2physics_add_dpl_workflow(tag-and-probe SOURCES TagAndProbe.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2::EMCALBase O2::EMCALCalib O2Physics::AnalysisCore O2Physics::PWGEMPhotonMesonCore O2Physics::MLCore @@ -205,3 +210,7 @@ o2physics_add_dpl_workflow(pair-qc-task SOURCES pairQcTask.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::PWGEMPhotonMesonCore COMPONENT_NAME Analysis) +o2physics_add_dpl_workflow(two-photon-mc-task + SOURCES twoPhotonMcTask.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::PWGEMPhotonMesonCore + COMPONENT_NAME Analysis) diff --git a/PWGEM/PhotonMeson/Tasks/Pi0EtaToGammaGammaMCPCMEMC.cxx b/PWGEM/PhotonMeson/Tasks/Pi0EtaToGammaGammaMCPCMEMC.cxx new file mode 100644 index 00000000000..6bd33481ce2 --- /dev/null +++ b/PWGEM/PhotonMeson/Tasks/Pi0EtaToGammaGammaMCPCMEMC.cxx @@ -0,0 +1,39 @@ +// 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 Pi0EtaToGammaGammaMCPCMEMC.cxx +/// \brief This code loops over photons and makes pairs for neutral mesons analyses in MC for PCM-EMC. +/// \author M. Hemmer, marvin.hemmer@cern.ch + +#include "PWGEM/PhotonMeson/Core/Pi0EtaToGammaGammaMC.h" +#include "PWGEM/PhotonMeson/DataModel/GammaTablesRedux.h" +#include "PWGEM/PhotonMeson/DataModel/gammaTables.h" +#include "PWGEM/PhotonMeson/Utils/PairUtilities.h" + +#include +#include +#include + +using namespace o2; +using namespace o2::aod; +using namespace o2::framework; +using namespace o2::aod::pwgem::photonmeson::photonpair; + +using MyV0Photons = o2::soa::Filtered>; +using MyMCV0Legs = soa::Join; +using MyEMCClusters = soa::Join; + +WorkflowSpec defineDataProcessing(ConfigContext const& context) +{ + return WorkflowSpec{ + adaptAnalysisTask>(context, TaskName{"pi0eta-to-gammagamma-mc-pcmemc"}), + }; +} diff --git a/PWGEM/PhotonMeson/Tasks/taskPi0FlowEMC.cxx b/PWGEM/PhotonMeson/Tasks/taskPi0FlowEMC.cxx index e824436c401..230ad975aec 100644 --- a/PWGEM/PhotonMeson/Tasks/taskPi0FlowEMC.cxx +++ b/PWGEM/PhotonMeson/Tasks/taskPi0FlowEMC.cxx @@ -133,7 +133,6 @@ struct TaskPi0FlowEMC { ConfigurableAxis thnConfigAxisCosDeltaPhi{"thnConfigAxisCosDeltaPhi", {8, -1., 1.}, "cos(delta phi) axis for the current event"}; ConfigurableAxis thnConfigAxisM02{"thnConfigAxisM02", {200, 0., 5.}, "M02 axis for the EMCal cluster"}; ConfigurableAxis thnConfigAxisEnergyCalib{"thnConfigAxisEnergyCalib", {200, 0., 20.}, "energy axis for the emcal clusters for the calibration process"}; - ConfigurableAxis thConfigAxisEtaPhiAngle{"thConfigAxisEtaPhiAngle", {180, -o2::constants::math::PI, o2::constants::math::PI}, "atan2(#Delta#eta,#Delta#varphi) axis (deg)"}; EMPhotonEventCut fEMEventCut; struct : ConfigurableGroup { @@ -222,6 +221,11 @@ struct TaskPi0FlowEMC { Configurable minTanThetadPhi{"minTanThetadPhi", 4., "apply min opening angle in delta theta delta phi to cut on late conversion"}; Configurable cfgEnableQA{"cfgEnableQA", false, "flag to turn QA plots on/off"}; ConfigurableAxis thConfigAxisTanThetaPhi{"thConfigAxisTanThetaPhi", {180, -90.f, 90.f}, ""}; + Configurable cfgPairBkgMeson{"cfgPairBkgMeson", 0, "meson for the mass binning of the pair background QA THn (0: pi0, 1: eta)"}; + ConfigurableAxis thnConfigAxisPairBkgMassPi0{"thnConfigAxisPairBkgMassPi0", {VARIABLE_WIDTH, 0.0f, 0.07f, 0.10f, 0.17f, 0.30f, 0.50f}, "invariant mass axis (sidebands and peak) of the pair background QA THn for the pi0"}; + ConfigurableAxis thnConfigAxisPairBkgMassEta{"thnConfigAxisPairBkgMassEta", {VARIABLE_WIDTH, 0.30f, 0.45f, 0.50f, 0.60f, 0.70f, 0.80f}, "invariant mass axis (sidebands and peak) of the pair background QA THn for the eta"}; + ConfigurableAxis thnConfigAxisPairBkgPt{"thnConfigAxisPairBkgPt", {VARIABLE_WIDTH, 0.f, 1.f, 2.f, 3.f, 5.f, 8.f, 20.f}, "pT axis of the pair background QA THn"}; + ConfigurableAxis thnConfigAxisPairBkgCent{"thnConfigAxisPairBkgCent", {VARIABLE_WIDTH, 0.f, 10.f, 30.f, 50.f, 90.f}, "centrality axis of the pair background QA THn"}; } mesonConfig; struct : ConfigurableGroup { @@ -417,7 +421,6 @@ struct TaskPi0FlowEMC { const AxisSpec thnAxisMixingVtx{mixingConfig.cfgVtxBins, "#it{z} (cm)"}; const AxisSpec thnAxisMixingCent{mixingConfig.cfgCentBins, "Centrality (%)"}; const AxisSpec thnAxisMixingEP{mixingConfig.cfgEPBins, Form("cos(%d#varphi)", harmonic.value)}; - const AxisSpec thAxisEtaPhiAngle{thConfigAxisEtaPhiAngle, "atan2(#Delta#eta,#Delta#varphi) (rad)"}; const AxisSpec thnQaAxisMinv{16, 0, 0.8, "#it{M}_{#gamma#gamma} (GeV/#it{c}^{2})"}; // bin witdh 0.05 GeV/c^2 const AxisSpec thnQaAxisPt{10, 0, 20., "#it{p}_{T} (GeV/#it{c})"}; // bin width 2 GeV/c @@ -477,8 +480,14 @@ struct TaskPi0FlowEMC { registry.add("mesonQA/hInvMassPtMixed", "Histo for inv pair mass vs pt for mixed event", HistType::kTH2D, {thnAxisInvMass, thnAxisPt}); registry.add("mesonQA/hTanThetaPhiMixed", "Histo for identification of conversion cluster for mixed event", HistType::kTH2D, {thnAxisInvMass, thAxisTanThetaPhi}); registry.add("mesonQA/hAlphaPtMixed", "Histo of meson asymmetry vs pT for mixed event", HistType::kTH2D, {thAxisAlpha, thnAxisPt}); - registry.add("mesonQA/hEtaPhiAngleMassCent", "atan2(#Delta#eta,#Delta#varphi) vs m_{inv} vs cent, same event", HistType::kTH3D, {thAxisEtaPhiAngle, thnAxisInvMass, thnAxisCent}); - registry.add("mesonQA/hEtaPhiAngleMassCentMixed", "atan2(#Delta#eta,#Delta#varphi) vs m_{inv} vs cent, mixed event", HistType::kTH3D, {thAxisEtaPhiAngle, thnAxisInvMass, thnAxisCent}); + const bool isPCMEMC = doprocessEMCalPCMC.value || doprocessEMCalPCMMixed.value; + const AxisSpec thnAxisPairBkgMass{mesonConfig.cfgPairBkgMeson.value == 0 ? mesonConfig.thnConfigAxisPairBkgMassPi0 : mesonConfig.thnConfigAxisPairBkgMassEta, "#it{M}_{#gamma#gamma} (GeV/#it{c}^{2})"}; + const AxisSpec thnAxisPairBkgPt{mesonConfig.thnConfigAxisPairBkgPt, "#it{p}_{T} (GeV/#it{c})"}; + const AxisSpec thnAxisPairBkgCent{mesonConfig.thnConfigAxisPairBkgCent, "Centrality (%)"}; + const AxisSpec thnAxisPairBkgAlpha = isPCMEMC ? AxisSpec{40, -1.f, 1.f, "#alpha = (#it{E}_{PCM} - #it{E}_{EMC}) / (#it{E}_{PCM} + #it{E}_{EMC})"} : AxisSpec{20, 0.f, 1.f, "#alpha = |#it{E}_{1} - #it{E}_{2}| / (#it{E}_{1} + #it{E}_{2})"}; + const AxisSpec thnAxisPairBkgDirection{36, -o2::constants::math::PIHalf, o2::constants::math::PIHalf, "atan2(#Delta#eta, #Delta#varphi) (rad)"}; + registry.add("mesonQA/hPairBkgSE", "#it{m}_{inv} vs #it{p}_{T} vs cent vs #alpha vs pair direction for same event", HistType::kTHnF, {thnAxisPairBkgMass, thnAxisPairBkgPt, thnAxisPairBkgCent, thnAxisPairBkgAlpha, thnAxisPairBkgDirection}); + registry.add("mesonQA/hPairBkgME", "#it{m}_{inv} vs #it{p}_{T} vs cent vs #alpha vs pair direction for mixed event", HistType::kTHnF, {thnAxisPairBkgMass, thnAxisPairBkgPt, thnAxisPairBkgCent, thnAxisPairBkgAlpha, thnAxisPairBkgDirection}); registry.add("mesonQA/hDeltaEtaDeltaPhiMinvPtCentSE", "#it{m}_{inv} vs #it{p}_{T} vs #Delta#eta vs #Delta#varphi vs cent for same event", HistType::kTHnF, {thnQaAxisMinv, thnQaAxisPt, thnQaAxisDeltaEta, thnQaAxisDeltaPhi, thnQaAxisCent}); registry.add("mesonQA/hDeltaEtaDeltaPhiMinvPtCentME", "#it{m}_{inv} vs #it{p}_{T} vs #Delta#eta vs #Delta#varphi vs cent for mixed event", HistType::kTHnF, {thnQaAxisMinv, thnQaAxisPt, thnQaAxisDeltaEta, thnQaAxisDeltaPhi, thnQaAxisCent}); } @@ -541,6 +550,38 @@ struct TaskPi0FlowEMC { registry.fill(HIST(HistTypes[histType]), mass, pt, cent); } + /// \brief Fill the pair background QA THn (m_inv, pT, cent, alpha, direction) + /// The direction is folded such that dPhi is in [0, pi] and dEta keeps its sign, identical for same and mixed event. + /// \tparam isMixed fill the mixed instead of the same event histogram + /// \tparam isSignedAlpha use the signed asymmetry (E1 - E2) / (E1 + E2), for PCM-EMC photon1 has to be the PCM photon + /// \param photon1 first photon + /// \param photon2 second photon + /// \param meson pair of photon1 and photon2 + /// \param cent centrality of the (first) collision + template + void fillPairBkg(ROOT::Math::PtEtaPhiMVector const& photon1, ROOT::Math::PtEtaPhiMVector const& photon2, ROOT::Math::PtEtaPhiMVector const& meson, const float cent) + { + if (!mesonConfig.cfgEnableQA.value) { + return; + } + float alpha = (photon1.E() - photon2.E()) / (photon1.E() + photon2.E()); + if constexpr (!isSignedAlpha) { + alpha = std::fabs(alpha); + } + float dPhi = RecoDecay::constrainAngle(photon1.Phi() - photon2.Phi(), -o2::constants::math::PI); // wrap to [-pi, pi) + float dEta = photon1.Eta() - photon2.Eta(); + if (dPhi < 0.f) { // swap the labels: flip BOTH signs + dPhi = -dPhi; + dEta = -dEta; + } + const float direction = std::atan2(dEta, dPhi); // in [-pi/2, pi/2] since dPhi >= 0 + if constexpr (isMixed) { + registry.fill(HIST("mesonQA/hPairBkgME"), meson.M(), meson.Pt(), cent, alpha, direction); + } else { + registry.fill(HIST("mesonQA/hPairBkgSE"), meson.M(), meson.Pt(), cent, alpha, direction); + } + } + /// \brief eta a cluster would have if its own vertex vzOld is swapped for vzNew, /// assuming a nominal cylindrical EMCal surface at transverse radius emcalR (cm). /// phi is untouched: transverse vertex spread is negligible next to emcalR. @@ -1077,7 +1118,12 @@ struct TaskPi0FlowEMC { float dTheta = v1.Theta() - v2.Theta(); float dPhi = v1.Phi() - v2.Phi(); + float dPhiWrapped = RecoDecay::constrainAngle(v1.Phi() - v2.Phi(), -o2::constants::math::PI); // wrap to [-π, π) float dEta = v1.Eta() - v2.Eta(); + if (dPhiWrapped < 0.) { // swap the labels: flip BOTH signs + dPhiWrapped = -dPhiWrapped; + dEta = -dEta; + } float openingAngle = std::acos(v1.Vect().Dot(v2.Vect()) / (v1.P() * v2.P())); registry.fill(HIST("hMesonCuts"), 1); @@ -1100,14 +1146,14 @@ struct TaskPi0FlowEMC { registry.fill(HIST("mesonQA/hInvMassPt"), vMeson.M(), vMeson.Pt()); registry.fill(HIST("mesonQA/hTanThetaPhi"), vMeson.M(), getAngleDegree(std::atan(dTheta / dPhi))); registry.fill(HIST("mesonQA/hAlphaPt"), (v1.E() - v2.E()) / (v1.E() + v2.E()), vMeson.Pt()); - registry.fill(HIST("mesonQA/hEtaPhiAngleMassCent"), RecoDecay::constrainAngle(std::atan2(dEta, dPhi), -o2::constants::math::PI), vMeson.M(), getCentrality(collision)); - registry.fill(HIST("mesonQA/hDeltaEtaDeltaPhiMinvPtCentSE"), vMeson.M(), vMeson.Pt(), dEta, std::fabs(RecoDecay::constrainAngle(dPhi, -o2::constants::math::PI)), getCentrality(collision)); + registry.fill(HIST("mesonQA/hDeltaEtaDeltaPhiMinvPtCentSE"), vMeson.M(), vMeson.Pt(), dEta, dPhiWrapped, getCentrality(collision)); } if (mesonConfig.enableTanThetadPhi.value && mesonConfig.minTanThetadPhi > std::fabs(getAngleDegree(std::atan(dTheta / dPhi)))) { registry.fill(HIST("hMesonCuts"), 5); continue; } registry.fill(HIST("hMesonCuts"), 6); + fillPairBkg(v1, v2, vMeson, getCentrality(collision)); runFlowAnalysis<0>(collision, vMeson); } } @@ -1241,7 +1287,12 @@ struct TaskPi0FlowEMC { float dTheta = v1.Theta() - v2.Theta(); float dPhi = v1.Phi() - v2.Phi(); - float dEta = v1.Eta() - eta2; + float dPhiWrapped = RecoDecay::constrainAngle(v1.Phi() - v2.Phi(), -o2::constants::math::PI); // wrap to [-π, π) + float dEta = v1.Eta() - v2.Eta(); + if (dPhiWrapped < 0.) { // swap the labels: flip BOTH signs + dPhiWrapped = -dPhiWrapped; + dEta = -dEta; + } float openingAngle = std::acos(v1.Vect().Dot(v2.Vect()) / (v1.P() * v2.P())); registry.fill(HIST("hMesonCutsMixed"), 1); @@ -1261,14 +1312,14 @@ struct TaskPi0FlowEMC { registry.fill(HIST("mesonQA/hInvMassPtMixed"), vMeson.M(), vMeson.Pt()); registry.fill(HIST("mesonQA/hTanThetaPhiMixed"), vMeson.M(), getAngleDegree(std::atan(dTheta / dPhi))); registry.fill(HIST("mesonQA/hAlphaPtMixed"), (v1.E() - v2.E()) / (v1.E() + v2.E()), vMeson.Pt()); - registry.fill(HIST("mesonQA/hEtaPhiAngleMassCentMixed"), RecoDecay::constrainAngle(std::atan2(dEta, dPhi), -o2::constants::math::PI), vMeson.M(), getCentrality(c1)); - registry.fill(HIST("mesonQA/hDeltaEtaDeltaPhiMinvPtCentME"), vMeson.M(), vMeson.Pt(), dEta, std::fabs(RecoDecay::constrainAngle(dPhi, -o2::constants::math::PI)), getCentrality(c1)); + registry.fill(HIST("mesonQA/hDeltaEtaDeltaPhiMinvPtCentME"), vMeson.M(), vMeson.Pt(), dEta, dPhiWrapped, getCentrality(c1)); } if (mesonConfig.enableTanThetadPhi.value && mesonConfig.minTanThetadPhi > std::fabs(getAngleDegree(std::atan(dTheta / dPhi)))) { registry.fill(HIST("hMesonCutsMixed"), 5); continue; } registry.fill(HIST("hMesonCutsMixed"), 6); + fillPairBkg(v1, v2, vMeson, getCentrality(c1)); runFlowAnalysis<2>(c1, vMeson); } } @@ -1335,7 +1386,12 @@ struct TaskPi0FlowEMC { float dTheta = v1.Theta() - v2.Theta(); float dPhi = v1.Phi() - v2.Phi(); + float dPhiWrapped = RecoDecay::constrainAngle(v1.Phi() - v2.Phi(), -o2::constants::math::PI); // wrap to [-π, π) float dEta = v1.Eta() - v2.Eta(); + if (dPhiWrapped < 0.) { // swap the labels: flip BOTH signs + dPhiWrapped = -dPhiWrapped; + dEta = -dEta; + } float openingAngle = std::acos(v1.Vect().Dot(v2.Vect()) / (v1.P() * v2.P())); registry.fill(HIST("hMesonCuts"), 1); if (openingAngle <= mesonConfig.minOpenAngle) { @@ -1357,13 +1413,13 @@ struct TaskPi0FlowEMC { registry.fill(HIST("mesonQA/hInvMassPt"), vMeson.M(), vMeson.Pt()); registry.fill(HIST("mesonQA/hTanThetaPhi"), vMeson.M(), getAngleDegree(std::atan(dTheta / dPhi))); registry.fill(HIST("mesonQA/hAlphaPt"), (v1.E() - v2.E()) / (v1.E() + v2.E()), vMeson.Pt()); - registry.fill(HIST("mesonQA/hEtaPhiAngleMassCent"), RecoDecay::constrainAngle(std::atan2(dEta, dPhi), -o2::constants::math::PI), vMeson.M(), getCentrality(collision)); - registry.fill(HIST("mesonQA/hDeltaEtaDeltaPhiMinvPtCentSE"), vMeson.M(), vMeson.Pt(), dEta, std::fabs(RecoDecay::constrainAngle(dPhi, -o2::constants::math::PI)), getCentrality(collision)); + registry.fill(HIST("mesonQA/hDeltaEtaDeltaPhiMinvPtCentSE"), vMeson.M(), vMeson.Pt(), dEta, dPhiWrapped, getCentrality(collision)); } if (mesonConfig.enableTanThetadPhi.value && mesonConfig.minTanThetadPhi > std::fabs(getAngleDegree(std::atan(dTheta / dPhi)))) { registry.fill(HIST("hMesonCuts"), 5); continue; } + fillPairBkg(v2, v1, vMeson, getCentrality(collision)); // v2 is the PCM photon runFlowAnalysis<0>(collision, vMeson); } } @@ -1434,7 +1490,12 @@ struct TaskPi0FlowEMC { float dTheta = v1.Theta() - v2.Theta(); float dPhi = v1.Phi() - v2.Phi(); + float dPhiWrapped = RecoDecay::constrainAngle(v1.Phi() - v2.Phi(), -o2::constants::math::PI); // wrap to [-π, π) float dEta = v1.Eta() - v2.Eta(); + if (dPhiWrapped < 0.) { // swap the labels: flip BOTH signs + dPhiWrapped = -dPhiWrapped; + dEta = -dEta; + } float openingAngle = std::acos(v1.Vect().Dot(v2.Vect()) / (v1.P() * v2.P())); registry.fill(HIST("hMesonCutsMixed"), 1); @@ -1454,14 +1515,14 @@ struct TaskPi0FlowEMC { registry.fill(HIST("mesonQA/hInvMassPtMixed"), vMeson.M(), vMeson.Pt()); registry.fill(HIST("mesonQA/hTanThetaPhiMixed"), vMeson.M(), getAngleDegree(std::atan(dTheta / dPhi))); registry.fill(HIST("mesonQA/hAlphaPtMixed"), (v1.E() - v2.E()) / (v1.E() + v2.E()), vMeson.Pt()); - registry.fill(HIST("mesonQA/hEtaPhiAngleMassCentMixed"), RecoDecay::constrainAngle(std::atan2(dEta, dPhi), -o2::constants::math::PI), vMeson.M(), getCentrality(c1)); - registry.fill(HIST("mesonQA/hDeltaEtaDeltaPhiMinvPtCentME"), vMeson.M(), vMeson.Pt(), dEta, std::fabs(RecoDecay::constrainAngle(dPhi, -o2::constants::math::PI)), getCentrality(c1)); + registry.fill(HIST("mesonQA/hDeltaEtaDeltaPhiMinvPtCentME"), vMeson.M(), vMeson.Pt(), dEta, dPhiWrapped, getCentrality(c1)); } if (mesonConfig.enableTanThetadPhi.value && mesonConfig.minTanThetadPhi > std::fabs(getAngleDegree(std::atan(dTheta / dPhi)))) { registry.fill(HIST("hMesonCutsMixed"), 5); continue; } registry.fill(HIST("hMesonCutsMixed"), 6); + fillPairBkg(v2, v1, vMeson, getCentrality(c1)); // v2 is the PCM photon runFlowAnalysis<2>(c1, vMeson); } } @@ -1555,8 +1616,12 @@ struct TaskPi0FlowEMC { ROOT::Math::PtEtaPhiMVector v2(g2.corrPt(), g2.eta(), g2.phi(), 0.); ROOT::Math::PtEtaPhiMVector vMeson = v1 + v2; - float dPhi = v1.Phi() - v2.Phi(); + float dPhiWrapped = RecoDecay::constrainAngle(v1.Phi() - v2.Phi(), -o2::constants::math::PI); // wrap to [-π, π) float dEta = v1.Eta() - v2.Eta(); + if (dPhiWrapped < 0.) { // swap the labels: flip BOTH signs + dPhiWrapped = -dPhiWrapped; + dEta = -dEta; + } float openingAngle = std::acos(v1.Vect().Dot(v2.Vect()) / (v1.P() * v2.P())); registry.fill(HIST("hMesonCuts"), 1); @@ -1575,10 +1640,10 @@ struct TaskPi0FlowEMC { if (mesonConfig.cfgEnableQA.value) { registry.fill(HIST("mesonQA/hInvMassPt"), vMeson.M(), vMeson.Pt()); registry.fill(HIST("mesonQA/hAlphaPt"), (v1.E() - v2.E()) / (v1.E() + v2.E()), vMeson.Pt()); - registry.fill(HIST("mesonQA/hEtaPhiAngleMassCent"), RecoDecay::constrainAngle(std::atan2(dEta, dPhi), -o2::constants::math::PI), vMeson.M(), getCentrality(collision)); - registry.fill(HIST("mesonQA/hDeltaEtaDeltaPhiMinvPtCentSE"), vMeson.M(), vMeson.Pt(), dEta, std::fabs(RecoDecay::constrainAngle(dPhi, -o2::constants::math::PI)), getCentrality(collision)); + registry.fill(HIST("mesonQA/hDeltaEtaDeltaPhiMinvPtCentSE"), vMeson.M(), vMeson.Pt(), dEta, dPhiWrapped, getCentrality(collision)); } registry.fill(HIST("hMesonCuts"), 6); + fillPairBkg(v1, v2, vMeson, getCentrality(collision)); runFlowAnalysis<0>(collision, vMeson); } } // end of loop over collisions @@ -1630,8 +1695,12 @@ struct TaskPi0FlowEMC { ROOT::Math::PtEtaPhiMVector v2(g2.corrPt(), g2.eta(), g2.phi(), 0.); ROOT::Math::PtEtaPhiMVector vMeson = v1 + v2; - float dPhi = v1.Phi() - v2.Phi(); + float dPhiWrapped = RecoDecay::constrainAngle(v1.Phi() - v2.Phi(), -o2::constants::math::PI); // wrap to [-π, π) float dEta = v1.Eta() - v2.Eta(); + if (dPhiWrapped < 0.) { // swap the labels: flip BOTH signs + dPhiWrapped = -dPhiWrapped; + dEta = -dEta; + } float openingAngle = std::acos(v1.Vect().Dot(v2.Vect()) / (v1.P() * v2.P())); registry.fill(HIST("hMesonCutsMixed"), 1); @@ -1650,10 +1719,10 @@ struct TaskPi0FlowEMC { if (mesonConfig.cfgEnableQA.value) { registry.fill(HIST("mesonQA/hInvMassPtMixed"), vMeson.M(), vMeson.Pt()); registry.fill(HIST("mesonQA/hAlphaPtMixed"), (v1.E() - v2.E()) / (v1.E() + v2.E()), vMeson.Pt()); - registry.fill(HIST("mesonQA/hEtaPhiAngleMassCent"), RecoDecay::constrainAngle(std::atan2(dEta, dPhi), -o2::constants::math::PI), vMeson.M(), getCentrality(c1)); - registry.fill(HIST("mesonQA/hDeltaEtaDeltaPhiMinvPtCentME"), vMeson.M(), vMeson.Pt(), dEta, std::fabs(RecoDecay::constrainAngle(dPhi, -o2::constants::math::PI)), getCentrality(c1)); + registry.fill(HIST("mesonQA/hDeltaEtaDeltaPhiMinvPtCentME"), vMeson.M(), vMeson.Pt(), dEta, dPhiWrapped, getCentrality(c1)); } registry.fill(HIST("hMesonCutsMixed"), 6); + fillPairBkg(v1, v2, vMeson, getCentrality(c1)); runFlowAnalysis<2>(c1, vMeson); } } diff --git a/PWGEM/PhotonMeson/Tasks/twoPhotonMcTask.cxx b/PWGEM/PhotonMeson/Tasks/twoPhotonMcTask.cxx new file mode 100644 index 00000000000..6b2e5fd2aa7 --- /dev/null +++ b/PWGEM/PhotonMeson/Tasks/twoPhotonMcTask.cxx @@ -0,0 +1,1278 @@ +// 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 twoPhotonMcTask.cxx +/// \brief Analysis task to obtain the invariant mass distributions for different MC truth classes +/// \author M. Hemmer, marvin.hemmer@cern.ch + +#include "PWGEM/Dilepton/Utils/MCUtilities.h" +#include "PWGEM/PhotonMeson/Core/EMBitFlags.h" +#include "PWGEM/PhotonMeson/Core/EMCPhotonCut.h" +#include "PWGEM/PhotonMeson/Core/EMPhotonEventCut.h" +#include "PWGEM/PhotonMeson/Core/V0PhotonCut.h" +#include "PWGEM/PhotonMeson/DataModel/EventTables.h" +#include "PWGEM/PhotonMeson/DataModel/GammaTablesRedux.h" +#include "PWGEM/PhotonMeson/DataModel/gammaTables.h" +#include "PWGEM/PhotonMeson/Utils/EventHistograms.h" +#include "PWGEM/PhotonMeson/Utils/MCUtilities.h" + +#include "Common/DataModel/Centrality.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include // IWYU pragma: keep (do not replace with Math/Vector4Dfwd.h) +#include +#include +#include +#include +#include +#include +#include +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace o2; +using namespace o2::aod; +using namespace o2::framework; +using namespace o2::framework::expressions; +using namespace o2::soa; +using namespace o2::aod::pwgem::photon; +using namespace o2::aod::pwgem::photonmeson::utils::mcutil; + +static constexpr int kMaxChainLength = 16; +using Chain = std::array; // filled with -1 after the last ancestor + +/// \brief MC truth of a single photon candidate, ordered from best to worst such that std::max gives the worse flag of a pair +enum class PhotonTruth : uint8_t { + kTruePhoton = 0, // PCM: both legs from the same photon, EMC: leading contributor is a photon (or e+-, converted in front of the EMCal) + kNonPhoton, // PCM: both legs have the same direct mother which is not a photon (K0S, Lambda, pi0 Dalitz), EMC: leading contributor is not a photon + kFake, // PCM: leg without MC label or legs without common direct mother, EMC: cluster without MC label + kNTruth +}; + +/// \brief Class of the closest common ancestor of the two photon candidates of a pair +enum class PairBucket : uint8_t { + kNone = 0, // no common ancestor + kSameParticle, // one particle seen twice: the common ancestor is one of the two candidates themselves (EMC conversion split, Dalitz), + // or a photon/electron whose electromagnetic shower made both candidates (e.g. bremsstrahlung of both conversion legs) + kSignalPi0, // pi0 is the direct mother of both + kSignalEta, // eta is the direct mother of both + kCrossPi0, // pi0 further up in the chain: the pi0 has no hadronic decays, so these are pi0 decays with an EM interaction in between (conversion, bremsstrahlung, shower in front of the EMCal) + kCrossEta, // eta further up in the chain: hadronic decays (eta -> 3 pi0, photons of different pi0) and eta decays with an EM interaction in between + kCrossOmega, // omega + kCrossK0S, // K0S + kCrossEtaPrime, // eta' + kCrossOtherHadron, // any other hadron (including the rare Sigma0) + kPartonString, // quark, gluon, diquark or string/cluster: jet-like + kNBuckets +}; + +/// \brief What the leading contributor of the EMC cluster is in a PCM-EMC same particle pair +enum class SameParticleOrigin : uint8_t { + kLegElectron = 0, // the e+- that was reconstructed as one of the V0 legs, should have been removed by track matching + kBremsFromLeg, // a photon radiated by one of the V0 leg e+-, cannot be removed by track matching + kOtherElectron, // any other e+- + kOther, // anything else + kNOrigins +}; + +/// \brief Why the V0 leg did not veto the EMC cluster, from the tracks matched to the cluster +enum class LegMatchStatus : uint8_t { + kNoMatchStored = 0, // no matched track stored for this cluster at all + kNoCompatibleMatch, // matched tracks stored, but none with the momentum of the leg + kOutsideWindow, // leg found, but outside the track matching window of the analysis + kInsideWindowEoverP, // leg found inside the window, but E/p above the cut (primaries only) + kInsideWindow, // leg found inside the window: the cluster should have been vetoed + kNStatus +}; + +/// \brief Species of the closest common ancestor of an "other hadron" pair +enum class OtherHadronSpecies : uint8_t { + kChargedPion = 0, + kChargedKaon, + kK0L, + kRho, + kKstar, + kPhi, + kA1, // a1(1260) + kExcitedK, // K1(1270), K1(1400), K0*(1430), K2*(1430) + kOpenCharm, // D, D*, Ds, ... (one c quark) + kCharmonium, // J/psi, psi(2S), chi_c, ... + kOpenBeauty, // B, B*, Bs, Bc, ... (one b quark) + kBottomonium, // Upsilon, chi_b, ... + kOtherMeson, // incl. f0, f2, b, h resonances + kProton, + kNeutron, + kDelta, + kLambda, + kSigma, + kXi, + kOtherBaryon, + kOther, + kNSpecies +}; + +/// \brief How the particles directly below the common ancestor were produced, i.e. how the ancestor split into the two branches. +/// Ordered such that std::max of the two branches gives the material interaction if one branch comes from it. +enum class SplitProduction : uint8_t { + kGeneratorDecay = 0, // produced by the event generator: decay in the generator (resonances, Sigma0, ...) + kTransportDecay, // produced by the transport code in a decay (weak decays of K+-, K0L, Lambda, Xi, ...) + kHadronicInteraction, // produced by the transport code in a hadronic interaction in the detector material + kOtherProcess, // produced by the transport code in any other process + kNProduction +}; + +enum CentralityEstimator { + None = 0, + CFT0A = 1, + CFT0C = 2, + CFT0M = 3, + NCentralityEstimators = 4 +}; + +struct TwoPhotonMcTask { + // configurable axis + ConfigurableAxis thnConfigAxisInvMass{"thnConfigAxisInvMass", {400, 0.0, 0.8}, "invariant mass axis for the neutral meson"}; + ConfigurableAxis thnConfigAxisPt{"thnConfigAxisPt", {100, 0., 20.}, "pT axis for the neutral meson"}; + ConfigurableAxis thnConfigAxisERelative{"thnConfigAxisERelative", {600, -1., 5.}, "(E rec - E true) / E true axis"}; + ConfigurableAxis thnConfigAxisPRelative{"thnConfigAxisPRelative", {300, -1., 2.}, "(P rec - P true) / P true axis"}; + ConfigurableAxis thnConfigAxisEtaRelative{"thnConfigAxisEtaRelative", {300, -1., 2.}, "(eta rec - eta true) / eta true axis"}; + ConfigurableAxis thnConfigAxisPhiRelative{"thnConfigAxisPhiRelative", {300, -1., 2.}, "(phi rec - phi true) / phi true axis"}; + ConfigurableAxis thnConfigAxisCent{"thnConfigAxisCent", {10, 0., 100.}, "centrality axis for the current event"}; + ConfigurableAxis thnConfigAxisMult{"thnConfigAxisMult", {60, 0., 60000.}, "multiplicity axis for the current event"}; + Configurable useCent{"useCent", false, "flag to enable usage of centrality instead of multiplicity as axis."}; + + EMPhotonEventCut fEMEventCut; + struct : ConfigurableGroup { + std::string prefix = "eventcuts"; + Configurable cfgZvtxMax{"cfgZvtxMax", 10.f, "max. Zvtx"}; + Configurable cfgRequireSel8{"cfgRequireSel8", true, "require sel8 in event cut"}; + Configurable cfgRequireFT0AND{"cfgRequireFT0AND", true, "require FT0AND in event cut"}; + Configurable cfgRequireNoTFB{"cfgRequireNoTFB", false, "require No time frame border in event cut"}; + Configurable cfgRequireNoITSROFB{"cfgRequireNoITSROFB", false, "require no ITS readout frame border in event cut"}; + Configurable cfgRequireNoSameBunchPileup{"cfgRequireNoSameBunchPileup", false, "require no same bunch pileup in event cut"}; + Configurable cfgRequireVertexITSTPC{"cfgRequireVertexITSTPC", false, "require Vertex ITSTPC in event cut"}; // ITS-TPC matched track contributes PV. + Configurable cfgRequireGoodZvtxFT0vsPV{"cfgRequireGoodZvtxFT0vsPV", false, "require good Zvtx between FT0 vs. PV in event cut"}; + Configurable cfgRequireEMCReadoutInMB{"cfgRequireEMCReadoutInMB", true, "require the EMC to be read out in an MB collision (kTVXinEMC)"}; + Configurable cfgRequireEMCHardwareTriggered{"cfgRequireEMCHardwareTriggered", false, "require the EMC to be hardware triggered (kEMC7 or kDMC7)"}; + Configurable cfgFT0COccupancyMin{"cfgFT0COccupancyMin", -1, "min. FT0C occupancy"}; + Configurable cfgFT0COccupancyMax{"cfgFT0COccupancyMax", 1000000000, "max. FT0C occupancy"}; + Configurable cfgMinCent{"cfgMinCent", 0, "min. centrality (%)"}; + Configurable cfgMaxCent{"cfgMaxCent", 90, "max. centrality (%)"}; + Configurable centEstimator{"centEstimator", 2, "Centrality estimation (FT0A: 1, FT0C: 2, FT0M: 3)"}; + Configurable onlyKeepWeightedEvents{"onlyKeepWeightedEvents", false, "flag to keep only weighted events (for JJ MCs) and remove all MB events (with weight = 1)"}; + } eventcuts; + + EMCPhotonCut fEMCCut; + struct : ConfigurableGroup { + std::string prefix = "emccuts"; + Configurable clusterDefinition{"clusterDefinition", "kV3MostSplitSmallestTimeDiff", "Clusterizer to be selected, e.g. V3Default"}; + Configurable cfgEMCminTime{"cfgEMCminTime", -25., "Minimum cluster time for EMCal time cut"}; + Configurable cfgEMCmaxTime{"cfgEMCmaxTime", +30., "Maximum cluster time for EMCal time cut"}; + Configurable cfgEMCminM02{"cfgEMCminM02", 0.1, "Minimum M02 for EMCal M02 cut"}; + Configurable cfgEMCmaxM02{"cfgEMCmaxM02", 0.7, "Maximum M02 for EMCal M02 cut"}; + Configurable cfgEMCminE{"cfgEMCminE", 0.7, "Minimum cluster energy for EMCal energy cut"}; + Configurable cfgEMCminNCell{"cfgEMCminNCell", 1, "Minimum number of cells per cluster for EMCal NCell cut"}; + Configurable> cfgEMCTMEta{"cfgEMCTMEta", {0.01f, 4.07f, -2.5f}, "|eta| <= [0]+(pT+[1])^[2] for EMCal track matching"}; + Configurable> cfgEMCTMPhi{"cfgEMCTMPhi", {0.015f, 3.65f, -2.f}, "|phi| <= [0]+(pT+[1])^[2] for EMCal track matching"}; + Configurable> emcSecTMEta{"emcSecTMEta", {0.01f, 4.07f, -2.5f}, "|eta| <= [0]+(pT+[1])^[2] for EMCal track matching"}; + Configurable> emcSecTMPhi{"emcSecTMPhi", {0.015f, 3.65f, -2.f}, "|phi| <= [0]+(pT+[1])^[2] for EMCal track matching"}; + Configurable cfgEMCEoverp{"cfgEMCEoverp", 1.75, "Minimum cluster energy over track momentum for EMCal track matching"}; + Configurable cfgEMCUseExoticCut{"cfgEMCUseExoticCut", true, "FLag to use the EMCal exotic cluster cut"}; + Configurable cfgEMCUseTM{"cfgEMCUseTM", false, "flag to use EMCal track matching cut or not"}; + Configurable emcUseSecondaryTM{"emcUseSecondaryTM", false, "flag to use EMCal secondary track matching cut or not"}; + Configurable cfgEnableQA{"cfgEnableQA", false, "flag to turn QA plots on/off"}; + } emccuts; + + V0PhotonCut fV0PhotonCut; + struct : o2::framework::ConfigurableGroup { + std::string prefix = "PCMcuts"; + o2::framework::Configurable requireV0WithITSTPC{"requireV0WithITSTPC", false, "flag to enforce V0s have ITS and TPC"}; + o2::framework::Configurable requireV0WithITSOnly{"requireV0WithITSOnly", false, "flag to select V0s with ITSonly tracks"}; + o2::framework::Configurable requireV0WithTPCOnly{"requireV0WithTPCOnly", false, "flag to select V0s with TPConly tracks"}; + o2::framework::Configurable minPtV0{"minPtV0", 0.1, "min pT for v0 photons at PV"}; + o2::framework::Configurable maxPtV0{"maxPtV0", 1e+10, "max pT for v0 photons at PV"}; + o2::framework::Configurable minEtaV0{"minEtaV0", -0.8, "min eta for v0 photons at PV"}; + o2::framework::Configurable maxEtaV0{"maxEtaV0", 0.8, "max eta for v0 photons at PV"}; + o2::framework::Configurable minRV0{"minRV0", 4.0, "min v0 radius"}; + o2::framework::Configurable maxRV0{"maxRV0", 90.0, "max v0 radius"}; + o2::framework::Configurable maxAlphaAP{"maxAlphaAP", 0.95, "max alpha for AP cut"}; + o2::framework::Configurable maxQtAP{"maxQtAP", 0.01, "max qT for AP cut"}; + o2::framework::Configurable minCosPA{"minCosPA", 0.999, "min V0 CosPA"}; + o2::framework::Configurable maxPCA{"maxPCA", 1.5, "max distance btween 2 legs"}; + o2::framework::Configurable maxChi2KF{"maxChi2KF", 1.e+10f, "max chi2/ndf with KF"}; + o2::framework::Configurable rejectV0onITSib{"rejectV0onITSib", true, "flag to reject V0s on ITSib"}; + o2::framework::Configurable applyPrefilter{"applyPrefilter", false, "flag to apply prefilter to V0"}; + + o2::framework::Configurable minNClusterTPC{"minNClusterTPC", 0, "min NCluster TPC"}; + o2::framework::Configurable minNCrossedRowsTPC{"minNCrossedRowsTPC", 40, "min ncrossed rows in TPC"}; + o2::framework::Configurable minNCrossedRowsOverFindableClustersTPC{"minNCrossedRowsOverFindableClustersTPC", 0.8f, "min fraction of crossed rows over findable clusters in TPC"}; + o2::framework::Configurable maxFracSharedClustersTPC{"maxFracSharedClustersTPC", 999.f, "max fraction of shared clusters in TPC"}; + o2::framework::Configurable minNClusterITS{"minNClusterITS", 0, "min NCluster ITS"}; + o2::framework::Configurable minMeanClusterSizeITSob{"minMeanClusterSizeITSob", 0.f, "min ITSob"}; + o2::framework::Configurable maxMeanClusterSizeITSob{"maxMeanClusterSizeITSob", 16.f, "max ITSob"}; + o2::framework::Configurable macChi2TPC{"macChi2TPC", 4.f, "max chi2/NclsTPC"}; + o2::framework::Configurable macChi2ITS{"macChi2ITS", 36.f, "max chi2/NclsITS"}; + o2::framework::Configurable minTPCNSigmaEl{"minTPCNSigmaEl", -3.0f, "min. TPC n sigma for electron"}; + o2::framework::Configurable maxTPCNSigmaEl{"maxTPCNSigmaEl", +3.0f, "max. TPC n sigma for electron"}; + o2::framework::Configurable disableITSOnly{"disableITSOnly", false, "flag to disable ITSonly tracks"}; + o2::framework::Configurable disableTPCOnly{"disableTPCOnly", false, "flag to disable TPConly tracks"}; + o2::framework::Configurable doQA{"doQA", false, "flag to set QA flag."}; + } pcmcuts; + + struct : ConfigurableGroup { + std::string prefix = "mesonConfig"; + Configurable minOpenAngle{"minOpenAngle", 0.0202, "apply min opening angle. Default value one EMCal cell"}; + Configurable cfgEnableQA{"cfgEnableQA", false, "flag to turn QA plots on/off"}; + } mesonConfig; + + // same mixing as taskPi0FlowEMC, but without event plane bins: the derived MC has no Q vectors + struct : ConfigurableGroup { + std::string prefix = "mixingConfig"; + ConfigurableAxis cfgVtxBins{"cfgVtxBins", {VARIABLE_WIDTH, -10.0f, -8.f, -6.f, -4.f, -2.f, 0.f, 2.f, 4.f, 6.f, 8.f, 10.f}, "Mixing bins - z-vertex"}; + ConfigurableAxis cfgCentBins{"cfgCentBins", {VARIABLE_WIDTH, 0.0f, 5.0f, 10.0f, 20.0f, 30.0f, 40.0f, 50.0f, 60.0f, 70.0f, 80.0f, 90.0f, 100.f}, "Mixing bins - centrality"}; + Configurable cfgMixingDepth{"cfgMixingDepth", 2, "Mixing depth"}; + Configurable cfgCorrectMixedVtxEta{"cfgCorrectMixedVtxEta", true, "EMC-EMC: re-project the eta of cluster 2 onto the vertex of collision 1 in the mixed event"}; + Configurable cfgEmcalEffRadius{"cfgEmcalEffRadius", 430.f, "effective EMCal radius (cm) used for the mixed event vertex correction"}; + } mixingConfig; + + using BinningType = ColumnBinningPolicy; + + SliceCache cache; + + using EMCalPhotons = soa::Join; + using PcmPhotons = soa::Join; + + using PcmMcLegs = soa::Join; + + using Colls = soa::Join; + + using McColls = o2::soa::Join; + using McParticles = EMMCParticles; + + PresliceOptional perCollisionEMC = o2::aod::emccluster::pmeventId; + PresliceOptional perCollisionPCM = aod::v0photonkf::pmeventId; + PresliceOptional perEMCClusterMT = o2::aod::mintm::minClusterId; + PresliceOptional perEMCClusterMS = o2::aod::mintm::minClusterId; + + HistogramRegistry registry{"registry", {}, OutputObjHandlingPolicy::AnalysisObject, false, false}; + + int mRunNumber{-1}; + float dBz{0.f}; + + static constexpr int PhotonPCM = 0; + static constexpr int PhotonEMC = 1; + static constexpr int PairPCMPCM = 0; + static constexpr int PairEMCEMC = 1; + static constexpr int PairPCMEMC = 2; + + // ancestor chains and truth flags, filled once per DF, index = globalIndex of the V0 photon / cluster + std::vector mChainsPCM; + std::vector mChainsEMC; + std::vector mTruthPCM; + std::vector mTruthEMC; + + void defineEMEventCut() + { + fEMEventCut = EMPhotonEventCut("fEMEventCut", "fEMEventCut"); + fEMEventCut.SetRequireSel8(eventcuts.cfgRequireSel8); + fEMEventCut.SetRequireFT0AND(eventcuts.cfgRequireFT0AND); + fEMEventCut.SetZvtxRange(-eventcuts.cfgZvtxMax, +eventcuts.cfgZvtxMax); + fEMEventCut.SetRequireNoTFB(eventcuts.cfgRequireNoTFB); + fEMEventCut.SetRequireNoITSROFB(eventcuts.cfgRequireNoITSROFB); + fEMEventCut.SetRequireNoSameBunchPileup(eventcuts.cfgRequireNoSameBunchPileup); + fEMEventCut.SetRequireVertexITSTPC(eventcuts.cfgRequireVertexITSTPC); + fEMEventCut.SetRequireGoodZvtxFT0vsPV(eventcuts.cfgRequireGoodZvtxFT0vsPV); + fEMEventCut.SetRequireEMCReadoutInMB(eventcuts.cfgRequireEMCReadoutInMB); + fEMEventCut.SetRequireEMCHardwareTriggered(eventcuts.cfgRequireEMCHardwareTriggered); + } + + void defineEMCCut() + { + fEMCCut = EMCPhotonCut("fEMCCut", "fEMCCut"); + + fEMCCut.SetTrackMatchingEtaParams(emccuts.cfgEMCTMEta->at(0), emccuts.cfgEMCTMEta->at(1), emccuts.cfgEMCTMEta->at(2)); + fEMCCut.SetTrackMatchingPhiParams(emccuts.cfgEMCTMPhi->at(0), emccuts.cfgEMCTMPhi->at(1), emccuts.cfgEMCTMPhi->at(2)); + + fEMCCut.SetSecTrackMatchingEtaParams(emccuts.emcSecTMEta->at(0), emccuts.emcSecTMEta->at(1), emccuts.emcSecTMEta->at(2)); + fEMCCut.SetSecTrackMatchingPhiParams(emccuts.emcSecTMPhi->at(0), emccuts.emcSecTMPhi->at(1), emccuts.emcSecTMPhi->at(2)); + fEMCCut.SetMinEoverP(emccuts.cfgEMCEoverp); + + fEMCCut.SetMinE(emccuts.cfgEMCminE); + fEMCCut.SetMinNCell(emccuts.cfgEMCminNCell); + fEMCCut.SetM02Range(emccuts.cfgEMCminM02, emccuts.cfgEMCmaxM02); + fEMCCut.SetTimeRange(emccuts.cfgEMCminTime, emccuts.cfgEMCmaxTime); + fEMCCut.SetUseExoticCut(emccuts.cfgEMCUseExoticCut); + fEMCCut.SetClusterizer(emccuts.clusterDefinition); + fEMCCut.SetUseTM(emccuts.cfgEMCUseTM.value); // disables or enables TM + fEMCCut.SetUseSecondaryTM(emccuts.emcUseSecondaryTM.value); // disables or enables secondary TM + fEMCCut.SetDoQA(emccuts.cfgEnableQA.value); + } + + void definePCMCut() + { + fV0PhotonCut = V0PhotonCut("fV0PhotonCut", "fV0PhotonCut"); + + // for v0 + fV0PhotonCut.SetV0PtRange(pcmcuts.minPtV0, pcmcuts.maxPtV0); + fV0PhotonCut.SetV0EtaRange(pcmcuts.minEtaV0, pcmcuts.maxEtaV0); + fV0PhotonCut.SetMinCosPA(pcmcuts.minCosPA); + fV0PhotonCut.SetMaxPCA(pcmcuts.maxPCA); + fV0PhotonCut.SetMaxChi2KF(pcmcuts.maxChi2KF); + fV0PhotonCut.SetRxyRange(pcmcuts.minRV0, pcmcuts.maxRV0); + fV0PhotonCut.SetAPRange(pcmcuts.maxAlphaAP, pcmcuts.maxQtAP); + fV0PhotonCut.RejectITSib(pcmcuts.rejectV0onITSib); + + // for track + fV0PhotonCut.SetMinNClustersTPC(pcmcuts.minNClusterTPC); + fV0PhotonCut.SetMinNCrossedRowsTPC(pcmcuts.minNCrossedRowsTPC); + fV0PhotonCut.SetMinNCrossedRowsOverFindableClustersTPC(pcmcuts.minNCrossedRowsOverFindableClustersTPC); + fV0PhotonCut.SetMaxFracSharedClustersTPC(pcmcuts.maxFracSharedClustersTPC); + fV0PhotonCut.SetChi2PerClusterTPC(0.f, pcmcuts.macChi2TPC); + fV0PhotonCut.SetTPCNsigmaElRange(pcmcuts.minTPCNSigmaEl, pcmcuts.maxTPCNSigmaEl); + fV0PhotonCut.SetChi2PerClusterITS(0.f, pcmcuts.macChi2ITS); + fV0PhotonCut.SetNClustersITS(pcmcuts.minNClusterITS, 7); + fV0PhotonCut.SetMeanClusterSizeITSob(pcmcuts.minMeanClusterSizeITSob, pcmcuts.maxMeanClusterSizeITSob); + fV0PhotonCut.SetDisableITSonly(pcmcuts.disableITSOnly); + fV0PhotonCut.SetDisableTPConly(pcmcuts.disableTPCOnly); + fV0PhotonCut.SetRequireITSTPC(pcmcuts.requireV0WithITSTPC); + fV0PhotonCut.SetRequireITSonly(pcmcuts.requireV0WithITSOnly); + fV0PhotonCut.SetRequireTPConly(pcmcuts.requireV0WithTPCOnly); + + fV0PhotonCut.setDoQA(pcmcuts.doQA.value); + } + + void init(InitContext&) + { + mRunNumber = 0; + dBz = 0; + + defineEMEventCut(); + defineEMCCut(); + fEMCCut.addQAHistograms(®istry); + definePCMCut(); + fV0PhotonCut.addQAHistograms(®istry); + o2::aod::pwgem::photonmeson::utils::eventhistogram::addEventHistograms(®istry); + + const AxisSpec thnAxisPtGen{thnConfigAxisPt, "#it{p}_{T,Gen} (GeV/#it{c})"}; + const AxisSpec thnAxisPtRec{thnConfigAxisPt, "#it{p}_{T,Rec} (GeV/#it{c})"}; + const AxisSpec thnAxisPGen{thnConfigAxisPt, "#it{p}_{Gen} (GeV/#it{c})"}; + const AxisSpec thnAxisPRec{thnConfigAxisPt, "#it{p}_{Rec} (GeV/#it{c})"}; + const AxisSpec thnAxisPRelative{thnConfigAxisPRelative, "#it{p}_{Rec} - #it{p}_{Gen} / #it{p}_{Gen}"}; + const AxisSpec thnAxisEGen{thnConfigAxisPt, "#it{E}_{Gen} (GeV)"}; + const AxisSpec thnAxisERec{thnConfigAxisPt, "#it{E}_{Rec} (GeV)"}; + const AxisSpec thnAxisERelative{thnConfigAxisERelative, "#it{E}_{Rec} - #it{E}_{Gen} / #it{E}_{Gen}"}; + const AxisSpec thnAxisInvMass{thnConfigAxisInvMass, "#it{M}_{#gamma#gamma} (GeV/#it{c}^{2})"}; + + const AxisSpec thnAxisEtaRelative{thnConfigAxisEtaRelative, "#it{#eta}_{Rec} - #it{#eta}_{Gen} / #it{#eta}_{Gen}"}; + const AxisSpec thnAxisPhiRelative{thnConfigAxisPhiRelative, "#it{#varphi}_{Rec} - #it{#varphi}_{Gen} / #it{#varphi}_{Gen}"}; + + const AxisSpec thnAxisEtaGen{280, -0.7, 0.7, "#it{#eta}_{Gen}"}; + const AxisSpec thnAxisPhiGen{360, 0., o2::constants::math::TwoPI, "#it{#varphi}_{Gen} (rad)"}; + + const AxisSpec thnAxisPairTruth{static_cast(PhotonTruth::kNTruth), -0.5, static_cast(PhotonTruth::kNTruth) - 0.5, "pair truth (worse photon)"}; + const AxisSpec thnAxisPairBucket{static_cast(PairBucket::kNBuckets), -0.5, static_cast(PairBucket::kNBuckets) - 0.5, "closest common ancestor"}; + + AxisSpec thnAxisCentOrMult{1, 0., 1., "Centrality/Multiplicity"}; // placeholder, overwritten in init + if (useCent.value) { + // PbPb: use centrality + thnAxisCentOrMult = {thnConfigAxisCent, "Centrality (%)"}; + } else { + // pp: use multiplicity + thnAxisCentOrMult = {thnConfigAxisMult, "FT0C Multiplicity"}; + } + + auto hPairAncestor = registry.add("PCM/hPairAncestor", "minv vs pT vs cent vs pair truth vs closest common ancestor", HistType::kTHnF, {thnAxisInvMass, thnAxisPtRec, thnAxisCentOrMult, thnAxisPairTruth, thnAxisPairBucket}); + static constexpr std::array(PhotonTruth::kNTruth)> TruthLabels = {"true #gamma", "non #gamma", "fake"}; + for (std::size_t i = 0; i < TruthLabels.size(); ++i) { + hPairAncestor->GetAxis(3)->SetBinLabel(static_cast(i) + 1, TruthLabels[i].data()); + } + static constexpr std::array(PairBucket::kNBuckets)> BucketLabels = {"none", "same particle", "#pi^{0} direct", "#eta direct", "#pi^{0} cross", "#eta cross", "#omega", "K^{0}_{S}", "#eta'", "other hadron", "parton/string"}; + for (std::size_t i = 0; i < BucketLabels.size(); ++i) { + hPairAncestor->GetAxis(4)->SetBinLabel(static_cast(i) + 1, BucketLabels[i].data()); + } + // "other hadron" pairs: which hadron is the common ancestor, and did it split in a decay or a material interaction + const AxisSpec axisOHMass{160, 0., 0.8, "#it{M}_{#gamma#gamma} (GeV/#it{c}^{2})"}; + const AxisSpec axisOHPt{20, 0., 20., "#it{p}_{T} (GeV/#it{c})"}; + const AxisSpec axisOHSpecies{static_cast(OtherHadronSpecies::kNSpecies), -0.5, static_cast(OtherHadronSpecies::kNSpecies) - 0.5, "common ancestor"}; + const AxisSpec axisOHProduction{static_cast(SplitProduction::kNProduction), -0.5, static_cast(SplitProduction::kNProduction) - 0.5, "split of the common ancestor"}; + auto hOtherHadron = registry.add("PCM/OtherHadron/hAncestor", "other hadron pairs: minv vs pT vs ancestor species vs split production vs pair truth", HistType::kTHnF, {axisOHMass, axisOHPt, axisOHSpecies, axisOHProduction, thnAxisPairTruth}); + static constexpr std::array(OtherHadronSpecies::kNSpecies)> SpeciesLabels = {"#pi^{#pm}", "K^{#pm}", "K^{0}_{L}", "#rho", "K^{*}", "#phi", "a_{1}", "K_{1}, K^{*}_{0,2}", "open charm", "c#bar{c}", "open beauty", "b#bar{b}", "other meson", "p", "n", "#Delta", "#Lambda", "#Sigma", "#Xi", "other baryon", "other"}; + for (std::size_t i = 0; i < SpeciesLabels.size(); ++i) { + hOtherHadron->GetAxis(2)->SetBinLabel(static_cast(i) + 1, SpeciesLabels[i].data()); + } + static constexpr std::array(SplitProduction::kNProduction)> ProductionLabels = {"generator decay", "transport decay", "hadronic interaction", "other process"}; + for (std::size_t i = 0; i < ProductionLabels.size(); ++i) { + hOtherHadron->GetAxis(3)->SetBinLabel(static_cast(i) + 1, ProductionLabels[i].data()); + } + for (std::size_t i = 0; i < TruthLabels.size(); ++i) { + hOtherHadron->GetAxis(4)->SetBinLabel(static_cast(i) + 1, TruthLabels[i].data()); + } + registry.add("PCM/OtherHadron/hAncestorPdg", "other hadron pairs: |PDG| of the common ancestor;|PDG| (last bin: #geq 9999);pairs", HistType::kTH1D, {{10000, -0.5, 9999.5}}); + + // mixed event, same axes as hPairAncestor without the bucket (pairs from different events share no ancestor) + auto hMixed = registry.add("PCM/hMixed", "minv vs pT vs cent vs pair truth, mixed event", HistType::kTHnF, {thnAxisInvMass, thnAxisPtRec, thnAxisCentOrMult, thnAxisPairTruth}); + for (std::size_t i = 0; i < TruthLabels.size(); ++i) { + hMixed->GetAxis(3)->SetBinLabel(static_cast(i) + 1, TruthLabels[i].data()); + } + registry.addClone("PCM/", "EMC/"); + registry.addClone("PCM/", "PCMEMC/"); + + // ancestor chain QA + registry.add("QA/PCM/hChainLength", "length of the ancestor chain;chain length;counts", HistType::kTH1D, {{kMaxChainLength, 0.5, kMaxChainLength + 0.5}}); + auto hChainEnd = registry.add("QA/PCM/hChainEnd", "reason the ancestor chain ended;;counts", HistType::kTH1D, {{4, -0.5, 3.5}}); + hChainEnd->GetXaxis()->SetBinLabel(1, "parton-like"); + hChainEnd->GetXaxis()->SetBinLabel(2, "beam nucleus"); + hChainEnd->GetXaxis()->SetBinLabel(3, "no mothers"); + hChainEnd->GetXaxis()->SetBinLabel(4, "truncated"); + auto hChainEndKind = registry.add("QA/PCM/hChainEndKind", "kind of the last chain element;;counts", HistType::kTH1D, {{static_cast(AncestorKind::NClasses), -0.5, static_cast(AncestorKind::NClasses) - 0.5}}); + static constexpr std::array(AncestorKind::NClasses)> KindLabels = {"none", "string/parton", "beam/nucleus", "#pi^{0}", "#eta", "other meson", "baryon", "#gamma", "e^{#pm}", "other"}; + for (std::size_t i = 0; i < KindLabels.size(); ++i) { + hChainEndKind->GetXaxis()->SetBinLabel(static_cast(i) + 1, KindLabels[i].data()); + } + registry.add("QA/PCM/hChainEndPdg", "|PDG| of the last chain element;|PDG| (last bin: #geq 9999, incl. nuclei);counts", HistType::kTH1D, {{10000, -0.5, 9999.5}}); + registry.addClone("QA/PCM/", "QA/EMC/"); + + // origin of PCM-EMC same particle pairs and why track matching did not remove them + if (doprocessPcmEmcal.value) { + const AxisSpec axisOrigin{static_cast(SameParticleOrigin::kNOrigins), -0.5, static_cast(SameParticleOrigin::kNOrigins) - 0.5, "leading contributor of the cluster"}; + const AxisSpec axisStatus{static_cast(LegMatchStatus::kNStatus), -0.5, static_cast(LegMatchStatus::kNStatus) - 0.5, "V0 leg among the matched tracks"}; + const AxisSpec axisPtLeg{40, 0., 20., "#it{p}_{T,leg} (GeV/#it{c})"}; + const AxisSpec axisDeltaEta{100, -0.1, 0.1, "#Delta#eta (cluster - track)"}; + const AxisSpec axisDeltaPhi{100, -0.1, 0.1, "#Delta#varphi (cluster - track)"}; + const AxisSpec axisDp{200, -1., 1., "(#it{p}_{match} - #it{p}_{leg}) / #it{p}_{leg}"}; + static constexpr std::array(SameParticleOrigin::kNOrigins)> OriginLabels = {"leg e^{#pm}", "brems #gamma of leg", "other e^{#pm}", "other"}; + static constexpr std::array(LegMatchStatus::kNStatus)> StatusLabels = {"no match stored", "no match with leg #it{p}", "outside window", "in window, #it{E}/#it{p} #geq cut", "in window"}; + + auto hOrigin = registry.add("PCMEMC/SameParticle/hOrigin", "same particle pairs: leading contributor of the cluster", HistType::kTH3F, {thnAxisInvMass, thnAxisPtRec, axisOrigin}); + for (std::size_t i = 0; i < OriginLabels.size(); ++i) { + hOrigin->GetZaxis()->SetBinLabel(static_cast(i) + 1, OriginLabels[i].data()); + } + registry.add("PCMEMC/SameParticle/hLegEoverP", "cluster from a V0 leg: #it{E}_{cls}/#it{p}_{leg}", HistType::kTH2F, {{200, 0., 4., "#it{E}_{cls}/#it{p}_{leg}"}, axisPtLeg}); + for (const auto& [type, title] : {std::pair{"Prim", "primary"}, {"Sec", "secondary"}}) { + auto hStatus = registry.add(std::format("PCMEMC/SameParticle/hLeg{}Status", type.data()).data(), std::format("cluster from a V0 leg: leg among the matched {} tracks", title.data()).data(), HistType::kTH2F, {axisStatus, axisPtLeg}); + for (std::size_t i = 0; i < StatusLabels.size(); ++i) { + hStatus->GetXaxis()->SetBinLabel(static_cast(i) + 1, StatusLabels[i].data()); + } + registry.add(std::format("PCMEMC/SameParticle/hLeg{}Dp", type.data()).data(), std::format("cluster from a V0 leg: closest {} match in #it{{p}}", title.data()).data(), HistType::kTH1F, {axisDp}); + registry.add(std::format("PCMEMC/SameParticle/hLeg{}Residual", type.data()).data(), std::format("cluster from a V0 leg: residual of the {} match", title.data()).data(), HistType::kTH3F, {axisDeltaEta, axisDeltaPhi, axisPtLeg}); + } + } + + auto hMesonCuts = registry.add("hMesonCuts", "hMesonCuts;;Counts", kTH1D, {{6, 0.5, 6.5}}, false); + hMesonCuts->GetXaxis()->SetBinLabel(1, "in"); + hMesonCuts->GetXaxis()->SetBinLabel(2, "opening angle"); + hMesonCuts->GetXaxis()->SetBinLabel(3, "#it{M}_{#gamma#gamma}"); + hMesonCuts->GetXaxis()->SetBinLabel(4, "#it{p}_{T}"); + hMesonCuts->GetXaxis()->SetBinLabel(5, "conversion cut"); + hMesonCuts->GetXaxis()->SetBinLabel(6, "out"); + if (mesonConfig.cfgEnableQA.value) { + registry.add("mesonQA/hInvMassPt", "Histo for inv pair mass vs pt", HistType::kTH2D, {thnAxisInvMass, thnAxisPtRec}); + } + }; // end init + + template + void initCCDB(TCollision const& collision) + { + if (mRunNumber == collision.runNumber()) { + return; + } + mRunNumber = collision.runNumber(); + + dBz = collision.grpMagField().getNominalL3Field(); + LOG(info) << "Retrieved GRP for timestamp " << collision.timestamp() << " with magnetic field of " << dBz << " kZG"; + fV0PhotonCut.SetD_Bz(dBz); + } + + template + float getCentralityOrMultiplicity(TCollision const& collision) + { + if (useCent.value) { + return getCentrality(collision); + } + // pp: use raw FT0C multiplicity + return collision.multFT0C(); + } + + /// Get the centrality + /// \param collision is the collision with the centrality information + template + float getCentrality(TCollision const& collision) + { + float cent = -999.; + switch (eventcuts.centEstimator) { + case CentralityEstimator::CFT0M: + cent = collision.centFT0M(); + break; + case CentralityEstimator::CFT0A: + cent = collision.centFT0A(); + break; + case CentralityEstimator::CFT0C: + cent = collision.centFT0C(); + break; + default: + LOG(warning) << "Centrality estimator not valid. Possible values are T0M, T0A, T0C. Fallback to T0C"; + cent = collision.centFT0C(); + break; + } + return cent; + } + + /// \brief check if standard event cuts + FT0 occupancy + centrality + QVec good is + /// \param collision collision that will be checked + /// \return true if collision survives all checks, otherwise false + template + bool isFullEventSelected(TCollision const& collision, bool fillHisto = false) + { + if (fillHisto) { + o2::aod::pwgem::photonmeson::utils::eventhistogram::fillEventInfo<0>(®istry, collision); + } + if (!(fEMEventCut.IsSelected(collision))) { + // general event selection + return false; + } + if (!(eventcuts.cfgFT0COccupancyMin <= collision.ft0cOccupancyInTimeRange() && collision.ft0cOccupancyInTimeRange() < eventcuts.cfgFT0COccupancyMax)) { + // occupancy selection + return false; + } + float centOrMult = getCentralityOrMultiplicity(collision); + if (useCent && (centOrMult < eventcuts.cfgMinCent || centOrMult > eventcuts.cfgMaxCent)) { + // event selection + return false; + } + if (fillHisto) { + o2::aod::pwgem::photonmeson::utils::eventhistogram::fillEventInfo<1>(®istry, collision); + registry.fill(HIST("Event/before/hCollisionCounter"), 12.0); // accepted + registry.fill(HIST("Event/after/hCollisionCounter"), 12.0); // accepted + } + return true; + } + + /// \brief Reset the per DF cache of one photon type + /// \param chains ancestor chains, set to -1 (not {}, which gives 0, a valid index) + /// \param truth truth flags, set to kFake + /// \param size size of the full photon table + static void resetCache(std::vector& chains, std::vector& truth, const std::size_t size) + { + Chain empty; + empty.fill(-1); + chains.assign(size, empty); + truth.assign(size, PhotonTruth::kFake); + } + + /// \brief Fill the ancestor chain QA + /// \param mcIter MC particle iterator, has to sit where getAncestorChain left it (last chain element or the nucleus that ended the chain) + /// \param length length of the chain as returned by getAncestorChain + template + void fillChainQA(TMCIter const& mcIter, const int length) + { + static constexpr std::array HistLength = {"QA/PCM/hChainLength", "QA/EMC/hChainLength"}; + static constexpr std::array HistEnd = {"QA/PCM/hChainEnd", "QA/EMC/hChainEnd"}; + static constexpr std::array HistEndKind = {"QA/PCM/hChainEndKind", "QA/EMC/hChainEndKind"}; + static constexpr std::array HistEndPdg = {"QA/PCM/hChainEndPdg", "QA/EMC/hChainEndPdg"}; + if (length <= 0) { + return; + } + static constexpr int MaxPdgBin = 9999; // nuclei and other large codes go to the last bin + const int lastPdg = mcIter.pdgCode(); + int endReason = 3; // chain is full without reaching a parton, a nucleus or a particle without mothers + if (isPartonLike(lastPdg)) { + endReason = 0; + } else if (isNucleus(lastPdg)) { + endReason = 1; + } else if (!mcIter.has_mothers()) { + endReason = 2; + } + registry.fill(HIST(HistLength[photonType]), length); + registry.fill(HIST(HistEnd[photonType]), endReason); + registry.fill(HIST(HistEndKind[photonType]), static_cast(ancestorKindOfPdg(lastPdg))); + registry.fill(HIST(HistEndPdg[photonType]), std::min(std::abs(lastPdg), MaxPdgBin)); + } + + /// \brief Fill ancestor chains and truth flags of all EMC clusters of the DF + /// \param clusters full (unfiltered) cluster table + /// \param mcIter iterator of the full MC particle table, WILL BE CHANGED + /// \param doChainQA fill the ancestor chain QA (off for the mixed event, so chains are not counted twice) + template + void fillEmcCache(TClusters const& clusters, TMCIter& mcIter, const bool doChainQA = true) + { + for (const auto& cluster : clusters) { + const auto mcIds = cluster.emmcparticleIds(); + if (mcIds.empty() || mcIds[0] < 0) { + continue; // noise cluster without MC label: kFake, chain stays -1 + } + auto& chain = mChainsEMC[cluster.globalIndex()]; + const int length = getAncestorChain(mcIter, mcIds[0], chain); // leading contributor + if (doChainQA) { + fillChainQA(mcIter, length); + } + if (length <= 0) { + continue; // leading contributor is a nucleus: kFake, chain stays -1 + } + + // e+- whose mother is a photon (conversion in front of the EMCal): start the chain at the photon, + // otherwise the direct mother check of the pair classification fails for converted photons + mcIter.setCursor(chain[0]); + if (length > 1 && std::abs(mcIter.pdgCode()) == PDG_t::kElectron) { + mcIter.setCursor(chain[1]); + if (mcIter.pdgCode() == PDG_t::kGamma) { + std::shift_left(chain.begin(), chain.end(), 1); + chain.back() = -1; + } + } + mcIter.setCursor(chain[0]); + mTruthEMC[cluster.globalIndex()] = mcIter.pdgCode() == PDG_t::kGamma ? PhotonTruth::kTruePhoton : PhotonTruth::kNonPhoton; + } + } + + /// \brief Fill ancestor chains and truth flags of all V0 photons of the DF + /// \param photons full (unfiltered) V0 photon table + /// \param legs full V0 leg table with MC labels + /// \param mcParticles full MC particle table the leg labels point into + /// \param doChainQA fill the ancestor chain QA (off for the mixed event, so chains are not counted twice) + template + void fillPcmCache(TPhotons const& photons, TLegs const& legs, TMCParticles const& mcParticles, const bool doChainQA = true) + { + auto mcIter = mcParticles.begin(); + auto posLeg = legs.begin(); + auto negLeg = legs.begin(); + auto posMc = mcParticles.begin(); + auto negMc = mcParticles.begin(); + + Chain chainPos{}; + Chain chainNeg{}; + + for (const auto& photon : photons) { + posLeg.setCursor(photon.posTrackId()); + negLeg.setCursor(photon.negTrackId()); + if (posLeg.emmcparticleId() < 0 || negLeg.emmcparticleId() < 0) { + continue; // leg without MC label: kFake, chain stays -1 + } + posMc.setCursor(posLeg.emmcparticleId()); + negMc.setCursor(negLeg.emmcparticleId()); + + auto& chain = mChainsPCM[photon.globalIndex()]; + auto& truth = mTruthPCM[photon.globalIndex()]; + const int photonId = o2::aod::pwgem::dilepton::utils::mcutil::FindCommonMotherFrom2Prongs(posMc, negMc, PDG_t::kPositron, PDG_t::kElectron, PDG_t::kGamma, mcParticles); + if (photonId >= 0) { + truth = PhotonTruth::kTruePhoton; + const int length = getAncestorChain(mcIter, photonId, chain); + if (doChainQA) { + fillChainQA(mcIter, length); + } + continue; + } + + // not a true conversion: K0S/Lambda V0, Dalitz pair or combinatorial fake + const int lengthPos = getAncestorChain(mcIter, posLeg.emmcparticleId(), chainPos); + if (doChainQA) { + fillChainQA(mcIter, lengthPos); // the kept chain is a tail of this one, so it ends the same way + } + getAncestorChain(mcIter, negLeg.emmcparticleId(), chainNeg); + const int ancestor = getFirstCommonAncestor(chainPos, chainNeg); + if (ancestor < 0) { + continue; // no common ancestor: kFake, chain stays -1 + } + const bool isDirectMother = chainPos[1] == ancestor && chainNeg[1] == ancestor; + truth = isDirectMother ? PhotonTruth::kNonPhoton : PhotonTruth::kFake; + const auto start = std::ranges::find(chainPos, ancestor); + std::ranges::fill(std::copy(start, chainPos.end(), chain.begin()), chain.end(), -1); + } + } + + /// \brief Classify a pair by the closest common ancestor of the two candidates + /// \param chain1 ancestor chain of the first candidate + /// \param chain2 ancestor chain of the second candidate + /// \param mcIter iterator of the full MC particle table, WILL BE CHANGED + template + PairBucket classifyPair(Chain const& chain1, Chain const& chain2, TMCIter& mcIter) + { + const int ancestor = getFirstCommonAncestor(chain1, chain2); + if (ancestor < 0) { + return PairBucket::kNone; + } + if (ancestor == chain1[0] || ancestor == chain2[0]) { + return PairBucket::kSameParticle; + } + mcIter.setCursor(ancestor); + const int absPdg = std::abs(mcIter.pdgCode()); + if (absPdg == PDG_t::kGamma || absPdg == PDG_t::kElectron) { + return PairBucket::kSameParticle; // both candidates come from the electromagnetic shower of one photon or electron + } + const bool isDirectMother = chain1[1] == ancestor && chain2[1] == ancestor; + if (isDirectMother && absPdg == PDG_t::kPi0) { + return PairBucket::kSignalPi0; + } + if (isDirectMother && absPdg == o2::constants::physics::Pdg::kEta) { + return PairBucket::kSignalEta; + } + switch (absPdg) { + case PDG_t::kPi0: + return PairBucket::kCrossPi0; + case o2::constants::physics::Pdg::kEta: + return PairBucket::kCrossEta; + case o2::constants::physics::Pdg::kOmega: + return PairBucket::kCrossOmega; + case PDG_t::kK0Short: + return PairBucket::kCrossK0S; + case o2::constants::physics::Pdg::kEtaPrime: + return PairBucket::kCrossEtaPrime; + default: + break; + } + if (isPartonLike(absPdg)) { + return PairBucket::kPartonString; + } + return PairBucket::kCrossOtherHadron; + } + + /// \brief Look for the V0 leg among the tracks matched to an EMC cluster and classify why it did not veto the cluster. + /// The matched track tables store no track index, so the leg is identified by its momentum. + /// \tparam isSecondary matches are secondaries (V0 legs, no E/p condition) instead of primaries + /// \param matches tracks matched to this cluster + /// \param eCluster cluster energy + /// \param pLeg momentum of the V0 leg + /// \param ptLeg transverse momentum of the V0 leg + template + void fillLegMatch(TMatches const& matches, const float eCluster, const float pLeg, const float ptLeg) + { + static constexpr std::array HistStatus = {"PCMEMC/SameParticle/hLegPrimStatus", "PCMEMC/SameParticle/hLegSecStatus"}; + static constexpr std::array HistDp = {"PCMEMC/SameParticle/hLegPrimDp", "PCMEMC/SameParticle/hLegSecDp"}; + static constexpr std::array HistResidual = {"PCMEMC/SameParticle/hLegPrimResidual", "PCMEMC/SameParticle/hLegSecResidual"}; + static constexpr float MaxRelDp = 0.1f; // max |p_match - p_leg| / p_leg to identify the match as the leg + + LegMatchStatus status = LegMatchStatus::kNoMatchStored; + if (matches.size() > 0) { + float bestDp = 1e9f; + float dEta = 0.f, dPhi = 0.f, trackPt = 0.f, trackP = 0.f; + for (const auto& match : matches) { + const float dp = (match.trackP() - pLeg) / pLeg; + if (std::fabs(dp) < std::fabs(bestDp)) { + bestDp = dp; + dEta = match.deltaEta(); + dPhi = match.deltaPhi(); + trackPt = match.trackPt(); + trackP = match.trackP(); + } + } + registry.fill(HIST(HistDp[isSecondary]), bestDp); + if (std::fabs(bestDp) > MaxRelDp) { + status = LegMatchStatus::kNoCompatibleMatch; + } else { + registry.fill(HIST(HistResidual[isSecondary]), dEta, dPhi, ptLeg); + // same windows as EMCPhotonCut::checkTrackMatching + bool isInWindow = false; + if constexpr (isSecondary) { + isInWindow = std::fabs(dEta) <= fEMCCut.GetSecTrackMatchingEta(trackPt) && std::fabs(dPhi) <= fEMCCut.GetSecTrackMatchingPhi(trackPt); + } else { + isInWindow = std::fabs(dEta) <= fEMCCut.GetTrackMatchingEta(trackPt) && std::fabs(dPhi) <= fEMCCut.GetTrackMatchingPhi(trackPt); + } + if (!isInWindow) { + status = LegMatchStatus::kOutsideWindow; + } else if (!isSecondary && eCluster / trackP >= emccuts.cfgEMCEoverp) { + status = LegMatchStatus::kInsideWindowEoverP; + } else { + status = LegMatchStatus::kInsideWindow; + } + } + } + registry.fill(HIST(HistStatus[isSecondary]), static_cast(status), ptLeg); + } + + /// \brief Find out where a PCM-EMC same particle pair comes from: what made the EMC cluster and why track matching did not remove it + /// \param pcm V0 photon + /// \param emc EMC cluster + /// \param legs full V0 leg table with MC labels + /// \param matchedPrims full table of primary tracks matched to clusters + /// \param matchedSeconds full table of secondary tracks (V0 legs) matched to clusters + /// \param mcIter iterator of the full MC particle table, WILL BE CHANGED + template + void diagnoseSameParticle(TPCM const& pcm, TEMC const& emc, TLegs const& legs, TMatchedPrims const& matchedPrims, TMatchedSecs const& matchedSeconds, TMCIter& mcIter) + { + const auto mcIds = emc.emmcparticleIds(); + if (mcIds.empty() || mcIds[0] < 0) { + return; + } + const int lead = mcIds[0]; // leading contributor, before the conversion shift of the chain + const auto posLeg = legs.iteratorAt(pcm.posTrackId()); + const auto negLeg = legs.iteratorAt(pcm.negTrackId()); + const std::array legMcIds = {posLeg.emmcparticleId(), negLeg.emmcparticleId()}; + + SameParticleOrigin origin = SameParticleOrigin::kOther; + int legIndex = -1; + if (lead == legMcIds[0] || lead == legMcIds[1]) { + origin = SameParticleOrigin::kLegElectron; + legIndex = lead == legMcIds[0] ? 0 : 1; + } else { + mcIter.setCursor(lead); + const int absPdg = std::abs(mcIter.pdgCode()); + const int motherId = mcIter.has_mothers() ? mcIter.mothersIds()[0] : -1; + if (absPdg == PDG_t::kGamma && motherId >= 0 && (motherId == legMcIds[0] || motherId == legMcIds[1])) { + origin = SameParticleOrigin::kBremsFromLeg; + } else if (absPdg == PDG_t::kElectron) { + origin = SameParticleOrigin::kOtherElectron; + } + } + + ROOT::Math::PtEtaPhiMVector v1(pcm.pt(), pcm.eta(), pcm.phi(), 0.); + ROOT::Math::PtEtaPhiMVector v2(emc.pt(), emc.eta(), emc.phi(), 0.); + ROOT::Math::PtEtaPhiMVector vMeson = v1 + v2; + registry.fill(HIST("PCMEMC/SameParticle/hOrigin"), vMeson.M(), vMeson.Pt(), static_cast(origin)); + if (origin != SameParticleOrigin::kLegElectron) { + return; + } + + const auto& leg = legIndex == 0 ? posLeg : negLeg; + const float pLeg = std::hypot(leg.px(), leg.py(), leg.pz()); + const float ptLeg = std::hypot(leg.px(), leg.py()); + if (pLeg <= 0.f) { + return; + } + registry.fill(HIST("PCMEMC/SameParticle/hLegEoverP"), emc.e() / pLeg, ptLeg); + fillLegMatch(matchedPrims.sliceBy(perEMCClusterMT, emc.globalIndex()), emc.e(), pLeg, ptLeg); + fillLegMatch(matchedSeconds.sliceBy(perEMCClusterMS, emc.globalIndex()), emc.e(), pLeg, ptLeg); + } + + /// \brief Species class of the common ancestor of an "other hadron" pair + /// \param pdg PDG code of the ancestor + static OtherHadronSpecies otherHadronSpecies(const int pdg) + { + const int absPdg = std::abs(pdg); + switch (absPdg) { + case PDG_t::kPiPlus: + return OtherHadronSpecies::kChargedPion; + case PDG_t::kKPlus: + return OtherHadronSpecies::kChargedKaon; + case PDG_t::kK0Long: + return OtherHadronSpecies::kK0L; + case PDG_t::kRho770_0: + case PDG_t::kRho770Plus: + return OtherHadronSpecies::kRho; + case o2::constants::physics::Pdg::kK0Star892: + case o2::constants::physics::Pdg::kKPlusStar892: + return OtherHadronSpecies::kKstar; + case o2::constants::physics::Pdg::kPhi: + return OtherHadronSpecies::kPhi; + case PDG_t::kProton: + return OtherHadronSpecies::kProton; + case PDG_t::kNeutron: + return OtherHadronSpecies::kNeutron; + case 1114: // o2-linter: disable=pdg/explicit-code (Delta-, no named constant) + case 2114: // o2-linter: disable=pdg/explicit-code (Delta0, no named constant) + case 2214: // o2-linter: disable=pdg/explicit-code (Delta+, no named constant) + case 2224: // o2-linter: disable=pdg/explicit-code (Delta++, no named constant) + return OtherHadronSpecies::kDelta; + case PDG_t::kLambda0: + return OtherHadronSpecies::kLambda; + case PDG_t::kSigmaMinus: + case PDG_t::kSigma0: + case PDG_t::kSigmaPlus: + return OtherHadronSpecies::kSigma; + case PDG_t::kXiMinus: + case 3322: // o2-linter: disable=pdg/explicit-code (Xi0, no named constant) + return OtherHadronSpecies::kXi; + default: + break; + } + constexpr int MinHadronPdg = 100; + constexpr int MaxHadronPdg = 10000000; + if (absPdg < MinHadronPdg || absPdg >= MaxHadronPdg) { + return OtherHadronSpecies::kOther; + } + // PDG numbering: mesons have no third quark, i.e. the thousands digit is 0 (also for 9000000 type resonances) + if ((absPdg / 1000) % 10 != 0) { // o2-linter: disable=magic-number (PDG digit) + return OtherHadronSpecies::kOtherBaryon; + } + // mesons: hundreds digit = heavier quark, tens digit = lighter quark, last digit = 2J+1 + // (also covers excited states with a radial/orbital prefix, e.g. 10441, 20443, 100553) + const int heavyQuark = (absPdg / 100) % 10; // o2-linter: disable=magic-number (PDG digit) + const int lightQuark = (absPdg / 10) % 10; // o2-linter: disable=magic-number (PDG digit) + constexpr int Strange = 3, Charm = 4, Beauty = 5; + constexpr int ExcitedPrefix = 10000; + constexpr int SpinTwo = 5; // 2J+1 for J = 2 + if (heavyQuark == Charm) { + return lightQuark == Charm ? OtherHadronSpecies::kCharmonium : OtherHadronSpecies::kOpenCharm; + } + if (heavyQuark == Beauty) { + return lightQuark == Beauty ? OtherHadronSpecies::kBottomonium : OtherHadronSpecies::kOpenBeauty; + } + if (absPdg == 20213 || absPdg == 20113) { // o2-linter: disable=pdg/explicit-code (a1(1260), no named constant) + return OtherHadronSpecies::kA1; + } + if (heavyQuark == Strange && lightQuark < Strange && (absPdg >= ExcitedPrefix || absPdg % 10 == SpinTwo)) { // o2-linter: disable=magic-number (PDG digit) + return OtherHadronSpecies::kExcitedK; + } + return OtherHadronSpecies::kOtherMeson; + } + + /// \brief How an MC particle was produced + /// \param mcIter iterator of the full MC particle table, WILL BE CHANGED + /// \param id index of the particle + template + static SplitProduction productionOf(TMCIter& mcIter, const int id) + { + mcIter.setCursor(id); + if (mcIter.producedByGenerator()) { + return SplitProduction::kGeneratorDecay; + } + switch (mcIter.getProcess()) { + case TMCProcess::kPDecay: + return SplitProduction::kTransportDecay; + case TMCProcess::kPHadronic: + case TMCProcess::kPHInhelastic: + case TMCProcess::kPNuclearAbsorption: + case TMCProcess::kPNCapture: + case TMCProcess::kPPbarAnnihilation: + case TMCProcess::kPNbarAnnihilation: + case TMCProcess::kPEvaporation: + case TMCProcess::kPNuclearFission: + case TMCProcess::kPPhotoNuclear: + case TMCProcess::kPPhotonInhelastic: + return SplitProduction::kHadronicInteraction; + default: + return SplitProduction::kOtherProcess; + } + } + + /// \brief Fill which hadron is the common ancestor of an "other hadron" pair and how it split into the two branches + /// \param chain1 ancestor chain of the first candidate + /// \param chain2 ancestor chain of the second candidate + /// \param pairTruth truth flag of the pair + /// \param vMeson pair four vector + /// \param mcIter iterator of the full MC particle table, WILL BE CHANGED + template + void diagnoseOtherHadron(Chain const& chain1, Chain const& chain2, const PhotonTruth pairTruth, ROOT::Math::PtEtaPhiMVector const& vMeson, TMCIter& mcIter) + { + static constexpr std::array HistAncestor = {"PCM/OtherHadron/hAncestor", "EMC/OtherHadron/hAncestor", "PCMEMC/OtherHadron/hAncestor"}; + static constexpr std::array HistAncestorPdg = {"PCM/OtherHadron/hAncestorPdg", "EMC/OtherHadron/hAncestorPdg", "PCMEMC/OtherHadron/hAncestorPdg"}; + static constexpr int MaxPdgBin = 9999; + + const int ancestor = getFirstCommonAncestor(chain1, chain2); + const auto pos1 = std::distance(chain1.begin(), std::ranges::find(chain1, ancestor)); + const auto pos2 = std::distance(chain2.begin(), std::ranges::find(chain2, ancestor)); + if (ancestor < 0 || pos1 < 1 || pos2 < 1 || pos1 >= std::ssize(chain1) || pos2 >= std::ssize(chain2)) { + return; // not possible for "other hadron" (the ancestor is neither of the candidates), protection only + } + // the particle directly below the ancestor in each branch tells how the ancestor split + const SplitProduction production = std::max(productionOf(mcIter, chain1[pos1 - 1]), productionOf(mcIter, chain2[pos2 - 1])); + mcIter.setCursor(ancestor); + const int pdg = mcIter.pdgCode(); + registry.fill(HIST(HistAncestor[pairType]), vMeson.M(), vMeson.Pt(), static_cast(otherHadronSpecies(pdg)), static_cast(production), static_cast(pairTruth)); + registry.fill(HIST(HistAncestorPdg[pairType]), std::min(std::abs(pdg), MaxPdgBin)); + } + + /// \brief Apply the pair cuts and fill the pair ancestor histogram + /// \return bucket of the pair, std::nullopt if the pair was rejected + template + std::optional fillPair(TPhoton1 const& g1, TPhoton2 const& g2, Chain const& chain1, Chain const& chain2, const PhotonTruth truth1, const PhotonTruth truth2, const float centOrMult, TMCIter& mcIter) + { + static constexpr std::array HistPairAncestor = {"PCM/hPairAncestor", "EMC/hPairAncestor", "PCMEMC/hPairAncestor"}; + + ROOT::Math::PtEtaPhiMVector v1(g1.pt(), g1.eta(), g1.phi(), 0.); + ROOT::Math::PtEtaPhiMVector v2(g2.pt(), g2.eta(), g2.phi(), 0.); + ROOT::Math::PtEtaPhiMVector vMeson = v1 + v2; + float openingAngle = std::acos(v1.Vect().Dot(v2.Vect()) / (v1.P() * v2.P())); + + registry.fill(HIST("hMesonCuts"), 1); + if (openingAngle <= mesonConfig.minOpenAngle) { + registry.fill(HIST("hMesonCuts"), 2); + return std::nullopt; + } + registry.fill(HIST("hMesonCuts"), 6); + + const PhotonTruth pairTruth = std::max(truth1, truth2); + const PairBucket bucket = classifyPair(chain1, chain2, mcIter); + registry.fill(HIST(HistPairAncestor[pairType]), vMeson.M(), vMeson.Pt(), centOrMult, static_cast(pairTruth), static_cast(bucket)); + if (bucket == PairBucket::kCrossOtherHadron) { + diagnoseOtherHadron(chain1, chain2, pairTruth, vMeson, mcIter); + } + if (mesonConfig.cfgEnableQA.value) { + registry.fill(HIST("mesonQA/hInvMassPt"), vMeson.M(), vMeson.Pt()); + } + return bucket; + } + + /// \brief eta a cluster would have if its own vertex vzOld is swapped for vzNew, same as in taskPi0FlowEMC + /// \param etaOld old eta value + /// \param vzOld old primary vertex z position + /// \param vzNew new primary vertex z position + /// \param emcalR radius of the EMCal + static float correctEtaForVertexShift(float etaOld, float vzOld, float vzNew, float emcalR) + { + return std::asinh(std::sinh(etaOld) + (vzOld - vzNew) / emcalR); + } + + /// \brief Apply the pair cuts to a mixed event pair and fill the mixed event histogram + /// \param eta2 eta of the second photon, possibly re-projected onto the vertex of the first collision + template + void fillMixedPair(TPhoton1 const& g1, TPhoton2 const& g2, const float eta2, const PhotonTruth truth1, const PhotonTruth truth2, const float centOrMult) + { + static constexpr std::array HistMixed = {"PCM/hMixed", "EMC/hMixed", "PCMEMC/hMixed"}; + + ROOT::Math::PtEtaPhiMVector v1(g1.pt(), g1.eta(), g1.phi(), 0.); + ROOT::Math::PtEtaPhiMVector v2(g2.pt(), eta2, g2.phi(), 0.); + ROOT::Math::PtEtaPhiMVector vMeson = v1 + v2; + float openingAngle = std::acos(v1.Vect().Dot(v2.Vect()) / (v1.P() * v2.P())); + if (openingAngle <= mesonConfig.minOpenAngle) { + return; + } + registry.fill(HIST(HistMixed[pairType]), vMeson.M(), vMeson.Pt(), centOrMult, static_cast(std::max(truth1, truth2))); + } + + // PCM-PCM same event + void processPcm(Colls const& collisions, PcmPhotons const& photons, PcmMcLegs const& legs, McParticles const& mcParticles) + { + if (photons.size() <= 0) { + LOG(info) << "Skipping DF because there are not photons!"; + return; + } + // none of the photon tables is filtered, so its size is the size of the full table the globalIndex points into + resetCache(mChainsPCM, mTruthPCM, photons.size()); + fillPcmCache(photons, legs, mcParticles); + + EMBitFlags v0flags(photons.size()); + fV0PhotonCut.AreSelectedRunning(v0flags, photons, ®istry); + + auto mcIter = mcParticles.begin(); + for (const auto& collision : collisions) { + initCCDB(collision); + if (!isFullEventSelected(collision, true)) { + continue; + } + const float centOrMult = getCentralityOrMultiplicity(collision); + auto photonsPerCollision = photons.sliceBy(perCollisionPCM, collision.globalIndex()); + + for (const auto& [g1, g2] : combinations(CombinationsStrictlyUpperIndexPolicy(photonsPerCollision, photonsPerCollision))) { + if (!(v0flags.test(g1.globalIndex())) || !(v0flags.test(g2.globalIndex()))) { + continue; + } + fillPair(g1, g2, mChainsPCM[g1.globalIndex()], mChainsPCM[g2.globalIndex()], mTruthPCM[g1.globalIndex()], mTruthPCM[g2.globalIndex()], centOrMult, mcIter); + } + } + } + PROCESS_SWITCH(TwoPhotonMcTask, processPcm, "Process for pcm photons", false); + + // EMCal-EMCal same event + void processEmcal(Colls const& collisions, EMCalPhotons const& clusters, MinMTracks const& matchedPrims, MinMSTracks const& matchedSeconds, McParticles const& mcParticles) + { + if (clusters.size() <= 0) { + LOG(info) << "Skipping DF because there are not photons!"; + return; + } + auto mcIter = mcParticles.begin(); + // none of the photon tables is filtered, so its size is the size of the full table the globalIndex points into + resetCache(mChainsEMC, mTruthEMC, clusters.size()); + fillEmcCache(clusters, mcIter); + + EMBitFlags emcFlags(clusters.size()); + fEMCCut.AreSelectedRunning(emcFlags, clusters, matchedPrims, matchedSeconds, ®istry); + + for (const auto& collision : collisions) { + initCCDB(collision); + if (!isFullEventSelected(collision, true)) { + continue; + } + const float centOrMult = getCentralityOrMultiplicity(collision); + auto photonsPerCollision = clusters.sliceBy(perCollisionEMC, collision.globalIndex()); + + for (const auto& [g1, g2] : combinations(CombinationsStrictlyUpperIndexPolicy(photonsPerCollision, photonsPerCollision))) { + if (!(emcFlags.test(g1.globalIndex())) || !(emcFlags.test(g2.globalIndex()))) { + continue; + } + fillPair(g1, g2, mChainsEMC[g1.globalIndex()], mChainsEMC[g2.globalIndex()], mTruthEMC[g1.globalIndex()], mTruthEMC[g2.globalIndex()], centOrMult, mcIter); + } + } + } + PROCESS_SWITCH(TwoPhotonMcTask, processEmcal, "Process for emcal photons", false); + + // PCM-EMCal same event + void processPcmEmcal(Colls const& collisions, EMCalPhotons const& clusters, PcmPhotons const& photons, PcmMcLegs const& legs, MinMTracks const& matchedPrims, MinMSTracks const& matchedSeconds, McParticles const& mcParticles) + { + if (clusters.size() <= 0 || photons.size() <= 0) { + LOG(info) << "Skipping DF because there are not photons!"; + return; + } + auto mcIter = mcParticles.begin(); + // none of the photon tables is filtered, so its size is the size of the full table the globalIndex points into + resetCache(mChainsPCM, mTruthPCM, photons.size()); + resetCache(mChainsEMC, mTruthEMC, clusters.size()); + fillPcmCache(photons, legs, mcParticles); + fillEmcCache(clusters, mcIter); + + EMBitFlags emcFlags(clusters.size()); + fEMCCut.AreSelectedRunning(emcFlags, clusters, matchedPrims, matchedSeconds, ®istry); + EMBitFlags v0flags(photons.size()); + fV0PhotonCut.AreSelectedRunning(v0flags, photons, ®istry); + + for (const auto& collision : collisions) { + initCCDB(collision); + if (!isFullEventSelected(collision, true)) { + continue; + } + const float centOrMult = getCentralityOrMultiplicity(collision); + auto photonsEMCPerCollision = clusters.sliceBy(perCollisionEMC, collision.globalIndex()); + auto photonsPCMPerCollision = photons.sliceBy(perCollisionPCM, collision.globalIndex()); + + // PCM - EMC loop, two different tables: all combinations + for (const auto& [g1, g2] : combinations(CombinationsFullIndexPolicy(photonsPCMPerCollision, photonsEMCPerCollision))) { + if (!(v0flags.test(g1.globalIndex())) || !(emcFlags.test(g2.globalIndex()))) { + continue; + } + const auto bucket = fillPair(g1, g2, mChainsPCM[g1.globalIndex()], mChainsEMC[g2.globalIndex()], mTruthPCM[g1.globalIndex()], mTruthEMC[g2.globalIndex()], centOrMult, mcIter); + if (bucket == PairBucket::kSameParticle) { + diagnoseSameParticle(g1, g2, legs, matchedPrims, matchedSeconds, mcIter); + } + } + } + } + PROCESS_SWITCH(TwoPhotonMcTask, processPcmEmcal, "Process for pcm and emcal photons", true); + + // PCM-PCM mixed event + void processPcmMixed(Colls const& collisions, PcmPhotons const& photons, PcmMcLegs const& legs, McParticles const& mcParticles) + { + if (photons.size() <= 0) { + LOG(info) << "Skipping DF because there are not photons!"; + return; + } + // the mixed pairs are built inside the DF, so the per DF cache is valid for them + resetCache(mChainsPCM, mTruthPCM, photons.size()); + fillPcmCache(photons, legs, mcParticles, false); + + EMBitFlags v0flags(photons.size()); + fV0PhotonCut.AreSelectedRunning(v0flags, photons); + + BinningType binningMixedEvent{{mixingConfig.cfgVtxBins, mixingConfig.cfgCentBins}}; + auto photonTuple = std::make_tuple(photons); + SameKindPair pair{binningMixedEvent, mixingConfig.cfgMixingDepth, -1, collisions, photonTuple, &cache}; + + for (const auto& [c1, photons1, c2, photons2] : pair) { + if (!isFullEventSelected(c1) || !isFullEventSelected(c2)) { + continue; + } + const float centOrMult = getCentralityOrMultiplicity(c1); + for (const auto& [g1, g2] : combinations(CombinationsFullIndexPolicy(photons1, photons2))) { + if (!(v0flags.test(g1.globalIndex())) || !(v0flags.test(g2.globalIndex()))) { + continue; + } + fillMixedPair(g1, g2, g2.eta(), mTruthPCM[g1.globalIndex()], mTruthPCM[g2.globalIndex()], centOrMult); + } + } + } + PROCESS_SWITCH(TwoPhotonMcTask, processPcmMixed, "Process for pcm photons, mixed event", false); + + // EMCal-EMCal mixed event + void processEmcalMixed(Colls const& collisions, EMCalPhotons const& clusters, MinMTracks const& matchedPrims, MinMSTracks const& matchedSeconds, McParticles const& mcParticles) + { + if (clusters.size() <= 0) { + LOG(info) << "Skipping DF because there are not photons!"; + return; + } + auto mcIter = mcParticles.begin(); + // the mixed pairs are built inside the DF, so the per DF cache is valid for them + resetCache(mChainsEMC, mTruthEMC, clusters.size()); + fillEmcCache(clusters, mcIter, false); + + EMBitFlags emcFlags(clusters.size()); + fEMCCut.AreSelectedRunning(emcFlags, clusters, matchedPrims, matchedSeconds); + + BinningType binningMixedEvent{{mixingConfig.cfgVtxBins, mixingConfig.cfgCentBins}}; + auto clusterTuple = std::make_tuple(clusters); + SameKindPair pair{binningMixedEvent, mixingConfig.cfgMixingDepth, -1, collisions, clusterTuple, &cache}; + + for (const auto& [c1, clusters1, c2, clusters2] : pair) { + if (!isFullEventSelected(c1) || !isFullEventSelected(c2)) { + continue; + } + const float centOrMult = getCentralityOrMultiplicity(c1); + for (const auto& [g1, g2] : combinations(CombinationsFullIndexPolicy(clusters1, clusters2))) { + if (!(emcFlags.test(g1.globalIndex())) || !(emcFlags.test(g2.globalIndex()))) { + continue; + } + // cluster 2 seen from the vertex of collision 1, as in taskPi0FlowEMC + const float eta2 = mixingConfig.cfgCorrectMixedVtxEta.value ? correctEtaForVertexShift(g2.eta(), c2.posZ(), c1.posZ(), mixingConfig.cfgEmcalEffRadius.value) : g2.eta(); + fillMixedPair(g1, g2, eta2, mTruthEMC[g1.globalIndex()], mTruthEMC[g2.globalIndex()], centOrMult); + } + } + } + PROCESS_SWITCH(TwoPhotonMcTask, processEmcalMixed, "Process for emcal photons, mixed event", false); + + // PCM-EMCal mixed event + void processPcmEmcalMixed(Colls const& collisions, EMCalPhotons const& clusters, PcmPhotons const& photons, PcmMcLegs const& legs, MinMTracks const& matchedPrims, MinMSTracks const& matchedSeconds, McParticles const& mcParticles) + { + if (clusters.size() <= 0 || photons.size() <= 0) { + LOG(info) << "Skipping DF because there are not photons!"; + return; + } + auto mcIter = mcParticles.begin(); + // the mixed pairs are built inside the DF, so the per DF cache is valid for them + resetCache(mChainsPCM, mTruthPCM, photons.size()); + resetCache(mChainsEMC, mTruthEMC, clusters.size()); + fillPcmCache(photons, legs, mcParticles, false); + fillEmcCache(clusters, mcIter, false); + + EMBitFlags emcFlags(clusters.size()); + fEMCCut.AreSelectedRunning(emcFlags, clusters, matchedPrims, matchedSeconds); + EMBitFlags v0flags(photons.size()); + fV0PhotonCut.AreSelectedRunning(v0flags, photons); + + BinningType binningMixedEvent{{mixingConfig.cfgVtxBins, mixingConfig.cfgCentBins}}; + auto associatedTables = std::make_tuple(photons, clusters); + Pair pair{binningMixedEvent, mixingConfig.cfgMixingDepth, -1, collisions, associatedTables, &cache}; + + for (const auto& [c1, photonsPCM, c2, photonsEMC] : pair) { + if (!isFullEventSelected(c1) || !isFullEventSelected(c2)) { + continue; + } + const float centOrMult = getCentralityOrMultiplicity(c1); + for (const auto& [g1, g2] : combinations(CombinationsFullIndexPolicy(photonsPCM, photonsEMC))) { + if (!(v0flags.test(g1.globalIndex())) || !(emcFlags.test(g2.globalIndex()))) { + continue; + } + fillMixedPair(g1, g2, g2.eta(), mTruthPCM[g1.globalIndex()], mTruthEMC[g2.globalIndex()], centOrMult); + } + } + } + PROCESS_SWITCH(TwoPhotonMcTask, processPcmEmcalMixed, "Process for pcm and emcal photons, mixed event", false); + +}; // End struct TwoPhotonMcTask + +WorkflowSpec defineDataProcessing(ConfigContext const& context) +{ + return WorkflowSpec{adaptAnalysisTask(context)}; +} diff --git a/PWGEM/PhotonMeson/Utils/MCUtilities.h b/PWGEM/PhotonMeson/Utils/MCUtilities.h index 270d6d2ead9..1326f513aef 100644 --- a/PWGEM/PhotonMeson/Utils/MCUtilities.h +++ b/PWGEM/PhotonMeson/Utils/MCUtilities.h @@ -32,12 +32,16 @@ #include #include #include +#include #include //_______________________________________________________________________ namespace o2::aod::pwgem::photonmeson::utils::mcutil { +template +concept is_unfiltered_iterator = o2::soa::is_iterator && !o2::soa::has_filtered_policy; + constexpr float kVertexEps = 1e-4f; // cm template @@ -1012,6 +1016,66 @@ float trueOpeningAngle(TMCParticles const& mcParticles, const int id1, const int return std::acos(std::clamp(c, -1.f, 1.f)); } +//_______________________________________________________________________ +/// \brief Whether the PDG code is a quark, gluon, diquark or string/cluster, i.e. where the decay chain ends +[[nodiscard]] inline bool isPartonLike(int pdg) +{ + const int absPdg = std::abs(pdg); + const bool isDiquark = absPdg > 1000 && absPdg < 6000 && (absPdg / 10) % 10 == 0; // o2-linter: disable=magic-number + return (absPdg >= PDG_t::kDown && absPdg <= PDG_t::kTop) || absPdg == PDG_t::kGluon || absPdg == 91 || absPdg == 92 || isDiquark; // o2-linter: disable=magic-number +} + +//_______________________________________________________________________ +/// \brief Whether the PDG code is a nucleus (10LZZZAAAI), e.g. one of the colliding beam nuclei +[[nodiscard]] inline bool isNucleus(int pdg) +{ + constexpr int NucleusOffset = 1000000000; + return std::abs(pdg) >= NucleusOffset; +} + +//_______________________________________________________________________ +/// \brief Fill the ancestor chain of a particle, starting with the particle itself. +/// Ends with the first parton-like ancestor (included), a particle without mothers, or when chain is full. +/// A nucleus (the colliding beam nuclei) is never part of the chain, otherwise all particles of an event would share it. +/// Unused entries of chain are set to -1. +/// \param mcParticle iterator of the full (unsliced) MC particle table. !WILL BE CHANGED! +/// Afterwards it points to the last chain element, or to the nucleus that ended the chain. +/// \param startId global index of the particle of interest +/// \param chain output, global indices, closest first +/// \return number of filled entries +template +int getAncestorChain(TMCParticle& mcParticle, int startId, std::span chain) +{ + int length = 0; + int id = startId; + while (id >= 0 && length < std::ssize(chain)) { + mcParticle.setCursor(id); + if (isNucleus(mcParticle.pdgCode())) { + break; // colliding nuclei are not an ancestor + } + chain[length++] = id; + if (isPartonLike(mcParticle.pdgCode()) || !mcParticle.has_mothers()) { + break; + } + id = mcParticle.mothersIds()[0]; + } + std::ranges::fill(chain.subspan(length), -1); + return length; +} + +//_______________________________________________________________________ +/// \brief Find the closest common ancestor of two ancestor chains +/// \param chain1 ancestor chain of first particle (particle itself first, may be padded with -1) +/// \param chain2 ancestor chain of second particle (particle itself first, may be padded with -1) +/// \return global index of the closest common ancestor, -1 if none found +[[nodiscard]] inline int getFirstCommonAncestor(std::span chain1, std::span chain2) +{ + chain1 = chain1.first(std::ranges::distance(chain1.begin(), std::ranges::find(chain1, -1))); + chain2 = chain2.first(std::ranges::distance(chain2.begin(), std::ranges::find(chain2, -1))); + const auto it = std::ranges::find_first_of(chain1, chain2); + return it != chain1.end() ? *it : -1; +} + } // namespace o2::aod::pwgem::photonmeson::utils::mcutil #endif // PWGEM_PHOTONMESON_UTILS_MCUTILITIES_H_ diff --git a/PWGEM/PhotonMeson/Utils/NMHistograms.h b/PWGEM/PhotonMeson/Utils/NMHistograms.h index ab69f458636..448e94c3007 100644 --- a/PWGEM/PhotonMeson/Utils/NMHistograms.h +++ b/PWGEM/PhotonMeson/Utils/NMHistograms.h @@ -33,12 +33,19 @@ namespace o2::aod::pwgem::photonmeson::utils::nmhistogram { +/// \brief Add the pair histograms +/// \param isMC book the MC truth histograms (only mass and pT axes) instead of the same/mixed event histograms +/// \param pairname name of the pair used in the axis titles +/// \param sparseFullAxes data only: add EP, occupancy and z-vertex axes to the same/mixed event histograms +/// \param occBins data only: occupancy bin edges +/// \param epBins data only: event plane bin edges +/// \param vtxBins data only: z-vertex bin edges inline void addNMHistograms(o2::framework::HistogramRegistry* fRegistry, bool isMC, + const char* pairname = "#gamma#gamma", bool sparseFullAxes = false, std::vector const& occBins = {0, 500, 1000, 2000, 3000, 5000, 10000}, std::vector const& epBins = {-o2::constants::math::PIHalf, -o2::constants::math::PIQuarter, 0, o2::constants::math::PIQuarter, o2::constants::math::PIHalf}, - std::vector const& vtxBins = {-10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10}, - const char* pairname = "#gamma#gamma") + std::vector const& vtxBins = {-10, -8, -6, -4, -2, 0, 2, 4, 6, 8, 10}) { // !!Don't change pt,eta,y binning. These binnings have to be consistent with binned data at skimming.!! std::vector ptbins; @@ -57,9 +64,6 @@ inline void addNMHistograms(o2::framework::HistogramRegistry* fRegistry, bool is } const o2::framework::AxisSpec axis_pt{ptbins, Form("p_{T,%s} (GeV/c)", pairname)}; const o2::framework::AxisSpec axis_mass{400, 0, 0.8, Form("m_{%s} (GeV/c^{2})", pairname)}; - const o2::framework::AxisSpec axis_occ{occBins, "occupancy"}; - const o2::framework::AxisSpec axis_ep{epBins, "#Psi_{2} (rad)"}; - const o2::framework::AxisSpec axis_vtx{vtxBins, "Z_{vtx} (cm)"}; if (isMC) { fRegistry->add("Pair/Pi0/hs_Primary", "rec. true pi0", o2::framework::HistType::kTHnSparseD, {axis_mass, axis_pt}, true); @@ -84,6 +88,9 @@ inline void addNMHistograms(o2::framework::HistogramRegistry* fRegistry, bool is fRegistry->get(HIST("Generated/Eta/hPtY"))->SetYTitle("rapidity |y|"); } else { if (sparseFullAxes) { + const o2::framework::AxisSpec axis_occ{occBins, "occupancy"}; + const o2::framework::AxisSpec axis_ep{epBins, "#Psi_{2} (rad)"}; + const o2::framework::AxisSpec axis_vtx{vtxBins, "Z_{vtx} (cm)"}; fRegistry->add("Pair/same/hs", "diphoton", o2::framework::HistType::kTHnSparseD, {axis_mass, axis_pt, axis_ep, axis_occ, axis_vtx}, true);