#!/usr/bin/env python
import sys
import math
from ROOT import gROOT, gStyle, TFile, TTree, TChain, TMVA, TCut, TCanvas, gDirectory, TH1, TGraph, gPad, TF1
import getopt

def print_usage():
    print "\nUsage: {0} <output file basename> <input ROOT file name>".format(sys.argv[0])
    print "Arguments: "
    print '\t-e: use this beam energy'
    print '\t-h: this help message'
    print

ebeam=1.056

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

# Parse the command line arguments
for opt, arg in options:
        if opt=='-e':
            ebeam=float(arg)
        if opt=='-h':
            print_usage()
            sys.exit(0)

molmass=math.sqrt(2.0*5.11e-4*ebeam)

if len(remainder)!=2:
    print_usage()
    sys.exit(0)

gROOT.SetBatch(True)
gStyle.SetOptStat(0)
gStyle.SetOptFit(1)
inFile = TFile(remainder[1])
#outFile = TFile(remainder[0]+".root","RECREATE")
events = inFile.Get("ntuple")


firstenergy=round(ebeam)*0.3
energystep=round(ebeam)*0.1

c = TCanvas("c","c",1200,900);
c.Print(remainder[0]+".pdf[")

events.Draw("topPY/topP:topPX/topP>>hnew(100,-0.05,0.05,100,0,0.1)","max(topMatchChisq,botMatchChisq)<5&&abs(topClT-botClT)<2&&abs(uncPX/uncP)<0.005&&abs(uncPY/uncP)<0.005&&(topHasL1&&botHasL1)&&bscChisq<5&&minIso>0.5&&max(topTrkChisq,botTrkChisq)<30&&abs(topP-0.4)<0.05","colz")
gDirectory.Get("hnew").SetTitle("Moller electron angles, 350-450 MeV")
gDirectory.Get("hnew").GetXaxis().SetTitle("theta_x relative to beam axis [rad]")
gDirectory.Get("hnew").GetYaxis().SetTitle("theta_y [rad]")
c.Print(remainder[0]+".pdf","Title:top_etheta")

events.Draw("topPY/topP:topPX/topP>>hnew(100,-0.05,0.05,100,0,0.1)","max(topMatchChisq,botMatchChisq)<5&&abs(topClT-botClT)<2&&abs(uncPX/uncP)<0.005&&abs(uncPY/uncP)<0.005&&(topHasL1&&botHasL1)&&bscChisq<5&&minIso>0.5&&max(topTrkChisq,botTrkChisq)<30&&abs(topP-0.4)<0.025","colz")
gDirectory.Get("hnew").SetTitle("Moller electron angles, 375-425 MeV")
gDirectory.Get("hnew").GetXaxis().SetTitle("theta_x relative to beam axis [rad]")
gDirectory.Get("hnew").GetYaxis().SetTitle("theta_y [rad]")
c.Print(remainder[0]+".pdf","Title:top_etheta")
#c.Print(remainder[0]+".pdf]")
#sys.exit(1)

c.Clear()
c.Divide(4,2)
for i in xrange(1,5):
    energy=firstenergy+(i-1)*energystep
    events.Draw("acos(topPZ/topP)-acos(1-0.511e-3*(1/topP-1/{0})):atan2(topPX,topPY)>>hnew(30,-1,1,50,-0.01,0.01)".format(ebeam),"abs(vzcM*{0}/vzcP-{1})<0.003&&abs(topTrkT-botTrkT)<3&&abs(uncPY/uncP)<0.005&&abs(uncPX/uncP)<0.005&&abs(topP-{2})<{3}".format(ebeam,molmass,energy,energystep/2),"goff,colz")
    hnew = gDirectory.Get("hnew")
    c.cd(i)
    hnew.DrawCopy("colz")
    hnew.FitSlicesY()
    hnew_1 = gDirectory.Get("hnew_1")
    c.cd(i+4)
    hnew_1.GetYaxis().SetRangeUser(-0.005,0.005)
    hnew_1.Fit("pol1","","",-0.4,0.4)
    hnew_1.DrawCopy()
c.Print(remainder[0]+".pdf","Title:top_etheta")

