Belle II Software  release-05-01-25
testCosmicT0Finder.py
1 #!/usr/bin/env python
2 # -*- coding: utf-8 -*-
3 
4 # ---------------------------------------------------------------------------------
5 # Test of performance and validation of TOPCosmicT0Finder with simulation
6 # Input file can be prepared with top/examples/simGCRData.py
7 # Output ntuple to cosmicT0FinderNtuple.root
8 # ---------------------------------------------------------------------------------
9 
10 from basf2 import *
11 from ROOT import Belle2
12 from ROOT import TH1F, TFile
13 from ROOT import gROOT, AddressOf
14 import math
15 import ROOT
16 
17 timeShift = 100 # ns (can be arbitrary chosen)
18 
19 
20 class ShiftDigits(Module):
21  ''' shift digits in time by timeShift'''
22 
23  def event(self):
24  ''' event processing '''
25 
26  digits = Belle2.PyStoreArray('TOPDigits')
27  for digit in digits:
28  digit.setTime(digit.getTime() + timeShift)
29 
30 
31 gROOT.ProcessLine('struct TreeStruct {\
32  int slot; /* slot number */ \
33  int nfot; /* number of photons in the slot */ \
34  int pdg; /* PDG code from barHit */ \
35  float x; /* local x from barHit */ \
36  float y; /* local y from barHit */ \
37  float z; /* local z from barHit */ \
38  float t; /* time from barHit */ \
39  float p; /* momentum from barHit */ \
40  float theta; /* momentum theta [deg] from barHit */ \
41  float phi; /* momentum local phi [deg] from barHit */ \
42  float dx; /* extHit - barHit difference in local x */ \
43  float dy; /* extHit - barHit difference in local y */ \
44  float dz; /* extHit - barHit difference in local z */ \
45  float dt; /* extHit - barHit difference in time */ \
46  float dp; /* extHit - barHit difference in momentum */ \
47  float dtheta; /* extHit - barHit difference in theta [deg] */ \
48  float dphi; /* extHit - barHit difference in local phi [deg] */ \
49  float t0; /* t0 determined by TOPCosmicT0Finder */ \
50  float t0err; /* error on t0 */ \
51  float timeShift; /* time shift applied to TOPDigits */ \
52  float dt0; /* difference: dt + t0 - timeShift */ \
53 };')
54 
55 from ROOT import TreeStruct
56 
57 
58 class Ntuple(Module):
59  ''' makes ntuple for MC studies of TOPCosmicT0Finder performance '''
60 
61  def initialize(self):
62  ''' initialize '''
63 
64 
65  self.file = ROOT.TFile('cosmicT0FinderNtuple.root', 'recreate')
66 
67  self.tree = ROOT.TTree('tree', '')
68 
69  self.data = TreeStruct()
70 
71  for key in TreeStruct.__dict__.keys():
72  if '__' not in key:
73  formstring = '/F'
74  if isinstance(self.data.__getattribute__(key), int):
75  formstring = '/I'
76  self.tree.Branch(key, AddressOf(self.data, key), key + formstring)
77 
78 
79  self.reso = TH1F('reso', 'T0 residuals', 200, -1.0, 1.0)
80  self.reso.SetXTitle('#Delta t0 [ns]')
81 
82  self.pull = TH1F('pull', 'T0 pulls', 200, -10.0, 10.0)
83 
84  def event(self):
85  ''' event processing: fill ntuple and histograms '''
86 
87  mcParticles = Belle2.PyStoreArray('MCParticles')
88  timeZeros = Belle2.PyStoreArray('TOPTimeZeros')
89  for timeZero in timeZeros:
90  extHit = timeZero.getRelated('ExtHits')
91  if not extHit:
92  B2WARNING('no related extHits')
93  continue
94  barHit = None
95  for barhit in mcParticles[0].getRelationsWith('TOPBarHits'):
96  if barhit.getModuleID() == extHit.getCopyID():
97  barHit = barhit
98  if not barHit:
99  B2WARNING('no corresponding TOPBarHit')
100  continue
101 
102  moduleID = timeZero.getModuleID()
103  if extHit.getCopyID() != moduleID:
104  B2ERROR('moduleID in extHit differs: ' + str(extHit.getCopyID()) +
105  ' ' + str(moduleID))
106  if barHit.getModuleID() != moduleID:
107  B2ERROR('moduleID in barHit differs: ' + str(barHit.getModuleID()) +
108  ' ' + str(moduleID))
109 
110  self.data.slot = moduleID
111  self.data.nfot = timeZero.getNumPhotons()
112  self.data.pdg = barHit.getPDG()
113  pos = barHit.getLocalPosition()
114  self.data.x = pos.x()
115  self.data.y = pos.y()
116  self.data.z = pos.z()
117  self.data.t = barHit.getTime()
118  phiBar = math.pi / 2 - math.pi / 8 * (moduleID - 0.5)
119  mom = barHit.getMomentum()
120  mom.RotateZ(phiBar)
121  self.data.p = mom.Mag()
122  self.data.theta = math.degrees(mom.Theta())
123  self.data.phi = math.degrees(mom.Phi())
124  dr = extHit.getPosition() - barHit.getPosition()
125  dr.RotateZ(phiBar)
126  self.data.dx = dr.x()
127  self.data.dy = dr.y()
128  self.data.dz = dr.z()
129  self.data.dt = extHit.getTOF() - self.data.t
130  momExt = extHit.getMomentum()
131  momExt.RotateZ(phiBar)
132  self.data.dp = momExt.Mag() - self.data.p
133  self.data.dtheta = math.degrees(momExt.Theta()) - self.data.theta
134  self.data.dphi = math.degrees(momExt.Phi()) - self.data.phi
135  self.data.t0 = timeZero.getTime()
136  self.data.t0err = timeZero.getError()
137  self.data.timeShift = timeShift
138  self.data.dt0 = self.data.dt + self.data.t0 - self.data.timeShift
139  self.file.cd()
140  self.tree.Fill()
141  self.reso.Fill(self.data.dt0)
142  self.pull.Fill(self.data.dt0/self.data.t0err)
143  if timeZero.getChi2().GetEntries() > 0:
144  timeZero.getChi2().Write()
145  timeZero.getPDF().Write()
146  timeZero.getHits().Write()
147 
148  def terminate(self):
149  ''' terminate: write and close root file '''
150 
151  self.file.cd()
152  self.file.Write()
153  self.file.Close()
154 
155 
156 # Suppress messages and warnings during processing:
157 set_log_level(LogLevel.WARNING)
158 
159 # Define a global tag (note: the one given bellow will become out-dated!)
160 use_central_database('data_reprocessing_proc8')
161 
162 # Create path
163 main = create_path()
164 
165 # input
166 roinput = register_module('RootInput')
167 main.add_module(roinput)
168 
169 # Initialize TOP geometry parameters (creation of Geant geometry is not needed)
170 main.add_module('TOPGeometryParInitializer')
171 
172 # shift digits in time
173 main.add_module(ShiftDigits())
174 
175 # t0 finder
176 finder = register_module('TOPCosmicT0Finder')
177 finder.param('useIncomingTrack', False) # or True
178 # finder.param('saveHistograms', True)
179 main.add_module(finder)
180 
181 # ntuple
182 main.add_module(Ntuple())
183 
184 # Print progress
185 main.add_module('Progress')
186 
187 # Process events
188 process(main)
189 
190 # Print statistics
191 print(statistics)
testCosmicT0Finder.Ntuple.file
file
root file
Definition: testCosmicT0Finder.py:65
testCosmicT0Finder.Ntuple
Definition: testCosmicT0Finder.py:58
testCosmicT0Finder.Ntuple.terminate
def terminate(self)
Definition: testCosmicT0Finder.py:148
testCosmicT0Finder.ShiftDigits
Definition: testCosmicT0Finder.py:20
testCosmicT0Finder.Ntuple.event
def event(self)
Definition: testCosmicT0Finder.py:84
testCosmicT0Finder.Ntuple.initialize
def initialize(self)
Definition: testCosmicT0Finder.py:61
testCosmicT0Finder.Ntuple.reso
reso
histogram of t0 residuals
Definition: testCosmicT0Finder.py:79
testCosmicT0Finder.ShiftDigits.event
def event(self)
Definition: testCosmicT0Finder.py:23
Belle2::PyStoreArray
a (simplified) python wrapper for StoreArray.
Definition: PyStoreArray.h:58
testCosmicT0Finder.Ntuple.data
data
data structure to be written to tree
Definition: testCosmicT0Finder.py:69
testCosmicT0Finder.Ntuple.tree
tree
root tree
Definition: testCosmicT0Finder.py:67
testCosmicT0Finder.Ntuple.pull
pull
histogram of t0 pulls
Definition: testCosmicT0Finder.py:82