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 // 26 // 27 /// \brief This class is a slightly modified v 27 /// \brief This class is a slightly modified version of G4Transportation 28 /// initially written by John Apostolak 28 /// initially written by John Apostolakis and colleagues 29 /// But it should use the exact same al 29 /// But it should use the exact same algorithm 30 // 30 // 31 // Contact : Mathieu Karamitros (kara (AT) cen 31 // Contact : Mathieu Karamitros (kara (AT) cenbg . in2p3 . fr) 32 // 32 // 33 // History : 33 // History : 34 // ----------- 34 // ----------- 35 // =========================================== 35 // ======================================================================= 36 // Modified: 36 // Modified: 37 // 28 Oct 2011, P.Gumpl./J.Ap: Detect gravi 37 // 28 Oct 2011, P.Gumpl./J.Ap: Detect gravity field, use magnetic moment 38 // 20 Nov 2008, J.Apostolakis: Push safety 38 // 20 Nov 2008, J.Apostolakis: Push safety to helper - after ComputeSafety 39 // 9 Nov 2007, J.Apostolakis: Flag for sho 39 // 9 Nov 2007, J.Apostolakis: Flag for short steps, push safety to helper 40 // 19 Jan 2006, P.MoraDeFreitas: Fix for su 40 // 19 Jan 2006, P.MoraDeFreitas: Fix for suspended tracks (StartTracking) 41 // 11 Aug 2004, M.Asai: Add G4VSensitiveDet 41 // 11 Aug 2004, M.Asai: Add G4VSensitiveDetector* for updating stepPoint. 42 // 21 June 2003, J.Apostolakis: Calling fiel 42 // 21 June 2003, J.Apostolakis: Calling field manager with 43 // track, to enable it to 43 // track, to enable it to configure its accuracy 44 // 13 May 2003, J.Apostolakis: Zero field a 44 // 13 May 2003, J.Apostolakis: Zero field areas now taken into 45 // account correclty in al 45 // account correclty in all cases (thanks to W Pokorski). 46 // 29 June 2001, J.Apostolakis, D.Cote-Ahern 46 // 29 June 2001, J.Apostolakis, D.Cote-Ahern, P.Gumplinger: 47 // correction for spin tra 47 // correction for spin tracking 48 // 20 Febr 2001, J.Apostolakis: update for 48 // 20 Febr 2001, J.Apostolakis: update for new FieldTrack 49 // 22 Sept 2000, V.Grichine: update of K 49 // 22 Sept 2000, V.Grichine: update of Kinetic Energy 50 // ------------------------------------------ 50 // --------------------------------------------------- 51 // 10 Oct 2011, M.Karamitros: G4ITTranspo 51 // 10 Oct 2011, M.Karamitros: G4ITTransportation created 52 // Created: 19 March 1997, J. Apostolakis 52 // Created: 19 March 1997, J. Apostolakis 53 // =========================================== 53 // ======================================================================= 54 // 54 // 55 // ------------------------------------------- 55 // ------------------------------------------------------------------- 56 56 57 #include "G4ITTransportation.hh" 57 #include "G4ITTransportation.hh" 58 << 59 #include "G4ChordFinder.hh" << 60 #include "G4FieldManager.hh" << 61 #include "G4FieldManagerStore.hh" << 62 #include "G4IT.hh" 58 #include "G4IT.hh" 63 #include "G4ITNavigator.hh" << 59 #include "G4TrackingInformation.hh" 64 #include "G4ITSafetyHelper.hh" << 60 #include "G4SystemOfUnits.hh" >> 61 #include "G4TransportationManager.hh" 65 #include "G4ITTransportationManager.hh" 62 #include "G4ITTransportationManager.hh" 66 #include "G4LowEnergyEmProcessSubType.hh" << 67 #include "G4ParticleTable.hh" << 68 #include "G4ProductionCutsTable.hh" 63 #include "G4ProductionCutsTable.hh" >> 64 #include "G4ParticleTable.hh" >> 65 #include "G4ITNavigator.hh" 69 #include "G4PropagatorInField.hh" 66 #include "G4PropagatorInField.hh" 70 #include "G4ReferenceCast.hh" << 67 #include "G4FieldManager.hh" 71 #include "G4SystemOfUnits.hh" << 68 #include "G4ChordFinder.hh" 72 #include "G4TrackingInformation.hh" << 69 #include "G4ITSafetyHelper.hh" 73 #include "G4TransportationManager.hh" << 70 #include "G4FieldManagerStore.hh" 74 #include "G4UnitsTable.hh" << 75 71 76 #include <memory> << 72 #include "G4UnitsTable.hh" >> 73 #include "G4ReferenceCast.hh" 77 74 78 class G4VSensitiveDetector; 75 class G4VSensitiveDetector; 79 76 80 #ifndef PrepareState 77 #ifndef PrepareState 81 # define PrepareState() << 78 #define PrepareState() G4ITTransportationState* __state = this->GetState<G4ITTransportationState>(); 82 G4ITTransportationState* __state = this->G << 83 #endif 79 #endif 84 80 85 #ifndef State 81 #ifndef State 86 #define State(theXInfo) (__state->theXInfo) 82 #define State(theXInfo) (__state->theXInfo) 87 #endif 83 #endif 88 84 89 //#define DEBUG_MEM 85 //#define DEBUG_MEM 90 86 91 #ifdef DEBUG_MEM 87 #ifdef DEBUG_MEM 92 #include "G4MemStat.hh" 88 #include "G4MemStat.hh" 93 using namespace G4MemStat; 89 using namespace G4MemStat; 94 using G4MemStat::MemStat; 90 using G4MemStat::MemStat; 95 #endif 91 #endif 96 92 97 //#define G4DEBUG_TRANSPORT 1 93 //#define G4DEBUG_TRANSPORT 1 98 94 99 G4ITTransportation::G4ITTransportation(const G 95 G4ITTransportation::G4ITTransportation(const G4String& aName, int verbose) : 100 G4VITProcess(aName, fTransportation), 96 G4VITProcess(aName, fTransportation), 101 fThreshold_Warning_Energy(100 * MeV), 97 fThreshold_Warning_Energy(100 * MeV), 102 fThreshold_Important_Energy(250 * MeV), 98 fThreshold_Important_Energy(250 * MeV), 103 << 99 fThresholdTrials(10), 104 fUnimportant_Energy(1 * MeV), // Old defa << 100 fUnimportant_Energy(1 * MeV), // Not used >> 101 fSumEnergyKilled(0.0), >> 102 fMaxEnergyKilled(0.0), >> 103 fShortStepOptimisation(false), // Old default: true (=fast short steps) 105 fVerboseLevel(verbose) 104 fVerboseLevel(verbose) 106 { 105 { 107 pParticleChange = &fParticleChange; 106 pParticleChange = &fParticleChange; 108 G4TransportationManager* transportMgr; 107 G4TransportationManager* transportMgr; 109 transportMgr = G4TransportationManager::GetT 108 transportMgr = G4TransportationManager::GetTransportationManager(); 110 G4ITTransportationManager* ITtransportMgr; 109 G4ITTransportationManager* ITtransportMgr; 111 ITtransportMgr = G4ITTransportationManager:: 110 ITtransportMgr = G4ITTransportationManager::GetTransportationManager(); 112 fLinearNavigator = ITtransportMgr->GetNaviga 111 fLinearNavigator = ITtransportMgr->GetNavigatorForTracking(); 113 fFieldPropagator = transportMgr->GetPropagat 112 fFieldPropagator = transportMgr->GetPropagatorInField(); 114 fpSafetyHelper = ITtransportMgr->GetSafetyHe 113 fpSafetyHelper = ITtransportMgr->GetSafetyHelper(); // New 115 114 116 // Cannot determine whether a field exists h 115 // Cannot determine whether a field exists here, as it would 117 // depend on the relative order of creating 116 // depend on the relative order of creating the detector's 118 // field and this process. That order is no 117 // field and this process. That order is not guaranted. 119 // Instead later the method DoesGlobalFieldE 118 // Instead later the method DoesGlobalFieldExist() is called 120 119 121 enableAtRestDoIt = false; 120 enableAtRestDoIt = false; 122 enableAlongStepDoIt = true; 121 enableAlongStepDoIt = true; 123 enablePostStepDoIt = true; 122 enablePostStepDoIt = true; 124 SetProcessSubType(fLowEnergyTransportation); << 123 SetProcessSubType(60); 125 SetInstantiateProcessState(true); 124 SetInstantiateProcessState(true); 126 G4VITProcess::SetInstantiateProcessState(fal 125 G4VITProcess::SetInstantiateProcessState(false); 127 fInstantiateProcessState = true; 126 fInstantiateProcessState = true; 128 127 129 G4VITProcess::fpState = std::make_shared<G4I << 128 G4VITProcess::fpState.reset(new G4ITTransportationState()); 130 /* 129 /* 131 if(fTransportationState == 0) 130 if(fTransportationState == 0) 132 { 131 { 133 G4cout << "KILL in G4ITTransportation::G4IT 132 G4cout << "KILL in G4ITTransportation::G4ITTransportation" << G4endl; 134 abort(); 133 abort(); 135 } 134 } 136 */ 135 */ 137 } 136 } 138 137 139 G4ITTransportation::G4ITTransportation(const G 138 G4ITTransportation::G4ITTransportation(const G4ITTransportation& right) : 140 G4VITProcess(right) 139 G4VITProcess(right) 141 { 140 { 142 // Copy attributes 141 // Copy attributes 143 fVerboseLevel = right.fVerboseLevel; 142 fVerboseLevel = right.fVerboseLevel; 144 fThreshold_Warning_Energy = right.fThreshold 143 fThreshold_Warning_Energy = right.fThreshold_Warning_Energy; 145 fThreshold_Important_Energy = right.fThresho 144 fThreshold_Important_Energy = right.fThreshold_Important_Energy; 146 fThresholdTrials = right.fThresholdTrials; 145 fThresholdTrials = right.fThresholdTrials; 147 fUnimportant_Energy = right.fUnimportant_Ene 146 fUnimportant_Energy = right.fUnimportant_Energy; 148 fSumEnergyKilled = right.fSumEnergyKilled; 147 fSumEnergyKilled = right.fSumEnergyKilled; 149 fMaxEnergyKilled = right.fMaxEnergyKilled; 148 fMaxEnergyKilled = right.fMaxEnergyKilled; 150 fShortStepOptimisation = right.fShortStepOpt 149 fShortStepOptimisation = right.fShortStepOptimisation; 151 150 152 // Setup Navigators 151 // Setup Navigators 153 G4TransportationManager* transportMgr; 152 G4TransportationManager* transportMgr; 154 transportMgr = G4TransportationManager::GetT 153 transportMgr = G4TransportationManager::GetTransportationManager(); 155 G4ITTransportationManager* ITtransportMgr; 154 G4ITTransportationManager* ITtransportMgr; 156 ITtransportMgr = G4ITTransportationManager:: 155 ITtransportMgr = G4ITTransportationManager::GetTransportationManager(); 157 fLinearNavigator = ITtransportMgr->GetNaviga 156 fLinearNavigator = ITtransportMgr->GetNavigatorForTracking(); 158 fFieldPropagator = transportMgr->GetPropagat 157 fFieldPropagator = transportMgr->GetPropagatorInField(); 159 fpSafetyHelper = ITtransportMgr->GetSafetyHe 158 fpSafetyHelper = ITtransportMgr->GetSafetyHelper(); // New 160 159 161 // Cannot determine whether a field exists h 160 // Cannot determine whether a field exists here, as it would 162 // depend on the relative order of creating 161 // depend on the relative order of creating the detector's 163 // field and this process. That order is no 162 // field and this process. That order is not guaranted. 164 // Instead later the method DoesGlobalFieldE 163 // Instead later the method DoesGlobalFieldExist() is called 165 164 166 enableAtRestDoIt = false; 165 enableAtRestDoIt = false; 167 enableAlongStepDoIt = true; 166 enableAlongStepDoIt = true; 168 enablePostStepDoIt = true; 167 enablePostStepDoIt = true; 169 168 170 pParticleChange = &fParticleChange; 169 pParticleChange = &fParticleChange; 171 SetInstantiateProcessState(true); 170 SetInstantiateProcessState(true); 172 G4VITProcess::SetInstantiateProcessState(fal 171 G4VITProcess::SetInstantiateProcessState(false); 173 fInstantiateProcessState = right.fInstantiat 172 fInstantiateProcessState = right.fInstantiateProcessState; 174 } 173 } 175 174 176 G4ITTransportation& G4ITTransportation::operat 175 G4ITTransportation& G4ITTransportation::operator=(const G4ITTransportation& /*right*/) 177 { 176 { 178 // if (this == &right) return *this; 177 // if (this == &right) return *this; 179 return *this; 178 return *this; 180 } 179 } 181 180 182 ////////////////////////////////////////////// 181 ////////////////////////////////////////////////////////////////////////////// 183 /// Process State 182 /// Process State 184 ////////////////////////////////////////////// 183 ////////////////////////////////////////////////////////////////////////////// 185 G4ITTransportation::G4ITTransportationState::G 184 G4ITTransportation::G4ITTransportationState::G4ITTransportationState() : 186 fCurrentTouchableHandle(nullptr) << 185 G4ProcessState(), fCurrentTouchableHandle(0) 187 { 186 { 188 fTransportEndPosition = G4ThreeVector(0, 0, 187 fTransportEndPosition = G4ThreeVector(0, 0, 0); 189 fTransportEndMomentumDir = G4ThreeVector(0, 188 fTransportEndMomentumDir = G4ThreeVector(0, 0, 0); 190 fTransportEndKineticEnergy = -1; 189 fTransportEndKineticEnergy = -1; 191 fTransportEndSpin = G4ThreeVector(0, 0, 0); 190 fTransportEndSpin = G4ThreeVector(0, 0, 0); 192 fMomentumChanged = false; 191 fMomentumChanged = false; 193 fEnergyChanged = false; 192 fEnergyChanged = false; 194 fEndGlobalTimeComputed = false; 193 fEndGlobalTimeComputed = false; 195 fCandidateEndGlobalTime = -1; 194 fCandidateEndGlobalTime = -1; 196 fParticleIsLooping = false; 195 fParticleIsLooping = false; 197 196 198 static G4ThreadLocal G4TouchableHandle *null << 197 static G4ThreadLocal G4TouchableHandle *nullTouchableHandle = 0; 199 if (nullTouchableHandle == nullptr) nullTouc << 198 if (!nullTouchableHandle) nullTouchableHandle = new G4TouchableHandle; 200 // Points to (G4VTouchable*) 0 199 // Points to (G4VTouchable*) 0 201 200 202 fCurrentTouchableHandle = *nullTouchableHand 201 fCurrentTouchableHandle = *nullTouchableHandle; 203 fGeometryLimitedStep = false; 202 fGeometryLimitedStep = false; 204 fPreviousSftOrigin = G4ThreeVector(0, 0, 0); 203 fPreviousSftOrigin = G4ThreeVector(0, 0, 0); 205 fPreviousSafety = 0.0; 204 fPreviousSafety = 0.0; 206 fNoLooperTrials = 0; << 205 fNoLooperTrials = false; 207 fEndPointDistance = -1; 206 fEndPointDistance = -1; 208 } 207 } 209 208 210 G4ITTransportation::G4ITTransportationState::~ 209 G4ITTransportation::G4ITTransportationState::~G4ITTransportationState() 211 { 210 { 212 ; 211 ; 213 } 212 } 214 213 215 G4ITTransportation::~G4ITTransportation() 214 G4ITTransportation::~G4ITTransportation() 216 { 215 { 217 #ifdef G4VERBOSE 216 #ifdef G4VERBOSE 218 if ((fVerboseLevel > 0) && (fSumEnergyKilled 217 if ((fVerboseLevel > 0) && (fSumEnergyKilled > 0.0)) 219 { 218 { 220 G4cout << " G4ITTransportation: Statistics 219 G4cout << " G4ITTransportation: Statistics for looping particles " 221 << G4endl; 220 << G4endl; 222 G4cout << " Sum of energy of loopers kil 221 G4cout << " Sum of energy of loopers killed: " 223 << fSumEnergyKilled << G4endl; 222 << fSumEnergyKilled << G4endl; 224 G4cout << " Max energy of loopers killed 223 G4cout << " Max energy of loopers killed: " 225 << fMaxEnergyKilled << G4endl; 224 << fMaxEnergyKilled << G4endl; 226 } 225 } 227 #endif 226 #endif 228 } 227 } 229 228 230 void G4ITTransportation::BuildPhysicsTable(con 229 void G4ITTransportation::BuildPhysicsTable(const G4ParticleDefinition&) 231 { 230 { 232 fpSafetyHelper->InitialiseHelper(); 231 fpSafetyHelper->InitialiseHelper(); 233 } 232 } 234 233 235 G4bool G4ITTransportation::DoesGlobalFieldExis 234 G4bool G4ITTransportation::DoesGlobalFieldExist() 236 { 235 { 237 G4TransportationManager* transportMgr; 236 G4TransportationManager* transportMgr; 238 transportMgr = G4TransportationManager::GetT 237 transportMgr = G4TransportationManager::GetTransportationManager(); 239 238 240 // fFieldExists= transportMgr->GetFieldManag 239 // fFieldExists= transportMgr->GetFieldManager()->DoesFieldExist(); 241 // return fFieldExists; 240 // return fFieldExists; 242 return transportMgr->GetFieldManager()->Does 241 return transportMgr->GetFieldManager()->DoesFieldExist(); 243 } 242 } 244 243 245 ////////////////////////////////////////////// 244 ////////////////////////////////////////////////////////////////////////// 246 // 245 // 247 // Responsibilities: 246 // Responsibilities: 248 // Find whether the geometry limits the Ste 247 // Find whether the geometry limits the Step, and to what length 249 // Calculate the new value of the safety an 248 // Calculate the new value of the safety and return it. 250 // Store the final time, position and momen 249 // Store the final time, position and momentum. 251 G4double 250 G4double 252 G4ITTransportation:: 251 G4ITTransportation:: 253 AlongStepGetPhysicalInteractionLength(const G4 252 AlongStepGetPhysicalInteractionLength(const G4Track& track, 254 G4double 253 G4double, 255 G4double 254 G4double currentMinimumStep, 256 G4double 255 G4double& currentSafety, 257 G4GPILSe 256 G4GPILSelection* selection) 258 { 257 { 259 PrepareState(); << 258 PrepareState() 260 G4double geometryStepLength(-1.0), newSafety 259 G4double geometryStepLength(-1.0), newSafety(-1.0); 261 260 262 State(fParticleIsLooping) = false; 261 State(fParticleIsLooping) = false; 263 State(fEndGlobalTimeComputed) = false; 262 State(fEndGlobalTimeComputed) = false; 264 State(fGeometryLimitedStep) = false; 263 State(fGeometryLimitedStep) = false; 265 264 266 // Initial actions moved to StartTrack() 265 // Initial actions moved to StartTrack() 267 // -------------------------------------- 266 // -------------------------------------- 268 // Note: in case another process changes tou 267 // Note: in case another process changes touchable handle 269 // it will be necessary to add here (for 268 // it will be necessary to add here (for all steps) 270 // State(fCurrentTouchableHandle) = track.Ge 269 // State(fCurrentTouchableHandle) = track.GetTouchableHandle(); 271 270 272 // GPILSelection is set to defaule value of 271 // GPILSelection is set to defaule value of CandidateForSelection 273 // It is a return value 272 // It is a return value 274 // 273 // 275 *selection = CandidateForSelection; 274 *selection = CandidateForSelection; 276 275 277 // Get initial Energy/Momentum of the track 276 // Get initial Energy/Momentum of the track 278 // 277 // 279 const G4DynamicParticle* pParticle = track.G 278 const G4DynamicParticle* pParticle = track.GetDynamicParticle(); 280 // const G4ParticleDefinition* pParticleDef 279 // const G4ParticleDefinition* pParticleDef = pParticle->GetDefinition() ; 281 G4ThreeVector startMomentumDir = pParticle-> 280 G4ThreeVector startMomentumDir = pParticle->GetMomentumDirection(); 282 G4ThreeVector startPosition = track.GetPosit 281 G4ThreeVector startPosition = track.GetPosition(); 283 282 284 // G4double theTime = track.GetGlob 283 // G4double theTime = track.GetGlobalTime() ; 285 284 286 // The Step Point safety can be limited by o 285 // The Step Point safety can be limited by other geometries and/or the 287 // assumptions of any process - it's not alw 286 // assumptions of any process - it's not always the geometrical safety. 288 // We calculate the starting point's isotrop 287 // We calculate the starting point's isotropic safety here. 289 // 288 // 290 G4ThreeVector OriginShift = startPosition - 289 G4ThreeVector OriginShift = startPosition - State(fPreviousSftOrigin); 291 G4double MagSqShift = OriginShift.mag2(); 290 G4double MagSqShift = OriginShift.mag2(); 292 if (MagSqShift >= sqr(State(fPreviousSafety) 291 if (MagSqShift >= sqr(State(fPreviousSafety))) 293 { 292 { 294 currentSafety = 0.0; 293 currentSafety = 0.0; 295 } 294 } 296 else 295 else 297 { 296 { 298 currentSafety = State(fPreviousSafety) - s 297 currentSafety = State(fPreviousSafety) - std::sqrt(MagSqShift); 299 } 298 } 300 299 301 // Is the particle charged ? 300 // Is the particle charged ? 302 // 301 // 303 G4double particleCharge = pParticle->GetChar 302 G4double particleCharge = pParticle->GetCharge(); 304 303 305 // There is no need to locate the current vo 304 // There is no need to locate the current volume. It is Done elsewhere: 306 // On track construction 305 // On track construction 307 // By the tracking, after all AlongStepDoI 306 // By the tracking, after all AlongStepDoIts, in "Relocation" 308 307 309 // Check whether the particle have an (EM) f 308 // Check whether the particle have an (EM) field force exerting upon it 310 // 309 // 311 G4FieldManager* fieldMgr = nullptr; << 310 G4FieldManager* fieldMgr = 0; 312 G4bool fieldExertsForce = false; 311 G4bool fieldExertsForce = false; 313 if ((particleCharge != 0.0)) 312 if ((particleCharge != 0.0)) 314 { 313 { 315 fieldMgr = fFieldPropagator->FindAndSetFie 314 fieldMgr = fFieldPropagator->FindAndSetFieldManager(track.GetVolume()); 316 if (fieldMgr != nullptr) << 315 if (fieldMgr != 0) 317 { 316 { 318 // Message the field Manager, to configu 317 // Message the field Manager, to configure it for this track 319 fieldMgr->ConfigureForTrack(&track); 318 fieldMgr->ConfigureForTrack(&track); 320 // Moved here, in order to allow a trans 319 // Moved here, in order to allow a transition 321 // from a zero-field status (with fie 320 // from a zero-field status (with fieldMgr->(field)0 322 // to a finite field status 321 // to a finite field status 323 322 324 // If the field manager has no field, th 323 // If the field manager has no field, there is no field ! 325 fieldExertsForce = (fieldMgr->GetDetecto << 324 fieldExertsForce = (fieldMgr->GetDetectorField() != 0); 326 } 325 } 327 } 326 } 328 327 329 // G4cout << " G4Transport: field exerts fo 328 // G4cout << " G4Transport: field exerts force= " << fieldExertsForce 330 // << " fieldMgr= " << fieldMgr << G4endl 329 // << " fieldMgr= " << fieldMgr << G4endl; 331 330 332 // Choose the calculation of the transportat 331 // Choose the calculation of the transportation: Field or not 333 // 332 // 334 if (!fieldExertsForce) 333 if (!fieldExertsForce) 335 { 334 { 336 G4double linearStepLength; 335 G4double linearStepLength; 337 if (fShortStepOptimisation && (currentMini 336 if (fShortStepOptimisation && (currentMinimumStep <= currentSafety)) 338 { 337 { 339 // The Step is guaranteed to be taken 338 // The Step is guaranteed to be taken 340 // 339 // 341 geometryStepLength = currentMinimumStep; 340 geometryStepLength = currentMinimumStep; 342 State(fGeometryLimitedStep) = false; 341 State(fGeometryLimitedStep) = false; 343 } 342 } 344 else 343 else 345 { 344 { 346 // Find whether the straight path inter 345 // Find whether the straight path intersects a volume 347 // 346 // 348 // fLinearNavigator->SetNavigatorState(G 347 // fLinearNavigator->SetNavigatorState(GetIT(track)->GetTrackingInfo()->GetNavigatorState()); 349 linearStepLength = fLinearNavigator->Com 348 linearStepLength = fLinearNavigator->ComputeStep(startPosition, 350 349 startMomentumDir, 351 350 currentMinimumStep, 352 351 newSafety); 353 352 354 // G4cout << "linearStepLength : " << G4 353 // G4cout << "linearStepLength : " << G4BestUnit(linearStepLength,"Length") 355 // << " | currentMinimumStep: " << cu 354 // << " | currentMinimumStep: " << currentMinimumStep 356 // << " | trackID: " << track.GetTrac 355 // << " | trackID: " << track.GetTrackID() << G4endl; 357 356 358 // Remember last safety origin & value. 357 // Remember last safety origin & value. 359 // 358 // 360 State(fPreviousSftOrigin) = startPositio 359 State(fPreviousSftOrigin) = startPosition; 361 State(fPreviousSafety) = newSafety; 360 State(fPreviousSafety) = newSafety; 362 361 363 G4TrackStateManager& trackStateMan = Get 362 G4TrackStateManager& trackStateMan = GetIT(track)->GetTrackingInfo() 364 ->GetTrackStateManager(); 363 ->GetTrackStateManager(); 365 fpSafetyHelper->LoadTrackState(trackStat 364 fpSafetyHelper->LoadTrackState(trackStateMan); 366 // fpSafetyHelper->SetTrackState(state); 365 // fpSafetyHelper->SetTrackState(state); 367 fpSafetyHelper->SetCurrentSafety(newSafe 366 fpSafetyHelper->SetCurrentSafety(newSafety, 368 State(f 367 State(fTransportEndPosition)); 369 fpSafetyHelper->ResetTrackState(); 368 fpSafetyHelper->ResetTrackState(); 370 369 371 // The safety at the initial point has b 370 // The safety at the initial point has been re-calculated: 372 // 371 // 373 currentSafety = newSafety; 372 currentSafety = newSafety; 374 373 375 State(fGeometryLimitedStep) = (linearSte 374 State(fGeometryLimitedStep) = (linearStepLength <= currentMinimumStep); 376 if (State(fGeometryLimitedStep)) 375 if (State(fGeometryLimitedStep)) 377 { 376 { 378 // The geometry limits the Step size ( 377 // The geometry limits the Step size (an intersection was found.) 379 geometryStepLength = linearStepLength; 378 geometryStepLength = linearStepLength; 380 } 379 } 381 else 380 else 382 { 381 { 383 // The full Step is taken. 382 // The full Step is taken. 384 geometryStepLength = currentMinimumSte 383 geometryStepLength = currentMinimumStep; 385 } 384 } 386 } 385 } 387 State(fEndPointDistance) = geometryStepLen 386 State(fEndPointDistance) = geometryStepLength; 388 387 389 // Calculate final position 388 // Calculate final position 390 // 389 // 391 State(fTransportEndPosition) = startPositi 390 State(fTransportEndPosition) = startPosition 392 + geometryStepLength * startMomentumDi 391 + geometryStepLength * startMomentumDir; 393 392 394 // Momentum direction, energy and polarisa 393 // Momentum direction, energy and polarisation are unchanged by transport 395 // 394 // 396 State(fTransportEndMomentumDir) = startMom 395 State(fTransportEndMomentumDir) = startMomentumDir; 397 State(fTransportEndKineticEnergy) = track. 396 State(fTransportEndKineticEnergy) = track.GetKineticEnergy(); 398 State(fTransportEndSpin) = track.GetPolari 397 State(fTransportEndSpin) = track.GetPolarization(); 399 State(fParticleIsLooping) = false; 398 State(fParticleIsLooping) = false; 400 State(fMomentumChanged) = false; 399 State(fMomentumChanged) = false; 401 State(fEndGlobalTimeComputed) = true; 400 State(fEndGlobalTimeComputed) = true; 402 State(theInteractionTimeLeft) = State(fEnd 401 State(theInteractionTimeLeft) = State(fEndPointDistance) 403 / track.GetVelocity(); 402 / track.GetVelocity(); 404 State(fCandidateEndGlobalTime) = State(the 403 State(fCandidateEndGlobalTime) = State(theInteractionTimeLeft) 405 + track.GetGlobalTime(); 404 + track.GetGlobalTime(); 406 /* 405 /* 407 G4cout << "track.GetVelocity() : " 406 G4cout << "track.GetVelocity() : " 408 << track.GetVelocity() << G4endl; 407 << track.GetVelocity() << G4endl; 409 G4cout << "State(endpointDistance) : " 408 G4cout << "State(endpointDistance) : " 410 << G4BestUnit(State(endpointDistanc 409 << G4BestUnit(State(endpointDistance),"Length") << G4endl; 411 G4cout << "State(theInteractionTimeLeft) : 410 G4cout << "State(theInteractionTimeLeft) : " 412 << G4BestUnit(State(theInteractionT 411 << G4BestUnit(State(theInteractionTimeLeft),"Time") << G4endl; 413 G4cout << "track.GetGlobalTime() : " 412 G4cout << "track.GetGlobalTime() : " 414 << G4BestUnit(track.GetGlobalTime() 413 << G4BestUnit(track.GetGlobalTime(),"Time") << G4endl; 415 */ 414 */ 416 } 415 } 417 else // A field exerts force 416 else // A field exerts force 418 { 417 { 419 418 420 G4ExceptionDescription exceptionDescriptio 419 G4ExceptionDescription exceptionDescription; 421 exceptionDescription 420 exceptionDescription 422 << "ITTransportation does not support 421 << "ITTransportation does not support external fields."; 423 exceptionDescription 422 exceptionDescription 424 << " If you are dealing with a tradiat 423 << " If you are dealing with a tradiational MC simulation, "; 425 exceptionDescription << "please use G4Tran 424 exceptionDescription << "please use G4Transportation."; 426 425 427 G4Exception("G4ITTransportation::AlongStep 426 G4Exception("G4ITTransportation::AlongStepGetPhysicalInteractionLength", 428 "NoExternalFieldSupport", Fata 427 "NoExternalFieldSupport", FatalException, exceptionDescription); 429 /* 428 /* 430 G4double momentumMagnitude = pParti 429 G4double momentumMagnitude = pParticle->GetTotalMomentum() ; 431 // G4ThreeVector EndUnitMomentum 430 // G4ThreeVector EndUnitMomentum ; 432 G4double lengthAlongCurve ; 431 G4double lengthAlongCurve ; 433 G4double restMass = pParticleDef->G 432 G4double restMass = pParticleDef->GetPDGMass() ; 434 433 435 fFieldPropagator->SetChargeMomentumMass( 434 fFieldPropagator->SetChargeMomentumMass( particleCharge, // in e+ units 436 momentumMagnitude, // in Mev/c 435 momentumMagnitude, // in Mev/c 437 restMass ) ; 436 restMass ) ; 438 437 439 G4ThreeVector spin = track.GetPola 438 G4ThreeVector spin = track.GetPolarization() ; 440 G4FieldTrack aFieldTrack = G4FieldTrack( 439 G4FieldTrack aFieldTrack = G4FieldTrack( startPosition, 441 track.GetMomentumDirection(), 440 track.GetMomentumDirection(), 442 0.0, 441 0.0, 443 track.GetKineticEnergy(), 442 track.GetKineticEnergy(), 444 restMass, 443 restMass, 445 track.GetVelocity(), 444 track.GetVelocity(), 446 track.GetGlobalTime(), // Lab. 445 track.GetGlobalTime(), // Lab. 447 track.GetProperTime(), // Part. 446 track.GetProperTime(), // Part. 448 &spin ) ; 447 &spin ) ; 449 if( currentMinimumStep > 0 ) 448 if( currentMinimumStep > 0 ) 450 { 449 { 451 // Do the Transport in the field (non rec 450 // Do the Transport in the field (non recti-linear) 452 // 451 // 453 lengthAlongCurve = fFieldPropagator->Comp 452 lengthAlongCurve = fFieldPropagator->ComputeStep( aFieldTrack, 454 currentMinimumStep, 453 currentMinimumStep, 455 currentSafety, 454 currentSafety, 456 track.GetVolume() ) ; 455 track.GetVolume() ) ; 457 State(fGeometryLimitedStep)= lengthAlongC 456 State(fGeometryLimitedStep)= lengthAlongCurve < currentMinimumStep; 458 if( State(fGeometryLimitedStep) ) 457 if( State(fGeometryLimitedStep) ) 459 { 458 { 460 geometryStepLength = lengthAlongCurve ; 459 geometryStepLength = lengthAlongCurve ; 461 } 460 } 462 else 461 else 463 { 462 { 464 geometryStepLength = currentMinimumStep 463 geometryStepLength = currentMinimumStep ; 465 } 464 } 466 465 467 // Remember last safety origin & value. 466 // Remember last safety origin & value. 468 // 467 // 469 State(fPreviousSftOrigin) = startPosition 468 State(fPreviousSftOrigin) = startPosition ; 470 State(fPreviousSafety) = currentSafety 469 State(fPreviousSafety) = currentSafety ; 471 fpSafetyHelper->SetCurrentSafety( newSafe 470 fpSafetyHelper->SetCurrentSafety( newSafety, startPosition); 472 } 471 } 473 else 472 else 474 { 473 { 475 geometryStepLength = lengthAlongCurve= 474 geometryStepLength = lengthAlongCurve= 0.0 ; 476 State(fGeometryLimitedStep) = false ; 475 State(fGeometryLimitedStep) = false ; 477 } 476 } 478 477 479 // Get the End-Position and End-Momentum 478 // Get the End-Position and End-Momentum (Dir-ection) 480 // 479 // 481 State(fTransportEndPosition) = aFieldTrac 480 State(fTransportEndPosition) = aFieldTrack.GetPosition() ; 482 481 483 // Momentum: Magnitude and direction can 482 // Momentum: Magnitude and direction can be changed too now ... 484 // 483 // 485 State(fMomentumChanged) = true ; 484 State(fMomentumChanged) = true ; 486 State(fTransportEndMomentumDir) = aFieldT 485 State(fTransportEndMomentumDir) = aFieldTrack.GetMomentumDir() ; 487 486 488 State(fTransportEndKineticEnergy) = aFie 487 State(fTransportEndKineticEnergy) = aFieldTrack.GetKineticEnergy() ; 489 488 490 if( fFieldPropagator->GetCurrentFieldMana 489 if( fFieldPropagator->GetCurrentFieldManager()->DoesFieldChangeEnergy() ) 491 { 490 { 492 // If the field can change energy, then t 491 // If the field can change energy, then the time must be integrated 493 // - so this should have been updated 492 // - so this should have been updated 494 // 493 // 495 State(fCandidateEndGlobalTime) = aField 494 State(fCandidateEndGlobalTime) = aFieldTrack.GetLabTimeOfFlight(); 496 State(fEndGlobalTimeComputed) = true; 495 State(fEndGlobalTimeComputed) = true; 497 496 498 State(theInteractionTimeLeft) = State(fCa 497 State(theInteractionTimeLeft) = State(fCandidateEndGlobalTime) - 499 track 498 track.GetGlobalTime() ; 500 499 501 // was ( State(fCandidateEndGlobalTime) ! 500 // was ( State(fCandidateEndGlobalTime) != track.GetGlobalTime() ); 502 // a cleaner way is to have FieldTrack kn 501 // a cleaner way is to have FieldTrack knowing whether time is updated. 503 } 502 } 504 else 503 else 505 { 504 { 506 // The energy should be unchanged by fiel 505 // The energy should be unchanged by field transport, 507 // - so the time changed will be calcu 506 // - so the time changed will be calculated elsewhere 508 // 507 // 509 State(fEndGlobalTimeComputed) = false; 508 State(fEndGlobalTimeComputed) = false; 510 509 511 // Check that the integration preserved t 510 // Check that the integration preserved the energy 512 // - and if not correct this! 511 // - and if not correct this! 513 G4double startEnergy= track.GetKineticEn 512 G4double startEnergy= track.GetKineticEnergy(); 514 G4double endEnergy= State(fTransportEndK 513 G4double endEnergy= State(fTransportEndKineticEnergy); 515 514 516 static G4int no_inexact_steps=0, no_large 515 static G4int no_inexact_steps=0, no_large_ediff; 517 G4double absEdiff = std::fabs(startEnergy 516 G4double absEdiff = std::fabs(startEnergy- endEnergy); 518 if( absEdiff > perMillion * endEnergy ) 517 if( absEdiff > perMillion * endEnergy ) 519 { 518 { 520 no_inexact_steps++; 519 no_inexact_steps++; 521 // Possible statistics keeping here ... 520 // Possible statistics keeping here ... 522 } 521 } 523 #ifdef G4VERBOSE 522 #ifdef G4VERBOSE 524 if( fVerboseLevel > 1 ) 523 if( fVerboseLevel > 1 ) 525 { 524 { 526 if( std::fabs(startEnergy- endEnergy) > p 525 if( std::fabs(startEnergy- endEnergy) > perThousand * endEnergy ) 527 { 526 { 528 static G4int no_warnings= 0, warnModulo=1 527 static G4int no_warnings= 0, warnModulo=1, moduloFactor= 10; 529 no_large_ediff ++; 528 no_large_ediff ++; 530 if( (no_large_ediff% warnModulo) == 0 ) 529 if( (no_large_ediff% warnModulo) == 0 ) 531 { 530 { 532 no_warnings++; 531 no_warnings++; 533 G4cout << "WARNING - G4Transportation::Al 532 G4cout << "WARNING - G4Transportation::AlongStepGetPIL() " 534 << " Energy change in Step is above 1^- 533 << " Energy change in Step is above 1^-3 relative value. " << G4endl 535 << " Relative change in 'tracking' step 534 << " Relative change in 'tracking' step = " 536 << std::setw(15) << (endEnergy-startEnerg 535 << std::setw(15) << (endEnergy-startEnergy)/startEnergy << G4endl 537 << " Starting E= " << std::setw(12) < 536 << " Starting E= " << std::setw(12) << startEnergy / MeV << " MeV " 538 << G4endl 537 << G4endl 539 << " Ending E= " << std::setw(12) < 538 << " Ending E= " << std::setw(12) << endEnergy / MeV << " MeV " 540 << G4endl; 539 << G4endl; 541 G4cout << " Energy has been corrected -- 540 G4cout << " Energy has been corrected -- however, review" 542 << " field propagation parameters for acc 541 << " field propagation parameters for accuracy." << G4endl; 543 if( (fVerboseLevel > 2 ) || (no_warnings< 542 if( (fVerboseLevel > 2 ) || (no_warnings<4) || 544 (no_large_ediff == warnModulo * modul 543 (no_large_ediff == warnModulo * moduloFactor) ) 545 { 544 { 546 G4cout << " These include EpsilonStepMax( 545 G4cout << " These include EpsilonStepMax(/Min) in G4FieldManager " 547 << " which determine fractional error per 546 << " which determine fractional error per step for integrated quantities. " 548 << G4endl 547 << G4endl 549 << " Note also the influence of the permi 548 << " Note also the influence of the permitted number of integration steps." 550 << G4endl; 549 << G4endl; 551 } 550 } 552 G4cerr << "ERROR - G4Transportation::Alon 551 G4cerr << "ERROR - G4Transportation::AlongStepGetPIL()" << G4endl 553 << " Bad 'endpoint'. Energy change 552 << " Bad 'endpoint'. Energy change detected" 554 << " and corrected. " 553 << " and corrected. " 555 << " Has occurred already " 554 << " Has occurred already " 556 << no_large_ediff << " times." << G4endl; 555 << no_large_ediff << " times." << G4endl; 557 if( no_large_ediff == warnModulo * modulo 556 if( no_large_ediff == warnModulo * moduloFactor ) 558 { 557 { 559 warnModulo *= moduloFactor; 558 warnModulo *= moduloFactor; 560 } 559 } 561 } 560 } 562 } 561 } 563 } // end of if (fVerboseLevel) 562 } // end of if (fVerboseLevel) 564 #endif 563 #endif 565 // Correct the energy for fields that con 564 // Correct the energy for fields that conserve it 566 // This - hides the integration error 565 // This - hides the integration error 567 // - but gives a better physical an 566 // - but gives a better physical answer 568 State(fTransportEndKineticEnergy)= track. 567 State(fTransportEndKineticEnergy)= track.GetKineticEnergy(); 569 } 568 } 570 569 571 State(fTransportEndSpin) = aFieldTrack.Ge 570 State(fTransportEndSpin) = aFieldTrack.GetSpin(); 572 State(fParticleIsLooping) = fFieldPropaga 571 State(fParticleIsLooping) = fFieldPropagator->IsParticleLooping() ; 573 State(endpointDistance) = (State(fTrans 572 State(endpointDistance) = (State(fTransportEndPosition) - 574 startPositio 573 startPosition).mag() ; 575 // State(theInteractionTimeLeft) = 574 // State(theInteractionTimeLeft) = 576 track.GetVelocity()/Sta 575 track.GetVelocity()/State(endpointDistance) ; 577 */ 576 */ 578 } 577 } 579 578 580 // If we are asked to go a step length of 0, 579 // If we are asked to go a step length of 0, and we are on a boundary 581 // then a boundary will also limit the step 580 // then a boundary will also limit the step -> we must flag this. 582 // 581 // 583 if (currentMinimumStep == 0.0) 582 if (currentMinimumStep == 0.0) 584 { 583 { 585 if (currentSafety == 0.0) 584 if (currentSafety == 0.0) 586 { 585 { 587 State(fGeometryLimitedStep) = true; 586 State(fGeometryLimitedStep) = true; 588 // G4cout << "!!!! Safety is NULL, 587 // G4cout << "!!!! Safety is NULL, on the Boundary !!!!!" << G4endl; 589 // G4cout << " Track position : " < 588 // G4cout << " Track position : " << track.GetPosition() /nanometer 590 // << G4endl; 589 // << G4endl; 591 } 590 } 592 } 591 } 593 592 594 // Update the safety starting from the end-p 593 // Update the safety starting from the end-point, 595 // if it will become negative at the end-poi 594 // if it will become negative at the end-point. 596 // 595 // 597 if (currentSafety < State(fEndPointDistance) 596 if (currentSafety < State(fEndPointDistance)) 598 { 597 { 599 // if( particleCharge == 0.0 ) 598 // if( particleCharge == 0.0 ) 600 // G4cout << " Avoiding call to Compu 599 // G4cout << " Avoiding call to ComputeSafety : charge = 0.0 " << G4endl; 601 600 602 if (particleCharge != 0.0) 601 if (particleCharge != 0.0) 603 { 602 { 604 603 605 G4double endSafety = fLinearNavigator->C 604 G4double endSafety = fLinearNavigator->ComputeSafety( 606 State(fTransportEndPosition)); 605 State(fTransportEndPosition)); 607 currentSafety = endSafety; 606 currentSafety = endSafety; 608 State(fPreviousSftOrigin) = State(fTrans 607 State(fPreviousSftOrigin) = State(fTransportEndPosition); 609 State(fPreviousSafety) = currentSafety; 608 State(fPreviousSafety) = currentSafety; 610 609 611 /* 610 /* 612 G4VTrackStateHandle state = 611 G4VTrackStateHandle state = 613 GetIT(track)->GetTrackingInfo()->GetTra 612 GetIT(track)->GetTrackingInfo()->GetTrackState(fpSafetyHelper); 614 */ 613 */ 615 G4TrackStateManager& trackStateMan = Get 614 G4TrackStateManager& trackStateMan = GetIT(track)->GetTrackingInfo() 616 ->GetTrackStateManager(); 615 ->GetTrackStateManager(); 617 fpSafetyHelper->LoadTrackState(trackStat 616 fpSafetyHelper->LoadTrackState(trackStateMan); 618 // fpSafetyHelper->SetTrackState(state); 617 // fpSafetyHelper->SetTrackState(state); 619 fpSafetyHelper->SetCurrentSafety(current 618 fpSafetyHelper->SetCurrentSafety(currentSafety, 620 State(f 619 State(fTransportEndPosition)); 621 fpSafetyHelper->ResetTrackState(); 620 fpSafetyHelper->ResetTrackState(); 622 621 623 // Because the Stepping Manager assumes 622 // Because the Stepping Manager assumes it is from the start point, 624 // add the StepLength 623 // add the StepLength 625 // 624 // 626 currentSafety += State(fEndPointDistance 625 currentSafety += State(fEndPointDistance); 627 626 628 #ifdef G4DEBUG_TRANSPORT 627 #ifdef G4DEBUG_TRANSPORT 629 G4cout.precision(12); 628 G4cout.precision(12); 630 G4cout << "***G4Transportation::AlongSte 629 G4cout << "***G4Transportation::AlongStepGPIL ** " << G4endl; 631 G4cout << " Called Navigator->ComputeSa 630 G4cout << " Called Navigator->ComputeSafety at " 632 << State(fTransportEndPosition) 631 << State(fTransportEndPosition) 633 << " and it returned safety= " << 632 << " and it returned safety= " << endSafety << G4endl; 634 G4cout << " Adding endpoint distance " 633 G4cout << " Adding endpoint distance " << State(fEndPointDistance) 635 << " to obtain pseudo-safety= " 634 << " to obtain pseudo-safety= " << currentSafety << G4endl; 636 #endif 635 #endif 637 } 636 } 638 } 637 } 639 638 640 // fParticleChange.ProposeTrueStepLength(geo 639 // fParticleChange.ProposeTrueStepLength(geometryStepLength) ; 641 640 642 // G4cout << "G4ITTransportation::AlongStepGe 641 // G4cout << "G4ITTransportation::AlongStepGetPhysicalInteractionLength = " 643 // << G4BestUnit(geometryStepLength,"L 642 // << G4BestUnit(geometryStepLength,"Length") << G4endl; 644 643 645 return geometryStepLength; 644 return geometryStepLength; 646 } 645 } 647 646 648 void G4ITTransportation::ComputeStep(const G4T 647 void G4ITTransportation::ComputeStep(const G4Track& track, 649 const G4S 648 const G4Step& /*step*/, 650 const dou 649 const double timeStep, 651 double& o 650 double& oPhysicalStep) 652 { 651 { 653 PrepareState(); << 652 PrepareState() 654 const G4DynamicParticle* pParticle = track.G 653 const G4DynamicParticle* pParticle = track.GetDynamicParticle(); 655 G4ThreeVector startMomentumDir = pParticle-> 654 G4ThreeVector startMomentumDir = pParticle->GetMomentumDirection(); 656 G4ThreeVector startPosition = track.GetPosit 655 G4ThreeVector startPosition = track.GetPosition(); 657 656 658 track.CalculateVelocity(); 657 track.CalculateVelocity(); 659 G4double initialVelocity = track.GetVelocity 658 G4double initialVelocity = track.GetVelocity(); 660 659 661 State(fGeometryLimitedStep) = false; 660 State(fGeometryLimitedStep) = false; 662 661 663 ///////////////////////// 662 ///////////////////////// 664 // !!! CASE NO FIELD !!! 663 // !!! CASE NO FIELD !!! 665 ///////////////////////// 664 ///////////////////////// 666 State(fCandidateEndGlobalTime) = timeStep + 665 State(fCandidateEndGlobalTime) = timeStep + track.GetGlobalTime(); 667 State(fEndGlobalTimeComputed) = true; 666 State(fEndGlobalTimeComputed) = true; 668 667 669 // Choose the calculation of the transportat 668 // Choose the calculation of the transportation: Field or not 670 // 669 // 671 if (!State(fMomentumChanged)) 670 if (!State(fMomentumChanged)) 672 { 671 { 673 // G4cout << "Momentum has not chan 672 // G4cout << "Momentum has not changed" << G4endl; 674 fParticleChange.ProposeVelocity(initialVel 673 fParticleChange.ProposeVelocity(initialVelocity); 675 oPhysicalStep = initialVelocity * timeStep 674 oPhysicalStep = initialVelocity * timeStep; 676 675 677 // Calculate final position 676 // Calculate final position 678 // 677 // 679 State(fTransportEndPosition) = startPositi 678 State(fTransportEndPosition) = startPosition 680 + oPhysicalStep * startMomentumDir; 679 + oPhysicalStep * startMomentumDir; 681 } 680 } 682 } 681 } 683 682 684 ////////////////////////////////////////////// 683 ////////////////////////////////////////////////////////////////////////// 685 // 684 // 686 // Initialize ParticleChange (by setting al 685 // Initialize ParticleChange (by setting all its members equal 687 // to correspond 686 // to corresponding members in G4Track) 688 #include "G4ParticleTable.hh" 687 #include "G4ParticleTable.hh" 689 G4VParticleChange* G4ITTransportation::AlongSt 688 G4VParticleChange* G4ITTransportation::AlongStepDoIt(const G4Track& track, 690 689 const G4Step& stepData) 691 { 690 { 692 691 693 #if defined (DEBUG_MEM) 692 #if defined (DEBUG_MEM) 694 MemStat mem_first, mem_second, mem_diff; 693 MemStat mem_first, mem_second, mem_diff; 695 #endif 694 #endif 696 695 697 #if defined (DEBUG_MEM) 696 #if defined (DEBUG_MEM) 698 mem_first = MemoryUsage(); 697 mem_first = MemoryUsage(); 699 #endif 698 #endif 700 699 701 PrepareState(); << 700 PrepareState() 702 << 701 703 // G4cout << "G4ITTransportation::AlongStepD 702 // G4cout << "G4ITTransportation::AlongStepDoIt" << G4endl; 704 // set pdefOpticalPhoton 703 // set pdefOpticalPhoton 705 // Andrea Dotti: the following statement sho 704 // Andrea Dotti: the following statement should be in a single line: 706 // G4-MT transformation tools get confused i 705 // G4-MT transformation tools get confused if statement spans two lines 707 // If needed contact: adotti@slac.stanford.e 706 // If needed contact: adotti@slac.stanford.edu 708 static G4ThreadLocal G4ParticleDefinition* p << 707 static G4ThreadLocal G4ParticleDefinition* pdefOpticalPhoton = 0; 709 if (pdefOpticalPhoton == nullptr) pdefOptica << 708 if (!pdefOpticalPhoton) pdefOpticalPhoton = 710 G4ParticleTable::GetParticleTable()->Fin 709 G4ParticleTable::GetParticleTable()->FindParticle("opticalphoton"); 711 710 712 static G4ThreadLocal G4int noCalls = 0; 711 static G4ThreadLocal G4int noCalls = 0; 713 noCalls++; 712 noCalls++; 714 713 715 fParticleChange.Initialize(track); 714 fParticleChange.Initialize(track); 716 715 717 // Code for specific process 716 // Code for specific process 718 // 717 // 719 fParticleChange.ProposePosition(State(fTrans 718 fParticleChange.ProposePosition(State(fTransportEndPosition)); 720 fParticleChange.ProposeMomentumDirection(Sta 719 fParticleChange.ProposeMomentumDirection(State(fTransportEndMomentumDir)); 721 fParticleChange.ProposeEnergy(State(fTranspo 720 fParticleChange.ProposeEnergy(State(fTransportEndKineticEnergy)); 722 fParticleChange.SetMomentumChanged(State(fMo 721 fParticleChange.SetMomentumChanged(State(fMomentumChanged)); 723 722 724 fParticleChange.ProposePolarization(State(fT 723 fParticleChange.ProposePolarization(State(fTransportEndSpin)); 725 724 726 G4double deltaTime = 0.0; 725 G4double deltaTime = 0.0; 727 726 728 // Calculate Lab Time of Flight (ONLY if fi 727 // Calculate Lab Time of Flight (ONLY if field Equations used it!) 729 // G4double endTime = State(fCandidateEndG 728 // G4double endTime = State(fCandidateEndGlobalTime); 730 // G4double delta_time = endTime - startTime 729 // G4double delta_time = endTime - startTime; 731 730 732 G4double startTime = track.GetGlobalTime(); 731 G4double startTime = track.GetGlobalTime(); 733 ///_________________________________________ 732 ///___________________________________________________________________________ 734 /// !!!!!!! 733 /// !!!!!!! 735 /// A REVOIR !!!! 734 /// A REVOIR !!!! 736 if (State(fEndGlobalTimeComputed) == 0) << 735 if (State(fEndGlobalTimeComputed) == false) 737 { 736 { 738 // The time was not integrated .. make the 737 // The time was not integrated .. make the best estimate possible 739 // 738 // 740 G4double initialVelocity = stepData.GetPre 739 G4double initialVelocity = stepData.GetPreStepPoint()->GetVelocity(); 741 G4double stepLength = track.GetStepLength( 740 G4double stepLength = track.GetStepLength(); 742 741 743 deltaTime = 0.0; // in case initialVelocit 742 deltaTime = 0.0; // in case initialVelocity = 0 744 if (track.GetParticleDefinition() == pdefO 743 if (track.GetParticleDefinition() == pdefOpticalPhoton) 745 { 744 { 746 // For only Optical Photon, final veloci 745 // For only Optical Photon, final velocity is used 747 double finalVelocity = track.CalculateVe 746 double finalVelocity = track.CalculateVelocityForOpticalPhoton(); 748 fParticleChange.ProposeVelocity(finalVel 747 fParticleChange.ProposeVelocity(finalVelocity); 749 deltaTime = stepLength / finalVelocity; 748 deltaTime = stepLength / finalVelocity; 750 } 749 } 751 else if (initialVelocity > 0.0) 750 else if (initialVelocity > 0.0) 752 { 751 { 753 deltaTime = stepLength / initialVelocity 752 deltaTime = stepLength / initialVelocity; 754 } 753 } 755 754 756 State(fCandidateEndGlobalTime) = startTime 755 State(fCandidateEndGlobalTime) = startTime + deltaTime; 757 } 756 } 758 else 757 else 759 { 758 { 760 deltaTime = State(fCandidateEndGlobalTime) 759 deltaTime = State(fCandidateEndGlobalTime) - startTime; 761 } 760 } 762 761 763 fParticleChange.ProposeGlobalTime(State(fCan 762 fParticleChange.ProposeGlobalTime(State(fCandidateEndGlobalTime)); 764 fParticleChange.ProposeLocalTime(track.GetLo 763 fParticleChange.ProposeLocalTime(track.GetLocalTime() + deltaTime); 765 /* 764 /* 766 // Now Correct by Lorentz factor to get del 765 // Now Correct by Lorentz factor to get delta "proper" Time 767 766 768 G4double restMass = track.GetDynamic 767 G4double restMass = track.GetDynamicParticle()->GetMass() ; 769 G4double deltaProperTime = deltaTime*( rest 768 G4double deltaProperTime = deltaTime*( restMass/track.GetTotalEnergy() ) ; 770 769 771 fParticleChange.ProposeProperTime(track.Get 770 fParticleChange.ProposeProperTime(track.GetProperTime() + deltaProperTime) ; 772 */ 771 */ 773 772 774 fParticleChange.ProposeTrueStepLength(track. 773 fParticleChange.ProposeTrueStepLength(track.GetStepLength()); 775 774 776 ///_________________________________________ 775 ///___________________________________________________________________________ 777 /// 776 /// 778 777 779 // If the particle is caught looping or is s 778 // If the particle is caught looping or is stuck (in very difficult 780 // boundaries) in a magnetic field (doing ma 779 // boundaries) in a magnetic field (doing many steps) 781 // THEN this kills it ... 780 // THEN this kills it ... 782 // 781 // 783 if (State(fParticleIsLooping)) 782 if (State(fParticleIsLooping)) 784 { 783 { 785 G4double endEnergy = State(fTransportEndKi 784 G4double endEnergy = State(fTransportEndKineticEnergy); 786 785 787 if ((endEnergy < fThreshold_Important_Ener 786 if ((endEnergy < fThreshold_Important_Energy) || (State(fNoLooperTrials) 788 >= fThresholdTrials)) 787 >= fThresholdTrials)) 789 { 788 { 790 // Kill the looping particle 789 // Kill the looping particle 791 // 790 // 792 // G4cout << "G4ITTransportation will ki 791 // G4cout << "G4ITTransportation will killed the molecule"<< G4endl; 793 fParticleChange.ProposeTrackStatus(fStop 792 fParticleChange.ProposeTrackStatus(fStopAndKill); 794 793 795 // 'Bare' statistics 794 // 'Bare' statistics 796 fSumEnergyKilled += endEnergy; 795 fSumEnergyKilled += endEnergy; 797 if (endEnergy > fMaxEnergyKilled) 796 if (endEnergy > fMaxEnergyKilled) 798 { 797 { 799 fMaxEnergyKilled = endEnergy; 798 fMaxEnergyKilled = endEnergy; 800 } 799 } 801 800 802 #ifdef G4VERBOSE 801 #ifdef G4VERBOSE 803 if ((fVerboseLevel > 1) || (endEnergy > 802 if ((fVerboseLevel > 1) || (endEnergy > fThreshold_Warning_Energy)) 804 { 803 { 805 G4cout 804 G4cout 806 << " G4ITTransportation is killing 805 << " G4ITTransportation is killing track that is looping or stuck " 807 << G4endl<< " This track has " < 806 << G4endl<< " This track has " << track.GetKineticEnergy() / MeV 808 << " MeV energy." << G4endl; 807 << " MeV energy." << G4endl; 809 G4cout << " Number of trials = " << 808 G4cout << " Number of trials = " << State(fNoLooperTrials) 810 << " No of calls to AlongStepDoIt = 809 << " No of calls to AlongStepDoIt = " << noCalls 811 << G4endl; 810 << G4endl; 812 } 811 } 813 #endif 812 #endif 814 State(fNoLooperTrials) = 0; 813 State(fNoLooperTrials) = 0; 815 } 814 } 816 else 815 else 817 { 816 { 818 State(fNoLooperTrials)++; 817 State(fNoLooperTrials)++; 819 #ifdef G4VERBOSE 818 #ifdef G4VERBOSE 820 if ((fVerboseLevel > 2)) 819 if ((fVerboseLevel > 2)) 821 { 820 { 822 G4cout << " G4ITTransportation::Alon 821 G4cout << " G4ITTransportation::AlongStepDoIt(): Particle looping - " 823 << " Number of trials = " << 822 << " Number of trials = " << State(fNoLooperTrials) 824 << " No of calls to = " << n 823 << " No of calls to = " << noCalls << G4endl; 825 } 824 } 826 #endif 825 #endif 827 } 826 } 828 } 827 } 829 else 828 else 830 { 829 { 831 State(fNoLooperTrials)=0; 830 State(fNoLooperTrials)=0; 832 } 831 } 833 832 834 // Another (sometimes better way) is to use 833 // Another (sometimes better way) is to use a user-limit maximum Step size 835 // to alleviate this problem .. 834 // to alleviate this problem .. 836 835 837 // Introduce smooth curved trajectories to p 836 // Introduce smooth curved trajectories to particle-change 838 // 837 // 839 fParticleChange.SetPointerToVectorOfAuxiliar 838 fParticleChange.SetPointerToVectorOfAuxiliaryPoints( 840 fFieldPropagator->GimmeTrajectoryVectorA 839 fFieldPropagator->GimmeTrajectoryVectorAndForgetIt()); 841 840 842 #if defined (DEBUG_MEM) 841 #if defined (DEBUG_MEM) 843 mem_second = MemoryUsage(); 842 mem_second = MemoryUsage(); 844 mem_diff = mem_second-mem_first; 843 mem_diff = mem_second-mem_first; 845 G4cout << "\t || MEM || End of G4ITTransport 844 G4cout << "\t || MEM || End of G4ITTransportation::AlongStepDoIt, diff is: " 846 << mem_diff << G4endl; 845 << mem_diff << G4endl; 847 #endif 846 #endif 848 847 849 return &fParticleChange; 848 return &fParticleChange; 850 } 849 } 851 850 852 ////////////////////////////////////////////// 851 ////////////////////////////////////////////////////////////////////////// 853 // 852 // 854 // This ensures that the PostStep action is a 853 // This ensures that the PostStep action is always called, 855 // so that it can do the relocation if it is 854 // so that it can do the relocation if it is needed. 856 // 855 // 857 856 858 G4double 857 G4double 859 G4ITTransportation:: 858 G4ITTransportation:: 860 PostStepGetPhysicalInteractionLength(const G4T 859 PostStepGetPhysicalInteractionLength(const G4Track&, // track 861 G4double, 860 G4double, // previousStepSize 862 G4ForceCo 861 G4ForceCondition* pForceCond) 863 { 862 { 864 *pForceCond = Forced; 863 *pForceCond = Forced; 865 return DBL_MAX; // was kInfinity ; but conve 864 return DBL_MAX; // was kInfinity ; but convention now is DBL_MAX 866 } 865 } 867 866 868 ////////////////////////////////////////////// 867 ///////////////////////////////////////////////////////////////////////////// 869 // 868 // 870 869 871 G4VParticleChange* G4ITTransportation::PostSte 870 G4VParticleChange* G4ITTransportation::PostStepDoIt(const G4Track& track, 872 871 const G4Step&) 873 { 872 { 874 // G4cout << "G4ITTransportation::PostSte 873 // G4cout << "G4ITTransportation::PostStepDoIt" << G4endl; 875 << 874 876 PrepareState(); << 875 PrepareState() 877 G4TouchableHandle retCurrentTouchable; // Th 876 G4TouchableHandle retCurrentTouchable; // The one to return 878 G4bool isLastStep = false; 877 G4bool isLastStep = false; 879 878 880 // Initialize ParticleChange (by setting al 879 // Initialize ParticleChange (by setting all its members equal 881 // to correspond 880 // to corresponding members in G4Track) 882 fParticleChange.Initialize(track); // To ini 881 fParticleChange.Initialize(track); // To initialise TouchableChange 883 882 884 fParticleChange.ProposeTrackStatus(track.Get 883 fParticleChange.ProposeTrackStatus(track.GetTrackStatus()); 885 884 886 // If the Step was determined by the volume 885 // If the Step was determined by the volume boundary, 887 // logically relocate the particle 886 // logically relocate the particle 888 887 889 if (State(fGeometryLimitedStep)) 888 if (State(fGeometryLimitedStep)) 890 { 889 { 891 890 892 if(fVerboseLevel != 0) << 891 if(fVerboseLevel) 893 { 892 { 894 G4cout << "Step is limited by geometry " 893 G4cout << "Step is limited by geometry " 895 << "track ID : " << track.GetTrac 894 << "track ID : " << track.GetTrackID() << G4endl; 896 } 895 } 897 896 898 // fCurrentTouchable will now become the p 897 // fCurrentTouchable will now become the previous touchable, 899 // and what was the previous will be freed 898 // and what was the previous will be freed. 900 // (Needed because the preStepPoint can po 899 // (Needed because the preStepPoint can point to the previous touchable) 901 900 902 if ( State(fCurrentTouchableHandle)->GetVo << 901 if ( State(fCurrentTouchableHandle)->GetVolume() == 0) 903 { 902 { 904 G4ExceptionDescription exceptionDescript 903 G4ExceptionDescription exceptionDescription; 905 exceptionDescription << "No current touc 904 exceptionDescription << "No current touchable found "; 906 G4Exception(" G4ITTransportation::PostSt 905 G4Exception(" G4ITTransportation::PostStepDoIt", "G4ITTransportation001", 907 FatalErrorInArgument, except 906 FatalErrorInArgument, exceptionDescription); 908 } 907 } 909 908 910 fLinearNavigator->SetGeometricallyLimitedS 909 fLinearNavigator->SetGeometricallyLimitedStep(); 911 fLinearNavigator->LocateGlobalPointAndUpda 910 fLinearNavigator->LocateGlobalPointAndUpdateTouchableHandle( 912 track.GetPosition(), track.GetMomentum 911 track.GetPosition(), track.GetMomentumDirection(), 913 State(fCurrentTouchableHandle), true); 912 State(fCurrentTouchableHandle), true); 914 // Check whether the particle is out of th 913 // Check whether the particle is out of the world volume 915 // If so it has exited and must be killed. 914 // If so it has exited and must be killed. 916 // 915 // 917 if ( State(fCurrentTouchableHandle)->GetVo << 916 if ( State(fCurrentTouchableHandle)->GetVolume() == 0) 918 { 917 { 919 // abort(); 918 // abort(); 920 #ifdef G4VERBOSE 919 #ifdef G4VERBOSE 921 if (fVerboseLevel > 0) 920 if (fVerboseLevel > 0) 922 { 921 { 923 G4cout << "Track position : " << track 922 G4cout << "Track position : " << track.GetPosition() / nanometer 924 << " [nm]" << " Track ID : " << 923 << " [nm]" << " Track ID : " << track.GetTrackID() << G4endl; 925 G4cout << "G4ITTransportation will kil 924 G4cout << "G4ITTransportation will killed the track because " 926 "State(fCurrentTouchableHandle)->G 925 "State(fCurrentTouchableHandle)->GetVolume() == 0"<< G4endl; 927 } 926 } 928 #endif 927 #endif 929 fParticleChange.ProposeTrackStatus( fSto 928 fParticleChange.ProposeTrackStatus( fStopAndKill ); 930 } 929 } 931 930 932 retCurrentTouchable = State(fCurrentToucha 931 retCurrentTouchable = State(fCurrentTouchableHandle); 933 932 934 // G4cout << "Current volume : " << tra 933 // G4cout << "Current volume : " << track.GetVolume()->GetName() 935 // << " Next volume : " 934 // << " Next volume : " 936 // << (State(fCurrentTouchableHa 935 // << (State(fCurrentTouchableHandle)->GetVolume() ? 937 // State(fCurrentTouchableHandle)->GetVolum 936 // State(fCurrentTouchableHandle)->GetVolume()->GetName():"OutWorld") 938 // << " Position : " << track.Ge 937 // << " Position : " << track.GetPosition() / nanometer 939 // << " track ID : " << track.Ge 938 // << " track ID : " << track.GetTrackID() 940 // << G4endl; 939 // << G4endl; 941 940 942 fParticleChange.SetTouchableHandle(State(f 941 fParticleChange.SetTouchableHandle(State(fCurrentTouchableHandle)); 943 942 944 // Update the Step flag which identifies t 943 // Update the Step flag which identifies the Last Step in a volume 945 isLastStep = fLinearNavigator->ExitedMothe 944 isLastStep = fLinearNavigator->ExitedMotherVolume() 946 || fLinearNavigator->EnteredDau << 945 | fLinearNavigator->EnteredDaughterVolume(); 947 946 948 #ifdef G4DEBUG_TRANSPORT 947 #ifdef G4DEBUG_TRANSPORT 949 // Checking first implementation of flagg 948 // Checking first implementation of flagging Last Step in Volume 950 G4bool exiting = fLinearNavigator->ExitedM 949 G4bool exiting = fLinearNavigator->ExitedMotherVolume(); 951 G4bool entering = fLinearNavigator->Entere 950 G4bool entering = fLinearNavigator->EnteredDaughterVolume(); 952 951 953 if( ! (exiting || entering) ) 952 if( ! (exiting || entering) ) 954 { 953 { 955 G4cout << " Transport> : Proposed isLas 954 G4cout << " Transport> : Proposed isLastStep= " << isLastStep 956 << " Exiting " << fLinearNavigator->Exit 955 << " Exiting " << fLinearNavigator->ExitedMotherVolume() 957 << " Entering " << fLinearNavigator->Ent 956 << " Entering " << fLinearNavigator->EnteredDaughterVolume() 958 << " Track position : " << track.GetPosi 957 << " Track position : " << track.GetPosition() /nanometer << " [nm]" 959 << G4endl; 958 << G4endl; 960 G4cout << " Track position : " << track. 959 G4cout << " Track position : " << track.GetPosition() /nanometer 961 << G4endl; 960 << G4endl; 962 } 961 } 963 #endif 962 #endif 964 } 963 } 965 else // fGeometryLimitedStep is false 964 else // fGeometryLimitedStep is false 966 { 965 { 967 // This serves only to move the Navigator' 966 // This serves only to move the Navigator's location 968 // 967 // 969 // abort(); 968 // abort(); 970 fLinearNavigator->LocateGlobalPointWithinV 969 fLinearNavigator->LocateGlobalPointWithinVolume(track.GetPosition()); 971 970 972 // The value of the track's current Toucha 971 // The value of the track's current Touchable is retained. 973 // (and it must be correct because we must 972 // (and it must be correct because we must use it below to 974 // overwrite the (unset) one in particle c 973 // overwrite the (unset) one in particle change) 975 // It must be fCurrentTouchable too ?? 974 // It must be fCurrentTouchable too ?? 976 // 975 // 977 fParticleChange.SetTouchableHandle(track.G 976 fParticleChange.SetTouchableHandle(track.GetTouchableHandle()); 978 retCurrentTouchable = track.GetTouchableHa 977 retCurrentTouchable = track.GetTouchableHandle(); 979 978 980 isLastStep = false; 979 isLastStep = false; 981 #ifdef G4DEBUG_TRANSPORT 980 #ifdef G4DEBUG_TRANSPORT 982 // Checking first implementation of flagg 981 // Checking first implementation of flagging Last Step in Volume 983 G4cout << " Transport> Proposed isLastStep 982 G4cout << " Transport> Proposed isLastStep= " << isLastStep 984 << " Geometry did not limit step. Position 983 << " Geometry did not limit step. Position : " 985 << track.GetPosition()/ nanometer << G4end 984 << track.GetPosition()/ nanometer << G4endl; 986 #endif 985 #endif 987 } // endif ( fGeometryLimitedStep ) 986 } // endif ( fGeometryLimitedStep ) 988 987 989 fParticleChange.ProposeLastStepInVolume(isLa 988 fParticleChange.ProposeLastStepInVolume(isLastStep); 990 989 991 const G4VPhysicalVolume* pNewVol = retCurren 990 const G4VPhysicalVolume* pNewVol = retCurrentTouchable->GetVolume(); 992 const G4Material* pNewMaterial = nullptr; << 991 const G4Material* pNewMaterial = 0; 993 G4VSensitiveDetector* pNewSensitiveDetector << 992 const G4VSensitiveDetector* pNewSensitiveDetector = 0; 994 993 995 if (pNewVol != nullptr) << 994 if (pNewVol != 0) 996 { 995 { 997 pNewMaterial = pNewVol->GetLogicalVolume() 996 pNewMaterial = pNewVol->GetLogicalVolume()->GetMaterial(); 998 pNewSensitiveDetector = pNewVol->GetLogica 997 pNewSensitiveDetector = pNewVol->GetLogicalVolume()->GetSensitiveDetector(); 999 } 998 } 1000 999 1001 // ( <const_cast> pNewMaterial ) ; 1000 // ( <const_cast> pNewMaterial ) ; >> 1001 // ( <const_cast> pNewSensitiveDetector) ; 1002 1002 1003 fParticleChange.SetMaterialInTouchable((G4M 1003 fParticleChange.SetMaterialInTouchable((G4Material *) pNewMaterial); 1004 fParticleChange.SetSensitiveDetectorInTouch << 1004 fParticleChange.SetSensitiveDetectorInTouchable( >> 1005 (G4VSensitiveDetector *) pNewSensitiveDetector); 1005 1006 1006 const G4MaterialCutsCouple* pNewMaterialCut << 1007 const G4MaterialCutsCouple* pNewMaterialCutsCouple = 0; 1007 if (pNewVol != nullptr) << 1008 if (pNewVol != 0) 1008 { 1009 { 1009 pNewMaterialCutsCouple = 1010 pNewMaterialCutsCouple = 1010 pNewVol->GetLogicalVolume()->GetMater 1011 pNewVol->GetLogicalVolume()->GetMaterialCutsCouple(); 1011 } 1012 } 1012 1013 1013 if (pNewVol != nullptr && pNewMaterialCutsC << 1014 if (pNewVol != 0 && pNewMaterialCutsCouple != 0 1014 && pNewMaterialCutsCouple->GetMaterial( 1015 && pNewMaterialCutsCouple->GetMaterial() != pNewMaterial) 1015 { 1016 { 1016 // for parametrized volume 1017 // for parametrized volume 1017 // 1018 // 1018 pNewMaterialCutsCouple = G4ProductionCuts 1019 pNewMaterialCutsCouple = G4ProductionCutsTable::GetProductionCutsTable() 1019 ->GetMaterialCutsCouple(pNewMaterial, 1020 ->GetMaterialCutsCouple(pNewMaterial, 1020 pNewMaterialC 1021 pNewMaterialCutsCouple->GetProductionCuts()); 1021 } 1022 } 1022 fParticleChange.SetMaterialCutsCoupleInTouc 1023 fParticleChange.SetMaterialCutsCoupleInTouchable(pNewMaterialCutsCouple); 1023 1024 1024 // temporarily until Get/Set Material of Pa 1025 // temporarily until Get/Set Material of ParticleChange, 1025 // and StepPoint can be made const. 1026 // and StepPoint can be made const. 1026 // Set the touchable in ParticleChange 1027 // Set the touchable in ParticleChange 1027 // this must always be done because the par 1028 // this must always be done because the particle change always 1028 // uses this value to overwrite the current 1029 // uses this value to overwrite the current touchable pointer. 1029 // 1030 // 1030 fParticleChange.SetTouchableHandle(retCurre 1031 fParticleChange.SetTouchableHandle(retCurrentTouchable); 1031 1032 1032 return &fParticleChange; 1033 return &fParticleChange; 1033 } 1034 } 1034 1035 1035 // New method takes over the responsibility t 1036 // New method takes over the responsibility to reset the state of 1036 // G4Transportation object at the start of a 1037 // G4Transportation object at the start of a new track or the resumption of 1037 // a suspended track. 1038 // a suspended track. 1038 1039 1039 void G4ITTransportation::StartTracking(G4Trac 1040 void G4ITTransportation::StartTracking(G4Track* track) 1040 { 1041 { 1041 G4VProcess::StartTracking(track); 1042 G4VProcess::StartTracking(track); 1042 if (fInstantiateProcessState) 1043 if (fInstantiateProcessState) 1043 { 1044 { 1044 // G4VITProcess::fpState = new G4ITTra 1045 // G4VITProcess::fpState = new G4ITTransportationState(); 1045 G4VITProcess::fpState = std::make_shared< << 1046 G4VITProcess::fpState.reset(new G4ITTransportationState()); 1046 // Will set in the same time fTransportat 1047 // Will set in the same time fTransportationState 1047 } 1048 } 1048 1049 1049 fpSafetyHelper->NewTrackState(); 1050 fpSafetyHelper->NewTrackState(); 1050 fpSafetyHelper->SaveTrackState( 1051 fpSafetyHelper->SaveTrackState( 1051 GetIT(track)->GetTrackingInfo()->GetTra 1052 GetIT(track)->GetTrackingInfo()->GetTrackStateManager()); 1052 1053 1053 // The actions here are those that were tak 1054 // The actions here are those that were taken in AlongStepGPIL 1054 // when track.GetCurrentStepNumber()==1 1055 // when track.GetCurrentStepNumber()==1 1055 1056 1056 // reset safety value and center 1057 // reset safety value and center 1057 // 1058 // 1058 // State(fPreviousSafety) = 0.0 ; 1059 // State(fPreviousSafety) = 0.0 ; 1059 // State(fPreviousSftOrigin) = G4ThreeVe 1060 // State(fPreviousSftOrigin) = G4ThreeVector(0.,0.,0.) ; 1060 1061 1061 // reset looping counter -- for motion in f 1062 // reset looping counter -- for motion in field 1062 // State(fNoLooperTrials)= 0; 1063 // State(fNoLooperTrials)= 0; 1063 // Must clear this state .. else it depends 1064 // Must clear this state .. else it depends on last track's value 1064 // --> a better solution would set this fr 1065 // --> a better solution would set this from state of suspended track TODO ? 1065 // Was if( aTrack->GetCurrentStepNumber()== 1066 // Was if( aTrack->GetCurrentStepNumber()==1 ) { .. } 1066 1067 1067 // ChordFinder reset internal state 1068 // ChordFinder reset internal state 1068 // 1069 // 1069 if (DoesGlobalFieldExist()) 1070 if (DoesGlobalFieldExist()) 1070 { 1071 { 1071 fFieldPropagator->ClearPropagatorState(); 1072 fFieldPropagator->ClearPropagatorState(); 1072 // Resets all state of field propagator c 1073 // Resets all state of field propagator class (ONLY) 1073 // including safety values (in case of ov 1074 // including safety values (in case of overlaps and to wipe for first track). 1074 1075 1075 // G4ChordFinder* chordF= fFieldPropagato 1076 // G4ChordFinder* chordF= fFieldPropagator->GetChordFinder(); 1076 // if( chordF ) chordF->ResetStepEstimate 1077 // if( chordF ) chordF->ResetStepEstimate(); 1077 } 1078 } 1078 1079 1079 // Make sure to clear the chord finders of 1080 // Make sure to clear the chord finders of all fields (ie managers) 1080 static G4ThreadLocal G4FieldManagerStore* f << 1081 static G4ThreadLocal G4FieldManagerStore* fieldMgrStore = 0; 1081 if (fieldMgrStore == nullptr) fieldMgrStore << 1082 if (!fieldMgrStore) fieldMgrStore = G4FieldManagerStore::GetInstance(); 1082 fieldMgrStore->ClearAllChordFindersState(); 1083 fieldMgrStore->ClearAllChordFindersState(); 1083 1084 1084 // Update the current touchable handle (fr 1085 // Update the current touchable handle (from the track's) 1085 // 1086 // 1086 PrepareState(); << 1087 PrepareState() 1087 State(fCurrentTouchableHandle) = track->Get 1088 State(fCurrentTouchableHandle) = track->GetTouchableHandle(); 1088 1089 1089 G4VITProcess::StartTracking(track); 1090 G4VITProcess::StartTracking(track); 1090 } 1091 } 1091 1092 1092 #undef State 1093 #undef State 1093 #undef PrepareState 1094 #undef PrepareState 1094 1095