c.Clear()
c.Divide(4,2)
for i in xrange(1,5):
    energy=firstenergy+(i-1)*energystep
    events.Draw("acos(botPZ/botP)-acos(1-0.511e-3*(1/botP-1/{0})):atan2(botPX,-botPY)>>hnew(30,-1,1,50,-0.01,0.01)".format(ebeam),"abs(vzcM*{0}/vzcP-{1})<0.003&&abs(topTrkT-botTrkT)<3&&abs(uncPY/uncP)<0.005&&abs(uncPX/uncP)<0.005&&abs(botP-{2})<{3}".format(ebeam,molmass,energy,energystep/2),"goff,colz")
    hnew = gDirectory.Get("hnew")
    c.cd(i)
    hnew.DrawCopy("colz")
    hnew.FitSlicesY()
    hnew_1 = gDirectory.Get("hnew_1")
    c.cd(i+4)
    hnew_1.GetYaxis().SetRangeUser(-0.005,0.005)
    hnew_1.Fit("pol1","","",-0.4,0.4)
    hnew_1.DrawCopy()
c.Print(remainder[0]+".pdf","Title:bot_etheta")

c.Clear()
c.Divide(1,2)
c.cd(1)
events.Draw("topTrkZ0:topPY/topP>>hnew(100,0,0.1,100,-1,1)","abs(vzcM*{0}/vzcP-{1})<0.003&&abs(topTrkT-botTrkT)<3&&abs(uncPY/uncP)<0.005&&abs(uncPX/uncP)<0.005&&abs(topPX/topP)<0.01".format(ebeam,molmass),"colz")
hnew = gDirectory.Get("hnew")
hnew.FitSlicesY()
hnew_1 = gDirectory.Get("hnew_1")
c.cd(2)
hnew_1.Draw()
hnew_1.Fit("pol1","","",0.025,0.04)
hnew_1.GetYaxis().SetRangeUser(-0.4,0.4)
c.Print(remainder[0]+".pdf","Title:top_yz")

c.Clear()
c.Divide(1,2)
c.cd(1)
events.Draw("botTrkZ0:botPY/botP>>hnew(100,-0.1,0,100,-1,1)","abs(vzcM*{0}/vzcP-{1})<0.003&&abs(topTrkT-botTrkT)<3&&abs(uncPY/uncP)<0.005&&abs(uncPX/uncP)<0.005&&abs(botPX/botP)<0.01".format(ebeam,molmass),"colz")
hnew = gDirectory.Get("hnew")
hnew.FitSlicesY()
hnew_1 = gDirectory.Get("hnew_1")
c.cd(2)
hnew_1.Draw()
hnew_1.Fit("pol1","","",-0.04,-0.025)
hnew_1.GetYaxis().SetRangeUser(-0.4,0.4)
c.Print(remainder[0]+".pdf","Title:bot_yz")

c.Clear()
c.Divide(4,2)
for i in xrange(1,5):
    energy=firstenergy+(i-1)*energystep
    events.Draw("topTrkD0:topPX/topP>>hnew(100,-0.05,0.05,100,-2,2)","abs(vzcM*{0}/vzcP-{1})<0.003&&abs(topTrkT-botTrkT)<3&&abs(uncPY/uncP)<0.005&&abs(uncPX/uncP)<0.005&&abs(topP-{2})<{3}".format(ebeam,molmass,energy,energystep/2),"goff,colz")
    hnew = gDirectory.Get("hnew")
    c.cd(i)
    hnew.DrawCopy("colz")
    hnew.FitSlicesY()
    hnew_1 = gDirectory.Get("hnew_1")
    c.cd(i+4)
    hnew_1.GetYaxis().SetRangeUser(-0.5,0.5)
    hnew_1.Fit("pol1","","",-0.01,0.01)
    hnew_1.DrawCopy()
c.Print(remainder[0]+".pdf","Title:top_xz")

c.Clear()
c.Divide(4,2)
for i in xrange(1,5):
    energy=firstenergy+(i-1)*energystep
    events.Draw("botTrkD0:botPX/botP>>hnew(100,-0.05,0.05,100,-2,2)","abs(vzcM*{0}/vzcP-{1})<0.003&&abs(topTrkT-botTrkT)<3&&abs(uncPY/uncP)<0.005&&abs(uncPX/uncP)<0.005&&abs(botP-{2})<{3}".format(ebeam,molmass,energy,energystep/2),"goff,colz")
    hnew = gDirectory.Get("hnew")
    c.cd(i)
    hnew.DrawCopy("colz")
    hnew.FitSlicesY()
    hnew_1 = gDirectory.Get("hnew_1")
    c.cd(i+4)
    hnew_1.GetYaxis().SetRangeUser(-0.5,0.5)
    hnew_1.Fit("pol1","","",-0.01,0.01)
    hnew_1.DrawCopy()
c.Print(remainder[0]+".pdf","Title:bot_xz")

