Functions | Variables
dumpTree Namespace Reference

Functions

def printPrimaryParticle (depth, primaryParticle)
 
def printPrimaryVertex (depth, primaryVertex)
 
def printTrajectoryPoint (depth, trajectoryPoint)
 
def printTrajectory (depth, trajectory)
 
def printHitSegment (depth, hitSegment)
 
def printSegmentContainer (depth, containerName, hitSegments)
 
def dump (input_file, output_file)
 
def loop (evt, tgeo, tout)
 

Variables

 fileName = sys.argv[1]
 
 tFile = ROOT.TFile(fileName,"OLD")
 
 events = tFile.Get("EDepSimEvents")
 
 event = ROOT.TG4Event()
 
 hit = det.second[0]
 
 parser = OptionParser()
 
 help
 
 default
 
 type
 
 args
 
 dummy
 
 fout = ROOT.TFile( args.outfile, "RECREATE" )
 
 tout = ROOT.TTree( "tree","tree" )
 
 t_ievt = array('i',[0])
 
 t_Ev = array('f', [0.])
 
 t_p3lep = array('f',3*[0.0])
 
 t_vtx = array('f',3*[0.0])
 
 t_lepDeath = array('f',3*[0.0])
 
 t_lepPdg = array('i',[0])
 
 t_lepKE = array('f',[0])
 
 t_muonExitPt = array('f',3*[0.0])
 
 t_muonExitMom = array('f',3*[0.0])
 
 t_muonReco = array('i',[0])
 
 t_muon_endVolName = ROOT.std.string()
 
 t_muGArLen = array('f',[0])
 
 t_muECalLen = array('f',[0])
 
 t_hadTot = array('f', [0.] )
 
 t_hadP = array('f', [0.] )
 
 t_hadN = array('f', [0.] )
 
 t_hadPip = array('f', [0.] )
 
 t_hadPim = array('f', [0.] )
 
 t_hadPi0 = array('f', [0.] )
 
 t_hadOther = array('f', [0.] )
 
 t_hadCollar = array('f', [0.] )
 
 t_nFS = array('i',[0])
 
 t_fsPdg = array('i',100*[0])
 
 t_fsPx = array('f',100*[0.])
 
 t_fsPy = array('f',100*[0.])
 
 t_fsPz = array('f',100*[0.])
 
 t_fsE = array('f',100*[0.])
 
 t_fsTrkLen = array('f',100*[0.])
 
 t_fsTrkFrontdEdX = array('f',100*[0.])
 
 t_fsTrkEnddEdX = array('f',100*[0.])
 
 t_fsTrkEndpointBall = array('f', 100*[0.])
 
 t_fsGamma1 = array('f',100*[0.])
 
 t_fsGamma2 = array('f',100*[0.])
 
 t_fsTrkCalo = array('f',100*[0.])
 
 t_geoEffThrowResults = ROOT.std.vector('std::vector< std::vector < uint64_t > >')()
 
 tGeoEfficiencyThrowsOut = ROOT.TTree( "geoEffThrows","geoEffThrows")
 
 t_geoEffThrowsY = ROOT.std.vector('float')()
 
 t_geoEffThrowsZ = ROOT.std.vector('float')()
 
 t_geoEffThrowsPhi = ROOT.std.vector('float')()
 
 tf = ROOT.TFile( args.infile )
 
 tgeo = tf.Get("EDepSimGeometry")
 

Function Documentation

def dumpTree.dump (   input_file,
  output_file 
)

Definition at line 102 of file dumpTree.py.

