drawUncertainties.py

Go to the documentation of this file.

import ROOT,array
from optparse import OptionParser
# Helper function that makes a error histogram out of a bunch of universes
# by calculating rms of universes in each bin
"""
usage = "usage: %prog[options]"
parser = OptionParser(usage=usage)

parser.add_option("-f","--filename",dest="filename",help="filename for histogram creation",default="NA")

(options,args)=parser.parse_args()
pathname =options.filename
"""

#Use different approach instead
import os, sys, shutil

if len(sys.argv)<2:
    print " drawUncertainties.py <filename1> [filename2]"
    sys.exit(1)
pathname = str(sys.argv[1])
pathname1=""
if len(sys.argv)>2:
    pathname1 = str(sys.argv[2])
def MakeErrorHistogram(universes,asFrac=True):
    #try the rebinning as well
    counter=0
    newbins=[]
    for num in  range(0,121): #using the 1 GeV uniform binning here
        newbins.append(num*1.0)
        counter+=1
    n_universes = len(universes)
    avg_hist = universes[0].Clone(universes[0].GetName()+"_avg")
    for i in range(0,avg_hist.GetNbinsX()+1):
        avg_hist.SetBinContent(i,0.0)
        avg_hist.SetBinError(i,0.0)
    avg_hist.Scale(1,"width")
    error_hist = universes[0].Clone(universes[0].GetName()+"_error")
    for i in range(0,error_hist.GetNbinsX()+1):
        error_hist.SetBinContent(i,0.0)
        error_hist.SetBinError(i,0.0)
    for i in range(n_universes):
        #universes[i].Print()
        universes[i].Scale(1,"width")
        avg_hist.Add(universes[i])
    avg_hist.Scale(1/float(n_universes))

    # Calculate Error = RMS of universes
    for j in range(error_hist.GetNbinsX()+2):
        error = 0
        for i in range(n_universes):
            error += 1/float(n_universes)*ROOT.TMath.Power(universes[i].GetBinContent(j)-avg_hist.GetBinContent(j),2)
        error_hist.SetBinContent(j,ROOT.TMath.Sqrt(error))
    if asFrac:
        error_hist.Divide(avg_hist)
    return error_hist
#Additional function to plot things easily

def AddUniverses(universes):
    n_universes=len(universes)
    avg_hist=universes[0].Clone(universes[0].GetName()+"_avg")
    for i in range(n_universes):
        avg_hist.Add(universes[i])
        avg_hist.Scale(1/float(n_universes))
    return avg_hist


def PrintPlots(rootfile,in_dir,in_files,universes,frac,urange):
    histo_universes=[]
    for i in range(universes):
        histo_universes.append(rootfile.Get(in_dir+"/"+in_files+"_"+str(i)))
    error_universes = MakeErrorHistogram(histo_universes,frac)
    c1 = ROOT.TCanvas("c1","c1",750,750)
    error_universes.SetLineWidth(3)
    if "htotal" in error_universes.GetName():
        print error_universes.GetName()[:-8],type(error_universes)
    if urange:
          (error_universes.GetXaxis()).SetRangeUser(0,20)
          error_universes.Draw("hist")
          c1.Print(str(in_files)+".pdf")
    else:
        error_universes.Draw("hist")
        c1.Print(str(in_files)+".pdf")
    del rootfile


def DrawFluxTogether(rootfile,in_dir,in_file,universes):
    histo_universes=[]
    for i in range(universes):
        histo_universes.append(rootfile.Get(in_dir+"/"+in_file+"_"+str(i)))
    c1 = ROOT.TCanvas("c1","c1",750,750)
    c1.cd()
    (histo_universes[0].GetXaxis()).SetRangeUser(0,20)
    maxm= histo_universes[0].GetMaximum()
    histo_universes[0].SetMaximum(maxm+0.1*maxm)
    histo_universes[0].Draw("hist")
    for i in range(0,universes):
        histo_universes[i].Draw("same")
    raw_input()



