Project

General

Profile

Feature #1983 » tracker_purity_analysis.py

MC-driven purity analysis - Hunt, Christopher, 05 November 2018 17:34

 
1

    
2
import MAUS
3

    
4
# Generic Python imports
5
import sys
6
import os
7
import argparse
8
import math
9
from math import sqrt
10
import array
11

    
12
# Third Party library import statements
13
import json
14
import event_loader
15
import analysis
16
from analysis import tools
17
from analysis import covariances
18
from analysis import hit_types
19
import ROOT
20

    
21

    
22
RECON_STATION = 1
23
RECON_PLANE = 0
24
SEED_STATION = 1
25
SEED_PLANE = 0
26
EXPECTED_STRAIGHT_TRACKPOINTS = 9
27
EXPECTED_HELIX_TRACKPOINTS = 12
28
REQUIRE_DATA = True
29
P_VALUE_CUT = 0.0
30
MUON_PID = [13, -13]
31

    
32
RECON_TRACKERS = [0, 1]
33

    
34
REQUIRE_ALL_PLANES = True
35

    
36
ANALYSE_REFITS = False
37

    
38
P_MIN = 0.0
39
P_MAX = 1000.0
40

    
41
PT_MIN = 0.0
42
PT_MAX = 100.0
43
PT_BIN = 10
44
PT_BIN_WIDTH = 10.0
45

    
46
PZ_MIN = 120.0
47
PZ_MAX = 260.0
48
PZ_BIN = 14
49
PZ_BIN_WIDTH = 10.0
50

    
51
ALIGNMENT_TOLERANCE = 0.1
52
RESOLUTION_BINS = 10
53
EFFICIENCY_BINS = 10
54

    
55
TRACK_ALGORITHM = 1
56

    
57
VIRTUAL_PLANE_DICT = None
58
INVERSE_PLANE_DICT = {}
59
TRACKER_PLANE_RADIUS = 150.0
60

    
61
SELECT_EVENTS = False
62
GOOD_EVENTS = None
63

    
64

    
65
def init_plots_data() :
66

    
67
  PZ_BIN = int(((PZ_MAX-PZ_MIN) / PZ_BIN_WIDTH) + 0.5)
68
  PT_BIN = int(((PT_MAX-PT_MIN) / PT_BIN_WIDTH) + 0.5)
69
  plot_dict = {'upstream' : {}, 'downstream' : {}}
70

    
71
  for tracker in [ 'upstream', 'downstream' ] :
72
    tracker_dict = {}
73
    tracker_dict['trackpoint_efficiency'] = ROOT.TEfficiency( tracker+'_trackpoint_efficiency', "Track Point Efficiency in P_{z} and P_{#perp}", PZ_BIN, PZ_MIN, PZ_MAX, PT_BIN, PT_MIN, PT_MAX )
74
    tracker_dict['trackpoint_efficiency_pt'] = ROOT.TEfficiency( tracker+'_trackpoint_efficiency_pt', "Track Point Efficiency in P_{#perp}", PT_BIN, PT_MIN, PT_MAX )
75
    tracker_dict['trackpoint_efficiency_pz'] = ROOT.TEfficiency( tracker+'_trackpoint_efficiency_pz', "Track Point Efficiency in P_z", PZ_BIN, PZ_MIN, PZ_MAX )
76

    
77
    tracker_dict['spacepoint_efficiency'] = ROOT.TEfficiency( tracker+'_spacepoint_efficiency', "Track Point Efficiency in P_{z} and P_{#perp}", PZ_BIN, PZ_MIN, PZ_MAX, PT_BIN, PT_MIN, PT_MAX )
78
    tracker_dict['spacepoint_efficiency_pt'] = ROOT.TEfficiency( tracker+'_spacepoint_efficiency_pt', "Track Point Efficiency in P_{#perp}", PT_BIN, PT_MIN, PT_MAX )
79
    tracker_dict['spacepoint_efficiency_pz'] = ROOT.TEfficiency( tracker+'_spacepoint_efficiency_pz', "Track Point Efficiency in P_z", PZ_BIN, PZ_MIN, PZ_MAX )
80

    
81
    tracker_dict['track_efficiency'] = ROOT.TEfficiency( tracker+'_track_efficiency', "Track Efficiency in P_z and P_{#perp}",  PZ_BIN, PZ_MIN, PZ_MAX, PT_BIN, PT_MIN, PT_MAX )
