First version of cascadeflow

PWGLF/DataModel/cascqaanalysis.h

@@ -129,6 +129,10 @@ DECLARE_SOA_TABLE(MyCascades, "AOD", "MYCASCADES", o2::soa::Index<>,
                   mycascades::IsINELgt0<mycascades::EventSelFilterBitMask>,
                   mycascades::IsINELgt1<mycascades::EventSelFilterBitMask>);
 
+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
 {
 

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

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<aod::DauTrackExtras, aod::DauTrackTPCPIDs>;
+
+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<std::string> massSigmaParameterNames{"p0", "p1", "p2", "p3"};
+static const std::vector<std::string> speciesNames{"Xi", "Omega"};
+} // namespace cascadev2
+
+struct cascadeFlow {
+
+  Configurable<double> sideBandStart{"sideBandStart", 5, "Start of the sideband region in number of sigmas"};
+  Configurable<double> sideBandEnd{"sideBandEnd", 7, "End of the sideband region in number of sigmas"};
+  Configurable<double> downsample{"downsample", 1., "Downsample training output tree"};
+  Configurable<bool> doNTPCSigmaCut{"doNTPCSigmaCut", 1, "doNtpcSigmaCut"};
+  Configurable<float> nsigmatpcPr{"nsigmatpcPr", 5, "nsigmatpcPr"};
+  Configurable<float> nsigmatpcPi{"nsigmatpcPi", 5, "nsigmatpcPi"};
+  Configurable<float> mintpccrrows{"mintpccrrows", 3, "mintpccrrows"};
+
+  template <typename TCascade, typename TDaughter>
+  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<aod::CascTraining> trainingSample;
+  Configurable<LabeledArray<double>> 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 <class collision_t, class cascade_t>
+  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<aod::StraCollisions, aod::StraCents, aod::StraEvSels>::iterator const& coll, soa::Join<aod::CascCollRefs, aod::CascCores, aod::CascExtras, aod::CascBBs> 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<soa::Join<aod::DauTrackExtras, aod::DauTrackTPCPIDs>>();
+      auto posExtra = casc.posTrackExtra_as<soa::Join<aod::DauTrackExtras, aod::DauTrackTPCPIDs>>();
+      auto bachExtra = casc.bachTrackExtra_as<soa::Join<aod::DauTrackExtras, aod::DauTrackTPCPIDs>>();
+
+      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<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&)
+  {
+
+    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<soa::Join<aod::DauTrackExtras, aod::DauTrackTPCPIDs>>();
+      auto posExtra = casc.posTrackExtra_as<soa::Join<aod::DauTrackExtras, aod::DauTrackTPCPIDs>>();
+      auto bachExtra = casc.bachTrackExtra_as<soa::Join<aod::DauTrackExtras, aod::DauTrackTPCPIDs>>();
+
+      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<cascadeFlow>(cfgc, TaskName{"lf-cascade-flow"})};
+}