diff --git a/PWGLF/Tasks/Resonances/CMakeLists.txt b/PWGLF/Tasks/Resonances/CMakeLists.txt index 06d226e887c..f8a6669861e 100644 --- a/PWGLF/Tasks/Resonances/CMakeLists.txt +++ b/PWGLF/Tasks/Resonances/CMakeLists.txt @@ -316,7 +316,7 @@ o2physics_add_dpl_workflow(k892hadronphoton o2physics_add_dpl_workflow(k892hadronphotonbkg SOURCES k892hadronphotonBkg.cxx - PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2Physics::AnalysisCCDB + PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2Physics::AnalysisCCDB O2Physics::MLCore COMPONENT_NAME Analysis) o2physics_add_dpl_workflow(chk892li diff --git a/PWGLF/Tasks/Resonances/k892hadronphotonBkg.cxx b/PWGLF/Tasks/Resonances/k892hadronphotonBkg.cxx index 07e7490a3a4..08ffc2d1bbe 100644 --- a/PWGLF/Tasks/Resonances/k892hadronphotonBkg.cxx +++ b/PWGLF/Tasks/Resonances/k892hadronphotonBkg.cxx @@ -16,11 +16,13 @@ #include "PWGLF/DataModel/LFStrangenessMLTables.h" #include "PWGLF/DataModel/LFStrangenessPIDTables.h" #include "PWGLF/DataModel/LFStrangenessTables.h" +#include "PWGLF/Utils/ResonanceMlResponse.h" #include "Common/CCDB/EventSelectionParams.h" #include "Common/CCDB/ctpRateFetcher.h" #include "Common/Core/RecoDecay.h" #include "Common/DataModel/Centrality.h" +#include "Tools/ML/MlResponse.h" #include #include @@ -29,6 +31,7 @@ #include #include #include +#include #include #include #include @@ -37,19 +40,23 @@ #include #include +#include #include // IWYU pragma: keep (do not replace with Math/Vector4Dfwd.h) #include #include #include +#include #include #include #include #include #include +#include #include using namespace o2; +using namespace o2::ml; using namespace o2::framework; using namespace o2::framework::expressions; using std::array; @@ -73,7 +80,10 @@ enum BkgResonance { struct k892hadronphotonBkg { Service ccdb{}; + o2::ccdb::CcdbApi ccdbApi; ctpRateFetcher rateFetcher; + o2::analysis::ResonanceMlResponse mlResponse; + TRandom3 rotRng{12345}; // struct member; fixed seed for reproducibility across grid jobs // Histogram registry @@ -119,6 +129,22 @@ struct k892hadronphotonBkg { Configurable rotGamma{"rotGamma", false, "Flag to rotate the photon direction"}; } lstarBkgConfig; + struct : ConfigurableGroup { + std::string prefix = "bdt"; // JSON group name + Configurable ccdbUrl{"ccdbUrl", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; + Configurable> onnxFileNames{"onnxFileNames", std::vector{"BDTModel.onnx"}, "Local .onnx file names, one per pT bin"}; + Configurable> modelPathsCCDB{"modelPathsCCDB", std::vector{"Users/o/obenchik/MLModels/BDT"}, "Model paths on CCDB, one per pT bin (each model needs its own folder)"}; + Configurable timestampCCDB{"timestampCCDB", 1695750420200, "timestamp of the ONNX file for ML model used to query in CCDB. Please use 1695750420200"}; + 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)"}; + Configurable numThreads{"numThreads", 1, "ONNX intra-op threads. 0 lets ONNX Runtime default to one thread per physical core"}; + Configurable enableML{"enableML", false, "Enables bdt model"}; + Configurable> ptBinEdges{"ptBinEdges", {0., 30.}, "Pair-pT bin edges of the BDT models, one model per bin (pairs outside are rejected)"}; + Configurable> scoreCuts{"scoreCuts", {std::array{0., 0.}.data(), 1, 2, {"pT bin 0"}, {"Background score", "Signal score"}}, "BDT score cuts, one row per pT bin"}; + Configurable> cutDir{"cutDir", std::vector{o2::cuts_ml::CutNot, o2::cuts_ml::CutNot}, "Cut direction per class: 0 = keep score < cut, 1 = keep score >= cut, 2 = no cut"}; + Configurable> namesInputFeatures{"namesInputFeatures", std::vector{"lambdaDCADau", "lambdaAlpha", "lambdaDCANegPV", "lambdaDCAPosPV", "lambdaQt", "photonAlpha", "photonCosPA", "photonDCADau", "photonDCANegPV", "photonDCAPosPV", "photonQt", "photonRadius", "opAngle"}, "Names and order of the BDT input features (see ResonanceMlResponse.h): must match FeaturesToTrain"}; + } bdt; + ConfigurableAxis axisVertexMixBkg{"axisVertexMixBkg", {VARIABLE_WIDTH, -10.f, -8.f, -6.f, -4.f, -2.f, 0.f, 2.f, 4.f, 6.f, 8.f, 10.f}, "z-vertex bins for mixing"}; ConfigurableAxis axisCentralityMixBkg{"axisCentralityMixBkg", {VARIABLE_WIDTH, 0.0f, 1.0f, 5.0f, 10.0f, 20.0f, 30.0f, 40.0f, 50.0f, 60.0f, 70.0f, 80.0f, 90.0f, 100.0f, 110.0f}, "centrality bins for mixing"}; @@ -250,6 +276,14 @@ struct k892hadronphotonBkg { ConfigurableAxis axisAPQt{"axisAPQt", {220, 0.0f, 1.1f}, "Resonance AP q_{T} (GeV/c)"}; ConfigurableAxis axisCandSel{"axisCandSel", {15, 0.5f, +15.5f}, "Candidate Selection"}; ConfigurableAxis axisOPAngle{"axisOPAngle", {140, 0.0f, 7.0f}, "Opening angle (rad)"}; + // BDT QA axes + ConfigurableAxis mlProb{"mlProb", {100, 0.0f, 1.0f}, "BDT signal score"}; + ConfigurableAxis axisCosPA{"axisCosPA", {200, 0.5f, 1.0f}, "Cosine of pointing angle"}; + ConfigurableAxis axisDCAdau{"axisDCAdau", {50, 0.0f, 5.0f}, "DCA (cm)"}; + ConfigurableAxis axisSignedDCAtoPV{"axisSignedDCAtoPV", {1000, -50.0f, 50.0f}, "signed DCA (cm)"}; + ConfigurableAxis axisSignedDCAtoPVLambda{"axisSignedDCAtoPVLambda", {500, -10.0f, 10.0f}, "signed DCA (cm)"}; + ConfigurableAxis axisV0APQt{"axisV0APQt", {220, 0.0f, 0.5f}, "V0 AP q_{T} (GeV/c)"}; + ConfigurableAxis axisV0Radius{"axisV0Radius", {240, 0.0f, 120.0f}, "V0 radius (cm)"}; } axisConfig; void init(InitContext const&) @@ -259,6 +293,58 @@ struct k892hadronphotonBkg { ccdb->setCaching(true); ccdb->setFatalWhenNull(false); + if (bdt.enableML) { + ccdb->setURL(bdt.ccdbUrl.value); + + // One model per pair-pT bin. MlResponse checks the model files and cutDir, not the rows of scoreCuts + constexpr uint8_t NClassesML = 2; // background, signal + if (bdt.scoreCuts.value.rows() != bdt.ptBinEdges.value.size() - 1 || bdt.scoreCuts.value.cols() != NClassesML) { + LOG(fatal) << "bdt.scoreCuts needs one row per pT bin and " << static_cast(NClassesML) << " columns"; + } + mlResponse.configure(bdt.ptBinEdges.value, bdt.scoreCuts.value, bdt.cutDir.value, NClassesML); + mlResponse.cacheInputFeaturesIndices(bdt.namesInputFeatures); + + if (bdt.loadModelsFromCCDB) { + ccdbApi.init(bdt.ccdbUrl); + LOG(info) << "Fetching models for timestamp: " << bdt.timestampCCDB.value; + mlResponse.setModelPathsCCDB(bdt.onnxFileNames.value, ccdbApi, bdt.modelPathsCCDB.value, bdt.timestampCCDB.value); + } else { + mlResponse.setModelPathsLocal(bdt.onnxFileNames.value); + } + mlResponse.init(bdt.enableOptimizations.value, bdt.numThreads.value); + + // It is applied to the Lambda(1520) MIXED background only (for now!!) + if (!lstarBkgConfig.doSameEvtRotation && !lstarBkgConfig.doEvtMixing) { + LOG(warning) << "bdt.enableML is set but no Lambda(1520) background is requested: the BDT will not be applied."; + } + if (kstarBkgConfig.doSameEvtRotation || kstarBkgConfig.doEvtMixing) { + LOG(info) << "The BDT (gamma + Lambda features) is not applied to the K*(892) background."; + } + + histos.add("BDT/hScoreSignal", "hScoreSignal", kTH1D, {axisConfig.mlProb}); + histos.add("BDT/hScoreBackground", "hScoreBackground", kTH1D, {axisConfig.mlProb}); + histos.add("BDT/h2dScoreVsMassSignal", "h2dScoreVsMassSignal", kTH2D, {axisConfig.axisLambdaStarMass, axisConfig.mlProb}); + histos.add("BDT/h2dScoreVsPtSignal", "h2dScoreVsPtSignal", kTH2D, {axisConfig.axisPt, axisConfig.mlProb}); + histos.add("BDT/h3dScoreSignal", "h3dScoreSignal", kTH3D, {axisConfig.axisPt, axisConfig.axisLambdaStarMass, axisConfig.mlProb}); + histos.add("BDT/h2dScoreVsMassBackground", "h2dScoreVsMassBackground", kTH2D, {axisConfig.axisLambdaStarMass, axisConfig.mlProb}); + histos.add("BDT/h2dScoreVsPtBackground", "h2dScoreVsPtBackground", kTH2D, {axisConfig.axisPt, axisConfig.mlProb}); + histos.add("BDT/h3dScoreBackground", "h3dScoreBackground", kTH3D, {axisConfig.axisPt, axisConfig.axisLambdaStarMass, axisConfig.mlProb}); + histos.add("BDT/h2dLambdaDCADaughters", "h2dLambdaDCADaughters", kTH2D, {axisConfig.mlProb, axisConfig.axisDCAdau}); + + histos.add("BDT/h2dLambdaAlpha", "h2dLambdaAlpha", kTH2D, {axisConfig.mlProb, axisConfig.axisAPAlpha}); + histos.add("BDT/h2dLambdaDCANegPV", "h2dLambdaDCANegPV", kTH2D, {axisConfig.mlProb, axisConfig.axisSignedDCAtoPVLambda}); + histos.add("BDT/h2dLambdaDCAPosPV", "h2dLambdaDCAPosPV", kTH2D, {axisConfig.mlProb, axisConfig.axisSignedDCAtoPVLambda}); + histos.add("BDT/h2dLambdaQt", "h2dLambdaQt", kTH2D, {axisConfig.mlProb, axisConfig.axisV0APQt}); + histos.add("BDT/h2dPhotonAlpha", "h2dPhotonAlpha", kTH2D, {axisConfig.mlProb, axisConfig.axisAPAlpha}); + histos.add("BDT/h2dPhotonCosPA", "h2dPhotonCosPA", kTH2D, {axisConfig.mlProb, axisConfig.axisCosPA}); + histos.add("BDT/h2dPhotonDCADau", "h2dPhotonDCADau", kTH2D, {axisConfig.mlProb, axisConfig.axisDCAdau}); + histos.add("BDT/h2dPhotonDCANegPV", "h2dPhotonDCANegPV", kTH2D, {axisConfig.mlProb, axisConfig.axisSignedDCAtoPV}); + histos.add("BDT/h2dPhotonDCAPosPV", "h2dPhotonDCAPosPV", kTH2D, {axisConfig.mlProb, axisConfig.axisSignedDCAtoPV}); + histos.add("BDT/h2dPhotonQt", "h2dPhotonQt", kTH2D, {axisConfig.mlProb, axisConfig.axisV0APQt}); + histos.add("BDT/h2dPhotonRadius", "h2dPhotonRadius", kTH2D, {axisConfig.mlProb, axisConfig.axisV0Radius}); + histos.add("BDT/h2dOPAngle", "h2dOPAngle", kTH2D, {axisConfig.mlProb, axisConfig.axisOPAngle}); + } + histos.add("hEventCentrality", "hEventCentrality", kTH1D, {axisConfig.axisCentrality}); if (eventSelections.fUseEventSelection) { @@ -351,6 +437,8 @@ struct k892hadronphotonBkg { histos.add("LambdaStarBkg/h4dMixedLambdaStarPtVsAPAlphaVsAPQt", "h4dMixedLambdaStarPtVsAPAlphaVsAPQt", kTHnD, {axisConfig.axisAPAlpha, axisConfig.axisAPQt, axisConfig.axisPt, axisConfig.axisLambdaStarMass}); } } + + histos.print(); } //_______________________________________________ @@ -793,6 +881,73 @@ struct k892hadronphotonBkg { return std::sqrt(RecoDecay::p2(hadronP) - dp * dp / momTot2); } + //_______________________________________________ + // The two V0s share a daughter track (or are the same V0): rejected by the builder + template + static bool shareDaughters(TV0Object const& photon, TV0Object const& hadron) + { + return photon.globalIndex() == hadron.globalIndex() || + photon.posTrackExtraId() == hadron.posTrackExtraId() || + photon.negTrackExtraId() == hadron.negTrackExtraId() || + photon.posTrackExtraId() == hadron.negTrackExtraId() || + photon.negTrackExtraId() == hadron.posTrackExtraId(); + } + + //_______________________________________________ + // Fill BDT performance QA + template + void fillBDTPerformance(TV0Object const& lambda, TV0Object const& photon, float openAngle, float score, float pt, float mass) + { + float bkgScore = 1.0f - score; + + // Signal-probability output + histos.fill(HIST("BDT/hScoreSignal"), score); + histos.fill(HIST("BDT/h2dScoreVsMassSignal"), mass, score); + histos.fill(HIST("BDT/h2dScoreVsPtSignal"), pt, score); + histos.fill(HIST("BDT/h3dScoreSignal"), pt, mass, score); + + // Background-probability output + histos.fill(HIST("BDT/hScoreBackground"), bkgScore); + histos.fill(HIST("BDT/h2dScoreVsMassBackground"), mass, bkgScore); + histos.fill(HIST("BDT/h2dScoreVsPtBackground"), pt, bkgScore); + histos.fill(HIST("BDT/h3dScoreBackground"), pt, mass, bkgScore); + + // Signal score vs the main topological variables + histos.fill(HIST("BDT/h2dLambdaDCADaughters"), score, lambda.dcaV0daughters()); + histos.fill(HIST("BDT/h2dLambdaAlpha"), score, lambda.alpha()); + histos.fill(HIST("BDT/h2dLambdaDCANegPV"), score, lambda.dcanegtopv()); + histos.fill(HIST("BDT/h2dLambdaDCAPosPV"), score, lambda.dcapostopv()); + histos.fill(HIST("BDT/h2dLambdaQt"), score, lambda.qtarm()); + histos.fill(HIST("BDT/h2dPhotonAlpha"), score, photon.alpha()); + histos.fill(HIST("BDT/h2dPhotonCosPA"), score, photon.v0cosPA()); + histos.fill(HIST("BDT/h2dPhotonDCADau"), score, photon.dcaV0daughters()); + histos.fill(HIST("BDT/h2dPhotonDCANegPV"), score, photon.dcanegtopv()); + histos.fill(HIST("BDT/h2dPhotonDCAPosPV"), score, photon.dcapostopv()); + histos.fill(HIST("BDT/h2dPhotonQt"), score, photon.qtarm()); + histos.fill(HIST("BDT/h2dPhotonRadius"), score, photon.v0radius()); + histos.fill(HIST("BDT/h2dOPAngle"), score, openAngle); + } + + //_______________________________________________ + // BDT selection of a Lambda + photon pair + template + bool selectML(TV0Object const& lambda, TV0Object const& photon, + float openAngle, float pt, float mass) + { + // No model outside the bdt.ptBinEdges range + if (pt < bdt.ptBinEdges.value.front() || pt >= bdt.ptBinEdges.value.back()) + return false; + + // Features in the order of bdt.namesInputFeatures + auto inputFeatures = mlResponse.getInputFeatures(lambda, photon, openAngle); + std::vector outputMl; + const bool isSelected = mlResponse.isSelectedMl(inputFeatures, pt, outputMl); // model and cut of the pT bin + + fillBDTPerformance(lambda, photon, openAngle, outputMl[1], pt, mass); + + return isSelected; + } + //_______________________________________________ // Compute same-event rotational background within a single collision. template @@ -820,6 +975,10 @@ struct k892hadronphotonBkg { for (const int& pIdx : photonIndices) { const auto& photon = fullV0s.rawIteratorAt(pIdx); + // Same pair rejection + if (shareDaughters(photon, hadron)) + continue; + // photon as a massless 4-vector ROOT::Math::PtEtaPhiMVector pGamma(photon.pt(), photon.eta(), @@ -853,6 +1012,16 @@ struct k892hadronphotonBkg { // Opening angle between photon and hadron double cosOA = gammaLeg.Vect().Dot(hadronLeg.Vect()) / (gammaLeg.P() * hadronLeg.P()); double openAngle = std::acos(cosOA); + double pt = reso.Pt(); + double mass = reso.M(); + + // // To:Do BDT selection (Lambda(1520)) + // if constexpr (resonance == kResoLambdaStar) { + // if (bdt.enableML) { + // if (!selectML(hadron, photon, openAngle, pt, mass)) + // continue; + // } + // } // Armenteros-Podolanski of the rotated pair const std::array gammaMom{static_cast(gammaLeg.Px()), static_cast(gammaLeg.Py()), static_cast(gammaLeg.Pz())}; @@ -882,9 +1051,9 @@ struct k892hadronphotonBkg { //_______________________________________________ // Mixed-event pairing: hadrons and photons come from two different collisions. - // Centrality is taken from the reference collision (the first of the pair) - template - void calculateMixedBackground(TCollision const& refColl, + // Centrality is taken from the reference collision (the first of the pair). + template + void calculateMixedBackground(TRefColl const& refColl, std::vector const& hadronIndices, std::vector const& photonIndices, TV0s const& fullV0s) @@ -909,6 +1078,11 @@ struct k892hadronphotonBkg { for (const int& pIdx : photonIndices) { const auto& photon = fullV0s.rawIteratorAt(pIdx); + + // Same pair rejection as the builder + if (shareDaughters(photon, hadron)) + continue; + float pP = std::hypot(photon.px(), photon.py(), photon.pz()); ROOT::Math::PxPyPzEVector fourMomPhoton( photon.px(), photon.py(), photon.pz(), pP); @@ -931,6 +1105,14 @@ struct k892hadronphotonBkg { if (std::abs(rapidity) > maxRap) continue; + // BDT selection (Lambda(1520) only) + if constexpr (resonance == kResoLambdaStar) { + if (bdt.enableML) { + if (!selectML(hadron, photon, openAngle, pt, mass)) + continue; + } + } + // Armenteros-Podolanski of the mixed pair const std::array gammaMom{photon.px(), photon.py(), photon.pz()}; const std::array hadronMom{hadron.px(), hadron.py(), hadron.pz()}; diff --git a/PWGLF/Tasks/Strangeness/CMakeLists.txt b/PWGLF/Tasks/Strangeness/CMakeLists.txt index 7a757d05e88..5d3711f586a 100644 --- a/PWGLF/Tasks/Strangeness/CMakeLists.txt +++ b/PWGLF/Tasks/Strangeness/CMakeLists.txt @@ -86,7 +86,7 @@ o2physics_add_dpl_workflow(hstrangecorrelation o2physics_add_dpl_workflow(sigmaanalysis SOURCES sigmaanalysis.cxx - PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2Physics::AnalysisCCDB + PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2Physics::AnalysisCCDB O2Physics::MLCore COMPONENT_NAME Analysis) o2physics_add_dpl_workflow(phik0shortanalysis diff --git a/PWGLF/Tasks/Strangeness/sigmaanalysis.cxx b/PWGLF/Tasks/Strangeness/sigmaanalysis.cxx index fae9617f013..645bf6fcc3f 100644 --- a/PWGLF/Tasks/Strangeness/sigmaanalysis.cxx +++ b/PWGLF/Tasks/Strangeness/sigmaanalysis.cxx @@ -21,17 +21,21 @@ #include "PWGLF/DataModel/LFSigmaTables.h" #include "PWGLF/DataModel/LFStrangenessTables.h" +#include "PWGLF/Utils/ResonanceMlResponse.h" #include "Common/CCDB/EventSelectionParams.h" #include "Common/CCDB/ctpRateFetcher.h" +#include "Tools/ML/MlResponse.h" #include +// #include #include #include #include #include #include #include +#include #include #include #include @@ -47,6 +51,7 @@ #include #include +#include #include #include #include @@ -54,6 +59,7 @@ #include using namespace o2; +using namespace o2::ml; using namespace o2::framework; using namespace o2::framework::expressions; @@ -84,7 +90,9 @@ enum CentEstimator { struct sigmaanalysis { Service ccdb; + o2::ccdb::CcdbApi ccdbApi; ctpRateFetcher rateFetcher; + o2::analysis::ResonanceMlResponse mlResponse; //__________________________________________________ HistogramRegistry histos{"Histos", {}, OutputObjHandlingPolicy::AnalysisObject}; @@ -100,6 +108,24 @@ struct sigmaanalysis { Configurable fIRCrashOnNull{"fIRCrashOnNull", false, "Flag to avoid CTP RateFetcher crash."}; Configurable irSource{"irSource", "T0VTX", "Estimator of the interaction rate (Recommended: pp --> T0VTX, Pb-Pb --> ZNC hadronic)"}; + struct : ConfigurableGroup { + std::string prefix = "bdt"; // JSON group name + Configurable ccdbUrl{"ccdbUrl", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; + Configurable> onnxFileNames{"onnxFileNames", std::vector{"BDTModel.onnx"}, "Local .onnx file names, one per pT bin"}; + Configurable> modelPathsCCDB{"modelPathsCCDB", std::vector{"Users/o/obenchik/MLModels/BDT"}, "Model paths on CCDB, one per pT bin (each model needs its own folder)"}; + Configurable timestampCCDB{"timestampCCDB", 1695750420200, "timestamp of the ONNX file for ML model used to query in CCDB. Please use 1695750420200"}; + 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)"}; + Configurable numThreads{"numThreads", 1, "ONNX intra-op threads. 0 lets ONNX Runtime default to one thread per physical core"}; + Configurable enableML{"enableML", false, "Enables bdt model"}; + Configurable> ptBinEdges{"ptBinEdges", {0., 30.}, "Candidate-pT bin edges of the BDT models, one model per bin (candidates outside are rejected)"}; + Configurable> scoreCuts{"scoreCuts", {std::array{0., 0.}.data(), 1, 2, {"pT bin 0"}, {"Background score", "Signal score"}}, "BDT score cuts, one row per pT bin"}; + Configurable> cutDir{"cutDir", std::vector{o2::cuts_ml::CutNot, o2::cuts_ml::CutNot}, "Cut direction per class: 0 = keep score < cut, 1 = keep score >= cut, 2 = no cut"}; + // Configurable mlCorrThreshold{"mlCorrThreshold", 0.5, "Threshold for correlated K* score"}; + Configurable> namesInputFeatures{"namesInputFeatures", std::vector{"lambdaDCADau", "lambdaAlpha", "lambdaDCANegPV", "lambdaDCAPosPV", "lambdaQt", "photonAlpha", "photonCosPA", "photonDCADau", "photonDCANegPV", "photonDCAPosPV", "photonQt", "photonRadius", "opAngle"}, "Names and order of the BDT input features (see ResonanceMlResponse.h): must match FeaturesToTrain"}; + + } bdt; + struct : ConfigurableGroup { std::string prefix = "eventSelections"; // JSON group name Configurable requireSel8{"requireSel8", true, "require sel8 event selection"}; @@ -287,6 +313,8 @@ struct sigmaanalysis { ConfigurableAxis axisV0Radius{"axisV0Radius", {240, 0.0f, 120.0f}, "V0 radius (cm)"}; ConfigurableAxis axisV0PairRadius{"axisV0PairRadius", {200, 0.0f, 20.0f}, "V0Pair radius (cm)"}; ConfigurableAxis axisDCAtoPV{"axisDCAtoPV", {500, 0.0f, 50.0f}, "DCA (cm)"}; + ConfigurableAxis axisSignedDCAtoPV{"axisSignedDCAtoPV", {1000, -50.0f, 50.0f}, "signed DCA (cm)"}; + ConfigurableAxis axisSignedDCAtoPVLambda{"axisSignedDCAtoPVLambda", {500, -10.0f, 10.0f}, "signed DCA (cm)"}; ConfigurableAxis axisDCAdau{"axisDCAdau", {50, 0.0f, 5.0f}, "DCA (cm)"}; ConfigurableAxis axisCosPA{"axisCosPA", {200, 0.5f, 1.0f}, "Cosine of pointing angle"}; ConfigurableAxis axisPA{"axisPA", {100, 0.0f, 1}, "Pointing angle"}; @@ -305,7 +333,7 @@ struct sigmaanalysis { ConfigurableAxis axisCandSel{"axisCandSel", {20, 0.5f, +20.5f}, "Candidate Selection"}; // ML - ConfigurableAxis MLProb{"MLOutput", {100, 0.0f, 1.0f}, ""}; + ConfigurableAxis mlProb{"MLOutput", {100, 0.0f, 1.0f}, ""}; int NSigma0Cand = 0; void init(InitContext const&) @@ -321,6 +349,65 @@ struct sigmaanalysis { ccdb->setCaching(true); ccdb->setFatalWhenNull(false); + if (bdt.enableML) { + ccdb->setURL(bdt.ccdbUrl.value); + + // One model per candidate-pT bin + constexpr uint8_t NClassesML = 2; // background, signal + if (bdt.scoreCuts.value.rows() != bdt.ptBinEdges.value.size() - 1 || bdt.scoreCuts.value.cols() != NClassesML) { + LOG(fatal) << "bdt.scoreCuts needs one row per pT bin and " << static_cast(NClassesML) << " columns"; + } + mlResponse.configure(bdt.ptBinEdges.value, bdt.scoreCuts.value, bdt.cutDir.value, NClassesML); + mlResponse.cacheInputFeaturesIndices(bdt.namesInputFeatures); + + if (bdt.loadModelsFromCCDB) { + ccdbApi.init(bdt.ccdbUrl); + LOG(info) << "Fetching models for timestamp: " << bdt.timestampCCDB.value; + mlResponse.setModelPathsCCDB(bdt.onnxFileNames.value, ccdbApi, bdt.modelPathsCCDB.value, bdt.timestampCCDB.value); + } else { + mlResponse.setModelPathsLocal(bdt.onnxFileNames.value); + } + mlResponse.init(bdt.enableOptimizations.value, bdt.numThreads.value); + + // The model is trained on PCM (V0) photon features, which EMCal clusters do not have + if (doprocessRealDataWithEMCal || doprocessMonteCarloWithEMCal) { + LOG(fatal) << "BDT selection is only available for PCM photons: disable bdt.enableML or use processRealData/processMonteCarlo."; + } + + // BDT performance QA + histos.add("BDT/hScoreSignal", "hScoreSignal", kTH1D, {mlProb}); + histos.add("BDT/hScoreBackground", "hScoreBackground", kTH1D, {mlProb}); + histos.add("BDT/h2dScoreVsMassSignal", "h2dScoreVsMassSignal", kTH2D, {axisSigmaMass, mlProb}); + histos.add("BDT/h2dScoreVsPtSignal", "h2dScoreVsPtSignal", kTH2D, {axisPt, mlProb}); + histos.add("BDT/h3dScoreSignal", "h3dScoreSignal", kTH3D, {axisPt, axisSigmaMass, mlProb}); + histos.add("BDT/h2dScoreVsMassBackground", "h2dScoreVsMassBackground", kTH2D, {axisSigmaMass, mlProb}); + histos.add("BDT/h2dScoreVsPtBackground", "h2dScoreVsPtBackground", kTH2D, {axisPt, mlProb}); + histos.add("BDT/h3dScoreBackground", "h3dScoreBackground", kTH3D, {axisPt, axisSigmaMass, mlProb}); + + // Signal score vs the main topological variables + histos.add("BDT/h2dDCADaughters", "h2dDCADaughters", kTH2D, {mlProb, axisDCAdau}); + histos.add("BDT/h2dLambdaAlpha", "h2dLambdaAlpha", kTH2D, {mlProb, axisAPAlpha}); + histos.add("BDT/h2dLambdaDCANegPV", "h2dLambdaDCANegPV", kTH2D, {mlProb, axisSignedDCAtoPVLambda}); + histos.add("BDT/h2dLambdaDCAPosPV", "h2dLambdaDCAPosPV", kTH2D, {mlProb, axisSignedDCAtoPVLambda}); + histos.add("BDT/h2dLambdaQt", "h2dLambdaQt", kTH2D, {mlProb, axisAPQt}); + histos.add("BDT/h2dPhotonAlpha", "h2dPhotonAlpha", kTH2D, {mlProb, axisAPAlpha}); + histos.add("BDT/h2dPhotonCosPA", "h2dPhotonCosPA", kTH2D, {mlProb, axisCosPA}); + histos.add("BDT/h2dPhotonDCADau", "h2dPhotonDCADau", kTH2D, {mlProb, axisDCAdau}); + histos.add("BDT/h2dPhotonDCANegPV", "h2dPhotonDCANegPV", kTH2D, {mlProb, axisSignedDCAtoPV}); + histos.add("BDT/h2dPhotonDCAPosPV", "h2dPhotonDCAPosPV", kTH2D, {mlProb, axisSignedDCAtoPV}); + histos.add("BDT/h2dPhotonQt", "h2dPhotonQt", kTH2D, {mlProb, axisAPQt}); + histos.add("BDT/h2dPhotonRadius", "h2dPhotonRadius", kTH2D, {mlProb, axisV0Radius}); + histos.add("BDT/h2dOPAngle", "h2dOPAngle", kTH2D, {mlProb, axisOPAngle}); + + // MC-truth-based score + if (doprocessMonteCarlo) { + histos.add("BDT/hScoreTrueSignal", "hScoreTrueSignal", kTH1D, {mlProb}); + histos.add("BDT/hScoreTrueBackground", "hScoreTrueBackground", kTH1D, {mlProb}); + histos.add("BDT/h2dScoreVsPtTrueSignal", "h2dScoreVsPtTrueSignal", kTH2D, {axisPt, mlProb}); + histos.add("BDT/h2dScoreVsPtTrueBackground", "h2dScoreVsPtTrueBackground", kTH2D, {axisPt, mlProb}); + } + } + // Event Counters histos.add("hEventCentrality", "hEventCentrality", kTH1D, {axisCentrality}); histos.add("hCentralityVsNch", "hCentralityVsNch", kTH2D, {{101, 0.0f, 101.0f}, axisNch}); @@ -1652,10 +1739,82 @@ struct sigmaanalysis { return true; } + // Fill BDT performance QA + template + void fillBDTPerformance(TSigma0Object const& cand, float score) + { + float pt = cand.pt(); + float mass = cand.sigma0Mass(); + float bkgScore = 1.0f - score; + + // Signal-probability output + histos.fill(HIST("BDT/hScoreSignal"), score); + histos.fill(HIST("BDT/h2dScoreVsMassSignal"), mass, score); + histos.fill(HIST("BDT/h2dScoreVsPtSignal"), pt, score); + histos.fill(HIST("BDT/h3dScoreSignal"), pt, mass, score); + + // Background-probability output + histos.fill(HIST("BDT/hScoreBackground"), bkgScore); + histos.fill(HIST("BDT/h2dScoreVsMassBackground"), mass, bkgScore); + histos.fill(HIST("BDT/h2dScoreVsPtBackground"), pt, bkgScore); + histos.fill(HIST("BDT/h3dScoreBackground"), pt, mass, bkgScore); + + // Signal score vs the main topological variables + histos.fill(HIST("BDT/h2dDCADaughters"), score, cand.lambdaDCADau()); + histos.fill(HIST("BDT/h2dLambdaAlpha"), score, cand.lambdaAlpha()); + histos.fill(HIST("BDT/h2dLambdaDCANegPV"), score, cand.lambdaDCANegPV()); + histos.fill(HIST("BDT/h2dLambdaDCAPosPV"), score, cand.lambdaDCAPosPV()); + histos.fill(HIST("BDT/h2dLambdaQt"), score, cand.lambdaQt()); + histos.fill(HIST("BDT/h2dPhotonAlpha"), score, cand.photonAlpha()); + histos.fill(HIST("BDT/h2dPhotonCosPA"), score, cand.photonCosPA()); + histos.fill(HIST("BDT/h2dPhotonDCADau"), score, cand.photonDCADau()); + histos.fill(HIST("BDT/h2dPhotonDCANegPV"), score, cand.photonDCANegPV()); + histos.fill(HIST("BDT/h2dPhotonDCAPosPV"), score, cand.photonDCAPosPV()); + histos.fill(HIST("BDT/h2dPhotonQt"), score, cand.photonQt()); + histos.fill(HIST("BDT/h2dPhotonRadius"), score, cand.photonRadius()); + histos.fill(HIST("BDT/h2dOPAngle"), score, cand.opAngle()); + + // MC-truth-based separation (signal = particle + antiparticle) + if constexpr (requires { cand.isSigma0(); cand.isLambdaStar(); }) { + bool isTrueSignal = doLambdaStar ? (cand.isLambdaStar() || cand.isAntiLambdaStar()) : (cand.isSigma0() || cand.isAntiSigma0()); + if (isTrueSignal) { + histos.fill(HIST("BDT/hScoreTrueSignal"), score); + histos.fill(HIST("BDT/h2dScoreVsPtTrueSignal"), pt, score); + } else { + histos.fill(HIST("BDT/hScoreTrueBackground"), score); + histos.fill(HIST("BDT/h2dScoreVsPtTrueBackground"), pt, score); + } + } + } + + template + bool selectML(TSigma0Object const& cand) + { + // No model outside the bdt.ptBinEdges range + const float pt = cand.pt(); + if (pt < bdt.ptBinEdges.value.front() || pt >= bdt.ptBinEdges.value.back()) + return false; + + // Features in the order of bdt.namesInputFeatures + auto inputFeatures = mlResponse.getInputFeatures(cand, cand, cand.opAngle()); + std::vector outputMl; // [background, signal] + const bool isSelected = mlResponse.isSelectedMl(inputFeatures, pt, outputMl); // model and cut of the pT bin + + fillBDTPerformance(cand, outputMl[1]); + + return isSelected; + } + // Apply selections in sigma0 candidates template bool processSigma0Candidate(TSigma0Object const& cand) { + // BDT selection: model trained on PCM photon features only + if constexpr (requires { cand.photonV0Type(); }) { + if (bdt.enableML && !selectML(cand)) + return false; + } + // Photon specific selections if constexpr (requires { cand.photonV0Type(); }) { // Processing PCM photon if (!selectPhoton(cand)) diff --git a/PWGLF/Utils/ResonanceMlResponse.h b/PWGLF/Utils/ResonanceMlResponse.h new file mode 100644 index 00000000000..8a7932bc42f --- /dev/null +++ b/PWGLF/Utils/ResonanceMlResponse.h @@ -0,0 +1,316 @@ +// 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 ResonanceMlResponse.h +/// \brief Class to compute the ML response for Lambda + photon resonance selections (Sigma0, Lambda(1520)), +/// either from sigma0builder candidates or from a Lambda and a photon V0 of the strangeness derived data +/// \author Oussama Benchikhi + +#ifndef PWGLF_UTILS_RESONANCEMLRESPONSE_H_ +#define PWGLF_UTILS_RESONANCEMLRESPONSE_H_ + +#include "PWGLF/DataModel/LFStrangenessPIDTables.h" +#include "PWGLF/DataModel/LFStrangenessTables.h" + +#include "Tools/ML/MlResponse.h" + +#include + +#include +#include + +namespace o2::analysis +{ +// list of input features that can be requested via the namesInputFeatures configurable (sigma0builder getter names) +enum class InputFeaturesResonance : uint8_t { + // Lambda + lambdaQt = 0, + lambdaAlpha, + lambdaRadius, + lambdaCosPA, + lambdaDCADau, + lambdaDCANegPV, + lambdaDCAPosPV, + lambdaPosEta, + lambdaNegEta, + lambdaPosPrTPCNSigma, + lambdaPosPiTPCNSigma, + lambdaNegPrTPCNSigma, + lambdaNegPiTPCNSigma, + // Photon + photonQt, + photonAlpha, + photonCosPA, + photonDCADau, + photonDCANegPV, + photonDCAPosPV, + photonRadius, + photonZconv, + photonPsiPair, + photonPosEta, + photonNegEta, + photonPosTPCNSigmaEl, + photonNegTPCNSigmaEl, + // Photon-Lambda pair + opAngle +}; + +template +class ResonanceMlResponse : public MlResponse +{ + public: + ResonanceMlResponse() = default; + ~ResonanceMlResponse() override = default; + + /// Input features of a photon-Lambda pair + template , typename TLambda, typename TPhoton> + std::vector getInputFeatures(TLambda const& lambda, TPhoton const& photon, float opAngle) + { + // sigma0builder candidates carry prefixed columns (lambdaQt(), photonQt(), ...), derived V0s the plain V0 getters + constexpr bool LambdaFromSigma0 = requires(TLambda const& cand) { cand.lambdaQt(); }; + constexpr bool PhotonFromSigma0 = requires(TPhoton const& cand) { cand.photonQt(); }; + + std::vector inputFeatures; + inputFeatures.reserve(MlResponse::mCachedIndices.size()); + + for (const auto& idx : MlResponse::mCachedIndices) { + switch (static_cast(idx)) { + // Lambda + case InputFeaturesResonance::lambdaQt: + if constexpr (LambdaFromSigma0) { + inputFeatures.emplace_back(lambda.lambdaQt()); + } else { + inputFeatures.emplace_back(lambda.qtarm()); + } + break; + case InputFeaturesResonance::lambdaAlpha: + if constexpr (LambdaFromSigma0) { + inputFeatures.emplace_back(lambda.lambdaAlpha()); + } else { + inputFeatures.emplace_back(lambda.alpha()); + } + break; + case InputFeaturesResonance::lambdaRadius: + if constexpr (LambdaFromSigma0) { + inputFeatures.emplace_back(lambda.lambdaRadius()); + } else { + inputFeatures.emplace_back(lambda.v0radius()); + } + break; + case InputFeaturesResonance::lambdaCosPA: + if constexpr (LambdaFromSigma0) { + inputFeatures.emplace_back(lambda.lambdaCosPA()); + } else { + inputFeatures.emplace_back(lambda.v0cosPA()); + } + break; + case InputFeaturesResonance::lambdaDCADau: + if constexpr (LambdaFromSigma0) { + inputFeatures.emplace_back(lambda.lambdaDCADau()); + } else { + inputFeatures.emplace_back(lambda.dcaV0daughters()); + } + break; + case InputFeaturesResonance::lambdaDCANegPV: + if constexpr (LambdaFromSigma0) { + inputFeatures.emplace_back(lambda.lambdaDCANegPV()); + } else { + inputFeatures.emplace_back(lambda.dcanegtopv()); + } + break; + case InputFeaturesResonance::lambdaDCAPosPV: + if constexpr (LambdaFromSigma0) { + inputFeatures.emplace_back(lambda.lambdaDCAPosPV()); + } else { + inputFeatures.emplace_back(lambda.dcapostopv()); + } + break; + case InputFeaturesResonance::lambdaPosEta: + if constexpr (LambdaFromSigma0) { + inputFeatures.emplace_back(lambda.lambdaPosEta()); + } else { + inputFeatures.emplace_back(lambda.positiveeta()); + } + break; + case InputFeaturesResonance::lambdaNegEta: + if constexpr (LambdaFromSigma0) { + inputFeatures.emplace_back(lambda.lambdaNegEta()); + } else { + inputFeatures.emplace_back(lambda.negativeeta()); + } + break; + case InputFeaturesResonance::lambdaPosPrTPCNSigma: + if constexpr (LambdaFromSigma0) { + inputFeatures.emplace_back(lambda.lambdaPosPrTPCNSigma()); + } else { + inputFeatures.emplace_back(lambda.template posTrackExtra_as().tpcNSigmaPr()); + } + break; + case InputFeaturesResonance::lambdaPosPiTPCNSigma: + if constexpr (LambdaFromSigma0) { + inputFeatures.emplace_back(lambda.lambdaPosPiTPCNSigma()); + } else { + inputFeatures.emplace_back(lambda.template posTrackExtra_as().tpcNSigmaPi()); + } + break; + case InputFeaturesResonance::lambdaNegPrTPCNSigma: + if constexpr (LambdaFromSigma0) { + inputFeatures.emplace_back(lambda.lambdaNegPrTPCNSigma()); + } else { + inputFeatures.emplace_back(lambda.template negTrackExtra_as().tpcNSigmaPr()); + } + break; + case InputFeaturesResonance::lambdaNegPiTPCNSigma: + if constexpr (LambdaFromSigma0) { + inputFeatures.emplace_back(lambda.lambdaNegPiTPCNSigma()); + } else { + inputFeatures.emplace_back(lambda.template negTrackExtra_as().tpcNSigmaPi()); + } + break; + // Photon + case InputFeaturesResonance::photonQt: + if constexpr (PhotonFromSigma0) { + inputFeatures.emplace_back(photon.photonQt()); + } else { + inputFeatures.emplace_back(photon.qtarm()); + } + break; + case InputFeaturesResonance::photonAlpha: + if constexpr (PhotonFromSigma0) { + inputFeatures.emplace_back(photon.photonAlpha()); + } else { + inputFeatures.emplace_back(photon.alpha()); + } + break; + case InputFeaturesResonance::photonCosPA: + if constexpr (PhotonFromSigma0) { + inputFeatures.emplace_back(photon.photonCosPA()); + } else { + inputFeatures.emplace_back(photon.v0cosPA()); + } + break; + case InputFeaturesResonance::photonDCADau: + if constexpr (PhotonFromSigma0) { + inputFeatures.emplace_back(photon.photonDCADau()); + } else { + inputFeatures.emplace_back(photon.dcaV0daughters()); + } + break; + case InputFeaturesResonance::photonDCANegPV: + if constexpr (PhotonFromSigma0) { + inputFeatures.emplace_back(photon.photonDCANegPV()); + } else { + inputFeatures.emplace_back(photon.dcanegtopv()); + } + break; + case InputFeaturesResonance::photonDCAPosPV: + if constexpr (PhotonFromSigma0) { + inputFeatures.emplace_back(photon.photonDCAPosPV()); + } else { + inputFeatures.emplace_back(photon.dcapostopv()); + } + break; + case InputFeaturesResonance::photonRadius: + if constexpr (PhotonFromSigma0) { + inputFeatures.emplace_back(photon.photonRadius()); + } else { + inputFeatures.emplace_back(photon.v0radius()); + } + break; + case InputFeaturesResonance::photonZconv: + if constexpr (PhotonFromSigma0) { + inputFeatures.emplace_back(photon.photonZconv()); + } else { + inputFeatures.emplace_back(photon.z()); + } + break; + case InputFeaturesResonance::photonPsiPair: + if constexpr (PhotonFromSigma0) { + inputFeatures.emplace_back(photon.photonPsiPair()); + } else { + inputFeatures.emplace_back(photon.psipair()); + } + break; + case InputFeaturesResonance::photonPosEta: + if constexpr (PhotonFromSigma0) { + inputFeatures.emplace_back(photon.photonPosEta()); + } else { + inputFeatures.emplace_back(photon.positiveeta()); + } + break; + case InputFeaturesResonance::photonNegEta: + if constexpr (PhotonFromSigma0) { + inputFeatures.emplace_back(photon.photonNegEta()); + } else { + inputFeatures.emplace_back(photon.negativeeta()); + } + break; + case InputFeaturesResonance::photonPosTPCNSigmaEl: + if constexpr (PhotonFromSigma0) { + inputFeatures.emplace_back(photon.photonPosTPCNSigmaEl()); + } else { + inputFeatures.emplace_back(photon.template posTrackExtra_as().tpcNSigmaEl()); + } + break; + case InputFeaturesResonance::photonNegTPCNSigmaEl: + if constexpr (PhotonFromSigma0) { + inputFeatures.emplace_back(photon.photonNegTPCNSigmaEl()); + } else { + inputFeatures.emplace_back(photon.template negTrackExtra_as().tpcNSigmaEl()); + } + break; + // Photon-Lambda pair + case InputFeaturesResonance::opAngle: + inputFeatures.emplace_back(opAngle); + break; + } + } + return inputFeatures; + } + + protected: + /// Method to fill the map of available input features + void setAvailableInputFeatures() override + { + MlResponse::mAvailableInputFeatures = { + {"lambdaQt", static_cast(InputFeaturesResonance::lambdaQt)}, + {"lambdaAlpha", static_cast(InputFeaturesResonance::lambdaAlpha)}, + {"lambdaRadius", static_cast(InputFeaturesResonance::lambdaRadius)}, + {"lambdaCosPA", static_cast(InputFeaturesResonance::lambdaCosPA)}, + {"lambdaDCADau", static_cast(InputFeaturesResonance::lambdaDCADau)}, + {"lambdaDCANegPV", static_cast(InputFeaturesResonance::lambdaDCANegPV)}, + {"lambdaDCAPosPV", static_cast(InputFeaturesResonance::lambdaDCAPosPV)}, + {"lambdaPosEta", static_cast(InputFeaturesResonance::lambdaPosEta)}, + {"lambdaNegEta", static_cast(InputFeaturesResonance::lambdaNegEta)}, + {"lambdaPosPrTPCNSigma", static_cast(InputFeaturesResonance::lambdaPosPrTPCNSigma)}, + {"lambdaPosPiTPCNSigma", static_cast(InputFeaturesResonance::lambdaPosPiTPCNSigma)}, + {"lambdaNegPrTPCNSigma", static_cast(InputFeaturesResonance::lambdaNegPrTPCNSigma)}, + {"lambdaNegPiTPCNSigma", static_cast(InputFeaturesResonance::lambdaNegPiTPCNSigma)}, + {"photonQt", static_cast(InputFeaturesResonance::photonQt)}, + {"photonAlpha", static_cast(InputFeaturesResonance::photonAlpha)}, + {"photonCosPA", static_cast(InputFeaturesResonance::photonCosPA)}, + {"photonDCADau", static_cast(InputFeaturesResonance::photonDCADau)}, + {"photonDCANegPV", static_cast(InputFeaturesResonance::photonDCANegPV)}, + {"photonDCAPosPV", static_cast(InputFeaturesResonance::photonDCAPosPV)}, + {"photonRadius", static_cast(InputFeaturesResonance::photonRadius)}, + {"photonZconv", static_cast(InputFeaturesResonance::photonZconv)}, + {"photonPsiPair", static_cast(InputFeaturesResonance::photonPsiPair)}, + {"photonPosEta", static_cast(InputFeaturesResonance::photonPosEta)}, + {"photonNegEta", static_cast(InputFeaturesResonance::photonNegEta)}, + {"photonPosTPCNSigmaEl", static_cast(InputFeaturesResonance::photonPosTPCNSigmaEl)}, + {"photonNegTPCNSigmaEl", static_cast(InputFeaturesResonance::photonNegTPCNSigmaEl)}, + {"opAngle", static_cast(InputFeaturesResonance::opAngle)}}; + } +}; + +} // namespace o2::analysis + +#endif // PWGLF_UTILS_RESONANCEMLRESPONSE_H_