Project

General

Profile

Feature #1983 » efficiency_analysis.py

Data-driven efficiency analysis - Hunt, Christopher, 05 November 2018 17:34

 
1
#!/bin/env /vols/build/mice/ch308/maus-bzr/maus_v3.2.0/third_party/install/bin/python
2

    
3
import MAUS
4

    
5
import sys
6
import os
7
import argparse
8

    
9
import json
10
import event_loader
11
import analysis
12
from analysis import tools
13
from analysis import covariances
14
from analysis import hit_types
15
import ROOT
16

    
17
NPE_CUT = 5.0
18

    
19
NEIGHBOUR_PT_CUT = 2000.0
20

    
21
TOF1_CUT_LOW = 0.0
22
TOF1_CUT_HIGH = 40.0
23

    
24
TOF2_CUT_LOW = 20.0
25
TOF2_CUT_HIGH = 40.0
26

    
27
TOF_PLOT1 = ROOT.TH1F('tof_plot_2', '', 400, 0.0, 40.0 )
28
TOF_PLOT2 = ROOT.TH1F('tof_plot_1', '', 400, 20.0, 60.0 )
29

    
30

    
31
def passes_cuts( tof_event ) :
32
  event_spacepoints = tof_event.GetTOFEventSpacePoint()
33

    
34
  tof0_sp_size = event_spacepoints.GetTOF0SpacePointArraySize()
35
  tof1_sp_size = event_spacepoints.GetTOF1SpacePointArraySize()
36
  tof2_sp_size = event_spacepoints.GetTOF2SpacePointArraySize()
37

    
38
  if tof0_sp_size < 1 :
39
    return False
40
  if tof1_sp_size != 1 or tof2_sp_size != 1 :
41
    return False
42

    
43
  tof0_sp = event_spacepoints.GetTOF0SpacePointArrayElement(0)
44
  tof1_sp = event_spacepoints.GetTOF1SpacePointArrayElement(0)
45
  tof2_sp = event_spacepoints.GetTOF2SpacePointArrayElement(0)
46

    
47
  diff_0_1 = tof1_sp.GetTime() - tof0_sp.GetTime()
48
  diff_1_2 = tof2_sp.GetTime() - tof1_sp.GetTime()
49

    
50
  TOF_PLOT1.Fill(diff_0_1)
51
  if diff_0_1 < TOF1_CUT_LOW or diff_0_1 > TOF1_CUT_HIGH :
52
    return False
53

    
54
  TOF_PLOT2.Fill(diff_1_2)
55
  if diff_1_2 < TOF2_CUT_LOW or diff_1_2 > TOF2_CUT_HIGH :
56
    return False
57

    
58
  return True
59

    
60

    
61

    
62
class Tracker_Analyser(object) :
63

    
64
  def __init__(self, tracker_no) :
65
    self.__tracker = tracker_no
66

    
67
    self.__event_counter = 0
68
    self.__good_event_counter = 0
69
    self.__bad_event_counter = 0
70
    self.__digit_counter = 0
71
    self.__cluster_counter = 0
72
    self.__spacepoint_counter = 0
73
    self.__spacepoint_CUT_counter = 0
74
    self.__used_spacepoint_counter = 0
75
    self.__helical_counter = 0
76
    self.__straight_counter = 0
77
    self.__helical_events = 0
78
    self.__straight_events = 0
79

    
80
    self.__expected_track_events = 0
81
    self.__no_expected_track_events = 0
82
    self.__expected_missing_track_events = 0
83

    
84
    self.__new_select_missed = 0
85
    self.__old_select_missed = 0
86

    
87
    self.__new_missed = 0
88
    self.__old_missed = 0
89
    self.__new_found = 0
90
    self.__old_found = 0
91

    
92
    self.__neighbour_found = 0
93
    self.__neighbour_missed = 0
94
    self.__neighbour_expected_missed = 0
95

    
96
    self.__plot_cluster_npe = ROOT.TH1F(str(self.__tracker)+"_cluster_npe", "", 100, 0.0, 100.0)
97
    self.__plot_spacepoint_npe = ROOT.TH1F(str(self.__tracker)+"_spacepoint_npe", "", 100, 0.0, 100.0)
98
    self.__plot_spacepoints = ROOT.TH1F(str(self.__tracker)+"_spacepoint", "", 101, -0.5, 100.5)