102 def dump(input_file, output_file):
103 
104  # The input file is generated in a previous test (100TestTree.sh).
105  inputFile = TFile(input_file)
106 
107  # Get the input tree out of the file.
108  inputTree = inputFile.Get("EDepSimEvents")
109  print("Class:", inputTree.ClassName())
110 
111  # Attach a brach to the events.
112  event = TG4Event()
113  inputTree.SetBranchAddress("Event",event)
114 
115  # Read all of the events.
116  entries = inputTree.GetEntriesFast()
117 
118  segments_dtype = np.dtype([('eventID', 'u4'), ('z_end', 'f4'),
119  ('trackID', 'u4'), ('tran_diff', 'f4'),
120  ('z_start', 'f4'), ('x_end', 'f4'),
121  ('y_end', 'f4'), ('n_electrons', 'u4'),
122  ('pdgId', 'i4'), ('x_start', 'f4'),
123  ('y_start', 'f4'), ('t_start', 'f4'),
124  ('dx', 'f4'), ('long_diff', 'f4'),
125  ('pixel_plane', 'i4'), ('t_end', 'f4'),
126  ('dEdx', 'f4'), ('dE', 'f4'), ('t', 'f4'),
127  ('y', 'f4'), ('x', 'f4'), ('z', 'f4'),
128  ('n_photons','f4')])
129 
130  trajectories_dtype = np.dtype([('eventID', 'u4'), ('trackID', 'u4'),
131  ('parentID', 'i4'),
132  ('pxyz_start', 'f4', (3,)),
133  ('xyz_start', 'f4', (3,)), ('t_start', 'f4'),
134  ('pxyz_end', 'f4', (3,)),
135  ('xyz_end', 'f4', (3,)), ('t_end', 'f4'),
136  ('pdgId', 'i4'), ('start_process', 'u4'),
137  ('start_subprocess', 'u4'),
138  ('end_process', 'u4'),
139  ('end_subprocess', 'u4')])
140 
141  segments_list = []
142  trajectories_list = []
143 
144  for jentry in range(entries):
145  print(jentry)
146  nb = inputTree.GetEntry(jentry)
147  if nb <= 0:
148  continue
149 
150  print("Class: ", event.ClassName())
151  print("Event number:", event.EventId)
152 
153  # Dump the primary vertices
154  # for primaryVertex in event.Primaries:
155  # printPrimaryVertex("PP", primaryVertex)
156 
157  # Dump the trajectories
158  print("Number of trajectories ", len(event.Trajectories))
159  trajectories = np.empty(len(event.Trajectories), dtype=trajectories_dtype)
160  for iTraj, trajectory in enumerate(event.Trajectories):
161  start_pt, end_pt = trajectory.Points[0], trajectory.Points[-1]
162  trajectories[iTraj]['eventID'] = jentry
163  trajectories[iTraj]['trackID'] = trajectory.GetTrackId()
164  trajectories[iTraj]['parentID'] = trajectory.GetParentId()
165  trajectories[iTraj]['pxyz_start'] = (start_pt.GetMomentum().X(), start_pt.GetMomentum().Y(), start_pt.GetMomentum().Z())
166  trajectories[iTraj]['pxyz_end'] = (end_pt.GetMomentum().X(), end_pt.GetMomentum().Y(), end_pt.GetMomentum().Z())
167  trajectories[iTraj]['xyz_start'] = (start_pt.GetPosition().X(), start_pt.GetPosition().Y(), start_pt.GetPosition().Z())
168  trajectories[iTraj]['xyz_end'] = (end_pt.GetPosition().X(), end_pt.GetPosition().Y(), end_pt.GetPosition().Z())
169  trajectories[iTraj]['t_start'] = start_pt.GetPosition().T()
170  trajectories[iTraj]['t_end'] = end_pt.GetPosition().T()
171  trajectories[iTraj]['start_process'] = start_pt.GetProcess()
172  trajectories[iTraj]['start_subprocess'] = start_pt.GetSubprocess()
173  trajectories[iTraj]['end_process'] = end_pt.GetProcess()
174  trajectories[iTraj]['end_subprocess'] = end_pt.GetSubprocess()
175  trajectories[iTraj]['pdgId'] = trajectory.GetPDGCode()
176  trajectories_list.append(trajectories)
177 
178  # Dump the segment containers
179  print("Number of segment containers:", event.SegmentDetectors.size())
180 
181  for containerName, hitSegments in event.SegmentDetectors:
182 
183  segment = np.empty(len(hitSegments), dtype=segments_dtype)
184  for iHit, hitSegment in enumerate(hitSegments):
185  segment[iHit]['eventID'] = jentry
186  segment[iHit]['trackID'] = trajectories[hitSegment.Contrib[0]]['trackID']
187  segment[iHit]['x_start'] = hitSegment.GetStart().X() / 10
188  segment[iHit]['y_start'] = hitSegment.GetStart().Y() / 10
189  segment[iHit]['z_start'] = hitSegment.GetStart().Z() / 10
190  segment[iHit]['x_end'] = hitSegment.GetStop().X() / 10
191  segment[iHit]['y_end'] = hitSegment.GetStop().Y() / 10
192  segment[iHit]['z_end'] = hitSegment.GetStop().Z() / 10
193  segment[iHit]['dE'] = hitSegment.GetEnergyDeposit()
194  segment[iHit]['t'] = 0
195  segment[iHit]['t_start'] = 0
196  segment[iHit]['t_end'] = 0
197  xd = segment[iHit]['x_end'] - segment[iHit]['x_start']
198  yd = segment[iHit]['y_end'] - segment[iHit]['y_start']
199  zd = segment[iHit]['z_end'] - segment[iHit]['z_start']
200  dx = sqrt(xd**2 + yd**2 + zd**2)
201  segment[iHit]['dx'] = dx
202  segment[iHit]['x'] = (segment[iHit]['x_start'] + segment[iHit]['x_end']) / 2.
203  segment[iHit]['y'] = (segment[iHit]['y_start'] + segment[iHit]['y_end']) / 2.
204  segment[iHit]['z'] = (segment[iHit]['z_start'] + segment[iHit]['z_end']) / 2.
205  segment[iHit]['dEdx'] = hitSegment.GetEnergyDeposit() / dx if dx > 0 else 0
206  segment[iHit]['pdgId'] = trajectories[hitSegment.Contrib[0]]['pdgId']
207  segment[iHit]['n_electrons'] = 0
208  segment[iHit]['long_diff'] = 0
209  segment[iHit]['tran_diff'] = 0
210  segment[iHit]['pixel_plane'] = 0
211  segment[iHit]['n_photons'] = 0
212  segments_list.append(segment)
213  trajectories_list = np.concatenate(trajectories_list, axis=0)
214  segments_list = np.concatenate(segments_list, axis=0)
215 
216  with h5py.File(output_file, 'w') as f:
217  f.create_dataset("trajectories", data=trajectories_list)
218  f.create_dataset("segments", data=segments_list)
219 
220 
util::enumerate
auto enumerate(Iterables &&...iterables)
Range-for loop helper tracking the number of iteration.
Definition: enumerate.h:69
dumpTree.dump
def dump(input_file, output_file)
Definition: dumpTree.py:102
def dumpTree.loop (   evt,
  tgeo,
  tout 
)

Definition at line 14 of file dumpTree.py.

