« back to all changes in this revision
Viewing changes to src/common_py/analysis/scifi_efficiency.py
- browse files at revision 724
- compare with another revision
- download diff
- download tarball
- view history from revision 724
- Committer: Adam Dobbs
- Date: 2016-12-19 15:18:10 UTC
- mfrom: (659.2.53 release-candidate)
- Revision ID: phuccj@gmail.com-20161219151810-d4d1v2fas3tsjw1t
Tags: MAUS-v2.7, MAUS-v2.7.0
MAUS-v2.7.0
- tests/py_unit/test_analysis
- files removed:
- third_party/source/gtest-1.5.0.tar.gz.md5
- files modified:
- README
added
removed
src/common_py/analysis/scifi_efficiency.py
1
#!/usr/bin/env python
2
3
""" Check tracker efficiency """
4
5
import os
6
import math
7
# import abc
8
import ROOT
9
import libMausCpp #pylint: disable = W0611
10
import analysis.scifitools as tools
11
12
#pylint: disable = R0902
13
#pylint: disable = R0912
14
#pylint: disable = R0914
15
#pylint: disable = R0915
16
#pylint: disable = C0103
17
#pylint: disable = W0612
18
#pylint: disable = E1101
19
20
class EfficiencyDataReal():
21
""" Class to store pattern recognition efficiency data from real data """
22
def __init__(self):
23
""" Initialise member variables """
24
self.n_events_total = 0
25
self.n_passed_tof_timing_events = 0 # Events which pass the TOF cut
26
self.n_passed_tof_spoint_events = 0 # Events with 1 sp in TOF1 & TOF2
27
self.n_10spoint_events = 0 # Events with 5 spoints in both TkUS & TkDS
28
self.n_5spoint_tkus_events = 0 # Events with 5 spoints in TkUS
29
self.n_5spoint_tkds_events = 0 # Events with 5 spoints in TkDS
30
self.n_passed_tkus_events = 0 # Events which pass all the set cuts
31
self.n_passed_tkds_events = 0 # Events which pass all the set cuts
32
33
self.n_10spoint_tracks = 0
34
self.n_5spoint_tkus_tracks = 0
35
self.n_5spoint_tkds_tracks = 0
36
self.n_3to5spoint_tkus_tracks = 0
37
self.n_3to5spoint_tkds_tracks = 0
38
39
# Efficiency numbers
40
self.eff_tkus_5pt = 0.0
41
self.eff_tkus_5pt_err = 0.0
42
self.eff_tkus_3_5pt = 0.0
43
self.eff_tkus_3_5pt_err = 0.0
44
self.eff_tkds_5pt = 0.0
45
self.eff_tkds_5pt_err = 0.0
46
self.eff_tkds_3_5pt = 0.0
47
self.eff_tkds_3_5pt_err = 0.0
48
49
def __iadd__(self, other):
50
""" Add another instance's data to that here,
51
just the totals, not the final efficiency data"""
52
self.n_events_total += other.n_events_total
53
self.n_passed_tof_timing_events += other.n_passed_tof_timing_events
54
self.n_passed_tof_spoint_events += other.n_passed_tof_spoint_events
55
self.n_10spoint_events += other.n_10spoint_events
56
self.n_5spoint_tkus_events += other.n_5spoint_tkus_events
57
self.n_5spoint_tkds_events += other.n_5spoint_tkds_events
58
self.n_passed_tkus_events += other.n_passed_tkus_events
59
self.n_passed_tkds_events += other.n_passed_tkds_events
60
61
self.n_10spoint_tracks += other.n_10spoint_tracks
62
self.n_5spoint_tkus_tracks += other.n_5spoint_tkus_tracks
63
self.n_5spoint_tkds_tracks += other.n_5spoint_tkds_tracks
64
self.n_3to5spoint_tkus_tracks += other.n_3to5spoint_tkus_tracks
65
self.n_3to5spoint_tkds_tracks += other.n_3to5spoint_tkds_tracks
66
return self
67
68
def calculate_efficiency(self):
69
""" Calculate the final efficiency figures """
70
self.eff_tkus_5pt = 0.0
71
self.eff_tkus_5pt_err = 0.0
72
self.eff_tkus_3_5pt = 0.0
73
self.eff_tkus_3_5pt_err = 0.0
74
self.eff_tkds_5pt = 0.0
75
self.eff_tkds_5pt_err = 0.0
76
self.eff_tkds_3_5pt = 0.0
77
self.eff_tkds_3_5pt_err = 0.0
78
79
# Upstream tracker
80
try:
81
self.eff_tkus_5pt = float(self.n_5spoint_tkus_tracks) \
82
/float(self.n_passed_tkus_events)
83
self.eff_tkus_5pt_err = \
84
(self.eff_tkus_5pt * (1-self.eff_tkus_5pt)) \
85
/(float((self.n_passed_tkus_events)) ** (0.5))
86
except (ZeroDivisionError, ValueError):
87
self.eff_tkus_5pt = 0.0
88
self.eff_tkus_5pt_err = 0.0
89
try:
90
self.eff_tkus_3_5pt = float(self.n_3to5spoint_tkus_tracks) \
91
/float(self.n_passed_tkus_events)
92
self.eff_tkus_3_5pt_err = \
93
(self.eff_tkus_3_5pt * (1-self.eff_tkus_3_5pt)) \
94
/(float((self.n_passed_tkus_events)) ** (0.5))
95
self.eff_tkus_3_5pt_err = 0.001
96
except (ZeroDivisionError, ValueError):
97
self.eff_tkus_3_5pt = 0.0
98
self.eff_tkus_3_5pt_err = 0.0
99
100
# Downstream tracker
101
try:
102
self.eff_tkds_5pt = float(self.n_5spoint_tkds_tracks) \
103
/ float(self.n_passed_tkds_events)
104
self.eff_tkds_5pt_err = \
105
(self.eff_tkds_5pt * (1-self.eff_tkds_5pt)) \
106
/(float((self.n_passed_tkds_events)) ** (0.5))
107
except (ZeroDivisionError, ValueError):
108
self.eff_tkds_5pt = 0.0
109
self.eff_tkds_5pt_err = 0.0
110
try:
111
self.eff_tkds_3_5pt = float(self.n_3to5spoint_tkds_tracks) \
112
/float(self.n_passed_tkds_events)
113
self.eff_tkds_3_5pt_err = \
114
(self.eff_tkds_3_5pt * (1-self.eff_tkds_3_5pt)) \
115
/(float((self.n_passed_tkds_events)) ** (0.5))
116
#self.eff_tkds_3_5pt_err = 0.001
117
except (ZeroDivisionError, ValueError):
118
self.eff_tkds_3_5pt = 0.0
119
self.eff_tkds_3_5pt_err = 0.0
120
121
def clear(self):
122
""" Clear all data """
123
self.n_events_total = 0
124
self.n_passed_tof_timing_events = 0
125
self.n_passed_tof_spoint_events = 0
126
self.n_10spoint_events = 0
127
self.n_5spoint_tkus_events = 0
128
self.n_5spoint_tkds_events = 0
129
self.n_passed_tkus_events = 0
130
self.n_passed_tkds_events = 0
131
132
self.n_10spoint_tracks = 0
133
self.n_5spoint_tkus_tracks = 0
134
self.n_5spoint_tkds_tracks = 0
135
self.n_3to5spoint_tkus_tracks = 0
136
self.n_3to5spoint_tkds_tracks = 0
137
138
self.eff_tkus_5pt = 0.0
139
self.eff_tkus_5pt_err = 0.0
140
self.eff_tkus_3_5pt = 0.0
141
self.eff_tkus_3_5pt_err = 0.0
142
self.eff_tkds_5pt = 0.0
143
self.eff_tkds_5pt_err = 0.0
144
self.eff_tkds_3_5pt = 0.0
145
self.eff_tkds_3_5pt_err = 0.0
146
147
class EfficiencyDataMC():
148
""" Class to store pattern recognition efficiency data from MC data """
149
def __init__(self):
150
""" Initialise member variables """
151
self.n_events_total = 0 # Total number of events
152
self.n_mc_tku_tracks_total = 0 # Total number of tkus tracks
153
self.n_mc_tkd_tracks_total = 0 # Total number of tkus tracks
154
155
# MC tracks which created a viable # of spoints
156
self.n_mc_tku_tracks_5pt = 0
157
self.n_mc_tku_tracks_4to5pt = 0
158
self.n_mc_tku_tracks_3to5pt = 0
159
self.n_mc_tkd_tracks_5pt = 0
160
self.n_mc_tkd_tracks_4to5pt = 0
161
self.n_mc_tkd_tracks_3to5pt = 0
162
163
self.n_rec_tku_tracks_5pt_good = 0
164
self.n_rec_tkd_tracks_5pt_good = 0
165
self.n_rec_tku_tracks_4to5pt_good = 0
166
self.n_rec_tkd_tracks_4to5pt_good = 0
167
self.n_rec_tku_tracks_3to5pt_good = 0
168
self.n_rec_tkd_tracks_3to5pt_good = 0
169
170
self.eff_tkus_tracks = 0.0 # Efficiency of finding good tracks in tkus
171
self.eff_tkds_tracks = 0.0 # Efficiency of finding good tracks in tkds
172
173
def __iadd__(self, other):
174
""" Add another instance's data to that here,
175
just the totals, not the final efficiency data"""
176
self.n_events_total += other.n_events_total
177
self.n_mc_tku_tracks_total += other.n_mc_tku_tracks_total
178
self.n_mc_tkd_tracks_total += other.n_mc_tkd_tracks_total
179
180
self.n_mc_tku_tracks_5pt += other.n_mc_tku_tracks_5pt
181
self.n_mc_tku_tracks_4to5pt += other.n_mc_tku_tracks_4to5pt
182
self.n_mc_tku_tracks_3to5pt += other.n_mc_tku_tracks_3to5pt
183
184
self.n_mc_tkd_tracks_5pt += other.n_mc_tkd_tracks_5pt
185
self.n_mc_tkd_tracks_4to5pt += other.n_mc_tkd_tracks_4to5pt
186
self.n_mc_tkd_tracks_3to5pt += other.n_mc_tkd_tracks_3to5pt
187
188
self.n_rec_tku_tracks_5pt_good += other.n_rec_tku_tracks_5pt_good
189
self.n_rec_tkd_tracks_5pt_good += other.n_rec_tkd_tracks_5pt_good
190
self.n_rec_tku_tracks_4to5pt_good += other.n_rec_tku_tracks_4to5pt_good
191
self.n_rec_tkd_tracks_4to5pt_good += other.n_rec_tkd_tracks_4to5pt_good
192
self.n_rec_tku_tracks_3to5pt_good += other.n_rec_tku_tracks_3to5pt_good
193
self.n_rec_tkd_tracks_3to5pt_good += other.n_rec_tkd_tracks_3to5pt_good
194
195
return self
196
197
def calculate_efficiency(self):
198
""" Calculate the final efficiency figures """
199
self.eff_tkus_tracks = 0.0
200
self.eff_tkds_tracks = 0.0
201
202
# Upstream tracker
203
try:
204
self.eff_tkus_tracks = float(self.n_rec_tku_tracks_5pt_good) \
205
/ float(self.n_mc_tku_tracks_5pt)
206
except (ZeroDivisionError, ValueError):
207
self.eff_tkus_tracks = 0.0
208
# Downstream tracker
209
try:
210
self.eff_tkds_tracks = float(self.n_rec_tkd_tracks_5pt_good) \
211
/ float(self.n_mc_tkd_tracks_5pt)
212
except (ZeroDivisionError, ValueError):
213
self.eff_tkds_tracks = 0.0
214
215
def clear(self):
216
""" Clear all data """
217
self.n_events_total = 0
218
self.n_mc_tku_tracks_total = 0
219
self.n_mc_tkd_tracks_total = 0
220
221
self.n_mc_tku_tracks_5pt = 0
222
self.n_mc_tku_tracks_4to5pt = 0
223
self.n_mc_tku_tracks_3to5pt = 0
224
self.n_mc_tkd_tracks_5pt = 0
225
self.n_mc_tkd_tracks_4to5pt = 0
226
self.n_mc_tkd_tracks_3to5pt = 0
227
228
self.n_rec_tku_tracks_5pt_good = 0
229
self.n_rec_tkd_tracks_5pt_good = 0
230
self.n_rec_tku_tracks_4to5pt_good = 0
231
self.n_rec_tkd_tracks_4to5pt_good = 0
232
self.n_rec_tku_tracks_3to5pt_good = 0
233
self.n_rec_tkd_tracks_3to5pt_good = 0
234
235
self.eff_tkus_tracks = 0.0
236
self.eff_tkds_tracks = 0.0
237
238
class EfficiencyBase(object):
239
""" Base class for other Pattern Recognition Efficiency classes """
240
#__metaclass__ = abc.ABCMeta
241
242
def __init__(self):
243
""" Initialise member variables """
244
print 'Running init base ',
245
self.root_files = []
246
247
# Options for the user
248
self.cut_on_tof = True
249
self.cut_on_tof_time = True
250
self.tof_upper_cut = 50.0
251
self.tof_lower_cut = 27.0
252
print self.cut_on_tof
253
254
# Results
255
self.bool_2tof_timing_event = False # Tof timing ok
256
self.bool_2tof_spoint_event = False # 1 spacepoint only in each tof
257
258
def check_tof(self, tof_evt):
259
""" Analyse tof data. Return boolean indicating if tof cuts pass"""
260
try:
261
tof1 = tof_evt.GetTOFEventSpacePoint().GetTOF1SpacePointArray()
262
tof2 = tof_evt.GetTOFEventSpacePoint().GetTOF2SpacePointArray()
263
except ReferenceError:
264
print "Bad TOF data"
265
if self.cut_on_tof or self.cut_on_tof_time:
266
return False
267
# The data might be bad, but we aren't cutting on TOF so passes
268
return True
269
270
# Require 1 and only 1 sp in both TOF1 and TOF2
271
self.bool_2tof_spoint_event = True
272
if ((tof1.size() != 1) or (tof2.size() != 1)):
273
self.bool_2tof_spoint_event = False
274
275
# Require timing coincidence between TOF1 and TOF2
276
self.bool_2tof_timing_event = False
277
if self.bool_2tof_spoint_event:
278
for j in range(tof1.size()):
279
if (tof1.size() != 1) or (tof2.size() != 1):
280
continue
281
dt = tof2[j].GetTime() - tof1[j].GetTime()
282
283
if (dt < self.tof_upper_cut) and (dt > self.tof_lower_cut):
284
self.bool_2tof_timing_event = True
285
286
# Are the cuts choosen passed?
287
tof_good = True
288
if (not self.bool_2tof_spoint_event) and self.cut_on_tof:
289
tof_good = False
290
291
if (not self.bool_2tof_timing_event) and self.cut_on_tof_time:
292
tof_good = False
293
294
return tof_good
295
296
def run(self, files):
297
""" Loop over input root file, send each file for processing, print
298
summary information and produce any plots"""
299
self.print_welcome()
300
301
self.root_files = tools.load_data(files)
302
if len(self.root_files) < 1:
303
return False
304
305
counter = 0
306
for root_file_name in self.root_files:
307
self.process_file(root_file_name)
308
self.print_info(root_file_name, 'file')
309
self.accumulate_job_data()
310
counter += 1
311
self.calculate_job_efficiency()
312
self.print_info('Final job results', 'job')
313
return True
314
315
def process_file(self, root_file_name):
316
""" Process one root file of data"""
317
self.clear() # Start clean each time
318
319
# Load the ROOT file
320
root_file = ROOT.TFile(root_file_name, "READ") # pylint: disable = E1101
321
tree = root_file.Get("Spill")
322
data = ROOT.MAUS.Data() # pylint: disable = E1101
323
tree.SetBranchAddress("data", data)
324
325
# Loop over spills
326
for i in range(tree.GetEntries()):
327
tree.GetEntry(i)
328
spill = data.GetSpill()
329
self.process_spill(spill) # Analyse the spill
330
331
# Perform final efficiency calculations
332
self.calculate_file_efficiency()
333
334
# Close the ROOT file
335
root_file.Close()
336
337
#@abc.abstractmethod
338
def accumulate_job_data(self):
339
""" Accumulate the job data """
340
pass
341
342
#@abc.abstractmethod
343
def calculate_file_efficiency(self):
344
""" Calculate the file efficiency """
345
pass
346
347
#@abc.abstractmethod
348
def calculate_job_efficiency(self):
349
""" Calculate the job efficiency """
350
pass
351
352
#@abc.abstractmethod
353
def clear(self):
354
""" Set the internal data counters to zero & booleans to false """
355
pass
356
357
#@abc.abstractmethod
358
def print_info(self, data_name, info_type='file'):
359
""" Print the results """
360
pass
361
362
#@abc.abstractmethod
363
def print_welcome(self):
364
""" Message to be printed at programme start """
365
pass
366
367
#@abc.abstractmethod
368
def process_spill(self, spill):
369
""" Process one spill of data """
370
pass
371
372
class PatternRecognitionEfficiencyReal(EfficiencyBase):
373
""" Class to check Pattern Recognition efficiency with real data """
374
375
def __init__(self):
376
""" Initialise member variables """
377
EfficiencyBase.__init__(self)
378
379
# Options for the user
380
self.cut_on_tof = True
381
self.cut_on_tof_time = True
382
self.tof_upper_cut = 50.0
383
self.tof_lower_cut = 27.0
384
self.check_helical = True
385
self.check_straight = True
386
self.cut_on_tracker_10spnt = True
387
self.cut_on_tracker_5spnt = True # Should always be true
388
389
# Results
390
self.bool_2tof_timing_event = False # Tof timing ok
391
self.bool_2tof_spoint_event = False # 1 spacepoint only in each tof
392
self.bool_10spoint_event = False # Expect one 5pt track in each trcker
393
self.bool_tkus_5spoint_event = False # Expect one 5pt track in TkUS
394
self.bool_tkds_5spoint_event = False # Expect one 5pt track in TkDS
395
self.bool_tkus_1track = False # Found one track in TkUS
396
self.bool_tkds_1track = False # Found one track in TkDS
397
self.bool_tkus_5spoint_track = False # Found one 5pt track in TkUS
398
self.bool_tkds_5spoint_track = False # Found one 5pt track in TkDS
399
self.bool_passed_tkus_event = False # Have all the choosen cuts passed?
400
self.bool_passed_tkds_event = False # Have all the choosen cuts passed?
401
402
self.fdata = EfficiencyDataReal() # Per file data
403
self.jdata = EfficiencyDataReal() # Per job data
404
405
# Other objects not reset by the clear() function
406
self.n_print_calls = 0
407
408
def accumulate_job_data(self):
409
""" Accumulate the job data """
410
self.jdata += self.fdata
411
412
def calculate_file_efficiency(self):
413
""" Calculate the file efficiency """
414
self.fdata.calculate_efficiency()
415
416
def calculate_job_efficiency(self):
417
""" Calculate the job efficiency """
418
self.jdata.calculate_efficiency()
419
420
def check_tof(self, tof_evt):
421
""" Analyse tof data. Return boolean indicating if tof cuts pass"""
422
try:
423
tof1 = tof_evt.GetTOFEventSpacePoint().GetTOF1SpacePointArray()
424
tof2 = tof_evt.GetTOFEventSpacePoint().GetTOF2SpacePointArray()
425
except ReferenceError:
426
print "Bad TOF data"
427
if self.cut_on_tof or self.cut_on_tof_time:
428
return False
429
# The data might be bad, but we aren't cutting on TOF so passes
430
self.fdata.n_passed_tof_spoint_events += 1
431
self.fdata.n_passed_tof_timing_events += 1
432
return True
433
434
# Require 1 and only 1 sp in both TOF1 and TOF2
435
self.bool_2tof_spoint_event = True
436
if ((tof1.size() != 1) or (tof2.size() != 1)):
437
self.bool_2tof_spoint_event = False
438
439
# Require timing coincidence between TOF1 and TOF2
440
self.bool_2tof_timing_event = False
441
if self.bool_2tof_spoint_event:
442
for j in range(tof1.size()):
443
if (tof1.size() != 1) or (tof2.size() != 1):
444
continue
445
dt = tof2[j].GetTime() - tof1[j].GetTime()
446
447
if (dt < self.tof_upper_cut) and \
448
(dt > self.tof_lower_cut):
449
self.bool_2tof_timing_event = True
450
451
# Are the cuts choosen passed?
452
tof_good = True
453
if self.bool_2tof_spoint_event:
454
self.fdata.n_passed_tof_spoint_events += 1
455
elif self.cut_on_tof:
456
tof_good = False
457
else:
458
self.fdata.n_passed_tof_spoint_events += 1
459
460
if self.bool_2tof_timing_event:
461
self.fdata.n_passed_tof_timing_events += 1
462
elif self.cut_on_tof_time:
463
tof_good = False
464
else:
465
self.fdata.n_passed_tof_timing_events += 1
466
467
return tof_good
468
469
def print_welcome(self):
470
""" Message to be printed at programme start """
471
print '\nPattern Recognition Efficiency Calculator (Real Data)'
472
print '*******************************************************\n'
473
474
print 'Parameters:'
475
print 'Check helical\t' + str(self.check_helical)
476
print 'Check straight\t' + str(self.check_straight)
477
print 'Cut on TOF\t' + str(self.cut_on_tof)
478
print 'Cut on TOF Time\t' + str(self.cut_on_tof_time)
479
print 'Cut on trackers\t' + str(self.cut_on_tracker_10spnt) + '\n'
480
481
def process_spill(self, spill):
482
""" Process one spill of data """
483
#print "Spill number " + str(spill.GetSpillNumber())
484
if spill.GetDaqEventType() != "physics_event":
485
return False # remove event from consideration
486
487
# Loop over recon events
488
for i in range(spill.GetReconEvents().size()):
489
self.fdata.n_events_total += 1
490
491
# Pull out tof data and check if cuts pass
492
tof_evt = spill.GetReconEvents()[i].GetTOFEvent()
493
tof_good = self.check_tof(tof_evt)
494
if not tof_good:
495
continue # remove event from consideration
496
497
# Pull out tracker data and examine
498
tk_evt = spill.GetReconEvents()[i].GetSciFiEvent()
499
500
# Look at tracker spacepoint data if we expect tracks in either
501
# tracker
502
self.check_tracker_spacepoints(tk_evt.spacepoints())
503
504
# If we expect good tracks in neither tracker, skip this event
505
if ( not self.bool_passed_tkus_event ) and \
506
( not self.bool_passed_tkds_event ) :
507
continue # remove event from consideration
508
509
# Now switch from calculating expected tracks to what
510
# was actually reconstructed
511
512
# Extract the desired pat rec tracks (straight, helical or both)
513
prtracks = self.extract_prtracks(tk_evt)
514
515
# Analyse the pat rec tracks and see what we actually found
516
self.analyse_prtracks(prtracks)
517
518
return True
519
520
def analyse_prtracks(self, prtracks):
521
""" Analyse the reconstructed tracks, see how many were found when
522
they were expected """
523
524
# Require 1 and only 1 track in a tracker
525
n_tkus_tracks = 0
526
n_tkds_tracks = 0
527
self.bool_tkus_1track = False
528
self.bool_tkds_1track = False
529
for trk in prtracks:
530
if trk.get_tracker() == 0:
531
n_tkus_tracks += 1
532
elif trk.get_tracker() == 1:
533
n_tkds_tracks += 1
534
if n_tkus_tracks == 1:
535
self.bool_tkus_1track = True
536
if n_tkds_tracks == 1:
537
self.bool_tkds_1track = True
538
539
# If there were 5 spacepoints in a tracker and 1 track only was
540
# found in it, increment the 3 to 5 spoint track counter
541
if self.bool_tkus_5spoint_event and self.bool_tkus_1track:
542
self.fdata.n_3to5spoint_tkus_tracks += 1
543
if self.bool_tkds_5spoint_event and self.bool_tkds_1track:
544
self.fdata.n_3to5spoint_tkds_tracks += 1
545
546
# Now check the tracks found had 5 spoints, 1 from each station
547
self.bool_tkus_5spoint_track = False
548
self.bool_tkds_5spoint_track = False
549
for trk in prtracks:
550
if trk.get_tracker() == 0:
551
if (trk.get_spacepoints_pointers().size() == 5) \
552
and self.bool_tkus_5spoint_event and self.bool_tkus_1track:
553
self.fdata.n_5spoint_tkus_tracks += 1
554
self.bool_tkus_5spoint_track = True
555
elif trk.get_tracker() == 1:
556
if (trk.get_spacepoints_pointers().size() == 5) \
557
and self.bool_tkds_5spoint_event and self.bool_tkds_1track:
558
self.fdata.n_5spoint_tkds_tracks += 1
559
self.bool_tkds_5spoint_track = True
560
561
# Lastly see if we found only one 5 spoint track in BOTH trackers
562
if (len(prtracks) == 2 and self.bool_tkus_5spoint_track and \
563
self.bool_tkds_5spoint_track):
564
self.fdata.n_10spoint_tracks += 1
565
return True
566
567
568
def check_tracker_spacepoints(self, spoints):
569
""" Look for 5 spacepoint events in each tracker, increment the internal
570
good event counters, and return if the tracker spacepoint data
571
passes the selected cuts """
572
self.bool_10spoint_event = True
573
self.bool_tkus_5spoint_event = True
574
self.bool_tkds_5spoint_event = True
575
576
# Loop over trackers
577
for i in range(2):
578
tracker = [sp for sp in spoints if sp.get_tracker() == i]
579
# Loop over spacepoints
580
for j in range(1, 6):
581
station = \
582
[sp for sp in tracker if sp.get_station() == j]
583
# Does each station have one and only one sp
584
if len(station) != 1:
585
self.bool_10spoint_event = False
586
if i == 0:
587
self.bool_tkus_5spoint_event = False
588
elif i == 1:
589
self.bool_tkds_5spoint_event = False
590
591
# Update the internal expected tracks counters
592
if self.bool_10spoint_event:
593
self.fdata.n_10spoint_events += 1
594
if self.bool_tkus_5spoint_event:
595
self.fdata.n_5spoint_tkus_events += 1
596
if self.bool_tkds_5spoint_event:
597
self.fdata.n_5spoint_tkds_events += 1
598
599
# Are the complete set of tracker spacepoint cuts passed?
600
self.bool_passed_tkus_event = False
601
self.bool_passed_tkds_event = False
602
# Do we require 5 spoints in BOTH trackers?
603
if self.cut_on_tracker_10spnt:
604
if self.bool_10spoint_event:
605
self.bool_passed_tkus_event = True
606
self.fdata.n_passed_tkus_events += 1
607
self.bool_passed_tkds_event = True
608
self.fdata.n_passed_tkds_events += 1
609
else:
610
self.bool_passed_tkus_event = False
611
self.bool_passed_tkds_event = False
612
# If not, just see if have 5 spoints in each tracker individually
613
elif self.cut_on_tracker_5spnt:
614
if self.bool_tkus_5spoint_event:
615
self.bool_passed_tkus_event = True
616
self.fdata.n_passed_tkus_events += 1
617
else:
618
self.bool_passed_tkus_event = False
619
if self.bool_tkds_5spoint_event:
620
self.bool_passed_tkds_event = True
621
self.fdata.n_passed_tkds_events += 1
622
else:
623
self.bool_passed_tkds_event = False
624
# Not cutting on tracker spacepoints so everything passes
625
else:
626
self.bool_passed_tkus_event = True
627
self.fdata.n_passed_tkus_events += 1
628
self.bool_passed_tkds_event = True
629
self.fdata.n_passed_tkds_events += 1
630
631
def clear(self):
632
""" Set the internal file data counters to zero & booleans to false """
633
self.bool_2tof_timing_event = False
634
self.bool_2tof_spoint_event = False
635
self.bool_10spoint_event = False
636
self.bool_tkus_5spoint_event = False
637
self.bool_tkds_5spoint_event = False
638
self.bool_tkus_1track = False
639
self.bool_tkds_1track = False
640
self.bool_tkus_5spoint_track = False
641
self.bool_tkds_5spoint_track = False
642
self.bool_passed_tkus_event = False
643
self.bool_passed_tkds_event = False
644
645
self.fdata.clear()
646
647
def extract_prtracks(self, tk_evt):
648
""" Pull out the pattern recognition tracks selected for analysis -
649
either straight, helical or both """
650
prtracks = []
651
if (self.check_helical) and (not self.check_straight):
652
prtracks = tools.vector_to_list(tk_evt.helicalprtracks())
653
# print 'Looking at helical tracks'
654
elif (not self.check_helical) and (self.check_straight):
655
prtracks = tools.vector_to_list(tk_evt.straightprtracks())
656
# print 'Looking at straight tracks'
657
elif (self.check_helical) and (self.check_straight):
658
htracks = tools.vector_to_list(tk_evt.helicalprtracks())
659
stracks = tools.vector_to_list(tk_evt.straightprtracks())
660
prtracks = htracks + stracks
661
else:
662
print 'Warning: Both track type options not set'
663
return prtracks
664
665
#def make_plots(self, file_name_prefix):
666
#""" Draw the histogramms and save """
667
#c1 = ROOT.TCanvas()
668
#c1.Divide(2)
669
#c1.cd(1)
670
#self.pt_hist.Draw()
671
#c1.cd(2)
672
#self.pz_hist.Draw()
673
#c1.SaveAs(file_name_prefix + "_momentum_failed_events.pdf")
674
675
def print_info(self, data_name, info_type='file'):
676
""" Print the results """
677
if info_type == 'file':
678
edata = self.fdata
679
elif info_type == 'job':
680
edata = self.jdata
681
682
if self.n_print_calls == 0:
683
print '\nFile\t\t\tRecon_evt\tTOF\tTOF_SP\tTU_5pt' + \
684
'\tTD_5pt\tT_10pt\tGood_TU\tGood_TD\tTU_5trk' + \
685
'\tTU_5trk_err\tTU_3-5trk\tTU_3-5trk_err\t',
686
print 'TD_5trk\tTD_5trk_err\tTD_3-5trk\tTD_3-5trk_err'
687
688
print os.path.basename(data_name) + '\t',
689
print str(edata.n_events_total) + '\t\t',
690
print str(edata.n_passed_tof_timing_events) + '\t' + \
691
str(edata.n_passed_tof_spoint_events) + '\t' + \
692
str(edata.n_5spoint_tkus_events) + '\t' + \
693
str(edata.n_5spoint_tkds_events) + '\t' + \
694
str(edata.n_10spoint_events) + '\t' + \
695
str(edata.n_passed_tkus_events) + '\t' + \
696
str(edata.n_passed_tkds_events) + '\t',
697
698
f = '%.4f \t%.4f \t%.4f \t%.4f \t%.4f \t%.4f \t%.4f \t%.4f'
699
print f % (edata.eff_tkus_5pt, edata.eff_tkus_5pt_err, \
700
edata.eff_tkus_3_5pt, edata.eff_tkus_3_5pt_err, \
701
edata.eff_tkds_5pt, edata.eff_tkds_5pt_err, \
702
edata.eff_tkds_3_5pt, edata.eff_tkds_3_5pt_err)
703
704
self.n_print_calls += 1
705
706
class PatternRecognitionEfficiencyMC(EfficiencyBase):
707
""" Class to check Pattern Recognition efficiency with MC data """
708
709
def __init__(self):
710
""" Initialise member variables """
711
EfficiencyBase.__init__(self)
712
713
# Parameters
714
self.n_hits_cut = 5
715
716
# Internal counters
717
718
# Data storage
719
self.fdata = EfficiencyDataMC() # Per file data
720
self.jdata = EfficiencyDataMC() # Per job data
721
722
# Other objects not reset by the clear() function
723
self.n_print_calls = 0
724
725
self.hpt_m_tku = \
726
ROOT.TH1F('hpt_m_tku', 'pt of missed tracks TKU', 100, 0.0, 100.0)
727
self.hpt_m_tkd = \
728
ROOT.TH1F('hpt_m_tkd', 'pt of missed tracks TKD', 100, 0.0, 100.0)
729
self.hpz_m_tku = \
730
ROOT.TH1F('hpz_m_tku', 'pz of missed tracks TKU', 100, 0.0, 260.0)
731
self.hpz_m_tkd = \
732
ROOT.TH1F('hpz_m_tkd', 'pz of missed tracks TKD', 100, 0.0, 260.0)
733
734
def process_spill(self, spill):
735
""" Process one spill of data """
736
if spill.GetDaqEventType() != "physics_event":
737
return False # remove event from consideration
738
739
# print 'Processing spill ' + str(spill.GetSpillNumber())
740
741
# Loop over events
742
for i in range(spill.GetReconEvents().size()):
743
# print ' Processing event ' + str(i)
744
revent = spill.GetReconEvents()[i]
745
mcevent = spill.GetMCEvents()[i]
746
self.process_event(revent, mcevent)
747
748
def process_event(self, revent, mcevent):
749
""" Process the data for one event """
750
751
# Calculate the expected number of tracks by calculating how many
752
# MC tracks produced spacepoints in every station
753
# ***************************************************************
754
755
# Pull out tof data and check if cuts pass
756
tof_evt = revent.GetTOFEvent()
757
tof_good = self.check_tof(tof_evt)
758
if not tof_good:
759
return False # remove event from consideration
760
761
# Now the tracker data
762
bool_tkus_5spoint_event = False
763
bool_tkds_5spoint_event = False
764
765
sfevent = revent.GetSciFiEvent()
766
spoints_tku = \
767
[sp for sp in sfevent.spacepoints() if sp.get_tracker() == 0]
768
spoints_tkd = \
769
[sp for sp in sfevent.spacepoints() if sp.get_tracker() == 1]
770
#print ' Found ' +str(len(spoints_tku)) + ' spoints in TKU'
771
#print ' Found ' +str(len(spoints_tkd)) + ' spoints in TKD'
772
lkup = tools.SciFiLookup(mcevent)
773
774
tracks_tku_5pt = tools.find_mc_tracks_from_spoints(lkup, \
775
spoints_tku, 5, self.n_hits_cut)
776
tracks_tkd_5pt = tools.find_mc_tracks_from_spoints(lkup, \
777
spoints_tkd, 5, self.n_hits_cut)
778
tracks_tku_4to5pt = tools.find_mc_tracks_from_spoints(lkup, \
779
spoints_tku, 4, self.n_hits_cut)
780
tracks_tkd_4to5pt = tools.find_mc_tracks_from_spoints(lkup, \
781
spoints_tkd, 4, self.n_hits_cut)
782
tracks_tku_3to5pt = tools.find_mc_tracks_from_spoints(lkup, \
783
spoints_tku, 3, self.n_hits_cut)
784
tracks_tkd_3to5pt = tools.find_mc_tracks_from_spoints(lkup, \
785
spoints_tkd, 3, self.n_hits_cut)
786
787
# Only use the event for each detector if we have one 5pt track only
788
if len(tracks_tku_5pt) == 1 and len(tracks_tku_4to5pt) == 1 and \
789
len(tracks_tku_3to5pt) == 1:
790
bool_tkus_5spoint_event = True
791
else:
792
bool_tkus_5spoint_event = False
793
794
if len(tracks_tkd_5pt) == 1 and len(tracks_tkd_4to5pt) == 1 and \
795
len(tracks_tkd_3to5pt) == 1:
796
bool_tkds_5spoint_event = True
797
else:
798
bool_tkds_5spoint_event = False
799
800
# If we do not have good events in either tracker, save time and stop
801
if (not bool_tkus_5spoint_event) and (not bool_tkds_5spoint_event):
802
return False
803
804
#print 'Found ' + str(len(tracks_tku_5pt)) + ' 5pt MC tracks in TKU'
805
#print 'Found ' + str(len(tracks_tkd_5pt)) + ' 5pt MC tracks in TKD'
806
#print 'Found ' + str(len(tracks_tku_4to5pt)) + ' 4pt+ MC tracks in TKU'
807
#print 'Found ' + str(len(tracks_tkd_4to5pt)) + ' 4pt+ MC tracks in TKD'
808
#print 'Found ' + str(len(tracks_tku_3to5pt)) + ' 3pt+ MC tracks in TKU'
809
#print 'Found ' + str(len(tracks_tkd_3to5pt)) + ' 3pt+ MC tracks in TKD'
810
811
# Update internal counters
812
if bool_tkus_5spoint_event:
813
self.fdata.n_mc_tku_tracks_5pt += len(tracks_tku_5pt)
814
self.fdata.n_mc_tku_tracks_4to5pt += len(tracks_tku_4to5pt)
815
self.fdata.n_mc_tku_tracks_3to5pt += len(tracks_tku_3to5pt)
816
if bool_tkds_5spoint_event:
817
self.fdata.n_mc_tkd_tracks_5pt += len(tracks_tkd_5pt)
818
self.fdata.n_mc_tkd_tracks_4to5pt += len(tracks_tkd_4to5pt)
819
self.fdata.n_mc_tkd_tracks_3to5pt += len(tracks_tkd_3to5pt)
820
821
# Find the number of tracks actually reconstructed by pat rec
822
# ***********************************************************
823
824
stracks_tku = [trk for trk in sfevent.straightprtracks() \
825
if trk.get_tracker() == 0]
826
stracks_tkd = [trk for trk in sfevent.straightprtracks() \
827
if trk.get_tracker() == 1]
828
htracks_tku = [trk for trk in sfevent.helicalprtracks() \
829
if trk.get_tracker() == 0]
830
htracks_tkd = [trk for trk in sfevent.helicalprtracks() \
831
if trk.get_tracker() == 1]
832
833
# Do these tracks come from 1 MC track each? Are there mixed tracks?
834
tracks_tku = htracks_tku + stracks_tku
835
tracks_tkd = htracks_tkd + stracks_tkd
836
837
results_tku_5pt = self.check_tracks(lkup, 5, tracks_tku)
838
results_tkd_5pt = self.check_tracks(lkup, 5, tracks_tkd)
839
results_tku_4pt = self.check_tracks(lkup, 4, tracks_tku)
840
results_tkd_4pt = self.check_tracks(lkup, 4, tracks_tkd)
841
results_tku_3pt = self.check_tracks(lkup, 3, tracks_tku)
842
results_tkd_3pt = self.check_tracks(lkup, 3, tracks_tkd)
843
844
# Update internal counters
845
if bool_tkus_5spoint_event:
846
self.fdata.n_rec_tku_tracks_5pt_good += results_tku_5pt[0]
847
self.fdata.n_rec_tku_tracks_4to5pt_good += results_tku_4pt[0]
848
self.fdata.n_rec_tku_tracks_3to5pt_good += results_tku_3pt[0]
849
if bool_tkds_5spoint_event:
850
self.fdata.n_rec_tkd_tracks_5pt_good += results_tkd_5pt[0]
851
self.fdata.n_rec_tkd_tracks_4to5pt_good += results_tkd_4pt[0]
852
self.fdata.n_rec_tkd_tracks_3to5pt_good += results_tkd_3pt[0]
853
854
# Update histos, require that only 1 MC track is present
855
if bool_tkus_5spoint_event:
856
# Need the 2nd index to get track id rather than particle event id
857
track_id = tracks_tku_5pt[0]
858
mom = \
859
tools.find_mc_momentum_sfhits(lkup, spoints_tku, track_id, 0)
860
if (results_tku_3pt[0] < 1): # Fill if no rec track found
861
self.hpt_m_tku.Fill(math.sqrt(mom[0]**2 + mom[1]**2))
862
self.hpz_m_tku.Fill(mom[2])
863
864
if bool_tkds_5spoint_event:
865
track_id = tracks_tkd_5pt[0]
866
mom = \
867
tools.find_mc_momentum_sfhits(lkup, spoints_tkd, track_id, 1)
868
if (results_tkd_3pt[0] < 1): # Fill if no rec track found
869
self.hpt_m_tkd.Fill(math.sqrt(mom[0]**2 + mom[1]**2))
870
self.hpz_m_tkd.Fill(mom[2])
871
872
return True
873
874
def accumulate_job_data(self):
875
""" Accumulate the job data """
876
self.jdata += self.fdata
877
878
def calculate_file_efficiency(self):
879
""" Calculate the file efficiency """
880
self.fdata.calculate_efficiency()
881
882
def calculate_job_efficiency(self):
883
""" Calculate the job efficiency """
884
self.jdata.calculate_efficiency()
885
886
def check_tracks(self, lkup, nstations, recon_tracks):
887
""" Are the recon tracks found good tracks wrt the MC data?
888
nstations is the number of seed spacepoints originating from
889
the same MC track required for a recon track to be classed as
890
good. Hence if nstations = the actual number of tracker stations,
891
only perfect recon tracks will count as good """
892
n_good_tracks = 0 # tracks which can be associated with an MC track
893
# tracks generated from spoints from different MC tracks
894
n_mixed_tracks = 0
895
n_bad_tracks = 0 # tracks with mixed spoints or < expected spoints
896
for trk in recon_tracks:
897
spoints = trk.get_spacepoints_pointers()
898
mc_tids = tools.find_mc_tracks_from_spoints(lkup, spoints, \
899
nstations, self.n_hits_cut)
900
if len(mc_tids) == 1:
901
n_good_tracks = n_good_tracks + 1
902
elif len(mc_tids) > 1:
903
n_mixed_tracks = n_mixed_tracks + 1
904
elif len(mc_tids) < 1:
905
n_bad_tracks = n_bad_tracks + 1
906
return n_good_tracks, n_mixed_tracks, n_bad_tracks
907
908
def clear(self):
909
""" Set the internal file data counters to zero & booleans to false """
910
self.fdata.clear()
911
912
def print_info(self, data_name, info_type='file'):
913
""" Print the results """
914
if info_type == 'file':
915
edata = self.fdata
916
elif info_type == 'job':
917
edata = self.jdata
918
919
if self.n_print_calls == 0:
920
# print '\nRecon_evt\tTUMC5\tTURG5\tTUMC45\tTURG45\tTUMC35\tTURG35',
921
# print '\tTDMC5\tTDRG5\tTDMC45\tTDRG45\tTDMC35\tTDRG35\tFile'
922
print '\nRecon_evt\tTUMC5\tTURG5\tTURG45\tTURG35',
923
print '\tTDMC5\tTDRG5\tTDRG45\tTDRG35\tFile'
924
925
print str(edata.n_events_total) + '\t\t',
926
print str(edata.n_mc_tku_tracks_5pt) + '\t' + \
927
str(edata.n_rec_tku_tracks_5pt_good) + '\t' + \
928
str(edata.n_rec_tku_tracks_4to5pt_good) + '\t' + \
929
str(edata.n_rec_tku_tracks_3to5pt_good) + '\t' + \
930
str(edata.n_mc_tkd_tracks_5pt) + '\t' + \
931
str(edata.n_rec_tkd_tracks_5pt_good) + '\t' + \
932
str(edata.n_rec_tkd_tracks_4to5pt_good) + '\t' + \
933
str(edata.n_rec_tkd_tracks_3to5pt_good) + '\t' + \
934
os.path.basename(data_name)
935
self.n_print_calls += 1
936
937
def print_welcome(self):
938
""" Message to be printed at programme start """
939
print '\nPattern Recognition Efficiency Calculator (MC)'
940
print '************************************************\n'
941
942
# print 'Parameters:'
Loggerhead is a web-based interface for Breezy
Version: 2.0.1