def RebinHistogram(input_histos,uniform=False):
    output_histos=[]
    counter=0
    newbins=array.array('d',[])
    if uniform:
        for num in range(0,17): #upto 8 GeV, 0.5 bins size
            newbins.append(num*0.5)
            counter+=1
        for num in range(5,11): #From 8 to 20 GeV, we will use 2 GeV bins size
            newbins.append(num*2)
            counter+=1
    if not uniform:
        for num in range(0,41):
            newbins.append(num*1.0)
            counter+=1

    for hist in input_histos:
        temphistname=hist.GetName()+"_rebinned"
        temphist=hist.Rebin(counter-1,temphistname,newbins)
        #if not drawing the fractional uncertainties, need to scale
        output_histos.append(temphist)
    return output_histos

def DrawTogether(rootfile,nutype,universes,frac,urange):
    errhists=[]
    directory=["total","others","targetpcpi","targetpck","targetncpi",
               "targetpcnu","targetmesonInc","targetnuA","totabsorption","tgtabsorption"]
    files=["htotal","hothers","htargetpcpi","htargetpcK","htargetncpi",
               "htargetpcnu","htargetmesonInc","htargetnucleonA","habsorption","htargetabsorption"]
    legends=["Total HP","others","pC-->#pi X","pC-->KX","nC-->#pi X","pC-->nucleonX","meson inc.","nucloen-A",
             "other abs.","Target Absorption"]
    color=[38,38,4,4,46,46,43,43,3,3]
    #style=[0,1,2,1,2,1,2,1,2,1,2]
    style=[0,1,1,2,1,2,1,2,1,2,1]
    for j in range(len(directory)):
        histo_universes=[]
        in_dir=nutype+"_"+directory[j]
        in_files=files[j]+"_"+nutype
        print in_dir+"/"+in_files
        for i in range(universes):
            histo_universes.append(rootfile.Get(str(in_dir)+"/"+str(in_files)+"_"+str(i)))
        rebinned_histos=RebinHistogram(histo_universes,False)
        error_universes = MakeErrorHistogram(rebinned_histos,frac)
        errhists.append(error_universes)
    #for i in range(0,errhists[4].GetNbinsX()+1):
        #print errhists[4].GetBinContent(i)
    #print len(errhists)   
    c1 = ROOT.TCanvas("c1","c1",750,750)
    leg = ROOT.TLegend(0.3,0.5,0.5,0.9)
    c1.cd()
    errhists[0].SetLineColor(1)
    errhists[0].SetLineWidth(3)
    errhists[0].SetMinimum(0)
    (errhists[0].GetXaxis()).SetRangeUser(0,20)
    leg.AddEntry(errhists[0],str(legends[0]),"lp")
    errhists[0].SetMaximum(0.20)
    errhists[0].SetMinimum(0)
    errhists[0].Draw("hist")
    #print len(errhists),len(color)
    for i in range(1,len(errhists)):
        errhists[i].SetLineColor(color[i])
        errhists[i].SetLineStyle(style[i])
        errhists[i].SetLineWidth(3)
        errhists[i].Draw("histsame")
        leg.AddEntry(errhists[i],str(legends[i]),"lp")
    leg.Draw("")
    #c1.Print("Interaction_info.pdf")
    raw_input()

 #Additional function only useful for comparing two uncertainty plots
 
def PrintStackPlots(rfile1,rfile2,in_dir,in_files,universes,frac,urange):
    if len(sys.argv)<2:
        print "need at least two files to print stack plots...."
        sys.exit(1)
    histo1_universes=[]
    histo2_universes=[]
    for i in range(universes):
        histo1_universes.append(rfile1.Get(in_dir+"/"+in_files+"_"+str(i)))
        histo2_universes.append(rfile2.Get(in_dir+"/"+in_files+"_"+str(i)))
    error_universes1 = MakeErrorHistogram(histo1_universes,frac)
    error_universes2 = MakeErrorHistogram(histo2_universes,frac)
    error_universes1.SetLineColor(1)
    error_universes2.SetLineColor(2)
    error_universes1.SetLineWidth(3)
    error_universes2.SetLineWidth(3)
    if urange:
        (error_universes1.GetXaxis()).SetRangeUser(0,20)
        (error_universes2.GetXaxis()).SetRangeUser(0,20)
        
    max1 = error_universes1.GetMaximum()
    max2 = error_universes2.GetMaximum()

    c1 = ROOT.TCanvas("c1","c1",750,750)
    if max1>max2:
        error_universes1.SetMinimum(0)
        error_universes1.Draw("hist")
        error_universes2.Draw("histSame")
    else:
        error_universes2.SetMinimum(0)
        error_universes2.Draw("hist")
        error_universes1.Draw("histSame")
    leg1 = ROOT.TLegend(0.5,0.7,0.7,0.9)
    #Assuming that error_universes1 will be for the reference beam here.
    leg1.AddEntry(error_universes1,"FTFP","lp")
    leg1.AddEntry(error_universes2,"QGSP","lp")
    leg1.Draw("")
    current=os.getcwd()
    newpath=current+"/comp_plots"
    if not os.path.exists(newpath):
        os.makedirs(newpath)
    c1.Print(str(newpath)+"/"+str(in_files)+"_comp.pdf")

    return



