From 816d9f4670f4df55bc99c69bf1ecf7a44b2c8af7 Mon Sep 17 00:00:00 2001 From: Chiara De Martin Date: Thu, 7 Mar 2024 17:54:16 +0100 Subject: [PATCH 01/11] retrieve ML model and configure it + compute v2 + process data --- PWGLF/TableProducer/cascadeflow.cxx | 154 ++++++++++++++++++++++++++-- 1 file changed, 145 insertions(+), 9 deletions(-) diff --git a/PWGLF/TableProducer/cascadeflow.cxx b/PWGLF/TableProducer/cascadeflow.cxx index c62fa6a5972..bdb9177bfe9 100644 --- a/PWGLF/TableProducer/cascadeflow.cxx +++ b/PWGLF/TableProducer/cascadeflow.cxx @@ -12,19 +12,22 @@ /// \brief Task to create derived data for cascade flow analyses /// \authors: Chiara De Martin (chiara.de.martin@cern.ch), Maximiliano Puccio (maximiliano.puccio@cern.ch) +#include "Math/Vector3D.h" +#include "TRandom3.h" #include "Common/DataModel/Centrality.h" #include "Common/DataModel/EventSelection.h" #include "Common/DataModel/Multiplicity.h" #include "Common/DataModel/PIDResponse.h" #include "Common/DataModel/TrackSelectionTables.h" #include "Framework/AnalysisTask.h" +#include "Framework/ASoAHelpers.h" #include "Framework/O2DatabasePDGPlugin.h" #include "Framework/runDataProcessing.h" #include "PWGLF/DataModel/cascqaanalysis.h" #include "PWGLF/DataModel/LFStrangenessTables.h" #include "PWGLF/DataModel/LFStrangenessPIDTables.h" -#include "TRandom3.h" -#include "Framework/ASoAHelpers.h" +#include "Tools/ML/MlResponse.h" +#include "PWGHF/Core/HfHelper.h" using namespace o2; using namespace o2::framework; @@ -32,6 +35,9 @@ using namespace o2::framework::expressions; using std::array; using DauTracks = soa::Join; +using CollEventPlane = soa::Join::iterator; + +ROOT::Math::XYZVector eventplaneVecT0A, eventplaneVecT0C, cascphiVec; namespace cascadev2 { @@ -46,6 +52,8 @@ static const std::vector massSigmaParameterNames{"p0", "p1", "p2", static const std::vector speciesNames{"Xi", "Omega"}; } // namespace cascadev2 +static constexpr double defaultCutsMl[1][2] = {{0.5, 0.5}}; + struct cascadeFlow { Configurable sideBandStart{"sideBandStart", 5, "Start of the sideband region in number of sigmas"}; @@ -55,6 +63,21 @@ struct cascadeFlow { Configurable nsigmatpcPr{"nsigmatpcPr", 5, "nsigmatpcPr"}; Configurable nsigmatpcPi{"nsigmatpcPi", 5, "nsigmatpcPi"}; Configurable mintpccrrows{"mintpccrrows", 3, "mintpccrrows"}; + Configurable loadModelsFromCCDB{"loadModelsFromCCDB", false, "Flag to enable or disable the loading of models from CCDB"}; + Configurable ccdbUrl{"ccdbUrl", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; + Configurable modelPathsCCDB{"modelPathsCCDB", "Users/c/chdemart/CascadesFlow", "Path on CCDB"}; + Configurable timestampCCDB{"timestampCCDB", -1, "timestamp of the ONNX file for ML model used to query in CCDB"}; + Configurable> binsPtMl{"binsPtMl", std::vector{0.6, 10.}, "pT bin limits for ML application"}; + Configurable> cutDirMl{"cutDirMl", std::vector{cuts_ml::CutSmaller, cuts_ml::CutNot}, "Whether to reject score values greater or smaller than the threshold"}; + Configurable> cutsMl{"cutsMl", {defaultCutsMl[0], 1, 3, {"pT bin 0"}, {"score signal", "score bkg"}}, "ML selections per pT bin"}; + Configurable nClassesMl{"nClassesMl", (int8_t)2, "Number of classes in ML model"}; + + HfHelper hfHelper; + o2::ccdb::CcdbApi ccdbApi; + + // Add objects needed for ML inference + o2::analysis::MlResponse mlResponse; + std::vector outputMl = {}; template bool IsCascAccepted(TCascade casc, TDaughter negExtra, TDaughter posExtra, TDaughter bachExtra, int& counter) // loose cuts on topological selections of cascades @@ -78,10 +101,23 @@ struct cascadeFlow { return true; } + double GetPhiInRange(double phi) + { + double result = phi; + while (result < 0) { + result = result + 2. * TMath::Pi() / 2; + } + while (result > 2. * TMath::Pi() / 2) { + result = result - 2. * TMath::Pi() / 2; + } + return result; + } + HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject}; // Tables to produce Produces trainingSample; + Produces analysisSample; Configurable> parSigmaMass{ "parSigmaMass", {cascadev2::massSigmaParameters[0], 4, 2, @@ -119,13 +155,50 @@ struct cascadeFlow { pdgCode); } + template + // void fillAnalisedTable(collision_t coll, cascade_t casc, float BDTresponse) + void fillAnalysedTable(collision_t coll, cascade_t casc) + { + analysisSample(coll.centFT0M(), + casc.sign(), + casc.pt(), + casc.eta(), + casc.mXi(), + casc.mOmega(), + casc.mLambda()); + } + void init(InitContext const&) { + ConfigurableAxis vertexZ{"vertexZ", {20, -10, 10}, "vertex axis for histograms"}; histos.add("hEventVertexZ", "hEventVertexZ", kTH1F, {vertexZ}); histos.add("hEventCentrality", "hEventCentrality", kTH1F, {{101, 0, 101}}); histos.add("hCandidate", "hCandidate", HistType::kTH1D, {{22, -0.5, 21.5}}); histos.add("hCascadeSignal", "hCascadeSignal", HistType::kTH1D, {{6, -0.5, 5.5}}); + histos.add("hCascade", "hCascade", HistType::kTH1D, {{6, -0.5, 5.5}}); + histos.add("hCascadePhi", "hCascadePhi", HistType::kTH1D, {{100, 0, 2*TMath::Pi()}}); + histos.add("hcascminuspsiT0A", "hcascminuspsiT0A", HistType::kTH1D, {{100, 0, 2*TMath::Pi()}}); + histos.add("hcascminuspsiT0C", "hcascminuspsiT0C", HistType::kTH1D, {{100, 0, 2*TMath::Pi()}}); + histos.add("hPsiT0A", "hPsiT0A", HistType::kTH1D, {{100, 0, 2*TMath::Pi()}}); + histos.add("hPsiT0C", "hPsiT0C", HistType::kTH1D, {{100, 0, 2*TMath::Pi()}}); + histos.add("hFT0ARe", "hFT0ARe", HistType::kTH1D, {{100, -1, 1}}); + histos.add("hFT0AIm", "hFT0AIm", HistType::kTH1D, {{100, -1, 1}}); + histos.add("hFT0A", "hFT0A", HistType::kTH2D, {{100, -1, 1}, {100, -1, 1}}); + histos.add("hv2Norm", "hv2Norm", HistType::kTH1D, {{100, 0, 1}}); + + // Configure and initialise the ML class + mlResponse.configure(binsPtMl, cutsMl, cutDirMl, nClassesMl); + + // Bonus: retrieve the model from CCDB (needed for ML application on the GRID) + if (loadModelsFromCCDB) { + ccdbApi.init(ccdbUrl); + mlResponse.setModelPathsCCDB(onnxFileNames, ccdbApi, modelPathsCCDB.value, timestampCCDB); + } else { + mlResponse.setModelPathsLocal(onnxFileNames); + } + + mlResponse.init(); } void processTrainingBackground(soa::Join::iterator const& coll, soa::Join const& Cascades, DauTracks const&) @@ -153,9 +226,9 @@ struct cascadeFlow { } /// Add some minimal cuts for single track variables (min number of TPC clusters) - auto negExtra = casc.negTrackExtra_as>(); - auto posExtra = casc.posTrackExtra_as>(); - auto bachExtra = casc.bachTrackExtra_as>(); + auto negExtra = casc.negTrackExtra_as(); + auto posExtra = casc.posTrackExtra_as(); + auto bachExtra = casc.bachTrackExtra_as(); if (doNTPCSigmaCut) { if (casc.sign() < 0) { @@ -178,7 +251,7 @@ struct cascadeFlow { } } - void processTrainingSignal(soa::Join::iterator const& coll, soa::Join const& Cascades, soa::Join const&) + void processTrainingSignal(soa::Join::iterator const& coll, soa::Join const& Cascades, DauTracks const&) { if (!coll.sel8() || std::abs(coll.posZ()) > 10.) { @@ -191,9 +264,9 @@ struct cascadeFlow { && !(std::abs(pdgCode) == 3334 && std::abs(casc.pdgCodeV0()) == 3122 && std::abs(casc.pdgCodeBachelor()) == 321)) // Omega continue; - auto negExtra = casc.negTrackExtra_as>(); - auto posExtra = casc.posTrackExtra_as>(); - auto bachExtra = casc.bachTrackExtra_as>(); + auto negExtra = casc.negTrackExtra_as(); + auto posExtra = casc.posTrackExtra_as(); + auto bachExtra = casc.bachTrackExtra_as(); int counter = 0; IsCascAccepted(casc, negExtra, posExtra, bachExtra, counter); @@ -204,8 +277,71 @@ struct cascadeFlow { } } + void processAnalyseData(CollEventPlane const& coll, soa::Join const& Cascades, DauTracks const&) + { + + if (!coll.sel8() || std::abs(coll.posZ()) > 10.) { + return; + } + + histos.fill(HIST("hEventCentrality"), coll.centFT0M()); + histos.fill(HIST("hEventVertexZ"), coll.posZ()); + + eventplaneVecT0A = ROOT::Math::XYZVector(coll.qvecFT0ARe(), coll.qvecFT0AIm(), 0); + eventplaneVecT0C = ROOT::Math::XYZVector(coll.qvecFT0CRe(), coll.qvecFT0CIm(), 0); + //float eventplaneVecT0ANorm = sqrt(eventplaneVecT0A.Dot(eventplaneVecT0A)); + //float eventplaneVecT0CNorm = sqrt(eventplaneVecT0C.Dot(eventplaneVecT0C)); + // float v2Norm = eventplaneVecT0A.Dot(eventplaneVecT0C)/coll.sumAmplFT0A()/coll.sumAmplFT0C(); + float v2Norm = eventplaneVecT0A.Dot(eventplaneVecT0C); + //std::cout << "eventplaneVecT0ANorm: " << eventplaneVecT0ANorm << " = " << sqrt(pow(coll.qvecFT0ARe(), 2) + pow(coll.qvecFT0AIm(), 2)) << std::endl; //these two are equal + //std::cout << "stored norm value: " << coll.sumAmplFT0A() << std::endl; + histos.fill(HIST("hv2Norm"), v2Norm); + + float PsiT0A = TMath::ACos(coll.qvecFT0ARe()) / 2; + float PsiT0C = TMath::ACos(coll.qvecFT0CRe()) / 2; + if (coll.qvecFT0AIm() < 0) PsiT0A = -PsiT0A + TMath::Pi()/2; //to get dstribution between 0 and Pi + if (coll.qvecFT0CIm() < 0) PsiT0C = -PsiT0C + TMath::Pi()/2; //to get dstribution between 0 and Pi + std::cout << "PsiT0A : " << PsiT0A << std::endl; + std::cout << "PsiT0C : " << PsiT0C << std::endl; + + histos.fill(HIST("hFT0ARe"), coll.qvecFT0ARe()); + histos.fill(HIST("hFT0AIm"), coll.qvecFT0AIm()); + histos.fill(HIST("hFT0A"), coll.qvecFT0ARe(), coll.qvecFT0AIm()); + histos.fill(HIST("hPsiT0A"), PsiT0A); + histos.fill(HIST("hPsiT0C"), PsiT0C); + + for (auto& casc : Cascades) { + + /// Add some minimal cuts for single track variables (min number of TPC clusters) + auto negExtra = casc.negTrackExtra_as(); + auto posExtra = casc.posTrackExtra_as(); + auto bachExtra = casc.bachTrackExtra_as(); + + int counter = 0; + IsCascAccepted(casc, negExtra, posExtra, bachExtra, counter); + histos.fill(HIST("hCascade"), counter); + + cascphiVec = ROOT::Math::XYZVector(cos(2*casc.phi()), sin(2*casc.phi()), 0); + auto v2A = cascphiVec.Dot(eventplaneVecT0A); + auto v2C = cascphiVec.Dot(eventplaneVecT0C); + auto cascminuspsiT0A = GetPhiInRange(casc.phi() - PsiT0A); + auto cascminuspsiT0C = GetPhiInRange(casc.phi() - PsiT0C); + float v2A_diff = TMath::Cos(2.0 * cascminuspsiT0A); + float v2C_diff = TMath::Cos(2.0 * cascminuspsiT0C); + std::cout << "v2A: " << v2A << " v2C " << v2C << std::endl; + std::cout << "v2A_diff: " << v2A_diff << " v2C_diff " << v2C_diff << std::endl; + histos.fill(HIST("hCascadePhi"), casc.phi()); + histos.fill(HIST("hcascminuspsiT0A"), cascminuspsiT0A); + histos.fill(HIST("hcascminuspsiT0C"), cascminuspsiT0C); + fillAnalysedTable(coll, casc); + } + + } + + PROCESS_SWITCH(cascadeFlow, processTrainingBackground, "Process to create the training dataset for the background", true); PROCESS_SWITCH(cascadeFlow, processTrainingSignal, "Process to create the training dataset for the signal", false); + PROCESS_SWITCH(cascadeFlow, processAnalyseData, "Process to apply ML model to the data", false); }; WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) From a76d7f9b280408666bcd870f788608bc5e474832 Mon Sep 17 00:00:00 2001 From: Chiara De Martin Date: Fri, 8 Mar 2024 14:26:35 +0100 Subject: [PATCH 02/11] apply ML response --- PWGLF/TableProducer/CMakeLists.txt | 2 +- PWGLF/TableProducer/cascadeflow.cxx | 131 ++++++++++++++++++++++++++-- 2 files changed, 123 insertions(+), 10 deletions(-) diff --git a/PWGLF/TableProducer/CMakeLists.txt b/PWGLF/TableProducer/CMakeLists.txt index 0f2d9924d1a..bcd780ccb60 100644 --- a/PWGLF/TableProducer/CMakeLists.txt +++ b/PWGLF/TableProducer/CMakeLists.txt @@ -190,7 +190,7 @@ o2physics_add_dpl_workflow(cascqaanalysis o2physics_add_dpl_workflow(cascadeflow SOURCES cascadeflow.cxx - PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore + PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore O2Physics::MLCore COMPONENT_NAME Analysis) # Spectra diff --git a/PWGLF/TableProducer/cascadeflow.cxx b/PWGLF/TableProducer/cascadeflow.cxx index bdb9177bfe9..730e1efd3d1 100644 --- a/PWGLF/TableProducer/cascadeflow.cxx +++ b/PWGLF/TableProducer/cascadeflow.cxx @@ -27,7 +27,6 @@ #include "PWGLF/DataModel/LFStrangenessTables.h" #include "PWGLF/DataModel/LFStrangenessPIDTables.h" #include "Tools/ML/MlResponse.h" -#include "PWGHF/Core/HfHelper.h" using namespace o2; using namespace o2::framework; @@ -52,10 +51,70 @@ static const std::vector massSigmaParameterNames{"p0", "p1", "p2", static const std::vector speciesNames{"Xi", "Omega"}; } // namespace cascadev2 +namespace cascade_flow_cuts_ml +{ + // direction of the cut + enum CutDirection { + CutGreater = 0, // require score < cut value + CutSmaller, // require score > cut value + CutNot // do not cut on score + }; + + static constexpr int nBinsPt = 8; + static constexpr int nCutScores = 2; + // default values for the pT bin edges, offset by 1 from the bin numbers in cuts array + constexpr double binsPt[nBinsPt + 1] = { + 0.6, + 1., + 2., + 3., + 4., + 5., + 6., + 8., + 10.}; + auto vecBinsPt = std::vector{binsPt, binsPt + nBinsPt + 1}; + + // default values for the ML model paths, one model per pT bin + static const std::vector modelPaths = { + ""}; + + // default values for the cut directions + constexpr int cutDir[nCutScores] = {CutGreater, CutNot}; + auto vecCutDir = std::vector{cutDir, cutDir + nCutScores}; + + // default values for the cuts + constexpr double cuts[nBinsPt][nCutScores] = { + {0.5, 0.5}, + {0.5, 0.5}, + {0.5, 0.5}, + {0.5, 0.5}, + {0.5, 0.5}, + {0.5, 0.5}, + {0.5, 0.5}, + {0.5, 0.5}}; + + // row labels + static const std::vector labelsPt = { + "pT bin 0", + "pT bin 1", + "pT bin 2", + "pT bin 3", + "pT bin 4", + "pT bin 5", + "pT bin 6", + "pT bin 7"}; + + // column labels + static const std::vector labelsCutScore = {"score class 1", "score class 2"}; + static const std::vector labelsDmesCutScore = {"ML score signal", "ML score bkg"}; +} // namespace cascade_flow_cuts_ml + static constexpr double defaultCutsMl[1][2] = {{0.5, 0.5}}; struct cascadeFlow { + Configurable isOmega{"isOmega", 0, "Xi or Omega"}; Configurable sideBandStart{"sideBandStart", 5, "Start of the sideband region in number of sigmas"}; Configurable sideBandEnd{"sideBandEnd", 7, "End of the sideband region in number of sigmas"}; Configurable downsample{"downsample", 1., "Downsample training output tree"}; @@ -63,16 +122,21 @@ struct cascadeFlow { Configurable nsigmatpcPr{"nsigmatpcPr", 5, "nsigmatpcPr"}; Configurable nsigmatpcPi{"nsigmatpcPi", 5, "nsigmatpcPi"}; Configurable mintpccrrows{"mintpccrrows", 3, "mintpccrrows"}; - Configurable loadModelsFromCCDB{"loadModelsFromCCDB", false, "Flag to enable or disable the loading of models from CCDB"}; + Configurable ccdbUrl{"ccdbUrl", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; - Configurable modelPathsCCDB{"modelPathsCCDB", "Users/c/chdemart/CascadesFlow", "Path on CCDB"}; + Configurable> modelPathsCCDB{"modelPathsCCDB", std::vector{"Users/c/chdemart/CascadesFlow"}, "Paths of models on CCDB"}; + Configurable> onnxFileNames{"onnxFileNames", std::vector{"model_onnx.onnx"}, "ONNX file names for each pT bin (if not from CCDB full path)"}; Configurable timestampCCDB{"timestampCCDB", -1, "timestamp of the ONNX file for ML model used to query in CCDB"}; - Configurable> binsPtMl{"binsPtMl", std::vector{0.6, 10.}, "pT bin limits for ML application"}; - Configurable> cutDirMl{"cutDirMl", std::vector{cuts_ml::CutSmaller, cuts_ml::CutNot}, "Whether to reject score values greater or smaller than the threshold"}; - Configurable> cutsMl{"cutsMl", {defaultCutsMl[0], 1, 3, {"pT bin 0"}, {"score signal", "score bkg"}}, "ML selections per pT bin"}; - Configurable nClassesMl{"nClassesMl", (int8_t)2, "Number of classes in ML model"}; + Configurable loadModelsFromCCDB{"loadModelsFromCCDB", false, "Flag to enable or disable the loading of models from CCDB"}; + + // ML inference + Configurable applyMl{"applyMl", true, "Flag to apply ML selections"}; + Configurable> binsPtMl{"binsPtMl", std::vector{cascade_flow_cuts_ml::vecBinsPt}, "pT bin limits for ML application"}; + Configurable> cutDirMl{"cutDirMl", std::vector{cascade_flow_cuts_ml::vecCutDir}, "Whether to reject score values greater or smaller than the threshold"}; + Configurable> cutsMl{"cutsMl", {cascade_flow_cuts_ml::cuts[0], cascade_flow_cuts_ml::nBinsPt, cascade_flow_cuts_ml::nCutScores, cascade_flow_cuts_ml::labelsPt, cascade_flow_cuts_ml::labelsCutScore}, "ML selections per pT bin"}; + Configurable nClassesMl{"nClassesMl", (int8_t)cascade_flow_cuts_ml::nCutScores, "Number of classes in ML model"}; + // Configurable> namesInputFeatures{"namesInputFeatures", std::vector{"feature1", "feature2"}, "Names of ML model input features"}; - HfHelper hfHelper; o2::ccdb::CcdbApi ccdbApi; // Add objects needed for ML inference @@ -171,7 +235,12 @@ struct cascadeFlow { void init(InitContext const&) { + float MinMass = 1.28; + if (isOmega) MinMass = 1.6; + float MaxMass = 1.36; + if (isOmega) MaxMass = 1.73; ConfigurableAxis vertexZ{"vertexZ", {20, -10, 10}, "vertex axis for histograms"}; + ConfigurableAxis massCascAxis{"massCascAxis", {100, MinMass, MaxMass}, ""}; histos.add("hEventVertexZ", "hEventVertexZ", kTH1F, {vertexZ}); histos.add("hEventCentrality", "hEventCentrality", kTH1F, {{101, 0, 101}}); histos.add("hCandidate", "hCandidate", HistType::kTH1D, {{22, -0.5, 21.5}}); @@ -186,6 +255,12 @@ struct cascadeFlow { histos.add("hFT0AIm", "hFT0AIm", HistType::kTH1D, {{100, -1, 1}}); histos.add("hFT0A", "hFT0A", HistType::kTH2D, {{100, -1, 1}, {100, -1, 1}}); histos.add("hv2Norm", "hv2Norm", HistType::kTH1D, {{100, 0, 1}}); + histos.add("hMassBeforeSel", "hMassBeforeSel", HistType::kTH1F, {massCascAxis}); + histos.add("hMassAfterSel", "hMassAfterSel", HistType::kTH1F, {massCascAxis}); + histos.add("hSignalScoreBeforeSel", "Signal score before selection;BDT first score;entries", HistType::kTH1F, {{100, 0., 1.}}); + histos.add("hBkgScoreBeforeSel", "Bkg score before selection;BDT first score;entries", HistType::kTH1F, {{100, 0., 1.}}); + histos.add("hSignalScoreAfterSel", "Signal score after selection;BDT first score;entries", HistType::kTH1F, {{100, 0., 1.}}); + histos.add("hBkgScoreAfterSel", "Bkg score after selection;BDT first score;entries", HistType::kTH1F, {{100, 0., 1.}}); // Configure and initialise the ML class mlResponse.configure(binsPtMl, cutsMl, cutDirMl, nClassesMl); @@ -321,6 +396,43 @@ struct cascadeFlow { IsCascAccepted(casc, negExtra, posExtra, bachExtra, counter); histos.fill(HIST("hCascade"), counter); + // ML selections + bool isSelectedCasc = false; + + std::vector inputFeaturesCasc{casc.cascradius(), + casc.v0radius(), + casc.casccosPA(coll.posX(), coll.posY(), coll.posZ()), + casc.v0cosPA(coll.posX(), coll.posY(), coll.posZ()), + casc.dcapostopv(), + casc.dcanegtopv(), + casc.dcabachtopv(), + casc.dcacascdaughters(), + casc.dcaV0daughters(), + casc.dcav0topv(coll.posX(), coll.posY(), coll.posZ()), + casc.bachBaryonCosPA(), + casc.bachBaryonDCAxyToPV()}; + + // Retrieve model output and selection outcome + isSelectedCasc = mlResponse.isSelectedMl(inputFeaturesCasc, casc.pt(), outputMl); + + // Fill BDT score histograms before selection + histos.fill(HIST("hSignalScoreBeforeSel"), outputMl[0]); + histos.fill(HIST("hBkgScoreBeforeSel"), outputMl[1]); + + // Fill histograms for selected candidates + float massCasc = casc.mXi(); + if (isOmega) massCasc = casc.mOmega(); + if (isSelectedCasc) { + histos.fill(HIST("hMassAfterSel"), massCasc); + histos.fill(HIST("hSignalScoreAfterSel"), outputMl[0]); + histos.fill(HIST("hBkgScoreAfterSel"), outputMl[1]); + histos.fill(HIST("hMassAfterSelVsPt"), massCasc, casc.pt()); + histos.fill(HIST("hPromptScoreAfterSelVsPt"), outputMl[0], casc.pt()); + histos.fill(HIST("hBkgScoreAfterSelVsPt"), outputMl[1], casc.pt()); + } + + outputMl.clear(); // not necessary in this case but for good measure + cascphiVec = ROOT::Math::XYZVector(cos(2*casc.phi()), sin(2*casc.phi()), 0); auto v2A = cascphiVec.Dot(eventplaneVecT0A); auto v2C = cascphiVec.Dot(eventplaneVecT0C); @@ -333,7 +445,8 @@ struct cascadeFlow { histos.fill(HIST("hCascadePhi"), casc.phi()); histos.fill(HIST("hcascminuspsiT0A"), cascminuspsiT0A); histos.fill(HIST("hcascminuspsiT0C"), cascminuspsiT0C); - fillAnalysedTable(coll, casc); + + if (isSelectedCasc) fillAnalysedTable(coll, casc); } } From cfb4f5fdffb1918eac69b4a3f9f9a71c4425cb1b Mon Sep 17 00:00:00 2001 From: Chiara De Martin Date: Fri, 8 Mar 2024 14:42:56 +0100 Subject: [PATCH 03/11] introduce pt dependent onnx file --- PWGLF/TableProducer/cascadeflow.cxx | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/PWGLF/TableProducer/cascadeflow.cxx b/PWGLF/TableProducer/cascadeflow.cxx index 730e1efd3d1..a486e52a793 100644 --- a/PWGLF/TableProducer/cascadeflow.cxx +++ b/PWGLF/TableProducer/cascadeflow.cxx @@ -75,6 +75,16 @@ namespace cascade_flow_cuts_ml 10.}; auto vecBinsPt = std::vector{binsPt, binsPt + nBinsPt + 1}; + static const std::vector onnxBinsPt = { + "model_onnx0.6.onnx", + "model_onnx1.onnx", + "model_onnx2.onnx", + "model_onnx3.onnx", + "model_onnx4.onnx", + "model_onnx5.onnx", + "model_onnx6.onnx", + "model_onnx8.onnx"}; + // default values for the ML model paths, one model per pT bin static const std::vector modelPaths = { ""}; @@ -125,7 +135,7 @@ struct cascadeFlow { Configurable ccdbUrl{"ccdbUrl", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; Configurable> modelPathsCCDB{"modelPathsCCDB", std::vector{"Users/c/chdemart/CascadesFlow"}, "Paths of models on CCDB"}; - Configurable> onnxFileNames{"onnxFileNames", std::vector{"model_onnx.onnx"}, "ONNX file names for each pT bin (if not from CCDB full path)"}; + Configurable> onnxFileNames{"onnxFileNames", std::vector{cascade_flow_cuts_ml::onnxBinsPt}, "ONNX file names for each pT bin (if not from CCDB full path)"}; Configurable timestampCCDB{"timestampCCDB", -1, "timestamp of the ONNX file for ML model used to query in CCDB"}; Configurable loadModelsFromCCDB{"loadModelsFromCCDB", false, "Flag to enable or disable the loading of models from CCDB"}; From 7d29d06c4062e632c68a57d3242516c04b34d99f Mon Sep 17 00:00:00 2001 From: Chiara De Martin Date: Fri, 8 Mar 2024 14:47:57 +0100 Subject: [PATCH 04/11] boolean to chose whether to apply ML model or not --- PWGLF/TableProducer/cascadeflow.cxx | 65 ++++++++++++++++------------- 1 file changed, 36 insertions(+), 29 deletions(-) diff --git a/PWGLF/TableProducer/cascadeflow.cxx b/PWGLF/TableProducer/cascadeflow.cxx index a486e52a793..89fcfc1011b 100644 --- a/PWGLF/TableProducer/cascadeflow.cxx +++ b/PWGLF/TableProducer/cascadeflow.cxx @@ -125,6 +125,7 @@ static constexpr double defaultCutsMl[1][2] = {{0.5, 0.5}}; struct cascadeFlow { Configurable isOmega{"isOmega", 0, "Xi or Omega"}; + Configurable isApplyML{"isApplyML", 1, ""}; Configurable sideBandStart{"sideBandStart", 5, "Start of the sideband region in number of sigmas"}; Configurable sideBandEnd{"sideBandEnd", 7, "End of the sideband region in number of sigmas"}; Configurable downsample{"downsample", 1., "Downsample training output tree"}; @@ -272,18 +273,20 @@ struct cascadeFlow { histos.add("hSignalScoreAfterSel", "Signal score after selection;BDT first score;entries", HistType::kTH1F, {{100, 0., 1.}}); histos.add("hBkgScoreAfterSel", "Bkg score after selection;BDT first score;entries", HistType::kTH1F, {{100, 0., 1.}}); - // Configure and initialise the ML class - mlResponse.configure(binsPtMl, cutsMl, cutDirMl, nClassesMl); + if (isApplyML){ + // Configure and initialise the ML class + mlResponse.configure(binsPtMl, cutsMl, cutDirMl, nClassesMl); - // Bonus: retrieve the model from CCDB (needed for ML application on the GRID) - if (loadModelsFromCCDB) { - ccdbApi.init(ccdbUrl); - mlResponse.setModelPathsCCDB(onnxFileNames, ccdbApi, modelPathsCCDB.value, timestampCCDB); - } else { - mlResponse.setModelPathsLocal(onnxFileNames); - } + // Bonus: retrieve the model from CCDB (needed for ML application on the GRID) + if (loadModelsFromCCDB) { + ccdbApi.init(ccdbUrl); + mlResponse.setModelPathsCCDB(onnxFileNames, ccdbApi, modelPathsCCDB.value, timestampCCDB); + } else { + mlResponse.setModelPathsLocal(onnxFileNames); + } - mlResponse.init(); + mlResponse.init(); + } } void processTrainingBackground(soa::Join::iterator const& coll, soa::Join const& Cascades, DauTracks const&) @@ -422,27 +425,31 @@ struct cascadeFlow { casc.bachBaryonCosPA(), casc.bachBaryonDCAxyToPV()}; - // Retrieve model output and selection outcome - isSelectedCasc = mlResponse.isSelectedMl(inputFeaturesCasc, casc.pt(), outputMl); - - // Fill BDT score histograms before selection - histos.fill(HIST("hSignalScoreBeforeSel"), outputMl[0]); - histos.fill(HIST("hBkgScoreBeforeSel"), outputMl[1]); - - // Fill histograms for selected candidates - float massCasc = casc.mXi(); - if (isOmega) massCasc = casc.mOmega(); - if (isSelectedCasc) { - histos.fill(HIST("hMassAfterSel"), massCasc); - histos.fill(HIST("hSignalScoreAfterSel"), outputMl[0]); - histos.fill(HIST("hBkgScoreAfterSel"), outputMl[1]); - histos.fill(HIST("hMassAfterSelVsPt"), massCasc, casc.pt()); - histos.fill(HIST("hPromptScoreAfterSelVsPt"), outputMl[0], casc.pt()); - histos.fill(HIST("hBkgScoreAfterSelVsPt"), outputMl[1], casc.pt()); + if (isApplyML){ + // Retrieve model output and selection outcome + isSelectedCasc = mlResponse.isSelectedMl(inputFeaturesCasc, casc.pt(), outputMl); + + // Fill BDT score histograms before selection + histos.fill(HIST("hSignalScoreBeforeSel"), outputMl[0]); + histos.fill(HIST("hBkgScoreBeforeSel"), outputMl[1]); + + // Fill histograms for selected candidates + float massCasc = casc.mXi(); + if (isOmega) massCasc = casc.mOmega(); + if (isSelectedCasc) { + histos.fill(HIST("hMassAfterSel"), massCasc); + histos.fill(HIST("hSignalScoreAfterSel"), outputMl[0]); + histos.fill(HIST("hBkgScoreAfterSel"), outputMl[1]); + histos.fill(HIST("hMassAfterSelVsPt"), massCasc, casc.pt()); + histos.fill(HIST("hPromptScoreAfterSelVsPt"), outputMl[0], casc.pt()); + histos.fill(HIST("hBkgScoreAfterSelVsPt"), outputMl[1], casc.pt()); + } + + outputMl.clear(); // not necessary in this case but for good measure + } else { + isSelectedCasc = true; } - outputMl.clear(); // not necessary in this case but for good measure - cascphiVec = ROOT::Math::XYZVector(cos(2*casc.phi()), sin(2*casc.phi()), 0); auto v2A = cascphiVec.Dot(eventplaneVecT0A); auto v2C = cascphiVec.Dot(eventplaneVecT0C); From f28046197193bca68296fc1c71123551d5ab1bdf Mon Sep 17 00:00:00 2001 From: Chiara De Martin Date: Tue, 12 Mar 2024 10:47:54 +0100 Subject: [PATCH 05/11] fix v2 calculation --- PWGLF/TableProducer/cascadeflow.cxx | 38 ++++++++++++++++++++--------- 1 file changed, 26 insertions(+), 12 deletions(-) diff --git a/PWGLF/TableProducer/cascadeflow.cxx b/PWGLF/TableProducer/cascadeflow.cxx index 89fcfc1011b..5dbff8b81ea 100644 --- a/PWGLF/TableProducer/cascadeflow.cxx +++ b/PWGLF/TableProducer/cascadeflow.cxx @@ -232,7 +232,7 @@ struct cascadeFlow { template // void fillAnalisedTable(collision_t coll, cascade_t casc, float BDTresponse) - void fillAnalysedTable(collision_t coll, cascade_t casc) + void fillAnalysedTable(collision_t coll, cascade_t casc, float v2A, float v2C, float BDTresponse) { analysisSample(coll.centFT0M(), casc.sign(), @@ -240,7 +240,9 @@ struct cascadeFlow { casc.eta(), casc.mXi(), casc.mOmega(), - casc.mLambda()); + v2A, + v2C, + BDTresponse); } void init(InitContext const&) @@ -252,6 +254,7 @@ struct cascadeFlow { if (isOmega) MaxMass = 1.73; ConfigurableAxis vertexZ{"vertexZ", {20, -10, 10}, "vertex axis for histograms"}; ConfigurableAxis massCascAxis{"massCascAxis", {100, MinMass, MaxMass}, ""}; + ConfigurableAxis ptAxis{"ptAxis", {100, 0, 10}, ""}; histos.add("hEventVertexZ", "hEventVertexZ", kTH1F, {vertexZ}); histos.add("hEventCentrality", "hEventCentrality", kTH1F, {{101, 0, 101}}); histos.add("hCandidate", "hCandidate", HistType::kTH1D, {{22, -0.5, 21.5}}); @@ -268,6 +271,8 @@ struct cascadeFlow { histos.add("hv2Norm", "hv2Norm", HistType::kTH1D, {{100, 0, 1}}); histos.add("hMassBeforeSel", "hMassBeforeSel", HistType::kTH1F, {massCascAxis}); histos.add("hMassAfterSel", "hMassAfterSel", HistType::kTH1F, {massCascAxis}); + histos.add("hMassBeforeSelVsPt", "hMassBeforeSelVsPt", HistType::kTH2F, {{massCascAxis}, {ptAxis}}); + histos.add("hMassAfterSelVsPt", "hMassAfterSelVsPt", HistType::kTH2F, {{massCascAxis}, {ptAxis}}); histos.add("hSignalScoreBeforeSel", "Signal score before selection;BDT first score;entries", HistType::kTH1F, {{100, 0., 1.}}); histos.add("hBkgScoreBeforeSel", "Bkg score before selection;BDT first score;entries", HistType::kTH1F, {{100, 0., 1.}}); histos.add("hSignalScoreAfterSel", "Signal score after selection;BDT first score;entries", HistType::kTH1F, {{100, 0., 1.}}); @@ -377,10 +382,10 @@ struct cascadeFlow { eventplaneVecT0A = ROOT::Math::XYZVector(coll.qvecFT0ARe(), coll.qvecFT0AIm(), 0); eventplaneVecT0C = ROOT::Math::XYZVector(coll.qvecFT0CRe(), coll.qvecFT0CIm(), 0); - //float eventplaneVecT0ANorm = sqrt(eventplaneVecT0A.Dot(eventplaneVecT0A)); - //float eventplaneVecT0CNorm = sqrt(eventplaneVecT0C.Dot(eventplaneVecT0C)); + float eventplaneVecT0ANorm = sqrt(eventplaneVecT0A.Dot(eventplaneVecT0A)); + float eventplaneVecT0CNorm = sqrt(eventplaneVecT0C.Dot(eventplaneVecT0C)); // float v2Norm = eventplaneVecT0A.Dot(eventplaneVecT0C)/coll.sumAmplFT0A()/coll.sumAmplFT0C(); - float v2Norm = eventplaneVecT0A.Dot(eventplaneVecT0C); + float v2Norm = eventplaneVecT0A.Dot(eventplaneVecT0C) / eventplaneVecT0ANorm / eventplaneVecT0CNorm; //std::cout << "eventplaneVecT0ANorm: " << eventplaneVecT0ANorm << " = " << sqrt(pow(coll.qvecFT0ARe(), 2) + pow(coll.qvecFT0AIm(), 2)) << std::endl; //these two are equal //std::cout << "stored norm value: " << coll.sumAmplFT0A() << std::endl; histos.fill(HIST("hv2Norm"), v2Norm); @@ -425,6 +430,12 @@ struct cascadeFlow { casc.bachBaryonCosPA(), casc.bachBaryonDCAxyToPV()}; + float massCasc = casc.mXi(); + if (isOmega) massCasc = casc.mOmega(); + + histos.fill(HIST("hMassAfterSel"), massCasc); + histos.fill(HIST("hMassAfterSelVsPt"), massCasc, casc.pt()); + if (isApplyML){ // Retrieve model output and selection outcome isSelectedCasc = mlResponse.isSelectedMl(inputFeaturesCasc, casc.pt(), outputMl); @@ -434,13 +445,9 @@ struct cascadeFlow { histos.fill(HIST("hBkgScoreBeforeSel"), outputMl[1]); // Fill histograms for selected candidates - float massCasc = casc.mXi(); - if (isOmega) massCasc = casc.mOmega(); if (isSelectedCasc) { - histos.fill(HIST("hMassAfterSel"), massCasc); histos.fill(HIST("hSignalScoreAfterSel"), outputMl[0]); histos.fill(HIST("hBkgScoreAfterSel"), outputMl[1]); - histos.fill(HIST("hMassAfterSelVsPt"), massCasc, casc.pt()); histos.fill(HIST("hPromptScoreAfterSelVsPt"), outputMl[0], casc.pt()); histos.fill(HIST("hBkgScoreAfterSelVsPt"), outputMl[1], casc.pt()); } @@ -450,9 +457,14 @@ struct cascadeFlow { isSelectedCasc = true; } + if (isSelectedCasc) { + histos.fill(HIST("hMassAfterSel"), massCasc); + histos.fill(HIST("hMassAfterSelVsPt"), massCasc, casc.pt()); + } + cascphiVec = ROOT::Math::XYZVector(cos(2*casc.phi()), sin(2*casc.phi()), 0); - auto v2A = cascphiVec.Dot(eventplaneVecT0A); - auto v2C = cascphiVec.Dot(eventplaneVecT0C); + auto v2A = cascphiVec.Dot(eventplaneVecT0A)/sqrt(eventplaneVecT0A.Dot(eventplaneVecT0A)); + auto v2C = cascphiVec.Dot(eventplaneVecT0C)/sqrt(eventplaneVecT0C.Dot(eventplaneVecT0C)); auto cascminuspsiT0A = GetPhiInRange(casc.phi() - PsiT0A); auto cascminuspsiT0C = GetPhiInRange(casc.phi() - PsiT0C); float v2A_diff = TMath::Cos(2.0 * cascminuspsiT0A); @@ -463,7 +475,9 @@ struct cascadeFlow { histos.fill(HIST("hcascminuspsiT0A"), cascminuspsiT0A); histos.fill(HIST("hcascminuspsiT0C"), cascminuspsiT0C); - if (isSelectedCasc) fillAnalysedTable(coll, casc); + float BDTresponse = 0; + if (isApplyML) BDTresponse = outputMl[0]; + if (isSelectedCasc) fillAnalysedTable(coll, casc, v2A, v2C, BDTresponse); } } From 3ba0882ee25675a80280adad90a7c25cb5da24e1 Mon Sep 17 00:00:00 2001 From: Chiara De Martin Date: Tue, 12 Mar 2024 10:48:40 +0100 Subject: [PATCH 06/11] define new table for flow --- PWGLF/DataModel/cascqaanalysis.h | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) diff --git a/PWGLF/DataModel/cascqaanalysis.h b/PWGLF/DataModel/cascqaanalysis.h index 9efe2d9976a..370c4c7bebb 100644 --- a/PWGLF/DataModel/cascqaanalysis.h +++ b/PWGLF/DataModel/cascqaanalysis.h @@ -109,6 +109,21 @@ DECLARE_SOA_DYNAMIC_COLUMN(IsINELgt1, isINELgt1, //! True if the Event belongs t [](uint8_t flags) -> bool { return (flags & EvFlags::EvINELgt1) == EvFlags::EvINELgt1; }); } // namespace mycascades +namespace cascadesflow +{ + +DECLARE_SOA_COLUMN(MultFT0M, multFT0M, float); +DECLARE_SOA_COLUMN(Sign, sign, int); +DECLARE_SOA_COLUMN(Pt, pt, float); +DECLARE_SOA_COLUMN(Eta, eta, float); +DECLARE_SOA_COLUMN(MassXi, massxi, float); +DECLARE_SOA_COLUMN(MassOmega, massomega, float); +DECLARE_SOA_COLUMN(V2A, v2A, float); +DECLARE_SOA_COLUMN(V2C, v2C, float); +DECLARE_SOA_COLUMN(BDTResponse, bdtResponse, float); + +} // namespace cascadesflow + DECLARE_SOA_TABLE(MyCascades, "AOD", "MYCASCADES", o2::soa::Index<>, mycascades::CollisionZ, mycascades::MultFT0M, mycascades::MultFV0A, mycascades::Sign, mycascades::Pt, mycascades::RapXi, mycascades::RapOmega, mycascades::Eta, mycascades::MassXi, mycascades::MassOmega, mycascades::MassLambdaDau, mycascades::CascRadius, mycascades::V0Radius, mycascades::CascCosPA, mycascades::V0CosPA, mycascades::DCAPosToPV, mycascades::DCANegToPV, @@ -133,6 +148,9 @@ DECLARE_SOA_TABLE(CascTraining, "AOD", "CascTraining", o2::soa::Index<>, mycascades::MultFT0M, mycascades::Sign, mycascades::Pt, mycascades::Eta, mycascades::MassXi, mycascades::MassOmega, mycascades::MassLambdaDau, mycascades::CascRadius, mycascades::V0Radius, mycascades::CascCosPA, mycascades::V0CosPA, mycascades::DCAPosToPV, mycascades::DCANegToPV, mycascades::DCABachToPV, mycascades::DCACascDaughters, mycascades::DCAV0Daughters, mycascades::DCAV0ToPV, mycascades::BachBaryonCosPA, mycascades::BachBaryonDCAxyToPV, mycascades::McPdgCode); +DECLARE_SOA_TABLE(CascAnalysis, "AOD", "CascAnalysis", o2::soa::Index<>, + cascadesflow::MultFT0M, cascadesflow::Sign, cascadesflow::Pt, cascadesflow::Eta, cascadesflow::MassXi, cascadesflow::MassOmega, cascadesflow::V2A, cascadesflow::V2C, cascadesflow::BDTResponse); + namespace myMCcascades { From 51c97b648e06be3ed98fc86f3b40611c0faf19b5 Mon Sep 17 00:00:00 2001 From: Chiara De Martin Date: Tue, 12 Mar 2024 11:19:48 +0100 Subject: [PATCH 07/11] fix bug in filling histogram --- PWGLF/TableProducer/cascadeflow.cxx | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/PWGLF/TableProducer/cascadeflow.cxx b/PWGLF/TableProducer/cascadeflow.cxx index 5dbff8b81ea..99e9ce5bc6a 100644 --- a/PWGLF/TableProducer/cascadeflow.cxx +++ b/PWGLF/TableProducer/cascadeflow.cxx @@ -433,8 +433,8 @@ struct cascadeFlow { float massCasc = casc.mXi(); if (isOmega) massCasc = casc.mOmega(); - histos.fill(HIST("hMassAfterSel"), massCasc); - histos.fill(HIST("hMassAfterSelVsPt"), massCasc, casc.pt()); + histos.fill(HIST("hMassBeforeSel"), massCasc); + histos.fill(HIST("hMassBeforeSelVsPt"), massCasc, casc.pt()); if (isApplyML){ // Retrieve model output and selection outcome From 94eade89ed87318f7591c87c9a575d5de074c985 Mon Sep 17 00:00:00 2001 From: Chiara De Martin Date: Wed, 13 Mar 2024 09:40:12 +0100 Subject: [PATCH 08/11] fix configurables + histograms --- PWGLF/TableProducer/cascadeflow.cxx | 46 +++++++++++++---------------- 1 file changed, 21 insertions(+), 25 deletions(-) diff --git a/PWGLF/TableProducer/cascadeflow.cxx b/PWGLF/TableProducer/cascadeflow.cxx index 99e9ce5bc6a..6e22b06b5cd 100644 --- a/PWGLF/TableProducer/cascadeflow.cxx +++ b/PWGLF/TableProducer/cascadeflow.cxx @@ -30,6 +30,7 @@ using namespace o2; using namespace o2::framework; +using namespace o2::analysis; using namespace o2::framework::expressions; using std::array; @@ -75,19 +76,8 @@ namespace cascade_flow_cuts_ml 10.}; auto vecBinsPt = std::vector{binsPt, binsPt + nBinsPt + 1}; - static const std::vector onnxBinsPt = { - "model_onnx0.6.onnx", - "model_onnx1.onnx", - "model_onnx2.onnx", - "model_onnx3.onnx", - "model_onnx4.onnx", - "model_onnx5.onnx", - "model_onnx6.onnx", - "model_onnx8.onnx"}; - // default values for the ML model paths, one model per pT bin - static const std::vector modelPaths = { - ""}; + static const std::vector modelPaths = {""}; // default values for the cut directions constexpr int cutDir[nCutScores] = {CutGreater, CutNot}; @@ -125,6 +115,8 @@ static constexpr double defaultCutsMl[1][2] = {{0.5, 0.5}}; struct cascadeFlow { Configurable isOmega{"isOmega", 0, "Xi or Omega"}; + Configurable MinPt{"MinPt", 0.6, "Min pt of cascade"}; + Configurable MaxPt{"MaxPt", 10, "Max pt of cascade"}; Configurable isApplyML{"isApplyML", 1, ""}; Configurable sideBandStart{"sideBandStart", 5, "Start of the sideband region in number of sigmas"}; Configurable sideBandEnd{"sideBandEnd", 7, "End of the sideband region in number of sigmas"}; @@ -136,7 +128,7 @@ struct cascadeFlow { Configurable ccdbUrl{"ccdbUrl", "http://alice-ccdb.cern.ch", "url of the ccdb repository"}; Configurable> modelPathsCCDB{"modelPathsCCDB", std::vector{"Users/c/chdemart/CascadesFlow"}, "Paths of models on CCDB"}; - Configurable> onnxFileNames{"onnxFileNames", std::vector{cascade_flow_cuts_ml::onnxBinsPt}, "ONNX file names for each pT bin (if not from CCDB full path)"}; + Configurable> onnxFileNames{"onnxFileNames", std::vector{"model_onnx.onnx"}, "ONNX file names for each pT bin (if not from CCDB full path)"}; Configurable timestampCCDB{"timestampCCDB", -1, "timestamp of the ONNX file for ML model used to query in CCDB"}; Configurable loadModelsFromCCDB{"loadModelsFromCCDB", false, "Flag to enable or disable the loading of models from CCDB"}; @@ -281,15 +273,13 @@ struct cascadeFlow { if (isApplyML){ // Configure and initialise the ML class mlResponse.configure(binsPtMl, cutsMl, cutDirMl, nClassesMl); - // Bonus: retrieve the model from CCDB (needed for ML application on the GRID) if (loadModelsFromCCDB) { ccdbApi.init(ccdbUrl); - mlResponse.setModelPathsCCDB(onnxFileNames, ccdbApi, modelPathsCCDB.value, timestampCCDB); + mlResponse.setModelPathsCCDB(onnxFileNames, ccdbApi, modelPathsCCDB, timestampCCDB); } else { mlResponse.setModelPathsLocal(onnxFileNames); } - mlResponse.init(); } } @@ -373,6 +363,8 @@ struct cascadeFlow { void processAnalyseData(CollEventPlane const& coll, soa::Join const& Cascades, DauTracks const&) { + //std::cout << "onnxFileNames " << onnxFileNames.value[0] << std::endl; + if (!coll.sel8() || std::abs(coll.posZ()) > 10.) { return; } @@ -394,8 +386,6 @@ struct cascadeFlow { float PsiT0C = TMath::ACos(coll.qvecFT0CRe()) / 2; if (coll.qvecFT0AIm() < 0) PsiT0A = -PsiT0A + TMath::Pi()/2; //to get dstribution between 0 and Pi if (coll.qvecFT0CIm() < 0) PsiT0C = -PsiT0C + TMath::Pi()/2; //to get dstribution between 0 and Pi - std::cout << "PsiT0A : " << PsiT0A << std::endl; - std::cout << "PsiT0C : " << PsiT0C << std::endl; histos.fill(HIST("hFT0ARe"), coll.qvecFT0ARe()); histos.fill(HIST("hFT0AIm"), coll.qvecFT0AIm()); @@ -433,13 +423,21 @@ struct cascadeFlow { float massCasc = casc.mXi(); if (isOmega) massCasc = casc.mOmega(); + //inv mass loose cut + if (casc.pt() < MinPt) continue; + if (casc.pt() > MaxPt) continue; + if (isOmega) { + if (massCasc < 1.6 || massCasc > 1.73) continue; + } else { + if (massCasc < 1.28 || massCasc > 1.36) continue; + } + histos.fill(HIST("hMassBeforeSel"), massCasc); histos.fill(HIST("hMassBeforeSelVsPt"), massCasc, casc.pt()); if (isApplyML){ // Retrieve model output and selection outcome isSelectedCasc = mlResponse.isSelectedMl(inputFeaturesCasc, casc.pt(), outputMl); - // Fill BDT score histograms before selection histos.fill(HIST("hSignalScoreBeforeSel"), outputMl[0]); histos.fill(HIST("hBkgScoreBeforeSel"), outputMl[1]); @@ -448,8 +446,6 @@ struct cascadeFlow { if (isSelectedCasc) { histos.fill(HIST("hSignalScoreAfterSel"), outputMl[0]); histos.fill(HIST("hBkgScoreAfterSel"), outputMl[1]); - histos.fill(HIST("hPromptScoreAfterSelVsPt"), outputMl[0], casc.pt()); - histos.fill(HIST("hBkgScoreAfterSelVsPt"), outputMl[1], casc.pt()); } outputMl.clear(); // not necessary in this case but for good measure @@ -467,10 +463,10 @@ struct cascadeFlow { auto v2C = cascphiVec.Dot(eventplaneVecT0C)/sqrt(eventplaneVecT0C.Dot(eventplaneVecT0C)); auto cascminuspsiT0A = GetPhiInRange(casc.phi() - PsiT0A); auto cascminuspsiT0C = GetPhiInRange(casc.phi() - PsiT0C); - float v2A_diff = TMath::Cos(2.0 * cascminuspsiT0A); - float v2C_diff = TMath::Cos(2.0 * cascminuspsiT0C); - std::cout << "v2A: " << v2A << " v2C " << v2C << std::endl; - std::cout << "v2A_diff: " << v2A_diff << " v2C_diff " << v2C_diff << std::endl; + //float v2A_diff = TMath::Cos(2.0 * cascminuspsiT0A); + //float v2C_diff = TMath::Cos(2.0 * cascminuspsiT0C); + //std::cout << "v2A: " << v2A << " v2C " << v2C << std::endl; + //std::cout << "v2A_diff: " << v2A_diff << " v2C_diff " << v2C_diff << std::endl; histos.fill(HIST("hCascadePhi"), casc.phi()); histos.fill(HIST("hcascminuspsiT0A"), cascminuspsiT0A); histos.fill(HIST("hcascminuspsiT0C"), cascminuspsiT0C); From 380bef420a4eae1b539974234ef1058e9ea2cd32 Mon Sep 17 00:00:00 2001 From: ALICE Action Bot Date: Wed, 13 Mar 2024 11:05:09 +0000 Subject: [PATCH 09/11] Please consider the following formatting changes --- PWGLF/TableProducer/cascadeflow.cxx | 243 ++++++++++++++-------------- 1 file changed, 126 insertions(+), 117 deletions(-) diff --git a/PWGLF/TableProducer/cascadeflow.cxx b/PWGLF/TableProducer/cascadeflow.cxx index 6e22b06b5cd..20391cbbf28 100644 --- a/PWGLF/TableProducer/cascadeflow.cxx +++ b/PWGLF/TableProducer/cascadeflow.cxx @@ -54,60 +54,60 @@ static const std::vector speciesNames{"Xi", "Omega"}; namespace cascade_flow_cuts_ml { - // direction of the cut - enum CutDirection { - CutGreater = 0, // require score < cut value - CutSmaller, // require score > cut value - CutNot // do not cut on score - }; - - static constexpr int nBinsPt = 8; - static constexpr int nCutScores = 2; - // default values for the pT bin edges, offset by 1 from the bin numbers in cuts array - constexpr double binsPt[nBinsPt + 1] = { - 0.6, - 1., - 2., - 3., - 4., - 5., - 6., - 8., - 10.}; - auto vecBinsPt = std::vector{binsPt, binsPt + nBinsPt + 1}; - - // default values for the ML model paths, one model per pT bin - static const std::vector modelPaths = {""}; - - // default values for the cut directions - constexpr int cutDir[nCutScores] = {CutGreater, CutNot}; - auto vecCutDir = std::vector{cutDir, cutDir + nCutScores}; - - // default values for the cuts - constexpr double cuts[nBinsPt][nCutScores] = { - {0.5, 0.5}, - {0.5, 0.5}, - {0.5, 0.5}, - {0.5, 0.5}, - {0.5, 0.5}, - {0.5, 0.5}, - {0.5, 0.5}, - {0.5, 0.5}}; - - // row labels - static const std::vector labelsPt = { - "pT bin 0", - "pT bin 1", - "pT bin 2", - "pT bin 3", - "pT bin 4", - "pT bin 5", - "pT bin 6", - "pT bin 7"}; - - // column labels - static const std::vector labelsCutScore = {"score class 1", "score class 2"}; - static const std::vector labelsDmesCutScore = {"ML score signal", "ML score bkg"}; +// direction of the cut +enum CutDirection { + CutGreater = 0, // require score < cut value + CutSmaller, // require score > cut value + CutNot // do not cut on score +}; + +static constexpr int nBinsPt = 8; +static constexpr int nCutScores = 2; +// default values for the pT bin edges, offset by 1 from the bin numbers in cuts array +constexpr double binsPt[nBinsPt + 1] = { + 0.6, + 1., + 2., + 3., + 4., + 5., + 6., + 8., + 10.}; +auto vecBinsPt = std::vector{binsPt, binsPt + nBinsPt + 1}; + +// default values for the ML model paths, one model per pT bin +static const std::vector modelPaths = {""}; + +// default values for the cut directions +constexpr int cutDir[nCutScores] = {CutGreater, CutNot}; +auto vecCutDir = std::vector{cutDir, cutDir + nCutScores}; + +// default values for the cuts +constexpr double cuts[nBinsPt][nCutScores] = { + {0.5, 0.5}, + {0.5, 0.5}, + {0.5, 0.5}, + {0.5, 0.5}, + {0.5, 0.5}, + {0.5, 0.5}, + {0.5, 0.5}, + {0.5, 0.5}}; + +// row labels +static const std::vector labelsPt = { + "pT bin 0", + "pT bin 1", + "pT bin 2", + "pT bin 3", + "pT bin 4", + "pT bin 5", + "pT bin 6", + "pT bin 7"}; + +// column labels +static const std::vector labelsCutScore = {"score class 1", "score class 2"}; +static const std::vector labelsDmesCutScore = {"ML score signal", "ML score bkg"}; } // namespace cascade_flow_cuts_ml static constexpr double defaultCutsMl[1][2] = {{0.5, 0.5}}; @@ -233,17 +233,19 @@ struct cascadeFlow { casc.mXi(), casc.mOmega(), v2A, - v2C, - BDTresponse); + v2C, + BDTresponse); } void init(InitContext const&) { float MinMass = 1.28; - if (isOmega) MinMass = 1.6; + if (isOmega) + MinMass = 1.6; float MaxMass = 1.36; - if (isOmega) MaxMass = 1.73; + if (isOmega) + MaxMass = 1.73; ConfigurableAxis vertexZ{"vertexZ", {20, -10, 10}, "vertex axis for histograms"}; ConfigurableAxis massCascAxis{"massCascAxis", {100, MinMass, MaxMass}, ""}; ConfigurableAxis ptAxis{"ptAxis", {100, 0, 10}, ""}; @@ -252,11 +254,11 @@ struct cascadeFlow { histos.add("hCandidate", "hCandidate", HistType::kTH1D, {{22, -0.5, 21.5}}); histos.add("hCascadeSignal", "hCascadeSignal", HistType::kTH1D, {{6, -0.5, 5.5}}); histos.add("hCascade", "hCascade", HistType::kTH1D, {{6, -0.5, 5.5}}); - histos.add("hCascadePhi", "hCascadePhi", HistType::kTH1D, {{100, 0, 2*TMath::Pi()}}); - histos.add("hcascminuspsiT0A", "hcascminuspsiT0A", HistType::kTH1D, {{100, 0, 2*TMath::Pi()}}); - histos.add("hcascminuspsiT0C", "hcascminuspsiT0C", HistType::kTH1D, {{100, 0, 2*TMath::Pi()}}); - histos.add("hPsiT0A", "hPsiT0A", HistType::kTH1D, {{100, 0, 2*TMath::Pi()}}); - histos.add("hPsiT0C", "hPsiT0C", HistType::kTH1D, {{100, 0, 2*TMath::Pi()}}); + histos.add("hCascadePhi", "hCascadePhi", HistType::kTH1D, {{100, 0, 2 * TMath::Pi()}}); + histos.add("hcascminuspsiT0A", "hcascminuspsiT0A", HistType::kTH1D, {{100, 0, 2 * TMath::Pi()}}); + histos.add("hcascminuspsiT0C", "hcascminuspsiT0C", HistType::kTH1D, {{100, 0, 2 * TMath::Pi()}}); + histos.add("hPsiT0A", "hPsiT0A", HistType::kTH1D, {{100, 0, 2 * TMath::Pi()}}); + histos.add("hPsiT0C", "hPsiT0C", HistType::kTH1D, {{100, 0, 2 * TMath::Pi()}}); histos.add("hFT0ARe", "hFT0ARe", HistType::kTH1D, {{100, -1, 1}}); histos.add("hFT0AIm", "hFT0AIm", HistType::kTH1D, {{100, -1, 1}}); histos.add("hFT0A", "hFT0A", HistType::kTH2D, {{100, -1, 1}, {100, -1, 1}}); @@ -270,15 +272,15 @@ struct cascadeFlow { histos.add("hSignalScoreAfterSel", "Signal score after selection;BDT first score;entries", HistType::kTH1F, {{100, 0., 1.}}); histos.add("hBkgScoreAfterSel", "Bkg score after selection;BDT first score;entries", HistType::kTH1F, {{100, 0., 1.}}); - if (isApplyML){ + if (isApplyML) { // Configure and initialise the ML class mlResponse.configure(binsPtMl, cutsMl, cutDirMl, nClassesMl); // Bonus: retrieve the model from CCDB (needed for ML application on the GRID) if (loadModelsFromCCDB) { - ccdbApi.init(ccdbUrl); - mlResponse.setModelPathsCCDB(onnxFileNames, ccdbApi, modelPathsCCDB, timestampCCDB); + ccdbApi.init(ccdbUrl); + mlResponse.setModelPathsCCDB(onnxFileNames, ccdbApi, modelPathsCCDB, timestampCCDB); } else { - mlResponse.setModelPathsLocal(onnxFileNames); + mlResponse.setModelPathsLocal(onnxFileNames); } mlResponse.init(); } @@ -363,7 +365,7 @@ struct cascadeFlow { void processAnalyseData(CollEventPlane const& coll, soa::Join const& Cascades, DauTracks const&) { - //std::cout << "onnxFileNames " << onnxFileNames.value[0] << std::endl; + // std::cout << "onnxFileNames " << onnxFileNames.value[0] << std::endl; if (!coll.sel8() || std::abs(coll.posZ()) > 10.) { return; @@ -378,14 +380,16 @@ struct cascadeFlow { float eventplaneVecT0CNorm = sqrt(eventplaneVecT0C.Dot(eventplaneVecT0C)); // float v2Norm = eventplaneVecT0A.Dot(eventplaneVecT0C)/coll.sumAmplFT0A()/coll.sumAmplFT0C(); float v2Norm = eventplaneVecT0A.Dot(eventplaneVecT0C) / eventplaneVecT0ANorm / eventplaneVecT0CNorm; - //std::cout << "eventplaneVecT0ANorm: " << eventplaneVecT0ANorm << " = " << sqrt(pow(coll.qvecFT0ARe(), 2) + pow(coll.qvecFT0AIm(), 2)) << std::endl; //these two are equal - //std::cout << "stored norm value: " << coll.sumAmplFT0A() << std::endl; + // std::cout << "eventplaneVecT0ANorm: " << eventplaneVecT0ANorm << " = " << sqrt(pow(coll.qvecFT0ARe(), 2) + pow(coll.qvecFT0AIm(), 2)) << std::endl; //these two are equal + // std::cout << "stored norm value: " << coll.sumAmplFT0A() << std::endl; histos.fill(HIST("hv2Norm"), v2Norm); - float PsiT0A = TMath::ACos(coll.qvecFT0ARe()) / 2; + float PsiT0A = TMath::ACos(coll.qvecFT0ARe()) / 2; float PsiT0C = TMath::ACos(coll.qvecFT0CRe()) / 2; - if (coll.qvecFT0AIm() < 0) PsiT0A = -PsiT0A + TMath::Pi()/2; //to get dstribution between 0 and Pi - if (coll.qvecFT0CIm() < 0) PsiT0C = -PsiT0C + TMath::Pi()/2; //to get dstribution between 0 and Pi + if (coll.qvecFT0AIm() < 0) + PsiT0A = -PsiT0A + TMath::Pi() / 2; // to get dstribution between 0 and Pi + if (coll.qvecFT0CIm() < 0) + PsiT0C = -PsiT0C + TMath::Pi() / 2; // to get dstribution between 0 and Pi histos.fill(HIST("hFT0ARe"), coll.qvecFT0ARe()); histos.fill(HIST("hFT0AIm"), coll.qvecFT0AIm()); @@ -408,77 +412,82 @@ struct cascadeFlow { bool isSelectedCasc = false; std::vector inputFeaturesCasc{casc.cascradius(), - casc.v0radius(), - casc.casccosPA(coll.posX(), coll.posY(), coll.posZ()), - casc.v0cosPA(coll.posX(), coll.posY(), coll.posZ()), - casc.dcapostopv(), - casc.dcanegtopv(), - casc.dcabachtopv(), - casc.dcacascdaughters(), - casc.dcaV0daughters(), - casc.dcav0topv(coll.posX(), coll.posY(), coll.posZ()), - casc.bachBaryonCosPA(), - casc.bachBaryonDCAxyToPV()}; + casc.v0radius(), + casc.casccosPA(coll.posX(), coll.posY(), coll.posZ()), + casc.v0cosPA(coll.posX(), coll.posY(), coll.posZ()), + casc.dcapostopv(), + casc.dcanegtopv(), + casc.dcabachtopv(), + casc.dcacascdaughters(), + casc.dcaV0daughters(), + casc.dcav0topv(coll.posX(), coll.posY(), coll.posZ()), + casc.bachBaryonCosPA(), + casc.bachBaryonDCAxyToPV()}; float massCasc = casc.mXi(); - if (isOmega) massCasc = casc.mOmega(); + if (isOmega) + massCasc = casc.mOmega(); - //inv mass loose cut - if (casc.pt() < MinPt) continue; - if (casc.pt() > MaxPt) continue; + // inv mass loose cut + if (casc.pt() < MinPt) + continue; + if (casc.pt() > MaxPt) + continue; if (isOmega) { - if (massCasc < 1.6 || massCasc > 1.73) continue; + if (massCasc < 1.6 || massCasc > 1.73) + continue; } else { - if (massCasc < 1.28 || massCasc > 1.36) continue; + if (massCasc < 1.28 || massCasc > 1.36) + continue; } histos.fill(HIST("hMassBeforeSel"), massCasc); histos.fill(HIST("hMassBeforeSelVsPt"), massCasc, casc.pt()); - if (isApplyML){ - // Retrieve model output and selection outcome - isSelectedCasc = mlResponse.isSelectedMl(inputFeaturesCasc, casc.pt(), outputMl); - // Fill BDT score histograms before selection - histos.fill(HIST("hSignalScoreBeforeSel"), outputMl[0]); - histos.fill(HIST("hBkgScoreBeforeSel"), outputMl[1]); - - // Fill histograms for selected candidates - if (isSelectedCasc) { - histos.fill(HIST("hSignalScoreAfterSel"), outputMl[0]); - histos.fill(HIST("hBkgScoreAfterSel"), outputMl[1]); - } + if (isApplyML) { + // Retrieve model output and selection outcome + isSelectedCasc = mlResponse.isSelectedMl(inputFeaturesCasc, casc.pt(), outputMl); + // Fill BDT score histograms before selection + histos.fill(HIST("hSignalScoreBeforeSel"), outputMl[0]); + histos.fill(HIST("hBkgScoreBeforeSel"), outputMl[1]); + + // Fill histograms for selected candidates + if (isSelectedCasc) { + histos.fill(HIST("hSignalScoreAfterSel"), outputMl[0]); + histos.fill(HIST("hBkgScoreAfterSel"), outputMl[1]); + } - outputMl.clear(); // not necessary in this case but for good measure + outputMl.clear(); // not necessary in this case but for good measure } else { - isSelectedCasc = true; + isSelectedCasc = true; } if (isSelectedCasc) { - histos.fill(HIST("hMassAfterSel"), massCasc); - histos.fill(HIST("hMassAfterSelVsPt"), massCasc, casc.pt()); + histos.fill(HIST("hMassAfterSel"), massCasc); + histos.fill(HIST("hMassAfterSelVsPt"), massCasc, casc.pt()); } - cascphiVec = ROOT::Math::XYZVector(cos(2*casc.phi()), sin(2*casc.phi()), 0); - auto v2A = cascphiVec.Dot(eventplaneVecT0A)/sqrt(eventplaneVecT0A.Dot(eventplaneVecT0A)); - auto v2C = cascphiVec.Dot(eventplaneVecT0C)/sqrt(eventplaneVecT0C.Dot(eventplaneVecT0C)); + cascphiVec = ROOT::Math::XYZVector(cos(2 * casc.phi()), sin(2 * casc.phi()), 0); + auto v2A = cascphiVec.Dot(eventplaneVecT0A) / sqrt(eventplaneVecT0A.Dot(eventplaneVecT0A)); + auto v2C = cascphiVec.Dot(eventplaneVecT0C) / sqrt(eventplaneVecT0C.Dot(eventplaneVecT0C)); auto cascminuspsiT0A = GetPhiInRange(casc.phi() - PsiT0A); auto cascminuspsiT0C = GetPhiInRange(casc.phi() - PsiT0C); - //float v2A_diff = TMath::Cos(2.0 * cascminuspsiT0A); - //float v2C_diff = TMath::Cos(2.0 * cascminuspsiT0C); - //std::cout << "v2A: " << v2A << " v2C " << v2C << std::endl; - //std::cout << "v2A_diff: " << v2A_diff << " v2C_diff " << v2C_diff << std::endl; + // float v2A_diff = TMath::Cos(2.0 * cascminuspsiT0A); + // float v2C_diff = TMath::Cos(2.0 * cascminuspsiT0C); + // std::cout << "v2A: " << v2A << " v2C " << v2C << std::endl; + // std::cout << "v2A_diff: " << v2A_diff << " v2C_diff " << v2C_diff << std::endl; histos.fill(HIST("hCascadePhi"), casc.phi()); histos.fill(HIST("hcascminuspsiT0A"), cascminuspsiT0A); histos.fill(HIST("hcascminuspsiT0C"), cascminuspsiT0C); float BDTresponse = 0; - if (isApplyML) BDTresponse = outputMl[0]; - if (isSelectedCasc) fillAnalysedTable(coll, casc, v2A, v2C, BDTresponse); + if (isApplyML) + BDTresponse = outputMl[0]; + if (isSelectedCasc) + fillAnalysedTable(coll, casc, v2A, v2C, BDTresponse); } - } - PROCESS_SWITCH(cascadeFlow, processTrainingBackground, "Process to create the training dataset for the background", true); PROCESS_SWITCH(cascadeFlow, processTrainingSignal, "Process to create the training dataset for the signal", false); PROCESS_SWITCH(cascadeFlow, processAnalyseData, "Process to apply ML model to the data", false); From b7976c7d3f28042272ec39c02c6243e7bdaf4307 Mon Sep 17 00:00:00 2001 From: Chiara De Martin Date: Wed, 13 Mar 2024 18:27:01 +0100 Subject: [PATCH 10/11] remove unused variable --- PWGLF/TableProducer/cascadeflow.cxx | 1 - 1 file changed, 1 deletion(-) diff --git a/PWGLF/TableProducer/cascadeflow.cxx b/PWGLF/TableProducer/cascadeflow.cxx index 6e22b06b5cd..f3ff85f26f9 100644 --- a/PWGLF/TableProducer/cascadeflow.cxx +++ b/PWGLF/TableProducer/cascadeflow.cxx @@ -110,7 +110,6 @@ namespace cascade_flow_cuts_ml static const std::vector labelsDmesCutScore = {"ML score signal", "ML score bkg"}; } // namespace cascade_flow_cuts_ml -static constexpr double defaultCutsMl[1][2] = {{0.5, 0.5}}; struct cascadeFlow { From eaa30cfe04996f5e4a1a0e1343e737f37103bf88 Mon Sep 17 00:00:00 2001 From: ALICE Action Bot Date: Wed, 13 Mar 2024 17:28:22 +0000 Subject: [PATCH 11/11] Please consider the following formatting changes --- PWGLF/TableProducer/cascadeflow.cxx | 1 - 1 file changed, 1 deletion(-) diff --git a/PWGLF/TableProducer/cascadeflow.cxx b/PWGLF/TableProducer/cascadeflow.cxx index 075313899c9..993edee88c4 100644 --- a/PWGLF/TableProducer/cascadeflow.cxx +++ b/PWGLF/TableProducer/cascadeflow.cxx @@ -110,7 +110,6 @@ static const std::vector labelsCutScore = {"score class 1", "score static const std::vector labelsDmesCutScore = {"ML score signal", "ML score bkg"}; } // namespace cascade_flow_cuts_ml - struct cascadeFlow { Configurable isOmega{"isOmega", 0, "Xi or Omega"};