99
    self.__plot_spacepoints_CUT = ROOT.TH1F(str(self.__tracker)+"_spacepoint_CUT", "", 101, -0.5, 100.5)
100

    
101
    self.__plot_found_spacepoint_npe = ROOT.TH1F(str(self.__tracker)+"_found_spacepoint_npe", "", 100, 0.0, 100.0)
102
    self.__plot_missed_spacepoint_npe = ROOT.TH1F(str(self.__tracker)+"_missed_spacepoint_npe", "", 100, 0.0, 100.0)
103

    
104
    self.__plot_spacepoints_CUT_new = ROOT.TH1F(str(self.__tracker)+"_spacepoint_CUT_new", "", 101, -0.5, 100.5)
105
    self.__plot_spacepoints_CUT_old = ROOT.TH1F(str(self.__tracker)+"_spacepoint_CUT_old", "", 101, -0.5, 100.5)
106

    
107
    self.__plot_spacepoints_station_good = ROOT.TH2F(str(self.__tracker)+"_spacepoint_station_good", "", 5, 0.5, 5.5, 10, -0.5, 9.5)
108
    self.__plot_spacepoints_station_old = ROOT.TH2F(str(self.__tracker)+"_spacepoint_station_old", "", 5, 0.5, 5.5, 10, -0.5, 9.5)
109
    self.__plot_spacepoints_station_new = ROOT.TH2F(str(self.__tracker)+"_spacepoint_station_new", "", 5, 0.5, 5.5, 10, -0.5, 9.5)
110

    
111
    self.__plot_chisq_ndf = ROOT.TH1F(str(self.__tracker)+"_chisq_ndf", "", 500, 0.0, 100.0 )
112

    
113

    
114
  def analyse_event(self, scifi_event, passed) :
115
    digits = scifi_event.digits()
116
    clusters = scifi_event.clusters()
117
    spacepoints = scifi_event.spacepoints()
118

    
119
    helical_tracks = scifi_event.helicalprtracks()
120
    straight_tracks = scifi_event.straightprtracks()
121

    
122
    kalman_tracks = scifi_event.scifitracks()
123

    
124
    has_neighbours = False
125

    
126
    if self.__tracker == 0 :
127
      tracks_expected = scifi_event.get_expected_track_upstream()
128
    else :
129
      tracks_expected = scifi_event.get_expected_track_downstream()
130

    
131
    num_digits = 0
132
    num_clusters = 0
133
    num_spacepoints = 0
134
    num_spacepoints_CUT = 0
135
    num_used_spacepoints = 0
136
    num_helicals = 0
137
    num_straights = 0
138

    
139
    spacepoints_station = { 1: 0.0, 2: 0.0, 3: 0.0, 4: 0.0, 5: 0.0 }
140

    
141

    
142
    for helical in helical_tracks :
143
      if helical.get_tracker() == self.__tracker :
144
        num_helicals += 1
145
      elif helical.get_reference_momentum().Pt() < NEIGHBOUR_PT_CUT :
146
        has_neighbours = True
147

    
148
    for straight in straight_tracks :
149
      if straight.get_tracker() == self.__tracker :
150
        num_straights += 1
151

    
152
    total_tracks = num_straights + num_helicals
153

    
154
    for digit in digits :
155
      if digit.get_tracker() == self.__tracker :
156
        num_digits += 1
157

    
158
    for cluster in clusters :
159
      if cluster.get_tracker() == self.__tracker :
160
        num_clusters += 1
161

    
162
    for spacepoint in spacepoints :
163
      if spacepoint.get_tracker() == self.__tracker :
164
        num_spacepoints += 1
165
        spacepoints_station[ spacepoint.get_station() ] += 1.0
166

    
167
        if spacepoint.get_npe() > NPE_CUT :
168
          num_spacepoints_CUT += 1
169

    
170
        if spacepoint.get_used() :
171
          num_used_spacepoints += 1
172

    
173
    if num_spacepoints_CUT > 3 :
174
      self.__plot_spacepoints_CUT_new.Fill( num_spacepoints_CUT )
175
      for key, value in spacepoints_station.iteritems() :
176
        self.__plot_spacepoints_station_new.Fill( key, value )
177
    if total_tracks > 0 :
178
      self.__plot_spacepoints_CUT_old.Fill( num_spacepoints_CUT )
