Project

General

Profile

Feature #1152 » Plot.py

Carlisle, Timothy, 08 October 2012 14:24

 
1
from xboa  import *
2
from xboa.Hit   import Hit
3
from xboa.Bunch import Bunch
4
import sys
5
import os
6
import subprocess
7
from array import array
8
import math
9
import operator
10

    
11
import ROOT 
12

    
13
#some input data
14
filename = sys.argv[1]
15
filetype = "maus_root_virtual_hit"
16

    
17
print "First loading the data... "
18
bunch_list = Bunch.new_list_from_read_builtin(filetype, filename)
19

    
20
# now load the ROOT file
21

    
22
root_file = ROOT.TFile(filename, "READ") # pylint: disable = E1101
23

    
24
# and set up the data tree to be filled by ROOT IO
25

    
26
data = ROOT.MAUS.Data() # pylint: disable = E1101
27
tree = root_file.Get("Spill")
28
tree.SetBranchAddress("data", data)
29

    
30
vhit_z = ROOT.TH1D("z[m]",  #pylint: disable = E1101
31
                                  " z [m]; ",
32
                                  50, -6, 6)
33
vhit_Pz = ROOT.TH1D("Pz[MeV/c]",  #pylint: disable = E1101
34
                                   " Pz [MeV/c]; ",
35
                                   50, 170, 210)
36

    
37
bunch_list = Bunch.new_list_from_read_builtin(filetype, filename)
38

    
39
print " no. Virtual Planes: " + str( len(bunch_list) )
40

    
41
emittance4, beta4, Pt = array( 'd' ), array( 'd' ), array('d')
42
Bz, x, y, z, Pz, Px, Energy, Amplitude = array( 'd' ), array( 'd' ), array( 'd' ), array( 'd' ), array( 'd' ), array( 'd' ), array( 'd'), array( 'd')
43

    
44
Nmu = array( 'd' )
45

    
46
bunch_list[100].cut({'y':-300.}, operator.le, global_cut=True)
47
bunch_list[100].cut({'y':300.}, operator.ge, global_cut=True)
48
bunch_list[100].root_histogram('y','mm')
49
raw_input('Press Enter to exit')
50

    
51
bunch_list[100].clear_weights()
52

    
53
bunch_list[100].cut({'y':-300.}, operator.le, global_cut=False)
54
bunch_list[100].cut({'y':300.}, operator.ge, global_cut=False)
55
bunch_list[100].root_histogram('y','mm')
56
raw_input('Press Enter to exit')
57

    
58
for Plane in range( len(bunch_list) ):
59
    #    bunch_list[Plane].cut({'px':0.})#, operator.le, global_cut=True)
60
    #print str(Plane) + " loop 2"
61
    my_bunch = bunch_list[Plane]
62
    emittance4.append( my_bunch.get_emittance(['x','y'])  )
63
    beta4.append( my_bunch.get_beta(['x','y']) /10 )
64

    
65
    x.append(  my_bunch.get('mean', ['x']  ))
66
    y.append(  my_bunch.get('mean', ['y']  ))
67
    z.append(  my_bunch.get('mean', ['z']  ))
68

    
69
    Pz.append( my_bunch.get('mean',['pz']  ))
70
    Px.append( my_bunch.get('mean',['px']  ))
71
    Energy.append(  my_bunch.get('mean', ['energy']  ) ) 
72
    Pt.append( (  (my_bunch.get( 'mean',['px'] )  )**2 +  (my_bunch.get( 'mean',['py']) )**2   ) ** 0.5 )
73
    Bz.append(my_bunch.get('mean',['bz']) * 1000 )
74
    #Amplitude.append(my_bunch.get('mean',['amplitude x y'] ))
75
    Nmu.append( my_bunch.bunch_weight() )
76

    
77

    
78
cX = ROOT.TCanvas("cX", # pylint: disable = E1101
79
                        "cX")
80
cX.Divide(2)
81
cX.cd(1)                      
82
g27 =  ROOT.TGraph(len(z), z, Nmu)
83
g27.Draw("AL")
84

    
85
cX.cd(2)  
86
g28 =  ROOT.TGraph(len(z), z, x )
87
g28.Draw("AL")
88

    
89

    
90
c2 = ROOT.TCanvas("c2", # pylint: disable = E1101
91
                        "c2")
92
c2.Divide(2)
93
c2.cd(1)
94
g21 =  ROOT.TGraph(len(z), z, emittance4)
95
g21.SetTitle("#epsilon 4d [mm]")
96
g21.GetXaxis().SetTitle("z [m]")
97
g21.GetYaxis().SetTitle("#e [mm]")
98
g21.Draw("AL")
99

    
100
c2.cd(2)
101
g22 = ROOT.TGraph(len(z), z, Energy)
102
g22.SetTitle("Energy [MeV]")
103
g22.GetXaxis().SetTitle("z [m]")
104
g22.GetYaxis().SetTitle("Energy [MeV]")
105
g22.Draw("AL")
106

    
107
c3 = ROOT.TCanvas("c3", # pylint: disable = E1101
108
                        "c3")
109
c3.Divide(2)
110
c3.cd(1)
111
g31 = ROOT.TGraph(len(z), z, Pt)
112
g31.SetTitle("Pt [MeV/c]")
113
g31.GetXaxis().SetTitle("z [m]")
114
g31.GetYaxis().SetTitle("")
115
g31.Draw("AL")
116

    
117
c3.cd(2)
118
g32 = ROOT.TGraph(len(z), z, Bz)
119
g32.SetTitle("Bz [T]")
120
g32.GetXaxis().SetTitle("z [m]")
121
g32.GetYaxis().SetTitle("Bz [T]")
122
g32.Draw("AL")
123

    
124
c1 = ROOT.TCanvas("c1", # pylint: disable = E1101
125
                        "c1")
126
c1.Divide(2)
127
c1.cd(1)
128

    
129
g = ROOT.TGraph(len(z), z, Pz)
130
g.SetTitle("Pz [MeV/c]")
131
g.GetXaxis().SetTitle("z [m]")
132
#g.GetYaxis().SetTitle("Pz [MeV/c]")
133
g.Draw("AL")
134

    
135
c1.cd(2)
136
g1 = ROOT.TGraph(len(z), z, beta4)
137
g1.SetTitle("#beta [cm]")
138
g1.GetXaxis().SetTitle("z [m]")
139
#g1.GetYaxis().SetTitle("#beta [cm]")
140
g1.Draw("AL")
141
print " change in emittance: " + str( abs( emittance4[0] - emittance4[ len(emittance4) - 1 ]) )
142

    
143
raw_input('Press Enter to exit')
144

    
145
    
146
raise SystemExit
147

    
148

    
149

    
150

    
151
                  
152

    
153

    
154
            
155

    
156

    
157

    
158

    
159
                    
(1-1/4)