Reference Manual

Go to the documentation of this file.
00001 
00002 ####### 
00003 
00004 #  automatized plots generator for b-tagging performances
00005 #  Adrien Caudron, 2013, UCL
00006 
00007 #######
00008 
00009 #do all import
00010 import os, sys
00011 
00012 try:
00013     import ROOT
00014 except:
00015     print "\nCannot load PYROOT, make sure you have setup ROOT in the path"
00016     print "and pyroot library is also defined in the variable PYTHONPATH, try:\n"
00017     if (os.getenv("PYTHONPATH")):
00018         print " setenv PYTHONPATH ${PYTHONPATH}:$ROOTSYS/lib\n"
00019     else:
00020         print " setenv PYTHONPATH $ROOTSYS/lib\n"
00021         sys.exit()
00022 
00023 from ROOT import TFile
00024 from ROOT import TCanvas
00025 from ROOT import TPad
00026 from ROOT import TLegend
00027 from ROOT import TLatex
00028 from ROOT import TH1F
00029 from ROOT import TF1
00030 from ROOT import TVectorD
00031 from ROOT import TGraphErrors
00032 from ROOT import Double
00033 
00034 import Style
00035 from listHistos import *
00036 
00037 #define the input root files
00038 #fileVal = TFile("BTagRelVal_TTbar_Startup_14_612SLHC2.root","READ")
00039 #fileRef = TFile("BTagRelVal_TTbar_Startup_612SLHC1_14.root","READ")
00040 fileNameVal = "BTagRelVal_TTbar_Startup_14_612SLHC2.root"
00041 fileNameRef = "BTagRelVal_TTbar_Startup_612SLHC1_14.root"
00042 #define the val/ref labels
00043 ValRel = "612SLHC2_14"
00044 RefRel = "612SLHC1_14"
00045 #define the sample labels
00046 ValSample = "TTbar_FullSim"
00047 RefSample = "TTbar_FullSim"
00048 #define different settings
00049 batch = False #run on batch mode ?
00050 weight = 1. #rescale the histos according to this weight
00051 drawLegend = False #draw legend ?
00052 printBanner = False #draw text Banner on top of the histos
00053 Banner = "CMS Preliminary"
00054 doRatio = True #plot the ratios
00055 drawOption = "" # "" or "HIST"
00056 #path in the file
00057 pathInFile = "/DQMData/Run 1/Btag/Run summary/"
00058 #ETA/PT bins, GLOBAL ?
00059 EtaPtBin =[
00060     "GLOBAL",
00061     #"ETA_0-1v4",
00062     #"ETA_1v4-2v4",
00063     #"PT_50-80",
00064     #"PT_80-120",
00065     ]
00066 #list of taggers to look at
00067 listTag = [
00068     "CSV",
00069     "CSVMVA",
00070     "JP",
00071     "JBP",
00072     "TCHE",
00073     "TCHP",
00074     "SSVHE",
00075     "SSVHP",
00076     "SMT",
00077     #"SMTIP3d",
00078     #"SMTPt",
00079     "SET",
00080     ]
00081 #list of flavors to look at
00082 listFlavors = [
00083         #"ALL",
00084         "B",
00085         "C",
00086         #"G",
00087         #"DUS",
00088         "DUSG",
00089         #"NI",
00090         ]
00091 #map for marker color for flav-col and tag-col
00092 mapColor = {
00093     "ALL"  : 4 ,
00094     "B"    : 3 ,
00095     "C"    : 1 ,
00096     "G"    : 2 ,
00097     "DUS"  : 2 ,
00098     "DUSG" : 2 ,
00099     "NI"   : 5 ,
00100     "CSV"       : 5 ,
00101     "CSVMVA"   : 6 ,
00102     "JP"        : 3 ,
00103     "JBP"       : 9 ,
00104     "TCHE"      : 1,
00105     "TCHP"      : 2,
00106     "SSVHE"     : 4,
00107     "SSVHP"     : 7,
00108     "SMT"       : 8 ,
00109     "SMTIP3d" : 11 ,
00110     "SMTPt"   : 12
00111     }
00112 #marker style map for Val/Ref
00113 mapMarker = {
00114     "Val" : 22,
00115     "Ref" :  8
00116     }
00117 mapLineWidth = {
00118     "Val" : 3,
00119     "Ref" : 2
00120     }
00121 mapLineStyle = {
00122     "Val" : 2,
00123     "Ref" : 1
00124     }
00125 #choose the formats to save the plots 
00126 listFromats = [
00127     "gif",
00128     ]
00129 #unity function
00130 unity = TF1("unity","1",-1000,1000)
00131 unity.SetLineColor(8)
00132 unity.SetLineWidth(1)
00133 unity.SetLineStyle(1)
00134 #list of histos to plots
00135 listHistos = [
00136     jetPt,
00137     jetEta,
00138     discr,
00139     effVsDiscrCut_discr,
00140     FlavEffVsBEff_discr,
00141     performance,
00142     #performanceC,
00143 
00144     #IP,
00145     #IPe,
00146     #IPs,
00147     #NTracks,
00148     #decayLength,
00149     #distToJetAxis,
00150     #NHits,
00151     #NPixelHits,
00152     #NormChi2,
00153     #trackPt,
00154 
00155 
00156     #flightDist3Dval,
00157     #flightDist3Dsig,
00158     #jetNSecondaryVertices,
00159     
00160     #vertexMass,
00161     #vertexNTracks,
00162     #vertexJetDeltaR,
00163     #vertexEnergyRatio,
00164 
00165     #vertexCategory,
00166     #trackSip3dVal,
00167     #trackSip3dSig,
00168     #trackSip3dSigAboveCharm,
00169     #trackDeltaR,
00170     #trackEtaRel,
00171     #trackDecayLenVal,
00172     #trackSumJetDeltaR,
00173     #trackJetDist,
00174     #trackSumJetEtRatio,
00175     #trackPtRel,
00176     #trackPtRatio,
00177     #trackMomentum,
00178     #trackPPar,
00179     #trackPParRatio,
00180     ]
00181 
00182 #methode to do a plot from histos       
00183 def histoProducer(plot,histos,keys,isVal=True):
00184     if histos is None : return
00185     if isVal : sample = "Val"
00186     else : sample = "Ref"
00187     outhistos = []
00188     minY=9999.
00189     maxY=0.
00190     for k in keys :
00191         #Binning
00192         if plot.binning and len(plot.binning)==3 :
00193             histos[k].SetBins(plot.binning[0],plot.binning[1],plot.binning[2])
00194         elif plot.binning and len(plot.binning)==2 :
00195             nbins=plot.binning[1]+1-plot.binning[0]
00196             xmin=histos[k].GetBinLowEdge(plot.binning[0])
00197             xmax=histos[k].GetBinLowEdge(plot.binning[1]+1)
00198             valtmp=TH1F(histos[k].GetName(),histos[k].GetTitle(),nbins,xmin,xmax)
00199             i=1
00200             for bin in range(plot.binning[0],plot.binning[1]+1) :
00201                 valtmp.SetBinContent(i,histos[k].GetBinContent(bin))
00202                 i+=1
00203             histos[k]=valtmp
00204         if plot.Rebin and plot.Rebin > 0 :
00205             histos[k].Rebin(plot.Rebin)
00206         #Style
00207         histos[k].SetLineColor(mapColor[k])
00208         histos[k].SetMarkerColor(mapColor[k])
00209         histos[k].SetMarkerStyle(mapMarker[sample])
00210         if drawOption == "HIST" :
00211             histos[k].SetLineWidth(mapLineWidth[sample])
00212             histos[k].SetLineStyle(mapLineStyle[sample])
00213         #compute errors
00214         histos[k].Sumw2()
00215         #do the norm
00216         if plot.doNormalization :
00217             histos[k].Scale(1./histos[k].Integral())
00218         elif weight!=1 :
00219             histos[k].Scale(weight)
00220         #get Y min
00221         if histos[k].GetMinimum(0.) < minY :
00222             minY = histos[k].GetMinimum(0.)
00223         #get Y max
00224         if histos[k].GetBinContent(histos[k].GetMaximumBin()) > maxY :
00225             maxY = histos[k].GetBinContent(histos[k].GetMaximumBin())+histos[k].GetBinError(histos[k].GetMaximumBin())
00226         #Axis
00227         if plot.Xlabel :
00228             histos[k].SetXTitle(plot.Xlabel)
00229         if plot.Ylabel :
00230             histos[k].SetYTitle(plot.Ylabel)
00231         outhistos.append(histos[k])    
00232     #Range
00233     if not plot.logY : outhistos[0].GetYaxis().SetRangeUser(0,1.1*maxY)
00234     #else : outhistos[0].GetYaxis().SetRangeUser(0.0001,1.05)
00235     else : outhistos[0].GetYaxis().SetRangeUser(max(0.0001,0.5*minY),1.1*maxY)
00236     return outhistos        
00237 
00238 #method to do a plot from a graph
00239 def graphProducer(plot,histos,tagFlav="B",mistagFlav=["C","DUSG"],isVal=True):
00240     if histos is None : return
00241     if isVal : sample = "Val"
00242     else : sample = "Ref"
00243     #define graphs
00244     g = {}
00245     g_out = []
00246     if tagFlav not in listFlavors :
00247         return
00248     if plot.tagFlavor and plot.mistagFlavor :
00249         tagFlav = plot.tagFlavor
00250         mistagFlav = plot.mistagFlavor
00251     for f in listFlavors :
00252         #compute errors, in case not already done
00253         histos[f].Sumw2()
00254     #efficiency lists
00255     Eff = {}
00256     EffErr = {}
00257     for f in listFlavors :
00258         Eff[f] = []
00259         EffErr[f] = []
00260     #define mapping points for the histos
00261     maxnpoints = histos[tagFlav].GetNbinsX()
00262     for f in listFlavors :
00263         Eff[f].append(histos[f].GetBinContent(1))
00264         EffErr[f].append(histos[f].GetBinError(1))
00265     for bin in range(2,maxnpoints+1) :
00266         #check if we add the point to the graph for Val sample
00267         if len(Eff[tagFlav])>0 :
00268             delta = Eff[tagFlav][-1]-histos[tagFlav].GetBinContent(bin)
00269             if delta>max(0.005,EffErr[tagFlav][-1]) :
00270                 #get efficiencies
00271                 for f in listFlavors :
00272                     Eff[f].append(histos[f].GetBinContent(bin))
00273                     EffErr[f].append(histos[f].GetBinError(bin))
00274     #create TVector
00275     len_ = len(Eff[tagFlav])
00276     TVec_Eff = {}
00277     TVec_EffErr = {}
00278     for f in listFlavors :
00279         TVec_Eff[f] = TVectorD(len_)
00280         TVec_EffErr[f] = TVectorD(len_)
00281     #Fill the vector
00282     for j in range(0,len_) :
00283         for f in listFlavors :
00284             TVec_Eff[f][j] = Eff[f][j]
00285             TVec_EffErr[f][j] = EffErr[f][j]
00286     #fill TGraph
00287     for mis in mistagFlav :
00288         g[tagFlav+mis]=TGraphErrors(TVec_Eff[tagFlav],TVec_Eff[mis],TVec_EffErr[tagFlav],TVec_EffErr[mis])
00289     #style
00290     for f in listFlavors :
00291         if f not in mistagFlav : continue
00292         g[tagFlav+f].SetLineColor(mapColor[f])
00293         g[tagFlav+f].SetMarkerStyle(mapMarker[sample])
00294         g[tagFlav+f].SetMarkerColor(mapColor[f])
00295         g_out.append(g[tagFlav+f])
00296     index = -1     
00297     for g_i in g_out :
00298         index+=1
00299         if g_i is not None : break
00300     #Axis
00301     g_out[index].GetXaxis().SetRangeUser(0,1)
00302     g_out[index].GetYaxis().SetRangeUser(0.0001,1)
00303     if plot.Xlabel :
00304         g_out[index].GetXaxis().SetTitle(plot.Xlabel)
00305     if plot.Ylabel :
00306         g_out[index].GetYaxis().SetTitle(plot.Ylabel)
00307     #add in the list None for element in listFlavors for which no TGraph is computed
00308     for index,f in enumerate(listFlavors) :
00309         if f not in mistagFlav : g_out.insert(index,None)
00310     return g_out   
00311 
00312 #method to draw the plot and save it
00313 def savePlots(title,saveName,listFromats,plot,Histos,keyHisto,listLegend,options,ratios=None,legendName="") :
00314     #create canvas
00315     c = {}
00316     pads = {}
00317     if options.doRatio :
00318         c[keyHisto] = TCanvas(saveName,keyHisto+plot.title,700,700+24*len(listFlavors))
00319         pads["hist"] = TPad("hist", saveName+plot.title,0,0.11*len(listFlavors),1.0,1.0)    
00320     else :
00321         c[keyHisto] = TCanvas(keyHisto,saveName+plot.title,700,700)
00322         pads["hist"] = TPad("hist", saveName+plot.title,0,0.,1.0,1.0)
00323     pads["hist"].Draw()
00324     if ratios :
00325         for r in range(0,len(ratios)) :
00326             pads["ratio_"+str(r)] = TPad("ratio_"+str(r), saveName+plot.title+str(r),0,0.11*r,1.0,0.11*(r+1))
00327             pads["ratio_"+str(r)].Draw()
00328     pads["hist"].cd()
00329     #canvas style                                                                                                                                                                          
00330     if plot.logY : pads["hist"].SetLogy()
00331     if plot.grid : pads["hist"].SetGrid()
00332     #legend                                                                                                                                                                                
00333     leg = TLegend(0.6,0.4,0.8,0.6)
00334     leg.SetMargin(0.12)
00335     leg.SetTextSize(0.035)
00336     leg.SetFillColor(10)
00337     leg.SetBorderSize(0)
00338     #draw histos                                                                                                                                                                           
00339     first = True
00340     option = drawOption
00341     optionSame = drawOption+"same"
00342     if plot.doPerformance :
00343         option = "AP"
00344         optionSame = "sameP"
00345     for i in range(0,len(Histos)) :
00346         if Histos[i] is None : continue
00347         if first :
00348             if not plot.doPerformance : Histos[i].GetPainter().PaintStat(ROOT.gStyle.GetOptStat(),0)
00349             Histos[i].SetTitle(title)
00350             Histos[i].Draw(option)
00351             first = False
00352         else : Histos[i].Draw(optionSame)
00353         #Fill legend                                                                                                                                                                       
00354         if plot.legend and len(Histos)%len(listLegend)==0:
00355             r=len(Histos)/len(listLegend)
00356             index=i-r*len(listLegend)
00357             while(index<0):
00358                 index+=len(listLegend)
00359             legName = legendName.replace("KEY",listLegend[index])
00360             if i<len(listLegend) : legName = legName.replace("isVAL",options.ValRel)
00361             else : legName = legName.replace("isVAL",options.RefRel)
00362             if drawOption=="HIST" :
00363                 leg.AddEntry(Histos[i],legName,"L")
00364             else :
00365                 leg.AddEntry(Histos[i],legName,"P")
00366     #Draw legend                                                                                                                                                                           
00367     if plot.legend and options.drawLegend : leg.Draw("same")
00368     tex = None
00369     if options.printBanner :
00370         print type(options.printBanner)
00371         tex = TLatex(0.55,0.75,options.Banner)
00372         tex.SetNDC()
00373         tex.SetTextSize(0.05)
00374         tex.Draw()
00375     #save canvas
00376     if ratios :
00377         for r in range(0,len(ratios)) :
00378             pads["ratio_"+str(r)].cd()
00379             if ratios[r] is None : continue
00380             pads["ratio_"+str(r)].SetGrid()
00381             ratios[r].GetYaxis().SetTitle(listLegend[r]+"-jets")
00382             ratios[r].GetYaxis().SetTitleSize(0.15)
00383             ratios[r].GetYaxis().SetTitleOffset(0.2)
00384             ratios[r].GetYaxis().SetNdivisions(3,3,2)
00385             ratios[r].Draw("")
00386             unity.Draw("same")
00387     for format in listFromats :
00388         save = saveName+"."+format
00389         c[keyHisto].Print(save)
00390     return [c,leg,tex,pads]    
00391 #to create ratio plots from histograms
00392 def createRatio(hVal,hRef):
00393     ratio = []
00394     for h_i in range(0,len(hVal)): 
00395         if hVal[h_i] is None : continue
00396         r = TH1F(hVal[h_i].GetName()+"ratio","ratio "+hVal[h_i].GetTitle(),hVal[h_i].GetNbinsX(),hVal[h_i].GetXaxis().GetXmin(),hVal[h_i].GetXaxis().GetXmax())
00397         r.Add(hVal[h_i])
00398         r.Divide(hRef[h_i])
00399         r.GetYaxis().SetRangeUser(0.25,1.75)
00400         r.SetMarkerColor(hVal[h_i].GetMarkerColor())
00401         r.SetLineColor(hVal[h_i].GetLineColor())
00402         r.GetYaxis().SetLabelSize(0.15)
00403         r.GetXaxis().SetLabelSize(0.15)
00404         ratio.append(r)
00405     return ratio
00406 #to create ratio plots from TGraphErrors
00407 def createRatioFromGraph(hVal,hRef):
00408     ratio = []
00409     for g_i in range(0,len(hVal)):
00410         if hVal[g_i] is None :
00411             ratio.append(None)
00412             continue
00413         tmp = hVal[g_i].GetHistogram()
00414         histVal = TH1F(hVal[g_i].GetName()+"_ratio",hVal[g_i].GetTitle()+"_ratio",tmp.GetNbinsX(),tmp.GetXaxis().GetXmin(),tmp.GetXaxis().GetXmax())
00415         histRef = TH1F(hRef[g_i].GetName()+"_ratio",hRef[g_i].GetTitle()+"_ratio",histVal.GetNbinsX(),histVal.GetXaxis().GetXmin(),histVal.GetXaxis().GetXmax())
00416         #loop over the N points
00417         for p in range(0,hVal[g_i].GetN()-1):
00418             #get point p
00419             x = Double(0)
00420             y = Double(0)
00421             hVal[g_i].GetPoint(p,x,y)
00422             xerr = hVal[g_i].GetErrorX(p)
00423             yerr = hVal[g_i].GetErrorY(p)
00424             bin_p = histVal.FindBin(x)
00425             xHist = histVal.GetBinCenter(bin_p)
00426             #get the other point as xHist in [x,xbis]
00427             xbis = Double(0)
00428             ybis = Double(0)
00429             #points are odered from high x to low x
00430             if xHist>x : 
00431                 if p==0 : continue
00432                 xbiserr = hVal[g_i].GetErrorX(p-1)
00433                 ybiserr = hVal[g_i].GetErrorY(p-1)
00434                 hVal[g_i].GetPoint(p-1,xbis,ybis)
00435             else :
00436                 xbiserr = hVal[g_i].GetErrorX(p+1)
00437                 ybiserr = hVal[g_i].GetErrorY(p+1)
00438                 hVal[g_i].GetPoint(p+1,xbis,ybis)
00439             if ybis==y : 
00440                 #just take y at x
00441                 bin_p_valContent = y
00442                 bin_p_valContent_errP = y+yerr
00443                 bin_p_valContent_errM = y-yerr
00444             else :
00445                 #do a linear extrapolation (equivalent to do Eval(xHist))
00446                 a=(ybis-y)/(xbis-x)
00447                 b=y-a*x
00448                 bin_p_valContent = a*xHist+b
00449                 #extrapolate the error
00450                 aerrP = ( (ybis+ybiserr)-(y+yerr) ) / (xbis-x)
00451                 berrP = (y+yerr)-aerrP*x
00452                 bin_p_valContent_errP = aerrP*xHist+berrP
00453                 aerrM = ( (ybis-ybiserr)-(y-yerr) ) / (xbis-x)
00454                 berrM = (y-yerr)-aerrM*x
00455                 bin_p_valContent_errM = aerrM*xHist+berrM
00456             #fill val hist
00457             histVal.SetBinContent(bin_p,bin_p_valContent)
00458             histVal.SetBinError(bin_p,(bin_p_valContent_errP-bin_p_valContent_errM)/2)
00459             #loop over the reference TGraph to get the corresponding point
00460             for pRef in range(0,hRef[g_i].GetN()):
00461                 #get point pRef
00462                 xRef = Double(0)
00463                 yRef = Double(0)
00464                 hRef[g_i].GetPoint(pRef,xRef,yRef)
00465                 #take the first point as xRef < xHist
00466                 if xRef > xHist : continue
00467                 xReferr = hRef[g_i].GetErrorX(pRef)
00468                 yReferr = hRef[g_i].GetErrorY(pRef)
00469                 #get the other point as xHist in [xRef,xRefbis]
00470                 xRefbis = Double(0)
00471                 yRefbis = Double(0)
00472                 xRefbiserr = hRef[g_i].GetErrorX(pRef+1)
00473                 yRefbiserr = hRef[g_i].GetErrorY(pRef+1)
00474                 hRef[g_i].GetPoint(pRef+1,xRefbis,yRefbis)
00475                 if yRefbis==yRef :
00476                     #just take yRef at xRef
00477                     bin_p_refContent = yRef
00478                     bin_p_refContent_errP = yRef+yReferr
00479                     bin_p_refContent_errM = yRef-yReferr
00480                 else :
00481                     #do a linear extrapolation (equivalent to do Eval(xHist))
00482                     aRef=(ybis-y)/(xbis-x)
00483                     bRef=yRef-aRef*xRef
00484                     bin_p_refContent = aRef*xHist+bRef
00485                     #extrapolate the error
00486                     aReferrP = ((yRefbis+yRefbiserr)-(yRef+yReferr))/((xRefbis)-(xRef))
00487                     bReferrP = (yRef+yReferr)-aReferrP*(xRef-xReferr)
00488                     bin_p_refContent_errP = aReferrP*xHist+bReferrP
00489                     aReferrM = ((yRefbis-yRefbiserr)-(yRef-yReferr))/((xRefbis)-(xRef))
00490                     bReferrM = (yRef-yReferr)-aReferrM*(xRef+xReferr)
00491                     bin_p_refContent_errM = aReferrM*xHist+bReferrM
00492                 break
00493             #fill ref hist
00494             histRef.SetBinContent(bin_p,bin_p_refContent)
00495             histRef.SetBinError(bin_p,(bin_p_refContent_errP-bin_p_refContent_errM)/2)
00496         #do the ratio
00497         histVal.Sumw2()
00498         histRef.Sumw2()
00499         histVal.Divide(histRef)
00500         #ratio style
00501         histVal.GetXaxis().SetRangeUser(0.,1.)
00502         #histRef.GetXaxis().SetRangeUser(0.,1.)
00503         histVal.GetYaxis().SetRangeUser(0.25,1.75)
00504         histVal.SetMarkerColor(hVal[g_i].GetMarkerColor())
00505         histVal.SetLineColor(hVal[g_i].GetLineColor())
00506         histVal.GetYaxis().SetLabelSize(0.15)
00507         histVal.GetXaxis().SetLabelSize(0.15)
00508         ratio.append(histVal)
00509     return ratio
/afs/cern.ch/work/a/aaltunda/public/www/CMSSW_6_2_5/src/Validation/RecoB/scripts/histoStyle.py
