Back to home page

Project CMSSW displayed by LXR

 
 

    


File indexing completed on 2024-11-25 02:29:53

0001 #! /usr/bin/env python3
0002 #
0003 #
0004 # Francisco Yumiceva
0005 # yumiceva@fnal.gov
0006 #
0007 # Fermilab, 2010
0008 #
0009 #____________________________________________________________
0010 
0011 from builtins import range
0012 import sys
0013 import math
0014 import subprocess
0015 from array import array
0016 from ROOT import TCanvas, TFile, TGraphErrors
0017 
0018 def main():
0019 
0020     executable = "getbtagPerformance"
0021     rootFile = "rootFile=btag_performance_CSVM.root"
0022     payload  = "payload=BTAGCSVM" #"payload=MISTAGCSVM"
0023     OP = "CSVM"
0024     flavor   = "flavor=b" #l
0025     typeSF   = "type=SF"
0026     typeEff  = "type=eff"
0027     jetpt    = "pt="
0028     jeteta   = "eta="
0029 
0030     sp = " "
0031 
0032     ptmin = 25.
0033     ptmax = 240. #240.
0034     etamin = 0.
0035     etamax = 2.4
0036     ptNbins = 30
0037     etaNbins = 3
0038 
0039     ptbinwidth = (ptmax - ptmin)/ptNbins
0040     etabinwidth = (etamax - etamin)/etaNbins
0041     ptarray = array('d')
0042     pterrarray = array('d')
0043     etaarray = array('d')
0044     etaerrarray = array('d')
0045     SFarray = array('d')
0046     SFerrarray = array('d')
0047     effarray = array('d')
0048     efferrarray = array('d')
0049     MCeffarray = array('d')
0050     MCefferrarray = array('d')    
0051     SFarray_alleta = array('d')
0052     SFerrarray_alleta = array('d')
0053     effarray_alleta = array('d')
0054     efferrarray_alleta = array('d')
0055     MCeffarray_alleta = array('d')
0056     MCefferrarray_alleta = array('d')
0057         
0058 
0059     for ipt in range(0,ptNbins):
0060         SF_alleta = 0.
0061         SFerr_alleta = 0.
0062         Eff_alleta = 0.
0063         Efferr_alleta = 0.
0064         MCEff_alleta = 0.
0065         MCEfferr_alleta = 0.
0066                 
0067         apt = ptmin + ipt*ptbinwidth
0068         ptarray.append( apt )
0069         pterrarray.append( ptbinwidth*0.5 )
0070                     
0071         for ieta in range(0,etaNbins):
0072 
0073             aeta = etamin = ieta*etabinwidth
0074 
0075             print("pt = "+str(apt) + " eta = "+str(aeta))
0076             tmpjetpt = jetpt+str(apt)
0077             tmpjeteta = jeteta+str(aeta)
0078     
0079             allcmd = executable+ sp +rootFile+ sp +payload+ sp +flavor+ sp +typeSF+ sp +tmpjetpt+ sp +tmpjeteta
0080             print(allcmd)
0081             output = subprocess.getstatusoutput(allcmd)
0082             print(output[1])
0083             
0084             if output[0]!=0:
0085                 print(" Error retrieving data from file DB.")
0086 
0087             aSF = float( output[1].split()[2] )
0088             aSFerr = float( output[1].split()[4] )
0089 
0090             SFarray.append( aSF )
0091             SFerrarray.append( aSFerr )
0092             
0093             allcmd = executable+ sp +rootFile+ sp +payload+ sp +flavor+ sp +typeEff+ sp +tmpjetpt+ sp +tmpjeteta
0094             #print allcmd
0095             output = subprocess.getstatusoutput(allcmd)
0096             #print output[1]
0097             
0098             if output[0]!=0:
0099                 print(" Error retrieving data from file DB.")
0100                 
0101             aEff = float( output[1].split()[2] )
0102             aEfferr = float( output[1].split()[4] )
0103 
0104             effarray.append( aEff )
0105             efferrarray.append( aEfferr )
0106 
0107             MCeffarray.append( aEff/aSF )
0108             aMCEfferr = math.sqrt( (aSF*aSF*aEfferr*aEfferr + aEff*aEff*aSFerr*aSFerr)/(aSF*aSF*aSF*aSF) )
0109             MCefferrarray.append( aMCEfferr )
0110                         
0111             SF_alleta += aSF
0112             SFerr_alleta += aSF*aSF
0113             Eff_alleta += aEff
0114             Efferr_alleta += aEfferr*aEfferr
0115             MCEff_alleta += aEff/aSF
0116             MCEfferr_alleta += aMCEfferr*aMCEfferr
0117                                   
0118         SF_alleta = SF_alleta/etaNbins
0119         SFerr_alleta = math.sqrt(SFerr_alleta)/etaNbins
0120         Eff_alleta = Eff_alleta/etaNbins
0121         Efferr_alleta = math.sqrt(Efferr_alleta)/etaNbins
0122         MCEff_alleta = MCEff_alleta/etaNbins
0123         MCEfferr_alleta = math.sqrt(MCEfferr_alleta)/etaNbins
0124                 
0125         SFarray_alleta.append( SF_alleta )
0126         SFerrarray_alleta.append( SFerr_alleta )
0127         effarray_alleta.append( Eff_alleta )
0128         efferrarray_alleta.append( Efferr_alleta )
0129         MCeffarray_alleta.append( MCEff_alleta )
0130         MCefferrarray_alleta.append( MCEfferr_alleta )
0131                 
0132     histogram = {}
0133     histogram["SF_jet_pt"] = TGraphErrors(len(SFarray_alleta), ptarray, SFarray_alleta, pterrarray, SFerrarray_alleta)
0134     histogram["Eff_jet_pt"] = TGraphErrors(len(effarray_alleta), ptarray, effarray_alleta, pterrarray, efferrarray_alleta)
0135     histogram["MC_Eff_jet_pt"] = TGraphErrors(len(MCeffarray_alleta), ptarray, MCeffarray_alleta, pterrarray, MCefferrarray_alleta)
0136     
0137     cv1 = TCanvas("SF_jet_pt","SF_jet_pt",700,700)
0138     histogram["SF_jet_pt"].Draw("ap")
0139     histogram["SF_jet_pt"].SetName("SF_jet_pt")
0140     histogram["SF_jet_pt"].GetHistogram().SetXTitle("Jet p_{T} [GeV/c]")
0141     histogram["SF_jet_pt"].GetHistogram().SetYTitle("SF_{b}=#epsilon^{"+OP+"}_{data}/#epsilon^{"+OP+"}_{MC}")
0142     
0143     cv2 = TCanvas("Eff_jet_pt","Eff_jet_pt",700,700)
0144     histogram["Eff_jet_pt"].Draw("ap")
0145     histogram["Eff_jet_pt"].SetName("Eff_jet_pt")
0146     histogram["Eff_jet_pt"].GetHistogram().SetXTitle("Jet p_{T} [GeV/c]")
0147     histogram["Eff_jet_pt"].GetHistogram().SetYTitle("b-tag efficiency #epsilon^{"+OP+"}_{data}")
0148 
0149     cv3 = TCanvas("MC_Eff_jet_pt","MC_Eff_jet_pt",700,700)
0150     histogram["MC_Eff_jet_pt"].Draw("ap")
0151     histogram["MC_Eff_jet_pt"].SetName("MC_Eff_jet_pt")
0152     histogram["MC_Eff_jet_pt"].GetHistogram().SetXTitle("Jet p_{T} [GeV/c]")
0153     histogram["MC_Eff_jet_pt"].GetHistogram().SetYTitle("b-tag efficiency #epsilon^{"+OP+"}_{MC}")
0154 
0155                     
0156     outfilename = "plotPerformanceDB.root"
0157     outputroot = TFile( outfilename, "RECREATE")
0158         
0159     for key in histogram.keys():
0160         histogram[key].Write()
0161 
0162     outputroot.Close()
0163                     
0164     
0165     raw_input ("Enter to quit:")
0166      
0167 
0168 if __name__ == '__main__':
0169     main()
0170