#!/usr/bin/env python """ This is a complete analysis of K40 signals recorded in the first ORCA line. We use c++ for the hit-processing, leading to many histograms of the delta-t distibutions for many different pairs of PMTs. The c++ code is compiled on-the-fly by the cxx function. The estimated random-coincidence rate is computed from the SummarySlices. Next, the offsets, efficiencies and resolutions are fitted to the histograms using Minuit. The histograms, the result of the fit, and the fitted functions are written to html file for inspection in a web browser. The default run (2711) is special since the trigger used a 20ns-wide coindidence window (rather than the more usual 10ns). options: -w : www dir, default=/afs/in2p3.fr/home/h/heijboer/www/k40 -f : input file, default=/in2p3/km3net/data/raw/sea/KM3NeT_00000029/2/KM3NeT_00000029_00002711.root -d : det file, default=/sps/km3net/users/zaborov/km3net/det/KM3NeT_00000029_20170908_DU2.detx -r : force re-filling of histograms default=False -h : this help message and exit default=False -i : Python interative mode (prompt when done) default=False """ import copy, collections, itertools, sys, os import aa import ROOT gStyle.SetOptStat(0); # aa.py contains helper to deal with cmd line options options = aa.Options( __doc__, sys.argv[1:] ) html_dir = options.w # Compile and load the c++ code for this analysis cxx ( open("k40_histfill.cc").read() ) os.system("mkdir "+ options.w ) det = Det( options.d ) # Create the histograms unless the file exists # or when 'refill' is on the cmd line histo_file_name = "histograms.root" fill_histograms = not os.path.exists( histo_file_name ) or options.r #-------------------------------------------------------------- # Part 1 : Open input file and process the hits, save to file #-------------------------------------------------------------- if fill_histograms : EventFile.read_timeslices = True # timeslices will be read into the evt.hits vector. f = EventFile( options.f ) for evt in f: # Load the rates from the summarytimeslice and set the pmt.rate fields # in the detector object: with Timer("summary_slices"): f.get_rates( det ) # Note that the c++ function process_coincidences2 returns the same, # updated, collection of histogreams each time. with Timer("hit processing") : histograms = process_coincidences( evt.hits , det ) if f.index%10000 == 0 : print f.index,"/",f.roottree().GetEntries(), " " , f.report() Timer.report() # Write the map object that holds the histogram to file TFile(histo_file_name,"RECREATE").WriteObject( histograms, "histograms") else : # load the histograms from root file histograms = MakeNullPointer( "map" ) TFile( histo_file_name ).GetObject( "histograms", histograms ) #------------------------------------------------------------------------------------------- # part 2: Post-process and the histograms #------------------------------------------------------------------------------------------- def k40_rate( cos_angle ): "function derived from MC, describing the rate(Hz) of k40 coincidences vs the angle" coefficients = (-1.07061, 3.17173, -1.35769, 1.6885) return exp( sum( coefficients[i] * (cos_angle**i) for i in range(4) ) ) F = TF1("F","gaus(0)+[3]",-20,20); F.SetParLimits( 0, 0, 10000 ) # gaus norm F.SetParLimits( 1, -5, 5 ) # mean F.SetParLimits( 2, 1, 10 ) # sigma F.SetParameter( 1, 0 ) # mean = 0 F.SetParameter( 2, 3 ) # sigma = 3 histos = collections.defaultdict( list ) for domid, vector_of_hists in histograms : for i,j in itertools.combinations( range(31), 2 ) : h = vector_of_hists[i*31+j] h.domid = domid # on the python-side, we can just add attributes (1) h.ch1, h.ch2 = i,j h.cos_angle = det.doms[domid].pmts[i].dir.dot( det.doms[domid].pmts[j].dir ) if h.cos_angle < 0 : continue # >90 deg was not triggered on in this run h.random_bg = h.GetBinContent(0) h.nslices = h.GetBinContent( h.GetNbinsX()+1 ) h.nk40 = h.nslices * 0.10 * k40_rate( h.cos_angle ) # expected k40 rate h.SetTitle("%d %d %d angle = %f" % ( domid, i,j, 57.2 * acos ( h.cos_angle ) ) ) # make a copy of F, normalisation so that integral is nk40 h.F = copy.copy(F) h.F.SetParameter( 0, h.GetBinWidth(0) * h.nk40 / ( (sqrt(2*pi) * F.GetParameter(2)) ) ) h.F.FixParameter( 3, h.random_bg ) h.GetListOfFunctions().Add( h.F ) # so that it's plotted automatically h.SetLineWidth(2) h.SetFillColor(5) # We have to store the h objects on the python side now # because of (1), so a new list histos[domid].append( h ) def write_html( domid, results_table, histos, height=100, width=100, imgdir = "imgs" ): oldbatch = gROOT.IsBatch(); gROOT.SetBatch( True ) os.system("mkdir -p "+html_dir +"/"+imgdir ); f = open ( html_dir + "/" +str(domid) + ".html" ,'w') f.write ("\n") f.write ("

