diff --git a/PWGEM/PhotonMeson/TableProducer/skimmerPrimaryElectronFromDalitzEE.cxx b/PWGEM/PhotonMeson/TableProducer/skimmerPrimaryElectronFromDalitzEE.cxx index 9c0c52b24de..cdedb5d5f4a 100644 --- a/PWGEM/PhotonMeson/TableProducer/skimmerPrimaryElectronFromDalitzEE.cxx +++ b/PWGEM/PhotonMeson/TableProducer/skimmerPrimaryElectronFromDalitzEE.cxx @@ -16,6 +16,7 @@ #include "PWGEM/Dilepton/Utils/PairUtilities.h" #include "PWGEM/PhotonMeson/DataModel/gammaTables.h" +#include "Common/DataModel/CollisionAssociationTables.h" #include "Common/DataModel/EventSelection.h" #include "Common/DataModel/PIDResponseTOF.h" #include "Common/DataModel/PIDResponseTPC.h" @@ -45,8 +46,10 @@ #include #include #include +#include #include #include +#include #include #include @@ -60,17 +63,40 @@ using MyCollisions = soa::Join; using MyCollisionsWithSWT = soa::Join; using MyCollisionsMC = soa::Join; -using MyTracks = soa::Join; +using MyTracks = soa::Join; using MyTrack = MyTracks::iterator; using MyTracksMC = soa::Join; using MyTrackMC = MyTracksMC::iterator; -struct SkimmerPrimaryElectronFromDalitzEE { +namespace o2::aod +{ +namespace pwgem::pm::recalculatedtofpid +{ +DECLARE_SOA_COLUMN(BetaRecalculated, betaRecalculated, float); +DECLARE_SOA_COLUMN(TOFNSigmaElRecalculated, tofNSigmaElRecalculated, float); +} // namespace pwgem::pm::recalculatedtofpid + +DECLARE_SOA_TABLE(EMTOFNSigmas, "AOD", "EMTOFNSIGMA", // make std::map in your tasks later. // Don't store this table in the derived data. + o2::aod::emprimaryelectron::CollisionId, o2::aod::emprimaryelectron::TrackId, + o2::aod::pwgem::pm::recalculatedtofpid::BetaRecalculated, o2::aod::pwgem::pm::recalculatedtofpid::TOFNSigmaElRecalculated); + +using EMTOFNSigma = EMTOFNSigmas::iterator; +} // namespace o2::aod + +struct skimmerPrimaryElectronFromDalitzEE { SliceCache cache; Preslice perCol = o2::aod::track::collisionId; - Preslice perColPcm = o2::aod::v0photonkf::collisionId; + PresliceOptional perTracksCollision = aod::track::collisionId; + Preslice perCol_pcm = o2::aod::v0photonkf::collisionId; + Preslice trackIndicesPerCollision = aod::track_association::collisionId; Produces emprimaryelectrons; Produces emprimaryelectronsDeDxMC; + Service mTOFResponse; + + Produces emtofs; // Configurables Configurable ccdburl{"ccdb-url", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; @@ -82,11 +108,13 @@ struct SkimmerPrimaryElectronFromDalitzEE { Configurable dBzInput{"dBzInput", -999, "bz field in kG, -999 is automatic"}; Configurable min_ncluster_tpc{"min_ncluster_tpc", 0, "min ncluster tpc"}; // o2-linter: disable=name/function-variable (renaming configs would mess up hyperloop) Configurable mincrossedrows{"mincrossedrows", 70, "min. crossed rows"}; - Configurable min_tpc_cr_findable_ratio{"min_tpc_cr_findable_ratio", 0.8, "min. TPC Ncr/Nf ratio"}; // o2-linter: disable=name/function-variable (renaming configs would mess up hyperloop) - Configurable max_frac_shared_clusters_tpc{"max_frac_shared_clusters_tpc", 999.f, "max fraction of shared clusters in TPC"}; // o2-linter: disable=name/function-variable (renaming configs would mess up hyperloop) - Configurable min_ncluster_its{"min_ncluster_its", 4, "min ncluster its"}; // o2-linter: disable=name/function-variable (renaming configs would mess up hyperloop) - Configurable min_ncluster_itsib{"min_ncluster_itsib", 1, "min ncluster itsib"}; // o2-linter: disable=name/function-variable (renaming configs would mess up hyperloop) + Configurable min_tpc_cr_findable_ratio{"min_tpc_cr_findable_ratio", 0.8, "min. TPC Ncr/Nf ratio"}; + Configurable max_frac_shared_clusters_tpc{"max_frac_shared_clusters_tpc", 999.f, "max fraction of shared clusters in TPC"}; + Configurable min_ncluster_its{"min_ncluster_its", 4, "min ncluster its"}; + Configurable min_ncluster_itsib{"min_ncluster_itsib", 1, "min ncluster itsib"}; + Configurable minchi2tpc{"minchi2tpc", 0.0, "min. chi2/NclsTPC"}; Configurable maxchi2tpc{"maxchi2tpc", 5.0, "max. chi2/NclsTPC"}; + Configurable minchi2its{"minchi2its", 0.0, "min. chi2/NclsITS"}; Configurable maxchi2its{"maxchi2its", 36.0, "max. chi2/NclsITS"}; Configurable minpt{"minpt", 0.05, "min pt for ITS-TPC track"}; Configurable maxeta{"maxeta", 2.0, "max eta acceptance"}; @@ -99,13 +127,23 @@ struct SkimmerPrimaryElectronFromDalitzEE { Configurable minTPCNsigmaPi{"minTPCNsigmaPi", 0.0, "min. TPC n sigma for pion exclusion"}; Configurable minTOFNsigmaEl{"minTOFNsigmaEl", -3.5, "min. TOF n sigma for electron inclusion"}; Configurable maxTOFNsigmaEl{"maxTOFNsigmaEl", +3.5, "max. TOF n sigma for electron inclusion"}; + Configurable minTPCNsigmaKa{"minTPCNsigmaKa", -2.5, "min. TPC n sigma for kaon exclusion"}; + Configurable maxTPCNsigmaKa{"maxTPCNsigmaKa", 2.5, "max. TPC n sigma for kaon exclusion"}; + Configurable minTPCNsigmaPr{"minTPCNsigmaPr", -2.5, "min. TPC n sigma for proton exclusion"}; + Configurable maxTPCNsigmaPr{"maxTPCNsigmaPr", 2.5, "max. TPC n sigma for proton exclusion"}; + Configurable requireTOF{"requireTOF", false, "require TOF hit"}; + Configurable min_pin_for_pion_rejection{"min_pin_for_pion_rejection", 0.0, "pion rejection is applied above this pin"}; // this is used only in TOFreq + Configurable max_pin_for_pion_rejection{"max_pin_for_pion_rejection", 0.5, "pion rejection is applied below this pin"}; Configurable maxMee{"maxMee", 0.04, "max. mee to store dalitz ee pairs"}; Configurable fillLS{"fillLS", true, "flag to fill LS histograms for QA"}; + Configurable fillWithPairs{"fillWithPairs", false, "flag to fill table based on pair information"}; Configurable includeITSsa{"includeITSsa", false, "Flag to include ITSsa tracks"}; Configurable maxpt_itssa{"maxpt_itssa", 0.15, "max pt for ITSsa track"}; // o2-linter: disable=name/function-variable (renaming configs would mess up hyperloop) Configurable maxMeanITSClusterSize{"maxMeanITSClusterSize", 16, "max x cos(lambda)"}; Configurable slope{"slope", 0.0185, "slope for m vs. phiv"}; 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"}; HistogramRegistry fRegistry{"output", {}, OutputObjHandlingPolicy::AnalysisObject, false, false}; static constexpr std::array DileptonSigns = {"uls/", "lspp/", "lsmm/"}; @@ -114,8 +152,9 @@ struct SkimmerPrimaryElectronFromDalitzEE { float dBz = 0.; Service ccdb{}; o2::base::Propagator::MatCorrType matCorr = o2::base::Propagator::MatCorrType::USEMatCorrNONE; + o2::dataformats::VertexBase mVtx; - void init(InitContext&) + void init(InitContext& initContext) { // if (doprocessRec && doprocessRec_SWT) { // LOGF(fatal, "Cannot enable doprocessRec and doprocessRec_SWT at the same time. Please choose one."); @@ -129,6 +168,10 @@ struct SkimmerPrimaryElectronFromDalitzEE { ccdb->setLocalObjectValidityChecking(); ccdb->setFatalWhenNull(false); + LOGF(info, "before TOF initSetup"); + mTOFResponse->initSetup(ccdb, initContext); + LOGF(info, "after TOF initSetup"); + fRegistry.add("Track/hPt", "pT;p_{T} (GeV/c)", kTH1F, {{1000, 0.0f, 10}}, false); fRegistry.add("Track/hEtaPhi", "#eta vs. #varphi;#varphi (rad.);#eta", kTH2F, {{180, 0, o2::constants::math::TwoPI}, {400, -2.0f, 2.0f}}, false); fRegistry.add("Track/hQoverPt", "q/pT;q/p_{T} (GeV/c)^{-1}", kTH1F, {{400, -20, 20}}, false); @@ -167,8 +210,18 @@ struct SkimmerPrimaryElectronFromDalitzEE { fRegistry.add("Track/hTOFNsigmaPi", "TOF n sigma pi;p_{pv} (GeV/c);n #sigma_{#pi}^{TOF}", kTH2F, {{1000, 0, 10}, {100, -5, +5}}, false); // 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); fRegistry.add("Pair/uls/hMvsPt", "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/hMvsPhiV", "m_{ee} vs. #varphi_{V};#varphi_{V} (rad.);m_{ee} (GeV/c^{2})", kTH2F, {{180, 0, o2::constants::math::PI}, {100, 0, 0.1}}, false); + fRegistry.add("Pair/uls/hMCutMvsPt", "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/hMPhiCutMvsPt", "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/hTrackMvsPhiV", "m_{ee} vs. #varphi_{V};#varphi_{V} (rad.);m_{ee} (GeV/c^{2})", kTH2F, {{180, 0, M_PI}, {100, 0, 0.1}}, false); + fRegistry.add("Pair/uls/hCheckEMvsPhiV", "m_{ee} vs. #varphi_{V};#varphi_{V} (rad.);m_{ee} (GeV/c^{2})", kTH2F, {{180, 0, M_PI}, {100, 0, 0.1}}, false); + fRegistry.add("Pair/uls/hMvsPhiV", "m_{ee} vs. #varphi_{V};#varphi_{V} (rad.);m_{ee} (GeV/c^{2})", kTH2F, {{180, 0, M_PI}, {100, 0, 0.1}}, false); + fRegistry.add("Pair/uls/hMCutMvsPhiV", "m_{ee} vs. #varphi_{V};#varphi_{V} (rad.);m_{ee} (GeV/c^{2})", kTH2F, {{180, 0, M_PI}, {100, 0, 0.1}}, false); + fRegistry.add("Pair/uls/hMPhiCutMvsPhiV", "m_{ee} vs. #varphi_{V};#varphi_{V} (rad.);m_{ee} (GeV/c^{2})", kTH2F, {{180, 0, M_PI}, {100, 0, 0.1}}, false); + fRegistry.addClone("Pair/uls/", "Pair/lspp/"); fRegistry.addClone("Pair/uls/", "Pair/lsmm/"); } @@ -215,20 +268,105 @@ struct SkimmerPrimaryElectronFromDalitzEE { mRunNumber = bc.runNumber(); } - template - bool checkTrack(TCollision const&, TTrack const& track) + template + void calculateTOFNSigmaWithReassociation(TCollisions const& collisions, TBCs const&, TTracks const& tracks, TTrackAssoc const& trackIndices) + { + if (useTOFNSigmaDeltaBC) { + if constexpr (withTTCA) { + for (const auto& collision : collisions) { + if (mapCollisionTime.find(collision.globalIndex()) == mapCollisionTime.end()) { + continue; + } + auto bcCollision = collision.template bc_as(); + auto trackIdsThisCollision = trackIndices.sliceBy(trackIndicesPerCollision, collision.globalIndex()); + for (const auto& trackId : trackIdsThisCollision) { + auto track = trackId.template track_as(); + if (!track.hasITS() || !track.hasTPC()) { // apply only minimal cut + continue; + } + + if (track.hasTOF() && track.has_collision()) { // TTCA may use orphan tracks. + auto bcTrack = track.template collision_as().template bc_as(); + float tofNSigmaEl = mTOFResponse->nSigma(track.tofSignalInAnotherBC(bcTrack.globalBC(), bcCollision.globalBC()), track.tofExpMom(), track.length(), track.p(), track.eta(), mapCollisionTime[collision.globalIndex()], mapCollisionTimeError[collision.globalIndex()]); + float beta = track.length() / (track.tofSignalInAnotherBC(bcTrack.globalBC(), bcCollision.globalBC()) - mapCollisionTime[collision.globalIndex()]) / (TMath::C() * 1e+2 * 1e-12); + mapTOFNsigmaReassociated[std::make_pair(collision.globalIndex(), track.globalIndex())] = tofNSigmaEl; + mapTOFBetaReassociated[std::make_pair(collision.globalIndex(), track.globalIndex())] = beta; + emtofs(collision.globalIndex(), track.globalIndex(), mapTOFBetaReassociated[std::make_pair(collision.globalIndex(), track.globalIndex())], mapTOFNsigmaReassociated[std::make_pair(collision.globalIndex(), track.globalIndex())]); + } else { + mapTOFNsigmaReassociated[std::make_pair(collision.globalIndex(), track.globalIndex())] = track.tofNSigmaEl(); + mapTOFBetaReassociated[std::make_pair(collision.globalIndex(), track.globalIndex())] = track.beta(); + emtofs(collision.globalIndex(), track.globalIndex(), mapTOFBetaReassociated[std::make_pair(collision.globalIndex(), track.globalIndex())], mapTOFNsigmaReassociated[std::make_pair(collision.globalIndex(), track.globalIndex())]); + } + } // end of track loop + } // end of collision loop + } else { + for (const auto& collision : collisions) { + auto tracks_per_coll = tracks.sliceBy(perCol, collision.globalIndex()); + for (const auto& track : tracks_per_coll) { + if (!track.hasITS() || !track.hasTPC()) { // apply only minimal cut + continue; + } + mapTOFNsigmaReassociated[std::make_pair(collision.globalIndex(), track.globalIndex())] = track.tofNSigmaEl(); + mapTOFBetaReassociated[std::make_pair(collision.globalIndex(), track.globalIndex())] = track.beta(); + emtofs(collision.globalIndex(), track.globalIndex(), mapTOFBetaReassociated[std::make_pair(collision.globalIndex(), track.globalIndex())], mapTOFNsigmaReassociated[std::make_pair(collision.globalIndex(), track.globalIndex())]); + } + } // end of track loop + } // end of collision loop + } else { + if constexpr (withTTCA) { + for (const auto& collision : collisions) { + auto trackIdsThisCollision = trackIndices.sliceBy(trackIndicesPerCollision, collision.globalIndex()); + for (const auto& trackId : trackIdsThisCollision) { + auto track = trackId.template track_as(); + if (!track.hasITS() || !track.hasTPC()) { // apply only minimal cut + continue; + } + mapTOFNsigmaReassociated[std::make_pair(collision.globalIndex(), track.globalIndex())] = track.tofNSigmaEl(); + mapTOFBetaReassociated[std::make_pair(collision.globalIndex(), track.globalIndex())] = track.beta(); + emtofs(collision.globalIndex(), track.globalIndex(), mapTOFBetaReassociated[std::make_pair(collision.globalIndex(), track.globalIndex())], mapTOFNsigmaReassociated[std::make_pair(collision.globalIndex(), track.globalIndex())]); + } // end of track loop + } // end of collision loop + } else { + for (const auto& collision : collisions) { + auto tracks_per_coll = tracks.sliceBy(perCol, collision.globalIndex()); + for (const auto& track : tracks_per_coll) { + if (!track.hasITS() || !track.hasTPC()) { // apply only minimal cut + continue; + } + mapTOFNsigmaReassociated[std::make_pair(collision.globalIndex(), track.globalIndex())] = track.tofNSigmaEl(); + mapTOFBetaReassociated[std::make_pair(collision.globalIndex(), track.globalIndex())] = track.beta(); + emtofs(collision.globalIndex(), track.globalIndex(), mapTOFBetaReassociated[std::make_pair(collision.globalIndex(), track.globalIndex())], mapTOFNsigmaReassociated[std::make_pair(collision.globalIndex(), track.globalIndex())]); + } + } // end of track loop + } // end of collision loop + } + } + + template + bool checkTrack(TCollision const& collision, TTrack const& track) { if constexpr (isMC) { if (!track.has_mcParticle()) { return false; } + if (storeOnlyTrueElectronMC) { + const auto& mcParticle = track.template mcParticle_as(); + if (std::abs(mcParticle.pdgCode()) != 11) { + return false; + } + } + } + + float tofNSigmaEl = mapTOFNsigmaReassociated[std::make_pair(collision.globalIndex(), track.globalIndex())]; + if (requireTOF && !(track.hasTOF() && std::fabs(tofNSigmaEl) < maxTOFNsigmaEl)) { + return false; } if (!track.hasITS()) { return false; } - if (track.itsChi2NCl() < 0.f || maxchi2its < track.itsChi2NCl()) { + if (track.itsChi2NCl() < minchi2its || maxchi2its < track.itsChi2NCl()) { return false; } @@ -239,12 +377,12 @@ struct SkimmerPrimaryElectronFromDalitzEE { return false; } - if (!includeITSsa && !track.hasTPC()) { + if (!includeITSsa && (!track.hasITS() || !track.hasTPC())) { return false; } if (track.hasTPC()) { - if (track.tpcChi2NCl() < 0.f || maxchi2tpc < track.tpcChi2NCl()) { + if (track.tpcChi2NCl() < minchi2tpc || maxchi2tpc < track.tpcChi2NCl()) { return false; } @@ -265,80 +403,213 @@ struct SkimmerPrimaryElectronFromDalitzEE { } } - if (std::fabs(track.dcaXY()) > dca_xy_max || std::fabs(track.dcaZ()) > dca_z_max) { + o2::dataformats::DCA mDcaInfoCov; // for Track association + mDcaInfoCov.set(999, 999, 999, 999, 999); + auto trackParCov = getTrackParCov(track); + trackParCov.setPID(o2::track::PID::Electron); + mVtx.setPos({collision.posX(), collision.posY(), collision.posZ()}); + 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; + } + float dcaXY = mDcaInfoCov.getY(); + float dcaZ = mDcaInfoCov.getZ(); + + if (std::fabs(dcaXY) > dca_xy_max || std::fabs(dcaZ) > dca_z_max) { return false; } float dca3D = 999.f; - float det = track.cYY() * track.cZZ() - track.cZY() * track.cZY(); + float det = trackParCov.getSigmaY2() * trackParCov.getSigmaZ2() - trackParCov.getSigmaZY() * trackParCov.getSigmaZY(); if (det < 0) { dca3D = 999.f; } else { - float chi2 = (track.dcaXY() * track.dcaXY() * track.cZZ() + track.dcaZ() * track.dcaZ() * track.cYY() - 2. * track.dcaXY() * track.dcaZ() * track.cZY()) / det; + float chi2 = (dcaXY * dcaXY * trackParCov.getSigmaZ2() + dcaZ * dcaZ * trackParCov.getSigmaY2() - 2. * dcaXY * dcaZ * trackParCov.getSigmaZY()) / det; dca3D = std::sqrt(std::fabs(chi2) / 2.); } if (dca3D > dca_3d_sigma_max) { return false; } - if (std::fabs(track.eta()) > maxeta) { + if (trackParCov.getPt() < minpt || std::fabs(trackParCov.getEta()) > maxeta) { return false; } - if ((track.hasITS() && track.hasTPC()) && track.pt() < minpt) { + + int total_cluster_size = 0, nl = 0; + for (unsigned int layer = 0; layer < 7; layer++) { + int cluster_size_per_layer = track.itsClsSizeInLayer(layer); + if (cluster_size_per_layer > 0) { + nl++; + } + total_cluster_size += cluster_size_per_layer; + } + + if (maxMeanITSClusterSize < static_cast(total_cluster_size) / static_cast(nl) * std::cos(std::atan(trackParCov.getTgl()))) { return false; } + if ((track.hasITS() && !track.hasTPC() && !track.hasTRD() && !track.hasTOF()) && maxpt_itssa < track.pt()) { return false; } + // if ((track.hasITS() && !track.hasTPC() && !track.hasTRD() && !track.hasTOF()) && maxpt_itssa < track.pt()) { + // return false; + // } + return true; } - template - bool isElectron(TTrack const& track) + template + bool isElectron(TCollision const& collision, TTrack const& track) { - if (includeITSsa && (track.hasITS() && !track.hasTPC() && !track.hasTRD() && !track.hasTOF())) { - int totalClusterSize = 0, nl = 0; - for (unsigned int layer = 0; layer < 7; layer++) { // o2-linter: disable=magic-number (number of its layers) - int clusterSizePerLayer = track.itsClsSizeInLayer(layer); - if (clusterSizePerLayer > 0) { + if (includeITSsa && (track.hasITS() && !track.hasTPC() && !track.hasTRD() && !track.hasTOF())) [[unlikely]] { + int total_cluster_size = 0, nl = 0; + for (unsigned int layer = 0; layer < 7; layer++) { + int cluster_size_per_layer = track.itsClsSizeInLayer(layer); + if (cluster_size_per_layer > 0) { nl++; } - totalClusterSize += clusterSizePerLayer; + total_cluster_size += cluster_size_per_layer; } - return (maxMeanITSClusterSize > static_cast(totalClusterSize) / static_cast(nl) * std::cos(std::atan(track.tgl()))); + return (maxMeanITSClusterSize > static_cast(total_cluster_size) / static_cast(nl) * std::cos(std::atan(track.tgl()))); } + return isElectron_TPChadrej(collision, track) || isElectron_TOFreq(collision, track); + } + + template + bool isElectron_TPChadrej(TCollision const& collision, TTrack const& track) + { + float tofNSigmaEl = mapTOFNsigmaReassociated[std::make_pair(collision.globalIndex(), track.globalIndex())]; if (track.tpcNSigmaEl() < minTPCNsigmaEl || maxTPCNsigmaEl < track.tpcNSigmaEl()) { return false; } - if (minTPCNsigmaPi < track.tpcNSigmaPi() && track.tpcNSigmaPi() < maxTPCNsigmaPi) { + if (minTPCNsigmaPi < track.tpcNSigmaPi() && track.tpcNSigmaPi() < maxTPCNsigmaPi && track.tpcInnerParam() < max_pin_for_pion_rejection) { + return false; + } + if (minTPCNsigmaKa < track.tpcNSigmaKa() && track.tpcNSigmaKa() < maxTPCNsigmaKa) { return false; } - if (track.hasTOF() && (track.tofNSigmaEl() < minTOFNsigmaEl || maxTOFNsigmaEl < track.tofNSigmaEl())) { // TOFif + if (minTPCNsigmaPr < track.tpcNSigmaPr() && track.tpcNSigmaPr() < maxTPCNsigmaPr) { + return false; + } + if (track.hasTOF() && (maxTOFNsigmaEl < std::fabs(tofNSigmaEl))) { return false; } return true; } - template - void fillTrackTable(TCollision const& collision, TTrack const& track) + template + bool isElectron_TOFreq(TCollision const& collision, TTrack const& track) { - emprimaryelectrons(collision.globalIndex(), track.globalIndex(), track.sign(), - track.pt(), track.eta(), track.phi(), track.dcaXY(), track.dcaZ(), track.cYY(), track.cZY(), track.cZZ(), - track.tpcNClsFindable(), track.tpcNClsFindableMinusFound(), track.tpcNClsFindableMinusCrossedRows(), track.tpcNClsShared(), - track.tpcChi2NCl(), track.tpcInnerParam(), - track.tpcSignal(), track.tpcNSigmaEl(), track.tpcNSigmaPi(), - track.beta(), track.tofNSigmaEl(), - track.itsClusterSizes(), track.itsChi2NCl(), track.tofChi2(), track.detectorMap()); + float tofNSigmaEl = mapTOFNsigmaReassociated[std::make_pair(collision.globalIndex(), track.globalIndex())]; - if constexpr (isMC) { - emprimaryelectronsDeDxMC(track.mcTunedTPCSignal()); + if (minTPCNsigmaPi < track.tpcNSigmaPi() && track.tpcNSigmaPi() < maxTPCNsigmaPi && (min_pin_for_pion_rejection < track.tpcInnerParam() && track.tpcInnerParam() < max_pin_for_pion_rejection)) { + return false; + } + return minTPCNsigmaEl < track.tpcNSigmaEl() && track.tpcNSigmaEl() < maxTPCNsigmaEl && std::fabs(tofNSigmaEl) < maxTOFNsigmaEl; + } + + template + void fillTrackInfo(TCollision const& collision, TTracks const& tracks) + { + + for (const auto& track : tracks) { + if (!checkTrack(collision, track)) { + continue; + } + + if (!isElectron(collision, track)) { + continue; + } + + fillTrackHistograms(track); + fillTrackTable(collision, track); + } + } + + template + void fillPairInfo(TCollision const& collision, TTracks1 const& tracks1, TTracks2 const& tracks2) + { + if constexpr (pairtype == 0) { // ULS + for (const auto& [t1, t2] : combinations(CombinationsFullIndexPolicy(tracks1, tracks2))) { + + if (!checkTrack(collision, t1) || !checkTrack(collision, t2)) { + continue; + } + + if (!isElectron(collision, t1) || !isElectron(collision, t2)) { + continue; + } + + ROOT::Math::PtEtaPhiMVector v1(t1.pt(), t1.eta(), t1.phi(), o2::constants::physics::MassElectron); + ROOT::Math::PtEtaPhiMVector v2(t2.pt(), t2.eta(), t2.phi(), o2::constants::physics::MassElectron); + ROOT::Math::PtEtaPhiMVector v12 = v1 + v2; + float phiv = o2::aod::pwgem::dilepton::utils::pairutil::getPhivPair(t1.px(), t1.py(), t1.pz(), t2.px(), t2.py(), t2.pz(), t1.sign(), t2.sign(), dBz); + + fRegistry.fill(HIST("Pair/") + HIST(DileptonSigns[pairtype]) + HIST("hMvsPt"), v12.M(), v12.Pt()); + fRegistry.fill(HIST("Pair/") + HIST(DileptonSigns[pairtype]) + HIST("hMvsPhiV"), phiv, v12.M()); + + if (v12.M() > maxMee) { // don't store + continue; + } + + fRegistry.fill(HIST("Pair/") + HIST(DileptonSigns[pairtype]) + HIST("hMCutMvsPt"), v12.M(), v12.Pt()); + fRegistry.fill(HIST("Pair/") + HIST(DileptonSigns[pairtype]) + HIST("hMCutMvsPhiV"), phiv, v12.M()); + + if (v12.M() < slope * phiv + intercept) { + continue; + } + + fRegistry.fill(HIST("Pair/") + HIST(DileptonSigns[pairtype]) + HIST("hMPhiCutMvsPt"), v12.M(), v12.Pt()); + fRegistry.fill(HIST("Pair/") + HIST(DileptonSigns[pairtype]) + HIST("hMPhiCutMvsPhiV"), phiv, v12.M()); + + if (t1.sign() > 0) { // for positron + if (std::find(acceptedPosTrackIds_per_collision.begin(), acceptedPosTrackIds_per_collision.end(), t1.globalIndex()) == acceptedPosTrackIds_per_collision.end()) { + fillTrackHistograms(t1); + acceptedPosTrackIds_per_collision.emplace_back(t1.globalIndex()); + } + } else { // for electron + if (std::find(acceptedNegTrackIds_per_collision.begin(), acceptedNegTrackIds_per_collision.end(), t1.globalIndex()) == acceptedNegTrackIds_per_collision.end()) { + fillTrackHistograms(t1); + acceptedNegTrackIds_per_collision.emplace_back(t1.globalIndex()); + } + } + + if (t2.sign() > 0) { // for positron + if (std::find(acceptedPosTrackIds_per_collision.begin(), acceptedPosTrackIds_per_collision.end(), t2.globalIndex()) == acceptedPosTrackIds_per_collision.end()) { + fillTrackHistograms(t2); + acceptedPosTrackIds_per_collision.emplace_back(t2.globalIndex()); + } + } else { // for electron + if (std::find(acceptedNegTrackIds_per_collision.begin(), acceptedNegTrackIds_per_collision.end(), t2.globalIndex()) == acceptedNegTrackIds_per_collision.end()) { + fillTrackHistograms(t2); + acceptedNegTrackIds_per_collision.emplace_back(t2.globalIndex()); + } + } + } // end of ULS pairing + } else { // LS + for (auto& [t1, t2] : combinations(CombinationsStrictlyUpperIndexPolicy(tracks1, tracks2))) { + if (!checkTrack(collision, t1) || !checkTrack(collision, t2)) { + continue; + } + if (!isElectron(collision, t1) || !isElectron(collision, t2)) { + continue; + } + + ROOT::Math::PtEtaPhiMVector v1(t1.pt(), t1.eta(), t1.phi(), o2::constants::physics::MassElectron); + ROOT::Math::PtEtaPhiMVector v2(t2.pt(), t2.eta(), t2.phi(), o2::constants::physics::MassElectron); + ROOT::Math::PtEtaPhiMVector v12 = v1 + v2; + float phiv = o2::aod::pwgem::dilepton::utils::pairutil::getPhivPair(t1.px(), t1.py(), t1.pz(), t2.px(), t2.py(), t2.pz(), t1.sign(), t2.sign(), dBz); + fRegistry.fill(HIST("Pair/") + HIST(DileptonSigns[pairtype]) + HIST("hMvsPt"), v12.M(), v12.Pt()); + fRegistry.fill(HIST("Pair/") + HIST(DileptonSigns[pairtype]) + HIST("hMvsPhiV"), phiv, v12.M()); + } // end of LS pairing } } - template + template void fillTrackHistograms(TTrack const& track) { float mcTunedTPCSignal = 0.f; @@ -408,235 +679,237 @@ struct SkimmerPrimaryElectronFromDalitzEE { fRegistry.fill(HIST("Track/hMeanClusterSizeITSob"), track.p(), static_cast(totalClusterSizeOb) / static_cast(nlOb) * std::cos(std::atan(track.tgl()))); } - template - void fillPairInfo(TCollision const& collision, TTracks1 const& tracks1, TTracks2 const& tracks2) + template + void fillTrackTable(TCollision const& collision, TTrack const& track) { - if constexpr (pairtype == 0) { // ULS - for (const auto& [t1, t2] : combinations(CombinationsFullIndexPolicy(tracks1, tracks2))) { - if (!checkTrack(collision, t1) || !checkTrack(collision, t2)) { - continue; - } - if (!isElectron(t1) || !isElectron(t2)) { - continue; - } - - ROOT::Math::PtEtaPhiMVector v1(t1.pt(), t1.eta(), t1.phi(), o2::constants::physics::MassElectron); - ROOT::Math::PtEtaPhiMVector v2(t2.pt(), t2.eta(), t2.phi(), o2::constants::physics::MassElectron); - ROOT::Math::PtEtaPhiMVector v12 = v1 + v2; - float phiv = o2::aod::pwgem::dilepton::utils::pairutil::getPhivPair(t1.px(), t1.py(), t1.pz(), t2.px(), t2.py(), t2.pz(), t1.sign(), t2.sign(), dBz); - fRegistry.fill(HIST("Pair/") + HIST(DileptonSigns[pairtype]) + HIST("hMvsPt"), v12.M(), v12.Pt()); - fRegistry.fill(HIST("Pair/") + HIST(DileptonSigns[pairtype]) + HIST("hMvsPhiV"), phiv, v12.M()); - - if (v12.M() > maxMee) { // don't store - continue; - } - - if (v12.M() < slope * phiv + intercept) { - continue; - } - - if (t1.sign() > 0) { // for positron - if (std::find(acceptedPosTrackIdsPerCollision.begin(), acceptedPosTrackIdsPerCollision.end(), t1.globalIndex()) == acceptedPosTrackIdsPerCollision.end()) { - fillTrackHistograms(t1); - acceptedPosTrackIdsPerCollision.emplace_back(t1.globalIndex()); - } - } else { // for electron - if (std::find(acceptedNegTrackIdsPerCollision.begin(), acceptedNegTrackIdsPerCollision.end(), t1.globalIndex()) == acceptedNegTrackIdsPerCollision.end()) { - fillTrackHistograms(t1); - acceptedNegTrackIdsPerCollision.emplace_back(t1.globalIndex()); - } - } - - if (t2.sign() > 0) { // for positron - if (std::find(acceptedPosTrackIdsPerCollision.begin(), acceptedPosTrackIdsPerCollision.end(), t2.globalIndex()) == acceptedPosTrackIdsPerCollision.end()) { - fillTrackHistograms(t2); - acceptedPosTrackIdsPerCollision.emplace_back(t2.globalIndex()); - } - } else { // for electron - if (std::find(acceptedNegTrackIdsPerCollision.begin(), acceptedNegTrackIdsPerCollision.end(), t2.globalIndex()) == acceptedNegTrackIdsPerCollision.end()) { - fillTrackHistograms(t2); - acceptedNegTrackIdsPerCollision.emplace_back(t2.globalIndex()); - } - } - } // end of ULS pairing - } else { // LS - for (const auto& [t1, t2] : combinations(CombinationsStrictlyUpperIndexPolicy(tracks1, tracks2))) { - if (!checkTrack(collision, t1) || !checkTrack(collision, t2)) { - continue; - } - if (!isElectron(t1) || !isElectron(t2)) { - continue; - } + emprimaryelectrons(collision.globalIndex(), track.globalIndex(), track.sign(), + track.pt(), track.eta(), track.phi(), track.dcaXY(), track.dcaZ(), track.cYY(), track.cZY(), track.cZZ(), + track.tpcNClsFindable(), track.tpcNClsFindableMinusFound(), track.tpcNClsFindableMinusCrossedRows(), track.tpcNClsShared(), + track.tpcChi2NCl(), track.tpcInnerParam(), + track.tpcSignal(), track.tpcNSigmaEl(), track.tpcNSigmaPi(), + track.beta(), track.tofNSigmaEl(), + track.itsClusterSizes(), track.itsChi2NCl(), track.tofChi2(), track.detectorMap()); - ROOT::Math::PtEtaPhiMVector v1(t1.pt(), t1.eta(), t1.phi(), o2::constants::physics::MassElectron); - ROOT::Math::PtEtaPhiMVector v2(t2.pt(), t2.eta(), t2.phi(), o2::constants::physics::MassElectron); - ROOT::Math::PtEtaPhiMVector v12 = v1 + v2; - float phiv = o2::aod::pwgem::dilepton::utils::pairutil::getPhivPair(t1.px(), t1.py(), t1.pz(), t2.px(), t2.py(), t2.pz(), t1.sign(), t2.sign(), dBz); - fRegistry.fill(HIST("Pair/") + HIST(DileptonSigns[pairtype]) + HIST("hMvsPt"), v12.M(), v12.Pt()); - fRegistry.fill(HIST("Pair/") + HIST(DileptonSigns[pairtype]) + HIST("hMvsPhiV"), phiv, v12.M()); - } // end of LS pairing + if constexpr (isMC) { + emprimaryelectronsDeDxMC(track.mcTunedTPCSignal()); } } - std::vector acceptedPosTrackIdsPerCollision; - std::vector acceptedNegTrackIdsPerCollision; - std::vector> storedTrackIds; - Filter trackFilter = minpt < o2::aod::track::pt && nabs(o2::aod::track::eta) < maxeta && o2::aod::track::itsChi2NCl < maxchi2its && ncheckbit(aod::track::v001::detectorMap, (uint8_t)o2::aod::track::ITS) == true && nabs(o2::aod::track::dcaXY) < dca_xy_max && nabs(o2::aod::track::dcaZ) < dca_z_max; - using MyFilteredTracks = soa::Filtered; - Partition posTracks = o2::aod::track::signed1Pt > 0.f; - Partition negTracks = o2::aod::track::signed1Pt < 0.f; + std::vector acceptedPosTrackIds_per_collision; + std::vector acceptedNegTrackIds_per_collision; + std::vector acceptedTrackIds_per_collision; + std::vector> stored_trackIds; + // Filter trackFilter = minpt < o2::aod::track::pt && nabs(o2::aod::track::eta) < maxeta && o2::aod::track::itsChi2NCl < maxchi2its && ncheckbit(aod::track::v001::detectorMap, (uint8_t)o2::aod::track::ITS) == true && nabs(o2::aod::track::dcaXY) < dca_xy_max && nabs(o2::aod::track::dcaZ) < dca_z_max; + // Filter trackFilter + // using MyFilteredTracks = soa::Filtered; + Partition posTracks = o2::aod::track::signed1Pt > 0.f; + Partition negTracks = o2::aod::track::signed1Pt < 0.f; + + std::unordered_map mapCollisionTime; + std::unordered_map mapCollisionTimeError; + + std::map, float> mapTOFNsigmaReassociated; // map pair(collisionId, trackId) -> tof n sigma + std::map, float> mapTOFBetaReassociated; // map pair(collisionId, trackId) -> tof beta // ---------- for data ---------- - void processRec(MyCollisions const& collisions, aod::BCsWithTimestamps const&, MyFilteredTracks const& tracks, aod::V0PhotonsKF const& v0photons) + void processRec(MyCollisions const& collisions, aod::BCsWithTimestamps const& bcs, MyTracks const& tracks, aod::V0PhotonsKF const& v0photons, aod::TrackAssoc const& trackIndices) { - storedTrackIds.reserve(tracks.size()); + initCCDB(bcs.iteratorAt(0)); + mTOFResponse->processSetup(bcs.iteratorAt(0)); + + for (const auto& track : tracks) { + if (mapCollisionTime.find(track.collisionId()) == mapCollisionTime.end()) { + // LOGF(info, "track.collisionId() = %d, track.tofEvTime() = %f, track.tofEvTimeErr() = %f", track.collisionId(), track.tofEvTime(), track.tofEvTimeErr()); + mapCollisionTime[track.collisionId()] = track.tofEvTime(); + mapCollisionTimeError[track.collisionId()] = track.tofEvTimeErr(); + } + } + calculateTOFNSigmaWithReassociation(collisions, bcs, tracks, trackIndices); for (const auto& collision : collisions) { - auto bc = collision.template foundBC_as(); - initCCDB(bc); if (!collision.isSelected()) { continue; } - const auto& v0photonsPerColl = v0photons.sliceBy(perColPcm, collision.globalIndex()); - const auto& posTracksPerColl = posTracks->sliceByCached(o2::aod::track::collisionId, collision.globalIndex(), cache); - const auto& negTracksPerColl = negTracks->sliceByCached(o2::aod::track::collisionId, collision.globalIndex(), cache); - acceptedPosTrackIdsPerCollision.reserve(posTracksPerColl.size()); - acceptedNegTrackIdsPerCollision.reserve(negTracksPerColl.size()); - - fillPairInfo(collision, posTracksPerColl, negTracksPerColl); // ULS - if (fillLS) { - fillPairInfo(collision, posTracksPerColl, posTracksPerColl); // LS++ - fillPairInfo(collision, negTracksPerColl, negTracksPerColl); // LS-- - } - - if ((v0photonsPerColl.size() >= 1 && !acceptedPosTrackIdsPerCollision.empty() && !acceptedNegTrackIdsPerCollision.empty()) || (acceptedPosTrackIdsPerCollision.size() >= 2 && acceptedNegTrackIdsPerCollision.size() >= 2)) { // o2-linter: disable=magic-number (check to see if we have enough particles that it make sense to fill the tables) - // LOGF(info, "v0photonsPerColl.size() = %d, acceptedPosTrackIdsPerCollision.size() = %d, acceptedNegTrackIdsPerCollision.size() = %d", v0photonsPerColl.size(), acceptedPosTrackIdsPerCollision.size(), acceptedNegTrackIdsPerCollision.size()); - for (const auto& posId : acceptedPosTrackIdsPerCollision) { - const auto& pos = tracks.rawIteratorAt(posId); - fillTrackTable(collision, pos); + const auto& v0photons_per_coll = v0photons.sliceBy(perCol_pcm, collision.globalIndex()); + const auto& posTracks_per_coll = posTracks->sliceByCached(o2::aod::track::collisionId, collision.globalIndex(), cache); + const auto& negTracks_per_coll = negTracks->sliceByCached(o2::aod::track::collisionId, collision.globalIndex(), cache); + const auto& slicedTracks = tracks.sliceBy(perTracksCollision, collision.globalIndex()); + acceptedPosTrackIds_per_collision.reserve(posTracks_per_coll.size()); + acceptedNegTrackIds_per_collision.reserve(negTracks_per_coll.size()); + + if (!fillWithPairs) { + fillTrackInfo(collision, slicedTracks); + } else { + fillPairInfo(collision, posTracks_per_coll, negTracks_per_coll); // ULS + if (fillLS) { + fillPairInfo(collision, posTracks_per_coll, posTracks_per_coll); // LS++ + fillPairInfo(collision, negTracks_per_coll, negTracks_per_coll); // LS-- } - for (const auto& eleId : acceptedNegTrackIdsPerCollision) { - const auto& ele = tracks.rawIteratorAt(eleId); - fillTrackTable(collision, ele); + + if ((v0photons_per_coll.size() >= 1 && !acceptedPosTrackIds_per_collision.empty() && !acceptedNegTrackIds_per_collision.empty()) || (acceptedPosTrackIds_per_collision.size() >= 2 && acceptedNegTrackIds_per_collision.size() >= 2)) { + for (const auto& posId : acceptedPosTrackIds_per_collision) { + const auto& pos = tracks.rawIteratorAt(posId); + fillTrackTable(collision, pos); + } + for (const auto& eleId : acceptedNegTrackIds_per_collision) { + const auto& ele = tracks.rawIteratorAt(eleId); + fillTrackTable(collision, ele); + } } } - acceptedPosTrackIdsPerCollision.clear(); - acceptedPosTrackIdsPerCollision.shrink_to_fit(); - acceptedNegTrackIdsPerCollision.clear(); - acceptedNegTrackIdsPerCollision.shrink_to_fit(); - } // end of collision loop + acceptedPosTrackIds_per_collision.clear(); + acceptedNegTrackIds_per_collision.clear(); - storedTrackIds.clear(); - storedTrackIds.shrink_to_fit(); + } // end of collision loop } - PROCESS_SWITCH(SkimmerPrimaryElectronFromDalitzEE, processRec, "process reconstructed info only", true); // standalone + PROCESS_SWITCH(skimmerPrimaryElectronFromDalitzEE, processRec, "process reconstructed info only", false); // standalone - // void processRec_SWT(MyCollisionsWithSWT const& collisions, aod::BCsWithTimestamps const&, MyFilteredTracks const& tracks, aod::V0PhotonsKF const& v0photons) + // void processRec_SWT(MyCollisionsWithSWT const& collisions, aod::BCsWithTimestamps const& bcs, MyTracks const& tracks, aod::V0PhotonsKF const& v0photons, aod::TrackAssoc const& trackIndices) // { - // storedTrackIds.reserve(tracks.size()); + // initCCDB(bcs.iteratorAt(0)); + + // mTOFResponse->processSetup(bcs.iteratorAt(0)); + // for (const auto& track : tracks) { + // if (mapCollisionTime.find(track.collisionId()) == mapCollisionTime.end()) { + // // LOGF(info, "track.collisionId() = %d, track.tofEvTime() = %f, track.tofEvTimeErr() = %f", track.collisionId(), track.tofEvTime(), track.tofEvTimeErr()); + // mapCollisionTime[track.collisionId()] = track.tofEvTime(); + // mapCollisionTimeError[track.collisionId()] = track.tofEvTimeErr(); + // } + // } + // calculateTOFNSigmaWithReassociation(collisions, bcs, tracks, trackIndices); // for (const auto& collision : collisions) { - // auto bc = collision.template foundBC_as(); - // initCCDB(bc); // if (!collision.isSelected()) { // continue; // } - // if (collision.triggerMask_raw() == 0) { + // if (collision.swtaliastmp_raw() == 0) { // continue; // } - // const auto& v0photonsPerColl = v0photons.sliceBy(perColPcm, collision.globalIndex()); - // const auto& posTracksPerColl = posTracks->sliceByCached(o2::aod::track::collisionId, collision.globalIndex(), cache); - // const auto& negTracksPerColl = negTracks->sliceByCached(o2::aod::track::collisionId, collision.globalIndex(), cache); - // acceptedPosTrackIdsPerCollision.reserve(posTracksPerColl.size()); - // acceptedNegTrackIdsPerCollision.reserve(negTracksPerColl.size()); - - // fillPairInfo(collision, posTracksPerColl, negTracksPerColl); // ULS - // if (fillLS) { - // fillPairInfo(collision, posTracksPerColl, posTracksPerColl); // LS++ - // fillPairInfo(collision, negTracksPerColl, negTracksPerColl); // LS-- - // } - - // if ((v0photonsPerColl.size() >= 1 && acceptedPosTrackIdsPerCollision.size() >= 1 && acceptedNegTrackIdsPerCollision.size() >= 1) || (acceptedPosTrackIdsPerCollision.size() >= 2 && acceptedNegTrackIdsPerCollision.size() >= 2)) { - // // LOGF(info, "v0photonsPerColl.size() = %d, acceptedPosTrackIdsPerCollision.size() = %d, acceptedNegTrackIdsPerCollision.size() = %d", v0photonsPerColl.size(), acceptedPosTrackIdsPerCollision.size(), acceptedNegTrackIdsPerCollision.size()); - // for (const auto& posId : acceptedPosTrackIdsPerCollision) { - // const auto& pos = tracks.rawIteratorAt(posId); - // fillTrackTable(collision, pos); + // const auto& v0photons_per_coll = v0photons.sliceBy(perCol_pcm, collision.globalIndex()); + // const auto& posTracks_per_coll = posTracks->sliceByCached(o2::aod::track::collisionId, collision.globalIndex(), cache); + // const auto& negTracks_per_coll = negTracks->sliceByCached(o2::aod::track::collisionId, collision.globalIndex(), cache); + // const auto& slicedTracks = tracks.sliceBy(perTracksCollision, collision.globalIndex()); + // acceptedPosTrackIds_per_collision.reserve(posTracks_per_coll.size()); + // acceptedNegTrackIds_per_collision.reserve(negTracks_per_coll.size()); + + // if(!fillWithPairs){ + // fillTrackInfo(collision, slicedTracks); + // }else{ + // fillPairInfo(collision, posTracks_per_coll, negTracks_per_coll); // ULS + // if (fillLS) { + // fillPairInfo(collision, posTracks_per_coll, posTracks_per_coll); // LS++ + // fillPairInfo(collision, negTracks_per_coll, negTracks_per_coll); // LS-- // } - // for (const auto& eleId : acceptedNegTrackIdsPerCollision) { - // const auto& ele = tracks.rawIteratorAt(eleId); - // fillTrackTable(collision, ele); + + // if ((v0photons_per_coll.size() >= 1 && !acceptedPosTrackIds_per_collision.empty() && !acceptedNegTrackIds_per_collision.empty()) || (acceptedPosTrackIds_per_collision.size() >= 2 && acceptedNegTrackIds_per_collision.size() >= 2)) { + // // LOGF(info, "v0photons_per_coll.size() = %d, acceptedPosTrackIds_per_collision.size() = %d, acceptedNegTrackIds_per_collision.size() = %d", v0photons_per_coll.size(), acceptedPosTrackIds_per_collision.size(), acceptedNegTrackIds_per_collision.size()); + // for (const auto& posId : acceptedPosTrackIds_per_collision) { + // const auto& pos = tracks.rawIteratorAt(posId); + // fillTrackTable(collision, pos); + // } + // for (const auto& eleId : acceptedNegTrackIds_per_collision) { + // const auto& ele = tracks.rawIteratorAt(eleId); + // fillTrackTable(collision, ele); + // } // } // } - // acceptedPosTrackIdsPerCollision.clear(); - // acceptedPosTrackIdsPerCollision.shrink_to_fit(); - // acceptedNegTrackIdsPerCollision.clear(); - // acceptedNegTrackIdsPerCollision.shrink_to_fit(); - // } // end of collision loop + // acceptedPosTrackIds_per_collision.clear(); + // acceptedNegTrackIds_per_collision.clear(); - // storedTrackIds.clear(); - // storedTrackIds.shrink_to_fit(); + // } // end of collision loop // } - // PROCESS_SWITCH(SkimmerPrimaryElectronFromDalitzEE, processRec_SWT, "process reconstructed info with CEFP", false); // with cefp + // PROCESS_SWITCH(skimmerPrimaryElectronFromDalitzEE, processRec_SWT, "process reconstructed info with CEFP", false); // with cefp - using MyFilteredTracksMC = soa::Filtered; - Partition posTracksMC = o2::aod::track::signed1Pt > 0.f; - Partition negTracksMC = o2::aod::track::signed1Pt < 0.f; + // using MyFilteredTracksMC = soa::Filtered; + Partition posTracksMC = o2::aod::track::signed1Pt > 0.f; + Partition negTracksMC = o2::aod::track::signed1Pt < 0.f; // ---------- for MC ---------- - void processMC(MyCollisionsMC const& collisions, aod::McCollisions const&, aod::BCsWithTimestamps const&, MyFilteredTracksMC const& tracks, aod::V0PhotonsKF const& v0photons) + void processMC(MyCollisionsMC const& collisions, aod::McCollisions const&, aod::BCsWithTimestamps const& bcs, MyTracksMC const& tracks, aod::V0PhotonsKF const& v0photons, aod::TrackAssoc const& trackIndices) { - storedTrackIds.reserve(tracks.size()); + uint64_t nCollisSel = 0; + uint64_t nNoMcColl = 0; + uint64_t nProcessedCollisions = 0; + uint64_t nCheckTrack = 0; + + initCCDB(bcs.iteratorAt(0)); + + mTOFResponse->processSetup(bcs.iteratorAt(0)); + for (const auto& track : tracks) { + if (mapCollisionTime.find(track.collisionId()) == mapCollisionTime.end()) { + // LOGF(info, "track.collisionId() = %d, track.tofEvTime() = %f, track.tofEvTimeErr() = %f", track.collisionId(), track.tofEvTime(), track.tofEvTimeErr()); + mapCollisionTime[track.collisionId()] = track.tofEvTime(); + mapCollisionTimeError[track.collisionId()] = track.tofEvTimeErr(); + } + } + calculateTOFNSigmaWithReassociation(collisions, bcs, tracks, trackIndices); for (const auto& collision : collisions) { - auto bc = collision.template foundBC_as(); - initCCDB(bc); if (!collision.has_mcCollision()) { + nNoMcColl++; continue; } + if (!collision.isSelected()) { + nCollisSel++; continue; } - - const auto& v0photonsPerColl = v0photons.sliceBy(perColPcm, collision.globalIndex()); - const auto& posTracksPerColl = posTracksMC->sliceByCached(o2::aod::track::collisionId, collision.globalIndex(), cache); - const auto& negTracksPerColl = negTracksMC->sliceByCached(o2::aod::track::collisionId, collision.globalIndex(), cache); - acceptedPosTrackIdsPerCollision.reserve(posTracksPerColl.size()); - acceptedNegTrackIdsPerCollision.reserve(negTracksPerColl.size()); - - fillPairInfo(collision, posTracksPerColl, negTracksPerColl); // ULS - if (fillLS) { - fillPairInfo(collision, posTracksPerColl, posTracksPerColl); // LS++ - fillPairInfo(collision, negTracksPerColl, negTracksPerColl); // LS-- - } - if ((v0photonsPerColl.size() >= 1 && !acceptedPosTrackIdsPerCollision.empty() && !acceptedNegTrackIdsPerCollision.empty()) || (acceptedPosTrackIdsPerCollision.size() >= 2 && acceptedNegTrackIdsPerCollision.size() >= 2)) { // o2-linter: disable=magic-number (check to see if we have enough particles that it make sense to fill the tables) - // LOGF(info, "v0photonsPerColl.size() = %d, acceptedPosTrackIdsPerCollision.size() = %d, acceptedNegTrackIdsPerCollision.size() = %d", v0photonsPerColl.size(), acceptedPosTrackIdsPerCollision.size(), acceptedNegTrackIdsPerCollision.size()); - for (const auto& posId : acceptedPosTrackIdsPerCollision) { - const auto& pos = tracks.rawIteratorAt(posId); - fillTrackTable(collision, pos); + nProcessedCollisions++; + + const auto& v0photons_per_coll = v0photons.sliceBy(perCol_pcm, collision.globalIndex()); + const auto& posTracks_per_coll = posTracksMC->sliceByCached(o2::aod::track::collisionId, collision.globalIndex(), cache); + const auto& negTracks_per_coll = negTracksMC->sliceByCached(o2::aod::track::collisionId, collision.globalIndex(), cache); + const auto& slicedTracks = tracks.sliceBy(perTracksCollision, collision.globalIndex()); + acceptedPosTrackIds_per_collision.reserve(posTracks_per_coll.size()); + acceptedNegTrackIds_per_collision.reserve(negTracks_per_coll.size()); + acceptedTrackIds_per_collision.reserve(2 * (negTracks_per_coll.size())); + + if (!fillWithPairs) { + fillTrackInfo(collision, slicedTracks); + } else { + fillPairInfo(collision, posTracks_per_coll, negTracks_per_coll); // ULS + if (fillLS) { + fillPairInfo(collision, posTracks_per_coll, posTracks_per_coll); // LS++ + fillPairInfo(collision, negTracks_per_coll, negTracks_per_coll); // LS-- } - for (const auto& eleId : acceptedNegTrackIdsPerCollision) { - const auto& ele = tracks.rawIteratorAt(eleId); - fillTrackTable(collision, ele); + if ((acceptedPosTrackIds_per_collision.empty() && acceptedNegTrackIds_per_collision.empty()) || (acceptedPosTrackIds_per_collision.size() >= 2 && acceptedNegTrackIds_per_collision.size() >= 2)) { // v0photons_per_coll.size() >= 1 && + for (const auto& posId : acceptedPosTrackIds_per_collision) { + const auto& pos = tracks.rawIteratorAt(posId); + fillTrackTable(collision, pos); + } + for (const auto& eleId : acceptedNegTrackIds_per_collision) { + const auto& ele = tracks.rawIteratorAt(eleId); + fillTrackTable(collision, ele); + } } - } + } // end of fill loop + + acceptedPosTrackIds_per_collision.clear(); + acceptedNegTrackIds_per_collision.clear(); - acceptedPosTrackIdsPerCollision.clear(); - acceptedPosTrackIdsPerCollision.shrink_to_fit(); - acceptedNegTrackIdsPerCollision.clear(); - acceptedNegTrackIdsPerCollision.shrink_to_fit(); } // end of collision loop - storedTrackIds.clear(); - storedTrackIds.shrink_to_fit(); + LOG(info) << "Total collisions: " << collisions.size(); + LOG(info) << "No MC collision: " << nNoMcColl; + LOG(info) << "Rejected by selection: " << nCollisSel; + LOG(info) << "Processed collisions: " << nProcessedCollisions; + LOG(info) << "Invalid track: " << nCheckTrack; } - PROCESS_SWITCH(SkimmerPrimaryElectronFromDalitzEE, processMC, "process reconstructed and MC info ", false); + PROCESS_SWITCH(skimmerPrimaryElectronFromDalitzEE, processMC, "process reconstructed and MC info ", true); }; +// WorkflowSpec defineDataProcessing(ConfigContext const& context) +// { +// return WorkflowSpec{adaptAnalysisTask(context, TaskName{"skimmer-primary-electron-from-dalitzee"})}; +// } WorkflowSpec defineDataProcessing(ConfigContext const& context) { - return WorkflowSpec{adaptAnalysisTask(context, TaskName{"skimmer-primary-electron-from-dalitzee"})}; + o2::pid::tof::TOFResponseImpl::metadataInfo.initMetadata(context); + + return WorkflowSpec{ + adaptAnalysisTask(context, TaskName{"skimmer-primary-electron-from-dalitzee"})}; }