#!/usr/bin/env python
import sys, array
import getopt
from ROOT import gROOT, TCanvas, TF1, TFile, gStyle, TFormula, TGraph, TProfile, TGraphErrors, TH1D, TCutG, TH2D, TH2F, gDirectory, TLegend

options, remainder = getopt.gnu_getopt(sys.argv[1:], 'h', ['help',])

for opt, arg in options:
    if opt in ('-h', '--help'):
        print "\nUsage: "+sys.argv[0]+" <output basename> <root files>"
        print "Arguments: "
        print "\n"
        sys.exit(0)


#if (len(remainder)!=4):
#        print sys.argv[0]+' <output basename> <RAD> <tritrig> <WAB>'
#        sys.exit()

gStyle.SetOptStat(0)
c = TCanvas("c","c",800,600);
#print "\nUsage: {0} <output basename> <data ROOT file> <MC ROOT file>".format(sys.argv[0])
radFile_2 = TFile(remainder[1])
radFile_3 = TFile(remainder[2])
radFile_4 = TFile(remainder[3])
radFile_5 = TFile(remainder[4])
radFile_6 = TFile(remainder[5])
radFile_7 = TFile(remainder[6])
radFile_8 = TFile(remainder[7])
radFile_9 = TFile(remainder[8])
radFile_10 = TFile(remainder[9])
#tritrig
triFile_2 = TFile(remainder[10])
triFile_3 = TFile(remainder[11])
triFile_4 = TFile(remainder[12])
triFile_5 = TFile(remainder[13])
triFile_6 = TFile(remainder[14])
triFile_7 = TFile(remainder[15])
triFile_8 = TFile(remainder[16])
triFile_9 = TFile(remainder[17])
triFile_10 = TFile(remainder[18])
# wab
wabFile_2 = TFile(remainder[19])
wabFile_3 = TFile(remainder[20])
wabFile_4 = TFile(remainder[21])
wabFile_5 = TFile(remainder[22])
wabFile_6 = TFile(remainder[23])
wabFile_7 = TFile(remainder[24])
wabFile_8 = TFile(remainder[25])
wabFile_9 = TFile(remainder[26])
wabFile_10 = TFile(remainder[27])

radEvents_2 = radFile_2.Get("cut")
triEvents_2 = triFile_2.Get("cut")
wabEvents_2 = wabFile_2.Get("cut")
radEvents_3 = radFile_3.Get("cut")
triEvents_3 = triFile_3.Get("cut")
wabEvents_3 = wabFile_3.Get("cut")
radEvents_4 = radFile_4.Get("cut")
triEvents_4 = triFile_4.Get("cut")
wabEvents_4 = wabFile_4.Get("cut")
radEvents_5 = radFile_5.Get("cut")
triEvents_5 = triFile_5.Get("cut")
wabEvents_5 = wabFile_5.Get("cut")
radEvents_6 = radFile_6.Get("cut")
triEvents_6 = triFile_6.Get("cut")
wabEvents_6 = wabFile_6.Get("cut")
radEvents_7 = radFile_7.Get("cut")
triEvents_7 = triFile_7.Get("cut")
wabEvents_7 = wabFile_7.Get("cut")
radEvents_8 = radFile_8.Get("cut")
triEvents_8 = triFile_8.Get("cut")
wabEvents_8 = wabFile_8.Get("cut")
radEvents_9 = radFile_9.Get("cut")
triEvents_9 = triFile_9.Get("cut")
wabEvents_9 = wabFile_9.Get("cut")
radEvents_10 = radFile_10.Get("cut")
triEvents_10 = triFile_10.Get("cut")
wabEvents_10 = wabFile_10.Get("cut")

#nbins = 50
#width = 0.1
#scaling = nbins/width#500

# weight= (*/(n_gen/XS)*n_files)
#radEvents.SetWeight(500e-6/((10e4/0.154e6)*10))
#triEvents.SetWeight(500e-6/(0.0682*375))
#wabEvents.SetWeight(500e-6/((9.5e4/0.7e9)*995))
#1.056GeV
#radEvents.SetWeight(500e-6/((10e4/0.05565822e6)*99))

# 2016 rotationFix
#radEvents.SetWeight(500e-6/((10e4/0.023850e6/0.75962615429)*99))
#triEvents.SetWeight(500e-6/((10e4/0.495271980e6/0.75962615429)*100))
#wabEvents.SetWeight(500e-6/((15.785e4/0.163739790600e9/0.81367305033)*498))

# 2016 MG_alphaFix
# RAD
radEvents_2.SetWeight(500e-6/((10e4/0.2236e6)*10*100))
radEvents_3.SetWeight(500e-6/((10e4/0.1617e6)*10*99))
radEvents_4.SetWeight(500e-6/((10e4/0.1146e6)*10*99))
radEvents_5.SetWeight(500e-6/((10e4/0.08193e6)*10*99))
radEvents_6.SetWeight(500e-6/((10e4/0.06049e6)*10*98))
radEvents_7.SetWeight(500e-6/((10e4/0.04626e6)*10*98))
radEvents_8.SetWeight(500e-6/((10e4/0.03628e6)*10*99))
radEvents_9.SetWeight(500e-6/((10e4/0.02913e6)*10*100))
radEvents_10.SetWeight(500e-6/((10e4/0.02381e6)*10*100))
# tritrig
triEvents_2.SetWeight(500e-6/((10e4/4.690e6)*10*98))
triEvents_3.SetWeight(500e-6/((10e4/2.985e6)*10*98))
triEvents_4.SetWeight(500e-6/((10e4/2.018e6)*10*99))
triEvents_5.SetWeight(500e-6/((10e4/1.416e6)*10*98))
triEvents_6.SetWeight(500e-6/((10e4/1.027e6)*10*99))
triEvents_7.SetWeight(500e-6/((10e4/0.7652e6)*10*100))
triEvents_8.SetWeight(500e-6/((10e4/0.5832e6)*10*99))
triEvents_9.SetWeight(500e-6/((10e4/0.4533e6)*10*99))
triEvents_10.SetWeight(500e-6/((10e4/0.3592e6)*10*98))
# wab
wabEvents_2.SetWeight(500e-6/((10e4/1315.0e6)*10*994))
wabEvents_3.SetWeight(500e-6/((10e4/647.8e6)*10*989))
wabEvents_4.SetWeight(500e-6/((10e4/355.5e6)*10*989))
wabEvents_5.SetWeight(500e-6/((10e4/220.4e6)*10*990))
wabEvents_6.SetWeight(500e-6/((10e4/148.1e6)*10*989))
wabEvents_7.SetWeight(500e-6/((10e4/104.8e6)*10*988))
wabEvents_8.SetWeight(500e-6/((10e4/77.32e6)*10*992))
wabEvents_9.SetWeight(500e-6/((10e4/58.73e6)*10*989))
wabEvents_10.SetWeight(500e-6/((10e4/45.61e6)*10*984))