82
    tracker_dict['track_efficiency_pt'] = ROOT.TEfficiency(  tracker+'_track_efficiency_pt', "Track Efficiency in P_{#perp}", PT_BIN, PT_MIN, PT_MAX )
83
    tracker_dict['track_efficiency_pz'] = ROOT.TEfficiency( tracker+'_track_efficiency_pz', "Track Efficiency in P_z", PZ_BIN, PZ_MIN, PZ_MAX )
84

    
85
    tracker_dict['trackpoint_purity'] = ROOT.TEfficiency( tracker+'_trackpoint_purity', "Track Efficiency in P_z and P_{#perp}",  PZ_BIN, PZ_MIN, PZ_MAX, PT_BIN, PT_MIN, PT_MAX )
86
    tracker_dict['trackpoint_purity_pt'] = ROOT.TEfficiency(  tracker+'_trackpoint_purity_pt', "Track Efficiency in P_{#perp}", PT_BIN, PT_MIN, PT_MAX )
87
    tracker_dict['trackpoint_purity_pz'] = ROOT.TEfficiency( tracker+'_trackpoint_purity_pz', "Track Efficiency in P_z", PZ_BIN, PZ_MIN, PZ_MAX )
88

    
89
    tracker_dict['spacepoint_purity'] = ROOT.TEfficiency( tracker+'_spacepoint_purity', "Track Efficiency in P_z and P_{#perp}",  PZ_BIN, PZ_MIN, PZ_MAX, PT_BIN, PT_MIN, PT_MAX )
90
    tracker_dict['spacepoint_purity_pt'] = ROOT.TEfficiency(  tracker+'_spacepoint_purity_pt', "Track Efficiency in P_{#perp}", PT_BIN, PT_MIN, PT_MAX )
91
    tracker_dict['spacepoint_purity_pz'] = ROOT.TEfficiency( tracker+'_spacepoint_purity_pz', "Track Efficiency in P_z", PZ_BIN, PZ_MIN, PZ_MAX )
92

    
93
    plot_dict[tracker] = tracker_dict
94

    
95
  data_dict = { 'counters' : {'upstream' : {}, 'downstream' : {} }, \
96
                                                                  'data' : {} }
97
  data_dict['counters']['number_events'] = 0
98

    
99
  for tracker in ['upstream', 'downstream'] :
100
    data_dict['counters'][tracker]['number_virtual'] = 0
101
    data_dict['counters'][tracker]['missing_virtuals'] = 0
102

    
103
    data_dict['counters'][tracker]['number_tracks'] = 0
104
    data_dict['counters'][tracker]['number_candidates'] = 0
105
    data_dict['counters'][tracker]['found_tracks'] = 0
106
    data_dict['counters'][tracker]['wrong_track_type'] = 0
107
    data_dict['counters'][tracker]['p_value_cut'] = 0
108
    data_dict['counters'][tracker]['superfluous_track_events'] = 0
109
    data_dict['counters'][tracker]['missing_tracks'] = 0
110
    data_dict['counters'][tracker]['missing_reference_hits'] = 0
111

    
112
    data_dict['counters'][tracker]['momentum_cut'] = 0
113
    data_dict['counters'][tracker]['gradient_cut'] = 0
114

    
115
    data_dict['counters'][tracker]['found_pairs'] = 0
116

    
117
  return plot_dict, data_dict
118

    
119

    
120
def create_virtual_plane_dict(file_reader) :
121
  """
122
    Matches up scifitrackpoints to virtual planes to make a lookup dictionary 
123
  """
124
  virtual_plane_dict = {}
125
  for num in range( -15, 0, 1 ) :
126
    virtual_plane_dict[ num ] = ( -1, (ALIGNMENT_TOLERANCE * 100.0) )
127
  for num in range( 1, 16, 1 ) :
128
    virtual_plane_dict[ num ] = ( -1, (ALIGNMENT_TOLERANCE * 100.0) )
129

    
130
  while file_reader.next_event() :
131
    scifi_event = file_reader.get_event( 'scifi' )
132
    mc_event = file_reader.get_event( 'mc' )
133

    
134
    tracks = scifi_event.scifitracks()
135
    for track in tracks :
136
      if track.tracker() not in RECON_TRACKERS :
137
        continue
138
      trackpoints = track.scifitrackpoints()
139
      for trkpt in trackpoints :
