diff --git a/PWGEM/Dilepton/Core/DielectronCut.cxx b/PWGEM/Dilepton/Core/DielectronCut.cxx index db2514f9dbe..8516f174ca9 100644 --- a/PWGEM/Dilepton/Core/DielectronCut.cxx +++ b/PWGEM/Dilepton/Core/DielectronCut.cxx @@ -103,10 +103,11 @@ void DielectronCut::SetChi2PerClusterITS(float min, float max) mMaxChi2PerClusterITS = max; LOG(info) << "Dielectron Cut, set chi2 per cluster ITS range: " << mMinChi2PerClusterITS << " - " << mMaxChi2PerClusterITS; } -void DielectronCut::SetMeanClusterSizeITSob(float min, float max) +void DielectronCut::SetMeanClusterSizeITS(float min, float max, float maxP) { mMinMeanClusterSizeITS = min; mMaxMeanClusterSizeITS = max; + mMaxP_ITSClusterSize = maxP; LOG(info) << "Dielectron Cut, set mean cluster size ITS range: " << mMinMeanClusterSizeITS << " - " << mMaxMeanClusterSizeITS; } void DielectronCut::SetDca3DRange(float min, float max) diff --git a/PWGEM/Dilepton/Core/DielectronCut.h b/PWGEM/Dilepton/Core/DielectronCut.h index ab4b443b68f..ff7155b778a 100644 --- a/PWGEM/Dilepton/Core/DielectronCut.h +++ b/PWGEM/Dilepton/Core/DielectronCut.h @@ -225,7 +225,7 @@ class DielectronCut : public TNamed std::vector inputFeatures{static_cast(collision.numContrib()), track.p(), track.tgl(), track.tpcNSigmaEl(), /*track.tpcNSigmaMu(),*/ track.tpcNSigmaPi(), track.tpcNSigmaKa(), track.tpcNSigmaPr(), track.tofNSigmaEl(), /*track.tofNSigmaMu(),*/ track.tofNSigmaPi(), track.tofNSigmaKa(), track.tofNSigmaPr(), - track.meanClusterSizeITSob() * std::cos(std::atan(track.tgl()))}; + track.meanClusterSizeITS() * std::cos(std::atan(track.tgl()))}; // calculate classifier float prob_ele = mPIDModel->evalModel(inputFeatures)[0]; @@ -341,7 +341,7 @@ class DielectronCut : public TNamed return mMinChi2PerClusterITS < track.itsChi2NCl() && track.itsChi2NCl() < mMaxChi2PerClusterITS; case DielectronCuts::kITSCluserSize: - return mMinMeanClusterSizeITS < track.meanClusterSizeITSob() * std::cos(std::atan(track.tgl())) && track.meanClusterSizeITSob() * std::cos(std::atan(track.tgl())) < mMaxMeanClusterSizeITS; + return track.p() < mMaxP_ITSClusterSize ? mMinMeanClusterSizeITS < track.meanClusterSizeITS() * std::cos(std::atan(track.tgl())) && track.meanClusterSizeITS() * std::cos(std::atan(track.tgl())) < mMaxMeanClusterSizeITS : true; case DielectronCuts::kPrefilter: return track.pfb() <= 0; @@ -367,7 +367,7 @@ class DielectronCut : public TNamed void SetChi2PerClusterTPC(float min, float max); void SetNClustersITS(int min, int max); void SetChi2PerClusterITS(float min, float max); - void SetMeanClusterSizeITSob(float min, float max); + void SetMeanClusterSizeITS(float min, float max, float maxP = 0.f); void SetPIDScheme(int scheme); void SetMinPinTOF(float min); @@ -440,6 +440,7 @@ class DielectronCut : public TNamed bool mApplyPhiV{true}; bool mApplyPF{false}; float mMinMeanClusterSizeITS{-1e10f}, mMaxMeanClusterSizeITS{1e10f}; // max x cos(Lmabda) + float mMaxP_ITSClusterSize{0.0}; // pid cuts int mPIDScheme{-1}; diff --git a/PWGEM/Dilepton/Core/Dilepton.h b/PWGEM/Dilepton/Core/Dilepton.h index 55926c19882..0dfae159071 100644 --- a/PWGEM/Dilepton/Core/Dilepton.h +++ b/PWGEM/Dilepton/Core/Dilepton.h @@ -530,7 +530,7 @@ struct Dilepton { 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.SetMeanClusterSizeITS(0, 16); fDielectronCut.SetMaxDcaXY(dielectroncuts.cfg_max_dcaxy); fDielectronCut.SetMaxDcaZ(dielectroncuts.cfg_max_dcaz); diff --git a/PWGEM/Dilepton/Core/DileptonMC.h b/PWGEM/Dilepton/Core/DileptonMC.h index d9e975d9b78..0f836a8bf71 100644 --- a/PWGEM/Dilepton/Core/DileptonMC.h +++ b/PWGEM/Dilepton/Core/DileptonMC.h @@ -490,7 +490,7 @@ struct DileptonMC { 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.SetMeanClusterSizeITS(0, 16); fDielectronCut.SetMaxDcaXY(dielectroncuts.cfg_max_dcaxy); fDielectronCut.SetMaxDcaZ(dielectroncuts.cfg_max_dcaz); diff --git a/PWGEM/Dilepton/Core/PhotonHBT.h b/PWGEM/Dilepton/Core/PhotonHBT.h index 6facd2ab1c4..2ca433b8bd5 100644 --- a/PWGEM/Dilepton/Core/PhotonHBT.h +++ b/PWGEM/Dilepton/Core/PhotonHBT.h @@ -188,6 +188,9 @@ struct PhotonHBT { Configurable cfg_max_dcaz{"cfg_max_dcaz", 1.0, "max dca Z for single track in cm"}; Configurable cfg_min_TOFbeta{"cfg_min_TOFbeta", 0.985, "min TOF beta for single track"}; //|beta - 1| < 0.015 corresponds to 3 sigma in pp Configurable cfg_max_TOFbeta{"cfg_max_TOFbeta", 1.015, "max TOF beta for single track"}; //|beta - 1| < 0.015 corresponds to 3 sigma in pp + Configurable cfg_min_its_cluster_size{"cfg_min_its_cluster_size", 0.f, "min ITS cluster size"}; + Configurable cfg_max_its_cluster_size{"cfg_max_its_cluster_size", 16.f, "max ITS cluster size"}; + Configurable cfg_max_p_its_cluster_size{"cfg_max_p_its_cluster_size", 0.2, "max p to apply ITS cluster size cut"}; 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"}; @@ -463,7 +466,7 @@ struct PhotonHBT { 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.SetMeanClusterSizeITS(dielectroncuts.cfg_min_its_cluster_size, dielectroncuts.cfg_max_its_cluster_size, dielectroncuts.cfg_max_p_its_cluster_size); fDielectronCut.SetMaxDcaXY(dielectroncuts.cfg_max_dcaxy); fDielectronCut.SetMaxDcaZ(dielectroncuts.cfg_max_dcaz); diff --git a/PWGEM/Dilepton/Core/SingleTrackQC.h b/PWGEM/Dilepton/Core/SingleTrackQC.h index a7754c3d50a..7f45154dc64 100644 --- a/PWGEM/Dilepton/Core/SingleTrackQC.h +++ b/PWGEM/Dilepton/Core/SingleTrackQC.h @@ -307,7 +307,7 @@ struct SingleTrackQC { 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.SetMeanClusterSizeITS(0, 16); fDielectronCut.SetMaxDcaXY(dielectroncuts.cfg_max_dcaxy); fDielectronCut.SetMaxDcaZ(dielectroncuts.cfg_max_dcaz); fDielectronCut.RequireITSibAny(dielectroncuts.cfg_require_itsib_any); diff --git a/PWGEM/Dilepton/Core/SingleTrackQCMC.h b/PWGEM/Dilepton/Core/SingleTrackQCMC.h index 997a9e3c4c4..e87949f7dd4 100644 --- a/PWGEM/Dilepton/Core/SingleTrackQCMC.h +++ b/PWGEM/Dilepton/Core/SingleTrackQCMC.h @@ -351,7 +351,7 @@ struct SingleTrackQCMC { 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.SetMeanClusterSizeITS(0, 16); fDielectronCut.SetMaxDcaXY(dielectroncuts.cfg_max_dcaxy); fDielectronCut.SetMaxDcaZ(dielectroncuts.cfg_max_dcaz); fDielectronCut.RequireITSibAny(dielectroncuts.cfg_require_itsib_any); diff --git a/PWGEM/Dilepton/Tasks/CMakeLists.txt b/PWGEM/Dilepton/Tasks/CMakeLists.txt index 644490bbddf..3cae03b233b 100644 --- a/PWGEM/Dilepton/Tasks/CMakeLists.txt +++ b/PWGEM/Dilepton/Tasks/CMakeLists.txt @@ -55,6 +55,16 @@ o2physics_add_dpl_workflow(single-electron-qc-mc PUBLIC_LINK_LIBRARIES O2::Framework O2::DetectorsBase O2Physics::AnalysisCore O2Physics::MLCore O2Physics::PWGEMDileptonCore COMPONENT_NAME Analysis) +o2physics_add_dpl_workflow(vp-pair-qc + SOURCES vpPairQC.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::MLCore O2Physics::PWGEMDileptonCore + COMPONENT_NAME Analysis) + +o2physics_add_dpl_workflow(vp-pair-qc-mc + SOURCES vpPairQCMC.cxx + PUBLIC_LINK_LIBRARIES O2::Framework 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 diff --git a/PWGEM/Dilepton/Tasks/vpPairQC.cxx b/PWGEM/Dilepton/Tasks/vpPairQC.cxx new file mode 100644 index 00000000000..75563e9d737 --- /dev/null +++ b/PWGEM/Dilepton/Tasks/vpPairQC.cxx @@ -0,0 +1,478 @@ +// 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 dalitz ee table for dalitz 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 "Common/Core/RecoDecay.h" + +#include "DetectorsBase/GeometryManager.h" +#include "DataFormatsParameters/GRPObject.h" +#include "DataFormatsParameters/GRPMagField.h" +#include "CCDB/BasicCCDBManager.h" + +#include "PWGEM/Dilepton/DataModel/dileptonTables.h" +#include "PWGEM/Dilepton/Core/DielectronCut.h" +#include "PWGEM/Dilepton/Core/EMEventCut.h" +#include "PWGEM/Dilepton/Utils/EMTrack.h" +#include "PWGEM/Dilepton/Utils/EventHistograms.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 vpPairQC { + // 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 cfgCentEstimator{"cfgCentEstimator", 2, "FT0M:0, FT0A:1, FT0C:2"}; + Configurable cfgCentMin{"cfgCentMin", 0, "min. centrality"}; + Configurable cfgCentMax{"cfgCentMax", 999.f, "max. centrality"}; + Configurable maxY{"maxY", 0.9, "maximum rapidity for reconstructed particles"}; + ConfigurableAxis ConfPtlBins{"ConfPtlBins", {VARIABLE_WIDTH, 0.00, 0.10, 0.11, 0.12, 0.13, 0.14, 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.50, 3.00, 3.50, 4.00, 4.50, 5.00, 6.00, 7.00, 8.00, 9.00, 10.00}, "pTl 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"}; + Configurable cfgRequireNoCollInTimeRangeStandard{"cfgRequireNoCollInTimeRangeStandard", false, "require no collision in time range standard"}; + } eventcuts; + + DielectronCut fDielectronCut; + struct : ConfigurableGroup { + std::string prefix = "dielectroncut_group"; + Configurable cfg_min_mass{"cfg_min_mass", 0.0, "min mass"}; + Configurable cfg_max_mass{"cfg_max_mass", 0.01, "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_y{"cfg_min_pair_y", -0.9, "min pair rapidity"}; + Configurable cfg_max_pair_y{"cfg_max_pair_y", +0.9, "max pair rapidity"}; + Configurable cfg_min_pair_dca3d{"cfg_min_pair_dca3d", 0.0, "min pair dca3d in sigma"}; + Configurable cfg_max_pair_dca3d{"cfg_max_pair_dca3d", 1e+10, "max pair dca3d in sigma"}; + Configurable cfg_apply_phiv{"cfg_apply_phiv", true, "flag to apply phiv cut"}; + Configurable cfg_apply_pf{"cfg_apply_pf", false, "flag to apply phiv prefilter"}; + Configurable cfg_require_itsib_any{"cfg_require_itsib_any", false, "flag to require ITS ib any hits"}; + Configurable cfg_require_itsib_1st{"cfg_require_itsib_1st", true, "flag to require ITS ib 1st hit"}; + Configurable cfg_phiv_slope{"cfg_phiv_slope", 0.0185, "slope for m vs. phiv"}; + Configurable cfg_phiv_intercept{"cfg_phiv_intercept", -0.0280, "intercept for m vs. phiv"}; + + 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", -0.9, "min eta for single track"}; + Configurable cfg_max_eta_track{"cfg_max_eta_track", +0.9, "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_min_its_cluster_size{"cfg_min_its_cluster_size", 0.f, "min ITS cluster size"}; + Configurable cfg_max_its_cluster_size{"cfg_max_its_cluster_size", 16.f, "max ITS cluster size"}; + Configurable cfg_max_p_its_cluster_size{"cfg_max_p_its_cluster_size", 0.2, "max p to apply ITS cluster size cut"}; + + Configurable cfg_pid_scheme{"cfg_pid_scheme", static_cast(DielectronCut::PIDSchemes::kTOFif), "pid scheme [kTOFreq : 0, kTPChadrej : 1, kTPChadrejORTOFreq : 2, kTPConly : 3, kTOFif : 4, kPIDML : 5]"}; + 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", -0.0, "min. TPC n sigma for pion exclusion"}; + Configurable cfg_max_TPCNsigmaPi{"cfg_max_TPCNsigmaPi", +0.0, "max. TPC n sigma for pion exclusion"}; + Configurable cfg_min_TPCNsigmaKa{"cfg_min_TPCNsigmaKa", -0.0, "min. TPC n sigma for kaon exclusion"}; + Configurable cfg_max_TPCNsigmaKa{"cfg_max_TPCNsigmaKa", +0.0, "max. TPC n sigma for kaon exclusion"}; + Configurable cfg_min_TPCNsigmaPr{"cfg_min_TPCNsigmaPr", -0.0, "min. TPC n sigma for proton exclusion"}; + Configurable cfg_max_TPCNsigmaPr{"cfg_max_TPCNsigmaPr", +0.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", true, "Flag to enable or disable TTCA"}; + + // 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; + int mRunNumber; + float d_bz; + + HistogramRegistry fRegistry{"output", {}, OutputObjHandlingPolicy::AnalysisObject, false, false}; + static constexpr std::string_view event_cut_types[2] = {"before/", "after/"}; + + void init(InitContext& /*context*/) + { + DefineEMEventCut(); + DefineDielectronCut(); + 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(); + } + + ~vpPairQC() + { + if (eid_bdt) { + delete eid_bdt; + } + } + + void addhistograms() + { + // event info + o2::aod::pwgem::dilepton::utils::eventhistogram::addEventHistograms(&fRegistry); + + const AxisSpec axis_pt{ConfPtlBins, "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_dca{{0.0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9, 1.0, 1.1, 1.2, 1.3, 1.4, 1.5, 1.6, 1.7, 1.8, 1.9, 2.0, 2.5, 3.0, 3.5, 4.0, 4.5, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0}, "DCA_{e}^{3D} (#sigma)"}; + fRegistry.add("Track/positive/hs", "rec. single electron", kTHnSparseD, {axis_pt, axis_eta, axis_phi, axis_dca}, true); + fRegistry.add("Track/positive/hQoverPt", "q/pT;q/p_{T} (GeV/c)^{-1}", kTH1F, {{400, -20, 20}}, false); + fRegistry.add("Track/positive/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/positive/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/positive/hDCAxyRes_Pt", "DCA_{xy} resolution vs. pT;p_{T} (GeV/c);DCA_{xy} resolution (#mum)", kTH2F, {{200, 0, 10}, {200, 0., 400}}, false); + fRegistry.add("Track/positive/hDCAzRes_Pt", "DCA_{z} resolution vs. pT;p_{T} (GeV/c);DCA_{z} resolution (#mum)", kTH2F, {{200, 0, 10}, {200, 0., 400}}, false); + fRegistry.add("Track/positive/hNclsTPC", "number of TPC clusters", kTH1F, {{161, -0.5, 160.5}}, false); + fRegistry.add("Track/positive/hNcrTPC", "number of TPC crossed rows", kTH1F, {{161, -0.5, 160.5}}, false); + fRegistry.add("Track/positive/hChi2TPC", "chi2/number of TPC clusters", kTH1F, {{100, 0, 10}}, false); + fRegistry.add("Track/positive/hTPCNcr2Nf", "TPC Ncr/Nfindable", kTH1F, {{200, 0, 2}}, false); + fRegistry.add("Track/positive/hTPCNcls2Nf", "TPC Ncls/Nfindable", kTH1F, {{200, 0, 2}}, false); + fRegistry.add("Track/positive/hNclsITS", "number of ITS clusters", kTH1F, {{8, -0.5, 7.5}}, false); + fRegistry.add("Track/positive/hChi2ITS", "chi2/number of ITS clusters", kTH1F, {{100, 0, 10}}, false); + fRegistry.add("Track/positive/hITSClusterMap", "ITS cluster map", kTH1F, {{128, -0.5, 127.5}}, false); + fRegistry.add("Track/positive/hTPCdEdx", "TPC dE/dx;p_{in} (GeV/c);TPC dE/dx (a.u.)", kTH2F, {{1000, 0, 10}, {200, 0, 200}}, false); + fRegistry.add("Track/positive/hTOFbeta", "TOF #beta;p_{pv} (GeV/c);#beta", kTH2F, {{1000, 0, 10}, {240, 0, 1.2}}, false); + fRegistry.add("Track/positive/hMeanClusterSizeITS", "mean cluster size ITS;p_{pv} (GeV/c); on ITS #times cos(#lambda);", kTH2F, {{1000, 0.f, 10.f}, {160, 0, 16}}, false); + fRegistry.add("Track/positive/hTPCNsigmaEl", "TPC n sigma el;p_{in} (GeV/c);n #sigma_{e}^{TPC}", kTH2F, {{1000, 0, 10}, {100, -5, +5}}, false); + fRegistry.add("Track/positive/hTPCNsigmaMu", "TPC n sigma mu;p_{in} (GeV/c);n #sigma_{#mu}^{TPC}", kTH2F, {{1000, 0, 10}, {100, -5, +5}}, false); + fRegistry.add("Track/positive/hTPCNsigmaPi", "TPC n sigma pi;p_{in} (GeV/c);n #sigma_{#pi}^{TPC}", kTH2F, {{1000, 0, 10}, {100, -5, +5}}, false); + fRegistry.add("Track/positive/hTPCNsigmaKa", "TPC n sigma ka;p_{in} (GeV/c);n #sigma_{K}^{TPC}", kTH2F, {{1000, 0, 10}, {100, -5, +5}}, false); + fRegistry.add("Track/positive/hTPCNsigmaPr", "TPC n sigma pr;p_{in} (GeV/c);n #sigma_{p}^{TPC}", kTH2F, {{1000, 0, 10}, {100, -5, +5}}, false); + fRegistry.add("Track/positive/hTOFNsigmaEl", "TOF n sigma el;p_{pv} (GeV/c);n #sigma_{e}^{TOF}", kTH2F, {{1000, 0, 10}, {100, -5, +5}}, false); + fRegistry.add("Track/positive/hTOFNsigmaMu", "TOF n sigma mu;p_{pv} (GeV/c);n #sigma_{#mu}^{TOF}", kTH2F, {{1000, 0, 10}, {100, -5, +5}}, false); + fRegistry.add("Track/positive/hTOFNsigmaPi", "TOF n sigma pi;p_{pv} (GeV/c);n #sigma_{#pi}^{TOF}", kTH2F, {{1000, 0, 10}, {100, -5, +5}}, false); + fRegistry.add("Track/positive/hTOFNsigmaKa", "TOF n sigma ka;p_{pv} (GeV/c);n #sigma_{K}^{TOF}", kTH2F, {{1000, 0, 10}, {100, -5, +5}}, false); + fRegistry.add("Track/positive/hTOFNsigmaPr", "TOF n sigma pr;p_{pv} (GeV/c);n #sigma_{p}^{TOF}", kTH2F, {{1000, 0, 10}, {100, -5, +5}}, false); + fRegistry.addClone("Track/positive/", "Track/negative/"); + + // for pair + fRegistry.add("Pair/hMvsPt", "m_{ee} vs. p_{T,ee};m_{ee} (GeV/c^{2});p_{T,ee} (GeV/c)", kTH2F, {{100, 0, 0.1}, {100, 0, 1}}, true); + fRegistry.add("Pair/hMvsPhiV", "m_{ee} vs. #varphi_{V};#varphi (rad.);m_{ee} (GeV/c^{2})", kTH2F, {{90, 0, M_PI}, {100, 0.0f, 0.1f}}, true); + } + + 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); + fEMEventCut.SetRequireNoCollInTimeRangeStandard(eventcuts.cfgRequireNoCollInTimeRangeStandard); + } + + o2::ml::OnnxModel* eid_bdt = nullptr; + void DefineDielectronCut() + { + fDielectronCut = DielectronCut("fDielectronCut", "fDielectronCut"); + + // for pair + fDielectronCut.SetMeeRange(dielectroncuts.cfg_min_mass, dielectroncuts.cfg_max_mass); + fDielectronCut.SetPairPtRange(dielectroncuts.cfg_min_pair_pt, dielectroncuts.cfg_max_pair_pt); + fDielectronCut.SetPairYRange(dielectroncuts.cfg_min_pair_y, dielectroncuts.cfg_max_pair_y); + fDielectronCut.SetPairDCARange(dielectroncuts.cfg_min_pair_dca3d, dielectroncuts.cfg_max_pair_dca3d); // in sigma + fDielectronCut.SetMaxPhivPairMeeDep([&](float mll) { return (mll - dielectroncuts.cfg_phiv_intercept) / dielectroncuts.cfg_phiv_slope; }); + 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 + fDielectronCut.SetTrackPtRange(dielectroncuts.cfg_min_pt_track, 1e+10f); + fDielectronCut.SetTrackEtaRange(dielectroncuts.cfg_min_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.SetMeanClusterSizeITS(dielectroncuts.cfg_min_its_cluster_size, dielectroncuts.cfg_max_its_cluster_size, dielectroncuts.cfg_max_p_its_cluster_size); + fDielectronCut.SetMaxDcaXY(dielectroncuts.cfg_max_dcaxy); + fDielectronCut.SetMaxDcaZ(dielectroncuts.cfg_max_dcaz); + + // for eID + 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(); + 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); + } + + fDielectronCut.SetPIDModel(eid_bdt); + } // end of PID ML + } + + template + bool fillPairInfo(TCollision const&, TTrack1 const& t1, TTrack2 const& t2) + { + // if (t1.trackId() == t2.trackId()) { // this is protection against pairing identical 2 tracks. This happens, when TTCA is used. TTCA can assign a track to several possible collisions. + // return false; + // } + + if (!fDielectronCut.IsSelectedTrack(t1) || !fDielectronCut.IsSelectedTrack(t2)) { + return false; + } + + if (!fDielectronCut.IsSelectedPair(t1, t2, d_bz)) { + return false; + } + + ROOT::Math::PtEtaPhiMVector v1(t1.pt(), t1.eta(), t1.phi(), o2::constants::physics::MassElectron); + ROOT::Math::PtEtaPhiMVector v2(t2.pt(), t2.eta(), t2.phi(), o2::constants::physics::MassElectron); + ROOT::Math::PtEtaPhiMVector v12 = v1 + v2; + if (abs(v12.Rapidity()) > maxY) { + return false; + } + float phiv = o2::aod::pwgem::dilepton::utils::pairutil::getPhivPair(t1.px(), t1.py(), t1.pz(), t2.px(), t2.py(), t2.pz(), t1.sign(), t2.sign(), d_bz); + + if (t1.sign() * t2.sign() < 0) { // ULS + fRegistry.fill(HIST("Pair/hMvsPt"), v12.M(), v12.Pt()); + fRegistry.fill(HIST("Pair/hMvsPhiV"), phiv, v12.M()); + } + + if (std::find(used_trackIds.begin(), used_trackIds.end(), t1.globalIndex()) == used_trackIds.end()) { + used_trackIds.emplace_back(t1.globalIndex()); + fillTrackInfo(t1); + } + if (std::find(used_trackIds.begin(), used_trackIds.end(), t2.globalIndex()) == used_trackIds.end()) { + used_trackIds.emplace_back(t2.globalIndex()); + fillTrackInfo(t2); + } + + return true; + } + + template + void fillTrackInfo(TTrack const& track) + { + float weight = 1.f; + float dca_3d = dca3DinSigma(track); + if (track.sign() > 0) { + fRegistry.fill(HIST("Track/positive/hs"), track.pt(), track.eta(), track.phi(), dca_3d, weight); + fRegistry.fill(HIST("Track/positive/hQoverPt"), track.sign() / track.pt()); + fRegistry.fill(HIST("Track/positive/hDCAxyz"), track.dcaXY(), track.dcaZ()); + fRegistry.fill(HIST("Track/positive/hDCAxyzSigma"), track.dcaXY() / sqrt(track.cYY()), track.dcaZ() / sqrt(track.cZZ())); + fRegistry.fill(HIST("Track/positive/hDCAxyRes_Pt"), track.pt(), sqrt(track.cYY()) * 1e+4); // convert cm to um + fRegistry.fill(HIST("Track/positive/hDCAzRes_Pt"), track.pt(), sqrt(track.cZZ()) * 1e+4); // convert cm to um + fRegistry.fill(HIST("Track/positive/hNclsITS"), track.itsNCls()); + fRegistry.fill(HIST("Track/positive/hNclsTPC"), track.tpcNClsFound()); + fRegistry.fill(HIST("Track/positive/hNcrTPC"), track.tpcNClsCrossedRows()); + fRegistry.fill(HIST("Track/positive/hTPCNcr2Nf"), track.tpcCrossedRowsOverFindableCls()); + fRegistry.fill(HIST("Track/positive/hTPCNcls2Nf"), track.tpcFoundOverFindableCls()); + fRegistry.fill(HIST("Track/positive/hChi2TPC"), track.tpcChi2NCl()); + fRegistry.fill(HIST("Track/positive/hChi2ITS"), track.itsChi2NCl()); + fRegistry.fill(HIST("Track/positive/hITSClusterMap"), track.itsClusterMap()); + + fRegistry.fill(HIST("Track/positive/hTPCdEdx"), track.tpcInnerParam(), track.tpcSignal()); + fRegistry.fill(HIST("Track/positive/hTOFbeta"), track.p(), track.beta()); + fRegistry.fill(HIST("Track/positive/hMeanClusterSizeITS"), track.p(), track.meanClusterSizeITS() * std::cos(std::atan(track.tgl()))); + fRegistry.fill(HIST("Track/positive/hTPCNsigmaEl"), track.tpcInnerParam(), track.tpcNSigmaEl()); + fRegistry.fill(HIST("Track/positive/hTPCNsigmaMu"), track.tpcInnerParam(), track.tpcNSigmaMu()); + fRegistry.fill(HIST("Track/positive/hTPCNsigmaPi"), track.tpcInnerParam(), track.tpcNSigmaPi()); + fRegistry.fill(HIST("Track/positive/hTPCNsigmaKa"), track.tpcInnerParam(), track.tpcNSigmaKa()); + fRegistry.fill(HIST("Track/positive/hTPCNsigmaPr"), track.tpcInnerParam(), track.tpcNSigmaPr()); + fRegistry.fill(HIST("Track/positive/hTOFNsigmaEl"), track.p(), track.tofNSigmaEl()); + fRegistry.fill(HIST("Track/positive/hTOFNsigmaMu"), track.p(), track.tofNSigmaMu()); + fRegistry.fill(HIST("Track/positive/hTOFNsigmaPi"), track.p(), track.tofNSigmaPi()); + fRegistry.fill(HIST("Track/positive/hTOFNsigmaKa"), track.p(), track.tofNSigmaKa()); + fRegistry.fill(HIST("Track/positive/hTOFNsigmaPr"), track.p(), track.tofNSigmaPr()); + } else { + fRegistry.fill(HIST("Track/negative/hs"), track.pt(), track.eta(), track.phi(), dca_3d, weight); + fRegistry.fill(HIST("Track/negative/hQoverPt"), track.sign() / track.pt()); + fRegistry.fill(HIST("Track/negative/hDCAxyz"), track.dcaXY(), track.dcaZ()); + fRegistry.fill(HIST("Track/negative/hDCAxyzSigma"), track.dcaXY() / sqrt(track.cYY()), track.dcaZ() / sqrt(track.cZZ())); + fRegistry.fill(HIST("Track/negative/hDCAxyRes_Pt"), track.pt(), sqrt(track.cYY()) * 1e+4); // convert cm to um + fRegistry.fill(HIST("Track/negative/hDCAzRes_Pt"), track.pt(), sqrt(track.cZZ()) * 1e+4); // convert cm to um + fRegistry.fill(HIST("Track/negative/hNclsITS"), track.itsNCls()); + fRegistry.fill(HIST("Track/negative/hNclsTPC"), track.tpcNClsFound()); + fRegistry.fill(HIST("Track/negative/hNcrTPC"), track.tpcNClsCrossedRows()); + fRegistry.fill(HIST("Track/negative/hTPCNcr2Nf"), track.tpcCrossedRowsOverFindableCls()); + fRegistry.fill(HIST("Track/negative/hTPCNcls2Nf"), track.tpcFoundOverFindableCls()); + fRegistry.fill(HIST("Track/negative/hChi2TPC"), track.tpcChi2NCl()); + fRegistry.fill(HIST("Track/negative/hChi2ITS"), track.itsChi2NCl()); + fRegistry.fill(HIST("Track/negative/hITSClusterMap"), track.itsClusterMap()); + + fRegistry.fill(HIST("Track/negative/hTPCdEdx"), track.tpcInnerParam(), track.tpcSignal()); + fRegistry.fill(HIST("Track/negative/hTOFbeta"), track.p(), track.beta()); + fRegistry.fill(HIST("Track/negative/hMeanClusterSizeITS"), track.p(), track.meanClusterSizeITS() * std::cos(std::atan(track.tgl()))); + fRegistry.fill(HIST("Track/negative/hTPCNsigmaEl"), track.tpcInnerParam(), track.tpcNSigmaEl()); + fRegistry.fill(HIST("Track/negative/hTPCNsigmaMu"), track.tpcInnerParam(), track.tpcNSigmaMu()); + fRegistry.fill(HIST("Track/negative/hTPCNsigmaPi"), track.tpcInnerParam(), track.tpcNSigmaPi()); + fRegistry.fill(HIST("Track/negative/hTPCNsigmaKa"), track.tpcInnerParam(), track.tpcNSigmaKa()); + fRegistry.fill(HIST("Track/negative/hTPCNsigmaPr"), track.tpcInnerParam(), track.tpcNSigmaPr()); + fRegistry.fill(HIST("Track/negative/hTOFNsigmaEl"), track.p(), track.tofNSigmaEl()); + fRegistry.fill(HIST("Track/negative/hTOFNsigmaMu"), track.p(), track.tofNSigmaMu()); + fRegistry.fill(HIST("Track/negative/hTOFNsigmaPi"), track.p(), track.tofNSigmaPi()); + fRegistry.fill(HIST("Track/negative/hTOFNsigmaKa"), track.p(), track.tofNSigmaKa()); + fRegistry.fill(HIST("Track/negative/hTOFNsigmaPr"), track.p(), 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; + + SliceCache cache; + Preslice perCollision_track = aod::emprimaryelectron::emeventId; + Filter trackFilter_electron = dielectroncuts.cfg_min_pt_track < o2::aod::track::pt && dielectroncuts.cfg_min_eta_track < o2::aod::track::eta && 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_electron = dielectroncuts.cfg_min_TPCNsigmaEl < o2::aod::pidtpc::tpcNSigmaEl && o2::aod::pidtpc::tpcNSigmaEl < dielectroncuts.cfg_max_TPCNsigmaEl; + Filter ttcaFilter_electron = ifnode(dielectroncuts.enableTTCA.node(), o2::aod::emprimaryelectron::isAssociatedToMPC == true || o2::aod::emprimaryelectron::isAssociatedToMPC == false, o2::aod::emprimaryelectron::isAssociatedToMPC == true); + using FilteredMyTracks = soa::Filtered; + + Partition posTracks = o2::aod::emprimaryelectron::sign > int8_t(0); + Partition negTracks = o2::aod::emprimaryelectron::sign < int8_t(0); + + std::vector used_trackIds; + void processQC(FilteredMyCollisions const& collisions, FilteredMyTracks const& /*tracks*/) + { + for (auto& collision : collisions) { + initCCDB(collision); + const 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>(&fRegistry, collision); + if (!fEMEventCut.IsSelected(collision)) { + continue; + } + o2::aod::pwgem::dilepton::utils::eventhistogram::fillEventInfo<1>(&fRegistry, collision); + fRegistry.fill(HIST("Event/before/hCollisionCounter"), 12.0); // accepted + fRegistry.fill(HIST("Event/after/hCollisionCounter"), 12.0); // accepted + + auto posTracks_per_coll = posTracks->sliceByCached(o2::aod::emprimaryelectron::emeventId, collision.globalIndex(), cache); + auto negTracks_per_coll = negTracks->sliceByCached(o2::aod::emprimaryelectron::emeventId, collision.globalIndex(), cache); + // LOGF(info, "centrality = %f , posTracks_per_coll.size() = %d, negTracks_per_coll.size() = %d", centralities[cfgCentEstimator], posTracks_per_coll.size(), negTracks_per_coll.size()); + + for (auto& [pos, ele] : combinations(CombinationsFullIndexPolicy(posTracks_per_coll, negTracks_per_coll))) { // ULS + fillPairInfo(collision, pos, ele); + } + } // end of collision loop + + used_trackIds.clear(); + used_trackIds.shrink_to_fit(); + + } // end of process + PROCESS_SWITCH(vpPairQC, processQC, "run vp pair QC", true); + + void processDummy(MyCollisions const&) {} + PROCESS_SWITCH(vpPairQC, processDummy, "Dummy function", false); +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{ + adaptAnalysisTask(cfgc, TaskName{"vp-pair-qc"})}; +} diff --git a/PWGEM/Dilepton/Tasks/vpPairQCMC.cxx b/PWGEM/Dilepton/Tasks/vpPairQCMC.cxx new file mode 100644 index 00000000000..07c0cabd5d4 --- /dev/null +++ b/PWGEM/Dilepton/Tasks/vpPairQCMC.cxx @@ -0,0 +1,621 @@ +// 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 dalitz ee table for dalitz 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 "Common/Core/RecoDecay.h" + +#include "PWGEM/Dilepton/DataModel/dileptonTables.h" +#include "PWGEM/Dilepton/Core/DielectronCut.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 MyCollisions = soa::Join; +using MyCollision = MyCollisions::iterator; + +using MyMCTracks = soa::Join; +using MyMCTrack = MyMCTracks::iterator; + +struct vpPairQCMC { + + // 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 cfgCentEstimator{"cfgCentEstimator", 2, "FT0M:0, FT0A:1, FT0C:2"}; + Configurable cfgCentMin{"cfgCentMin", 0, "min. centrality"}; + Configurable cfgCentMax{"cfgCentMax", 999.f, "max. centrality"}; + Configurable maxY{"maxY", 0.9, "maximum rapidity for reconstructed particles"}; + + 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"}; + Configurable cfgRequireNoCollInTimeRangeStandard{"cfgRequireNoCollInTimeRangeStandard", false, "require no collision in time range standard"}; + } eventcuts; + + DielectronCut fDielectronCut; + struct : ConfigurableGroup { + std::string prefix = "dielectroncut_group"; + Configurable cfg_min_mass{"cfg_min_mass", 0.0, "min mass"}; + Configurable cfg_max_mass{"cfg_max_mass", 0.01, "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_y{"cfg_min_pair_y", -0.9, "min pair rapidity"}; + Configurable cfg_max_pair_y{"cfg_max_pair_y", +0.9, "max pair rapidity"}; + Configurable cfg_min_pair_dca3d{"cfg_min_pair_dca3d", 0.0, "min pair dca3d in sigma"}; + Configurable cfg_max_pair_dca3d{"cfg_max_pair_dca3d", 1e+10, "max pair dca3d in sigma"}; + Configurable cfg_apply_phiv{"cfg_apply_phiv", true, "flag to apply phiv cut"}; + Configurable cfg_apply_pf{"cfg_apply_pf", false, "flag to apply phiv prefilter"}; + Configurable cfg_require_itsib_any{"cfg_require_itsib_any", false, "flag to require ITS ib any hits"}; + Configurable cfg_require_itsib_1st{"cfg_require_itsib_1st", true, "flag to require ITS ib 1st hit"}; + Configurable cfg_phiv_slope{"cfg_phiv_slope", 0.0185, "slope for m vs. phiv"}; + Configurable cfg_phiv_intercept{"cfg_phiv_intercept", -0.0280, "intercept for m vs. phiv"}; + + 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", -0.9, "min eta for single track"}; + Configurable cfg_max_eta_track{"cfg_max_eta_track", +0.9, "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_pid_scheme{"cfg_pid_scheme", static_cast(DielectronCut::PIDSchemes::kTOFif), "pid scheme [kTOFreq : 0, kTPChadrej : 1, kTPChadrejORTOFreq : 2, kTPConly : 3, kTOFif : 4, kPIDML : 5]"}; + 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", -0.0, "min. TPC n sigma for pion exclusion"}; + Configurable cfg_max_TPCNsigmaPi{"cfg_max_TPCNsigmaPi", +0.0, "max. TPC n sigma for pion exclusion"}; + Configurable cfg_min_TPCNsigmaKa{"cfg_min_TPCNsigmaKa", -0.0, "min. TPC n sigma for kaon exclusion"}; + Configurable cfg_max_TPCNsigmaKa{"cfg_max_TPCNsigmaKa", +0.0, "max. TPC n sigma for kaon exclusion"}; + Configurable cfg_min_TPCNsigmaPr{"cfg_min_TPCNsigmaPr", -0.0, "min. TPC n sigma for proton exclusion"}; + Configurable cfg_max_TPCNsigmaPr{"cfg_max_TPCNsigmaPr", +0.0, "max. TPC n sigma for proton exclusion"}; + Configurable cfg_min_TOFNsigmaEl{"cfg_min_TOFNsigmaEl", -0.0, "min. TOF n sigma for electron inclusion"}; + Configurable cfg_max_TOFNsigmaEl{"cfg_max_TOFNsigmaEl", +0.0, "max. TOF n sigma for electron inclusion"}; + Configurable enableTTCA{"enableTTCA", true, "Flag to enable or disable TTCA"}; + + // 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; + int mRunNumber; + float d_bz; + + 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 max_mcEta{"max_mcEta", 0.9, "max. MC eta"}; + } mctrackcuts; + + HistogramRegistry fRegistry{"output", {}, OutputObjHandlingPolicy::AnalysisObject, false, false}; + static constexpr std::string_view event_cut_types[2] = {"before/", "after/"}; + + ~vpPairQCMC() {} + + void addhistograms() + { + // event info + o2::aod::pwgem::dilepton::utils::eventhistogram::addEventHistograms(&fRegistry); + + std::vector ptbins; + std::vector massbins; + + for (int i = 0; i < 51; i++) { + massbins.emplace_back(0.01 * (i - 0) + 0.0); // every 0.01 GeV/c2 from 0.0 to 0.5 GeV/c2 + } + const AxisSpec axis_mass{massbins, "m_{ee} (GeV/c^{2})"}; + + for (int i = 0; i < 50; i++) { + ptbins.emplace_back(0.1 * (i - 0) + 0.0); // every 0.1 GeV/c from 0.0 to 5.0 GeV/c + } + for (int i = 50; i < 61; i++) { + ptbins.emplace_back(0.5 * (i - 50) + 5.0); // every 0.5 GeV/c from 5.0 to 10 GeV/c + } + const AxisSpec axis_pt{ptbins, "p_{T,ee} (GeV/c)"}; + + // generated info + fRegistry.add("Generated/sm/Pi0/hMvsPt", "m_{ee} vs. p_{T,ee} ULS", kTH2F, {axis_mass, axis_pt}, true); + fRegistry.addClone("Generated/sm/Pi0/", "Generated/sm/Eta/"); + fRegistry.addClone("Generated/sm/Pi0/", "Generated/sm/EtaPrime/"); + fRegistry.addClone("Generated/sm/Pi0/", "Generated/sm/Rho/"); + fRegistry.addClone("Generated/sm/Pi0/", "Generated/sm/Omega/"); + fRegistry.addClone("Generated/sm/Pi0/", "Generated/sm/Phi/"); + + // reconstructed pair info + fRegistry.add("Pair/sm/Photon/hMvsPt", "m_{ee} vs. p_{T,ee} ULS", kTH2F, {axis_mass, axis_pt}, true); + fRegistry.add("Pair/sm/Photon/hMvsPhiV", "m_{ee} vs. #varphi_{V};#varphi (rad.);m_{ee} (GeV/c^{2})", kTH2F, {{90, 0, M_PI}, {100, 0.0f, 0.1f}}, false); + fRegistry.addClone("Pair/sm/Photon/", "Pair/sm/Pi0/"); + fRegistry.addClone("Pair/sm/Photon/", "Pair/sm/Eta/"); + fRegistry.addClone("Pair/sm/Photon/", "Pair/sm/EtaPrime/"); + fRegistry.addClone("Pair/sm/Photon/", "Pair/sm/Rho/"); + fRegistry.addClone("Pair/sm/Photon/", "Pair/sm/Omega/"); + fRegistry.addClone("Pair/sm/Photon/", "Pair/sm/Phi/"); + + // 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/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/hTPCNsigmaPi", "TPC n sigma pi;p_{in} (GeV/c);n #sigma_{#pi}^{TPC}", 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&) + { + DefineEMEventCut(); + DefineDielectronCut(); + 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"); + 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); + fEMEventCut.SetRequireNoCollInTimeRangeStandard(eventcuts.cfgRequireNoCollInTimeRangeStandard); + } + + o2::ml::OnnxModel* eid_bdt = nullptr; + void DefineDielectronCut() + { + fDielectronCut = DielectronCut("fDielectronCut", "fDielectronCut"); + + // for pair + fDielectronCut.SetMeeRange(dielectroncuts.cfg_min_mass, dielectroncuts.cfg_max_mass); + fDielectronCut.SetPairPtRange(dielectroncuts.cfg_min_pair_pt, dielectroncuts.cfg_max_pair_pt); + fDielectronCut.SetPairYRange(dielectroncuts.cfg_min_pair_y, dielectroncuts.cfg_max_pair_y); + fDielectronCut.SetPairDCARange(dielectroncuts.cfg_min_pair_dca3d, dielectroncuts.cfg_max_pair_dca3d); // in sigma + fDielectronCut.SetMaxPhivPairMeeDep([&](float mll) { return (mll - dielectroncuts.cfg_phiv_intercept) / dielectroncuts.cfg_phiv_slope; }); + 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 + 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.SetMeanClusterSizeITS(0, 16); + fDielectronCut.SetMaxDcaXY(dielectroncuts.cfg_max_dcaxy); + fDielectronCut.SetMaxDcaZ(dielectroncuts.cfg_max_dcaz); + + // for eID + 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 DefineDielectronCut + // o2::ml::OnnxModel* eid_bdt = new o2::ml::OnnxModel(); + 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); + } + + fDielectronCut.SetPIDModel(eid_bdt); + } // end of PID ML + } + + template + int FindLF(TTrack const& posmc, TTrack const& elemc, TMCParticles const& mcparticles) + { + int arr[] = { + FindCommonMotherFrom2Prongs(posmc, elemc, -11, 11, 22, mcparticles), + FindCommonMotherFrom2Prongs(posmc, elemc, -11, 11, 111, mcparticles), + FindCommonMotherFrom2Prongs(posmc, elemc, -11, 11, 221, mcparticles), + FindCommonMotherFrom2Prongs(posmc, elemc, -11, 11, 331, mcparticles), + FindCommonMotherFrom2Prongs(posmc, elemc, -11, 11, 113, mcparticles), + FindCommonMotherFrom2Prongs(posmc, elemc, -11, 11, 223, mcparticles), + FindCommonMotherFrom2Prongs(posmc, elemc, -11, 11, 333, mcparticles), + FindCommonMotherFrom2Prongs(posmc, elemc, -11, 11, 443, mcparticles), + FindCommonMotherFrom2Prongs(posmc, elemc, -11, 11, 100443, mcparticles)}; + int size = sizeof(arr) / sizeof(*arr); + int max = *std::max_element(arr, arr + size); + return max; + } + + template + bool isInAcceptance(T const& t1) + { + if ((mctrackcuts.min_mcPt < t1.pt() && t1.pt() < mctrackcuts.max_mcPt) && abs(t1.eta()) < mctrackcuts.max_mcEta) { + return true; + } else { + return false; + } + } + + template + bool fillTruePairInfo(TCollision const&, TTrack1 const& t1, TTrack2 const& t2, TMCParticles const& mcparticles) + { + if (!fDielectronCut.IsSelectedTrack(t1) || !fDielectronCut.IsSelectedTrack(t2)) { + return false; + } + + if (!fDielectronCut.IsSelectedPair(t1, t2, d_bz)) { + return false; + } + + auto t1mc = t1.template emmcparticle_as(); + auto t2mc = t2.template emmcparticle_as(); + + int mother_id = FindLF(t1mc, t2mc, mcparticles); + int hfee_type = IsHF(t1mc, t2mc, mcparticles); + if (mother_id < 0 && hfee_type < 0) { + return false; + } + ROOT::Math::PtEtaPhiMVector v1(t1.pt(), t1.eta(), t1.phi(), o2::constants::physics::MassElectron); + ROOT::Math::PtEtaPhiMVector v2(t2.pt(), t2.eta(), t2.phi(), o2::constants::physics::MassElectron); + ROOT::Math::PtEtaPhiMVector v12 = v1 + v2; + if (abs(v12.Rapidity()) > maxY) { + return false; + } + float phiv = o2::aod::pwgem::dilepton::utils::pairutil::getPhivPair(t1.px(), t1.py(), t1.pz(), t2.px(), t2.py(), t2.pz(), t1.sign(), t2.sign(), d_bz); + + if (mother_id > -1 && t1mc.pdgCode() * t2mc.pdgCode() < 0) { + auto mcmother = mcparticles.iteratorAt(mother_id); + if (mcmother.isPhysicalPrimary() || mcmother.producedByGenerator()) { + if ((t1mc.isPhysicalPrimary() || t1mc.producedByGenerator()) && (t2mc.isPhysicalPrimary() || t2mc.producedByGenerator())) { + switch (abs(mcmother.pdgCode())) { + case 111: + fRegistry.fill(HIST("Pair/sm/Pi0/hMvsPt"), v12.M(), v12.Pt()); + fRegistry.fill(HIST("Pair/sm/Pi0/hMvsPhiV"), phiv, v12.M()); + break; + case 221: + fRegistry.fill(HIST("Pair/sm/Eta/hMvsPt"), v12.M(), v12.Pt()); + fRegistry.fill(HIST("Pair/sm/Eta/hMvsPhiV"), phiv, v12.M()); + break; + case 331: + fRegistry.fill(HIST("Pair/sm/EtaPrime/hMvsPt"), v12.M(), v12.Pt()); + fRegistry.fill(HIST("Pair/sm/EtaPrime/hMvsPhiV"), phiv, v12.M()); + break; + case 113: + fRegistry.fill(HIST("Pair/sm/Rho/hMvsPt"), v12.M(), v12.Pt()); + fRegistry.fill(HIST("Pair/sm/Rho/hMvsPhiV"), phiv, v12.M()); + break; + case 223: + fRegistry.fill(HIST("Pair/sm/Omega/hMvsPt"), v12.M(), v12.Pt()); + fRegistry.fill(HIST("Pair/sm/Omega/hMvsPhiV"), phiv, v12.M()); + break; + case 333: + fRegistry.fill(HIST("Pair/sm/Phi/hMvsPt"), v12.M(), v12.Pt()); + fRegistry.fill(HIST("Pair/sm/Phi/hMvsPhiV"), phiv, v12.M()); + break; + default: + break; + } + } else if (!(t1mc.isPhysicalPrimary() || t1mc.producedByGenerator()) && !(t2mc.isPhysicalPrimary() || t2mc.producedByGenerator())) { + switch (abs(mcmother.pdgCode())) { + case 22: + fRegistry.fill(HIST("Pair/sm/Photon/hMvsPt"), v12.M(), v12.Pt()); + fRegistry.fill(HIST("Pair/sm/Photon/hMvsPhiV"), phiv, v12.M()); + break; + default: + break; + } + } // end of primary/secondary selection + } // end of primary selection for same mother + } + + // fill track info that belong to true pairs. + if (t1.sign() > 0) { + if (std::find(used_trackIds.begin(), used_trackIds.end(), t1.globalIndex()) == used_trackIds.end()) { + used_trackIds.emplace_back(t1.globalIndex()); + fillTrackInfo(t1); + } + } else { + if (std::find(used_trackIds.begin(), used_trackIds.end(), t1.globalIndex()) == used_trackIds.end()) { + used_trackIds.emplace_back(t1.globalIndex()); + fillTrackInfo(t1); + } + } + if (t2.sign() > 0) { + if (std::find(used_trackIds.begin(), used_trackIds.end(), t1.globalIndex()) == used_trackIds.end()) { + used_trackIds.emplace_back(t1.globalIndex()); + fillTrackInfo(t2); + } + } else { + if (std::find(used_trackIds.begin(), used_trackIds.end(), t2.globalIndex()) == used_trackIds.end()) { + used_trackIds.emplace_back(t2.globalIndex()); + fillTrackInfo(t2); + } + } + + return true; + } + + template + void fillTrackInfo(TTrack const& 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/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/hTPCNsigmaPi"), track.tpcInnerParam(), track.tpcNSigmaPi()); + } + + std::vector used_trackIds; + SliceCache cache; + Preslice perCollision_track = aod::emprimaryelectron::emeventId; + Filter trackFilter_electron = dielectroncuts.cfg_min_pt_track < o2::aod::track::pt && dielectroncuts.cfg_min_eta_track < o2::aod::track::eta && 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_electron = dielectroncuts.cfg_min_TPCNsigmaEl < o2::aod::pidtpc::tpcNSigmaEl && o2::aod::pidtpc::tpcNSigmaEl < dielectroncuts.cfg_max_TPCNsigmaEl; + Filter ttcaFilter_electron = ifnode(dielectroncuts.enableTTCA.node(), o2::aod::emprimaryelectron::isAssociatedToMPC == true || o2::aod::emprimaryelectron::isAssociatedToMPC == false, o2::aod::emprimaryelectron::isAssociatedToMPC == true); + + using FilteredMyMCTracks = soa::Filtered; + Partition posTracks = o2::aod::emprimaryelectron::sign > int8_t(0); + Partition negTracks = o2::aod::emprimaryelectron::sign < int8_t(0); + + 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&) + { + 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>(&fRegistry, collision); + if (!fEMEventCut.IsSelected(collision)) { + continue; + } + o2::aod::pwgem::dilepton::utils::eventhistogram::fillEventInfo<1>(&fRegistry, collision); + fRegistry.fill(HIST("Event/before/hCollisionCounter"), 12.0); // accepted + fRegistry.fill(HIST("Event/after/hCollisionCounter"), 12.0); // accepted + + auto posTracks_per_coll = posTracks->sliceByCached(o2::aod::emprimaryelectron::emeventId, collision.globalIndex(), cache); + auto negTracks_per_coll = negTracks->sliceByCached(o2::aod::emprimaryelectron::emeventId, collision.globalIndex(), cache); + // LOGF(info, "centrality = %f , posTracks_per_coll.size() = %d, negTracks_per_coll.size() = %d", centralities[cfgCentEstimator], posTracks_per_coll.size(), negTracks_per_coll.size()); + + for (auto& [pos, ele] : combinations(CombinationsFullIndexPolicy(posTracks_per_coll, negTracks_per_coll))) { // ULS + fillTruePairInfo(collision, pos, ele, mcparticles); + } // end of ULS pair loop + + } // end of collision loop + + used_trackIds.clear(); + used_trackIds.shrink_to_fit(); + } // end of process + PROCESS_SWITCH(vpPairQCMC, processQCMC, "run Dalitz QC", true); + + Partition posTracksMC = o2::aod::mcparticle::pdgCode == -11; // e+ + Partition negTracksMC = o2::aod::mcparticle::pdgCode == +11; // e- + 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(); + + auto posTracks_per_coll = posTracksMC->sliceByCachedUnsorted(o2::aod::emmcparticle::emmceventId, mccollision.globalIndex(), cache); + auto negTracks_per_coll = negTracksMC->sliceByCachedUnsorted(o2::aod::emmcparticle::emmceventId, mccollision.globalIndex(), cache); + + for (auto& [t1, t2] : combinations(CombinationsFullIndexPolicy(posTracks_per_coll, negTracks_per_coll))) { // ULS + // LOGF(info, "pdg1 = %d, pdg2 = %d", t1.pdgCode(), t2.pdgCode()); + + if (!isInAcceptance(t1) || !isInAcceptance(t2)) { + continue; + } + + if (!t1.isPhysicalPrimary() && !t1.producedByGenerator()) { + continue; + } + if (!t2.isPhysicalPrimary() && !t2.producedByGenerator()) { + continue; + } + + int mother_id = FindLF(t1, t2, mcparticles); + int hfee_type = IsHF(t1, t2, mcparticles); + if (mother_id < 0 && hfee_type < 0) { + continue; + } + ROOT::Math::PtEtaPhiMVector v1(t1.pt(), t1.eta(), t1.phi(), o2::constants::physics::MassElectron); + ROOT::Math::PtEtaPhiMVector v2(t2.pt(), t2.eta(), t2.phi(), o2::constants::physics::MassElectron); + ROOT::Math::PtEtaPhiMVector v12 = v1 + v2; + + if (abs(v12.Rapidity()) > maxY) { + continue; + } + + if (mother_id > -1) { + auto mcmother = mcparticles.iteratorAt(mother_id); + if (mcmother.isPhysicalPrimary() || mcmother.producedByGenerator()) { + + switch (abs(mcmother.pdgCode())) { + case 111: + fRegistry.fill(HIST("Generated/sm/Pi0/hMvsPt"), v12.M(), v12.Pt()); + break; + case 221: + fRegistry.fill(HIST("Generated/sm/Eta/hMvsPt"), v12.M(), v12.Pt()); + break; + case 331: + fRegistry.fill(HIST("Generated/sm/EtaPrime/hMvsPt"), v12.M(), v12.Pt()); + break; + case 113: + fRegistry.fill(HIST("Generated/sm/Rho/hMvsPt"), v12.M(), v12.Pt()); + break; + case 223: + fRegistry.fill(HIST("Generated/sm/Omega/hMvsPt"), v12.M(), v12.Pt()); + break; + case 333: + fRegistry.fill(HIST("Generated/sm/Phi/hMvsPt"), v12.M(), v12.Pt()); + break; + default: + break; + } + } + } + } // end of true ULS pair loop + } // end of collision loop + } + PROCESS_SWITCH(vpPairQCMC, processGen, "run genrated info", true); + + void processDummy(MyCollisions const&) {} + PROCESS_SWITCH(vpPairQCMC, processDummy, "Dummy function", false); +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{ + adaptAnalysisTask(cfgc, TaskName{"vp-pair-qc-mc"})}; +} diff --git a/PWGEM/PhotonMeson/Core/V0PhotonCut.h b/PWGEM/PhotonMeson/Core/V0PhotonCut.h index 3f0e2ede17d..55e9718e156 100644 --- a/PWGEM/PhotonMeson/Core/V0PhotonCut.h +++ b/PWGEM/PhotonMeson/Core/V0PhotonCut.h @@ -61,7 +61,7 @@ class V0PhotonCut : public TNamed kDCAz, kITSNCls, kITSChi2NDF, - kITSCluserSize, + kITSClusterSize, kIsWithinBeamPipe, kRequireITSTPC, kRequireITSonly, @@ -208,7 +208,7 @@ class V0PhotonCut : public TNamed if (!IsSelectedTrack(track, V0PhotonCuts::kITSChi2NDF)) { return false; } - if (!IsSelectedTrack(track, V0PhotonCuts::kITSCluserSize)) { + if (!IsSelectedTrack(track, V0PhotonCuts::kITSClusterSize)) { return false; } return true; @@ -395,7 +395,7 @@ class V0PhotonCut : public TNamed case V0PhotonCuts::kITSChi2NDF: return mMinChi2PerClusterITS < track.itsChi2NCl() && track.itsChi2NCl() < mMaxChi2PerClusterITS; - case V0PhotonCuts::kITSCluserSize: { + case V0PhotonCuts::kITSClusterSize: { if (!isITSonlyTrack(track)) { return true; }