c.Print(remainder[0]+".pdf[")
outfile = TFile(remainder[0]+".root","RECREATE")

# 2mrad plots
radEvents_2.Draw("uncM>>rad_2mrad(50,0.0,0.2)","triPair1P>0.8*2.3","E1")
c.Print(remainder[0]+".pdf","Title:data_1d")
triEvents_2.Draw("uncM>>tri_2mrad(50,0.0,0.2)","","E1")
c.Print(remainder[0]+".pdf","Title:data_1d")
wabEvents_2.Draw("uncM>>wab_2mrad(50,0.0,0.2)","","E1")
c.Print(remainder[0]+".pdf","Title:data_1d")

radHist_2 = gDirectory.Get("rad_2mrad")
triHist_2 = gDirectory.Get("tri_2mrad")
wabHist_2 = gDirectory.Get("wab_2mrad")
totalHist_2 = triHist_2.Clone("total_2mrad")
totalHist_2.Add(wabHist_2)
totalHist_2.Draw()
c.Print(remainder[0]+".pdf","Title:data_1d")

totalHist_2.Draw()
totalHist_2.SetLineColor(1)
totalHist_2.SetTitle("Pair cross-sections 2mrad")
totalHist_2.GetXaxis().SetTitle("pair invariant mass [GeV]")
totalHist_2.GetYaxis().SetTitle("cross-section [mb/GeV]")
radHist_2.Draw("same")
radHist_2.SetLineColor(2)
wabHist_2.Draw("same")
wabHist_2.SetLineColor(3)
othertriHist_2 = triHist_2.Clone("othertri 2mrad")
othertriHist_2.Add(radHist_2,-1)
othertriHist_2.Draw("same")
othertriHist_2.SetLineColor(4)
leg = TLegend(0.5,0.75,0.9,0.9)
leg.AddEntry(totalHist_2,"Combined rate 2mrad")
leg.AddEntry(radHist_2,"Radiative tridents 2mrad")
leg.AddEntry(wabHist_2,"Wide-angle bremsstrahlung conversions 2mrad")
leg.AddEntry(othertriHist_2,"Non-radiative tridents 2mrad")
leg.Draw()
c.Print(remainder[0]+".pdf","Title:data_1d")

radfracHist_2 = radHist_2.Clone("radfrac_2mrad")
radfracHist_2.SetLineColor(1)
radfracHist_2.SetTitle("Radiative fraction (2mrad)")
radfracHist_2.GetXaxis().SetTitle("pair invariant mass [GeV]")
radfracHist_2.GetYaxis().SetTitle("radiative fraction")
radfracHist_2.Divide(totalHist_2)
radfracHist_2.Draw()
radfracHist_2.GetXaxis().SetRangeUser(0,0.17)
radfracHist_2.GetYaxis().SetRangeUser(0,0.5)
radfracHist_2.Fit("pol3","","",0.027,0.18)
c.Print(remainder[0]+".pdf","Title:data_1d")

# 3mrad plots
radEvents_3.Draw("uncM>>rad_3mrad(50,0.0,0.2)","triPair1P>0.8*2.3","E1")
c.Print(remainder[0]+".pdf","Title:data_1d")
triEvents_3.Draw("uncM>>tri_3mrad(50,0.0,0.2)","","E1")
c.Print(remainder[0]+".pdf","Title:data_1d")
wabEvents_3.Draw("uncM>>wab_3mrad(50,0.0,0.2)","","E1")
c.Print(remainder[0]+".pdf","Title:data_1d")

radHist_3 = gDirectory.Get("rad_3mrad")
triHist_3 = gDirectory.Get("tri_3mrad")
wabHist_3 = gDirectory.Get("wab_3mrad")
totalHist_3 = triHist_3.Clone("total_3mrad")
totalHist_3.Add(wabHist_3)
totalHist_3.Draw()
c.Print(remainder[0]+".pdf","Title:data_1d")

totalHist_3.Draw()
totalHist_3.SetLineColor(1)
totalHist_3.SetTitle("Pair cross-sections (3mrad)")
totalHist_3.GetXaxis().SetTitle("pair invariant mass [GeV]")
totalHist_3.GetYaxis().SetTitle("cross-section [mb/GeV]")
radHist_3.Draw("same")
radHist_3.SetLineColor(2)
wabHist_3.Draw("same")
wabHist_3.SetLineColor(3)
othertriHist_3 = triHist_3.Clone("othertri (3mrad)")
othertriHist_3.Add(radHist_3,-1)
othertriHist_3.Draw("same")
othertriHist_3.SetLineColor(4)
leg = TLegend(0.5,0.75,0.9,0.9)
leg.AddEntry(totalHist_3,"Combined rate (3mrad)")
leg.AddEntry(radHist_3,"Radiative tridents (3mrad)")
leg.AddEntry(wabHist_3,"Wide-angle bremsstrahlung conversions (3mrad)")
leg.AddEntry(othertriHist_3,"Non-radiative tridents (3mrad)")
leg.Draw()
c.Print(remainder[0]+".pdf","Title:data_1d")

