import ROOT,os,sys,argparse from math import ceil, sqrt pandapath=os.environ.get('PANDAPATH') sys.path.append(pandapath+'/macro/tpc/FOPI/mberger') from functions import * from anaFile import anaFile sys.path.append('./images/macros/python') class histColl(): def __init__(self,name): self.name=name self.anafile=None self.fileinfo=None self.hists={} self.proj={} self.cproj={} self.fitdata={} self.graphs={} self.step=8 self.zmin=0 self.zmax=72 self.stepwidth=8 self.nstep=0 def SetAnaFile(self,str): self.anafile=anaFile(str) self.fileinfo=self.anafile.getInfos() def SetHistNames(self,names): for name in names: self.hists[name]=None def AddHistName(self,name): self.hist[name]=None def GetHists(self): for name in self.hists: self.hists[name]=self.anafile.getHist(name) self.zmin=self.hists[name].GetXaxis().GetXmin() self.zmax=self.hists[name].GetXaxis().GetXmax() self.nstep=int(self.hists[name].GetNbinsX()/self.stepwidth) self.stepwidth_z=self.hists[name].GetBinCenter(self.stepwidth)-self.hists[name].GetBinCenter(0) def MakeProj(self): for hname in self.hists: self.proj[hname]=[] start=0 stop=start+self.stepwidth for i in range(self.nstep): start_z=self.hists[hname].GetBinCenter(start) stop_z =self.hists[hname].GetBinCenter(stop) self.proj[hname].append(self.hists[hname].ProjectionY("{0}_py_{1}({2:.0f}-{3:.0f})".format(hname,i,start_z,stop_z),start,stop)) if args.rebin>1: self.proj[hname][-1].Rebin(args.rebin) start=stop stop+=self.stepwidth def DrawProj(self): for hname in self.proj: self.cproj[hname]=ROOT.TCanvas("c_proj_"+hname+self.name,self.name+" "+hname+" Projection") ndiv=int(ceil(sqrt(self.nstep))) self.cproj[hname].Divide(ndiv,ndiv) for iproj in range(len(self.proj[hname])): self.cproj[hname].cd(iproj+1) self.proj[hname][iproj].Draw() if self.fitdata.get(hname,None)!=None: drawGaussians(self.fitdata[hname][iproj]) def FitProj(self): for hname in self.proj: self.fitdata[hname]=[] for iproj in range(len(self.proj[hname])): self.fitdata[hname].append(fitres(self.proj[hname][iproj],0,[self.proj[hname][iproj].GetMaximum(),self.proj[hname][iproj].GetMean(),self.proj[hname][iproj].GetRMS()],[-2000,2000])) def MakeGraph(self): #print self.zmin,self.stepwidth for hname in self.proj: self.graphs[hname]={} self.graphs[hname]["csig"]=ROOT.TGraphErrors() for iproj in range(len(self.proj[hname])): self.graphs[hname]["csig"].SetPoint(iproj,self.zmin+iproj*self.stepwidth_z+self.stepwidth_z/2.,self.fitdata[hname][iproj]["csig"]) self.graphs[hname]["csig"].SetPointError(iproj,self.stepwidth_z/2.,self.fitdata[hname][iproj]["csigerr"]) def SetGraphStyle(self,icol): for hname in self.proj: for gname in self.graphs[hname]: self.graphs[hname][gname].SetTitle(hname) self.graphs[hname][gname].SetMarkerSize(1.5) self.graphs[hname][gname].SetLineWidth(2) self.graphs[hname][gname].SetMarkerColor(colors[icol]) self.graphs[hname][gname].SetLineColor(colors[icol]) self.graphs[hname][gname].GetYaxis().SetTitleOffset(1.2) self.graphs[hname][gname].GetXaxis().SetTitle("z (cm)") self.graphs[hname][gname].GetYaxis().SetTitle("#sigma (#mum)") def getMaxReslDiff(gres,first='recl clcorr',second='recl clcorr cut'): x1,y1=ROOT.Double(0),ROOT.Double(0) x2,y2=ROOT.Double(0),ROOT.Double(0) maxreldiff=0 for ip in range(gres[first].GetN()): gres[first].GetPoint(ip,x1,y1) gres[second].GetPoint(ip,x2,y2) reldiff=(y1-y2)/y1 maxreldiff=max(maxreldiff,reldiff) return maxreldiff def plotProj(coll,hres,hresProj,canv,gres,data): start=0 stop=8 for i in range(9): coll.AddProj() hresProj[corr].append(hres[corr].ProjectionY(hres[corr].GetName()+corr+"_py"+str(i),start,stop)) start=stop stop+=8 opt="ALP" counter=-1 for corr in hres: counter+=1 canv[corr]=ROOT.TCanvas(hres[corr].GetName()+corr,hres[corr].GetTitle()+corr,1000,1000) canv[corr].Divide(3,3) data[corr]=[] gres[corr]=ROOT.TGraphErrors() for i in range(len(hresXProj[corr])): canv[corr].cd(i+1) histo_style.setCanvasMargins(canv[corr].cd(i+1), 0.18,0.05,0.12,0.05) data[corr].append(fitres(hresProj[corr][i],0,[hresProj[corr][i].GetMaximum(),hresProj[corr][i].GetMean(),hresProj[corr][i].GetRMS()],[-2000,2000])) hresProj[corr][i].SetMaximum(hresProj[corr][i].GetMaximum()*1.5) hresProj[corr][i].SetStats(0) hresProj[corr][i].SetTitle("") hresProj[corr][i].GetYaxis().SetTitleOffset(1.7) hresProj[corr][i].GetXaxis().SetNdivisions(506) hresProj[corr][i].SetLineWidth(2) hresProj[corr][i].SetLineColor(ROOT.kAzure+7) histo_style.setAxisLabels(hresProj[corr][i], "Residual (#mum)", "#", "") histo_style.setHistText(hresProj[corr][i],0.06) hresProj[corr][i].Draw() drawGaussians(data[corr][i]) drawresdata(data[corr][i],"sigma ngaus bgaus cgaus title={0}cm-{1}cm times tsize=0.06".format(i*8,(i+1)*8),"#mum",1,0.90,0.055) gres[corr].SetPoint(i,i*8+4, data[corr][i]["csig"]) gres[corr].SetPointError(i,0,data[corr][i]["csigerr"]) canv[corr].Update() if filetype==1: canv[corr].SaveAs("images/chap_Results/sec_Resolution/fitslices_"+corr.replace(" ","_")+"_"+hres[corr].GetName()+".pdf") histo_style.setCanvasMargins(canv['graph'], 0.15, -1 , 0.15, -1) canv["graph"].cd() gres[corr].SetMarkerStyle(20) gres[corr].SetMarkerSize(1.5) gres[corr].SetLineWidth(2) gres[corr].SetMarkerColor(colors[counter]) gres[corr].SetLineColor(colors[counter]) histo_style.setHistText(gres[corr], 0.06, opt) histo_style.setAxisLabels(gres[corr], "z (cm)", "#sigma (#mum)") gres[corr].GetYaxis().SetTitleOffset(1.2) gres[corr].Draw(opt) opt="LP" canv["graph"].Update() parser=argparse.ArgumentParser(description="plot residuals as a function of Z") parser.add_argument('files',help='the input files (output from analyse_cosmics2)',type=str,nargs='+') parser.add_argument('--titles',help='titles to be used for the legend (given in same order as input files)',type=str,nargs='+',default=None) parser.add_argument('--proj',help='draw also the projections',action='store_true') parser.add_argument('--dists',help='draw the distributions themselves',action='store_true') parser.add_argument('--rebin',help='rebin projection histos (%(default)s)',type=int,default=1) parser.add_argument('--hlp',help='help',action='store_true') args=parser.parse_args() if args.hlp: parser.print_help() exit(0) colors=getColors2() if args.titles==None: args.titles=[] for i in range(len(args.files)): args.titles.append(str(i)) leg={} collection=[] canv={} gopt="ALP" for ifile in range(len(args.files)): print 'at file: ',args.files[ifile] collection.append(histColl(args.titles[ifile])) collection[-1].SetAnaFile(args.files[ifile]) collection[-1].SetHistNames(["hClusterResXvsZ","hClusterResYvsZ","hClusterResZvsZ","hClusterResXvsZcut10","hClusterResYvsZcut10","hClusterResZvsZcut10","hClusterResXvsZcut20","hClusterResYvsZcut20","hClusterResZvsZcut20","hClusterResUvsZ","hClusterResVvsZ","hClusterResUvsZcut","hClusterResVvsZcut"]) collection[-1].GetHists() for hname in collection[-1].hists: if ifile==0: leg[hname]=ROOT.TLegend(0.15,0.7,0.5,0.9) leg[hname].SetFillStyle(0) leg[hname].SetBorderSize(0) leg[hname].SetTextFont(132) leg[hname].SetTextSize(0.045) opt='same' if 0: if canv.get(hname,None)==None: canv[hname]=ROOT.TCanvas("c"+hname+args.titles[ifile],hname+" "+args.titles[ifile]) opt='' canv[hname].cd() collection[-1].hists[hname].Draw('colz') collection[-1].MakeProj() collection[-1].FitProj() if args.proj: collection[-1].DrawProj() collection[-1].MakeGraph() collection[-1].SetGraphStyle(ifile) for hname in collection[-1].hists: leg[hname].AddEntry(collection[-1].graphs[hname]["csig"],args.titles[ifile],'lpe') if canv.get("g"+hname,None)==None: cx=700 cy=500 if args.dists: cx=700 cy=1000 canv["g"+hname]=ROOT.TCanvas("c_g_"+hname+args.titles[ifile],"graph "+hname+" "+args.titles[ifile],cx,cy) if args.dists: canv["g"+hname].Divide(1,2) canv["g"+hname].cd(1) collection[-1].hists[hname].Draw("colz") canv["g"+hname].cd(2) else: canv["g"+hname].cd() collection[-1].graphs[hname]["csig"].Draw(gopt) leg[hname].Draw() canv["g"+hname].Update() gopt="LP" thisIsTheEnd()