Extend analysis tools and LUT creation routine (#86)

- Configure analysis workflows from config file
- Update workflows
- Add option for nuclei when creating AODs
- Extend LUT creator, add possibility to set path, particles and outputpath

Author: njacazio

examples/aod/createO2tables.C

@@ -57,6 +57,7 @@ const char* inputFileAccMuonPID = "muonAccEffPID.root";
 
 // Simulation parameters
 const bool do_vertexing = true;
+const bool enable_nuclei = false;
 
 int createO2tables(const char* inputFile = "delphes.root",
                    const char* outputFile = "AODRun5.root",
@@ -98,6 +99,10 @@ int createO2tables(const char* inputFile = "delphes.root",
   smearer.loadTable(211, "lutCovm.pi.dat");
   smearer.loadTable(321, "lutCovm.ka.dat");
   smearer.loadTable(2212, "lutCovm.pr.dat");
+  if (enable_nuclei) {
+    smearer.loadTable(1000010020, "lutCovm.de.dat");
+    smearer.loadTable(1000020030, "lutCovm.he3.dat");
+  }
 
   // TOF layer
   o2::delphes::TOFLayer tof_layer;
@@ -356,7 +361,7 @@ int createO2tables(const char* inputFile = "delphes.root",
     std::array<float, 2> ftzero;
 
     forward_tof_layer.eventTime(ftof_tracks, ftzero);
-    for (int i = 0; i < ftof_tracks.size(); i++) {
+    for (unsigned int i = 0; i < ftof_tracks.size(); i++) {
       auto track = ftof_tracks[i];
       ftof.fIndexCollisions = ftof_tracks_indices[i].first;
       ftof.fIndexTracks = ftof_tracks_indices[i].second; // Index in the Track table

examples/scripts/common.py

@@ -7,6 +7,7 @@
 import argparse
 import multiprocessing
 import sys
+import os
 try:
     import tqdm
 except ImportError as e:
@@ -79,3 +80,26 @@ def run_in_parallel(processes, job_runner, job_arguments, job_message):
             r = list(tqdm.tqdm(pool.imap(job_runner, job_arguments),
                                total=len(job_arguments),
                                bar_format='{l_bar}{bar:10}{r_bar}{bar:-10b}'))
+
+
+def run_cmd(cmd, comment="", check_status=True):
+    """
+    Function to run a command in bash, allows to check the status of the command and to log the command output
+    """
+    verbose_msg("Running", f"'{cmd}'", bcolors.BOKBLUE + comment)
+    try:
+        to_run = cmd
+        if check_status:
+            to_run = f"{cmd} && echo OK"
+        content = os.popen(to_run).read()
+        if content:
+            content = content.strip()
+            for i in content.strip().split("\n"):
+                verbose_msg("++", i)
+        if "Encountered error" in content:
+            warning_msg("Error encountered runtime error in", cmd)
+        if check_status:
+            if "OK" not in content and "root" not in cmd:
+                fatal_msg("Command", cmd, "does not have the OK tag", content)
+    except:
+        fatal_msg("Error while running", f"'{cmd}'")

examples/scripts/createO2tables.py

@@ -11,30 +11,7 @@
 import glob
 import random
 from datetime import datetime
-from common import bcolors, msg, fatal_msg, verbose_msg, run_in_parallel, set_verbose_mode, get_default_parser, warning_msg
-
-
-def run_cmd(cmd, comment="", check_status=True):
-    """
-    Function to run a command in bash, allows to check the status of the command and to log the command output
-    """
-    verbose_msg("Running", f"'{cmd}'", bcolors.BOKBLUE + comment)
-    try:
-        to_run = cmd
-        if check_status:
-            to_run = f"{cmd} && echo OK"
-        content = os.popen(to_run).read()
-        if content:
-            content = content.strip()
-            for i in content.strip().split("\n"):
-                verbose_msg("++", i)
-        if "Encountered error" in content:
-            warning_msg("Error encountered runtime error in", cmd)
-        if check_status:
-            if "OK" not in content and "root" not in cmd:
-                fatal_msg("Command", cmd, "does not have the OK tag", content)
-    except:
-        fatal_msg("Error while running", f"'{cmd}'")
+from common import bcolors, msg, fatal_msg, verbose_msg, run_in_parallel, set_verbose_mode, get_default_parser, run_cmd
 
 
 def process_run(run_number):
@@ -60,7 +37,8 @@ def main(configuration_file,
          clean_delphes_files,
          create_luts,
          turn_off_vertexing,
-         append_production):
+         append_production,
+         use_nuclei):
     arguments = locals()  # List of arguments to put into the log
     parser = configparser.RawConfigParser()
     parser.read(configuration_file)
@@ -125,6 +103,9 @@ def do_copy(in_file, out_file=".", in_path=None):
 
     lut_path = opt("lut_path")
     lut_tag = opt("lut_tag")
+    lut_particles = ["el", "mu", "pi", "ka", "pr"]
+    if use_nuclei:
+        lut_particles += ["de", "he3"]
     if create_luts:
         # Creating LUTs
         minimum_track_radius = opt("minimum_track_radius")
@@ -135,12 +116,13 @@ def do_copy(in_file, out_file=".", in_path=None):
     else:
         # Fetching LUTs
         verbose_msg(f"Fetching LUTs with tag {lut_tag} from path {lut_path}")
-        for i in ["el", "mu", "pi", "ka", "pr"]:
+        for i in lut_particles:
             lut_bg = "{}kG".format(bField).replace(".", "")
-            do_copy(f"lutCovm.{i}.{lut_bg}.{lut_tag}.dat", f"lutCovm.{i}.dat", in_path=lut_path)
+            do_copy(f"lutCovm.{i}.{lut_bg}.{lut_tag}.dat",
+                    f"lutCovm.{i}.dat", in_path=lut_path)
 
     # Checking that we actually have LUTs
-    for i in ["el", "mu", "pi", "ka", "pr"]:
+    for i in lut_particles:
         i = f"lutCovm.{i}.dat"
         if not os.path.isfile(i):
             fatal_msg("Did not find LUT file", i)
@@ -236,6 +218,9 @@ def set_config(config_file, config, value):
     if turn_off_vertexing:
         set_config("createO2tables.C",
                    "const bool do_vertexing = ", "false\;/")
+    if use_nuclei:
+        set_config("createO2tables.C",
+                   "const bool enable_nuclei = ", "true\;/")
     else:  # Check that the geometry file for the vertexing is there
         if not os.path.isfile("o2sim_grp.root") or not os.path.isfile("o2sim_geometry.root"):
             run_cmd("mkdir tmpo2sim && cd tmpo2sim && o2-sim -m PIPE ITS MFT -g boxgen -n 1 -j 1 --configKeyValues 'BoxGun.number=1' && cp o2sim_grp.root .. && cp o2sim_geometry.root .. && cd .. && rm -r tmpo2sim")
@@ -359,7 +344,10 @@ def copy_and_link(file_name):
                             log_file=aod_log_file,
                             check_status=True)
             # Check that there were no O2 errors
-            write_to_runner(f"if grep -q \"\[ERROR\]\" {aod_log_file}; then echo \": got some errors in '{aod_log_file}'\" && exit 1; fi")
+            write_to_runner(
+                f"if grep -q \"\[ERROR\]\" {aod_log_file}; then echo \": got some errors in '{aod_log_file}'\" && exit 1; fi")
+            write_to_runner(
+                f"if grep -q \"\[FATAL\]\" {aod_log_file}; then echo \": got some fatals in '{aod_log_file}'\" && exit 1; fi")
             # Rename the temporary AODs to standard AODs
             write_to_runner(f"mv tmp_{aod_file} {aod_file}", check_status=True)
             if not clean_delphes_files:
@@ -481,6 +469,9 @@ def write_config(entry, prefix=""):
     parser.add_argument("--append", "-a",
                         action="store_true",
                         help="Option to append the results instead of starting over by shifting the AOD indexing. N.B. the user is responsible of the compatibility between appended AODs. Only works in conjuction by specifying an output path (option '-o')")
+    parser.add_argument("--nuclei",
+                        action="store_true",
+                        help="Option use nuclei LUTs")
     parser.add_argument("--use-preexisting-luts", "-l",
                         action="store_true",
                         help="Option to use preexisting LUTs instead of creating new ones, in this case LUTs with the requested tag are fetched from the LUT path. By default new LUTs are created at each run.")
@@ -501,4 +492,5 @@ def write_config(entry, prefix=""):
          qa=args.qa,
          create_luts=not args.use_preexisting_luts,
          turn_off_vertexing=args.no_vertexing,
-         append_production=args.append)
+         append_production=args.append,
+         use_nuclei=args.nuclei)

examples/scripts/create_luts.sh

@@ -5,27 +5,36 @@ if [[ $* == *"-h"* ]]; then
     echo "1) tag of the LUT writer [default]"
     echo "2) Magnetic field in T [0.5]"
     echo "3) Minimum radius of the track in cm [100]"
+    echo "4) Path where the LUT writers are located [\$DELPHESO2_ROOT/lut/]"
+    echo "5) Output path where to write the LUTs [.]"
+    echo "6) Particles to consider [\"0 1 2 3 4\"]"
     exit 0
 fi
 
 WHAT=default
 FIELD=0.5
 RMIN=100.
+OUT_PATH=.
+WRITER_PATH=$DELPHESO2_ROOT/lut/
+PARTICLES="0 1 2 3 4"
 
 [ -z "$1" ] || WHAT=$1
 [ -z "$2" ] || FIELD=$2
 [ -z "$3" ] || RMIN=$3
+[ -z "$4" ] || WRITER_PATH=$4
+[ -z "$5" ] || OUT_PATH=$5
+[ -z "$6" ] || PARTICLES=$6
 
-cp    "$DELPHESO2_ROOT/lut/lutWrite.$WHAT.cc" . || { echo "cannot find lut writer: $DELPHESO2_ROOT/lut/lutWrite.$WHAT.cc" ; exit 1; }
-cp    "$DELPHESO2_ROOT/lut/DetectorK/DetectorK.cxx" .
-cp    "$DELPHESO2_ROOT/lut/DetectorK/DetectorK.h" .
-cp -r "$DELPHESO2_ROOT/lut/fwdRes" .
-cp    "$DELPHESO2_ROOT/lut/lutWrite.cc" .
-cp    "$DELPHESO2_ROOT/lut/lutCovm.hh" .
+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" .
+cp    "${WRITER_PATH}/lutCovm.hh" .
 
 echo " --- creating LUTs: config = $WHAT, field = $FIELD T, min tracking radius = $RMIN cm"
 
-for i in 0 1 2 3 4; do
+for i in $PARTICLES; do
     root -l -b <<EOF
     .L DetectorK.cxx+
     .L lutWrite.${WHAT}.cc
@@ -39,19 +48,22 @@ for i in 0 1 2 3 4; do
     const TString pn[7] = {"el", "mu", "pi", "ka", "pr", "de", "he3"};
     const int pc[7] = {11, 13, 211, 321, 2212, 1000010020, 1000020030 };
     const float field = ${FIELD}f;
-    const float rmin = ${RMIN}f;
+    const float rmin = ${RMIN};
     const int i = ${i};
-    lutWrite_${WHAT}("lutCovm." + pn[i] + ".dat", pc[i], field, rmin);
+    lutWrite_${WHAT}("${OUT_PATH}/lutCovm." + pn[i] + ".dat", pc[i], field, rmin);
 
 EOF
 done
 
 # Checking that the output LUTs are OK
 NullSize=""
-for i in el mu pi ka pr; do
-    if [[ ! -s lutCovm.$i.dat ]]; then
+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"
         NullSize="${NullSize} ${i}"
+    else
+        echo "lutCovm.${P[$i]}.dat is ok"
     fi
 done