radfracHist_3 = radHist_3.Clone("radfrac_3mrad")
radfracHist_3.SetLineColor(1)
radfracHist_3.SetTitle("Radiative fraction (3mrad)")
radfracHist_3.GetXaxis().SetTitle("pair invariant mass [GeV]")
radfracHist_3.GetYaxis().SetTitle("radiative fraction")
radfracHist_3.Divide(totalHist_3)
radfracHist_3.Draw()
radfracHist_3.GetXaxis().SetRangeUser(0,0.17)
radfracHist_3.GetYaxis().SetRangeUser(0,0.5)
radfracHist_3.Fit("pol3","","",0.027,0.18)


c.Print(remainder[0]+".pdf","Title:data_1d")


# 4mrad plots
radEvents_4.Draw("uncM>>rad_4mrad(50,0.0,0.2)","triPair1P>0.8*2.3","E1")
c.Print(remainder[0]+".pdf","Title:data_1d")
triEvents_4.Draw("uncM>>tri_4mrad(50,0.0,0.2)","","E1")
c.Print(remainder[0]+".pdf","Title:data_1d")
wabEvents_4.Draw("uncM>>wab_4mrad(50,0.0,0.2)","","E1")
c.Print(remainder[0]+".pdf","Title:data_1d")

radHist_4 = gDirectory.Get("rad_4mrad")
triHist_4 = gDirectory.Get("tri_4mrad")
wabHist_4 = gDirectory.Get("wab_4mrad")
totalHist_4 = triHist_4.Clone("total_4mrad")
totalHist_4.Add(wabHist_4)
totalHist_4.Draw()
c.Print(remainder[0]+".pdf","Title:data_1d")

totalHist_4.Draw()
totalHist_4.SetLineColor(1)
totalHist_4.SetTitle("Pair cross-sections (4mrad)")
totalHist_4.GetXaxis().SetTitle("pair invariant mass [GeV]")
totalHist_4.GetYaxis().SetTitle("cross-section [mb/GeV]")
radHist_4.Draw("same")
radHist_4.SetLineColor(2)
wabHist_4.Draw("same")
wabHist_4.SetLineColor(3)
othertriHist_4 = triHist_4.Clone("othertri (4mrad)")
othertriHist_4.Add(radHist_4,-1)
othertriHist_4.Draw("same")
othertriHist_4.SetLineColor(4)
leg = TLegend(0.5,0.75,0.9,0.9)
leg.AddEntry(totalHist_4,"Combined rate (4mrad)")
leg.AddEntry(radHist_4,"Radiative tridents (4mrad)")
leg.AddEntry(wabHist_4,"Wide-angle bremsstrahlung conversions (4mrad)")
leg.AddEntry(othertriHist_4,"Non-radiative tridents (4mrad)")
leg.Draw()
c.Print(remainder[0]+".pdf","Title:data_1d")

radfracHist_4 = radHist_4.Clone("radfrac_4mrad")
radfracHist_4.SetLineColor(1)
radfracHist_4.SetTitle("Radiative fraction (4mrad)")
radfracHist_4.GetXaxis().SetTitle("pair invariant mass [GeV]")
radfracHist_4.GetYaxis().SetTitle("radiative fraction")
radfracHist_4.Divide(totalHist_4)
radfracHist_4.Draw()
radfracHist_4.GetXaxis().SetRangeUser(0,0.17)
radfracHist_4.GetYaxis().SetRangeUser(0,0.5)
radfracHist_4.Fit("pol3","","",0.027,0.18)
c.Print(remainder[0]+".pdf","Title:data_1d")


# 5mrad plots
radEvents_5.Draw("uncM>>rad_5mrad(50,0.0,0.2)","triPair1P>0.8*2.3","E1")
c.Print(remainder[0]+".pdf","Title:data_1d")
triEvents_5.Draw("uncM>>tri_5mrad(50,0.0,0.2)","","E1")
c.Print(remainder[0]+".pdf","Title:data_1d")
wabEvents_5.Draw("uncM>>wab_5mrad(50,0.0,0.2)","","E1")
c.Print(remainder[0]+".pdf","Title:data_1d")

radHist_5 = gDirectory.Get("rad_5mrad")
triHist_5 = gDirectory.Get("tri_5mrad")
wabHist_5 = gDirectory.Get("wab_5mrad")
totalHist_5 = triHist_5.Clone("total_5mrad")
totalHist_5.Add(wabHist_5)
totalHist_5.Draw()
c.Print(remainder[0]+".pdf","Title:data_1d")

totalHist_5.Draw()
totalHist_5.SetLineColor(1)
totalHist_5.SetTitle("Pair cross-sections (5mrad)")
totalHist_5.GetXaxis().SetTitle("pair invariant mass [GeV]")
totalHist_5.GetYaxis().SetTitle("cross-section [mb/GeV]")
radHist_5.Draw("same")
radHist_5.SetLineColor(2)
wabHist_5.Draw("same")
wabHist_5.SetLineColor(3)
othertriHist_5 = triHist_5.Clone("othertri (5mrad)")
othertriHist_5.Add(radHist_5,-1)
othertriHist_5.Draw("same")
othertriHist_5.SetLineColor(4)
leg = TLegend(0.5,0.75,0.9,0.9)
leg.AddEntry(totalHist_5,"Combined rate (5mrad)")
leg.AddEntry(radHist_5,"Radiative tridents (5mrad)")
leg.AddEntry(wabHist_5,"Wide-angle bremsstrahlung conversions (5mrad)")
leg.AddEntry(othertriHist_5,"Non-radiative tridents (5mrad)")
leg.Draw()
c.Print(remainder[0]+".pdf","Title:data_1d")