140
        z_pos = trkpt.pos().z()
141
        plane_id = analysis.tools.calculate_plane_id(\
142
                               trkpt.tracker(), trkpt.station(), trkpt.plane())
143

    
144
        for vhit_num in xrange(mc_event.GetVirtualHitsSize()) :
145
          vhit = mc_event.GetAVirtualHit(vhit_num)
146
          diff = math.fabs(vhit.GetPosition().z() - z_pos)
147

    
148
          if diff < virtual_plane_dict[ plane_id ][1] :
149
            virtual_plane_dict[ plane_id ] = ( vhit.GetStationId(), diff )
150

    
151
    done = True
152
    for tracker in RECON_TRACKERS :
153
      for station in [1, 2, 3, 4, 5] :
154
        for plane in [0, 1, 2] :
155
          plane_id = analysis.tools.calculate_plane_id( \
156
                                                      tracker, station, plane )
157
          if virtual_plane_dict[plane_id][1] > ALIGNMENT_TOLERANCE :
158
#            print plane_id, virtual_plane_dict[plane]
159
            done = False
160
    if done :
161
      break
162
  else :
163
    if REQUIRE_ALL_PLANES :
164
      print
165
      print virtual_plane_dict
166
      raise ValueError("Could not locate all virtuals planes")
167

    
168
  file_reader.reset()
169
  return virtual_plane_dict
170

    
171

    
172
def inverse_virtual_plane_dict(virtual_plane_dict) :
173
  """
174
    Create the inverse lookup.
175
  """
176
  inverse_dict = {}
177
  for num in range( -15, 0, 1 ) :
178
    inverse_dict[virtual_plane_dict[num][0]] = num
179
  for num in range( 1, 16, 1 ) :
180
    inverse_dict[virtual_plane_dict[num][0]] = num
181

    
182
  return inverse_dict
183

    
184

    
185
def get_expected_tracks(mc_event, virtual_plane_dict) :
186
  upstream_planes = [ virtual_plane_dict[i][0] for i in range(-15, 0)]
187
  downstream_planes = [ virtual_plane_dict[i][0] for i in range(1, 16)]
188

    
189
  upstream_track = None
190
  downstream_track = None
191

    
192
  upstream_hits = {}
193
  downstream_hits = {}
194

    
195
  for vhit_num in xrange(mc_event.GetVirtualHitsSize()) :
196
    vhit = mc_event.GetAVirtualHit(vhit_num)
197
    if vhit.GetParticleId() not in MUON_PID :
198
      continue
199

    
200
    station_id = vhit.GetStationId()
201
    radius = math.sqrt( vhit.GetPosition().x()**2 + vhit.GetPosition().y()**2 )
202
    if radius > TRACKER_PLANE_RADIUS : 
203
      continue
204

    
205
    if station_id in upstream_planes :
206
      plane_id = INVERSE_PLANE_DICT[station_id]
207
      upstream_hits[plane_id] = vhit
208
    if station_id in downstream_planes :
209
      plane_id = INVERSE_PLANE_DICT[station_id]
210
      downstream_hits[plane_id] = vhit
211

    
212
  if TRACK_ALGORITHM == 1 :
213
    if len(upstream_hits) > EXPECTED_HELIX_TRACKPOINTS :
214
      upstream_track = upstream_hits
215
    if len(downstream_hits) > EXPECTED_HELIX_TRACKPOINTS :
216
      downstream_track = downstream_hits
217
  elif TRACK_ALGORITHM == 0 :
218
    if len(upstream_hits) > EXPECTED_STRAIGHT_TRACKPOINTS :
219
      upstream_track = upstream_hits
220
    if len(downstream_hits) > EXPECTED_STRAIGHT_TRACKPOINTS :
221
      downstream_track = downstream_hits
222
  else:
223
    raise ValueError("Unknown track algorithm found!")
224

    
225
  return upstream_track, downstream_track
226

    
227

    
228
def get_found_tracks(scifi_event, plot_dict, data_dict) :
229
  """
230
    Find all the single tracks that pass the cuts.
231
  """
232
  upstream_tracks = []
233
  downstream_tracks = []
234

    
235
  tracks = scifi_event.scifitracks()
236
  for track in tracks :
237
    if track.tracker() == 0 :
238
      tracker = "upstream"
239
    else :
240
      tracker = "downstream"
241

    
242
    data_dict['counters'][tracker]['number_tracks'] += 1
