From 0ef0af10eef986686d6a49e5477d82b8dfaa98ba Mon Sep 17 00:00:00 2001 From: Daiki Sekihata Date: Tue, 21 Oct 2025 14:58:59 +0200 Subject: [PATCH 1/2] PWGEM/Dilepton: add a task to evaluate pair acc. in MC --- PWGEM/Dilepton/Core/Dilepton.h | 6 +- PWGEM/Dilepton/Core/DileptonMC.h | 68 ++- .../TableProducer/associateMCinfoDilepton.cxx | 8 +- PWGEM/Dilepton/Tasks/CMakeLists.txt | 4 + PWGEM/Dilepton/Tasks/evaluateAcceptance.cxx | 407 ++++++++++++++++++ PWGEM/Dilepton/Utils/MCUtilities.h | 36 +- 6 files changed, 485 insertions(+), 44 deletions(-) create mode 100644 PWGEM/Dilepton/Tasks/evaluateAcceptance.cxx diff --git a/PWGEM/Dilepton/Core/Dilepton.h b/PWGEM/Dilepton/Core/Dilepton.h index 4b8c4658b2d..ac355b6c698 100644 --- a/PWGEM/Dilepton/Core/Dilepton.h +++ b/PWGEM/Dilepton/Core/Dilepton.h @@ -131,8 +131,8 @@ struct Dilepton { // ConfigurableAxis ConfMmumuBins{"ConfMmumuBins", {VARIABLE_WIDTH, 0.20, 0.21, 0.22, 0.23, 0.24, 0.25, 0.26, 0.27, 0.28, 0.29, 0.30, 0.31, 0.32, 0.33, 0.34, 0.35, 0.36, 0.37, 0.38, 0.39, 0.40, 0.41, 0.42, 0.43, 0.44, 0.45, 0.46, 0.47, 0.48, 0.49, 0.50, 0.51, 0.52, 0.53, 0.54, 0.55, 0.56, 0.57, 0.58, 0.59, 0.60, 0.61, 0.62, 0.63, 0.64, 0.65, 0.66, 0.67, 0.68, 0.69, 0.70, 0.71, 0.72, 0.73, 0.74, 0.75, 0.76, 0.77, 0.78, 0.79, 0.80, 0.81, 0.82, 0.83, 0.84, 0.85, 0.86, 0.87, 0.88, 0.89, 0.90, 0.91, 0.92, 0.93, 0.94, 0.95, 0.96, 0.97, 0.98, 0.99, 1.00, 1.01, 1.02, 1.03, 1.04, 1.05, 1.06, 1.07, 1.08, 1.09, 1.10, 1.11,1.12,1.13,1.14,1.15,1.16,1.17,1.18,1.19, 1.20, 1.30, 1.40, 1.50, 1.60, 1.70, 1.80, 1.90, 2.00, 2.10, 2.20, 2.30, 2.40, 2.50, 2.60, 2.70, 2.75, 2.80, 2.85, 2.90, 2.95, 3.00, 3.05, 3.10, 3.15, 3.20, 3.25, 3.30, 3.35, 3.40, 3.45, 3.50, 3.55, 3.60, 3.65, 3.70, 3.75, 3.80, 3.85, 3.90, 3.95, 4.00, 4.10, 4.20, 4.30, 4.40, 4.50, 4.60, 4.70, 4.80, 4.90, 5.00, 5.10, 5.20, 5.30, 5.40, 5.50, 5.60, 5.70, 5.80, 5.90, 6.00, 6.10, 6.20, 6.30, 6.40, 6.50, 6.60, 6.70, 6.80, 6.90, 7.00, 7.10, 7.20, 7.30, 7.40, 7.50, 7.60, 7.70, 7.80, 7.90, 8.00, 8.10, 8.20, 8.30, 8.40, 8.50, 8.60, 8.70, 8.80, 8.90, 9.00, 9.10, 9.20, 9.30, 9.40, 9.50, 9.60, 9.70, 9.80, 9.90, 10.00, 10.10, 10.20, 10.30, 10.40, 10.50, 10.60, 10.70, 10.80, 10.90, 11.00, 11.50, 12.00}, "mmumu bins for output histograms"}; // for dimuon. one can copy bins here to hyperloop page. ConfigurableAxis ConfSPBins{"ConfSPBins", {200, -5, 5}, "SP bins for flow analysis"}; - ConfigurableAxis ConfPolarizationPhiBins{"ConfPolarizationPhiBins", {1, 0.f, 2 * M_PI}, "phi bins for polarization analysis"}; ConfigurableAxis ConfPolarizationCosThetaBins{"ConfPolarizationCosThetaBins", {20, -1.f, 1.f}, "cos(theta) bins for polarization analysis"}; + ConfigurableAxis ConfPolarizationPhiBins{"ConfPolarizationPhiBins", {1, -M_PI, M_PI}, "phi bins for polarization analysis"}; ConfigurableAxis ConfPolarizationQuadMomBins{"ConfPolarizationQuadMomBins", {15, -0.5, 1}, "quadrupole moment bins for polarization analysis"}; // quardrupole moment <(3 x cos^2(theta) -1)/2> EMEventCut fEMEventCut; @@ -973,7 +973,7 @@ struct Dilepton { } else if (cfgPolarizationFrame == 1) { o2::aod::pwgem::dilepton::utils::pairutil::getAngleHX(std::array{t1.px(), t1.py(), t1.pz(), leptonM1}, std::array{t2.px(), t2.py(), t2.pz(), leptonM2}, beamE1, beamE2, beamP1, beamP2, t1.sign(), cos_thetaPol, phiPol); } - o2::math_utils::bringTo02Pi(phiPol); + o2::math_utils::bringToPMPi(phiPol); if (t1.sign() * t2.sign() < 0) { // ULS fRegistry.fill(HIST("Pair/") + HIST(event_pair_types[ev_id]) + HIST("uls/hs"), v12.M(), v12.Pt(), pair_dca, v12.Rapidity(), aco, asym, std::fabs(dphi_e_ee), cos_thetaPol, weight); @@ -1030,7 +1030,7 @@ struct Dilepton { } else if (cfgPolarizationFrame == 1) { o2::aod::pwgem::dilepton::utils::pairutil::getAngleHX(std::array{t1.px(), t1.py(), t1.pz(), leptonM1}, std::array{t2.px(), t2.py(), t2.pz(), leptonM2}, beamE1, beamE2, beamP1, beamP2, t1.sign(), cos_thetaPol, phiPol); } - o2::math_utils::bringTo02Pi(phiPol); + o2::math_utils::bringToPMPi(phiPol); float quadmom = (3.f * std::pow(cos_thetaPol, 2) - 1.f) / 2.f; if (t1.sign() * t2.sign() < 0) { // ULS diff --git a/PWGEM/Dilepton/Core/DileptonMC.h b/PWGEM/Dilepton/Core/DileptonMC.h index 51de312162f..1c01522934e 100644 --- a/PWGEM/Dilepton/Core/DileptonMC.h +++ b/PWGEM/Dilepton/Core/DileptonMC.h @@ -127,8 +127,8 @@ struct DileptonMC { Configurable cfg_nbin_aco{"cfg_nbin_aco", 1, "number of bins for acoplanarity"}; // 10 Configurable cfg_nbin_asym_pt{"cfg_nbin_asym_pt", 1, "number of bins for pt asymmetry"}; // 10 Configurable cfg_nbin_dphi_e_ee{"cfg_nbin_dphi_e_ee", 1, "number of bins for dphi_ee_e"}; // 18 - ConfigurableAxis ConfPolarizationPhiBins{"ConfPolarizationPhiBins", {1, 0.f, 2 * M_PI}, "phi bins for polarization analysis"}; ConfigurableAxis ConfPolarizationCosThetaBins{"ConfPolarizationCosThetaBins", {1, -1.f, 1.f}, "cos(theta) bins for polarization analysis"}; + ConfigurableAxis ConfPolarizationPhiBins{"ConfPolarizationPhiBins", {1, -M_PI, M_PI}, "phi bins for polarization analysis"}; Configurable cfgPolarizationFrame{"cfgPolarizationFrame", 0, "frame of polarization. 0:CS, 1:HX, else:FATAL"}; ConfigurableAxis ConfPolarizationQuadMomBins{"ConfPolarizationQuadMomBins", {1, -0.5, 1}, "quadrupole moment bins for polarization analysis"}; // quardrupole moment <(3 x cos^2(theta) -1)/2> @@ -803,27 +803,57 @@ struct DileptonMC { } template - int FindLF(TTrack const& posmc, TTrack const& negmc, TMCParticles const& mcparticles) + int FindSMULS(TTrack const& t1mc, TTrack const& t2mc, TMCParticles const& mcparticles) { int arr[] = { - FindCommonMotherFrom2Prongs(posmc, negmc, -pdg_lepton, pdg_lepton, 22, mcparticles), - FindCommonMotherFrom2Prongs(posmc, negmc, -pdg_lepton, pdg_lepton, 111, mcparticles), - FindCommonMotherFrom2Prongs(posmc, negmc, -pdg_lepton, pdg_lepton, 221, mcparticles), - FindCommonMotherFrom2Prongs(posmc, negmc, -pdg_lepton, pdg_lepton, 331, mcparticles), - FindCommonMotherFrom2Prongs(posmc, negmc, -pdg_lepton, pdg_lepton, 113, mcparticles), - FindCommonMotherFrom2Prongs(posmc, negmc, -pdg_lepton, pdg_lepton, 223, mcparticles), - FindCommonMotherFrom2Prongs(posmc, negmc, -pdg_lepton, pdg_lepton, 333, mcparticles), - FindCommonMotherFrom2Prongs(posmc, negmc, -pdg_lepton, pdg_lepton, 443, mcparticles), - FindCommonMotherFrom2Prongs(posmc, negmc, -pdg_lepton, pdg_lepton, 100443, mcparticles) - // FindCommonMotherFrom2Prongs(posmc, negmc, -pdg_lepton, pdg_lepton, 553, mcparticles), - // FindCommonMotherFrom2Prongs(posmc, negmc, -pdg_lepton, pdg_lepton, 100553, mcparticles), - // FindCommonMotherFrom2Prongs(posmc, negmc, -pdg_lepton, pdg_lepton, 200553, mcparticles) + FindCommonMotherFrom2Prongs(t1mc, t2mc, -pdg_lepton, pdg_lepton, 22, mcparticles), + FindCommonMotherFrom2Prongs(t1mc, t2mc, -pdg_lepton, pdg_lepton, 111, mcparticles), + FindCommonMotherFrom2Prongs(t1mc, t2mc, -pdg_lepton, pdg_lepton, 221, mcparticles), + FindCommonMotherFrom2Prongs(t1mc, t2mc, -pdg_lepton, pdg_lepton, 331, mcparticles), + FindCommonMotherFrom2Prongs(t1mc, t2mc, -pdg_lepton, pdg_lepton, 113, mcparticles), + FindCommonMotherFrom2Prongs(t1mc, t2mc, -pdg_lepton, pdg_lepton, 223, mcparticles), + FindCommonMotherFrom2Prongs(t1mc, t2mc, -pdg_lepton, pdg_lepton, 333, mcparticles), + FindCommonMotherFrom2Prongs(t1mc, t2mc, -pdg_lepton, pdg_lepton, 443, mcparticles), + FindCommonMotherFrom2Prongs(t1mc, t2mc, -pdg_lepton, pdg_lepton, 100443, mcparticles) + // FindCommonMotherFrom2Prongs(t1mc, t2mc, -pdg_lepton, pdg_lepton, 553, mcparticles), + // FindCommonMotherFrom2Prongs(t1mc, t2mc, -pdg_lepton, pdg_lepton, 100553, mcparticles), + // FindCommonMotherFrom2Prongs(t1mc, t2mc, -pdg_lepton, pdg_lepton, 200553, mcparticles) }; int size = sizeof(arr) / sizeof(*arr); int max = *std::max_element(arr, arr + size); return max; } + template + int FindSMLSPP(TTrack const& t1mc, TTrack const& t2mc, TMCParticles const& mcparticles) + { + int arr[] = { + FindCommonMotherFrom2Prongs(t1mc, t2mc, -pdg_lepton, -pdg_lepton, 111, mcparticles), + FindCommonMotherFrom2Prongs(t1mc, t2mc, -pdg_lepton, -pdg_lepton, 221, mcparticles), + FindCommonMotherFrom2Prongs(t1mc, t2mc, -pdg_lepton, -pdg_lepton, 331, mcparticles), + FindCommonMotherFrom2Prongs(t1mc, t2mc, -pdg_lepton, -pdg_lepton, 113, mcparticles), + FindCommonMotherFrom2Prongs(t1mc, t2mc, -pdg_lepton, -pdg_lepton, 223, mcparticles), + FindCommonMotherFrom2Prongs(t1mc, t2mc, -pdg_lepton, -pdg_lepton, 333, mcparticles)}; + int size = sizeof(arr) / sizeof(*arr); + int max = *std::max_element(arr, arr + size); + return max; + } + + template + int FindSMLSMM(TTrack const& t1mc, TTrack const& t2mc, TMCParticles const& mcparticles) + { + int arr[] = { + FindCommonMotherFrom2Prongs(t1mc, t2mc, pdg_lepton, pdg_lepton, 111, mcparticles), + FindCommonMotherFrom2Prongs(t1mc, t2mc, pdg_lepton, pdg_lepton, 221, mcparticles), + FindCommonMotherFrom2Prongs(t1mc, t2mc, pdg_lepton, pdg_lepton, 331, mcparticles), + FindCommonMotherFrom2Prongs(t1mc, t2mc, pdg_lepton, pdg_lepton, 113, mcparticles), + FindCommonMotherFrom2Prongs(t1mc, t2mc, pdg_lepton, pdg_lepton, 223, mcparticles), + FindCommonMotherFrom2Prongs(t1mc, t2mc, pdg_lepton, pdg_lepton, 333, mcparticles)}; + int size = sizeof(arr) / sizeof(*arr); + int max = *std::max_element(arr, arr + size); + return max; + } + template bool isInAcceptance(T const& lepton) { @@ -1490,7 +1520,7 @@ struct DileptonMC { } else if (cfgPolarizationFrame == 1) { o2::aod::pwgem::dilepton::utils::pairutil::getAngleHX(arrP1, arrP2, beamE1, beamE2, beamP1, beamP2, t1.sign(), cos_thetaPol, phiPol); } - o2::math_utils::bringTo02Pi(phiPol); + o2::math_utils::bringToPMPi(phiPol); float quadmom = (3.f * std::pow(cos_thetaPol, 2) - 1.f) / 2.f; if ((FindCommonMotherFrom2ProngsWithoutPDG(t1mc, t2mc) > 0 || IsHF(t1mc, t2mc, mcparticles) > 0) && is_pair_from_same_mcevent) { // for bkg study @@ -1533,7 +1563,7 @@ struct DileptonMC { if (cfgRequireTrueAssociation && (t1mc.emmceventId() != collision.emmceventId() || t2mc.emmceventId() != collision.emmceventId())) { return false; } - int mother_id = FindLF(t1mc, t2mc, mcparticles); + int mother_id = std::max({FindSMULS(t1mc, t2mc, mcparticles), FindSMULS(t2mc, t1mc, mcparticles), FindSMLSPP(t1mc, t2mc, mcparticles), FindSMLSMM(t1mc, t2mc, mcparticles)}); int hfee_type = IsHF(t1mc, t2mc, mcparticles); if (mother_id < 0 && hfee_type < 0) { return false; @@ -1708,7 +1738,7 @@ struct DileptonMC { return false; } - int mother_id = FindLF(t1, t2, mcparticles); + int mother_id = std::max({FindSMULS(t1, t2, mcparticles), FindSMULS(t2, t1, mcparticles), FindSMLSPP(t1, t2, mcparticles), FindSMLSMM(t1, t2, mcparticles)}); int hfee_type = IsHF(t1, t2, mcparticles); if (mother_id < 0 && hfee_type < 0) { return false; @@ -1797,7 +1827,7 @@ struct DileptonMC { } else if (cfgPolarizationFrame == 1) { o2::aod::pwgem::dilepton::utils::pairutil::getAngleHX(arrP1, arrP2, beamE1, beamE2, beamP1, beamP2, -t1.pdgCode() / pdg_lepton, cos_thetaPol, phiPol); } - o2::math_utils::bringTo02Pi(phiPol); + o2::math_utils::bringToPMPi(phiPol); float quadmom = (3.f * std::pow(cos_thetaPol, 2) - 1.f) / 2.f; // bool isInAcc = isInAcceptance(t1) && isInAcceptance(t2); @@ -2226,7 +2256,7 @@ struct DileptonMC { if (!((t1mc.isPhysicalPrimary() || t1mc.producedByGenerator()) && (t2mc.isPhysicalPrimary() || t2mc.producedByGenerator()))) { return false; } - int mother_id = FindLF(t1mc, t2mc, mcparticles); + int mother_id = std::max({FindSMULS(t1mc, t2mc, mcparticles), FindSMULS(t2mc, t1mc, mcparticles), FindSMLSPP(t1mc, t2mc, mcparticles), FindSMLSMM(t1mc, t2mc, mcparticles)}); int hfee_type = IsHF(t1mc, t2mc, mcparticles); if (mother_id < 0 && hfee_type < 0) { return false; diff --git a/PWGEM/Dilepton/TableProducer/associateMCinfoDilepton.cxx b/PWGEM/Dilepton/TableProducer/associateMCinfoDilepton.cxx index 0d62389cdfc..92078356fd6 100644 --- a/PWGEM/Dilepton/TableProducer/associateMCinfoDilepton.cxx +++ b/PWGEM/Dilepton/TableProducer/associateMCinfoDilepton.cxx @@ -60,10 +60,10 @@ struct AssociateMCInfoDilepton { Configurable n_dummy_loop{"n_dummy_loop", 0, "for loop runs over n times"}; Configurable down_scaling_omega{"down_scaling_omega", 1.1, "down scaling factor to store omega"}; Configurable down_scaling_phi{"down_scaling_phi", 1.1, "down scaling factor to store phi"}; - Configurable min_eta_gen_primary{"min_eta_gen_primary", -1.5, "min rapidity Y to store generated information"}; // smearing is applied at analysis stage. set wider value. - Configurable max_eta_gen_primary{"max_eta_gen_primary", +1.5, "max rapidity Y to store generated information"}; // smearing is applied at analysis stage. set wider value. - Configurable min_eta_gen_primary_fwd{"min_eta_gen_primary_fwd", -5.0, "min eta to store generated information"}; // smearing is applied at analysis stage. set wider value. - Configurable max_eta_gen_primary_fwd{"max_eta_gen_primary_fwd", -1.5, "max eta to store generated information"}; // smearing is applied at analysis stage. set wider value. + Configurable min_eta_gen_primary{"min_eta_gen_primary", -1.5, "min eta to store generated information"}; // smearing is applied at analysis stage. set wider value. + Configurable max_eta_gen_primary{"max_eta_gen_primary", +1.5, "max eta to store generated information"}; // smearing is applied at analysis stage. set wider value. + Configurable min_eta_gen_primary_fwd{"min_eta_gen_primary_fwd", -6.0, "min eta to store generated information"}; // smearing is applied at analysis stage. set wider value. + Configurable max_eta_gen_primary_fwd{"max_eta_gen_primary_fwd", -1.0, "max eta to store generated information"}; // smearing is applied at analysis stage. set wider value. HistogramRegistry registry{"EMMCEvent"}; std::mt19937 engine; diff --git a/PWGEM/Dilepton/Tasks/CMakeLists.txt b/PWGEM/Dilepton/Tasks/CMakeLists.txt index 4e75843e48d..24cf04752b3 100644 --- a/PWGEM/Dilepton/Tasks/CMakeLists.txt +++ b/PWGEM/Dilepton/Tasks/CMakeLists.txt @@ -171,3 +171,7 @@ o2physics_add_dpl_workflow(study-dcafitter PUBLIC_LINK_LIBRARIES O2::Framework O2::DCAFitter O2Physics::AnalysisCore COMPONENT_NAME Analysis) +o2physics_add_dpl_workflow(evaluate-acceptance + SOURCES evaluateAcceptance.cxx + PUBLIC_LINK_LIBRARIES O2::Framework O2Physics::AnalysisCore + COMPONENT_NAME Analysis) diff --git a/PWGEM/Dilepton/Tasks/evaluateAcceptance.cxx b/PWGEM/Dilepton/Tasks/evaluateAcceptance.cxx new file mode 100644 index 00000000000..b11b73ec5c9 --- /dev/null +++ b/PWGEM/Dilepton/Tasks/evaluateAcceptance.cxx @@ -0,0 +1,407 @@ +// 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. + +/// \file evaluateAcceptance.cxx +/// \brief a task to evaluate pair acceptance in MC +/// \author daiki.sekihata@cern.ch + +#include "PWGEM/Dilepton/Utils/MCUtilities.h" +#include "PWGEM/Dilepton/Utils/PairUtilities.h" + +#include "Common/Core/RecoDecay.h" +#include "Common/Core/TableHelper.h" + +// #include "Common/Core/trackUtilities.h" +// #include "Common/DataModel/Centrality.h" +// #include "Common/DataModel/CollisionAssociationTables.h" +// #include "Common/DataModel/EventSelection.h" +// #include "Common/DataModel/Multiplicity.h" +// #include "Common/DataModel/PIDResponse.h" +// #include "Common/DataModel/TrackSelectionTables.h" + +#include "CCDB/BasicCCDBManager.h" +#include "CommonConstants/LHCConstants.h" +#include "CommonConstants/PhysicsConstants.h" +#include "DataFormatsParameters/GRPLHCIFData.h" + +// #include "DataFormatsCalibration/MeanVertexObject.h" +// #include "DataFormatsParameters/GRPMagField.h" +// #include "DataFormatsParameters/GRPObject.h" +// #include "DetectorsBase/GeometryManager.h" +// #include "DetectorsBase/Propagator.h" + +#include "Framework/AnalysisDataModel.h" +#include "Framework/AnalysisTask.h" +#include "Framework/runDataProcessing.h" + +#include "Math/Vector4D.h" + +#include +#include +#include +#include +#include +#include + +using namespace o2; +using namespace o2::soa; +using namespace o2::framework; +using namespace o2::framework::expressions; +using namespace o2::constants::physics; +using namespace o2::aod::pwgem::dilepton::utils::mcutil; +using namespace o2::aod::pwgem::dilepton::utils::pairutil; + +struct evaluateAcceptance { + // Configurables + Configurable ccdburl{"ccdb-url", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; + + Configurable cfgEventGeneratorType{"cfgEventGeneratorType", -1, "if positive, select event generator type. i.e. gap or signal"}; + Configurable cfgPdgLepton{"cfgPdgLepton", 11, "pdg code 11 or 13"}; + ConfigurableAxis ConfMllBins{"ConfMllBins", {400, 0, 4}, "mll bins"}; + ConfigurableAxis ConfPtllBins{"ConfPtllBins", {100, 0, 10}, "pTll bins"}; + ConfigurableAxis ConfYllBins{"ConfYllBins", {400, -10, +10}, "yll bins"}; + ConfigurableAxis ConfCosThetaBins{"ConfCosThetaBins", {40, -1, +1}, "cos theta bins for polarization"}; + ConfigurableAxis ConfPhiBins{"ConfPhiBins", {72, -M_PI, M_PI}, "phi bins for polarization"}; + ConfigurableAxis ConfQuadMomBins{"ConfQuadMomBins", {150, -0.5, 1}, "quadrupole moment bins for polarization"}; + ConfigurableAxis ConfPtlBins{"ConfPtlBins", {200, 0, 10}, "pTl bins"}; + ConfigurableAxis ConfEtalBins{"ConfEtalBins", {200, -10, 10}, "etal bins"}; + + HistogramRegistry fRegistry{"fRegistry"}; + Service ccdb; + int mRunNumber = 0; + + float beamM1 = o2::constants::physics::MassProton; // mass of beam + float beamM2 = o2::constants::physics::MassProton; // mass of beam + float beamE1 = 0.f; // beam energy + float beamE2 = 0.f; // beam energy + float beamP1 = 0.f; // beam momentum + float beamP2 = 0.f; // beam momentum + + float leptonM1 = 0.f; + float leptonM2 = 0.f; + void init(o2::framework::InitContext&) + { + ccdb->setURL(ccdburl); + ccdb->setCaching(true); + ccdb->setLocalObjectValidityChecking(); + ccdb->setFatalWhenNull(false); + + if (cfgPdgLepton.value == 11) { + leptonM1 = o2::constants::physics::MassElectron; + leptonM2 = o2::constants::physics::MassElectron; + } else if (cfgPdgLepton.value == 13) { + leptonM1 = o2::constants::physics::MassMuon; + leptonM2 = o2::constants::physics::MassMuon; + } else { + LOGF(fatal, "pdg code must be 11 or 13."); + } + + addHistograms(); + } + + template + void initCCDB(TBC const& bc) + { + if (mRunNumber == bc.runNumber()) { + return; + } + + auto grplhcif = ccdb->getForTimeStamp("GLO/Config/GRPLHCIF", bc.timestamp()); + int beamZ1 = grplhcif->getBeamZ(o2::constants::lhc::BeamC); + int beamZ2 = grplhcif->getBeamZ(o2::constants::lhc::BeamA); + int beamA1 = grplhcif->getBeamA(o2::constants::lhc::BeamC); + int beamA2 = grplhcif->getBeamA(o2::constants::lhc::BeamA); + beamE1 = grplhcif->getBeamEnergyPerNucleonInGeV(o2::constants::lhc::BeamC); + beamE2 = grplhcif->getBeamEnergyPerNucleonInGeV(o2::constants::lhc::BeamA); + beamM1 = o2::constants::physics::MassProton * beamA1; + beamM2 = o2::constants::physics::MassProton * beamA2; + beamP1 = std::sqrt(std::pow(beamE1, 2) - std::pow(beamM1, 2)); + beamP2 = std::sqrt(std::pow(beamE2, 2) - std::pow(beamM2, 2)); + LOGF(info, "beamZ1 = %d, beamZ2 = %d, beamA1 = %d, beamA2 = %d, beamE1 = %f (GeV), beamE2 = %f (GeV), beamM1 = %f (GeV), beamM2 = %f (GeV), beamP1 = %f (GeV), beamP2 = %f (GeV)", beamZ1, beamZ2, beamA1, beamA2, beamE1, beamE2, beamM1, beamM2, beamP1, beamP2); + mRunNumber = bc.runNumber(); + } + + static constexpr std::string_view pair_sign_types[3] = {"uls/", "lspp/", "lsmm/"}; + static constexpr std::string_view dilepton_source_types[20] = { + "sm/Pi0/", // 0 + "sm/Eta/", // 1 + "sm/EtaPrime/", // 2 + "sm/Rho/", // 3 + "sm/Omega/", // 4 + "sm/Omega2ll/", // 5 + "sm/Phi/", // 6 + "sm/Phi2ll/", // 7 + "sm/PromptJPsi/", // 8 + "sm/NonPromptJPsi/", // 9 + "sm/PromptPsi2S/", // 10 + "sm/NonPromptPsi2S/", // 11 + "sm/Upsilon1S/", // 12 + "sm/Upsilon2S/", // 13 + "sm/Upsilon3S/", // 14 + "ccbar/c2l_c2l/", // 15 + "bbbar/b2l_b2l/", // 16 + "bbbar/b2c2l_b2c2l/", // 17 + "bbbar/b2c2l_b2l_sameb/", // 18 + "bbbar/b2c2l_b2l_diffb/" // 19 + }; // unordered_map is better, but cannot be constexpr. + + void addHistograms() + { + auto hCollisionCounter = fRegistry.add("Event/hCollisionCounter", "collision counter", kTH1D, {{2, -0.5f, 1.5f}}, false); + hCollisionCounter->GetXaxis()->SetBinLabel(1, "all"); + hCollisionCounter->GetXaxis()->SetBinLabel(2, "accepted"); + + const AxisSpec axisMll{ConfMllBins, "m_{ll} (GeV/c^{2})"}; + const AxisSpec axisPtll{ConfPtllBins, "p_{T,ll} (GeV/c)"}; + const AxisSpec axisYll{ConfYllBins, "y_{ll}"}; + const AxisSpec axisCosThetaCS{ConfCosThetaBins, "cos(#theta^{CS})"}; + const AxisSpec axisPhiCS{ConfPhiBins, "#varphi^{CS} (rad.)"}; + const AxisSpec axisQuadMomCS{ConfQuadMomBins, "#frac{3 cos^{2}(#theta^{CS}) #minus 1}{2}"}; + const AxisSpec axisCosThetaHX{ConfCosThetaBins, "cos(#theta^{HX})"}; + const AxisSpec axisPhiHX{ConfPhiBins, "#varphi^{HX} (rad.)"}; + const AxisSpec axisQuadMomHX{ConfQuadMomBins, "#frac{3 cos^{2}(#theta^{HX}) #minus 1}{2}"}; + + const AxisSpec axisPtl1{ConfPtlBins, "p_{T,l1} (GeV/c)"}; + const AxisSpec axisPtl2{ConfPtlBins, "p_{T,l2} (GeV/c)"}; + const AxisSpec axisEtal1{ConfEtalBins, "#eta_{l1}"}; + const AxisSpec axisEtal2{ConfEtalBins, "#eta_{l2}"}; + + // for pairs + fRegistry.add("Generated/sm/Pi0/uls/hs", "gen. dilepton", kTHnSparseD, {axisMll, axisPtll, axisYll, axisCosThetaCS, axisPhiCS, axisQuadMomCS, axisCosThetaHX, axisPhiHX, axisQuadMomHX, axisPtl1, axisPtl2, axisEtal1, axisEtal2}, true); + fRegistry.addClone("Generated/sm/Pi0/uls/", "Generated/sm/Pi0/lspp/"); + fRegistry.addClone("Generated/sm/Pi0/uls/", "Generated/sm/Pi0/lsmm/"); + 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/Omega2ll/"); + fRegistry.addClone("Generated/sm/Pi0/", "Generated/sm/Phi/"); + fRegistry.addClone("Generated/sm/Pi0/", "Generated/sm/Phi2ll/"); + fRegistry.addClone("Generated/sm/Pi0/", "Generated/sm/PromptJPsi/"); + fRegistry.addClone("Generated/sm/Pi0/", "Generated/sm/NonPromptJPsi/"); + fRegistry.addClone("Generated/sm/Pi0/", "Generated/sm/PromptPsi2S/"); + fRegistry.addClone("Generated/sm/Pi0/", "Generated/sm/NonPromptPsi2S/"); + fRegistry.addClone("Generated/sm/Pi0/", "Generated/sm/Upsilon1S/"); + fRegistry.addClone("Generated/sm/Pi0/", "Generated/sm/Upsilon2S/"); + fRegistry.addClone("Generated/sm/Pi0/", "Generated/sm/Upsilon3S/"); + fRegistry.addClone("Generated/sm/Pi0/", "Generated/ccbar/c2l_c2l/"); + fRegistry.addClone("Generated/sm/Pi0/", "Generated/bbbar/b2l_b2l/"); + fRegistry.addClone("Generated/sm/Pi0/", "Generated/bbbar/b2c2l_b2c2l/"); + fRegistry.addClone("Generated/sm/Pi0/", "Generated/bbbar/b2c2l_b2l_sameb/"); // ULS + fRegistry.addClone("Generated/sm/Pi0/", "Generated/bbbar/b2c2l_b2l_diffb/"); // LS + } + + template + void fillGenPairInfo(TLepton const& t1, TLepton const& t2, TMCParticles const& mcParticles) + { + if (!t1.isPhysicalPrimary() && !t1.producedByGenerator()) { + return; + } + if (!t2.isPhysicalPrimary() && !t2.producedByGenerator()) { + return; + } + + int mother_id = std::max({FindSMULS(t1, t2, mcParticles), FindSMULS(t2, t1, mcParticles), FindSMLSPP(t1, t2, mcParticles), FindSMLSMM(t1, t2, mcParticles)}); + int hfee_type = IsHF(t1, t2, mcParticles); + if (mother_id < 0 && hfee_type < 0) { + return; + } + + ROOT::Math::PtEtaPhiMVector v1(t1.pt(), t1.eta(), t1.phi(), leptonM1); + ROOT::Math::PtEtaPhiMVector v2(t2.pt(), t2.eta(), t2.phi(), leptonM2); + ROOT::Math::PtEtaPhiMVector v12 = v1 + v2; + + int sign1 = t1.pdgCode() > 0 ? -1 : 1; + std::array arrP1 = {t1.px(), t1.py(), t1.pz(), leptonM1}; + std::array arrP2 = {t2.px(), t2.py(), t2.pz(), leptonM2}; + float cosThetaCS = 999, phiCS = 999.f; + float cosThetaHX = 999, phiHX = 999.f; + o2::aod::pwgem::dilepton::utils::pairutil::getAngleCS(arrP1, arrP2, beamE1, beamE2, beamP1, beamP2, sign1, cosThetaCS, phiCS); + o2::aod::pwgem::dilepton::utils::pairutil::getAngleHX(arrP1, arrP2, beamE1, beamE2, beamP1, beamP2, sign1, cosThetaHX, phiHX); + o2::math_utils::bringToPMPi(phiCS); + o2::math_utils::bringToPMPi(phiHX); + float quadmomCS = (3.f * std::pow(cosThetaCS, 2) - 1.f) / 2.f; + float quadmomHX = (3.f * std::pow(cosThetaHX, 2) - 1.f) / 2.f; + + if (mother_id > -1) { + auto mcmother = mcParticles.iteratorAt(mother_id); + int nd = mcmother.daughtersIds()[1] - mcmother.daughtersIds()[0] + 1; // number of daughters + switch (std::abs(mcmother.pdgCode())) { + case 111: + fRegistry.fill(HIST("Generated/sm/Pi0/") + HIST(pair_sign_types[pairSignId]) + HIST("hs"), v12.M(), v12.Pt(), v12.Rapidity(), cosThetaCS, phiCS, quadmomCS, cosThetaHX, phiHX, quadmomHX, t1.pt(), t2.pt(), t1.eta(), t2.eta()); + break; + case 221: + fRegistry.fill(HIST("Generated/sm/Eta/") + HIST(pair_sign_types[pairSignId]) + HIST("hs"), v12.M(), v12.Pt(), v12.Rapidity(), cosThetaCS, phiCS, quadmomCS, cosThetaHX, phiHX, quadmomHX, t1.pt(), t2.pt(), t1.eta(), t2.eta()); + break; + case 331: + fRegistry.fill(HIST("Generated/sm/EtaPrime/") + HIST(pair_sign_types[pairSignId]) + HIST("hs"), v12.M(), v12.Pt(), v12.Rapidity(), cosThetaCS, phiCS, quadmomCS, cosThetaHX, phiHX, quadmomHX, t1.pt(), t2.pt(), t1.eta(), t2.eta()); + break; + case 113: + fRegistry.fill(HIST("Generated/sm/Rho/") + HIST(pair_sign_types[pairSignId]) + HIST("hs"), v12.M(), v12.Pt(), v12.Rapidity(), cosThetaCS, phiCS, quadmomCS, cosThetaHX, phiHX, quadmomHX, t1.pt(), t2.pt(), t1.eta(), t2.eta()); + break; + case 223: + fRegistry.fill(HIST("Generated/sm/Omega/") + HIST(pair_sign_types[pairSignId]) + HIST("hs"), v12.M(), v12.Pt(), v12.Rapidity(), cosThetaCS, phiCS, quadmomCS, cosThetaHX, phiHX, quadmomHX, t1.pt(), t2.pt(), t1.eta(), t2.eta()); + if (nd == 2) { + fRegistry.fill(HIST("Generated/sm/Omega2ll/") + HIST(pair_sign_types[pairSignId]) + HIST("hs"), v12.M(), v12.Pt(), v12.Rapidity(), cosThetaCS, phiCS, quadmomCS, cosThetaHX, phiHX, quadmomHX, t1.pt(), t2.pt(), t1.eta(), t2.eta()); + } + break; + case 333: + fRegistry.fill(HIST("Generated/sm/Phi/") + HIST(pair_sign_types[pairSignId]) + HIST("hs"), v12.M(), v12.Pt(), v12.Rapidity(), cosThetaCS, phiCS, quadmomCS, cosThetaHX, phiHX, quadmomHX, t1.pt(), t2.pt(), t1.eta(), t2.eta()); + if (nd == 2) { + fRegistry.fill(HIST("Generated/sm/Phi2ll/") + HIST(pair_sign_types[pairSignId]) + HIST("hs"), v12.M(), v12.Pt(), v12.Rapidity(), cosThetaCS, phiCS, quadmomCS, cosThetaHX, phiHX, quadmomHX, t1.pt(), t2.pt(), t1.eta(), t2.eta()); + } + break; + case 443: + if (IsFromBeauty(mcmother, mcParticles) > 0) { + fRegistry.fill(HIST("Generated/sm/NonPromptJPsi/") + HIST(pair_sign_types[pairSignId]) + HIST("hs"), v12.M(), v12.Pt(), v12.Rapidity(), cosThetaCS, phiCS, quadmomCS, cosThetaHX, phiHX, quadmomHX, t1.pt(), t2.pt(), t1.eta(), t2.eta()); + } else { + fRegistry.fill(HIST("Generated/sm/PromptJPsi/") + HIST(pair_sign_types[pairSignId]) + HIST("hs"), v12.M(), v12.Pt(), v12.Rapidity(), cosThetaCS, phiCS, quadmomCS, cosThetaHX, phiHX, quadmomHX, t1.pt(), t2.pt(), t1.eta(), t2.eta()); + } + break; + case 100443: + if (IsFromBeauty(mcmother, mcParticles) > 0) { + fRegistry.fill(HIST("Generated/sm/NonPromptPsi2S/") + HIST(pair_sign_types[pairSignId]) + HIST("hs"), v12.M(), v12.Pt(), v12.Rapidity(), cosThetaCS, phiCS, quadmomCS, cosThetaHX, phiHX, quadmomHX, t1.pt(), t2.pt(), t1.eta(), t2.eta()); + } else { + fRegistry.fill(HIST("Generated/sm/PromptPsi2S/") + HIST(pair_sign_types[pairSignId]) + HIST("hs"), v12.M(), v12.Pt(), v12.Rapidity(), cosThetaCS, phiCS, quadmomCS, cosThetaHX, phiHX, quadmomHX, t1.pt(), t2.pt(), t1.eta(), t2.eta()); + } + break; + case 553: + fRegistry.fill(HIST("Generated/sm/Upsilon1S/") + HIST(pair_sign_types[pairSignId]) + HIST("hs"), v12.M(), v12.Pt(), v12.Rapidity(), cosThetaCS, phiCS, quadmomCS, cosThetaHX, phiHX, quadmomHX, t1.pt(), t2.pt(), t1.eta(), t2.eta()); + break; + case 100553: + fRegistry.fill(HIST("Generated/sm/Upsilon2S/") + HIST(pair_sign_types[pairSignId]) + HIST("hs"), v12.M(), v12.Pt(), v12.Rapidity(), cosThetaCS, phiCS, quadmomCS, cosThetaHX, phiHX, quadmomHX, t1.pt(), t2.pt(), t1.eta(), t2.eta()); + break; + case 200553: + fRegistry.fill(HIST("Generated/sm/Upsilon3S/") + HIST(pair_sign_types[pairSignId]) + HIST("hs"), v12.M(), v12.Pt(), v12.Rapidity(), cosThetaCS, phiCS, quadmomCS, cosThetaHX, phiHX, quadmomHX, t1.pt(), t2.pt(), t1.eta(), t2.eta()); + break; + default: + break; + } + } else if (hfee_type > -1) { + switch (hfee_type) { + case static_cast(EM_HFeeType::kCe_Ce): + fRegistry.fill(HIST("Generated/ccbar/c2l_c2l/") + HIST(pair_sign_types[pairSignId]) + HIST("hs"), v12.M(), v12.Pt(), v12.Rapidity(), cosThetaCS, phiCS, quadmomCS, cosThetaHX, phiHX, quadmomHX, t1.pt(), t2.pt(), t1.eta(), t2.eta()); + break; + case static_cast(EM_HFeeType::kBe_Be): + fRegistry.fill(HIST("Generated/bbbar/b2l_b2l/") + HIST(pair_sign_types[pairSignId]) + HIST("hs"), v12.M(), v12.Pt(), v12.Rapidity(), cosThetaCS, phiCS, quadmomCS, cosThetaHX, phiHX, quadmomHX, t1.pt(), t2.pt(), t1.eta(), t2.eta()); + break; + case static_cast(EM_HFeeType::kBCe_BCe): + fRegistry.fill(HIST("Generated/bbbar/b2c2l_b2c2l/") + HIST(pair_sign_types[pairSignId]) + HIST("hs"), v12.M(), v12.Pt(), v12.Rapidity(), cosThetaCS, phiCS, quadmomCS, cosThetaHX, phiHX, quadmomHX, t1.pt(), t2.pt(), t1.eta(), t2.eta()); + break; + case static_cast(EM_HFeeType::kBCe_Be_SameB): + fRegistry.fill(HIST("Generated/bbbar/b2c2l_b2l_sameb/") + HIST(pair_sign_types[pairSignId]) + HIST("hs"), v12.M(), v12.Pt(), v12.Rapidity(), cosThetaCS, phiCS, quadmomCS, cosThetaHX, phiHX, quadmomHX, t1.pt(), t2.pt(), t1.eta(), t2.eta()); + break; + case static_cast(EM_HFeeType::kBCe_Be_DiffB): + fRegistry.fill(HIST("Generated/bbbar/b2c2l_b2l_diffb/") + HIST(pair_sign_types[pairSignId]) + HIST("hs"), v12.M(), v12.Pt(), v12.Rapidity(), cosThetaCS, phiCS, quadmomCS, cosThetaHX, phiHX, quadmomHX, t1.pt(), t2.pt(), t1.eta(), t2.eta()); + break; + default: + break; + } + } + } + + template + int FindSMULS(TTrack const& t1mc, TTrack const& t2mc, TMCParticles const& mcParticles) + { + int arr[] = { + FindCommonMotherFrom2Prongs(t1mc, t2mc, -cfgPdgLepton, cfgPdgLepton, 111, mcParticles), + FindCommonMotherFrom2Prongs(t1mc, t2mc, -cfgPdgLepton, cfgPdgLepton, 221, mcParticles), + FindCommonMotherFrom2Prongs(t1mc, t2mc, -cfgPdgLepton, cfgPdgLepton, 331, mcParticles), + FindCommonMotherFrom2Prongs(t1mc, t2mc, -cfgPdgLepton, cfgPdgLepton, 113, mcParticles), + FindCommonMotherFrom2Prongs(t1mc, t2mc, -cfgPdgLepton, cfgPdgLepton, 223, mcParticles), + FindCommonMotherFrom2Prongs(t1mc, t2mc, -cfgPdgLepton, cfgPdgLepton, 333, mcParticles), + FindCommonMotherFrom2Prongs(t1mc, t2mc, -cfgPdgLepton, cfgPdgLepton, 443, mcParticles), + FindCommonMotherFrom2Prongs(t1mc, t2mc, -cfgPdgLepton, cfgPdgLepton, 100443, mcParticles), + FindCommonMotherFrom2Prongs(t1mc, t2mc, -cfgPdgLepton, cfgPdgLepton, 553, mcParticles), + FindCommonMotherFrom2Prongs(t1mc, t2mc, -cfgPdgLepton, cfgPdgLepton, 100553, mcParticles), + FindCommonMotherFrom2Prongs(t1mc, t2mc, -cfgPdgLepton, cfgPdgLepton, 200553, mcParticles)}; + int size = sizeof(arr) / sizeof(*arr); + int max = *std::max_element(arr, arr + size); + return max; + } + + template + int FindSMLSPP(TTrack const& t1mc, TTrack const& t2mc, TMCParticles const& mcParticles) + { + int arr[] = { + FindCommonMotherFrom2Prongs(t1mc, t2mc, -cfgPdgLepton, -cfgPdgLepton, 221, mcParticles), + FindCommonMotherFrom2Prongs(t1mc, t2mc, -cfgPdgLepton, -cfgPdgLepton, 331, mcParticles), + FindCommonMotherFrom2Prongs(t1mc, t2mc, -cfgPdgLepton, -cfgPdgLepton, 113, mcParticles), + FindCommonMotherFrom2Prongs(t1mc, t2mc, -cfgPdgLepton, -cfgPdgLepton, 223, mcParticles), + FindCommonMotherFrom2Prongs(t1mc, t2mc, -cfgPdgLepton, -cfgPdgLepton, 333, mcParticles)}; + int size = sizeof(arr) / sizeof(*arr); + int max = *std::max_element(arr, arr + size); + return max; + } + + template + int FindSMLSMM(TTrack const& t1mc, TTrack const& t2mc, TMCParticles const& mcParticles) + { + int arr[] = { + FindCommonMotherFrom2Prongs(t1mc, t2mc, cfgPdgLepton, cfgPdgLepton, 221, mcParticles), + FindCommonMotherFrom2Prongs(t1mc, t2mc, cfgPdgLepton, cfgPdgLepton, 331, mcParticles), + FindCommonMotherFrom2Prongs(t1mc, t2mc, cfgPdgLepton, cfgPdgLepton, 113, mcParticles), + FindCommonMotherFrom2Prongs(t1mc, t2mc, cfgPdgLepton, cfgPdgLepton, 223, mcParticles), + FindCommonMotherFrom2Prongs(t1mc, t2mc, cfgPdgLepton, cfgPdgLepton, 333, mcParticles)}; + int size = sizeof(arr) / sizeof(*arr); + int max = *std::max_element(arr, arr + size); + return max; + } + + SliceCache cache; + Preslice perMCCollision = o2::aod::mcparticle::mcCollisionId; + Partition posLeptons = o2::aod::mcparticle::pdgCode == -cfgPdgLepton; // l+ + Partition negLeptons = o2::aod::mcparticle::pdgCode == cfgPdgLepton; // l- + + void process(aod::McCollisions const& mcCollisions, aod::McParticles const& mcParticles, aod::BCsWithTimestamps const&) + { + for (const auto& mcCollision : mcCollisions) { + auto bc = mcCollision.template bc_as(); + initCCDB(bc); + + fRegistry.fill(HIST("Event/hCollisionCounter"), 0); + if (cfgEventGeneratorType >= 0 && mcCollision.getSubGeneratorId() != cfgEventGeneratorType) { + continue; + } + fRegistry.fill(HIST("Event/hCollisionCounter"), 1); + + auto posLeptons_per_coll = posLeptons->sliceByCached(o2::aod::mcparticle::mcCollisionId, mcCollision.globalIndex(), cache); + auto negLeptons_per_coll = negLeptons->sliceByCached(o2::aod::mcparticle::mcCollisionId, mcCollision.globalIndex(), cache); + // LOGF(info, "mcCollision.globalIndex() = %d, posLeptons_per_coll.size() = %d, negLeptons_per_coll.size() = %d", mcCollision.globalIndex(), posLeptons_per_coll.size(), negLeptons_per_coll.size()); + + for (const auto& [t1, t2] : combinations(CombinationsFullIndexPolicy(posLeptons_per_coll, negLeptons_per_coll))) { // ULS + if (!(t1.isPhysicalPrimary() || t1.producedByGenerator()) || !(t2.isPhysicalPrimary() || t2.producedByGenerator())) { + continue; + } + fillGenPairInfo<0>(t1, t2, mcParticles); + } // end of ULS pairing + + for (const auto& [t1, t2] : combinations(CombinationsStrictlyUpperIndexPolicy(posLeptons_per_coll, posLeptons_per_coll))) { // LS++ + if (!(t1.isPhysicalPrimary() || t1.producedByGenerator()) || !(t2.isPhysicalPrimary() || t2.producedByGenerator())) { + continue; + } + fillGenPairInfo<1>(t1, t2, mcParticles); + } // end of LS++ pairing + + for (const auto& [t1, t2] : combinations(CombinationsStrictlyUpperIndexPolicy(negLeptons_per_coll, negLeptons_per_coll))) { // LS-- + if (!(t1.isPhysicalPrimary() || t1.producedByGenerator()) || !(t2.isPhysicalPrimary() || t2.producedByGenerator())) { + continue; + } + fillGenPairInfo<2>(t1, t2, mcParticles); + } // end of LS++ pairing + + } // end of mc collision loop + } +}; +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{adaptAnalysisTask(cfgc, TaskName{"evaluate-acceptance"})}; +} diff --git a/PWGEM/Dilepton/Utils/MCUtilities.h b/PWGEM/Dilepton/Utils/MCUtilities.h index dd12d2576f2..27ee530e4f1 100644 --- a/PWGEM/Dilepton/Utils/MCUtilities.h +++ b/PWGEM/Dilepton/Utils/MCUtilities.h @@ -241,14 +241,14 @@ int IsFromBeauty(TMCParticle const& p, TMCParticles const& mcparticles) int motherid = p.mothersIds()[0]; // first mother index auto mp_tmp = mcparticles.iteratorAt(motherid); - if (abs(mp_tmp.pdgCode()) < 1e+9 && (std::to_string(abs(mp_tmp.pdgCode()))[std::to_string(abs(mp_tmp.pdgCode())).length() - 2] == '5' && std::to_string(abs(mp_tmp.pdgCode()))[std::to_string(abs(mp_tmp.pdgCode())).length() - 3] == '5') && abs(mp_tmp.pdgCode()) % 2 == 1) { + if (std::abs(mp_tmp.pdgCode()) < 1e+9 && (std::to_string(std::abs(mp_tmp.pdgCode()))[std::to_string(std::abs(mp_tmp.pdgCode())).length() - 2] == '5' && std::to_string(std::abs(mp_tmp.pdgCode()))[std::to_string(std::abs(mp_tmp.pdgCode())).length() - 3] == '5') && std::abs(mp_tmp.pdgCode()) % 2 == 1) { return -999; // reject bottomonia } while (motherid > -1) { if (motherid < mcparticles.size()) { // protect against bad mother indices. why is this needed? auto mp = mcparticles.iteratorAt(motherid); - if (abs(mp.pdgCode()) < 1e+9 && (std::to_string(abs(mp.pdgCode()))[std::to_string(abs(mp.pdgCode())).length() - 3] == '5' || std::to_string(abs(mp.pdgCode()))[std::to_string(abs(mp.pdgCode())).length() - 4] == '5')) { + if (std::abs(mp.pdgCode()) < 1e+9 && (std::to_string(std::abs(mp.pdgCode()))[std::to_string(std::abs(mp.pdgCode())).length() - 3] == '5' || std::to_string(std::abs(mp.pdgCode()))[std::to_string(std::abs(mp.pdgCode())).length() - 4] == '5')) { return motherid; } if (mp.has_mothers()) { @@ -274,13 +274,13 @@ int IsFromCharm(TMCParticle const& p, TMCParticles const& mcparticles) int motherid = p.mothersIds()[0]; // first mother index auto mp_tmp = mcparticles.iteratorAt(motherid); - if (abs(mp_tmp.pdgCode()) < 1e+9 && (std::to_string(abs(mp_tmp.pdgCode()))[std::to_string(abs(mp_tmp.pdgCode())).length() - 2] == '4' && std::to_string(abs(mp_tmp.pdgCode()))[std::to_string(abs(mp_tmp.pdgCode())).length() - 3] == '4') && abs(mp_tmp.pdgCode()) % 2 == 1) { + if (std::abs(mp_tmp.pdgCode()) < 1e+9 && (std::to_string(std::abs(mp_tmp.pdgCode()))[std::to_string(std::abs(mp_tmp.pdgCode())).length() - 2] == '4' && std::to_string(std::abs(mp_tmp.pdgCode()))[std::to_string(std::abs(mp_tmp.pdgCode())).length() - 3] == '4') && std::abs(mp_tmp.pdgCode()) % 2 == 1) { return -999; // reject bottomonia } while (motherid > -1) { if (motherid < mcparticles.size()) { // protect against bad mother indices. why is this needed? auto mp = mcparticles.iteratorAt(motherid); - if (abs(mp.pdgCode()) < 1e+9 && (std::to_string(abs(mp.pdgCode()))[std::to_string(abs(mp.pdgCode())).length() - 3] == '4' || std::to_string(abs(mp.pdgCode()))[std::to_string(abs(mp.pdgCode())).length() - 4] == '4')) { + if (std::abs(mp.pdgCode()) < 1e+9 && (std::to_string(std::abs(mp.pdgCode()))[std::to_string(std::abs(mp.pdgCode())).length() - 3] == '4' || std::to_string(std::abs(mp.pdgCode()))[std::to_string(std::abs(mp.pdgCode())).length() - 4] == '4')) { return motherid; } if (mp.has_mothers()) { @@ -395,8 +395,8 @@ int IsHF(TMCParticle1 const& p1, TMCParticle2 const& p2, TMCParticles const& mcp if (mid1 == mid2) { auto common_mp = mcparticles.iteratorAt(mid1); int mp_pdg = common_mp.pdgCode(); - bool is_mp_diquark = (1100 < abs(mp_pdg) && abs(mp_pdg) < 5600) && std::to_string(mp_pdg)[std::to_string(mp_pdg).length() - 2] == '0'; - if (!is_mp_diquark && abs(mp_pdg) < 1e+9 && (std::to_string(abs(mp_pdg))[std::to_string(abs(mp_pdg)).length() - 3] == '5' || std::to_string(abs(mp_pdg))[std::to_string(abs(mp_pdg)).length() - 4] == '5')) { + bool is_mp_diquark = (1100 < std::abs(mp_pdg) && std::abs(mp_pdg) < 5600) && std::to_string(mp_pdg)[std::to_string(mp_pdg).length() - 2] == '0'; + if (!is_mp_diquark && std::abs(mp_pdg) < 1e+9 && (std::to_string(std::abs(mp_pdg))[std::to_string(std::abs(mp_pdg)).length() - 3] == '5' || std::to_string(std::abs(mp_pdg))[std::to_string(std::abs(mp_pdg)).length() - 4] == '5')) { mothers_id1.clear(); mothers_pdg1.clear(); mothers_id2.clear(); @@ -417,8 +417,8 @@ int IsHF(TMCParticle1 const& p1, TMCParticle2 const& p2, TMCParticles const& mcp if (mid1 == mid2) { auto common_mp = mcparticles.iteratorAt(mid1); int mp_pdg = common_mp.pdgCode(); - bool is_mp_diquark = (1100 < abs(mp_pdg) && abs(mp_pdg) < 5600) && std::to_string(mp_pdg)[std::to_string(mp_pdg).length() - 2] == '0'; - if (!is_mp_diquark && abs(mp_pdg) < 1e+9 && (std::to_string(abs(mp_pdg))[std::to_string(abs(mp_pdg)).length() - 3] == '5' || std::to_string(abs(mp_pdg))[std::to_string(abs(mp_pdg)).length() - 4] == '5')) { + bool is_mp_diquark = (1100 < std::abs(mp_pdg) && std::abs(mp_pdg) < 5600) && std::to_string(mp_pdg)[std::to_string(mp_pdg).length() - 2] == '0'; + if (!is_mp_diquark && std::abs(mp_pdg) < 1e+9 && (std::to_string(std::abs(mp_pdg))[std::to_string(std::abs(mp_pdg)).length() - 3] == '5' || std::to_string(std::abs(mp_pdg))[std::to_string(std::abs(mp_pdg)).length() - 4] == '5')) { is_same_mother_found = true; } } @@ -466,7 +466,7 @@ int searchMothers(T& p, U& mcParticles, int pdg, bool equal) } else if (mothersids[1] < mothersids[0]) { allmothersids.push_back(mothersids[0]); allmothersids.push_back(mothersids[1]); - } else if ((80 < abs(o2::mcgenstatus::getGenStatusCode(p.statusCode())) && abs(o2::mcgenstatus::getGenStatusCode(p.statusCode())) < 90) || (100 < abs(o2::mcgenstatus::getGenStatusCode(p.statusCode())) && abs(o2::mcgenstatus::getGenStatusCode(p.statusCode())) < 110)) { // NOTE: THIS IS GENERATOR DEPENDENT AND WORKS ONLY FOR PYTHIA! + } else if ((80 < std::abs(o2::mcgenstatus::getGenStatusCode(p.statusCode())) && std::abs(o2::mcgenstatus::getGenStatusCode(p.statusCode())) < 90) || (100 < std::abs(o2::mcgenstatus::getGenStatusCode(p.statusCode())) && std::abs(o2::mcgenstatus::getGenStatusCode(p.statusCode())) < 110)) { // NOTE: THIS IS GENERATOR DEPENDENT AND WORKS ONLY FOR PYTHIA! for (int i = mothersids[0]; i <= mothersids[1]; i++) { allmothersids.push_back(i); } @@ -484,16 +484,16 @@ int searchMothers(T& p, U& mcParticles, int pdg, bool equal) for (int i : allmothersids) { auto mother = mcParticles.iteratorAt(i); int mpdg = mother.pdgCode(); - // if (abs(mpdg) == pdg && mpdg * p.pdgCode() > 0) { // check for quark - if (abs(mpdg) == pdg) { // check for quark to allow for beauty and charm + oscillation + // if (std::abs(mpdg) == pdg && mpdg * p.pdgCode() > 0) { // check for quark + if (std::abs(mpdg) == pdg) { // check for quark to allow for beauty and charm + oscillation if (quark_id > -1 || next_mother_id > -1) { // we already found a possible candidate in the list of mothers, so now we have (at least) two // LOG(warning) << "Flavour tracking is ambiguous. Stopping here."; return -1; } quark_id = i; - //} else if ((static_cast(abs(mpdg) / 100) == pdg || static_cast(abs(mpdg) / 1000) == pdg) && mpdg * p.pdgCode() > 0) { // check for other mothers with flavour content - } else if ((static_cast(abs(mpdg) / 100) == pdg || static_cast(abs(mpdg) / 1000) == pdg)) { // check for other mothers with flavour content to allow for beauty and charm - if (quark_id > -1 || next_mother_id > -1) { // we already found a possible candidate in the list of mothers, so now we have (at least) two + //} else if ((static_cast(std::abs(mpdg) / 100) == pdg || static_cast(std::abs(mpdg) / 1000) == pdg) && mpdg * p.pdgCode() > 0) { // check for other mothers with flavour content + } else if ((static_cast(std::abs(mpdg) / 100) == pdg || static_cast(std::abs(mpdg) / 1000) == pdg)) { // check for other mothers with flavour content to allow for beauty and charm + if (quark_id > -1 || next_mother_id > -1) { // we already found a possible candidate in the list of mothers, so now we have (at least) two // LOG(warning) << "Flavour tracking is ambiguous. Stopping here."; return -1; } @@ -555,7 +555,7 @@ bool checkFromSameQuarkPair(T& p1, T& p2, U& mcParticles, int pdg) template bool isCharmMeson(T const& track) { - if (400 < abs(track.pdgCode()) && abs(track.pdgCode()) < 500) { + if (400 < std::abs(track.pdgCode()) && std::abs(track.pdgCode()) < 500) { return true; } else { return false; @@ -565,7 +565,7 @@ bool isCharmMeson(T const& track) template bool isCharmBaryon(T const& track) { - if (4000 < abs(track.pdgCode()) && abs(track.pdgCode()) < 5000) { + if (4000 < std::abs(track.pdgCode()) && std::abs(track.pdgCode()) < 5000) { return true; } else { return false; @@ -575,7 +575,7 @@ bool isCharmBaryon(T const& track) template bool isBeautyMeson(T const& track) { - if (500 < abs(track.pdgCode()) && abs(track.pdgCode()) < 600) { + if (500 < std::abs(track.pdgCode()) && std::abs(track.pdgCode()) < 600) { return true; } else { return false; @@ -585,7 +585,7 @@ bool isBeautyMeson(T const& track) template bool isBeautyBaryon(T const& track) { - if (5000 < abs(track.pdgCode()) && abs(track.pdgCode()) < 6000) { + if (5000 < std::abs(track.pdgCode()) && std::abs(track.pdgCode()) < 6000) { return true; } else { return false; From d8ac3d75fc270fa3c5a66712aa3612fb7cddb79f Mon Sep 17 00:00:00 2001 From: Daiki Sekihata Date: Tue, 21 Oct 2025 16:29:03 +0200 Subject: [PATCH 2/2] Add new includes for algorithm, array, format, map, and string --- PWGEM/Dilepton/Core/DileptonMC.h | 1 + 1 file changed, 1 insertion(+) diff --git a/PWGEM/Dilepton/Core/DileptonMC.h b/PWGEM/Dilepton/Core/DileptonMC.h index 1c01522934e..bbbdaaca444 100644 --- a/PWGEM/Dilepton/Core/DileptonMC.h +++ b/PWGEM/Dilepton/Core/DileptonMC.h @@ -47,6 +47,7 @@ #include "Math/Vector4D.h" #include "TString.h" +#include #include #include #include