radfracHist_5 = radHist_5.Clone("radfrac_5mrad")
radfracHist_5.SetLineColor(1)
radfracHist_5.SetTitle("Radiative fraction (5mrad)")
radfracHist_5.GetXaxis().SetTitle("pair invariant mass [GeV]")
radfracHist_5.GetYaxis().SetTitle("radiative fraction")
radfracHist_5.Divide(totalHist_5)
radfracHist_5.Draw()
radfracHist_5.GetXaxis().SetRangeUser(0,0.17)
radfracHist_5.GetYaxis().SetRangeUser(0,0.5)
radfracHist_5.Fit("pol3","","",0.027,0.18)
c.Print(remainder[0]+".pdf","Title:data_1d")

# 6mrad plots
radEvents_6.Draw("uncM>>rad_6mrad(50,0.0,0.2)","triPair1P>0.8*2.3","E1")
c.Print(remainder[0]+".pdf","Title:data_1d")
triEvents_6.Draw("uncM>>tri_6mrad(50,0.0,0.2)","","E1")
c.Print(remainder[0]+".pdf","Title:data_1d")
wabEvents_6.Draw("uncM>>wab_6mrad(50,0.0,0.2)","","E1")
c.Print(remainder[0]+".pdf","Title:data_1d")

radHist_6 = gDirectory.Get("rad_6mrad")
triHist_6 = gDirectory.Get("tri_6mrad")
wabHist_6 = gDirectory.Get("wab_6mrad")
totalHist_6 = triHist_6.Clone("total_6mrad")
totalHist_6.Add(wabHist_6)
totalHist_6.Draw()
c.Print(remainder[0]+".pdf","Title:data_1d")

totalHist_6.Draw()
totalHist_6.SetLineColor(1)
totalHist_6.SetTitle("Pair cross-sections (6mrad)")
totalHist_6.GetXaxis().SetTitle("pair invariant mass [GeV]")
totalHist_6.GetYaxis().SetTitle("cross-section [mb/GeV]")
radHist_6.Draw("same")
radHist_6.SetLineColor(2)
wabHist_6.Draw("same")
wabHist_6.SetLineColor(3)
othertriHist_6 = triHist_6.Clone("othertri (6mrad)")
othertriHist_6.Add(radHist_6,-1)
othertriHist_6.Draw("same")
othertriHist_6.SetLineColor(4)
leg = TLegend(0.5,0.75,0.9,0.9)
leg.AddEntry(totalHist_6,"Combined rate (6mrad)")
leg.AddEntry(radHist_6,"Radiative tridents (6mrad)")
leg.AddEntry(wabHist_6,"Wide-angle bremsstrahlung conversions (6mrad)")
leg.AddEntry(othertriHist_6,"Non-radiative tridents (6mrad)")
leg.Draw()
c.Print(remainder[0]+".pdf","Title:data_1d")

radfracHist_6 = radHist_6.Clone("radfrac_6mrad")
radfracHist_6.SetLineColor(1)
radfracHist_6.SetTitle("Radiative fraction (6mrad)")
radfracHist_6.GetXaxis().SetTitle("pair invariant mass [GeV]")
radfracHist_6.GetYaxis().SetTitle("radiative fraction")
radfracHist_6.Divide(totalHist_6)
radfracHist_6.Draw()
radfracHist_6.GetXaxis().SetRangeUser(0,0.17)
radfracHist_6.GetYaxis().SetRangeUser(0,0.5)
radfracHist_6.Fit("pol3","","",0.027,0.18)
c.Print(remainder[0]+".pdf","Title:data_1d")


# 7mrad plots
radEvents_7.Draw("uncM>>rad_7mrad(50,0.0,0.2)","triPair1P>0.8*2.3","E1")
c.Print(remainder[0]+".pdf","Title:data_1d")
triEvents_7.Draw("uncM>>tri_7mrad(50,0.0,0.2)","","E1")
c.Print(remainder[0]+".pdf","Title:data_1d")
wabEvents_7.Draw("uncM>>wab_7mrad(50,0.0,0.2)","","E1")
c.Print(remainder[0]+".pdf","Title:data_1d")

radHist_7 = gDirectory.Get("rad_7mrad")
triHist_7 = gDirectory.Get("tri_7mrad")
wabHist_7 = gDirectory.Get("wab_7mrad")
totalHist_7 = triHist_7.Clone("total_7mrad")
totalHist_7.Add(wabHist_7)
totalHist_7.Draw()
c.Print(remainder[0]+".pdf","Title:data_1d")

totalHist_7.Draw()
totalHist_7.SetLineColor(1)
totalHist_7.SetTitle("Pair cross-sections (7mrad)")
totalHist_7.GetXaxis().SetTitle("pair invariant mass [GeV]")
totalHist_7.GetYaxis().SetTitle("cross-section [mb/GeV]")
radHist_7.Draw("same")
radHist_7.SetLineColor(2)
wabHist_7.Draw("same")
wabHist_7.SetLineColor(3)
othertriHist_7 = triHist_7.Clone("othertri (7mrad)")
othertriHist_7.Add(radHist_7,-1)
othertriHist_7.Draw("same")
othertriHist_7.SetLineColor(4)
leg = TLegend(0.5,0.75,0.9,0.9)
leg.AddEntry(totalHist_7,"Combined rate (7mrad)")
leg.AddEntry(radHist_7,"Radiative tridents (7mrad)")
leg.AddEntry(wabHist_7,"Wide-angle bremsstrahlung conversions (7mrad)")
leg.AddEntry(othertriHist_7,"Non-radiative tridents (7mrad)")
leg.Draw()
c.Print(remainder[0]+".pdf","Title:data_1d")