14 def loop( evt, tgeo, tout ):
15 
16  offset = [ 0., 5.5, 411. ]
17  collarLo = [ -320., -120., 30. ]
18  collarHi = [ 320., 120., 470. ]
19 
20  # Initialize geometric efficiency module.
21  geoEff = pyGeoEff.geoEff(args.seed)
22  # Multiple of 64 doesn't waste bits
23  geoEff.setNthrows(4992)
24  # Use neutrino decay position, rather than fixed neutrino direction as symmetry axis
25  geoEff.setUseFixedBeamDir(False)
26 
27  # Average neutrino decay position in beam coordinates as a function of vertex x (from Luke): Will be used to set the decay position event-by-event.
28  OffAxisPoints = array('f', [-2, 0.5, 3, 5.5, 8, 10.5, 13, 15.5, 18, 20.5, 23, 25.5, 28, 30.5])
29  meanPDPZ = array('f', [ 93.6072, 93.362, 90.346, 85.6266, 81.1443, 76.6664, 73.0865, 69.8348, 67.5822, 65.005, 62.4821, 60.8336, 59.1433, 57.7352])
30  gDecayZ = ROOT.TGraph(14, OffAxisPoints, meanPDPZ)
31 
32  # Get beam parameters
33  try :
34  GNuMIFluxTree = ET.parse(os.environ["GNUMIXML"])
35  beamLineRotation = float(GNuMIFluxTree.getroot()[0].findall('beamdir')[0].findall('rotation')[0].text)
36  posText = GNuMIFluxTree.getroot()[0].findall('beampos')[0].text
37  posText = posText.replace("(", "")
38  posText = posText.replace(")", "")
39  posText = posText.replace("=", ",")
40  posText = posText.split(",")
41  beamRefDetCoord = [ float(posText[i]) for i in range(3)]
42  detRefBeamCoord = [ float(posText[i]) for i in range(3,6)]
43  except Exception, e:
44  print(str(e))
45  print("dumpTree.py WARNING!!! Error reading GNUMIXML file. Set GNUMIXML enviornment variable to point at GNuMI flux configuration xml file. Exitting!")
46  exit(-1)
47 
48  # 30 cm veto
49  geoEff.setVetoSizes([30])
50  # 20, 30 and 40 MeV threshold
51  geoEff.setVetoEnergyThresholds([20, 30, 40])
52  # Active detector dimensions
53  geoEff.setActiveX(collarLo[0]-30, collarHi[0]+30)
54  geoEff.setActiveY(collarLo[1]-30, collarHi[1]+30)
55  geoEff.setActiveZ(collarLo[2]-30, collarHi[2]+30)
56  # Range for translation throws. Use full active volume but fix X.
57  geoEff.setRangeX(-1, -1)
58  geoEff.setRandomizeX(False)
59  geoEff.setRangeY(collarLo[1]-30, collarHi[1]+30)
60  geoEff.setRangeZ(collarLo[2]-30, collarHi[2]+30)
61  # Set offset between MC coordinate system and volumes defined above.
62  geoEff.setOffsetX(offset[0])
63  geoEff.setOffsetY(offset[1])
64  geoEff.setOffsetZ(offset[2])
65 
66  event = ROOT.TG4Event()
67  events.SetBranchAddress("Event",ROOT.AddressOf(event))
68 
69  N = events.GetEntries()
70 
71  print "Starting loop over %d entries" % N
72  ient = 0 # This is unnecessary
73  iwritten = 0
74  for ient in range(N):
75 
76  if ient % 100 == 0:
77  print "Event %d of %d..." % (ient,N)
78  events.GetEntry(ient)
79 
80  for ivtx,vertex in enumerate(event.Primaries):
81 
82  ## initialize output variables
83  t_ievt[0] = ient;
84  t_vtx[0]=0.0; t_vtx[1]=0.0; t_vtx[2]=0.0;
85  t_p3lep[0]=0.0; t_p3lep[1]=0.0; t_p3lep[2]=0.0;
86  t_lepDeath[0]=0.0; t_lepDeath[1]=0.0; t_lepDeath[2]=0.0;
87  t_lepPdg[0] = 0
88  t_lepKE[0] = 0.
89  t_muonExitPt[0] = 0.0; t_muonExitPt[1] = 0.0; t_muonExitPt[2] = 0.0;
90  t_muonExitMom[0] = 0.0; t_muonExitMom[1] = 0.0; t_muonExitMom[2] = 0.0;
91  t_muonReco[0] = -1;
92  t_muon_endVolName.replace(0, ROOT.std.string.npos, "")
93  t_muGArLen[0]=0.0;
94  t_hadTot[0] = 0.
95  t_hadP[0] = 0.
96  t_hadN[0] = 0.
97  t_hadPip[0] = 0.
98  t_hadPim[0] = 0.
99  t_hadPi0[0] = 0.
100  t_hadOther[0] = 0.
101  t_hadCollar[0] = 0.
102  t_nFS[0] = 0
103  ## done
104 
105  # now ID numucc
106  reaction=vertex.Reaction
107 
108  # set the vertex location for output
109  for i in range(3):
110  t_vtx[i] = vertex.Position[i] / 10. - offset[i] # cm
111 
112  # fiducial vertex pre-cut
113  #if abs(t_vtx[0]) > 310. or abs(t_vtx[1]) > 110. or t_vtx[2] < 40. or t_vtx[2] > 360.:
114  # continue
115 
116  # Use vertex to determine mean decay point
117  decayZbeamCoord = gDecayZ.Eval(vertex.Position[0] / 1000 - detRefBeamCoord[0])*100 # in cm
118  decayZdetCoord = (-1*detRefBeamCoord[2]*100+decayZbeamCoord)*cos(beamLineRotation) - (-1*detRefBeamCoord[1]*100*sin(beamLineRotation)) + beamRefDetCoord[2]*100
119  decayYdetCoord = (-1*detRefBeamCoord[1]*100*cos(beamLineRotation)) + (-1*detRefBeamCoord[2]*100+decayZbeamCoord)*sin(beamLineRotation) + beamRefDetCoord[1]*100
120  decayXdetCoord = -1*detRefBeamCoord[0]*100 + beamRefDetCoord[0]*100
121 
122  geoEff.setDecayPos(decayXdetCoord, decayYdetCoord, decayZdetCoord)
123  geoEff.setVertex(vertex.Position[0] / 10., vertex.Position[1] / 10.,vertex.Position[2] / 10.)
124 
125  # Renew throws every 100th event written to the output file.
126  if (iwritten % 100) == 0 :
127  geoEff.throwTransforms()
128  t_geoEffThrowsY.clear()
129  for i in geoEff.getCurrentThrowTranslationsY() :
130  t_geoEffThrowsY.push_back(i)
131  t_geoEffThrowsZ.clear()
132  for i in geoEff.getCurrentThrowTranslationsZ() :
133  t_geoEffThrowsZ.push_back(i)
134  t_geoEffThrowsPhi.clear()
135  for i in geoEff.getCurrentThrowRotations() :
136  t_geoEffThrowsPhi.push_back(i)
137  tGeoEfficiencyThrowsOut.Fill()
138 
139  ileptraj = -1
140  nfsp = 0
141  nHadrons = 0
142  # get the lepton kinematics from the edepsim file
143  fsParticleIdx = {}
144  for ipart,particle in enumerate(vertex.Particles):
145  e = particle.Momentum[3]
146  p = (particle.Momentum[0]**2 + particle.Momentum[1]**2 + particle.Momentum[2]**2)**0.5
147  m = (e**2 - p**2)**0.5
148  t_fsPdg[nfsp] = particle.PDGCode
149  t_fsPx[nfsp] = particle.Momentum[0]
150  t_fsPy[nfsp] = particle.Momentum[1]
151  t_fsPz[nfsp] = particle.Momentum[2]
152  t_fsE[nfsp] = e
153  fsParticleIdx[particle.TrackId] = nfsp
154  nfsp += 1
155  pdg = particle.PDGCode
156  if abs(pdg) in [11,12,13,14]:
157  ileptraj = particle.TrackId
158  t_lepPdg[0] = pdg
159  # set the muon momentum for output
160  for i in range(3): t_p3lep[i] = particle.Momentum[i]
161  t_lepKE[0] = e - m
162 
163  assert ileptraj != -1, "There isn't a lepton??"
164  t_nFS[0] = nfsp
165 
166  # If there is a muon, determine how to reconstruct its momentum and charge
167  if abs(t_lepPdg[0]) == 13:
168  leptraj = event.Trajectories[ileptraj]
169  for p in leptraj.Points:
170  pt = p.Position
171  node = tgeo.FindNode( pt.X(), pt.Y(), pt.Z() )
172  volName = node.GetName()
173  active = False
174  if "LAr" in volName or "PixelPlane" in volName or "sPlane" in volName: # in active volume, update exit points
175  t_muonExitPt[0] = pt.X() / 10. - offset[0]
176  t_muonExitPt[1] = pt.Y() / 10. - offset[1]
177  t_muonExitPt[2] = pt.Z() / 10. - offset[2]
178  t_muonExitMom[0] = p.Momentum.x()
179  t_muonExitMom[1] = p.Momentum.y()
180  t_muonExitMom[2] = p.Momentum.z()
181  else:
182  t_muonExitPt[0] = pt.X() / 10. - offset[0]
183  t_muonExitPt[1] = pt.Y() / 10. - offset[1]
184  t_muonExitPt[2] = pt.Z() / 10. - offset[2]
185  break
186 
187  endpt = leptraj.Points[-1].Position
188 
189  node = tgeo.FindNode( endpt.X(), endpt.Y(), endpt.Z() )
190 
191  t_lepDeath[0] = endpt.X()/10. - offset[0]
192  t_lepDeath[1] = endpt.Y()/10. - offset[1]
193  t_lepDeath[2] = endpt.Z()/10. - offset[2]
194 
195  endVolName = node.GetName()
196  t_muon_endVolName.replace(0, ROOT.std.string.npos, endVolName)
197  if "LArActive" in endVolName: t_muonReco[0] = 1 # contained
198  elif "ECal" in endVolName: t_muonReco[0] = 3 # ECAL stopper
199  else: t_muonReco[0] = 0 # endpoint not in active material, but might still be reconstructed by curvature if GAr length > 0
200 
201  # look for muon hits in the gas TPC
202  hits = []
203  for key in event.SegmentDetectors:
204  if key.first in ["TPC_Drift1", "TPC_Drift2"]:
205  hits += key.second
206 
207  tot_length = 0.0
208  for hit in hits:
209  if hit.PrimaryId == ileptraj: # hit is due to the muon
210  # TG4HitSegment::TrackLength includes all delta-rays, which spiral in gas TPC and give ridiculously long tracks
211  hStart = ROOT.TVector3( hit.Start[0]/10.-offset[0], hit.Start[1]/10.-offset[1], hit.Start[2]/10.-offset[2] )
212  hStop = ROOT.TVector3( hit.Stop[0]/10.-offset[0], hit.Stop[1]/10.-offset[1], hit.Stop[2]/10.-offset[2] )
213  tot_length += (hStop-hStart).Mag()
214 
215  # look for muon hits in the ECAL
216  ehits = []
217  for key in event.SegmentDetectors:
218  if key.first in ["BarrelECal_vol", "EndcapECal_vol"]:
219  ehits += key.second
220 
221  etot_length = 0.0
222  for hit in ehits:
223  if hit.PrimaryId == ileptraj: # hit is due to the muon
224  # TG4HitSegment::TrackLength includes all delta-rays, which spiral in gas TPC and give ridiculously long tracks
225  hStart = ROOT.TVector3( hit.Start[0]/10.-offset[0], hit.Start[1]/10.-offset[1], hit.Start[2]/10.-offset[2] )
226  hStop = ROOT.TVector3( hit.Stop[0]/10.-offset[0], hit.Stop[1]/10.-offset[1], hit.Stop[2]/10.-offset[2] )
227  etot_length += (hStop-hStart).Mag()
228 
229  t_muGArLen[0] = tot_length
230  t_muECalLen[0] = etot_length
231 
232  if tot_length > 50.: t_muonReco[0] = 2 # tracker
233 
234  # hadronic containment -- find hits in ArgonCube
235  hits = []
236  for key in event.SegmentDetectors:
237  if key.first == "ArgonCube":
238  hits += key.second
239 
240  # Truth-matching energy -- make dictionary of trajectory --> primary pdg
241  traj_to_pdg = {}
242  # For pi0s, stop at the photons to determine the visible energy due to gamma1/gamma2
243  tid_to_gamma = {}
244  gamma_tids = []
245  for traj in event.Trajectories:
246  mom = traj.ParentId
247  tid = traj.TrackId
248  if event.Trajectories[mom].PDGCode == 111 and event.Trajectories[tid].PDGCode == 22 and event.Trajectories[mom].ParentId == -1:
249  gamma_tids.append(tid)
250  while mom != -1:
251  tid = mom
252  mom = event.Trajectories[mom].ParentId
253  if mom in gamma_tids:
254  tid_to_gamma[tid] = mom
255  traj_to_pdg[traj] = event.Trajectories[tid].PDGCode
256 
257  collar_energy = 0.
258  total_energy = 0.
259 
260  geoEff_EDepPosition = []
261  geoEff_EDepEnergy = []
262 
263  track_length = [0. for i in range(nfsp)]
264  dEdX = [[] for i in range(nfsp)]
265  this_step = [[0.,0.] for i in range(nfsp)]
266  end_point = [None for i in range(nfsp)]
267  int_energy = [0. for i in range(nfsp)]
268  trk_calo = [0. for i in range(nfsp)]
269  gamma_energy = {}
270  for g in gamma_tids:
271  gamma_energy[g] = 0.
272  for hit in hits:
273  hStart = ROOT.TVector3( hit.Start[0]/10.-offset[0], hit.Start[1]/10.-offset[1], hit.Start[2]/10.-offset[2] )
274  hStop = ROOT.TVector3( hit.Stop[0]/10.-offset[0], hit.Stop[1]/10.-offset[1], hit.Stop[2]/10.-offset[2] )
275 
276  # Don't use edep-sim's PrimaryId, which thinks you want to associate absoltely everything with the primary
277  # Instead, get the actual contributors (usually only one) and take the biggest
278  tid = hit.Contrib[0]
279 
280  traj = event.Trajectories[tid]
281  if traj.ParentId == -1: # primary particle
282  idx = fsParticleIdx[hit.PrimaryId]
283  trk_calo[idx] += hit.EnergyDeposit
284  end_point[idx] = hStop
285  dx = (hStop-hStart).Mag()
286  track_length[idx] += dx
287  this_step[idx][1] += dx
288  this_step[idx][0] += hit.EnergyDeposit
289  if this_step[idx][1] > 0.5:
290  dEdX[idx].append( (this_step[idx][0], this_step[idx][1], hStart) ) # MeV/cm
291  this_step[idx] = [0., 0.]
292  else: # non-primary energy
293  for k,ep in enumerate(end_point):
294  if ep is None: continue
295  if (hStart-ep).Mag() < 10.:
296  int_energy[k] += hit.EnergyDeposit
297 
298  if tid in tid_to_gamma:
299  gamma_energy[tid_to_gamma[tid]] += hit.EnergyDeposit
300 
301  if hit.PrimaryId != ileptraj: # here we do want to associate stuff to the lepton
302  hStart = ROOT.TVector3( hit.Start[0]/10.-offset[0], hit.Start[1]/10.-offset[1], hit.Start[2]/10.-offset[2] )
303  total_energy += hit.EnergyDeposit
304 
305  # check if hit is in collar region
306  if hStart.x() < collarLo[0] or hStart.x() > collarHi[0] or hStart.y() < collarLo[1] or hStart.y() > collarHi[1] or hStart.z() < collarLo[2] or hStart.z() > collarHi[2]:
307  collar_energy += hit.EnergyDeposit
308 
309  # Set up arrays for geometric efficiency
310  for dim in range(3) :
311  geoEff_EDepPosition.append((hit.Start[dim] + hit.Stop[dim])/2./10.)
312  geoEff_EDepEnergy.append(hit.EnergyDeposit)
313 
314  # Determine primary particle
315  pdg = traj_to_pdg[traj]
316  if pdg in [11, -11, 13, -13]: continue # lepton
317  elif pdg == 2212: t_hadP[0] += hit.EnergyDeposit
318  elif pdg == 2112: t_hadN[0] += hit.EnergyDeposit
319  elif pdg == 211: t_hadPip[0] += hit.EnergyDeposit
320  elif pdg == -211: t_hadPim[0] += hit.EnergyDeposit
321  elif pdg == 111: t_hadPi0[0] += hit.EnergyDeposit
322  else: t_hadOther[0] += hit.EnergyDeposit
323 
324  t_hadTot[0] = total_energy
325  t_hadCollar[0] = collar_energy
326 
327  geoEff.setHitSegEdeps(geoEff_EDepEnergy)
328  geoEff.setHitSegPoss(geoEff_EDepPosition)
329 
330  geoEffThrowResultsList = geoEff.getHadronContainmentThrows()
331 
332  t_geoEffThrowResults.clear()
333  for i in range(len(geoEffThrowResultsList)) :
334  iVec = ROOT.std.vector('vector< uint64_t >')()
335  for j in range (len(geoEffThrowResultsList[i])) :
336  jVec = ROOT.std.vector('uint64_t')()
337  for k in range(len(geoEffThrowResultsList[i][j])) :
338  jVec.push_back(geoEffThrowResultsList[i][j][k])
339  iVec.push_back(jVec)
340  t_geoEffThrowResults.push_back(iVec)
341 
342  for i in range(nfsp):
343  t_fsTrkLen[i] = track_length[i]
344  # Average of anything in the last 3cm of track
345  t_fsTrkEnddEdX[i] = 0.
346  t_fsTrkFrontdEdX[i] = 0.
347  totlen = 0.
348  j = len(dEdX[i])-1
349  while j >= 0:
350  totlen += dEdX[i][j][1]
351  t_fsTrkEnddEdX[i] += dEdX[i][j][0]
352  j -= 1
353  if totlen > 3.: break
354  if totlen > 0.: t_fsTrkEnddEdX[i] /= totlen
355 
356  totlen = 0.
357  for k in range(j+1):
358  totlen += dEdX[i][k][1]
359  t_fsTrkFrontdEdX[i] += dEdX[i][k][0]
360  if totlen > 0.: t_fsTrkFrontdEdX[i] /= totlen
361 
362  t_fsTrkEndpointBall[i] = int_energy[i]
363  t_fsTrkCalo[i] = trk_calo[i]
364 
365  t_fsGamma1[i] = 0.
366  t_fsGamma2[i] = 0.
367 
368  # Photons
369  for t in gamma_tids:
370  mom = event.Trajectories[t].ParentId
371  if t_fsGamma1[mom] == 0.: t_fsGamma1[mom] = gamma_energy[t]
372  elif t_fsGamma2[mom] == 0.: t_fsGamma2[mom] = gamma_energy[t]
373  else:
374  print "Pi0 has more than two photons wtf"
375 
376 
377  tout.Fill()
378  iwritten += 1
379  ient += 1 # This is unnecessary
380 
381 
reco_momentum_tuples.float
float
Definition: reco_momentum_tuples.py:12
util::enumerate
auto enumerate(Iterables &&...iterables)
Range-for loop helper tracking the number of iteration.
Definition: enumerate.h:69
dumpTree.loop
def loop(evt, tgeo, tout)
Definition: dumpTree.py:14
abs
T abs(T value)
Definition: sparse_vector_test.cc:30
gar::dump::array
auto array(Array const &a)
Returns a manipulator which will print the specified array.
Definition: DumpUtils.h:228
str
static QCString str
Definition: fortrancode.cpp:27098
def dumpTree.printHitSegment (   depth,
  hitSegment 
)