179
      for key, value in spacepoints_station.iteritems() :
180
        self.__plot_spacepoints_station_old.Fill( key, value )
181

    
182

    
183
    if passed :
184
      self.__event_counter += 1
185

    
186
      self.__helical_counter += num_helicals
187
      self.__straight_counter += num_straights
188
      self.__digit_counter += num_digits
189
      self.__cluster_counter += num_clusters
190
      self.__spacepoint_counter += num_spacepoints
191

    
192
      if has_neighbours :
193
        self.__neighbour_found += 1
194

    
195
      if total_tracks == 0 and tracks_expected :
196
        self.__expected_missing_track_events += 1
197

    
198
      if tracks_expected :
199
        self.__expected_track_events += 1
200
      else :
201
        self.__no_expected_track_events += 1
202

    
203
      if total_tracks == 0 and has_neighbours :
204
        self.__neighbour_missed += 1
205

    
206
      if total_tracks == 0 and has_neighbours and tracks_expected :
207
        self.__neighbour_expected_missed += 1
208

    
209

    
210
      for track in kalman_tracks :
211
        self.__plot_chisq_ndf.Fill(track.chi2()/track.ndf())
212

    
213
      for key, value in spacepoints_station.iteritems() :
214
        self.__plot_spacepoints_station_good.Fill( key, value )
215

    
216
      for cluster in clusters :
217
        if cluster.get_tracker() == self.__tracker :
218
          self.__plot_cluster_npe.Fill(cluster.get_npe())
219

    
220
      for spacepoint in spacepoints :
221
        if spacepoint.get_tracker() == self.__tracker :
222
          self.__plot_spacepoint_npe.Fill(spacepoint.get_npe())
223

    
224
          if total_tracks == 1 :
225
            if spacepoint.is_used() :
226
              self.__plot_found_spacepoint_npe.Fill( spacepoint.get_npe() )
227
            else :
228
              self.__plot_missed_spacepoint_npe.Fill( spacepoint.get_npe() )
229

    
230
      self.__used_spacepoint_counter += num_used_spacepoints
231

    
232
      if num_spacepoints_CUT > 3 :
233
        self.__good_event_counter += 1
234
        if total_tracks < 1 :
235
          self.__old_select_missed += 1
236
      else :
237
        if total_tracks > 0 :
238
          self.__new_select_missed += 1
239

    
240
      if num_helicals > 0 :
241
        self.__helical_events += 1
242
      if num_straights > 0 :
243
        self.__straight_events += 1
244

    
245
      self.__plot_spacepoints.Fill(num_spacepoints)
246
      self.__plot_spacepoints_CUT.Fill(num_spacepoints_CUT)
247

    
248
    else :
249
      self.__bad_event_counter += 1
250

    
251
      if num_spacepoints_CUT > 2 :
252
        if total_tracks < 1 :
253
          self.__old_missed += 1
254
        else :
255
          self.__old_found += 1
256
      else :
257
        if total_tracks > 0 :
258
          self.__new_missed += 1
259
        else :
260
          self.__new_found += 1
261

    
262

    
263

    
264
  def conclude(self) :
265
    data = ""
266
    data += "\nTracker Number : {0}".format(self.__tracker)
267
    data += "\n"
268
    data += "\nAnalysed {0} Events".format(self.__event_counter)
269
    data += "\n"
270
    data += "\nMean Numbers :"
271
    data += "\n  Digits      : {0}".format( float(self.__digit_counter)/float(self.__event_counter))
272
    data += "\n  Clusters    : {0}".format( float(self.__cluster_counter)/float(self.__event_counter))
273
    data += "\n  Spacepoints : {0}".format( float(self.__spacepoint_counter)/float(self.__event_counter))
274
    data += "\n  Helicals    : {0}".format( float(self.__helical_counter)/float(self.__event_counter))
275
    data += "\n  Straights   : {0}".format( float(self.__straight_counter)/float(self.__event_counter))
276
    data += "\n"
277
    data += "\n  Spacepoints > {0:3d}        : {1:2.4f}".format( int(NPE_CUT), float(self.__spacepoint_CUT_counter)/float(self.__event_counter) )
278
    data += "\n  Spacepoints > {0:3d} (Good) : {1:2.4f}".format( int(NPE_CUT), float(self.__spacepoint_CUT_counter)/float(self.__good_event_counter) )