radfracHist_7 = radHist_7.Clone("radfrac_7mrad")
radfracHist_7.SetLineColor(1)
radfracHist_7.SetTitle("Radiative fraction (7mrad)")
radfracHist_7.GetXaxis().SetTitle("pair invariant mass [GeV]")
radfracHist_7.GetYaxis().SetTitle("radiative fraction")
radfracHist_7.Divide(totalHist_7)
radfracHist_7.Draw()
radfracHist_7.GetXaxis().SetRangeUser(0,0.17)
radfracHist_7.GetYaxis().SetRangeUser(0,0.5)
radfracHist_7.Fit("pol3","","",0.027,0.18)
c.Print(remainder[0]+".pdf","Title:data_1d")

# 8mrad plots
radEvents_8.Draw("uncM>>rad_8mrad(50,0.0,0.2)","triPair1P>0.8*2.3","E1")
c.Print(remainder[0]+".pdf","Title:data_1d")
triEvents_8.Draw("uncM>>tri_8mrad(50,0.0,0.2)","","E1")
c.Print(remainder[0]+".pdf","Title:data_1d")
wabEvents_8.Draw("uncM>>wab_8mrad(50,0.0,0.2)","","E1")
c.Print(remainder[0]+".pdf","Title:data_1d")

radHist_8 = gDirectory.Get("rad_8mrad")
triHist_8 = gDirectory.Get("tri_8mrad")
wabHist_8 = gDirectory.Get("wab_8mrad")
totalHist_8 = triHist_8.Clone("total_8mrad")
totalHist_8.Add(wabHist_8)
totalHist_8.Draw()
c.Print(remainder[0]+".pdf","Title:data_1d")

totalHist_8.Draw()
totalHist_8.SetLineColor(1)
totalHist_8.SetTitle("Pair cross-sections (8mrad)")
totalHist_8.GetXaxis().SetTitle("pair invariant mass [GeV]")
totalHist_8.GetYaxis().SetTitle("cross-section [mb/GeV]")
radHist_8.Draw("same")
radHist_8.SetLineColor(2)
wabHist_8.Draw("same")
wabHist_8.SetLineColor(3)
othertriHist_8 = triHist_8.Clone("othertri (8mrad)")
othertriHist_8.Add(radHist_8,-1)
othertriHist_8.Draw("same")
othertriHist_8.SetLineColor(4)
leg = TLegend(0.5,0.75,0.9,0.9)
leg.AddEntry(totalHist_8,"Combined rate (8mrad)")
leg.AddEntry(radHist_8,"Radiative tridents (8mrad)")
leg.AddEntry(wabHist_8,"Wide-angle bremsstrahlung conversions (8mrad)")
leg.AddEntry(othertriHist_8,"Non-radiative tridents (8mrad)")
leg.Draw()
c.Print(remainder[0]+".pdf","Title:data_1d")

radfracHist_8 = radHist_8.Clone("radfrac_8mrad")
radfracHist_8.SetLineColor(1)
radfracHist_8.SetTitle("Radiative fraction (8mrad)")
radfracHist_8.GetXaxis().SetTitle("pair invariant mass [GeV]")
radfracHist_8.GetYaxis().SetTitle("radiative fraction")
radfracHist_8.Divide(totalHist_8)
radfracHist_8.Draw()
radfracHist_8.GetXaxis().SetRangeUser(0,0.17)
radfracHist_8.GetYaxis().SetRangeUser(0,0.5)
radfracHist_8.Fit("pol3","","",0.027,0.18)
c.Print(remainder[0]+".pdf","Title:data_1d")


# 9mrad plots
radEvents_9.Draw("uncM>>rad_9mrad(50,0.0,0.2)","triPair1P>0.8*2.3","E1")
c.Print(remainder[0]+".pdf","Title:data_1d")
triEvents_9.Draw("uncM>>tri_9mrad(50,0.0,0.2)","","E1")
c.Print(remainder[0]+".pdf","Title:data_1d")
wabEvents_9.Draw("uncM>>wab_9mrad(50,0.0,0.2)","","E1")
c.Print(remainder[0]+".pdf","Title:data_1d")

radHist_9 = gDirectory.Get("rad_9mrad")
triHist_9 = gDirectory.Get("tri_9mrad")
wabHist_9 = gDirectory.Get("wab_9mrad")
totalHist_9 = triHist_9.Clone("total_9mrad")
totalHist_9.Add(wabHist_9)
totalHist_9.Draw()
c.Print(remainder[0]+".pdf","Title:data_1d")

totalHist_9.Draw()
totalHist_9.SetLineColor(1)
totalHist_9.SetTitle("Pair cross-sections (9mrad)")
totalHist_9.GetXaxis().SetTitle("pair invariant mass [GeV]")
totalHist_9.GetYaxis().SetTitle("cross-section [mb/GeV]")
radHist_9.Draw("same")
radHist_9.SetLineColor(2)
wabHist_9.Draw("same")
wabHist_9.SetLineColor(3)
othertriHist_9 = triHist_9.Clone("othertri (9mrad)")
othertriHist_9.Add(radHist_9,-1)
othertriHist_9.Draw("same")
othertriHist_9.SetLineColor(4)
leg = TLegend(0.5,0.75,0.9,0.9)
leg.AddEntry(totalHist_9,"Combined rate (9mrad)")
leg.AddEntry(radHist_9,"Radiative tridents (9mrad)")
leg.AddEntry(wabHist_9,"Wide-angle bremsstrahlung conversions (9mrad)")
leg.AddEntry(othertriHist_9,"Non-radiative tridents (9mrad)")
leg.Draw()
c.Print(remainder[0]+".pdf","Title:data_1d")

radfracHist_9 = radHist_9.Clone("radfrac_9mrad")
radfracHist_9.SetLineColor(1)
radfracHist_9.SetTitle("Radiative fraction (9mrad)")
radfracHist_9.GetXaxis().SetTitle("pair invariant mass [GeV]")
radfracHist_9.GetYaxis().SetTitle("radiative fraction")
radfracHist_9.Divide(totalHist_9)
radfracHist_9.Draw()
radfracHist_9.GetXaxis().SetRangeUser(0,0.17)
radfracHist_9.GetYaxis().SetRangeUser(0,0.5)
radfracHist_9.Fit("pol3","","",0.027,0.18)
c.Print(remainder[0]+".pdf","Title:data_1d")