243

    
244
    if track.GetAlgorithmUsed() != TRACK_ALGORITHM :
245
      data_dict['counters'][tracker]['wrong_track_type'] += 1
246
      continue
247

    
248
    if track.P_value() < P_VALUE_CUT :
249
      data_dict['counters'][tracker]['p_value_cut'] += 1
250
      continue
251

    
252
    data_dict['counters'][tracker]['number_candidates'] += 1
253

    
254

    
255
    if track.tracker() == 0 :
256
      upstream_tracks.append(track)
257
    if track.tracker() == 1 :
258
      downstream_tracks.append(track)
259

    
260
  if len(upstream_tracks) > 1 :
261
    data_dict['counters']['upstream']['superfluous_track_events'] += 1
262
  if len(downstream_tracks) > 1 :
263
    data_dict['counters']['downstream']['superfluous_track_events'] += 1
264

    
265
  if len(upstream_tracks) == 1 :
266
    upstream_track = upstream_tracks[0]
267
    data_dict['counters']['upstream']['found_tracks'] += 1
268
  else :
269
    upstream_track = None
270

    
271
  if len(downstream_tracks) == 1 :
272
    downstream_track = downstream_tracks[0]
273
    data_dict['counters']['downstream']['found_tracks'] += 1
274
  else :
275
    downstream_track = None
276

    
277
  return upstream_track, downstream_track
278

    
279

    
280
def make_scifi_mc_pairs(plot_dict, data_dict, virtual_plane_dict, scif_event, mc_event) :
281
  """
282
    Make pairs of SciFiTrackpoints and MC VirtualHits
283
  """
284
  paired_hits = []
285
  paired_seeds = []
286

    
287
  expected_up, expected_down = get_expected_tracks(mc_event, virtual_plane_dict)
288
  found_up, found_down = get_found_tracks(scifi_event, plot_dict, data_dict)
289

    
290
  downstream_pt = 0.0
291
  downstream_pz = 0.0
292

    
293
  data_dict['counters']['number_events'] += 1
294

    
295
  for tracker_num, tracker, scifi_track, virtual_track in \
296
                            [ (0, "upstream", found_up, expected_up), \
297
                              (1, "downstream", found_down, expected_down) ] :
298
    if virtual_track is None :
299
      continue
300
    if ANALYSE_REFITS and (scifi_track is None or scifi_track.GetWasRefit() == 0) :
301
      continue
302

    
303
    ref_plane = tools.calculate_plane_id(tracker_num,
304
                                                    RECON_STATION, RECON_PLANE)
305
    seed_plane = tools.calculate_plane_id(tracker_num,
306
                                                    SEED_STATION, SEED_PLANE)
307
    virtual_pt = 0.0
308
    virtual_pz = 0.0
309
    virtual_radius = 0.0
310
    virtual_hits = 0
311
    scifi_hits = 0
312

    
313
    seed_virt = None
314
    reference_virt = None
315
    reference_scifi = None
316
    virtual_spacepoints = set()
317

    
318
    for plane in virtual_track :
319
      if virtual_track[plane] is not None :
320
        hit = virtual_track[plane]
321
        px = hit.GetMomentum().x()
322
        py = hit.GetMomentum().y()
323
        pt = math.sqrt( px**2 + py**2)
324
        virtual_pt += pt
325
        virtual_pz += hit.GetMomentum().z()
326

    
327
        virtual_hits += 1
328
        if plane == ref_plane :
329
          reference_virt = virtual_track[plane]
330
        if plane == seed_plane :
331
          seed_virt = virtual_track[plane]
332

    
333
        virtual_spacepoints.add(int((abs(plane)-1) / 3))
334

    
335

    
336
    virtual_pt /= virtual_hits
337
    virtual_pz /= virtual_hits
338
    virtual_p = math.sqrt( virtual_pt**2 + virtual_pz**2 )
339
    virtual_radius /= virtual_hits
340

    
341
    virtual_spacepoints = len(virtual_spacepoints)
342

    
343
    if virtual_p > P_MAX or virtual_p < P_MIN :
344
      data_dict['counters'][tracker]['momentum_cut'] += 1
345
      continue
346
    else :
347
      data_dict['counters'][tracker]['number_virtual'] += 1
348

    
349
    if scifi_track is None :
350
      plot_dict[tracker]['track_efficiency'].Fill(False, virtual_pz, virtual_pt)
