Project

General

Profile

Feature #1667 » ckov_momentum_scan_v2.py

Updated analysis script using surveyed TOF1 position - Bayes, Ryan, 13 April 2015 11:24

 
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
    
(7-7/7)