# 10mrad plots
radEvents_10.Draw("uncM>>rad_10mrad(50,0.0,0.2)","triPair1P>0.8*2.3","E1")
c.Print(remainder[0]+".pdf","Title:data_1d")
triEvents_10.Draw("uncM>>tri_10mrad(50,0.0,0.2)","","E1")
c.Print(remainder[0]+".pdf","Title:data_1d")
wabEvents_10.Draw("uncM>>wab_10mrad(50,0.0,0.2)","","E1")
c.Print(remainder[0]+".pdf","Title:data_1d")

radHist_10 = gDirectory.Get("rad_10mrad")
triHist_10 = gDirectory.Get("tri_10mrad")
wabHist_10 = gDirectory.Get("wab_10mrad")
totalHist_10 = triHist_9.Clone("total_10mrad")
totalHist_10.Add(wabHist_9)
totalHist_10.Draw()
c.Print(remainder[0]+".pdf","Title:data_1d")

totalHist_10.Draw()
totalHist_10.SetLineColor(1)
totalHist_10.SetTitle("Pair cross-sections (10mrad)")
totalHist_10.GetXaxis().SetTitle("pair invariant mass [GeV]")
totalHist_10.GetYaxis().SetTitle("cross-section [mb/GeV]")
radHist_10.Draw("same")
radHist_10.SetLineColor(2)
wabHist_10.Draw("same")
wabHist_10.SetLineColor(3)
othertriHist_10 = triHist_10.Clone("othertri (10mrad)")
othertriHist_10.Add(radHist_10,-1)
othertriHist_10.Draw("same")
othertriHist_10.SetLineColor(4)
leg = TLegend(0.5,0.75,0.9,0.9)
leg.AddEntry(totalHist_10,"Combined rate (10mrad)")
leg.AddEntry(radHist_10,"Radiative tridents (10mrad)")
leg.AddEntry(wabHist_10,"Wide-angle bremsstrahlung conversions (10mrad)")
leg.AddEntry(othertriHist_10,"Non-radiative tridents (10mrad)")
leg.Draw()
c.Print(remainder[0]+".pdf","Title:data_1d")

radfracHist_10 = radHist_10.Clone("radfrac_10mrad")
radfracHist_10.SetLineColor(1)
radfracHist_10.SetTitle("Radiative fraction (10mrad)")
radfracHist_10.GetXaxis().SetTitle("pair invariant mass [GeV]")
radfracHist_10.GetYaxis().SetTitle("radiative fraction")
radfracHist_10.Divide(totalHist_10)
radfracHist_10.Draw()
radfracHist_10.GetXaxis().SetRangeUser(0,0.17)
radfracHist_10.GetYaxis().SetRangeUser(0,0.5)
radfracHist_10.Fit("pol3","","",0.027,0.18)
c.Print(remainder[0]+".pdf","Title:data_1d")

radFracAt40 = TH2F("radFracAt40","RadFrac vs. theta cut (40 MeV)",9,1,10,10,0.065,0.1)
radFracAt50 = TH2F("radFracAt50","RadFrac vs. theta cut (50 MeV)",9,1,10,10,0.065,0.1)
radFracAt60 = TH2F("radFracAt60","RadFrac vs. theta cut (60 MeV)",9,1,10,10,0.065,0.1)
radFracAt70 = TH2F("radFracAt70","RadFrac vs. theta cut (70 MeV)",9,1,10,10,0.065,0.1)
radFracAt80 = TH2F("radFracAt80","RadFrac vs. theta cut (80 MeV)",9,1,10,10,0.065,0.1)
radFracAt90 = TH2F("radFracAt90","RadFrac vs. theta cut (90 MeV)",9,1,10,10,0.065,0.1)
radFracAt100 = TH2F("radFracAt100","RadFrac vs. theta cut (100 MeV)",9,1,10,10,0.065,0.1)
radFracAt110 = TH2F("radFracAt110","RadFrac vs. theta cut (110 MeV)",9,1,10,10,0.065,0.1)
radFracAt120 = TH2F("radFracAt120","RadFrac vs. theta cut (120 MeV)",9,1,10,10,0.065,0.1)

fit_40_2 = radfracHist_2.GetFunction("pol3").Eval(0.040)
fit_40_3 = radfracHist_3.GetFunction("pol3").Eval(0.040)
fit_40_4 = radfracHist_4.GetFunction("pol3").Eval(0.040)
fit_40_5 = radfracHist_5.GetFunction("pol3").Eval(0.040)
fit_40_6 = radfracHist_6.GetFunction("pol3").Eval(0.040)
fit_40_7 = radfracHist_7.GetFunction("pol3").Eval(0.040)
fit_40_8 = radfracHist_8.GetFunction("pol3").Eval(0.040)
fit_40_9 = radfracHist_9.GetFunction("pol3").Eval(0.040)
fit_40_10 = radfracHist_10.GetFunction("pol3").Eval(0.040)

fit_50_2 = radfracHist_2.GetFunction("pol3").Eval(0.050)
fit_50_3 = radfracHist_3.GetFunction("pol3").Eval(0.050)
fit_50_4 = radfracHist_4.GetFunction("pol3").Eval(0.050)
fit_50_5 = radfracHist_5.GetFunction("pol3").Eval(0.050)
fit_50_6 = radfracHist_6.GetFunction("pol3").Eval(0.050)
fit_50_7 = radfracHist_7.GetFunction("pol3").Eval(0.050)
fit_50_8 = radfracHist_8.GetFunction("pol3").Eval(0.050)
fit_50_9 = radfracHist_9.GetFunction("pol3").Eval(0.050)
fit_50_10 = radfracHist_10.GetFunction("pol3").Eval(0.050)

