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()