c.Clear()
c.Divide(1,2)
c.cd(1)
events.Draw("topTrkD0:topP>>hnew(100,0,{0},100,-1,1)".format(ebeam),"abs(vzcM*{0}/vzcP-{1})<0.003&&abs(topTrkT-botTrkT)<3&&abs(uncPY/uncP)<0.005&&abs(uncPX/uncP)<0.005&&abs(topPX/topP)<0.01".format(ebeam,molmass),"colz")
hnew = gDirectory.Get("hnew")
hnew.FitSlicesY()
hnew_1 = gDirectory.Get("hnew_1")
c.cd(2)
hnew_1.Draw()
hnew_1.Fit("pol1","","",ebeam/3,ebeam*2/3)
hnew_1.GetYaxis().SetRangeUser(-0.5,0.5)
c.Print(remainder[0]+".pdf","Title:top_xp")

c.Clear()
c.Divide(1,2)
c.cd(1)
events.Draw("botTrkD0:botP>>hnew(100,0,{0},100,-1,1)".format(ebeam),"abs(vzcM*{0}/vzcP-{1})<0.003&&abs(topTrkT-botTrkT)<3&&abs(uncPY/uncP)<0.005&&abs(uncPX/uncP)<0.005&&abs(botPX/botP)<0.01".format(ebeam,molmass),"colz")
hnew = gDirectory.Get("hnew")
hnew.FitSlicesY()
hnew_1 = gDirectory.Get("hnew_1")
c.cd(2)
hnew_1.Draw()
hnew_1.Fit("pol1","","",ebeam/3,ebeam*2/3)
hnew_1.GetYaxis().SetRangeUser(-0.5,0.5)
c.Print(remainder[0]+".pdf","Title:bot_xp")

c.Clear()
c.Divide(1,2)
c.cd(1)
events.Draw("(vzcM*{0}/vzcP):atan2(topPY-botPY,topPX-botPX)>>hnew(100,0,3,100,{1},{2})".format(ebeam,0.7*molmass,1.3*molmass),"abs(topTrkT-botTrkT)<3&&abs(uncPY/uncP)<0.005&&abs(uncPX/uncP)<0.005&&abs(topP-botP)<0.1","colz")
hnew = gDirectory.Get("hnew")
hnew.FitSlicesY()
hnew_1 = gDirectory.Get("hnew_1")
c.cd(2)
hnew_1.Draw()
hnew_1.Fit("pol1","","",1,2)
hnew_1.GetYaxis().SetRangeUser(0.9*molmass,1.1*molmass)
c.Print(remainder[0]+".pdf","Title:m_phi")

c.Clear()
c.Divide(1,2)
c.cd(1)
events.Draw("(vzcM*{0}/vzcP):topP-botP>>hnew(100,-1,1,100,{1},{2})".format(ebeam,0.7*molmass,1.3*molmass),"abs(topTrkT-botTrkT)<3&&abs(uncPY/uncP)<0.005&&abs(uncPX/uncP)<0.005&&abs(atan2(topPY-botPY,topPX-botPX)-TMath::Pi()/2)<0.1","colz")
hnew = gDirectory.Get("hnew")
hnew.FitSlicesY()
hnew_1 = gDirectory.Get("hnew_1")
c.cd(2)
hnew_1.Draw()
hnew_1.Fit("pol1","","",-0.4,0.4)
hnew_1.GetYaxis().SetRangeUser(0.9*molmass,1.1*molmass)
c.Print(remainder[0]+".pdf","Title:m_theta")

c.Clear()
c.Divide(1,2)
c.cd(1)
events.Draw("uncPY/uncP:uncPX/uncP>>hnew(100,-0.01,0.01,100,-0.01,0.01)","abs(topTrkT-botTrkT)<3&&abs(vzcM*{0}/vzcP-{1})<0.003".format(ebeam,molmass),"colz")
c.cd(2)
events.Draw("uncM:uncP>>(100,{0},{1},100,{2},{3})".format(0.9*ebeam,1.1*ebeam,0.8*molmass,1.2*molmass),"","colz")
c.Print(remainder[0]+".pdf","Title:vtx")

cb = TF1("cb","crystalball")

