| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 /// \file SteppingAction.cc
27 /// \brief Implementation of the SteppingActio
28 //
29 //
30
31 //....oooOO0OOooo........oooOO0OOooo........oo
32 //....oooOO0OOooo........oooOO0OOooo........oo
33
34 #include "SteppingAction.hh"
35
36 #include "Run.hh"
37
38 #include "G4DecayProducts.hh"
39 #include "G4DecayTable.hh"
40 #include "G4LossTableManager.hh"
41 #include "G4ParticleDefinition.hh"
42 #include "G4ParticleTypes.hh"
43 #include "G4Step.hh"
44 #include "G4StepPoint.hh"
45 #include "G4SystemOfUnits.hh"
46 #include "G4TouchableHistory.hh"
47 #include "G4Track.hh"
48 #include "G4VDecayChannel.hh"
49 #include "G4VPhysicalVolume.hh"
50 #include "G4VTouchable.hh"
51
52 //....oooOO0OOooo........oooOO0OOooo........oo
53
54 SteppingAction::SteppingAction() : G4UserStepp
55 {
56 Initialize();
57 }
58
59 //....oooOO0OOooo........oooOO0OOooo........oo
60
61 SteppingAction::~SteppingAction() {}
62
63 //....oooOO0OOooo........oooOO0OOooo........oo
64
65 void SteppingAction::Initialize()
66 {
67 // Initialization needed at the beginning of
68 fRunPtr = nullptr;
69 fToleranceEPviolations = 1.0 * CLHEP::eV; /
70 fPrimaryParticleId = 0;
71 fPrimaryParticleInitialKineticEnergy = 0.0;
72 fPrimaryParticleInitialTotalEnergy = 0.0;
73 fPrimaryParticleInitialMomentum = 0.0;
74 fPrimaryParticleInitialBeta = 1.0;
75 fPrimaryParticleInitialGamma = 1.0;
76 fPrimaryParticleInitial3Momentum = G4ThreeVe
77 fPrimaryParticleInitialPosition = G4ThreeVec
78 fMaxEkin_deltaMax = 0.0;
79 fMaxEtot_deltaMax = 0.0;
80 fMaxP_deltaMax = 0.0;
81 fMaxPdir_deltaMax = 0.0;
82 fMaxMass_deltaMax1 = 0.0;
83 fMaxMass_deltaMax2 = 0.0;
84 fMaxMass_deltaMax3 = 0.0;
85 fMeanMass_deltaMax3 = 0.0;
86 fMaxBeta_deltaMax1 = 0.0;
87 fMaxBeta_deltaMax2 = 0.0;
88 fMaxGamma_deltaMax1 = 0.0;
89 fMaxGamma_deltaMax2 = 0.0;
90 fMaxGamma_deltaMax3 = 0.0;
91 fMaxT_proper_deltaMax = 0.0;
92 fMaxT_lab_deltaMax = 0.0;
93 fMaxMc_truth_rPos_deltaMax = 0.0;
94 fMeanMc_truth_rPos_deltaMax = 0.0;
95 fMeanDeltaR_primaryDecay = 0.0;
96 fMinDeltaR_primaryDecay = 9999999.9;
97 fMaxDeltaR_primaryDecay = -9999999.9;
98 fMeanR_primaryDecay = 0.0;
99 fMinR_primaryDecay = 9999999.9;
100 fMaxR_primaryDecay = -9999999.9;
101 fMeanX_primaryDecay = 0.0;
102 fMinX_primaryDecay = 9999999.9;
103 fMaxX_primaryDecay = -9999999.9;
104 fMeanY_primaryDecay = 0.0;
105 fMinY_primaryDecay = 9999999.9;
106 fMaxY_primaryDecay = -9999999.9;
107 fMeanZ_primaryDecay = 0.0;
108 fMinZ_primaryDecay = 9999999.9;
109 fMaxZ_primaryDecay = -9999999.9;
110 fMeanDeltaAngle_primaryDecay = 0.0;
111 fMinDeltaAngle_primaryDecay = 9999999.9;
112 fMaxDeltaAngle_primaryDecay = -9999999.9;
113 fMeanDeltaEkin_primaryDecay = 0.0;
114 fMinDeltaEkin_primaryDecay = 9999999.9;
115 fMaxDeltaEkin_primaryDecay = -9999999.9;
116 fMeanEkin_primaryDecay = 0.0;
117 fMinEkin_primaryDecay = 9999999.9;
118 fMaxEkin_primaryDecay = -9999999.9;
119 fMeanPx_primaryDecay = 0.0;
120 fMinPx_primaryDecay = 9999999.9;
121 fMaxPx_primaryDecay = -9999999.9;
122 fMeanPy_primaryDecay = 0.0;
123 fMinPy_primaryDecay = 9999999.9;
124 fMaxPy_primaryDecay = -9999999.9;
125 fMeanPz_primaryDecay = 0.0;
126 fMinPz_primaryDecay = 9999999.9;
127 fMaxPz_primaryDecay = -9999999.9;
128 fMinUnderestimated_mc_truth_rPos_delta = 999
129 fMaxOverestimated_mc_truth_rPos_delta = -999
130 fMeanUnderestimated_mc_truth_rPos_delta = 0.
131 fMeanOverestimated_mc_truth_rPos_delta = 0.0
132 fMinUnderestimated_rDeltaPos = 9999999.9;
133 fMaxOverestimated_rDeltaPos = -9999999.9;
134 fMeanUnderestimated_rDeltaPos = 0.0;
135 fMeanOverestimated_rDeltaPos = 0.0;
136 fMaxFloat_rDeltaPos_deltaMax = -9999999.9;
137 fMeanViolationE_primaryDecay = 0.0;
138 fMinViolationE_primaryDecay = 9999999.9;
139 fMaxViolationE_primaryDecay = -9999999.9;
140 fMeanViolationPx_primaryDecay = 0.0;
141 fMinViolationPx_primaryDecay = 9999999.9;
142 fMaxViolationPx_primaryDecay = -9999999.9;
143 fMeanViolationPy_primaryDecay = 0.0;
144 fMinViolationPy_primaryDecay = 9999999.9;
145 fMaxViolationPy_primaryDecay = -9999999.9;
146 fMeanViolationPz_primaryDecay = 0.0;
147 fMinViolationPz_primaryDecay = 9999999.9;
148 fMaxViolationPz_primaryDecay = -9999999.9;
149 }
150
151 //....oooOO0OOooo........oooOO0OOooo........oo
152
153 void SteppingAction::UserSteppingAction(const
154 {
155 // Store the information about the ID and th
156 // at the first step of the first event.
157 // Note that for the kinetic energy, we are
158 if (theStep->GetTrack()->GetParentID() == 0
159 fPrimaryParticleId = theStep->GetTrack()->
160 fPrimaryParticleInitialKineticEnergy = the
161 fPrimaryParticleInitialTotalEnergy = theSt
162 fPrimaryParticleInitial3Momentum = theStep
163 fPrimaryParticleInitialMomentum = fPrimary
164 fPrimaryParticleInitialPosition = theStep-
165 fPrimaryParticleInitialBeta = theStep->Get
166 fPrimaryParticleInitialGamma = theStep->Ge
167 // As tolerance for EP violations, conside
168 // and 1 billionth of the initial, primary
169 if (fToleranceEPviolations < fPrimaryParti
170 fToleranceEPviolations = fPrimaryParticl
171 }
172 // Set the values of this run to the Run o
173 if (fRunPtr) {
174 fRunPtr->SetPrimaryParticleId(fPrimaryPa
175 fRunPtr->SetPrimaryParticleInitialKineti
176 fRunPtr->SetPrimaryParticleInitialTotalE
177 fRunPtr->SetPrimaryParticleInitialMoment
178 fRunPtr->SetPrimaryParticleInitialBeta(f
179 fRunPtr->SetPrimaryParticleInitialGamma(
180 fRunPtr->SetPrimaryParticleInitial3Momen
181 fRunPtr->SetPrimaryParticleInitialPositi
182 fRunPtr->SetToleranceEPviolations(Tolera
183 fRunPtr->SetToleranceDeltaDecayRadius(To
184 fRunPtr->SetIsPreassignedDecayEnabled(Is
185 fRunPtr->SetIsBoostToLabEnabled(IsBoostT
186 }
187 // Use the preassigned decay is enabled
188 if (IsPreassignedDecayEnabled() && (!theSt
189 G4DynamicParticle* dynamicParent =
190 const_cast<G4DynamicParticle*>(theStep
191 if (dynamicParent != nullptr) {
192 G4DecayProducts* decayProducts =
193 (G4DecayProducts*)(dynamicParent->Ge
194 if (decayProducts == nullptr) {
195 G4ParticleDefinition* parentDef = th
196 G4DecayTable* decayTable = (parentDe
197 if (decayTable != nullptr) {
198 G4double parentMass = dynamicParen
199 G4VDecayChannel* decayChannel = de
200 if (decayChannel != nullptr) {
201 decayProducts = decayChannel->De
202 if (!decayProducts->IsChecked())
203 if (IsBoostToLabEnabled()) {
204 // boost all decay products to
205 decayProducts->Boost(dynamicPa
206 dynamicPa
207 }
208 }
209 else {
210 decayProducts = new G4DecayProdu
211 }
212 dynamicParent->SetPreAssignedDecay
213 }
214 }
215 else {
216 G4cout << "WARNING : already present
217 }
218 }
219 }
220 }
221
222 // G4cout << theStep->GetPostStepPoint()->Ge
223
224 // If the primary decays somewhere inside th
225 if (theStep->GetTrack()->GetParentID() == 0
226 && theStep->GetPostStepPoint()->GetProce
227 && theStep->GetPostStepPoint()->GetProce
228 != std::string::npos)
229 {
230 // Get properties of the primary particle
231
232 //--- Get values in different ways and che
233 // Kinetic energy of the primary particle
234 const G4double ekin_dynamicParticle =
235 theStep->GetTrack()->GetDynamicParticle(
236 const G4double ekin_track = theStep->GetTr
237 const G4double ekin_postStepPoint = theSte
238 const G4double ekin_deltaMax = std::max(st
239 st
240 // G4cout << "\t ekin_deltaMax [eV] = " <<
241 const G4double ekin_val = ekin_dynamicPart
242 // Total energy of the primary particle at
243 const G4double etot_dynamicParticle =
244 theStep->GetTrack()->GetDynamicParticle(
245 const G4double etot_track = theStep->GetTr
246 const G4double etot_postStepPoint = theSte
247 const G4double etot_deltaMax = std::max(st
248 st
249 // G4cout << "\t etot_deltaMax [eV] = " <<
250 const G4double etot_val = etot_dynamicPart
251 // Module of the 3-momentum of the primary
252 const G4double p_dynamicParticle =
253 theStep->GetTrack()->GetDynamicParticle(
254 const G4double p_track = theStep->GetTrack
255 const G4double p_postStepPoint = theStep->
256 const G4double p_deltaMax = std::max(std::
257 std::
258 // G4cout << "\t p_deltaMax [eV] = " << p_
259 const G4double p_val = p_dynamicParticle;
260 // 3-momentum direction (adimensional) of
261 const G4ThreeVector pdir_dynamicParticle =
262 theStep->GetTrack()->GetDynamicParticle(
263 const G4ThreeVector pdir_track = theStep->
264 const G4ThreeVector pdir_postStepPoint = t
265 const G4double pdir_x_deltaMax =
266 std::max(std::abs(pdir_dynamicParticle.x
267 std::abs(pdir_dynamicParticle.x
268 const G4double pdir_y_deltaMax =
269 std::max(std::abs(pdir_dynamicParticle.y
270 std::abs(pdir_dynamicParticle.y
271 const G4double pdir_z_deltaMax =
272 std::max(std::abs(pdir_dynamicParticle.z
273 std::abs(pdir_dynamicParticle.z
274 const G4double pdir_deltaMax =
275 std::max(std::max(pdir_x_deltaMax, pdir_
276 // G4cout << "\t pdir_deltaMax = " << pdir
277 // Mass of the primary particle at the de
278 const G4double mass_dynamicParticle = theS
279 const G4double mass_preStepPoint = theStep
280 const G4double mass_postStepPoint = theSte
281 const G4double mass_from_etot_ekin = etot_
282 const G4double mass_from4mom = std::sqrt(e
283 G4double mass_deltaMax1 = std::max(std::ab
284 std::ab
285 G4double mass_deltaMax2 = std::abs(mass_dy
286 G4double mass_deltaMax3 = std::abs(mass_dy
287 fMeanMass_deltaMax3 += mass_deltaMax3;
288 // G4cout << "\t mass_deltaMax{1,2,3} [eV]
289 // << mass_deltaMax2 / CLHEP::eV <<
290 const G4double mass_val = mass_dynamicPart
291 // Lorentz beta of the primary particle at
292 // The following line works only for G4 ve
293 const G4double beta_dynamicParticle = theS
294 const G4double beta_postStepPoint = theSte
295 // Before-10.7 const G4double beta_dynami
296 const G4double beta_velocity_track = theSt
297 const G4double beta_velocity_postStepPoint
298 theStep->GetPostStepPoint()->GetVelocity
299 const G4double beta_p_over_etot = p_val /
300 G4double beta_deltaMax1 = std::max(std::ab
301 std::ab
302 beta_deltaMax1 =
303 std::max(beta_deltaMax1, std::abs(beta_d
304 const G4double beta_deltaMax2 = std::abs(b
305 // G4cout << "\t beta_deltaMax{1,2} = " <<
306 const G4double beta_val = beta_dynamicPart
307 // Lorentz gamma of the primary particle a
308 const G4double gamma_postStepPoint = theSt
309 const G4double gamma_from_e_over_m = etot_
310 const G4double gamma_deltaMax1 = std::abs(
311 G4double gamma_from_beta = 0.0;
312 G4double gamma_deltaMax2 = 0.0;
313 G4double gamma_deltaMax3 = 0.0;
314 if (beta_val < 1.0) {
315 gamma_from_beta = 1.0 / std::sqrt(1.0 -
316 gamma_deltaMax2 = std::abs(gamma_postSte
317 gamma_deltaMax3 = std::abs(gamma_from_e_
318 }
319 const G4double gamma_val = gamma_postStepP
320 // G4cout << "\t gamma_deltaMax{1,2,3} = "
321 // << " , " << gamma_deltaMax3 <<
322 // << " ; gamma_from_e_over_m = " <
323 // << gamma_from_beta << G4endl;
324 // Proper time of the primary particle at
325 const G4double t_proper_track = theStep->G
326 const G4double t_proper_postStepPoint = th
327 const G4double t_proper_deltaMax = std::ab
328 // G4cout << "\t t_proper_deltaMax [fs] =
329 const G4double t_proper_val = t_proper_tra
330 // Lab time of the primary particle at the
331 // (Note: it would be wrong to trying to c
332 // above proper time via the simple
333 // const G4double t_lab_from_gamm
334 // because the gamma value of the p
335 // during its lifetime.)
336 const G4double t_local_track = theStep->Ge
337 const G4double t_local_postStepPoint = the
338 const G4double t_global_track = theStep->G
339 const G4double t_global_postStepPoint = th
340 G4double t_lab_deltaMax = std::max(std::ab
341 std::ab
342 t_lab_deltaMax = std::max(t_lab_deltaMax,
343 // G4cout << "\t t_lab_deltaMax [fs] = " <
344 const G4double t_lab_val = t_local_track;
345 // "MC-truth" decay radius of the primary
346 // (defined as the one that would happen i
347 // nor interactions with matter).
348 const G4double primaryBeta =
349 fPrimaryParticleInitialMomentum / fPrima
350 const G4double mc_truth_rPos1 = t_lab_val
351 const G4double mc_truth_rPos2 = t_lab_val
352 const G4double mc_truth_rPos_deltaMax = st
353 fMeanMc_truth_rPos_deltaMax += mc_truth_rP
354 // G4cout << "\t mc_truth_rPos_deltaMax [m
355 // << mc_truth_rPos_deltaMax / CLHE
356 if (mc_truth_rPos_deltaMax > ToleranceDelt
357 // G4cout << std::setprecision(6)
358 // << " Large : mc_truth_rPos_del
359 // << mc_truth_rPos_deltaMax / CL
360 // << " ; " << mc_truth_rPos1 <<
361 if (fRunPtr) fRunPtr->IncrementNumber_mc
362 }
363 const G4double mc_truth_rPos_val = mc_trut
364 // Keep note of the biggest discrepancies
365 fMaxEkin_deltaMax = std::max(fMaxEkin_delt
366 fMaxEtot_deltaMax = std::max(fMaxEtot_delt
367 fMaxP_deltaMax = std::max(fMaxP_deltaMax,
368 fMaxPdir_deltaMax = std::max(fMaxPdir_delt
369 fMaxMass_deltaMax1 = std::max(fMaxMass_del
370 fMaxMass_deltaMax2 = std::max(fMaxMass_del
371 fMaxMass_deltaMax3 = std::max(fMaxMass_del
372 fMaxBeta_deltaMax1 = std::max(fMaxBeta_del
373 fMaxBeta_deltaMax2 = std::max(fMaxBeta_del
374 fMaxGamma_deltaMax1 = std::max(fMaxGamma_d
375 fMaxGamma_deltaMax2 = std::max(fMaxGamma_d
376 fMaxGamma_deltaMax3 = std::max(fMaxGamma_d
377 fMaxT_lab_deltaMax = std::max(fMaxT_lab_de
378 fMaxT_proper_deltaMax = std::max(fMaxT_pro
379 fMaxMc_truth_rPos_deltaMax = std::max(fMax
380 //--- End consistency checks ---
381
382 // Global position
383 const G4double xPos = theStep->GetPostStep
384 const G4double yPos = theStep->GetPostStep
385 const G4double zPos = theStep->GetPostStep
386 const G4double rPos = std::sqrt(xPos * xPo
387 // I have verified that for this case in w
388 // "GetGlobalTime()" is the same as "GetLo
389 // Moreover, this value is also the same a
390 G4double tPos = theStep->GetPostStepPoint(
391 // The "MC-truth" decay radius is defined
392 // neither magnetic field effects nor inte
393 const G4double mc_truth_rPos = tPos * fPri
394 const G4double rDeltaPos = mc_truth_rPos -
395 const G4double eKin = theStep->GetPostStep
396 const G4double xMom = theStep->GetPostStep
397 const G4double yMom = theStep->GetPostStep
398 const G4double zMom = theStep->GetPostStep
399 // The compute here the angular deflection
400 // the primary particle - which is along t
401 G4double xDirection = std::min(theStep->Ge
402 if (xDirection < -1.0) xDirection = -1.0;
403 const G4double deflection_angle_in_degrees
404 const G4double delta_ekin = fPrimaryPartic
405 // G4cout << std::setprecision(6)
406 // << " Decay: tPos[ns]=" << tPos <
407 // << rDeltaPos /CLHEP::micrometer
408 // << " ; deltaAngle(deg)=" << defl
409 // If the absolute difference between the
410 // given threshold, then we notify this s
411 // for post-processing evaluation. Moreov
412 // smaller than the real one, then we cou
413 // this special situation in the output w
414 // evaluation.
415 if (std::abs(rDeltaPos) > ToleranceDeltaDe
416 // G4cout << "\t LARGE_DELTA_R : mc_trut
417 // << " ; rPos[mm]=" << rPos;
418 if (rDeltaPos < 0.0) {
419 // G4cout << "\t ***UNEXPECTED***";
420 if (fRunPtr) fRunPtr->IncrementNumberU
421 }
422 // G4cout << G4endl;
423 }
424 fMeanDeltaR_primaryDecay += rDeltaPos;
425 fMinDeltaR_primaryDecay = std::min(fMinDel
426 fMaxDeltaR_primaryDecay = std::max(fMaxDel
427 fMeanR_primaryDecay += rPos;
428 fMinR_primaryDecay = std::min(fMinR_primar
429 fMaxR_primaryDecay = std::max(fMaxR_primar
430 fMeanX_primaryDecay += xPos;
431 fMinX_primaryDecay = std::min(fMinX_primar
432 fMaxX_primaryDecay = std::max(fMaxX_primar
433 fMeanY_primaryDecay += yPos;
434 fMinY_primaryDecay = std::min(fMinY_primar
435 fMaxY_primaryDecay = std::max(fMaxY_primar
436 fMeanZ_primaryDecay += zPos;
437 fMinZ_primaryDecay = std::min(fMinZ_primar
438 fMaxZ_primaryDecay = std::max(fMaxZ_primar
439 fMeanDeltaAngle_primaryDecay += deflection
440 fMinDeltaAngle_primaryDecay =
441 std::min(fMinDeltaAngle_primaryDecay, de
442 fMaxDeltaAngle_primaryDecay =
443 std::max(fMaxDeltaAngle_primaryDecay, de
444 fMeanDeltaEkin_primaryDecay += delta_ekin;
445 fMinDeltaEkin_primaryDecay = std::min(fMin
446 fMaxDeltaEkin_primaryDecay = std::max(fMax
447 fMeanEkin_primaryDecay += eKin;
448 fMinEkin_primaryDecay = std::min(fMinEkin_
449 fMaxEkin_primaryDecay = std::max(fMaxEkin_
450 fMeanPx_primaryDecay += xMom;
451 fMinPx_primaryDecay = std::min(fMinPx_prim
452 fMaxPx_primaryDecay = std::max(fMaxPx_prim
453 fMeanPy_primaryDecay += yMom;
454 fMinPy_primaryDecay = std::min(fMinPy_prim
455 fMaxPy_primaryDecay = std::max(fMaxPy_prim
456 fMeanPz_primaryDecay += zMom;
457 fMinPz_primaryDecay = std::min(fMinPz_prim
458 fMaxPz_primaryDecay = std::max(fMaxPz_prim
459
460 //--- Extra checks ---
461 // Compute the "MC-truth" decay radius usi
462 // it decays.
463 // To do this, we would need an "effective
464 // particle during its lifetime, whereas i
465 // values at the beginning and at the end
466 // overestimate of the "MC-truth" decay ra
467 // or an underestimate of it - by using th
468 // We want to check the average values and
469 const G4double underestimated_mc_truth_rPo
470 t_proper_val * gamma_val * beta_val * CL
471 const G4double overestimated_mc_truth_rPos
472 t_proper_val * fPrimaryParticleInitialGa
473 const G4double underestimated_mc_truth_rPo
474 underestimated_mc_truth_rPos - mc_truth_
475 const G4double overestimated_mc_truth_rPos
476 overestimated_mc_truth_rPos - mc_truth_r
477 fMeanUnderestimated_mc_truth_rPos_delta +=
478 fMeanOverestimated_mc_truth_rPos_delta +=
479 // G4cout << "\t underestimated_mc_truth_r
480 // << underestimated_mc_truth_rPos_
481 // << " ; overestimated_mc_truth_rP
482 // << overestimated_mc_truth_rPos_d
483 if (-underestimated_mc_truth_rPos_delta >
484 // G4cout << std::setprecision(6)
485 // << " Large : underestimated_mc
486 // << underestimated_mc_truth_rPo
487 // << " ; " << underestimated_mc_
488 // << mc_truth_rPos_val << " mm"
489 if (fRunPtr) fRunPtr->IncrementNumber_un
490 }
491 if (overestimated_mc_truth_rPos_delta > To
492 // G4cout << std::setprecision(6)
493 // << " Large : overestimated_mc_
494 // << overestimated_mc_truth_rPos
495 // << " ; " << overestimated_mc_t
496 // << mc_truth_rPos_val << " mm"
497 if (fRunPtr) fRunPtr->IncrementNumber_ov
498 }
499 const G4double underestimated_rDeltaPos =
500 const G4double overestimated_rDeltaPos = o
501 fMeanUnderestimated_rDeltaPos += underesti
502 fMeanOverestimated_rDeltaPos += overestima
503 // G4cout << std::setprecision(6)
504 // << "\t underestimated_rDeltaPos=
505 // << " ; overestimated_rDeltaPos="
506 // << " mum" << G4endl;
507 if (-underestimated_rDeltaPos > ToleranceD
508 if (fRunPtr) fRunPtr->IncrementNumberLar
509 }
510 if (overestimated_rDeltaPos > ToleranceDel
511 if (fRunPtr) fRunPtr->IncrementNumberLar
512 }
513 // Keep note of the biggest discrepancies
514 fMinUnderestimated_mc_truth_rPos_delta =
515 std::min(fMinUnderestimated_mc_truth_rPo
516 fMaxOverestimated_mc_truth_rPos_delta =
517 std::max(fMaxOverestimated_mc_truth_rPos
518 fMinUnderestimated_rDeltaPos = std::min(fM
519 fMaxOverestimated_rDeltaPos = std::max(fMa
520 //--- End extra checks ---
521
522 // Check numerical errors due to the use o
523 // several, equivalent computations, takin
524 const G4float float_xPos = static_cast<G4f
525 const G4float float_yPos = static_cast<G4f
526 const G4float float_zPos = static_cast<G4f
527 const G4float float_rPos =
528 std::sqrt(float_xPos * float_xPos + floa
529 const G4float float_tPos = static_cast<G4f
530 const G4float float_initialBeta1 = static_
531 const G4float float_initialBeta2 = static_
532 / stati
533 const G4float float_initialGamma = static_
534 const G4float float_initialBeta3 =
535 std::sqrt(float_initialGamma * float_ini
536 const G4float float_c_light = static_cast<
537 const G4float float_mc_truth_rPos1 = float
538 const G4float float_mc_truth_rPos2 = float
539 const G4float float_mc_truth_rPos3 = float
540 const G4float float_rDeltaPos_0 = static_c
541 const G4float float_rDeltaPos_1 = float_mc
542 const G4float float_rDeltaPos_2 = float_mc
543 const G4float float_rDeltaPos_3 = float_mc
544 const G4float float_rDeltaPos_4 = static_c
545 const G4float float_rDeltaPos_5 = float_mc
546 const G4float float_rDeltaPos_6 = float_mc
547 const G4float float_rDeltaPos_7 = float_mc
548 G4double rDeltaPos_deltaMax =
549 std::max(std::abs(float_rDeltaPos_0 - rD
550 rDeltaPos_deltaMax = std::max(rDeltaPos_de
551 rDeltaPos_deltaMax = std::max(rDeltaPos_de
552 rDeltaPos_deltaMax = std::max(rDeltaPos_de
553 rDeltaPos_deltaMax = std::max(rDeltaPos_de
554 rDeltaPos_deltaMax = std::max(rDeltaPos_de
555 rDeltaPos_deltaMax = std::max(rDeltaPos_de
556 // G4cout << std::setprecision(6) << " rDe
557 // << rDeltaPos_deltaMax / CLHEP::m
558 fMaxFloat_rDeltaPos_deltaMax = std::max(fM
559
560 // Get properties of the decay products an
561 // decay
562 std::size_t nSec = theStep->GetNumberOfSec
563 const std::vector<const G4Track*>* ptrVecS
564 G4double deltaE = 0.0, deltaPx = 0.0, delt
565 if (nSec > 0 && ptrVecSecondaries != nullp
566 G4double sumEsecondaries = 0.0;
567 G4ThreeVector sumPsecondaries(0.0, 0.0,
568 for (std::size_t i = 0; i < nSec; ++i) {
569 if ((*ptrVecSecondaries)[i]) {
570 sumEsecondaries += (*ptrVecSecondari
571 sumPsecondaries += (*ptrVecSecondari
572 }
573 }
574 deltaE = sumEsecondaries - theStep->GetP
575 fMeanViolationE_primaryDecay += deltaE;
576 fMinViolationE_primaryDecay = std::min(f
577 fMaxViolationE_primaryDecay = std::max(f
578 if (std::abs(deltaE) > ToleranceEPviolat
579 if (fRunPtr) fRunPtr->IncrementNumberE
580 }
581 deltaPx = sumPsecondaries.x() - xMom;
582 fMeanViolationPx_primaryDecay += deltaPx
583 fMinViolationPx_primaryDecay = std::min(
584 fMaxViolationPx_primaryDecay = std::max(
585 deltaPy = sumPsecondaries.y() - yMom;
586 fMeanViolationPy_primaryDecay += deltaPy
587 fMinViolationPy_primaryDecay = std::min(
588 fMaxViolationPy_primaryDecay = std::max(
589 deltaPz = sumPsecondaries.z() - zMom;
590 fMeanViolationPz_primaryDecay += deltaPz
591 fMinViolationPz_primaryDecay = std::min(
592 fMaxViolationPz_primaryDecay = std::max(
593 if (std::abs(deltaPx) > ToleranceEPviola
594 || std::abs(deltaPz) > ToleranceEPvi
595 {
596 if (fRunPtr) fRunPtr->IncrementNumberP
597 }
598 }
599 else {
600 if (fRunPtr) fRunPtr->IncrementNumberBad
601 }
602
603 if (fRunPtr) {
604 fRunPtr->IncrementNumberDecays();
605 fRunPtr->SetDecayT(tPos);
606 fRunPtr->SetDecayR_mc_truth(mc_truth_rPo
607 fRunPtr->SetDecayR(rPos);
608 fRunPtr->SetDecayX(xPos);
609 fRunPtr->SetDecayY(yPos);
610 fRunPtr->SetDecayZ(zPos);
611 fRunPtr->SetDeltaDecayR(rDeltaPos);
612 fRunPtr->SetDeflectionAngle(deflection_a
613 fRunPtr->SetDeltaEkin(delta_ekin);
614 fRunPtr->SetDecayEkin(eKin);
615 fRunPtr->SetDecayPx(xMom);
616 fRunPtr->SetDecayPy(yMom);
617 fRunPtr->SetDecayPz(zMom);
618 fRunPtr->SetDecayEtotViolation(deltaE);
619 fRunPtr->SetDecayPxViolation(deltaPx);
620 fRunPtr->SetDecayPyViolation(deltaPy);
621 fRunPtr->SetDecayPzViolation(deltaPz);
622 fRunPtr->SetMaxEkin_deltaMax(ekin_deltaM
623 fRunPtr->SetMaxEtot_deltaMax(etot_deltaM
624 fRunPtr->SetMaxP_deltaMax(p_deltaMax);
625 fRunPtr->SetMaxPdir_deltaMax(pdir_deltaM
626 fRunPtr->SetMaxMass_deltaMax1(mass_delta
627 fRunPtr->SetMaxMass_deltaMax2(mass_delta
628 fRunPtr->SetMaxMass_deltaMax3(mass_delta
629 fRunPtr->SetMaxBeta_deltaMax1(beta_delta
630 fRunPtr->SetMaxBeta_deltaMax2(beta_delta
631 fRunPtr->SetMaxGamma_deltaMax1(gamma_del
632 fRunPtr->SetMaxGamma_deltaMax2(gamma_del
633 fRunPtr->SetMaxGamma_deltaMax3(gamma_del
634 fRunPtr->SetMaxT_proper_deltaMax(t_prope
635 fRunPtr->SetMaxT_lab_deltaMax(t_lab_delt
636 fRunPtr->SetMaxMc_truth_rPos_deltaMax(mc
637 fRunPtr->SetMinUnderestimated_mc_truth_r
638 fRunPtr->SetMaxOverestimated_mc_truth_rP
639 fRunPtr->SetMinUnderestimated_rDeltaPos(
640 fRunPtr->SetMaxOverestimated_rDeltaPos(o
641 fRunPtr->SetMaxFloat_rDeltaPos_deltaMax(
642 }
643 }
644 }
645
646 //....oooOO0OOooo........oooOO0OOooo........oo
647