← Back to branch summary

~maus-detector-integration/maus/ckov_merge

« back to all changes in this revision

Viewing changes to src/reduce/ReducePyCherenkov/ReducePyCherenkov.py

  • Committer: Chris Rogers
  • Date: 2012-03-05 22:08:21 UTC
  • mfrom: (340.51.8 maus)
  • Revision ID: chris.rogers@stfc.ac.uk-20120305220821-u428cxmzu537cdoi
Toy test for gene

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/reduce/ReducePyCherenkov/ReducePyCherenkov.py
 
1
import ROOT
 
2
from ReducePyROOTHistogram import ReducePyROOTHistogram
 
3
from math import sqrt
 
4
class ReducePyCherenkov(ReducePyROOTHistogram): # pylint: disable=R0902
 
5
 
 
6
    def __init__(self):
 
7
                
 
8
        ReducePyROOTHistogram.__init__(self)
 
9
    
 
10
        self.refresh_rate = 4
 
11
        self.root_batch_mode = 0
 
12
       
 
13
        self._hcharge = None
 
14
        self._htime = None
 
15
        self._htof = None
 
16
        self._htof_A = None
 
17
        self._htof_B = None
 
18
        self._hPMu = None
 
19
        self._hPPi = None
 
20
        self._hPMuA = None
 
21
        self._hPPiA = None
 
22
        self._hlight_Mu = None
 
23
        self._hEstCkov = None
 
24
 
 
25
        self.canvas_charge = None
 
26
        self.canvas_time = None
 
27
        self.canvas_tof = None
 
28
        self.canvas_tof_A = None
 
29
        self.canvas_tof_B = None
 
30
        self.canvas_hPMu = None
 
31
        self.canvas_hPPi = None
 
32
        self.canvas_hlight_Mu = None
 
33
        self.canvas_hPMuA = None
 
34
        self.canvas_hPPiA = None
 
35
        
 
36
    def _configure_at_birth(self, config_doc):
 
37
 
 
38
        if 'refresh_rate' in config_doc:
 
39
            self.refresh_rate = int(config_doc["refresh_rate"])
 
40
        self.__init_histos()
 
41
        return True
 
42
    
 
43
    def _update_histograms(self, spill):
 
44
 
 
45
        if "END_OF_RUN" in spill:
 
46
            self.update_histos()
 
47
            return self.get_histogram_images()
 
48
        
 
49
        if not self.get_space_points(spill):
 
50
            raise ValueError("space points not in spill")
 
51
 
 
52
        # Get charge points & fill histograms.
 
53
        if not self.get_charge(spill):
 
54
            raise ValueError("charge not in spill")
 
55
        
 
56
        # Refresh canvases at requested frequency.
 
57
    
 
58
        if self.spill_count % self.refresh_rate == 0:
 
59
            self.update_histos()
 
60
            return self.get_histogram_images()
 
61
        else:
 
62
            return [spill]
 
63
    
 
64
    def get_space_points(self, spill):
 
65
 
 
66
        if 'space_points' not in spill:
 
67
            return False
 
68
        
 
69
        # check for NoneType --there are events with no reconstructed SP
 
70
        if spill['space_points'] is None:
 
71
            return False
 
72
        
 
73
        space_points = spill['space_points']
 
74
        
 
75
        if 'tof0' not in space_points:
 
76
            return False
 
77
        sp_tof0 = space_points['tof0']
 
78
        
 
79
        if 'tof1' not in space_points:
 
80
            return False
 
81
        sp_tof1 = space_points['tof1']
 
82
        
 
83
        if 'tof2' not in space_points:
 
84
            return False
 
85
        sp_tof2 = space_points['tof2']
 
86
 
 
87
        if 'digits' not in spill:
 
88
            return False
 
89
 
 
90
        digits = spill['digits']
 
91
 
 
92
        #The cuts on e and mu TOFs are made first by 'eyeing' the 1D TOF distributions.
 
93
        
 
94
        tof_cut_e = 27.00
 
95
        tof_cut_mu = 31.00
 
96
 
 
97
        d = 7.8241 #Distance between TOF1 and TOF0
 
98
        
 
99
        MASS_MU = 105.658
 
100
        MASS_PI = 139.570
 
101
 
 
102
        CKOV_TURN_ON_107 = 2.6270345
 
103
        CKOV_TURN_ON_112 = 1.9826290
 
104
        
 
105
        m = 0
 
106
        p = 0
 