351
      plot_dict[tracker]['track_efficiency_pt'].Fill(False, virtual_pt)
352
      plot_dict[tracker]['track_efficiency_pz'].Fill(False, virtual_pz)
353
      data_dict['counters'][tracker]['missing_tracks'] += 1
354

    
355
      continue # Can't do anything else without a scifi track
356

    
357

    
358
    scifi_spacepoints = set()
359
    for scifi_hit in scifi_track.scifitrackpoints() :
360
      if scifi_hit.has_data() :
361
        scifi_hits += 1
362
        scifi_spacepoints.add( scifi_hit.station()-1 )
363

    
364
        pl_id = analysis.tools.calculate_plane_id(scifi_hit.tracker(), scifi_hit.station(), scifi_hit.plane())
365

    
366
        if scifi_hit.station() == RECON_STATION and scifi_hit.plane() == RECON_PLANE :
367
          reference_scifi = scifi_hit
368

    
369
    scifi_spacepoints = len(scifi_spacepoints)
370

    
371

    
372
    plot_dict[tracker]['track_efficiency'].Fill(True, virtual_pz, virtual_pt)
373
    plot_dict[tracker]['track_efficiency_pt'].Fill(True, virtual_pt)
374
    plot_dict[tracker]['track_efficiency_pz'].Fill(True, virtual_pz)
375

    
376
    if scifi_hits >= virtual_hits :
377
      for i in range(virtual_hits) :
378
        plot_dict[tracker]['trackpoint_efficiency'].Fill(True, virtual_pz, virtual_pt)
379
        plot_dict[tracker]['trackpoint_efficiency_pt'].Fill(True, virtual_pt)
380
        plot_dict[tracker]['trackpoint_efficiency_pz'].Fill(True, virtual_pz)
381
    else :
382
      for i in range( virtual_hits - scifi_hits ) :
383
        plot_dict[tracker]['trackpoint_efficiency'].Fill(False, virtual_pz, virtual_pt)
384
        plot_dict[tracker]['trackpoint_efficiency_pt'].Fill(False, virtual_pt)
385
        plot_dict[tracker]['trackpoint_efficiency_pz'].Fill(False, virtual_pz)
386
      for i in range( scifi_hits ) :
387
        plot_dict[tracker]['trackpoint_efficiency'].Fill(True, virtual_pz, virtual_pt)
388
        plot_dict[tracker]['trackpoint_efficiency_pt'].Fill(True, virtual_pt)
389
        plot_dict[tracker]['trackpoint_efficiency_pz'].Fill(True, virtual_pz)
390

    
391
    if scifi_spacepoints >= virtual_spacepoints :
392
      for i in range(virtual_hits) :
393
        plot_dict[tracker]['spacepoint_efficiency'].Fill(True, virtual_pz, virtual_pt)
394
        plot_dict[tracker]['spacepoint_efficiency_pt'].Fill(True, virtual_pt)
395
        plot_dict[tracker]['spacepoint_efficiency_pz'].Fill(True, virtual_pz)
396
    else :
397
      for i in range( virtual_spacepoints - scifi_spacepoints ) :
398
        plot_dict[tracker]['spacepoint_efficiency'].Fill(False, virtual_pz, virtual_pt)
399
        plot_dict[tracker]['spacepoint_efficiency_pt'].Fill(False, virtual_pt)
400
        plot_dict[tracker]['spacepoint_efficiency_pz'].Fill(False, virtual_pz)
401
      for i in range( scifi_spacepoints ) :
402
        plot_dict[tracker]['spacepoint_efficiency'].Fill(True, virtual_pz, virtual_pt)
403
        plot_dict[tracker]['spacepoint_efficiency_pt'].Fill(True, virtual_pt)
404
        plot_dict[tracker]['spacepoint_efficiency_pz'].Fill(True, virtual_pz)
405

    
406
    spacepoint_set = set()
407
    for scifi_hit in scifi_track.scifitrackpoints() :
408
      if scifi_hit.has_data() :
409
        scifi_plane = analysis.tools.calculate_plane_id(scifi_hit.tracker(), scifi_hit.station(), scifi_hit.plane())
410

    
411
        if scifi_plane in virtual_track :
412
          virtual_hit = virtual_track[scifi_plane]
413
          vx = virtual_hit.GetPosition().X()
414
          vy = virtual_hit.GetPosition().Y()
