diff --git a/PWGEM/PhotonMeson/Tasks/Pi0EtaToGammaGammaEMCEMC.cxx b/PWGEM/PhotonMeson/Tasks/Pi0EtaToGammaGammaEMCEMC.cxx index 1e0567d6ee6..a57661b5805 100644 --- a/PWGEM/PhotonMeson/Tasks/Pi0EtaToGammaGammaEMCEMC.cxx +++ b/PWGEM/PhotonMeson/Tasks/Pi0EtaToGammaGammaEMCEMC.cxx @@ -14,6 +14,7 @@ /// \author D. Sekihata, daiki.sekihata@cern.ch #include "PWGEM/PhotonMeson/Core/Pi0EtaToGammaGamma.h" +#include "PWGEM/PhotonMeson/DataModel/GammaTablesRedux.h" #include "PWGEM/PhotonMeson/DataModel/gammaTables.h" #include "PWGEM/PhotonMeson/Utils/PairUtilities.h" @@ -26,11 +27,11 @@ using namespace o2::aod; using namespace o2::framework; using namespace o2::aod::pwgem::photonmeson::photonpair; -using MyEMCClusters = soa::Join; +using MyEMCClusters = soa::Join; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { return WorkflowSpec{ - adaptAnalysisTask>(cfgc, TaskName{"pi0eta-to-gammagamma-emcemc"}), + adaptAnalysisTask>(cfgc, TaskName{"pi0eta-to-gammagamma-emcemc"}), }; } diff --git a/PWGEM/PhotonMeson/Tasks/Pi0EtaToGammaGammaMCEMCEMC.cxx b/PWGEM/PhotonMeson/Tasks/Pi0EtaToGammaGammaMCEMCEMC.cxx index 36fa6c724b9..da93889b9db 100644 --- a/PWGEM/PhotonMeson/Tasks/Pi0EtaToGammaGammaMCEMCEMC.cxx +++ b/PWGEM/PhotonMeson/Tasks/Pi0EtaToGammaGammaMCEMCEMC.cxx @@ -14,6 +14,7 @@ /// \author D. Sekihata, daiki.sekihata@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" @@ -26,11 +27,11 @@ using namespace o2::aod; using namespace o2::framework; using namespace o2::aod::pwgem::photonmeson::photonpair; -using MyEMCClusters = soa::Join; +using MyEMCClusters = soa::Join; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) { return WorkflowSpec{ - adaptAnalysisTask>(cfgc, TaskName{"pi0eta-to-gammagamma-mc-emcemc"}), + adaptAnalysisTask>(cfgc, TaskName{"pi0eta-to-gammagamma-mc-emcemc"}), }; } diff --git a/PWGEM/PhotonMeson/Tasks/photonResoTask.cxx b/PWGEM/PhotonMeson/Tasks/photonResoTask.cxx index 575952f9abc..a9b4d5ed81d 100644 --- a/PWGEM/PhotonMeson/Tasks/photonResoTask.cxx +++ b/PWGEM/PhotonMeson/Tasks/photonResoTask.cxx @@ -86,8 +86,10 @@ struct PhotonResoTask { // configurable axis ConfigurableAxis thnConfigAxisInvMass{"thnConfigAxisInvMass", {400, 0.0, 0.8}, "invariant mass axis for the neutral meson"}; ConfigurableAxis thnConfigAxisPt{"thnConfigAxisPt", {400, 0., 20.}, "pT axis for the neutral meson"}; - ConfigurableAxis thnConfigAxisERelative{"thnConfigAxisERelative", {400, -1., 19.}, "(E rec - E true) / E true axis"}; - ConfigurableAxis thnConfigAxisPRelative{"thnConfigAxisPRelative", {60, -1., 2.}, "(P rec - P true) / P true axis"}; + 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", {20, 0., 100.}, "centrality axis for the current event"}; ConfigurableAxis thnConfigAxisMult{"thnConfigAxisMult", {60, 0., 60000.}, "multiplicity axis for the current event"}; Configurable useCent{"useCent", 0, "flag to enable usage of centrality instead of multiplicity as axis."}; @@ -296,11 +298,11 @@ struct PhotonResoTask { 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 thnAxisEtaGen{280, -0.7, 0.7, "#it{#eta}_{Gen}"}; - const AxisSpec thnAxisEtaRec{280, -0.7, 0.7, "#it{#eta}_{Rec}"}; + 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 thnAxisPhiRec{360, 0., o2::constants::math::TwoPI, "#it{#varphi}_{Rec} (rad)"}; AxisSpec thnAxisCentOrMult{1, 0., 1., "Centrality/Multiplicity"}; // placeholder, overwritten in init if (useCent.value) { @@ -313,24 +315,28 @@ struct PhotonResoTask { registry.add("EMCal/hPhotonReso", "EMCal photon rec pT vs true pT vs cent", HistType::kTH3D, {thnAxisPtRec, thnAxisPtGen, thnAxisCentOrMult}); registry.add("EMCal/hConvPhotonReso", "EMCal conversion photon rec pT vs true pT vs cent ", HistType::kTH3D, {thnAxisPtRec, thnAxisPtGen, thnAxisCentOrMult}); + registry.add("EMCal/hFullConvPhotonReso", "full EMCal conversion photon rec pT vs true pT vs cent ", HistType::kTH3D, {thnAxisPtRec, thnAxisPtGen, thnAxisCentOrMult}); registry.add("EMCal/hErecEmcPhotons", "EMCal photon rec E - true E vs true E vs cent", HistType::kTH3D, {thnAxisERelative, thnAxisEGen, thnAxisCentOrMult}); registry.add("EMCal/hErecEmcConvPhotons", "EMCal conversion photon rec E - true E vs true E vs cent ", HistType::kTH3D, {thnAxisERelative, thnAxisEGen, thnAxisCentOrMult}); + registry.add("EMCal/hErecEmcFullConvPhotons", "full EMCal conversion photon rec E - true E vs true E vs cent ", HistType::kTH3D, {thnAxisERelative, thnAxisEGen, thnAxisCentOrMult}); registry.add("EMCal/hPi0Reso", "EMCal pi0 rec pT vs true pT vs min vs cent ", HistType::kTHnSparseF, {thnAxisPtRec, thnAxisPtGen, thnConfigAxisInvMass, thnAxisCentOrMult}); registry.add("EMCal/hEtaReso", "EMCal eta rec pT vs true pT vs min vs cent ", HistType::kTHnSparseF, {thnAxisPtRec, thnAxisPtGen, thnConfigAxisInvMass, thnAxisCentOrMult}); - registry.add("EMCal/hPhotonResoEta", "EMCal photon rec eta vs true eta vs cent", HistType::kTH3D, {thnAxisEtaRec, thnAxisEtaGen, thnAxisCentOrMult}); - registry.add("EMCal/hConvPhotonResoEta", "EMCal conversion photon rec eta vs true eta vs cent ", HistType::kTH3D, {thnAxisEtaRec, thnAxisEtaGen, thnAxisCentOrMult}); + registry.add("EMCal/hPhotonResoEta", "EMCal photon rec eta vs true eta vs cent", HistType::kTH3D, {thnAxisEtaGen, thnAxisEtaRelative, thnAxisCentOrMult}); + registry.add("EMCal/hConvPhotonResoEta", "EMCal conversion photon rec eta vs true eta vs cent ", HistType::kTH3D, {thnAxisEtaGen, thnAxisEtaRelative, thnAxisCentOrMult}); + registry.add("EMCal/hFullConvPhotonResoEta", "full EMCal conversion photon rec eta vs true eta vs cent ", HistType::kTH3D, {thnAxisEtaGen, thnAxisEtaRelative, thnAxisCentOrMult}); - registry.add("EMCal/hPi0ResoEta", "EMCal pi0 rec eta vs true eta vs min vs cent ", HistType::kTHnSparseF, {thnAxisEtaRec, thnAxisEtaGen, thnConfigAxisInvMass, thnAxisCentOrMult}); - registry.add("EMCal/hEtaResoEta", "EMCal eta rec eta vs true eta vs min vs cent ", HistType::kTHnSparseF, {thnAxisEtaRec, thnAxisEtaGen, thnConfigAxisInvMass, thnAxisCentOrMult}); + registry.add("EMCal/hPi0ResoEta", "EMCal pi0 rec eta vs true eta vs min vs cent ", HistType::kTHnSparseF, {thnAxisEtaGen, thnAxisEtaRelative, thnConfigAxisInvMass, thnAxisCentOrMult}); + registry.add("EMCal/hEtaResoEta", "EMCal eta rec eta vs true eta vs min vs cent ", HistType::kTHnSparseF, {thnAxisEtaGen, thnAxisEtaRelative, thnConfigAxisInvMass, thnAxisCentOrMult}); - registry.add("EMCal/hPhotonResoPhi", "EMCal photon rec phi vs true phi vs cent", HistType::kTH3D, {thnAxisPhiRec, thnAxisPhiGen, thnAxisCentOrMult}); - registry.add("EMCal/hConvPhotonResoPhi", "EMCal conversion photon rec phi vs true phi vs cent ", HistType::kTH3D, {thnAxisPhiRec, thnAxisPhiGen, thnAxisCentOrMult}); + registry.add("EMCal/hPhotonResoPhi", "EMCal photon rec phi vs true phi vs cent", HistType::kTH3D, {thnAxisPhiGen, thnAxisPhiRelative, thnAxisCentOrMult}); + registry.add("EMCal/hConvPhotonResoPhi", "EMCal conversion photon rec phi vs true phi vs cent ", HistType::kTH3D, {thnAxisPhiGen, thnAxisPhiRelative, thnAxisCentOrMult}); + registry.add("EMCal/hFullConvPhotonResoPhi", "full EMCal conversion photon rec phi vs true phi vs cent ", HistType::kTH3D, {thnAxisPhiGen, thnAxisPhiRelative, thnAxisCentOrMult}); - registry.add("EMCal/hPi0ResoPhi", "EMCal pi0 rec phi vs true phi vs min vs cent ", HistType::kTHnSparseF, {thnAxisPhiRec, thnAxisPhiGen, thnConfigAxisInvMass, thnAxisCentOrMult}); - registry.add("EMCal/hEtaResoPhi", "EMCal eta rec phi vs true phi vs min vs cent ", HistType::kTHnSparseF, {thnAxisPhiRec, thnAxisPhiGen, thnConfigAxisInvMass, thnAxisCentOrMult}); + registry.add("EMCal/hPi0ResoPhi", "EMCal pi0 rec phi vs true phi vs min vs cent ", HistType::kTHnSparseF, {thnAxisPhiGen, thnAxisPhiRelative, thnConfigAxisInvMass, thnAxisCentOrMult}); + registry.add("EMCal/hEtaResoPhi", "EMCal eta rec phi vs true phi vs min vs cent ", HistType::kTHnSparseF, {thnAxisPhiGen, thnAxisPhiRelative, thnConfigAxisInvMass, thnAxisCentOrMult}); registry.add("PCM/hPhotonReso", "PCM photon rec pT vs true pT vs cent", HistType::kTH3D, {thnAxisPtRec, thnAxisPtGen, thnAxisCentOrMult}); @@ -460,6 +466,7 @@ struct PhotonResoTask { // create iterators for photon mc particles auto mcPhoton1 = mcParticles.begin(); auto mcPhoton2 = mcParticles.begin(); + auto mcConvLeg = mcParticles.begin(); // leg iterators for PCM auto pos1 = legs.begin(); @@ -489,21 +496,49 @@ struct PhotonResoTask { // we only want to look at the largest contribution mcPhoton1.setCursor(photonEMC.emmcparticleIds()[0]); + // if the largest contribution is photon, just fill the photon histograms if (std::abs(mcPhoton1.pdgCode()) == PDG_t::kGamma) { registry.fill(HIST("EMCal/hPhotonReso"), photonEMC.pt(), mcPhoton1.pt(), centOrMult); - registry.fill(HIST("EMCal/hPhotonResoEta"), photonEMC.eta(), mcPhoton1.eta(), centOrMult); - registry.fill(HIST("EMCal/hPhotonResoPhi"), photonEMC.phi(), mcPhoton1.phi(), centOrMult); + registry.fill(HIST("EMCal/hPhotonResoEta"), mcPhoton1.eta(), (photonEMC.eta() - mcPhoton1.eta()) / mcPhoton1.eta(), centOrMult); + registry.fill(HIST("EMCal/hPhotonResoPhi"), mcPhoton1.phi(), (photonEMC.phi() - mcPhoton1.phi()) / mcPhoton1.phi(), centOrMult); registry.fill(HIST("EMCal/hErecEmcPhotons"), (photonEMC.e() - mcPhoton1.e()) / mcPhoton1.e(), mcPhoton1.e(), centOrMult); } else if (std::abs(mcPhoton1.pdgCode()) == PDG_t::kElectron) { - if (!o2::aod::pwgem::photonmeson::utils::mcutil::isMotherPDG(mcPhoton1, PDG_t::kGamma)) { + // if largest contribution is e+ or e-, check if its from a photon conversion + int32_t mcPhoton1MotherId = o2::aod::pwgem::photonmeson::utils::mcutil::getMotherIndexFromChain(mcPhoton1, PDG_t::kGamma); + int32_t mcConvLegMotherId = -1; + bool hasBothLegs = false; + if (mcPhoton1MotherId == -1) { continue; } - registry.fill(HIST("EMCal/hConvPhotonReso"), photonEMC.pt(), mcPhoton1.pt(), centOrMult); - registry.fill(HIST("EMCal/hConvPhotonResoEta"), photonEMC.eta(), mcPhoton1.eta(), centOrMult); - registry.fill(HIST("EMCal/hConvPhotonResoPhi"), photonEMC.phi(), mcPhoton1.phi(), centOrMult); - registry.fill(HIST("EMCal/hErecEmcConvPhotons"), (photonEMC.e() - mcPhoton1.e()) / mcPhoton1.e(), mcPhoton1.e(), centOrMult); - } - } + // mcPhoton1 now points to a photon that produced the e+/e- that was the largest contributor in the cluster + // check if both legs of the photon hit the cluster to decide which conversion histogram to fill + for (int32_t emmcparticleId = 1; emmcparticleId < static_cast(photonEMC.emmcparticleIds().size()); ++emmcparticleId) { + mcConvLeg.setCursor(photonEMC.emmcparticleIds()[emmcparticleId]); + const int convLegPDG = mcConvLeg.pdgCode(); // store the pdg value, because getMotherIndexFromChain will alter the cursor! + mcConvLegMotherId = o2::aod::pwgem::photonmeson::utils::mcutil::getMotherIndexFromChain(mcConvLeg, PDG_t::kGamma); + if (!(std::abs(convLegPDG) == PDG_t::kElectron) || mcConvLegMotherId == -1) { + continue; + } + if (mcPhoton1MotherId == mcConvLegMotherId) { + hasBothLegs = true; + break; + } + } + if (hasBothLegs) { + // if its a conversion with both legs, fill the full conversion histograms + registry.fill(HIST("EMCal/hFullConvPhotonReso"), photonEMC.pt(), mcPhoton1.pt(), centOrMult); + registry.fill(HIST("EMCal/hFullConvPhotonResoEta"), mcPhoton1.eta(), (photonEMC.eta() - mcPhoton1.eta()) / mcPhoton1.eta(), centOrMult); + registry.fill(HIST("EMCal/hFullConvPhotonResoPhi"), mcPhoton1.phi(), (photonEMC.phi() - mcPhoton1.phi()) / mcPhoton1.phi(), centOrMult); + registry.fill(HIST("EMCal/hErecEmcFullConvPhotons"), (photonEMC.e() - mcPhoton1.e()) / mcPhoton1.e(), mcPhoton1.e(), centOrMult); + } else { + // if its a conversion with only one leg, fill the conversion histograms + registry.fill(HIST("EMCal/hConvPhotonReso"), photonEMC.pt(), mcPhoton1.pt(), centOrMult); + registry.fill(HIST("EMCal/hConvPhotonResoEta"), mcPhoton1.eta(), (photonEMC.eta() - mcPhoton1.eta()) / mcPhoton1.eta(), centOrMult); + registry.fill(HIST("EMCal/hConvPhotonResoPhi"), mcPhoton1.phi(), (photonEMC.phi() - mcPhoton1.phi()) / mcPhoton1.phi(), centOrMult); + registry.fill(HIST("EMCal/hErecEmcConvPhotons"), (photonEMC.e() - mcPhoton1.e()) / mcPhoton1.e(), mcPhoton1.e(), centOrMult); + } + } // else if (std::abs(mcPhoton1.pdgCode()) == PDG_t::kElectron) + } // for (const auto& photonEMC : photonsEMCPerCollision) { for (const auto& photonPCM : photonsPCMPerCollision) { if (!(v0flags.test(photonPCM.globalIndex()))) { @@ -543,6 +578,8 @@ struct PhotonResoTask { ROOT::Math::PtEtaPhiMVector v2(g2.pt(), g2.eta(), g2.phi(), 0.); ROOT::Math::PtEtaPhiMVector vMeson = v1 + v2; + const float mesonPhi = RecoDecay::constrainAngle(vMeson.Phi()); + float openingAngle = std::acos(v1.Vect().Dot(v2.Vect()) / (v1.P() * v2.P())); registry.fill(HIST("hMesonCuts"), 1); @@ -586,14 +623,14 @@ struct PhotonResoTask { if (pi0id >= 0) { const auto pi0mc = mcParticles.iteratorAt(pi0id); registry.fill(HIST("EMCal/hPi0Reso"), vMeson.Pt(), pi0mc.pt(), vMeson.M(), centOrMult); - registry.fill(HIST("EMCal/hPi0ResoEta"), vMeson.Eta(), pi0mc.eta(), vMeson.M(), centOrMult); - registry.fill(HIST("EMCal/hPi0ResoPhi"), RecoDecay::constrainAngle(vMeson.Phi()), pi0mc.phi(), vMeson.M(), centOrMult); + registry.fill(HIST("EMCal/hPi0ResoEta"), pi0mc.eta(), (vMeson.Eta() - pi0mc.eta()) / pi0mc.eta(), vMeson.M(), centOrMult); + registry.fill(HIST("EMCal/hPi0ResoPhi"), pi0mc.phi(), (mesonPhi - pi0mc.phi()) / pi0mc.phi(), vMeson.M(), centOrMult); } if (etaid >= 0) { const auto etamc = mcParticles.iteratorAt(etaid); registry.fill(HIST("EMCal/hEtaReso"), vMeson.Pt(), etamc.pt(), vMeson.M(), centOrMult); - registry.fill(HIST("EMCal/hEtaResoEta"), vMeson.Eta(), etamc.eta(), vMeson.M(), centOrMult); - registry.fill(HIST("EMCal/hEtaResoPhi"), RecoDecay::constrainAngle(vMeson.Phi()), etamc.phi(), vMeson.M(), centOrMult); + registry.fill(HIST("EMCal/hEtaResoEta"), etamc.eta(), (vMeson.Eta() - etamc.eta()) / etamc.eta(), vMeson.M(), centOrMult); + registry.fill(HIST("EMCal/hEtaResoPhi"), etamc.phi(), (mesonPhi - etamc.phi()) / etamc.phi(), vMeson.M(), centOrMult); } } } diff --git a/PWGEM/PhotonMeson/Utils/MCUtilities.h b/PWGEM/PhotonMeson/Utils/MCUtilities.h index 9d766466a00..90568398c63 100644 --- a/PWGEM/PhotonMeson/Utils/MCUtilities.h +++ b/PWGEM/PhotonMeson/Utils/MCUtilities.h @@ -282,21 +282,46 @@ bool isGammaGammaDecay(TMCParticle const& mcParticle, TMCParticles const& mcPart } //_______________________________________________________________________ -// Go up the decay chain of a mcparticle looking for a mother with the given pdg codes, if found return true else false -// E.g. if electron cluster is coming from a photon return true, if primary electron return false +/// \brief Go up the decay chain of a mcparticle looking for a mother with the given pdg codes, if found return true else false +/// E.g. if electron cluster is coming from a photon return true, if primary electron return false +/// \param mcparticle iterator of mxparticle, WILL BE CHANGED by this function! +/// \param motherPDG target mother PDG value +/// \param depth how many steps in the chain this check should go maximum before failing template -bool isMotherPDG(T& mcparticle, const int motherPDG, const int Depth = 10) // o2-linter: disable=pdg/explicit-code (false positive) +bool isMotherPDG(T& mcparticle, const int motherPDG, const int depth = 10) // o2-linter: disable=pdg/explicit-code (false positive) { - if (!mcparticle.has_mothers() || Depth < 1) { + if (!mcparticle.has_mothers() || depth < 1) { return false; } int motherid = mcparticle.mothersIds()[0]; mcparticle.setCursor(motherid); - if (mcparticle.pdgCode() != motherPDG) { + if (mcparticle.pdgCode() == motherPDG) { return true; // The mother has the required pdg code, so return its daughters global mc particle code. } else { - return isMotherPDG(mcparticle, motherPDG, Depth - 1); + return isMotherPDG(mcparticle, motherPDG, depth - 1); + } +} + +//_______________________________________________________________________ +/// \brief Go up the decay chain of a mcparticle looking for a mother with the given pdg codes, if found return id else -1 +/// E.g. if electron cluster is coming from a photon return true, if primary electron return false +/// \param mcparticle iterator of mxparticle, WILL BE CHANGED by this function! +/// \param motherPDG target mother PDG value +/// \param depth how many steps in the chain this check should go maximum before failing +template +int32_t getMotherIndexFromChain(T& mcparticle, const int motherPDG, const int depth = 10) // o2-linter: disable=pdg/explicit-code (false positive) +{ + if (!mcparticle.has_mothers() || depth < 1) { + return -1; + } + + int32_t motherid = mcparticle.mothersIds()[0]; + mcparticle.setCursor(motherid); + if (mcparticle.pdgCode() == motherPDG) { + return motherid; // The mother has the required pdg code, so return its daughters global mc particle code. + } else { + return getMotherIndexFromChain(mcparticle, motherPDG, depth - 1); } }