Geant4 Cross Reference |
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 const G4double deflection_angle_in_degrees = 402 if (xDirection < -1.0) xDirection = -1.0; << 401 57.29*std::acos( theStep->GetPostStepPoint()->GetMomentumDirection().x() ); 403 const G4double deflection_angle_in_degrees << 404 const G4double delta_ekin = fPrimaryPartic 402 const G4double delta_ekin = fPrimaryParticleInitialKineticEnergy - eKin; 405 // G4cout << std::setprecision(6) << 403 //G4cout << std::setprecision(6) 406 // << " Decay: tPos[ns]=" << tPos < << 404 // << " Decay: tPos[ns]=" << tPos << " ; rPos[mm]=" << rPos << " ; deltaR[mum]=" 407 // << rDeltaPos /CLHEP::micrometer << 405 // << rDeltaPos /CLHEP::micrometer << " ; deltaEkin[MeV]=" << delta_ekin 408 // << " ; deltaAngle(deg)=" << defl << 406 // << " ; deltaAngle(deg)=" << deflection_angle_in_degrees << G4endl; 409 // If the absolute difference between the << 407 // 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 << 408 // given threshold, then we notify this special situation in the output, with "LARGE_DELTA_R" 411 // for post-processing evaluation. Moreov << 409 // for post-processing evaluation. Moreover, in this case, if the "MC-truth" decay radius is 412 // smaller than the real one, then we cou << 410 // smaller than the real one, then we count this unexpected occurrence and we further notify 413 // this special situation in the output w << 411 // this special situation in the output with "***UNEXPECTED***" for post-processing 414 // evaluation. << 412 // evaluation. 415 if (std::abs(rDeltaPos) > ToleranceDeltaDe << 413 if ( std::abs( rDeltaPos ) > ToleranceDeltaDecayRadius() ) { 416 // G4cout << "\t LARGE_DELTA_R : mc_trut << 414 //G4cout << "\t LARGE_DELTA_R : mc_truth_rPos[mm]=" << mc_truth_rPos 417 // << " ; rPos[mm]=" << rPos; << 415 // << " ; rPos[mm]=" << rPos; 418 if (rDeltaPos < 0.0) { << 416 if ( rDeltaPos < 0.0 ) { 419 // G4cout << "\t ***UNEXPECTED***"; << 417 //G4cout << "\t ***UNEXPECTED***"; 420 if (fRunPtr) fRunPtr->IncrementNumberU << 418 if ( fRunPtr ) fRunPtr->IncrementNumberUnexpectedDecays(); 421 } 419 } 422 // G4cout << G4endl; << 420 //G4cout << G4endl; 423 } 421 } 424 fMeanDeltaR_primaryDecay += rDeltaPos; 422 fMeanDeltaR_primaryDecay += rDeltaPos; 425 fMinDeltaR_primaryDecay = std::min(fMinDel << 423 fMinDeltaR_primaryDecay = std::min( fMinDeltaR_primaryDecay, rDeltaPos ); 426 fMaxDeltaR_primaryDecay = std::max(fMaxDel << 424 fMaxDeltaR_primaryDecay = std::max( fMaxDeltaR_primaryDecay, rDeltaPos ); 427 fMeanR_primaryDecay += rPos; 425 fMeanR_primaryDecay += rPos; 428 fMinR_primaryDecay = std::min(fMinR_primar << 426 fMinR_primaryDecay = std::min( fMinR_primaryDecay, rPos ); 429 fMaxR_primaryDecay = std::max(fMaxR_primar << 427 fMaxR_primaryDecay = std::max( fMaxR_primaryDecay, rPos ); 430 fMeanX_primaryDecay += xPos; 428 fMeanX_primaryDecay += xPos; 431 fMinX_primaryDecay = std::min(fMinX_primar << 429 fMinX_primaryDecay = std::min( fMinX_primaryDecay, xPos ); 432 fMaxX_primaryDecay = std::max(fMaxX_primar << 430 fMaxX_primaryDecay = std::max( fMaxX_primaryDecay, xPos ); 433 fMeanY_primaryDecay += yPos; 431 fMeanY_primaryDecay += yPos; 434 fMinY_primaryDecay = std::min(fMinY_primar << 432 fMinY_primaryDecay = std::min( fMinY_primaryDecay, yPos ); 435 fMaxY_primaryDecay = std::max(fMaxY_primar << 433 fMaxY_primaryDecay = std::max( fMaxY_primaryDecay, yPos ); 436 fMeanZ_primaryDecay += zPos; 434 fMeanZ_primaryDecay += zPos; 437 fMinZ_primaryDecay = std::min(fMinZ_primar << 435 fMinZ_primaryDecay = std::min( fMinZ_primaryDecay, zPos ); 438 fMaxZ_primaryDecay = std::max(fMaxZ_primar << 436 fMaxZ_primaryDecay = std::max( fMaxZ_primaryDecay, zPos ); 439 fMeanDeltaAngle_primaryDecay += deflection 437 fMeanDeltaAngle_primaryDecay += deflection_angle_in_degrees; 440 fMinDeltaAngle_primaryDecay = << 438 fMinDeltaAngle_primaryDecay = std::min( fMinDeltaAngle_primaryDecay, 441 std::min(fMinDeltaAngle_primaryDecay, de << 439 deflection_angle_in_degrees ); 442 fMaxDeltaAngle_primaryDecay = << 440 fMaxDeltaAngle_primaryDecay = std::max( fMaxDeltaAngle_primaryDecay, 443 std::max(fMaxDeltaAngle_primaryDecay, de << 441 deflection_angle_in_degrees ); 444 fMeanDeltaEkin_primaryDecay += delta_ekin; 442 fMeanDeltaEkin_primaryDecay += delta_ekin; 445 fMinDeltaEkin_primaryDecay = std::min(fMin << 443 fMinDeltaEkin_primaryDecay = std::min( fMinDeltaEkin_primaryDecay, delta_ekin ); 446 fMaxDeltaEkin_primaryDecay = std::max(fMax << 444 fMaxDeltaEkin_primaryDecay = std::max( fMaxDeltaEkin_primaryDecay, delta_ekin ); 447 fMeanEkin_primaryDecay += eKin; 445 fMeanEkin_primaryDecay += eKin; 448 fMinEkin_primaryDecay = std::min(fMinEkin_ << 446 fMinEkin_primaryDecay = std::min( fMinEkin_primaryDecay, eKin ); 449 fMaxEkin_primaryDecay = std::max(fMaxEkin_ << 447 fMaxEkin_primaryDecay = std::max( fMaxEkin_primaryDecay, eKin ); 450 fMeanPx_primaryDecay += xMom; 448 fMeanPx_primaryDecay += xMom; 451 fMinPx_primaryDecay = std::min(fMinPx_prim << 449 fMinPx_primaryDecay = std::min( fMinPx_primaryDecay, xMom ); 452 fMaxPx_primaryDecay = std::max(fMaxPx_prim << 450 fMaxPx_primaryDecay = std::max( fMaxPx_primaryDecay, xMom ); 453 fMeanPy_primaryDecay += yMom; 451 fMeanPy_primaryDecay += yMom; 454 fMinPy_primaryDecay = std::min(fMinPy_prim << 452 fMinPy_primaryDecay = std::min( fMinPy_primaryDecay, yMom ); 455 fMaxPy_primaryDecay = std::max(fMaxPy_prim << 453 fMaxPy_primaryDecay = std::max( fMaxPy_primaryDecay, yMom ); 456 fMeanPz_primaryDecay += zMom; 454 fMeanPz_primaryDecay += zMom; 457 fMinPz_primaryDecay = std::min(fMinPz_prim << 455 fMinPz_primaryDecay = std::min( fMinPz_primaryDecay, zMom ); 458 fMaxPz_primaryDecay = std::max(fMaxPz_prim << 456 fMaxPz_primaryDecay = std::max( fMaxPz_primaryDecay, zMom ); 459 457 460 //--- Extra checks --- 458 //--- Extra checks --- 461 // Compute the "MC-truth" decay radius usi 459 // Compute the "MC-truth" decay radius using the proper time of the primary particle when 462 // it decays. 460 // it decays. 463 // To do this, we would need an "effective 461 // To do this, we would need an "effective" or "average" Lorentz beta and gamma of the primary 464 // particle during its lifetime, whereas i 462 // particle during its lifetime, whereas in practice we have only the Lorentz beta and gamma 465 // values at the beginning and at the end 463 // 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 464 // overestimate of the "MC-truth" decay radius - by using the initial Lorentz beta and gamma - 467 // or an underestimate of it - by using th 465 // or an underestimate of it - by using the Lorentz beta and gamma at the decay. 468 // We want to check the average values and 466 // We want to check the average values and the largest values of these wrong estimates. 469 const G4double underestimated_mc_truth_rPo 467 const G4double underestimated_mc_truth_rPos = 470 t_proper_val * gamma_val * beta_val * CL 468 t_proper_val * gamma_val * beta_val * CLHEP::c_light; 471 const G4double overestimated_mc_truth_rPos 469 const G4double overestimated_mc_truth_rPos = 472 t_proper_val * fPrimaryParticleInitialGa 470 t_proper_val * fPrimaryParticleInitialGamma * fPrimaryParticleInitialBeta * CLHEP::c_light; 473 const G4double underestimated_mc_truth_rPo 471 const G4double underestimated_mc_truth_rPos_delta = 474 underestimated_mc_truth_rPos - mc_truth_ 472 underestimated_mc_truth_rPos - mc_truth_rPos_val; 475 const G4double overestimated_mc_truth_rPos 473 const G4double overestimated_mc_truth_rPos_delta = 476 overestimated_mc_truth_rPos - mc_truth_r << 474 overestimated_mc_truth_rPos - mc_truth_rPos_val; 477 fMeanUnderestimated_mc_truth_rPos_delta += 475 fMeanUnderestimated_mc_truth_rPos_delta += underestimated_mc_truth_rPos_delta; 478 fMeanOverestimated_mc_truth_rPos_delta += << 476 fMeanOverestimated_mc_truth_rPos_delta += overestimated_mc_truth_rPos_delta; 479 // G4cout << "\t underestimated_mc_truth_r << 477 //G4cout << "\t underestimated_mc_truth_rPos_delta [mum] = " 480 // << underestimated_mc_truth_rPos_ << 478 // << underestimated_mc_truth_rPos_delta / CLHEP::micrometer 481 // << " ; overestimated_mc_truth_rP << 479 // << " ; overestimated_mc_truth_rPos_delta [mum] = " 482 // << overestimated_mc_truth_rPos_d << 480 // << overestimated_mc_truth_rPos_delta / CLHEP::micrometer << G4endl; 483 if (-underestimated_mc_truth_rPos_delta > << 481 if ( -underestimated_mc_truth_rPos_delta > ToleranceDeltaDecayRadius() ) { 484 // G4cout << std::setprecision(6) << 482 //G4cout << std::setprecision(6) 485 // << " Large : underestimated_mc << 483 // << " Large : underestimated_mc_truth_rPos_delta [mum]=" 486 // << underestimated_mc_truth_rPo << 484 // << underestimated_mc_truth_rPos_delta / CLHEP::micrometer 487 // << " ; " << underestimated_mc_ << 485 // << " ; " << underestimated_mc_truth_rPos << " , " 488 // << mc_truth_rPos_val << " mm" << 486 // << mc_truth_rPos_val << " mm" << G4endl; 489 if (fRunPtr) fRunPtr->IncrementNumber_un << 487 if ( fRunPtr ) fRunPtr->IncrementNumber_underestimated_mc_truth_rPos_delta_above(); 490 } << 488 } 491 if (overestimated_mc_truth_rPos_delta > To << 489 if ( overestimated_mc_truth_rPos_delta > ToleranceDeltaDecayRadius() ) { 492 // G4cout << std::setprecision(6) << 490 //G4cout << std::setprecision(6) 493 // << " Large : overestimated_mc_ << 491 // << " Large : overestimated_mc_truth_rPos_delta [mum]=" 494 // << overestimated_mc_truth_rPos << 492 // << overestimated_mc_truth_rPos_delta / CLHEP::micrometer 495 // << " ; " << overestimated_mc_t << 493 // << " ; " << overestimated_mc_truth_rPos << " , " 496 // << mc_truth_rPos_val << " mm" << 494 // << mc_truth_rPos_val << " mm" << G4endl; 497 if (fRunPtr) fRunPtr->IncrementNumber_ov << 495 if ( fRunPtr ) fRunPtr->IncrementNumber_overestimated_mc_truth_rPos_delta_above(); 498 } << 496 } 499 const G4double underestimated_rDeltaPos = 497 const G4double underestimated_rDeltaPos = underestimated_mc_truth_rPos - rPos; 500 const G4double overestimated_rDeltaPos = o << 498 const G4double overestimated_rDeltaPos = overestimated_mc_truth_rPos - rPos; 501 fMeanUnderestimated_rDeltaPos += underesti 499 fMeanUnderestimated_rDeltaPos += underestimated_rDeltaPos; 502 fMeanOverestimated_rDeltaPos += overestima << 500 fMeanOverestimated_rDeltaPos += overestimated_rDeltaPos; 503 // G4cout << std::setprecision(6) << 501 //G4cout << std::setprecision(6) 504 // << "\t underestimated_rDeltaPos= << 502 // << "\t underestimated_rDeltaPos=" << underestimated_rDeltaPos/CLHEP::micrometer 505 // << " ; overestimated_rDeltaPos=" << 503 // << " ; overestimated_rDeltaPos=" << overestimated_rDeltaPos/CLHEP::micrometer 506 // << " mum" << G4endl; << 504 // << " mum" << G4endl; 507 if (-underestimated_rDeltaPos > ToleranceD << 505 if ( -underestimated_rDeltaPos > ToleranceDeltaDecayRadius() ) { 508 if (fRunPtr) fRunPtr->IncrementNumberLar << 506 if ( fRunPtr ) fRunPtr->IncrementNumberLargeUnderestimates(); 509 } 507 } 510 if (overestimated_rDeltaPos > ToleranceDel << 508 if ( overestimated_rDeltaPos > ToleranceDeltaDecayRadius() ) { 511 if (fRunPtr) fRunPtr->IncrementNumberLar << 509 if ( fRunPtr ) fRunPtr->IncrementNumberLargeOverestimates(); 512 } 510 } 513 // Keep note of the biggest discrepancies 511 // Keep note of the biggest discrepancies 514 fMinUnderestimated_mc_truth_rPos_delta = << 512 fMinUnderestimated_mc_truth_rPos_delta = std::min( fMinUnderestimated_mc_truth_rPos_delta, 515 std::min(fMinUnderestimated_mc_truth_rPo << 513 underestimated_mc_truth_rPos_delta ); 516 fMaxOverestimated_mc_truth_rPos_delta = << 514 fMaxOverestimated_mc_truth_rPos_delta = std::max( fMaxOverestimated_mc_truth_rPos_delta, 517 std::max(fMaxOverestimated_mc_truth_rPos << 515 overestimated_mc_truth_rPos_delta ); 518 fMinUnderestimated_rDeltaPos = std::min(fM << 516 fMinUnderestimated_rDeltaPos = std::min( fMinUnderestimated_rDeltaPos, 519 fMaxOverestimated_rDeltaPos = std::max(fMa << 517 underestimated_rDeltaPos ); >> 518 fMaxOverestimated_rDeltaPos = std::max( fMaxOverestimated_rDeltaPos, >> 519 overestimated_rDeltaPos ); 520 //--- End extra checks --- 520 //--- End extra checks --- 521 521 522 // Check numerical errors due to the use o 522 // Check numerical errors due to the use of float instead of double : try out 523 // several, equivalent computations, takin 523 // several, equivalent computations, taking the one with the largest numerical error. 524 const G4float float_xPos = static_cast<G4f << 524 const G4float float_xPos = 525 const G4float float_yPos = static_cast<G4f << 525 static_cast< G4float >( theStep->GetPostStepPoint()->GetPosition().x() ); 526 const G4float float_zPos = static_cast<G4f << 526 const G4float float_yPos = >> 527 static_cast< G4float >( theStep->GetPostStepPoint()->GetPosition().y() ); >> 528 const G4float float_zPos = >> 529 static_cast< G4float >( theStep->GetPostStepPoint()->GetPosition().z() ); 527 const G4float float_rPos = 530 const G4float float_rPos = 528 std::sqrt(float_xPos * float_xPos + floa << 531 std::sqrt( float_xPos*float_xPos + float_yPos*float_yPos + float_zPos*float_zPos ); 529 const G4float float_tPos = static_cast<G4f << 532 const G4float float_tPos = 530 const G4float float_initialBeta1 = static_ << 533 static_cast< G4float >( theStep->GetPostStepPoint()->GetLocalTime() ); 531 const G4float float_initialBeta2 = static_ << 534 const G4float float_initialBeta1 = static_cast< G4float >( fPrimaryParticleInitialBeta ); 532 / stati << 535 const G4float float_initialBeta2 = 533 const G4float float_initialGamma = static_ << 536 static_cast< G4float >( fPrimaryParticleInitialMomentum ) / >> 537 static_cast< G4float >( fPrimaryParticleInitialTotalEnergy ); >> 538 const G4float float_initialGamma = static_cast< G4float >( fPrimaryParticleInitialGamma ); 534 const G4float float_initialBeta3 = 539 const G4float float_initialBeta3 = 535 std::sqrt(float_initialGamma * float_ini << 540 std::sqrt( float_initialGamma*float_initialGamma - 1.0 ) / float_initialGamma; 536 const G4float float_c_light = static_cast< << 541 const G4float float_c_light = static_cast< G4float >( CLHEP::c_light ); 537 const G4float float_mc_truth_rPos1 = float 542 const G4float float_mc_truth_rPos1 = float_tPos * float_initialBeta1 * float_c_light; 538 const G4float float_mc_truth_rPos2 = float 543 const G4float float_mc_truth_rPos2 = float_tPos * float_initialBeta2 * float_c_light; 539 const G4float float_mc_truth_rPos3 = float 544 const G4float float_mc_truth_rPos3 = float_tPos * float_initialBeta3 * float_c_light; 540 const G4float float_rDeltaPos_0 = static_c << 545 const G4float float_rDeltaPos_0 = static_cast< G4float >( rDeltaPos ); 541 const G4float float_rDeltaPos_1 = float_mc 546 const G4float float_rDeltaPos_1 = float_mc_truth_rPos1 - float_rPos; 542 const G4float float_rDeltaPos_2 = float_mc 547 const G4float float_rDeltaPos_2 = float_mc_truth_rPos2 - float_rPos; 543 const G4float float_rDeltaPos_3 = float_mc 548 const G4float float_rDeltaPos_3 = float_mc_truth_rPos3 - float_rPos; 544 const G4float float_rDeltaPos_4 = static_c << 549 const G4float float_rDeltaPos_4 = static_cast< G4float >( mc_truth_rPos ) - float_rPos; 545 const G4float float_rDeltaPos_5 = float_mc << 550 const G4float float_rDeltaPos_5 = float_mc_truth_rPos1 - static_cast< G4float >( rPos ); 546 const G4float float_rDeltaPos_6 = float_mc << 551 const G4float float_rDeltaPos_6 = float_mc_truth_rPos2 - static_cast< G4float >( rPos ); 547 const G4float float_rDeltaPos_7 = float_mc << 552 const G4float float_rDeltaPos_7 = float_mc_truth_rPos3 - static_cast< G4float >( rPos ); 548 G4double rDeltaPos_deltaMax = << 553 G4double rDeltaPos_deltaMax = std::max( std::abs( float_rDeltaPos_0 - rDeltaPos ), 549 std::max(std::abs(float_rDeltaPos_0 - rD << 554 std::abs( float_rDeltaPos_1 - rDeltaPos ) ); 550 rDeltaPos_deltaMax = std::max(rDeltaPos_de << 555 rDeltaPos_deltaMax = std::max( rDeltaPos_deltaMax, 551 rDeltaPos_deltaMax = std::max(rDeltaPos_de << 556 std::abs( float_rDeltaPos_2 - rDeltaPos ) ); 552 rDeltaPos_deltaMax = std::max(rDeltaPos_de << 557 rDeltaPos_deltaMax = std::max( rDeltaPos_deltaMax, 553 rDeltaPos_deltaMax = std::max(rDeltaPos_de << 558 std::abs( float_rDeltaPos_3 - rDeltaPos ) ); 554 rDeltaPos_deltaMax = std::max(rDeltaPos_de << 559 rDeltaPos_deltaMax = std::max( rDeltaPos_deltaMax, 555 rDeltaPos_deltaMax = std::max(rDeltaPos_de << 560 std::abs( float_rDeltaPos_4 - rDeltaPos ) ); 556 // G4cout << std::setprecision(6) << " rDe << 561 rDeltaPos_deltaMax = std::max( rDeltaPos_deltaMax, 557 // << rDeltaPos_deltaMax / CLHEP::m << 562 std::abs( float_rDeltaPos_5 - rDeltaPos ) ); 558 fMaxFloat_rDeltaPos_deltaMax = std::max(fM << 563 rDeltaPos_deltaMax = std::max( rDeltaPos_deltaMax, 559 << 564 std::abs( float_rDeltaPos_6 - rDeltaPos ) ); >> 565 rDeltaPos_deltaMax = std::max( rDeltaPos_deltaMax, >> 566 std::abs( float_rDeltaPos_7 - rDeltaPos ) ); >> 567 //G4cout << std::setprecision(6) << " rDeltaPos_deltaMax[mum]=" >> 568 // << rDeltaPos_deltaMax / CLHEP::micrometer << G4endl; >> 569 fMaxFloat_rDeltaPos_deltaMax = std::max( fMaxFloat_rDeltaPos_deltaMax, rDeltaPos_deltaMax ); >> 570 560 // Get properties of the decay products an 571 // Get properties of the decay products and check the energy-momentum conservation of the 561 // decay 572 // decay 562 std::size_t nSec = theStep->GetNumberOfSec 573 std::size_t nSec = theStep->GetNumberOfSecondariesInCurrentStep(); 563 const std::vector<const G4Track*>* ptrVecS << 574 const std::vector< const G4Track* >* ptrVecSecondaries = theStep->GetSecondaryInCurrentStep(); 564 G4double deltaE = 0.0, deltaPx = 0.0, delt 575 G4double deltaE = 0.0, deltaPx = 0.0, deltaPy = 0.0, deltaPz = 0.0; 565 if (nSec > 0 && ptrVecSecondaries != nullp << 576 if ( nSec > 0 && ptrVecSecondaries != nullptr ) { 566 G4double sumEsecondaries = 0.0; 577 G4double sumEsecondaries = 0.0; 567 G4ThreeVector sumPsecondaries(0.0, 0.0, << 578 G4ThreeVector sumPsecondaries( 0.0, 0.0, 0.0 ); 568 for (std::size_t i = 0; i < nSec; ++i) { << 579 for ( std::size_t i = 0; i < nSec; ++i ) { 569 if ((*ptrVecSecondaries)[i]) { << 580 if ( (*ptrVecSecondaries)[i] ) { 570 sumEsecondaries += (*ptrVecSecondari 581 sumEsecondaries += (*ptrVecSecondaries)[i]->GetTotalEnergy(); 571 sumPsecondaries += (*ptrVecSecondari 582 sumPsecondaries += (*ptrVecSecondaries)[i]->GetMomentum(); 572 } 583 } 573 } 584 } 574 deltaE = sumEsecondaries - theStep->GetP 585 deltaE = sumEsecondaries - theStep->GetPostStepPoint()->GetTotalEnergy(); 575 fMeanViolationE_primaryDecay += deltaE; 586 fMeanViolationE_primaryDecay += deltaE; 576 fMinViolationE_primaryDecay = std::min(f << 587 fMinViolationE_primaryDecay = std::min( fMinViolationE_primaryDecay, deltaE ); 577 fMaxViolationE_primaryDecay = std::max(f << 588 fMaxViolationE_primaryDecay = std::max( fMaxViolationE_primaryDecay, deltaE ); 578 if (std::abs(deltaE) > ToleranceEPviolat << 589 if ( std::abs( deltaE ) > ToleranceEPviolations() ) { 579 if (fRunPtr) fRunPtr->IncrementNumberE << 590 if ( fRunPtr ) fRunPtr->IncrementNumberEviolations(); 580 } 591 } 581 deltaPx = sumPsecondaries.x() - xMom; 592 deltaPx = sumPsecondaries.x() - xMom; 582 fMeanViolationPx_primaryDecay += deltaPx 593 fMeanViolationPx_primaryDecay += deltaPx; 583 fMinViolationPx_primaryDecay = std::min( << 594 fMinViolationPx_primaryDecay = std::min( fMinViolationPx_primaryDecay, deltaPx ); 584 fMaxViolationPx_primaryDecay = std::max( << 595 fMaxViolationPx_primaryDecay = std::max( fMaxViolationPx_primaryDecay, deltaPx ); 585 deltaPy = sumPsecondaries.y() - yMom; 596 deltaPy = sumPsecondaries.y() - yMom; 586 fMeanViolationPy_primaryDecay += deltaPy 597 fMeanViolationPy_primaryDecay += deltaPy; 587 fMinViolationPy_primaryDecay = std::min( << 598 fMinViolationPy_primaryDecay = std::min( fMinViolationPy_primaryDecay, deltaPy ); 588 fMaxViolationPy_primaryDecay = std::max( << 599 fMaxViolationPy_primaryDecay = std::max( fMaxViolationPy_primaryDecay, deltaPy ); 589 deltaPz = sumPsecondaries.z() - zMom; 600 deltaPz = sumPsecondaries.z() - zMom; 590 fMeanViolationPz_primaryDecay += deltaPz 601 fMeanViolationPz_primaryDecay += deltaPz; 591 fMinViolationPz_primaryDecay = std::min( << 602 fMinViolationPz_primaryDecay = std::min( fMinViolationPz_primaryDecay, deltaPz ); 592 fMaxViolationPz_primaryDecay = std::max( << 603 fMaxViolationPz_primaryDecay = std::max( fMaxViolationPz_primaryDecay, deltaPz ); 593 if (std::abs(deltaPx) > ToleranceEPviola << 604 if ( std::abs( deltaPx ) > ToleranceEPviolations() || 594 || std::abs(deltaPz) > ToleranceEPvi << 605 std::abs( deltaPy ) > ToleranceEPviolations() || 595 { << 606 std::abs( deltaPz ) > ToleranceEPviolations() ) { 596 if (fRunPtr) fRunPtr->IncrementNumberP << 607 if ( fRunPtr ) fRunPtr->IncrementNumberPviolations(); 597 } 608 } 598 } << 609 } else { 599 else { << 610 if ( fRunPtr ) fRunPtr->IncrementNumberBadPrimaryDecays(); 600 if (fRunPtr) fRunPtr->IncrementNumberBad << 601 } 611 } 602 612 603 if (fRunPtr) { << 613 if ( fRunPtr ) { 604 fRunPtr->IncrementNumberDecays(); 614 fRunPtr->IncrementNumberDecays(); 605 fRunPtr->SetDecayT(tPos); << 615 fRunPtr->SetDecayT( tPos ); 606 fRunPtr->SetDecayR_mc_truth(mc_truth_rPo << 616 fRunPtr->SetDecayR_mc_truth( mc_truth_rPos ); 607 fRunPtr->SetDecayR(rPos); << 617 fRunPtr->SetDecayR( rPos ); 608 fRunPtr->SetDecayX(xPos); << 618 fRunPtr->SetDecayX( xPos ); 609 fRunPtr->SetDecayY(yPos); << 619 fRunPtr->SetDecayY( yPos ); 610 fRunPtr->SetDecayZ(zPos); << 620 fRunPtr->SetDecayZ( zPos ); 611 fRunPtr->SetDeltaDecayR(rDeltaPos); << 621 fRunPtr->SetDeltaDecayR( rDeltaPos ); 612 fRunPtr->SetDeflectionAngle(deflection_a << 622 fRunPtr->SetDeflectionAngle( deflection_angle_in_degrees ); 613 fRunPtr->SetDeltaEkin(delta_ekin); << 623 fRunPtr->SetDeltaEkin( delta_ekin ); 614 fRunPtr->SetDecayEkin(eKin); << 624 fRunPtr->SetDecayEkin( eKin ); 615 fRunPtr->SetDecayPx(xMom); << 625 fRunPtr->SetDecayPx( xMom ); 616 fRunPtr->SetDecayPy(yMom); << 626 fRunPtr->SetDecayPy( yMom ); 617 fRunPtr->SetDecayPz(zMom); << 627 fRunPtr->SetDecayPz( zMom ); 618 fRunPtr->SetDecayEtotViolation(deltaE); << 628 fRunPtr->SetDecayEtotViolation( deltaE ); 619 fRunPtr->SetDecayPxViolation(deltaPx); << 629 fRunPtr->SetDecayPxViolation( deltaPx ); 620 fRunPtr->SetDecayPyViolation(deltaPy); << 630 fRunPtr->SetDecayPyViolation( deltaPy ); 621 fRunPtr->SetDecayPzViolation(deltaPz); << 631 fRunPtr->SetDecayPzViolation( deltaPz ); 622 fRunPtr->SetMaxEkin_deltaMax(ekin_deltaM << 632 fRunPtr->SetMaxEkin_deltaMax( ekin_deltaMax ); 623 fRunPtr->SetMaxEtot_deltaMax(etot_deltaM << 633 fRunPtr->SetMaxEtot_deltaMax( etot_deltaMax ); 624 fRunPtr->SetMaxP_deltaMax(p_deltaMax); << 634 fRunPtr->SetMaxP_deltaMax( p_deltaMax ); 625 fRunPtr->SetMaxPdir_deltaMax(pdir_deltaM << 635 fRunPtr->SetMaxPdir_deltaMax( pdir_deltaMax ); 626 fRunPtr->SetMaxMass_deltaMax1(mass_delta << 636 fRunPtr->SetMaxMass_deltaMax1( mass_deltaMax1 ); 627 fRunPtr->SetMaxMass_deltaMax2(mass_delta << 637 fRunPtr->SetMaxMass_deltaMax2( mass_deltaMax2 ); 628 fRunPtr->SetMaxMass_deltaMax3(mass_delta << 638 fRunPtr->SetMaxMass_deltaMax3( mass_deltaMax3 ); 629 fRunPtr->SetMaxBeta_deltaMax1(beta_delta << 639 fRunPtr->SetMaxBeta_deltaMax1( beta_deltaMax1 ); 630 fRunPtr->SetMaxBeta_deltaMax2(beta_delta << 640 fRunPtr->SetMaxBeta_deltaMax2( beta_deltaMax2 ); 631 fRunPtr->SetMaxGamma_deltaMax1(gamma_del << 641 fRunPtr->SetMaxGamma_deltaMax1( gamma_deltaMax1 ); 632 fRunPtr->SetMaxGamma_deltaMax2(gamma_del << 642 fRunPtr->SetMaxGamma_deltaMax2( gamma_deltaMax2 ); 633 fRunPtr->SetMaxGamma_deltaMax3(gamma_del << 643 fRunPtr->SetMaxGamma_deltaMax3( gamma_deltaMax3 ); 634 fRunPtr->SetMaxT_proper_deltaMax(t_prope << 644 fRunPtr->SetMaxT_proper_deltaMax( t_proper_deltaMax ); 635 fRunPtr->SetMaxT_lab_deltaMax(t_lab_delt << 645 fRunPtr->SetMaxT_lab_deltaMax( t_lab_deltaMax ); 636 fRunPtr->SetMaxMc_truth_rPos_deltaMax(mc << 646 fRunPtr->SetMaxMc_truth_rPos_deltaMax( mc_truth_rPos_deltaMax ); 637 fRunPtr->SetMinUnderestimated_mc_truth_r << 647 fRunPtr->SetMinUnderestimated_mc_truth_rPos_delta( underestimated_mc_truth_rPos_delta ); 638 fRunPtr->SetMaxOverestimated_mc_truth_rP << 648 fRunPtr->SetMaxOverestimated_mc_truth_rPos_delta( overestimated_mc_truth_rPos_delta ); 639 fRunPtr->SetMinUnderestimated_rDeltaPos( << 649 fRunPtr->SetMinUnderestimated_rDeltaPos( underestimated_rDeltaPos ); 640 fRunPtr->SetMaxOverestimated_rDeltaPos(o << 650 fRunPtr->SetMaxOverestimated_rDeltaPos( overestimated_rDeltaPos ); 641 fRunPtr->SetMaxFloat_rDeltaPos_deltaMax( << 651 fRunPtr->SetMaxFloat_rDeltaPos_deltaMax( fMaxFloat_rDeltaPos_deltaMax ); 642 } 652 } >> 653 643 } 654 } 644 } 655 } 645 656 646 //....oooOO0OOooo........oooOO0OOooo........oo 657 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 647 658