Burg jeff E1 E2.py

import mytools2 as mt2
import mytools as mt
import ROOT
import numpy as np
from TH1Wrapper import TH2F,TH3F

target = "D2O"
forward = 1
particle = "neutrons"
binwidth=1.25
delta_thetaperp = 45
cutmode = 0



if forward:
    treeSP, n_pulses_SP = getattr(mt2.NCorrRun("SP",target,Forward = forward),"{}_doubles_tree".format(particle))
    treeDP,n_pulses_DP = getattr(mt2.NCorrRun("DP",target, Forward = forward),"{}_doubles_tree".format(particle))
else:
    treeSP, n_pulses_SP = getattr(mt2.NCorrRun("SP", target), "{}_doubles_tree".format(particle))
    treeDP, n_pulses_DP = getattr(mt2.NCorrRun("DP", target), "{}_doubles_tree".format(particle))


delta_thetapar= 90- delta_thetaperp
cutPerp = mt.cut_rangeAND([(3.1415/180)*(90-delta_thetaperp), (3.1415/180)*(90+delta_thetaperp)],"neutrons.coinc_hits.theta_abs[0][1]")
cutPar = mt.cut_range_outer([(3.1415/180)*(delta_thetapar), (3.1415/180)*(180-delta_thetapar)],"neutrons.coinc_hits.theta_abs[0][1]")


histSP = TH2F(mins=0,maxs=6,binwidths=1.5)
histDP = TH2F(mins=0,maxs=6,binwidths=1.5)


if forward:
    # cut = "{0}.coinc_hits[].ForwardDet[0]==0 && {0}.coinc_hits[].ForwardDet[1]==0 ".format(particle)
    cut = ["",cutPar,cutPerp][cutmode]
    drw_str = "neutrons.coinc_hits[0].erg[0]:neutrons.coinc_hits[0].erg[1]"
else:
    drw_str = "neutrons.hits[0].erg:neutrons.hits[1].erg"
    cut = ""


histSP.Project(treeSP,drw_str, cut, weight=1.0/n_pulses_SP)
histDP.Project(treeDP,drw_str, cut, weight=1.0/n_pulses_DP)

c1 = ROOT.TCanvas()

c1.Divide(2)
c1.cd(1)




histSP -= 0.5*histDP

histSP.normalize()
histDP.normalize()

histSP /= histDP

histSP.GetXaxis().SetTitle("E1 [MeV]")
histSP.GetYaxis().SetTitle("E2 [MeV]")
histSP.GetZaxis().SetTitle("normalized correlation [arb. units]")

histSP.GetYaxis().SetTitleOffset(1.75)
histSP.GetXaxis().SetTitleOffset(1.75)
histSP.GetZaxis().SetTitleOffset(1.4)
# print histSP.GetTitleOffset(0)
# print histSP.GetTitleOffset(1)

histSPTr = histSP.__copy__()
histSPTr.transpose()
histSP += histSPTr
histSP *=0.5

histSP.SetName("E1vsE2")
histSP.SetTitle(["No #theta_{{abs}} cut", "parallel to beam  (#parallel #pm {par}#circ) ", "perpendicular to beam (#perp #pm {perp})"][cutmode].format(par=delta_thetapar, perp = delta_thetaperp))
histSP.SetStats(0)

histSP.Draw("Lego2 E",False)

c1.cd(2)
ROOT.TGaxis.SetMaxDigits(2)

_ = histSP.__copy__()
_.SetTitle("same data")
_.SetStats(0)
# _.Draw("colz",False)

errors = histSP.getErrorsHistogram(relative_error=False if target=="D2O" else False)
errors.SetStats(0)
errors.GetZaxis().SetTitleOffset(1.65)

t = "relative error" if target!="D2O" else "Absolute error"

errors.GetZaxis().SetTitle(t)

errors.SetTitle(t)
errors.Draw("Lego",False)




if __name__ == "__main__":
    import ROOT as ROOT
    from multiprocessing import Process, Queue
    import time, sys, os


    def input_thread(q, stdin):
        while True:
            print 'ROOT: '
            cmd = stdin.readline()
            q.put(cmd)


    def root(char):
        assert isinstance(char, str), "Argument must be string!"
        ROOT.gROOT.ProcessLine(char)


    if __name__ == '__main__':
        ___queue___ = Queue()
        ___newstdin___ = os.fdopen(os.dup(sys.stdin.fileno()))
        ___input_p___ = Process(target=input_thread,
                                args=(___queue___, ___newstdin___))
        ___input_p___.daemon = True
        ___input_p___.start()
        ___g___ = ROOT.gSystem.ProcessEvents
        try:
            while 1:
                if not ___queue___.empty():
                    ___cmd___ = ___queue___.get()
                    try:
                        exec (___cmd___, globals())
                    except:
                        print sys.exc_info()
                time.sleep(0.01)
                ___g___()
        except(KeyboardInterrupt, SystemExit):
            ___input_p___.terminate()