c.Clear()
c.Divide(2,2)
c.cd(1)
gPad.SetLogy()
events.Draw("uncM>>hnew(200,{0},{1})".format(0.5*molmass,2*molmass),"max(topMatchChisq,botMatchChisq)<5&&abs(topClT-botClT)<2&&abs(uncPX/uncP)<0.005&&abs(uncPY/uncP)<0.005","colz")
hnew = gDirectory.Get("hnew")
hnew.Fit("gaus","","",0.9*molmass,1.1*molmass)
gaus = hnew.GetFunction("gaus")
cb.SetParameters(gaus.GetParameter(0),gaus.GetParameter(1),gaus.GetParameter(2),-1.5,2.5)
hnew.Fit(cb,"L","",0.9*molmass,1.5*molmass)
gDirectory.Get("hnew").DrawCopy()
c.cd(2)
gPad.SetLogy()
events.Draw("vzcM>>hnew(200,{0},{1})".format(0.5*molmass,2*molmass),"max(topMatchChisq,botMatchChisq)<5&&abs(topClT-botClT)<2&&abs(uncPX/uncP)<0.005&&abs(uncPY/uncP)<0.005","colz")
hnew = gDirectory.Get("hnew")
hnew.Fit("gaus","","",0.9*molmass,1.1*molmass)
gaus = hnew.GetFunction("gaus")
cb.SetParameters(gaus.GetParameter(0),gaus.GetParameter(1),gaus.GetParameter(2),-1.5,2.5)
hnew.Fit(cb,"","",0.9*molmass,1.5*molmass)
gDirectory.Get("hnew").DrawCopy()
c.cd(3)
gPad.SetLogy()
events.Draw("uncM*{0}/uncP>>hnew(200,{1},{2})".format(ebeam,0.5*molmass,2*molmass),"max(topMatchChisq,botMatchChisq)<5&&abs(topClT-botClT)<2&&abs(uncPX/uncP)<0.005&&abs(uncPY/uncP)<0.005","colz")
hnew = gDirectory.Get("hnew")
hnew.Fit("gaus","","",0.9*molmass,1.1*molmass)
gaus = hnew.GetFunction("gaus")
cb.SetParameters(gaus.GetParameter(0),gaus.GetParameter(1),gaus.GetParameter(2),-1.5,2.5)
hnew.Fit(cb,"L","",0.9*molmass,1.5*molmass)
gDirectory.Get("hnew").DrawCopy()
c.cd(4)
gPad.SetLogy()
events.Draw("vzcM*{0}/vzcP>>hnew(200,{1},{2})".format(ebeam,0.5*molmass,2*molmass),"max(topMatchChisq,botMatchChisq)<5&&abs(topClT-botClT)<2&&abs(uncPX/uncP)<0.005&&abs(uncPY/uncP)<0.005","colz")
hnew = gDirectory.Get("hnew")
hnew.Fit("gaus","","",0.9*molmass,1.1*molmass)
gaus = hnew.GetFunction("gaus")
cb.SetParameters(gaus.GetParameter(0),gaus.GetParameter(1),gaus.GetParameter(2),-1.5,2.5)
hnew.Fit(cb,"L","",0.95*molmass,1.5*molmass)
gDirectory.Get("hnew").DrawCopy()
c.Print(remainder[0]+".pdf","Title:mass_ecalmatch")

