#!/usr/bin/env python """ Usage any_to_aa.py -f input_file -o output_file Produce a KM3NeT summary file. The summary files use the offline, with hits removed, and only the best tracks for each reconstruction algorihtm. Additional summary information is recorded in a friend TTree. The Monte-carlo header informtion is checked for consistency and merged. Aliases are added to the main TTree for easy access to often-used variables and to contants from the header options: -f : input file default=../../data/mc5.1.numuCC.ps10.aa.root -o : output file default=dst.root -d : detector files default=get_from_header -p : print 1-line summaries of the events default=False -h : this help message and exit default=False -i : Python interative mode (prompt when done) default=False -r : ids of reco algorithms default=jppmuon,aashower -H : print header of output file default=False -m : record the 'index' from filename as mc_run_id default=True -v : report every default=100 -s : sample/prescale: process every nth event default=1 -O : oscprob directory default= -C : store cerenkov hits default=False -T : apply time slewing correction default=False -n : process this many input events (before prescale) default = -1 -N2020 : store neutrino2020 variables default=False -RBR : Run by run approach, info for DB query default=ORCA,-1 -RBR_livetime : RBR livetime, bypass DB query default=-1.0 -casc : store cascade variables default=False -keepMC : keep the full mc_trks information default=False """ import os, sys, aa sys.path += [ aa.aadir + "/ana/dst"] from aliases import * from dst_util import * from math import * from copy import copy from ROOT import EventFile, best_reco_track, Trk, Timer import ROOT options = aa.Options( __doc__, sys.argv[1:] ) infiles, outfile = options.f.strip() , options.o # if multiple input files, put them in list if " " in infiles : infiles = infiles.split() elif "," in infiles : infiles = infiles.split(",") else : infiles = [infiles] # make list of rectypes to put in output file rec_types = list( options.r.split(",") ) for rec in rec_types: if rec in ['jppmuon','jppshower','aashower']: continue raise Exception(f"Unknown reconstruction type: {rec}") #===================================================================== # we are going to use oscprob... todo: allow running without oscprob oscprobdir = options.O if oscprobdir == "": print ("OscProb path not provided, look for environment variable OSCPROBDIR.") if "OSCPROBDIR" in os.environ: oscprobdir = os.environ["OSCPROBDIR"] print ("OSCPROBDIR found.") else: raise Exception("OSCPROBDIR not found. Use either OSCPROBDIR environment variable or set it manually using the -O command line argument") else: print (f"OscProb directory set with command line argument: {oscprobdir}") ROOT.gInterpreter.ProcessLine(".L "+oscprobdir + "/libOscProb.so" ) os.chdir(aa.aadir + "/ana/dst") aa.include("summarize_hits.hh") aa.include("summarize_mc_trks.hh") aa.include("summarize_rec_trks.hh") aa.include("summarize_jppshower.hh") aa.include("summarize_crkv_hits.hh") aa.include("celestial_coords.hh") aa.include("nu_summary.hh") aa.include("summarize_mc_evts.hh") aa.include("time_slewing.hh") aa.include("features_ORCA_Neutrino2020.hh") aa.include("cascade_summary.hh") # ---- definitions from km3net-dataformat sys.path += [ aa.aadir + "/externals/km3net-dataformat/definitions"] import w2list_gseagen, trigger, fitparameters # do the flavour - challenge if len ( { flavour(_) for _ in infiles } ) > 1 : print ("ERROR: make_dst.py called with input files with >1 flavours:") print ( "\n".join(infiles) ) sys.exit() flav = flavour( infiles[0] ) is_mc = False is_cc = None is_mu = None is_nu = None is_noise = None if flav == 'data' : pass else : is_mc = True if flav is None: is_noise = True flav = 'pure-noise' is_mu = 'mupage' in flav is_nu = not is_mu and not is_noise is_cc = ( flav.endswith("CC") or flav.endswith("CCHEDIS") ) print (" flav ", flav ) print (" is_mc ", is_mc ) print (" is_mu ", is_mu ) print (" is_nu ", is_nu ) print (" is_cc ", is_cc ) print (" is_noise ", is_noise ) def update_detector( evtfile , verbose = False ) : global det if evtfile.current_filename == update_detector.filename : return update_detector.filename = evtfile.current_filename run = evtfile.evt.run_id print ("current run ", run) if options.d == "get_from_header" : det = get_detector( f.current_file_header ) if not det : print (" failed to get detector from the event-file header") import sys sys.exit() else : detfile = options.d.replace( "RUN",str(run) ) if hasattr( det, 'filename' ) and detfile == det.filename : return else: print ("new detector :", detfile) ok = det.read( detfile ) if not ok : print ("failed to load detector") import sys sys.exit(44) det.filename = detfile update_detector.filename = None #check we can open and read all the files for fn in infiles[:] : try : f = EventFile( fn ) except : print ("File open error", sys.exc_info()[0] ) print ("BADFILE : ", fn) infiles.remove( fn ) #========================================================================== f = EventFile( infiles ) det = ROOT.Det() f.set_output( outfile ) # this opens an aanet output file # 'meta' info and how the file was produced f.autosave = False f.output_file.mkdir("dst_history") f.output_file.cd("dst_history") write_as_tnamed("input_files", " ".join(infiles) ) write_as_tnamed("command_line", " ".join(sys.argv) ) # Compute the exposure # For pure MC files, exposure is set to 1 year # For data or RBR MC files, exposure is set to the # the run ones. 3 solutions: # - Manually provided livetime # - Provided by a DB call # - Taken from header (default, JPP > 17.0.0) DAQ_livetime = 365.25*24.*3600. detector, run_number = options.RBR.split(',') run_number = int(run_number) RBR_info = None if options.RBR_livetime > 0.: # Manually provided livetime DAQ_livetime = options.RBR_livetime RBR_info = {'livetime_s':DAQ_livetime, 'runNumber':run_number} print ('RBR approach, manually set livetime. Call to DB being bypassed.') print (f' \tlivetime: {DAQ_livetime}s') elif run_number > 0: # DB request print ('RBR approach, connect to DB to retrieve livetime_s (and more ...)') print (' \tDetector: {}, run number: {}'.format(detector, run_number)) RBR_info = getRBRinfo(detector, run_number) DAQ_livetime = RBR_info['livetime_s'] print ('Done') elif "DAQ" in dict(f.header): # Use header value h = dict(f.header) DAQ_livetime = float(h['DAQ']) RBR_info = {"livetime_s" : DAQ_livetime} f.roottree().GetEntry(0) RBR_info['runNumber'] = int(f.roottree().GetLeaf("run_id").GetValue()) else : # It's a pure MC file RBR_info = {'runNumber':1, 'livetime_s':DAQ_livetime} print ('Weight be expressed by year') print ('DAQ livetime : {0} s'.format(DAQ_livetime)) # summary tree if RBR_info is not None: sumTree = summary_tree(f) for key, value in RBR_info.items() : sumTree.addBranch(key, value) sumTree.close() # Add another tree as "friend" f.output_file.cd() T = ROOT.TTree("T","summary variables") mc_trks_summary = ROOT.MC_trks_summary() nu_summary = ROOT.nu_summary() hits_summary = ROOT.hits_summary() trig_hits_summary = ROOT.hits_summary() mc_hits_summary = ROOT.hits_summary() cel_coords = ROOT.Celestial_coordinates() crkv_hits = ROOT.vector('crkv_hits')() hand_Neutrino2020 = ROOT.features_ORCA_Neutrino2020() feat_Neutrino2020 = ROOT.ORCA_Neutrino2020() mc_evt_summary = ROOT.MC_evts_summary() cascade_summary = ROOT.cascade_summary() jppshower_summary = ROOT.jppshower_summary() jppmuon_summary = ROOT.rec_trks_summary() if is_mc : if is_nu or is_mu: T.Branch( "sum_mc_trks" , mc_trks_summary , 32000, 4 ) T.Branch( "sum_mc_hits" , mc_hits_summary , 32000, 4 ) if is_nu: T.Branch( "sum_mc_nu" , nu_summary , 32000, 4 ) T.Branch( "sum_mc_evt" , mc_evt_summary , 32000, 4 ) T.Branch( "sum_hits" , hits_summary , 32000, 4 ) T.Branch( "sum_trig_hits" , trig_hits_summary , 32000, 4 ) T.Branch( "coords" , cel_coords , 32000, 4 ) if options.C : T.Branch( "crkv_hits" , crkv_hits , 32000, 4 ) if options.N2020 : T.Branch( "feat_Neutrino2020" ,feat_Neutrino2020 , 32000, 4 ) if options.casc : T.Branch( "sum_casc" ,cascade_summary , 32000, 4 ) if "jppmuon" in rec_types: T.Branch("sum_jppmuon", jppmuon_summary, 32000, 4) if "jppshower" in rec_types: T.Branch("sum_jppshower", jppshower_summary, 32000, 4) f.output_tree.AddFriend(T) write_header_dir( f.header, f.output_file ) E = f.output_tree if is_mc : deepAlias(E, ROOT.Trk, "nu", "mc_trks[0]") deepAlias(E, ROOT.Trk, "mu", "mc_trks[1]") E.SetAlias("angle_nu_mu", angle("nu","mu") ) E.SetAlias("angle_nu_trackfit", angle("nu","trackfit" , True ) ) E.SetAlias("angle_nu_showerfit", angle("nu","showerfit", True ) ) E.SetAlias("angle_mu_trackfit", angle("mu","trackfit" , True ) ) E.SetAlias("angle_mu_showerfit", angle("mu","showerfit", True ) ) deepAlias(E, ROOT.Trk, "trackfit", "trks[0]") deepAlias(E, ROOT.Trk, "showerfit", "trks[1]") if is_mc : define_vector_aliases( E, "w2list", w2list_gseagen.data ) define_vector_aliases( E, "trackfit.fitinf", fitparameters.data ) define_bits_aliases( E, trigger.data ) if is_mc : # initialize oscoprob (in summarize_mc_trks.hh) dm21, dm32, th12, th23, th13, deltacp = 7.400e-5, 2.494e-3, 5.868e-1,8.238e-1,1.491e-1,234.*pi/180. ROOT.oscprob_pmns().SetDeltaMsqrs( dm21, dm32 ) ROOT.oscprob_pmns().SetMix( th12, th23, th13, deltacp ) for n in "dm21, dm32, th12, th23, th13, deltacp".split() : E.SetAlias("osc_"+n, '('+str(eval(n))+')' ) # make an alias for every header entry for the moment, only 'known' tags generate_header_aliases( f.header, f.output_tree ) NNN =1 maxhits, maxmchits = 1,1 print(rec_types) rectypesv = ROOT.vector(str)() for x in rec_types : rectypesv.push_back( x ) for evt in f : with Timer("event") : # prescale (use global index to correctly prescale many small input files) if f.global_index() % options.s > 0 : continue if options.n > 0 and f.global_index() >= options.n : print("requested number of events (",options.n,") reached") q = f.global_index() / float( f.total_events ) print ('extrapolated total time:', f.wallwatch.RealTime() / q ) break if options.p : print( f.global_index() , f.index , evt) if options.m and is_mc : # record the 'index' number found in the filename as mc_run_id set_mc_run_id_to_file_index( f ) # ensure we have a valid (run-dependent) det for the file we are processing update_detector( f ) # record summary information / derived quantities in out friend Tree if (f.index % options.v == 0 ) : Timer.report() with Timer("summarize") : def call ( func, *args ) : with Timer(func.__name__ ) : ok = func( *args ) if not ok : print ( func, "failed" ) print ( "file ", f.current_filename() ) print ("------------------------------------------------------------") import sys sys.exit(1) if is_mc : if is_nu: call( ROOT.summarize_nu , evt, nu_summary ) call ( ROOT.summarize_mc_trks , evt, mc_trks_summary , is_nu, is_cc, options.keepMC) call ( ROOT.summarize_hits , det, evt, mc_hits_summary , "mc_hits" ) if DAQ_livetime > 0 : call ( ROOT.summarize_mc_evt , is_nu, is_mu, mc_evt_summary, evt, f.header, nu_summary, DAQ_livetime ) if options.T : det.apply(evt) call ( ROOT.apply_time_slewing, evt ) if options.N2020 : call (hand_Neutrino2020.fill_Neutrino2020 , evt, det , feat_Neutrino2020 ) call ( ROOT.summarize_hits , det, evt, hits_summary , "hits" ) call ( ROOT.summarize_hits , det, evt, trig_hits_summary , "trig_hits" ) if "jppshower" in rec_types: call ( ROOT.summarize_jppshower, evt, jppshower_summary) if "jppmuon" in rec_types: call ( ROOT.summarize_rec_trks, evt, jppmuon_summary, ROOT.JPP_RECONSTRUCTION_TYPE , ROOT.JMUONGANDALF) call ( ROOT.keep_best_trks, evt, rectypesv ) call ( ROOT.fill_celestial_coordinates, evt, det , cel_coords ) if options.C : crkv_hits.clear() for trk in evt.trks : ch = ROOT.crkv_hits() call ( ROOT.summarize_crkv_hits, evt, trk, ch ) crkv_hits.push_back(ch) if options.casc: call ( ROOT.summarize_cascade, evt, cascade_summary ) # this clears the hits and mc hits and removes the hit-refs from the tracks evt.clear_hits() T.Fill() f.write() hitsbranch = T.GetBranch("hits") T.GetListOfBranches().Remove( hitsbranch ) print("-------------------------------") v= outfile.rsplit("/",1) if len(v) == 1 : outdir = "./www/" else : outdir = v[0] # for safety, write the file (eventhout eventfile dtor should do it too) f.output_file.Write() del f