def PrintRatioErrorPlots(rfile1,rfile2,in_dir,in_files,universes,frac,urange,RPlots=False,rebin=True,uniform=True):
    if len(sys.argv)<2:
        print "need at least 2 files to compare the ratio error plots...exiting "
        sys.exit(1)
    print "Ratio and Error info for "+rfile1.GetName()+" / "+rfile2.GetName()
    
    histo1_universes=[]
    histo2_universes=[]
    ratio_universes=[]
    for i in range(universes):
        histo1_universes.append(rfile1.Get(in_dir+"/"+in_files+"_"+str(i)))
        histo2_universes.append(rfile2.Get(in_dir+"/"+in_files+"_"+str(i)))
    if not rebin:      
        for i in range(universes):
            histo1_universes[i].Divide(histo2_universes[i])
        error_universes = MakeErrorHistogram(histo1_universes,frac)
        ratio_universe=AddUniverses(histo1_universes)
    if rebin:
        unibin=True
        if not uniform:
            unibin=False
        rebin_histo1=RebinHistogram(histo1_universes,unibin)
        rebin_histo2=RebinHistogram(histo2_universes,unibin)
        for i in range(universes):
            rebin_histo1[i].Divide(rebin_histo2[i])
        error_universes=MakeErrorHistogram(rebin_histo1,frac)
        ratio_universe=AddUniverses(rebin_histo1)
        
    if urange:
      (error_universes.GetXaxis()).SetRangeUser(0,20)
      (ratio_universe.GetXaxis()).SetRangeUser(0,20)
    c1 = ROOT.TCanvas("c1","c1",750,750)
    leg1 = ROOT.TLegend(0.5,0.7,0.7,0.9)
    c1.cd()
    if RPlots:
        ratio_universe.SetLineWidth(3)
        ratio_universe.Draw("histe")
        raw_input()
        #c1.Print(str(in_files)+"_ratio.pdf")
        del c1
    else:
        
        error_universes.SetMinimum(0)
        error_universes.SetLineWidth(3)
        error_universes.Draw("hist")
        leg1.AddEntry(error_universes,"Error for Far/Near ")
        leg1.Draw("")
        raw_input()
        #c1.Print(str(in_files)+"_ratioerror.pdf")
# Open ppfx output file
in_file = ROOT.TFile(str(pathname))
rfile = ROOT.TFile(str(pathname1))
#PrintStackPlots(in_file,rfile,"numu_targetnuA","htargetnucleonA_numu",100,True,True)
#DrawFluxTogether(in_file,"numu_targetnuA","htargetnucleonA_numu",100)
#PrintPlots(in_file,"numu_targetncpi","htargetncpi_numu",100,True,True)
#PrintStackPlots(in_file,rfile,"numu_targetnuA","htargetnucleonA_numu",100,True,True)
#PrintPlots(in_file,"numu_attenuation","hatt_numu",100,True,True)
#PrintPlots(in_file,"numu_totabsorption","habsorption_numu",100,True,True)
#PrintPlots(in_file,"numu_targetncpi","htargetncpi_numu",100,True,True)
#PrintRatioErrorPlots(in_file,rfile,"numu_total","htotal_numu",100,True,True,False)