c.Clear()
c.Divide(2,2)
c.cd(1)
gPad.SetLogy()
events.Draw("uncM>>hnew(200,{0},{1})".format(0.5*molmass,2*molmass),"abs(topTrkT-botTrkT)<3&&abs(uncPX/uncP)<0.005&&abs(uncPY/uncP)<0.005","colz")
hnew = gDirectory.Get("hnew")
hnew.Fit("gaus","","",0.9*molmass,1.1*molmass)
gaus = hnew.GetFunction("gaus")
cb.SetParameters(gaus.GetParameter(0),gaus.GetParameter(1),gaus.GetParameter(2),-1.5,2.5)
hnew.Fit(cb,"L","",0.9*molmass,1.5*molmass)
gDirectory.Get("hnew").DrawCopy()
c.cd(2)
gPad.SetLogy()
events.Draw("vzcM>>hnew(200,{0},{1})".format(0.5*molmass,2*molmass),"abs(topTrkT-botTrkT)<3&&abs(uncPX/uncP)<0.005&&abs(uncPY/uncP)<0.005","colz")
hnew = gDirectory.Get("hnew")
hnew.Fit("gaus","","",0.9*molmass,1.1*molmass)
gaus = hnew.GetFunction("gaus")
cb.SetParameters(gaus.GetParameter(0),gaus.GetParameter(1),gaus.GetParameter(2),-1.5,2.5)
hnew.Fit(cb,"L","",0.9*molmass,1.5*molmass)
gDirectory.Get("hnew").DrawCopy()
c.cd(3)
gPad.SetLogy()
events.Draw("uncM*{0}/uncP>>hnew(200,{1},{2})".format(ebeam,0.5*molmass,2*molmass),"abs(topTrkT-botTrkT)<3&&abs(uncPX/uncP)<0.005&&abs(uncPY/uncP)<0.005","colz")
hnew = gDirectory.Get("hnew")
hnew.Fit("gaus","","",0.9*molmass,1.1*molmass)
gaus = hnew.GetFunction("gaus")
cb.SetParameters(gaus.GetParameter(0),gaus.GetParameter(1),gaus.GetParameter(2),-1.5,2.5)
hnew.Fit(cb,"L","",0.9*molmass,1.5*molmass)
gDirectory.Get("hnew").DrawCopy()
c.cd(4)
gPad.SetLogy()
events.Draw("vzcM*{0}/vzcP>>hnew(200,{1},{2})".format(ebeam,0.5*molmass,2*molmass),"abs(topTrkT-botTrkT)<3&&abs(uncPX/uncP)<0.005&&abs(uncPY/uncP)<0.005","colz")
hnew = gDirectory.Get("hnew")
hnew.Fit("gaus","","",0.9*molmass,1.1*molmass)
gaus = hnew.GetFunction("gaus")
cb.SetParameters(gaus.GetParameter(0),gaus.GetParameter(1),gaus.GetParameter(2),-1.5,2.5)
hnew.Fit(cb,"L","",0.95*molmass,1.5*molmass)
gDirectory.Get("hnew").DrawCopy()
c.Print(remainder[0]+".pdf","Title:mass_noecalmatch")

c.Clear()
c.Divide(2,2)
c.cd(1)
gPad.SetLogy()
events.Draw("uncM>>hnew(200,{0},{1})".format(0.5*molmass,2*molmass),"abs(topTrkT-botTrkT)<3&&abs(uncPX/uncP)<0.005&&abs(uncPY/uncP)<0.005&&(topHasL1&&botHasL1)","colz")
hnew = gDirectory.Get("hnew")
hnew.Fit("gaus","","",0.9*molmass,1.1*molmass)
gaus = hnew.GetFunction("gaus")
cb.SetParameters(gaus.GetParameter(0),gaus.GetParameter(1),gaus.GetParameter(2),-1.5,2.5)
hnew.Fit(cb,"L","",0.9*molmass,1.5*molmass)
gDirectory.Get("hnew").DrawCopy()
c.cd(2)
gPad.SetLogy()
events.Draw("vzcM>>hnew(200,{0},{1})".format(0.5*molmass,2*molmass),"abs(topTrkT-botTrkT)<3&&abs(uncPX/uncP)<0.005&&abs(uncPY/uncP)<0.005&&(topHasL1&&botHasL1)","colz")
hnew = gDirectory.Get("hnew")
hnew.Fit("gaus","","",0.9*molmass,1.1*molmass)
gaus = hnew.GetFunction("gaus")
cb.SetParameters(gaus.GetParameter(0),gaus.GetParameter(1),gaus.GetParameter(2),-1.5,2.5)
hnew.Fit(cb,"L","",0.9*molmass,1.5*molmass)
gDirectory.Get("hnew").DrawCopy()
c.cd(3)
gPad.SetLogy()
events.Draw("uncM*{0}/uncP>>hnew(200,{1},{2})".format(ebeam,0.5*molmass,2*molmass),"abs(topTrkT-botTrkT)<3&&abs(uncPX/uncP)<0.005&&abs(uncPY/uncP)<0.005&&(topHasL1&&botHasL1)","colz")
hnew = gDirectory.Get("hnew")
hnew.Fit("gaus","","",0.9*molmass,1.1*molmass)
gaus = hnew.GetFunction("gaus")
cb.SetParameters(gaus.GetParameter(0),gaus.GetParameter(1),gaus.GetParameter(2),-1.5,2.5)
hnew.Fit(cb,"L","",0.9*molmass,1.5*molmass)
gDirectory.Get("hnew").DrawCopy()
c.cd(4)
gPad.SetLogy()
events.Draw("vzcM*{0}/vzcP>>hnew(200,{1},{2})".format(ebeam,0.5*molmass,2*molmass),"abs(topTrkT-botTrkT)<3&&abs(uncPX/uncP)<0.005&&abs(uncPY/uncP)<0.005&&(topHasL1&&botHasL1)","colz")
hnew = gDirectory.Get("hnew")
hnew.Fit("gaus","","",0.9*molmass,1.1*molmass)
gaus = hnew.GetFunction("gaus")
cb.SetParameters(gaus.GetParameter(0),gaus.GetParameter(1),gaus.GetParameter(2),-1.5,2.5)
hnew.Fit(cb,"L","",0.95*molmass,1.5*molmass)
gDirectory.Get("hnew").DrawCopy()
c.Print(remainder[0]+".pdf","Title:mass_noecalmatch_hasL1")

