Add nuclei LUTs in tracking with checks in AO2D creation (#93)

* Add nuclei LUTs in tracking

* Extend analysis tools

* Using dNdEta to smear tracks

* Update LUT maker script

Author: njacazio

examples/aod/createO2tables.C

@@ -94,14 +94,21 @@ int createO2tables(const char* inputFile = "delphes.root",
 
   // smearer
   o2::delphes::TrackSmearer smearer;
-  smearer.loadTable(11, "lutCovm.el.dat");
-  smearer.loadTable(13, "lutCovm.mu.dat");
-  smearer.loadTable(211, "lutCovm.pi.dat");
-  smearer.loadTable(321, "lutCovm.ka.dat");
-  smearer.loadTable(2212, "lutCovm.pr.dat");
+  std::map<int, const char*> mapPdgLut;
+  mapPdgLut.insert(std::make_pair(11, "lutCovm.el.dat"));
+  mapPdgLut.insert(std::make_pair(13, "lutCovm.mu.dat"));
+  mapPdgLut.insert(std::make_pair(211, "lutCovm.pi.dat"));
+  mapPdgLut.insert(std::make_pair(321, "lutCovm.ka.dat"));
+  mapPdgLut.insert(std::make_pair(2212, "lutCovm.pr.dat"));
   if (enable_nuclei) {
-    smearer.loadTable(1000010020, "lutCovm.de.dat");
-    smearer.loadTable(1000020030, "lutCovm.he3.dat");
+    mapPdgLut.insert(std::make_pair(1000010020, "lutCovm.de.dat"));
+    mapPdgLut.insert(std::make_pair(1000020030, "lutCovm.he.dat"));
+  }
+  for (auto e : mapPdgLut) {
+    if (!smearer.loadTable(e.first, e.second)) {
+      Printf("Having issue with loading the LUT %i '%s'", e.first, e.second);
+      return 1;
+    }
   }
 
   // TOF layer
@@ -170,6 +177,8 @@ int createO2tables(const char* inputFile = "delphes.root",
 
     // Load selected branches with data from specified event
     treeReader->ReadEntry(ientry);
+    const float multEtaRange = 2.f; // Range in eta to count the charged particles
+    float dNdEta = 0.f;             // Charged particle multiplicity to use in the efficiency evaluation
 
     for (Int_t iparticle = 0; iparticle < particles->GetEntries(); ++iparticle) { // Loop over particles
       auto particle = (GenParticle*)particles->At(iparticle);
@@ -207,8 +216,12 @@ int createO2tables(const char* inputFile = "delphes.root",
       mcparticle.fVz = particle->Z * 0.1;
       mcparticle.fVt = particle->T;
 
+      if (TMath::Abs(particle->Eta) <= multEtaRange && particle->D1 < 0 && particle->D2 < 0 && particle->Charge != 0) {
+        dNdEta += 1.f;
+      }
       FillTree(kMcParticle);
     }
+    dNdEta = 0.5f * dNdEta / multEtaRange;
     fOffsetLabel += particles->GetEntries();
 
     // loop over tracks
@@ -246,7 +259,7 @@ int createO2tables(const char* inputFile = "delphes.root",
 
       O2Track o2track; // tracks in internal O2 format
       o2::delphes::TrackUtils::convertTrackToO2Track(*track, o2track, true);
-      if (!smearer.smearTrack(o2track, track->PID)) { // Skipping inefficient/not correctly smeared tracks
+      if (!smearer.smearTrack(o2track, track->PID, dNdEta)) { // Skipping inefficient/not correctly smeared tracks
         continue;
       }
       o2::delphes::TrackUtils::convertO2TrackToTrack(o2track, *track, true);

examples/scripts/create_luts.sh

@@ -60,11 +60,19 @@ while getopts ":t:B:R:p:o:T:P:h" option; do
     esac
 done
 
-cp    "${WRITER_PATH}/lutWrite.$WHAT.cc" . || { echo "cannot find lut writer: ${WRITER_PATH}/lutWrite.$WHAT.cc" ; exit 1; }
-cp    "${WRITER_PATH}/DetectorK/DetectorK.cxx" .
-cp    "${WRITER_PATH}/DetectorK/DetectorK.h" .
-cp -r "${WRITER_PATH}/fwdRes" .
-cp    "${WRITER_PATH}/lutWrite.cc" .
+function do_copy() {
+    cp "${1}" . || {
+        echo "Cannot find $2: ${1}"
+        exit 1
+    }
+}
+
+do_copy "${WRITER_PATH}/lutWrite.$WHAT.cc" "lut writer"
+do_copy "${WRITER_PATH}/DetectorK/DetectorK.cxx"
+do_copy "${WRITER_PATH}/DetectorK/DetectorK.h"
+do_copy "${WRITER_PATH}/lutWrite.cc"
+do_copy "${WRITER_PATH}/lutCovm.hh"
+cp -r   "${WRITER_PATH}/fwdRes" .
 
 echo " --- creating LUTs: config = ${WHAT}, field = ${FIELD} T, min tracking radius = ${RMIN} cm"
 
@@ -93,11 +101,12 @@ done
 NullSize=""
 P=(el mu pi ka pr de he3)
 for i in $PARTICLES; do
-    if [[ ! -s lutCovm.${P[$i]}.dat ]]; then
-        echo "${i} has zero size"
+    LUT_FILE=${OUT_PATH}/lutCovm.${P[$i]}${OUT_TAG}.dat
+    if [[ ! -s ${LUT_FILE} ]]; then
+        echo "LUT file ${LUT_FILE} for particle ${i} has zero size"
         NullSize="${NullSize} ${i}"
     else
-        echo "lutCovm.${P[$i]}.dat is ok"
+        echo "${LUT_FILE} is ok"
     fi
 done
 

examples/scripts/diagnostic_tools/doanalysis.py

@@ -124,6 +124,7 @@ def main(mode,
          merge_output=False,
          merge_only=False,
          shm_mem_size=16000000000,
+         rate_lim=1000000000,
          readers=1,
          avoid_overwriting_merge=False,
          extra_arguments=""):
@@ -139,7 +140,7 @@ def main(mode,
     else:
         msg("Merging output of", f"'{mode}'",
             "analysis", color=bcolors.BOKBLUE)
-    o2_arguments = f"-b --shm-segment-size {shm_mem_size} --readers {readers}"
+    o2_arguments = f"-b --shm-segment-size {shm_mem_size} --aod-memory-rate-limit {rate_lim} --readers {readers}"
     o2_arguments += extra_arguments
     if mode not in analyses:
         raise ValueError("Did not find analyses matching mode",

src/TrackSmearer.hh

@@ -34,21 +34,23 @@ public:
 
   int getIndexPDG(int pdg) {
     switch(abs(pdg)) {
-    case 11: return 0;
-    case 13: return 1;
-    case 211: return 2;
-    case 321: return 3;
-    case 2212: return 4;
-    default: return 2;
+    case 11: return 0; // Electron
+    case 13: return 1; // Muon
+    case 211: return 2; // Pion
+    case 321: return 3; // Kaon
+    case 2212: return 4; // Proton
+    case 1000010020: return 5; // Deuteron
+    case 1000020030: return 6; // Helium3
+    default: return 2; // Default: pion
     };
   };
 
   void setdNdEta(float val) { mdNdEta = val; };
 
 protected:
-
-  lutHeader_t *mLUTHeader[5] = {nullptr};
-  lutEntry_t *****mLUTEntry[5] = {nullptr};
+  static constexpr unsigned int nLUTs = 7;
+  lutHeader_t *mLUTHeader[nLUTs] = {nullptr};
+  lutEntry_t *****mLUTEntry[nLUTs] = {nullptr};
   bool mUseEfficiency = true;
   float mdNdEta =  1600.;