Extend AOD creation, extend analysis tools, Improve FTOF writing to AODs (#71)

* Computing T0 before filling fTOF tree with all FTOF tracks

* Add more arguments to the analysis runner

* Add possibility to append to an existing production

Author: njacazio

examples/aod/createO2tables.C

@@ -206,7 +206,10 @@ int createO2tables(const char* inputFile = "delphes.root",
     // loop over tracks
     std::vector<TrackAlice3> tracks_for_vertexing;
     std::vector<o2::InteractionRecord> bcData;
+    // Tracks used for the T0 evaluation
     std::vector<Track*> tof_tracks;
+    std::vector<Track*> ftof_tracks;
+    std::vector<std::pair<int, int>> ftof_tracks_indices;
     const int multiplicity = tracks->GetEntries();
 
     // Build index array of tracks to randomize track writing order
@@ -321,25 +324,8 @@ int createO2tables(const char* inputFile = "delphes.root",
 
       // check if has Forward TOF
       if (forward_tof_layer.hasTOF(*track)) {
-        ftof.fIndexCollisions = ientry + eventOffset;
-        ftof.fIndexTracks = fTrackCounter; // Index in the Track table
-
-        ftof.fFTOFLength = track->L * 0.1;        // [cm]
-        ftof.fFTOFSignal = track->TOuter * 1.e12; // [ps]
-
-        std::array<float, 5> deltat, nsigma;
-        forward_tof_layer.makePID(*track, deltat, nsigma);
-        ftof.fFTOFDeltaEl = deltat[0];
-        ftof.fFTOFDeltaMu = deltat[1];
-        ftof.fFTOFDeltaPi = deltat[2];
-        ftof.fFTOFDeltaKa = deltat[3];
-        ftof.fFTOFDeltaPr = deltat[4];
-        ftof.fFTOFNsigmaEl = nsigma[0];
-        ftof.fFTOFNsigmaMu = nsigma[1];
-        ftof.fFTOFNsigmaPi = nsigma[2];
-        ftof.fFTOFNsigmaKa = nsigma[3];
-        ftof.fFTOFNsigmaPr = nsigma[4];
-        FillTree(kFTOF);
+        ftof_tracks.push_back(track);
+        ftof_tracks_indices.push_back(std::pair<int, int>{ientry + eventOffset, fTrackCounter});
       }
 
       // check if it is within the acceptance of the MID
@@ -362,6 +348,34 @@ int createO2tables(const char* inputFile = "delphes.root",
       fTrackCounter++;
       // fill histograms
     }
+
+    // Filling the fTOF tree after computing its T0
+    std::array<float, 2> ftzero;
+
+    forward_tof_layer.eventTime(ftof_tracks, ftzero);
+    for (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
+
+      ftof.fFTOFLength = track->L * 0.1;        // [cm]
+      ftof.fFTOFSignal = track->TOuter * 1.e12; // [ps]
+
+      std::array<float, 5> deltat, nsigma;
+      forward_tof_layer.makePID(*track, deltat, nsigma);
+      ftof.fFTOFDeltaEl = deltat[0];
+      ftof.fFTOFDeltaMu = deltat[1];
+      ftof.fFTOFDeltaPi = deltat[2];
+      ftof.fFTOFDeltaKa = deltat[3];
+      ftof.fFTOFDeltaPr = deltat[4];
+      ftof.fFTOFNsigmaEl = nsigma[0];
+      ftof.fFTOFNsigmaMu = nsigma[1];
+      ftof.fFTOFNsigmaPi = nsigma[2];
+      ftof.fFTOFNsigmaKa = nsigma[3];
+      ftof.fFTOFNsigmaPr = nsigma[4];
+      FillTree(kFTOF);
+    }
+
     if (eventextra.fNentries[kTracks] != eventextra.fNentries[kTracksCov] || eventextra.fNentries[kTracks] != eventextra.fNentries[kTracksExtra]) {
       Printf("Issue with the counters");
       return 1;

examples/scripts/createO2tables.py

@@ -98,26 +98,19 @@ def main(configuration_file,
          output_path,
          clean_delphes_files,
          create_luts,
-         turn_off_vertexing):
+         turn_off_vertexing,
+         append_production):
+    arguments = locals()
     global verbose_mode
     verbose_mode = verbose
     parser = configparser.RawConfigParser()
     parser.read(configuration_file)
 
     run_cmd("./clean.sh > /dev/null 2>&1", check_status=False)
     # Dictionary of fetched options
-    running_options = {
-        "ARG configuration_file": configuration_file,
-        "ARG config_entry": config_entry,
-        "ARG njobs": njobs,
-        "ARG nruns": nruns,
-        "ARG nevents": nevents,
-        "ARG verbose": verbose,
-        "ARG qa": qa,
-        "ARG output_path": output_path,
-        "ARG clean_delphes_files": clean_delphes_files,
-        "ARG create_luts": create_luts
-    }
+    running_options = {}
+    for i in arguments:
+        running_options["ARG "+i] = arguments[i]
 
     def opt(entry, require=True):
         try:
@@ -211,18 +204,19 @@ def check_duplicate(option_name):
 
     # Printing configuration
     msg(" --- running createO2tables.py", color=bcolors.HEADER)
-    msg("  njobs    =", njobs)
-    msg("  nruns    =", nruns)
-    msg("  nevents  =", nevents)
-    msg("  lut path =", lut_path)
+    msg("  n. jobs        =", njobs)
+    msg("  n. runs        =", nruns)
+    msg("  events per run =", nevents)
+    msg("  tot. events    =", "{:.0e}".format(nevents*nruns))
+    msg("  LUT path       =", f"'{lut_path}'")
     msg(" --- with detector configuration", color=bcolors.HEADER)
-    msg("  bField               =", bField, "[kG]")
+    msg("  B field              =", bField, "[kG]")
     msg("  sigmaT               =", sigmaT, "[ns]")
     msg("  sigmaT0              =", sigmaT0, "[ns]")
-    msg("  barrel_radius        =", barrel_radius, "[cm]")
-    msg("  barrel_half_length   =", barrel_half_length, "[cm]")
+    msg("  Barrel radius        =", barrel_radius, "[cm]")
+    msg("  Barrel half length   =", barrel_half_length, "[cm]")
     if create_luts:
-        msg("  minimum_track_radius =", minimum_track_radius, "[cm]")
+        msg("  Minimum track radius =", minimum_track_radius, "[cm]")
     msg("  LUT                  =", lut_tag)
     msg("  etaMax               =", etaMax)
 
@@ -283,8 +277,20 @@ def set_config(config_file, config, value):
                "const double tof_sigmat =", f"{sigmaT}""\;/")
     set_config("createO2tables.C",
                "const double tof_sigmat0 =", f"{sigmaT0}""\;/")
-
     run_list = range(nruns)
+    if append_production:
+        if output_path is None:
+            fatal_msg("Output path is not defined, cannot append")
+        last_preexisting_aod = [each for each in os.listdir(output_path)
+                                if each.endswith('.root') and "AODRun5" in each]
+        if len(last_preexisting_aod) == 0:
+            fatal_msg("Appending to a non existing production")
+        last_preexisting_aod = sorted([int(each.replace("AODRun5.", "").replace(".root", ""))
+                                       for each in last_preexisting_aod])[-1] + 1
+        msg(f" Appending to production with {last_preexisting_aod} AODs",
+            color=bcolors.BWARNING)
+        run_list = range(last_preexisting_aod,
+                         last_preexisting_aod + nruns)
 
     def configure_run(run_number):
         # Create executable that runs Generation, Delphes and analysis
@@ -422,11 +428,28 @@ def write_config(entry, prefix=""):
             if "ARG" not in i:
                 write_config(i, prefix=" * ")
 
+        output_size = sum(os.path.getsize(os.path.join(output_path, f))
+                          for f in os.listdir(output_path)
+                          if os.path.isfile(os.path.join(output_path, f)))
+        f.write("\n##  Size of the ouput ##\n")
+        f.write(f" - {output_size} bytes\n")
+        f.write(f" - {output_size/1e6} MB\n")
+        f.write(f" - {output_size/1e9} GB\n")
     run_cmd("echo  >> " + summaryfile)
     run_cmd("echo + DelphesO2 Version + >> " + summaryfile)
     run_cmd("git rev-parse HEAD >> " + summaryfile, check_status=False)
+
     if os.path.normpath(output_path) != os.getcwd():
-        run_cmd(f"mv {summaryfile} {output_path}")
+        if append_production:
+            s = os.path.join(output_path, summaryfile)
+            run_cmd(f"echo '' >> {s}")
+            run_cmd(f"echo '  **' >> {s}")
+            run_cmd(f"echo 'Appended production' >> {s}")
+            run_cmd(f"echo '  **' >> {s}")
+            run_cmd(f"echo '' >> {s}")
+            run_cmd(f"cat {summaryfile} >> {s}")
+        else:
+            run_cmd(f"mv {summaryfile} {output_path}")
 
     if qa:
         msg(" --- running test analysis", color=bcolors.HEADER)
@@ -441,7 +464,7 @@ def write_config(entry, prefix=""):
     parser.add_argument("--entry", "-e", type=str,
                         default="DEFAULT",
                         help="Entry in the configuration file, e.g. the INEL or CCBAR entries in the configuration file.")
-    parser.add_argument("--output-path", "-o", type=str,
+    parser.add_argument("--output-path", "--output_path", "-o", type=str,
                         default=None,
                         help="Output path, by default the current path is used as output.")
     parser.add_argument("--njobs", "-j", type=int,
@@ -463,10 +486,18 @@ def write_config(entry, prefix=""):
     parser.add_argument("--no-vertexing",
                         action="store_true",
                         help="Option turning off the vertexing.")
+    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("--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.")
     args = parser.parse_args()
+    # Check arguments
+    if args.append and args.output_path is None:
+        fatal_msg(
+            "Asked to append production but did not specify output path (option '-o')")
+
     main(configuration_file=args.configuration_file,
          config_entry=args.entry,
          njobs=args.njobs,
@@ -477,4 +508,5 @@ def write_config(entry, prefix=""):
          clean_delphes_files=args.clean_delphes,
          qa=args.qa,
          create_luts=not args.use_preexisting_luts,
-         turn_off_vertexing=args.no_vertexing)
+         turn_off_vertexing=args.no_vertexing,
+         append_production=args.append)

examples/scripts/diagnostic_tools/doanalysis.py

@@ -168,16 +168,24 @@ def main(mode,
          dpl_configuration_file=None,
          njobs=1,
          merge_output=False,
-         only_merge=False):
+         merge_only=False,
+         shm_mem_size=16000000000,
+         readers=1,
+         extra_arguments=""):
     if len(input_file) == 1:
         input_file = input_file[0]
     else:
         input_file = input_file[0:n_max_files]
-    msg("Running", f"'{mode}'", "analysis on",
-        f"'{input_file}'", color=bcolors.BOKBLUE)
-    msg("Maximum", n_max_files, "files with batch size",
-        batch_size, "and", njobs, "jobs" if njobs > 1 else "job", color=bcolors.BOKBLUE)
-    args = f"-b --shm-segment-size 16000000000 --readers 4"
+    if not merge_only:
+        msg("Running", f"'{mode}'", "analysis on",
+            f"'{input_file}'", color=bcolors.BOKBLUE)
+        msg("Maximum", n_max_files, "files with batch size",
+            batch_size, "and", njobs, "jobs" if njobs > 1 else "job", color=bcolors.BOKBLUE)
+    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 += extra_arguments
     if mode not in analyses:
         raise ValueError("Did not find analyses matching mode",
                          mode, ", please choose in", ", ".join(analyses.keys()))
@@ -229,16 +237,16 @@ def build_list_of_files(file_list):
     run_list = []
     for i, j in enumerate(input_file_list):
         run_list.append(set_o2_analysis(an,
-                                        o2_arguments=args,
+                                        o2_arguments=o2_arguments,
                                         input_file=j,
                                         tag=tag,
                                         dpl_configuration_file=dpl_configuration_file))
-    if not only_merge:
+    if not merge_only:
         with multiprocessing.Pool(processes=njobs) as pool:
             pool.map(run_o2_analysis, run_list)
-        msg("Analysis completed", color=bcolors.BOKGREEN)
 
-    if merge_output or only_merge:
+    if merge_output or merge_only:
+        msg("Merging results", color=bcolors.BOKBLUE)
         files_to_merge = []
         for i in input_file_list:
             p = os.path.dirname(os.path.abspath(i))
@@ -304,10 +312,17 @@ def build_list_of_files(file_list):
                         help="Name of the dpl configuration file e.g. dpl-config_std.json")
     parser.add_argument("--merge_output", "--merge-output", "--merge",
                         action="store_true", help="Flag to merge the output files into one")
+    parser.add_argument("--readers", "-r",
+                        default=1, type=int,
+                        help="Number of parallel readers")
+    parser.add_argument("--mem", "-m",
+                        default=16000000000, type=int,
+                        help="Size of the shared memory to allocate")
+    parser.add_argument("--extra_arguments", "-e",
+                        default="", type=str,
+                        help="Extra arguments to feed to the workflow")
     parser.add_argument("--merge_only", "--merge-only", "--mergeonly",
                         action="store_true", help="Flag avoid running the analysis and to merge the output files into one")
-    parser.add_argument("-b",
-                        action="store_true", help="Background mode")
     args = parser.parse_args()
     if args.verbose:
         verbose_mode = False,
@@ -322,4 +337,7 @@ def build_list_of_files(file_list):
              out_tag=args.tag,
              merge_output=args.merge_output,
              out_path=args.out_path,
-             only_merge=args.merge_only)
+             merge_only=args.merge_only,
+             readers=args.readers,
+             extra_arguments=args.extra_arguments,
+             shm_mem_size=args.mem)