107
 
 
108
        #TOF vs. Total Charge
 
109
        
 
110
        for i in range(len(sp_tof1)):
 
111
            if sp_tof0[i] and sp_tof1[i] :
 
112
                if len(sp_tof0[i])==1 and len(sp_tof1[i])==1:
 
113
                    if sp_tof0[i][0]["part_event_number"]==digits[i]['B']['part_event_number']:
 
114
                        t_0 = sp_tof0[i][0]["time"]
 
115
                        t_1 = sp_tof1[i][0]["time"]
 
116
                        self._htof.Fill(t_1-t_0)
 
117
                        
 
118
                        charge_B = digits[i]["B"]["total_charge"] #This should be done pe and done in the mapper.
 
119
                        charge_A = digits[i]["A"]["total_charge"]
 
120
                        PE_B = charge_B/20.00 #Photoelectrons Normalized (can be done in mapper)
 
121
                        PE_A = charge_A/20.00 #Photoelectrons Normalized
 
122
                        if PE_B > 0:
 
123
                            self._htof_B.Fill(PE_B, (t_1 - t_0))
 
124
                        
 
125
                        if PE_A > 0:
 
126
                            self._htof_A.Fill(PE_A, (t_1-t_0))
 
127
                        
 
128
                        #Convert to momenta from TOF using m and d while making the appropriate cuts; fill momenta histograms
 
129
                        if ((t_1-t_0) < tof_cut_e):
 
130
                            p = -1.0
 
131
                        else:
 
132
                            if ((t_1-t_0) < tof_cut_mu):
 
133
                                m = MASS_MU
 
134
                            else:
 
135
                                m = MASS_PI
 
136
                            if (((t_1-t_0)*0.3) > d):
 
137
                                p = m / sqrt(((t_1-t_0)*0.3/d)*((t_1-t_0)*0.3/d)-1)
 
138
                            else:
 
139
                                p = -1.0
 
140
 
 
141
                        if (m==MASS_MU and p > 0.0 and PE_A > 1):
 
142
                            self._hPMu.Fill(p)
 
143
                            self._hlight_Mu.Fill(p, PE_A)
 
144
 
 
145
                        if (m==MASS_PI and p > 0.0):
 
146
                            self._hPPi.Fill(p)                                
 
147
 
 
148
                       # if (p > CKOV_TURN_ON_107*MASS_MU):
 
149
                       #     n_photons_est_ckov107 = 18*(1 - (CKOV_TURN_ON_107*MASS_MU/p)*(CKOV_TURN_ON_107*MASS_MU/p))
 
150
                       #     self._hEstCkov.Fill(p, n_phtons_est_ckov107)
 
151
                       # else:
 
152
                       #     n_photons_est_ckov107 = 0.0
 
153
                       #     self._hEstCkov.Fill(p, n_photons_est_ckov107)
 
154
                        
 
155
        return True
 
156
    
 
157
    def get_charge(self, spill):
 
158
 
 
159
        if 'digits' not in spill:
 
160
            return False
 
161
 
 
162
        digits = spill['digits']
 
163
                
 
164
        for i in range(len(digits)):
 
165
 
 
166
            for pmt in range(1,9):
 
167
                pulse = "pulse_%d" % (pmt)
 
168
                arrival_time = "arrival_time_%d" % (pmt)
 
169
 
 
170
                #CkovB
 
171
                if pmt < 5:
 
172
                    charge = digits[i]['B'][pulse]
 
173
                    if charge > -1:
 
174
                        self._hcharge[i-1].Fill(charge)
 
175
 
 
176
                    time = digits[i]['B'][arrival_time]
 
177
                    if time < 255:
 
178
                        self._htime[i-1].Fill(time)
 
179
                #CkovA
 
180
                if pmt > 4:
 
181
                    charge = digits[i]['A'][pulse]
 
182
                    if charge > -1:
 
183
                        self._hcharge[i-1].Fill(charge)
 
184
 
 
185
                    time = digits[i]['A'][arrival_time]
 
186
                    if time < 255:
 
187
                        self._htime[1-i].Fill(time)
 
188
               
 
189
        return True
 
190
    
 
191
    def __init_histos(self): #pylint: disable=R0201, R0914
 
192
        # have root run quietly without verbose informationals
 
193
        ROOT.gErrorIgnoreLevel = 1001
 
194
 
 
195
        # white canvas
 
196
        ROOT.gROOT.SetStyle("Plain")
 