415
          sx = scifi_hit.pos().X()
416
          sy = scifi_hit.pos().Y()
417
          delta = math.sqrt((sx-vx)**2 + (sy-vy)**2)
418
          
419
          if delta < 5.0 :
420
            plot_dict[tracker]["trackpoint_purity"].Fill(True, virtual_pz, virtual_pt)
421
            plot_dict[tracker]["trackpoint_purity_pt"].Fill(True, virtual_pt)
422
            plot_dict[tracker]["trackpoint_purity_pz"].Fill(True, virtual_pz)
423
          else :
424
            plot_dict[tracker]["trackpoint_purity"].Fill(False, virtual_pz, virtual_pt)
425
            plot_dict[tracker]["trackpoint_purity_pt"].Fill(False, virtual_pt)
426
            plot_dict[tracker]["trackpoint_purity_pz"].Fill(False, virtual_pz)
427

    
428
          if not (scifi_hit.station() in spacepoint_set) :
429
            spacepoint_set.add(scifi_hit.station())
430
            if delta < 5.0 :
431
              plot_dict[tracker]["spacepoint_purity"].Fill(True, virtual_pz, virtual_pt)
432
              plot_dict[tracker]["spacepoint_purity_pt"].Fill(True, virtual_pt)
433
              plot_dict[tracker]["spacepoint_purity_pz"].Fill(True, virtual_pz)
434
            else :
435
              plot_dict[tracker]["spacepoint_purity"].Fill(False, virtual_pz, virtual_pt)
436
              plot_dict[tracker]["spacepoint_purity_pt"].Fill(False, virtual_pt)
437
              plot_dict[tracker]["spacepoint_purity_pz"].Fill(False, virtual_pz)
438

    
439

    
440

    
441
    if reference_virt is None :
442
      data_dict['counters'][tracker]['missing_virtuals'] += 1
443

    
444
    if reference_scifi is None :
445
      data_dict['counters'][tracker]['missing_reference_hits'] += 1
446

    
447
    if reference_virt is not None and reference_scifi is not None :
448
      paired_hits.append( (reference_scifi, reference_virt) )
449
      data_dict['counters'][tracker]['found_pairs'] += 1
450

    
451
    if seed_virt is not None and scifi_track is not None :
452
      paired_seeds.append( (scifi_track, seed_virt))
453

    
454
  return paired_hits, paired_seeds
455

    
456

    
457
def analyse_plots(plot_dict, data_dict) :
458

    
459
  for tracker in [ 'upstream', 'downstream' ] :
460
    plot_dict[tracker]['trackpoint_efficiency'] = plot_dict[tracker]['trackpoint_efficiency'].CreateHistogram()
461
    plot_dict[tracker]['trackpoint_efficiency_pt'] = plot_dict[tracker]['trackpoint_efficiency_pt'].CreateGraph()
462
    plot_dict[tracker]['trackpoint_efficiency_pz'] = plot_dict[tracker]['trackpoint_efficiency_pz'].CreateGraph()
463

    
464
    plot_dict[tracker]['track_efficiency'] = plot_dict[tracker]['track_efficiency'].CreateHistogram()
465
    plot_dict[tracker]['track_efficiency_pt'] = plot_dict[tracker]['track_efficiency_pt'].CreateGraph()
466
    plot_dict[tracker]['track_efficiency_pz'] = plot_dict[tracker]['track_efficiency_pz'].CreateGraph()
467
                                                                                                   
468
    plot_dict[tracker]['trackpoint_purity'] = plot_dict[tracker]['trackpoint_purity'] .CreateHistogram()
469
    plot_dict[tracker]['trackpoint_purity_pt'] = plot_dict[tracker]['trackpoint_purity_pt'].CreateGraph()
470
    plot_dict[tracker]['trackpoint_purity_pz'] = plot_dict[tracker]['trackpoint_purity_pz'].CreateGraph()
471

    
472

    
473
    plot_dict[tracker]['spacepoint_efficiency'] = plot_dict[tracker]['spacepoint_efficiency'].CreateHistogram()
474
    plot_dict[tracker]['spacepoint_efficiency_pt'] = plot_dict[tracker]['spacepoint_efficiency_pt'].CreateGraph()
475
    plot_dict[tracker]['spacepoint_efficiency_pz'] = plot_dict[tracker]['spacepoint_efficiency_pz'].CreateGraph()
