RecoB/scripts/histoStyle.py

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