#c.Clear()
#c.Divide(2,2)
#c.cd(1)
#gPad.SetLogy()
#events.Draw("uncM>>hnew(200,{0},{1})".format(0.5*molmass,2*molmass),"abs(topTrkT-botTrkT)<3&&abs(uncPX/uncP)<0.005&&abs(uncPY/uncP)<0.005&&!(topHasL1&&botHasL1)","colz")
#hnew = gDirectory.Get("hnew")
#hnew.Fit("gaus","","",0.9*molmass,1.1*molmass)
#gaus = hnew.GetFunction("gaus")
#cb.SetParameters(gaus.GetParameter(0),gaus.GetParameter(1),gaus.GetParameter(2),-1.5,2.5)
#hnew.Fit(cb,"L","",0.9*molmass,1.5*molmass)
#gDirectory.Get("hnew").DrawCopy()
#c.cd(2)
#gPad.SetLogy()
#events.Draw("vzcM>>hnew(200,{0},{1})".format(0.5*molmass,2*molmass),"abs(topTrkT-botTrkT)<3&&abs(uncPX/uncP)<0.005&&abs(uncPY/uncP)<0.005&&!(topHasL1&&botHasL1)","colz")
#hnew = gDirectory.Get("hnew")
#hnew.Fit("gaus","","",0.9*molmass,1.1*molmass)
#gaus = hnew.GetFunction("gaus")
#cb.SetParameters(gaus.GetParameter(0),gaus.GetParameter(1),gaus.GetParameter(2),-1.5,2.5)
#hnew.Fit(cb,"L","",0.9*molmass,1.5*molmass)
#gDirectory.Get("hnew").DrawCopy()
#c.cd(3)
#gPad.SetLogy()
#events.Draw("uncM*{0}/uncP>>hnew(200,{1},{2})".format(ebeam,0.5*molmass,2*molmass),"abs(topTrkT-botTrkT)<3&&abs(uncPX/uncP)<0.005&&abs(uncPY/uncP)<0.005&&!(topHasL1&&botHasL1)","colz")
#hnew = gDirectory.Get("hnew")
#hnew.Fit("gaus","","",0.9*molmass,1.1*molmass)
#gaus = hnew.GetFunction("gaus")
#cb.SetParameters(gaus.GetParameter(0),gaus.GetParameter(1),gaus.GetParameter(2),-1.5,2.5)
#hnew.Fit(cb,"L","",0.9*molmass,1.5*molmass)
#gDirectory.Get("hnew").DrawCopy()
#c.cd(4)
#gPad.SetLogy()
#events.Draw("vzcM*{0}/vzcP>>hnew(200,{1},{2})".format(ebeam,0.5*molmass,2*molmass),"abs(topTrkT-botTrkT)<3&&abs(uncPX/uncP)<0.005&&abs(uncPY/uncP)<0.005&&!(topHasL1&&botHasL1)","colz")
#hnew = gDirectory.Get("hnew")
#hnew.Fit("gaus","","",0.9*molmass,1.1*molmass)
#gaus = hnew.GetFunction("gaus")
#cb.SetParameters(gaus.GetParameter(0),gaus.GetParameter(1),gaus.GetParameter(2),-1.5,2.5)
#hnew.Fit(cb,"L","",0.95*molmass,1.5*molmass)
#gDirectory.Get("hnew").DrawCopy()
#c.Print(remainder[0]+".pdf","Title:mass_noecalmatch_noL1")

