diff --git a/PWGEM/Dilepton/Core/DielectronCut.h b/PWGEM/Dilepton/Core/DielectronCut.h index 16f50e79179..ef5e47fbe60 100644 --- a/PWGEM/Dilepton/Core/DielectronCut.h +++ b/PWGEM/Dilepton/Core/DielectronCut.h @@ -40,6 +40,7 @@ class DielectronCut : public TNamed public: DielectronCut() = default; DielectronCut(const char* name, const char* title) : TNamed(name, title) {} + ~DielectronCut() {} enum class DielectronCuts : int { // pair cut diff --git a/PWGEM/Dilepton/Core/DimuonCut.h b/PWGEM/Dilepton/Core/DimuonCut.h index 07b9cb5d6f3..8af3161d670 100644 --- a/PWGEM/Dilepton/Core/DimuonCut.h +++ b/PWGEM/Dilepton/Core/DimuonCut.h @@ -27,6 +27,9 @@ #include "Framework/Logger.h" #include "Framework/DataTypes.h" #include "CommonConstants/PhysicsConstants.h" +#include "PWGEM/Dilepton/Utils/EMTrackUtilities.h" + +using namespace o2::aod::pwgem::dilepton::utils::emtrackutil; class DimuonCut : public TNamed { @@ -34,6 +37,8 @@ class DimuonCut : public TNamed DimuonCut() = default; DimuonCut(const char* name, const char* title) : TNamed(name, title) {} + ~DimuonCut() {} + enum class DimuonCuts : int { // pair cut kMass = 0, @@ -79,36 +84,9 @@ class DimuonCut : public TNamed ROOT::Math::PtEtaPhiMVector v2(t2.pt(), t2.eta(), t2.phi(), o2::constants::physics::MassMuon); ROOT::Math::PtEtaPhiMVector v12 = v1 + v2; - float dcaX1 = t1.fwdDcaX(); - float dcaY1 = t1.fwdDcaY(); - float cXX1 = t1.cXX(); - float cYY1 = t1.cYY(); - float cXY1 = t1.cXY(); - float dcaX2 = t2.fwdDcaX(); - float dcaY2 = t2.fwdDcaY(); - float cXX2 = t2.cXX(); - float cYY2 = t2.cYY(); - float cXY2 = t2.cXY(); - - float dca_xy1 = 999.f; - float det1 = cXX1 * cYY1 - cXY1 * cXY1; - if (det1 < 0) { - dca_xy1 = 999.f; - } else { - float chi2 = (dcaX1 * dcaX1 * cYY1 + dcaY1 * dcaY1 * cXX1 - 2. * dcaX1 * dcaY1 * cXY1) / det1; - dca_xy1 = std::sqrt(std::abs(chi2) / 2.); // in sigma - } - - float dca_xy2 = 999.f; - float det2 = cXX2 * cYY2 - cXY2 * cXY2; - if (det2 < 0) { - dca_xy2 = 999.f; - } else { - float chi2 = (dcaX2 * dcaX2 * cYY2 + dcaY2 * dcaY2 * cXX2 - 2. * dcaX2 * dcaY2 * cXY2) / det2; - dca_xy2 = std::sqrt(std::abs(chi2) / 2.); // in sigma - } - - float pair_dca_xy = std::sqrt((std::pow(dca_xy1, 2) + std::pow(dca_xy2, 2)) / 2.); + float dca_xy_t1 = fwdDcaXYinSigma(t1); + float dca_xy_t2 = fwdDcaXYinSigma(t2); + float pair_dca_xy = std::sqrt((dca_xy_t1 * dca_xy_t1 + dca_xy_t2 * dca_xy_t2) / 2.); if (v12.M() < mMinMass || mMaxMass < v12.M()) { return false; diff --git a/PWGEM/Dilepton/Core/EMEventCut.h b/PWGEM/Dilepton/Core/EMEventCut.h index db2d8d803c6..1950414c0c4 100644 --- a/PWGEM/Dilepton/Core/EMEventCut.h +++ b/PWGEM/Dilepton/Core/EMEventCut.h @@ -27,6 +27,7 @@ class EMEventCut : public TNamed public: EMEventCut() = default; EMEventCut(const char* name, const char* title) : TNamed(name, title) {} + ~EMEventCut() {} enum class EMEventCuts : int { kSel8 = 0, diff --git a/PWGEM/Dilepton/Core/PhotonHBT.h b/PWGEM/Dilepton/Core/PhotonHBT.h index d25bacb045b..81cd48c0b44 100644 --- a/PWGEM/Dilepton/Core/PhotonHBT.h +++ b/PWGEM/Dilepton/Core/PhotonHBT.h @@ -207,6 +207,10 @@ struct PhotonHBT { used_photonIds.shrink_to_fit(); used_dileptonIds.clear(); used_dileptonIds.shrink_to_fit(); + + if (eid_bdt) { + delete eid_bdt; + } } HistogramRegistry fRegistry{"output", {}, OutputObjHandlingPolicy::AnalysisObject, false, false}; @@ -298,12 +302,13 @@ struct PhotonHBT { const AxisSpec axis_kt{ConfKtBins, "k_{T} (GeV/c)"}; const AxisSpec axis_qinv{30, 0.0, +0.3, "q_{inv} (GeV/c)"}; - const AxisSpec axis_qout_cms{60, -0.3, +0.3, "q_{out}^{CMS} (GeV/c)"}; - const AxisSpec axis_qside_cms{60, -0.3, +0.3, "q_{side}^{CMS} (GeV/c)"}; - const AxisSpec axis_qlong_cms{60, -0.3, +0.3, "q_{long}^{CMS} (GeV/c)"}; - const AxisSpec axis_qlong_lcms{60, -0.3, +0.3, "q_{long}^{LCMS} (GeV/c)"}; + const AxisSpec axis_qabs_lcms{30, 0.0, +0.3, "|q|^{LCMS} (GeV/c)"}; + const AxisSpec axis_qout{60, -0.3, +0.3, "q_{out} (GeV/c)"}; // qout does not change between LAB and LCMS frame + const AxisSpec axis_qside{60, -0.3, +0.3, "q_{side} (GeV/c)"}; // qside does not change between LAB and LCMS frame + const AxisSpec axis_qlong{60, -0.3, +0.3, "q_{long} (GeV/c)"}; - fRegistry.add("Pair/same/hs", "diphoton correlation", kTHnSparseD, {axis_kt, axis_qinv, axis_qout_cms, axis_qside_cms, axis_qlong_cms, axis_qlong_lcms}, true); + fRegistry.add("Pair/same/hs_1d", "diphoton correlation 1D", kTHnSparseD, {axis_kt, axis_qinv, axis_qabs_lcms}, true); + fRegistry.add("Pair/same/hs_3d", "diphoton correlation 3D LCMS", kTHnSparseD, {axis_kt, axis_qout, axis_qside, axis_qlong}, true); fRegistry.addClone("Pair/same/", "Pair/mix/"); } @@ -376,6 +381,7 @@ struct PhotonHBT { } } + o2::ml::OnnxModel* eid_bdt = nullptr; void DefineDileptonCut() { fDielectronCut = DielectronCut("fDielectronCut", "fDielectronCut"); @@ -414,7 +420,7 @@ struct PhotonHBT { 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(); + eid_bdt = new o2::ml::OnnxModel(); if (dielectroncuts.loadModelsFromCCDB) { ccdbApi.init(ccdburl); std::map metadata; @@ -453,14 +459,12 @@ struct PhotonHBT { ROOT::Math::PtEtaPhiMVector k12 = 0.5 * (v1 + v2); float qinv = -q12.M(); float kt = k12.Pt(); - float qlong_cms = q12.Pz(); - ROOT::Math::XYZVector q_3d = q12.Vect(); // 3D q vector + // ROOT::Math::XYZVector q_3d = q12.Vect(); // 3D q vector ROOT::Math::XYZVector uv_out(k12.Px() / k12.Pt(), k12.Py() / k12.Pt(), 0); // unit vector for out. i.e. parallel to kt ROOT::Math::XYZVector uv_long(0, 0, 1); // unit vector for long, beam axis ROOT::Math::XYZVector uv_side = uv_out.Cross(uv_long); // unit vector for side - float qout_cms = q_3d.Dot(uv_out); - float qside_cms = q_3d.Dot(uv_side); + // float qlong_lab = q_3d.Dot(uv_long); // longitudinally co-moving system (LCMS) ROOT::Math::PxPyPzEVector v1_cartesian(v1.Px(), v1.Py(), v1.Pz(), v1.E()); @@ -469,17 +473,28 @@ struct PhotonHBT { float beta_z = (v1 + v2).Pz() / (v1 + v2).E(); ROOT::Math::Boost bst_z(0, 0, -beta_z); // Boost supports only PxPyPzEVector ROOT::Math::PxPyPzEVector q12_lcms = bst_z(q12_cartesian); - float qlong_lcms = q12_lcms.Pz(); + ROOT::Math::XYZVector q_3d_lcms = q12_lcms.Vect(); // 3D q vector in LCMS + float qout_lcms = q_3d_lcms.Dot(uv_out); + float qside_lcms = q_3d_lcms.Dot(uv_side); + float qlong_lcms = q_3d_lcms.Dot(uv_long); + float qabs_lcms = q_3d_lcms.R(); // ROOT::Math::PxPyPzEVector v1_lcms_cartesian = bst_z(v1_cartesian); // ROOT::Math::PxPyPzEVector v2_lcms_cartesian = bst_z(v2_cartesian); // ROOT::Math::PxPyPzEVector q12_lcms_cartesian = bst_z(q12_cartesian); - // LOGF(info, "q12.Pz() = %f, q12_cartesian.Pz() = %f",q12.Pz(), q12_cartesian.Pz()); - // LOGF(info, "v1.Pz() = %f, v2.Pz() = %f",v1.Pz(), v2.Pz()); - // LOGF(info, "v1_lcms_cartesian.Pz() = %f, v2_lcms_cartesian.Pz() = %f",v1_lcms_cartesian.Pz(), v2_lcms_cartesian.Pz()); + // LOGF(info, "q12.Pz() = %f, q12_cartesian.Pz() = %f", q12.Pz(), q12_cartesian.Pz()); + // LOGF(info, "v1.Pz() = %f, v2.Pz() = %f", v1.Pz(), v2.Pz()); + // LOGF(info, "v1_lcms_cartesian.Pz() = %f, v2_lcms_cartesian.Pz() = %f", v1_lcms_cartesian.Pz(), v2_lcms_cartesian.Pz()); // LOGF(info, "q12_lcms_cartesian.Pz() = %f", q12_lcms_cartesian.Pz()); - - fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("hs"), kt, qinv, qout_cms, qside_cms, qlong_cms, qlong_lcms); + // LOGF(info, "q_3d_lcms.Dot(uv_out) = %f, q_3d_lcms.Dot(uv_side) = %f, q_3d.Dot(uv_out) = %f, q_3d.Dot(uv_side) = %f", q_3d_lcms.Dot(uv_out), q_3d_lcms.Dot(uv_side), q_3d.Dot(uv_out), q_3d.Dot(uv_side)); + // LOGF(info, "q12_lcms.Pz() = %f, q_3d_lcms.Dot(uv_long) = %f", q12_lcms.Pz(), q_3d_lcms.Dot(uv_long)); + // ROOT::Math::PxPyPzEVector q12_lcms_tmp = bst_z(v1_cartesian) - bst_z(v2_cartesian); + // LOGF(info, "q12_lcms.Px() = %f, q12_lcms.Py() = %f, q12_lcms.Pz() = %f, q12_lcms_tmp.Px() = %f, q12_lcms_tmp.Py() = %f, q12_lcms_tmp.Pz() = %f", q12_lcms.Px(), q12_lcms.Py(), q12_lcms.Pz(), q12_lcms_tmp.Px(), q12_lcms_tmp.Py(), q12_lcms_tmp.Pz()); + // float qabs_lcms_tmp = q12_lcms.P(); + // LOGF(info, "qabs_lcms = %f, qabs_lcms_tmp = %f", qabs_lcms, qabs_lcms_tmp); + + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("hs_1d"), kt, qinv, qabs_lcms); + fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("hs_3d"), kt, qout_lcms, qside_lcms, qlong_lcms); } template diff --git a/PWGEM/Dilepton/DataModel/dileptonTables.h b/PWGEM/Dilepton/DataModel/dileptonTables.h index 392fdd1df49..47071b770a3 100644 --- a/PWGEM/Dilepton/DataModel/dileptonTables.h +++ b/PWGEM/Dilepton/DataModel/dileptonTables.h @@ -9,6 +9,8 @@ // granted to it by virtue of its status as an Intergovernmental Organization // or submit itself to any jurisdiction. +#include +#include #include "Common/Core/RecoDecay.h" #include "Framework/AnalysisDataModel.h" #include "Common/DataModel/PIDResponse.h" @@ -25,9 +27,39 @@ namespace o2::aod { +namespace pwgem::dilepton::swt +{ +enum class swtAliases : int { // software trigger aliases for EM + kUnDef = 0, + kHighTrackMult = 1, + kHighFt0Mult, + kSingleE, + kLMeeIMR, + kLMeeHMR, + kDiElectron, + kSingleMuLow, + kSingleMuHigh, + kDiMuon, + kNaliases +}; + +const std::unordered_map aliasLabels = { + {"fHighTrackMult", static_cast(o2::aod::pwgem::dilepton::swt::swtAliases::kHighTrackMult)}, + {"fHighFt0Mult", static_cast(o2::aod::pwgem::dilepton::swt::swtAliases::kHighFt0Mult)}, + {"fSingleE", static_cast(o2::aod::pwgem::dilepton::swt::swtAliases::kSingleE)}, + {"fLMeeIMR", static_cast(o2::aod::pwgem::dilepton::swt::swtAliases::kLMeeIMR)}, + {"fLMeeHMR", static_cast(o2::aod::pwgem::dilepton::swt::swtAliases::kLMeeHMR)}, + {"fDiElectron", static_cast(o2::aod::pwgem::dilepton::swt::swtAliases::kDiElectron)}, + {"fSingleMuLow", static_cast(o2::aod::pwgem::dilepton::swt::swtAliases::kSingleMuLow)}, + {"fSingleMuHigh", static_cast(o2::aod::pwgem::dilepton::swt::swtAliases::kSingleMuHigh)}, + {"fDiMuon", static_cast(o2::aod::pwgem::dilepton::swt::swtAliases::kDiMuon)}, +}; +} // namespace pwgem::dilepton::swt + namespace emevent { DECLARE_SOA_COLUMN(CollisionId, collisionId, int); +DECLARE_SOA_BITMAP_COLUMN(SWTAlias, swtalias, 16); //! Bitmask of fired trigger aliases (see above for definitions) DECLARE_SOA_COLUMN(NeeULS, neeuls, int); DECLARE_SOA_COLUMN(NeeLSpp, neelspp, int); DECLARE_SOA_COLUMN(NeeLSmm, neelsmm, int); @@ -134,6 +166,10 @@ DECLARE_SOA_TABLE(EMEventsQvec, "AOD", "EMEVENTQVEC", //! event q vector table emevent::EP3BTot); using EMEventQvec = EMEventsQvec::iterator; +DECLARE_SOA_TABLE(EMSWTriggerBits, "AOD", "EMSWTRIGGERBIT", //! + emevent::SWTAlias); +using EMSWTriggerBit = EMSWTriggerBits::iterator; + DECLARE_SOA_TABLE(EMEventsNee, "AOD", "EMEVENTNEE", emevent::NeeULS, emevent::NeeLSpp, emevent::NeeLSmm); // joinable to EMEvents using EMEventNee = EMEventsNee::iterator; diff --git a/PWGEM/Dilepton/TableProducer/CMakeLists.txt b/PWGEM/Dilepton/TableProducer/CMakeLists.txt index 14c74a69c03..9ca8e020b29 100644 --- a/PWGEM/Dilepton/TableProducer/CMakeLists.txt +++ b/PWGEM/Dilepton/TableProducer/CMakeLists.txt @@ -47,7 +47,7 @@ o2physics_add_dpl_workflow(skimmer-secondary-electron o2physics_add_dpl_workflow(create-emevent-dilepton SOURCES createEMEventDilepton.cxx - PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore O2Physics::EventFilteringUtils COMPONENT_NAME Analysis) o2physics_add_dpl_workflow(associate-mc-info-dilepton diff --git a/PWGEM/Dilepton/TableProducer/createEMEventDilepton.cxx b/PWGEM/Dilepton/TableProducer/createEMEventDilepton.cxx index 8d49ee134f5..2df68a64baf 100644 --- a/PWGEM/Dilepton/TableProducer/createEMEventDilepton.cxx +++ b/PWGEM/Dilepton/TableProducer/createEMEventDilepton.cxx @@ -24,6 +24,7 @@ #include "DataFormatsParameters/GRPObject.h" #include "DataFormatsParameters/GRPMagField.h" #include "CCDB/BasicCCDBManager.h" +#include "EventFiltering/Zorro.h" #include "PWGEM/Dilepton/DataModel/dileptonTables.h" #include "PWGEM/PhotonMeson/DataModel/gammaTables.h" @@ -50,6 +51,7 @@ struct CreateEMEventDilepton { Produces event_mult; Produces event_cent; Produces event_qvec; + Produces emswtbit; enum class EMEventType : int { kEvent = 0, @@ -64,6 +66,8 @@ struct CreateEMEventDilepton { Configurable skipGRPOquery{"skipGRPOquery", true, "skip grpo query"}; Configurable d_bz_input{"d_bz", -999, "bz field, -999 is automatic"}; Configurable applyEveSel_at_skimming{"applyEveSel_at_skimming", false, "flag to apply minimal event selection at the skimming level"}; + Configurable enable_swt{"enable_swt", false, "flag to process skimmed data (swt triggered)"}; + Configurable cfg_swt_names{"cfg_swt_names", "fHighTrackMult,fHighFt0Mult", "comma-separated software trigger names"}; // !trigger names have to be pre-registered in dileptonTable.h for bit operation! HistogramRegistry registry{"registry"}; void init(o2::framework::InitContext&) @@ -73,6 +77,14 @@ struct CreateEMEventDilepton { hEventCounter->GetXaxis()->SetBinLabel(2, "sel8"); } + ~CreateEMEventDilepton() + { + swt_names.clear(); + swt_names.shrink_to_fit(); + } + + Zorro zorro; + std::vector swt_names; int mRunNumber; float d_bz; Service ccdb; @@ -84,6 +96,16 @@ struct CreateEMEventDilepton { return; } + if (enable_swt) { + LOGF(info, "enable software triggers : %s", cfg_swt_names.value.data()); + zorro.initCCDB(ccdb.service, bc.runNumber(), bc.timestamp(), cfg_swt_names.value); + std::stringstream tokenizer(cfg_swt_names.value); + std::string token; + while (std::getline(tokenizer, token, ',')) { + swt_names.emplace_back(token); + } + } + // In case override, don't proceed, please - no CCDB access required if (d_bz_input > -990) { d_bz = d_bz_input; @@ -137,13 +159,27 @@ struct CreateEMEventDilepton { continue; } + if (enable_swt) { + if (zorro.isSelected(bc.globalBC())) { // triggered event + uint16_t trigger_bitmap = 0; + for (auto& swtname : swt_names) { + LOGF(info, "swtname = %s , swt index = %d", swtname.data(), o2::aod::pwgem::dilepton::swt::aliasLabels.at(swtname)); + trigger_bitmap |= BIT(o2::aod::pwgem::dilepton::swt::aliasLabels.at(swtname)); + } // end of desired trigger names loop + LOGF(info, "trigger_bitmap = %d", trigger_bitmap); + emswtbit(trigger_bitmap); + } else { // rejected + continue; + } + } + // LOGF(info, "collision.multNTracksPV() = %d, collision.multFT0A() = %f, collision.multFT0C() = %f", collision.multNTracksPV(), collision.multFT0A(), collision.multFT0C()); + registry.fill(HIST("hEventCounter"), 1); if (collision.sel8()) { registry.fill(HIST("hEventCounter"), 2); } - // uint64_t tag = collision.selection_raw(); event(collision.globalIndex(), bc.runNumber(), bc.globalBC(), collision.alias_raw(), collision.selection_raw(), bc.timestamp(), collision.posX(), collision.posY(), collision.posZ(), collision.numContrib(), collision.trackOccupancyInTimeRange()); diff --git a/PWGEM/Dilepton/Tasks/dielectronQC.cxx b/PWGEM/Dilepton/Tasks/dielectronQC.cxx index 8f2a89fdc3c..73562dacf36 100644 --- a/PWGEM/Dilepton/Tasks/dielectronQC.cxx +++ b/PWGEM/Dilepton/Tasks/dielectronQC.cxx @@ -277,6 +277,10 @@ struct dielectronQC { used_trackIds.clear(); used_trackIds.shrink_to_fit(); + + if (eid_bdt) { + delete eid_bdt; + } } void addhistograms() @@ -412,6 +416,7 @@ struct dielectronQC { fEMEventCut.SetOccupancyRange(eventcuts.cfgOccupancyMin, eventcuts.cfgOccupancyMax); } + o2::ml::OnnxModel* eid_bdt = nullptr; void DefineDileptonCut() { fDielectonCut = DielectronCut("fDielectonCut", "fDielectonCut"); @@ -448,7 +453,7 @@ struct dielectronQC { 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 - o2::ml::OnnxModel* eid_bdt = new o2::ml::OnnxModel(); + eid_bdt = new o2::ml::OnnxModel(); if (dielectroncuts.loadModelsFromCCDB) { ccdbApi.init(ccdburl); std::map metadata; diff --git a/PWGEM/Dilepton/Tasks/dielectronQCMC.cxx b/PWGEM/Dilepton/Tasks/dielectronQCMC.cxx index 67b381b21cb..0104ffcbc53 100644 --- a/PWGEM/Dilepton/Tasks/dielectronQCMC.cxx +++ b/PWGEM/Dilepton/Tasks/dielectronQCMC.cxx @@ -157,7 +157,12 @@ struct dielectronQCMC { 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() {} + ~dielectronQCMC() + { + if (eid_bdt) { + delete eid_bdt; + } + } void addhistograms() { @@ -288,7 +293,6 @@ struct dielectronQCMC { float beamP1 = 0.f; // beam momentum float beamP2 = 0.f; // beam momentum - bool cfgDoFlow = false; void init(InitContext&) { DefineEMEventCut(); @@ -382,6 +386,7 @@ struct dielectronQCMC { fEMEventCut.SetOccupancyRange(eventcuts.cfgOccupancyMin, eventcuts.cfgOccupancyMax); } + o2::ml::OnnxModel* eid_bdt = nullptr; void DefineDileptonCut() { fDielectonCut = DielectronCut("fDielectonCut", "fDielectonCut"); @@ -418,7 +423,8 @@ struct dielectronQCMC { 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 - o2::ml::OnnxModel* eid_bdt = new o2::ml::OnnxModel(); + // o2::ml::OnnxModel* eid_bdt = new o2::ml::OnnxModel(); + eid_bdt = new o2::ml::OnnxModel(); if (dielectroncuts.loadModelsFromCCDB) { ccdbApi.init(ccdburl); std::map metadata; @@ -645,7 +651,7 @@ struct dielectronQCMC { fRegistry.fill(HIST("Pair/sm/Phi/hMvsPhiV"), phiv, v12.M()); fillTrackInfo<0, TMCParticles>(t1); fillTrackInfo<0, TMCParticles>(t2); - if (mcmother.daughtersIds().size() == 2) { // omeag->ee + 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()); } @@ -841,13 +847,13 @@ struct dielectronQCMC { continue; } - o2::aod::pwgem::dilepton::utils::eventhistogram::fillEventInfo<0>(&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>(&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 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); @@ -955,9 +961,6 @@ struct dielectronQCMC { fRegistry.fill(HIST("Generated/sm/Omega/hs"), v12.M(), v12.Pt(), abs(dphi), abs(cos_thetaCS), abs(phiCS), aco, asym, abs(dphi_e_ee)); if (mcmother.daughtersIds().size() == 2) { // omega->ee fRegistry.fill(HIST("Generated/sm/Omega2ee/hs"), v12.M(), v12.Pt(), abs(dphi), abs(cos_thetaCS), abs(phiCS), aco, asym, abs(dphi_e_ee)); - // float mt = std::sqrt(std::pow(v12.M(), 2) + std::pow(v12.Pt(),2)); - // float cos_thetaCS_byhand = 2.f * (v1.E() * v2.Pz() - v2.E() * v1.Pz()) / (v12.M() * mt); - // LOGF(info, "cos_thetaCS = %f, cos_thetaCS_byhand = %f", cos_thetaCS, cos_thetaCS_byhand); } break; case 333: @@ -1193,6 +1196,11 @@ struct dielectronQCMC { auto mctracks_per_coll = mcparticles.sliceBy(perMcCollision, mccollision.globalIndex()); for (auto& mctrack : mctracks_per_coll) { + + if (!(mctrack.isPhysicalPrimary() || mctrack.producedByGenerator())) { + continue; + } + if (!(mctrack.isPhysicalPrimary() || mctrack.producedByGenerator()) || abs(mctrack.y()) > maxY) { continue; } diff --git a/PWGEM/Dilepton/Tasks/dimuonQCMC.cxx b/PWGEM/Dilepton/Tasks/dimuonQCMC.cxx index d1a037b1741..375314f7c56 100644 --- a/PWGEM/Dilepton/Tasks/dimuonQCMC.cxx +++ b/PWGEM/Dilepton/Tasks/dimuonQCMC.cxx @@ -142,6 +142,8 @@ struct dimuonQCMC { const AxisSpec axis_pt{ConfPtmumuBins, "p_{T,#mu#mu} (GeV/c)"}; const AxisSpec axis_dca{ConfDCAmumuBins, "DCA_{#mu#mu}^{xy} (#sigma)"}; // const AxisSpec axis_pca{ConfPCAmumuBins, "PCA (mm)"}; // particle closest approach + const AxisSpec axis_pt_meson{ConfPtmumuBins, "p_{T} (GeV/c)"}; // for omega, phi meson pT spectra + const AxisSpec axis_y_meson{25, -4.5, -2.0, "y"}; // rapidity of meson const AxisSpec axis_dphi_ee{18, 0, M_PI, "#Delta#varphi = #varphi_{#mu1} - #varphi_{#mu2} (rad.)"}; // for kHFll const AxisSpec axis_cos_theta_cs{10, 0.f, 1.f, "|cos(#theta_{CS})|"}; // for kPolarization @@ -162,6 +164,10 @@ struct dimuonQCMC { fRegistry.addClone("Generated/sm/Eta/", "Generated/sm/NonPromptJPsi/"); fRegistry.addClone("Generated/sm/Eta/", "Generated/sm/PromptPsi2S/"); fRegistry.addClone("Generated/sm/Eta/", "Generated/sm/NonPromptPsi2S/"); + fRegistry.add("Generated/sm/Omega2mumu/hPt", "pT of #omega meson", kTH1F, {axis_pt_meson}, true); + fRegistry.add("Generated/sm/Omega2mumu/hY", "rapidity of #omega meson", kTH1F, {axis_y_meson}, true); + fRegistry.add("Generated/sm/Phi2mumu/hPt", "pT of #phi meson", kTH1F, {axis_pt_meson}, true); + fRegistry.add("Generated/sm/Phi2mumu/hY", "rapidity of #phi meson", kTH1F, {axis_y_meson}, true); fRegistry.add("Generated/ccbar/c2mu_c2mu/hadron_hadron/hs", "m_{#mu#mu} vs. p_{T,#mu#mu}", kTHnSparseF, {axis_mass, axis_pt, axis_dphi_ee, axis_cos_theta_cs, axis_phi_cs, axis_aco, axis_asym_pt, axis_dphi_e_ee}, true); fRegistry.addClone("Generated/ccbar/c2mu_c2mu/hadron_hadron/", "Generated/ccbar/c2mu_c2mu/meson_meson/"); @@ -230,7 +236,6 @@ struct dimuonQCMC { fRegistry.addClone("Track/lf/", "Track/b2c2mu/"); } - bool cfgDoFlow = false; void init(InitContext&) { DefineEMEventCut(); @@ -706,13 +711,13 @@ struct dimuonQCMC { continue; } - o2::aod::pwgem::dilepton::utils::eventhistogram::fillEventInfo<0>(&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>(&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 posTracks_per_coll = posTracks->sliceByCached(o2::aod::emprimarymuon::emeventId, collision.globalIndex(), cache); auto negTracks_per_coll = negTracks->sliceByCached(o2::aod::emprimarymuon::emeventId, collision.globalIndex(), cache); @@ -815,9 +820,6 @@ struct dimuonQCMC { fRegistry.fill(HIST("Generated/sm/Omega/hs"), v12.M(), v12.Pt(), abs(dphi), abs(cos_thetaCS), abs(phiCS), aco, asym, abs(dphi_e_ee)); if (mcmother.daughtersIds().size() == 2) { // omega->ee fRegistry.fill(HIST("Generated/sm/Omega2mumu/hs"), v12.M(), v12.Pt(), abs(dphi), abs(cos_thetaCS), abs(phiCS), aco, asym, abs(dphi_e_ee)); - // float mt = std::sqrt(std::pow(v12.M(), 2) + std::pow(v12.Pt(),2)); - // float cos_thetaCS_byhand = 2.f * (v1.E() * v2.Pz() - v2.E() * v1.Pz()) / (v12.M() * mt); - // LOGF(info, "cos_thetaCS = %f, cos_thetaCS_byhand = %f", cos_thetaCS, cos_thetaCS_byhand); } break; case 333: @@ -1038,6 +1040,46 @@ struct dimuonQCMC { } // end of true LS++ pair loop } // end of collision loop + + // for oemga, phi efficiency + 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 mctracks_per_coll = mcparticles.sliceBy(perMcCollision, mccollision.globalIndex()); + + for (auto& mctrack : mctracks_per_coll) { + + if (!(mctrack.isPhysicalPrimary() || mctrack.producedByGenerator())) { + continue; + } + + if (!(mctrack.isPhysicalPrimary() || mctrack.producedByGenerator()) || (mctrack.y() < minY || maxY < mctrack.y())) { + continue; + } + switch (abs(mctrack.pdgCode())) { + case 223: + fRegistry.fill(HIST("Generated/sm/Omega2mumu/hPt"), mctrack.pt()); + fRegistry.fill(HIST("Generated/sm/Omega2mumu/hY"), mctrack.y()); + break; + case 333: + fRegistry.fill(HIST("Generated/sm/Phi2mumu/hPt"), mctrack.pt()); + fRegistry.fill(HIST("Generated/sm/Phi2mumu/hY"), mctrack.y()); + break; + default: + break; + } + + } // end of mctracks per mccollision + + } // end of collision loop } PROCESS_SWITCH(dimuonQCMC, processGen, "run genrated info", true); diff --git a/PWGEM/Dilepton/Tasks/singleElectronQCMC.cxx b/PWGEM/Dilepton/Tasks/singleElectronQCMC.cxx index b5e7fddacb8..32ee54a7c69 100644 --- a/PWGEM/Dilepton/Tasks/singleElectronQCMC.cxx +++ b/PWGEM/Dilepton/Tasks/singleElectronQCMC.cxx @@ -136,7 +136,12 @@ struct singleElectronQCMC { 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/"}; - ~singleElectronQCMC() {} + ~singleElectronQCMC() + { + if (eid_bdt) { + delete eid_bdt; + } + } void addhistograms() { @@ -271,6 +276,7 @@ struct singleElectronQCMC { fEMEventCut.SetOccupancyRange(eventcuts.cfgOccupancyMin, eventcuts.cfgOccupancyMax); } + o2::ml::OnnxModel* eid_bdt = nullptr; void DefineDileptonCut() { fDielectonCut = DielectronCut("fDielectonCut", "fDielectonCut"); @@ -300,7 +306,7 @@ struct singleElectronQCMC { 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 - o2::ml::OnnxModel* eid_bdt = new o2::ml::OnnxModel(); + eid_bdt = new o2::ml::OnnxModel(); if (dielectroncuts.loadModelsFromCCDB) { ccdbApi.init(ccdburl); std::map metadata;