FairShip
Loading...
Searching...
No Matches
convertToACTS.py
Go to the documentation of this file.
1# SPDX-License-Identifier: LGPL-3.0-or-later
2# SPDX-FileCopyrightText: Copyright CERN for the benefit of the SHiP Collaboration
3
4import array
5from collections import defaultdict
6
7import global_variables
8import numpy as np
9import ROOT
10import shipunit as u
11
12PDG = ROOT.TDatabasePDG.Instance()
13import acts
14import shipDet_conf
15from ShipGeoConfig import load_from_root_file
16
17
18def main():
19 inputFile = global_variables.inputFile
20 sTree = ROOT.TChain("cbmsim")
21 sTree.Add(inputFile)
22
23 geoFile = global_variables.geoFile
24 fgeo = ROOT.TFile.Open(geoFile, "READ")
25
26 outFile = global_variables.inputFile.replace(".root", "_ACTS.root")
27
28 # load Shipgeo dictionary
29 ShipGeo = load_from_root_file(fgeo, "ShipGeo")
30 run = ROOT.FairRunSim()
31 run.SetName("TGeant4") # Transport engine
32 run.SetSink(ROOT.FairRootFileSink(ROOT.TMemFile("output", "recreate"))) # Dummy output file
33 run.SetUserConfig("g4Config_basic.C") # geant4 transport not used, only needed for creating VMC field
34 run.GetRuntimeDb()
35 shipDet_conf.configure(run, ShipGeo)
36
37 with ROOT.TFile.Open(outFile, "RECREATE") as f:
38 # Trees to hold converted data: particles, vertices, SiliconTarget hits, Strawtube hits.
39 particleTree = ROOT.TTree("particles", "particles")
40 vertexTree = ROOT.TTree("vertices", "vertices")
41 siHitTree = ROOT.TTree("siHits", "siHits")
42 mtcHitTree = ROOT.TTree("mtcHits", "mtcHits")
43 strawHitTree = ROOT.TTree("strawHits", "strawHits")
44
45 # Silicon target hit branches
46 # Not all branches are used by ROOTSimHitReader, but must exist nevertheless
47 event_id = array.array("i", [0])
48 siHitTree.Branch("event_id", event_id, "event_id/i")
49 volume_id = array.array("i", [0])
50 siHitTree.Branch("volume_id", volume_id, "volume_id/i")
51 boundary_id = array.array("i", [0])
52 siHitTree.Branch("boundary_id", boundary_id, "boundary_id/i")
53 layer_id = array.array("i", [0])
54 siHitTree.Branch("layer_id", layer_id, "layer_id/i")
55 approach_id = array.array("i", [0])
56 siHitTree.Branch("approach_id", approach_id, "approach_id/i")
57 sensitive_id = array.array("i", [0])
58 siHitTree.Branch("sensitive_id", sensitive_id, "sensitive_id/i")
59 geometry_id = array.array("l", [0])
60 siHitTree.Branch("geometry_id", geometry_id, "geometry_id/l")
61 particle_id = array.array("l", [0])
62 siHitTree.Branch("particle_id", particle_id, "particle_id/l")
63 index_si = array.array("i", [0])
64 siHitTree.Branch("index", index_si, "index/I")
65 tx = array.array("f", [0])
66 siHitTree.Branch("tx", tx, "tx/F")
67 ty = array.array("f", [0])
68 siHitTree.Branch("ty", ty, "ty/F")
69 tz = array.array("f", [0])
70 siHitTree.Branch("tz", tz, "tz/F")
71 tt = array.array("f", [0])
72 siHitTree.Branch("tt", tt, "tt/F")
73 tpx = array.array("f", [0])
74 siHitTree.Branch("tpx", tpx, "tpx/F")
75 tpy = array.array("f", [0])
76 siHitTree.Branch("tpy", tpy, "tpy/F")
77 tpz = array.array("f", [0])
78 siHitTree.Branch("tpz", tpz, "tpz/F")
79 tenergy = array.array("f", [0])
80 siHitTree.Branch("te", tenergy, "te/F")
81 deltapx = array.array("f", [0])
82 siHitTree.Branch("deltapx", deltapx, "deltapx/F")
83 deltapy = array.array("f", [0])
84 siHitTree.Branch("deltapy", deltapy, "deltapy/F")
85 deltapz = array.array("f", [0])
86 siHitTree.Branch("deltapz", deltapz, "deltapz/F")
87 deltae = array.array("f", [0])
88 siHitTree.Branch("deltae", deltae, "deltae/F")
89
90 # MTC hit branches
91 event_id_mtc = array.array("i", [0])
92 mtcHitTree.Branch("event_id", event_id_mtc, "event_id/i")
93 volume_id_mtc = array.array("i", [0])
94 mtcHitTree.Branch("volume_id", volume_id_mtc, "volume_id/i")
95 boundary_id_mtc = array.array("i", [0])
96 mtcHitTree.Branch("boundary_id", boundary_id_mtc, "boundary_id/i")
97 layer_id_mtc = array.array("i", [0])
98 mtcHitTree.Branch("layer_id", layer_id_mtc, "layer_id/i")
99 approach_id_mtc = array.array("i", [0])
100 mtcHitTree.Branch("approach_id", approach_id_mtc, "approach_id/i")
101 sensitive_id_mtc = array.array("i", [0])
102 mtcHitTree.Branch("sensitive_id", sensitive_id_mtc, "sensitive_id/i")
103 geometry_id_mtc = array.array("l", [0])
104 mtcHitTree.Branch("geometry_id", geometry_id_mtc, "geometry_id/l")
105 particle_id_mtc = array.array("l", [0])
106 mtcHitTree.Branch("particle_id", particle_id_mtc, "particle_id/l")
107 index_mtc = array.array("i", [0])
108 mtcHitTree.Branch("index", index_mtc, "index/I")
109 tx_mtc = array.array("f", [0])
110 mtcHitTree.Branch("tx", tx_mtc, "tx/F")
111 ty_mtc = array.array("f", [0])
112 mtcHitTree.Branch("ty", ty_mtc, "ty/F")
113 tz_mtc = array.array("f", [0])
114 mtcHitTree.Branch("tz", tz_mtc, "tz/F")
115 tt_mtc = array.array("f", [0])
116 mtcHitTree.Branch("tt", tt_mtc, "tt/F")
117 tpx_mtc = array.array("f", [0])
118 mtcHitTree.Branch("tpx", tpx_mtc, "tpx/F")
119 tpy_mtc = array.array("f", [0])
120 mtcHitTree.Branch("tpy", tpy_mtc, "tpy/F")
121 tpz_mtc = array.array("f", [0])
122 mtcHitTree.Branch("tpz", tpz_mtc, "tpz/F")
123 te_mtc = array.array("f", [0])
124 mtcHitTree.Branch("te", te_mtc, "te/F")
125 deltapx_mtc = array.array("f", [0])
126 mtcHitTree.Branch("deltapx", deltapx_mtc, "deltapx/F")
127 deltapy_mtc = array.array("f", [0])
128 mtcHitTree.Branch("deltapy", deltapy_mtc, "deltapy/F")
129 deltapz_mtc = array.array("f", [0])
130 mtcHitTree.Branch("deltapz", deltapz_mtc, "deltapz/F")
131 deltae_mtc = array.array("f", [0])
132 mtcHitTree.Branch("deltae", deltae_mtc, "deltae/F")
133
134 # Straw tracker hit branches
135 event_id_straw = array.array("i", [0])
136 strawHitTree.Branch("event_id", event_id_straw, "event_id/i")
137 volume_id_straw = array.array("i", [0])
138 strawHitTree.Branch("volume_id", volume_id_straw, "volume_id/i")
139 boundary_id_straw = array.array("i", [0])
140 strawHitTree.Branch("boundary_id", boundary_id_straw, "boundary_id/i")
141 layer_id_straw = array.array("i", [0])
142 strawHitTree.Branch("layer_id", layer_id_straw, "layer_id/i")
143 approach_id_straw = array.array("i", [0])
144 strawHitTree.Branch("approach_id", approach_id_straw, "approach_id/i")
145 sensitive_id_straw = array.array("i", [0])
146 strawHitTree.Branch("sensitive_id", sensitive_id_straw, "sensitive_id/i")
147 geometry_id_straw = array.array("l", [0])
148 strawHitTree.Branch("geometry_id", geometry_id_straw, "geometry_id/l")
149 particle_id_straw = array.array("l", [0])
150 strawHitTree.Branch("particle_id", particle_id_straw, "particle_id/l")
151 index_straw = array.array("i", [0])
152 strawHitTree.Branch("index", index_straw, "index/I")
153 tx_straw = array.array("f", [0])
154 strawHitTree.Branch("tx", tx_straw, "tx/F")
155 ty_straw = array.array("f", [0])
156 strawHitTree.Branch("ty", ty_straw, "ty/F")
157 tz_straw = array.array("f", [0])
158 strawHitTree.Branch("tz", tz_straw, "tz/F")
159 tt_straw = array.array("f", [0])
160 strawHitTree.Branch("tt", tt_straw, "tt/F")
161 tpx_straw = array.array("f", [0])
162 strawHitTree.Branch("tpx", tpx_straw, "tpx/F")
163 tpy_straw = array.array("f", [0])
164 strawHitTree.Branch("tpy", tpy_straw, "tpy/F")
165 tpz_straw = array.array("f", [0])
166 strawHitTree.Branch("tpz", tpz_straw, "tpz/F")
167 te_straw = array.array("f", [0])
168 strawHitTree.Branch("te", te_straw, "te/F")
169 deltapx_straw = array.array("f", [0])
170 strawHitTree.Branch("deltapx", deltapx_straw, "deltapx/F")
171 deltapy_straw = array.array("f", [0])
172 strawHitTree.Branch("deltapy", deltapy_straw, "deltapy/F")
173 deltapz_straw = array.array("f", [0])
174 strawHitTree.Branch("deltapz", deltapz_straw, "deltapz/F")
175 deltae_straw = array.array("f", [0])
176 strawHitTree.Branch("deltae", deltae_straw, "deltae/F")
177
178 # Vertex branches
179 event_id_vertex = array.array("i", [0])
180 vertexTree.Branch("event_id", event_id_vertex, "event_id/i")
181 vertex_id = ROOT.std.vector("unsigned long")([0])
182 vertexTree.Branch("vertex_id", vertex_id)
183 process_vtx = ROOT.std.vector("int")([0])
184 vertexTree.Branch("process", process_vtx)
185 vx_vtx = ROOT.std.vector("float")([0])
186 vertexTree.Branch("vx", vx_vtx)
187 vy_vtx = ROOT.std.vector("float")([0])
188 vertexTree.Branch("vy", vy_vtx)
189 vz_vtx = ROOT.std.vector("float")([0])
190 vertexTree.Branch("vz", vz_vtx)
191 vt_vtx = ROOT.std.vector("float")([0])
192 vertexTree.Branch("vt", vt_vtx)
193 outgoing_particlesVec = ROOT.std.vector("unsigned long")()
194 outgoing_particles = ROOT.std.vector(ROOT.std.vector("unsigned long"))()
195 vertexTree.Branch("outgoing_particles", outgoing_particles)
196 vertex_primary_vtx = ROOT.std.vector("int")([0])
197 vertexTree.Branch("vertex_primary", vertex_primary_vtx)
198 vertex_secondary_vtx = ROOT.std.vector("int")([0])
199 vertexTree.Branch("vertex_secondary", vertex_secondary_vtx)
200 generation_vtx = ROOT.std.vector("int")([0])
201 vertexTree.Branch("generation", generation_vtx)
202
203 # Particle branches
204 event_id_particle = array.array("i", [0])
205 particleTree.Branch("event_id", event_id_particle, "event_id/i")
206 particle_id_particle = ROOT.std.vector("unsigned long")()
207 particleTree.Branch("particle_id", particle_id_particle)
208 particle_type = ROOT.std.vector("int")()
209 particleTree.Branch("particle_type", particle_type)
210 process = ROOT.std.vector("int")()
211 particleTree.Branch("process", process)
212 vx = ROOT.std.vector("float")()
213 particleTree.Branch("vx", vx)
214 vy = ROOT.std.vector("float")()
215 particleTree.Branch("vy", vy)
216 vz = ROOT.std.vector("float")()
217 particleTree.Branch("vz", vz)
218 vt = ROOT.std.vector("float")()
219 particleTree.Branch("vt", vt)
220 pp = ROOT.std.vector("float")()
221 particleTree.Branch("p", pp)
222 px = ROOT.std.vector("float")()
223 particleTree.Branch("px", px)
224 py = ROOT.std.vector("float")()
225 particleTree.Branch("py", py)
226 pz = ROOT.std.vector("float")()
227 particleTree.Branch("pz", pz)
228 m = ROOT.std.vector("float")()
229 particleTree.Branch("m", m)
230 q = ROOT.std.vector("float")()
231 particleTree.Branch("q", q)
232 eta = ROOT.std.vector("float")()
233 particleTree.Branch("eta", eta)
234 phi = ROOT.std.vector("float")()
235 particleTree.Branch("phi", phi)
236 pt = ROOT.std.vector("float")()
237 particleTree.Branch("pt", pt)
238 vertex_primary = ROOT.std.vector("int")()
239 particleTree.Branch("vertex_primary", vertex_primary)
240 vertex_secondary = ROOT.std.vector("int")()
241 particleTree.Branch("vertex_secondary", vertex_secondary)
242 particle = ROOT.std.vector("int")()
243 particleTree.Branch("particle", particle)
244 generation_particle = ROOT.std.vector("int")()
245 particleTree.Branch("generation", generation_particle)
246 sub_particle = ROOT.std.vector("float")()
247 particleTree.Branch("sub_particle", sub_particle)
248 e_loss = ROOT.std.vector("float")()
249 particleTree.Branch("e_loss", e_loss)
250 total_x0 = ROOT.std.vector("float")()
251 particleTree.Branch("total_x0", total_x0)
252 total_l0 = ROOT.std.vector("float")()
253 particleTree.Branch("total_l0", total_l0)
254 number_of_hits = ROOT.std.vector("int")()
255 particleTree.Branch("number_of_hits", number_of_hits)
256 outcome = ROOT.std.vector("float")()
257 particleTree.Branch("outcome", outcome)
258
259 for ievent, event in enumerate(sTree):
260 motherMap = defaultdict(list)
261 motherMapVal = defaultdict(list)
262 detHitMap = defaultdict(list)
263 detHitArray = []
264 strawHitArr = []
265 trID = []
266
267 if global_variables.detector == "SiliconTarget" or global_variables.detector == "SND":
268 detHitArray.clear()
269 trID.clear()
270 detHitMap.clear()
271 pointsDict = defaultdict(list)
272 for index, point in enumerate(event.SiliconTargetPoint):
273 detID = point.GetDetectorID()
274 pointsDict[detID].append(point)
275
276 for index, detID in enumerate(pointsDict):
277 allPoints = ROOT.std.vector("SiliconTargetPoint*")()
278 trID.clear()
279 for p in pointsDict[detID]:
280 allPoints.push_back(p)
281 trID.append(p.GetTrackID())
282 detHitArray.append([ROOT.SiliconTargetHit(detID, allPoints), trID[0]])
283 for index, hit in enumerate(detHitArray):
284 trackID = hit[1]
285 detHitMap[trackID].append([hit[0].GetDetectorID(), hit[0]])
286 for i in detHitMap:
287 detHitMap[i].sort(key=lambda x: x[0])
288
289 for iTrack in detHitMap:
290 for index, hit in enumerate(detHitMap[iTrack]):
291 if iTrack < -1: # If track id == -2 barcode will fail
292 continue
293 hit = hit[1]
294 layerID = (
295 6 * hit.GetLayer() + 2 * hit.GetPlane() + 4
296 ) # Layers correspond to planes, including W planes and navigation layers in between
297 event_id[0] = ievent
298 volume_id[0] = 1
299 boundary_id[0] = 0
300 layer_id[0] = layerID
301 approach_id[0] = 0
302 sens = hit.GetModule()
303 sensitive_id[0] = sens
304 geometry_id[0] = acts.GeometryIdentifier(volume=1, layer=layerID, sensitive=sens).value
305 barVal = acts.Barcode(primaryVertex=1, secondaryVertex=0, part=iTrack, gen=0, subpart=0).value
306 particle_id[0] = barVal
307 index_si[0] = index
308
309 # SHiP coordinate system translated to reconstruction coord system
310 # z -> x
311 # y -> y
312 # x -> -z
313 tx[0] = hit.GetZ() * 10 # Units cm -> mm
314 ty[0] = hit.GetY() * 10
315 tz[0] = -hit.GetX() * 10
316 tt[0] = 0
317 tpx[0] = 0.0
318 tpy[0] = 0.0
319 tpz[0] = 0.0
320 tenergy[0] = 0.0
321 deltapx[0] = 0.0
322 deltapy[0] = 0.0
323 deltapz[0] = 0.0
324 deltae[0] = -hit.GetSignal()
325
326 siHitTree.Fill()
327
328 if global_variables.detector == "MTC" or global_variables.detector == "SND":
329 detHitArray.clear()
330 trID.clear()
331 detHitMap.clear()
332 pointsDict = defaultdict(list)
333 for index, point in enumerate(event.MTCPoint):
334 detID = point.GetDetectorID()
335 pointsDict[detID].append(point)
336
337 for index, detID in enumerate(pointsDict):
338 allPoints = ROOT.std.vector("MTCPoint*")()
339 trID.clear()
340 for p in pointsDict[detID]:
341 allPoints.push_back(p)
342 trID.append(p.GetTrackID())
343 detHitArray.append([ROOT.MTCHit(detID, allPoints), trID[0]])
344 for index, hit in enumerate(detHitArray):
345 trackID = hit[1]
346 detHitMap[trackID].append([hit[0].GetDetectorID(), hit[0]])
347 for i in detHitMap:
348 detHitMap[i].sort(key=lambda x: x[0])
349
350 for iTrack in detHitMap:
351 for index, hit in enumerate(detHitMap[iTrack]):
352 if iTrack < -1: # If track id == -2 barcode will fail
353 continue
354 hit = hit[1]
355 layerID = (
356 6 * hit.GetLayer() + 2 * hit.GetPlane() + 4
357 ) # Layers correspond to SciFi planes, including Fe planes and navigation layers in between
358 event_id[0] = ievent
359 volume_id[0] = 1
360 boundary_id[0] = 0
361 layer_id[0] = layerID
362 approach_id[0] = 0
363 sens = hit.GetChannelID() # Does this work properly?
364 sensitive_id[0] = sens
365 if global_variables.detector == "SND":
366 volumeVar = 2
367 else:
368 volumeVar = 1
369 geometry_id[0] = acts.GeometryIdentifier(volume=volumeVar, layer=layerID, sensitive=sens).value
370 barVal = acts.Barcode(primaryVertex=1, secondaryVertex=0, part=iTrack, gen=0, subpart=0).value
371 particle_id[0] = barVal
372 index_mtc[0] = index
373
374 # SHiP coordinate system translated to reconstruction coord system
375 # z -> x
376 # y -> y
377 # x -> -z
378 tx_mtc[0] = hit.GetZ() * 10 # Units cm -> mm
379 ty_mtc[0] = hit.GetY() * 10
380 tz_mtc[0] = -hit.GetX() * 10
381 tt_mtc[0] = 0
382 tpx_mtc[0] = 0.0
383 tpy_mtc[0] = 0.0
384 tpz_mtc[0] = 0.0
385 te_mtc[0] = 0.0
386 deltapx_mtc[0] = 0.0
387 deltapy_mtc[0] = 0.0
388 deltapz_mtc[0] = 0.0
389 deltae_mtc[0] = -hit.GetSignal()
390
391 mtcHitTree.Fill()
392
393 if global_variables.detector == "StrawTracker":
394 # Follow shipDigiReco method of only selecting the point with the earliest time
395 ROOT.gRandom.SetSeed(13)
396 t0 = ROOT.TRandom().Rndm() * 1 * u.microsecond
397 strawHitArr.clear()
398 earliestHitPerDet = {}
399 for index, point in enumerate(event.strawtubesPoint):
400 hit = ROOT.strawtubesHit(point, t0)
401 strawHitArr.append([hit, point, point.GetTrackID()])
402 if hit.isValid():
403 detID = hit.GetDetectorID()
404 if detID in earliestHitPerDet:
405 earliest = earliestHitPerDet[detID]
406 if strawHitArr[earliest][0].GetTDC() > hit.GetTDC():
407 strawHitArr[earliest][0].setInvalid()
408 earliestHitPerDet[detID] = index
409 else:
410 strawHitArr[index][0].setInvalid()
411 else:
412 earliestHitPerDet[detID] = index
413
414 # Loop through earliest hit in each straw
415 for index, hit in enumerate(strawHitArr):
416 if hit[0].isValid():
417
418 detHitMap[hit[2]].append([hit[0].GetDetectorID(), hit[0], hit[1]])
419
420 for i in detHitMap:
421
422 detHitMap[i].sort(key=lambda x: x[0])
423
424 for iTrack in detHitMap:
425 for index, hit in enumerate(detHitMap[iTrack]):
426 if iTrack < -1: # If track id == -2 barcode will fail
427 continue
428 point = hit[2]
429 hit = hit[1]
430 layerID = 8 * hit.GetStationNumber() + 2 * hit.GetViewNumber() - 6
431
432 event_id_straw[0] = ievent
433 volume_id_straw[0] = 1 # Will only change when we have multiple tracking volumes built in geom
434 boundary_id_straw[0] = 0
435 layer_id_straw[0] = layerID
436 approach_id_straw[0] = 0
437 # Sensitive surface depends on view and layer
438 if hit.GetViewNumber() == 0 or hit.GetViewNumber() == 3:
439 sens = hit.GetStrawNumber() + hit.GetLayerNumber() * 300
440 else:
441 sens = hit.GetStrawNumber() + hit.GetLayerNumber() * 315
442
443 sensitive_id_straw[0] = sens
444
445 geometry_id_straw[0] = acts.GeometryIdentifier(volume=1, layer=layerID, sensitive=sens).value
446 particle_id_straw[0] = acts.Barcode(
447 primaryVertex=1, secondaryVertex=0, part=iTrack, gen=0, subpart=0
448 ).value
449 index_straw[0] = index
450
451 # SHiP coordinate system translated to reconstruction coord system
452 # z -> x
453 # y -> y
454 # x -> -z
455 tx_straw[0] = point.GetZ() * 10 # Units cm -> mm
456 ty_straw[0] = point.GetY() * 10
457 tz_straw[0] = -point.GetX() * 10
458 tt_straw[0] = hit.GetTDC() # In microseconds
459 tpx_straw[0] = point.GetPz()
460 tpy_straw[0] = point.GetPy()
461 tpz_straw[0] = -point.GetPx()
462 te_straw[0] = 0.0
463 deltapx_straw[0] = 0.0
464 deltapy_straw[0] = 0.0
465 deltapz_straw[0] = 0.0
466 deltae_straw[0] = -point.GetEnergyLoss() # Energy for only a single point
467
468 strawHitTree.Fill()
469
470
471 for count, part in enumerate(event.MCTrack):
472 # He3, alpha, deutron, triton not in database, $ROOTSYS/etc/pdg_table.txt
473 if part.GetPdgCode() >= 1000010020:
474 continue
475 event_id_particle[0] = ievent
476 generation_particle.push_back(0)
477 barVal = acts.Barcode(primaryVertex=1, secondaryVertex=0, part=count, gen=0, subpart=0).value
478 particle_id_particle.push_back(barVal)
479 particle_type.push_back(part.GetPdgCode())
480 process.push_back(0) # Process not recorded, fill with zeros
481 # SHiP coordinate system translated to reconstruction coord system
482 # z -> x
483 # y -> y
484 # x -> -z
485 fourVec = ROOT.Math.LorentzVector("ROOT::Math::PxPyPzE4D<double>")(
486 part.GetPz(), part.GetPy(), -part.GetPx(), part.GetEnergy()
487 )
488 vx.push_back(part.GetStartZ() * 10) # Units mm
489 vy.push_back(part.GetStartY() * 10)
490 vz.push_back(-part.GetStartX() * 10)
491 vt.push_back(part.GetStartT())
492 px.push_back(part.GetPz())
493 py.push_back(part.GetPy())
494 pz.push_back(-part.GetPx())
495 pp.push_back(part.GetP())
496 m.push_back(part.GetMass())
497 charge = PDG.GetParticle(part.GetPdgCode()).Charge()
498 q.push_back(np.sign(charge))
499 eta.push_back(fourVec.Eta())
500 phi.push_back(fourVec.Phi())
501 pt.push_back(fourVec.Pt())
502 particle.push_back(count)
503 sub_particle.push_back(0)
504 vertex_primary.push_back(1)
505 vertex_secondary.push_back(0)
506 number_of_hits.push_back(len(detHitMap[count]))
507 total_x0.push_back(0)
508 total_l0.push_back(0)
509 outcome.push_back(0)
510
511
512 motherID = part.GetMotherId()
513 motherMap[str(motherID)].append(count) # Vector of particles grouped by mother
514 motherMapVal[str(motherID)].append(barVal) # Vector of particles grouped by mother
515
516 particleTree.Fill()
517 particle_id_particle.clear()
518 particle_type.clear()
519 process.clear()
520 vx.clear()
521 vy.clear()
522 vz.clear()
523 vt.clear()
524 px.clear()
525 py.clear()
526 pz.clear()
527 pp.clear()
528 m.clear()
529 q.clear()
530 eta.clear()
531 phi.clear()
532 pt.clear()
533 particle.clear()
534 sub_particle.clear()
535 generation_particle.clear()
536 vertex_primary.clear()
537 vertex_secondary.clear()
538 number_of_hits.clear()
539 total_x0.clear()
540 total_l0.clear()
541 outcome.clear()
542
543 for c, i in enumerate(motherMap):
544 particleCodes = motherMapVal[str(i)]
545 # If there are not at least 2 tracks at decay vtx, do not write
546 if len(particleCodes) < 2:
547 continue
548 event_id_vertex[0] = ievent
549 vertexVal = acts.Barcode(primaryVertex=1, secondaryVertex=0, part=0, gen=0, subpart=0).value
550 vertex_id.push_back(vertexVal)
551 particle_0ID = motherMap[str(i)]
552 particle_0 = event.MCTrack[particle_0ID[0]]
553 # SHiP coordinate system translated to reconstruction coord system
554 # z -> x
555 # y -> y
556 # x -> -z
557 vx_vtx.push_back(particle_0.GetStartZ() * 10)
558 vy_vtx.push_back(particle_0.GetStartY() * 10)
559 vz_vtx.push_back(-particle_0.GetStartX() * 10)
560 vt_vtx.push_back(0)
561 vertex_primary_vtx.push_back(1)
562 vertex_secondary_vtx.push_back(0)
563 generation_vtx.push_back(0)
564 process_vtx.push_back(0)
565 for j in particleCodes:
566 outgoing_particlesVec.push_back(j)
567 outgoing_particles.push_back(outgoing_particlesVec)
568
569 vertexTree.Fill()
570 outgoing_particles.clear()
571 outgoing_particlesVec.clear()
572 vertex_id.clear()
573 vx_vtx.clear()
574 vy_vtx.clear()
575 vz_vtx.clear()
576 vt_vtx.clear()
577 vertex_primary_vtx.clear()
578 vertex_secondary_vtx.clear()
579 generation_vtx.clear()
580 process_vtx.clear()
581 motherMap.clear()
582 motherMapVal.clear()
583
584 f.WriteObject(siHitTree, "siHits")
585 f.WriteObject(mtcHitTree, "mtcHits")
586 f.WriteObject(strawHitTree, "strawHits")
587 f.WriteObject(particleTree, "particles")
588 f.WriteObject(vertexTree, "vertices")
589
590
591if __name__ == "__main__":
592 ROOT.gErrorIgnoreLevel = ROOT.kWarning
593 ROOT.gROOT.SetBatch(True)
594 main()
shipDet_conf.configure
def configure(run, ship_geo)
Definition: shipDet_conf.py:265