c.Clear()
c.Divide(2,2)
c.cd(1)
gPad.SetLogy()
events.Draw("uncM>>hnew(200,{0},{1})".format(0.5*molmass,2*molmass),"max(topMatchChisq,botMatchChisq)<5&&abs(topClT-botClT)<2&&abs(uncPX/uncP)<0.005&&abs(uncPY/uncP)<0.005&&(topHasL1&&botHasL1)","colz")
hnew = gDirectory.Get("hnew")
hnew.Fit("gaus","","",0.9*molmass,1.1*molmass)
gaus = hnew.GetFunction("gaus")
cb.SetParameters(gaus.GetParameter(0),gaus.GetParameter(1),gaus.GetParameter(2),-1.5,2.5)
hnew.Fit(cb,"L","",0.9*molmass,1.5*molmass)
gDirectory.Get("hnew").DrawCopy()
c.cd(2)
gPad.SetLogy()
events.Draw("vzcM>>hnew(200,{0},{1})".format(0.5*molmass,2*molmass),"max(topMatchChisq,botMatchChisq)<5&&abs(topClT-botClT)<2&&abs(uncPX/uncP)<0.005&&abs(uncPY/uncP)<0.005&&(topHasL1&&botHasL1)","colz")
hnew = gDirectory.Get("hnew")
hnew.Fit("gaus","","",0.9*molmass,1.1*molmass)
gaus = hnew.GetFunction("gaus")
cb.SetParameters(gaus.GetParameter(0),gaus.GetParameter(1),gaus.GetParameter(2),-1.5,2.5)
hnew.Fit(cb,"","",0.9*molmass,1.5*molmass)
gDirectory.Get("hnew").DrawCopy()
c.cd(3)
gPad.SetLogy()
events.Draw("uncM*{0}/uncP>>hnew(200,{1},{2})".format(ebeam,0.5*molmass,2*molmass),"max(topMatchChisq,botMatchChisq)<5&&abs(topClT-botClT)<2&&abs(uncPX/uncP)<0.005&&abs(uncPY/uncP)<0.005&&(topHasL1&&botHasL1)","colz")
hnew = gDirectory.Get("hnew")
hnew.Fit("gaus","","",0.9*molmass,1.1*molmass)
gaus = hnew.GetFunction("gaus")
cb.SetParameters(gaus.GetParameter(0),gaus.GetParameter(1),gaus.GetParameter(2),-1.5,2.5)
hnew.Fit(cb,"L","",0.9*molmass,1.5*molmass)
gDirectory.Get("hnew").DrawCopy()
c.cd(4)
gPad.SetLogy()
events.Draw("vzcM*{0}/vzcP>>hnew(200,{1},{2})".format(ebeam,0.5*molmass,2*molmass),"max(topMatchChisq,botMatchChisq)<5&&abs(topClT-botClT)<2&&abs(uncPX/uncP)<0.005&&abs(uncPY/uncP)<0.005&&(topHasL1&&botHasL1)","colz")
hnew = gDirectory.Get("hnew")
hnew.Fit("gaus","","",0.9*molmass,1.1*molmass)
gaus = hnew.GetFunction("gaus")
cb.SetParameters(gaus.GetParameter(0),gaus.GetParameter(1),gaus.GetParameter(2),-1.5,2.5)
hnew.Fit(cb,"L","",0.95*molmass,1.5*molmass)
gDirectory.Get("hnew").DrawCopy()
c.Print(remainder[0]+".pdf","Title:mass_ecalmatch_hasl1")

c.Clear()
gPad.SetLogy()
events.Draw("uncM>>hnew(200,{0},{1})".format(0.5*molmass,2*molmass),"max(topMatchChisq,botMatchChisq)<5&&abs(topClT-botClT)<2&&abs(uncPX/uncP)<0.005&&abs(uncPY/uncP)<0.005&&(topHasL1&&botHasL1)&&bscChisq<5&&minIso>0.5&&abs(topP-botP)<0.4&&max(topTrkChisq,botTrkChisq)<30","colz")
hnew = gDirectory.Get("hnew")
hnew.SetTitle("Moller mass resolution")
hnew.GetXaxis().SetTitle("mass [GeV]")
hnew.GetYaxis().SetTitle("events/({0} GeV)".format(1.5*molmass/200))
hnew.Fit("gaus","","",0.9*molmass,1.1*molmass)
gaus = hnew.GetFunction("gaus")
cb.SetParameters(gaus.GetParameter(0),gaus.GetParameter(1),gaus.GetParameter(2),-1.5,2.5)
hnew.Fit(cb,"L","",0.9*molmass,1.5*molmass)
c.Print(remainder[0]+".pdf","Title:mollermass")

