RecoB/scripts/histoStyle.py

0001 from __future__ import print_function
0002
0003 #######
0004
0005 #  automatized plots generator for b-tagging performances
0006 #  Adrien Caudron, 2013, UCL
0007
0008 #######
0009
0010 #do all import
0011 import os, sys
0012
0013 try:
0014     import ROOT
0015 except:
0016     print("\nCannot load PYROOT, make sure you have setup ROOT in the path")
0017     print("and pyroot library is also defined in the variable PYTHONPATH, try:\n")
0018     if (os.getenv("PYTHONPATH")):
0019         print(" setenv PYTHONPATH ${PYTHONPATH}:$ROOTSYS/lib\n")
0020     else:
0021         print(" setenv PYTHONPATH $ROOTSYS/lib\n")
0022         sys.exit()
0023
0024 from ROOT import TFile
0025 from ROOT import TCanvas
0026 from ROOT import TPad
0027 from ROOT import TLegend
0028 from ROOT import TLatex
0029 from ROOT import TH1F
0030 from ROOT import TF1
0031 from ROOT import TVectorD
0032 from ROOT import TGraphErrors
0033 from ROOT import Double
0034
0035 import Style
0036 from listHistos import *
0037
0038 #define the input root files
0039 #fileVal = TFile("BTagRelVal_TTbar_Startup_14_612SLHC2.root","READ")
0040 #fileRef = TFile("BTagRelVal_TTbar_Startup_612SLHC1_14.root","READ")
0041 fileNameVal = "BTagRelVal_TTbar_Startup_14_612SLHC2.root"
0042 fileNameRef = "BTagRelVal_TTbar_Startup_612SLHC1_14.root"
0043 #define the val/ref labels
0044 ValRel = "612SLHC2_14"
0045 RefRel = "612SLHC1_14"
0046 #define the sample labels
0047 ValSample = "TTbar_FullSim"
0048 RefSample = "TTbar_FullSim"
0049 #define different settings
0050 batch = False #run on batch mode ?
0051 weight = 1. #rescale the histos according to this weight
0052 drawLegend = False #draw legend ?
0053 printBanner = False #draw text Banner on top of the histos
0054 Banner = "CMS Preliminary"
0055 doRatio = True #plot the ratios
0056 drawOption = "" # "" or "HIST"
0057 #path in the file
0058 pathInFile = "/DQMData/Run 1/Btag/Run summary/"
0059 #ETA/PT bins, GLOBAL ?
0060 EtaPtBin =[
0061     "GLOBAL",
0062     #"ETA_0-1v4",
0063     #"ETA_1v4-2v4",
0064     #"PT_50-80",
0065     #"PT_80-120",
0066     ]
0067 #list of taggers to look at
0068 listTag = [
0069     "CSV",
0070     "CSVMVA",
0071     "JP",
0072     "JBP",
0073     "TCHE",
0074     "TCHP",
0075     "SSVHE",
0076     "SSVHP",
0077     "SMT",
0078     #"SMTIP3d",
0079     #"SMTPt",
0080     "SET",
0081     ]
0082 #list of flavors to look at
0083 listFlavors = [
0084         #"ALL",
0085         "B",
0086         "C",
0087         #"G",
0088         #"DUS",
0089         "DUSG",
0090         #"NI",
0091         ]
0092 #map for marker color for flav-col and tag-col
0093 mapColor = {
0094     "ALL"  : 4 ,
0095     "B"    : 3 ,
0096     "C"    : 1 ,
0097     "G"    : 2 ,
0098     "DUS"  : 2 ,
0099     "DUSG" : 2 ,
0100     "NI"   : 5 ,
0101     "CSV"       : 5 ,
0102     "CSVMVA"   : 6 ,
0103     "JP"        : 3 ,
0104     "JBP"       : 9 ,
0105     "TCHE"      : 1,
0106     "TCHP"      : 2,
0107     "SSVHE"     : 4,
0108     "SSVHP"     : 7,
0109     "SMT"       : 8 ,
0110     "SMTIP3d" : 11 ,
0111     "SMTPt"   : 12
0112     }
0113 #marker style map for Val/Ref
0114 mapMarker = {
0115     "Val" : 22,
0116     "Ref" :  8
0117     }
0118 mapLineWidth = {
0119     "Val" : 3,
0120     "Ref" : 2
0121     }
0122 mapLineStyle = {
0123     "Val" : 2,
0124     "Ref" : 1
0125     }
0126 #choose the formats to save the plots
0127 listFromats = [
0128     "gif",
0129     ]
0130 #unity function
0131 unity = TF1("unity","1",-1000,1000)
0132 unity.SetLineColor(8)
0133 unity.SetLineWidth(1)
0134 unity.SetLineStyle(1)
0135 #list of histos to plots
0136 listHistos = [
0137     jetPt,
0138     jetEta,
0139     discr,
0140     effVsDiscrCut_discr,
0141     FlavEffVsBEff_discr,
0142     performance,
0143     #performanceC,
0144
0145     #IP,
0146     #IPe,
0147     #IPs,
0148     #NTracks,
0149     #decayLength,
0150     #distToJetAxis,
0151     #NHits,
0152     #NPixelHits,
0153     #NormChi2,
0154     #trackPt,
0155
0156
0157     #flightDist3Dval,
0158     #flightDist3Dsig,
0159     #jetNSecondaryVertices,
0160
0161     #vertexMass,
0162     #vertexNTracks,
0163     #vertexJetDeltaR,
0164     #vertexEnergyRatio,
0165
0166     #vertexCategory,
0167     #trackSip3dVal,
0168     #trackSip3dSig,
0169     #trackSip3dSigAboveCharm,
0170     #trackDeltaR,
0171     #trackEtaRel,
0172     #trackDecayLenVal,
0173     #trackSumJetDeltaR,
0174     #trackJetDist,
0175     #trackSumJetEtRatio,
0176     #trackPtRel,
0177     #trackPtRatio,
0178     #trackMomentum,
0179     #trackPPar,
0180     #trackPParRatio,
0181     ]
0182
0183 #methode to do a plot from histos
0184 def histoProducer(plot,histos,keys,isVal=True):
0185     if histos is None : return
0186     if isVal : sample = "Val"
0187     else : sample = "Ref"
0188     outhistos = []
0189     minY=9999.
0190     maxY=0.
0191     for k in keys :
0192         #Binning
0193         if plot.binning and len(plot.binning)==3 :
0194             histos[k].SetBins(plot.binning[0],plot.binning[1],plot.binning[2])
0195         elif plot.binning and len(plot.binning)==2 :
0196             nbins=plot.binning[1]+1-plot.binning[0]
0197             xmin=histos[k].GetBinLowEdge(plot.binning[0])
0198             xmax=histos[k].GetBinLowEdge(plot.binning[1]+1)
0199             valtmp=TH1F(histos[k].GetName(),histos[k].GetTitle(),nbins,xmin,xmax)
0200             i=1
0201             for bin in range(plot.binning[0],plot.binning[1]+1) :
0202                 valtmp.SetBinContent(i,histos[k].GetBinContent(bin))
0203                 i+=1
0204             histos[k]=valtmp
0205         if plot.Rebin and plot.Rebin > 0 :
0206             histos[k].Rebin(plot.Rebin)
0207         #Style
0208         histos[k].SetLineColor(mapColor[k])
0209         histos[k].SetMarkerColor(mapColor[k])
0210         histos[k].SetMarkerStyle(mapMarker[sample])
0211         if drawOption == "HIST" :
0212             histos[k].SetLineWidth(mapLineWidth[sample])
0213             histos[k].SetLineStyle(mapLineStyle[sample])
0214         #compute errors
0215         histos[k].Sumw2()
0216         #do the norm
0217         if plot.doNormalization :
0218             histos[k].Scale(1./histos[k].Integral())
0219         elif weight!=1 :
0220             histos[k].Scale(weight)
0221         #get Y min
0222         if histos[k].GetMinimum(0.) < minY :
0223             minY = histos[k].GetMinimum(0.)
0224         #get Y max
0225         if histos[k].GetBinContent(histos[k].GetMaximumBin()) > maxY :
0226             maxY = histos[k].GetBinContent(histos[k].GetMaximumBin())+histos[k].GetBinError(histos[k].GetMaximumBin())
0227         #Axis
0228         if plot.Xlabel :
0229             histos[k].SetXTitle(plot.Xlabel)
0230         if plot.Ylabel :
0231             histos[k].SetYTitle(plot.Ylabel)
0232         outhistos.append(histos[k])
0233     #Range
0234     if not plot.logY : outhistos[0].GetYaxis().SetRangeUser(0,1.1*maxY)
0235     #else : outhistos[0].GetYaxis().SetRangeUser(0.0001,1.05)
0236     else : outhistos[0].GetYaxis().SetRangeUser(max(0.0001,0.5*minY),1.1*maxY)
0237     return outhistos
0238
0239 #method to do a plot from a graph
0240 def graphProducer(plot,histos,tagFlav="B",mistagFlav=["C","DUSG"],isVal=True):
0241     if histos is None : return
0242     if isVal : sample = "Val"
0243     else : sample = "Ref"
0244     #define graphs
0245     g = {}
0246     g_out = []
0247     if tagFlav not in listFlavors :
0248         return
0249     if plot.tagFlavor and plot.mistagFlavor :
0250         tagFlav = plot.tagFlavor
0251         mistagFlav = plot.mistagFlavor
0252     for f in listFlavors :
0253         #compute errors, in case not already done
0254         histos[f].Sumw2()
0255     #efficiency lists
0256     Eff = {}
0257     EffErr = {}
0258     for f in listFlavors :
0259         Eff[f] = []
0260         EffErr[f] = []
0261     #define mapping points for the histos
0262     maxnpoints = histos[tagFlav].GetNbinsX()
0263     for f in listFlavors :
0264         Eff[f].append(histos[f].GetBinContent(1))
0265         EffErr[f].append(histos[f].GetBinError(1))
0266     for bin in range(2,maxnpoints+1) :
0267         #check if we add the point to the graph for Val sample
0268         if len(Eff[tagFlav])>0 :
0269             delta = Eff[tagFlav][-1]-histos[tagFlav].GetBinContent(bin)
0270             if delta>max(0.005,EffErr[tagFlav][-1]) :
0271                 #get efficiencies
0272                 for f in listFlavors :
0273                     Eff[f].append(histos[f].GetBinContent(bin))
0274                     EffErr[f].append(histos[f].GetBinError(bin))
0275     #create TVector
0276     len_ = len(Eff[tagFlav])
0277     TVec_Eff = {}
0278     TVec_EffErr = {}
0279     for f in listFlavors :
0280         TVec_Eff[f] = TVectorD(len_)
0281         TVec_EffErr[f] = TVectorD(len_)
0282     #Fill the vector
0283     for j in range(0,len_) :
0284         for f in listFlavors :
0285             TVec_Eff[f][j] = Eff[f][j]
0286             TVec_EffErr[f][j] = EffErr[f][j]
0287     #fill TGraph
0288     for mis in mistagFlav :
0289         g[tagFlav+mis]=TGraphErrors(TVec_Eff[tagFlav],TVec_Eff[mis],TVec_EffErr[tagFlav],TVec_EffErr[mis])
0290     #style
0291     for f in listFlavors :
0292         if f not in mistagFlav : continue
0293         g[tagFlav+f].SetLineColor(mapColor[f])
0294         g[tagFlav+f].SetMarkerStyle(mapMarker[sample])
0295         g[tagFlav+f].SetMarkerColor(mapColor[f])
0296         g_out.append(g[tagFlav+f])
0297     index = -1
0298     for g_i in g_out :
0299         index+=1
0300         if g_i is not None : break
0301     #Axis
0302     g_out[index].GetXaxis().SetRangeUser(0,1)
0303     g_out[index].GetYaxis().SetRangeUser(0.0001,1)
0304     if plot.Xlabel :
0305         g_out[index].GetXaxis().SetTitle(plot.Xlabel)
0306     if plot.Ylabel :
0307         g_out[index].GetYaxis().SetTitle(plot.Ylabel)
0308     #add in the list None for element in listFlavors for which no TGraph is computed
0309     for index,f in enumerate(listFlavors) :
0310         if f not in mistagFlav : g_out.insert(index,None)
0311     return g_out
0312
0313 #method to draw the plot and save it
0314 def savePlots(title,saveName,listFromats,plot,Histos,keyHisto,listLegend,options,ratios=None,legendName="") :
0315     #create canvas
0316     c = {}
0317     pads = {}
0318     if options.doRatio :
0319         c[keyHisto] = TCanvas(saveName,keyHisto+plot.title,700,700+24*len(listFlavors))
0320         pads["hist"] = TPad("hist", saveName+plot.title,0,0.11*len(listFlavors),1.0,1.0)
0321     else :
0322         c[keyHisto] = TCanvas(keyHisto,saveName+plot.title,700,700)
0323         pads["hist"] = TPad("hist", saveName+plot.title,0,0.,1.0,1.0)
0324     pads["hist"].Draw()
0325     if ratios :
0326         for r in range(0,len(ratios)) :
0327             pads["ratio_"+str(r)] = TPad("ratio_"+str(r), saveName+plot.title+str(r),0,0.11*r,1.0,0.11*(r+1))
0328             pads["ratio_"+str(r)].Draw()
0329     pads["hist"].cd()
0330     #canvas style
0331     if plot.logY : pads["hist"].SetLogy()
0332     if plot.grid : pads["hist"].SetGrid()
0333     #legend
0334     leg = TLegend(0.6,0.4,0.8,0.6)
0335     leg.SetMargin(0.12)
0336     leg.SetTextSize(0.035)
0337     leg.SetFillColor(10)
0338     leg.SetBorderSize(0)
0339     #draw histos
0340     first = True
0341     option = drawOption
0342     optionSame = drawOption+"same"
0343     if plot.doPerformance :
0344         option = "AP"
0345         optionSame = "sameP"
0346     for i in range(0,len(Histos)) :
0347         if Histos[i] is None : continue
0348         if first :
0349             if not plot.doPerformance : Histos[i].GetPainter().PaintStat(ROOT.gStyle.GetOptStat(),0)
0350             Histos[i].SetTitle(title)
0351             Histos[i].Draw(option)
0352             first = False
0353         else : Histos[i].Draw(optionSame)
0354         #Fill legend
0355         if plot.legend and len(Histos)%len(listLegend)==0:
0356             r=len(Histos)/len(listLegend)
0357             index=i-r*len(listLegend)
0358             while(index<0):
0359                 index+=len(listLegend)
0360             legName = legendName.replace("KEY",listLegend[index])
0361             if i<len(listLegend) : legName = legName.replace("isVAL",options.ValRel)
0362             else : legName = legName.replace("isVAL",options.RefRel)
0363             if drawOption=="HIST" :
0364                 leg.AddEntry(Histos[i],legName,"L")
0365             else :
0366                 leg.AddEntry(Histos[i],legName,"P")
0367     #Draw legend
0368     if plot.legend and options.drawLegend : leg.Draw("same")
0369     tex = None
0370     if options.printBanner :
0371         print(type(options.printBanner))
0372         tex = TLatex(0.55,0.75,options.Banner)
0373         tex.SetNDC()
0374         tex.SetTextSize(0.05)
0375         tex.Draw()
0376     #save canvas
0377     if ratios :
0378         for r in range(0,len(ratios)) :
0379             pads["ratio_"+str(r)].cd()
0380             if ratios[r] is None : continue
0381             pads["ratio_"+str(r)].SetGrid()
0382             ratios[r].GetYaxis().SetTitle(listLegend[r]+"-jets")
0383             ratios[r].GetYaxis().SetTitleSize(0.15)
0384             ratios[r].GetYaxis().SetTitleOffset(0.2)
0385             ratios[r].GetYaxis().SetNdivisions(3,3,2)
0386             ratios[r].Draw("")
0387             unity.Draw("same")
0388     for format in listFromats :
0389         save = saveName+"."+format
0390         c[keyHisto].Print(save)
0391     return [c,leg,tex,pads]
0392 #to create ratio plots from histograms
0393 def createRatio(hVal,hRef):
0394     ratio = []
0395     for h_i in range(0,len(hVal)):
0396         if hVal[h_i] is None : continue
0397         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())
0398         r.Add(hVal[h_i])
0399         r.Divide(hRef[h_i])
0400         r.GetYaxis().SetRangeUser(0.25,1.75)
0401         r.SetMarkerColor(hVal[h_i].GetMarkerColor())
0402         r.SetLineColor(hVal[h_i].GetLineColor())
0403         r.GetYaxis().SetLabelSize(0.15)
0404         r.GetXaxis().SetLabelSize(0.15)
0405         ratio.append(r)
0406     return ratio
0407 #to create ratio plots from TGraphErrors
0408 def createRatioFromGraph(hVal,hRef):
0409     ratio = []
0410     for g_i in range(0,len(hVal)):
0411         if hVal[g_i] is None :
0412             ratio.append(None)
0413             continue
0414         tmp = hVal[g_i].GetHistogram()
0415         histVal = TH1F(hVal[g_i].GetName()+"_ratio",hVal[g_i].GetTitle()+"_ratio",tmp.GetNbinsX(),tmp.GetXaxis().GetXmin(),tmp.GetXaxis().GetXmax())
0416         histRef = TH1F(hRef[g_i].GetName()+"_ratio",hRef[g_i].GetTitle()+"_ratio",histVal.GetNbinsX(),histVal.GetXaxis().GetXmin(),histVal.GetXaxis().GetXmax())
0417         #loop over the N points
0418         for p in range(0,hVal[g_i].GetN()-1):
0419             #get point p
0420             x = Double(0)
0421             y = Double(0)
0422             hVal[g_i].GetPoint(p,x,y)
0423             xerr = hVal[g_i].GetErrorX(p)
0424             yerr = hVal[g_i].GetErrorY(p)
0425             bin_p = histVal.FindBin(x)
0426             xHist = histVal.GetBinCenter(bin_p)
0427             #get the other point as xHist in [x,xbis]
0428             xbis = Double(0)
0429             ybis = Double(0)
0430             #points are odered from high x to low x
0431             if xHist>x :
0432                 if p==0 : continue
0433                 xbiserr = hVal[g_i].GetErrorX(p-1)
0434                 ybiserr = hVal[g_i].GetErrorY(p-1)
0435                 hVal[g_i].GetPoint(p-1,xbis,ybis)
0436             else :
0437                 xbiserr = hVal[g_i].GetErrorX(p+1)
0438                 ybiserr = hVal[g_i].GetErrorY(p+1)
0439                 hVal[g_i].GetPoint(p+1,xbis,ybis)
0440             if ybis==y :
0441                 #just take y at x
0442                 bin_p_valContent = y
0443                 bin_p_valContent_errP = y+yerr
0444                 bin_p_valContent_errM = y-yerr
0445             else :
0446                 #do a linear extrapolation (equivalent to do Eval(xHist))
0447                 a=(ybis-y)/(xbis-x)
0448                 b=y-a*x
0449                 bin_p_valContent = a*xHist+b
0450                 #extrapolate the error
0451                 aerrP = ( (ybis+ybiserr)-(y+yerr) ) / (xbis-x)
0452                 berrP = (y+yerr)-aerrP*x
0453                 bin_p_valContent_errP = aerrP*xHist+berrP
0454                 aerrM = ( (ybis-ybiserr)-(y-yerr) ) / (xbis-x)
0455                 berrM = (y-yerr)-aerrM*x
0456                 bin_p_valContent_errM = aerrM*xHist+berrM
0457             #fill val hist
0458             histVal.SetBinContent(bin_p,bin_p_valContent)
0459             histVal.SetBinError(bin_p,(bin_p_valContent_errP-bin_p_valContent_errM)/2)
0460             #loop over the reference TGraph to get the corresponding point
0461             for pRef in range(0,hRef[g_i].GetN()):
0462                 #get point pRef
0463                 xRef = Double(0)
0464                 yRef = Double(0)
0465                 hRef[g_i].GetPoint(pRef,xRef,yRef)
0466                 #take the first point as xRef < xHist
0467                 if xRef > xHist : continue
0468                 xReferr = hRef[g_i].GetErrorX(pRef)
0469                 yReferr = hRef[g_i].GetErrorY(pRef)
0470                 #get the other point as xHist in [xRef,xRefbis]
0471                 xRefbis = Double(0)
0472                 yRefbis = Double(0)
0473                 xRefbiserr = hRef[g_i].GetErrorX(pRef+1)
0474                 yRefbiserr = hRef[g_i].GetErrorY(pRef+1)
0475                 hRef[g_i].GetPoint(pRef+1,xRefbis,yRefbis)
0476                 if yRefbis==yRef :
0477                     #just take yRef at xRef
0478                     bin_p_refContent = yRef
0479                     bin_p_refContent_errP = yRef+yReferr
0480                     bin_p_refContent_errM = yRef-yReferr
0481                 else :
0482                     #do a linear extrapolation (equivalent to do Eval(xHist))
0483                     aRef=(ybis-y)/(xbis-x)
0484                     bRef=yRef-aRef*xRef
0485                     bin_p_refContent = aRef*xHist+bRef
0486                     #extrapolate the error
0487                     aReferrP = ((yRefbis+yRefbiserr)-(yRef+yReferr))/((xRefbis)-(xRef))
0488                     bReferrP = (yRef+yReferr)-aReferrP*(xRef-xReferr)
0489                     bin_p_refContent_errP = aReferrP*xHist+bReferrP
0490                     aReferrM = ((yRefbis-yRefbiserr)-(yRef-yReferr))/((xRefbis)-(xRef))
0491                     bReferrM = (yRef-yReferr)-aReferrM*(xRef+xReferr)
0492                     bin_p_refContent_errM = aReferrM*xHist+bReferrM
0493                 break
0494             #fill ref hist
0495             histRef.SetBinContent(bin_p,bin_p_refContent)
0496             histRef.SetBinError(bin_p,(bin_p_refContent_errP-bin_p_refContent_errM)/2)
0497         #do the ratio
0498         histVal.Sumw2()
0499         histRef.Sumw2()
0500         histVal.Divide(histRef)
0501         #ratio style
0502         histVal.GetXaxis().SetRangeUser(0.,1.)
0503         #histRef.GetXaxis().SetRangeUser(0.,1.)
0504         histVal.GetYaxis().SetRangeUser(0.25,1.75)
0505         histVal.SetMarkerColor(hVal[g_i].GetMarkerColor())
0506         histVal.SetLineColor(hVal[g_i].GetLineColor())
0507         histVal.GetYaxis().SetLabelSize(0.15)
0508         histVal.GetXaxis().SetLabelSize(0.15)
0509         ratio.append(histVal)
0510     return ratio