Definition at line 77 of file dumpTree.py.

77 def printHitSegment(depth, hitSegment):
78  print(depth,"Class: ", hitSegment.ClassName())
79  print(depth,"Primary Id:", hitSegment.GetPrimaryId())
80  print(depth,"Energy Deposit:",hitSegment.GetEnergyDeposit())
81  print(depth,"Secondary Deposit:", hitSegment.GetSecondaryDeposit())
82  print(depth,"Track Length:",hitSegment.GetTrackLength())
83  print(depth,"Start:", hitSegment.GetStart().X(),
84  hitSegment.GetStart().Y(),
85  hitSegment.GetStart().Z(),
86  hitSegment.GetStart().T())
87  print(depth,"Stop:", hitSegment.GetStop().X(),
88  hitSegment.GetStop().Y(),
89  hitSegment.GetStop().Z(),
90  hitSegment.GetStop().T())
91  print(depth,"Contributor:", [contributor for contributor in hitSegment.Contrib])
92 
93 # Print the fields in a single element of the SegmentDetectors map.
94 # The container name is the key, and the hitSegments is the value (a
95 # vector of TG4HitSegment objects).
dumpTree.printHitSegment
def printHitSegment(depth, hitSegment)
Definition: dumpTree.py:77
def dumpTree.printPrimaryParticle (   depth,
  primaryParticle 
)

