MCUtils: functionality to query ALICE properties on MCtracks/particles

Author: sawenzel

DataFormats/simulation/CMakeLists.txt

@@ -18,6 +18,7 @@ o2_add_library(SimulationDataFormat
                        src/StackParam.cxx
                        src/MCEventHeader.cxx
                        src/CustomStreamers.cxx
+                       src/MCUtils.cxx
                PUBLIC_LINK_LIBRARIES Microsoft.GSL::GSL
                                      O2::DetectorsCommonDataFormats
                                      O2::GPUCommon O2::DetectorsBase
@@ -40,6 +41,7 @@ o2_target_root_dictionary(
           include/SimulationDataFormat/MCEventHeader.h
           include/SimulationDataFormat/MCEventStats.h
           include/SimulationDataFormat/IOMCTruthContainerView.h
+          include/SimulationDataFormat/MCUtils.h
   LINKDEF src/SimulationDataLinkDef.h)
 # note the explicit LINKDEF as the linkdef in src is
 #

DataFormats/simulation/include/SimulationDataFormat/MCUtils.h

@@ -0,0 +1,78 @@
+// 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.
+
+/// \file MCUtils.h
+/// \brief Utility functions for MC particles
+/// \author S. Wenzel - November 2021
+
+#ifndef O2_MCUTILS_H
+#define O2_MCUTILS_H
+
+#include <SimulationDataFormat/MCTrack.h>
+#include "TPDGCode.h"
+
+namespace o2
+{
+namespace mcutils
+{
+/// Function to determine if a MC track/particle p is a primary according to physics criteria.
+/// Needs the particle as input, as well as the whole navigable container of particles
+/// (of which p needs to be a part itself). The container can be fetched via MCKinematicsReader.
+bool isPhysicalPrimary(o2::MCTrack const& p, std::vector<o2::MCTrack> const& pcontainer);
+
+/// Determine if a particle (identified by pdg) is stable
+inline bool isStable(int pdg)
+{
+  //
+  // Decide whether particle (pdg) is stable
+  //
+
+  // All ions/nucleons are considered as stable
+  // Nuclear code is 10LZZZAAAI
+  if (pdg > 1000000000) {
+    return true;
+  }
+
+  const int kNstable = 18;
+  int pdgStable[kNstable] = {
+    kGamma,      // Photon
+    kElectron,   // Electron
+    kMuonPlus,   // Muon
+    kPiPlus,     // Pion
+    kKPlus,      // Kaon
+    kK0Short,    // K0s
+    kK0Long,     // K0l
+    kProton,     // Proton
+    kNeutron,    // Neutron
+    kLambda0,    // Lambda_0
+    kSigmaMinus, // Sigma Minus
+    kSigmaPlus,  // Sigma Plus
+    3312,        // Xsi Minus
+    3322,        // Xsi
+    3334,        // Omega
+    kNuE,        // Electron Neutrino
+    kNuMu,       // Muon Neutrino
+    kNuTau       // Tau Neutrino
+  };
+
+  // this is linear search ---> a hash map or binary search should be more appropriate??
+  for (int i = 0; i < kNstable; ++i) {
+    if (pdg == std::abs(pdgStable[i])) {
+      return true;
+    }
+  }
+  return false;
+}
+
+} // namespace mcutils
+} // namespace o2
+
+#endif //O2_MCUTILS_H

DataFormats/simulation/src/MCUtils.cxx

@@ -0,0 +1,121 @@
+// 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.
+
+//
+// Created by Sandro Wenzel on 03.11.21.
+//
+
+#include <SimulationDataFormat/MCUtils.h>
+
+using namespace o2::mcutils;
+
+// taken from AliRoot and adapted to use of o2::MCTrack class
+bool isPhysicalPrimary(o2::MCTrack const& p, std::vector<o2::MCTrack> const& pcontainer)
+{
+  // Test if a particle is a physical primary according to the following definition:
+  // Particles produced in the collision including products of strong and
+  // electromagnetic decay and excluding feed-down from weak decays of strange
+  // particles.
+  //
+
+  const int ist = p.getStatusCode(); // the generator status code
+  const int pdg = std::abs(p.GetPdgCode());
+  //
+  // Initial state particle
+  // Solution for K0L decayed by Pythia6
+  // ->
+  // ist > 1 --> essentially means unstable
+  if ((ist > 1) && (pdg != 130) && p.isPrimary()) {
+    return false;
+  }
+  if ((ist > 1) && p.isSecondary()) {
+    return false;
+  }
+  // <-
+
+  if (!isStable(pdg)) {
+    return false;
+  }
+
+  if (p.isPrimary()) {
+    //
+    // Particle produced by generator
+    // Solution for K0L decayed by Pythia6
+    // ->
+    const int ipm = p.getMotherTrackId();
+    if (ipm > -1) {
+      const auto& ppm = pcontainer[ipm];
+      if (std::abs(ppm.GetPdgCode()) == 130) {
+        return false;
+      }
+    }
+    // <-
+    // check for direct photon in parton shower
+    // ->
+    const int ipd = p.getFirstDaughterTrackId();
+    if (pdg == 22 && ipd > -1) {
+      const auto& ppd = pcontainer[ipd];
+      if (ppd.GetPdgCode() == 22) {
+        return false;
+      }
+    }
+    // <-
+    return true;
+  } else {
+    //
+    // Particle produced during transport
+    //
+
+    int imo = p.getMotherTrackId();
+    auto pm = &(pcontainer[imo]);
+    int mpdg = std::abs(pm->GetPdgCode());
+    // Check for Sigma0
+    if ((mpdg == 3212) && pm->isPrimary()) {
+      return true;
+    }
+
+    //
+    // Check if it comes from a pi0 decay
+    //
+    if ((mpdg == kPi0) && pm->isPrimary()) {
+      return true;
+    }
+
+    // Check if this is a heavy flavor decay product
+    int mfl = int(mpdg / std::pow(10, int(std::log10(mpdg))));
+    //
+    // Light hadron
+    if (mfl < 4) {
+      return false;
+    }
+
+    //
+    // Heavy flavor hadron produced by generator
+    if (pm->isPrimary()) {
+      return true;
+    }
+
+    // To be sure that heavy flavor has not been produced in a secondary interaction
+    // Loop back to the generated mother
+    while (!pm->isPrimary()) {
+      imo = pm->getMotherTrackId();
+      pm = &(pcontainer[imo]);
+    }
+    mpdg = std::abs(pm->GetPdgCode());
+    mfl = int(mpdg / std::pow(10, int(std::log10(mpdg))));
+
+    if (mfl < 4) {
+      return false;
+    } else {
+      return true;
+    }
+  } // end else branch produced by generator
+}