197
        
 
198
        #turn off stat box
 
199
        ROOT.gStyle.SetOptStat(0)
 
200
        
 
201
        #sensible color palette
 
202
        ROOT.gStyle.SetPalette(1)
 
203
        
 
204
        # xy grid on canvas
 
205
        ROOT.gStyle.SetPadGridX(1)
 
206
        ROOT.gStyle.SetPadGridY(1)
 
207
 
 
208
        self._hcharge = []
 
209
        
 
210
        #c_x = 800
 
211
 
 
212
        #Make Canvas for TOF
 
213
        
 
214
        self.canvas_tof_A = ROOT.TCanvas("tof_A", "tof_A", 600, 600)
 
215
        self._htof_A = ROOT.TH2F("tof_A", "tof_A", 100, 0, 50, 200, 15, 35)
 
216
        self._htof_A.GetXaxis().SetTitle("#PE")
 
217
        self._htof_A.GetYaxis().SetTitle("TOF(CKOVA) (ns)")
 
218
        self.canvas_tof_A.cd()
 
219
        self._htof_A.Draw()
 
220
 
 
221
        self.canvas_tof_B = ROOT.TCanvas("tof_B", "tof_B", 600, 600)
 
222
        self._htof_B = ROOT.TH2F("tof_B", "tof_B", 100, 0, 50, 200, 15, 35)
 
223
        self._htof_B.GetXaxis().SetTitle("#PE")
 
224
        self._htof_B.GetYaxis().SetTitle("TOF(CKOVB) (ns)")
 
225
        self.canvas_tof_B.cd()
 
226
        self._htof_B.Draw()
 
227
 
 
228
        self.canvas_tof = ROOT.TCanvas("tof", "tof", 600, 600)
 
229
        self._htof =  ROOT.TH1F("htof", "htof", 100, 10, 40)
 
230
        self._htof.GetXaxis().SetTitle("TOF(ns)")
 
231
        self.canvas_tof.cd()
 
232
        self._htof.Draw()
 
233
                                        
 
234
        
 
235
        #Make Canvas for PMT Charges
 
236
        self.canvas_charge = ROOT.TCanvas("charge", "charge", 1200, 800)
 
237
        self.canvas_charge.Divide(4,2)
 
238
              
 
239
        for i in range(0,8):
 
240
            #ROOT.gPad.SetLogy()
 
241
 
 
242
            hname = "hPMT%d" % (i)
 
243
            htitle = "PMT%d" % (i)
 
244
            nbins = 200
 
245
            x_lo = 0
 
246
            x_hi = 200
 
247
 
 
248
            self._hcharge.append(ROOT.TH1F(hname, htitle, nbins, x_lo, x_hi))
 
249
            self._hcharge[i].GetXaxis().SetTitle("Charge")
 
250
            self.canvas_charge.cd(i+1)
 
251
            ROOT.gPad.SetLogy()
 
252
            self._hcharge[i].Draw()
 
253
 
 
254
        #Make Canvas for Arrival Times
 
255
 
 
256
        self._htime = []
 
257
        
 
258
        self.canvas_time = ROOT.TCanvas("time", "time", 1200, 800)
 
259
        self.canvas_time.Divide(4,2)
 
260
            
 
261
        for i in range(0,8):
 
262
            #ROOT.gPad.SetLogy()
 
263
            
 
264
            hname = "harr_time%d" % (i)
 
265
            htitle = "Arrival Times%d" % (i)
 
266
            nbins = 200
 
267
            x_lo = 0
 
268
            x_hi = 200
 
269
            
 
270
            self._htime.append(ROOT.TH1F(hname, htitle, nbins, x_lo, x_hi))
 
271
            self._htime[i].GetXaxis().SetTitle("Time (ns)")
 
272
            self.canvas_time.cd(i+1)
 
273
            self._htime[i].Draw()
 
274
 
 
275
        self.canvas_hPMu = ROOT.TCanvas("hPMu", "hPMu", 600, 600)
 
276
        self._hPMu = ROOT.TH1F("hPMu", "hPMu", 500, 0, 500)
 
277
        self._hPMu.GetXaxis().SetTitle("Momentum (MeV/c)")
 
278
        self.canvas_hPMu.cd()
 
279
        self._hPMu.Draw()
 
280
 
 
281
        self.canvas_hPPi = ROOT.TCanvas("hPPi", "hPPi", 600, 600)
 
