15from ROOT
import Belle2
16from ROOT.Belle2
import TrackFindingCDC
as TFCDC
33 return logging.getLogger(__name__)
36CONTACT =
"oliver.frost@desy.de"
39def getNaiveBetheEnergyLoss(p, pdg_code, s):
41 eloss = s * 1 * 10**-4
45def getBetheStoppingPower(p, pdg_code):
47 energy = np.sqrt(mass * mass + p * p)
52 gamma2 = gamma * gamma
55 eMass = 0.511 * Belle2.Unit.MeV
57 m_eDensity = 0.000515726
59 Wmax = 2 * eMass * beta2 * gamma2 / (1 + 2 * gamma * eMass / mass)
62 K = 0.307075 * Belle2.Unit.MeV * Belle2.Unit.cm2
64 K * m_eDensity / beta2 *
65 (1.0 / 2.0 * math.log(2 * eMass * beta2 * gamma2 * Wmax / I2) - beta2)
72 factor = 178.6 * 0.00015e-04 / 2
74 dEdx = factor / beta2 * (math.log(2 * eMass * beta2 * gamma2 * Wmax / I2) - beta2)
83def getBetheEnergyLoss(p, pdg_code, path):
84 dEdx = getBetheStoppingPower(p, pdg_code)
89def getMomentumLossFactor(p, pdg_code, eloss):
90 p_factor = (p - eloss) / p
94def DeltaR(path, particleID, P):
95 eloss = getBetheEnergyLoss(P, particleID, path)
96 return getMomentumLossFactor(P, particleID, eloss)
100 """Harvester to generate, postprocess and inspect MC events for energy-loss evaluation"""
104 generator_module =
"eloss_gun"
106 detector_setup =
"TrackingDetectorConstB"
110 """Get the output ROOT filename"""
114 """Harvest and post-process the generated events"""
117 path.add_module(harvesting_module)
118 return harvesting_module
121 """Convert command-line arguments to basf2 argument list"""
123 return argument_parser
125 def create_path(self):
127 Sets up a path that plays back pregenerated events or generates events
128 based on the properties
in the base
class.
130 path = super().create_path()
132 path.add_module("TFCDC_WireHitPreparer",
134 flightTimeEstimation=
"outwards",
137 path.add_module(
'TFCDC_AxialTrackCreatorMCTruth',
139 useOnlyBeforeTOP=
True,
141 reconstructedDriftLength=
True,
142 reconstructedPositions=
True)
149 """Module to collect information about the generated segments and
150 compose validation plots on terminate."""
154 super().__init__(foreach='CDCTrackVector',
155 output_file_name=output_file_name)
160 origin_track_fitter = TFCDC.CDCRiemannFitter()
161 origin_track_fitter.setOriginConstrained()
166 """Receive signal at the start of event processing"""
178 """Initialize the MC-hit lookup method"""
179 TFCDC.CDCMCHitLookUp.getInstance().fill()
182 """Select tracks with at least 4 segments and associated primary MC particle with pt >= 0.25 GeV/c"""
187 mc_particle = mc_track_lookup.getMCParticle(track)
190 if mc_particle
is None:
193 if mc_particle.getMomentum().Rho() < 0.250:
196 return is_primary(mc_particle)
199 """Aggregate the track and MC information for dE/dx analysis"""
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()
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()
226 first_sim_hit = first_hit.getRelated(
"CDCSimHits")
227 if first_sim_hit
is None:
230 if first_sim_hit.getWireID().getICLayer() != 0:
234 last_sim_hit = last_hit.getRelated(
"CDCSimHits")
235 if last_sim_hit
is None:
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()
248 first_sim_mom2D = first_sim_mom3D.xy()
251 first_sim_trajectory2D = TFCDC.CDCTrajectory2D(first_sim_pos3D.xy(), first_sim_tof, first_sim_mom2D, charge)
256 for reco_hit3D
in track:
257 sim_hit = self.
mc_hit_lookup.getSimHit(reco_hit3D.getWireHit().getHit())
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()
266 mc_eloss_truth = mc_energy - sim_energy
267 first_eloss_truth = first_sim_energy - sim_energy
269 sim_pos3D = TFCDC.Vector3D(sim_hit.getPosTrack())
270 sim_pos2D = sim_pos3D.xy()
272 layer_cid = reco_hit3D.getWire().getICLayer()
273 bz = self.
bfield.getBFieldZ(sim_pos3D)
274 r = sim_pos3D.cylindricalR()
282 mc_s2D = mc_trajectory2D.calcArcLength2D(sim_pos2D)
283 first_s2D = first_sim_trajectory2D.calcArcLength2D(sim_pos2D)
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)
293 sasha_eloss_estimate = getBetheEnergyLoss(first_sim_pt, pdg_code, first_s2D)
294 sasha_ploss_factor = DeltaR(first_s2D, pdg_code, first_sim_pt)
297 first_residual2D = first_sim_trajectory2D.getDist2D(sim_pos2D)
298 first_disp2D = charge * first_residual2D
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
304 mc_residual2D = mc_trajectory2D.getDist2D(sim_pos2D)
305 mc_disp2D = charge * mc_residual2D
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
330 if abs(mc_residual2D) > 6:
333 if abs(first_residual2D) > 6:
341 pdg_code=abs(pdg_code),
345 first_sim_pt=first_sim_pt,
348 diff_pt=first_sim_pt - sim_pt,
349 diff_pz=first_sim_pz - sim_pz,
351 mc_eloss_truth=mc_eloss_truth,
352 mc_eloss_estimate=mc_eloss_estimate,
354 first_eloss_truth=first_eloss_truth,
356 first_eloss_estimate=first_eloss_estimate,
357 sasha_eloss_estimate=sasha_eloss_estimate,
359 first_ploss_factor=first_ploss_factor,
360 sasha_ploss_factor=sasha_ploss_factor,
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,
371 mc_residual2D=mc_residual2D,
373 mc_loss_disp2D_estimate=mc_loss_disp2D_estimate,
374 mc_delossed_disp2D=mc_delossed_disp2D,
386 save_tree = refiners.save_tree()
388 save_histograms = refiners.save_histograms(
394 save_histograms_stackby_charge = refiners.save_histograms(
408 save_scatter = refiners.save_scatters(
414 groupby=[
None,
"charge"],
415 filter=
lambda x: x == 211,
416 filter_on=
"pdg_code",
420 save_profiles = refiners.save_profiles(
427 groupby=[
None,
"charge"],
431 save_bz_profiles = refiners.save_profiles(
440 save_cid_histogram = refiners.save_histograms(
445 'first_delossed_disp2D',
450 groupby=[
"layer_cid"],
455 save_cid_profiles = refiners.save_profiles(
461 'first_delossed_disp2D',
464 groupby=[
"layer_cid"],
469 save_cid_scatters = refiners.save_scatters(
475 'first_delossed_disp2D',
478 groupby=[
"layer_cid"],
485 save_energy_cid_histogram = refiners.save_histograms(
488 'first_eloss_estimate',
492 groupby=[
"layer_cid"],
494 folder_name=
'energy/{groupby_addition}',
498 save_energy_cid_profiles = refiners.save_profiles(
502 'first_eloss_estimate',
505 groupby=[
"layer_cid"],
507 folder_name=
'energy/{groupby_addition}',
511 save_energy_cid_scatters = refiners.save_scatters(
515 'first_eloss_estimate',
518 groupby=[
"layer_cid"],
520 folder_name=
'energy/{groupby_addition}',
526 run.configure_and_execute_from_commandline()
529if __name__ ==
"__main__":
530 logging.basicConfig(stream=sys.stdout, level=logging.INFO, format=
'%(levelname)s:%(message)s')
mc_track_lookup
by default, there is no method to find matching MC tracks
mc_hit_lookup
Method to find matching MC hits.
track_fitter
Use the CDCReimannFitter with a constrained origin for track fitting.
bfield
Method to interrogate the magnetic field values.
eloss_estimator
Method to estimate dE/dx in the CDC.
def __init__(self, output_file_name)
def create_argument_parser(self, **kwds)
def harvesting_module(self, path=None)
def output_file_name(self)
None output_file_name
Disable the writing of an output ROOT file.
None output_file_name
There is no default for the name of the output TFile.