476

    
477
    plot_dict[tracker]['spacepoint_purity'] = plot_dict[tracker]['spacepoint_purity'] .CreateHistogram()
478
    plot_dict[tracker]['spacepoint_purity_pt'] = plot_dict[tracker]['spacepoint_purity_pt'].CreateGraph()
479
    plot_dict[tracker]['spacepoint_purity_pz'] = plot_dict[tracker]['spacepoint_purity_pz'].CreateGraph()
480

    
481

    
482
if __name__ == "__main__" : 
483
  ROOT.gROOT.SetBatch( True )
484
  ROOT.gErrorIgnoreLevel = ROOT.kError
485

    
486
  parser = argparse.ArgumentParser( description='An example script showing '+\
487
      'some basic data extraction and analysis routines' )
488

    
489
  parser.add_argument( 'maus_root_files', nargs='+', help='List of MAUS '+\
490
                  'output root files containing reconstructed straight tracks')
491

    
492
  parser.add_argument( '-N', '--max_num_events', type=int, \
493
                                   help='Maximum number of events to analyse.')
494

    
495
  parser.add_argument( '-O', '--output_filename', \
496
            default='tracker_purity_plots', help='Set the output filename')
497

    
498
  parser.add_argument( '-D', '--output_directory', \
499
                                 default='./', help='Set the output directory')
500

    
501
  parser.add_argument( '-V', '--virtual_plane_dictionary', default=None, \
502
                   help='Specify a json file containing a dictionary of the '+\
503
                                                       'virtual plane lookup' )
504

    
505
  parser.add_argument( '-P', '--print_plots', action='store_true', \
506
                        help="Flag to save the plots as individual pdf files" )
507

    
508
  parser.add_argument( '--cut_p_value', type=float, default=0.0, \
509
  help="Specify the P-Value below which tracks are removed from the analysis" )
510

    
511
  parser.add_argument( '--track_algorithm', type=int, default=1, \
512
                          help="Specify the track reconstruction algorithm. "+\
513
                             "1 for Helical Tracks and 0 for Straight Tracks" )
514

    
515
  parser.add_argument( '--pz_bin', type=float, default=PZ_BIN_WIDTH, \
516
             help="Specify the size of the Pz bins which are used to select "+\
517
                     "particles for the reconstruction of optical functions." )
518
  parser.add_argument( '--pz_window', type=float, nargs=2, \
519
        default=[PZ_MIN, PZ_MAX], help="Specify the range of Pz to consider "+\
520
                               "for the reconstruction of optical functions." )
521

    
522
  parser.add_argument( '--pt_bin', type=float, default=PT_BIN_WIDTH, \
523
             help="Specify the size of the Pt bins which are used to select "+\
524
                     "particles for the reconstruction of optical functions." )
525
  parser.add_argument( '--pt_window', type=float, nargs=2, \
526
        default=[PT_MIN, PT_MAX], help="Specify the range of Pt to consider "+\
527
                               "for the reconstruction of optical functions." )
528

    
529
  parser.add_argument( '--trackers', type=int, default=RECON_TRACKERS, \
530
                          nargs='+', help="Specifies the trackers to analyse" )
531

    
532
  parser.add_argument( '--p_window', type=float, nargs=2, \
533
             default=[P_MIN, P_MAX], help="Specify the range of the total " + \
534
                                         "momentum to consider for analysis." )
535

    
536
  parser.add_argument( '--selection_file', default=None, \
537
                 help='Name of a JSON file containing the events to analyses' )
538

    
539

    
540
  parser.add_argument( '--analyse_refits', action='store_true', \
541
                       help='Only reconstruct the tracks flagged as refitted' )
542

    
543

    
544
  try :
545
    namespace = parser.parse_args()
546

    
547
    P_VALUE_CUT = namespace.cut_p_value
548
    TRACK_ALGORITHM = namespace.track_algorithm
549

    
550
    RECON_TRACKERS = namespace.trackers
551

    
552
    ANALYSE_REFITS = namespace.analyse_refits
553

    
554
    P_MIN = namespace.p_window[0]
555
    P_MAX = namespace.p_window[1]
556

    
557
    PZ_MIN = namespace.pz_window[0]
558
    PZ_MAX = namespace.pz_window[1]
559
    PZ_BIN_WIDTH = namespace.pz_bin
560

    
561
    PT_MIN = namespace.pt_window[0]
