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 ]

Diff markup

Differences between /examples/extended/hadronic/Hadr10/src/SteppingAction.cc (Version 11.3.0) and /examples/extended/hadronic/Hadr10/src/SteppingAction.cc (Version 11.1.2)


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