Results for DOM "+str(domid)+"

\n\n") f.write( results_table.html() ) f.write( "
\n\n"*2 ) c = TCanvas("c","c",500,500) for h in histos : fn = imgdir + "/" + h.GetName()+".png" h.Draw() latex = TLatex() latex.SetNDC(True); latex.SetTextSize( 0.1 ); latex.DrawLatex( 0.2 , 0.8 , str(h.ch1)+"-"+str(h.ch2)); c.SaveAs( html_dir+"/"+fn ) f.write(' \n') f.write("") gROOT.SetBatch ( oldbatch ) #------------------------------------------------------------------------------------------- # part 3: 93-parameter Fit the time offsets and other variables #------------------------------------------------------------------------------------------- def score_function ( histograms, offsets, efficiencies, resolutions ) : for h in histograms : F = h.F F.SetParameter( 1, offsets[h.ch1]-offsets[h.ch2] ) F.SetParameter( 2, (resolutions[h.ch1]**2 + resolutions[h.ch2] ) ** 0.5 ) F.SetParameter( 0, h.nk40 / ( sqrt(2*pi) * F.GetParameter(2) ) *\ efficiencies[h.ch1] * efficiencies[h.ch2] ) F.SetParameter( 3, h.random_bg ) return sum( h.Chisquare( h.F ) for h in histograms ) def fit( histograms ) : # wrap the score_function into something that minuit can call def fcn( npar, gin, f, par__, iflag ): par = [par__[i] for i in range(npar[0])] # the par__ object does support slicing, so make alist f[0] = score_function( histograms, par[0:31], par[31:62], par[62:93] ) pars = "offset efficiency resolution".split() # define the start-value, step_size and min and max for each type of parameter descmap = { 'offset' : (0, 0.1, -5, 5 ), 'efficiency' : (1, 0.1, 0.5, 1.5), 'resolution' : (3, 0.1, 2, 5 ) } M = TMinuit( 93 ) M.SetFCN( fcn ) ierflg = Long() for ivar, vartype in enumerate( pars ) : start, step,min_,max_ = descmap[vartype] for ipmt in range(31) : M.mnparm( ivar*31+ipmt , vartype+"_of_pmt"+str(ipmt) , start, step,min_,max_, ierflg ) # call minuit maxcalls, tollerance = 1000000, 1e-6 import array M.mnexcm( "MIGRAD", array.array( 'd', [maxcalls, tollerance] ) , 2 , ierflg ) # build a table with results T = Table(* ["pmt"]+pars ) for ipmt in range(31): T.add(ipmt) for ivar,vartype in enumerate( pars ): cv,err = double(), double() M.GetParameter( ipmt + ivar*31 , cv, err ) T.add("%5.3f +/- %5.3f" % ( float(cv), float(err) )) return T for domid in histos.keys() : histos_for_this_dom = histos[ domid ] results_table = fit ( histos_for_this_dom ) results_table.title = "fit for DOM " +str(domid) write_html( domid, results_table, histos_for_this_dom ) break # only one dom in this example