import sys, os, copy from math import sqrt,cos,sin pandapath = os.environ.get('PANDAPATH') sys.path.append(pandapath + '/macro/tpc/FOPI/python/argparse-1.2.1') sys.path.append(pandapath + '/macro/tpc/FOPI/mberger') import argparse, ROOT from functions import set_palette, openTree, thisIsTheEnd, getSlices import random def writeStepToFile(offxy, offz, diffXY, diffZ, amp, counter, deltaAmp): foldername = "offXY{0}_offZ{1}_diffXY{2}_diffZ{3}_amp{4}".format(offxy, offz, diffXY, diffXY, amp) print 'at ', foldername, counter if(testfile.GetDirectory(foldername)): print "folder exists already, skipping" # amp+=deltaAmp return deltaAmp, counter counter += 1 chi2 = clusterCalibrator.getChi2p(offxy, diffXY, offz, diffXY, amp) print 'found chi2=', chi2 tchi2 = ROOT.TNamed("chi2", str(chi2)) testfile.mkdir(foldername) testfile.cd(foldername) tchi2.Write() hpulls1 = (clusterCalibrator.getCurrentPullSet1()) tmpgraphs1 = getSlices(hpulls1, None, 0, 0) tmpc = ROOT.TCanvas("tmppulls", 'pulls', 1500, 1000) tmpc.Divide(2, 3) tmpc.cd(1) hpulls1.Draw('colz') tmpc.cd(3) tmpgraphs1['Mean'].Draw('AP') tmpc.cd(5) tmpgraphs1['RMS'].Draw('AP') for g in tmpgraphs1: tmpgraphs1[g].Write() hpulls2 = (clusterCalibrator.getCurrentPullSet2()) tmpgraphs2 = getSlices(hpulls1, None, 0, 0) tmpc.cd(2) hpulls2.Draw('colz') tmpc.cd(4) tmpgraphs2['Mean'].Draw('AP') tmpc.cd(6) tmpgraphs2['RMS'].Draw('AP') tmpc.Write() hpulls1.Write() hpulls2.Write() for g in tmpgraphs2: tmpgraphs2[g].Write() del tmpc del hpulls1 del hpulls2 del tmpgraphs1['Mean'] del tmpgraphs1['RMS'] del tmpgraphs2['Mean'] del tmpgraphs2['RMS'] testfile.cd("/") # pargraph.SetPoint(pargraph.GetN()+1,offxy+offz*100+diff*1000+amp*10000000,chi2) return deltaAmp , counter parser = argparse.ArgumentParser(description='calibrate parametrization of cluster error') parser.add_argument('recofile', help='the reco file, or path to the folder with the runfolders', type=str,nargs='+') parser.add_argument('paramfile', help='the param file', type=str) parser.add_argument('--list',help='the recofile is actually a list of files',action='store_true') parser.add_argument('--runs', help='if data, the runs to use (%(default)s)', type=str, default='3888') parser.add_argument('--pattern', help='if data, the file name patters ('')', type=str, default='') parser.add_argument('--events', type=int, default=1000,help="number of events to process ( %(default)i )") parser.add_argument('--ncluster', help='how many cluster to analyze', type=int, default=-1) parser.add_argument('--batPdf', help='name of the bat pdf', type=str, default='bat.pdf') parser.add_argument('--fake', help='use fake distributions', action='store_true') parser.add_argument('--nIter', help='MaxIterations PreRun and MaxIterations Run (%(default)s)', type=int, nargs=2, default=[6000, 3000]) parser.add_argument('--chi2Gaus', help='use Gaus chi2 calc', action='store_true') parser.add_argument('--chi2Slope', help='use Slope chi2 calc', action='store_true') parser.add_argument('--chi2Mean', help='use Mean chi2 calc', action='store_true') parser.add_argument('--chi2RMS', help='use RMS chi2 calc', action='store_true') parser.add_argument('--ltestGaus', help='use the likelihood test', action='store_true') parser.add_argument('--fixGausMean', help='fix the mean of the gaus in chi2 gaus calc', action='store_true') parser.add_argument('--fitGaus', help='fit the gaus or calc by hand', action='store_true') parser.add_argument('--diffPar', help='fix the diffusion parameter in bat model', type=float, default=-1) parser.add_argument('--ampPar', help='fix the amplitude parameter in bat model', type=float, default=-1) parser.add_argument('--offxyPar', help='fix the offsetXY parameter in bat model', type=float, default=-1) parser.add_argument('--offzPar', help='fix the offsetZ parameter in bat model', type=float, default=-1) parser.add_argument('--batch', help='batchmode', action='store_true') parser.add_argument('--parscan', help='perform a parameter scan', action='store_true') parser.add_argument('--noMarg',help='dont run marginalization',action='store_true') parser.add_argument('--scanXY',help='scan over parameter offsetXY (start,stop,step)',type=float,nargs=3,default=[0,0,1]) parser.add_argument('--scanZ',help='scan over parameter offsetZ (start,stop,step)',type=float,nargs=3,default=[0,0,1]) parser.add_argument('--scanAmp',help='scan over parameter Amp (start,stop,step)',type=float,nargs=3,default=[0,0,1]) parser.add_argument('--scanDiffXY',help='scan over parameter diffXY (start,stop,step)',type=float,nargs=3,default=[0,0,1]) parser.add_argument('--scanDiffZ',help='scan over parameter diffZ (start,stop,step)',type=float,nargs=3,default=[0,0,1]) parser.add_argument('--scanFile',help='name of the scan filename',type=str,default="bat_parameter_scan.root") parser.add_argument('--noprf', help='disable the prf weighting', action='store_true') parser.add_argument('--dataCuts', help='enable cuts for real cosmics data', action='store_true') parser.add_argument('--badZCut',help='enable cuts for steep tracks at small drift',action='store_true') parser.add_argument('--wait',help='wait after filling of histograms',action='store_true') #parser.add_argument('--errParFile',help='path to error par file to use',type=str,default='None') parser.add_argument('--useChi2Theta',help='enable chi2 calculation as function of theta',action='store_true') parser.add_argument('--useChi2ZTheta',help='enable chi2 calculation in bins of Z and theta',action='store_true') parser.add_argument('--fopi',help='use fopi data',action='store_true') parser.add_argument('--tpcID',help="detector ID of the tpc (default=%(default)s)",default=22,type=int) args = parser.parse_args() ROOT.gROOT.ProcessLine(".x rootlogon.C") ROOT.gROOT.ProcessLine('gROOT->SetStyle("Plain")') ROOT.gROOT.LoadMacro("stlPYROOT.h+") set_palette() if args.batch: ROOT.gROOT.SetBatch(True) if args.ncluster != -1: args.events = 1000000000 # silence roofit ROOT.RooMsgService.instance().setGlobalKillBelow(ROOT.RooFit.ERROR) ROOT.RooMsgService.instance().setSilentMode(True) pcanv = ROOT.TCanvas("cpulls", 'pulls', 1500, 1000) pcanv.Divide(4, 3) pcanv2 = ROOT.TCanvas('cangles', 'Angles', 1500, 1000) pcanv2.Divide(4, 2) pcanv3 = ROOT.TCanvas("ctimes", "times", 1000, 500) pcanv3.Divide(2, 1) pcanv4=ROOT.TCanvas("cpullsTheta",'pulls vs theta',1500,1000) pcanv4.Divide(4,3) rcanv = ROOT.TCanvas("cres","residuals",1500,1000) rcanv.Divide(4,3) scanv=ROOT.TCanvas('csig','errors',1500,1000) scanv.Divide(4,3) posCanv=ROOT.TCanvas("cpos",'pos',1500,1000) # exit() # first create the new tree to work with force_stop = False if force_stop: exit() ev = -1 global clusterCalibrator clusterCalibrator = ROOT.TpcClusterErrorCalib() if not args.fake: trees=[] files=[] if args.list: for line in open(args.recofile[0]): line=line.replace('\n','') words=line.split(';') files.append(words[0]) else: files=args.recofile for itree in files: trees.append(openTree(itree,args.runs,args.pattern)) trees[-1].SetBranchStatus("*", 0) if args.fopi: trees[-1].SetBranchStatus("Particles.*",1) else: trees[-1].SetBranchStatus("CutCosmics.*", 1) #tree = openTree(args.recofile, args.runs, args.pattern) #tree.SetBranchStatus("*", 0) #tree.SetBranchStatus("CutCosmics.*", 1) # initialize the singletons if os.path.isfile(args.paramfile): parfile = ROOT.TFile.Open(args.paramfile) digiPar = parfile.Get('TpcDigiPar') errScaler = ROOT.TpcClusterErrorScaler.getInstance() errScaler.SetNoFairRun(True) errScaler.init(digiPar.getClusterErrorFile(),61) #if args.errParFile!='None': # errScaler.init(args.errParFile,61) #else: # if args.pattern=='': # errScaler.init("./tpc/FOPI/par/tpc.errorParam.86.par", 61) # else: # errScaler.init("./tpc/FOPI/par/tpc.errorParam.86.data.par", 61) digiMap = ROOT.TpcDigiMapper.getInstance() digiMap.init(digiPar) padplane=digiPar.getPadPlane() wx=padplane.GetPad(0).GetWeight(0.99,1.,0); wy=padplane.GetPad(0).GetWeight(0.99,0,1.); print "weights:",wx,wy Dl = digiPar.getGas().Dl() Dt = digiPar.getGas().Dt() print 'found diffusion: Dl=', Dl, 'Dt=', Dt dx, dy = ROOT.Double(0), ROOT.Double(0) digiMap.padsize(0, dx, dy) dummyID = ROOT.McId(1, 1, 0) dummyColl = ROOT.McIdCollection() dummyColl.AddID(dummyID) zDiff1 = ROOT.TVector3() zDiff2 = ROOT.TVector3() digi1 = ROOT.TpcDigi(1, 1, 1, dummyColl) digi2 = ROOT.TpcDigi(1, 2, 1, dummyColl) digiMap.map(digi1, zDiff1) digiMap.map(digi2, zDiff2) dz = zDiff2.Z() - zDiff1.Z() print 'dx=', dx, 'dy=', dy, '(',dy / 2 * (1 + ROOT.TMath.Sin(30 * ROOT.TMath.DegToRad())),')', 'dz=', dz, 'gain=', digiPar.getGain() globcounter = -1 treecounter=-1 lastglob=0 lastev=0 print 'collecting clusters' for tree in trees: treecounter+=1 print "at tree",treecounter,"from file",files[treecounter],"(",tree.GetEntries(),")" ev = -1 lastev=0 for e in tree: ev += 1 if ev > args.events: break if ev % 100 == 0 and args.ncluster==-1: print 'at event:', ev, "( nCluster="+str(globcounter)+")" elif args.ncluster!=-1 and globcounter % 1000 == 0 and lastglob!=globcounter: print 'at event:', ev, "( nCluster="+str(globcounter)+")" lastglob=globcounter if args.fopi: branch=e.Particles else: branch=e.CutCosmics for tr in branch: #if args.dataCuts: # if abs(tr.GetTheta() - 90) > 2: # break for icl in range(tr.nhits()): if args.fopi: if tr.getDetID()!=args.tpcID: continue if args.dataCuts: if tr.GetCogPos(icl).Perp() > 13 or tr.GetCogPos(icl).Perp() < 7: continue if tr.GetCogPos(icl).Z()>55: continue if args.badZCut: if tr.GetCogPos(icl).Z()<12 and abs(tr.GetTheta() - 90)<10: continue if args.ncluster != -1 and globcounter > args.ncluster: ev = args.events + 1 break globcounter += 1 calibCluster = ROOT.TpcCalibCluster() calibCluster.setRefPos(tr.GetProjectionPoint(icl)) calibCluster.setPlane(tr.GetMidPlaneO(icl), tr.GetMidPlaneU(icl), tr.GetMidPlaneV(icl)) calibCluster.setSize(tr.GetFullClusterSize(icl)) calibCluster.setShapeCov(ROOT.TMatrixDSym(2,tr.GetShapeCovOnPlane(icl).GetMatrixArray())) calibCluster.setPos(tr.GetHitPosistionXYZ(icl)) calibCluster.setID(globcounter) calibCluster.setDAFWeight(tr.GetDAFWeights(icl)) #offsets calibCluster.setParam(0,errScaler.getParameter(0, 0)) calibCluster.setParam(1,errScaler.getParameter(1, 0)) calibCluster.setParam(2,errScaler.getParameter(2, 0)) #diffusion calibCluster.setParam(3,errScaler.getParameter(0, 1)) calibCluster.setParam(4,errScaler.getParameter(1, 1)) calibCluster.setParam(5,errScaler.getParameter(2, 1)) #finite pad calibCluster.setParam(6,errScaler.getParameter(0, 2)) calibCluster.setParam(7,errScaler.getParameter(1, 2)) calibCluster.setParam(8,errScaler.getParameter(2, 2)) #angle calibCluster.setParam(9,errScaler.getParameter(0, 3)) calibCluster.setParam(10,errScaler.getParameter(1, 3)) calibCluster.setParam(11,errScaler.getParameter(2, 3)) #angle decay calibCluster.setParam(12,errScaler.getParameter(0, 4)) calibCluster.setParam(13,errScaler.getParameter(1, 4)) calibCluster.setParam(14,errScaler.getParameter(2, 4)) calibCluster.setParam(15,errScaler.getParameter(0, 5)) calibCluster.setParam(16,errScaler.getParameter(1, 5)) calibCluster.setParam(17,errScaler.getParameter(2, 5)) calibCluster.setParam(18,errScaler.getParameter(0, 6)) calibCluster.setParam(19,errScaler.getParameter(1, 6)) calibCluster.setParam(20,errScaler.getParameter(2, 6)) calibCluster.setParam(21,errScaler.getParameter(0, 7)) clusterCalibrator.addCalibCluster(calibCluster) lastev=ev else: print 'faking clusters' for i in range(args.events): calibCluster = ROOT.TpcCalibCluster() normal=ROOT.TVector3(1,0,0) normal.SetTheta(random.uniform(0,ROOT.TMath.Pi())) normal.SetMag(1) refpos=ROOT.TVector3(0,0,ROOT.gRandom.Uniform(0, 75)) calibCluster.setRefPos(refpos) calibCluster.setPlaneON(refpos,normal) diffVect=ROOT.TVector3(0,0,0) diffVect.SetX(random.gauss(0,0.01*sqrt(refpos.Z()) ))#+0.02*(cos(normal.Theta()))**2 ) ) diffVect.SetY(random.gauss(0,0.01*sqrt(refpos.Z()) ))#+0.02*(cos(normal.Theta()))**2 ) ) diffVect.SetZ(random.gauss(0,0.02*sqrt(refpos.Z()) ))#+0.02*(sin(normal.Theta())) ) ) clpos=ROOT.TVector3(refpos) clpos+=diffVect calibCluster.setPos(clpos) calibCluster.setDiffConst(1, 1) calibCluster.setGain(1) calibCluster.setPadGeom(1 , 1 , 1) calibCluster.setWeights(1,1) calibCluster.setFake(True) #calibCluster.setID(globcounter) #if i%10==0: # calibCluster.addDigi(ROOT.TVector3(random.gauss(0,8),random.gauss(0,8),random.gauss(0,8)),1) ''' x=random.gauss(0, 1) y=random.gauss(0, 1) z=random.gauss(0, 1) calibCluster.addDigi(ROOT.TVector3(x,y,z), 1) ''' #for i in range(20): for i in range(12): calibCluster.setParam(i,0) calibCluster.setParam(2,0.015**2) calibCluster.setParam(3,0.015**2) #calibCluster.setParams(0, 1, 0, 1) clusterCalibrator.addCalibCluster(calibCluster) print 'start calibration' clusterCalibrator.setMCMCIterations(args.nIter[0], args.nIter[1]) clusterCalibrator.setUseChi2Gaus(args.chi2Gaus) clusterCalibrator.setUseChi2Mean(args.chi2Mean) clusterCalibrator.setUseChi2Slope(args.chi2Slope) clusterCalibrator.setUseChi2RMS(args.chi2RMS) clusterCalibrator.setUseLTest(args.ltestGaus) clusterCalibrator.setFixGausMean(args.fixGausMean) clusterCalibrator.setFixDiffusion(args.diffPar) clusterCalibrator.setFixAmp(args.ampPar) clusterCalibrator.setFixOffsetXY(args.offxyPar) clusterCalibrator.setFixOffsetZ(args.offzPar) #if(args.fake): # clusterCalibrator.setFixOffsetZ(0) clusterCalibrator.setFitGaus(args.fitGaus) clusterCalibrator.setUseFake(args.fake) #clusterCalibrator.setDrawGauss(True) clusterCalibrator.setPdfName(args.batPdf) clusterCalibrator.fillHPulls(True,0) clusterCalibrator.fillAngleHists(0) clusterCalibrator.setUseThetaChi2(args.useChi2Theta) clusterCalibrator.setUseZThetaChi2(args.useChi2ZTheta) pulls = {} projs = {} graphs = {} angles = {} times = {} residuals={} errors={} set_palette() pcanv.cd(1) pulls['hP1vsZ_before'] = clusterCalibrator.getHist("P1vsZ before") pulls['hP1vsZ_before'].Draw('colz') projs['hP1vsZ_before_proj'] = [] graphs['gP1vsZ_before'] = copy.deepcopy(getSlices(pulls['hP1vsZ_before'], projs['hP1vsZ_before_proj'], 0, 0,'t')) pcanv.cd(5) graphs['gP1vsZ_before']['Mean'].Draw('AP') pcanv.cd(9) graphs['gP1vsZ_before']['RMS'].Draw('AP') pcanv.cd(3) pulls['hP2vsZ_before'] = clusterCalibrator.getHist("P2vsZ before") pulls['hP2vsZ_before'].Draw('colz') projs['hP2vsZ_before_proj'] = [] graphs['gP2vsZ_before'] = copy.deepcopy(getSlices(pulls['hP2vsZ_before'], projs['hP2vsZ_before_proj'], 0, 0,'t')) pcanv.cd(7) graphs['gP2vsZ_before']['Mean'].Draw('AP') pcanv.cd(11) graphs['gP2vsZ_before']['RMS'].Draw('AP') pcanv.Update() pcanv2.cd(1) angles['htheta1_before'] = copy.deepcopy(clusterCalibrator.getHist("ThetaEv1 before")) angles['htheta1_before'].Draw('colz') pcanv2.cd(3) angles['htheta2_before'] = copy.deepcopy(clusterCalibrator.getHist("ThetaEv2 before")) angles['htheta2_before'].Draw('colz') pcanv2.Update() rcanv.cd(1) residuals['R1vsZ before']=clusterCalibrator.getHist("R1vsZ before") residuals['R1vsZ before'].Draw('colz') rcanv.cd(3) residuals['R2vsZ before']=clusterCalibrator.getHist("R2vsZ before") residuals['R2vsZ before'].Draw('colz') scanv.cd(1) errors['E1vsZ before']=clusterCalibrator.getHist("E1vsZ before") errors['E1vsZ before'].Draw('colz') scanv.cd(3) errors['E2vsZ before']=clusterCalibrator.getHist("E2vsZ before") errors['E2vsZ before'].Draw('colz') # clusterCalibrator.getChi2p(0,1,0.1,0) scanv.cd(5) errors['E1vsT before']=clusterCalibrator.getHist("E1vsT before") errors['E1vsT before'].Draw('colz') scanv.cd(7) errors['E2vsT before']=clusterCalibrator.getHist("E2vsT before") errors['E2vsT before'].Draw('colz') pcanv4.cd(1) pulls['hP1vsTheta_before'] = clusterCalibrator.getHist("P1vsTheta before") pulls['hP1vsTheta_before'].Draw('colz') projs['hP1vsTheta_before_proj'] = [] graphs['gP1vsTheta_before'] = copy.deepcopy(getSlices(pulls['hP1vsTheta_before'], projs['hP1vsTheta_before_proj'], 0, 0,'t')) pcanv4.cd(5) graphs['gP1vsTheta_before']['Mean'].Draw('AP') pcanv4.cd(9) graphs['gP1vsTheta_before']['RMS'].Draw('AP') pcanv4.cd(3) pulls['hP2vsTheta_before'] = clusterCalibrator.getHist("P2vsTheta before") pulls['hP2vsTheta_before'].Draw('colz') projs['hP2vsTheta_before_proj'] = [] graphs['gP2vsTheta_before'] = copy.deepcopy(getSlices(pulls['hP2vsTheta_before'], projs['hP2vsTheta_before_proj'], 0, 0,'t')) pcanv4.cd(7) graphs['gP2vsTheta_before']['Mean'].Draw('AP') pcanv4.cd(11) graphs['gP2vsTheta_before']['RMS'].Draw('AP') pcanv4.Update() posCanv.cd() posHist=clusterCalibrator.getHist("hClPos") posHist.Draw() if args.parscan : global testfile testfile = ROOT.TFile(args.scanFile, "update") ROOT.gROOT.SetBatch(True) # print testfile.Get('pargraph') # new=False # if testfile.Get('pargraph')!=None: # print 'found pargraph' # pargraph=testfile.Get('pargraph') # new=False # else: # pargraph=ROOT.TGraph() # pargraph.SetName('pargraph') counter = -1 offxy = args.scanXY[0] print args.scanXY print args.scanZ print args.scanAmp print args.scanDiffXY print args.scanDiffZ deltaAmp = args.scanAmp[2] while offxy <= args.scanXY[1] and not force_stop: print 'at offxy:',offxy offz = args.scanZ[0] while offz <= args.scanZ[1] and not force_stop: print 'at offz:',offz diffXY = args.scanDiffXY[0] while diffXY <= args.scanDiffXY[1] and not force_stop: print 'at diffXY:',diffXY diffZ=args.scanDiffZ[0] while diffZ <=args.scanDiffZ[1] and not force_stop: print 'at diffZ:',diffZ amp = args.scanAmp[0] # deltaAmp #if amp==0: # amp=deltaAmp while amp <= args.scanAmp[1]: print 'at amp::',amp if counter > 200: print 'counter break' force_stop = True break da, counter = writeStepToFile(offxy, offz, diffXY, diffZ, amp, counter, deltaAmp) amp += da diffZ += args.scanDiffZ[2] diffXY += args.scanDiffXY[2] offz += args.scanZ[2] offxy += args.scanXY[2] ROOT.gROOT.SetBatch(False) #tmpc = ROOT.TCanvas() # setGraph(pargraph, "Chi2 vs parsetting", 20, 0.5, 1, "offXY+offZ*100+diff*10000+amp*10000000", "chi2") # pargraph.Draw('AP') testfile.cd('/') # if(new): # pargraph.Write() testfile.Close() elif not args.noMarg: u='w' if args.wait: u=raw_input('start') if u=='q': exit() clusterCalibrator.execute() parfile=open(args.batPdf.replace(".pdf","_param.txt"),'w') parfile_err=open(args.batPdf.replace(".pdf","_paramErr.txt"),'w') print "getting vest parameter" pars=clusterCalibrator.getBestParameter() print "getting best parameter errors" pars_err=clusterCalibrator.getBestParameterErrors() ipar=0 for i in range(3): for j in range(7): parfile.write("{0:.5e} ".format(pars[j*3+i])) parfile.write("{0:.5e}\n".format(pars[21])) parfile.close() parfile_err.close() else: pars=[] pars.append(errScaler.getParameter(0, 0)) #off x 0 pars.append(errScaler.getParameter(1, 0)) #off y 1 pars.append(errScaler.getParameter(2, 0)) #off z 2 pars.append(errScaler.getParameter(0, 1)) #diff x 3 pars.append(errScaler.getParameter(1, 1)) #diff y 4 pars.append(errScaler.getParameter(2, 1)) #diff z 5 pars.append(errScaler.getParameter(0, 2)) #pad decay X 6 pars.append(errScaler.getParameter(1, 2)) #pad decay Y 7 pars.append(errScaler.getParameter(2, 2)) #pad decay Z 8 pars.append(errScaler.getParameter(0, 3)) #pad decay amp x 9 pars.append(errScaler.getParameter(1, 3)) #pad decay amp y 10 pars.append(errScaler.getParameter(2, 3)) #pad decay amp z 11 pars.append(errScaler.getParameter(0, 4)) #angle x 12 pars.append(errScaler.getParameter(1, 4)) #angle 7 13 pars.append(errScaler.getParameter(2, 4)) #angle z 14 pars.append(errScaler.getParameter(0, 5)) #angle^2 x 15 pars.append(errScaler.getParameter(1, 5)) #angle^2 y 16 pars.append(errScaler.getParameter(2, 5)) #angle^2 z 17 pars.append(errScaler.getParameter(0, 6)) #angle mix x 18 pars.append(errScaler.getParameter(1, 6)) #angle mix y 19 pars.append(errScaler.getParameter(2, 6)) #angle mix z 20 pars.append(errScaler.getParameter(0, 7)) #fit contrib 21 for i in pars: clusterCalibrator.addFakePar(i) chi2_tmp=clusterCalibrator.getChi2p(True) print' chi2({0})={1}'.format(1,chi2_tmp) ''' graph=ROOT.TGraph() dpar_step=0.0005 #dpar=0.015-50*dpar_step dpar=0.001 for idiff in range(100): dpar+=dpar_step clusterCalibrator.setFakePar(2,dpar**2) clusterCalibrator.setFakePar(3,dpar**2) chi2_tmp=clusterCalibrator.getChi2p(True) print' chi2({0})={1}'.format(dpar,chi2_tmp) graph.SetPoint(graph.GetN(),dpar,chi2_tmp) ctmp=ROOT.TCanvas() graph.Draw("APL") ''' ''' chi2_tmp=clusterCalibrator.getChi2p(True) print' chi2=',chi2_tmp clusterCalibrator.setFakePar(2,0.0015**2) clusterCalibrator.setFakePar(3,0.0015**2) chi2_tmp=clusterCalibrator.getChi2p(True) print' chi2=',chi2_tmp ''' u='w' if args.wait: u=raw_input('start') if u=='q': exit() print "filling hpulls" clusterCalibrator.fillHPulls(False,1) print "filling angle" clusterCalibrator.fillAngleHists(1) set_palette() pcanv.cd(2) print "getting pulls1 after" pulls['hP1vsZ_after'] = clusterCalibrator.getHist("P1vsZ after") pulls['hP1vsZ_after'].Draw('colz') projs['hP1vsZ_after_proj'] = [] graphs['gP1vsZ_after'] = copy.deepcopy(getSlices(pulls['hP1vsZ_after'], projs['hP1vsZ_after_proj'], 0, 0,'t')) pcanv.cd(6) graphs['gP1vsZ_after']['Mean'].Draw('AP') pcanv.cd(10) graphs['gP1vsZ_after']['RMS'].Draw('AP') print "getting pulls2 after" pcanv.cd(4) pulls['hP2vsZ_after'] = clusterCalibrator.getHist("P2vsZ after") pulls['hP2vsZ_after'].Draw('colz') projs['hP2vsZ_after_proj'] = [] graphs['gP2vsZ_after'] = copy.deepcopy(getSlices(pulls['hP2vsZ_after'], projs['hP2vsZ_after_proj'], 0, 0,'t')) pcanv.cd(8) graphs['gP2vsZ_after']['Mean'].Draw('AP') pcanv.cd(12) graphs['gP2vsZ_after']['RMS'].Draw('AP') pcanv.Update() print "getting theta1 after" pcanv2.cd(2) angles['htheta1_after'] = copy.deepcopy(clusterCalibrator.getHist("ThetaEv1 after")) angles['htheta1_after'].Draw('colz') pcanv2.cd(4) print "getting theta2 after" angles['htheta2_after'] = copy.deepcopy(clusterCalibrator.getHist("ThetaEv2 after")) angles['htheta2_after'].Draw('colz') pcanv2.Update() rcanv.cd(2) print "getting res1 after" residuals['R1vsZ after']=clusterCalibrator.getHist("R1vsZ after") residuals['R1vsZ after'].Draw('colz') rcanv.cd(4) print "getting res2 after" residuals['R2vsZ after']=clusterCalibrator.getHist("R2vsZ after") residuals['R2vsZ after'].Draw('colz') rcanv.Update() scanv.cd(2) errors['E1vsZ after']=clusterCalibrator.getHist("E1vsZ after") errors['E1vsZ after'].Draw('colz') scanv.cd(4) errors['E2vsZ after']=clusterCalibrator.getHist("E2vsZ after") errors['E2vsZ after'].Draw('colz') scanv.cd(6) errors['E1vsT after']=clusterCalibrator.getHist("E1vsT after") errors['E1vsT after'].Draw('colz') scanv.cd(8) errors['E2vsT after']=clusterCalibrator.getHist("E2vsT after") errors['E2vsT after'].Draw('colz') scanv.Update() pcanv4.cd(2) print "getting pullsvsTheta1 after" pulls['hP1vsTheta after'] = clusterCalibrator.getHist("P1vsTheta after") pulls['hP1vsTheta after'].Draw('colz') projs['hP1vsTheta_after_proj'] = [] graphs['gP1vsTheta_after'] = copy.deepcopy(getSlices(pulls['hP1vsTheta after'], projs['hP1vsTheta_after_proj'], 0, 0,'t')) pcanv4.cd(6) graphs['gP1vsTheta_after']['Mean'].Draw('AP') pcanv4.cd(10) graphs['gP1vsTheta_after']['RMS'].Draw('AP') pcanv4.cd(4) print "getting pullsvsTheta2 after" pulls['hP2vsTheta after'] = clusterCalibrator.getHist("P2vsTheta after") pulls['hP2vsTheta after'].Draw('colz') projs['hP2vsTheta_after_proj'] = [] graphs['gP2vsTheta_after'] = copy.deepcopy(getSlices(pulls['hP2vsTheta after'], projs['hP2vsTheta_after_proj'], 0, 0,'t')) pcanv4.cd(8) graphs['gP2vsTheta_after']['Mean'].Draw('AP') pcanv4.cd(12) graphs['gP2vsTheta_after']['RMS'].Draw('AP') pcanv4.Update() pcanv3.cd(1) times['htimeChi2'] = copy.deepcopy(clusterCalibrator.getHist("TimeChi2")) times['htimeChi2'].Draw() pcanv3.cd(2) times['htimeData'] = copy.deepcopy(clusterCalibrator.getHist("TimeData")) times['htimeData'].Draw() pcanv3.Update() # calibModel=ROOT.BatClusterErrorModel() routfile = ROOT.TFile(args.batPdf.replace(".pdf", ".root"), "recreate") pcanv.Write() pcanv4.Write() scanv.Write() rcanv.Write() routfile.cd('/') routfile.mkdir('pulls') routfile.cd('pulls') print 'write pulls' for h in pulls: pulls[h].Write() routfile.cd('/') routfile.mkdir("projs") routfile.cd('projs') print 'write projs' for p in projs: routfile.mkdir('projs/' + p) routfile.cd('projs/' + p) for ip in projs[p]: ip.Write() routfile.cd('/') routfile.mkdir('graphs') print 'write graphs' for g in graphs: routfile.mkdir('graphs/' + g) routfile.cd('graphs/' + g) for ig in graphs[g]: graphs[g][ig].Write() routfile.cd('/') routfile.mkdir('angles') routfile.cd('angles') print 'write angles' for h in angles: angles[h].Write() routfile.cd('/') routfile.mkdir('times') routfile.cd('times') print 'writing times' for h in times: times[h].Write() routfile.cd('/') routfile.mkdir('res') routfile.cd('res') for h in residuals: residuals[h].Write() routfile.cd('/') routfile.mkdir('err') routfile.cd('err') for h in errors: errors[h].Write() if not args.batch: thisIsTheEnd() routfile.Close()