Definition at line 15 of file dumpTree.py.

15 def printPrimaryParticle(depth, primaryParticle):
16  print(depth,"Class: ", primaryParticle.ClassName())
17  print(depth,"Track Id:", primaryParticle.GetTrackId())
18  print(depth,"Name:", primaryParticle.GetName())
19  print(depth,"PDG Code:",primaryParticle.GetPDGCode())
20  print(depth,"Momentum:",primaryParticle.GetMomentum().X(),
21  primaryParticle.GetMomentum().Y(),
22  primaryParticle.GetMomentum().Z(),
23  primaryParticle.GetMomentum().E(),
24  primaryParticle.GetMomentum().P(),
25  primaryParticle.GetMomentum().M())
26 
27 # Print the fields in an TG4PrimaryVertex object
gar::P
std::pair< float, std::string > P
Definition: anatree_module.cc:75
dumpTree.printPrimaryParticle
def printPrimaryParticle(depth, primaryParticle)
Definition: dumpTree.py:15
def dumpTree.printPrimaryVertex (   depth,
  primaryVertex 
)

Definition at line 28 of file dumpTree.py.

28 def printPrimaryVertex(depth, primaryVertex):
29  print(depth,"Class: ", primaryVertex.ClassName())
30  print(depth,"Position:", primaryVertex.GetPosition().X(),
31  primaryVertex.GetPosition().Y(),
32  primaryVertex.GetPosition().Z(),
33  primaryVertex.GetPosition().T())
34  print(depth,"Generator:",primaryVertex.GetGeneratorName())
35  print(depth,"Reaction:",primaryVertex.GetReaction())
36  print(depth,"Filename:",primaryVertex.GetFilename())
37  print(depth,"InteractionNumber:",primaryVertex.GetInteractionNumber())
38  depth = depth + ".."
39  for infoVertex in primaryVertex.Informational:
40  printPrimaryVertex(depth,infoVertex)
41  for primaryParticle in primaryVertex.Particles:
42  printPrimaryParticle(depth,primaryParticle)
43 
44 # Print the fields in a TG4TrajectoryPoint object
dumpTree.printPrimaryVertex
def printPrimaryVertex(depth, primaryVertex)
Definition: dumpTree.py:28
dumpTree.printPrimaryParticle
def printPrimaryParticle(depth, primaryParticle)
Definition: dumpTree.py:15
def dumpTree.printSegmentContainer (   depth,
  containerName,
  hitSegments 
)

