diff --git a/PWGEM/PhotonMeson/TableProducer/skimmerPrimaryElectronFromDalitzEE.h b/PWGEM/PhotonMeson/TableProducer/skimmerPrimaryElectronFromDalitzEE.h index 559da2e0d8d..148f1667dd5 100644 --- a/PWGEM/PhotonMeson/TableProducer/skimmerPrimaryElectronFromDalitzEE.h +++ b/PWGEM/PhotonMeson/TableProducer/skimmerPrimaryElectronFromDalitzEE.h @@ -48,6 +48,8 @@ #include // IWYU pragma: keep (do not replace with Math/Vector4Dfwd.h) #include +#include +#include #include #include @@ -121,6 +123,33 @@ struct skimmerPrimaryElectronFromDalitzEE { EpEmPairsAndPhoton = 3 }; + enum RejectReason { + kAllTracks = 0, + kPassed, + kNoMcParticle, + kNotTrueElectron, + kTOFRequired, + kNoITS, + kITSChi2, + kITSNCls, + kITSNClsInnerBarrel, + kNoTPC, + kTPCChi2, + kTPCNClsFound, + kTPCCrossedRows, + kTPCCrossedRowsOverFindable, + kTPCSharedFraction, + kPropagation, + kDCAxyz, + kDCA3D, + kPtEta, + kITSClusterSize, + kITSsaPt, + kNotElectron, + kNRejectReasons + }; + static constexpr std::array RejectReasonLabels = {"all tracks", "accepted tracks", "no MC particle", "not true electron", "TOF required", "no ITS", "ITS chi2", "ITS Ncls", "ITS Ncls IB", "no TPC", "TPC chi2", "TPC Ncls found", "TPC crossed rows", "TPC Ncr/Nf", "TPC shared frac", "propagation", "DCAxy/z", "DCA 3D", "pt/eta", "ITS cluster size", "ITSsa pt", "not electron"}; + // Configurables Configurable ccdburl{"ccdb-url", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; Configurable grpPath{"grpPath", "GLO/GRP/GRP", "Path of the grp file"}; @@ -172,9 +201,11 @@ struct skimmerPrimaryElectronFromDalitzEE { Configurable intercept{"intercept", -0.0380, "intercept for m vs. phiv"}; Configurable useTOFNSigmaDeltaBC{"useTOFNSigmaDeltaBC", false, "Flag to shift delta BC for TOF n sigma (only with TTCA)"}; Configurable storeOnlyTrueElectronMC{"storeOnlyTrueElectronMC", false, "Flag to store only true electron in MC"}; + Configurable fillHistogramsBeforeSelection{"fillHistogramsBeforeSelection", false, "Flag to fill track histograms before track selections (only for fillingMode 1)"}; o2::framework::HistogramRegistry fRegistry{"output", {}, o2::framework::OutputObjHandlingPolicy::AnalysisObject, false, false}; static constexpr std::array DileptonSigns = {"uls/", "lspp/", "lsmm/"}; + static constexpr std::array TrackDirs = {"Track/", "TrackTrueElec/", "TrackNoSel/", "TrackNoSelTrueElec/"}; int mRunNumber = 0; float dBz = 0.; @@ -237,6 +268,24 @@ struct skimmerPrimaryElectronFromDalitzEE { fRegistry.add("Track/hTOFNsigmaEl", "TOF n sigma el;p_{pv} (GeV/c);n #sigma_{e}^{TOF}", kTH2F, {{1000, 0, 10}, {100, -5, +5}}, false); fRegistry.add("Track/hTOFNsigmaPi", "TOF n sigma pi;p_{pv} (GeV/c);n #sigma_{#pi}^{TOF}", kTH2F, {{1000, 0, 10}, {100, -5, +5}}, false); + fRegistry.add("Track/hTrackRejectReason", "track rejection reason vs. p_{T};;p_{T} (GeV/c)", kTH2F, {{static_cast(kNRejectReasons), -0.5, static_cast(kNRejectReasons) - 0.5}, {200, 0, 20}}, false); + auto hRejectReason = fRegistry.get(HIST("Track/hTrackRejectReason")); + for (int i = 0; i < kNRejectReasons; i++) { + hRejectReason->GetXaxis()->SetBinLabel(i + 1, RejectReasonLabels[i]); + } + + if (doprocessMC) { + LOG(info) << "adding TrackTrueElec folder"; + fRegistry.addClone("Track/", "TrackTrueElec/"); + } + + if (fillHistogramsBeforeSelection) { + fRegistry.addClone("Track/", "TrackNoSel/"); + if (doprocessMC) { + fRegistry.addClone("Track/", "TrackNoSelTrueElec/"); + } + } + // pair fRegistry.add("Pair/uls/hTrackMvsPt", "m_{ee} vs. p_{T,ee};m_{ee} (GeV/c^{2});p_{T,ee} (GeV/c)", kTH2F, {{100, 0, 0.1}, {200, 0, 2}}, false); fRegistry.add("Pair/uls/hCheckEMvsPt", "m_{ee} vs. p_{T,ee};m_{ee} (GeV/c^{2});p_{T,ee} (GeV/c)", kTH2F, {{100, 0, 0.1}, {200, 0, 2}}, false); @@ -354,63 +403,72 @@ struct skimmerPrimaryElectronFromDalitzEE { } template - bool checkTrack(TCollision const& collision, TTrack const& track) + bool checkTrack(TCollision const& collision, TTrack const& track, int* rejectReason = nullptr) { + // Lambda function: Sets the reject reason to the correct value and returns false. + // Used for every return statement in the checkTrack function where false is returned + auto reject = [rejectReason](RejectReason reason) { + if (rejectReason != nullptr) { + *rejectReason = reason; + } + return false; + }; + if constexpr (isMC) { if (!track.has_mcParticle()) { - return false; + return reject(kNoMcParticle); } if (storeOnlyTrueElectronMC) { const auto& mcParticle = track.template mcParticle_as(); - if (std::abs(mcParticle.pdgCode()) != 11) { - return false; + if (std::abs(mcParticle.pdgCode()) != PDG_t::kElectron) { + return reject(kNotTrueElectron); } } } float tofNSigmaEl = mapTOFNsigmaReassociated[std::make_pair(collision.globalIndex(), track.globalIndex())]; if (requireTOF && !(track.hasTOF() && std::fabs(tofNSigmaEl) < maxTOFNsigmaEl)) { - return false; + return reject(kTOFRequired); } if (!track.hasITS()) { - return false; + return reject(kNoITS); } if (track.itsChi2NCl() < minchi2its || maxchi2its < track.itsChi2NCl()) { - return false; + return reject(kITSChi2); } if (track.itsNCls() < min_ncluster_its) { - return false; + return reject(kITSNCls); } if (track.itsNClsInnerBarrel() < min_ncluster_itsib) { - return false; + return reject(kITSNClsInnerBarrel); } if (!includeITSsa && !track.hasTPC()) { - return false; + return reject(kNoTPC); } if (track.hasTPC()) { if (track.tpcChi2NCl() < minchi2tpc || maxchi2tpc < track.tpcChi2NCl()) { - return false; + return reject(kTPCChi2); } if (track.tpcNClsFound() < min_ncluster_tpc) { - return false; + return reject(kTPCNClsFound); } if (track.tpcNClsCrossedRows() < mincrossedrows) { - return false; + return reject(kTPCCrossedRows); } if (track.tpcCrossedRowsOverFindableCls() < min_tpc_cr_findable_ratio) { - return false; + return reject(kTPCCrossedRowsOverFindable); } if (track.tpcFractionSharedCls() > max_frac_shared_clusters_tpc) { - return false; + return reject(kTPCSharedFraction); } } @@ -422,13 +480,13 @@ struct skimmerPrimaryElectronFromDalitzEE { mVtx.setCov(collision.covXX(), collision.covXY(), collision.covYY(), collision.covXZ(), collision.covYZ(), collision.covZZ()); bool isPropOK = o2::base::Propagator::Instance()->propagateToDCABxByBz(mVtx, trackParCov, 2.f, matCorr, &mDcaInfoCov); if (!isPropOK) { - return false; + return reject(kPropagation); } float dcaXY = mDcaInfoCov.getY(); float dcaZ = mDcaInfoCov.getZ(); if (std::fabs(dcaXY) > dca_xy_max || std::fabs(dcaZ) > dca_z_max) { - return false; + return reject(kDCAxyz); } float dca3D = 999.f; @@ -440,11 +498,11 @@ struct skimmerPrimaryElectronFromDalitzEE { dca3D = std::sqrt(std::fabs(chi2) / 2.); } if (dca3D > dca_3d_sigma_max) { - return false; + return reject(kDCA3D); } if (trackParCov.getPt() < minpt || std::fabs(trackParCov.getEta()) > maxeta) { - return false; + return reject(kPtEta); } int total_cluster_size = 0, nl = 0; @@ -457,11 +515,11 @@ struct skimmerPrimaryElectronFromDalitzEE { } if (maxMeanITSClusterSize < static_cast(total_cluster_size) / static_cast(nl) * std::cos(std::atan(trackParCov.getTgl()))) { - return false; + return reject(kITSClusterSize); } if ((track.hasITS() && !track.hasTPC() && !track.hasTRD() && !track.hasTOF()) && maxpt_itssa < track.pt()) { - return false; + return reject(kITSsaPt); } // if ((track.hasITS() && !track.hasTPC() && !track.hasTRD() && !track.hasTOF()) && maxpt_itssa < track.pt()) { @@ -528,11 +586,24 @@ struct skimmerPrimaryElectronFromDalitzEE { { for (const auto& track : tracks) { - if (!checkTrack(collision, track)) { - continue; + + if (fillHistogramsBeforeSelection) { + fillTrackHistograms(track); } - if (!isElectron(collision, track)) { + fRegistry.fill(HIST("Track/hTrackRejectReason"), kAllTracks, track.pt()); + int rejectReason = kPassed; + if (checkTrack(collision, track, &rejectReason) && !isElectron(collision, track)) { + rejectReason = kNotElectron; + } + fRegistry.fill(HIST("Track/hTrackRejectReason"), rejectReason, track.pt()); + if constexpr (isMC) { + if (isTruePrimaryElectron(track)) { + fRegistry.fill(HIST("TrackTrueElec/hTrackRejectReason"), kAllTracks, track.pt()); + fRegistry.fill(HIST("TrackTrueElec/hTrackRejectReason"), rejectReason, track.pt()); + } + } + if (rejectReason != kPassed) { continue; } @@ -544,7 +615,16 @@ struct skimmerPrimaryElectronFromDalitzEE { template void fillPairInfo(TCollision const& collision, TTracks1 const& tracks1, TTracks2 const& tracks2) { - if constexpr (pairtype == 0) { // ULS + if constexpr (pairtype == 0) { // ULS + if (fillHistogramsBeforeSelection) { // we cannot fill the track histograms in the pair loop as they would get filled multiple times. Hence we loop over all tracks here + for (const auto& track : tracks1) { + fillTrackHistograms(track); + } + for (const auto& track : tracks2) { + fillTrackHistograms(track); + } + } + for (const auto& [t1, t2] : combinations(o2::soa::CombinationsFullIndexPolicy(tracks1, tracks2))) { if (!checkTrack(collision, t1) || !checkTrack(collision, t2)) { @@ -621,44 +701,66 @@ struct skimmerPrimaryElectronFromDalitzEE { } } - template + template + bool isTruePrimaryElectron(TTrack const& track) + { + if (!track.has_mcParticle()) { + return false; + } + const auto mcParticle = track.template mcParticle_as(); + return std::abs(mcParticle.pdgCode()) == PDG_t::kElectron && mcParticle.isPhysicalPrimary(); // o2-linter: disable=pdg/explicit-code (electron) + } + + template void fillTrackHistograms(TTrack const& track) + { + constexpr int kDirAllTracks = beforeSelection ? 2 : 0; + fillTrackHistogramsImpl(track); + if constexpr (isMC) { + if (isTruePrimaryElectron(track)) { + fillTrackHistogramsImpl(track); + } + } + } + + template + void fillTrackHistogramsImpl(TTrack const& track) { float mcTunedTPCSignal = 0.f; if constexpr (isMC) { mcTunedTPCSignal = track.mcTunedTPCSignal(); } - fRegistry.fill(HIST("Track/hPt"), track.pt()); - fRegistry.fill(HIST("Track/hEtaPhi"), track.phi(), track.eta()); - fRegistry.fill(HIST("Track/hQoverPt"), track.sign() / track.pt()); - fRegistry.fill(HIST("Track/hRelDeltaPt"), track.pt(), track.sigma1Pt() * track.pt()); - fRegistry.fill(HIST("Track/hDCAxyz"), track.dcaXY(), track.dcaZ()); - fRegistry.fill(HIST("Track/hDCAxy_Pt"), track.pt(), track.dcaXY()); - fRegistry.fill(HIST("Track/hDCAz_Pt"), track.pt(), track.dcaZ()); - fRegistry.fill(HIST("Track/hDCAxyzSigma"), track.dcaXY() / std::sqrt(track.cYY()), track.dcaZ() / std::sqrt(track.cZZ())); - fRegistry.fill(HIST("Track/hDCAxyRes_Pt"), track.pt(), std::sqrt(track.cYY()) * 1e+4); // convert cm to um - fRegistry.fill(HIST("Track/hDCAzRes_Pt"), track.pt(), std::sqrt(track.cZZ()) * 1e+4); // convert cm to um - - fRegistry.fill(HIST("Track/hNclsTPC"), track.tpcNClsFound()); - fRegistry.fill(HIST("Track/hNcrTPC"), track.tpcNClsCrossedRows()); - fRegistry.fill(HIST("Track/hChi2TPC"), track.tpcChi2NCl()); - fRegistry.fill(HIST("Track/hTPCNcr2Nf"), track.tpcCrossedRowsOverFindableCls()); - fRegistry.fill(HIST("Track/hTPCNcls2Nf"), track.tpcFoundOverFindableCls()); - fRegistry.fill(HIST("Track/hTPCNclsShared"), track.pt(), track.tpcFractionSharedCls()); - fRegistry.fill(HIST("Track/hTPCdEdx"), track.tpcInnerParam(), track.tpcSignal()); - fRegistry.fill(HIST("Track/hTPCdEdxMC"), track.tpcInnerParam(), mcTunedTPCSignal); - fRegistry.fill(HIST("Track/hTPCNsigmaEl"), track.tpcInnerParam(), track.tpcNSigmaEl()); - fRegistry.fill(HIST("Track/hTPCNsigmaPi"), track.tpcInnerParam(), track.tpcNSigmaPi()); - - fRegistry.fill(HIST("Track/hChi2TOF"), track.tofChi2()); - fRegistry.fill(HIST("Track/hTOFbeta"), track.p(), track.beta()); - fRegistry.fill(HIST("Track/hTOFNsigmaEl"), track.p(), track.tofNSigmaEl()); - fRegistry.fill(HIST("Track/hTOFNsigmaPi"), track.p(), track.tofNSigmaPi()); - - fRegistry.fill(HIST("Track/hNclsITS"), track.itsNCls()); - fRegistry.fill(HIST("Track/hChi2ITS"), track.itsChi2NCl()); - fRegistry.fill(HIST("Track/hITSClusterMap"), track.itsClusterMap()); + fRegistry.fill(HIST(TrackDirs[trackType]) + HIST("hPt"), track.pt()); + fRegistry.fill(HIST(TrackDirs[trackType]) + HIST("hEtaPhi"), track.phi(), track.eta()); + fRegistry.fill(HIST(TrackDirs[trackType]) + HIST("hQoverPt"), track.sign() / track.pt()); + fRegistry.fill(HIST(TrackDirs[trackType]) + HIST("hRelDeltaPt"), track.pt(), track.sigma1Pt() * track.pt()); + fRegistry.fill(HIST(TrackDirs[trackType]) + HIST("hDCAxyz"), track.dcaXY(), track.dcaZ()); + fRegistry.fill(HIST(TrackDirs[trackType]) + HIST("hDCAxy_Pt"), track.pt(), track.dcaXY()); + fRegistry.fill(HIST(TrackDirs[trackType]) + HIST("hDCAz_Pt"), track.pt(), track.dcaZ()); + fRegistry.fill(HIST(TrackDirs[trackType]) + HIST("hDCAxyzSigma"), track.dcaXY() / std::sqrt(track.cYY()), track.dcaZ() / std::sqrt(track.cZZ())); + fRegistry.fill(HIST(TrackDirs[trackType]) + HIST("hDCAxyRes_Pt"), track.pt(), std::sqrt(track.cYY()) * 1e+4); // convert cm to um + fRegistry.fill(HIST(TrackDirs[trackType]) + HIST("hDCAzRes_Pt"), track.pt(), std::sqrt(track.cZZ()) * 1e+4); // convert cm to um + + fRegistry.fill(HIST(TrackDirs[trackType]) + HIST("hNclsTPC"), track.tpcNClsFound()); + fRegistry.fill(HIST(TrackDirs[trackType]) + HIST("hNcrTPC"), track.tpcNClsCrossedRows()); + fRegistry.fill(HIST(TrackDirs[trackType]) + HIST("hChi2TPC"), track.tpcChi2NCl()); + fRegistry.fill(HIST(TrackDirs[trackType]) + HIST("hTPCNcr2Nf"), track.tpcCrossedRowsOverFindableCls()); + fRegistry.fill(HIST(TrackDirs[trackType]) + HIST("hTPCNcls2Nf"), track.tpcFoundOverFindableCls()); + fRegistry.fill(HIST(TrackDirs[trackType]) + HIST("hTPCNclsShared"), track.pt(), track.tpcFractionSharedCls()); + fRegistry.fill(HIST(TrackDirs[trackType]) + HIST("hTPCdEdx"), track.tpcInnerParam(), track.tpcSignal()); + fRegistry.fill(HIST(TrackDirs[trackType]) + HIST("hTPCdEdxMC"), track.tpcInnerParam(), mcTunedTPCSignal); + fRegistry.fill(HIST(TrackDirs[trackType]) + HIST("hTPCNsigmaEl"), track.tpcInnerParam(), track.tpcNSigmaEl()); + fRegistry.fill(HIST(TrackDirs[trackType]) + HIST("hTPCNsigmaPi"), track.tpcInnerParam(), track.tpcNSigmaPi()); + + fRegistry.fill(HIST(TrackDirs[trackType]) + HIST("hChi2TOF"), track.tofChi2()); + fRegistry.fill(HIST(TrackDirs[trackType]) + HIST("hTOFbeta"), track.p(), track.beta()); + fRegistry.fill(HIST(TrackDirs[trackType]) + HIST("hTOFNsigmaEl"), track.p(), track.tofNSigmaEl()); + fRegistry.fill(HIST(TrackDirs[trackType]) + HIST("hTOFNsigmaPi"), track.p(), track.tofNSigmaPi()); + + fRegistry.fill(HIST(TrackDirs[trackType]) + HIST("hNclsITS"), track.itsNCls()); + fRegistry.fill(HIST(TrackDirs[trackType]) + HIST("hChi2ITS"), track.itsChi2NCl()); + fRegistry.fill(HIST(TrackDirs[trackType]) + HIST("hITSClusterMap"), track.itsClusterMap()); int totalClusterSize = 0, nl = 0; for (unsigned int layer = 0; layer < 7; layer++) { // o2-linter: disable=magic-number (number of ITS layers) @@ -686,9 +788,9 @@ struct skimmerPrimaryElectronFromDalitzEE { } totalClusterSizeOb += clusterSizePerLayer; } - fRegistry.fill(HIST("Track/hMeanClusterSizeITS"), track.p(), static_cast(totalClusterSize) / static_cast(nl) * std::cos(std::atan(track.tgl()))); - fRegistry.fill(HIST("Track/hMeanClusterSizeITSib"), track.p(), static_cast(totalClusterSizeIb) / static_cast(nl_ib) * std::cos(std::atan(track.tgl()))); - fRegistry.fill(HIST("Track/hMeanClusterSizeITSob"), track.p(), static_cast(totalClusterSizeOb) / static_cast(nlOb) * std::cos(std::atan(track.tgl()))); + fRegistry.fill(HIST(TrackDirs[trackType]) + HIST("hMeanClusterSizeITS"), track.p(), static_cast(totalClusterSize) / static_cast(nl) * std::cos(std::atan(track.tgl()))); + fRegistry.fill(HIST(TrackDirs[trackType]) + HIST("hMeanClusterSizeITSib"), track.p(), static_cast(totalClusterSizeIb) / static_cast(nl_ib) * std::cos(std::atan(track.tgl()))); + fRegistry.fill(HIST(TrackDirs[trackType]) + HIST("hMeanClusterSizeITSob"), track.p(), static_cast(totalClusterSizeOb) / static_cast(nlOb) * std::cos(std::atan(track.tgl()))); } template