diff --git a/PWGEM/PhotonMeson/Core/Pi0EtaToGammaGammaMC.h b/PWGEM/PhotonMeson/Core/Pi0EtaToGammaGammaMC.h index 8e05897fd6d..48dacaddee8 100644 --- a/PWGEM/PhotonMeson/Core/Pi0EtaToGammaGammaMC.h +++ b/PWGEM/PhotonMeson/Core/Pi0EtaToGammaGammaMC.h @@ -634,8 +634,8 @@ struct Pi0EtaToGammaGammaMC { photonid1 = o2::aod::pwgem::dilepton::utils::mcutil::FindCommonMotherFrom2Prongs(pos1mc, ele1mc, -11, 11, 22, mcparticles); photonid2 = o2::aod::pwgem::dilepton::utils::mcutil::FindCommonMotherFrom2Prongs(pos2mc, ele2mc, -11, 11, 22, mcparticles); } else if constexpr (pairtype == o2::aod::pwgem::photonmeson::photonpair::PairType::kEMCEMC) { - auto cluster1mcparticle = mcparticles.iteratorAt(g1.emmcparticleId()); - auto cluster2mcparticle = mcparticles.iteratorAt(g2.emmcparticleId()); + auto cluster1mcparticle = mcparticles.iteratorAt(g1.emmcparticleIds()[0]); + auto cluster2mcparticle = mcparticles.iteratorAt(g2.emmcparticleIds()[0]); 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}); diff --git a/PWGEM/PhotonMeson/Core/V0PhotonCut.h b/PWGEM/PhotonMeson/Core/V0PhotonCut.h index 98db109534b..d1b2edc3949 100644 --- a/PWGEM/PhotonMeson/Core/V0PhotonCut.h +++ b/PWGEM/PhotonMeson/Core/V0PhotonCut.h @@ -897,23 +897,20 @@ class V0PhotonCut : public TNamed } template - bool IsConversionPointInAcceptance(TMCPhoton const& mcphoton) const + bool IsConversionPointInAcceptance(TMCPhoton const& mcphoton, float convRadius) const { - - float rGenXY = std::sqrt(std::pow(mcphoton.vx(), 2) + std::pow(mcphoton.vy(), 2)); - // eta cut - if (mcphoton.eta() >= mMinV0Eta && mcphoton.eta() <= mMaxV0Eta) { + if (mcphoton.eta() < mMinV0Eta || mcphoton.eta() > mMaxV0Eta) { return false; } // radius cut - if (rGenXY < mMinRxy || mMaxRxy < rGenXY) { + if (convRadius < mMinRxy || mMaxRxy < convRadius) { return false; } // line cut - if (rGenXY < std::abs(mcphoton.vz()) * std::tan(2 * std::atan(std::exp(-mMaxV0Eta))) - mMaxMarginZ) { + if (convRadius <= std::abs(mcphoton.vz()) * std::tan(2 * std::atan(std::exp(-mMaxV0Eta))) - mMaxMarginZ) { return false; } diff --git a/PWGEM/PhotonMeson/DataModel/gammaTables.h b/PWGEM/PhotonMeson/DataModel/gammaTables.h index 5fee7952cfd..0851203654e 100644 --- a/PWGEM/PhotonMeson/DataModel/gammaTables.h +++ b/PWGEM/PhotonMeson/DataModel/gammaTables.h @@ -98,25 +98,25 @@ using V0LegMCLabel = V0LegMCLabels::iterator; // * EMC cluster mc label tables: // 1. EMCALMCClusters in EMCalClusters.h: Vectors of global mc particle ids and energy fractions of the cluster // 2. EMCClusterMCLabels: Vector of global mc particle ids -// 3. EMEMCClusterMCLabels: EM MC particle ID of largest contributor to cluster +// 3. EMEMCClusterMCLabels: Vector of EM MC particle ID of largest contributor to cluster namespace emcclustermclabel { -DECLARE_SOA_ARRAY_INDEX_COLUMN(EMMCParticle, emmcparticle); //! +DECLARE_SOA_ARRAY_INDEX_COLUMN(McParticle, emmcparticle) //! } // namespace emcclustermclabel // NOTE: MC labels. This table has one vector of global mc particle ids for each reconstructed emc cluster (joinable with emccluster table) DECLARE_SOA_TABLE(EMCClusterMCLabels, "AOD", "EMCClsMCLABEL", //! - emcclustermclabel::EMMCParticleIds); + emcclustermclabel::McParticleIds); using EMCClusterMCLabel = EMCClusterMCLabels::iterator; namespace ememcclustermclabel { -DECLARE_SOA_INDEX_COLUMN(EMMCParticle, emmcparticle); //! +DECLARE_SOA_ARRAY_INDEX_COLUMN(EMMCParticle, emmcparticle); //! } // namespace ememcclustermclabel -// NOTE: MC labels. This table has one entry for each reconstructed emc cluster (joinable with emccluster table) +// NOTE: MC labels. This table has a vector of entries for each reconstructed emc cluster (joinable with emccluster table) DECLARE_SOA_TABLE(EMEMCClusterMCLabels, "AOD", "EMEMCClsMCLABEL", //! - ememcclustermclabel::EMMCParticleId); + ememcclustermclabel::EMMCParticleIds); using EMEMCClusterMCLabel = EMEMCClusterMCLabels::iterator; namespace v0leg diff --git a/PWGEM/PhotonMeson/TableProducer/associateMCinfoPhoton.cxx b/PWGEM/PhotonMeson/TableProducer/associateMCinfoPhoton.cxx index da560a82784..26f83737285 100644 --- a/PWGEM/PhotonMeson/TableProducer/associateMCinfoPhoton.cxx +++ b/PWGEM/PhotonMeson/TableProducer/associateMCinfoPhoton.cxx @@ -16,14 +16,18 @@ /// \author Daiki Sekihata (daiki.sekihata@cern.ch) /// +#include "PWGEM/PhotonMeson/DataModel/GammaTablesRedux.h" #include "PWGEM/PhotonMeson/DataModel/gammaTables.h" #include "PWGEM/PhotonMeson/Utils/MCUtilities.h" -#include "Framework/ASoAHelpers.h" -#include "Framework/AnalysisDataModel.h" -#include "Framework/AnalysisTask.h" -#include "Framework/runDataProcessing.h" -#include "ReconstructionDataFormats/Track.h" +#include +#include +#include +#include +#include +#include + +#include #include #include @@ -33,11 +37,12 @@ using namespace o2::framework; using namespace o2::framework::expressions; using namespace o2::soa; using namespace o2::aod::pwgem::photonmeson::utils::mcutil; +using namespace o2::constants::physics; using MyCollisionsMC = soa::Join; using TracksMC = soa::Join; using FwdTracksMC = soa::Join; -using MyEMCClusters = soa::Join; +using MyEMCClusters = soa::Join; struct AssociateMCInfoPhoton { enum SubSystem { @@ -71,26 +76,27 @@ struct AssociateMCInfoPhoton { // !!Don't change pt,eta,y binning. These binnings have to be consistent with binned data at analysis.!! std::vector ptbins; - for (int i = 0; i < 2; i++) { + for (int i = 0; i < 2; i++) { // o2-linter: disable=magic-number (just numbers for binning) ptbins.emplace_back(0.05 * (i - 0) + 0.0); // from 0 to 0.05 GeV/c, every 0.05 GeV/c } - for (int i = 2; i < 51; i++) { + for (int i = 2; i < 51; i++) { // o2-linter: disable=magic-number (just numbers for binning) ptbins.emplace_back(0.1 * (i - 2) + 0.1); // from 0.1 to 4.9 GeV/c, every 0.1 GeV/c } - for (int i = 51; i < 61; i++) { + for (int i = 51; i < 61; i++) { // o2-linter: disable=magic-number (just numbers for binning) ptbins.emplace_back(0.5 * (i - 51) + 5.0); // from 5 to 9.5 GeV/c, every 0.5 GeV/c } - for (int i = 61; i < 72; i++) { + for (int i = 61; i < 72; i++) { // o2-linter: disable=magic-number (just numbers for binning) ptbins.emplace_back(1.0 * (i - 61) + 10.0); // from 10 to 20 GeV/c, every 1 GeV/c } const AxisSpec axisPt{ptbins, "p_{T} (GeV/c)"}; const AxisSpec axisRapidity{{0.0, +0.8, +0.9}, "rapidity |y|"}; - static constexpr std::string_view ParticleNames[9] = { + static constexpr uint NParticleNames = 9; + static constexpr std::string_view ParticleNames[NParticleNames] = { "Gamma", "Pi0", "Eta", "Omega", "Phi", "ChargedPion", "ChargedKaon", "K0S", "Lambda"}; - for (int i = 0; i < 9; i++) { + for (uint i = 0; i < NParticleNames; i++) { registry.add(Form("Generated/h2PtY_%s", ParticleNames[i].data()), Form("Generated %s", ParticleNames[i].data()), kTH2F, {axisPt, axisRapidity}, true); } @@ -113,7 +119,7 @@ struct AssociateMCInfoPhoton { std::vector genEta; // primary, pt, y template - void skimmingMC(MyCollisionsMC const& collisions, aod::BCs const&, aod::McCollisions const&, aod::McParticles const& mcParticles, TTracks const& o2tracks, TFwdTracks const&, TPCMs const& v0photons, TPCMLegs const&, TPHOSs const&, TEMCs const& emcclusters, TEMPrimaryElectrons const& emprimaryelectrons) + void skimmingMC(MyCollisionsMC const& collisions, aod::BCs const&, aod::McCollisions const& mcCollisions, aod::McParticles const& mcParticles, TTracks const& o2tracks, TFwdTracks const&, TPCMs const& v0photons, TPCMLegs const& legs, TPHOSs const&, TEMCs const& emcclusters, TEMPrimaryElectrons const& emprimaryelectrons) { // temporary variables used for the indexing of the skimmed MC stack std::map fNewLabels; @@ -124,6 +130,21 @@ struct AssociateMCInfoPhoton { int fCounters[2] = {0, 0}; //! [0] - particle counter, [1] - event counter auto hBinFinder = registry.get(HIST("Generated/h2PtY_Gamma")); + // collision iterator from EMCal cluster + auto collisionIter = collisions.begin(); + + // mc particle iterator for photons + auto mcPhoton = mcParticles.begin(); + + // mc particles iterator for mother + auto motherParticle = mcParticles.begin(); + + // mc particles iterator for mother and other mc particles + auto mcParticleIter = mcParticles.begin(); + + // mc collision iter + auto mcCollisionIter = mcCollisions.begin(); + for (const auto& collision : collisions) { registry.fill(HIST("hEventCounter"), 1); @@ -142,25 +163,25 @@ struct AssociateMCInfoPhoton { std::fill(genPi0.begin(), genPi0.end(), 0); std::fill(genEta.begin(), genEta.end(), 0); - auto mcCollision = collision.mcCollision(); + mcCollisionIter.setCursor(collision.mcCollisionId()); // store mc particles - auto groupedMcParticles = mcParticles.sliceBy(perMcCollision, mcCollision.globalIndex()); + auto groupedMcParticles = mcParticles.sliceBy(perMcCollision, mcCollisionIter.globalIndex()); for (const auto& mcParticle : groupedMcParticles) { // store necessary information for denominator of efficiency if ((mcParticle.isPhysicalPrimary() || mcParticle.producedByGenerator()) && std::fabs(mcParticle.y()) < 0.9f && mcParticle.pt() < 20.f) { auto binNumber = hBinFinder->FindBin(mcParticle.pt(), std::fabs(mcParticle.y())); // caution: pack switch (std::abs(mcParticle.pdgCode())) { - case 22: + case PDG_t::kGamma: registry.fill(HIST("Generated/h2PtY_Gamma"), mcParticle.pt(), std::fabs(mcParticle.y())); genGamma[binNumber]++; break; - case 111: + case PDG_t::kPi0: if (requireGammaGammaDecay && !isGammaGammaDecay(mcParticle, mcParticles)) continue; registry.fill(HIST("Generated/h2PtY_Pi0"), mcParticle.pt(), std::fabs(mcParticle.y())); genPi0[binNumber]++; break; - case 221: + case Pdg::kEta: if (requireGammaGammaDecay && !isGammaGammaDecay(mcParticle, mcParticles)) continue; registry.fill(HIST("Generated/h2PtY_Eta"), mcParticle.pt(), std::fabs(mcParticle.y())); @@ -173,15 +194,15 @@ struct AssociateMCInfoPhoton { } // end of mc track loop // make an entry for this MC event only if it was not already added to the table - if (!(fEventLabels.find(mcCollision.globalIndex()) != fEventLabels.end())) { - mcevents(mcCollision.globalIndex(), mcCollision.generatorsID(), mcCollision.posX(), mcCollision.posY(), mcCollision.posZ(), mcCollision.impactParameter(), mcCollision.eventPlaneAngle()); - fEventLabels[mcCollision.globalIndex()] = fCounters[1]; + if (!(fEventLabels.find(mcCollisionIter.globalIndex()) != fEventLabels.end())) { + mcevents(mcCollisionIter.globalIndex(), mcCollisionIter.generatorsID(), mcCollisionIter.posX(), mcCollisionIter.posY(), mcCollisionIter.posZ(), mcCollisionIter.impactParameter(), mcCollisionIter.eventPlaneAngle()); + fEventLabels[mcCollisionIter.globalIndex()] = fCounters[1]; fCounters[1]++; binnedGenPt(genGamma, genPi0, genEta); } // LOGF(info, "collision.globalIndex() = %d , mceventlabels.lastIndex() = %d", collision.globalIndex(), mceventlabels.lastIndex()); - mceventlabels(fEventLabels.find(mcCollision.globalIndex())->second, collision.mcMask()); + mceventlabels(fEventLabels.find(mcCollisionIter.globalIndex())->second, collision.mcMask()); for (const auto& mcParticle : groupedMcParticles) { // store necessary information for denominator of efficiency if (mcParticle.pt() < 1e-3 || std::fabs(mcParticle.vz()) > 250 || std::sqrt(std::pow(mcParticle.vx(), 2) + std::pow(mcParticle.vy(), 2)) > max_rxy_gen) { @@ -195,30 +216,30 @@ struct AssociateMCInfoPhoton { // Note that pi0 from weak decay gives producedByGenerator() = false // LOGF(info,"index = %d , mc track pdg = %d , producedByGenerator = %d , isPhysicalPrimary = %d", mcParticle.index(), mcParticle.pdgCode(), mcParticle.producedByGenerator(), mcParticle.isPhysicalPrimary()); - if (std::abs(pdg) == 11 && mcParticle.has_mothers() && !(mcParticle.isPhysicalPrimary() || mcParticle.producedByGenerator())) { // secondary electrons. i.e. ele/pos from photon conversions. - int motherid = mcParticle.mothersIds()[0]; // first mother index - auto mp = mcParticles.iteratorAt(motherid); + if (std::abs(pdg) == PDG_t::kElectron && mcParticle.has_mothers() && !(mcParticle.isPhysicalPrimary() || mcParticle.producedByGenerator())) { // secondary electrons. i.e. ele/pos from photon conversions. + int motherid = mcParticle.mothersIds()[0]; // first mother index + motherParticle.setCursor(motherid); if (std::sqrt(std::pow(mcParticle.vx(), 2) + std::pow(mcParticle.vy(), 2)) < std::fabs(mcParticle.vz()) * std::tan(2 * std::atan(std::exp(-max_eta_gen_secondary))) - margin_z_gen) { continue; } - if (mp.pdgCode() == 22 && (mp.isPhysicalPrimary() || mp.producedByGenerator()) && std::fabs(mp.eta()) < max_eta_gen_secondary) { + if (motherParticle.pdgCode() == PDG_t::kGamma && (motherParticle.isPhysicalPrimary() || motherParticle.producedByGenerator()) && std::fabs(motherParticle.eta()) < max_eta_gen_secondary) { // if the MC truth particle corresponding to this reconstructed track which is not already written, add it to the skimmed MC stack if (!(fNewLabels.find(mcParticle.globalIndex()) != fNewLabels.end())) { // store electron information. !!Not photon!! fNewLabels[mcParticle.globalIndex()] = fCounters[0]; fNewLabelsReversed[fCounters[0]] = mcParticle.globalIndex(); // fMCFlags[mcParticle.globalIndex()] = mcflags; - fEventIdx[mcParticle.globalIndex()] = fEventLabels.find(mcCollision.globalIndex())->second; + fEventIdx[mcParticle.globalIndex()] = fEventLabels.find(mcCollisionIter.globalIndex())->second; fCounters[0]++; } // if the MC truth particle corresponding to this reconstructed track which is not already written, add it to the skimmed MC stack - if (!(fNewLabels.find(mp.globalIndex()) != fNewLabels.end())) { // store conversion photon - fNewLabels[mp.globalIndex()] = fCounters[0]; - fNewLabelsReversed[fCounters[0]] = mp.globalIndex(); - // fMCFlags[mp.globalIndex()] = mcflags; - fEventIdx[mp.globalIndex()] = fEventLabels.find(mcCollision.globalIndex())->second; + if (!(fNewLabels.find(motherParticle.globalIndex()) != fNewLabels.end())) { // store conversion photon + fNewLabels[motherParticle.globalIndex()] = fCounters[0]; + fNewLabelsReversed[fCounters[0]] = motherParticle.globalIndex(); + // fMCFlags[motherParticle.globalIndex()] = mcflags; + fEventIdx[motherParticle.globalIndex()] = fEventLabels.find(mcCollisionIter.globalIndex())->second; fCounters[0]++; } } @@ -227,59 +248,66 @@ struct AssociateMCInfoPhoton { } // end of rec. collision loop - if constexpr (static_cast(system & kPCM)) { + if constexpr (static_cast(system & kPCM) && !std::is_same_v) { + // electron and positron iterators from the TPCMLegs table as well as the o2Track table + auto ele = legs.begin(); + auto pos = legs.begin(); + + auto o2TrackEle = o2tracks.begin(); + auto o2TrackPos = o2tracks.begin(); + auto o2TrackIter = o2tracks.begin(); for (const auto& v0 : v0photons) { - auto collisionFromV0 = collisions.iteratorAt(v0.collisionId()); - if (!collisionFromV0.has_mcCollision()) { + collisionIter.setCursor(v0.collisionId()); + if (!collisionIter.has_mcCollision()) { continue; } - auto mcCollisionFromV0 = collisionFromV0.mcCollision(); + mcCollisionIter.setCursor(collisionIter.mcCollisionId()); - auto ele = v0.template negTrack_as(); - auto pos = v0.template posTrack_as(); + ele.setCursor(v0.negTrackId()); + pos.setCursor(v0.posTrackId()); - auto o2TrackEle = o2tracks.iteratorAt(pos.trackId()); - auto o2TrackPos = o2tracks.iteratorAt(ele.trackId()); + o2TrackEle.setCursor(ele.trackId()); + o2TrackPos.setCursor(pos.trackId()); if (!o2TrackEle.has_mcParticle() || !o2TrackPos.has_mcParticle()) { continue; // If no MC particle is found, skip the v0 } for (const auto& leg : {pos, ele}) { // be carefull of order {pos, ele}! - auto o2track = o2tracks.iteratorAt(leg.trackId()); - auto mcParticle = o2track.template mcParticle_as(); - // LOGF(info, "mcParticle.globalIndex() = %d, mcParticle.index() = %d", mcParticle.globalIndex(), mcParticle.index()); // these are exactly the same. + o2TrackIter.setCursor(leg.trackId()); + mcParticleIter.setCursor(o2TrackIter.mcParticleId()); + // LOGF(info, "mcParticleIter.globalIndex() = %d, mcParticleIter.index() = %d", mcParticleIter.globalIndex(), mcParticleIter.index()); // these are exactly the same. // if the MC truth particle corresponding to this reconstructed track which is not already written, add it to the skimmed MC stack - if (!(fNewLabels.find(mcParticle.globalIndex()) != fNewLabels.end())) { - fNewLabels[mcParticle.globalIndex()] = fCounters[0]; - fNewLabelsReversed[fCounters[0]] = mcParticle.globalIndex(); - // fMCFlags[mcParticle.globalIndex()] = mcflags; - fEventIdx[mcParticle.globalIndex()] = fEventLabels.find(mcCollisionFromV0.globalIndex())->second; + if (!(fNewLabels.find(mcParticleIter.globalIndex()) != fNewLabels.end())) { + fNewLabels[mcParticleIter.globalIndex()] = fCounters[0]; + fNewLabelsReversed[fCounters[0]] = mcParticleIter.globalIndex(); + // fMCFlags[mcParticleIter.globalIndex()] = mcflags; + fEventIdx[mcParticleIter.globalIndex()] = fEventLabels.find(mcCollisionIter.globalIndex())->second; fCounters[0]++; } - v0legmclabels(fNewLabels.find(mcParticle.index())->second, o2track.mcMask()); + v0legmclabels(fNewLabels.find(mcParticleIter.index())->second, o2TrackIter.mcMask()); // Next, store mother-chain of this reconstructed track. int motherid = -999; // first mother index - if (mcParticle.has_mothers()) { - motherid = mcParticle.mothersIds()[0]; // first mother index + if (mcParticleIter.has_mothers()) { + motherid = mcParticleIter.mothersIds()[0]; // first mother index } while (motherid > -1) { if (motherid < mcParticles.size()) { // protect against bad mother indices. why is this needed? - auto mp = mcParticles.iteratorAt(motherid); + motherParticle.setCursor(motherid); // if the MC truth particle corresponding to this reconstructed track which is not already written, add it to the skimmed MC stack - if (!(fNewLabels.find(mp.globalIndex()) != fNewLabels.end())) { - fNewLabels[mp.globalIndex()] = fCounters[0]; - fNewLabelsReversed[fCounters[0]] = mp.globalIndex(); - // fMCFlags[mp.globalIndex()] = mcflags; - fEventIdx[mp.globalIndex()] = fEventLabels.find(mcCollisionFromV0.globalIndex())->second; + if (!(fNewLabels.find(motherParticle.globalIndex()) != fNewLabels.end())) { + fNewLabels[motherParticle.globalIndex()] = fCounters[0]; + fNewLabelsReversed[fCounters[0]] = motherParticle.globalIndex(); + // fMCFlags[motherParticle.globalIndex()] = mcflags; + fEventIdx[motherParticle.globalIndex()] = fEventLabels.find(mcCollisionIter.globalIndex())->second; fCounters[0]++; } - if (mp.has_mothers()) { - motherid = mp.mothersIds()[0]; // first mother index + if (motherParticle.has_mothers()) { + motherid = motherParticle.mothersIds()[0]; // first mother index } else { motherid = -999; } @@ -292,50 +320,51 @@ struct AssociateMCInfoPhoton { } if constexpr (static_cast(system & kElectron)) { + auto o2TrackIter = o2tracks.begin(); // auto emprimaryelectrons_coll = emprimaryelectrons.sliceBy(perCollisionEl, collision.globalIndex()); for (const auto& emprimaryelectron : emprimaryelectrons) { - auto collisionFromEl = collisions.iteratorAt(emprimaryelectron.collisionId()); - if (!collisionFromEl.has_mcCollision()) { + collisionIter.setCursor(emprimaryelectron.collisionId()); + if (!collisionIter.has_mcCollision()) { continue; } - auto mcCollisionFromEl = collisionFromEl.mcCollision(); + mcCollisionIter.setCursor(collisionIter.mcCollisionId()); - auto o2track = o2tracks.iteratorAt(emprimaryelectron.trackId()); - if (!o2track.has_mcParticle()) { + o2TrackIter.setCursor(emprimaryelectron.trackId()); + if (!o2TrackIter.has_mcParticle()) { continue; // If no MC particle is found, skip the dilepton } - auto mcParticle = o2track.template mcParticle_as(); + mcParticleIter.setCursor(o2TrackIter.mcParticleId()); // if the MC truth particle corresponding to this reconstructed track which is not already written, add it to the skimmed MC stack - if (!(fNewLabels.find(mcParticle.globalIndex()) != fNewLabels.end())) { - fNewLabels[mcParticle.globalIndex()] = fCounters[0]; - fNewLabelsReversed[fCounters[0]] = mcParticle.globalIndex(); - // fMCFlags[mcParticle.globalIndex()] = mcflags; - fEventIdx[mcParticle.globalIndex()] = fEventLabels.find(mcCollisionFromEl.globalIndex())->second; + if (!(fNewLabels.find(mcParticleIter.globalIndex()) != fNewLabels.end())) { + fNewLabels[mcParticleIter.globalIndex()] = fCounters[0]; + fNewLabelsReversed[fCounters[0]] = mcParticleIter.globalIndex(); + // fMCFlags[mcParticleIter.globalIndex()] = mcflags; + fEventIdx[mcParticleIter.globalIndex()] = fEventLabels.find(mcCollisionIter.globalIndex())->second; fCounters[0]++; } - emprimaryelectronmclabels(fNewLabels.find(mcParticle.index())->second, o2track.mcMask()); + emprimaryelectronmclabels(fNewLabels.find(mcParticleIter.index())->second, o2TrackIter.mcMask()); // Next, store mother-chain of this reconstructed track. int motherid = -999; // first mother index - if (mcParticle.has_mothers()) { - motherid = mcParticle.mothersIds()[0]; // first mother index + if (mcParticleIter.has_mothers()) { + motherid = mcParticleIter.mothersIds()[0]; // first mother index } while (motherid > -1) { if (motherid < mcParticles.size()) { // protect against bad mother indices. why is this needed? - auto mp = mcParticles.iteratorAt(motherid); + motherParticle.setCursor(motherid); // if the MC truth particle corresponding to this reconstructed track which is not already written, add it to the skimmed MC stack - if (!(fNewLabels.find(mp.globalIndex()) != fNewLabels.end())) { - fNewLabels[mp.globalIndex()] = fCounters[0]; - fNewLabelsReversed[fCounters[0]] = mp.globalIndex(); - // fMCFlags[mp.globalIndex()] = mcflags; - fEventIdx[mp.globalIndex()] = fEventLabels.find(mcCollisionFromEl.globalIndex())->second; + if (!(fNewLabels.find(motherParticle.globalIndex()) != fNewLabels.end())) { + fNewLabels[motherParticle.globalIndex()] = fCounters[0]; + fNewLabelsReversed[fCounters[0]] = motherParticle.globalIndex(); + // fMCFlags[motherParticle.globalIndex()] = mcflags; + fEventIdx[motherParticle.globalIndex()] = fEventLabels.find(mcCollisionIter.globalIndex())->second; fCounters[0]++; } - if (mp.has_mothers()) { - motherid = mp.mothersIds()[0]; // first mother index + if (motherParticle.has_mothers()) { + motherid = motherParticle.mothersIds()[0]; // first mother index } else { motherid = -999; } @@ -350,62 +379,74 @@ struct AssociateMCInfoPhoton { if constexpr (static_cast(system & kEMC)) { // for emc photons // auto ememcclusters_coll = emcclusters.sliceBy(perCollisionEMC, collision.globalIndex()); for (const auto& emccluster : emcclusters) { - auto collisionFromEMC = collisions.iteratorAt(emccluster.collisionId()); - if (!collisionFromEMC.has_mcCollision()) { + collisionIter.setCursor(emccluster.collisionId()); + if (!collisionIter.has_mcCollision()) { continue; } - auto mcCollisionFromEMC = collisionFromEMC.mcCollision(); - - // TODO: check that this does not cause issues, since we have sometimes clusters that are only noise! - auto mcphoton = mcParticles.iteratorAt(emccluster.emmcparticleIds()[0]); + mcCollisionIter.setCursor(collisionIter.mcCollisionId()); - // if the MC truth particle corresponding to this reconstructed track which is not already written, add it to the skimmed MC stack - if (!(fNewLabels.find(mcphoton.globalIndex()) != fNewLabels.end())) { - fNewLabels[mcphoton.globalIndex()] = fCounters[0]; - fNewLabelsReversed[fCounters[0]] = mcphoton.globalIndex(); - fEventIdx[mcphoton.globalIndex()] = fEventLabels.find(mcCollisionFromEMC.globalIndex())->second; - fCounters[0]++; + // TODO: test + if (emccluster.emmcparticleIds().size() <= 0) { + continue; } - ememcclustermclabels(fNewLabels.find(mcphoton.index())->second); + std::vector vEmcMcParticleIds; - // Next, store mother-chain of this reconstructed track. - int motherid = -999; // first mother index - if (mcphoton.has_mothers()) { - motherid = mcphoton.mothersIds()[0]; // first mother index - } - while (motherid > -1) { - if (motherid < mcParticles.size()) { // protect against bad mother indices. why is this needed? - auto mp = mcParticles.iteratorAt(motherid); + vEmcMcParticleIds.reserve(emccluster.emmcparticleIds().size()); - // if the MC truth particle corresponding to this reconstructed track which is not already written, add it to the skimmed MC stack - if (!(fNewLabels.find(mp.globalIndex()) != fNewLabels.end())) { - fNewLabels[mp.globalIndex()] = fCounters[0]; - fNewLabelsReversed[fCounters[0]] = mp.globalIndex(); - fEventIdx[mp.globalIndex()] = fEventLabels.find(mcCollisionFromEMC.globalIndex())->second; - fCounters[0]++; - } + for (const auto& emcParticleId : emccluster.emmcparticleIds()) { + mcPhoton.setCursor(emcParticleId); - if (mp.has_mothers()) { - motherid = mp.mothersIds()[0]; // first mother index + // if the MC truth particle corresponding to this reconstructed track which is not already written, add it to the skimmed MC stack + if (!(fNewLabels.find(mcPhoton.globalIndex()) != fNewLabels.end())) { + fNewLabels[mcPhoton.globalIndex()] = fCounters[0]; + fNewLabelsReversed[fCounters[0]] = mcPhoton.globalIndex(); + fEventIdx[mcPhoton.globalIndex()] = fEventLabels.find(mcCollisionIter.globalIndex())->second; + fCounters[0]++; + } + vEmcMcParticleIds.emplace_back(fNewLabels.find(mcPhoton.index())->second); + // ememcclustermclabels(fNewLabels.find(mcPhoton.index())->second); + + // Next, store mother-chain of this reconstructed track. + int motherid = -999; // first mother index + if (mcPhoton.has_mothers()) { + motherid = mcPhoton.mothersIds()[0]; // first mother index + } + while (motherid > -1) { + if (motherid < mcParticles.size()) { // protect against bad mother indices. why is this needed? + motherParticle.setCursor(motherid); + + // if the MC truth particle corresponding to this reconstructed track which is not already written, add it to the skimmed MC stack + if (!(fNewLabels.find(motherParticle.globalIndex()) != fNewLabels.end())) { + fNewLabels[motherParticle.globalIndex()] = fCounters[0]; + fNewLabelsReversed[fCounters[0]] = motherParticle.globalIndex(); + fEventIdx[motherParticle.globalIndex()] = fEventLabels.find(mcCollisionIter.globalIndex())->second; + fCounters[0]++; + } + + if (motherParticle.has_mothers()) { + motherid = motherParticle.mothersIds()[0]; // first mother index + } else { + motherid = -999; + } } else { motherid = -999; } - } else { - motherid = -999; - } - } // end of mother chain loop + } // end of mother chain loop + + } // end of loop over mc particles of the current emc cluster + ememcclustermclabels(vEmcMcParticleIds); } // end of em emc cluster loop } // Loop over the label map, create the mother/daughter relationships if these exist and write the skimmed MC stack for (const auto& [newLabel, oldLabel] : fNewLabelsReversed) { - auto mcParticle = mcParticles.iteratorAt(oldLabel); + mcParticleIter.setCursor(oldLabel); // uint16_t mcflags = fMCFlags.find(oldLabel)->second; std::vector mothers; - if (mcParticle.has_mothers()) { - for (const auto& m : mcParticle.mothersIds()) { + if (mcParticleIter.has_mothers()) { + for (const auto& m : mcParticleIter.mothersIds()) { if (m < mcParticles.size()) { // protect against bad mother indices if (fNewLabels.find(m) != fNewLabels.end()) { mothers.push_back(fNewLabels.find(m)->second); @@ -419,9 +460,9 @@ struct AssociateMCInfoPhoton { // Note that not all daughters from the original table are preserved in the skimmed MC stack std::vector daughters; - if (mcParticle.has_daughters()) { - // LOGF(info, "daughter range in original MC stack pdg = %d | %d - %d , n dau = %d", mcParticle.pdgCode(), mcParticle.daughtersIds()[0], mcParticle.daughtersIds()[1], mcParticle.daughtersIds()[1] -mcParticle.daughtersIds()[0] +1); - for (int d = mcParticle.daughtersIds()[0]; d <= mcParticle.daughtersIds()[1]; ++d) { + if (mcParticleIter.has_daughters()) { + // LOGF(info, "daughter range in original MC stack pdg = %d | %d - %d , n dau = %d", mcParticleIter.pdgCode(), mcParticleIter.daughtersIds()[0], mcParticleIter.daughtersIds()[1], mcParticleIter.daughtersIds()[1] -mcParticleIter.daughtersIds()[0] +1); + for (int d = mcParticleIter.daughtersIds()[0]; d <= mcParticleIter.daughtersIds()[1]; ++d) { // TODO: remove this check as soon as issues with MC production are fixed if (d < mcParticles.size()) { // protect against bad daughter indices // auto dau_tmp = mcParticles.iteratorAt(d); @@ -430,16 +471,16 @@ struct AssociateMCInfoPhoton { daughters.push_back(fNewLabels.find(d)->second); } } else { - LOG(info) << "Daughter label (" << d << ") exceeds the McParticles size (" << mcParticles.size() << ")"; - LOG(info) << " Check the MC generator"; + LOG(error) << "Daughter label (" << d << ") exceeds the McParticles size (" << mcParticles.size() << ")"; + LOG(error) << " Check the MC generator"; } } } - emmcparticles(fEventIdx.find(oldLabel)->second, mcParticle.pdgCode(), mcParticle.flags(), mcParticle.statusCode(), + emmcparticles(fEventIdx.find(oldLabel)->second, mcParticleIter.pdgCode(), mcParticleIter.flags(), mcParticleIter.statusCode(), mothers, daughters, - mcParticle.px(), mcParticle.py(), mcParticle.pz(), mcParticle.e(), - mcParticle.vx(), mcParticle.vy(), mcParticle.vz()); + mcParticleIter.px(), mcParticleIter.py(), mcParticleIter.pz(), mcParticleIter.e(), + mcParticleIter.vx(), mcParticleIter.vy(), mcParticleIter.vz()); } // end loop over labels fNewLabels.clear(); diff --git a/PWGEM/PhotonMeson/TableProducer/bcWiseClusterSkimmer.cxx b/PWGEM/PhotonMeson/TableProducer/bcWiseClusterSkimmer.cxx index 72eb21b1ad4..dbc8c80c7d2 100644 --- a/PWGEM/PhotonMeson/TableProducer/bcWiseClusterSkimmer.cxx +++ b/PWGEM/PhotonMeson/TableProducer/bcWiseClusterSkimmer.cxx @@ -10,9 +10,7 @@ // or submit itself to any jurisdiction. /// /// \file bcWiseClusterSkimmer.cxx -/// /// \brief This task creates minimalistic skimmed tables containing EMC clusters and centrality information -/// /// \author Nicolas Strangmann (nicolas.strangmann@cern.ch) - Goethe University Frankfurt /// @@ -24,12 +22,15 @@ #include "Common/DataModel/Centrality.h" #include "Common/DataModel/EventSelection.h" -#include "CCDB/BasicCCDBManager.h" -#include "DataFormatsParameters/GRPLHCIFData.h" -#include "DetectorsBase/GeometryManager.h" -#include "EMCALBase/Geometry.h" -#include "Framework/AnalysisTask.h" -#include "Framework/runDataProcessing.h" +#include +#include +#include +#include +#include +#include +#include + +#include #include #include @@ -40,6 +41,7 @@ using namespace o2; using namespace o2::aod::emdownscaling; using namespace o2::framework; using namespace o2::framework::expressions; +using namespace o2::constants::physics; using MyCollisions = soa::Join; using MyMCCollisions = soa::Join; @@ -150,13 +152,13 @@ struct bcWiseClusterSkimmer { { for (const auto& cluster : clusters) { clusterTable(bcID, - convertForStorage(cluster.definition(), kDefinition), - convertForStorage(cluster.energy(), kEnergy), - convertForStorage(cluster.eta(), kEta), - convertForStorage(cluster.phi(), kPhi), - convertForStorage(cluster.nCells(), kNCells), - convertForStorage(cluster.m02(), kM02), - convertForStorage(cluster.time(), kTime), + convertForStorage(cluster.definition(), Observable::kDefinition), + convertForStorage(cluster.energy(), Observable::kEnergy), + convertForStorage(cluster.eta(), Observable::kEta), + convertForStorage(cluster.phi(), Observable::kPhi), + convertForStorage(cluster.nCells(), Observable::kNCells), + convertForStorage(cluster.m02(), Observable::kM02), + convertForStorage(cluster.time(), Observable::kTime), cluster.isExotic()); } } @@ -184,12 +186,12 @@ struct bcWiseClusterSkimmer { } bool isEta = false; if (mesonMCIndex >= 0) { - if (mcParticles.iteratorAt(mesonMCIndex).pdgCode() == 111) { + if (mcParticles.iteratorAt(mesonMCIndex).pdgCode() == PDG_t::kPi0) { if (fMapPi0Index.find(mesonMCIndex) != fMapPi0Index.end()) // Some pi0s might not be found (not gg decay or too large y) mesonMCIndex = fMapPi0Index[mesonMCIndex]; // If pi0 was stored in table, change index from the MC index to the pi0 index from this task else // If pi0 was not stored, treat photon as if not from pi0 mesonMCIndex = -1; - } else if (mcParticles.iteratorAt(mesonMCIndex).pdgCode() == 221) { + } else if (mcParticles.iteratorAt(mesonMCIndex).pdgCode() == Pdg::kEta) { isEta = true; if (fMapEtaIndex.find(mesonMCIndex) != fMapEtaIndex.end()) // Some etas might not be found (not gg decay or too large y) mesonMCIndex = fMapEtaIndex[mesonMCIndex]; // If eta was stored in table, change index from the MC index to the eta index from this task @@ -305,7 +307,7 @@ struct bcWiseClusterSkimmer { if (daughtersIds.size() != 2) return false; for (const auto& daughterId : daughtersIds) { - if (mcParticles.iteratorAt(daughterId).pdgCode() != 22) + if (mcParticles.iteratorAt(daughterId).pdgCode() != PDG_t::kGamma) return false; } return true; @@ -318,7 +320,7 @@ struct bcWiseClusterSkimmer { if (daughtersIds.size() != 2) return false; for (const auto& daughterId : daughtersIds) { - if (mcParticles.iteratorAt(daughterId).pdgCode() != 22) + if (mcParticles.iteratorAt(daughterId).pdgCode() != PDG_t::kGamma) return false; int iCellID = -1; try { @@ -380,10 +382,10 @@ struct bcWiseClusterSkimmer { bool isPrimary = mcParticle.isPhysicalPrimary() || mcParticle.producedByGenerator(); bool isFromWD = (aod::pwgem::photonmeson::utils::mcutil::IsFromWD(mcCollision, mcParticle, mcParticles)) > 0; - if (mcParticle.pdgCode() == 111) { + if (mcParticle.pdgCode() == PDG_t::kPi0) { mcpi0Table(bcTable.lastIndex(), convertForStorage(mcParticle.pt(), kpT), isAccepted(mcParticle, mcParticles), isPrimary, isFromWD); fMapPi0Index[mcParticle.globalIndex()] = static_cast(mcpi0Table.lastIndex()); - } else if (mcParticle.pdgCode() == 221) { + } else if (mcParticle.pdgCode() == Pdg::kEta) { mcetaTable(bcTable.lastIndex(), convertForStorage(mcParticle.pt(), kpT), isAccepted(mcParticle, mcParticles), isPrimary, isFromWD); fMapEtaIndex[mcParticle.globalIndex()] = static_cast(mcetaTable.lastIndex()); } diff --git a/PWGEM/PhotonMeson/Tasks/emcalBcWiseGammaGamma.cxx b/PWGEM/PhotonMeson/Tasks/emcalBcWiseGammaGamma.cxx index efbc38d799c..93e89d45c54 100644 --- a/PWGEM/PhotonMeson/Tasks/emcalBcWiseGammaGamma.cxx +++ b/PWGEM/PhotonMeson/Tasks/emcalBcWiseGammaGamma.cxx @@ -11,26 +11,23 @@ // /// /// \file emcalBcWiseGammaGamma.cxx -/// /// \brief Task that extracts pi0s and eta mesons from BC wise derived data of EMCal clusters -/// /// \author Nicolas Strangmann (nicolas.strangmann@cern.ch) Goethe University Frankfurt /// #include "PWGEM/PhotonMeson/DataModel/bcWiseTables.h" -#include "CommonConstants/MathConstants.h" -#include "EMCALBase/Geometry.h" -#include "Framework/AnalysisTask.h" -#include "Framework/HistogramRegistry.h" -#include "Framework/runDataProcessing.h" - -#include "Math/AxisAngle.h" -#include "Math/LorentzRotation.h" -#include "Math/Rotation3D.h" -#include "Math/Vector3D.h" -#include "Math/Vector4D.h" -#include "TString.h" +#include +#include +#include +#include +#include + +#include +#include +#include // IWYU pragma: keep +#include +#include using namespace o2; using namespace o2::framework; diff --git a/PWGEM/PhotonMeson/Tasks/photonResoTask.cxx b/PWGEM/PhotonMeson/Tasks/photonResoTask.cxx index cc1fe873798..2c6cecf0a27 100644 --- a/PWGEM/PhotonMeson/Tasks/photonResoTask.cxx +++ b/PWGEM/PhotonMeson/Tasks/photonResoTask.cxx @@ -23,24 +23,15 @@ #include "PWGEM/PhotonMeson/Utils/EventHistograms.h" #include "PWGEM/PhotonMeson/Utils/MCUtilities.h" -#include "Common/DataModel/EventSelection.h" - #include #include #include #include -#include -#include -#include #include #include -#include #include #include -#include #include -#include -#include #include #include #include @@ -48,14 +39,11 @@ #include #include -#include -#include #include // IWYU pragma: keep #include #include #include #include -#include #include #include @@ -179,11 +167,7 @@ struct PhotonResoTask { struct : ConfigurableGroup { std::string prefix = "mesonConfig"; Configurable minOpenAngle{"minOpenAngle", 0.0202, "apply min opening angle. Default value one EMCal cell"}; - Configurable enableTanThetadPhi{"enableTanThetadPhi", false, "flag to turn cut opening angle in delta theta delta phi on/off"}; - Configurable minTanThetadPhi{"minTanThetadPhi", 4., "apply min opening angle in delta theta delta phi to cut on late conversion"}; - Configurable maxEnergyAsymmetry{"maxEnergyAsymmetry", 1., "apply max energy asymmetry for meson candidate"}; Configurable cfgEnableQA{"cfgEnableQA", false, "flag to turn QA plots on/off"}; - ConfigurableAxis thConfigAxisTanThetaPhi{"thConfigAxisTanThetaPhi", {180, -90.f, 90.f}, ""}; } mesonConfig; struct : ConfigurableGroup { @@ -213,15 +197,12 @@ struct PhotonResoTask { SliceCache cache; - Filter collisionFilter = (nabs(aod::collision::posZ) <= eventcuts.cfgZvtxMax) && (aod::evsel::ft0cOccupancyInTimeRange <= eventcuts.cfgFT0COccupancyMax) && (aod::evsel::ft0cOccupancyInTimeRange >= eventcuts.cfgFT0COccupancyMin); - using EMCalPhotons = soa::Join; using PcmPhotons = soa::Join; using PcmMcLegs = soa::Join; using Colls = soa::Join; - using FilteredColls = soa::Filtered; using McColls = o2::soa::Join; using McParticles = EMMCParticles; @@ -446,7 +427,7 @@ struct PhotonResoTask { fEMCCut.AreSelectedRunning(emcFlags, clusters, matchedPrims, matchedSeconds, ®istry); EMBitFlags v0flags(photons.size()); - fV0PhotonCut.AreSelectedRunning(v0flags, photons, ®istry); + fV0PhotonCut.AreSelectedRunning(v0flags, photons, ®istry); // create iterators for photon mc particles auto mcPhoton1 = mcParticles.begin(); @@ -473,7 +454,12 @@ struct PhotonResoTask { if (!(emcFlags.test(photonEMC.globalIndex()))) { continue; } - mcPhoton1.setCursor(photonEMC.emmcparticleId()); + if (photonEMC.emmcparticleIds().size() <= 0) { + // this is a cluster with just noise, skip + continue; + } + // we only want to look at the largest contribution + mcPhoton1.setCursor(photonEMC.emmcparticleIds()[0]); if (std::abs(mcPhoton1.pdgCode()) == PDG_t::kGamma) { registry.fill(HIST("EMCal/hPhotonReso"), photonEMC.pt(), mcPhoton1.pt(), cent); @@ -502,9 +488,11 @@ struct PhotonResoTask { continue; } + float trueConvRadius = std::hypot(ele1mc.vx(), ele1mc.vy()); + mcPhoton1.setCursor(photonid1); - if (!fV0PhotonCut.IsConversionPointInAcceptance(mcPhoton1)) { + if (!fV0PhotonCut.IsConversionPointInAcceptance(mcPhoton1, trueConvRadius)) { continue; } @@ -547,8 +535,13 @@ struct PhotonResoTask { registry.fill(HIST("mesonQA/hInvMassPt"), vMeson.M(), vMeson.Pt()); } - mcPhoton1.setCursor(g1.emmcparticleId()); - mcPhoton2.setCursor(g2.emmcparticleId()); + if (g1.emmcparticleIds().size() <= 0 || g2.emmcparticleIds().size() <= 0) { + // there is a cluster which is just noise, skip + continue; + } + + mcPhoton1.setCursor(g1.emmcparticleIds()[0]); + mcPhoton2.setCursor(g2.emmcparticleIds()[0]); int photonid1 = -1, photonid2 = -1, pi0id = -1, etaid = -1; photonid1 = o2::aod::pwgem::photonmeson::utils::mcutil::FindMotherInChain(mcPhoton1, mcParticles, std::vector{111, 221});