562
    PT_MAX = namespace.pt_window[1]
563
    PT_BIN_WIDTH = namespace.pt_bin
564

    
565
    if namespace.selection_file is not None :
566
      SELECT_EVENTS = True
567
      with open(namespace.selection_file, 'r') as infile :
568
        GOOD_EVENTS = json.load(infile)
569
    else :
570
      SELECT_EVENTS = False
571

    
572
    if namespace.virtual_plane_dictionary is not None :
573
      VIRTUAL_PLANE_DICT = analysis.tools.load_virtual_plane_dict( \
574
                                           namespace.virtual_plane_dictionary )
575
  except BaseException as ex:
576
    raise
577
  else :
578
##### 1. Load MAUS globals and geometry. - NOT NECESSARY AT PRESENT
579
    # geom = load_tracker_geometry(namespace.configuration_file)
580

    
581
##### 2. Intialise plots ######################################################
582
    print
583
    sys.stdout.write( "\n- Initialising Plots : Running\r" )
584
    sys.stdout.flush()
585
    plot_dict, data_dict = init_plots_data()
586
    sys.stdout.write(   "- Initialising Plots : Done   \n" )
587

    
588
    file_reader = event_loader.maus_reader(namespace.maus_root_files)
589
    file_reader.set_max_num_events(10000)
590

    
591
##### 3. Initialise Plane Dictionary ##########################################
592
    if VIRTUAL_PLANE_DICT is None :
593
      sys.stdout.write( "\n- Finding Virtual Planes : Running\r" )
594
      sys.stdout.flush()
595
      virtual_plane_dictionary = create_virtual_plane_dict(file_reader)
596
      VIRTUAL_PLANE_DICT = virtual_plane_dictionary
597
      sys.stdout.write(   "- Finding Virtual Planes : Done   \n" )
598

    
599
    INVERSE_PLANE_DICT = inverse_virtual_plane_dict(VIRTUAL_PLANE_DICT)
600
    file_reader.select_events(GOOD_EVENTS)
601
    file_reader.set_max_num_events(namespace.max_num_events)
602
    file_reader.set_print_progress('spill')
603

    
604
##### 4. Load Events ##########################################################
605
    print "\n- Loading Spills...\n"
606
    try :
607
      while file_reader.next_selected_event() :
608

    
609
        try :
610
          scifi_event = file_reader.get_event( 'scifi' )
611
          mc_event = file_reader.get_event( 'mc' )
612

    
613
##### 5. Extract tracks and Fill Plots ########################################
614

    
615
          paired_hits, seed_pairs = make_scifi_mc_pairs(plot_dict, data_dict, \
616
                                     VIRTUAL_PLANE_DICT, scifi_event, mc_event)
617
#          fill_plots(plot_dict, data_dict, paired_hits)
618
#          fill_plots_seeds(plot_dict, data_dict, seed_pairs)
619

    
620
        except ValueError as ex :
621
          print "An Error Occured: " + str(ex)
622
          print "Skipping Event: " +\
623
                str(file_reader.get_current_event_number()) + " In Spill: " + \
624
                str(file_reader.get_current_spill_number()) + " In File: " + \
625
                str(file_reader.get_current_filenumber()) + "\n"
626
          continue
627

    
628
    except KeyboardInterrupt :
629
      print
630
      print " ###  Keyboard Interrupt  ###"
631
      print
632
    print "- {0:0.0f} Spills Loaded                                 ".format( \
633
                                            file_reader.get_total_num_spills())
634
##### 6. Analysing Plots ######################################################
635
    print"\n- Analysing Data...\n"
636

    
637
    analyse_plots(plot_dict, data_dict)
638

    
639
##### 7. Saving Plots and Data ################################################
640

    
641
    sys.stdout.write( "\n- Saving Plots and Data : Running\r" )
642
    sys.stdout.flush()
643
#    save_pretty(plot_dict, namespace.output_directory )
644

    
645
#    save_plots(plot_dict, namespace.output_directory, \
646
#                              namespace.output_filename, namespace.print_plots)
647
    filename = os.path.join(namespace.output_directory, \
648
                                                     namespace.output_filename)
649
    analysis.tools.save_plots(plot_dict, filename+'.root')
650
    sys.stdout.write(   "- Saving Plots and Data : Done   \n" )
651

    
652
  print 
653
  print "Complete."
654
  print
655

    
(2-2/2)