fit_60_2 = radfracHist_2.GetFunction("pol3").Eval(0.060)
fit_60_3 = radfracHist_3.GetFunction("pol3").Eval(0.060)
fit_60_4 = radfracHist_4.GetFunction("pol3").Eval(0.060)
fit_60_5 = radfracHist_5.GetFunction("pol3").Eval(0.060)
fit_60_6 = radfracHist_6.GetFunction("pol3").Eval(0.060)
fit_60_7 = radfracHist_7.GetFunction("pol3").Eval(0.060)
fit_60_8 = radfracHist_8.GetFunction("pol3").Eval(0.060)
fit_60_9 = radfracHist_9.GetFunction("pol3").Eval(0.060)
fit_60_10 = radfracHist_10.GetFunction("pol3").Eval(0.060)

fit_70_2 = radfracHist_2.GetFunction("pol3").Eval(0.070)
fit_70_3 = radfracHist_3.GetFunction("pol3").Eval(0.070)
fit_70_4 = radfracHist_4.GetFunction("pol3").Eval(0.070)
fit_70_5 = radfracHist_5.GetFunction("pol3").Eval(0.070)
fit_70_6 = radfracHist_6.GetFunction("pol3").Eval(0.070)
fit_70_7 = radfracHist_7.GetFunction("pol3").Eval(0.070)
fit_70_8 = radfracHist_8.GetFunction("pol3").Eval(0.070)
fit_70_9 = radfracHist_9.GetFunction("pol3").Eval(0.070)
fit_70_10 = radfracHist_10.GetFunction("pol3").Eval(0.070)

fit_80_2 = radfracHist_2.GetFunction("pol3").Eval(0.080)
fit_80_3 = radfracHist_3.GetFunction("pol3").Eval(0.080)
fit_80_4 = radfracHist_4.GetFunction("pol3").Eval(0.080)
fit_80_5 = radfracHist_5.GetFunction("pol3").Eval(0.080)
fit_80_6 = radfracHist_6.GetFunction("pol3").Eval(0.080)
fit_80_7 = radfracHist_7.GetFunction("pol3").Eval(0.080)
fit_80_8 = radfracHist_8.GetFunction("pol3").Eval(0.080)
fit_80_9 = radfracHist_9.GetFunction("pol3").Eval(0.080)
fit_80_10 = radfracHist_10.GetFunction("pol3").Eval(0.080)

fit_90_2 = radfracHist_2.GetFunction("pol3").Eval(0.090)
fit_90_3 = radfracHist_3.GetFunction("pol3").Eval(0.090)
fit_90_4 = radfracHist_4.GetFunction("pol3").Eval(0.090)
fit_90_5 = radfracHist_5.GetFunction("pol3").Eval(0.090)
fit_90_6 = radfracHist_6.GetFunction("pol3").Eval(0.090)
fit_90_7 = radfracHist_7.GetFunction("pol3").Eval(0.090)
fit_90_8 = radfracHist_8.GetFunction("pol3").Eval(0.090)
fit_90_9 = radfracHist_9.GetFunction("pol3").Eval(0.090)
fit_90_10 = radfracHist_10.GetFunction("pol3").Eval(0.090)

fit_100_2 = radfracHist_2.GetFunction("pol3").Eval(0.100)
fit_100_3 = radfracHist_3.GetFunction("pol3").Eval(0.100)
fit_100_4 = radfracHist_4.GetFunction("pol3").Eval(0.100)
fit_100_5 = radfracHist_5.GetFunction("pol3").Eval(0.100)
fit_100_6 = radfracHist_6.GetFunction("pol3").Eval(0.100)
fit_100_7 = radfracHist_7.GetFunction("pol3").Eval(0.100)
fit_100_8 = radfracHist_8.GetFunction("pol3").Eval(0.100)
fit_100_9 = radfracHist_9.GetFunction("pol3").Eval(0.100)
fit_100_10 = radfracHist_10.GetFunction("pol3").Eval(0.100)

fit_110_2 = radfracHist_2.GetFunction("pol3").Eval(0.110)
fit_110_3 = radfracHist_3.GetFunction("pol3").Eval(0.110)
fit_110_4 = radfracHist_4.GetFunction("pol3").Eval(0.110)
fit_110_5 = radfracHist_5.GetFunction("pol3").Eval(0.110)
fit_110_6 = radfracHist_6.GetFunction("pol3").Eval(0.110)
fit_110_7 = radfracHist_7.GetFunction("pol3").Eval(0.110)
fit_110_8 = radfracHist_8.GetFunction("pol3").Eval(0.110)
fit_110_9 = radfracHist_9.GetFunction("pol3").Eval(0.110)
fit_110_10 = radfracHist_10.GetFunction("pol3").Eval(0.110)

fit_120_2 = radfracHist_2.GetFunction("pol3").Eval(0.120)
fit_120_3 = radfracHist_3.GetFunction("pol3").Eval(0.120)
fit_120_4 = radfracHist_4.GetFunction("pol3").Eval(0.120)
fit_120_5 = radfracHist_5.GetFunction("pol3").Eval(0.120)
fit_120_6 = radfracHist_6.GetFunction("pol3").Eval(0.120)
fit_120_7 = radfracHist_7.GetFunction("pol3").Eval(0.120)
fit_120_8 = radfracHist_8.GetFunction("pol3").Eval(0.120)
fit_120_9 = radfracHist_9.GetFunction("pol3").Eval(0.120)
fit_120_10 = radfracHist_10.GetFunction("pol3").Eval(0.120)

radFracAt40.SetMarkerStyle(2)
radFracAt40.SetMarkerColor(2)
radFracAt40.Fill(2,fit_40_2)
radFracAt40.Fill(3,fit_40_3)
radFracAt40.Fill(4,fit_40_4)
radFracAt40.Fill(5,fit_40_5)
radFracAt40.Fill(6,fit_40_6)
radFracAt40.Fill(7,fit_40_7)
radFracAt40.Fill(8,fit_40_8)
radFracAt40.Fill(9,fit_40_9)
radFracAt40.Fill(10,fit_40_10)