Definition at line 96 of file dumpTree.py.

96 def printSegmentContainer(depth, containerName, hitSegments):
97  print(depth,"Detector: ", containerName, hitSegments.size())
98  depth = depth + ".."
99  for hitSegment in hitSegments: printHitSegment(depth, hitSegment)
100 
101 # Read a file and dump it.
dumpTree.printHitSegment
def printHitSegment(depth, hitSegment)
Definition: dumpTree.py:77
dumpTree.printSegmentContainer
def printSegmentContainer(depth, containerName, hitSegments)
Definition: dumpTree.py:96
def dumpTree.printTrajectory (   depth,
  trajectory 
)

Definition at line 59 of file dumpTree.py.

59 def printTrajectory(depth, trajectory):
60  print(depth,"Class: ", trajectory.ClassName())
61  depth = depth + ".."
62  print(depth,"Track Id/Parent Id:",
63  trajectory.GetTrackId(),
64  trajectory.GetParentId())
65  print(depth,"Name:",trajectory.GetName())
66  print(depth,"PDG Code",trajectory.GetPDGCode())
67  print(depth,"Initial Momentum:",trajectory.GetInitialMomentum().X(),
68  trajectory.GetInitialMomentum().Y(),
69  trajectory.GetInitialMomentum().Z(),
70  trajectory.GetInitialMomentum().E(),
71  trajectory.GetInitialMomentum().P(),
72  trajectory.GetInitialMomentum().M())
73  for trajectoryPoint in trajectory.Points:
74  printTrajectoryPoint(depth,trajectoryPoint)
75 
76 # Print the fields in a TG4HitSegment object
gar::P
std::pair< float, std::string > P
Definition: anatree_module.cc:75
dumpTree.printTrajectory
def printTrajectory(depth, trajectory)
Definition: dumpTree.py:59
dumpTree.printTrajectoryPoint
def printTrajectoryPoint(depth, trajectoryPoint)
Definition: dumpTree.py:45
def dumpTree.printTrajectoryPoint (   depth,
  trajectoryPoint 
)

