#!/usr/bin/env python """ R. Bayes An example of extracting information from the special virtual planes """ # pylint: disable = E1101, C0103, C0301, W0611, R0914, R0915 import os import sys import ROOT import math import subprocess import libMausCpp def generate_simulation(outfile): """ Run the simulation to make a data file for consideration """ analysis = os.path.join\ (os.environ["MAUS_ROOT_DIR"],"bin","simulate_mice.py") proc = subprocess.Popen(['python', analysis, \ '--output_root_file_name', outfile, '--simulation_geometry_file', 'Stage4.dat']) proc.wait() def main(args): """ extract information from the special virtual planes from the maus_output.root file in the current directory """ Bfield = [] zpos = [] h = ROOT.TH1D("hrate", \ ";Position along Z axis (m);Particles at Virtual Plane", \ 101,13.95,24.05) hf = ROOT.TH2D("hfield", \ ";Position along Z axis (m); Field Magnitude", \ 101,13.95,24.05,100,-5,5) hp = ROOT.TH2D("hp", \ ";Position along Z axis (m); Momentum (MeV/c)", \ 101,13.95,24.05,100,0,400) hpz = ROOT.TH2D("hpz", \ ";Position along Z axis (m); Longitudinal Momentum (MeV/c)", \ 101,13.95,24.05,100,0,400) hp200 = ROOT.TH2D("hp200", \ ";Position along Z axis (m); Momentum (MeV/c)", \ 101,13.95,24.05,100,0,400) hr = ROOT.TH2D("hr", \ ";Position along Z axis (m); Radial Position (m)", \ 101,13.95,24.05,100,0.0, 500) hpx = ROOT.TH2D("hpx", \ ";Position along Z axis (m); Momentum along X axis (MeV/c)", \ 101,13.95,24.05,100,-50,50) hpy = ROOT.TH2D("hpy", \ ";Position along Z axis (m); Momentum along Y axis (MeV/c)", \ 101,13.95,24.05,100,-50,50) hpx2 = ROOT.TH2D("hpx2", \ ";Position along Z axis (m); Momentum^{2} along X axis (MeV/c)", \ 101,13.95,24.05,100,-50,50) hpy2 = ROOT.TH2D("hpy2", \ ";Position along Z axis (m); Momentum^{2} along Y axis (MeV/c)", \ 101,13.95,24.05,100,-50,50) hx = ROOT.TH2D("hx", \ ";Position along Z axis (m); Position along X axis (m)", \ 101,13.95,24.05,100,-500, 500) hy = ROOT.TH2D("hy", \ ";Position along Z axis (m); Position along Y axis (m)", \ 101,13.95,24.05,100,-500, 500) hx2 = ROOT.TH2D("hx2", \ ";Position along Z axis (m); Position^{2} along X axis (m^{2})", \ 101,13.95,24.05,100,-500, 500) hy2 = ROOT.TH2D("hy2", \ ";Position along Z axis (m); Position^{2} along Y axis (m^{2})", \ 101,13.95,24.05,100,-500, 500) hypx = ROOT.TH2D("hypx", \ ";Position along Z axis (m); y - p_{x} correlation (m MeV)", \ 101,13.95,24.05,100,-50,50) hypy = ROOT.TH2D("hypy", \ ";Position along Z axis (m); y - p_{y} correlation (m MeV)", \ 101,13.95,24.05,100,-50,50) hxpx = ROOT.TH2D("hxpx", \ ";Position along Z axis (m); x - p_{x} correlation (m MeV)", \ 101,13.95,24.05,100,-50,50) hxpy = ROOT.TH2D("hxpy", \ ";Position along Z axis (m); x - p_{y} correlation (m MeV)", \ 101,13.95,24.05,100,-50,50) hxy = ROOT.TH2D("hxy", \ ";Position along Z axis (m); x - y correlation (m^{2})", \ 101,13.95,24.05,100,-500,500) hpxpy = ROOT.TH2D("hpxpy", \ ";Position along Z axis (m); p_{x} - p_{y} correlation (MeV^{2})", \ 101,13.95,24.05,100,-50,50) hbt = ROOT.TH2D("hbt", \ ";Position along Z axis (m); Transverse #beta", \ 101,13.95,24.05,100,0, 0.5) maxsize = -999 for fp in args: rfile = ROOT.TFile(str(fp),"R") if rfile.IsZombie(): continue data = ROOT.MAUS.Data() tree = rfile.Get("Spill") if not tree: continue tree.SetBranchAddress("data", data) print "Evaluating "+str(fp) for i in range(tree.GetEntries()): tree.GetEntry(i) spill = data.GetSpill() if spill.GetDaqEventType() == "physics_event": for j in range(spill.GetMCEvents().size()): # if this is not the largest set of virtual hit go to the # next particle if spill.GetMCEvents()[j].GetPrimary().GetParticleId() != -13: continue passTOF2 = False nvirthits = spill.GetMCEvents()[j].GetVirtualHits().size() if nvirthits >= 84: xpTOF2 = spill.GetMCEvents()[j].GetVirtualHits()[83]\ .GetPosition().x() ypTOF2 = spill.GetMCEvents()[j].GetVirtualHits()[83]\ .GetPosition().y() if math.fabs(xpTOF2) < 600 \ and math.fabs(ypTOF2) < 600: passTOF2 = True for k in range(nvirthits): h.Fill(spill.GetMCEvents()[j] \ .GetVirtualHits()[k].GetPosition().z() / 1000) r = math.sqrt((spill.GetMCEvents()[j] \ .GetVirtualHits()[k].GetPosition().x())**2 + (spill.GetMCEvents()[j] \ .GetVirtualHits()[k].GetPosition().y())**2 ) if r < 200.: hf.Fill( \ spill.GetMCEvents()[j] \ .GetVirtualHits()[k].GetPosition().z() / 1000, spill.GetMCEvents()[j] \ .GetVirtualHits()[k].GetBField().z() * 1000) hp200.Fill( \ spill.GetMCEvents()[j] \ .GetVirtualHits()[k].GetPosition().z() / 1000, spill.GetMCEvents()[j] \ .GetVirtualHits()[k].GetMomentum().mag()) if not passTOF2: continue z = spill.GetMCEvents()[j].GetVirtualHits()[k].GetPosition().z() / 1000. x = spill.GetMCEvents()[j].GetVirtualHits()[k].GetPosition().x() y = spill.GetMCEvents()[j].GetVirtualHits()[k].GetPosition().y() pz = spill.GetMCEvents()[j].GetVirtualHits()[k].GetMomentum().z() px = spill.GetMCEvents()[j].GetVirtualHits()[k].GetMomentum().x() py = spill.GetMCEvents()[j].GetVirtualHits()[k].GetMomentum().y() hp.Fill(z, spill.GetMCEvents()[j] \ .GetVirtualHits()[k].GetMomentum().mag()) hpx.Fill(z, px) hpy.Fill(z, py) hpx2.Fill(z, px*px) hpy2.Fill(z, py*py) hpz.Fill(z, pz) hr.Fill(z, math.sqrt(x*x + y*y)) hx.Fill(z, x) hy.Fill(z, y) hx2.Fill(z, x*x) hy2.Fill(z, y*y) hypx.Fill(z, y*px) hypy.Fill(z, y*py) hxpx.Fill(z, x*px) hxpy.Fill(z, x*py) hxy.Fill(z, x*y) hpxpy.Fill(z, px*py) hbt.Fill(spill.GetMCEvents()[j] \ .GetVirtualHits()[k].GetPosition().z() / 1000, math.sqrt((spill.GetMCEvents()[j] \ .GetVirtualHits()[k].GetMomentum().x())**2 + (spill.GetMCEvents()[j] \ .GetVirtualHits()[k].GetMomentum().y())**2 )/ math.sqrt((spill.GetMCEvents()[j] \ .GetVirtualHits()[k].GetMomentum().mag())**2 \ + 105.65**2) ) if maxsize < nvirthits: maxsize = spill.GetMCEvents()[j].GetVirtualHits().size() print maxsize # reinitialize the Bfiled and zpos objects Bfield = [] zpos = [] for k in range(maxsize): Bfield.append(spill.GetMCEvents()[j] \ .GetVirtualHits()[k].GetBField().z() * 1000) zpos.append(spill.GetMCEvents()[j] \ .GetVirtualHits()[k].GetPosition().z() / 1000) mz = [] mx = [] my = [] mpx = [] mpy = [] mpz = [] mxpx = [] mypx = [] mxpy = [] mypy = [] mxy = [] mpxpy = [] Ex = [] Ey = [] Et = [] bx = [] by = [] bt = [] ax = [] ay = [] at = [] mmu = 105.65 gEx = ROOT.TGraph(101) gEx.SetTitle(";Z Position (m);2D emittance (mm)") gEy = ROOT.TGraph(101) gEy.SetTitle(";Z Position (m);2D emittance (mm)") gEt = ROOT.TGraph(101) gEt.SetTitle(";Z Position (m);4D emittance (mm)") gbx = ROOT.TGraph(101) gbx.SetTitle(";Z Position (m);beta (mm)") gby = ROOT.TGraph(101) gby.SetTitle(";Z Position (m);beta (mm)") gbt = ROOT.TGraph(101) gbt.SetTitle(";Z Position (m);beta (mm)") gax = ROOT.TGraph(101) gax.SetTitle(";Z Position (m);alpha (mm)") gay = ROOT.TGraph(101) gay.SetTitle(";Z Position (m);alpha (mm)") gat = ROOT.TGraph(101) gat.SetTitle(";Z Position (m);alpha (mm)") for i in range(0,101): mz.append(hx.GetXaxis().GetBinCenter(i)) mx.append(math.fabs(hx2.ProfileX().GetBinContent(i) - \ hx.ProfileX().GetBinContent(i)*hx.ProfileX().GetBinContent(i))) my.append(math.fabs(hy2.ProfileX().GetBinContent(i) - \ hy.ProfileX().GetBinContent(i)*hy.ProfileX().GetBinContent(i))) mpx.append(math.fabs(hpx2.ProfileX().GetBinContent(i) - \ hpx.ProfileX().GetBinContent(i)*hpx.ProfileX().GetBinContent(i))) mpy.append(math.fabs(hpy2.ProfileX().GetBinContent(i) - \ hpy.ProfileX().GetBinContent(i)*hpy.ProfileX().GetBinContent(i))) mxpx.append((hxpx.ProfileX().GetBinContent(i) - \ hpx.ProfileX().GetBinContent(i)*hx.ProfileX().GetBinContent(i))) mypx.append((hypx.ProfileX().GetBinContent(i) - \ hpx.ProfileX().GetBinContent(i)*hy.ProfileX().GetBinContent(i))) mxpy.append((hxpy.ProfileX().GetBinContent(i) - \ hpy.ProfileX().GetBinContent(i)*hx.ProfileX().GetBinContent(i))) mypy.append((hypy.ProfileX().GetBinContent(i) - \ hpy.ProfileX().GetBinContent(i)*hy.ProfileX().GetBinContent(i))) mxy.append((hxy.ProfileX().GetBinContent(i) - \ hy.ProfileX().GetBinContent(i)*hx.ProfileX().GetBinContent(i))) mpxpy.append((hpxpy.ProfileX().GetBinContent(i) - \ hpx.ProfileX().GetBinContent(i)*hpy.ProfileX().GetBinContent(i))) mpz.append(hpz.ProfileX().GetBinContent(i)) Ex.append( math.sqrt(math.fabs( mx[-1]*mpx[-1] - \ mxpx[-1]*mxpx[-1] )) / mmu) Ey.append( math.sqrt(math.fabs( my[-1]*mpy[-1] - \ mypy[-1]*mypy[-1] )) / mmu) if mx[-1] > 0: # Et.append( math.sqrt( math.fabs( mx[-1]*mpx[-1] - mxpx[-1]*mxpx[-1]\ # -mxpy[-1]*mxpy[-1])) / mmu ) Et.append( math.pow( math.fabs(mx[-1] * ( mpx[-1] *(my[-1]*mpy[-1] - mypy[-1]*mypy[-1]) + mypx[-1] *(mypy[-1]*mpxpy[-1] - mpy[-1]*mpxpy[-1]) + mpxpy[-1]*(mypx[-1]*mypy[-1] - my[-1]*mpxpy[-1]) ) - mxpx[-1]*( mxpx[-1] *(my[-1]*mpy[-1] - mypy[-1]*mypy[-1]) + mypx[-1] *(mypy[-1]*mxpy[-1] - mpy[-1]*mxy[-1]) + mpxpy[-1]*(mxy[-1]*mypy[-1] - my[-1]*mxpy[-1]) ) + mxy[-1] *( mxpx[-1] *(mypx[-1]*mpy[-1] - mypy[-1]*mypy[-1]) + mpx[-1] *(mypy[-1]*mxpy[-1] - mpy[-1]*mxy[-1]) + mpxpy[-1]*(mxy[-1]*mpxpy[-1] - mypx[-1]*mxpy[-1]) ) - mxpy[-1]*( mxpx[-1] *(mypx[-1]*mypy[-1] - my[-1]*mpxpy[-1]) + mpx[-1] *(my[-1]*mxpy[-1] - mypy[-1]*mxy[-1]) + mypx[-1] *(mxy[-1]*mpxpy[-1] - mypx[-1]*mxpy[-1])) ), 0.25 ) / mmu ) else: Et.append( 0.0 ) if Ex[-1] > 0: bx.append( mx[-1]*mpz[-1] / mmu / Ex[-1] ) ax.append( -mxpx[-1]/ mmu / Ex[-1] ) else: bx.append( 0.0 ) ax.append( 0.0 ) if Ey[-1] > 0: by.append( my[-1]*mpz[-1] / mmu / Ey[-1] ) ay.append( -mypy[-1]/ mmu / Ey[-1] ) else: by.append( 0.0 ) ay.append( 0.0 ) if Et[-1] > 0: bt.append( (mx[-1] + my[-1])*mpz[-1] / mmu / Et[-1] ) at.append( -(mxpx[-1] + mypy[-1])/ mmu / Et[-1] ) else: bt.append( 0.0 ) at.append( 0.0 ) gEx.SetPoint(i, mz[-1], Ex[-1]) gEy.SetPoint(i, mz[-1], Ey[-1]) gEt.SetPoint(i, mz[-1], Et[-1]) gbx.SetPoint(i, mz[-1], bx[-1]) gby.SetPoint(i, mz[-1], by[-1]) gbt.SetPoint(i, mz[-1], bt[-1]) gax.SetPoint(i, mz[-1], ax[-1]) gay.SetPoint(i, mz[-1], ay[-1]) gat.SetPoint(i, mz[-1], at[-1]) xmin = 13400 xmax = 20000 g = ROOT.TGraph(maxsize) for i in range(maxsize): g.SetPoint(i, zpos[i], Bfield[i]) g.SetMarkerColor(2) g.SetMarkerStyle(21) g.SetTitle("; Position on Z axis (m); Field Magnitude (Tesla)") c = ROOT.TCanvas() g.GetHistogram().GetXaxis().SetRangeUser(xmin,xmax) g.Draw('al') c.Print("Bfield_v_z.eps") c0_1 = ROOT.TCanvas() gEx.Draw("al") gEx.GetHistogram().GetXaxis().SetRangeUser(xmin,xmax) c0_1.Print("emittance_X.eps") c0_2 = ROOT.TCanvas() gEy.Draw("al") gEy.GetHistogram().GetXaxis().SetRangeUser(xmin,xmax) c0_2.Print("emittance_Y.eps") c0_3 = ROOT.TCanvas() gEt.Draw("al") gEt.GetHistogram().GetXaxis().SetRangeUser(xmin,xmax) c0_3.Print("emittance_4D.eps") c0_4 = ROOT.TCanvas() gbx.Draw("al") gbx.GetHistogram().GetXaxis().SetRangeUser(xmin,xmax) c0_4.Print("beta_X.eps") c0_5 = ROOT.TCanvas() gby.Draw("al") gby.GetHistogram().GetXaxis().SetRangeUser(xmin,xmax) c0_5.Print("beta_Y.eps") c0_6 = ROOT.TCanvas() gbt.Draw("al") gbt.GetHistogram().GetXaxis().SetRangeUser(xmin,xmax) c0_6.Print("beta_4D.eps") c0_7 = ROOT.TCanvas() gax.Draw("al") gax.GetHistogram().GetXaxis().SetRangeUser(xmin,xmax) c0_7.Print("alpha_X.eps") c0_8 = ROOT.TCanvas() gay.Draw("al") gay.GetHistogram().GetXaxis().SetRangeUser(xmin,xmax) c0_8.Print("alpha_Y.eps") c0_9 = ROOT.TCanvas() gat.Draw("al") gat.GetHistogram().GetXaxis().SetRangeUser(xmin,xmax) c0_9.Print("alpha_4D.eps") c1 = ROOT.TCanvas() hf.GetXaxis().SetRangeUser(xmin,xmax) hf.Draw("colz") c1.Print("Bfield_v_z_hist.eps") c3 = ROOT.TCanvas() c3.SetLogz() hp.GetXaxis().SetRangeUser(xmin,xmax) hp.Draw('colz') c3.Print("Momentum_v_z_hist.eps") c3_1 = ROOT.TCanvas() hp_px = hp.ProjectionX() hp_px.GetXaxis().SetRangeUser(xmin,xmax) hp_px.Draw() c3_1.Print("Momentum_v_z_projx.eps") c3_2 = ROOT.TCanvas() hp_pfx = hp.ProfileX() hp_pfx.GetXaxis().SetRangeUser(xmin,xmax) hp_pfx.Draw() c3_2.Print("Momentum_v_z_prflx.eps") c3_3 = ROOT.TCanvas() c3_3.SetLogz() hp200.GetXaxis().SetRangeUser(xmin,xmax) hp200.Draw('colz') c3_3.Print("Momentum_v_z_hist_lt200mm.eps") c3_4 = ROOT.TCanvas() hp200_px = hp200.ProjectionX() hp200_px.GetXaxis().SetRangeUser(xmin,xmax) hp200_px.Draw() c3_4.Print("Momentum_v_z_projx_lt200mm.eps") c3_5 = ROOT.TCanvas() hp200_pfx = hp200.ProfileX() hp200_pfx.GetXaxis().SetRangeUser(xmin,xmax) hp200_pfx.Draw() c3_5.Print("Momentum_v_z_prflx_lt200mm.eps") c4 = ROOT.TCanvas() c4.SetLogz() hr.GetXaxis().SetRangeUser(xmin,xmax) hr.Draw('colz') c4.Print("BeamRadius_v_z_hist.eps") c5 = ROOT.TCanvas() c5.SetLogz() hbt.GetXaxis().SetRangeUser(xmin,xmax) hbt.Draw('colz') c5.Print("TransverseBeta_v_z_hist.eps") c6 = ROOT.TCanvas() c6.SetLogz() h.GetXaxis().SetRangeUser(xmin,xmax) h.Draw() c6.Print("MuonRate_v_z_hist.eps") c7_0 = ROOT.TCanvas() c7_0.SetLogz() hpx.GetXaxis().SetRangeUser(xmin,xmax) hpx.Draw('colz') c7_0.Print("XMomentum_v_z_hist.eps") c7_0a = ROOT.TCanvas() hpx_pfx = hpx.ProfileX() hpx_pfx.GetXaxis().SetRangeUser(xmin,xmax) hpx_pfx.Draw('colz') c7_0a.Print("XMomentum_v_z_prflx.eps") c7_1 = ROOT.TCanvas() c7_1.SetLogz() hpy.GetXaxis().SetRangeUser(xmin,xmax) hpy.Draw('colz') c7_1.Print("YMomentum_v_z_hist.eps") c7_1a = ROOT.TCanvas() hpy_pfx = hpy.ProfileX() hpy_pfx.GetXaxis().SetRangeUser(xmin,xmax) hpy_pfx.Draw('colz') c7_1a.Print("YMomentum_v_z_prflx.eps") c7_2 = ROOT.TCanvas() c7_2.SetLogz() hpz.GetXaxis().SetRangeUser(xmin,xmax) hpz.Draw('colz') c7_2.Print("ZMomentum_v_z_hist.eps") c7_2a = ROOT.TCanvas() hpz_pfx = hpz.ProfileX() hpz_pfx.GetXaxis().SetRangeUser(xmin,xmax) hpz_pfx.Draw('colz') c7_2a.Print("ZMomentum_v_z_prflx.eps") c8_0 = ROOT.TCanvas() c8_0.SetLogz() hx.GetXaxis().SetRangeUser(xmin,xmax) hx.Draw('colz') c8_0.Print("XPosition_v_z_hist.eps") c8_0a = ROOT.TCanvas() hx_pfx = hx.ProfileX() hx_pfx.GetXaxis().SetRangeUser(xmin,xmax) hx_pfx.Draw('colz') c8_0a.Print("XPosition_v_z_prflx.eps") c8_1 = ROOT.TCanvas() c8_1.SetLogz() hy.GetXaxis().SetRangeUser(xmin,xmax) hy.Draw('colz') c8_1.Print("YPosition_v_z_hist.eps") c8_1a = ROOT.TCanvas() hy_pfx = hy.ProfileX() hy_pfx.GetXaxis().SetRangeUser(xmin,xmax) hy_pfx.Draw('colz') c8_1a.Print("YPosition_v_z_prflx.eps") if __name__ == "__main__": if len(sys.argv) > 1: main(sys.argv[1:]) else: ofile = os.path.join\ (os.environ['MAUS_ROOT_DIR'], 'tmp','example_simulation_file.root') generate_simulation(ofile) main([ofile])