From 12f390360393e6acc1bfa1531772d60f396a853c Mon Sep 17 00:00:00 2001 From: Daiki Sekihata Date: Fri, 26 Jul 2024 03:23:48 +0200 Subject: [PATCH] PWGEM/Dilepton: add single electron QA task (#6997) --- PWGEM/Dilepton/Tasks/CMakeLists.txt | 15 + PWGEM/Dilepton/Tasks/dielectronQC.cxx | 148 ++----- PWGEM/Dilepton/Tasks/dielectronQCMC.cxx | 222 ++--------- PWGEM/Dilepton/Tasks/dimuonQC.cxx | 58 +-- PWGEM/Dilepton/Tasks/dimuonQCMC.cxx | 127 +----- PWGEM/Dilepton/Tasks/singleElectronQC.cxx | 335 ++++++++++++++++ PWGEM/Dilepton/Tasks/singleElectronQCMC.cxx | 115 ++---- PWGEM/Dilepton/Tasks/singleMuonQC.cxx | 252 ++++++++++++ PWGEM/Dilepton/Tasks/singleMuonQCMC.cxx | 409 ++++++++++++++++++++ 9 files changed, 1131 insertions(+), 550 deletions(-) create mode 100644 PWGEM/Dilepton/Tasks/singleElectronQC.cxx create mode 100644 PWGEM/Dilepton/Tasks/singleMuonQC.cxx create mode 100644 PWGEM/Dilepton/Tasks/singleMuonQCMC.cxx diff --git a/PWGEM/Dilepton/Tasks/CMakeLists.txt b/PWGEM/Dilepton/Tasks/CMakeLists.txt index 6cb75035f2e..114cc8676e9 100644 --- a/PWGEM/Dilepton/Tasks/CMakeLists.txt +++ b/PWGEM/Dilepton/Tasks/CMakeLists.txt @@ -40,11 +40,26 @@ o2physics_add_dpl_workflow(table-reader-barrel PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2::DetectorsBase O2Physics::AnalysisCore O2Physics::PWGDQCore COMPONENT_NAME Analysis) +o2physics_add_dpl_workflow(single-electron-qc + SOURCES singleElectronQC.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2::DetectorsBase O2Physics::AnalysisCore O2Physics::MLCore O2Physics::PWGEMDileptonCore + COMPONENT_NAME Analysis) + o2physics_add_dpl_workflow(single-electron-qc-mc SOURCES singleElectronQCMC.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2::DetectorsBase O2Physics::AnalysisCore O2Physics::MLCore O2Physics::PWGEMDileptonCore COMPONENT_NAME Analysis) +o2physics_add_dpl_workflow(single-muon-qc + SOURCES singleMuonQC.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2::DetectorsBase O2Physics::AnalysisCore O2Physics::PWGEMDileptonCore + COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(single-muon-qc-mc + SOURCES singleMuonQCMC.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2::DetectorsBase O2Physics::AnalysisCore O2Physics::PWGEMDileptonCore + COMPONENT_NAME Analysis) + o2physics_add_dpl_workflow(dielectron-qc SOURCES dielectronQC.cxx PUBLIC_LINK_LIBRARIES O2::Framework O2::DetectorsBase O2Physics::AnalysisCore O2Physics::MLCore O2::DCAFitter O2Physics::PWGEMDileptonCore diff --git a/PWGEM/Dilepton/Tasks/dielectronQC.cxx b/PWGEM/Dilepton/Tasks/dielectronQC.cxx index 9079da7cf46..4a524b961f6 100644 --- a/PWGEM/Dilepton/Tasks/dielectronQC.cxx +++ b/PWGEM/Dilepton/Tasks/dielectronQC.cxx @@ -102,7 +102,7 @@ struct dielectronQC { Configurable cfgOccupancyMax{"cfgOccupancyMax", 1000000000, "max. occupancy"}; } eventcuts; - DielectronCut fDielectonCut; + DielectronCut fDielectronCut; struct : ConfigurableGroup { std::string prefix = "dielectroncut_group"; Configurable cfg_min_mass{"cfg_min_mass", 0.0, "min mass"}; @@ -175,6 +175,14 @@ struct dielectronQC { void init(InitContext& /*context*/) { + mRunNumber = 0; + d_bz = 0; + + ccdb->setURL(ccdburl); + ccdb->setCaching(true); + ccdb->setLocalObjectValidityChecking(); + ccdb->setFatalWhenNull(false); + zvtx_bin_edges = std::vector(ConfVtxBins.value.begin(), ConfVtxBins.value.end()); zvtx_bin_edges.erase(zvtx_bin_edges.begin()); @@ -194,14 +202,6 @@ struct dielectronQC { DefineDileptonCut(); addhistograms(); - mRunNumber = 0; - d_bz = 0; - - ccdb->setURL(ccdburl); - ccdb->setCaching(true); - ccdb->setLocalObjectValidityChecking(); - ccdb->setFatalWhenNull(false); - // fitter.setPropagateToPCA(true); // fitter.setMaxR(5.f); // fitter.setMinParamChange(1e-3); @@ -359,38 +359,6 @@ struct dielectronQC { fRegistry.addClone("Pair/same/", "Pair/mix/"); } - // for track info - fRegistry.add("Track/hPt", "pT;p_{T} (GeV/c)", kTH1F, {{1000, 0.0f, 10}}, false); - fRegistry.add("Track/hQoverPt", "q/pT;q/p_{T} (GeV/c)^{-1}", kTH1F, {{400, -20, 20}}, false); - fRegistry.add("Track/hEtaPhi", "#eta vs. #varphi;#varphi (rad.);#eta", kTH2F, {{180, 0, 2 * M_PI}, {40, -2.0f, 2.0f}}, false); - fRegistry.add("Track/hDCAxyz", "DCA xy vs. z;DCA_{xy} (cm);DCA_{z} (cm)", kTH2F, {{200, -1.0f, 1.0f}, {200, -1.0f, 1.0f}}, false); - fRegistry.add("Track/hDCAxyzSigma", "DCA xy vs. z;DCA_{xy} (#sigma);DCA_{z} (#sigma)", kTH2F, {{200, -10.0f, 10.0f}, {200, -10.0f, 10.0f}}, false); - fRegistry.add("Track/hDCA3DSigma", "DCA 3D;DCA_{3D} (#sigma);", kTH1F, {{100, 0.0f, 10.0f}}, false); - fRegistry.add("Track/hDCAxyRes_Pt", "DCA_{xy} resolution vs. pT;p_{T} (GeV/c);DCA_{xy} resolution (#mum)", kTH2F, {{1000, 0, 10}, {500, 0., 500}}, false); - fRegistry.add("Track/hDCAzRes_Pt", "DCA_{z} resolution vs. pT;p_{T} (GeV/c);DCA_{z} resolution (#mum)", kTH2F, {{1000, 0, 10}, {500, 0., 500}}, false); - fRegistry.add("Track/hNclsTPC", "number of TPC clusters", kTH1F, {{161, -0.5, 160.5}}, false); - fRegistry.add("Track/hNcrTPC", "number of TPC crossed rows", kTH1F, {{161, -0.5, 160.5}}, false); - fRegistry.add("Track/hChi2TPC", "chi2/number of TPC clusters", kTH1F, {{100, 0, 10}}, false); - fRegistry.add("Track/hTPCdEdx", "TPC dE/dx;p_{in} (GeV/c);TPC dE/dx (a.u.)", kTH2F, {{1000, 0, 10}, {200, 0, 200}}, false); - fRegistry.add("Track/hTPCNsigmaEl", "TPC n sigma el;p_{in} (GeV/c);n #sigma_{e}^{TPC}", kTH2F, {{1000, 0, 10}, {100, -5, +5}}, false); - fRegistry.add("Track/hTPCNsigmaMu", "TPC n sigma mu;p_{in} (GeV/c);n #sigma_{#mu}^{TPC}", kTH2F, {{1000, 0, 10}, {100, -5, +5}}, false); - fRegistry.add("Track/hTPCNsigmaPi", "TPC n sigma pi;p_{in} (GeV/c);n #sigma_{#pi}^{TPC}", kTH2F, {{1000, 0, 10}, {100, -5, +5}}, false); - fRegistry.add("Track/hTPCNsigmaKa", "TPC n sigma ka;p_{in} (GeV/c);n #sigma_{K}^{TPC}", kTH2F, {{1000, 0, 10}, {100, -5, +5}}, false); - fRegistry.add("Track/hTPCNsigmaPr", "TPC n sigma pr;p_{in} (GeV/c);n #sigma_{p}^{TPC}", kTH2F, {{1000, 0, 10}, {100, -5, +5}}, false); - fRegistry.add("Track/hTOFbeta", "TOF beta;p_{in} (GeV/c);#beta", kTH2F, {{1000, 0, 10}, {600, 0, 1.2}}, false); - fRegistry.add("Track/h1overTOFbeta", "TOF beta;p_{in} (GeV/c);#beta", kTH2F, {{1000, 0, 10}, {1000, 0.8, 1.8}}, false); - fRegistry.add("Track/hTOFNsigmaEl", "TOF n sigma el;p_{in} (GeV/c);n #sigma_{e}^{TOF}", kTH2F, {{1000, 0, 10}, {100, -5, +5}}, false); - fRegistry.add("Track/hTOFNsigmaMu", "TOF n sigma mu;p_{in} (GeV/c);n #sigma_{#mu}^{TOF}", kTH2F, {{1000, 0, 10}, {100, -5, +5}}, false); - fRegistry.add("Track/hTOFNsigmaPi", "TOF n sigma pi;p_{in} (GeV/c);n #sigma_{#pi}^{TOF}", kTH2F, {{1000, 0, 10}, {100, -5, +5}}, false); - fRegistry.add("Track/hTOFNsigmaKa", "TOF n sigma ka;p_{in} (GeV/c);n #sigma_{K}^{TOF}", kTH2F, {{1000, 0, 10}, {100, -5, +5}}, false); - fRegistry.add("Track/hTOFNsigmaPr", "TOF n sigma pr;p_{in} (GeV/c);n #sigma_{p}^{TOF}", kTH2F, {{1000, 0, 10}, {100, -5, +5}}, false); - fRegistry.add("Track/hTPCNcr2Nf", "TPC Ncr/Nfindable", kTH1F, {{200, 0, 2}}, false); - fRegistry.add("Track/hTPCNcls2Nf", "TPC Ncls/Nfindable", kTH1F, {{200, 0, 2}}, false); - fRegistry.add("Track/hNclsITS", "number of ITS clusters", kTH1F, {{8, -0.5, 7.5}}, false); - fRegistry.add("Track/hChi2ITS", "chi2/number of ITS clusters", kTH1F, {{100, 0, 10}}, false); - fRegistry.add("Track/hITSClusterMap", "ITS cluster map", kTH1F, {{128, -0.5, 127.5}}, false); - fRegistry.add("Track/hMeanClusterSizeITS", "mean cluster size ITS; on ITS #times cos(#lambda)", kTH1F, {{32, 0, 16}}, false); - // event info if (nmod == 2) { o2::aod::pwgem::dilepton::utils::eventhistogram::addEventHistograms<2>(&fRegistry); @@ -422,38 +390,38 @@ struct dielectronQC { o2::ml::OnnxModel* eid_bdt = nullptr; void DefineDileptonCut() { - fDielectonCut = DielectronCut("fDielectonCut", "fDielectonCut"); + fDielectronCut = DielectronCut("fDielectronCut", "fDielectronCut"); // for pair - fDielectonCut.SetMeeRange(dielectroncuts.cfg_min_mass, dielectroncuts.cfg_max_mass); - fDielectonCut.SetMaxPhivPairMeeDep([&](float mll) { return (mll - dielectroncuts.cfg_phiv_intercept) / dielectroncuts.cfg_phiv_slope; }); - fDielectonCut.SetPairDCARange(dielectroncuts.cfg_min_pair_dca3d, dielectroncuts.cfg_max_pair_dca3d); // in sigma - fDielectonCut.ApplyPhiV(dielectroncuts.cfg_apply_phiv); - fDielectonCut.ApplyPrefilter(dielectroncuts.cfg_apply_pf); - fDielectonCut.RequireITSibAny(dielectroncuts.cfg_require_itsib_any); - fDielectonCut.RequireITSib1st(dielectroncuts.cfg_require_itsib_1st); + fDielectronCut.SetMeeRange(dielectroncuts.cfg_min_mass, dielectroncuts.cfg_max_mass); + fDielectronCut.SetMaxPhivPairMeeDep([&](float mll) { return (mll - dielectroncuts.cfg_phiv_intercept) / dielectroncuts.cfg_phiv_slope; }); + fDielectronCut.SetPairDCARange(dielectroncuts.cfg_min_pair_dca3d, dielectroncuts.cfg_max_pair_dca3d); // in sigma + fDielectronCut.ApplyPhiV(dielectroncuts.cfg_apply_phiv); + fDielectronCut.ApplyPrefilter(dielectroncuts.cfg_apply_pf); + fDielectronCut.RequireITSibAny(dielectroncuts.cfg_require_itsib_any); + fDielectronCut.RequireITSib1st(dielectroncuts.cfg_require_itsib_1st); // for track - fDielectonCut.SetTrackPtRange(dielectroncuts.cfg_min_pt_track, 1e+10f); - fDielectonCut.SetTrackEtaRange(-dielectroncuts.cfg_max_eta_track, +dielectroncuts.cfg_max_eta_track); - fDielectonCut.SetMinNClustersTPC(dielectroncuts.cfg_min_ncluster_tpc); - fDielectonCut.SetMinNCrossedRowsTPC(dielectroncuts.cfg_min_ncrossedrows); - fDielectonCut.SetMinNCrossedRowsOverFindableClustersTPC(0.8); - fDielectonCut.SetChi2PerClusterTPC(0.0, dielectroncuts.cfg_max_chi2tpc); - fDielectonCut.SetChi2PerClusterITS(0.0, dielectroncuts.cfg_max_chi2its); - fDielectonCut.SetNClustersITS(dielectroncuts.cfg_min_ncluster_its, 7); - fDielectonCut.SetMeanClusterSizeITSob(0, 16); - fDielectonCut.SetMaxDcaXY(dielectroncuts.cfg_max_dcaxy); - fDielectonCut.SetMaxDcaZ(dielectroncuts.cfg_max_dcaz); + fDielectronCut.SetTrackPtRange(dielectroncuts.cfg_min_pt_track, 1e+10f); + fDielectronCut.SetTrackEtaRange(-dielectroncuts.cfg_max_eta_track, +dielectroncuts.cfg_max_eta_track); + fDielectronCut.SetMinNClustersTPC(dielectroncuts.cfg_min_ncluster_tpc); + fDielectronCut.SetMinNCrossedRowsTPC(dielectroncuts.cfg_min_ncrossedrows); + fDielectronCut.SetMinNCrossedRowsOverFindableClustersTPC(0.8); + fDielectronCut.SetChi2PerClusterTPC(0.0, dielectroncuts.cfg_max_chi2tpc); + fDielectronCut.SetChi2PerClusterITS(0.0, dielectroncuts.cfg_max_chi2its); + fDielectronCut.SetNClustersITS(dielectroncuts.cfg_min_ncluster_its, 7); + fDielectronCut.SetMeanClusterSizeITSob(0, 16); + fDielectronCut.SetMaxDcaXY(dielectroncuts.cfg_max_dcaxy); + fDielectronCut.SetMaxDcaZ(dielectroncuts.cfg_max_dcaz); // for eID - fDielectonCut.SetPIDScheme(dielectroncuts.cfg_pid_scheme); - fDielectonCut.SetTPCNsigmaElRange(dielectroncuts.cfg_min_TPCNsigmaEl, dielectroncuts.cfg_max_TPCNsigmaEl); - fDielectonCut.SetTPCNsigmaMuRange(dielectroncuts.cfg_min_TPCNsigmaMu, dielectroncuts.cfg_max_TPCNsigmaMu); - fDielectonCut.SetTPCNsigmaPiRange(dielectroncuts.cfg_min_TPCNsigmaPi, dielectroncuts.cfg_max_TPCNsigmaPi); - fDielectonCut.SetTPCNsigmaKaRange(dielectroncuts.cfg_min_TPCNsigmaKa, dielectroncuts.cfg_max_TPCNsigmaKa); - fDielectonCut.SetTPCNsigmaPrRange(dielectroncuts.cfg_min_TPCNsigmaPr, dielectroncuts.cfg_max_TPCNsigmaPr); - fDielectonCut.SetTOFNsigmaElRange(dielectroncuts.cfg_min_TOFNsigmaEl, dielectroncuts.cfg_max_TOFNsigmaEl); + fDielectronCut.SetPIDScheme(dielectroncuts.cfg_pid_scheme); + fDielectronCut.SetTPCNsigmaElRange(dielectroncuts.cfg_min_TPCNsigmaEl, dielectroncuts.cfg_max_TPCNsigmaEl); + fDielectronCut.SetTPCNsigmaMuRange(dielectroncuts.cfg_min_TPCNsigmaMu, dielectroncuts.cfg_max_TPCNsigmaMu); + fDielectronCut.SetTPCNsigmaPiRange(dielectroncuts.cfg_min_TPCNsigmaPi, dielectroncuts.cfg_max_TPCNsigmaPi); + fDielectronCut.SetTPCNsigmaKaRange(dielectroncuts.cfg_min_TPCNsigmaKa, dielectroncuts.cfg_max_TPCNsigmaKa); + fDielectronCut.SetTPCNsigmaPrRange(dielectroncuts.cfg_min_TPCNsigmaPr, dielectroncuts.cfg_max_TPCNsigmaPr); + fDielectronCut.SetTOFNsigmaElRange(dielectroncuts.cfg_min_TOFNsigmaEl, dielectroncuts.cfg_max_TOFNsigmaEl); if (dielectroncuts.cfg_pid_scheme == static_cast(DielectronCut::PIDSchemes::kPIDML)) { // please call this at the end of DefineDileptonCut eid_bdt = new o2::ml::OnnxModel(); @@ -470,7 +438,7 @@ struct dielectronQC { eid_bdt->initModel(dielectroncuts.BDTLocalPathGamma.value, dielectroncuts.enableOptimizations.value); } - fDielectonCut.SetPIDModel(eid_bdt); + fDielectronCut.SetPIDModel(eid_bdt); } // end of PID ML } @@ -505,17 +473,17 @@ struct dielectronQC { if constexpr (ev_id == 0) { if (dielectroncuts.cfg_pid_scheme == static_cast(DielectronCut::PIDSchemes::kPIDML)) { - if (!fDielectonCut.IsSelectedTrack(t1, collision) || !fDielectonCut.IsSelectedTrack(t2, collision)) { + if (!fDielectronCut.IsSelectedTrack(t1, collision) || !fDielectronCut.IsSelectedTrack(t2, collision)) { return false; } } else { // cut-based - if (!fDielectonCut.IsSelectedTrack(t1) || !fDielectonCut.IsSelectedTrack(t2)) { + if (!fDielectronCut.IsSelectedTrack(t1) || !fDielectronCut.IsSelectedTrack(t2)) { return false; } } } - if (!fDielectonCut.IsSelectedPair(t1, t2, d_bz)) { + if (!fDielectronCut.IsSelectedPair(t1, t2, d_bz)) { return false; } @@ -719,7 +687,6 @@ struct dielectronQC { if (std::find(used_trackIds.begin(), used_trackIds.end(), pair_tmp_id1) == used_trackIds.end()) { used_trackIds.emplace_back(pair_tmp_id1); - fillTrackInfo(t1); if (cfgDoMix) { if (t1.sign() > 0) { emh_pos->AddTrackToEventPool(key_df_collision, EMTrackWithCov(t1.globalIndex(), collision.globalIndex(), t1.trackId(), t1.pt(), t1.eta(), t1.phi(), o2::constants::physics::MassElectron, t1.sign(), t1.dcaXY(), t1.dcaZ(), possibleIds1, @@ -734,7 +701,6 @@ struct dielectronQC { } if (std::find(used_trackIds.begin(), used_trackIds.end(), pair_tmp_id2) == used_trackIds.end()) { used_trackIds.emplace_back(pair_tmp_id2); - fillTrackInfo(t2); if (cfgDoMix) { if (t2.sign() > 0) { emh_pos->AddTrackToEventPool(key_df_collision, EMTrackWithCov(t2.globalIndex(), collision.globalIndex(), t2.trackId(), t2.pt(), t2.eta(), t2.phi(), o2::constants::physics::MassElectron, t2.sign(), t2.dcaXY(), t2.dcaZ(), possibleIds1, @@ -751,42 +717,6 @@ struct dielectronQC { return true; } - template - void fillTrackInfo(TTrack const& track) - { - float dca_3d = dca3DinSigma(track); - fRegistry.fill(HIST("Track/hPt"), track.pt()); - fRegistry.fill(HIST("Track/hQoverPt"), track.sign() / track.pt()); - fRegistry.fill(HIST("Track/hEtaPhi"), track.phi(), track.eta()); - fRegistry.fill(HIST("Track/hDCAxyz"), track.dcaXY(), track.dcaZ()); - fRegistry.fill(HIST("Track/hDCAxyzSigma"), track.dcaXY() / sqrt(track.cYY()), track.dcaZ() / sqrt(track.cZZ())); - fRegistry.fill(HIST("Track/hDCA3DSigma"), dca_3d); - fRegistry.fill(HIST("Track/hDCAxyRes_Pt"), track.pt(), sqrt(track.cYY()) * 1e+4); // convert cm to um - fRegistry.fill(HIST("Track/hDCAzRes_Pt"), track.pt(), sqrt(track.cZZ()) * 1e+4); // convert cm to um - fRegistry.fill(HIST("Track/hNclsITS"), track.itsNCls()); - fRegistry.fill(HIST("Track/hNclsTPC"), track.tpcNClsFound()); - fRegistry.fill(HIST("Track/hNcrTPC"), track.tpcNClsCrossedRows()); - fRegistry.fill(HIST("Track/hTPCNcr2Nf"), track.tpcCrossedRowsOverFindableCls()); - fRegistry.fill(HIST("Track/hTPCNcls2Nf"), track.tpcFoundOverFindableCls()); - fRegistry.fill(HIST("Track/hChi2TPC"), track.tpcChi2NCl()); - fRegistry.fill(HIST("Track/hChi2ITS"), track.itsChi2NCl()); - fRegistry.fill(HIST("Track/hITSClusterMap"), track.itsClusterMap()); - fRegistry.fill(HIST("Track/hMeanClusterSizeITS"), track.meanClusterSizeITS() * std::cos(std::atan(track.tgl()))); - fRegistry.fill(HIST("Track/hTPCdEdx"), track.tpcInnerParam(), track.tpcSignal()); - fRegistry.fill(HIST("Track/hTPCNsigmaEl"), track.tpcInnerParam(), track.tpcNSigmaEl()); - fRegistry.fill(HIST("Track/hTPCNsigmaMu"), track.tpcInnerParam(), track.tpcNSigmaMu()); - fRegistry.fill(HIST("Track/hTPCNsigmaPi"), track.tpcInnerParam(), track.tpcNSigmaPi()); - fRegistry.fill(HIST("Track/hTPCNsigmaKa"), track.tpcInnerParam(), track.tpcNSigmaKa()); - fRegistry.fill(HIST("Track/hTPCNsigmaPr"), track.tpcInnerParam(), track.tpcNSigmaPr()); - fRegistry.fill(HIST("Track/hTOFbeta"), track.tpcInnerParam(), track.beta()); - fRegistry.fill(HIST("Track/h1overTOFbeta"), track.tpcInnerParam(), 1. / track.beta()); - fRegistry.fill(HIST("Track/hTOFNsigmaEl"), track.tpcInnerParam(), track.tofNSigmaEl()); - fRegistry.fill(HIST("Track/hTOFNsigmaMu"), track.tpcInnerParam(), track.tofNSigmaMu()); - fRegistry.fill(HIST("Track/hTOFNsigmaPi"), track.tpcInnerParam(), track.tofNSigmaPi()); - fRegistry.fill(HIST("Track/hTOFNsigmaKa"), track.tpcInnerParam(), track.tofNSigmaKa()); - fRegistry.fill(HIST("Track/hTOFNsigmaPr"), track.tpcInnerParam(), track.tofNSigmaPr()); - } - Filter collisionFilter_centrality = (cfgCentMin < o2::aod::cent::centFT0M && o2::aod::cent::centFT0M < cfgCentMax) || (cfgCentMin < o2::aod::cent::centFT0A && o2::aod::cent::centFT0A < cfgCentMax) || (cfgCentMin < o2::aod::cent::centFT0C && o2::aod::cent::centFT0C < cfgCentMax); using FilteredMyCollisions = soa::Filtered; diff --git a/PWGEM/Dilepton/Tasks/dielectronQCMC.cxx b/PWGEM/Dilepton/Tasks/dielectronQCMC.cxx index 0104ffcbc53..0dc49fd8f86 100644 --- a/PWGEM/Dilepton/Tasks/dielectronQCMC.cxx +++ b/PWGEM/Dilepton/Tasks/dielectronQCMC.cxx @@ -50,7 +50,7 @@ using namespace o2::soa; using namespace o2::aod::pwgem::dilepton::utils::mcutil; using namespace o2::aod::pwgem::dilepton::utils::emtrackutil; -using MyCollisions = soa::Join; +using MyCollisions = soa::Join; using MyCollision = MyCollisions::iterator; using MyMCTracks = soa::Join; @@ -91,7 +91,7 @@ struct dielectronQCMC { Configurable cfgOccupancyMax{"cfgOccupancyMax", 1000000000, "max. occupancy"}; } eventcuts; - DielectronCut fDielectonCut; + DielectronCut fDielectronCut; struct : ConfigurableGroup { std::string prefix = "dielectroncut_group"; Configurable cfg_min_mass{"cfg_min_mass", 0.0, "min mass"}; @@ -155,7 +155,6 @@ struct dielectronQCMC { HistogramRegistry fRegistry{"output", {}, OutputObjHandlingPolicy::AnalysisObject, false, false}; static constexpr std::string_view event_cut_types[2] = {"before/", "after/"}; - static constexpr std::string_view ele_source_types[9] = {"lf/", "Photon/", "PromptJPsi/", "NonPromptJPsi/", "PromptPsi2S/", "NonPromptPsi2S/", "c2e/", "b2e/", "b2c2e/"}; ~dielectronQCMC() { @@ -241,49 +240,6 @@ struct dielectronQCMC { fRegistry.addClone("Pair/corr_bkg_eh/uls/", "Pair/corr_bkg_eh/lsmm/"); fRegistry.addClone("Pair/corr_bkg_eh/", "Pair/corr_bkg_hh/"); fRegistry.addClone("Pair/corr_bkg_eh/", "Pair/comb_bkg/"); - - // track info - fRegistry.add("Track/lf/hPt", "pT;p_{T} (GeV/c)", kTH1F, {{1000, 0.0f, 10}}, false); - fRegistry.add("Track/lf/hQoverPt", "q/pT;q/p_{T} (GeV/c)^{-1}", kTH1F, {{400, -20, 20}}, false); - fRegistry.add("Track/lf/hEtaPhi", "#eta vs. #varphi;#varphi (rad.);#eta", kTH2F, {{180, 0, 2 * M_PI}, {40, -2.0f, 2.0f}}, false); - fRegistry.add("Track/lf/hDCAxyz", "DCA xy vs. z;DCA_{xy} (cm);DCA_{z} (cm)", kTH2F, {{200, -1.0f, 1.0f}, {200, -1.0f, 1.0f}}, false); - fRegistry.add("Track/lf/hDCAxyzSigma", "DCA xy vs. z;DCA_{xy} (#sigma);DCA_{z} (#sigma)", kTH2F, {{200, -10.0f, 10.0f}, {200, -10.0f, 10.0f}}, false); - fRegistry.add("Track/lf/hDCA3DSigma", "DCA 3D;DCA_{3D} (#sigma);", kTH1F, {{100, 0.0f, 10.0f}}, false); - fRegistry.add("Track/lf/hDCAxyRes_Pt", "DCA_{xy} resolution vs. pT;p_{T} (GeV/c);DCA_{xy} resolution (#mum)", kTH2F, {{1000, 0, 10}, {500, 0., 500}}, false); - fRegistry.add("Track/lf/hDCAzRes_Pt", "DCA_{z} resolution vs. pT;p_{T} (GeV/c);DCA_{z} resolution (#mum)", kTH2F, {{1000, 0, 10}, {500, 0., 500}}, false); - fRegistry.add("Track/lf/hNclsTPC", "number of TPC clusters", kTH1F, {{161, -0.5, 160.5}}, false); - fRegistry.add("Track/lf/hNcrTPC", "number of TPC crossed rows", kTH1F, {{161, -0.5, 160.5}}, false); - fRegistry.add("Track/lf/hChi2TPC", "chi2/number of TPC clusters", kTH1F, {{100, 0, 10}}, false); - fRegistry.add("Track/lf/hTPCdEdx", "TPC dE/dx;p_{in} (GeV/c);TPC dE/dx (a.u.)", kTH2F, {{1000, 0, 10}, {200, 0, 200}}, false); - fRegistry.add("Track/lf/hTPCNsigmaEl", "TPC n sigma el;p_{in} (GeV/c);n #sigma_{e}^{TPC}", kTH2F, {{1000, 0, 10}, {100, -5, +5}}, false); - fRegistry.add("Track/lf/hTPCNsigmaMu", "TPC n sigma mu;p_{in} (GeV/c);n #sigma_{#mu}^{TPC}", kTH2F, {{1000, 0, 10}, {100, -5, +5}}, false); - fRegistry.add("Track/lf/hTPCNsigmaPi", "TPC n sigma pi;p_{in} (GeV/c);n #sigma_{#pi}^{TPC}", kTH2F, {{1000, 0, 10}, {100, -5, +5}}, false); - fRegistry.add("Track/lf/hTPCNsigmaKa", "TPC n sigma ka;p_{in} (GeV/c);n #sigma_{K}^{TPC}", kTH2F, {{1000, 0, 10}, {100, -5, +5}}, false); - fRegistry.add("Track/lf/hTPCNsigmaPr", "TPC n sigma pr;p_{in} (GeV/c);n #sigma_{p}^{TPC}", kTH2F, {{1000, 0, 10}, {100, -5, +5}}, false); - fRegistry.add("Track/lf/hTOFbeta", "TOF #beta;p_{in} (GeV/c);#beta", kTH2F, {{1000, 0, 10}, {600, 0, 1.2}}, false); - fRegistry.add("Track/lf/h1overTOFbeta", "TOF 1/#beta;p_{in} (GeV/c);1/#beta", kTH2F, {{1000, 0, 10}, {1000, 0.8, 1.8}}, false); - fRegistry.add("Track/lf/hTOFNsigmaEl", "TOF n sigma el;p_{in} (GeV/c);n #sigma_{e}^{TOF}", kTH2F, {{1000, 0, 10}, {100, -5, +5}}, false); - fRegistry.add("Track/lf/hTOFNsigmaMu", "TOF n sigma mu;p_{in} (GeV/c);n #sigma_{#mu}^{TOF}", kTH2F, {{1000, 0, 10}, {100, -5, +5}}, false); - fRegistry.add("Track/lf/hTOFNsigmaPi", "TOF n sigma pi;p_{in} (GeV/c);n #sigma_{#pi}^{TOF}", kTH2F, {{1000, 0, 10}, {100, -5, +5}}, false); - fRegistry.add("Track/lf/hTOFNsigmaKa", "TOF n sigma ka;p_{in} (GeV/c);n #sigma_{K}^{TOF}", kTH2F, {{1000, 0, 10}, {100, -5, +5}}, false); - fRegistry.add("Track/lf/hTOFNsigmaPr", "TOF n sigma pr;p_{in} (GeV/c);n #sigma_{p}^{TOF}", kTH2F, {{1000, 0, 10}, {100, -5, +5}}, false); - fRegistry.add("Track/lf/hTPCNcr2Nf", "TPC Ncr/Nfindable", kTH1F, {{200, 0, 2}}, false); - fRegistry.add("Track/lf/hTPCNcls2Nf", "TPC Ncls/Nfindable", kTH1F, {{200, 0, 2}}, false); - fRegistry.add("Track/lf/hNclsITS", "number of ITS clusters", kTH1F, {{8, -0.5, 7.5}}, false); - fRegistry.add("Track/lf/hChi2ITS", "chi2/number of ITS clusters", kTH1F, {{100, 0, 10}}, false); - fRegistry.add("Track/lf/hITSClusterMap", "ITS cluster map", kTH1F, {{128, -0.5, 127.5}}, false); - fRegistry.add("Track/lf/hMeanClusterSizeITS", "mean cluster size ITS; on ITS #times cos(#lambda)", kTH1F, {{32, 0, 16}}, false); - fRegistry.add("Track/lf/hPtGen_DeltaPtOverPtGen", "electron p_{T} resolution;p_{T}^{gen} (GeV/c);(p_{T}^{rec} - p_{T}^{gen})/p_{T}^{gen}", kTH2F, {{1000, 0, 10}, {400, -1.0f, 1.0f}}, true); - fRegistry.add("Track/lf/hPtGen_DeltaEta", "electron #eta resolution;p_{T}^{gen} (GeV/c);#eta^{rec} - #eta^{gen}", kTH2F, {{1000, 0, 10}, {400, -1.0f, 1.0f}}, true); - fRegistry.add("Track/lf/hPtGen_DeltaPhi", "electron #varphi resolution;p_{T}^{gen} (GeV/c);#varphi^{rec} - #varphi^{gen} (rad.)", kTH2F, {{1000, 0, 10}, {400, -1.0f, 1.0f}}, true); - fRegistry.addClone("Track/lf/", "Track/Photon/"); - fRegistry.addClone("Track/lf/", "Track/PromptJPsi/"); - fRegistry.addClone("Track/lf/", "Track/NonPromptJPsi/"); - fRegistry.addClone("Track/lf/", "Track/PromptPsi2S/"); - fRegistry.addClone("Track/lf/", "Track/NonPromptPsi2S/"); - fRegistry.addClone("Track/lf/", "Track/c2e/"); - fRegistry.addClone("Track/lf/", "Track/b2e/"); - fRegistry.addClone("Track/lf/", "Track/b2c2e/"); } float beamM1 = o2::constants::physics::MassProton; // mass of beam @@ -295,10 +251,6 @@ struct dielectronQCMC { void init(InitContext&) { - DefineEMEventCut(); - DefineDileptonCut(); - addhistograms(); - mRunNumber = 0; d_bz = 0; @@ -307,6 +259,10 @@ struct dielectronQCMC { ccdb->setLocalObjectValidityChecking(); ccdb->setFatalWhenNull(false); + DefineEMEventCut(); + DefineDileptonCut(); + addhistograms(); + // fitter.setPropagateToPCA(true); // fitter.setMaxR(5.f); // fitter.setMinParamChange(1e-3); @@ -389,38 +345,38 @@ struct dielectronQCMC { o2::ml::OnnxModel* eid_bdt = nullptr; void DefineDileptonCut() { - fDielectonCut = DielectronCut("fDielectonCut", "fDielectonCut"); + fDielectronCut = DielectronCut("fDielectronCut", "fDielectronCut"); // for pair - fDielectonCut.SetMeeRange(dielectroncuts.cfg_min_mass, dielectroncuts.cfg_max_mass); - fDielectonCut.SetMaxPhivPairMeeDep([&](float mll) { return (mll - dielectroncuts.cfg_phiv_intercept) / dielectroncuts.cfg_phiv_slope; }); - fDielectonCut.SetPairDCARange(dielectroncuts.cfg_min_pair_dca3d, dielectroncuts.cfg_max_pair_dca3d); // in sigma - fDielectonCut.ApplyPhiV(dielectroncuts.cfg_apply_phiv); - fDielectonCut.ApplyPrefilter(dielectroncuts.cfg_apply_pf); - fDielectonCut.RequireITSibAny(dielectroncuts.cfg_require_itsib_any); - fDielectonCut.RequireITSib1st(dielectroncuts.cfg_require_itsib_1st); + fDielectronCut.SetMeeRange(dielectroncuts.cfg_min_mass, dielectroncuts.cfg_max_mass); + fDielectronCut.SetMaxPhivPairMeeDep([&](float mll) { return (mll - dielectroncuts.cfg_phiv_intercept) / dielectroncuts.cfg_phiv_slope; }); + fDielectronCut.SetPairDCARange(dielectroncuts.cfg_min_pair_dca3d, dielectroncuts.cfg_max_pair_dca3d); // in sigma + fDielectronCut.ApplyPhiV(dielectroncuts.cfg_apply_phiv); + fDielectronCut.ApplyPrefilter(dielectroncuts.cfg_apply_pf); + fDielectronCut.RequireITSibAny(dielectroncuts.cfg_require_itsib_any); + fDielectronCut.RequireITSib1st(dielectroncuts.cfg_require_itsib_1st); // for track - fDielectonCut.SetTrackPtRange(dielectroncuts.cfg_min_pt_track, 1e+10f); - fDielectonCut.SetTrackEtaRange(-dielectroncuts.cfg_max_eta_track, +dielectroncuts.cfg_max_eta_track); - fDielectonCut.SetMinNClustersTPC(dielectroncuts.cfg_min_ncluster_tpc); - fDielectonCut.SetMinNCrossedRowsTPC(dielectroncuts.cfg_min_ncrossedrows); - fDielectonCut.SetMinNCrossedRowsOverFindableClustersTPC(0.8); - fDielectonCut.SetChi2PerClusterTPC(0.0, dielectroncuts.cfg_max_chi2tpc); - fDielectonCut.SetChi2PerClusterITS(0.0, dielectroncuts.cfg_max_chi2its); - fDielectonCut.SetNClustersITS(dielectroncuts.cfg_min_ncluster_its, 7); - fDielectonCut.SetMeanClusterSizeITSob(0, 16); - fDielectonCut.SetMaxDcaXY(dielectroncuts.cfg_max_dcaxy); - fDielectonCut.SetMaxDcaZ(dielectroncuts.cfg_max_dcaz); + fDielectronCut.SetTrackPtRange(dielectroncuts.cfg_min_pt_track, 1e+10f); + fDielectronCut.SetTrackEtaRange(-dielectroncuts.cfg_max_eta_track, +dielectroncuts.cfg_max_eta_track); + fDielectronCut.SetMinNClustersTPC(dielectroncuts.cfg_min_ncluster_tpc); + fDielectronCut.SetMinNCrossedRowsTPC(dielectroncuts.cfg_min_ncrossedrows); + fDielectronCut.SetMinNCrossedRowsOverFindableClustersTPC(0.8); + fDielectronCut.SetChi2PerClusterTPC(0.0, dielectroncuts.cfg_max_chi2tpc); + fDielectronCut.SetChi2PerClusterITS(0.0, dielectroncuts.cfg_max_chi2its); + fDielectronCut.SetNClustersITS(dielectroncuts.cfg_min_ncluster_its, 7); + fDielectronCut.SetMeanClusterSizeITSob(0, 16); + fDielectronCut.SetMaxDcaXY(dielectroncuts.cfg_max_dcaxy); + fDielectronCut.SetMaxDcaZ(dielectroncuts.cfg_max_dcaz); // for eID - fDielectonCut.SetPIDScheme(dielectroncuts.cfg_pid_scheme); - fDielectonCut.SetTPCNsigmaElRange(dielectroncuts.cfg_min_TPCNsigmaEl, dielectroncuts.cfg_max_TPCNsigmaEl); - fDielectonCut.SetTPCNsigmaMuRange(dielectroncuts.cfg_min_TPCNsigmaMu, dielectroncuts.cfg_max_TPCNsigmaMu); - fDielectonCut.SetTPCNsigmaPiRange(dielectroncuts.cfg_min_TPCNsigmaPi, dielectroncuts.cfg_max_TPCNsigmaPi); - fDielectonCut.SetTPCNsigmaKaRange(dielectroncuts.cfg_min_TPCNsigmaKa, dielectroncuts.cfg_max_TPCNsigmaKa); - fDielectonCut.SetTPCNsigmaPrRange(dielectroncuts.cfg_min_TPCNsigmaPr, dielectroncuts.cfg_max_TPCNsigmaPr); - fDielectonCut.SetTOFNsigmaElRange(dielectroncuts.cfg_min_TOFNsigmaEl, dielectroncuts.cfg_max_TOFNsigmaEl); + fDielectronCut.SetPIDScheme(dielectroncuts.cfg_pid_scheme); + fDielectronCut.SetTPCNsigmaElRange(dielectroncuts.cfg_min_TPCNsigmaEl, dielectroncuts.cfg_max_TPCNsigmaEl); + fDielectronCut.SetTPCNsigmaMuRange(dielectroncuts.cfg_min_TPCNsigmaMu, dielectroncuts.cfg_max_TPCNsigmaMu); + fDielectronCut.SetTPCNsigmaPiRange(dielectroncuts.cfg_min_TPCNsigmaPi, dielectroncuts.cfg_max_TPCNsigmaPi); + fDielectronCut.SetTPCNsigmaKaRange(dielectroncuts.cfg_min_TPCNsigmaKa, dielectroncuts.cfg_max_TPCNsigmaKa); + fDielectronCut.SetTPCNsigmaPrRange(dielectroncuts.cfg_min_TPCNsigmaPr, dielectroncuts.cfg_max_TPCNsigmaPr); + fDielectronCut.SetTOFNsigmaElRange(dielectroncuts.cfg_min_TOFNsigmaEl, dielectroncuts.cfg_max_TOFNsigmaEl); if (dielectroncuts.cfg_pid_scheme == static_cast(DielectronCut::PIDSchemes::kPIDML)) { // please call this at the end of DefineDileptonCut // o2::ml::OnnxModel* eid_bdt = new o2::ml::OnnxModel(); @@ -438,7 +394,7 @@ struct dielectronQCMC { eid_bdt->initModel(dielectroncuts.BDTLocalPathGamma.value, dielectroncuts.enableOptimizations.value); } - fDielectonCut.SetPIDModel(eid_bdt); + fDielectronCut.SetPIDModel(eid_bdt); } // end of PID ML } @@ -474,64 +430,20 @@ struct dielectronQCMC { } } - template - void fillTrackInfo(TTrack const& track) - { - // fill track info that belong to true pairs. - if (std::find(used_trackIds.begin(), used_trackIds.end(), track.globalIndex()) == used_trackIds.end()) { - auto mctrack = track.template emmcparticle_as(); - float dca_3d = dca3DinSigma(track); - fRegistry.fill(HIST("Track/") + HIST(ele_source_types[e_source_id]) + HIST("hPt"), track.pt()); - fRegistry.fill(HIST("Track/") + HIST(ele_source_types[e_source_id]) + HIST("hQoverPt"), track.sign() / track.pt()); - fRegistry.fill(HIST("Track/") + HIST(ele_source_types[e_source_id]) + HIST("hEtaPhi"), track.phi(), track.eta()); - fRegistry.fill(HIST("Track/") + HIST(ele_source_types[e_source_id]) + HIST("hDCAxyz"), track.dcaXY(), track.dcaZ()); - fRegistry.fill(HIST("Track/") + HIST(ele_source_types[e_source_id]) + HIST("hDCAxyzSigma"), track.dcaXY() / sqrt(track.cYY()), track.dcaZ() / sqrt(track.cZZ())); - fRegistry.fill(HIST("Track/") + HIST(ele_source_types[e_source_id]) + HIST("hDCA3DSigma"), dca_3d); - fRegistry.fill(HIST("Track/") + HIST(ele_source_types[e_source_id]) + HIST("hDCAxyRes_Pt"), track.pt(), sqrt(track.cYY()) * 1e+4); // convert cm to um - fRegistry.fill(HIST("Track/") + HIST(ele_source_types[e_source_id]) + HIST("hDCAzRes_Pt"), track.pt(), sqrt(track.cZZ()) * 1e+4); // convert cm to um - fRegistry.fill(HIST("Track/") + HIST(ele_source_types[e_source_id]) + HIST("hNclsITS"), track.itsNCls()); - fRegistry.fill(HIST("Track/") + HIST(ele_source_types[e_source_id]) + HIST("hNclsTPC"), track.tpcNClsFound()); - fRegistry.fill(HIST("Track/") + HIST(ele_source_types[e_source_id]) + HIST("hNcrTPC"), track.tpcNClsCrossedRows()); - fRegistry.fill(HIST("Track/") + HIST(ele_source_types[e_source_id]) + HIST("hTPCNcr2Nf"), track.tpcCrossedRowsOverFindableCls()); - fRegistry.fill(HIST("Track/") + HIST(ele_source_types[e_source_id]) + HIST("hTPCNcls2Nf"), track.tpcFoundOverFindableCls()); - fRegistry.fill(HIST("Track/") + HIST(ele_source_types[e_source_id]) + HIST("hChi2TPC"), track.tpcChi2NCl()); - fRegistry.fill(HIST("Track/") + HIST(ele_source_types[e_source_id]) + HIST("hChi2ITS"), track.itsChi2NCl()); - fRegistry.fill(HIST("Track/") + HIST(ele_source_types[e_source_id]) + HIST("hITSClusterMap"), track.itsClusterMap()); - fRegistry.fill(HIST("Track/") + HIST(ele_source_types[e_source_id]) + HIST("hMeanClusterSizeITS"), track.meanClusterSizeITS() * std::cos(std::atan(track.tgl()))); - fRegistry.fill(HIST("Track/") + HIST(ele_source_types[e_source_id]) + HIST("hTPCdEdx"), track.tpcInnerParam(), track.tpcSignal()); - fRegistry.fill(HIST("Track/") + HIST(ele_source_types[e_source_id]) + HIST("hTPCNsigmaEl"), track.tpcInnerParam(), track.tpcNSigmaEl()); - fRegistry.fill(HIST("Track/") + HIST(ele_source_types[e_source_id]) + HIST("hTPCNsigmaMu"), track.tpcInnerParam(), track.tpcNSigmaMu()); - fRegistry.fill(HIST("Track/") + HIST(ele_source_types[e_source_id]) + HIST("hTPCNsigmaPi"), track.tpcInnerParam(), track.tpcNSigmaPi()); - fRegistry.fill(HIST("Track/") + HIST(ele_source_types[e_source_id]) + HIST("hTPCNsigmaKa"), track.tpcInnerParam(), track.tpcNSigmaKa()); - fRegistry.fill(HIST("Track/") + HIST(ele_source_types[e_source_id]) + HIST("hTPCNsigmaPr"), track.tpcInnerParam(), track.tpcNSigmaPr()); - fRegistry.fill(HIST("Track/") + HIST(ele_source_types[e_source_id]) + HIST("hTOFbeta"), track.tpcInnerParam(), track.beta()); - fRegistry.fill(HIST("Track/") + HIST(ele_source_types[e_source_id]) + HIST("h1overTOFbeta"), track.tpcInnerParam(), 1. / track.beta()); - fRegistry.fill(HIST("Track/") + HIST(ele_source_types[e_source_id]) + HIST("hTOFNsigmaEl"), track.tpcInnerParam(), track.tofNSigmaEl()); - fRegistry.fill(HIST("Track/") + HIST(ele_source_types[e_source_id]) + HIST("hTOFNsigmaMu"), track.tpcInnerParam(), track.tofNSigmaMu()); - fRegistry.fill(HIST("Track/") + HIST(ele_source_types[e_source_id]) + HIST("hTOFNsigmaPi"), track.tpcInnerParam(), track.tofNSigmaPi()); - fRegistry.fill(HIST("Track/") + HIST(ele_source_types[e_source_id]) + HIST("hTOFNsigmaKa"), track.tpcInnerParam(), track.tofNSigmaKa()); - fRegistry.fill(HIST("Track/") + HIST(ele_source_types[e_source_id]) + HIST("hTOFNsigmaPr"), track.tpcInnerParam(), track.tofNSigmaPr()); - fRegistry.fill(HIST("Track/") + HIST(ele_source_types[e_source_id]) + HIST("hPtGen_DeltaPtOverPtGen"), mctrack.pt(), (track.pt() - mctrack.pt()) / mctrack.pt()); - fRegistry.fill(HIST("Track/") + HIST(ele_source_types[e_source_id]) + HIST("hPtGen_DeltaEta"), mctrack.pt(), track.eta() - mctrack.eta()); - fRegistry.fill(HIST("Track/") + HIST(ele_source_types[e_source_id]) + HIST("hPtGen_DeltaPhi"), mctrack.pt(), track.phi() - mctrack.phi()); - used_trackIds.emplace_back(track.globalIndex()); - } - } - template bool fillTruePairInfo(TCollision const& collision, TTrack1 const& t1, TTrack2 const& t2, TMCParticles const& mcparticles) { if (dielectroncuts.cfg_pid_scheme == static_cast(DielectronCut::PIDSchemes::kPIDML)) { - if (!fDielectonCut.IsSelectedTrack(t1, collision) || !fDielectonCut.IsSelectedTrack(t2, collision)) { + if (!fDielectronCut.IsSelectedTrack(t1, collision) || !fDielectronCut.IsSelectedTrack(t2, collision)) { return false; } } else { // cut-based - if (!fDielectonCut.IsSelectedTrack(t1) || !fDielectonCut.IsSelectedTrack(t2)) { + if (!fDielectronCut.IsSelectedTrack(t1) || !fDielectronCut.IsSelectedTrack(t2)) { return false; } } - if (!fDielectonCut.IsSelectedPair(t1, t2, d_bz)) { + if (!fDielectronCut.IsSelectedPair(t1, t2, d_bz)) { return false; } @@ -615,32 +527,22 @@ struct dielectronQCMC { case 111: fRegistry.fill(HIST("Pair/sm/Pi0/hs"), v12.M(), v12.Pt(), abs(dphi), abs(cos_thetaCS), abs(phiCS), aco, asym, abs(dphi_e_ee), dca_ee_3d); fRegistry.fill(HIST("Pair/sm/Pi0/hMvsPhiV"), phiv, v12.M()); - fillTrackInfo<0, TMCParticles>(t1); - fillTrackInfo<0, TMCParticles>(t2); break; case 221: fRegistry.fill(HIST("Pair/sm/Eta/hs"), v12.M(), v12.Pt(), abs(dphi), abs(cos_thetaCS), abs(phiCS), aco, asym, abs(dphi_e_ee), dca_ee_3d); fRegistry.fill(HIST("Pair/sm/Eta/hMvsPhiV"), phiv, v12.M()); - fillTrackInfo<0, TMCParticles>(t1); - fillTrackInfo<0, TMCParticles>(t2); break; case 331: fRegistry.fill(HIST("Pair/sm/EtaPrime/hs"), v12.M(), v12.Pt(), abs(dphi), abs(cos_thetaCS), abs(phiCS), aco, asym, abs(dphi_e_ee), dca_ee_3d); fRegistry.fill(HIST("Pair/sm/EtaPrime/hMvsPhiV"), phiv, v12.M()); - fillTrackInfo<0, TMCParticles>(t1); - fillTrackInfo<0, TMCParticles>(t2); break; case 113: fRegistry.fill(HIST("Pair/sm/Rho/hs"), v12.M(), v12.Pt(), abs(dphi), abs(cos_thetaCS), abs(phiCS), aco, asym, abs(dphi_e_ee), dca_ee_3d); fRegistry.fill(HIST("Pair/sm/Rho/hMvsPhiV"), phiv, v12.M()); - fillTrackInfo<0, TMCParticles>(t1); - fillTrackInfo<0, TMCParticles>(t2); break; case 223: fRegistry.fill(HIST("Pair/sm/Omega/hs"), v12.M(), v12.Pt(), abs(dphi), abs(cos_thetaCS), abs(phiCS), aco, asym, abs(dphi_e_ee), dca_ee_3d); fRegistry.fill(HIST("Pair/sm/Omega/hMvsPhiV"), phiv, v12.M()); - fillTrackInfo<0, TMCParticles>(t1); - fillTrackInfo<0, TMCParticles>(t2); if (mcmother.daughtersIds().size() == 2) { // omeag->ee fRegistry.fill(HIST("Pair/sm/Omega2ee/hs"), v12.M(), v12.Pt(), abs(dphi), abs(cos_thetaCS), abs(phiCS), aco, asym, abs(dphi_e_ee), dca_ee_3d); fRegistry.fill(HIST("Pair/sm/Omega2ee/hMvsPhiV"), phiv, v12.M()); @@ -649,8 +551,6 @@ struct dielectronQCMC { case 333: fRegistry.fill(HIST("Pair/sm/Phi/hs"), v12.M(), v12.Pt(), abs(dphi), abs(cos_thetaCS), abs(phiCS), aco, asym, abs(dphi_e_ee), dca_ee_3d); fRegistry.fill(HIST("Pair/sm/Phi/hMvsPhiV"), phiv, v12.M()); - fillTrackInfo<0, TMCParticles>(t1); - fillTrackInfo<0, TMCParticles>(t2); if (mcmother.daughtersIds().size() == 2) { // phi->ee fRegistry.fill(HIST("Pair/sm/Phi2ee/hs"), v12.M(), v12.Pt(), abs(dphi), abs(cos_thetaCS), abs(phiCS), aco, asym, abs(dphi_e_ee), dca_ee_3d); fRegistry.fill(HIST("Pair/sm/Phi2ee/hMvsPhiV"), phiv, v12.M()); @@ -660,13 +560,9 @@ struct dielectronQCMC { if (IsFromBeauty(mcmother, mcparticles) > 0) { fRegistry.fill(HIST("Pair/sm/NonPromptJPsi/hs"), v12.M(), v12.Pt(), abs(dphi), abs(cos_thetaCS), abs(phiCS), aco, asym, abs(dphi_e_ee), dca_ee_3d); fRegistry.fill(HIST("Pair/sm/NonPromptJPsi/hMvsPhiV"), phiv, v12.M()); - fillTrackInfo<3, TMCParticles>(t1); - fillTrackInfo<3, TMCParticles>(t2); } else { fRegistry.fill(HIST("Pair/sm/PromptJPsi/hs"), v12.M(), v12.Pt(), abs(dphi), abs(cos_thetaCS), abs(phiCS), aco, asym, abs(dphi_e_ee), dca_ee_3d); fRegistry.fill(HIST("Pair/sm/PromptJPsi/hMvsPhiV"), phiv, v12.M()); - fillTrackInfo<2, TMCParticles>(t1); - fillTrackInfo<2, TMCParticles>(t2); } break; } @@ -674,13 +570,9 @@ struct dielectronQCMC { if (IsFromBeauty(mcmother, mcparticles) > 0) { fRegistry.fill(HIST("Pair/sm/NonPromptPsi2S/hs"), v12.M(), v12.Pt(), abs(dphi), abs(cos_thetaCS), abs(phiCS), aco, asym, abs(dphi_e_ee), dca_ee_3d); fRegistry.fill(HIST("Pair/sm/NonPromptPsi2S/hMvsPhiV"), phiv, v12.M()); - fillTrackInfo<5, TMCParticles>(t1); - fillTrackInfo<5, TMCParticles>(t2); } else { fRegistry.fill(HIST("Pair/sm/PromptPsi2S/hs"), v12.M(), v12.Pt(), abs(dphi), abs(cos_thetaCS), abs(phiCS), aco, asym, abs(dphi_e_ee), dca_ee_3d); fRegistry.fill(HIST("Pair/sm/PromptPsi2S/hMvsPhiV"), phiv, v12.M()); - fillTrackInfo<4, TMCParticles>(t1); - fillTrackInfo<4, TMCParticles>(t2); } break; } @@ -693,8 +585,6 @@ struct dielectronQCMC { case 22: fRegistry.fill(HIST("Pair/sm/Photon/hs"), v12.M(), v12.Pt(), abs(dphi), abs(cos_thetaCS), abs(phiCS), aco, asym, abs(dphi_e_ee), dca_ee_3d); fRegistry.fill(HIST("Pair/sm/Photon/hMvsPhiV"), phiv, v12.M()); - fillTrackInfo<1, TMCParticles>(t1); - fillTrackInfo<1, TMCParticles>(t2); break; default: break; @@ -711,16 +601,10 @@ struct dielectronQCMC { fRegistry.fill(HIST("Pair/ccbar/c2e_c2e/hadron_hadron/hs"), v12.M(), v12.Pt(), abs(dphi), abs(cos_thetaCS), abs(phiCS), aco, asym, abs(dphi_e_ee), dca_ee_3d); if (isCharmMeson(mp1) && isCharmMeson(mp2)) { fRegistry.fill(HIST("Pair/ccbar/c2e_c2e/meson_meson/hs"), v12.M(), v12.Pt(), abs(dphi), abs(cos_thetaCS), abs(phiCS), aco, asym, abs(dphi_e_ee), dca_ee_3d); - fillTrackInfo<6, TMCParticles>(t1); - fillTrackInfo<6, TMCParticles>(t2); } else if (isCharmBaryon(mp1) && isCharmBaryon(mp2)) { fRegistry.fill(HIST("Pair/ccbar/c2e_c2e/baryon_baryon/hs"), v12.M(), v12.Pt(), abs(dphi), abs(cos_thetaCS), abs(phiCS), aco, asym, abs(dphi_e_ee), dca_ee_3d); - fillTrackInfo<6, TMCParticles>(t1); - fillTrackInfo<6, TMCParticles>(t2); } else { fRegistry.fill(HIST("Pair/ccbar/c2e_c2e/meson_baryon/hs"), v12.M(), v12.Pt(), abs(dphi), abs(cos_thetaCS), abs(phiCS), aco, asym, abs(dphi_e_ee), dca_ee_3d); - fillTrackInfo<6, TMCParticles>(t1); - fillTrackInfo<6, TMCParticles>(t2); } break; } @@ -728,16 +612,10 @@ struct dielectronQCMC { fRegistry.fill(HIST("Pair/bbbar/b2e_b2e/hadron_hadron/hs"), v12.M(), v12.Pt(), abs(dphi), abs(cos_thetaCS), abs(phiCS), aco, asym, abs(dphi_e_ee), dca_ee_3d); if (isBeautyMeson(mp1) && isBeautyMeson(mp2)) { fRegistry.fill(HIST("Pair/bbbar/b2e_b2e/meson_meson/hs"), v12.M(), v12.Pt(), abs(dphi), abs(cos_thetaCS), abs(phiCS), aco, asym, abs(dphi_e_ee), dca_ee_3d); - fillTrackInfo<7, TMCParticles>(t1); - fillTrackInfo<7, TMCParticles>(t2); } else if (isBeautyBaryon(mp1) && isBeautyBaryon(mp2)) { fRegistry.fill(HIST("Pair/bbbar/b2e_b2e/baryon_baryon/hs"), v12.M(), v12.Pt(), abs(dphi), abs(cos_thetaCS), abs(phiCS), aco, asym, abs(dphi_e_ee), dca_ee_3d); - fillTrackInfo<7, TMCParticles>(t1); - fillTrackInfo<7, TMCParticles>(t2); } else { fRegistry.fill(HIST("Pair/bbbar/b2e_b2e/meson_baryon/hs"), v12.M(), v12.Pt(), abs(dphi), abs(cos_thetaCS), abs(phiCS), aco, asym, abs(dphi_e_ee), dca_ee_3d); - fillTrackInfo<7, TMCParticles>(t1); - fillTrackInfo<7, TMCParticles>(t2); } break; } @@ -745,16 +623,10 @@ struct dielectronQCMC { fRegistry.fill(HIST("Pair/bbbar/b2c2e_b2c2e/hadron_hadron/hs"), v12.M(), v12.Pt(), abs(dphi), abs(cos_thetaCS), abs(phiCS), aco, asym, abs(dphi_e_ee), dca_ee_3d); if (isCharmMeson(mp1) && isCharmMeson(mp2)) { fRegistry.fill(HIST("Pair/bbbar/b2c2e_b2c2e/meson_meson/hs"), v12.M(), v12.Pt(), abs(dphi), abs(cos_thetaCS), abs(phiCS), aco, asym, abs(dphi_e_ee), dca_ee_3d); - fillTrackInfo<8, TMCParticles>(t1); - fillTrackInfo<8, TMCParticles>(t2); } else if (isCharmBaryon(mp1) && isCharmBaryon(mp2)) { fRegistry.fill(HIST("Pair/bbbar/b2c2e_b2c2e/baryon_baryon/hs"), v12.M(), v12.Pt(), abs(dphi), abs(cos_thetaCS), abs(phiCS), aco, asym, abs(dphi_e_ee), dca_ee_3d); - fillTrackInfo<8, TMCParticles>(t1); - fillTrackInfo<8, TMCParticles>(t2); } else { fRegistry.fill(HIST("Pair/bbbar/b2c2e_b2c2e/meson_baryon/hs"), v12.M(), v12.Pt(), abs(dphi), abs(cos_thetaCS), abs(phiCS), aco, asym, abs(dphi_e_ee), dca_ee_3d); - fillTrackInfo<8, TMCParticles>(t1); - fillTrackInfo<8, TMCParticles>(t2); } break; } @@ -767,13 +639,6 @@ struct dielectronQCMC { } else { fRegistry.fill(HIST("Pair/bbbar/b2c2e_b2e_sameb/meson_baryon/hs"), v12.M(), v12.Pt(), abs(dphi), abs(cos_thetaCS), abs(phiCS), aco, asym, abs(dphi_e_ee), dca_ee_3d); } - if ((isCharmMeson(mp1) || isCharmBaryon(mp1)) && (isBeautyMeson(mp2) || isBeautyBaryon(mp2))) { - fillTrackInfo<7, TMCParticles>(t1); - fillTrackInfo<8, TMCParticles>(t2); - } else { - fillTrackInfo<8, TMCParticles>(t1); - fillTrackInfo<7, TMCParticles>(t2); - } break; } case static_cast(EM_HFeeType::kBCe_Be_DiffB): // LS @@ -805,13 +670,6 @@ struct dielectronQCMC { } else { fRegistry.fill(HIST("Pair/bbbar/b2c2e_b2e_diffb/meson_baryon/hs"), v12.M(), v12.Pt(), abs(dphi), abs(cos_thetaCS), abs(phiCS), aco, asym, abs(dphi_e_ee), dca_ee_3d); } - if ((isCharmMeson(mp1) || isCharmBaryon(mp1)) && (isBeautyMeson(mp2) || isBeautyBaryon(mp2))) { - fillTrackInfo<7, TMCParticles>(t1); - fillTrackInfo<8, TMCParticles>(t2); - } else { - fillTrackInfo<8, TMCParticles>(t1); - fillTrackInfo<7, TMCParticles>(t2); - } break; } default: @@ -823,7 +681,6 @@ struct dielectronQCMC { return true; } - std::vector used_trackIds; SliceCache cache; Preslice perCollision_track = aod::emprimaryelectron::emeventId; Filter trackFilter = dielectroncuts.cfg_min_pt_track < o2::aod::track::pt && nabs(o2::aod::track::eta) < dielectroncuts.cfg_max_eta_track && o2::aod::track::tpcChi2NCl < dielectroncuts.cfg_max_chi2tpc && o2::aod::track::itsChi2NCl < dielectroncuts.cfg_max_chi2its && nabs(o2::aod::track::dcaXY) < dielectroncuts.cfg_max_dcaxy && nabs(o2::aod::track::dcaZ) < dielectroncuts.cfg_max_dcaz; @@ -836,10 +693,8 @@ struct dielectronQCMC { Filter collisionFilter_centrality = (cfgCentMin < o2::aod::cent::centFT0M && o2::aod::cent::centFT0M < cfgCentMax) || (cfgCentMin < o2::aod::cent::centFT0A && o2::aod::cent::centFT0A < cfgCentMax) || (cfgCentMin < o2::aod::cent::centFT0C && o2::aod::cent::centFT0C < cfgCentMax); using FilteredMyCollisions = soa::Filtered; - void processQCMC(FilteredMyCollisions const& collisions, FilteredMyMCTracks const& tracks, aod::EMMCParticles const& mcparticles, aod::EMMCEvents const&) + void processQCMC(FilteredMyCollisions const& collisions, FilteredMyMCTracks const&, aod::EMMCParticles const& mcparticles, aod::EMMCEvents const&) { - used_trackIds.reserve(tracks.size()); - for (auto& collision : collisions) { initCCDB(collision); float centralities[3] = {collision.centFT0M(), collision.centFT0A(), collision.centFT0C()}; @@ -873,9 +728,6 @@ struct dielectronQCMC { } // end of collision loop - used_trackIds.clear(); - used_trackIds.shrink_to_fit(); - } // end of process PROCESS_SWITCH(dielectronQCMC, processQCMC, "run dielectron QC MC", true); diff --git a/PWGEM/Dilepton/Tasks/dimuonQC.cxx b/PWGEM/Dilepton/Tasks/dimuonQC.cxx index 71ff6104497..1c972335501 100644 --- a/PWGEM/Dilepton/Tasks/dimuonQC.cxx +++ b/PWGEM/Dilepton/Tasks/dimuonQC.cxx @@ -153,6 +153,14 @@ struct dimuonQC { void init(InitContext& /*context*/) { + mRunNumber = 0; + d_bz = 0; + + ccdb->setURL(ccdburl); + ccdb->setCaching(true); + ccdb->setLocalObjectValidityChecking(); + ccdb->setFatalWhenNull(false); + zvtx_bin_edges = std::vector(ConfVtxBins.value.begin(), ConfVtxBins.value.end()); zvtx_bin_edges.erase(zvtx_bin_edges.begin()); @@ -172,14 +180,6 @@ struct dimuonQC { DefineDimuonCut(); addhistograms(); - mRunNumber = 0; - d_bz = 0; - - ccdb->setURL(ccdburl); - ccdb->setCaching(true); - ccdb->setLocalObjectValidityChecking(); - ccdb->setFatalWhenNull(false); - // fitter.setPropagateToPCA(true); // fitter.setMaxR(90.f); // fitter.setMinParamChange(1e-3); @@ -333,24 +333,6 @@ struct dimuonQC { fRegistry.addClone("Pair/same/", "Pair/mix/"); } - // for track info - fRegistry.add("Track/hPt", "pT;p_{T} (GeV/c)", kTH1F, {{1000, 0.0f, 10}}, false); - fRegistry.add("Track/hQoverPt", "q/pT;q/p_{T} (GeV/c)^{-1}", kTH1F, {{400, -20, 20}}, false); - fRegistry.add("Track/hEtaPhi", "#eta vs. #varphi;#varphi (rad.);#eta", kTH2F, {{180, 0, 2 * M_PI}, {25, -4.5f, -2.0f}}, false); - fRegistry.add("Track/hTrackType", "track type", kTH1F, {{6, -0.5f, 5.5}}, false); - fRegistry.add("Track/hDCAxy", "DCA x vs. y;DCA_{x} (cm);DCA_{y} (cm)", kTH2F, {{200, -1.0f, 1.0f}, {200, -1.0f, 1.0f}}, false); - fRegistry.add("Track/hDCAxySigma", "DCA x vs. y;DCA_{x} (#sigma);DCA_{y} (#sigma)", kTH2F, {{200, -10.0f, 10.0f}, {200, -10.0f, 10.0f}}, false); - fRegistry.add("Track/hDCA2DSigma", "DCA xy;DCA_{xy} (#sigma)", kTH1F, {{100, 0.0f, 10.0f}}, false); - fRegistry.add("Track/hDCAxRes_Pt", "DCA_{x} resolution vs. pT;p_{T} (GeV/c);DCA_{x} resolution (#mum)", kTH2F, {{1000, 0, 10}, {200, 0., 200}}, false); - fRegistry.add("Track/hDCAyRes_Pt", "DCA_{y} resolution vs. pT;p_{T} (GeV/c);DCA_{y} resolution (#mum)", kTH2F, {{1000, 0, 10}, {200, 0., 200}}, false); - fRegistry.add("Track/hNclsMCH", "number of MCH clusters", kTH1F, {{21, -0.5, 20.5}}, false); - fRegistry.add("Track/hNclsMFT", "number of MFT clusters", kTH1F, {{11, -0.5, 10.5}}, false); - fRegistry.add("Track/hPDCA", "pDCA;p_{T} at PV (GeV/c);p #times DCA (GeV/c #upoint cm)", kTH2F, {{100, 0, 10}, {100, 0.0f, 1000}}, false); - fRegistry.add("Track/hChi2", "chi2;chi2", kTH1F, {{100, 0.0f, 100}}, false); - fRegistry.add("Track/hChi2MatchMCHMID", "chi2 match MCH-MID;chi2", kTH1F, {{100, 0.0f, 100}}, false); - fRegistry.add("Track/hChi2MatchMCHMFT", "chi2 match MCH-MFT;chi2", kTH1F, {{100, 0.0f, 100}}, false); - fRegistry.add("Track/hMFTClusterMap", "MFT cluster map", kTH1F, {{1024, -0.5, 1023.5}}, false); - // event info if (nmod == 2) { o2::aod::pwgem::dilepton::utils::eventhistogram::addEventHistograms<2>(&fRegistry); @@ -623,7 +605,6 @@ struct dimuonQC { if (std::find(used_trackIds.begin(), used_trackIds.end(), pair_tmp_id1) == used_trackIds.end()) { used_trackIds.emplace_back(pair_tmp_id1); - fillTrackInfo(t1); if (cfgDoMix) { if (t1.sign() > 0) { emh_pos->AddTrackToEventPool(key_df_collision, EMFwdTrackWithCov(t1.globalIndex(), collision.globalIndex(), t1.fwdtrackId(), t1.pt(), t1.eta(), t1.phi(), o2::constants::physics::MassMuon, t1.sign(), t1.fwdDcaX(), t1.fwdDcaY(), possibleIds1, @@ -638,7 +619,6 @@ struct dimuonQC { } if (std::find(used_trackIds.begin(), used_trackIds.end(), pair_tmp_id2) == used_trackIds.end()) { used_trackIds.emplace_back(pair_tmp_id2); - fillTrackInfo(t2); if (cfgDoMix) { if (t2.sign() > 0) { emh_pos->AddTrackToEventPool(key_df_collision, EMFwdTrackWithCov(t2.globalIndex(), collision.globalIndex(), t2.fwdtrackId(), t2.pt(), t2.eta(), t2.phi(), o2::constants::physics::MassMuon, t2.sign(), t2.fwdDcaX(), t2.fwdDcaY(), possibleIds2, @@ -655,28 +635,6 @@ struct dimuonQC { return true; } - template - void fillTrackInfo(TTrack const& track) - { - float dca_xy = fwdDcaXYinSigma(track); - fRegistry.fill(HIST("Track/hPt"), track.pt()); - fRegistry.fill(HIST("Track/hQoverPt"), track.sign() / track.pt()); - fRegistry.fill(HIST("Track/hEtaPhi"), track.phi(), track.eta()); - fRegistry.fill(HIST("Track/hTrackType"), track.trackType()); - fRegistry.fill(HIST("Track/hDCAxy"), track.fwdDcaX(), track.fwdDcaY()); - fRegistry.fill(HIST("Track/hDCAxySigma"), track.fwdDcaX() / sqrt(track.cXX()), track.fwdDcaY() / sqrt(track.cYY())); - fRegistry.fill(HIST("Track/hDCA2DSigma"), dca_xy); - fRegistry.fill(HIST("Track/hDCAxRes_Pt"), track.pt(), sqrt(track.cXX()) * 1e+4); // convert cm to um - fRegistry.fill(HIST("Track/hDCAyRes_Pt"), track.pt(), sqrt(track.cYY()) * 1e+4); // convert cm to um - fRegistry.fill(HIST("Track/hNclsMFT"), track.nClustersMFT()); - fRegistry.fill(HIST("Track/hNclsMCH"), track.nClusters()); - fRegistry.fill(HIST("Track/hPDCA"), track.pt(), track.pDca()); - fRegistry.fill(HIST("Track/hChi2"), track.chi2()); - fRegistry.fill(HIST("Track/hChi2MatchMCHMFT"), track.chi2MatchMCHMFT()); - fRegistry.fill(HIST("Track/hChi2MatchMCHMID"), track.chi2MatchMCHMID()); - fRegistry.fill(HIST("Track/hMFTClusterMap"), track.mftClusterMap()); - } - Filter collisionFilter_centrality = (cfgCentMin < o2::aod::cent::centFT0M && o2::aod::cent::centFT0M < cfgCentMax) || (cfgCentMin < o2::aod::cent::centFT0A && o2::aod::cent::centFT0A < cfgCentMax) || (cfgCentMin < o2::aod::cent::centFT0C && o2::aod::cent::centFT0C < cfgCentMax); using FilteredMyCollisions = soa::Filtered; diff --git a/PWGEM/Dilepton/Tasks/dimuonQCMC.cxx b/PWGEM/Dilepton/Tasks/dimuonQCMC.cxx index 375314f7c56..5452af348b4 100644 --- a/PWGEM/Dilepton/Tasks/dimuonQCMC.cxx +++ b/PWGEM/Dilepton/Tasks/dimuonQCMC.cxx @@ -46,7 +46,7 @@ using namespace o2::soa; using namespace o2::aod::pwgem::dilepton::utils::mcutil; using namespace o2::aod::pwgem::dilepton::utils::emtrackutil; -using MyCollisions = soa::Join; +using MyCollisions = soa::Join; using MyCollision = MyCollisions::iterator; using MyMCTracks = soa::Join; @@ -129,7 +129,6 @@ struct dimuonQCMC { HistogramRegistry fRegistry{"output", {}, OutputObjHandlingPolicy::AnalysisObject, false, false}; static constexpr std::string_view event_cut_types[2] = {"before/", "after/"}; - static constexpr std::string_view ele_source_types[8] = {"lf/", "PromptJPsi/", "NonPromptJPsi/", "PromptPsi2S/", "NonPromptPsi2S/", "c2mu/", "b2mu/", "b2c2mu/"}; ~dimuonQCMC() {} @@ -206,42 +205,10 @@ struct dimuonQCMC { fRegistry.addClone("Pair/corr_bkg_muh/uls/", "Pair/corr_bkg_muh/lsmm/"); fRegistry.addClone("Pair/corr_bkg_muh/", "Pair/corr_bkg_hh/"); fRegistry.addClone("Pair/corr_bkg_muh/", "Pair/comb_bkg/"); - - // track info - fRegistry.add("Track/lf/hPt", "pT;p_{T} (GeV/c)", kTH1F, {{1000, 0.0f, 10}}, false); - fRegistry.add("Track/lf/hQoverPt", "q/pT;q/p_{T} (GeV/c)^{-1}", kTH1F, {{400, -20, 20}}, false); - fRegistry.add("Track/lf/hEtaPhi", "#eta vs. #varphi;#varphi (rad.);#eta", kTH2F, {{180, 0, 2 * M_PI}, {30, -4.5f, -2.0f}}, false); - fRegistry.add("Track/lf/hTrackType", "track type", kTH1F, {{6, -0.5f, 5.5}}, false); - fRegistry.add("Track/lf/hDCAxy", "DCA x vs. y;DCA_{x} (cm);DCA_{y} (cm)", kTH2F, {{200, -1.0f, 1.0f}, {200, -1.0f, 1.0f}}, false); - fRegistry.add("Track/lf/hDCAxySigma", "DCA x vs. y;DCA_{x} (#sigma);DCA_{y} (#sigma)", kTH2F, {{200, -10.0f, 10.0f}, {200, -10.0f, 10.0f}}, false); - fRegistry.add("Track/lf/hDCA2DSigma", "DCA xy;DCA_{xy} (#sigma)", kTH1F, {{100, 0.0f, 10.0f}}, false); - fRegistry.add("Track/lf/hDCAxRes_Pt", "DCA_{x} resolution vs. pT;p_{T} (GeV/c);DCA_{x} resolution (#mum)", kTH2F, {{1000, 0, 10}, {200, 0., 200}}, false); - fRegistry.add("Track/lf/hDCAyRes_Pt", "DCA_{y} resolution vs. pT;p_{T} (GeV/c);DCA_{y} resolution (#mum)", kTH2F, {{1000, 0, 10}, {200, 0., 200}}, false); - fRegistry.add("Track/lf/hNclsMCH", "number of MCH clusters", kTH1F, {{21, -0.5, 20.5}}, false); - fRegistry.add("Track/lf/hNclsMFT", "number of MFT clusters", kTH1F, {{11, -0.5, 10.5}}, false); - fRegistry.add("Track/lf/hPDCA", "pDCA;p_{T} at PV (GeV/c);p #times DCA (GeV/c #upoint cm)", kTH2F, {{100, 0, 10}, {100, 0.0f, 1000}}, false); - fRegistry.add("Track/lf/hChi2", "chi2;chi2", kTH1F, {{100, 0.0f, 100}}, false); - fRegistry.add("Track/lf/hChi2MatchMCHMID", "chi2 match MCH-MID;chi2", kTH1F, {{100, 0.0f, 100}}, false); - fRegistry.add("Track/lf/hChi2MatchMCHMFT", "chi2 match MCH-MFT;chi2", kTH1F, {{100, 0.0f, 100}}, false); - fRegistry.add("Track/lf/hMFTClusterMap", "MFT cluster map", kTH1F, {{1024, -0.5, 1023.5}}, false); - fRegistry.add("Track/lf/hPtGen_DeltaPtOverPtGen", "muon p_{T} resolution;p_{T}^{gen} (GeV/c);(p_{T}^{rec} - p_{T}^{gen})/p_{T}^{gen}", kTH2F, {{1000, 0, 10}, {400, -1.0f, 1.0f}}, true); - fRegistry.add("Track/lf/hPtGen_DeltaEta", "muon #eta resolution;p_{T}^{gen} (GeV/c);#eta^{rec} - #eta^{gen}", kTH2F, {{1000, 0, 10}, {400, -1.0f, 1.0f}}, true); - fRegistry.add("Track/lf/hPtGen_DeltaPhi", "muon #varphi resolution;p_{T}^{gen} (GeV/c);#varphi^{rec} - #varphi^{gen} (rad.)", kTH2F, {{1000, 0, 10}, {400, -1.0f, 1.0f}}, true); - fRegistry.addClone("Track/lf/", "Track/PromptJPsi/"); - fRegistry.addClone("Track/lf/", "Track/NonPromptJPsi/"); - fRegistry.addClone("Track/lf/", "Track/PromptPsi2S/"); - fRegistry.addClone("Track/lf/", "Track/NonPromptPsi2S/"); - fRegistry.addClone("Track/lf/", "Track/c2mu/"); - fRegistry.addClone("Track/lf/", "Track/b2mu/"); - fRegistry.addClone("Track/lf/", "Track/b2c2mu/"); } void init(InitContext&) { - DefineEMEventCut(); - DefineDimuonCut(); - addhistograms(); - mRunNumber = 0; d_bz = 0; @@ -250,6 +217,10 @@ struct dimuonQCMC { ccdb->setLocalObjectValidityChecking(); ccdb->setFatalWhenNull(false); + DefineEMEventCut(); + DefineDimuonCut(); + addhistograms(); + // fitter.setPropagateToPCA(true); // fitter.setMaxR(90.f); // fitter.setMinParamChange(1e-3); @@ -388,36 +359,6 @@ struct dimuonQCMC { } } - template - void fillTrackInfo(TTrack const& track) - { - // fill track info that belong to true pairs. - if (std::find(used_trackIds.begin(), used_trackIds.end(), track.globalIndex()) == used_trackIds.end()) { - auto mctrack = track.template emmcparticle_as(); - float dca_xy = fwdDcaXYinSigma(track); - fRegistry.fill(HIST("Track/") + HIST(ele_source_types[e_source_id]) + HIST("hPt"), track.pt()); - fRegistry.fill(HIST("Track/") + HIST(ele_source_types[e_source_id]) + HIST("hQoverPt"), track.sign() / track.pt()); - fRegistry.fill(HIST("Track/") + HIST(ele_source_types[e_source_id]) + HIST("hEtaPhi"), track.phi(), track.eta()); - fRegistry.fill(HIST("Track/") + HIST(ele_source_types[e_source_id]) + HIST("hTrackType"), track.trackType()); - fRegistry.fill(HIST("Track/") + HIST(ele_source_types[e_source_id]) + HIST("hDCAxy"), track.fwdDcaX(), track.fwdDcaY()); - fRegistry.fill(HIST("Track/") + HIST(ele_source_types[e_source_id]) + HIST("hDCAxySigma"), track.fwdDcaX() / std::sqrt(track.cXX()), track.fwdDcaY() / std::sqrt(track.cYY())); - fRegistry.fill(HIST("Track/") + HIST(ele_source_types[e_source_id]) + HIST("hDCA2DSigma"), dca_xy); - fRegistry.fill(HIST("Track/") + HIST(ele_source_types[e_source_id]) + HIST("hDCAxRes_Pt"), track.pt(), std::sqrt(track.cXX()) * 1e+4); - fRegistry.fill(HIST("Track/") + HIST(ele_source_types[e_source_id]) + HIST("hDCAyRes_Pt"), track.pt(), std::sqrt(track.cYY()) * 1e+4); - fRegistry.fill(HIST("Track/") + HIST(ele_source_types[e_source_id]) + HIST("hNclsMCH"), track.nClusters()); - fRegistry.fill(HIST("Track/") + HIST(ele_source_types[e_source_id]) + HIST("hNclsMFT"), track.nClustersMFT()); - fRegistry.fill(HIST("Track/") + HIST(ele_source_types[e_source_id]) + HIST("hPDCA"), track.pt(), track.pDca()); - fRegistry.fill(HIST("Track/") + HIST(ele_source_types[e_source_id]) + HIST("hChi2"), track.chi2()); - fRegistry.fill(HIST("Track/") + HIST(ele_source_types[e_source_id]) + HIST("hChi2MatchMCHMID"), track.chi2MatchMCHMID()); - fRegistry.fill(HIST("Track/") + HIST(ele_source_types[e_source_id]) + HIST("hChi2MatchMCHMFT"), track.chi2MatchMCHMFT()); - fRegistry.fill(HIST("Track/") + HIST(ele_source_types[e_source_id]) + HIST("hMFTClusterMap"), track.mftClusterMap()); - fRegistry.fill(HIST("Track/") + HIST(ele_source_types[e_source_id]) + HIST("hPtGen_DeltaPtOverPtGen"), mctrack.pt(), (track.pt() - mctrack.pt()) / mctrack.pt()); - fRegistry.fill(HIST("Track/") + HIST(ele_source_types[e_source_id]) + HIST("hPtGen_DeltaEta"), mctrack.pt(), track.eta() - mctrack.eta()); - fRegistry.fill(HIST("Track/") + HIST(ele_source_types[e_source_id]) + HIST("hPtGen_DeltaPhi"), mctrack.pt(), track.phi() - mctrack.phi()); - used_trackIds.emplace_back(track.globalIndex()); - } - } - template bool fillTruePairInfo(TCollision const& collision, TTrack1 const& t1, TTrack2 const& t2, TMCParticles const& mcparticles) { @@ -507,31 +448,21 @@ struct dimuonQCMC { switch (abs(mcmother.pdgCode())) { case 221: fRegistry.fill(HIST("Pair/sm/Eta/hs"), v12.M(), v12.Pt(), abs(dphi), abs(cos_thetaCS), abs(phiCS), aco, asym, abs(dphi_e_ee), dca_mumu_xy); - fillTrackInfo<0, TMCParticles>(t1); - fillTrackInfo<0, TMCParticles>(t2); break; case 331: fRegistry.fill(HIST("Pair/sm/EtaPrime/hs"), v12.M(), v12.Pt(), abs(dphi), abs(cos_thetaCS), abs(phiCS), aco, asym, abs(dphi_e_ee), dca_mumu_xy); - fillTrackInfo<0, TMCParticles>(t1); - fillTrackInfo<0, TMCParticles>(t2); break; case 113: fRegistry.fill(HIST("Pair/sm/Rho/hs"), v12.M(), v12.Pt(), abs(dphi), abs(cos_thetaCS), abs(phiCS), aco, asym, abs(dphi_e_ee), dca_mumu_xy); - fillTrackInfo<0, TMCParticles>(t1); - fillTrackInfo<0, TMCParticles>(t2); break; case 223: fRegistry.fill(HIST("Pair/sm/Omega/hs"), v12.M(), v12.Pt(), abs(dphi), abs(cos_thetaCS), abs(phiCS), aco, asym, abs(dphi_e_ee), dca_mumu_xy); - fillTrackInfo<0, TMCParticles>(t1); - fillTrackInfo<0, TMCParticles>(t2); if (mcmother.daughtersIds().size() == 2) { // omeag->mumu fRegistry.fill(HIST("Pair/sm/Omega2mumu/hs"), v12.M(), v12.Pt(), abs(dphi), abs(cos_thetaCS), abs(phiCS), aco, asym, abs(dphi_e_ee), dca_mumu_xy); } break; case 333: fRegistry.fill(HIST("Pair/sm/Phi/hs"), v12.M(), v12.Pt(), abs(dphi), abs(cos_thetaCS), abs(phiCS), aco, asym, abs(dphi_e_ee), dca_mumu_xy); - fillTrackInfo<0, TMCParticles>(t1); - fillTrackInfo<0, TMCParticles>(t2); if (mcmother.daughtersIds().size() == 2) { // omeag->mumu fRegistry.fill(HIST("Pair/sm/Phi2mumu/hs"), v12.M(), v12.Pt(), abs(dphi), abs(cos_thetaCS), abs(phiCS), aco, asym, abs(dphi_e_ee), dca_mumu_xy); } @@ -539,24 +470,16 @@ struct dimuonQCMC { case 443: { if (IsFromBeauty(mcmother, mcparticles) > 0) { fRegistry.fill(HIST("Pair/sm/NonPromptJPsi/hs"), v12.M(), v12.Pt(), abs(dphi), abs(cos_thetaCS), abs(phiCS), aco, asym, abs(dphi_e_ee), dca_mumu_xy); - fillTrackInfo<2, TMCParticles>(t1); - fillTrackInfo<2, TMCParticles>(t2); } else { fRegistry.fill(HIST("Pair/sm/PromptJPsi/hs"), v12.M(), v12.Pt(), abs(dphi), abs(cos_thetaCS), abs(phiCS), aco, asym, abs(dphi_e_ee), dca_mumu_xy); - fillTrackInfo<1, TMCParticles>(t1); - fillTrackInfo<1, TMCParticles>(t2); } break; } case 100443: { if (IsFromBeauty(mcmother, mcparticles) > 0) { fRegistry.fill(HIST("Pair/sm/NonPromptPsi2S/hs"), v12.M(), v12.Pt(), abs(dphi), abs(cos_thetaCS), abs(phiCS), aco, asym, abs(dphi_e_ee), dca_mumu_xy); - fillTrackInfo<4, TMCParticles>(t1); - fillTrackInfo<4, TMCParticles>(t2); } else { fRegistry.fill(HIST("Pair/sm/PromptPsi2S/hs"), v12.M(), v12.Pt(), abs(dphi), abs(cos_thetaCS), abs(phiCS), aco, asym, abs(dphi_e_ee), dca_mumu_xy); - fillTrackInfo<3, TMCParticles>(t1); - fillTrackInfo<3, TMCParticles>(t2); } break; } @@ -575,16 +498,10 @@ struct dimuonQCMC { fRegistry.fill(HIST("Pair/ccbar/c2mu_c2mu/hadron_hadron/hs"), v12.M(), v12.Pt(), abs(dphi), abs(cos_thetaCS), abs(phiCS), aco, asym, abs(dphi_e_ee), dca_mumu_xy); if (isCharmMeson(mp1) && isCharmMeson(mp2)) { fRegistry.fill(HIST("Pair/ccbar/c2mu_c2mu/meson_meson/hs"), v12.M(), v12.Pt(), abs(dphi), abs(cos_thetaCS), abs(phiCS), aco, asym, abs(dphi_e_ee), dca_mumu_xy); - fillTrackInfo<5, TMCParticles>(t1); - fillTrackInfo<5, TMCParticles>(t2); } else if (isCharmBaryon(mp1) && isCharmBaryon(mp2)) { fRegistry.fill(HIST("Pair/ccbar/c2mu_c2mu/baryon_baryon/hs"), v12.M(), v12.Pt(), abs(dphi), abs(cos_thetaCS), abs(phiCS), aco, asym, abs(dphi_e_ee), dca_mumu_xy); - fillTrackInfo<5, TMCParticles>(t1); - fillTrackInfo<5, TMCParticles>(t2); } else { fRegistry.fill(HIST("Pair/ccbar/c2mu_c2mu/meson_baryon/hs"), v12.M(), v12.Pt(), abs(dphi), abs(cos_thetaCS), abs(phiCS), aco, asym, abs(dphi_e_ee), dca_mumu_xy); - fillTrackInfo<5, TMCParticles>(t1); - fillTrackInfo<5, TMCParticles>(t2); } break; } @@ -592,16 +509,10 @@ struct dimuonQCMC { fRegistry.fill(HIST("Pair/bbbar/b2mu_b2mu/hadron_hadron/hs"), v12.M(), v12.Pt(), abs(dphi), abs(cos_thetaCS), abs(phiCS), aco, asym, abs(dphi_e_ee), dca_mumu_xy); if (isBeautyMeson(mp1) && isBeautyMeson(mp2)) { fRegistry.fill(HIST("Pair/bbbar/b2mu_b2mu/meson_meson/hs"), v12.M(), v12.Pt(), abs(dphi), abs(cos_thetaCS), abs(phiCS), aco, asym, abs(dphi_e_ee), dca_mumu_xy); - fillTrackInfo<6, TMCParticles>(t1); - fillTrackInfo<6, TMCParticles>(t2); } else if (isBeautyBaryon(mp1) && isBeautyBaryon(mp2)) { fRegistry.fill(HIST("Pair/bbbar/b2mu_b2mu/baryon_baryon/hs"), v12.M(), v12.Pt(), abs(dphi), abs(cos_thetaCS), abs(phiCS), aco, asym, abs(dphi_e_ee), dca_mumu_xy); - fillTrackInfo<6, TMCParticles>(t1); - fillTrackInfo<6, TMCParticles>(t2); } else { fRegistry.fill(HIST("Pair/bbbar/b2mu_b2mu/meson_baryon/hs"), v12.M(), v12.Pt(), abs(dphi), abs(cos_thetaCS), abs(phiCS), aco, asym, abs(dphi_e_ee), dca_mumu_xy); - fillTrackInfo<6, TMCParticles>(t1); - fillTrackInfo<6, TMCParticles>(t2); } break; } @@ -609,16 +520,10 @@ struct dimuonQCMC { fRegistry.fill(HIST("Pair/bbbar/b2c2mu_b2c2mu/hadron_hadron/hs"), v12.M(), v12.Pt(), abs(dphi), abs(cos_thetaCS), abs(phiCS), aco, asym, abs(dphi_e_ee), dca_mumu_xy); if (isCharmMeson(mp1) && isCharmMeson(mp2)) { fRegistry.fill(HIST("Pair/bbbar/b2c2mu_b2c2mu/meson_meson/hs"), v12.M(), v12.Pt(), abs(dphi), abs(cos_thetaCS), abs(phiCS), aco, asym, abs(dphi_e_ee), dca_mumu_xy); - fillTrackInfo<7, TMCParticles>(t1); - fillTrackInfo<7, TMCParticles>(t2); } else if (isCharmBaryon(mp1) && isCharmBaryon(mp2)) { fRegistry.fill(HIST("Pair/bbbar/b2c2mu_b2c2mu/baryon_baryon/hs"), v12.M(), v12.Pt(), abs(dphi), abs(cos_thetaCS), abs(phiCS), aco, asym, abs(dphi_e_ee), dca_mumu_xy); - fillTrackInfo<7, TMCParticles>(t1); - fillTrackInfo<7, TMCParticles>(t2); } else { fRegistry.fill(HIST("Pair/bbbar/b2c2mu_b2c2mu/meson_baryon/hs"), v12.M(), v12.Pt(), abs(dphi), abs(cos_thetaCS), abs(phiCS), aco, asym, abs(dphi_e_ee), dca_mumu_xy); - fillTrackInfo<7, TMCParticles>(t1); - fillTrackInfo<7, TMCParticles>(t2); } break; } @@ -631,13 +536,6 @@ struct dimuonQCMC { } else { fRegistry.fill(HIST("Pair/bbbar/b2c2mu_b2mu_sameb/meson_baryon/hs"), v12.M(), v12.Pt(), abs(dphi), abs(cos_thetaCS), abs(phiCS), aco, asym, abs(dphi_e_ee), dca_mumu_xy); } - if ((isCharmMeson(mp1) || isCharmBaryon(mp1)) && (isBeautyMeson(mp2) || isBeautyBaryon(mp2))) { - fillTrackInfo<6, TMCParticles>(t1); - fillTrackInfo<7, TMCParticles>(t2); - } else { - fillTrackInfo<7, TMCParticles>(t1); - fillTrackInfo<6, TMCParticles>(t2); - } break; } case static_cast(EM_HFeeType::kBCe_Be_DiffB): // LS @@ -669,13 +567,6 @@ struct dimuonQCMC { } else { fRegistry.fill(HIST("Pair/bbbar/b2c2mu_b2mu_diffb/meson_baryon/hs"), v12.M(), v12.Pt(), abs(dphi), abs(cos_thetaCS), abs(phiCS), aco, asym, abs(dphi_e_ee), dca_mumu_xy); } - if ((isCharmMeson(mp1) || isCharmBaryon(mp1)) && (isBeautyMeson(mp2) || isBeautyBaryon(mp2))) { - fillTrackInfo<6, TMCParticles>(t1); - fillTrackInfo<7, TMCParticles>(t2); - } else { - fillTrackInfo<7, TMCParticles>(t1); - fillTrackInfo<6, TMCParticles>(t2); - } break; } default: @@ -687,7 +578,6 @@ struct dimuonQCMC { return true; } - std::vector used_trackIds; SliceCache cache; Preslice perCollision_track = aod::emprimarymuon::emeventId; Filter trackFilter = o2::aod::fwdtrack::trackType == dimuoncuts.cfg_track_type && dimuoncuts.cfg_min_pt_track < o2::aod::fwdtrack::pt && dimuoncuts.cfg_min_eta_track < o2::aod::fwdtrack::eta && o2::aod::fwdtrack::eta < dimuoncuts.cfg_max_eta_track; @@ -700,10 +590,8 @@ struct dimuonQCMC { Filter collisionFilter_centrality = (cfgCentMin < o2::aod::cent::centFT0M && o2::aod::cent::centFT0M < cfgCentMax) || (cfgCentMin < o2::aod::cent::centFT0A && o2::aod::cent::centFT0A < cfgCentMax) || (cfgCentMin < o2::aod::cent::centFT0C && o2::aod::cent::centFT0C < cfgCentMax); using FilteredMyCollisions = soa::Filtered; - void processQCMC(FilteredMyCollisions const& collisions, FilteredMyMCTracks const& tracks, aod::EMMCParticles const& mcparticles, aod::EMMCEvents const&) + void processQCMC(FilteredMyCollisions const& collisions, FilteredMyMCTracks const&, aod::EMMCParticles const& mcparticles, aod::EMMCEvents const&) { - used_trackIds.reserve(tracks.size()); - for (auto& collision : collisions) { initCCDB(collision); float centralities[3] = {collision.centFT0M(), collision.centFT0A(), collision.centFT0C()}; @@ -737,9 +625,6 @@ struct dimuonQCMC { } // end of collision loop - used_trackIds.clear(); - used_trackIds.shrink_to_fit(); - } // end of process PROCESS_SWITCH(dimuonQCMC, processQCMC, "run dimuon QC MC", true); diff --git a/PWGEM/Dilepton/Tasks/singleElectronQC.cxx b/PWGEM/Dilepton/Tasks/singleElectronQC.cxx new file mode 100644 index 00000000000..aeb33889a17 --- /dev/null +++ b/PWGEM/Dilepton/Tasks/singleElectronQC.cxx @@ -0,0 +1,335 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. +// +// ======================== +// +// This code runs loop over electrons for QC. +// Please write to: daiki.sekihata@cern.ch + +#include "TString.h" +#include "Math/Vector4D.h" +#include "Framework/runDataProcessing.h" +#include "Framework/AnalysisTask.h" +#include "Framework/ASoAHelpers.h" + +#include "DetectorsBase/GeometryManager.h" +#include "DataFormatsParameters/GRPObject.h" +#include "DataFormatsParameters/GRPMagField.h" +#include "CCDB/BasicCCDBManager.h" +#include "Tools/ML/MlResponse.h" +#include "Tools/ML/model.h" + +#include "PWGEM/Dilepton/DataModel/dileptonTables.h" +#include "PWGEM/Dilepton/Core/DielectronCut.h" +#include "PWGEM/Dilepton/Core/EMEventCut.h" +#include "PWGEM/Dilepton/Utils/EventHistograms.h" +#include "PWGEM/Dilepton/Utils/EMTrackUtilities.h" +#include "PWGEM/Dilepton/Utils/PairUtilities.h" + +using namespace o2; +using namespace o2::aod; +using namespace o2::framework; +using namespace o2::framework::expressions; +using namespace o2::soa; +using namespace o2::aod::pwgem::dilepton::utils::emtrackutil; + +using MyCollisions = soa::Join; +using MyCollision = MyCollisions::iterator; + +using MyTracks = soa::Join; +using MyTrack = MyTracks::iterator; + +struct singleElectronQC { + + // Configurables + Configurable ccdburl{"ccdb-url", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; + + Configurable cfgCentEstimator{"cfgCentEstimator", 2, "FT0M:0, FT0A:1, FT0C:2"}; + Configurable cfgCentMin{"cfgCentMin", 0, "min. centrality"}; + Configurable cfgCentMax{"cfgCentMax", 999.f, "max. centrality"}; + + ConfigurableAxis ConfPteBins{"ConfPteBins", {VARIABLE_WIDTH, 0.00, 0.05, 0.10, 0.15, 0.20, 0.30, 0.40, 0.50, 0.60, 0.70, 0.80, 0.90, 1.00, 1.10, 1.20, 1.30, 1.40, 1.50, 1.60, 1.70, 1.80, 1.90, 2.00, 2.10, 2.20, 2.30, 2.40, 2.50, 2.60, 2.70, 2.80, 2.90, 3.00, 3.10, 3.20, 3.30, 3.40, 3.50, 3.60, 3.70, 3.80, 3.90, 4.00, 4.10, 4.20, 4.30, 4.40, 4.50, 4.60, 4.70, 4.80, 4.90, 5.00, 5.50, 6.00, 6.50, 7.00, 7.50, 8.00, 8.50, 9.00, 9.50, 10.00}, "pTe bins for output histograms"}; + + EMEventCut fEMEventCut; + struct : ConfigurableGroup { + std::string prefix = "eventcut_group"; + Configurable cfgZvtxMax{"cfgZvtxMax", 10.f, "max. Zvtx"}; + Configurable cfgRequireSel8{"cfgRequireSel8", true, "require sel8 in event cut"}; + Configurable cfgRequireFT0AND{"cfgRequireFT0AND", true, "require FT0AND in event cut"}; + Configurable cfgRequireNoTFB{"cfgRequireNoTFB", false, "require No time frame border in event cut"}; + Configurable cfgRequireNoITSROFB{"cfgRequireNoITSROFB", false, "require no ITS readout frame border in event cut"}; + Configurable cfgRequireNoSameBunchPileup{"cfgRequireNoSameBunchPileup", false, "require no same bunch pileup in event cut"}; + Configurable cfgRequireVertexITSTPC{"cfgRequireVertexITSTPC", false, "require Vertex ITSTPC in event cut"}; // ITS-TPC matched track contributes PV. + Configurable cfgRequireGoodZvtxFT0vsPV{"cfgRequireGoodZvtxFT0vsPV", false, "require good Zvtx between FT0 vs. PV in event cut"}; + Configurable cfgOccupancyMin{"cfgOccupancyMin", -1, "min. occupancy"}; + Configurable cfgOccupancyMax{"cfgOccupancyMax", 1000000000, "max. occupancy"}; + } eventcuts; + + DielectronCut fDielectonCut; + struct : ConfigurableGroup { + std::string prefix = "dielectroncut_group"; + + Configurable cfg_min_pt_track{"cfg_min_pt_track", 0.2, "min pT for single track"}; + Configurable cfg_max_eta_track{"cfg_max_eta_track", 0.8, "max eta for single track"}; + Configurable cfg_min_ncluster_tpc{"cfg_min_ncluster_tpc", 0, "min ncluster tpc"}; + Configurable cfg_min_ncluster_its{"cfg_min_ncluster_its", 5, "min ncluster its"}; + Configurable cfg_min_ncrossedrows{"cfg_min_ncrossedrows", 100, "min ncrossed rows"}; + Configurable cfg_max_chi2tpc{"cfg_max_chi2tpc", 4.0, "max chi2/NclsTPC"}; + Configurable cfg_max_chi2its{"cfg_max_chi2its", 5.0, "max chi2/NclsITS"}; + Configurable cfg_max_dcaxy{"cfg_max_dcaxy", 1.0, "max dca XY for single track in cm"}; + Configurable cfg_max_dcaz{"cfg_max_dcaz", 1.0, "max dca Z for single track in cm"}; + Configurable cfg_require_itsib_any{"cfg_require_itsib_any", true, "flag to require ITS ib any hits"}; + Configurable cfg_require_itsib_1st{"cfg_require_itsib_1st", false, "flag to require ITS ib 1st hit"}; + + Configurable cfg_pid_scheme{"cfg_pid_scheme", static_cast(DielectronCut::PIDSchemes::kTPChadrejORTOFreq), "pid scheme [kTOFreq : 0, kTPChadrej : 1, kTPChadrejORTOFreq : 2, kTPConly : 3]"}; + Configurable cfg_min_TPCNsigmaEl{"cfg_min_TPCNsigmaEl", -2.0, "min. TPC n sigma for electron inclusion"}; + Configurable cfg_max_TPCNsigmaEl{"cfg_max_TPCNsigmaEl", +3.0, "max. TPC n sigma for electron inclusion"}; + Configurable cfg_min_TPCNsigmaMu{"cfg_min_TPCNsigmaMu", -0.0, "min. TPC n sigma for muon exclusion"}; + Configurable cfg_max_TPCNsigmaMu{"cfg_max_TPCNsigmaMu", +0.0, "max. TPC n sigma for muon exclusion"}; + Configurable cfg_min_TPCNsigmaPi{"cfg_min_TPCNsigmaPi", -1e+10, "min. TPC n sigma for pion exclusion"}; + Configurable cfg_max_TPCNsigmaPi{"cfg_max_TPCNsigmaPi", +3.0, "max. TPC n sigma for pion exclusion"}; + Configurable cfg_min_TPCNsigmaKa{"cfg_min_TPCNsigmaKa", -3.0, "min. TPC n sigma for kaon exclusion"}; + Configurable cfg_max_TPCNsigmaKa{"cfg_max_TPCNsigmaKa", +3.0, "max. TPC n sigma for kaon exclusion"}; + Configurable cfg_min_TPCNsigmaPr{"cfg_min_TPCNsigmaPr", -3.0, "min. TPC n sigma for proton exclusion"}; + Configurable cfg_max_TPCNsigmaPr{"cfg_max_TPCNsigmaPr", +3.0, "max. TPC n sigma for proton exclusion"}; + Configurable cfg_min_TOFNsigmaEl{"cfg_min_TOFNsigmaEl", -3.0, "min. TOF n sigma for electron inclusion"}; + Configurable cfg_max_TOFNsigmaEl{"cfg_max_TOFNsigmaEl", +3.0, "max. TOF n sigma for electron inclusion"}; + Configurable enableTTCA{"enableTTCA", false, "Flag to enable or disable TTCA"}; // TTCA should not be used for single track to avoid double counting. + + // CCDB configuration for PID ML + Configurable BDTLocalPathGamma{"BDTLocalPathGamma", "pid_ml_xgboost.onnx", "Path to the local .onnx file"}; + Configurable BDTPathCCDB{"BDTPathCCDB", "Users/d/dsekihat/pwgem/pidml/", "Path on CCDB"}; + Configurable timestampCCDB{"timestampCCDB", -1, "timestamp of the ONNX file for ML model used to query in CCDB. Exceptions: > 0 for the specific timestamp, 0 gets the run dependent timestamp"}; + Configurable loadModelsFromCCDB{"loadModelsFromCCDB", false, "Flag to enable or disable the loading of models from CCDB"}; + Configurable enableOptimizations{"enableOptimizations", false, "Enables the ONNX extended model-optimization: sessionOptions.SetGraphOptimizationLevel(GraphOptimizationLevel::ORT_ENABLE_EXTENDED)"}; + } dielectroncuts; + + o2::ccdb::CcdbApi ccdbApi; + Service ccdb; + + HistogramRegistry fRegistry{"output", {}, OutputObjHandlingPolicy::AnalysisObject, false, false}; // 1 HistogramRegistry can keep up to 512 histograms + static constexpr std::string_view event_cut_types[2] = {"before/", "after/"}; + + ~singleElectronQC() + { + if (eid_bdt) { + delete eid_bdt; + } + } + + void addhistograms() + { + // event info + o2::aod::pwgem::dilepton::utils::eventhistogram::addEventHistograms<-1>(&fRegistry); + + const AxisSpec axis_pt{ConfPteBins, "p_{T,e} (GeV/c)"}; + const AxisSpec axis_eta{20, -1.0, +1.0, "#eta_{e}"}; + const AxisSpec axis_phi{36, 0.0, 2 * M_PI, "#varphi_{e} (rad.)"}; + const AxisSpec axis_charge_rec{3, -1.5, +1.5, "charge"}; + + // track info + fRegistry.add("Track/hs", "rec. single electron", kTHnSparseF, {axis_pt, axis_eta, axis_phi, axis_charge_rec}, true); + fRegistry.add("Track/hQoverPt", "q/pT;q/p_{T} (GeV/c)^{-1}", kTH1F, {{400, -20, 20}}, false); + fRegistry.add("Track/hDCAxyz", "DCA xy vs. z;DCA_{xy} (cm);DCA_{z} (cm)", kTH2F, {{200, -1.0f, 1.0f}, {200, -1.0f, 1.0f}}, false); + fRegistry.add("Track/hDCAxyzSigma", "DCA xy vs. z;DCA_{xy} (#sigma);DCA_{z} (#sigma)", kTH2F, {{200, -10.0f, 10.0f}, {200, -10.0f, 10.0f}}, false); + fRegistry.add("Track/hDCA3DSigma", "DCA 3D;DCA_{3D} (#sigma);", kTH1F, {{100, 0.0f, 10.0f}}, false); + fRegistry.add("Track/hDCAxyRes_Pt", "DCA_{xy} resolution vs. pT;p_{T} (GeV/c);DCA_{xy} resolution (#mum)", kTH2F, {{1000, 0, 10}, {500, 0., 500}}, false); + fRegistry.add("Track/hDCAzRes_Pt", "DCA_{z} resolution vs. pT;p_{T} (GeV/c);DCA_{z} resolution (#mum)", kTH2F, {{1000, 0, 10}, {500, 0., 500}}, false); + fRegistry.add("Track/hNclsTPC", "number of TPC clusters", kTH1F, {{161, -0.5, 160.5}}, false); + fRegistry.add("Track/hNcrTPC", "number of TPC crossed rows", kTH1F, {{161, -0.5, 160.5}}, false); + fRegistry.add("Track/hChi2TPC", "chi2/number of TPC clusters", kTH1F, {{100, 0, 10}}, false); + fRegistry.add("Track/hTPCdEdx", "TPC dE/dx;p_{in} (GeV/c);TPC dE/dx (a.u.)", kTH2F, {{1000, 0, 10}, {200, 0, 200}}, false); + fRegistry.add("Track/hTPCNsigmaEl", "TPC n sigma el;p_{in} (GeV/c);n #sigma_{e}^{TPC}", kTH2F, {{1000, 0, 10}, {100, -5, +5}}, false); + fRegistry.add("Track/hTPCNsigmaMu", "TPC n sigma mu;p_{in} (GeV/c);n #sigma_{#mu}^{TPC}", kTH2F, {{1000, 0, 10}, {100, -5, +5}}, false); + fRegistry.add("Track/hTPCNsigmaPi", "TPC n sigma pi;p_{in} (GeV/c);n #sigma_{#pi}^{TPC}", kTH2F, {{1000, 0, 10}, {100, -5, +5}}, false); + fRegistry.add("Track/hTPCNsigmaKa", "TPC n sigma ka;p_{in} (GeV/c);n #sigma_{K}^{TPC}", kTH2F, {{1000, 0, 10}, {100, -5, +5}}, false); + fRegistry.add("Track/hTPCNsigmaPr", "TPC n sigma pr;p_{in} (GeV/c);n #sigma_{p}^{TPC}", kTH2F, {{1000, 0, 10}, {100, -5, +5}}, false); + fRegistry.add("Track/hTOFbeta", "TOF #beta;p_{in} (GeV/c);#beta", kTH2F, {{1000, 0, 10}, {600, 0, 1.2}}, false); + fRegistry.add("Track/h1overTOFbeta", "TOF 1/#beta;p_{in} (GeV/c);1/#beta", kTH2F, {{1000, 0, 10}, {1000, 0.8, 1.8}}, false); + fRegistry.add("Track/hTOFNsigmaEl", "TOF n sigma el;p_{in} (GeV/c);n #sigma_{e}^{TOF}", kTH2F, {{1000, 0, 10}, {100, -5, +5}}, false); + fRegistry.add("Track/hTOFNsigmaMu", "TOF n sigma mu;p_{in} (GeV/c);n #sigma_{#mu}^{TOF}", kTH2F, {{1000, 0, 10}, {100, -5, +5}}, false); + fRegistry.add("Track/hTOFNsigmaPi", "TOF n sigma pi;p_{in} (GeV/c);n #sigma_{#pi}^{TOF}", kTH2F, {{1000, 0, 10}, {100, -5, +5}}, false); + fRegistry.add("Track/hTOFNsigmaKa", "TOF n sigma ka;p_{in} (GeV/c);n #sigma_{K}^{TOF}", kTH2F, {{1000, 0, 10}, {100, -5, +5}}, false); + fRegistry.add("Track/hTOFNsigmaPr", "TOF n sigma pr;p_{in} (GeV/c);n #sigma_{p}^{TOF}", kTH2F, {{1000, 0, 10}, {100, -5, +5}}, false); + fRegistry.add("Track/hTPCNcr2Nf", "TPC Ncr/Nfindable", kTH1F, {{200, 0, 2}}, false); + fRegistry.add("Track/hTPCNcls2Nf", "TPC Ncls/Nfindable", kTH1F, {{200, 0, 2}}, false); + fRegistry.add("Track/hNclsITS", "number of ITS clusters", kTH1F, {{8, -0.5, 7.5}}, false); + fRegistry.add("Track/hChi2ITS", "chi2/number of ITS clusters", kTH1F, {{100, 0, 10}}, false); + fRegistry.add("Track/hITSClusterMap", "ITS cluster map", kTH1F, {{128, -0.5, 127.5}}, false); + fRegistry.add("Track/hMeanClusterSizeITS", "mean cluster size ITS; on ITS #times cos(#lambda)", kTH1F, {{32, 0, 16}}, false); + } + + void init(InitContext&) + { + ccdb->setURL(ccdburl); + ccdb->setCaching(true); + ccdb->setLocalObjectValidityChecking(); + ccdb->setFatalWhenNull(false); + + DefineEMEventCut(); + DefineDileptonCut(); + addhistograms(); + } + + void DefineEMEventCut() + { + fEMEventCut = EMEventCut("fEMEventCut", "fEMEventCut"); + fEMEventCut.SetRequireSel8(eventcuts.cfgRequireSel8); + fEMEventCut.SetRequireFT0AND(eventcuts.cfgRequireFT0AND); + fEMEventCut.SetZvtxRange(-eventcuts.cfgZvtxMax, +eventcuts.cfgZvtxMax); + fEMEventCut.SetRequireNoTFB(eventcuts.cfgRequireNoTFB); + fEMEventCut.SetRequireNoITSROFB(eventcuts.cfgRequireNoITSROFB); + fEMEventCut.SetRequireNoSameBunchPileup(eventcuts.cfgRequireNoSameBunchPileup); + fEMEventCut.SetRequireVertexITSTPC(eventcuts.cfgRequireVertexITSTPC); + fEMEventCut.SetRequireGoodZvtxFT0vsPV(eventcuts.cfgRequireGoodZvtxFT0vsPV); + fEMEventCut.SetOccupancyRange(eventcuts.cfgOccupancyMin, eventcuts.cfgOccupancyMax); + } + + o2::ml::OnnxModel* eid_bdt = nullptr; + void DefineDileptonCut() + { + fDielectonCut = DielectronCut("fDielectonCut", "fDielectonCut"); + + // for track + fDielectonCut.SetTrackPtRange(dielectroncuts.cfg_min_pt_track, 1e+10f); + fDielectonCut.SetTrackEtaRange(-dielectroncuts.cfg_max_eta_track, +dielectroncuts.cfg_max_eta_track); + fDielectonCut.SetMinNClustersTPC(dielectroncuts.cfg_min_ncluster_tpc); + fDielectonCut.SetMinNCrossedRowsTPC(dielectroncuts.cfg_min_ncrossedrows); + fDielectonCut.SetMinNCrossedRowsOverFindableClustersTPC(0.8); + fDielectonCut.SetChi2PerClusterTPC(0.0, dielectroncuts.cfg_max_chi2tpc); + fDielectonCut.SetChi2PerClusterITS(0.0, dielectroncuts.cfg_max_chi2its); + fDielectonCut.SetNClustersITS(dielectroncuts.cfg_min_ncluster_its, 7); + fDielectonCut.SetMeanClusterSizeITSob(0, 16); + fDielectonCut.SetMaxDcaXY(dielectroncuts.cfg_max_dcaxy); + fDielectonCut.SetMaxDcaZ(dielectroncuts.cfg_max_dcaz); + fDielectonCut.RequireITSibAny(dielectroncuts.cfg_require_itsib_any); + fDielectonCut.RequireITSib1st(dielectroncuts.cfg_require_itsib_1st); + + // for eID + fDielectonCut.SetPIDScheme(dielectroncuts.cfg_pid_scheme); + fDielectonCut.SetTPCNsigmaElRange(dielectroncuts.cfg_min_TPCNsigmaEl, dielectroncuts.cfg_max_TPCNsigmaEl); + fDielectonCut.SetTPCNsigmaMuRange(dielectroncuts.cfg_min_TPCNsigmaMu, dielectroncuts.cfg_max_TPCNsigmaMu); + fDielectonCut.SetTPCNsigmaPiRange(dielectroncuts.cfg_min_TPCNsigmaPi, dielectroncuts.cfg_max_TPCNsigmaPi); + fDielectonCut.SetTPCNsigmaKaRange(dielectroncuts.cfg_min_TPCNsigmaKa, dielectroncuts.cfg_max_TPCNsigmaKa); + fDielectonCut.SetTPCNsigmaPrRange(dielectroncuts.cfg_min_TPCNsigmaPr, dielectroncuts.cfg_max_TPCNsigmaPr); + fDielectonCut.SetTOFNsigmaElRange(dielectroncuts.cfg_min_TOFNsigmaEl, dielectroncuts.cfg_max_TOFNsigmaEl); + + if (dielectroncuts.cfg_pid_scheme == static_cast(DielectronCut::PIDSchemes::kPIDML)) { // please call this at the end of DefineDileptonCut + eid_bdt = new o2::ml::OnnxModel(); + if (dielectroncuts.loadModelsFromCCDB) { + ccdbApi.init(ccdburl); + std::map metadata; + bool retrieveSuccessGamma = ccdbApi.retrieveBlob(dielectroncuts.BDTPathCCDB.value, ".", metadata, dielectroncuts.timestampCCDB.value, false, dielectroncuts.BDTLocalPathGamma.value); + if (retrieveSuccessGamma) { + eid_bdt->initModel(dielectroncuts.BDTLocalPathGamma.value, dielectroncuts.enableOptimizations.value); + } else { + LOG(fatal) << "Error encountered while fetching/loading the Gamma model from CCDB! Maybe the model doesn't exist yet for this runnumber/timestamp?"; + } + } else { + eid_bdt->initModel(dielectroncuts.BDTLocalPathGamma.value, dielectroncuts.enableOptimizations.value); + } + + fDielectonCut.SetPIDModel(eid_bdt); + } // end of PID ML + } + + template + void fillTrackInfo(TTrack const& track) + { + float dca_3d = dca3DinSigma(track); + fRegistry.fill(HIST("Track/hs"), track.pt(), track.eta(), track.phi(), track.sign()); + fRegistry.fill(HIST("Track/hQoverPt"), track.sign() / track.pt()); + fRegistry.fill(HIST("Track/hDCAxyz"), track.dcaXY(), track.dcaZ()); + fRegistry.fill(HIST("Track/hDCAxyzSigma"), track.dcaXY() / sqrt(track.cYY()), track.dcaZ() / sqrt(track.cZZ())); + fRegistry.fill(HIST("Track/hDCA3DSigma"), dca_3d); + fRegistry.fill(HIST("Track/hDCAxyRes_Pt"), track.pt(), sqrt(track.cYY()) * 1e+4); // convert cm to um + fRegistry.fill(HIST("Track/hDCAzRes_Pt"), track.pt(), sqrt(track.cZZ()) * 1e+4); // convert cm to um + fRegistry.fill(HIST("Track/hNclsITS"), track.itsNCls()); + fRegistry.fill(HIST("Track/hNclsTPC"), track.tpcNClsFound()); + fRegistry.fill(HIST("Track/hNcrTPC"), track.tpcNClsCrossedRows()); + fRegistry.fill(HIST("Track/hTPCNcr2Nf"), track.tpcCrossedRowsOverFindableCls()); + fRegistry.fill(HIST("Track/hTPCNcls2Nf"), track.tpcFoundOverFindableCls()); + fRegistry.fill(HIST("Track/hChi2TPC"), track.tpcChi2NCl()); + fRegistry.fill(HIST("Track/hChi2ITS"), track.itsChi2NCl()); + fRegistry.fill(HIST("Track/hITSClusterMap"), track.itsClusterMap()); + fRegistry.fill(HIST("Track/hMeanClusterSizeITS"), track.meanClusterSizeITS() * std::cos(std::atan(track.tgl()))); + fRegistry.fill(HIST("Track/hTPCdEdx"), track.tpcInnerParam(), track.tpcSignal()); + fRegistry.fill(HIST("Track/hTPCNsigmaEl"), track.tpcInnerParam(), track.tpcNSigmaEl()); + fRegistry.fill(HIST("Track/hTPCNsigmaMu"), track.tpcInnerParam(), track.tpcNSigmaMu()); + fRegistry.fill(HIST("Track/hTPCNsigmaPi"), track.tpcInnerParam(), track.tpcNSigmaPi()); + fRegistry.fill(HIST("Track/hTPCNsigmaKa"), track.tpcInnerParam(), track.tpcNSigmaKa()); + fRegistry.fill(HIST("Track/hTPCNsigmaPr"), track.tpcInnerParam(), track.tpcNSigmaPr()); + fRegistry.fill(HIST("Track/hTOFbeta"), track.tpcInnerParam(), track.beta()); + fRegistry.fill(HIST("Track/h1overTOFbeta"), track.tpcInnerParam(), 1. / track.beta()); + fRegistry.fill(HIST("Track/hTOFNsigmaEl"), track.tpcInnerParam(), track.tofNSigmaEl()); + fRegistry.fill(HIST("Track/hTOFNsigmaMu"), track.tpcInnerParam(), track.tofNSigmaMu()); + fRegistry.fill(HIST("Track/hTOFNsigmaPi"), track.tpcInnerParam(), track.tofNSigmaPi()); + fRegistry.fill(HIST("Track/hTOFNsigmaKa"), track.tpcInnerParam(), track.tofNSigmaKa()); + fRegistry.fill(HIST("Track/hTOFNsigmaPr"), track.tpcInnerParam(), track.tofNSigmaPr()); + } + + SliceCache cache; + Preslice perCollision_track = aod::emprimaryelectron::emeventId; + Filter trackFilter = dielectroncuts.cfg_min_pt_track < o2::aod::track::pt && nabs(o2::aod::track::eta) < dielectroncuts.cfg_max_eta_track && o2::aod::track::tpcChi2NCl < dielectroncuts.cfg_max_chi2tpc && o2::aod::track::itsChi2NCl < dielectroncuts.cfg_max_chi2its && nabs(o2::aod::track::dcaXY) < dielectroncuts.cfg_max_dcaxy && nabs(o2::aod::track::dcaZ) < dielectroncuts.cfg_max_dcaz; + Filter pidFilter = (dielectroncuts.cfg_min_TPCNsigmaEl < o2::aod::pidtpc::tpcNSigmaEl && o2::aod::pidtpc::tpcNSigmaEl < dielectroncuts.cfg_max_TPCNsigmaEl) && (o2::aod::pidtpc::tpcNSigmaPi < dielectroncuts.cfg_min_TPCNsigmaPi || dielectroncuts.cfg_max_TPCNsigmaPi < o2::aod::pidtpc::tpcNSigmaPi) && ((0.96f < o2::aod::pidtofbeta::beta && o2::aod::pidtofbeta::beta < 1.04f) || o2::aod::pidtofbeta::beta < 0.f); + Filter ttcaFilter = ifnode(dielectroncuts.enableTTCA.node(), o2::aod::emprimaryelectron::isAssociatedToMPC == true || o2::aod::emprimaryelectron::isAssociatedToMPC == false, o2::aod::emprimaryelectron::isAssociatedToMPC == true); + using FilteredMyTracks = soa::Filtered; + + Filter collisionFilter_centrality = (cfgCentMin < o2::aod::cent::centFT0M && o2::aod::cent::centFT0M < cfgCentMax) || (cfgCentMin < o2::aod::cent::centFT0A && o2::aod::cent::centFT0A < cfgCentMax) || (cfgCentMin < o2::aod::cent::centFT0C && o2::aod::cent::centFT0C < cfgCentMax); + using FilteredMyCollisions = soa::Filtered; + + void processQC(FilteredMyCollisions const& collisions, FilteredMyTracks const& tracks) + { + for (auto& collision : collisions) { + float centralities[3] = {collision.centFT0M(), collision.centFT0A(), collision.centFT0C()}; + if (centralities[cfgCentEstimator] < cfgCentMin || cfgCentMax < centralities[cfgCentEstimator]) { + continue; + } + + o2::aod::pwgem::dilepton::utils::eventhistogram::fillEventInfo<0, -1>(&fRegistry, collision); + if (!fEMEventCut.IsSelected(collision)) { + continue; + } + o2::aod::pwgem::dilepton::utils::eventhistogram::fillEventInfo<1, -1>(&fRegistry, collision); + fRegistry.fill(HIST("Event/before/hCollisionCounter"), o2::aod::pwgem::dilepton::utils::eventhistogram::nbin_ev); // accepted + fRegistry.fill(HIST("Event/after/hCollisionCounter"), o2::aod::pwgem::dilepton::utils::eventhistogram::nbin_ev); // accepted + + auto tracks_per_coll = tracks.sliceBy(perCollision_track, collision.globalIndex()); + + for (auto& track : tracks_per_coll) { + if (dielectroncuts.cfg_pid_scheme == static_cast(DielectronCut::PIDSchemes::kPIDML)) { + if (!fDielectonCut.IsSelectedTrack(track, collision)) { + continue; + } + } else { // cut-based + if (!fDielectonCut.IsSelectedTrack(track)) { + continue; + } + } + + fillTrackInfo(track); + + } // end of track loop + + } // end of collision loop + + } // end of process + PROCESS_SWITCH(singleElectronQC, processQC, "run dielectron QC", true); + + void processDummy(MyCollisions const&) {} + PROCESS_SWITCH(singleElectronQC, processDummy, "Dummy function", false); +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{ + adaptAnalysisTask(cfgc, TaskName{"single-electron-qc"})}; +} diff --git a/PWGEM/Dilepton/Tasks/singleElectronQCMC.cxx b/PWGEM/Dilepton/Tasks/singleElectronQCMC.cxx index 32ee54a7c69..72b48e72b01 100644 --- a/PWGEM/Dilepton/Tasks/singleElectronQCMC.cxx +++ b/PWGEM/Dilepton/Tasks/singleElectronQCMC.cxx @@ -43,7 +43,7 @@ using namespace o2::soa; using namespace o2::aod::pwgem::dilepton::utils::mcutil; using namespace o2::aod::pwgem::dilepton::utils::emtrackutil; -using MyCollisions = soa::Join; +using MyCollisions = soa::Join; using MyCollision = MyCollisions::iterator; using MyMCTracks = soa::Join; @@ -53,10 +53,6 @@ struct singleElectronQCMC { // 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"}; - Configurable grpmagPath{"grpmagPath", "GLO/Config/GRPMagField", "CCDB path of the GRPMagField object"}; - Configurable skipGRPOquery{"skipGRPOquery", true, "skip grpo query"}; - Configurable d_bz_input{"d_bz_input", -999, "bz field in kG, -999 is automatic"}; Configurable cfgEventGeneratorType{"cfgEventGeneratorType", -1, "if positive, select event generator type. i.e. gap or signal"}; Configurable cfgCentEstimator{"cfgCentEstimator", 2, "FT0M:0, FT0A:1, FT0C:2"}; @@ -80,7 +76,7 @@ struct singleElectronQCMC { Configurable cfgOccupancyMax{"cfgOccupancyMax", 1000000000, "max. occupancy"}; } eventcuts; - DielectronCut fDielectonCut; + DielectronCut fDielectronCut; struct : ConfigurableGroup { std::string prefix = "dielectroncut_group"; @@ -121,8 +117,6 @@ struct singleElectronQCMC { o2::ccdb::CcdbApi ccdbApi; Service ccdb; - int mRunNumber; - float d_bz; struct : ConfigurableGroup { std::string prefix = "mctrackcut_group"; @@ -207,61 +201,18 @@ struct singleElectronQCMC { fRegistry.addClone("Track/lf/", "Track/b2c2e/"); } - bool cfgDoFlow = false; void init(InitContext&) { DefineEMEventCut(); DefineDileptonCut(); addhistograms(); - mRunNumber = 0; - d_bz = 0; - ccdb->setURL(ccdburl); ccdb->setCaching(true); ccdb->setLocalObjectValidityChecking(); ccdb->setFatalWhenNull(false); } - template - void initCCDB(TCollision const& collision) - { - if (mRunNumber == collision.runNumber()) { - return; - } - - // In case override, don't proceed, please - no CCDB access required - if (d_bz_input > -990) { - d_bz = d_bz_input; - o2::parameters::GRPMagField grpmag; - if (fabs(d_bz) > 1e-5) { - grpmag.setL3Current(30000.f / (d_bz / 5.0f)); - } - mRunNumber = collision.runNumber(); - return; - } - - auto run3grp_timestamp = collision.timestamp(); - o2::parameters::GRPObject* grpo = 0x0; - o2::parameters::GRPMagField* grpmag = 0x0; - if (!skipGRPOquery) - grpo = ccdb->getForTimeStamp(grpPath, run3grp_timestamp); - if (grpo) { - // Fetch magnetic field from ccdb for current collision - d_bz = grpo->getNominalL3Field(); - LOG(info) << "Retrieved GRP for timestamp " << run3grp_timestamp << " with magnetic field of " << d_bz << " kZG"; - } else { - grpmag = ccdb->getForTimeStamp(grpmagPath, run3grp_timestamp); - if (!grpmag) { - LOG(fatal) << "Got nullptr from CCDB for path " << grpmagPath << " of object GRPMagField and " << grpPath << " of object GRPObject for timestamp " << run3grp_timestamp; - } - // Fetch magnetic field from ccdb for current collision - d_bz = std::lround(5.f * grpmag->getL3Current() / 30000.f); - LOG(info) << "Retrieved GRP for timestamp " << run3grp_timestamp << " with magnetic field of " << d_bz << " kZG"; - } - mRunNumber = collision.runNumber(); - } - void DefineEMEventCut() { fEMEventCut = EMEventCut("fEMEventCut", "fEMEventCut"); @@ -279,31 +230,31 @@ struct singleElectronQCMC { o2::ml::OnnxModel* eid_bdt = nullptr; void DefineDileptonCut() { - fDielectonCut = DielectronCut("fDielectonCut", "fDielectonCut"); + fDielectronCut = DielectronCut("fDielectronCut", "fDielectronCut"); // for track - fDielectonCut.SetTrackPtRange(dielectroncuts.cfg_min_pt_track, 1e+10f); - fDielectonCut.SetTrackEtaRange(-dielectroncuts.cfg_max_eta_track, +dielectroncuts.cfg_max_eta_track); - fDielectonCut.SetMinNClustersTPC(dielectroncuts.cfg_min_ncluster_tpc); - fDielectonCut.SetMinNCrossedRowsTPC(dielectroncuts.cfg_min_ncrossedrows); - fDielectonCut.SetMinNCrossedRowsOverFindableClustersTPC(0.8); - fDielectonCut.SetChi2PerClusterTPC(0.0, dielectroncuts.cfg_max_chi2tpc); - fDielectonCut.SetChi2PerClusterITS(0.0, dielectroncuts.cfg_max_chi2its); - fDielectonCut.SetNClustersITS(dielectroncuts.cfg_min_ncluster_its, 7); - fDielectonCut.SetMeanClusterSizeITSob(0, 16); - fDielectonCut.SetMaxDcaXY(dielectroncuts.cfg_max_dcaxy); - fDielectonCut.SetMaxDcaZ(dielectroncuts.cfg_max_dcaz); - fDielectonCut.RequireITSibAny(dielectroncuts.cfg_require_itsib_any); - fDielectonCut.RequireITSib1st(dielectroncuts.cfg_require_itsib_1st); + fDielectronCut.SetTrackPtRange(dielectroncuts.cfg_min_pt_track, 1e+10f); + fDielectronCut.SetTrackEtaRange(-dielectroncuts.cfg_max_eta_track, +dielectroncuts.cfg_max_eta_track); + fDielectronCut.SetMinNClustersTPC(dielectroncuts.cfg_min_ncluster_tpc); + fDielectronCut.SetMinNCrossedRowsTPC(dielectroncuts.cfg_min_ncrossedrows); + fDielectronCut.SetMinNCrossedRowsOverFindableClustersTPC(0.8); + fDielectronCut.SetChi2PerClusterTPC(0.0, dielectroncuts.cfg_max_chi2tpc); + fDielectronCut.SetChi2PerClusterITS(0.0, dielectroncuts.cfg_max_chi2its); + fDielectronCut.SetNClustersITS(dielectroncuts.cfg_min_ncluster_its, 7); + fDielectronCut.SetMeanClusterSizeITSob(0, 16); + fDielectronCut.SetMaxDcaXY(dielectroncuts.cfg_max_dcaxy); + fDielectronCut.SetMaxDcaZ(dielectroncuts.cfg_max_dcaz); + fDielectronCut.RequireITSibAny(dielectroncuts.cfg_require_itsib_any); + fDielectronCut.RequireITSib1st(dielectroncuts.cfg_require_itsib_1st); // for eID - fDielectonCut.SetPIDScheme(dielectroncuts.cfg_pid_scheme); - fDielectonCut.SetTPCNsigmaElRange(dielectroncuts.cfg_min_TPCNsigmaEl, dielectroncuts.cfg_max_TPCNsigmaEl); - fDielectonCut.SetTPCNsigmaMuRange(dielectroncuts.cfg_min_TPCNsigmaMu, dielectroncuts.cfg_max_TPCNsigmaMu); - fDielectonCut.SetTPCNsigmaPiRange(dielectroncuts.cfg_min_TPCNsigmaPi, dielectroncuts.cfg_max_TPCNsigmaPi); - fDielectonCut.SetTPCNsigmaKaRange(dielectroncuts.cfg_min_TPCNsigmaKa, dielectroncuts.cfg_max_TPCNsigmaKa); - fDielectonCut.SetTPCNsigmaPrRange(dielectroncuts.cfg_min_TPCNsigmaPr, dielectroncuts.cfg_max_TPCNsigmaPr); - fDielectonCut.SetTOFNsigmaElRange(dielectroncuts.cfg_min_TOFNsigmaEl, dielectroncuts.cfg_max_TOFNsigmaEl); + fDielectronCut.SetPIDScheme(dielectroncuts.cfg_pid_scheme); + fDielectronCut.SetTPCNsigmaElRange(dielectroncuts.cfg_min_TPCNsigmaEl, dielectroncuts.cfg_max_TPCNsigmaEl); + fDielectronCut.SetTPCNsigmaMuRange(dielectroncuts.cfg_min_TPCNsigmaMu, dielectroncuts.cfg_max_TPCNsigmaMu); + fDielectronCut.SetTPCNsigmaPiRange(dielectroncuts.cfg_min_TPCNsigmaPi, dielectroncuts.cfg_max_TPCNsigmaPi); + fDielectronCut.SetTPCNsigmaKaRange(dielectroncuts.cfg_min_TPCNsigmaKa, dielectroncuts.cfg_max_TPCNsigmaKa); + fDielectronCut.SetTPCNsigmaPrRange(dielectroncuts.cfg_min_TPCNsigmaPr, dielectroncuts.cfg_max_TPCNsigmaPr); + fDielectronCut.SetTOFNsigmaElRange(dielectroncuts.cfg_min_TOFNsigmaEl, dielectroncuts.cfg_max_TOFNsigmaEl); if (dielectroncuts.cfg_pid_scheme == static_cast(DielectronCut::PIDSchemes::kPIDML)) { // please call this at the end of DefineDileptonCut eid_bdt = new o2::ml::OnnxModel(); @@ -320,7 +271,7 @@ struct singleElectronQCMC { eid_bdt->initModel(dielectroncuts.BDTLocalPathGamma.value, dielectroncuts.enableOptimizations.value); } - fDielectonCut.SetPIDModel(eid_bdt); + fDielectronCut.SetPIDModel(eid_bdt); } // end of PID ML } @@ -373,7 +324,6 @@ struct singleElectronQCMC { fRegistry.fill(HIST("Track/") + HIST(ele_source_types[e_source_id]) + HIST("hPtGen_DeltaPhi"), mctrack.pt(), track.phi() - mctrack.phi()); } - std::vector used_trackIds; SliceCache cache; Preslice perCollision_track = aod::emprimaryelectron::emeventId; Filter trackFilter = dielectroncuts.cfg_min_pt_track < o2::aod::track::pt && nabs(o2::aod::track::eta) < dielectroncuts.cfg_max_eta_track && o2::aod::track::tpcChi2NCl < dielectroncuts.cfg_max_chi2tpc && o2::aod::track::itsChi2NCl < dielectroncuts.cfg_max_chi2its && nabs(o2::aod::track::dcaXY) < dielectroncuts.cfg_max_dcaxy && nabs(o2::aod::track::dcaZ) < dielectroncuts.cfg_max_dcaz; @@ -386,22 +336,20 @@ struct singleElectronQCMC { void processQCMC(FilteredMyCollisions const& collisions, FilteredMyMCTracks const& tracks, aod::EMMCParticles const& mcparticles, aod::EMMCEvents const&) { - used_trackIds.reserve(tracks.size()); for (auto& collision : collisions) { - initCCDB(collision); float centralities[3] = {collision.centFT0M(), collision.centFT0A(), collision.centFT0C()}; if (centralities[cfgCentEstimator] < cfgCentMin || cfgCentMax < centralities[cfgCentEstimator]) { continue; } - o2::aod::pwgem::dilepton::utils::eventhistogram::fillEventInfo<0, -1>(&fRegistry, collision, cfgDoFlow); + o2::aod::pwgem::dilepton::utils::eventhistogram::fillEventInfo<0, -1>(&fRegistry, collision); if (!fEMEventCut.IsSelected(collision)) { continue; } - o2::aod::pwgem::dilepton::utils::eventhistogram::fillEventInfo<1, -1>(&fRegistry, collision, cfgDoFlow); - fRegistry.fill(HIST("Event/before/hCollisionCounter"), 10.0); // accepted - fRegistry.fill(HIST("Event/after/hCollisionCounter"), 10.0); // accepted + o2::aod::pwgem::dilepton::utils::eventhistogram::fillEventInfo<1, -1>(&fRegistry, collision); + fRegistry.fill(HIST("Event/before/hCollisionCounter"), o2::aod::pwgem::dilepton::utils::eventhistogram::nbin_ev); // accepted + fRegistry.fill(HIST("Event/after/hCollisionCounter"), o2::aod::pwgem::dilepton::utils::eventhistogram::nbin_ev); // accepted auto mccollision = collision.emmcevent_as(); if (cfgEventGeneratorType >= 0 && mccollision.getSubGeneratorId() != cfgEventGeneratorType) { @@ -416,11 +364,11 @@ struct singleElectronQCMC { continue; } if (dielectroncuts.cfg_pid_scheme == static_cast(DielectronCut::PIDSchemes::kPIDML)) { - if (!fDielectonCut.IsSelectedTrack(track, collision)) { + if (!fDielectronCut.IsSelectedTrack(track, collision)) { continue; } } else { // cut-based - if (!fDielectonCut.IsSelectedTrack(track)) { + if (!fDielectronCut.IsSelectedTrack(track)) { continue; } } @@ -463,9 +411,6 @@ struct singleElectronQCMC { } // end of collision loop - used_trackIds.clear(); - used_trackIds.shrink_to_fit(); - } // end of process PROCESS_SWITCH(singleElectronQCMC, processQCMC, "run dielectron QC MC", true); diff --git a/PWGEM/Dilepton/Tasks/singleMuonQC.cxx b/PWGEM/Dilepton/Tasks/singleMuonQC.cxx new file mode 100644 index 00000000000..a4b49ee74e4 --- /dev/null +++ b/PWGEM/Dilepton/Tasks/singleMuonQC.cxx @@ -0,0 +1,252 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. +// +// ======================== +// +// Analysis task for single muon QC. +// Please write to: daiki.sekihata@cern.ch + +#include "TString.h" +#include "Math/Vector4D.h" +#include "Framework/runDataProcessing.h" +#include "Framework/AnalysisTask.h" +#include "Framework/ASoAHelpers.h" + +#include "CCDB/BasicCCDBManager.h" +#include "Common/Core/RecoDecay.h" + +#include "PWGEM/Dilepton/DataModel/dileptonTables.h" +#include "PWGEM/Dilepton/Core/DimuonCut.h" +#include "PWGEM/Dilepton/Core/EMEventCut.h" +#include "PWGEM/Dilepton/Utils/EventHistograms.h" +#include "PWGEM/Dilepton/Utils/EMTrackUtilities.h" + +using namespace o2; +using namespace o2::aod; +using namespace o2::framework; +using namespace o2::framework::expressions; +using namespace o2::soa; +using namespace o2::aod::pwgem::dilepton::utils::emtrackutil; + +using MyCollisions = soa::Join; +using MyCollision = MyCollisions::iterator; + +using MyTracks = soa::Join; +using MyTrack = MyTracks::iterator; + +struct singleMuonQC { + + // Configurables + Configurable ccdburl{"ccdb-url", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; + + Configurable cfgEventGeneratorType{"cfgEventGeneratorType", -1, "if positive, select event generator type. i.e. gap or signal"}; + Configurable cfgCentEstimator{"cfgCentEstimator", 2, "FT0M:0, FT0A:1, FT0C:2"}; + Configurable cfgCentMin{"cfgCentMin", 0, "min. centrality"}; + Configurable cfgCentMax{"cfgCentMax", 999.f, "max. centrality"}; + + ConfigurableAxis ConfPtmuBins{"ConfPtmuBins", {VARIABLE_WIDTH, 0.00, 0.05, 0.10, 0.15, 0.20, 0.30, 0.40, 0.50, 0.60, 0.70, 0.80, 0.90, 1.00, 1.10, 1.20, 1.30, 1.40, 1.50, 1.60, 1.70, 1.80, 1.90, 2.00, 2.10, 2.20, 2.30, 2.40, 2.50, 2.60, 2.70, 2.80, 2.90, 3.00, 3.10, 3.20, 3.30, 3.40, 3.50, 3.60, 3.70, 3.80, 3.90, 4.00, 4.10, 4.20, 4.30, 4.40, 4.50, 4.60, 4.70, 4.80, 4.90, 5.00, 5.50, 6.00, 6.50, 7.00, 7.50, 8.00, 8.50, 9.00, 9.50, 10.00}, "pTmu bins for output histograms"}; + + EMEventCut fEMEventCut; + struct : ConfigurableGroup { + std::string prefix = "eventcut_group"; + Configurable cfgZvtxMax{"cfgZvtxMax", 10.f, "max. Zvtx"}; + Configurable cfgRequireSel8{"cfgRequireSel8", true, "require sel8 in event cut"}; + Configurable cfgRequireFT0AND{"cfgRequireFT0AND", true, "require FT0AND in event cut"}; + Configurable cfgRequireNoTFB{"cfgRequireNoTFB", false, "require No time frame border in event cut"}; + Configurable cfgRequireNoITSROFB{"cfgRequireNoITSROFB", false, "require no ITS readout frame border in event cut"}; + Configurable cfgRequireNoSameBunchPileup{"cfgRequireNoSameBunchPileup", false, "require no same bunch pileup in event cut"}; + Configurable cfgRequireVertexITSTPC{"cfgRequireVertexITSTPC", false, "require Vertex ITSTPC in event cut"}; // ITS-TPC matched track contributes PV. + Configurable cfgRequireGoodZvtxFT0vsPV{"cfgRequireGoodZvtxFT0vsPV", false, "require good Zvtx between FT0 vs. PV in event cut"}; + Configurable cfgOccupancyMin{"cfgOccupancyMin", -1, "min. occupancy"}; + Configurable cfgOccupancyMax{"cfgOccupancyMax", 1000000000, "max. occupancy"}; + } eventcuts; + + DimuonCut fDimuonCut; + struct : ConfigurableGroup { + std::string prefix = "dimuoncut_group"; + Configurable cfg_min_mass{"cfg_min_mass", 0.0, "min mass"}; + Configurable cfg_max_mass{"cfg_max_mass", 1e+10, "max mass"}; + Configurable cfg_min_pair_pt{"cfg_min_pair_pt", 0.0, "min pair pt"}; + Configurable cfg_max_pair_pt{"cfg_max_pair_pt", 1e+10, "max pair pt"}; + Configurable cfg_min_pair_dcaxy{"cfg_min_pair_dcaxy", 0.0, "min pair dca3d in sigma"}; + Configurable cfg_max_pair_dcaxy{"cfg_max_pair_dcaxy", 1e+10, "max pair dca3d in sigma"}; + + Configurable cfg_track_type{"cfg_track_type", 3, "muon track type [0: MFT-MCH-MID, 3: MCH-MID]"}; + Configurable cfg_min_pt_track{"cfg_min_pt_track", 0.1, "min pT for single track"}; + Configurable cfg_min_eta_track{"cfg_min_eta_track", -4.0, "min eta for single track"}; + Configurable cfg_max_eta_track{"cfg_max_eta_track", -2.5, "max eta for single track"}; + Configurable cfg_min_ncluster_mft{"cfg_min_ncluster_mft", 5, "min ncluster MFT"}; + Configurable cfg_min_ncluster_mch{"cfg_min_ncluster_mch", 5, "min ncluster MCH"}; + Configurable cfg_max_chi2{"cfg_max_chi2", 1e+10, "max chi2/NclsTPC"}; + Configurable cfg_max_matching_chi2_mftmch{"cfg_max_matching_chi2_mftmch", 1e+10, "max chi2 for MFT-MCH matching"}; + Configurable cfg_max_matching_chi2_mchmid{"cfg_max_matching_chi2_mchmid", 1e+10, "max chi2 for MCH-MID matching"}; + Configurable cfg_max_dcaxy{"cfg_max_dcaxy", 1e+10, "max dca XY for single track in cm"}; + Configurable cfg_min_rabs{"cfg_min_rabs", 17.6, "min Radius at the absorber end"}; + Configurable cfg_max_rabs{"cfg_max_rabs", 89.5, "max Radius at the absorber end"}; + Configurable enableTTCA{"enableTTCA", true, "Flag to enable or disable TTCA"}; + } dimuoncuts; + + o2::ccdb::CcdbApi ccdbApi; + Service ccdb; + + HistogramRegistry fRegistry{"output", {}, OutputObjHandlingPolicy::AnalysisObject, false, false}; + static constexpr std::string_view event_cut_types[2] = {"before/", "after/"}; + + ~singleMuonQC() {} + + void addhistograms() + { + // event info + o2::aod::pwgem::dilepton::utils::eventhistogram::addEventHistograms(&fRegistry); + + const AxisSpec axis_pt{ConfPtmuBins, "p_{T,#mu} (GeV/c)"}; + const AxisSpec axis_eta{25, -4.5, -2.0, "#eta_{#mu}"}; + const AxisSpec axis_phi{36, 0.0, 2 * M_PI, "#varphi_{#mu} (rad.)"}; + const AxisSpec axis_charge_rec{3, -1.5, +1.5, "charge"}; + + // track info + fRegistry.add("Track/hs", "rec. single muon", kTHnSparseF, {axis_pt, axis_eta, axis_phi, axis_charge_rec}, true); + fRegistry.add("Track/hQoverPt", "q/pT;q/p_{T} (GeV/c)^{-1}", kTH1F, {{400, -20, 20}}, false); + fRegistry.add("Track/hTrackType", "track type", kTH1F, {{6, -0.5f, 5.5}}, false); + fRegistry.add("Track/hDCAxy", "DCA x vs. y;DCA_{x} (cm);DCA_{y} (cm)", kTH2F, {{200, -1.0f, 1.0f}, {200, -1.0f, 1.0f}}, false); + fRegistry.add("Track/hDCAxySigma", "DCA x vs. y;DCA_{x} (#sigma);DCA_{y} (#sigma)", kTH2F, {{200, -10.0f, 10.0f}, {200, -10.0f, 10.0f}}, false); + fRegistry.add("Track/hDCA2DSigma", "DCA xy;DCA_{xy} (#sigma)", kTH1F, {{100, 0.0f, 10.0f}}, false); + fRegistry.add("Track/hDCAxRes_Pt", "DCA_{x} resolution vs. pT;p_{T} (GeV/c);DCA_{x} resolution (#mum)", kTH2F, {{1000, 0, 10}, {200, 0., 200}}, false); + fRegistry.add("Track/hDCAyRes_Pt", "DCA_{y} resolution vs. pT;p_{T} (GeV/c);DCA_{y} resolution (#mum)", kTH2F, {{1000, 0, 10}, {200, 0., 200}}, false); + fRegistry.add("Track/hNclsMCH", "number of MCH clusters", kTH1F, {{21, -0.5, 20.5}}, false); + fRegistry.add("Track/hNclsMFT", "number of MFT clusters", kTH1F, {{11, -0.5, 10.5}}, false); + fRegistry.add("Track/hPDCA", "pDCA;p_{T} at PV (GeV/c);p #times DCA (GeV/c #upoint cm)", kTH2F, {{100, 0, 10}, {100, 0.0f, 1000}}, false); + fRegistry.add("Track/hChi2", "chi2;chi2", kTH1F, {{100, 0.0f, 100}}, false); + fRegistry.add("Track/hChi2MatchMCHMID", "chi2 match MCH-MID;chi2", kTH1F, {{100, 0.0f, 100}}, false); + fRegistry.add("Track/hChi2MatchMCHMFT", "chi2 match MCH-MFT;chi2", kTH1F, {{100, 0.0f, 100}}, false); + fRegistry.add("Track/hMFTClusterMap", "MFT cluster map", kTH1F, {{1024, -0.5, 1023.5}}, false); + } + + void init(InitContext&) + { + ccdb->setURL(ccdburl); + ccdb->setCaching(true); + ccdb->setLocalObjectValidityChecking(); + ccdb->setFatalWhenNull(false); + + DefineEMEventCut(); + DefineDimuonCut(); + addhistograms(); + } + + void DefineEMEventCut() + { + fEMEventCut = EMEventCut("fEMEventCut", "fEMEventCut"); + fEMEventCut.SetRequireSel8(eventcuts.cfgRequireSel8); + fEMEventCut.SetRequireFT0AND(eventcuts.cfgRequireFT0AND); + fEMEventCut.SetZvtxRange(-eventcuts.cfgZvtxMax, +eventcuts.cfgZvtxMax); + fEMEventCut.SetRequireNoTFB(eventcuts.cfgRequireNoTFB); + fEMEventCut.SetRequireNoITSROFB(eventcuts.cfgRequireNoITSROFB); + fEMEventCut.SetRequireNoSameBunchPileup(eventcuts.cfgRequireNoSameBunchPileup); + fEMEventCut.SetRequireVertexITSTPC(eventcuts.cfgRequireVertexITSTPC); + fEMEventCut.SetRequireGoodZvtxFT0vsPV(eventcuts.cfgRequireGoodZvtxFT0vsPV); + fEMEventCut.SetOccupancyRange(eventcuts.cfgOccupancyMin, eventcuts.cfgOccupancyMax); + } + + void DefineDimuonCut() + { + fDimuonCut = DimuonCut("fDimuonCut", "fDimuonCut"); + + // for pair + fDimuonCut.SetMassRange(dimuoncuts.cfg_min_mass, dimuoncuts.cfg_max_mass); + fDimuonCut.SetPairPtRange(dimuoncuts.cfg_min_pair_pt, dimuoncuts.cfg_max_pair_pt); + fDimuonCut.SetPairDCAxyRange(dimuoncuts.cfg_min_pair_dcaxy, dimuoncuts.cfg_max_pair_dcaxy); // DCAxy in cm + + // for track + fDimuonCut.SetTrackType(dimuoncuts.cfg_track_type); + fDimuonCut.SetTrackPtRange(dimuoncuts.cfg_min_pt_track, 1e10f); + fDimuonCut.SetTrackEtaRange(dimuoncuts.cfg_min_eta_track, dimuoncuts.cfg_max_eta_track); + fDimuonCut.SetNClustersMFT(dimuoncuts.cfg_min_ncluster_mft, 10); + fDimuonCut.SetNClustersMCHMID(dimuoncuts.cfg_min_ncluster_mch, 16); + fDimuonCut.SetChi2(0.f, dimuoncuts.cfg_max_chi2); + fDimuonCut.SetMatchingChi2MCHMFT(0.f, dimuoncuts.cfg_max_matching_chi2_mftmch); + fDimuonCut.SetMatchingChi2MCHMID(0.f, dimuoncuts.cfg_max_matching_chi2_mchmid); + fDimuonCut.SetDCAxy(0.f, dimuoncuts.cfg_max_dcaxy); + fDimuonCut.SetRabs(dimuoncuts.cfg_min_rabs, dimuoncuts.cfg_max_rabs); + fDimuonCut.SetMaxPDCARabsDep([&](float rabs) { return (rabs < 26.5 ? 594.f : 324.f); }); + } + + template + void fillTrackInfo(TTrack const& track) + { + float dca_xy = fwdDcaXYinSigma(track); + fRegistry.fill(HIST("Track/hs"), track.pt(), track.eta(), track.phi(), track.sign()); + fRegistry.fill(HIST("Track/hQoverPt"), track.sign() / track.pt()); + fRegistry.fill(HIST("Track/hTrackType"), track.trackType()); + fRegistry.fill(HIST("Track/hDCAxy"), track.fwdDcaX(), track.fwdDcaY()); + fRegistry.fill(HIST("Track/hDCAxySigma"), track.fwdDcaX() / std::sqrt(track.cXX()), track.fwdDcaY() / std::sqrt(track.cYY())); + fRegistry.fill(HIST("Track/hDCA2DSigma"), dca_xy); + fRegistry.fill(HIST("Track/hDCAxRes_Pt"), track.pt(), std::sqrt(track.cXX()) * 1e+4); + fRegistry.fill(HIST("Track/hDCAyRes_Pt"), track.pt(), std::sqrt(track.cYY()) * 1e+4); + fRegistry.fill(HIST("Track/hNclsMCH"), track.nClusters()); + fRegistry.fill(HIST("Track/hNclsMFT"), track.nClustersMFT()); + fRegistry.fill(HIST("Track/hPDCA"), track.pt(), track.pDca()); + fRegistry.fill(HIST("Track/hChi2"), track.chi2()); + fRegistry.fill(HIST("Track/hChi2MatchMCHMID"), track.chi2MatchMCHMID()); + fRegistry.fill(HIST("Track/hChi2MatchMCHMFT"), track.chi2MatchMCHMFT()); + fRegistry.fill(HIST("Track/hMFTClusterMap"), track.mftClusterMap()); + } + + SliceCache cache; + Preslice perCollision_track = aod::emprimarymuon::emeventId; + Filter trackFilter = o2::aod::fwdtrack::trackType == dimuoncuts.cfg_track_type && dimuoncuts.cfg_min_pt_track < o2::aod::fwdtrack::pt && dimuoncuts.cfg_min_eta_track < o2::aod::fwdtrack::eta && o2::aod::fwdtrack::eta < dimuoncuts.cfg_max_eta_track; + Filter ttcaFilter = ifnode(dimuoncuts.enableTTCA.node(), o2::aod::emprimarymuon::isAssociatedToMPC == true || o2::aod::emprimarymuon::isAssociatedToMPC == false, o2::aod::emprimarymuon::isAssociatedToMPC == true); + using FilteredMyTracks = soa::Filtered; + + Filter collisionFilter_centrality = (cfgCentMin < o2::aod::cent::centFT0M && o2::aod::cent::centFT0M < cfgCentMax) || (cfgCentMin < o2::aod::cent::centFT0A && o2::aod::cent::centFT0A < cfgCentMax) || (cfgCentMin < o2::aod::cent::centFT0C && o2::aod::cent::centFT0C < cfgCentMax); + using FilteredMyCollisions = soa::Filtered; + + void processQC(FilteredMyCollisions const& collisions, FilteredMyTracks const& tracks) + { + for (auto& collision : collisions) { + float centralities[3] = {collision.centFT0M(), collision.centFT0A(), collision.centFT0C()}; + if (centralities[cfgCentEstimator] < cfgCentMin || cfgCentMax < centralities[cfgCentEstimator]) { + continue; + } + + o2::aod::pwgem::dilepton::utils::eventhistogram::fillEventInfo<0, -1>(&fRegistry, collision); + if (!fEMEventCut.IsSelected(collision)) { + continue; + } + o2::aod::pwgem::dilepton::utils::eventhistogram::fillEventInfo<1, -1>(&fRegistry, collision); + fRegistry.fill(HIST("Event/before/hCollisionCounter"), o2::aod::pwgem::dilepton::utils::eventhistogram::nbin_ev); // accepted + fRegistry.fill(HIST("Event/after/hCollisionCounter"), o2::aod::pwgem::dilepton::utils::eventhistogram::nbin_ev); // accepted + + auto tracks_per_coll = tracks.sliceBy(perCollision_track, collision.globalIndex()); + + for (auto& track : tracks_per_coll) { + if (!fDimuonCut.IsSelectedTrack(track)) { + continue; + } + + fillTrackInfo(track); + + } // end of track loop + + } // end of collision loop + + } // end of process + PROCESS_SWITCH(singleMuonQC, processQC, "run QC", true); + + void processDummy(MyCollisions const&) {} + PROCESS_SWITCH(singleMuonQC, processDummy, "Dummy function", false); +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{ + adaptAnalysisTask(cfgc, TaskName{"single-muon-qc"})}; +} diff --git a/PWGEM/Dilepton/Tasks/singleMuonQCMC.cxx b/PWGEM/Dilepton/Tasks/singleMuonQCMC.cxx new file mode 100644 index 00000000000..54c30041a75 --- /dev/null +++ b/PWGEM/Dilepton/Tasks/singleMuonQCMC.cxx @@ -0,0 +1,409 @@ +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. +// All rights not expressly granted are reserved. +// +// This software is distributed under the terms of the GNU General Public +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". +// +// In applying this license CERN does not waive the privileges and immunities +// granted to it by virtue of its status as an Intergovernmental Organization +// or submit itself to any jurisdiction. +// +// ======================== +// +// Analysis task for dimuon in MC. +// Please write to: daiki.sekihata@cern.ch + +#include "TString.h" +#include "Math/Vector4D.h" +#include "Framework/runDataProcessing.h" +#include "Framework/AnalysisTask.h" +#include "Framework/ASoAHelpers.h" + +#include "CCDB/BasicCCDBManager.h" +#include "Common/Core/RecoDecay.h" + +#include "PWGEM/Dilepton/DataModel/dileptonTables.h" +#include "PWGEM/Dilepton/Core/DimuonCut.h" +#include "PWGEM/Dilepton/Core/EMEventCut.h" +#include "PWGEM/Dilepton/Utils/MCUtilities.h" +#include "PWGEM/Dilepton/Utils/EventHistograms.h" +#include "PWGEM/Dilepton/Utils/EMTrackUtilities.h" +#include "PWGEM/Dilepton/Utils/PairUtilities.h" + +using namespace o2; +using namespace o2::aod; +using namespace o2::framework; +using namespace o2::framework::expressions; +using namespace o2::soa; +using namespace o2::aod::pwgem::dilepton::utils::mcutil; +using namespace o2::aod::pwgem::dilepton::utils::emtrackutil; + +using MyCollisions = soa::Join; +using MyCollision = MyCollisions::iterator; + +using MyMCTracks = soa::Join; +using MyMCTrack = MyMCTracks::iterator; + +struct singleMuonQCMC { + + // Configurables + Configurable ccdburl{"ccdb-url", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; + + Configurable cfgEventGeneratorType{"cfgEventGeneratorType", -1, "if positive, select event generator type. i.e. gap or signal"}; + Configurable cfgCentEstimator{"cfgCentEstimator", 2, "FT0M:0, FT0A:1, FT0C:2"}; + Configurable cfgCentMin{"cfgCentMin", 0, "min. centrality"}; + Configurable cfgCentMax{"cfgCentMax", 999.f, "max. centrality"}; + + ConfigurableAxis ConfPtmuBins{"ConfPtmuBins", {VARIABLE_WIDTH, 0.00, 0.05, 0.10, 0.15, 0.20, 0.30, 0.40, 0.50, 0.60, 0.70, 0.80, 0.90, 1.00, 1.10, 1.20, 1.30, 1.40, 1.50, 1.60, 1.70, 1.80, 1.90, 2.00, 2.10, 2.20, 2.30, 2.40, 2.50, 2.60, 2.70, 2.80, 2.90, 3.00, 3.10, 3.20, 3.30, 3.40, 3.50, 3.60, 3.70, 3.80, 3.90, 4.00, 4.10, 4.20, 4.30, 4.40, 4.50, 4.60, 4.70, 4.80, 4.90, 5.00, 5.50, 6.00, 6.50, 7.00, 7.50, 8.00, 8.50, 9.00, 9.50, 10.00}, "pTmu bins for output histograms"}; + + EMEventCut fEMEventCut; + struct : ConfigurableGroup { + std::string prefix = "eventcut_group"; + Configurable cfgZvtxMax{"cfgZvtxMax", 10.f, "max. Zvtx"}; + Configurable cfgRequireSel8{"cfgRequireSel8", true, "require sel8 in event cut"}; + Configurable cfgRequireFT0AND{"cfgRequireFT0AND", true, "require FT0AND in event cut"}; + Configurable cfgRequireNoTFB{"cfgRequireNoTFB", false, "require No time frame border in event cut"}; + Configurable cfgRequireNoITSROFB{"cfgRequireNoITSROFB", false, "require no ITS readout frame border in event cut"}; + Configurable cfgRequireNoSameBunchPileup{"cfgRequireNoSameBunchPileup", false, "require no same bunch pileup in event cut"}; + Configurable cfgRequireVertexITSTPC{"cfgRequireVertexITSTPC", false, "require Vertex ITSTPC in event cut"}; // ITS-TPC matched track contributes PV. + Configurable cfgRequireGoodZvtxFT0vsPV{"cfgRequireGoodZvtxFT0vsPV", false, "require good Zvtx between FT0 vs. PV in event cut"}; + Configurable cfgOccupancyMin{"cfgOccupancyMin", -1, "min. occupancy"}; + Configurable cfgOccupancyMax{"cfgOccupancyMax", 1000000000, "max. occupancy"}; + } eventcuts; + + DimuonCut fDimuonCut; + struct : ConfigurableGroup { + std::string prefix = "dimuoncut_group"; + Configurable cfg_min_mass{"cfg_min_mass", 0.0, "min mass"}; + Configurable cfg_max_mass{"cfg_max_mass", 1e+10, "max mass"}; + Configurable cfg_min_pair_pt{"cfg_min_pair_pt", 0.0, "min pair pt"}; + Configurable cfg_max_pair_pt{"cfg_max_pair_pt", 1e+10, "max pair pt"}; + Configurable cfg_min_pair_dcaxy{"cfg_min_pair_dcaxy", 0.0, "min pair dca3d in sigma"}; + Configurable cfg_max_pair_dcaxy{"cfg_max_pair_dcaxy", 1e+10, "max pair dca3d in sigma"}; + + Configurable cfg_track_type{"cfg_track_type", 3, "muon track type [0: MFT-MCH-MID, 3: MCH-MID]"}; + Configurable cfg_min_pt_track{"cfg_min_pt_track", 0.1, "min pT for single track"}; + Configurable cfg_min_eta_track{"cfg_min_eta_track", -4.0, "min eta for single track"}; + Configurable cfg_max_eta_track{"cfg_max_eta_track", -2.5, "max eta for single track"}; + Configurable cfg_min_ncluster_mft{"cfg_min_ncluster_mft", 5, "min ncluster MFT"}; + Configurable cfg_min_ncluster_mch{"cfg_min_ncluster_mch", 5, "min ncluster MCH"}; + Configurable cfg_max_chi2{"cfg_max_chi2", 1e+10, "max chi2/NclsTPC"}; + Configurable cfg_max_matching_chi2_mftmch{"cfg_max_matching_chi2_mftmch", 1e+10, "max chi2 for MFT-MCH matching"}; + Configurable cfg_max_matching_chi2_mchmid{"cfg_max_matching_chi2_mchmid", 1e+10, "max chi2 for MCH-MID matching"}; + Configurable cfg_max_dcaxy{"cfg_max_dcaxy", 1e+10, "max dca XY for single track in cm"}; + Configurable cfg_min_rabs{"cfg_min_rabs", 17.6, "min Radius at the absorber end"}; + Configurable cfg_max_rabs{"cfg_max_rabs", 89.5, "max Radius at the absorber end"}; + Configurable enableTTCA{"enableTTCA", true, "Flag to enable or disable TTCA"}; + } dimuoncuts; + + o2::ccdb::CcdbApi ccdbApi; + Service ccdb; + + struct : ConfigurableGroup { + std::string prefix = "mctrackcut_group"; + Configurable min_mcPt{"min_mcPt", 0.05, "min. MC pT"}; + Configurable max_mcPt{"max_mcPt", 1e+10, "max. MC pT"}; + Configurable min_mcEta{"min_mcEta", -4.0, "min. MC eta"}; + Configurable max_mcEta{"max_mcEta", -2.5, "max. MC eta"}; + } mctrackcuts; + + HistogramRegistry fRegistry{"output", {}, OutputObjHandlingPolicy::AnalysisObject, false, false}; + static constexpr std::string_view event_cut_types[2] = {"before/", "after/"}; + static constexpr std::string_view muon_source_types[9] = {"lf/", "Photon/", "PromptJPsi/", "NonPromptJPsi/", "PromptPsi2S/", "NonPromptPsi2S/", "c2mu/", "b2mu/", "b2c2mu/"}; + + ~singleMuonQCMC() {} + + void addhistograms() + { + // event info + o2::aod::pwgem::dilepton::utils::eventhistogram::addEventHistograms(&fRegistry); + + const AxisSpec axis_pt{ConfPtmuBins, "p_{T,#mu} (GeV/c)"}; + const AxisSpec axis_eta{25, -4.5, -2.0, "#eta_{#mu}"}; + const AxisSpec axis_phi{36, 0.0, 2 * M_PI, "#varphi_{#mu} (rad.)"}; + const AxisSpec axis_charge_rec{3, -1.5, +1.5, "charge"}; + const AxisSpec axis_charge_gen{3, -1.5, +1.5, "true charge"}; + + // generated info + fRegistry.add("Generated/lf/hs", "gen. single muon", kTHnSparseF, {axis_pt, axis_eta, axis_phi, axis_charge_gen}, true); + fRegistry.addClone("Generated/lf/", "Generated/PromptJPsi/"); + fRegistry.addClone("Generated/lf/", "Generated/NonPromptJPsi/"); + fRegistry.addClone("Generated/lf/", "Generated/PromptPsi2S/"); + fRegistry.addClone("Generated/lf/", "Generated/NonPromptPsi2S/"); + fRegistry.addClone("Generated/lf/", "Generated/c2mu/"); + fRegistry.addClone("Generated/lf/", "Generated/b2mu/"); + fRegistry.addClone("Generated/lf/", "Generated/b2c2mu/"); + + // track info + fRegistry.add("Track/lf/hs", "rec. single muon", kTHnSparseF, {axis_pt, axis_eta, axis_phi, axis_charge_rec, axis_charge_gen}, true); + fRegistry.add("Track/lf/hQoverPt", "q/pT;q/p_{T} (GeV/c)^{-1}", kTH1F, {{400, -20, 20}}, false); + fRegistry.add("Track/lf/hTrackType", "track type", kTH1F, {{6, -0.5f, 5.5}}, false); + fRegistry.add("Track/lf/hDCAxy", "DCA x vs. y;DCA_{x} (cm);DCA_{y} (cm)", kTH2F, {{200, -1.0f, 1.0f}, {200, -1.0f, 1.0f}}, false); + fRegistry.add("Track/lf/hDCAxySigma", "DCA x vs. y;DCA_{x} (#sigma);DCA_{y} (#sigma)", kTH2F, {{200, -10.0f, 10.0f}, {200, -10.0f, 10.0f}}, false); + fRegistry.add("Track/lf/hDCA2DSigma", "DCA xy;DCA_{xy} (#sigma)", kTH1F, {{100, 0.0f, 10.0f}}, false); + fRegistry.add("Track/lf/hDCAxRes_Pt", "DCA_{x} resolution vs. pT;p_{T} (GeV/c);DCA_{x} resolution (#mum)", kTH2F, {{1000, 0, 10}, {200, 0., 200}}, false); + fRegistry.add("Track/lf/hDCAyRes_Pt", "DCA_{y} resolution vs. pT;p_{T} (GeV/c);DCA_{y} resolution (#mum)", kTH2F, {{1000, 0, 10}, {200, 0., 200}}, false); + fRegistry.add("Track/lf/hNclsMCH", "number of MCH clusters", kTH1F, {{21, -0.5, 20.5}}, false); + fRegistry.add("Track/lf/hNclsMFT", "number of MFT clusters", kTH1F, {{11, -0.5, 10.5}}, false); + fRegistry.add("Track/lf/hPDCA", "pDCA;p_{T} at PV (GeV/c);p #times DCA (GeV/c #upoint cm)", kTH2F, {{100, 0, 10}, {100, 0.0f, 1000}}, false); + fRegistry.add("Track/lf/hChi2", "chi2;chi2", kTH1F, {{100, 0.0f, 100}}, false); + fRegistry.add("Track/lf/hChi2MatchMCHMID", "chi2 match MCH-MID;chi2", kTH1F, {{100, 0.0f, 100}}, false); + fRegistry.add("Track/lf/hChi2MatchMCHMFT", "chi2 match MCH-MFT;chi2", kTH1F, {{100, 0.0f, 100}}, false); + fRegistry.add("Track/lf/hMFTClusterMap", "MFT cluster map", kTH1F, {{1024, -0.5, 1023.5}}, false); + fRegistry.add("Track/lf/hPtGen_DeltaPtOverPtGen", "muon p_{T} resolution;p_{T}^{gen} (GeV/c);(p_{T}^{rec} - p_{T}^{gen})/p_{T}^{gen}", kTH2F, {{1000, 0, 10}, {400, -1.0f, 1.0f}}, true); + fRegistry.add("Track/lf/hPtGen_DeltaEta", "muon #eta resolution;p_{T}^{gen} (GeV/c);#eta^{rec} - #eta^{gen}", kTH2F, {{1000, 0, 10}, {400, -1.0f, 1.0f}}, true); + fRegistry.add("Track/lf/hPtGen_DeltaPhi", "muon #varphi resolution;p_{T}^{gen} (GeV/c);#varphi^{rec} - #varphi^{gen} (rad.)", kTH2F, {{1000, 0, 10}, {400, -1.0f, 1.0f}}, true); + fRegistry.addClone("Track/lf/", "Track/Photon/"); // this is not for efficiency! only for contamination. We don't store generated photon conversions. + fRegistry.addClone("Track/lf/", "Track/PromptJPsi/"); + fRegistry.addClone("Track/lf/", "Track/NonPromptJPsi/"); + fRegistry.addClone("Track/lf/", "Track/PromptPsi2S/"); + fRegistry.addClone("Track/lf/", "Track/NonPromptPsi2S/"); + fRegistry.addClone("Track/lf/", "Track/c2mu/"); + fRegistry.addClone("Track/lf/", "Track/b2mu/"); + fRegistry.addClone("Track/lf/", "Track/b2c2mu/"); + } + + void init(InitContext&) + { + ccdb->setURL(ccdburl); + ccdb->setCaching(true); + ccdb->setLocalObjectValidityChecking(); + ccdb->setFatalWhenNull(false); + + DefineEMEventCut(); + DefineDimuonCut(); + addhistograms(); + } + + void DefineEMEventCut() + { + fEMEventCut = EMEventCut("fEMEventCut", "fEMEventCut"); + fEMEventCut.SetRequireSel8(eventcuts.cfgRequireSel8); + fEMEventCut.SetRequireFT0AND(eventcuts.cfgRequireFT0AND); + fEMEventCut.SetZvtxRange(-eventcuts.cfgZvtxMax, +eventcuts.cfgZvtxMax); + fEMEventCut.SetRequireNoTFB(eventcuts.cfgRequireNoTFB); + fEMEventCut.SetRequireNoITSROFB(eventcuts.cfgRequireNoITSROFB); + fEMEventCut.SetRequireNoSameBunchPileup(eventcuts.cfgRequireNoSameBunchPileup); + fEMEventCut.SetRequireVertexITSTPC(eventcuts.cfgRequireVertexITSTPC); + fEMEventCut.SetRequireGoodZvtxFT0vsPV(eventcuts.cfgRequireGoodZvtxFT0vsPV); + fEMEventCut.SetOccupancyRange(eventcuts.cfgOccupancyMin, eventcuts.cfgOccupancyMax); + } + + void DefineDimuonCut() + { + fDimuonCut = DimuonCut("fDimuonCut", "fDimuonCut"); + + // for pair + fDimuonCut.SetMassRange(dimuoncuts.cfg_min_mass, dimuoncuts.cfg_max_mass); + fDimuonCut.SetPairPtRange(dimuoncuts.cfg_min_pair_pt, dimuoncuts.cfg_max_pair_pt); + fDimuonCut.SetPairDCAxyRange(dimuoncuts.cfg_min_pair_dcaxy, dimuoncuts.cfg_max_pair_dcaxy); // DCAxy in cm + + // for track + fDimuonCut.SetTrackType(dimuoncuts.cfg_track_type); + fDimuonCut.SetTrackPtRange(dimuoncuts.cfg_min_pt_track, 1e10f); + fDimuonCut.SetTrackEtaRange(dimuoncuts.cfg_min_eta_track, dimuoncuts.cfg_max_eta_track); + fDimuonCut.SetNClustersMFT(dimuoncuts.cfg_min_ncluster_mft, 10); + fDimuonCut.SetNClustersMCHMID(dimuoncuts.cfg_min_ncluster_mch, 16); + fDimuonCut.SetChi2(0.f, dimuoncuts.cfg_max_chi2); + fDimuonCut.SetMatchingChi2MCHMFT(0.f, dimuoncuts.cfg_max_matching_chi2_mftmch); + fDimuonCut.SetMatchingChi2MCHMID(0.f, dimuoncuts.cfg_max_matching_chi2_mchmid); + fDimuonCut.SetDCAxy(0.f, dimuoncuts.cfg_max_dcaxy); + fDimuonCut.SetRabs(dimuoncuts.cfg_min_rabs, dimuoncuts.cfg_max_rabs); + fDimuonCut.SetMaxPDCARabsDep([&](float rabs) { return (rabs < 26.5 ? 594.f : 324.f); }); + } + + template + bool isInAcceptance(T const& t1) + { + if ((mctrackcuts.min_mcPt < t1.pt() && t1.pt() < mctrackcuts.max_mcPt) && (mctrackcuts.min_mcEta < t1.eta() && t1.eta() < mctrackcuts.max_mcEta)) { + return true; + } else { + return false; + } + } + + template + void fillTrackInfo(TTrack const& track) + { + auto mctrack = track.template emmcparticle_as(); + float dca_xy = fwdDcaXYinSigma(track); + fRegistry.fill(HIST("Track/") + HIST(muon_source_types[mu_source_id]) + HIST("hs"), track.pt(), track.eta(), track.phi(), track.sign(), -mctrack.pdgCode() / 11); + fRegistry.fill(HIST("Track/") + HIST(muon_source_types[mu_source_id]) + HIST("hQoverPt"), track.sign() / track.pt()); + fRegistry.fill(HIST("Track/") + HIST(muon_source_types[mu_source_id]) + HIST("hTrackType"), track.trackType()); + fRegistry.fill(HIST("Track/") + HIST(muon_source_types[mu_source_id]) + HIST("hDCAxy"), track.fwdDcaX(), track.fwdDcaY()); + fRegistry.fill(HIST("Track/") + HIST(muon_source_types[mu_source_id]) + HIST("hDCAxySigma"), track.fwdDcaX() / std::sqrt(track.cXX()), track.fwdDcaY() / std::sqrt(track.cYY())); + fRegistry.fill(HIST("Track/") + HIST(muon_source_types[mu_source_id]) + HIST("hDCA2DSigma"), dca_xy); + fRegistry.fill(HIST("Track/") + HIST(muon_source_types[mu_source_id]) + HIST("hDCAxRes_Pt"), track.pt(), std::sqrt(track.cXX()) * 1e+4); + fRegistry.fill(HIST("Track/") + HIST(muon_source_types[mu_source_id]) + HIST("hDCAyRes_Pt"), track.pt(), std::sqrt(track.cYY()) * 1e+4); + fRegistry.fill(HIST("Track/") + HIST(muon_source_types[mu_source_id]) + HIST("hNclsMCH"), track.nClusters()); + fRegistry.fill(HIST("Track/") + HIST(muon_source_types[mu_source_id]) + HIST("hNclsMFT"), track.nClustersMFT()); + fRegistry.fill(HIST("Track/") + HIST(muon_source_types[mu_source_id]) + HIST("hPDCA"), track.pt(), track.pDca()); + fRegistry.fill(HIST("Track/") + HIST(muon_source_types[mu_source_id]) + HIST("hChi2"), track.chi2()); + fRegistry.fill(HIST("Track/") + HIST(muon_source_types[mu_source_id]) + HIST("hChi2MatchMCHMID"), track.chi2MatchMCHMID()); + fRegistry.fill(HIST("Track/") + HIST(muon_source_types[mu_source_id]) + HIST("hChi2MatchMCHMFT"), track.chi2MatchMCHMFT()); + fRegistry.fill(HIST("Track/") + HIST(muon_source_types[mu_source_id]) + HIST("hMFTClusterMap"), track.mftClusterMap()); + fRegistry.fill(HIST("Track/") + HIST(muon_source_types[mu_source_id]) + HIST("hPtGen_DeltaPtOverPtGen"), mctrack.pt(), (track.pt() - mctrack.pt()) / mctrack.pt()); + fRegistry.fill(HIST("Track/") + HIST(muon_source_types[mu_source_id]) + HIST("hPtGen_DeltaEta"), mctrack.pt(), track.eta() - mctrack.eta()); + fRegistry.fill(HIST("Track/") + HIST(muon_source_types[mu_source_id]) + HIST("hPtGen_DeltaPhi"), mctrack.pt(), track.phi() - mctrack.phi()); + } + + SliceCache cache; + Preslice perCollision_track = aod::emprimarymuon::emeventId; + Filter trackFilter = o2::aod::fwdtrack::trackType == dimuoncuts.cfg_track_type && dimuoncuts.cfg_min_pt_track < o2::aod::fwdtrack::pt && dimuoncuts.cfg_min_eta_track < o2::aod::fwdtrack::eta && o2::aod::fwdtrack::eta < dimuoncuts.cfg_max_eta_track; + Filter ttcaFilter = ifnode(dimuoncuts.enableTTCA.node(), o2::aod::emprimarymuon::isAssociatedToMPC == true || o2::aod::emprimarymuon::isAssociatedToMPC == false, o2::aod::emprimarymuon::isAssociatedToMPC == true); + using FilteredMyMCTracks = soa::Filtered; + + Filter collisionFilter_centrality = (cfgCentMin < o2::aod::cent::centFT0M && o2::aod::cent::centFT0M < cfgCentMax) || (cfgCentMin < o2::aod::cent::centFT0A && o2::aod::cent::centFT0A < cfgCentMax) || (cfgCentMin < o2::aod::cent::centFT0C && o2::aod::cent::centFT0C < cfgCentMax); + using FilteredMyCollisions = soa::Filtered; + + void processQCMC(FilteredMyCollisions const& collisions, FilteredMyMCTracks const& tracks, aod::EMMCParticles const& mcparticles, aod::EMMCEvents const&) + { + for (auto& collision : collisions) { + float centralities[3] = {collision.centFT0M(), collision.centFT0A(), collision.centFT0C()}; + if (centralities[cfgCentEstimator] < cfgCentMin || cfgCentMax < centralities[cfgCentEstimator]) { + continue; + } + + o2::aod::pwgem::dilepton::utils::eventhistogram::fillEventInfo<0, -1>(&fRegistry, collision); + if (!fEMEventCut.IsSelected(collision)) { + continue; + } + o2::aod::pwgem::dilepton::utils::eventhistogram::fillEventInfo<1, -1>(&fRegistry, collision); + fRegistry.fill(HIST("Event/before/hCollisionCounter"), o2::aod::pwgem::dilepton::utils::eventhistogram::nbin_ev); // accepted + fRegistry.fill(HIST("Event/after/hCollisionCounter"), o2::aod::pwgem::dilepton::utils::eventhistogram::nbin_ev); // accepted + + auto mccollision = collision.emmcevent_as(); + if (cfgEventGeneratorType >= 0 && mccollision.getSubGeneratorId() != cfgEventGeneratorType) { + continue; + } + + auto tracks_per_coll = tracks.sliceBy(perCollision_track, collision.globalIndex()); + + for (auto& track : tracks_per_coll) { + auto mctrack = track.template emmcparticle_as(); + if (abs(mctrack.pdgCode()) != 13) { + continue; + } + if (!fDimuonCut.IsSelectedTrack(track)) { + continue; + } + + if (!mctrack.has_mothers()) { + continue; + } + auto mcmother = mcparticles.iteratorAt(mctrack.mothersIds()[0]); + int pdg_mother = abs(mcmother.pdgCode()); + + if (mctrack.isPhysicalPrimary() || mctrack.producedByGenerator()) { + if (pdg_mother == 111 || pdg_mother == 221 || pdg_mother == 331 || pdg_mother == 113 || pdg_mother == 223 || pdg_mother == 333) { + fillTrackInfo<0, aod::EMMCParticles>(track); + } else if (pdg_mother == 443) { + if (IsFromBeauty(mcmother, mcparticles) > 0) { // b is found in full decay chain. + fillTrackInfo<3, aod::EMMCParticles>(track); + } else { + fillTrackInfo<2, aod::EMMCParticles>(track); + } + } else if (pdg_mother == 100443) { + if (IsFromBeauty(mcmother, mcparticles) > 0) { // b is found in full decay chain. + fillTrackInfo<5, aod::EMMCParticles>(track); + } else { + fillTrackInfo<4, aod::EMMCParticles>(track); + } + } else if (IsFromBeauty(mctrack, mcparticles) > 0) { // b is found in full decay chain. + if (IsFromCharm(mctrack, mcparticles) > 0) { // c is found in full decay chain. + fillTrackInfo<8, aod::EMMCParticles>(track); + } else { + fillTrackInfo<7, aod::EMMCParticles>(track); + } + } else if (IsFromCharm(mctrack, mcparticles) > 0) { // c is found in full decay chain. Not from b. + fillTrackInfo<6, aod::EMMCParticles>(track); + } + } else { + fillTrackInfo<1, aod::EMMCParticles>(track); + } + + } // end of track loop + + } // end of collision loop + + } // end of process + PROCESS_SWITCH(singleMuonQCMC, processQCMC, "run QC MC", true); + + Partition muonsMC = nabs(o2::aod::mcparticle::pdgCode) == 13; // mu+, mu- + PresliceUnsorted perMcCollision = aod::emmcparticle::emmceventId; + void processGen(MyCollisions const& collisions, aod::EMMCEvents const&, aod::EMMCParticles const& mcparticles) + { + // loop over mc stack and fill histograms for pure MC truth signals + // all MC tracks which belong to the MC event corresponding to the current reconstructed event + + for (auto& collision : collisions) { + float centralities[3] = {collision.centFT0M(), collision.centFT0A(), collision.centFT0C()}; + if (centralities[cfgCentEstimator] < cfgCentMin || cfgCentMax < centralities[cfgCentEstimator]) { + continue; + } + + if (!fEMEventCut.IsSelected(collision)) { + continue; + } + + auto mccollision = collision.emmcevent_as(); + // LOGF(info, "mccollision.getGeneratorId() = %d", mccollision.getGeneratorId()); + // LOGF(info, "mccollision.getSubGeneratorId() = %d", mccollision.getSubGeneratorId()); + // LOGF(info, "mccollision.getSourceId() = %d", mccollision.getSourceId()); + if (cfgEventGeneratorType >= 0 && mccollision.getSubGeneratorId() != cfgEventGeneratorType) { + continue; + } + + auto muonsMC_per_coll = muonsMC->sliceByCachedUnsorted(o2::aod::emmcparticle::emmceventId, mccollision.globalIndex(), cache); + + for (auto& muon : muonsMC_per_coll) { + if (!(muon.isPhysicalPrimary() || muon.producedByGenerator())) { + continue; + } + if (!isInAcceptance(muon)) { + continue; + } + if (!muon.has_mothers()) { + continue; + } + auto mcmother = mcparticles.iteratorAt(muon.mothersIds()[0]); + int pdg_mother = abs(mcmother.pdgCode()); + if (pdg_mother == 111 || pdg_mother == 221 || pdg_mother == 331 || pdg_mother == 113 || pdg_mother == 223 || pdg_mother == 333) { + fRegistry.fill(HIST("Generated/lf/hs"), muon.pt(), muon.eta(), muon.phi(), -muon.pdgCode() / 13); + } else if (pdg_mother == 443) { + if (IsFromBeauty(mcmother, mcparticles) > 0) { // b is found in full decay chain. + fRegistry.fill(HIST("Generated/NonPromptJPsi/hs"), muon.pt(), muon.eta(), muon.phi(), -muon.pdgCode() / 13); + } else { + fRegistry.fill(HIST("Generated/PromptJPsi/hs"), muon.pt(), muon.eta(), muon.phi(), -muon.pdgCode() / 13); + } + } else if (pdg_mother == 100443) { + if (IsFromBeauty(mcmother, mcparticles) > 0) { // b is found in full decay chain. + fRegistry.fill(HIST("Generated/NonPromptPsi2S/hs"), muon.pt(), muon.eta(), muon.phi(), -muon.pdgCode() / 13); + } else { + fRegistry.fill(HIST("Generated/PromptPsi2S/hs"), muon.pt(), muon.eta(), muon.phi(), -muon.pdgCode() / 13); + } + } else if (IsFromBeauty(muon, mcparticles) > 0) { // b is found in full decay chain. + if (IsFromCharm(muon, mcparticles) > 0) { // c is found in full decay chain. + fRegistry.fill(HIST("Generated/b2c2mu/hs"), muon.pt(), muon.eta(), muon.phi(), -muon.pdgCode() / 13); + } else { + fRegistry.fill(HIST("Generated/b2mu/hs"), muon.pt(), muon.eta(), muon.phi(), -muon.pdgCode() / 13); + } + } else if (IsFromCharm(muon, mcparticles) > 0) { // c is found in full decay chain. Not from b. + fRegistry.fill(HIST("Generated/c2mu/hs"), muon.pt(), muon.eta(), muon.phi(), -muon.pdgCode() / 13); + } + } + + } // end of collision loop + } + PROCESS_SWITCH(singleMuonQCMC, processGen, "run genrated info", true); + + void processDummy(MyCollisions const&) {} + PROCESS_SWITCH(singleMuonQCMC, processDummy, "Dummy function", false); +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{ + adaptAnalysisTask(cfgc, TaskName{"single-muon-qc-mc"})}; +}