Definition at line 45 of file dumpTree.py.

45 def printTrajectoryPoint(depth, trajectoryPoint):
46  print(depth,"Class: ", trajectoryPoint.ClassName())
47  print(depth,"Position:", trajectoryPoint.GetPosition().X(),
48  trajectoryPoint.GetPosition().Y(),
49  trajectoryPoint.GetPosition().Z(),
50  trajectoryPoint.GetPosition().T())
51  print(depth,"Momentum:", trajectoryPoint.GetMomentum().X(),
52  trajectoryPoint.GetMomentum().Y(),
53  trajectoryPoint.GetMomentum().Z(),
54  trajectoryPoint.GetMomentum().Mag())
55  print(depth,"Process",trajectoryPoint.GetProcess())
56  print(depth,"Subprocess",trajectoryPoint.GetSubprocess())
57 
58 # Print the fields in a TG4Trajectory object
dumpTree.printTrajectoryPoint
def printTrajectoryPoint(depth, trajectoryPoint)
Definition: dumpTree.py:45

Variable Documentation

dumpTree.args

Definition at line 391 of file dumpTree.py.

dumpTree.default

Definition at line 387 of file dumpTree.py.

dumpTree.dummy

Definition at line 391 of file dumpTree.py.

dumpTree.event = ROOT.TG4Event()

Definition at line 32 of file dumpTree.py.

dumpTree.events = tFile.Get("EDepSimEvents")

Definition at line 23 of file dumpTree.py.

