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 main():
|
19
|
|
20
|
""" extract information from the special virtual planes
|
21
|
from the maus_output.root file in the current directory """
|
22
|
|
23
|
args = [["tmp/06270/06270_recon.root",300,28.,29,30,31.,25.,27.],
|
24
|
["tmp/06271/06271_recon.root",300,28.,29,30.,31.,25.,27.],
|
25
|
["tmp/06274/06274_recon.root",350,27.,28.3,28.9,30.,25.,27.],
|
26
|
["tmp/06275/06275_recon.root",375,27.,28.,28.5,29.5,25,26.7],
|
27
|
["tmp/06276/06276_recon.root",400,26.7,27.8,28.0,29.0,26.,26.7],
|
28
|
["tmp/06278/06278_recon.root",425,26.5,27.5,27.8,28.8,26,26.5]]
|
29
|
|
30
|
maxsize = -999
|
31
|
|
32
|
htof10 = [] # array of histograms
|
33
|
pmutof = []
|
34
|
ppitof = []
|
35
|
petof = []
|
36
|
ckovaMuPE = []
|
37
|
ckovbMuPE = []
|
38
|
ckovaPiPE = []
|
39
|
ckovbPiPE = []
|
40
|
ckovaEPE = []
|
41
|
ckovbEPE = []
|
42
|
CkovAMuPE = ROOT.TGraphErrors(len(args))
|
43
|
CkovBMuPE = ROOT.TGraphErrors(len(args))
|
44
|
CkovAPiPE = ROOT.TGraphErrors(len(args))
|
45
|
CkovBPiPE = ROOT.TGraphErrors(len(args))
|
46
|
CkovAEPE = ROOT.TGraphErrors(len(args))
|
47
|
CkovBEPE = ROOT.TGraphErrors(len(args))
|
48
|
|
49
|
ptofmuvckova = ROOT.TH2D("pmutofckova" ,"; Momentum from TOF1 - TOF0 (MeV/c); Number of Photo-electrons",
|
50
|
80, 100, 500, 500, 0, 100)
|
51
|
ptofmuvckovb = ROOT.TH2D("pmutofckovb","; Momentum from TOF1 - TOF0 (MeV/c); Number of Photo-electrons",
|
52
|
80, 100, 500, 500, 0, 100)
|
53
|
|
54
|
ptofpivckova = ROOT.TH2D("ppitofckova" ,"; Momentum from TOF1 - TOF0 (MeV/c); Number of Photo-electrons",
|
55
|
80, 100, 500, 500, 0, 100)
|
56
|
ptofpivckovb = ROOT.TH2D("ppitofckovb","; Momentum from TOF1 - TOF0 (MeV/c); Number of Photo-electrons",
|
57
|
80, 100, 500, 500, 0, 100)
|
58
|
|
59
|
ptofevckova = ROOT.TH2D("petofckova" ,"; Momentum from TOF1 - TOF0 (MeV/c); Number of Photo-electrons",
|
60
|
80, 100, 500, 500, 0, 100)
|
61
|
ptofevckovb = ROOT.TH2D("petofckovb","; Momentum from TOF1 - TOF0 (MeV/c); Number of Photo-electrons",
|
62
|
80, 100, 500, 500, 0, 100)
|
63
|
|
64
|
ifile = 0
|
65
|
for fp in args:
|
66
|
rfile = ROOT.TFile(str(fp[0]),"R")
|
67
|
if rfile.IsZombie():
|
68
|
continue
|
69
|
data = ROOT.MAUS.Data()
|
70
|
tree = rfile.Get("Spill")
|
71
|
if not tree:
|
72
|
continue
|
73
|
tree.SetBranchAddress("data", data)
|
74
|
|
75
|
htof10.append( ROOT.TH1D("htof_" + str(fp[1])+"_"+str(ifile),
|
76
|
"TOF distribution at "+
|
77
|
str(fp[1])+"; TOF1 - TOF0 (ns); Entries (per 0.1 ns)", 300, 10, 40) )
|
78
|
|
79
|
pmutof.append( ROOT.TH1D("pmutof_" + str(fp[1])+"_"+str(ifile),
|
80
|
"Momenta from TOF distribution at "+
|
81
|
str(fp[1])+" for muons; Momentum (MeV/c); Entries (per 0.1 MeV)", 400, 100, 500) )
|
82
|
|
83
|
ppitof.append( ROOT.TH1D("ppitof_" + str(fp[1])+"_"+str(ifile),
|
84
|
"Momenta from TOF distribution at "+
|
85
|
str(fp[1])+" for pions; Momentum (MeV/c); Entries (per 0.1 MeV)", 400, 100, 500) )
|
86
|
|
87
|
petof.append( ROOT.TH1D("petof_" + str(fp[1])+"_"+str(ifile),
|
88
|
"Momenta from TOF distribution at "+
|
89
|
str(fp[1])+" for electrons; Momentum (MeV/c); Entries (per 0.1 MeV)", 400, 100, 500) )
|
90
|
|
91
|
ckovaPiPE.append( ROOT.TH1D("ckovaPi_" + str(fp[1])+"_"+str(ifile),
|
92
|
"Ckov A Pion PE distribution at "+
|
93
|
str(fp[1])+"; Number of Photo-Electrons; Entries (per 0.2 PE)", 500, 0, 100) )
|
94
|
|
95
|
ckovbPiPE.append( ROOT.TH1D("ckovbPi_" + str(fp[1])+"_"+str(ifile),
|
96
|
"Ckov B Pion PE distribution at "+
|
97
|
str(fp[1])+"; Number of Photo-Electrons; Entries (per 0.2 PE)", 500, 0, 100) )
|
98
|
|
99
|
ckovaMuPE.append( ROOT.TH1D("ckovaMu_" + str(fp[1])+"_"+str(ifile),
|
100
|
"Ckov A Muon PE distribution at "+
|
101
|
str(fp[1])+"; Number of Photo-Electrons; Entries (per 0.2 PE)", 500, 0, 100) )
|
102
|
|
103
|
ckovbMuPE.append( ROOT.TH1D("ckovbMu_" + str(fp[1])+"_"+str(ifile),
|
104
|
"Ckov B Muon PE distribution at "+
|
105
|
str(fp[1])+"; Number of Photo-Electrons; Entries (per 0.2 PE)", 500, 0, 100) )
|
106
|
|
107
|
ckovaEPE.append( ROOT.TH1D("ckovaE_" + str(fp[1])+"_"+str(ifile),
|
108
|
"Ckov A Electron PE distribution at "+
|
109
|
str(fp[1])+"; Number of Photo-Electrons; Entries (per 0.2 PE)", 500, 0, 100) )
|
110
|
|
111
|
ckovbEPE.append( ROOT.TH1D("ckovbE_" + str(fp[1])+"_"+str(ifile),
|
112
|
"Ckov B Electron PE distribution at "+
|
113
|
str(fp[1])+"; Number of Photo-Electrons; Entries (per 0.2 PE)", 500, 0, 100) )
|
114
|
|
115
|
|
116
|
print "Evaluating "+str(fp[0])
|
117
|
for i in range(tree.GetEntries()):
|
118
|
tree.GetEntry(i)
|
119
|
spill = data.GetSpill()
|
120
|
|
121
|
if spill.GetDaqEventType() == "physics_event":
|
122
|
for j in range(spill.GetReconEvents().size()):
|
123
|
tof_event = spill.GetReconEvents()[j].GetTOFEvent()
|
124
|
if not tof_event:
|
125
|
continue
|
126
|
spacepoints = tof_event.GetTOFEventSpacePoint()
|
127
|
if spacepoints.GetTOF0SpacePointArraySize() != 1:
|
128
|
continue
|
129
|
if spacepoints.GetTOF1SpacePointArraySize() != 1:
|
130
|
continue
|
131
|
tof0time = spacepoints.GetTOF0SpacePointArrayElement(0).GetTime()
|
132
|
tof1time = spacepoints.GetTOF1SpacePointArrayElement(0).GetTime()
|
133
|
tof10 = tof1time - tof0time
|
134
|
htof10[-1].Fill(tof10)
|
135
|
|
136
|
ckov_event = spill.GetReconEvents()[j].GetCkovEvent()
|
137
|
ckovA = ckov_event.GetCkovDigitArrayElement(0).GetCkovA()
|
138
|
ckovB = ckov_event.GetCkovDigitArrayElement(0).GetCkovB()
|
139
|
tof10_0 = 7630.1 / 299.7924
|
140
|
if tof10 > fp[2] and tof10 < fp[3]:
|
141
|
pmutof[-1].Fill( 105.66 / math.sqrt(( tof10 / tof10_0 )**2 - 1) )
|
142
|
ckovaMuPE[-1].Fill(ckovA.GetNumberOfPes())
|
143
|
ckovbMuPE[-1].Fill(ckovB.GetNumberOfPes())
|
144
|
ptofmuvckova.Fill(105.66 / math.sqrt(( tof10 / tof10_0 )**2 - 1),
|
145
|
ckovA.GetNumberOfPes())
|
146
|
ptofmuvckovb.Fill(105.66 / math.sqrt(( tof10 / tof10_0 )**2 - 1),
|
147
|
ckovB.GetNumberOfPes())
|
148
|
elif tof10 > fp[4] and tof10 < fp[5]:
|
149
|
ppitof[-1].Fill( 139.57 / math.sqrt(( tof10 / tof10_0 )**2 - 1) )
|
150
|
ckovaPiPE[-1].Fill(ckovA.GetNumberOfPes())
|
151
|
ckovbPiPE[-1].Fill(ckovB.GetNumberOfPes())
|
152
|
|
153
|
ptofpivckova.Fill(139.57 / math.sqrt(( tof10 / tof10_0 )**2 - 1),
|
154
|
ckovA.GetNumberOfPes())
|
155
|
ptofpivckovb.Fill(139.57 / math.sqrt(( tof10 / tof10_0 )**2 - 1),
|
156
|
ckovB.GetNumberOfPes())
|
157
|
elif tof10 > tof10_0 and tof10 < fp[7]:
|
158
|
pmutof[-1].Fill( 0.511 / math.sqrt(( tof10 / tof10_0 )**2 - 1) )
|
159
|
ckovaEPE[-1].Fill(ckovA.GetNumberOfPes())
|
160
|
ckovbEPE[-1].Fill(ckovB.GetNumberOfPes())
|
161
|
ptofevckova.Fill(fp[1], ckovA.GetNumberOfPes())
|
162
|
ptofevckovb.Fill(fp[1], ckovB.GetNumberOfPes())
|
163
|
|
164
|
|
165
|
|
166
|
CkovAMuPE.SetPoint(ifile, pmutof[-1].GetMean(), ckovaMuPE[-1].GetMean())
|
167
|
CkovAMuPE.SetPointError(ifile, pmutof[-1].GetRMS()/math.sqrt(pmutof[-1].GetEntries()),
|
168
|
ckovaMuPE[-1].GetRMS()/math.sqrt(ckovaMuPE[-1].GetEntries()))
|
169
|
CkovBMuPE.SetPoint(ifile, pmutof[-1].GetMean(), ckovbMuPE[-1].GetMean())
|
170
|
CkovBMuPE.SetPointError(ifile, pmutof[-1].GetRMS()/math.sqrt(pmutof[-1].GetEntries()),
|
171
|
ckovbMuPE[-1].GetRMS()/math.sqrt(ckovbMuPE[-1].GetEntries()))
|
172
|
|
173
|
CkovAPiPE.SetPoint(ifile, ppitof[-1].GetMean(), ckovaPiPE[-1].GetMean())
|
174
|
CkovAPiPE.SetPointError(ifile, ppitof[-1].GetRMS()/math.sqrt(ppitof[-1].GetEntries()),
|
175
|
ckovaPiPE[-1].GetRMS()/math.sqrt(ckovaPiPE[-1].GetEntries()))
|
176
|
CkovBPiPE.SetPoint(ifile, ppitof[-1].GetMean(), ckovbPiPE[-1].GetMean())
|
177
|
CkovBPiPE.SetPointError(ifile, ppitof[-1].GetRMS()/math.sqrt(ppitof[-1].GetEntries()),
|
178
|
ckovbPiPE[-1].GetRMS()/math.sqrt(ckovbPiPE[-1].GetEntries()))
|
179
|
|
180
|
CkovAEPE.SetPoint(ifile, fp[1], ckovaEPE[-1].GetMean())
|
181
|
CkovAEPE.SetPointError(ifile, 0.0,
|
182
|
ckovaEPE[-1].GetRMS()/math.sqrt(ckovaEPE[-1].GetEntries()))
|
183
|
CkovBEPE.SetPoint(ifile, fp[1], ckovbEPE[-1].GetMean())
|
184
|
CkovBEPE.SetPointError(ifile, 0.0,
|
185
|
ckovbEPE[-1].GetRMS()/math.sqrt(ckovbEPE[-1].GetEntries()))
|
186
|
|
187
|
c1 = ROOT.TCanvas("tof_"+str(fp[1]),"tof_"+str(fp[1]))
|
188
|
htof10[-1].Draw()
|
189
|
c1.Print("ckov_momentum_scan_tof_"+str(fp[1])+"_"+str(ifile)+".eps")
|
190
|
c1.Clear()
|
191
|
pmutof[-1].Draw()
|
192
|
c1.Print("ckov_scan_tof_muon_momentum_"+str(fp[1])+"_"+str(ifile)+".eps")
|
193
|
c1.Clear()
|
194
|
ppitof[-1].Draw()
|
195
|
c1.Print("ckov_scan_tof_pion_momentum_"+str(fp[1])+"_"+str(ifile)+".eps")
|
196
|
c1.Clear()
|
197
|
petof[-1].Draw()
|
198
|
c1.Print("ckov_scan_tof_electron_momentum_"+str(fp[1])+"_"+str(ifile)+".eps")
|
199
|
c1.Clear()
|
200
|
ckovaMuPE[-1].Draw()
|
201
|
c1.Print("ckov_momentum_scan_ckova_"+str(fp[1])+"_"+str(ifile)+".eps")
|
202
|
c1.Clear()
|
203
|
ckovbMuPE[-1].Draw()
|
204
|
c1.Print("ckov_momentum_scan_ckovb_"+str(fp[1])+"_"+str(ifile)+".eps")
|
205
|
|
206
|
ifile += 1
|
207
|
|
208
|
CkovAMuPE.SetMarkerColor(2)
|
209
|
CkovAMuPE.SetMarkerStyle(20)
|
210
|
CkovAMuPE.SetTitle(";Momentum at target;#PE Ckov A")
|
211
|
CkovBMuPE.SetMarkerColor(2)
|
212
|
CkovBMuPE.SetMarkerStyle(20)
|
213
|
CkovBMuPE.SetTitle(";Momentum at target;#PE Ckov B")
|
214
|
CkovAPiPE.SetMarkerColor(4)
|
215
|
CkovAPiPE.SetMarkerStyle(21)
|
216
|
CkovBPiPE.SetMarkerColor(4)
|
217
|
CkovBPiPE.SetMarkerStyle(21)
|
218
|
CkovAEPE.SetMarkerColor(6)
|
219
|
CkovAEPE.SetMarkerStyle(22)
|
220
|
CkovBEPE.SetMarkerColor(6)
|
221
|
CkovBEPE.SetMarkerStyle(22)
|
222
|
leg = ROOT.TLegend(0.15,0.75,0.3,0.89)
|
223
|
leg.AddEntry(CkovAMuPE,"#mu","p")
|
224
|
leg.AddEntry(CkovAPiPE,"#pi","p")
|
225
|
leg.AddEntry(CkovAEPE,"e","p")
|
226
|
leg.SetFillColor(10)
|
227
|
c1 = ROOT.TCanvas("ckov_scan","ckov_scan",800,900)
|
228
|
c1.Divide(1,2)
|
229
|
c1.cd(1)
|
230
|
CkovAMuPE.Draw("ap")
|
231
|
CkovAMuPE.GetHistogram().GetYaxis().SetRangeUser(0,30)
|
232
|
CkovAPiPE.Draw("psame")
|
233
|
CkovAEPE.Draw("psame")
|
234
|
leg.Draw()
|
235
|
c1.cd(2)
|
236
|
CkovBMuPE.Draw("ap")
|
237
|
CkovBMuPE.GetHistogram().GetYaxis().SetRangeUser(0,30)
|
238
|
CkovBPiPE.Draw("psame")
|
239
|
CkovBEPE.Draw("psame")
|
240
|
|
241
|
c1.Print("ckov_scan_graph.eps")
|
242
|
|
243
|
c1.Clear()
|
244
|
c1.Divide(1,2)
|
245
|
c1.cd(1)
|
246
|
ptofmuvckova.Draw("colz")
|
247
|
c1.cd(2)
|
248
|
ptofmuvckovb.Draw("colz")
|
249
|
c1.Print("ckov_scan_mu_pvPE.eps")
|
250
|
|
251
|
c1.Clear()
|
252
|
c1.Divide(1,2)
|
253
|
c1.cd(1)
|
254
|
ptofpivckova.Draw("colz")
|
255
|
c1.cd(2)
|
256
|
ptofpivckovb.Draw("colz")
|
257
|
c1.Print("ckov_scan_pi_pvPE.eps")
|
258
|
|
259
|
pxmua = ptofmuvckova.ProfileX()
|
260
|
pxmua.SetMarkerColor(2)
|
261
|
pxmua.SetMarkerStyle(20)
|
262
|
pxmua.SetStats(0)
|
263
|
pxmua.SetTitle(";Momentum from TOF;#PE Ckov A")
|
264
|
pxmua.GetXaxis().SetLabelSize(0.05)
|
265
|
pxmua.GetXaxis().SetTitleSize(0.05)
|
266
|
pxmua.GetYaxis().SetLabelSize(0.05)
|
267
|
pxmua.GetYaxis().SetTitleSize(0.05)
|
268
|
pxmub = ptofmuvckovb.ProfileX()
|
269
|
pxmub.SetMarkerColor(2)
|
270
|
pxmub.SetMarkerStyle(20)
|
271
|
pxmub.SetStats(0)
|
272
|
pxmub.SetTitle(";Momentum from TOF;#PE Ckov B")
|
273
|
pxmub.GetXaxis().SetLabelSize(0.05)
|
274
|
pxmub.GetXaxis().SetTitleSize(0.05)
|
275
|
pxmub.GetYaxis().SetLabelSize(0.05)
|
276
|
pxmub.GetYaxis().SetTitleSize(0.05)
|
277
|
pxpia = ptofpivckova.ProfileX()
|
278
|
pxpia.SetMarkerColor(4)
|
279
|
pxpia.SetMarkerStyle(21)
|
280
|
pxpia.SetStats(0)
|
281
|
pxpia.SetTitle(";Momentum from TOF;#PE Ckov A")
|
282
|
pxpib = ptofpivckovb.ProfileX()
|
283
|
pxpib.SetMarkerColor(4)
|
284
|
pxpib.SetMarkerStyle(21)
|
285
|
pxpib.SetStats(0)
|
286
|
pxpib.SetTitle(";Momentum from TOF;#PE Ckov B")
|
287
|
c1.Clear()
|
288
|
c1.Divide(1,2)
|
289
|
c1.cd(1)
|
290
|
pxmua.Draw("p")
|
291
|
pxmua.GetYaxis().SetRangeUser(0,20)
|
292
|
pxpia.Draw("psame")
|
293
|
leg.Draw()
|
294
|
c1.cd(2)
|
295
|
pxmub.Draw("p")
|
296
|
pxmub.GetYaxis().SetRangeUser(0,20)
|
297
|
pxpib.Draw("psame")
|
298
|
c1.Print("ckov_scan_pvPE.eps")
|
299
|
|
300
|
if __name__ == "__main__":
|
301
|
main()
|
302
|
|