[PWGLF] Fix derived data production + adapt analysis task (#7373)

* Fix derived data production + adapt analysis task

* Fix duplicate entries + rework the processGenerated logic + add histo for event loss estimation

* Fix (unsorted) preslice + fix double MC particle association

* Enable all event selections in MC + add MC information on mother particle

* Add histograms for event loss estimation

* Add option to access casc. mom. info at prim. vtx

* Remove obsolote event selections

* Remove debugging "printf"

* Please consider the following formating changes

* Please consider the following formatting changes

* Please consider the following formatting changes

* Please consider the following formatting changes

* Please consider the following warning fixes

* Please consider the following warning fixes

Author: romainschotter

PWGLF/DataModel/LFStrangenessTables.h

@@ -20,6 +20,7 @@
 #include "Common/DataModel/Multiplicity.h"
 #include "Common/DataModel/Centrality.h"
 #include "Common/DataModel/Qvectors.h"
+#include "Common/DataModel/McCollisionExtra.h"
 #include "PWGLF/DataModel/EPCalibrationTables.h"
 
 namespace o2::aod
@@ -294,6 +295,9 @@ DECLARE_SOA_COLUMN(PzPosMC, pzPosMC, float);                    //! V0 positive
 DECLARE_SOA_COLUMN(PxNegMC, pxNegMC, float);                    //! V0 positive daughter px (GeV/c)
 DECLARE_SOA_COLUMN(PyNegMC, pyNegMC, float);                    //! V0 positive daughter py (GeV/c)
 DECLARE_SOA_COLUMN(PzNegMC, pzNegMC, float);                    //! V0 positive daughter pz (GeV/c)
+DECLARE_SOA_COLUMN(PxMC, pxMC, float);                          //! V0 px (GeV/c)
+DECLARE_SOA_COLUMN(PyMC, pyMC, float);                          //! V0 py (GeV/c)
+DECLARE_SOA_COLUMN(PzMC, pzMC, float);                          //! V0 pz (GeV/c)
 
 //______________________________________________________
 // Binned content for generated particles: derived data
@@ -464,6 +468,22 @@ DECLARE_SOA_DYNAMIC_COLUMN(IsPhotonTPConly, isPhotonTPConly, //! is tpc-only pho
                            [](uint8_t V0Type) -> bool { return V0Type & (1 << 1); });
 DECLARE_SOA_DYNAMIC_COLUMN(IsCollinear, isCollinear, //! is collinear V0
                            [](uint8_t V0Type) -> bool { return V0Type & (1 << 2); });
+
+DECLARE_SOA_DYNAMIC_COLUMN(RapidityMC, rapidityMC, //! rapidity (0:K0, 1:L, 2:Lbar)
+                           [](float PxMC, float PyMC, float PzMC, int value) -> float {
+                             if (value == 0)
+                               return RecoDecay::y(std::array{PxMC, PyMC, PzMC}, o2::constants::physics::MassKaonNeutral);
+                             if (value == 1 || value == 2)
+                               return RecoDecay::y(std::array{PxMC, PyMC, PzMC}, o2::constants::physics::MassLambda);
+                             return 0.0f;
+                           });
+
+DECLARE_SOA_DYNAMIC_COLUMN(NegativePtMC, negativeptMC, //! negative daughter pT
+                           [](float pxnegMC, float pynegMC) -> float { return RecoDecay::sqrtSumOfSquares(pxnegMC, pynegMC); });
+DECLARE_SOA_DYNAMIC_COLUMN(PositivePtMC, positiveptMC, //! positive daughter pT
+                           [](float pxposMC, float pyposMC) -> float { return RecoDecay::sqrtSumOfSquares(pxposMC, pyposMC); });
+DECLARE_SOA_DYNAMIC_COLUMN(PtMC, ptMC, //! V0 pT
+                           [](float pxMC, float pyMC) -> float { return RecoDecay::sqrtSumOfSquares(pxMC, pyMC); });
 } // namespace v0data
 
 DECLARE_SOA_TABLE(V0Indices, "AOD", "V0INDEX", //! index table when using AO2Ds
@@ -666,6 +686,19 @@ DECLARE_SOA_TABLE_VERSIONED(V0MCCores_001, "AOD", "V0MCCORE", 1, //! debug infor
                             v0data::PxPosMC, v0data::PyPosMC, v0data::PzPosMC,
                             v0data::PxNegMC, v0data::PyNegMC, v0data::PzNegMC);
 
+DECLARE_SOA_TABLE_VERSIONED(V0MCCores_002, "AOD", "V0MCCORE", 2, //! debug information
+                            v0data::ParticleIdMC,                //! MC properties of the V0 for posterior analysis
+                            v0data::PDGCode, v0data::PDGCodeMother,
+                            v0data::PDGCodePositive, v0data::PDGCodeNegative,
+                            v0data::IsPhysicalPrimary, v0data::XMC, v0data::YMC, v0data::ZMC,
+                            v0data::PxPosMC, v0data::PyPosMC, v0data::PzPosMC,
+                            v0data::PxNegMC, v0data::PyNegMC, v0data::PzNegMC,
+                            v0data::PxMC, v0data::PyMC, v0data::PzMC,
+                            v0data::RapidityMC<v0data::PxMC, v0data::PyMC, v0data::PzMC>,
+                            v0data::NegativePtMC<v0data::PxNegMC, v0data::PyNegMC>,
+                            v0data::PositivePtMC<v0data::PxPosMC, v0data::PyPosMC>,
+                            v0data::PtMC<v0data::PxMC, v0data::PyMC, v0data::PzMC>);
+
 DECLARE_SOA_TABLE(StoredV0MCCores_000, "AOD", "V0MCCORE", //! MC properties of the V0 for posterior analysis
                   v0data::PDGCode, v0data::PDGCodeMother,
                   v0data::PDGCodePositive, v0data::PDGCodeNegative,
@@ -683,6 +716,16 @@ DECLARE_SOA_TABLE_VERSIONED(StoredV0MCCores_001, "AOD", "V0MCCORE", 1, //! debug
                             v0data::PxNegMC, v0data::PyNegMC, v0data::PzNegMC,
                             o2::soa::Marker<1>);
 
+DECLARE_SOA_TABLE_VERSIONED(StoredV0MCCores_002, "AOD", "V0MCCORE", 2, //! debug information
+                            v0data::ParticleIdMC,                      //! MC properties of the V0 for posterior analysis
+                            v0data::PDGCode, v0data::PDGCodeMother,
+                            v0data::PDGCodePositive, v0data::PDGCodeNegative,
+                            v0data::IsPhysicalPrimary, v0data::XMC, v0data::YMC, v0data::ZMC,
+                            v0data::PxPosMC, v0data::PyPosMC, v0data::PzPosMC,
+                            v0data::PxNegMC, v0data::PyNegMC, v0data::PzNegMC,
+                            v0data::PxMC, v0data::PyMC, v0data::PzMC,
+                            o2::soa::Marker<1>);
+
 DECLARE_SOA_TABLE(V0MCCollRefs, "AOD", "V0MCCOLLREF", //! refers MC candidate back to proper MC Collision
                   o2::soa::Index<>, v0data::StraMCCollisionId, o2::soa::Marker<2>);
 
@@ -696,8 +739,8 @@ DECLARE_SOA_TABLE(V0MCMothers, "AOD", "V0MCMOTHER", //! optional table for MC mo
 DECLARE_SOA_TABLE(StoredV0MCMothers, "AOD1", "V0MCMOTHER", //! optional table for MC mothers
                   o2::soa::Index<>, v0data::MotherMCPartId, o2::soa::Marker<1>);
 
-using V0MCCores = V0MCCores_001;
-using StoredV0MCCores = StoredV0MCCores_001;
+using V0MCCores = V0MCCores_002;
+using StoredV0MCCores = StoredV0MCCores_002;
 
 using V0Index = V0Indices::iterator;
 using V0Core = V0Cores::iterator;
@@ -708,6 +751,7 @@ using V0fCDatas = soa::Join<V0fCIndices, V0fCTrackXs, V0fCCores>;
 using V0fCData = V0fCDatas::iterator;
 using V0MCDatas = soa::Join<V0MCCores, V0MCMothers>;
 using V0MCData = V0MCDatas::iterator;
+using V0MCCore = V0MCCores::iterator;
 
 // definitions of indices for interlink tables
 namespace v0data
@@ -997,6 +1041,24 @@ DECLARE_SOA_DYNAMIC_COLUMN(BachelorEta, bacheloreta, //! bachelor daughter eta
                            [](float PxPos, float PyPos, float PzPos) -> float { return RecoDecay::eta(std::array{PxPos, PyPos, PzPos}); });
 DECLARE_SOA_DYNAMIC_COLUMN(BachelorPhi, bachelorphi, //! bachelor daughter phi
                            [](float PxPos, float PyPos) -> float { return RecoDecay::phi(PxPos, PyPos); });
+
+DECLARE_SOA_DYNAMIC_COLUMN(RapidityMC, rapidityMC, //! rapidity (0, 1: Xi; 2, 3: Omega)
+                           [](float PxMC, float PyMC, float PzMC, int value) -> float {
+                             if (value == 0 || value == 1)
+                               return RecoDecay::y(std::array{PxMC, PyMC, PzMC}, o2::constants::physics::MassXiMinus);
+                             if (value == 2 || value == 3)
+                               return RecoDecay::y(std::array{PxMC, PyMC, PzMC}, o2::constants::physics::MassOmegaMinus);
+                             return 0.0f;
+                           });
+
+DECLARE_SOA_DYNAMIC_COLUMN(NegativePtMC, negativeptMC, //! negative daughter pT
+                           [](float pxNegMC, float pyNegMC) -> float { return RecoDecay::sqrtSumOfSquares(pxNegMC, pyNegMC); });
+DECLARE_SOA_DYNAMIC_COLUMN(PositivePtMC, positiveptMC, //! positive daughter pT
+                           [](float pxPosMC, float pyPosMC) -> float { return RecoDecay::sqrtSumOfSquares(pxPosMC, pyPosMC); });
+DECLARE_SOA_DYNAMIC_COLUMN(BachelorPtMC, bachelorptMC, //! bachelor daughter pT
+                           [](float pxBachMC, float pyBachMC) -> float { return RecoDecay::sqrtSumOfSquares(pxBachMC, pyBachMC); });
+DECLARE_SOA_DYNAMIC_COLUMN(PtMC, ptMC, //! cascade pT
+                           [](float pxMC, float pyMC) -> float { return RecoDecay::sqrtSumOfSquares(pxMC, pyMC); });
 } // namespace cascdata
 
 //______________________________________________________
@@ -1190,7 +1252,12 @@ DECLARE_SOA_TABLE(CascMCCores, "AOD", "CASCMCCORE", //! bachelor-baryon correlat
                   cascdata::PxPosMC, cascdata::PyPosMC, cascdata::PzPosMC,
                   cascdata::PxNegMC, cascdata::PyNegMC, cascdata::PzNegMC,
                   cascdata::PxBachMC, cascdata::PyBachMC, cascdata::PzBachMC,
-                  cascdata::PxMC, cascdata::PyMC, cascdata::PzMC);
+                  cascdata::PxMC, cascdata::PyMC, cascdata::PzMC,
+                  cascdata::RapidityMC<cascdata::PxMC, cascdata::PyMC, cascdata::PzMC>,
+                  cascdata::NegativePtMC<cascdata::PxNegMC, cascdata::PyNegMC>,
+                  cascdata::PositivePtMC<cascdata::PxPosMC, cascdata::PyPosMC>,
+                  cascdata::BachelorPtMC<cascdata::PxBachMC, cascdata::PyBachMC>,
+                  cascdata::PtMC<cascdata::PxMC, cascdata::PyMC>);
 
 namespace cascdata
 {

PWGLF/TableProducer/Strangeness/Converters/CMakeLists.txt

@@ -27,4 +27,9 @@ o2physics_add_dpl_workflow(strarawcentsconverter2v4
 o2physics_add_dpl_workflow(v0coresconverter
                     SOURCES v0coresconverter.cxx
                     PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore
+                    COMPONENT_NAME Analysis)
+
+o2physics_add_dpl_workflow(v0coresconverter2
+                    SOURCES v0coresconverter2.cxx
+                    PUBLIC_LINK_LIBRARIES O2Physics::AnalysisCore
                     COMPONENT_NAME Analysis)

PWGLF/TableProducer/Strangeness/Converters/v0coresconverter2.cxx

@@ -0,0 +1,52 @@
+// 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.
+#include "Framework/runDataProcessing.h"
+#include "Framework/AnalysisTask.h"
+#include "Framework/AnalysisDataModel.h"
+#include "PWGLF/DataModel/LFStrangenessTables.h"
+
+using namespace o2;
+using namespace o2::framework;
+
+// Converts V0 version 001 to 002
+struct v0coresconverter2 {
+  Produces<aod::V0MCCores_002> v0MCCores_002;
+
+  void process(aod::V0MCCores_001 const& v0MCCores_001)
+  {
+    for (auto& values : v0MCCores_001) {
+      v0MCCores_002(0,
+                    values.pdgCode(),
+                    values.pdgCodeMother(),
+                    values.pdgCodePositive(),
+                    values.pdgCodeNegative(),
+                    values.isPhysicalPrimary(),
+                    values.xMC(),
+                    values.yMC(),
+                    values.zMC(),
+                    values.pxPosMC(),
+                    values.pyPosMC(),
+                    values.pzPosMC(),
+                    values.pxNegMC(),
+                    values.pyNegMC(),
+                    values.pzNegMC(),
+                    values.pxPosMC() + values.pxNegMC(),
+                    values.pyPosMC() + values.pyNegMC(),
+                    values.pzPosMC() + values.pzNegMC());
+    }
+  }
+};
+
+WorkflowSpec defineDataProcessing(ConfigContext const& cfgc)
+{
+  return WorkflowSpec{
+    adaptAnalysisTask<v0coresconverter2>(cfgc)};
+}

PWGLF/TableProducer/Strangeness/cascadebuilder.cxx

@@ -194,6 +194,8 @@ struct cascadeBuilder {
   Configurable<bool> kfDoDCAFitterPreMinimV0{"kfDoDCAFitterPreMinimV0", true, "KF: do DCAFitter pre-optimization before KF fit to include material corrections for V0"};
   Configurable<bool> kfDoDCAFitterPreMinimCasc{"kfDoDCAFitterPreMinimCasc", true, "KF: do DCAFitter pre-optimization before KF fit to include material corrections for Xi"};
 
+  // for using cascade momentum at prim. vtx
+  Configurable<bool> useCascadeMomentumAtPrimVtx{"useCascadeMomentumAtPrimVtx", false, "if enabled, store cascade momentum at prim. vtx instead of decay point (= default)"};
   // for topo var QA
   struct : ConfigurableGroup {
     ConfigurableAxis axisTopoVarPointingAngle{"axisConfigurations.axisTopoVarPointingAngle", {50, 0.0, 1.0}, "pointing angle"};
@@ -1090,6 +1092,13 @@ struct cascadeBuilder {
     o2::base::Propagator::Instance()->propagateToDCABxByBz({collision.posX(), collision.posY(), collision.posZ()}, lCascadeTrack, 2.f, matCorrCascade, &dcaInfo);
     cascadecandidate.cascDCAxy = dcaInfo[0];
     cascadecandidate.cascDCAz = dcaInfo[1];
+    if (useCascadeMomentumAtPrimVtx) {
+      lCascadeTrack.getPxPyPzGlo(cascadecandidate.cascademom);
+    } else {
+      cascadecandidate.cascademom[0] = cascadecandidate.bachP[0] + cascadecandidate.v0mompos[0] + cascadecandidate.v0momneg[0];
+      cascadecandidate.cascademom[1] = cascadecandidate.bachP[1] + cascadecandidate.v0mompos[1] + cascadecandidate.v0momneg[1];
+      cascadecandidate.cascademom[2] = cascadecandidate.bachP[2] + cascadecandidate.v0mompos[2] + cascadecandidate.v0momneg[2];
+    }
 
     // Calculate masses a priori
     cascadecandidate.mXi = RecoDecay::m(array{array{cascadecandidate.bachP[0], cascadecandidate.bachP[1], cascadecandidate.bachP[2]}, array{v0.pxpos() + v0.pxneg(), v0.pypos() + v0.pyneg(), v0.pzpos() + v0.pzneg()}}, array{o2::constants::physics::MassPionCharged, o2::constants::physics::MassLambda});
@@ -1555,9 +1564,7 @@ struct cascadeBuilder {
              cascadecandidate.v0mompos[0], cascadecandidate.v0mompos[1], cascadecandidate.v0mompos[2],
              cascadecandidate.v0momneg[0], cascadecandidate.v0momneg[1], cascadecandidate.v0momneg[2],
              cascadecandidate.bachP[0], cascadecandidate.bachP[1], cascadecandidate.bachP[2],
-             cascadecandidate.bachP[0] + cascadecandidate.v0mompos[0] + cascadecandidate.v0momneg[0], // <--- redundant but ok
-             cascadecandidate.bachP[1] + cascadecandidate.v0mompos[1] + cascadecandidate.v0momneg[1], // <--- redundant but ok
-             cascadecandidate.bachP[2] + cascadecandidate.v0mompos[2] + cascadecandidate.v0momneg[2], // <--- redundant but ok
+             cascadecandidate.cascademom[0], cascadecandidate.cascademom[1], cascadecandidate.cascademom[2],
              cascadecandidate.v0dcadau, cascadecandidate.dcacascdau,
              cascadecandidate.v0dcapostopv, cascadecandidate.v0dcanegtopv,
              cascadecandidate.bachDCAxy, cascadecandidate.cascDCAxy, cascadecandidate.cascDCAz); // <--- no corresponding stratrack information available

PWGLF/TableProducer/Strangeness/cascademcbuilder.cxx

@@ -76,6 +76,9 @@ struct cascademcbuilder {
     int pdgCodeNegative;
     int pdgCodeBachelor;
     bool isPhysicalPrimary;
+    int processPositive = -1;
+    int processNegative = -1;
+    int processBachelor = -1;
     std::array<float, 3> xyz;
     std::array<float, 3> lxyz;
     std::array<float, 3> posP;
@@ -89,8 +92,6 @@ struct cascademcbuilder {
   mcCascinfo thisInfo;
   //*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*+-+*
 
-  void init(InitContext const&) {}
-
   template <typename TCascadeTable, typename TMCParticleTable>
   void generateCascadeMCinfo(TCascadeTable cascTable, TMCParticleTable mcParticles)
   {
@@ -142,6 +143,9 @@ struct cascademcbuilder {
         thisInfo.bachP[0] = lMCBachTrack.px();
         thisInfo.bachP[1] = lMCBachTrack.py();
         thisInfo.bachP[2] = lMCBachTrack.pz();
+        thisInfo.processPositive = lMCPosTrack.getProcess();
+        thisInfo.processNegative = lMCNegTrack.getProcess();
+        thisInfo.processBachelor = lMCBachTrack.getProcess();
 
         // Step 1: check if the mother is the same, go up a level
         if (lMCNegTrack.has_mothers() && lMCPosTrack.has_mothers()) {
@@ -194,8 +198,8 @@ struct cascademcbuilder {
         thisInfo.label, thisInfo.motherLabel);
 
       // Mark mcParticle as recoed (no searching necessary afterwards)
-      if (thisInfo.motherLabel > -1) {
-        mcParticleIsReco[thisInfo.motherLabel] = true;
+      if (thisInfo.label > -1) {
+        mcParticleIsReco[thisInfo.label] = true;
       }
 
       if (populateCascMCCoresSymmetric) {
@@ -216,10 +220,7 @@ struct cascademcbuilder {
         // step 1: check if this element is already provided in the table
         //         using the packedIndices variable calculated above
         for (uint32_t ii = 0; ii < mcCascinfos.size(); ii++) {
-          if (
-            thisInfo.mcParticlePositive == mcCascinfos[ii].mcParticlePositive && mcCascinfos[ii].mcParticlePositive > 0 &&
-            thisInfo.mcParticleNegative == mcCascinfos[ii].mcParticleNegative && mcCascinfos[ii].mcParticleNegative > 0 &&
-            thisInfo.mcParticleBachelor == mcCascinfos[ii].mcParticleBachelor && mcCascinfos[ii].mcParticleBachelor > 0) {
+          if (thisInfo.label == mcCascinfos[ii].label && mcCascinfos[ii].label > -1) {
             thisCascMCCoreIndex = ii;
             break; // this exists already in list
           }
@@ -272,6 +273,54 @@ struct cascademcbuilder {
           thisInfo.momentum[0] = mcParticle.px();
           thisInfo.momentum[1] = mcParticle.py();
           thisInfo.momentum[2] = mcParticle.pz();
+          thisInfo.label = mcParticle.globalIndex();
+
+          if (mcParticle.has_daughters()) {
+            auto const& daughters = mcParticle.template daughters_as<aod::McParticles>();
+            for (auto& dau : daughters) {
+              if (dau.getProcess() != 4) // check whether the daughter comes from a decay
+                continue;
+
+              if (TMath::Abs(dau.pdgCode()) == 211 || TMath::Abs(dau.pdgCode()) == 321) {
+                thisInfo.pdgCodeBachelor = dau.pdgCode();
+                thisInfo.bachP[0] = dau.px();
+                thisInfo.bachP[1] = dau.py();
+                thisInfo.bachP[2] = dau.pz();
+                thisInfo.xyz[0] = dau.vx();
+                thisInfo.xyz[1] = dau.vy();
+                thisInfo.xyz[2] = dau.vz();
+                thisInfo.mcParticleBachelor = dau.globalIndex();
+              }
+              if (TMath::Abs(dau.pdgCode()) == 2212) {
+                thisInfo.pdgCodeV0 = dau.pdgCode();
+
+                for (auto& v0Dau : dau.template daughters_as<aod::McParticles>()) {
+                  if (v0Dau.getProcess() != 4)
+                    continue;
+
+                  if (v0Dau.pdgCode() > 0) {
+                    thisInfo.pdgCodePositive = v0Dau.pdgCode();
+                    thisInfo.processPositive = v0Dau.getProcess();
+                    thisInfo.posP[0] = v0Dau.px();
+                    thisInfo.posP[1] = v0Dau.py();
+                    thisInfo.posP[2] = v0Dau.pz();
+                    thisInfo.lxyz[0] = v0Dau.vx();
+                    thisInfo.lxyz[1] = v0Dau.vy();
+                    thisInfo.lxyz[2] = v0Dau.vz();
+                    thisInfo.mcParticlePositive = v0Dau.globalIndex();
+                  }
+                  if (v0Dau.pdgCode() < 0) {
+                    thisInfo.pdgCodeNegative = v0Dau.pdgCode();
+                    thisInfo.processNegative = v0Dau.getProcess();
+                    thisInfo.negP[0] = v0Dau.px();
+                    thisInfo.negP[1] = v0Dau.py();
+                    thisInfo.negP[2] = v0Dau.pz();
+                    thisInfo.mcParticleNegative = v0Dau.globalIndex();
+                  }
+                }
+              }
+            }
+          }
 
           // if I got here, it means this MC particle was not recoed and is of interest. Add it please
           mcCascinfos.push_back(thisInfo);