← Back to branch summary

~vpec/maus/tof_calib_read

« back to all changes in this revision

Viewing changes to bin/examples/plot_virtuals.py

Preparing for release 0.9.3

Show diffs side-by-side

added added

removed removed

Lines of Context:
bin/examples/plot_virtuals.py
 
1
#!/usr/bin/env python
 
2
 
 
3
 
 
4
"""
 
5
R. Bayes
 
6
 
 
7
An example of extracting information from the special virtual planes
 
8
"""
 
9
 
 
10
# pylint: disable = E1101, C0103, C0301, W0611, R0914, R0915
 
11
import os
 
12
import sys
 
13
import ROOT
 
14
import math
 
15
import subprocess
 
16
import libMausCpp
 
17
 
 
18
def generate_simulation(outfile):
 
19
    """
 
20
    Run the simulation to make a data file for consideration
 
21
    """
 
22
    analysis = os.path.join\
 
23
               (os.environ["MAUS_ROOT_DIR"],"bin","simulate_mice.py")
 
24
    proc = subprocess.Popen(['python', analysis, \
 
25
                             '--output_root_file_name', outfile,
 
26
                             '--simulation_geometry_file', 'Stage4.dat'])
 
27
    proc.wait()
 
28
 
 
29
def main(args):
 
30
 
 
31
    """ extract information from the special virtual planes
 
32
    from the maus_output.root file in the current directory """
 
33
    
 
34
    Bfield = []
 
35
    zpos   = []
 
36
    h = ROOT.TH1D("hrate", \
 
37
                  ";Position along Z axis (m);Particles at Virtual Plane", \
 
38
                  240,0,24)
 
39
    hf = ROOT.TH2D("hfield", \
 
40
                   ";Position along Z axis (m); Field Magnitude", \
 
41
                   240,0,24,100,-5,5)
 
42
    hp = ROOT.TH2D("hp", \
 
43
                   ";Position along Z axis (m); Momentum (MeV/c)", \
 
44
                   240,0,24,100,0,400)
 
45
    hp200 = ROOT.TH2D("hp200", \
 
46
                      ";Position along Z axis (m); Momentum (MeV/c)", \
 
47
                      240,0,24,100,0,400)
 
48
    hr = ROOT.TH2D("hr", \
 
49
                   ";Position along Z axis (m); Radial Position (m)", \
 
50
                   240,0,24,100,0.0, 0.5)
 
51
    hbt = ROOT.TH2D("hbt", \
 
52
                    ";Position along Z axis (m); Transverse #beta", \
 
53
                    240,0,24,100,0, 0.5)
 
54
    maxsize = -999
 
55
 
 
56
    for fp in args:
 
57
        
 
58
        rfile = ROOT.TFile(str(fp),"R")
 
59
        if rfile.IsZombie():
 
60
            continue
 
61
        data = ROOT.MAUS.Data()
 
62
        tree = rfile.Get("Spill")
 
63
        if not tree:
 
64
            continue
 
65
        tree.SetBranchAddress("data", data) 
 
66
        print "Evaluating "+str(fp)
 
67
        for i in range(tree.GetEntries()):
 
68
            tree.GetEntry(i)
 
69
            spill = data.GetSpill()
 
70
            
 
71
            if spill.GetDaqEventType() == "physics_event":
 
72
                for j in range(spill.GetMCEvents().size()):
 
73
                    # if this is not the largest set of virtual hit go to the
 
74
                    # next particle
 
75
                    if spill.GetMCEvents()[j].GetPrimary().GetParticleId() != -13:
 
76
                        continue
 
77
                    nvirthits = spill.GetMCEvents()[j].GetVirtualHits().size()
 
78
                    for k in range(nvirthits):
 
79
                        h.Fill(spill.GetMCEvents()[j] \
 
80
                               .GetVirtualHits()[k].GetPosition().z() / 1000)
 
81
                        r = math.sqrt((spill.GetMCEvents()[j] \
 
82
                                       .GetVirtualHits()[k].GetPosition().x())**2 +
 
83
                                      (spill.GetMCEvents()[j] \
 
84
                                       .GetVirtualHits()[k].GetPosition().y())**2 )
 
85
                        if r < 200.:
 
86
                            hf.Fill( \
 
87
                                spill.GetMCEvents()[j] \
 
88
                                .GetVirtualHits()[k].GetPosition().z() / 1000,
 
89
                                spill.GetMCEvents()[j] \
 
90
                                .GetVirtualHits()[k].GetBField().z() * 1000)
 
91
                            hp200.Fill( \
 
92
                                spill.GetMCEvents()[j] \
 
93
                                .GetVirtualHits()[k].GetPosition().z() / 1000,
 
94
                                spill.GetMCEvents()[j] \
 
95
                                .GetVirtualHits()[k].GetMomentum().mag())
 
96
                        hp.Fill(spill.GetMCEvents()[j] \
 
97
                                .GetVirtualHits()[k].GetPosition().z() / 1000,
 
98
                                spill.GetMCEvents()[j] \
 
99
                                .GetVirtualHits()[k].GetMomentum().mag())
 
