|
| 1 | +// Copyright 2019-2020 CERN and copyright holders of ALICE O2. |
| 2 | +// See https://alice-o2.web.cern.ch/copyright for details of the copyright holders. |
| 3 | +// All rights not expressly granted are reserved. |
| 4 | +// |
| 5 | +// This software is distributed under the terms of the GNU General Public |
| 6 | +// License v3 (GPL Version 3), copied verbatim in the file "COPYING". |
| 7 | +// |
| 8 | +// In applying this license CERN does not waive the privileges and immunities |
| 9 | +// granted to it by virtue of its status as an Intergovernmental Organization |
| 10 | +// or submit itself to any jurisdiction. |
| 11 | +/// |
| 12 | +/// \brief Task to create derived data for cascade flow analyses |
| 13 | +/// \authors: Chiara De Martin (chiara.de.martin@cern.ch), Maximiliano Puccio (maximiliano.puccio@cern.ch) |
| 14 | + |
| 15 | +#include "Common/DataModel/Centrality.h" |
| 16 | +#include "Common/DataModel/EventSelection.h" |
| 17 | +#include "Common/DataModel/Multiplicity.h" |
| 18 | +#include "Common/DataModel/PIDResponse.h" |
| 19 | +#include "Common/DataModel/TrackSelectionTables.h" |
| 20 | +#include "Framework/AnalysisTask.h" |
| 21 | +#include "Framework/O2DatabasePDGPlugin.h" |
| 22 | +#include "Framework/runDataProcessing.h" |
| 23 | +#include "PWGLF/DataModel/cascqaanalysis.h" |
| 24 | +#include "PWGLF/DataModel/LFStrangenessTables.h" |
| 25 | +#include "PWGLF/DataModel/LFStrangenessPIDTables.h" |
| 26 | +#include "TRandom3.h" |
| 27 | +#include "Framework/ASoAHelpers.h" |
| 28 | + |
| 29 | +using namespace o2; |
| 30 | +using namespace o2::framework; |
| 31 | +using namespace o2::framework::expressions; |
| 32 | +using std::array; |
| 33 | + |
| 34 | +using DauTracks = soa::Join<aod::DauTrackExtras, aod::DauTrackTPCPIDs>; |
| 35 | + |
| 36 | +namespace cascadev2 |
| 37 | +{ |
| 38 | +enum species { Xi = 0, |
| 39 | + Omega = 1 }; |
| 40 | +constexpr double massSigmaParameters[4][2]{ |
| 41 | + {4.9736e-3, 0.006815}, |
| 42 | + {-2.39594, -2.257}, |
| 43 | + {1.8064e-3, 0.00138}, |
| 44 | + {1.03468e-1, 0.1898}}; |
| 45 | +static const std::vector<std::string> massSigmaParameterNames{"p0", "p1", "p2", "p3"}; |
| 46 | +static const std::vector<std::string> speciesNames{"Xi", "Omega"}; |
| 47 | +} // namespace cascadev2 |
| 48 | + |
| 49 | +struct cascadeFlow { |
| 50 | + |
| 51 | + Configurable<double> sideBandStart{"sideBandStart", 5, "Start of the sideband region in number of sigmas"}; |
| 52 | + Configurable<double> sideBandEnd{"sideBandEnd", 7, "End of the sideband region in number of sigmas"}; |
| 53 | + Configurable<double> downsample{"downsample", 1., "Downsample training output tree"}; |
| 54 | + Configurable<bool> doNTPCSigmaCut{"doNTPCSigmaCut", 1, "doNtpcSigmaCut"}; |
| 55 | + Configurable<float> nsigmatpcPr{"nsigmatpcPr", 5, "nsigmatpcPr"}; |
| 56 | + Configurable<float> nsigmatpcPi{"nsigmatpcPi", 5, "nsigmatpcPi"}; |
| 57 | + Configurable<float> mintpccrrows{"mintpccrrows", 3, "mintpccrrows"}; |
| 58 | + |
| 59 | + template <typename TCascade, typename TDaughter> |
| 60 | + bool IsCascAccepted(TCascade casc, TDaughter negExtra, TDaughter posExtra, TDaughter bachExtra, int& counter) // loose cuts on topological selections of cascades |
| 61 | + { |
| 62 | + // TPC cuts as those implemented for the training of the signal |
| 63 | + if (doNTPCSigmaCut) { |
| 64 | + if (casc.sign() < 0) { |
| 65 | + if (std::abs(posExtra.tpcNSigmaPr()) > nsigmatpcPr || std::abs(negExtra.tpcNSigmaPi()) > nsigmatpcPi) |
| 66 | + return false; |
| 67 | + } else if (casc.sign() > 0) { |
| 68 | + if (std::abs(posExtra.tpcNSigmaPi()) > nsigmatpcPi || std::abs(negExtra.tpcNSigmaPr()) > nsigmatpcPr) |
| 69 | + return false; |
| 70 | + } |
| 71 | + } |
| 72 | + counter++; |
| 73 | + |
| 74 | + if (posExtra.tpcCrossedRows() < mintpccrrows || negExtra.tpcCrossedRows() < mintpccrrows || bachExtra.tpcCrossedRows() < mintpccrrows) |
| 75 | + return false; |
| 76 | + |
| 77 | + counter++; |
| 78 | + return true; |
| 79 | + } |
| 80 | + |
| 81 | + HistogramRegistry histos{"histos", {}, OutputObjHandlingPolicy::AnalysisObject}; |
| 82 | + |
| 83 | + // Tables to produce |
| 84 | + Produces<aod::CascTraining> trainingSample; |
| 85 | + Configurable<LabeledArray<double>> parSigmaMass{ |
| 86 | + "parSigmaMass", |
| 87 | + {cascadev2::massSigmaParameters[0], 4, 2, |
| 88 | + cascadev2::massSigmaParameterNames, cascadev2::speciesNames}, |
| 89 | + "Mass resolution parameters: [0]*exp([1]*x)+[2]*exp([3]*x)"}; |
| 90 | + |
| 91 | + float getNsigmaMass(const cascadev2::species s, const float pt, const float nsigma = 6) |
| 92 | + { |
| 93 | + const auto sigma = parSigmaMass->get(0u, s) * exp(parSigmaMass->get(1, s) * pt) + parSigmaMass->get(2, s) * exp(parSigmaMass->get(3, s) * pt); |
| 94 | + return nsigma * sigma; |
| 95 | + } |
| 96 | + |
| 97 | + template <class collision_t, class cascade_t> |
| 98 | + void fillTrainingTable(collision_t coll, cascade_t casc, int pdgCode) |
| 99 | + { |
| 100 | + trainingSample(coll.centFT0M(), |
| 101 | + casc.sign(), |
| 102 | + casc.pt(), |
| 103 | + casc.eta(), |
| 104 | + casc.mXi(), |
| 105 | + casc.mOmega(), |
| 106 | + casc.mLambda(), |
| 107 | + casc.cascradius(), |
| 108 | + casc.v0radius(), |
| 109 | + casc.casccosPA(coll.posX(), coll.posY(), coll.posZ()), |
| 110 | + casc.v0cosPA(coll.posX(), coll.posY(), coll.posZ()), |
| 111 | + casc.dcapostopv(), |
| 112 | + casc.dcanegtopv(), |
| 113 | + casc.dcabachtopv(), |
| 114 | + casc.dcacascdaughters(), |
| 115 | + casc.dcaV0daughters(), |
| 116 | + casc.dcav0topv(coll.posX(), coll.posY(), coll.posZ()), |
| 117 | + casc.bachBaryonCosPA(), |
| 118 | + casc.bachBaryonDCAxyToPV(), |
| 119 | + pdgCode); |
| 120 | + } |
| 121 | + |
| 122 | + void init(InitContext const&) |
| 123 | + { |
| 124 | + ConfigurableAxis vertexZ{"vertexZ", {20, -10, 10}, "vertex axis for histograms"}; |
| 125 | + histos.add("hEventVertexZ", "hEventVertexZ", kTH1F, {vertexZ}); |
| 126 | + histos.add("hEventCentrality", "hEventCentrality", kTH1F, {{101, 0, 101}}); |
| 127 | + histos.add("hEventSelection", "hEventSelection", kTH1F, {{4, 0, 4}}); |
| 128 | + histos.add("hCandidate", "hCandidate", HistType::kTH1D, {{22, -0.5, 21.5}}); |
| 129 | + histos.add("hCascadeSignal", "hCascadeSignal", HistType::kTH1D, {{6, -0.5, 5.5}}); |
| 130 | + } |
| 131 | + |
| 132 | + void processTrainingBackground(soa::Join<aod::StraCollisions, aod::StraCents, aod::StraEvSels>::iterator const& coll, soa::Join<aod::CascCollRefs, aod::CascCores, aod::CascExtras, aod::CascBBs> const& Cascades, DauTracks const&) |
| 133 | + { |
| 134 | + |
| 135 | + int counter = 0; |
| 136 | + |
| 137 | + if (!coll.sel8() || std::abs(coll.posZ()) > 10.) { |
| 138 | + return; |
| 139 | + } |
| 140 | + |
| 141 | + for (auto& casc : Cascades) { |
| 142 | + if (gRandom->Uniform() > downsample) { |
| 143 | + continue; |
| 144 | + } |
| 145 | + |
| 146 | + float sigmaRangeXi[2]{getNsigmaMass(cascadev2::Xi, casc.pt(), sideBandStart), getNsigmaMass(cascadev2::Xi, casc.pt(), sideBandEnd)}; |
| 147 | + float sigmaRangeOmega[2]{getNsigmaMass(cascadev2::Omega, casc.pt(), sideBandStart), getNsigmaMass(cascadev2::Omega, casc.pt(), sideBandEnd)}; |
| 148 | + |
| 149 | + if ((std::abs(casc.mXi() - constants::physics::MassXiMinus) < sigmaRangeXi[0] || |
| 150 | + std::abs(casc.mXi() - constants::physics::MassXiMinus) > sigmaRangeXi[1]) && |
| 151 | + (std::abs(casc.mOmega() - constants::physics::MassOmegaMinus) < sigmaRangeOmega[0] || |
| 152 | + std::abs(casc.mOmega() - constants::physics::MassOmegaMinus) > sigmaRangeOmega[1])) { |
| 153 | + continue; |
| 154 | + } |
| 155 | + |
| 156 | + /// Add some minimal cuts for single track variables (min number of TPC clusters) |
| 157 | + auto negExtra = casc.negTrackExtra_as<soa::Join<aod::DauTrackExtras, aod::DauTrackTPCPIDs>>(); |
| 158 | + auto posExtra = casc.posTrackExtra_as<soa::Join<aod::DauTrackExtras, aod::DauTrackTPCPIDs>>(); |
| 159 | + auto bachExtra = casc.bachTrackExtra_as<soa::Join<aod::DauTrackExtras, aod::DauTrackTPCPIDs>>(); |
| 160 | + |
| 161 | + if (doNTPCSigmaCut) { |
| 162 | + if (casc.sign() < 0) { |
| 163 | + if (std::abs(posExtra.tpcNSigmaPr()) > nsigmatpcPr || std::abs(negExtra.tpcNSigmaPi()) > nsigmatpcPi) |
| 164 | + continue; |
| 165 | + histos.fill(HIST("hCandidate"), ++counter); |
| 166 | + } else if (casc.sign() > 0) { |
| 167 | + if (std::abs(posExtra.tpcNSigmaPi()) > nsigmatpcPi || std::abs(negExtra.tpcNSigmaPr()) > nsigmatpcPr) |
| 168 | + continue; |
| 169 | + histos.fill(HIST("hCandidate"), ++counter); |
| 170 | + } |
| 171 | + } else |
| 172 | + ++counter; |
| 173 | + |
| 174 | + if (posExtra.tpcCrossedRows() < mintpccrrows || negExtra.tpcCrossedRows() < mintpccrrows || bachExtra.tpcCrossedRows() < mintpccrrows) |
| 175 | + continue; |
| 176 | + histos.fill(HIST("hCandidate"), ++counter); |
| 177 | + |
| 178 | + fillTrainingTable(coll, casc, 0); |
| 179 | + } |
| 180 | + } |
| 181 | + |
| 182 | + void processTrainingSignal(soa::Join<aod::StraCollisions, aod::StraCents, aod::StraEvSels>::iterator const& coll, soa::Join<aod::CascCollRefs, aod::CascCores, aod::CascMCCores, aod::CascExtras, aod::CascBBs> const& Cascades, soa::Join<aod::DauTrackExtras, aod::DauTrackTPCPIDs> const&) |
| 183 | + { |
| 184 | + |
| 185 | + if (!coll.sel8() || std::abs(coll.posZ()) > 10.) { |
| 186 | + return; |
| 187 | + } |
| 188 | + |
| 189 | + for (auto& casc : Cascades) { |
| 190 | + int pdgCode{casc.pdgCode()}; |
| 191 | + if (!(std::abs(pdgCode) == 3312 && std::abs(casc.pdgCodeV0()) == 3122 && std::abs(casc.pdgCodeBachelor() == 211)) // Xi |
| 192 | + && !(std::abs(pdgCode) == 3334 && std::abs(casc.pdgCodeV0()) == 3122 && std::abs(casc.pdgCodeBachelor() == 321))) // Omega |
| 193 | + continue; |
| 194 | + |
| 195 | + auto negExtra = casc.negTrackExtra_as<soa::Join<aod::DauTrackExtras, aod::DauTrackTPCPIDs>>(); |
| 196 | + auto posExtra = casc.posTrackExtra_as<soa::Join<aod::DauTrackExtras, aod::DauTrackTPCPIDs>>(); |
| 197 | + auto bachExtra = casc.bachTrackExtra_as<soa::Join<aod::DauTrackExtras, aod::DauTrackTPCPIDs>>(); |
| 198 | + |
| 199 | + int counter = 0; |
| 200 | + IsCascAccepted(casc, negExtra, posExtra, bachExtra, counter); |
| 201 | + histos.fill(HIST("hCascadeSignal"), counter); |
| 202 | + |
| 203 | + // PDG cascades |
| 204 | + fillTrainingTable(coll, casc, pdgCode); |
| 205 | + } |
| 206 | + } |
| 207 | + |
| 208 | + PROCESS_SWITCH(cascadeFlow, processTrainingBackground, "Process to create the training dataset for the background", true); |
| 209 | + PROCESS_SWITCH(cascadeFlow, processTrainingSignal, "Process to create the training dataset for the signal", false); |
| 210 | +}; |
| 211 | + |
| 212 | +WorkflowSpec defineDataProcessing(ConfigContext const& cfgc) |
| 213 | +{ |
| 214 | + return WorkflowSpec{ |
| 215 | + adaptAnalysisTask<cascadeFlow>(cfgc, TaskName{"lf-cascade-flow"})}; |
| 216 | +} |
0 commit comments