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)

1. print histSP.GetTitleOffset(0)
2. print histSP.GetTitleOffset(1)

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

histSP.SetName("E1vsE2") histSP.SetTitle(["No #theta_Template: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)

1. _.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()