282
        self._hPPi = ROOT.TH1F("hPPi", "hPPi", 500, 0, 500)
 
283
        self._hPPi.GetXaxis().SetTitle("Momentum (MeV/c)")
 
284
        self.canvas_hPPi.cd()
 
285
        self._hPPi.Draw()
 
286
 
 
287
        self.canvas_hlight_Mu = ROOT.TCanvas("PE_MOM", "PE_MOM", 600, 600)
 
288
        self._hlight_Mu = ROOT.TH2F("hLight", "hLight", 500, 0, 500, 100, 0, 30)
 
289
        self._hEstCkov = ROOT.TH2F("hLight", "hLight", 500, 0, 500, 100, 0, 30)
 
290
        self._hlight_Mu.GetXaxis().SetTitle("Momentum (MeV/c)")
 
291
        self._hlight_Mu.GetYaxis().SetTitle("#PE")
 
292
        self.canvas_hlight_Mu.cd()
 
293
        self._hlight_Mu.Draw()
 
294
        self._hEstCkov.Draw("same")            
 
295
        return True
 
296
    
 
297
    def update_histos(self):
 
298
        
 
299
        self.canvas_charge.Update()
 
300
        self.canvas_time.Update()
 
301
        self.canvas_tof.Update()
 
302
        self.canvas_tof_A.Update()
 
303
        self.canvas_tof_B.Update()
 
304
        self.canvas_hPMu.Update()
 
305
        self.canvas_hPPi.Update()
 
306
        self.canvas_hlight_Mu.Update()
 
307
        
 
308
    def get_histogram_images(self):
 
309
        """
 
310
        Get histograms as JSON documents.
 
311
        @param self Object reference.
 
312
        @returns list of 3 JSON documents containing the images.
 
313
        """
 
314
        image_list = []
 
315
        
 
316
        # PMT Charge
 
317
        # file label = PTM1-8.eps
 
318
        tag = "PMT1-8"
 
319
        content = "PMTCharge"
 
320
        doc = ReducePyROOTHistogram.get_image_doc(self, content, tag, self.canvas_charge)
 
321
        image_list.append(doc)
 
322
 
 
323
        #ArrivalTimes
 
324
        #file label = ArrivalTime.eps
 
325
        tag = "ArrivalTime"
 
326
        content = "arrival_time"
 
327
        doc = ReducePyROOTHistogram.get_image_doc(self, content, tag, self.canvas_time)
 
328
        image_list.append(doc)
 
329
 
 
330
        #TOF
 
331
        #file label = "TOF"
 
332
        tag = "TOF"
 
333
        content = "TOF"
 
334
        doc = ReducePyROOTHistogram.get_image_doc(self, content, tag, self.canvas_tof)
 
335
        image_list.append(doc)                                
 
336
 
 
337
        #TOFA and Charge
 
338
        #file label = "TOF_A"
 
339
        tag = "TOF_A"
 
340
        content = "TOF_A"
 
341
        doc = ReducePyROOTHistogram.get_image_doc(self, content, tag, self.canvas_tof_A)
 
342
        image_list.append(doc)
 
343
 
 
344
        #TOFB and Charge
 
345
        #file label = "TOF_B"
 
346
        tag = "TOF_B"
 
347
        content = "TOF_B"
 
348
        doc = ReducePyROOTHistogram.get_image_doc(self, content, tag, self.canvas_tof_B)
 
349
        image_list.append(doc)
 
350
 
 
351
        #Muon Momentum
 
352
        #file label = hPMu
 
353
        tag = "hPMu"
 
354
        content = "hPMu"
 
355
        doc = ReducePyROOTHistogram.get_image_doc(self, content, tag, self.canvas_hPMu)
 
356
        image_list.append(doc)
 
357
 
 
358
        #Pion Momentum
 
359
        #file label = hPPi
 
360
        tag = "hPPi"
 
361
        content = "hPPi"
 
362
        doc = ReducePyROOTHistogram.get_image_doc(self, content, tag, self.canvas_hPPi)
 
363
        image_list.append(doc)
 
364
 
 
365
        #PE vs Momentum
 
366
        #file label = PE_MOM
 
367
        tag = "PE_MOM"
 
368
        content = "PE_MOM"
 
369
        doc = ReducePyROOTHistogram.get_image_doc(self, content, tag, self.canvas_hlight_Mu)
 
370
        image_list.append(doc)
 
371
                                
 
372
        return image_list
 
373