Project

General

Profile

Feature #1667 » ckov_momentum_scan.py

operation script - Bayes, Ryan, 07 April 2015 14:00

 
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
                    if tof10 > fp[2] and tof10 < fp[3]:
140
                        pmutof[-1].Fill( 105.66 / math.sqrt(( tof10 / 26.2 )**2 - 1) )
141
                        ckovaMuPE[-1].Fill(ckovA.GetNumberOfPes())
142
                        ckovbMuPE[-1].Fill(ckovB.GetNumberOfPes())
143
                        ptofmuvckova.Fill(105.66 / math.sqrt(( tof10 / 26.2 )**2 - 1),
144
                                          ckovA.GetNumberOfPes())
145
                        ptofmuvckovb.Fill(105.66 / math.sqrt(( tof10 / 26.2 )**2 - 1),
146
                                          ckovB.GetNumberOfPes())
147
                    elif tof10 > fp[4] and tof10 < fp[5]:
148
                        ppitof[-1].Fill( 139.57 / math.sqrt(( tof10 / 26.2 )**2 - 1) )
149
                        ckovaPiPE[-1].Fill(ckovA.GetNumberOfPes())
150
                        ckovbPiPE[-1].Fill(ckovB.GetNumberOfPes())
151
                        
152
                        ptofpivckova.Fill(139.57 / math.sqrt(( tof10 / 26.2 )**2 - 1),
153
                                          ckovA.GetNumberOfPes())
154
                        ptofpivckovb.Fill(139.57 / math.sqrt(( tof10 / 26.2 )**2 - 1),
155
                                          ckovB.GetNumberOfPes())
156
                    elif tof10 > 26.20666 and tof10 < fp[7]:
157
                        pmutof[-1].Fill( 0.511 / math.sqrt(( tof10 / 26.2 )**2 - 1) )
158
                        ckovaEPE[-1].Fill(ckovA.GetNumberOfPes())
159
                        ckovbEPE[-1].Fill(ckovB.GetNumberOfPes())
160
                        ptofevckova.Fill(fp[1], ckovA.GetNumberOfPes())
161
                        ptofevckovb.Fill(fp[1], ckovB.GetNumberOfPes())
162
                        
163

    
164

    
165
        CkovAMuPE.SetPoint(ifile, pmutof[-1].GetMean(), ckovaMuPE[-1].GetMean())
166
        CkovAMuPE.SetPointError(ifile, pmutof[-1].GetRMS()/math.sqrt(pmutof[-1].GetEntries()),
167
                                ckovaMuPE[-1].GetRMS()/math.sqrt(ckovaMuPE[-1].GetEntries()))
168
        CkovBMuPE.SetPoint(ifile, pmutof[-1].GetMean(), ckovbMuPE[-1].GetMean())
169
        CkovBMuPE.SetPointError(ifile, pmutof[-1].GetRMS()/math.sqrt(pmutof[-1].GetEntries()),
170
                                ckovbMuPE[-1].GetRMS()/math.sqrt(ckovbMuPE[-1].GetEntries()))
171

    
172
        CkovAPiPE.SetPoint(ifile, ppitof[-1].GetMean(), ckovaPiPE[-1].GetMean())
173
        CkovAPiPE.SetPointError(ifile, ppitof[-1].GetRMS()/math.sqrt(ppitof[-1].GetEntries()),
174
                                ckovaPiPE[-1].GetRMS()/math.sqrt(ckovaPiPE[-1].GetEntries()))
175
        CkovBPiPE.SetPoint(ifile, ppitof[-1].GetMean(), ckovbPiPE[-1].GetMean())
176
        CkovBPiPE.SetPointError(ifile, ppitof[-1].GetRMS()/math.sqrt(ppitof[-1].GetEntries()),
177
                                ckovbPiPE[-1].GetRMS()/math.sqrt(ckovbPiPE[-1].GetEntries()))
178

    
179
        CkovAEPE.SetPoint(ifile, fp[1], ckovaEPE[-1].GetMean())
180
        CkovAEPE.SetPointError(ifile, 0.0,
181
                               ckovaEPE[-1].GetRMS()/math.sqrt(ckovaEPE[-1].GetEntries()))
182
        CkovBEPE.SetPoint(ifile, fp[1], ckovbEPE[-1].GetMean())
183
        CkovBEPE.SetPointError(ifile, 0.0, 
184
                               ckovbEPE[-1].GetRMS()/math.sqrt(ckovbEPE[-1].GetEntries()))
185
        
186
        c1 = ROOT.TCanvas("tof_"+str(fp[1]),"tof_"+str(fp[1]))
187
        htof10[-1].Draw()
188
        c1.Print("ckov_momentum_scan_tof_"+str(fp[1])+"_"+str(ifile)+".eps")
189
        c1.Clear()
190
        pmutof[-1].Draw()
191
        c1.Print("ckov_scan_tof_muon_momentum_"+str(fp[1])+"_"+str(ifile)+".eps")
192
        c1.Clear()
193
        ppitof[-1].Draw()
194
        c1.Print("ckov_scan_tof_pion_momentum_"+str(fp[1])+"_"+str(ifile)+".eps")
195
        c1.Clear()
196
        petof[-1].Draw()
197
        c1.Print("ckov_scan_tof_electron_momentum_"+str(fp[1])+"_"+str(ifile)+".eps")
198
        c1.Clear()
