Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/examples/extended/hadronic/Hadr10/src/SteppingAction.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

  1 //
  2 // ********************************************************************
  3 // * License and Disclaimer                                           *
  4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.                             *
 10 // *                                                                  *
 11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                                                  *
 18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // ********************************************************************
 25 //
 26 /// \file SteppingAction.cc
 27 /// \brief Implementation of the SteppingAction class
 28 //
 29 //
 30 
 31 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 32 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 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........oooOO0OOooo........oooOO0OOooo......
 53 
 54 SteppingAction::SteppingAction() : G4UserSteppingAction()
 55 {
 56   Initialize();
 57 }
 58 
 59 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 60 
 61 SteppingAction::~SteppingAction() {}
 62 
 63 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 64 
 65 void SteppingAction::Initialize()
 66 {
 67   // Initialization needed at the beginning of each Run
 68   fRunPtr = nullptr;
 69   fToleranceEPviolations = 1.0 * CLHEP::eV;  //***LOOKHERE***
 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 = G4ThreeVector(0.0, 0.0, 0.0);
 77   fPrimaryParticleInitialPosition = G4ThreeVector(0.0, 0.0, 0.0);
 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 = 9999999.9;
129   fMaxOverestimated_mc_truth_rPos_delta = -9999999.9;
130   fMeanUnderestimated_mc_truth_rPos_delta = 0.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........oooOO0OOooo........oooOO0OOooo......
152 
153 void SteppingAction::UserSteppingAction(const G4Step* theStep)
154 {
155   // Store the information about the ID and the kinetic energy of the primary particle,
156   // at the first step of the first event.
157   // Note that for the kinetic energy, we are considering the "pre-step" point of such first step.
158   if (theStep->GetTrack()->GetParentID() == 0 && theStep->GetTrack()->GetCurrentStepNumber() == 1) {
159     fPrimaryParticleId = theStep->GetTrack()->GetDefinition()->GetPDGEncoding();
160     fPrimaryParticleInitialKineticEnergy = theStep->GetPreStepPoint()->GetKineticEnergy();
161     fPrimaryParticleInitialTotalEnergy = theStep->GetPreStepPoint()->GetTotalEnergy();
162     fPrimaryParticleInitial3Momentum = theStep->GetPreStepPoint()->GetMomentum();
163     fPrimaryParticleInitialMomentum = fPrimaryParticleInitial3Momentum.mag();
164     fPrimaryParticleInitialPosition = theStep->GetPreStepPoint()->GetPosition();
165     fPrimaryParticleInitialBeta = theStep->GetPreStepPoint()->GetBeta();
166     fPrimaryParticleInitialGamma = theStep->GetPreStepPoint()->GetGamma();
167     // As tolerance for EP violations, consider the max value between the default value
168     // and 1 billionth of the initial, primary particle kinetic energy.
169     if (fToleranceEPviolations < fPrimaryParticleInitialKineticEnergy * 1.0e-9) {
170       fToleranceEPviolations = fPrimaryParticleInitialKineticEnergy * 1.0e-9;
171     }
172     // Set the values of this run to the Run object
173     if (fRunPtr) {
174       fRunPtr->SetPrimaryParticleId(fPrimaryParticleId);
175       fRunPtr->SetPrimaryParticleInitialKineticEnergy(fPrimaryParticleInitialKineticEnergy);
176       fRunPtr->SetPrimaryParticleInitialTotalEnergy(fPrimaryParticleInitialTotalEnergy);
177       fRunPtr->SetPrimaryParticleInitialMomentum(fPrimaryParticleInitialMomentum);
178       fRunPtr->SetPrimaryParticleInitialBeta(fPrimaryParticleInitialBeta);
179       fRunPtr->SetPrimaryParticleInitialGamma(fPrimaryParticleInitialGamma);
180       fRunPtr->SetPrimaryParticleInitial3Momentum(fPrimaryParticleInitial3Momentum);
181       fRunPtr->SetPrimaryParticleInitialPosition(fPrimaryParticleInitialPosition);
182       fRunPtr->SetToleranceEPviolations(ToleranceEPviolations());
183       fRunPtr->SetToleranceDeltaDecayRadius(ToleranceDeltaDecayRadius());
184       fRunPtr->SetIsPreassignedDecayEnabled(IsPreassignedDecayEnabled());
185       fRunPtr->SetIsBoostToLabEnabled(IsBoostToLabEnabled());
186     }
187     // Use the preassigned decay is enabled
188     if (IsPreassignedDecayEnabled() && (!theStep->GetTrack()->GetDefinition()->GetPDGStable())) {
189       G4DynamicParticle* dynamicParent =
190         const_cast<G4DynamicParticle*>(theStep->GetTrack()->GetDynamicParticle());
191       if (dynamicParent != nullptr) {
192         G4DecayProducts* decayProducts =
193           (G4DecayProducts*)(dynamicParent->GetPreAssignedDecayProducts());
194         if (decayProducts == nullptr) {
195           G4ParticleDefinition* parentDef = theStep->GetTrack()->GetDefinition();
196           G4DecayTable* decayTable = (parentDef == nullptr ? nullptr : parentDef->GetDecayTable());
197           if (decayTable != nullptr) {
198             G4double parentMass = dynamicParent->GetMass();
199             G4VDecayChannel* decayChannel = decayTable->SelectADecayChannel(parentMass);
200             if (decayChannel != nullptr) {
201               decayProducts = decayChannel->DecayIt(parentMass);
202               if (!decayProducts->IsChecked()) decayProducts->DumpInfo();
203               if (IsBoostToLabEnabled()) {
204                 // boost all decay products to laboratory frame
205                 decayProducts->Boost(dynamicParent->GetTotalEnergy(),
206                                      dynamicParent->GetMomentumDirection());
207               }
208             }
209             else {
210               decayProducts = new G4DecayProducts(*dynamicParent);
211             }
212             dynamicParent->SetPreAssignedDecayProducts(decayProducts);
213           }
214         }
215         else {
216           G4cout << "WARNING : already present preassign decay !" << G4endl;
217         }
218       }
219     }
220   }
221 
222   // G4cout << theStep->GetPostStepPoint()->GetProcessDefinedStep()->GetProcessName() << G4endl;
223 
224   // If the primary decays somewhere inside the World volume, get the information about the decay
225   if (theStep->GetTrack()->GetParentID() == 0
226       && theStep->GetPostStepPoint()->GetProcessDefinedStep() != nullptr
227       && theStep->GetPostStepPoint()->GetProcessDefinedStep()->GetProcessName().find("Decay")
228            != std::string::npos)
229   {
230     // Get properties of the primary particle when it decays
231 
232     //--- Get values in different ways and check their consistency ---
233     // Kinetic energy of the primary particle at the decay
234     const G4double ekin_dynamicParticle =
235       theStep->GetTrack()->GetDynamicParticle()->GetKineticEnergy();
236     const G4double ekin_track = theStep->GetTrack()->GetKineticEnergy();
237     const G4double ekin_postStepPoint = theStep->GetPostStepPoint()->GetKineticEnergy();
238     const G4double ekin_deltaMax = std::max(std::abs(ekin_dynamicParticle - ekin_track),
239                                             std::abs(ekin_dynamicParticle - ekin_postStepPoint));
240     // G4cout << "\t ekin_deltaMax [eV] = " << ekin_deltaMax / CLHEP::eV << G4endl;
241     const G4double ekin_val = ekin_dynamicParticle;  // To be used later
242     // Total energy of the primary particle at the decay
243     const G4double etot_dynamicParticle =
244       theStep->GetTrack()->GetDynamicParticle()->GetTotalEnergy();
245     const G4double etot_track = theStep->GetTrack()->GetTotalEnergy();
246     const G4double etot_postStepPoint = theStep->GetPostStepPoint()->GetTotalEnergy();
247     const G4double etot_deltaMax = std::max(std::abs(etot_dynamicParticle - etot_track),
248                                             std::abs(etot_dynamicParticle - etot_postStepPoint));
249     // G4cout << "\t etot_deltaMax [eV] = " << etot_deltaMax / CLHEP::eV << G4endl;
250     const G4double etot_val = etot_dynamicParticle;  // To be used later
251     // Module of the 3-momentum of the primary particle at the decay
252     const G4double p_dynamicParticle =
253       theStep->GetTrack()->GetDynamicParticle()->GetMomentum().mag();
254     const G4double p_track = theStep->GetTrack()->GetMomentum().mag();
255     const G4double p_postStepPoint = theStep->GetPostStepPoint()->GetMomentum().mag();
256     const G4double p_deltaMax = std::max(std::abs(p_dynamicParticle - p_track),
257                                          std::abs(p_dynamicParticle - p_postStepPoint));
258     // G4cout << "\t p_deltaMax [eV] = " << p_deltaMax / CLHEP::eV << G4endl;
259     const G4double p_val = p_dynamicParticle;  // To be used later
260     // 3-momentum direction (adimensional) of the primary particle at the decay
261     const G4ThreeVector pdir_dynamicParticle =
262       theStep->GetTrack()->GetDynamicParticle()->GetMomentumDirection();
263     const G4ThreeVector pdir_track = theStep->GetTrack()->GetMomentumDirection();
264     const G4ThreeVector pdir_postStepPoint = theStep->GetPostStepPoint()->GetMomentumDirection();
265     const G4double pdir_x_deltaMax =
266       std::max(std::abs(pdir_dynamicParticle.x() - pdir_track.x()),
267                std::abs(pdir_dynamicParticle.x() - pdir_postStepPoint.x()));
268     const G4double pdir_y_deltaMax =
269       std::max(std::abs(pdir_dynamicParticle.y() - pdir_track.y()),
270                std::abs(pdir_dynamicParticle.y() - pdir_postStepPoint.y()));
271     const G4double pdir_z_deltaMax =
272       std::max(std::abs(pdir_dynamicParticle.z() - pdir_track.z()),
273                std::abs(pdir_dynamicParticle.z() - pdir_postStepPoint.z()));
274     const G4double pdir_deltaMax =
275       std::max(std::max(pdir_x_deltaMax, pdir_y_deltaMax), pdir_z_deltaMax);
276     // G4cout << "\t pdir_deltaMax = " << pdir_deltaMax << G4endl;
277     //  Mass of the primary particle at the decay
278     const G4double mass_dynamicParticle = theStep->GetTrack()->GetDynamicParticle()->GetMass();
279     const G4double mass_preStepPoint = theStep->GetPreStepPoint()->GetMass();
280     const G4double mass_postStepPoint = theStep->GetPostStepPoint()->GetMass();
281     const G4double mass_from_etot_ekin = etot_val - ekin_val;
282     const G4double mass_from4mom = std::sqrt(etot_val * etot_val - p_val * p_val);
283     G4double mass_deltaMax1 = std::max(std::abs(mass_dynamicParticle - mass_preStepPoint),
284                                        std::abs(mass_dynamicParticle - mass_postStepPoint));
285     G4double mass_deltaMax2 = std::abs(mass_dynamicParticle - mass_from_etot_ekin);
286     G4double mass_deltaMax3 = std::abs(mass_dynamicParticle - mass_from4mom);
287     fMeanMass_deltaMax3 += mass_deltaMax3;
288     // G4cout << "\t mass_deltaMax{1,2,3} [eV] = " << mass_deltaMax1 / CLHEP::eV << "\t"
289     //        << mass_deltaMax2 / CLHEP::eV << "\t" << mass_deltaMax3 / CLHEP::eV << G4endl;
290     const G4double mass_val = mass_dynamicParticle;  // To be used later
291     // Lorentz beta of the primary particle at the decay
292     // The following line works only for G4 versions >= 10.7
293     const G4double beta_dynamicParticle = theStep->GetTrack()->GetDynamicParticle()->GetBeta();
294     const G4double beta_postStepPoint = theStep->GetPostStepPoint()->GetBeta();
295     // Before-10.7  const G4double beta_dynamicParticle = beta_postStepPoint;
296     const G4double beta_velocity_track = theStep->GetTrack()->GetVelocity() / CLHEP::c_light;
297     const G4double beta_velocity_postStepPoint =
298       theStep->GetPostStepPoint()->GetVelocity() / CLHEP::c_light;
299     const G4double beta_p_over_etot = p_val / etot_val;
300     G4double beta_deltaMax1 = std::max(std::abs(beta_dynamicParticle - beta_postStepPoint),
301                                        std::abs(beta_dynamicParticle - beta_velocity_track));
302     beta_deltaMax1 =
303       std::max(beta_deltaMax1, std::abs(beta_dynamicParticle - beta_velocity_postStepPoint));
304     const G4double beta_deltaMax2 = std::abs(beta_dynamicParticle - beta_p_over_etot);
305     // G4cout << "\t beta_deltaMax{1,2} = " << beta_deltaMax1 << " , " << beta_deltaMax2 << G4endl;
306     const G4double beta_val = beta_dynamicParticle;  // To be used later
307     // Lorentz gamma of the primary particle at the decay
308     const G4double gamma_postStepPoint = theStep->GetPostStepPoint()->GetGamma();
309     const G4double gamma_from_e_over_m = etot_val / mass_val;
310     const G4double gamma_deltaMax1 = std::abs(gamma_postStepPoint - gamma_from_e_over_m);
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 - beta_val * beta_val);
316       gamma_deltaMax2 = std::abs(gamma_postStepPoint - gamma_from_beta);
317       gamma_deltaMax3 = std::abs(gamma_from_e_over_m - gamma_from_beta);
318     }
319     const G4double gamma_val = gamma_postStepPoint;  // To be used later;
320     // G4cout << "\t gamma_deltaMax{1,2,3} = " << gamma_deltaMax1 << " , " << gamma_deltaMax2
321     //        << " , " << gamma_deltaMax3  << " ; gamma_postStepPoint = " << gamma_postStepPoint
322     //        << " ; gamma_from_e_over_m = " << gamma_from_e_over_m << " ; gamma_from_beta = "
323     //        << gamma_from_beta << G4endl;
324     //  Proper time of the primary particle at the decay
325     const G4double t_proper_track = theStep->GetTrack()->GetProperTime();
326     const G4double t_proper_postStepPoint = theStep->GetPostStepPoint()->GetProperTime();
327     const G4double t_proper_deltaMax = std::abs(t_proper_track - t_proper_postStepPoint);
328     // G4cout << "\t t_proper_deltaMax [fs] = " << t_proper_deltaMax / femtosecond << G4endl;
329     const G4double t_proper_val = t_proper_track;  // To be used later
330     // Lab time of the primary particle at the decay
331     // (Note: it would be wrong to trying to compute this lab time from the
332     //        above proper time via the simple formula:
333     //          const G4double t_lab_from_gamma = t_proper_val * gamma_val;
334     //        because the gamma value of the primary particle has changed
335     //        during its lifetime.)
336     const G4double t_local_track = theStep->GetTrack()->GetLocalTime();
337     const G4double t_local_postStepPoint = theStep->GetPostStepPoint()->GetLocalTime();
338     const G4double t_global_track = theStep->GetTrack()->GetGlobalTime();
339     const G4double t_global_postStepPoint = theStep->GetPostStepPoint()->GetGlobalTime();
340     G4double t_lab_deltaMax = std::max(std::abs(t_local_track - t_local_postStepPoint),
341                                        std::abs(t_local_track - t_global_track));
342     t_lab_deltaMax = std::max(t_lab_deltaMax, std::abs(t_local_track - t_global_postStepPoint));
343     // G4cout << "\t t_lab_deltaMax [fs] = " << t_lab_deltaMax / femtosecond << G4endl;
344     const G4double t_lab_val = t_local_track;  // To be used later
345     // "MC-truth" decay radius of the primary particle at the decay
346     // (defined as the one that would happen if there are neither magnetic field effects
347     // nor interactions with matter).
348     const G4double primaryBeta =
349       fPrimaryParticleInitialMomentum / fPrimaryParticleInitialTotalEnergy;
350     const G4double mc_truth_rPos1 = t_lab_val * fPrimaryParticleInitialBeta * CLHEP::c_light;
351     const G4double mc_truth_rPos2 = t_lab_val * primaryBeta * CLHEP::c_light;
352     const G4double mc_truth_rPos_deltaMax = std::abs(mc_truth_rPos1 - mc_truth_rPos2);
353     fMeanMc_truth_rPos_deltaMax += mc_truth_rPos_deltaMax;
354     // G4cout << "\t mc_truth_rPos_deltaMax [mum] = "
355     //        << mc_truth_rPos_deltaMax / CLHEP::micrometer << G4endl;
356     if (mc_truth_rPos_deltaMax > ToleranceDeltaDecayRadius()) {
357       // G4cout << std::setprecision(6)
358       //        << " Large : mc_truth_rPos_deltaMax [mum]="
359       //        << mc_truth_rPos_deltaMax / CLHEP::micrometer
360       //        << " ; " << mc_truth_rPos1 << " , " << mc_truth_rPos2 << " mm" << G4endl;
361       if (fRunPtr) fRunPtr->IncrementNumber_mc_truth_rPos_deltaMax_above();
362     }
363     const G4double mc_truth_rPos_val = mc_truth_rPos1;  // To be used later
364     // Keep note of the biggest discrepancies
365     fMaxEkin_deltaMax = std::max(fMaxEkin_deltaMax, ekin_deltaMax);
366     fMaxEtot_deltaMax = std::max(fMaxEtot_deltaMax, etot_deltaMax);
367     fMaxP_deltaMax = std::max(fMaxP_deltaMax, p_deltaMax);
368     fMaxPdir_deltaMax = std::max(fMaxPdir_deltaMax, pdir_deltaMax);
369     fMaxMass_deltaMax1 = std::max(fMaxMass_deltaMax1, mass_deltaMax1);
370     fMaxMass_deltaMax2 = std::max(fMaxMass_deltaMax2, mass_deltaMax2);
371     fMaxMass_deltaMax3 = std::max(fMaxMass_deltaMax3, mass_deltaMax3);
372     fMaxBeta_deltaMax1 = std::max(fMaxBeta_deltaMax1, beta_deltaMax1);
373     fMaxBeta_deltaMax2 = std::max(fMaxBeta_deltaMax2, beta_deltaMax2);
374     fMaxGamma_deltaMax1 = std::max(fMaxGamma_deltaMax1, gamma_deltaMax1);
375     fMaxGamma_deltaMax2 = std::max(fMaxGamma_deltaMax2, gamma_deltaMax2);
376     fMaxGamma_deltaMax3 = std::max(fMaxGamma_deltaMax3, gamma_deltaMax3);
377     fMaxT_lab_deltaMax = std::max(fMaxT_lab_deltaMax, t_lab_deltaMax);
378     fMaxT_proper_deltaMax = std::max(fMaxT_proper_deltaMax, t_proper_deltaMax);
379     fMaxMc_truth_rPos_deltaMax = std::max(fMaxMc_truth_rPos_deltaMax, mc_truth_rPos_deltaMax);
380     //--- End consistency checks ---
381 
382     // Global position
383     const G4double xPos = theStep->GetPostStepPoint()->GetPosition().x();
384     const G4double yPos = theStep->GetPostStepPoint()->GetPosition().y();
385     const G4double zPos = theStep->GetPostStepPoint()->GetPosition().z();
386     const G4double rPos = std::sqrt(xPos * xPos + yPos * yPos + zPos * zPos);
387     // I have verified that for this case in which only primaries are considered, the
388     // "GetGlobalTime()" is the same as "GetLocalTime()" (the one we use).
389     // Moreover, this value is also the same as "GetProperTime()"*gamma .
390     G4double tPos = theStep->GetPostStepPoint()->GetLocalTime();
391     // The "MC-truth" decay radius is defined as the one that would happen if there are
392     // neither magnetic field effects nor interactions with matter.
393     const G4double mc_truth_rPos = tPos * fPrimaryParticleInitialBeta * CLHEP::c_light;
394     const G4double rDeltaPos = mc_truth_rPos - rPos;
395     const G4double eKin = theStep->GetPostStepPoint()->GetKineticEnergy();
396     const G4double xMom = theStep->GetPostStepPoint()->GetMomentum().x();
397     const G4double yMom = theStep->GetPostStepPoint()->GetMomentum().y();
398     const G4double zMom = theStep->GetPostStepPoint()->GetMomentum().z();
399     // The compute here the angular deflection, in degrees, between the initial direction of
400     // the primary particle - which is along the x-axis, and its direction when it decays.
401     G4double xDirection = std::min(theStep->GetPostStepPoint()->GetMomentumDirection().x(), 1.0);
402     if (xDirection < -1.0) xDirection = -1.0;
403     const G4double deflection_angle_in_degrees = 57.29 * std::acos(xDirection);
404     const G4double delta_ekin = fPrimaryParticleInitialKineticEnergy - eKin;
405     // G4cout << std::setprecision(6)
406     //        << " Decay: tPos[ns]=" << tPos << " ; rPos[mm]=" << rPos << " ; deltaR[mum]="
407     //        << rDeltaPos /CLHEP::micrometer << " ; deltaEkin[MeV]=" << delta_ekin
408     //        << " ; deltaAngle(deg)=" << deflection_angle_in_degrees << G4endl;
409     //  If the absolute difference between the "MC-truth" decay radius and the real one is above a
410     //  given threshold, then we notify this special situation in the output, with "LARGE_DELTA_R"
411     //  for post-processing evaluation. Moreover, in this case, if the "MC-truth" decay radius is
412     //  smaller than the real one, then we count this unexpected occurrence and we further notify
413     //  this special situation in the output with "***UNEXPECTED***" for post-processing
414     //  evaluation.
415     if (std::abs(rDeltaPos) > ToleranceDeltaDecayRadius()) {
416       // G4cout << "\t LARGE_DELTA_R : mc_truth_rPos[mm]=" << mc_truth_rPos
417       //        << " ; rPos[mm]=" << rPos;
418       if (rDeltaPos < 0.0) {
419         // G4cout << "\t ***UNEXPECTED***";
420         if (fRunPtr) fRunPtr->IncrementNumberUnexpectedDecays();
421       }
422       // G4cout << G4endl;
423     }
424     fMeanDeltaR_primaryDecay += rDeltaPos;
425     fMinDeltaR_primaryDecay = std::min(fMinDeltaR_primaryDecay, rDeltaPos);
426     fMaxDeltaR_primaryDecay = std::max(fMaxDeltaR_primaryDecay, rDeltaPos);
427     fMeanR_primaryDecay += rPos;
428     fMinR_primaryDecay = std::min(fMinR_primaryDecay, rPos);
429     fMaxR_primaryDecay = std::max(fMaxR_primaryDecay, rPos);
430     fMeanX_primaryDecay += xPos;
431     fMinX_primaryDecay = std::min(fMinX_primaryDecay, xPos);
432     fMaxX_primaryDecay = std::max(fMaxX_primaryDecay, xPos);
433     fMeanY_primaryDecay += yPos;
434     fMinY_primaryDecay = std::min(fMinY_primaryDecay, yPos);
435     fMaxY_primaryDecay = std::max(fMaxY_primaryDecay, yPos);
436     fMeanZ_primaryDecay += zPos;
437     fMinZ_primaryDecay = std::min(fMinZ_primaryDecay, zPos);
438     fMaxZ_primaryDecay = std::max(fMaxZ_primaryDecay, zPos);
439     fMeanDeltaAngle_primaryDecay += deflection_angle_in_degrees;
440     fMinDeltaAngle_primaryDecay =
441       std::min(fMinDeltaAngle_primaryDecay, deflection_angle_in_degrees);
442     fMaxDeltaAngle_primaryDecay =
443       std::max(fMaxDeltaAngle_primaryDecay, deflection_angle_in_degrees);
444     fMeanDeltaEkin_primaryDecay += delta_ekin;
445     fMinDeltaEkin_primaryDecay = std::min(fMinDeltaEkin_primaryDecay, delta_ekin);
446     fMaxDeltaEkin_primaryDecay = std::max(fMaxDeltaEkin_primaryDecay, delta_ekin);
447     fMeanEkin_primaryDecay += eKin;
448     fMinEkin_primaryDecay = std::min(fMinEkin_primaryDecay, eKin);
449     fMaxEkin_primaryDecay = std::max(fMaxEkin_primaryDecay, eKin);
450     fMeanPx_primaryDecay += xMom;
451     fMinPx_primaryDecay = std::min(fMinPx_primaryDecay, xMom);
452     fMaxPx_primaryDecay = std::max(fMaxPx_primaryDecay, xMom);
453     fMeanPy_primaryDecay += yMom;
454     fMinPy_primaryDecay = std::min(fMinPy_primaryDecay, yMom);
455     fMaxPy_primaryDecay = std::max(fMaxPy_primaryDecay, yMom);
456     fMeanPz_primaryDecay += zMom;
457     fMinPz_primaryDecay = std::min(fMinPz_primaryDecay, zMom);
458     fMaxPz_primaryDecay = std::max(fMaxPz_primaryDecay, zMom);
459 
460     //--- Extra checks ---
461     // Compute the "MC-truth" decay radius using the proper time of the primary particle when
462     // it decays.
463     // To do this, we would need an "effective" or "average" Lorentz beta and gamma of the primary
464     // particle during its lifetime, whereas in practice we have only the Lorentz beta and gamma
465     // values at the beginning and at the end when it decays. So, we can get only either an
466     // overestimate of the "MC-truth" decay radius - by using the initial Lorentz beta and gamma -
467     // or an underestimate of it - by using the Lorentz beta and gamma at the decay.
468     // We want to check the average values and the largest values of these wrong estimates.
469     const G4double underestimated_mc_truth_rPos =
470       t_proper_val * gamma_val * beta_val * CLHEP::c_light;
471     const G4double overestimated_mc_truth_rPos =
472       t_proper_val * fPrimaryParticleInitialGamma * fPrimaryParticleInitialBeta * CLHEP::c_light;
473     const G4double underestimated_mc_truth_rPos_delta =
474       underestimated_mc_truth_rPos - mc_truth_rPos_val;
475     const G4double overestimated_mc_truth_rPos_delta =
476       overestimated_mc_truth_rPos - mc_truth_rPos_val;
477     fMeanUnderestimated_mc_truth_rPos_delta += underestimated_mc_truth_rPos_delta;
478     fMeanOverestimated_mc_truth_rPos_delta += overestimated_mc_truth_rPos_delta;
479     // G4cout << "\t underestimated_mc_truth_rPos_delta [mum] = "
480     //        << underestimated_mc_truth_rPos_delta / CLHEP::micrometer
481     //        << " ; overestimated_mc_truth_rPos_delta [mum] = "
482     //        << overestimated_mc_truth_rPos_delta / CLHEP::micrometer << G4endl;
483     if (-underestimated_mc_truth_rPos_delta > ToleranceDeltaDecayRadius()) {
484       // G4cout << std::setprecision(6)
485       //        << " Large : underestimated_mc_truth_rPos_delta [mum]="
486       //        << underestimated_mc_truth_rPos_delta / CLHEP::micrometer
487       //        << " ; " << underestimated_mc_truth_rPos << " , "
488       //        << mc_truth_rPos_val << " mm" << G4endl;
489       if (fRunPtr) fRunPtr->IncrementNumber_underestimated_mc_truth_rPos_delta_above();
490     }
491     if (overestimated_mc_truth_rPos_delta > ToleranceDeltaDecayRadius()) {
492       // G4cout << std::setprecision(6)
493       //        << " Large : overestimated_mc_truth_rPos_delta [mum]="
494       //        << overestimated_mc_truth_rPos_delta / CLHEP::micrometer
495       //        << " ; " << overestimated_mc_truth_rPos << " , "
496       //        << mc_truth_rPos_val << " mm" << G4endl;
497       if (fRunPtr) fRunPtr->IncrementNumber_overestimated_mc_truth_rPos_delta_above();
498     }
499     const G4double underestimated_rDeltaPos = underestimated_mc_truth_rPos - rPos;
500     const G4double overestimated_rDeltaPos = overestimated_mc_truth_rPos - rPos;
501     fMeanUnderestimated_rDeltaPos += underestimated_rDeltaPos;
502     fMeanOverestimated_rDeltaPos += overestimated_rDeltaPos;
503     // G4cout << std::setprecision(6)
504     //        << "\t underestimated_rDeltaPos=" << underestimated_rDeltaPos/CLHEP::micrometer
505     //        << " ; overestimated_rDeltaPos=" << overestimated_rDeltaPos/CLHEP::micrometer
506     //        << " mum" << G4endl;
507     if (-underestimated_rDeltaPos > ToleranceDeltaDecayRadius()) {
508       if (fRunPtr) fRunPtr->IncrementNumberLargeUnderestimates();
509     }
510     if (overestimated_rDeltaPos > ToleranceDeltaDecayRadius()) {
511       if (fRunPtr) fRunPtr->IncrementNumberLargeOverestimates();
512     }
513     // Keep note of the biggest discrepancies
514     fMinUnderestimated_mc_truth_rPos_delta =
515       std::min(fMinUnderestimated_mc_truth_rPos_delta, underestimated_mc_truth_rPos_delta);
516     fMaxOverestimated_mc_truth_rPos_delta =
517       std::max(fMaxOverestimated_mc_truth_rPos_delta, overestimated_mc_truth_rPos_delta);
518     fMinUnderestimated_rDeltaPos = std::min(fMinUnderestimated_rDeltaPos, underestimated_rDeltaPos);
519     fMaxOverestimated_rDeltaPos = std::max(fMaxOverestimated_rDeltaPos, overestimated_rDeltaPos);
520     //--- End extra checks ---
521 
522     // Check numerical errors due to the use of  float  instead of  double : try out
523     // several, equivalent computations, taking the one with the largest numerical error.
524     const G4float float_xPos = static_cast<G4float>(theStep->GetPostStepPoint()->GetPosition().x());
525     const G4float float_yPos = static_cast<G4float>(theStep->GetPostStepPoint()->GetPosition().y());
526     const G4float float_zPos = static_cast<G4float>(theStep->GetPostStepPoint()->GetPosition().z());
527     const G4float float_rPos =
528       std::sqrt(float_xPos * float_xPos + float_yPos * float_yPos + float_zPos * float_zPos);
529     const G4float float_tPos = static_cast<G4float>(theStep->GetPostStepPoint()->GetLocalTime());
530     const G4float float_initialBeta1 = static_cast<G4float>(fPrimaryParticleInitialBeta);
531     const G4float float_initialBeta2 = static_cast<G4float>(fPrimaryParticleInitialMomentum)
532                                        / static_cast<G4float>(fPrimaryParticleInitialTotalEnergy);
533     const G4float float_initialGamma = static_cast<G4float>(fPrimaryParticleInitialGamma);
534     const G4float float_initialBeta3 =
535       std::sqrt(float_initialGamma * float_initialGamma - 1.0) / float_initialGamma;
536     const G4float float_c_light = static_cast<G4float>(CLHEP::c_light);
537     const G4float float_mc_truth_rPos1 = float_tPos * float_initialBeta1 * float_c_light;
538     const G4float float_mc_truth_rPos2 = float_tPos * float_initialBeta2 * float_c_light;
539     const G4float float_mc_truth_rPos3 = float_tPos * float_initialBeta3 * float_c_light;
540     const G4float float_rDeltaPos_0 = static_cast<G4float>(rDeltaPos);
541     const G4float float_rDeltaPos_1 = float_mc_truth_rPos1 - float_rPos;
542     const G4float float_rDeltaPos_2 = float_mc_truth_rPos2 - float_rPos;
543     const G4float float_rDeltaPos_3 = float_mc_truth_rPos3 - float_rPos;
544     const G4float float_rDeltaPos_4 = static_cast<G4float>(mc_truth_rPos) - float_rPos;
545     const G4float float_rDeltaPos_5 = float_mc_truth_rPos1 - static_cast<G4float>(rPos);
546     const G4float float_rDeltaPos_6 = float_mc_truth_rPos2 - static_cast<G4float>(rPos);
547     const G4float float_rDeltaPos_7 = float_mc_truth_rPos3 - static_cast<G4float>(rPos);
548     G4double rDeltaPos_deltaMax =
549       std::max(std::abs(float_rDeltaPos_0 - rDeltaPos), std::abs(float_rDeltaPos_1 - rDeltaPos));
550     rDeltaPos_deltaMax = std::max(rDeltaPos_deltaMax, std::abs(float_rDeltaPos_2 - rDeltaPos));
551     rDeltaPos_deltaMax = std::max(rDeltaPos_deltaMax, std::abs(float_rDeltaPos_3 - rDeltaPos));
552     rDeltaPos_deltaMax = std::max(rDeltaPos_deltaMax, std::abs(float_rDeltaPos_4 - rDeltaPos));
553     rDeltaPos_deltaMax = std::max(rDeltaPos_deltaMax, std::abs(float_rDeltaPos_5 - rDeltaPos));
554     rDeltaPos_deltaMax = std::max(rDeltaPos_deltaMax, std::abs(float_rDeltaPos_6 - rDeltaPos));
555     rDeltaPos_deltaMax = std::max(rDeltaPos_deltaMax, std::abs(float_rDeltaPos_7 - rDeltaPos));
556     // G4cout << std::setprecision(6) << " rDeltaPos_deltaMax[mum]="
557     //        << rDeltaPos_deltaMax / CLHEP::micrometer << G4endl;
558     fMaxFloat_rDeltaPos_deltaMax = std::max(fMaxFloat_rDeltaPos_deltaMax, rDeltaPos_deltaMax);
559 
560     // Get properties of the decay products and check the energy-momentum conservation of the
561     // decay
562     std::size_t nSec = theStep->GetNumberOfSecondariesInCurrentStep();
563     const std::vector<const G4Track*>* ptrVecSecondaries = theStep->GetSecondaryInCurrentStep();
564     G4double deltaE = 0.0, deltaPx = 0.0, deltaPy = 0.0, deltaPz = 0.0;
565     if (nSec > 0 && ptrVecSecondaries != nullptr) {
566       G4double sumEsecondaries = 0.0;
567       G4ThreeVector sumPsecondaries(0.0, 0.0, 0.0);
568       for (std::size_t i = 0; i < nSec; ++i) {
569         if ((*ptrVecSecondaries)[i]) {
570           sumEsecondaries += (*ptrVecSecondaries)[i]->GetTotalEnergy();
571           sumPsecondaries += (*ptrVecSecondaries)[i]->GetMomentum();
572         }
573       }
574       deltaE = sumEsecondaries - theStep->GetPostStepPoint()->GetTotalEnergy();
575       fMeanViolationE_primaryDecay += deltaE;
576       fMinViolationE_primaryDecay = std::min(fMinViolationE_primaryDecay, deltaE);
577       fMaxViolationE_primaryDecay = std::max(fMaxViolationE_primaryDecay, deltaE);
578       if (std::abs(deltaE) > ToleranceEPviolations()) {
579         if (fRunPtr) fRunPtr->IncrementNumberEviolations();
580       }
581       deltaPx = sumPsecondaries.x() - xMom;
582       fMeanViolationPx_primaryDecay += deltaPx;
583       fMinViolationPx_primaryDecay = std::min(fMinViolationPx_primaryDecay, deltaPx);
584       fMaxViolationPx_primaryDecay = std::max(fMaxViolationPx_primaryDecay, deltaPx);
585       deltaPy = sumPsecondaries.y() - yMom;
586       fMeanViolationPy_primaryDecay += deltaPy;
587       fMinViolationPy_primaryDecay = std::min(fMinViolationPy_primaryDecay, deltaPy);
588       fMaxViolationPy_primaryDecay = std::max(fMaxViolationPy_primaryDecay, deltaPy);
589       deltaPz = sumPsecondaries.z() - zMom;
590       fMeanViolationPz_primaryDecay += deltaPz;
591       fMinViolationPz_primaryDecay = std::min(fMinViolationPz_primaryDecay, deltaPz);
592       fMaxViolationPz_primaryDecay = std::max(fMaxViolationPz_primaryDecay, deltaPz);
593       if (std::abs(deltaPx) > ToleranceEPviolations() || std::abs(deltaPy) > ToleranceEPviolations()
594           || std::abs(deltaPz) > ToleranceEPviolations())
595       {
596         if (fRunPtr) fRunPtr->IncrementNumberPviolations();
597       }
598     }
599     else {
600       if (fRunPtr) fRunPtr->IncrementNumberBadPrimaryDecays();
601     }
602 
603     if (fRunPtr) {
604       fRunPtr->IncrementNumberDecays();
605       fRunPtr->SetDecayT(tPos);
606       fRunPtr->SetDecayR_mc_truth(mc_truth_rPos);
607       fRunPtr->SetDecayR(rPos);
608       fRunPtr->SetDecayX(xPos);
609       fRunPtr->SetDecayY(yPos);
610       fRunPtr->SetDecayZ(zPos);
611       fRunPtr->SetDeltaDecayR(rDeltaPos);
612       fRunPtr->SetDeflectionAngle(deflection_angle_in_degrees);
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_deltaMax);
623       fRunPtr->SetMaxEtot_deltaMax(etot_deltaMax);
624       fRunPtr->SetMaxP_deltaMax(p_deltaMax);
625       fRunPtr->SetMaxPdir_deltaMax(pdir_deltaMax);
626       fRunPtr->SetMaxMass_deltaMax1(mass_deltaMax1);
627       fRunPtr->SetMaxMass_deltaMax2(mass_deltaMax2);
628       fRunPtr->SetMaxMass_deltaMax3(mass_deltaMax3);
629       fRunPtr->SetMaxBeta_deltaMax1(beta_deltaMax1);
630       fRunPtr->SetMaxBeta_deltaMax2(beta_deltaMax2);
631       fRunPtr->SetMaxGamma_deltaMax1(gamma_deltaMax1);
632       fRunPtr->SetMaxGamma_deltaMax2(gamma_deltaMax2);
633       fRunPtr->SetMaxGamma_deltaMax3(gamma_deltaMax3);
634       fRunPtr->SetMaxT_proper_deltaMax(t_proper_deltaMax);
635       fRunPtr->SetMaxT_lab_deltaMax(t_lab_deltaMax);
636       fRunPtr->SetMaxMc_truth_rPos_deltaMax(mc_truth_rPos_deltaMax);
637       fRunPtr->SetMinUnderestimated_mc_truth_rPos_delta(underestimated_mc_truth_rPos_delta);
638       fRunPtr->SetMaxOverestimated_mc_truth_rPos_delta(overestimated_mc_truth_rPos_delta);
639       fRunPtr->SetMinUnderestimated_rDeltaPos(underestimated_rDeltaPos);
640       fRunPtr->SetMaxOverestimated_rDeltaPos(overestimated_rDeltaPos);
641       fRunPtr->SetMaxFloat_rDeltaPos_deltaMax(fMaxFloat_rDeltaPos_deltaMax);
642     }
643   }
644 }
645 
646 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
647