#!/usr/bin/python # import os,sys,numpy import os,sys nplanes, nwires = 44,80 offset_applied = 0. #to second histogram #offset_applied = -1.4 #to second histogram #offset_applied = -60. #to second histogram wire_cuts = (-1000,1000) #wire_cuts = (30,50) #wire_cuts = (25,75) stats_boxes = False # final_diff_plot_scale = 5. # final_diff_plot_scale = 3. final_diff_plot_scale = 1.0 # final_diff_plot_scale = 0.2 diff_1D_plot_scale = 5.0 ## __________ def check_usage(sa): if len(sa)!=3: print 'usage: compare_t0_files.py file1 file2' sys.exit(-1) for f in [sa[1],sa[2]]: if not os.path.isfile(f): print 'ERROR:',f,'does not exist. Exiting...' sys.exit(-1) return sa[1],sa[2] ## __________ def read_in_t0s(fn): f = file(fn,'r') fl = f.readlines() t0s = [] for i in range(nplanes): t0s.append([None]*nwires) for l in fl: if l.startswith('C') or l.startswith('#'): continue else: ls = l.split() ip, iw, t0, sigma, wtf = l.split() ip, iw = int(ip), int(iw) t0s[ip-1][iw-1] = float(t0) return t0s ## __________ def th2d_me(t0s,ti): true_min, true_max = 999999., -999999. h = TH2D(ti,ti, nplanes,0.5,float(nplanes+0.5), nwires,0.5,float(nwires+0.5)) for ip in range(nplanes): for iw in range(nwires): if ((iw+1)wire_cuts[1]): continue else: t0 = t0s[ip][iw] if t0>(-1000.): if t0true_max: true_max = t0 # print 'ip=',ip,'iw=',iw,'t0s=',t0s[ip][iw] h.Fill(ip+1,iw+1,t0) else: print 'DEAD WIRE: plane=',ip+1,'wire=',iw+1,' (skipping)' h.Fill(ip+1,iw+1,-999999.) return h, true_min, true_max ## __________ def make_histo_of_diffs(diff_t0s, file1, file2): h = TH1D('histo_of_diffs','Difference ('+file1+' - '+file2+');T0 [ns]',200,-1*diff_1D_plot_scale,diff_1D_plot_scale) for ip in range(nplanes): for iw in range(nwires): h.Fill(diff_t0s[ip][iw]) return h ## __________ def make_histo_of_diffs_US(diff_t0s): h = TH1D('histo_of_diffs_US','histo_of_diffs_US',200,-1.0*diff_1D_plot_scale,diff_1D_plot_scale) for ip in range(0,21): for iw in range(nwires): h.Fill(diff_t0s[ip][iw]) return h ## __________ def make_histo_of_diffs_DS(diff_t0s): h = TH1D('histo_of_diffs_DS','histo_of_diffs_DS',200,-1.0*diff_1D_plot_scale,diff_1D_plot_scale) for ip in range(22,44): for iw in range(nwires): h.Fill(diff_t0s[ip][iw]) return h ## __________ ## def mean_t0_me(t0s,offset=0.): ## if offset!=0.: print '*** mean_t0_me *** applying offset of',offset ## gr = TGraph() ## for ip in range(nplanes): ## alive_t0s = [] ## for iw in range(nwires): ## if t0s[ip][iw]>-1000.: ## if ((iw+1)wire_cuts[1]): ## continue ## else: ## alive_t0s.append(t0s[ip][iw]) ## gr.SetPoint(ip,ip+1,numpy.mean(alive_t0s)+offset) ## return gr ## __________ ##def th1d_me(l,ti): ## true_min, true_max = 999999., -999999. ## ## h = TH2D(ti,ti, ## nplanes,0.5,float(nplanes+0.5), ## nwires,0.5,float(nwires+0.5)) ## ## for ip in range(nplanes): ## for iw in range(nwires): ## t0 = t0s[ip][iw] ## if t0>(-1000.): ## if t0true_max: true_max = t0 ### print 'ip=',ip,'iw=',iw,'t0s=',t0s[ip][iw] ## h.Fill(ip+1,iw+1,t0) ## else: ## print 'DEAD WIRE: plane=',ip+1,'wire=',iw+1,' (skipping)' ## h.Fill(ip+1,iw+1,-999999.) ## return h, true_min, true_max ## __________ def nice_fn(ugly_fn): return ugly_fn.split('/')[-1] ## __________ def make_diff_t0s(l1,l2,offset=0.): ##NB offset is *added* to *second* file diff_t0s = [] for i in range(nplanes): diff_t0s.append([None]*nwires) for ip in range(nplanes): for iw in range(nwires): t0_1 = l1[ip][iw] t0_2 = l2[ip][iw] if t0_1>(-1000.) and t0_2>(-1000.): diff_t0s[ip][iw] = t0_1-(t0_2+offset) else: diff_t0s[ip][iw] = 0. return diff_t0s ## ____________________ from ROOT import * if __name__=='__main__': print 'checking usage...', f1, f2 = check_usage(sys.argv) print 'done' print 'reading in t0s from',f1,'and',f2,'...' t0_f1 = read_in_t0s(f1) t0_f2 = read_in_t0s(f2) print 'done' print 'importing ROOT ...' # ================ Display (Canvas, pads ...) =================== # # Canvas gStyle.SetCanvasColor ( 10 ) # Background color (white) # gStyle.SetCanvasDefH(500) # Default = 500 # gStyle.SetCanvasDefW(700) # Default = 700 # Output size # gStyle.SetPaperSize ( 20., 24. ) # Default = 20.,26. gStyle.SetPaperSize ( 12., 17. ) # Default = 20.,26. # Pads gStyle.SetPadBorderMode ( 0 ) gStyle.SetPadColor ( 10 ) # Background color (white) gStyle.SetPadLeftMargin (0.10) gStyle.SetPadBottomMargin (0.11) gStyle.SetPadRightMargin (0.16) gStyle.SetPadTopMargin (0.03) # Miscellaneous gStyle.SetFillColor ( 10 ) # Background color (white) gStyle.SetFrameFillColor ( 10 ) # Background color (white) # ================ Fonts variable ================= # Font = 132 # Times roman FontSize = 0.05 # ================ Plots =================== # gStyle.SetOptTitle ( 1 ) # Display histogram titles. gStyle.SetMarkerStyle ( 5 ) # Cross gStyle.SetPalette ( 1 ) # Nicer colour scale for 2D histograms. gStyle.SetHistLineWidth ( 2 ) # ================ Stat Box =================== # # Display name of the histogram, entries, mean, RMS, underflow, overflow # gStyle.SetOptStat ( 111111 ) # Good standard stats display. gStyle.SetOptStat ( 0 ) # Good standard stats display. gStyle.SetOptFit ( 1111 ) # Display fit info if available. gStyle.SetStatColor ( 10 ) # Background color (white) gStyle.SetStatX ( 0.95 ) # Statistic box X gStyle.SetStatY ( 0.9 ) # Statistic box Y gStyle.SetStatW ( 0.2 ) # Statistic box width # gStyle.SetStatH ( 0.2 ) # Statistic box height gStyle.SetStatBorderSize ( 1 ) gStyle.SetStatFont(Font) # No Stat box ! # gStyle.SetOptStat ( 0 ) # Good standard stats display. # gStyle.SetOptFit ( 0 ) # Display fit info if available. # ================ Title =================== # gStyle.SetTitleFillColor ( 10 ) # Background color (white) gStyle.SetTitleX ( 0.10 ) # Title box Y gStyle.SetTitleY ( 0.10 ) # Title box Y gStyle.SetTitleW ( 0.50 ) # Title box width gStyle.SetTitleH ( 0.09 ) # Title box height gStyle.SetTitleBorderSize ( 1 ) gStyle.SetTitleFont(Font) gStyle.SetTitleFont(Font, "X") gStyle.SetTitleFont(Font, "Y") gStyle.SetTitleFont(Font, "Z") gStyle.SetTitleFont(Font, "T") gStyle.SetTitleSize(FontSize ) gStyle.SetTitleSize(FontSize,"X") gStyle.SetTitleSize(FontSize,"Y") gStyle.SetTitleSize(FontSize,"Z") gStyle.SetTitleSize(FontSize,"T") # ================ Labels =================== # gStyle.SetLabelFont(Font) gStyle.SetLabelFont(Font, "X") gStyle.SetLabelFont(Font, "Y") gStyle.SetLabelFont(Font, "Z") gStyle.SetLabelFont(Font, "T") gStyle.SetLabelSize(FontSize) gStyle.SetLabelSize(FontSize,"X") gStyle.SetLabelSize(FontSize,"Y") gStyle.SetLabelSize(FontSize,"Z") gStyle.SetLabelSize(FontSize,"T") # =============== Text & Legends ================== # gStyle.SetTextFont(Font) gStyle.SetLegendBorderSize ( 1 ) # ============================================= # # Do not display the canvas in a window, plot in a postscript gROOT.SetBatch(); # Needed to make sure that all the gStyle parameters are passed to the plots gROOT.ForceStyle() print 'done' print 'making normal histograms ...' f1h,f1h_min,f1h_max = th2d_me(t0_f1,'t0s for '+nice_fn(f1)+';plane;wire') f2h,f2h_min,f2h_max = th2d_me(t0_f2,'t0s for '+nice_fn(f2)+';plane;wire') print 'done' print 'making diff histograms ...' if offset_applied!=0.: print '*** APPLYING OFFSET ***:',offset_applied,'to second histogram' diff_t0s = make_diff_t0s(t0_f1,t0_f2,offset_applied) # dh,dh_min,dh_max = th2d_me(diff_t0s, 'T0s difference ('+nice_fn(f2)+' - '+nice_fn(f1)+');plane;wire') dh,dh_min,dh_max = th2d_me(diff_t0s, ';plane number;wire number;Time offset difference [ns]') histo_diffs = make_histo_of_diffs(diff_t0s, nice_fn(f2), nice_fn(f1)) ## add separate diff histograms for US and DS hdiffs_US = make_histo_of_diffs_US(diff_t0s) hdiffs_DS = make_histo_of_diffs_DS(diff_t0s) ## ## __________ make simple diagnostic __________ ## diagnostic = TH1D('diagnostic','diagnostic',100,-10.0,10.0) ## for ip in range(int(0.5*nplanes)-1,nplanes): ## for iw in range(nwires): ## diagnostic.Fill(diff_t0s[ip][iw]) print 'done' print 'making mean t0 per plane plots' ## me1, me1_min, me1_max = mean_t0_me(t0_f1)#,'mean_t0s_for_'+nice_fn(f1)) ## me2, me2_min, me2_max = mean_t0_me(t0_f2)#,'mean_t0s_for_'+nice_fn(f2)) ## me1 = mean_t0_me(t0_f1) ## me2 = mean_t0_me(t0_f2,offset_applied) ## me_diff = mean_t0_me(diff_t0s) ## print 'making canvas and histos ...' ## gStyle = makegstyle(gStyle) ## gStyle = nice_colours(gStyle) # c1 = TCanvas() ## choose sensible limits real_min = min( [f1h_min, f2h_min] ) real_max = max( [f1h_max, f2h_max] ) for h in [f1h,f2h]: h.SetMinimum(real_min-0.1*(real_max-real_min)) h.SetMaximum(real_max+0.1*(real_max-real_min)) dh.SetMinimum(-3.) dh.SetMaximum(3.) ## #prettify f1h ## f1h.GetXaxis().CenterTitle() ## f1h.GetYaxis().CenterTitle() ## f1h.GetXaxis().SetTitleOffset(1.5) ## f1h.GetYaxis().SetTitleOffset(1.5) ## f2h.GetXaxis().CenterTitle() ## f2h.GetYaxis().CenterTitle() ## f2h.GetXaxis().SetTitleOffset(1.5) ## f2h.GetYaxis().SetTitleOffset(1.5) ## ## c2 = TCanvas() ## if stats_boxes: ## gStyle.SetOptStat("neMruoi") ## else: ## gStyle.SetOptStat(0) ## #style the diff histos ## histo_diffs.SetLineColor(1) ## histo_diffs.SetLineWidth(2) ## hdiffs_US.SetLineColor(2) ## hdiffs_DS.SetLineColor(4) ## hdiffs_US.SetLineStyle(7) ## hdiffs_DS.SetLineStyle(5) ## hdiffs_US.SetLineWidth(2) ## hdiffs_DS.SetLineWidth(2) ## #draw them ## histo_diffs.Draw() ## hdiffs_US.Draw('same') ## hdiffs_DS.Draw('same') ## #legend ## ldiff = TLegend(.8,.8,.9,.9) ## ldiff.SetBorderSize(1) ## ldiff.AddEntry(histo_diffs,'US+DS','l') ## ldiff.AddEntry(hdiffs_US,'US','l') ## ldiff.AddEntry(hdiffs_DS,'DS','l') ## ldiff.Draw() ## c2.Update() ## ## c3 = TCanvas() ## c3.Divide(1,2) ## # c1_n3_1.SetGridx() #AAH ## # c1_n3_1.SetGridy() #AAH ## # c1_n3_2.SetGridx() #AAH ## # c1_n3_2.SetGridy() #AAH ## c3.cd(1) ## me1.GetXaxis().SetLimits(-1,45) ## me_diff.GetXaxis().SetLimits(-1,45) ## me1.SetMarkerStyle(20) ## me2.SetMarkerStyle(21) ## me1.SetMarkerSize(1) ## me2.SetMarkerSize(1) ## me_diff.SetMarkerStyle(22) ## me_diff.SetMarkerSize(2) ## me1.SetMarkerColor(2) ## me2.SetMarkerColor(4) ## me_diff.SetMarkerColor(1) ## me1.Draw('AP') ## me2.Draw('P') ## me1.SetTitle("Absolute T0s;plane;T0s [ns]") ## #me1.GetHistogram().SetTitle("Absolute T0s;plane;T0s [ns]") ## l = TLegend(0.8,0.8,0.9,0.9) ## l.SetBorderSize(1) ## l.AddEntry(me1,nice_fn(f1),'p') ## l.AddEntry(me2,nice_fn(f2),'p') ## l.Draw() ## c3.cd(2) ## me_diff.Draw('AP') ## me_diff.SetTitle("T0s difference ("+nice_fn(f2)+" - "+nice_fn(f1)+");plane;T0s [ns]") ## c3.Update() #========================================== c4 = TCanvas() dh.SetMinimum(-1.*final_diff_plot_scale) dh.SetMaximum(1.*final_diff_plot_scale) dh.Draw('colz') dh.SetTitleOffset( 0.8, 'Y') c4.Update() #========================================== ## c5 = TCanvas() ## dh.Draw('lego2') ## c5.Update() ## ## ## diagnostic.Draw() ## ## c6 = TCanvas() ## f1h.Draw('lego2') ## c6.Update() ## c7 = TCanvas() ## f2h.Draw('lego2') ## c7.Update() ## ## au = raw_input('>') c4.Print('../T0_MC_Output-Input.eps')