199
        ckovaMuPE[-1].Draw()
200
        c1.Print("ckov_momentum_scan_ckova_"+str(fp[1])+"_"+str(ifile)+".eps")
201
        c1.Clear()
202
        ckovbMuPE[-1].Draw()
203
        c1.Print("ckov_momentum_scan_ckovb_"+str(fp[1])+"_"+str(ifile)+".eps")
204
        
205
        ifile += 1            
206

    
207
    CkovAMuPE.SetMarkerColor(2)
208
    CkovAMuPE.SetMarkerStyle(20)
209
    CkovAMuPE.SetTitle(";Momentum at target;#PE Ckov A")
210
    CkovBMuPE.SetMarkerColor(2)
211
    CkovBMuPE.SetMarkerStyle(20)
212
    CkovBMuPE.SetTitle(";Momentum at target;#PE Ckov B")
213
    CkovAPiPE.SetMarkerColor(4)
214
    CkovAPiPE.SetMarkerStyle(21)
215
    CkovBPiPE.SetMarkerColor(4)
216
    CkovBPiPE.SetMarkerStyle(21)
217
    CkovAEPE.SetMarkerColor(6)
218
    CkovAEPE.SetMarkerStyle(22)
219
    CkovBEPE.SetMarkerColor(6)
220
    CkovBEPE.SetMarkerStyle(22)
221
    leg = ROOT.TLegend(0.15,0.75,0.3,0.89)
222
    leg.AddEntry(CkovAMuPE,"#mu","p")
223
    leg.AddEntry(CkovAPiPE,"#pi","p")
224
    leg.AddEntry(CkovAEPE,"e","p")
225
    leg.SetFillColor(10)
226
    c1 = ROOT.TCanvas("ckov_scan","ckov_scan",800,900)
227
    c1.Divide(1,2)
228
    c1.cd(1)
229
    CkovAMuPE.Draw("ap")
230
    CkovAMuPE.GetHistogram().GetYaxis().SetRangeUser(0,30)
231
    CkovAPiPE.Draw("psame")
232
    CkovAEPE.Draw("psame")
233
    leg.Draw()
234
    c1.cd(2)
235
    CkovBMuPE.Draw("ap")
236
    CkovBMuPE.GetHistogram().GetYaxis().SetRangeUser(0,30)
237
    CkovBPiPE.Draw("psame")
238
    CkovBEPE.Draw("psame")
239

    
240
    c1.Print("ckov_scan_graph.eps")
241
    
242
    c1.Clear()
243
    c1.Divide(1,2)
244
    c1.cd(1)
245
    ptofmuvckova.Draw("colz")
246
    c1.cd(2)
247
    ptofmuvckovb.Draw("colz")
248
    c1.Print("ckov_scan_mu_pvPE.eps")
249
    
250
    c1.Clear()
251
    c1.Divide(1,2)
252
    c1.cd(1)
253
    ptofpivckova.Draw("colz")
254
    c1.cd(2)
255
    ptofpivckovb.Draw("colz")
256
    c1.Print("ckov_scan_pi_pvPE.eps")
257

    
258
    pxmua = ptofmuvckova.ProfileX()
259
    pxmua.SetMarkerColor(2)
260
    pxmua.SetMarkerStyle(20)
261
    pxmua.SetStats(0)
262
    pxmua.SetTitle(";Momentum from TOF;#PE Ckov A")
263
    pxmua.GetXaxis().SetLabelSize(0.05)
264
    pxmua.GetXaxis().SetTitleSize(0.05)
265
    pxmua.GetYaxis().SetLabelSize(0.05)
266
    pxmua.GetYaxis().SetTitleSize(0.05)
267
    pxmub = ptofmuvckovb.ProfileX()
268
    pxmub.SetMarkerColor(2)
269
    pxmub.SetMarkerStyle(20)
270
    pxmub.SetStats(0)
271
    pxmub.SetTitle(";Momentum from TOF;#PE Ckov B")
272
    pxmub.GetXaxis().SetLabelSize(0.05)
273
    pxmub.GetXaxis().SetTitleSize(0.05)
274
    pxmub.GetYaxis().SetLabelSize(0.05)
275
    pxmub.GetYaxis().SetTitleSize(0.05)
276
    pxpia = ptofpivckova.ProfileX()
277
    pxpia.SetMarkerColor(4)
278
    pxpia.SetMarkerStyle(21)
279
    pxpia.SetStats(0)
280
    pxpia.SetTitle(";Momentum from TOF;#PE Ckov A")
281
    pxpib = ptofpivckovb.ProfileX()
282
    pxpib.SetMarkerColor(4)
283
    pxpib.SetMarkerStyle(21)
284
    pxpib.SetStats(0)
285
    pxpib.SetTitle(";Momentum from TOF;#PE Ckov B")
286
    c1.Clear()
287
    c1.Divide(1,2)
288
    c1.cd(1)
289
    pxmua.Draw("p")
290
    pxmua.GetYaxis().SetRangeUser(0,20)
291
    pxpia.Draw("psame")
292
    leg.Draw()
293
    c1.cd(2)
294
    pxmub.Draw("p")
295
    pxmub.GetYaxis().SetRangeUser(0,20)
296
    pxpib.Draw("psame")
297
    c1.Print("ckov_scan_pvPE.eps")
298

    
299
if __name__ == "__main__":
300
    main()
301
    
(1-1/7)