100
                        hr.Fill(spill.GetMCEvents()[j]\
 
101
                                .GetVirtualHits()[k].GetPosition().z() / 1000,
 
102
                                math.sqrt((spill.GetMCEvents()[j] \
 
103
                                           .GetVirtualHits()[k].GetPosition().x())**2 +
 
104
                                          (spill.GetMCEvents()[j] \
 
105
                                           .GetVirtualHits()[k].GetPosition().y())**2 )/ 1000.)
 
106
                        hbt.Fill(spill.GetMCEvents()[j] \
 
107
                                 .GetVirtualHits()[k].GetPosition().z() / 1000,
 
108
                                 math.sqrt((spill.GetMCEvents()[j] \
 
109
                                            .GetVirtualHits()[k].GetMomentum().x())**2 +
 
110
                                           (spill.GetMCEvents()[j] \
 
111
                                            .GetVirtualHits()[k].GetMomentum().y())**2 )/
 
112
                                 math.sqrt((spill.GetMCEvents()[j] \
 
113
                                            .GetVirtualHits()[k].GetMomentum().mag())**2 \
 
114
                                           + 105.65**2) )
 
115
                    if maxsize < nvirthits:
 
116
                        maxsize = spill.GetMCEvents()[j].GetVirtualHits().size()
 
117
                        print maxsize
 
118
                        # reinitialize the Bfiled and zpos objects
 
119
                        Bfield = []
 
120
                        zpos   = []
 
121
                        for k in range(maxsize):
 
122
                            Bfield.append(spill.GetMCEvents()[j] \
 
123
                                          .GetVirtualHits()[k].GetBField().z() * 1000)
 
124
                            zpos.append(spill.GetMCEvents()[j] \
 
125
                                        .GetVirtualHits()[k].GetPosition().z() / 1000)
 
126
 
 
127
    g = ROOT.TGraph(maxsize)
 
128
    for i in range(maxsize):
 
129
        g.SetPoint(i, zpos[i], Bfield[i])
 
130
    g.SetMarkerColor(2)
 
131
    g.SetMarkerStyle(21)
 
132
    g.SetTitle("; Position on Z axis (m); Field Magnitude (Tesla)")
 
133
    
 
134
    c = ROOT.TCanvas()
 
135
    g.Draw('ap')
 
136
    c.Print("Bfield_v_z.eps")
 
137
 
 
138
    c1 = ROOT.TCanvas()
 
139
    hf.Draw("colz")
 
140
    c1.Print("Bfield_v_z_hist.eps")
 
141
 
 
142
 
 
143
    c3 = ROOT.TCanvas()
 
144
    c3.SetLogz()
 
145
    hp.Draw('colz')
 
146
    c3.Print("Momentum_v_z_hist.eps")
 
147
 
 
148
    c3_1 = ROOT.TCanvas()
 
149
    hp_px = hp.ProjectionX()
 
150
    hp_px.Draw()
 
151
    c3_1.Print("Momentum_v_z_projx.eps")
 
152
 
 
153
    c3_2 = ROOT.TCanvas()
 
154
    hp_pfx = hp.ProfileX()
 
155
    hp_pfx.Draw()
 
156
    c3_2.Print("Momentum_v_z_prflx.eps")
 
157
 
 
158
    c3_3 = ROOT.TCanvas()
 
159
    c3_3.SetLogz()
 
160
    hp200.Draw('colz')
 
161
    c3_3.Print("Momentum_v_z_hist_lt200mm.eps")
 
162
 
 
163
    c3_4 = ROOT.TCanvas()
 
164
    hp200_px = hp200.ProjectionX()
 
165
    hp200_px.Draw()
 
166
    c3_4.Print("Momentum_v_z_projx_lt200mm.eps")
 
167
 
 
168
    c3_5 = ROOT.TCanvas()
 
169
    hp200_pfx = hp200.ProfileX()
 
170
    hp200_pfx.Draw()
 
171
    c3_5.Print("Momentum_v_z_prflx_lt200mm.eps")
 
172
 
 
173
    
 
174
    c4 = ROOT.TCanvas()
 
175
    c4.SetLogz()
 
176
    hr.Draw('colz')
 
177
    c4.Print("BeamRadius_v_z_hist.eps")
 
178
    
 
179
    c5 = ROOT.TCanvas()
 
180
    c5.SetLogz()
 
181
    hbt.Draw('colz')
 
182
    c5.Print("TransverseBeta_v_z_hist.eps")
 
183
 
 
184
    
 
185
    c6 = ROOT.TCanvas()
 
186
    c6.SetLogz()
 
187
    h.Draw()
 
188
    c6.Print("MuonRate_v_z_hist.eps")
 
189
 
 
190
 
 
191
if __name__ == "__main__":
 
192
    if len(sys.argv) > 1:
 
193
        main(sys.argv[1:])
 
194
    else:
 
195
        ofile = os.path.join\
 
196
                  (os.environ['MAUS_ROOT_DIR'],
 
197
                   'tmp','example_simulation_file.root')
 
198
        generate_simulation(ofile)
 
199
        main([ofile])
 
200