279
    data += "\n"
280
    data += "\nTotal Numbers :"
281
    data += "\n  Digits      : {0}".format(self.__digit_counter)
282
    data += "\n  Clusters    : {0}".format(self.__cluster_counter)
283
    data += "\n  Spacepoints : {0}".format(self.__spacepoint_counter)
284
    data += "\n  Helicals    : {0}".format(self.__helical_counter)
285
    data += "\n  Straights   : {0}".format(self.__straight_counter)
286
    data += "\n"
287
    data += "\n  Events with 3xSpacepoints > {0:3d} : {1:3d}".format( int(NPE_CUT), self.__good_event_counter )
288
    data += "\n"
289
    data += "\n  Events with at least 1 helical  : {0}".format(self.__helical_events)
290
    data += "\n  Events with at least 1 straight : {0}".format(self.__straight_events)
291
    data += "\n  Efficiency = {0:2.2f}".format(100.0*(float(self.__helical_events) / float(self.__event_counter)))
292
    data += "\n"
293
    data += "\n Total Spacepoint Count = {0}".format(self.__used_spacepoint_counter)
294
    data += "\n Spacepoint Efficiency  = {0:2.2f}".format(100.0*self.__used_spacepoint_counter/float(5.0*self.__event_counter))
295
    data += "\n"
296
    data += "\nExpecting 3 Spacepoints:"
297
    data += "\n  New Selection Dropped : {0}".format(self.__new_select_missed)
298
    data += "\n  Old Selection Dropped : {0}".format(self.__old_select_missed)
299
    data += "\n"
300
    data += "\nEvents that didn't pass the cuts : {0}".format(self.__bad_event_counter)
301
    data += "\n"
302
    data += "\n New Selection Found/Dropped: {0}".format(self.__new_found, self.__new_missed)
303
    data += "\n Old Selection Found/Dropped: {0}".format(self.__old_found, self.__old_missed)
304
    data += "\n"
305
#    data += "\nEvents Expecting a track:         {0}".format(self.__expected_track_events)
306
#    data += "\nEvents Not Expecting a track:     {0}".format(self.__no_expected_track_events)
307
#    data += "\nEvents Missing an expected track: {0}".format(self.__expected_missing_track_events)
308
#    data += "\n"
309
    data += "\nEvents with a neighbour                 : {0}".format(self.__neighbour_found)
310
    data += "\nEvents with a neighbour, missing a track: {0}".format(self.__neighbour_missed)
311
#    data += "\nEvents with a neighbour, and expected, missing a track: {0}".format(self.__neighbour_expected_missed)
312
    data += "\nEfficiency = {0:2.2f}".format(100.0*float(self.__neighbour_found-self.__neighbour_missed) / float(self.__neighbour_found))
313
    data += "\n"
314
    data += "\n"
315
    data += "\n"
316
    return data
317

    
318

    
319
  def get_plots(self) :
320
    plot_dict = {}
321

    
322
    plot_dict['cluster_npe'] = self.__plot_cluster_npe
323
    plot_dict['spacepoint_npe'] = self.__plot_spacepoint_npe
324
    plot_dict['spacepoints'] = self.__plot_spacepoints
325
    plot_dict['spacepoints_CUT'] = self.__plot_spacepoints_CUT
326

    
327
    plot_dict['found_spacepoint_npe'] = self.__plot_found_spacepoint_npe
328
    plot_dict['missed_spacepoint_npe'] = self.__plot_missed_spacepoint_npe
329

    
330
    plot_dict['new_spacepoints_CUT'] = self.__plot_spacepoints_CUT_new
331
    plot_dict['old_spacepoints_CUT'] = self.__plot_spacepoints_CUT_old
332

    
333
    plot_dict['spacepoints_station_good'] = self.__plot_spacepoints_station_good
334
    plot_dict['spacepoints_station_new'] = self.__plot_spacepoints_station_new
335
    plot_dict['spacepoints_station_old'] = self.__plot_spacepoints_station_old
336

    
337
    return plot_dict
338

    
339

    
340

    
341
if __name__ == "__main__" : 
342
  ROOT.gROOT.SetBatch( True )
343
  ROOT.gErrorIgnoreLevel = ROOT.kError
344

    
345
  parser = argparse.ArgumentParser( description="" )