radFracAt50.SetMarkerStyle(3)
radFracAt50.SetMarkerColor(3)
radFracAt50.Fill(2,fit_50_2)
radFracAt50.Fill(3,fit_50_3)
radFracAt50.Fill(4,fit_50_4)
radFracAt50.Fill(5,fit_50_5)
radFracAt50.Fill(6,fit_50_6)
radFracAt50.Fill(7,fit_50_7)
radFracAt50.Fill(8,fit_50_8)
radFracAt50.Fill(9,fit_50_9)
radFracAt50.Fill(10,fit_50_10)

radFracAt60.SetMarkerStyle(4)
radFracAt60.SetMarkerColor(4)
radFracAt60.Fill(2,fit_60_2)
radFracAt60.Fill(3,fit_60_3)
radFracAt60.Fill(4,fit_60_4)
radFracAt60.Fill(5,fit_60_5)
radFracAt60.Fill(6,fit_60_6)
radFracAt60.Fill(7,fit_60_7)
radFracAt60.Fill(8,fit_60_8)
radFracAt60.Fill(9,fit_60_9)
radFracAt60.Fill(10,fit_60_10)

radFracAt70.SetMarkerStyle(21)
radFracAt70.SetMarkerColor(30)
radFracAt70.Fill(2,fit_70_2)
radFracAt70.Fill(3,fit_70_3)
radFracAt70.Fill(4,fit_70_4)
radFracAt70.Fill(5,fit_70_5)
radFracAt70.Fill(6,fit_70_6)
radFracAt70.Fill(7,fit_70_7)
radFracAt70.Fill(8,fit_70_8)
radFracAt70.Fill(9,fit_70_9)
radFracAt70.Fill(10,fit_70_10)

radFracAt80.SetMarkerStyle(29)
radFracAt80.SetMarkerColor(6)
radFracAt80.Fill(2,fit_80_2)
radFracAt80.Fill(3,fit_80_3)
radFracAt80.Fill(4,fit_80_4)
radFracAt80.Fill(5,fit_80_5)
radFracAt80.Fill(6,fit_80_6)
radFracAt80.Fill(7,fit_80_7)
radFracAt80.Fill(8,fit_80_8)
radFracAt80.Fill(9,fit_80_9)
radFracAt80.Fill(10,fit_80_10)

radFracAt90.SetMarkerStyle(33)
radFracAt90.SetMarkerColor(7)
radFracAt90.Fill(2,fit_90_2)
radFracAt90.Fill(3,fit_90_3)
radFracAt90.Fill(4,fit_90_4)
radFracAt90.Fill(5,fit_90_5)
radFracAt90.Fill(6,fit_90_6)
radFracAt90.Fill(7,fit_90_7)
radFracAt90.Fill(8,fit_90_8)
radFracAt90.Fill(9,fit_90_9)
radFracAt90.Fill(10,fit_90_10)

radFracAt100.SetMarkerStyle(8)
radFracAt100.SetMarkerColor(8)
radFracAt100.Fill(2,fit_100_2)
radFracAt100.Fill(3,fit_100_3)
radFracAt100.Fill(4,fit_100_4)
radFracAt100.Fill(5,fit_100_5)
radFracAt100.Fill(6,fit_100_6)
radFracAt100.Fill(7,fit_100_7)
radFracAt100.Fill(8,fit_100_8)
radFracAt100.Fill(9,fit_100_9)
radFracAt100.Fill(10,fit_100_10)

radFracAt110.SetMarkerStyle(34)
radFracAt110.SetMarkerColor(9)
radFracAt110.Fill(2,fit_110_2)
radFracAt110.Fill(3,fit_110_3)
radFracAt110.Fill(4,fit_110_4)
radFracAt110.Fill(5,fit_110_5)
radFracAt110.Fill(6,fit_110_6)
radFracAt110.Fill(7,fit_110_7)
radFracAt110.Fill(8,fit_110_8)
radFracAt110.Fill(9,fit_110_9)
radFracAt110.Fill(10,fit_110_10)

radFracAt120.SetMarkerStyle(22)
radFracAt120.SetMarkerColor(1)
radFracAt120.Fill(2,fit_120_2)
radFracAt120.Fill(3,fit_120_3)
radFracAt120.Fill(4,fit_120_4)
radFracAt120.Fill(5,fit_120_5)
radFracAt120.Fill(6,fit_120_6)
radFracAt120.Fill(7,fit_120_7)
radFracAt120.Fill(8,fit_120_8)
radFracAt120.Fill(9,fit_120_9)
radFracAt120.Fill(10,fit_120_10)

leg = TLegend(0.5,0.75,0.9,0.9)
leg.AddEntry(radFracAt40,"40 MeV")
leg.AddEntry(radFracAt50,"50 MeV")
leg.AddEntry(radFracAt60,"60 MeV")
leg.AddEntry(radFracAt70,"70 MeV")
leg.AddEntry(radFracAt80,"80 MeV")
leg.AddEntry(radFracAt90,"90 MeV")
leg.AddEntry(radFracAt100,"100 MeV")
leg.AddEntry(radFracAt110,"110 MeV")
leg.AddEntry(radFracAt120,"120 MeV")

radFracAt40.Draw()
radFracAt50.Draw("same")
radFracAt60.Draw("same")
radFracAt70.Draw("same")
radFracAt80.Draw("same")
radFracAt90.Draw("same")
radFracAt100.Draw("same")
radFracAt110.Draw("same")
radFracAt120.Draw("same")

leg.Draw()

c.Print(remainder[0]+".pdf]","Title:angleScan")

outfile.Write()
outfile.Close()

