diff --git a/PWGLF/DataModel/cascqaanalysis.h b/PWGLF/DataModel/cascqaanalysis.h index c5377b201fa..edbf9a42c32 100644 --- a/PWGLF/DataModel/cascqaanalysis.h +++ b/PWGLF/DataModel/cascqaanalysis.h @@ -129,6 +129,10 @@ DECLARE_SOA_TABLE(MyCascades, "AOD", "MYCASCADES", o2::soa::Index<>, mycascades::IsINELgt0, mycascades::IsINELgt1); +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); + namespace myMCcascades { diff --git a/PWGLF/TableProducer/CMakeLists.txt b/PWGLF/TableProducer/CMakeLists.txt index d3a58259240..0074e9b6b05 100644 --- a/PWGLF/TableProducer/CMakeLists.txt +++ b/PWGLF/TableProducer/CMakeLists.txt @@ -172,6 +172,11 @@ o2physics_add_dpl_workflow(cascqaanalysis PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore COMPONENT_NAME Analysis) +o2physics_add_dpl_workflow(cascadeflow + SOURCES cascadeflow.cxx + PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore + COMPONENT_NAME Analysis) + # Spectra o2physics_add_dpl_workflow(spectra-derived SOURCES spectraDerivedMaker.cxx diff --git a/PWGLF/TableProducer/cascadeflow.cxx b/PWGLF/TableProducer/cascadeflow.cxx new file mode 100644 index 00000000000..a9a16880adc --- /dev/null +++ b/PWGLF/TableProducer/cascadeflow.cxx @@ -0,0 +1,216 @@ +// 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. +/// +/// \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 "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/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" + +using namespace o2; +using namespace o2::framework; +using namespace o2::framework::expressions; +using std::array; + +using DauTracks = soa::Join; + +namespace cascadev2 +{ +enum species { Xi = 0, + Omega = 1 }; +constexpr double massSigmaParameters[4][2]{ + {4.9736e-3, 0.006815}, + {-2.39594, -2.257}, + {1.8064e-3, 0.00138}, + {1.03468e-1, 0.1898}}; +static const std::vector massSigmaParameterNames{"p0", "p1", "p2", "p3"}; +static const std::vector speciesNames{"Xi", "Omega"}; +} // namespace cascadev2 + +struct cascadeFlow { + + 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"}; + Configurable doNTPCSigmaCut{"doNTPCSigmaCut", 1, "doNtpcSigmaCut"}; + Configurable nsigmatpcPr{"nsigmatpcPr", 5, "nsigmatpcPr"}; + Configurable nsigmatpcPi{"nsigmatpcPi", 5, "nsigmatpcPi"}; + Configurable mintpccrrows{"mintpccrrows", 3, "mintpccrrows"}; + + template + bool IsCascAccepted(TCascade casc, TDaughter negExtra, TDaughter posExtra, TDaughter bachExtra, int& counter) // loose cuts on topological selections of cascades + { + // TPC cuts as those implemented for the training of the signal + if (doNTPCSigmaCut) { + if (casc.sign() < 0) { + if (std::abs(posExtra.tpcNSigmaPr()) > nsigmatpcPr || std::abs(negExtra.tpcNSigmaPi()) > nsigmatpcPi) + return false; + } else if (casc.sign() > 0) { + if (std::abs(posExtra.tpcNSigmaPi()) > nsigmatpcPi || std::abs(negExtra.tpcNSigmaPr()) > nsigmatpcPr) + return false; + } + } + counter++; + + if (posExtra.tpcCrossedRows() < mintpccrrows || negExtra.tpcCrossedRows() < mintpccrrows || bachExtra.tpcCrossedRows() < mintpccrrows) + return false; + + counter++; + return true; + } + + HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject}; + + // Tables to produce + Produces trainingSample; + Configurable> parSigmaMass{ + "parSigmaMass", + {cascadev2::massSigmaParameters[0], 4, 2, + cascadev2::massSigmaParameterNames, cascadev2::speciesNames}, + "Mass resolution parameters: [0]*exp([1]*x)+[2]*exp([3]*x)"}; + + float getNsigmaMass(const cascadev2::species s, const float pt, const float nsigma = 6) + { + const auto sigma = parSigmaMass->get(0u, s) * exp(parSigmaMass->get(1, s) * pt) + parSigmaMass->get(2, s) * exp(parSigmaMass->get(3, s) * pt); + return nsigma * sigma; + } + + template + void fillTrainingTable(collision_t coll, cascade_t casc, int pdgCode) + { + trainingSample(coll.centFT0M(), + casc.sign(), + casc.pt(), + casc.eta(), + casc.mXi(), + casc.mOmega(), + casc.mLambda(), + 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(), + pdgCode); + } + + 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("hEventSelection", "hEventSelection", kTH1F, {{4, 0, 4}}); + histos.add("hCandidate", "hCandidate", HistType::kTH1D, {{22, -0.5, 21.5}}); + histos.add("hCascadeSignal", "hCascadeSignal", HistType::kTH1D, {{6, -0.5, 5.5}}); + } + + void processTrainingBackground(soa::Join::iterator const& coll, soa::Join const& Cascades, DauTracks const&) + { + + int counter = 0; + + if (!coll.sel8() || std::abs(coll.posZ()) > 10.) { + return; + } + + for (auto& casc : Cascades) { + if (gRandom->Uniform() > downsample) { + continue; + } + + float sigmaRangeXi[2]{getNsigmaMass(cascadev2::Xi, casc.pt(), sideBandStart), getNsigmaMass(cascadev2::Xi, casc.pt(), sideBandEnd)}; + float sigmaRangeOmega[2]{getNsigmaMass(cascadev2::Omega, casc.pt(), sideBandStart), getNsigmaMass(cascadev2::Omega, casc.pt(), sideBandEnd)}; + + if ((std::abs(casc.mXi() - constants::physics::MassXiMinus) < sigmaRangeXi[0] || + std::abs(casc.mXi() - constants::physics::MassXiMinus) > sigmaRangeXi[1]) && + (std::abs(casc.mOmega() - constants::physics::MassOmegaMinus) < sigmaRangeOmega[0] || + std::abs(casc.mOmega() - constants::physics::MassOmegaMinus) > sigmaRangeOmega[1])) { + continue; + } + + /// 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>(); + + if (doNTPCSigmaCut) { + if (casc.sign() < 0) { + if (std::abs(posExtra.tpcNSigmaPr()) > nsigmatpcPr || std::abs(negExtra.tpcNSigmaPi()) > nsigmatpcPi) + continue; + histos.fill(HIST("hCandidate"), ++counter); + } else if (casc.sign() > 0) { + if (std::abs(posExtra.tpcNSigmaPi()) > nsigmatpcPi || std::abs(negExtra.tpcNSigmaPr()) > nsigmatpcPr) + continue; + histos.fill(HIST("hCandidate"), ++counter); + } + } else + ++counter; + + if (posExtra.tpcCrossedRows() < mintpccrrows || negExtra.tpcCrossedRows() < mintpccrrows || bachExtra.tpcCrossedRows() < mintpccrrows) + continue; + histos.fill(HIST("hCandidate"), ++counter); + + fillTrainingTable(coll, casc, 0); + } + } + + void processTrainingSignal(soa::Join::iterator const& coll, soa::Join const& Cascades, soa::Join const&) + { + + if (!coll.sel8() || std::abs(coll.posZ()) > 10.) { + return; + } + + for (auto& casc : Cascades) { + int pdgCode{casc.pdgCode()}; + if (!(std::abs(pdgCode) == 3312 && std::abs(casc.pdgCodeV0()) == 3122 && std::abs(casc.pdgCodeBachelor() == 211)) // Xi + && !(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>(); + + int counter = 0; + IsCascAccepted(casc, negExtra, posExtra, bachExtra, counter); + histos.fill(HIST("hCascadeSignal"), counter); + + // PDG cascades + fillTrainingTable(coll, casc, pdgCode); + } + } + + 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); +}; + +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) +{ + return WorkflowSpec{ + adaptAnalysisTask(cfgc, TaskName{"lf-cascade-flow"})}; +}