346

    
347
  parser.add_argument( 'maus_root_files', nargs='+', help='List of MAUS output root files containing reconstructed tracks')
348

    
349
  parser.add_argument( '--tof1_cut', nargs=2, type=float, default=[TOF1_CUT_LOW, TOF1_CUT_HIGH], help='')
350

    
351
  parser.add_argument( '--tof2_cut', nargs=2, type=float, default=[TOF2_CUT_LOW, TOF2_CUT_HIGH], help='')
352

    
353
  parser.add_argument( '-D', '--output_directory', default="./", help="Specify the directory in which to save the results.")
354
  parser.add_argument( '-O', '--output_filename', default="efficiency_analysis", help="Specify the file name to produce.")
355

    
356
  try :
357
    namespace = parser.parse_args()
358

    
359
    upstream_analyser = Tracker_Analyser(0)
360
    downstream_analyser = Tracker_Analyser(1)
361

    
362
    TOF1_CUT_LOW = namespace.tof1_cut[0]
363
    TOF1_CUT_HIGH = namespace.tof1_cut[1]
364

    
365
    TOF2_CUT_LOW = namespace.tof2_cut[0]
366
    TOF2_CUT_HIGH = namespace.tof2_cut[1]
367

    
368
  except BaseException as ex:
369
    raise
370
  else :
371

    
372

    
373
    file_reader = event_loader.maus_reader(namespace.maus_root_files)
374
#    file_reader.set_max_num_events(50000)
375

    
376
    try :
377
      while file_reader.next_selected_event() :
378

    
379
        try :
380
          scifi_event = file_reader.get_event( 'scifi' )
381
          tof_event = file_reader.get_event( 'tof' )
382

    
383
          if passes_cuts( tof_event ) :
384
            upstream_analyser.analyse_event( scifi_event, True )
385
            downstream_analyser.analyse_event( scifi_event, True )
386
          else :
387
            upstream_analyser.analyse_event( scifi_event, False )
388
            downstream_analyser.analyse_event( scifi_event, False )
389

    
390
        except ValueError as ex :
391
          print "An Error Occured: " + str(ex)
392
          print "Skipping Event: " +\
393
                str(file_reader.get_current_event_number()) + " In Spill: " + \
394
                str(file_reader.get_current_spill_number()) + " In File: " + \
395
                str(file_reader.get_current_filenumber()) + "\n"
396
          continue
397

    
398
    except KeyboardInterrupt :
399
      print
400
      print " ###  Keyboard Interrupt  ###"
401
      print
402
    print "- {0:0.0f} Spills Loaded                                 ".format( \
403
                                            file_reader.get_total_num_spills())
404

    
405

    
406
    print"\n- Analysing Data...\n"
407

    
408
    data = upstream_analyser.conclude()
409
    data += downstream_analyser.conclude()
410

    
411
    max_1 = TOF_PLOT1.GetMaximum()
412
    max_2 = TOF_PLOT2.GetMaximum()
413
    lower_line_1 = ROOT.TLine(TOF1_CUT_LOW, 0.0, TOF1_CUT_LOW, max_1)
414
    upper_line_1 = ROOT.TLine(TOF1_CUT_HIGH, 0.0, TOF1_CUT_HIGH, max_1)
415
    lower_line_2 = ROOT.TLine(TOF2_CUT_LOW, 0.0, TOF2_CUT_LOW, max_2)
416
    upper_line_2 = ROOT.TLine(TOF2_CUT_HIGH, 0.0, TOF2_CUT_HIGH, max_2)
417
    TOF_PLOT1.GetListOfFunctions().Add(lower_line_1)
418
    TOF_PLOT1.GetListOfFunctions().Add(upper_line_1)
419
    TOF_PLOT2.GetListOfFunctions().Add(lower_line_2)
420
    TOF_PLOT2.GetListOfFunctions().Add(upper_line_2)
421

    
422
    plots = { 'upstream' : upstream_analyser.get_plots(), 'downstream' : downstream_analyser.get_plots(), 'tof1' : TOF_PLOT1, 'tof2' : TOF_PLOT2 }
423

    
424
    filename = os.path.join(namespace.output_directory, namespace.output_filename)
425
    analysis.tools.save_plots(plots, filename+".root")
426
    with open(filename+".dat", "w") as outfile :
427
      outfile.write(data)
428
    
(1-1/2)