dumpTree.fileName = sys.argv[1]

Definition at line 9 of file dumpTree.py.

dumpTree.fout = ROOT.TFile( args.outfile, "RECREATE" )

Definition at line 394 of file dumpTree.py.

dumpTree.help

Definition at line 387 of file dumpTree.py.

dumpTree.hit = det.second[0]

Definition at line 50 of file dumpTree.py.

dumpTree.parser = OptionParser()

Definition at line 386 of file dumpTree.py.

dumpTree.t_Ev = array('f', [0.])

Definition at line 398 of file dumpTree.py.

dumpTree.t_fsE = array('f',100*[0.])

Definition at line 448 of file dumpTree.py.

dumpTree.t_fsGamma1 = array('f',100*[0.])

Definition at line 458 of file dumpTree.py.

dumpTree.t_fsGamma2 = array('f',100*[0.])

Definition at line 460 of file dumpTree.py.

dumpTree.t_fsPdg = array('i',100*[0])

Definition at line 440 of file dumpTree.py.

dumpTree.t_fsPx = array('f',100*[0.])

Definition at line 442 of file dumpTree.py.

dumpTree.t_fsPy = array('f',100*[0.])

Definition at line 444 of file dumpTree.py.

dumpTree.t_fsPz = array('f',100*[0.])

Definition at line 446 of file dumpTree.py.

dumpTree.t_fsTrkCalo = array('f',100*[0.])

Definition at line 462 of file dumpTree.py.

dumpTree.t_fsTrkEnddEdX = array('f',100*[0.])

Definition at line 454 of file dumpTree.py.

dumpTree.t_fsTrkEndpointBall = array('f', 100*[0.])

Definition at line 456 of file dumpTree.py.

dumpTree.t_fsTrkFrontdEdX = array('f',100*[0.])

Definition at line 452 of file dumpTree.py.

dumpTree.t_fsTrkLen = array('f',100*[0.])

Definition at line 450 of file dumpTree.py.

dumpTree.t_geoEffThrowResults = ROOT.std.vector('std::vector< std::vector < uint64_t > >')()

Definition at line 466 of file dumpTree.py.

dumpTree.t_geoEffThrowsPhi = ROOT.std.vector('float')()

Definition at line 475 of file dumpTree.py.

dumpTree.t_geoEffThrowsY = ROOT.std.vector('float')()

Definition at line 471 of file dumpTree.py.

dumpTree.t_geoEffThrowsZ = ROOT.std.vector('float')()

Definition at line 473 of file dumpTree.py.

dumpTree.t_hadCollar = array('f', [0.] )

Definition at line 436 of file dumpTree.py.

dumpTree.t_hadN = array('f', [0.] )

Definition at line 426 of file dumpTree.py.

dumpTree.t_hadOther = array('f', [0.] )

Definition at line 434 of file dumpTree.py.

dumpTree.t_hadP = array('f', [0.] )

Definition at line 424 of file dumpTree.py.

dumpTree.t_hadPi0 = array('f', [0.] )

Definition at line 432 of file dumpTree.py.

dumpTree.t_hadPim = array('f', [0.] )

Definition at line 430 of file dumpTree.py.

dumpTree.t_hadPip = array('f', [0.] )

Definition at line 428 of file dumpTree.py.

dumpTree.t_hadTot = array('f', [0.] )

Definition at line 422 of file dumpTree.py.

dumpTree.t_ievt = array('i',[0])

Definition at line 396 of file dumpTree.py.

dumpTree.t_lepDeath = array('f',3*[0.0])

Definition at line 404 of file dumpTree.py.

dumpTree.t_lepKE = array('f',[0])

Definition at line 408 of file dumpTree.py.

dumpTree.t_lepPdg = array('i',[0])

Definition at line 406 of file dumpTree.py.

dumpTree.t_muECalLen = array('f',[0])

Definition at line 420 of file dumpTree.py.

dumpTree.t_muGArLen = array('f',[0])

Definition at line 418 of file dumpTree.py.

dumpTree.t_muon_endVolName = ROOT.std.string()

Definition at line 416 of file dumpTree.py.

dumpTree.t_muonExitMom = array('f',3*[0.0])

Definition at line 412 of file dumpTree.py.

dumpTree.t_muonExitPt = array('f',3*[0.0])

Definition at line 410 of file dumpTree.py.

dumpTree.t_muonReco = array('i',[0])

Definition at line 414 of file dumpTree.py.

dumpTree.t_nFS = array('i',[0])

Definition at line 438 of file dumpTree.py.

dumpTree.t_p3lep = array('f',3*[0.0])

Definition at line 400 of file dumpTree.py.

dumpTree.t_vtx = array('f',3*[0.0])

Definition at line 402 of file dumpTree.py.

dumpTree.tf = ROOT.TFile( args.infile )

Definition at line 481 of file dumpTree.py.

dumpTree.tFile = ROOT.TFile(fileName,"OLD")

Definition at line 10 of file dumpTree.py.

dumpTree.tgeo = tf.Get("EDepSimGeometry")

Definition at line 485 of file dumpTree.py.

dumpTree.tGeoEfficiencyThrowsOut = ROOT.TTree( "geoEffThrows","geoEffThrows")

Definition at line 470 of file dumpTree.py.

dumpTree.tout = ROOT.TTree( "tree","tree" )

Definition at line 395 of file dumpTree.py.

dumpTree.type

Definition at line 389 of file dumpTree.py.

Generated by   doxygen 1.8.11