198 def peel(self, track):
199 """Aggregate the track and MC information for dE/dx analysis"""
200 mc_track_lookup = self.mc_track_lookup
201
202
203
204
205
206
207
208
209
210
211
212 mc_particle = mc_track_lookup.getMCParticle(track)
213 pdg_code = mc_particle.getPDG()
214 t_vertex = mc_particle.getVertex()
215 t_mom = mc_particle.getMomentum()
216 charge = mc_particle.getCharge()
217 mass = mc_particle.getMass()
218 mc_energy = mc_particle.getEnergy()
219
220 mc_vertex2D = TFCDC.Vector2D(t_vertex.XYvector())
221 mc_mom2D = TFCDC.Vector2D(t_mom.XYvector())
222 mc_trajectory2D = TFCDC.CDCTrajectory2D(mc_vertex2D, 0, mc_mom2D, charge)
223 mc_pt = mc_mom2D.norm()
224
225 first_hit = self.mc_track_lookup.getFirstHit(track)
226 first_sim_hit = first_hit.getRelated("CDCSimHits")
227 if first_sim_hit is None:
228 return
229
230 if first_sim_hit.getWireID().getICLayer() != 0:
231 return
232
233 last_hit = self.mc_track_lookup.getLastHit(track)
234 last_sim_hit = last_hit.getRelated("CDCSimHits")
235 if last_sim_hit is None:
236 return
237
238
239
240
241 first_sim_pos3D = TFCDC.Vector3D(first_sim_hit.getPosTrack())
242 first_sim_mom3D = TFCDC.Vector3D(first_sim_hit.getMomentum())
243 first_sim_tof = first_sim_hit.getFlightTime()
244 first_sim_energy = np.sqrt(first_sim_mom3D.norm() ** 2 + mass ** 2)
245 first_sim_pt = first_sim_mom3D.cylindricalR()
246 first_sim_pz = first_sim_mom3D.z()
247
248 first_sim_mom2D = first_sim_mom3D.xy()
249
250
251 first_sim_trajectory2D = TFCDC.CDCTrajectory2D(first_sim_pos3D.xy(), first_sim_tof, first_sim_mom2D, charge)
252
253
254
255
256 for reco_hit3D in track:
257 sim_hit = self.mc_hit_lookup.getSimHit(reco_hit3D.getWireHit().getHit())
258 if sim_hit is None:
259 continue
260
261 sim_mom3D = TFCDC.Vector3D(sim_hit.getMomentum())
262 sim_energy = np.sqrt(sim_mom3D.norm() ** 2 + mass ** 2)
263 sim_pt = sim_mom3D.cylindricalR()
264 sim_pz = sim_mom3D.z()
265
266 mc_eloss_truth = mc_energy - sim_energy
267 first_eloss_truth = first_sim_energy - sim_energy
268
269 sim_pos3D = TFCDC.Vector3D(sim_hit.getPosTrack())
270 sim_pos2D = sim_pos3D.xy()
271
272 layer_cid = reco_hit3D.getWire().getICLayer()
273 bz = self.bfield.getBFieldZ(sim_pos3D)
274 r = sim_pos3D.cylindricalR()
275
276
277
278
279
280
281
282 mc_s2D = mc_trajectory2D.calcArcLength2D(sim_pos2D)
283 first_s2D = first_sim_trajectory2D.calcArcLength2D(sim_pos2D)
284
285 if mc_s2D < 0:
286
287 continue
288
289 mc_eloss_estimate = self.eloss_estimator.getEnergyLoss(mc_pt, pdg_code, mc_s2D)
290 first_eloss_estimate = self.eloss_estimator.getEnergyLoss(first_sim_pt, pdg_code, first_s2D)
291 first_ploss_factor = self.eloss_estimator.getMomentumLossFactor(first_sim_pt, pdg_code, first_s2D)
292
293 sasha_eloss_estimate = getBetheEnergyLoss(first_sim_pt, pdg_code, first_s2D)
294 sasha_ploss_factor = DeltaR(first_s2D, pdg_code, first_sim_pt)
295
296
297 first_residual2D = first_sim_trajectory2D.getDist2D(sim_pos2D)
298 first_disp2D = charge * first_residual2D
299
300 first_loss_disp2D_estimate = abs(self.eloss_estimator.getLossDist2D(first_sim_pt, pdg_code, first_s2D))
301 first_delossed_disp2D = first_disp2D - first_loss_disp2D_estimate
302
303
304 mc_residual2D = mc_trajectory2D.getDist2D(sim_pos2D)
305 mc_disp2D = charge * mc_residual2D
306
307 mc_loss_disp2D_estimate = abs(self.eloss_estimator.getLossDist2D(mc_pt, pdg_code, mc_s2D))
308 mc_delossed_disp2D = mc_disp2D - mc_loss_disp2D_estimate
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330 if abs(mc_residual2D) > 6:
331 continue
332
333 if abs(first_residual2D) > 6:
334 continue
335
336 yield dict(
337 layer_cid=layer_cid,
338 r=r,
339 bz=bz,
340
341 pdg_code=abs(pdg_code),
342 mass=mass,
343 charge=charge,
344 mc_pt=mc_pt,
345 first_sim_pt=first_sim_pt,
346 sim_pt=sim_pt,
347
348 diff_pt=first_sim_pt - sim_pt,
349 diff_pz=first_sim_pz - sim_pz,
350
351 mc_eloss_truth=mc_eloss_truth,
352 mc_eloss_estimate=mc_eloss_estimate,
353
354 first_eloss_truth=first_eloss_truth,
355
356 first_eloss_estimate=first_eloss_estimate,
357 sasha_eloss_estimate=sasha_eloss_estimate,
358
359 first_ploss_factor=first_ploss_factor,
360 sasha_ploss_factor=sasha_ploss_factor,
361
362
363
364 first_s2D=first_s2D,
365 first_residual2D=first_residual2D,
366 first_disp2D=first_disp2D,
367 first_loss_disp2D_estimate=first_loss_disp2D_estimate,
368 first_delossed_disp2D=first_delossed_disp2D,
369
370 mc_s2D=mc_s2D,
371 mc_residual2D=mc_residual2D,
372 mc_disp2D=mc_disp2D,
373 mc_loss_disp2D_estimate=mc_loss_disp2D_estimate,
374 mc_delossed_disp2D=mc_delossed_disp2D,
375
376
377
378
379
380
381
382 )
383