diff --git a/PWGEM/PhotonMeson/Tasks/twoPhotonMcTask.cxx b/PWGEM/PhotonMeson/Tasks/twoPhotonMcTask.cxx index 6b2e5fd2aa7..e2b366714fe 100644 --- a/PWGEM/PhotonMeson/Tasks/twoPhotonMcTask.cxx +++ b/PWGEM/PhotonMeson/Tasks/twoPhotonMcTask.cxx @@ -119,6 +119,7 @@ enum class LegMatchStatus : uint8_t { 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 + kNotGlobalTrack, // primaries only: the leg is not a global track, so it is no candidate of the primary matching kNStatus }; @@ -478,12 +479,15 @@ struct TwoPhotonMcTask { 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"}; + 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", "not a global track"}; 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()); } + // labels make the axis extendable, and TH1Merger::LabelMerge in ROOT 6.36 reads the z labels from the y axis: + // merged TH3 content would end up in a new, unlabelled z bin. With a fixed axis the merge goes bin by bin. + hOrigin->GetZaxis()->SetCanExtend(false); 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}); @@ -762,6 +766,30 @@ struct TwoPhotonMcTask { return PairBucket::kCrossOtherHadron; } + /// \brief Whether the track of a V0 leg passes the global track selection used by the primary EMCal track matching. + /// Rebuilt from the leg columns of the derived data, following getGlobalTrackSelectionRun3ITSMatch(Run3ITSibAny): + /// the golden chi2 requirement cannot be checked here, and pT/eta are the ones of the leg at the V0 vertex. + /// \param leg V0 leg + template + static bool isLikelyGlobalTrack(TLeg const& leg) + { + constexpr float MinPt = 0.1f, MaxEta = 0.8f; + constexpr int MinCrossedRows = 70; + constexpr float MinCrossedRowsOverFindable = 0.8f, MaxChi2TPC = 4.f, MaxChi2ITS = 36.f, MaxDcaZ = 2.f; + const float pt = leg.pt(); + if (!leg.hasITS() || !leg.hasTPC() || pt < MinPt || std::fabs(leg.eta()) > MaxEta) { + return false; + } + if (leg.itsNClsInnerBarrel() < 1 || leg.itsChi2NCl() > MaxChi2ITS) { // Run3ITSibAny: at least one hit in ITS layers 0-2 + return false; + } + if (leg.tpcNClsCrossedRows() < MinCrossedRows || leg.tpcCrossedRowsOverFindableCls() < MinCrossedRowsOverFindable || leg.tpcChi2NCl() > MaxChi2TPC) { + return false; + } + const float maxDcaXY = 0.0105f + 0.035f / std::pow(pt, 1.1f); // o2-linter: disable=magic-number (global track selection) + return std::fabs(leg.dcaXY()) <= maxDcaXY && std::fabs(leg.dcaZ()) <= MaxDcaZ; + } + /// \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 @@ -865,7 +893,12 @@ struct TwoPhotonMcTask { return; } registry.fill(HIST("PCMEMC/SameParticle/hLegEoverP"), emc.e() / pLeg, ptLeg); - fillLegMatch(matchedPrims.sliceBy(perEMCClusterMT, emc.globalIndex()), emc.e(), pLeg, ptLeg); + // the primary matching only uses global tracks: other legs are no candidate there and get their own status + if (isLikelyGlobalTrack(leg)) { + fillLegMatch(matchedPrims.sliceBy(perEMCClusterMT, emc.globalIndex()), emc.e(), pLeg, ptLeg); + } else { + registry.fill(HIST("PCMEMC/SameParticle/hLegPrimStatus"), static_cast(LegMatchStatus::kNotGlobalTrack), ptLeg); + } fillLegMatch(matchedSeconds.sliceBy(perEMCClusterMS, emc.globalIndex()), emc.e(), pLeg, ptLeg); } diff --git a/PWGJE/Tasks/emcTmMonitor.cxx b/PWGJE/Tasks/emcTmMonitor.cxx index 3d6baefda3b..cef13ebaf16 100644 --- a/PWGJE/Tasks/emcTmMonitor.cxx +++ b/PWGJE/Tasks/emcTmMonitor.cxx @@ -22,7 +22,11 @@ /// Simple event selection using the flag doEventSel is provided, which selects INT7 events if set to 1 /// For pilot beam data, instead of relying on the event selection, one can veto specific BC IDS using the flag /// fDoVetoBCID. +/// The MC process (processMCTruth) follows the leading MC particle of each cluster through the matching chain: +/// is its track reconstructed, extrapolated to the EMCal, a matching candidate (global track or V0 leg), +/// and matched to the cluster as primary or secondary track. +#include "PWGEM/PhotonMeson/DataModel/gammaTables.h" // for the V0 legs used in the secondary matching #include "PWGJE/DataModel/EMCALClusters.h" #include "Common/CCDB/TriggerAliases.h" @@ -30,24 +34,40 @@ #include "Common/DataModel/PIDResponseTPC.h" #include "Common/DataModel/TrackSelectionTables.h" +#include #include #include #include #include +#include #include +#include #include #include #include #include #include +#include #include #include -#include +#include +#include +#include +#include +#include +#include + +#include #include +#include +#include #include #include #include +#include +#include +#include #include using namespace o2::framework; @@ -56,16 +76,55 @@ using CollisionEvSelIt = o2::soa::Join::it using BcEvSelIt = o2::soa::Join::iterator; using SelectedClusters = o2::soa::Filtered; using TracksPID = o2::soa::Join; +using SelectedMcClusters = o2::soa::Filtered>; +using TracksMC = o2::soa::Join; + +/// \brief Truth class of the leading MC particle of a cluster +enum class TmTruthClass : uint8_t { + kConversionElectron = 0, // e+- with a photon as mother, produced in the transport (photon conversion) + kPrimaryElectron, // e+- produced by the generator (Dalitz, heavy flavour, ...) + kOtherElectron, // any other e+- (delta rays, Compton, weak decays in the transport, ...) + kOtherCharged, // charged particle that is not an e+- (hadrons, muons) + kPhoton, + kNeutral, // neutral particle that is not a photon + kOther, // no MC particle or unknown PDG code + kNClasses +}; + +/// \brief What happened to the track of the leading MC particle of a cluster in the track matching +enum class TmTruthStatus : uint8_t { + kNoRecoTrack = 0, // no reconstructed track of this particle in this collision + kNotExtrapolated, // track(s) exist, but none was extrapolated to the EMCal (trackEta/PhiEmcal = -999) + kNoMatchingCandidate, // extrapolated, but neither a global track nor a V0 leg: input to no matching + kCandidateNotMatched, // matching candidate, but not matched to this cluster (dR > maxMatchingDistance or not among the closest tracks) + kMatchedPrimary, // matched to this cluster as primary (global) track + kMatchedSecondary, // matched to this cluster as secondary track (V0 leg) + kNStatus +}; + +/// \brief Type of the reconstructed track +enum class TmTrackType : uint8_t { + kNoTrack = 0, + kGlobal, // passes isGlobalTrack() + kItsTpc, // ITS and TPC, but not a global track + kTpcOnly, + kItsOnly, + kOtherType, + kNTypes +}; + struct EmcTmMonitor { HistogramRegistry mHistManager{"TrackMatchingMonitorHistograms", {}, OutputObjHandlingPolicy::AnalysisObject}; o2::emcal::Geometry* mGeometry = nullptr; + Service mPdgDatabase{}; Preslice perCluster = o2::aod::emcalclustercell::emcalclusterId; Preslice perClusterAmb = o2::aod::emcalclustercell::emcalambiguousclusterId; Preslice perClusterMatchedTracks = o2::aod::emcalclustercell::emcalclusterId; + Preslice perClusterMatchedSecondaries = o2::aod::emcalclustercell::emcalclusterId; // configurable parameters // TODO adapt doEventSel switch to also allow selection of other triggers (e.g. EMC7) - Configurable doEventSel{"doEventSel", 0, "demand kINT7"}; + Configurable doEventSel{"doEventSel", false, "demand kINT7"}; Configurable vetoBCID{"vetoBCID", "", "BC ID(s) to be excluded, this should be used as an alternative to the event selection"}; Configurable selectBCID{"selectBCID", "all", "BC ID(s) to be included, this should be used as an alternative to the event selection"}; Configurable vertexCut{"vertexCut", -1, "apply z-vertex cut with value in cm"}; @@ -85,6 +144,7 @@ struct EmcTmMonitor { Configurable minDPhi{"minDPhi", 0.01, "Minimum dPhi between track and cluster"}; Configurable> eOverPRange{"eOverPRange", {0.9, 1.2}, "E/p range where one would search for electrons (first <= E/p <= second)"}; Configurable hasTRD{"hasTRD", false, "demand track to have a hit in TRD."}; + Configurable emcalRadius{"emcalRadius", 440.f, "MC truth: radius (cm) used to move the cluster eta from the vertex to the origin frame of the track extrapolation"}; std::vector mVetoBCIDs; std::vector mSelectBCIDs; @@ -182,10 +242,67 @@ struct EmcTmMonitor { mHistManager.add("clusterTM_EoverP_hadron_ASide", "cluster E/hadron p in A-Side", HistType::kTH3F, {eoverpAxis, trackpAxis, nmatchedtrack}); // E over p vs track pT for the Nth closest hadron track in A-Side mHistManager.add("clusterTM_EoverP_hadron_CSide", "cluster E/hadron p in C-Side", HistType::kTH3F, {eoverpAxis, trackpAxis, nmatchedtrack}); // E over p vs track pT for the Nth closest hadron track in C-Side - if (vetoBCID->length()) { + // MC truth of the track matching, only booked if processMCTruth is enabled + if (doprocessMCTruth.value) { + static constexpr std::array(TmTruthClass::kNClasses)> ClassLabels = {"conversion e^{#pm}", "primary e^{#pm}", "other e^{#pm}", "charged (non e)", "#gamma", "neutral (non #gamma)", "other"}; + static constexpr std::array(TmTruthStatus::kNStatus)> StatusLabels = {"no reco track", "not extrapolated", "no matching candidate", "candidate not matched", "matched primary", "matched secondary"}; + static constexpr std::array(TmTrackType::kNTypes)> TypeLabels = {"no track", "global", "ITS-TPC non-global", "TPC only", "ITS only", "other"}; + const AxisSpec classAxis{static_cast(TmTruthClass::kNClasses), -0.5, static_cast(TmTruthClass::kNClasses) - 0.5, "leading MC particle"}; + const AxisSpec statusAxis{static_cast(TmTruthStatus::kNStatus), -0.5, static_cast(TmTruthStatus::kNStatus) - 0.5, "track matching status"}; + const AxisSpec typeAxis{static_cast(TmTrackType::kNTypes), -0.5, static_cast(TmTrackType::kNTypes) - 0.5, "track type"}; + const AxisSpec v0LegAxis{2, -0.5, 1.5, "V0 leg"}; + const AxisSpec mcPtAxis{std::vector{0., 0.25, 0.5, 0.75, 1., 1.5, 2., 3., 4., 6., 10., 20., 50.}, "#it{p}_{T,MC} (GeV/#it{c})"}; + // residuals: fine binning around 0, coarse binning up to MaxResidual, so the tails beyond the matching windows are visible + constexpr double MaxResidual = 0.4; // range of the residual axes + constexpr double FineRange = 0.04; // fine binning within +-FineRange + constexpr double FineWidth = 0.002; // bin width within +-FineRange + constexpr double CoarseWidth = 0.01; // bin width outside +-FineRange + const int nCoarse = static_cast(std::lround((MaxResidual - FineRange) / CoarseWidth)); + const int nFine = static_cast(std::lround((FineRange - (-FineRange)) / FineWidth)); + std::vector residualBins; + residualBins.reserve(nCoarse + nFine + nCoarse + 1); // negative coarse, fine, positive coarse edges incl. the upper end + for (int i = 0; i < nCoarse; ++i) { + residualBins.push_back(-MaxResidual + CoarseWidth * i); + } + for (int i = 0; i < nFine; ++i) { + residualBins.push_back(-FineRange + FineWidth * i); + } + for (int i = 0; i <= nCoarse; ++i) { + residualBins.push_back(FineRange + CoarseWidth * i); + } + const AxisSpec mcDEtaAxis{residualBins, "#Delta#eta (track - cluster)"}; + const AxisSpec mcDPhiAxis{residualBins, "#Delta#varphi (track - cluster) (rad)"}; + const AxisSpec radiusAxis{125, 0., 250., "#it{R}_{xy} of the conversion (cm)"}; + + auto setLabels = [](TAxis* axis, auto const& labels) { + for (std::size_t i = 0; i < labels.size(); ++i) { + axis->SetBinLabel(static_cast(i) + 1, labels[i].data()); + } + }; + // histograms with labelled axes are TH2 or THn: the TH3 label merge of ROOT 6.36 is broken + auto hLeadingClass = mHistManager.add("MCTruth/hLeadingClass", "leading MC particle of the cluster", HistType::kTH2F, {classAxis, energyAxis}); + setLabels(hLeadingClass->GetXaxis(), ClassLabels); + auto hStatus = mHistManager.add("MCTruth/hStatus", "track matching status of the leading MC particle", HistType::kTHnF, {statusAxis, classAxis, typeAxis, v0LegAxis, mcPtAxis}); + setLabels(hStatus->GetAxis(0), StatusLabels); + setLabels(hStatus->GetAxis(1), ClassLabels); + setLabels(hStatus->GetAxis(2), TypeLabels); + for (const auto& name : {"MCTruth/hResidualEta", "MCTruth/hResidualPhi"}) { + const bool isEta = std::string_view(name).ends_with("Eta"); + auto hRes = mHistManager.add(name, "residual of the extrapolated track of the leading MC particle (origin frame)", HistType::kTHnF, {isEta ? mcDEtaAxis : mcDPhiAxis, mcPtAxis, classAxis, typeAxis}); + setLabels(hRes->GetAxis(2), ClassLabels); + setLabels(hRes->GetAxis(3), TypeLabels); + } + mHistManager.add("MCTruth/hResidualConversion", "conversion e^{#pm}: residual of the extrapolated track", HistType::kTH3F, {mcDEtaAxis, mcDPhiAxis, mcPtAxis}); + mHistManager.add("MCTruth/hResidualVsStored", "matched tracks: own minus stored residual;#Delta#eta_{own} - #Delta#eta_{stored};#Delta#varphi_{own} - #Delta#varphi_{stored}", HistType::kTH2F, {{200, -0.02, 0.02}, {200, -0.02, 0.02}}); + auto hConvRadius = mHistManager.add("MCTruth/hConvRadius", "conversion e^{#pm}: matching status vs conversion radius", HistType::kTHnF, {statusAxis, radiusAxis, mcPtAxis}); + setLabels(hConvRadius->GetAxis(0), StatusLabels); + mHistManager.add("MCTruth/hEoverPConversion", "conversion e^{#pm}: #it{E}_{cls}/#it{p}_{track}", HistType::kTH2F, {eoverpAxis, trackptAxis}); + } + + if (!vetoBCID->empty()) { std::stringstream parser(vetoBCID.value); std::string token; - int bcid; + int bcid = -1; while (std::getline(parser, token, ',')) { bcid = std::stoi(token); LOG(info) << "Veto BCID " << bcid; @@ -195,7 +312,7 @@ struct EmcTmMonitor { if (selectBCID.value != "all") { std::stringstream parser(selectBCID.value); std::string token; - int bcid; + int bcid = -1; while (std::getline(parser, token, ',')) { bcid = std::stoi(token); LOG(info) << "Select BCID " << bcid; @@ -209,19 +326,22 @@ struct EmcTmMonitor { // an integer instead Filter clusterDefinitionSelection = (o2::aod::emcalcluster::definition == clusterDefinition) && (o2::aod::emcalcluster::time >= minTime) && (o2::aod::emcalcluster::time <= maxTime) && (o2::aod::emcalcluster::m02 > minM02) && (o2::aod::emcalcluster::m02 < maxM02); - /// \brief Process EMCAL clusters that are matched to a collisions - void processCollisions(CollisionEvSelIt const& theCollision, SelectedClusters const& clusters, o2::aod::EMCALClusterCells const&, o2::aod::Calos const&, o2::aod::EMCALMatchedTracks const& matchedtracks, TracksPID const& alltracks) + /// \brief Event selection of the monitor: trigger alias (if doEventSel) and z vertex + /// \param theCollision collision to check + /// \param fillHistos fill the event histograms (only once if several process functions run) + bool isEventSelected(CollisionEvSelIt const& theCollision, const bool fillHistos) { - mHistManager.fill(HIST("eventsAll"), 1); + if (fillHistos) { + mHistManager.fill(HIST("eventsAll"), 1); + } // do event selection if doEventSel is specified // currently the event selection is hard coded to kINT7 // but other selections are possible that are defined in TriggerAliases.h bool isSelected = true; const int beginningRun3 = 300000; - float lowP = 0.75f; - float midP = 1.25f; if (doEventSel) { + // needs aod::BCs in the process function, otherwise the bc() index is not bound if (theCollision.bc().runNumber() < beginningRun3) { if (!theCollision.alias_bit(kINT7)) { isSelected = false; @@ -234,20 +354,34 @@ struct EmcTmMonitor { } if (!isSelected) { LOG(debug) << "Event not selected because it is not kINT7 or does not have EMCAL in readout, skipping"; - return; + return false; + } + if (fillHistos) { + mHistManager.fill(HIST("eventVertexZAll"), theCollision.posZ()); } - mHistManager.fill(HIST("eventVertexZAll"), theCollision.posZ()); if (vertexCut > 0 && std::abs(theCollision.posZ()) > vertexCut) { LOG(debug) << "Event not selected because of z-vertex cut z= " << theCollision.posZ() << " > " << vertexCut << " cm, skipping"; + return false; + } + if (fillHistos) { + mHistManager.fill(HIST("eventsSelected"), 1); + mHistManager.fill(HIST("eventVertexZSelected"), theCollision.posZ()); + } + return true; + } + + /// \brief Process EMCAL clusters that are matched to a collisions + void processCollisions(CollisionEvSelIt const& theCollision, o2::aod::BCs const&, SelectedClusters const& clusters, o2::aod::EMCALClusterCells const&, o2::aod::Calos const&, o2::aod::EMCALMatchedTracks const& matchedtracks, TracksPID const& alltracks) + { + if (!isEventSelected(theCollision, true)) { return; } - mHistManager.fill(HIST("eventsSelected"), 1); - mHistManager.fill(HIST("eventVertexZSelected"), theCollision.posZ()); + float lowP = 0.75f; + float midP = 1.25f; for (const auto& alltrack : alltracks) { - double trackEta, trackPhi; - trackEta = alltrack.trackEtaEmcal(); - trackPhi = alltrack.trackPhiEmcal(); + double trackEta = alltrack.trackEtaEmcal(); + double trackPhi = alltrack.trackPhiEmcal(); if (alltrack.isGlobalTrack()) { // NSigma of all global tracks without matching mHistManager.fill(HIST("TrackTM_NSigma"), alltrack.tpcNSigmaEl(), alltrack.pt()); mHistManager.fill(HIST("TrackEtaPhi"), trackEta, trackPhi); @@ -284,11 +418,11 @@ struct EmcTmMonitor { // match.track_as() with // using globTracks = o2::soa::Join; // In this example the counter t is just used to only look at the closest match - float dEta, dPhi, pT, abs_p, trackEta, trackPhi, nSigmaEl; - int supermoduleID; + float dEta{0}, dPhi{0}, pT{0}, absP{0}, trackEta{0}, trackPhi{0}, nSigmaEl{0}; + int supermoduleID = -1; try { supermoduleID = mGeometry->SuperModuleNumberFromEtaPhi(cluster.eta(), cluster.phi()); - } catch (o2::emcal::InvalidPositionException& e) { + } catch (o2::emcal::InvalidPositionException const& e) { // Imprecision of the position at the sector boundaries, mostly due to // vertex imprecision. Skip these clusters for the now. continue; @@ -304,7 +438,7 @@ struct EmcTmMonitor { nSigmaEl = match.track_as().tpcNSigmaEl(); // exmple of how to access any property of the matched tracks (tracks are sorted by how close they are to cluster) LOG(debug) << "Pt of match" << match.track_as().pt(); - abs_p = std::abs(match.track_as().p()); + absP = std::abs(match.track_as().p()); // only consider closest match trackEta = match.track_as().trackEtaEmcal(); trackPhi = match.track_as().trackPhiEmcal(); @@ -330,14 +464,14 @@ struct EmcTmMonitor { mHistManager.fill(HIST("clusterTM_dEtadPhi"), dEta, dPhi, supermoduleID); mHistManager.fill(HIST("clusterTM_dEtaPt"), dEta, pT, supermoduleID); } - double eOverP = cluster.energy() / abs_p; + double eOverP = cluster.energy() / absP; mHistManager.fill(HIST("clusterTM_dEtaTN"), dEta, t); mHistManager.fill(HIST("clusterTM_dPhiTN"), dPhi, t); mHistManager.fill(HIST("clusterTM_dRTN"), std::sqrt(dPhi * dPhi + dEta * dEta), t); mHistManager.fill(HIST("clusterTM_EoverP_E"), eOverP, cluster.energy(), t); mHistManager.fill(HIST("clusterTM_dEtadPhi"), dEta, dPhi, t); mHistManager.fill(HIST("clusterEMatched"), cluster.energy(), t); - mHistManager.fill(HIST("clusterTM_EvsP"), cluster.energy(), abs_p, t); + mHistManager.fill(HIST("clusterTM_EvsP"), cluster.energy(), absP, t); mHistManager.fill(HIST("clusterTM_EoverP_Pt"), eOverP, match.track_as().pt(), t); mHistManager.fill(HIST("clusterTM_NSigma"), nSigmaEl, match.track_as().pt(), t); mHistManager.fill(HIST("MatchedTrackEtaPhi"), trackEta, trackPhi); @@ -357,9 +491,9 @@ struct EmcTmMonitor { if (t == 0) { mHistManager.fill(HIST("clusterTM_PosdPhiPt"), dPhi, pT, supermoduleID); mHistManager.fill(HIST("clusterTM_PosdEtadPhi"), dEta, dPhi, supermoduleID); - if (abs_p < lowP) { + if (absP < lowP) { mHistManager.fill(HIST("clusterTM_PosdEtadPhi_Pl0_75"), dEta, dPhi, supermoduleID); - } else if (abs_p >= midP) { + } else if (absP >= midP) { mHistManager.fill(HIST("clusterTM_PosdEtadPhi_Pgeq1_25"), dEta, dPhi, supermoduleID); } else { mHistManager.fill(HIST("clusterTM_PosdEtadPhi_0_75leqPl1_25"), dEta, dPhi, supermoduleID); @@ -371,9 +505,9 @@ struct EmcTmMonitor { if (t == 0) { mHistManager.fill(HIST("clusterTM_NegdPhiPt"), dPhi, pT, supermoduleID); mHistManager.fill(HIST("clusterTM_NegdEtadPhi"), dEta, dPhi, supermoduleID); - if (abs_p < lowP) { + if (absP < lowP) { mHistManager.fill(HIST("clusterTM_NegdEtadPhi_Pl0_75"), dEta, dPhi, supermoduleID); - } else if (abs_p >= midP) { + } else if (absP >= midP) { mHistManager.fill(HIST("clusterTM_NegdEtadPhi_Pgeq1_25"), dEta, dPhi, supermoduleID); } else { mHistManager.fill(HIST("clusterTM_NegdEtadPhi_0_75leqPl1_25"), dEta, dPhi, supermoduleID); @@ -446,6 +580,191 @@ struct EmcTmMonitor { } PROCESS_SWITCH(EmcTmMonitor, processCollisions, "Process clusters from collision", true); + /// \brief Information of a reconstructed track needed for the MC truth of the matching + struct McTrackInfo { + int64_t id{-1}; + bool isGlobal{false}; + bool hasITS{false}; + bool hasTPC{false}; + float etaEmcal{-999.f}; + float phiEmcal{-999.f}; + float p{0.f}; + float pt{0.f}; + }; + + static constexpr float NotOnEmcal = -900.f; // trackEta/PhiEmcal of tracks that were not extrapolated to the EMCal are -999 + + /// \brief Truth class of an MC particle + template + TmTruthClass truthClassOf(TMcParticle const& mc, TMcParticles const& mcParticles) + { + const int pdg = mc.pdgCode(); + if (std::abs(pdg) == PDG_t::kElectron) { + if (mc.producedByGenerator()) { + return TmTruthClass::kPrimaryElectron; + } + if (mc.getProcess() == TMCProcess::kPPair && mc.has_mothers() && mcParticles.iteratorAt(mc.mothersIds()[0]).pdgCode() == PDG_t::kGamma) { + return TmTruthClass::kConversionElectron; + } + return TmTruthClass::kOtherElectron; + } + if (pdg == PDG_t::kGamma) { + return TmTruthClass::kPhoton; + } + const TParticlePDG* info = mPdgDatabase->GetParticle(pdg); + if (info == nullptr) { + return TmTruthClass::kOther; + } + return info->Charge() != 0. ? TmTruthClass::kOtherCharged : TmTruthClass::kNeutral; + } + + /// \brief Type of a reconstructed track + static TmTrackType trackTypeOf(McTrackInfo const& info) + { + if (info.isGlobal) { + return TmTrackType::kGlobal; + } + if (info.hasITS && info.hasTPC) { + return TmTrackType::kItsTpc; + } + if (info.hasTPC) { + return TmTrackType::kTpcOnly; + } + return info.hasITS ? TmTrackType::kItsOnly : TmTrackType::kOtherType; + } + + /// \brief MC truth of the track matching: follow the leading MC particle of each cluster through the matching chain. + /// Needs emcal-correction-task with processMCWithSecondaries and the photon conversion builder (aod::V0Legs). + void processMCTruth(CollisionEvSelIt const& theCollision, o2::aod::BCs const&, SelectedMcClusters const& clusters, o2::aod::EMCALMatchedTracks const& matchedTracks, o2::aod::EMCMatchSecs const& matchedSecondaries, TracksMC const& tracks, o2::aod::V0Legs const& v0Legs, o2::aod::McParticles const& mcParticles) + { + if (!isEventSelected(theCollision, !doprocessCollisions.value)) { + return; + } + + // reconstructed tracks of each MC particle in this collision, and the tracks used as V0 legs (secondary matching) + std::unordered_map> tracksOfMcParticle; + for (const auto& track : tracks) { + if (!track.has_mcParticle()) { + continue; + } + tracksOfMcParticle[track.mcParticleId()].push_back({track.globalIndex(), track.isGlobalTrack(), track.hasITS(), track.hasTPC(), track.trackEtaEmcal(), track.trackPhiEmcal(), track.p(), track.pt()}); + } + std::unordered_set v0LegTracks; + for (const auto& leg : v0Legs) { + v0LegTracks.insert(leg.trackId()); + } + + for (const auto& cluster : clusters) { + try { + mGeometry->SuperModuleNumberFromEtaPhi(cluster.eta(), cluster.phi()); + } catch (o2::emcal::InvalidPositionException const& e) { + continue; // same as in processCollisions + } + const float clusterPt = cluster.energy() / std::cosh(cluster.eta()); + if (m02highPt > 0. && cluster.m02() >= maxM02HighPt && clusterPt >= m02highPt) { + continue; + } + + // leading contributor: the labels are ordered by energy fraction (ClusterLabel::orderLabels in the ClusterFactory) + const auto mcIds = cluster.mcParticleIds(); + const int64_t leadId = mcIds.empty() ? -1 : mcIds[0]; + TmTruthClass truthClass = TmTruthClass::kOther; + float mcPt = 0.f; + float convRadius = -1.f; + if (leadId >= 0) { + const auto mc = mcParticles.iteratorAt(leadId); + mcPt = mc.pt(); + truthClass = truthClassOf(mc, mcParticles); + convRadius = std::hypot(mc.vx(), mc.vy()); + } + mHistManager.fill(HIST("MCTruth/hLeadingClass"), static_cast(truthClass), cluster.energy()); + const bool isConversion = truthClass == TmTruthClass::kConversionElectron; + if (truthClass != TmTruthClass::kConversionElectron && truthClass != TmTruthClass::kPrimaryElectron && truthClass != TmTruthClass::kOtherElectron && truthClass != TmTruthClass::kOtherCharged) { + continue; // neutral leading particle: no track to match + } + + // best reconstructed track of the leading particle: extrapolated and matching candidate first, then highest pT + const McTrackInfo* best = nullptr; + int bestScore = -1; + if (auto it = tracksOfMcParticle.find(leadId); it != tracksOfMcParticle.end()) { + for (const auto& info : it->second) { + const bool isExtrapolated = info.etaEmcal > NotOnEmcal && info.phiEmcal > NotOnEmcal; + const bool isCandidate = info.isGlobal || v0LegTracks.contains(info.id); + const int score = 2 * (isExtrapolated ? 1 : 0) + (isCandidate ? 1 : 0); + if (score > bestScore || (score == bestScore && info.pt > best->pt)) { + best = &info; + bestScore = score; + } + } + } + + TmTruthStatus status = TmTruthStatus::kNoRecoTrack; + TmTrackType type = TmTrackType::kNoTrack; + bool isV0Leg = false; + if (best != nullptr) { + type = trackTypeOf(*best); + isV0Leg = v0LegTracks.contains(best->id); + const bool isExtrapolated = best->etaEmcal > NotOnEmcal && best->phiEmcal > NotOnEmcal; + + // is the track among the primary or secondary matches of this cluster + bool isMatchedPrimary = false, isMatchedSecondary = false; + float storedDEta = 0.f, storedDPhi = 0.f; + for (const auto& match : matchedTracks.sliceBy(perClusterMatchedTracks, cluster.globalIndex())) { + if (match.trackId() == best->id) { + isMatchedPrimary = true; + storedDEta = match.deltaEta(); + storedDPhi = match.deltaPhi(); + break; + } + } + if (!isMatchedPrimary) { + for (const auto& match : matchedSecondaries.sliceBy(perClusterMatchedSecondaries, cluster.globalIndex())) { + if (match.trackId() == best->id) { + isMatchedSecondary = true; + storedDEta = match.deltaEta(); + storedDPhi = match.deltaPhi(); + break; + } + } + } + if (isMatchedPrimary) { + status = TmTruthStatus::kMatchedPrimary; + } else if (isMatchedSecondary) { + status = TmTruthStatus::kMatchedSecondary; + } else if (!isExtrapolated) { + status = TmTruthStatus::kNotExtrapolated; + } else if (!best->isGlobal && !isV0Leg) { + status = TmTruthStatus::kNoMatchingCandidate; + } else { + status = TmTruthStatus::kCandidateNotMatched; + } + + if (isExtrapolated) { + // the matching compares to the cluster position seen from the origin, the cluster table has eta seen from the vertex + const float clusterEtaOrigin = std::asinh(std::sinh(cluster.eta()) + theCollision.posZ() / emcalRadius); + const float dEta = best->etaEmcal - clusterEtaOrigin; + const float dPhi = std::remainder(best->phiEmcal - cluster.phi(), o2::constants::math::TwoPI); + mHistManager.fill(HIST("MCTruth/hResidualEta"), dEta, mcPt, static_cast(truthClass), static_cast(type)); + mHistManager.fill(HIST("MCTruth/hResidualPhi"), dPhi, mcPt, static_cast(truthClass), static_cast(type)); + if (isConversion) { + mHistManager.fill(HIST("MCTruth/hResidualConversion"), dEta, dPhi, mcPt); + } + if (isMatchedPrimary || isMatchedSecondary) { + mHistManager.fill(HIST("MCTruth/hResidualVsStored"), dEta - storedDEta, dPhi - storedDPhi); + } + } + if (isConversion && best->p > 0.f) { + mHistManager.fill(HIST("MCTruth/hEoverPConversion"), cluster.energy() / best->p, best->pt); + } + } + mHistManager.fill(HIST("MCTruth/hStatus"), static_cast(status), static_cast(truthClass), static_cast(type), static_cast(isV0Leg), mcPt); + if (isConversion) { + mHistManager.fill(HIST("MCTruth/hConvRadius"), static_cast(status), convRadius, mcPt); + } + } + } + PROCESS_SWITCH(EmcTmMonitor, processMCTruth, "MC truth of the track matching for clusters from charged particles", false); + /// \brief Create binning for cluster energy axis (variable bin size) /// \return vector with bin limits std::vector makeEnergyBinning() const @@ -493,7 +812,7 @@ struct EmcTmMonitor { result.emplace_back(50. + 2.0 * (i - 190)); } else if (i < 235) { // o2-linter: disable=magic-number (just numbers for binning) result.emplace_back(100. + 5.0 * (i - 215)); - } else if (i < 245) { // o2-linter: disable=magic-number (just numbers for binning) + } else { result.emplace_back(200. + 10.0 * (i - 235)); } } @@ -529,9 +848,9 @@ struct EmcTmMonitor { } }; -WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +WorkflowSpec defineDataProcessing(ConfigContext const& context) { WorkflowSpec workflow{ - adaptAnalysisTask(cfgc)}; + adaptAnalysisTask(context)}; return workflow; }