Install root files for AOD making (#85)

- Add extra features and checks for AOD making
- Add helper for LUT creation tool
- Add check on created LUTs
- Linearize LUT maker script
- Fix empty LUTs for negative particles

Author: njacazio

examples/CMakeLists.txt

@@ -29,6 +29,7 @@ file(GLOB SMEARING
 file(GLOB AOD
   aod/*.C
   aod/*.h
+  aod/*.root
   )
 
 install(FILES    ${PYTHIA8}       DESTINATION examples/pythia8)

examples/scripts/createO2tables.py

@@ -109,15 +109,16 @@ def opt(entry, require=True):
     barrel_half_length = opt("barrel_half_length")
 
     # copy relevant files in the working directory
-    def do_copy(in_file, out_file):
+    def do_copy(in_file, out_file=".", in_path=None):
+        """Function to copy files"""
+        if in_path is not None:
+            in_file = os.path.join(in_path, in_file)
         in_file = os.path.expanduser(os.path.expandvars(in_file))
         verbose_msg("Copying", in_file, "to", out_file)
         shutil.copy2(in_file, out_file)
 
     # Fetching the propagation card
-    do_copy(os.path.join(opt("card_path"),
-                         opt("propagate_card")),
-            "propagate.tcl")
+    do_copy(opt("propagate_card"), "propagate.tcl", in_path=opt("card_path"))
 
     lut_path = opt("lut_path")
     lut_tag = opt("lut_tag")
@@ -133,8 +134,7 @@ def do_copy(in_file, out_file):
         verbose_msg(f"Fetching LUTs with tag {lut_tag} from path {lut_path}")
         for i in ["el", "mu", "pi", "ka", "pr"]:
             lut_bg = "{}kG".format(bField).replace(".", "")
-            do_copy(os.path.join(lut_path, f"lutCovm.{i}.{lut_bg}.{lut_tag}.dat"),
-                    f"lutCovm.{i}.dat")
+            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"]:
@@ -150,7 +150,7 @@ def do_copy(in_file, out_file):
                       config_entry, "in your configuration file", configuration_file)
         generators = opt("generators").split(" ")
         for i in generators:
-            do_copy(i, ".")
+            do_copy(i)
         msg("Using pythia with configuration", generators)
     else:
         def check_duplicate(option_name):
@@ -201,11 +201,11 @@ def check_duplicate(option_name):
     msg("  etaMax               =", etaMax)
 
     aod_path = opt("aod_path")
-    do_copy(os.path.join(aod_path, "createO2tables.h"), ".")
-    do_copy(os.path.join(aod_path, "createO2tables.C"), ".")
-    do_copy("../aod/muonAccEffPID.root", ".")
+    do_copy("createO2tables.h", in_path=aod_path)
+    do_copy("createO2tables.C", in_path=aod_path)
+    do_copy("muonAccEffPID.root", in_path=aod_path)
     if qa:
-        do_copy("diagnostic_tools/dpl-config_std.json", ".")
+        do_copy("diagnostic_tools/dpl-config_std.json")
 
     def set_config(config_file, config, value):
         config = config.strip()
@@ -355,8 +355,10 @@ def copy_and_link(file_name):
             write_to_runner(f"root -l -b -q 'createO2tables.C+(\"{delphes_file}\", \"tmp_{aod_file}\", 0)'",
                             log_file=aod_log_file,
                             check_status=True)
-            write_to_runner(f"mv tmp_{aod_file} {aod_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")
+            # 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:
                 copy_and_link(delphes_file)
                 if hepmc_file is not None:

examples/scripts/create_luts.sh

@@ -1,5 +1,13 @@
 #! /usr/bin/env bash
 
+if [[ $* == *"-h"* ]]; then
+    echo "Script to generate LUTs from LUT writer, arguments:"
+    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]"
+    exit 0
+fi
+
 WHAT=default
 FIELD=0.5
 RMIN=100.
@@ -8,37 +16,46 @@ RMIN=100.
 [ -z "$2" ] || FIELD=$2
 [ -z "$3" ] || RMIN=$3
 
-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 .
-
-echo " --- creating LUTs: config = $WHAT, field = $FIELD T, min tracking radius = $RMIN cm" 
-
-prepare_libs_macro() {
-    echo "void create_libs() {"
-    echo "   gROOT->ProcessLine(\".L DetectorK.cxx+\");"
-    echo "}"
-}
-
-prepare_luts_macro() {
-
-    echo "R__LOAD_LIBRARY(DetectorK_cxx.so)"
-    echo "#include \"lutWrite.$WHAT.cc\""
-    echo "void create_luts() {"
-    echo "   lutWrite_$WHAT(\"lutCovm.el.dat\",   11, $FIELD, $RMIN);"
-    echo "   lutWrite_$WHAT(\"lutCovm.mu.dat\",   13, $FIELD, $RMIN);"
-    echo "   lutWrite_$WHAT(\"lutCovm.pi.dat\",  211, $FIELD, $RMIN);"
-    echo "   lutWrite_$WHAT(\"lutCovm.ka.dat\",  321, $FIELD, $RMIN);"
-    echo "   lutWrite_$WHAT(\"lutCovm.pr.dat\", 2212, $FIELD, $RMIN);"
-    echo "}"
-
-}
-
-prepare_libs_macro > create_libs.C
-root -b -q -l create_libs.C
-
-prepare_luts_macro > create_luts.C
-root -b -q -l create_luts.C
-
+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" .
+
+echo " --- creating LUTs: config = $WHAT, field = $FIELD T, min tracking radius = $RMIN cm"
+
+for i in 0 1 2 3 4; do
+    root -l -b <<EOF
+    .L DetectorK.cxx+
+    .L lutWrite.${WHAT}.cc
+
+    TDatabasePDG::Instance()->AddParticle("deuteron", "deuteron", 1.8756134, kTRUE, 0.0, 3, "Nucleus", 1000010020);
+    TDatabasePDG::Instance()->AddAntiParticle("anti-deuteron", -1000010020);
+
+    TDatabasePDG::Instance()->AddParticle("helium3", "helium3", 2.80839160743, kTRUE, 0.0, 6, "Nucleus", 1000020030);
+    TDatabasePDG::Instance()->AddAntiParticle("anti-helium3", -1000020030);
+
+    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 int i = ${i};
+    lutWrite_${WHAT}("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
+        echo "${i} has zero size"
+        NullSize="${NullSize} ${i}"
+    fi
+done
+
+if [[ ! -z $NullSize ]]; then
+    echo "Created null sized LUTs!!"
+    exit 1
+fi

src/lutWrite.cc

@@ -60,9 +60,9 @@ lutWrite(const char *filename = "lutCovm.dat", int pdg = 211, float field = 0.2,
   // pid
   lutHeader.pdg = pdg;
   lutHeader.mass = TDatabasePDG::Instance()->GetParticle(pdg)->Mass();
-  const int q = TDatabasePDG::Instance()->GetParticle(pdg)->Charge() / 3;
+  const int q = std::abs(TDatabasePDG::Instance()->GetParticle(pdg)->Charge()) / 3;
   if (q <= 0) {
-    Printf("Negative or null charge for %i. Fix the charge!", pdg);
+    Printf("Negative or null charge (%f) for pdg code %i. Fix the charge!", TDatabasePDG::Instance()->GetParticle(pdg)->Charge(), pdg);
     return;
   }
   lutHeader.field = field;