#c.Clear()
#c.Divide(2,2)
#c.cd(1)
#gPad.SetLogy()
#events.Draw("uncM>>hnew(200,{0},{1})".format(0.5*molmass,2*molmass),"max(topMatchChisq,botMatchChisq)<5&&abs(topClT-botClT)<2&&abs(uncPX/uncP)<0.005&&abs(uncPY/uncP)<0.005&&!(topHasL1&&botHasL1)","colz")
#hnew = gDirectory.Get("hnew")
#hnew.Fit("gaus","","",0.9*molmass,1.1*molmass)
#gaus = hnew.GetFunction("gaus")
#cb.SetParameters(gaus.GetParameter(0),gaus.GetParameter(1),gaus.GetParameter(2),-1.5,2.5)
#hnew.Fit(cb,"L","",0.9*molmass,1.5*molmass)
#gDirectory.Get("hnew").DrawCopy()
#c.cd(2)
#gPad.SetLogy()
#events.Draw("vzcM>>hnew(200,{0},{1})".format(0.5*molmass,2*molmass),"max(topMatchChisq,botMatchChisq)<5&&abs(topClT-botClT)<2&&abs(uncPX/uncP)<0.005&&abs(uncPY/uncP)<0.005&&!(topHasL1&&botHasL1)","colz")
#hnew = gDirectory.Get("hnew")
#hnew.Fit("gaus","","",0.9*molmass,1.1*molmass)
#gaus = hnew.GetFunction("gaus")
#cb.SetParameters(gaus.GetParameter(0),gaus.GetParameter(1),gaus.GetParameter(2),-1.5,2.5)
#hnew.Fit(cb,"","",0.9*molmass,1.5*molmass)
#gDirectory.Get("hnew").DrawCopy()
#c.cd(3)
#gPad.SetLogy()
#events.Draw("uncM*{0}/uncP>>hnew(200,{1},{2})".format(ebeam,0.5*molmass,2*molmass),"max(topMatchChisq,botMatchChisq)<5&&abs(topClT-botClT)<2&&abs(uncPX/uncP)<0.005&&abs(uncPY/uncP)<0.005&&!(topHasL1&&botHasL1)","colz")
#hnew = gDirectory.Get("hnew")
#hnew.Fit("gaus","","",0.9*molmass,1.1*molmass)
#gaus = hnew.GetFunction("gaus")
#cb.SetParameters(gaus.GetParameter(0),gaus.GetParameter(1),gaus.GetParameter(2),-1.5,2.5)
#hnew.Fit(cb,"L","",0.9*molmass,1.5*molmass)
#gDirectory.Get("hnew").DrawCopy()
#c.cd(4)
#gPad.SetLogy()
#events.Draw("vzcM*{0}/vzcP>>hnew(200,{1},{2})".format(ebeam,0.5*molmass,2*molmass),"max(topMatchChisq,botMatchChisq)<5&&abs(topClT-botClT)<2&&abs(uncPX/uncP)<0.005&&abs(uncPY/uncP)<0.005&&!(topHasL1&&botHasL1)","colz")
#hnew = gDirectory.Get("hnew")
#hnew.Fit("gaus","","",0.9*molmass,1.1*molmass)
#gaus = hnew.GetFunction("gaus")
#cb.SetParameters(gaus.GetParameter(0),gaus.GetParameter(1),gaus.GetParameter(2),-1.5,2.5)
#hnew.Fit(cb,"L","",0.95*molmass,1.5*molmass)
#gDirectory.Get("hnew").DrawCopy()
#c.Print(remainder[0]+".pdf)","Title:mass_ecalmatch_noL1")

c.Print(remainder[0]+".pdf]")

outfile = TFile(sys.argv[1]+".root","RECREATE")
events.Draw("vzcM>>mollers_noecalmatch(200,{0},{1})".format(0.5*molmass,2*molmass),"abs(topTrkT-botTrkT)<3&&abs(uncPX/uncP)<0.005&&abs(uncPY/uncP)<0.005","colz")
events.Draw("vzcM>>mollers_ecalmatch(200,{0},{1})".format(0.5*molmass,2*molmass),"max(topMatchChisq,botMatchChisq)<5&&abs(topClT-botClT)<2&&abs(uncPX/uncP)<0.005&&abs(uncPY/uncP)<0.005","colz")
events.Draw("vzcM>>mollers_ecalmatch_L1(200,{0},{1})".format(0.5*molmass,2*molmass),"max(topMatchChisq,botMatchChisq)<5&&abs(topClT-botClT)<2&&abs(uncPX/uncP)<0.005&&abs(uncPY/uncP)<0.005&&(topHasL1&&botHasL1)","colz")
events.Draw("vzcM>>mollers_noecalmatch_L1(200,{0},{1})".format(0.5*molmass,2*molmass),"abs(topTrkT-botTrkT)<3&&abs(uncPX/uncP)<0.005&&abs(uncPY/uncP)<0.005&&(topHasL1&&botHasL1)","colz")
outfile.Write()
outfile.Close()
