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 exoticphysics/monopole/src/G4Monopol << 26 // $Id: G4MonopoleTransportation.cc,v 1.2 2010-11-29 15:14:17 vnivanch Exp $ 27 /// \brief Implementation of the G4MonopoleTra << 27 // GEANT4 tag $Name: not supported by cvs2svn $ 28 // << 29 // 28 // 30 //....oooOO0OOooo........oooOO0OOooo........oo 29 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 31 //....oooOO0OOooo........oooOO0OOooo........oo 30 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 32 // 31 // 33 // This class is a process responsible for the << 32 // This class is a process responsible for the transportation of 34 // magnetic monopoles, ie the geometrical prop << 33 // magnetic monopoles, ie the geometrical propagation that encounters the 35 // geometrical sub-volumes of the detectors. 34 // geometrical sub-volumes of the detectors. 36 // 35 // 37 // For monopoles, uses a different equation of 36 // For monopoles, uses a different equation of motion and ignores energy 38 // conservation. << 37 // conservation. 39 // 38 // 40 39 41 // =========================================== 40 // ======================================================================= 42 // Created: 3 May 2010, J. Apostolakis, B. Bo 41 // Created: 3 May 2010, J. Apostolakis, B. Bozsogi 43 // =========================================== 42 // ======================================================================= 44 43 45 #include "G4MonopoleTransportation.hh" 44 #include "G4MonopoleTransportation.hh" 46 << 45 #include "G4ProductionCutsTable.hh" 47 #include "DetectorConstruction.hh" << 46 #include "G4ParticleTable.hh" 48 << 49 #include "G4ChordFinder.hh" 47 #include "G4ChordFinder.hh" >> 48 #include "G4SafetyHelper.hh" 50 #include "G4FieldManagerStore.hh" 49 #include "G4FieldManagerStore.hh" 51 #include "G4Monopole.hh" 50 #include "G4Monopole.hh" 52 #include "G4ParticleTable.hh" << 53 #include "G4ProductionCutsTable.hh" << 54 #include "G4RunManager.hh" << 55 #include "G4SafetyHelper.hh" << 56 #include "G4SystemOfUnits.hh" << 57 #include "G4TransportationProcessType.hh" << 58 51 59 class G4VSensitiveDetector; 52 class G4VSensitiveDetector; 60 53 61 //....oooOO0OOooo........oooOO0OOooo........oo << 54 ////////////////////////////////////////////////////////////////////////// >> 55 // >> 56 // Constructor 62 57 63 G4MonopoleTransportation::G4MonopoleTransporta << 58 G4MonopoleTransportation::G4MonopoleTransportation( const G4Monopole* mpl, 64 : G4VProcess(G4String("MonopoleTransportatio << 59 G4int verb) 65 fParticleDef(mpl), << 60 : G4VProcess( G4String("MonopoleTransportation"), fTransportation ), 66 fMagSetup(0), << 61 fParticleIsLooping( false ), 67 fLinearNavigator(0), << 62 fPreviousSftOrigin (0.,0.,0.), 68 fFieldPropagator(0), << 63 fPreviousSafety ( 0.0 ), 69 fParticleIsLooping(false), << 64 fThreshold_Warning_Energy( 100 * MeV ), 70 fPreviousSftOrigin(0., 0., 0.), << 65 fThreshold_Important_Energy( 250 * MeV ), 71 fPreviousSafety(0.0), << 66 fThresholdTrials( 10 ), 72 fThreshold_Warning_Energy(100 * MeV), << 67 fUnimportant_Energy( 1 * MeV ), 73 fThreshold_Important_Energy(250 * MeV), << 74 fThresholdTrials(10), << 75 // fUnimportant_Energy( 1 * MeV ), << 76 fNoLooperTrials(0), 68 fNoLooperTrials(0), 77 fSumEnergyKilled(0.0), << 69 fSumEnergyKilled( 0.0 ), fMaxEnergyKilled( 0.0 ), 78 fMaxEnergyKilled(0.0), << 70 fShortStepOptimisation(false), // Old default: true (=fast short steps) 79 fShortStepOptimisation(false), // Old def << 71 fVerboseLevel( verb ) 80 fpSafetyHelper(0), << 81 noCalls(0) << 82 { 72 { 83 verboseLevel = verb; << 73 pParticleDef = mpl; >> 74 >> 75 magSetup = G4MonopoleFieldSetup::GetMonopoleFieldSetup(); >> 76 >> 77 G4TransportationManager* transportMgr ; 84 78 85 // set Process Sub Type << 79 transportMgr = G4TransportationManager::GetTransportationManager() ; 86 SetProcessSubType(TRANSPORTATION); << 87 80 88 // Do not finalize the G4MonopoleTransportat << 81 fLinearNavigator = transportMgr->GetNavigatorForTracking() ; 89 if (G4Threading::IsMasterThread() && G4Threa << 90 return; << 91 } << 92 << 93 const DetectorConstruction* detector = stati << 94 G4RunManager::GetRunManager()->GetUserDete << 95 << 96 fMagSetup = detector->GetMonopoleFieldSetup( << 97 82 98 G4TransportationManager* transportMgr = G4Tr << 83 // fGlobalFieldMgr = transportMgr->GetFieldManager() ; 99 84 100 fLinearNavigator = transportMgr->GetNavigato << 85 fFieldPropagator = transportMgr->GetPropagatorInField() ; 101 86 102 fFieldPropagator = transportMgr->GetPropagat << 87 fpSafetyHelper = transportMgr->GetSafetyHelper(); // New 103 fpSafetyHelper = transportMgr->GetSafetyHelp << 104 << 105 // New << 106 88 107 // Cannot determine whether a field exists h 89 // Cannot determine whether a field exists here, 108 // because it would only work if the field << 90 // because it would only work if the field manager has informed 109 // about the detector's field before this t << 91 // about the detector's field before this transportation process 110 // is constructed. 92 // is constructed. 111 // Instead later the method DoesGlobalFieldE 93 // Instead later the method DoesGlobalFieldExist() is called 112 fCurrentTouchableHandle = nullptr; << 113 94 114 fEndGlobalTimeComputed = false; << 95 static G4TouchableHandle nullTouchableHandle; // Points to (G4VTouchable*) 0 >> 96 fCurrentTouchableHandle = nullTouchableHandle; >> 97 >> 98 fEndGlobalTimeComputed = false; 115 fCandidateEndGlobalTime = 0; 99 fCandidateEndGlobalTime = 0; 116 } 100 } 117 101 118 //....oooOO0OOooo........oooOO0OOooo........oo << 102 ////////////////////////////////////////////////////////////////////////// 119 103 120 G4MonopoleTransportation::~G4MonopoleTransport 104 G4MonopoleTransportation::~G4MonopoleTransportation() 121 { 105 { 122 if ((verboseLevel > 0) && (fSumEnergyKilled << 106 if( (fVerboseLevel > 0) && (fSumEnergyKilled > 0.0 ) ){ 123 G4cout << " G4MonopoleTransportation: Stat 107 G4cout << " G4MonopoleTransportation: Statistics for looping particles " << G4endl; 124 G4cout << " Sum of energy of loopers kil << 108 G4cout << " Sum of energy of loopers killed: " << fSumEnergyKilled << G4endl; 125 G4cout << " Max energy of loopers killed << 109 G4cout << " Max energy of loopers killed: " << fMaxEnergyKilled << G4endl; 126 } << 110 } 127 } 111 } 128 112 129 //....oooOO0OOooo........oooOO0OOooo........oo << 113 ////////////////////////////////////////////////////////////////////////// 130 // 114 // 131 // Responsibilities: 115 // Responsibilities: 132 // Find whether the geometry limits the Ste 116 // Find whether the geometry limits the Step, and to what length 133 // Calculate the new value of the safety an 117 // Calculate the new value of the safety and return it. 134 // Store the final time, position and momen 118 // Store the final time, position and momentum. 135 119 136 G4double G4MonopoleTransportation::AlongStepGe << 120 G4double G4MonopoleTransportation:: 137 const G4Track& track, << 121 AlongStepGetPhysicalInteractionLength( const G4Track& track, 138 G4double, // previousStepSize << 122 G4double, // previousStepSize 139 G4double currentMinimumStep, G4double& curre << 123 G4double currentMinimumStep, 140 { << 124 G4double& currentSafety, 141 fMagSetup->SetStepperAndChordFinder(1); << 125 G4GPILSelection* selection ) >> 126 { >> 127 magSetup->SetStepperAndChordFinder(1); 142 // change to monopole equation 128 // change to monopole equation 143 129 144 G4double geometryStepLength, newSafety; << 130 G4double geometryStepLength, newSafety ; 145 fParticleIsLooping = false; << 131 fParticleIsLooping = false ; 146 132 147 // Initial actions moved to StartTrack() << 133 // Initial actions moved to StartTrack() 148 // -------------------------------------- 134 // -------------------------------------- 149 // Note: in case another process changes tou 135 // Note: in case another process changes touchable handle 150 // it will be necessary to add here (for << 136 // it will be necessary to add here (for all steps) 151 // fCurrentTouchableHandle = aTrack->GetTouc 137 // fCurrentTouchableHandle = aTrack->GetTouchableHandle(); 152 138 153 // GPILSelection is set to defaule value of 139 // GPILSelection is set to defaule value of CandidateForSelection 154 // It is a return value 140 // It is a return value 155 // 141 // 156 *selection = CandidateForSelection; << 142 *selection = CandidateForSelection ; 157 143 158 // Get initial Energy/Momentum of the track 144 // Get initial Energy/Momentum of the track 159 // 145 // 160 const G4DynamicParticle* pParticle = track.G << 146 const G4DynamicParticle* pParticle = track.GetDynamicParticle() ; 161 G4ThreeVector startMomentumDir = pParticle-> << 147 G4ThreeVector startMomentumDir = pParticle->GetMomentumDirection() ; 162 G4ThreeVector startPosition = track.GetPosit << 148 G4ThreeVector startPosition = track.GetPosition() ; 163 149 164 // G4double theTime = track.GetGlob 150 // G4double theTime = track.GetGlobalTime() ; 165 151 166 // The Step Point safety can be limited by o << 152 // The Step Point safety can be limited by other geometries and/or the 167 // assumptions of any process - it's not alw 153 // assumptions of any process - it's not always the geometrical safety. 168 // We calculate the starting point's isotrop 154 // We calculate the starting point's isotropic safety here. 169 // 155 // 170 G4ThreeVector OriginShift = startPosition - << 156 G4ThreeVector OriginShift = startPosition - fPreviousSftOrigin ; 171 G4double MagSqShift = OriginShift.mag2(); << 157 G4double MagSqShift = OriginShift.mag2() ; 172 if (MagSqShift >= sqr(fPreviousSafety)) { << 158 if( MagSqShift >= sqr(fPreviousSafety) ) 173 currentSafety = 0.0; << 159 { >> 160 currentSafety = 0.0 ; 174 } 161 } 175 else { << 162 else 176 currentSafety = fPreviousSafety - std::sqr << 163 { >> 164 currentSafety = fPreviousSafety - std::sqrt(MagSqShift) ; 177 } 165 } 178 166 179 // Is the monopole charged ? 167 // Is the monopole charged ? 180 // 168 // 181 G4double particleMagneticCharge = fParticleD << 169 G4double particleMagneticCharge = pParticleDef->MagneticCharge() ; 182 G4double particleElectricCharge = pParticle- 170 G4double particleElectricCharge = pParticle->GetCharge(); 183 171 184 fGeometryLimitedStep = false; << 172 fGeometryLimitedStep = false ; 185 // fEndGlobalTimeComputed = false ; 173 // fEndGlobalTimeComputed = false ; 186 174 187 // There is no need to locate the current vo 175 // There is no need to locate the current volume. It is Done elsewhere: 188 // On track construction << 176 // On track construction 189 // By the tracking, after all AlongStepDoI 177 // By the tracking, after all AlongStepDoIts, in "Relocation" 190 178 191 // Check whether the particle have an (EM) f 179 // Check whether the particle have an (EM) field force exerting upon it 192 // 180 // 193 G4FieldManager* fieldMgr = 0; << 181 G4FieldManager* fieldMgr=0; 194 G4bool fieldExertsForce = false; << 182 G4bool fieldExertsForce = false ; 195 << 183 196 if ((particleMagneticCharge != 0.0)) { << 184 if( (particleMagneticCharge != 0.0) ) 197 fieldMgr = fFieldPropagator->FindAndSetFie << 185 { 198 if (fieldMgr != 0) { << 186 fieldMgr= fFieldPropagator->FindAndSetFieldManager( track.GetVolume() ); 199 // Message the field Manager, to configu << 187 if (fieldMgr != 0) { 200 fieldMgr->ConfigureForTrack(&track); << 188 // Message the field Manager, to configure it for this track 201 // Moved here, in order to allow a trans << 189 fieldMgr->ConfigureForTrack( &track ); 202 // from a zero-field status (with fie << 190 // Moved here, in order to allow a transition 203 // to a finite field status << 191 // from a zero-field status (with fieldMgr->(field)0 204 << 192 // to a finite field status 205 // If the field manager has no field, th << 193 206 fieldExertsForce = (fieldMgr->GetDetecto << 194 // If the field manager has no field, there is no field ! 207 } << 195 fieldExertsForce = (fieldMgr->GetDetectorField() != 0); >> 196 } 208 } 197 } 209 198 210 // G4cout << " G4Transport: field exerts fo 199 // G4cout << " G4Transport: field exerts force= " << fieldExertsForce 211 // << " fieldMgr= " << fieldMgr << << 200 // << " fieldMgr= " << fieldMgr << G4endl; 212 201 213 // Choose the calculation of the transportat << 202 // Choose the calculation of the transportation: Field or not 214 // 203 // 215 if (!fieldExertsForce) { << 204 if( !fieldExertsForce ) 216 G4double linearStepLength; << 205 { 217 if (fShortStepOptimisation && (currentMini << 206 G4double linearStepLength ; 218 // The Step is guaranteed to be taken << 207 if( fShortStepOptimisation && (currentMinimumStep <= currentSafety) ) 219 // << 208 { 220 geometryStepLength = currentMinimumStep; << 209 // The Step is guaranteed to be taken 221 fGeometryLimitedStep = false; << 210 // 222 } << 211 geometryStepLength = currentMinimumStep ; 223 else { << 212 fGeometryLimitedStep = false ; 224 // Find whether the straight path inter << 213 } 225 // << 214 else 226 linearStepLength = fLinearNavigator->Com << 215 { 227 << 216 // Find whether the straight path intersects a volume 228 // Remember last safety origin & value. << 217 // 229 // << 218 linearStepLength = fLinearNavigator->ComputeStep( startPosition, 230 fPreviousSftOrigin = startPosition; << 219 startMomentumDir, 231 fPreviousSafety = newSafety; << 220 currentMinimumStep, 232 // fpSafetyHelper->SetCurrentSafety( new << 221 newSafety) ; 233 << 222 // Remember last safety origin & value. 234 // The safety at the initial point has b << 223 // 235 // << 224 fPreviousSftOrigin = startPosition ; 236 currentSafety = newSafety; << 225 fPreviousSafety = newSafety ; 237 << 226 // fpSafetyHelper->SetCurrentSafety( newSafety, startPosition); 238 fGeometryLimitedStep = (linearStepLength << 227 239 if (fGeometryLimitedStep) { << 228 // The safety at the initial point has been re-calculated: 240 // The geometry limits the Step size ( << 229 // 241 geometryStepLength = linearStepLength; << 230 currentSafety = newSafety ; 242 } << 231 243 else { << 232 fGeometryLimitedStep= (linearStepLength <= currentMinimumStep); 244 // The full Step is taken. << 233 if( fGeometryLimitedStep ) 245 geometryStepLength = currentMinimumSte << 234 { 246 } << 235 // The geometry limits the Step size (an intersection was found.) 247 } << 236 geometryStepLength = linearStepLength ; 248 endpointDistance = geometryStepLength; << 237 } 249 << 238 else 250 // Calculate final position << 239 { 251 // << 240 // The full Step is taken. 252 fTransportEndPosition = startPosition + ge << 241 geometryStepLength = currentMinimumStep ; 253 << 242 } 254 // Momentum direction, energy and polarisa << 243 } 255 // << 244 endpointDistance = geometryStepLength ; 256 fTransportEndMomentumDir = startMomentumDi << 245 257 fTransportEndKineticEnergy = track.GetKine << 246 // Calculate final position 258 fTransportEndSpin = track.GetPolarization( << 247 // 259 fParticleIsLooping = false; << 248 fTransportEndPosition = startPosition+geometryStepLength*startMomentumDir ; 260 fMomentumChanged = false; << 249 261 fEndGlobalTimeComputed = false; << 250 // Momentum direction, energy and polarisation are unchanged by transport 262 } << 251 // 263 else // A field exerts force << 252 fTransportEndMomentumDir = startMomentumDir ; 264 { << 253 fTransportEndKineticEnergy = track.GetKineticEnergy() ; 265 G4double momentumMagnitude = pParticle->Ge << 254 fTransportEndSpin = track.GetPolarization(); 266 G4ThreeVector EndUnitMomentum; << 255 fParticleIsLooping = false ; 267 G4double lengthAlongCurve; << 256 fMomentumChanged = false ; 268 G4double restMass = fParticleDef->GetPDGMa << 257 fEndGlobalTimeComputed = false ; 269 << 258 } 270 G4ChargeState chargeState(particleElectric << 259 else // A field exerts force 271 fParticleDef->Ge << 260 { 272 0, // Magnetic << 261 // G4double momentumMagnitude = pParticle->GetTotalMomentum() ; 273 0, // Electric << 262 G4ThreeVector EndUnitMomentum ; 274 particleMagnetic << 263 G4double lengthAlongCurve ; 275 << 264 G4double restMass = pParticleDef->GetPDGMass() ; 276 G4EquationOfMotion* equationOfMotion = << 265 277 fFieldPropagator->GetChordFinder()->GetI << 266 fFieldPropagator->SetChargeMomentumMass( particleMagneticCharge, // in Mev/c 278 << 267 particleElectricCharge, // in e+ units 279 equationOfMotion->SetChargeMomentumMass(ch << 268 restMass ) ; 280 // SetChargeMomentumMass now passes both t << 269 281 // charge in chargeState << 270 // SetChargeMomentumMass is _not_ used here as it would in everywhere else, 282 << 271 // it's just a workaround to pass the electric charge as well. 283 G4ThreeVector spin = track.GetPolarization << 272 284 G4FieldTrack aFieldTrack = G4FieldTrack(st << 273 285 tr << 274 G4ThreeVector spin = track.GetPolarization() ; 286 tr << 275 G4FieldTrack aFieldTrack = G4FieldTrack( startPosition, 287 tr << 276 track.GetMomentumDirection(), 288 &s << 277 0.0, 289 if (currentMinimumStep > 0) { << 278 track.GetKineticEnergy(), 290 // Do the Transport in the field (non re << 279 restMass, 291 // << 280 track.GetVelocity(), 292 lengthAlongCurve = fFieldPropagator->Com << 281 track.GetGlobalTime(), // Lab. 293 << 282 track.GetProperTime(), // Part. 294 fGeometryLimitedStep = lengthAlongCurve << 283 &spin ) ; 295 if (fGeometryLimitedStep) { << 284 if( currentMinimumStep > 0 ) 296 geometryStepLength = lengthAlongCurve; << 285 { 297 } << 286 // Do the Transport in the field (non recti-linear) 298 else { << 287 // 299 geometryStepLength = currentMinimumSte << 288 lengthAlongCurve = fFieldPropagator->ComputeStep( aFieldTrack, 300 } << 289 currentMinimumStep, 301 } << 290 currentSafety, 302 else { << 291 track.GetVolume() ) ; 303 geometryStepLength = lengthAlongCurve = << 292 fGeometryLimitedStep= lengthAlongCurve < currentMinimumStep; 304 fGeometryLimitedStep = false; << 293 if( fGeometryLimitedStep ) { 305 } << 294 geometryStepLength = lengthAlongCurve ; 306 << 295 } else { 307 // Remember last safety origin & value. << 296 geometryStepLength = currentMinimumStep ; 308 // << 297 } 309 fPreviousSftOrigin = startPosition; << 298 } 310 fPreviousSafety = currentSafety; << 299 else 311 // fpSafetyHelper->SetCurrentSafety( newSa << 300 { 312 << 301 geometryStepLength = lengthAlongCurve= 0.0 ; 313 // Get the End-Position and End-Momentum ( << 302 fGeometryLimitedStep = false ; 314 // << 303 } 315 fTransportEndPosition = aFieldTrack.GetPos << 304 316 << 305 // Remember last safety origin & value. 317 // Momentum: Magnitude and direction can << 306 // 318 // << 307 fPreviousSftOrigin = startPosition ; 319 fMomentumChanged = true; << 308 fPreviousSafety = currentSafety ; 320 fTransportEndMomentumDir = aFieldTrack.Get << 309 // fpSafetyHelper->SetCurrentSafety( newSafety, startPosition); 321 << 310 322 fTransportEndKineticEnergy = aFieldTrack.G << 311 // Get the End-Position and End-Momentum (Dir-ection) 323 << 312 // 324 fCandidateEndGlobalTime = aFieldTrack.GetL << 313 fTransportEndPosition = aFieldTrack.GetPosition() ; 325 fEndGlobalTimeComputed = true; << 314 326 << 315 // Momentum: Magnitude and direction can be changed too now ... 327 fTransportEndSpin = aFieldTrack.GetSpin(); << 316 // 328 fParticleIsLooping = fFieldPropagator->IsP << 317 fMomentumChanged = true ; 329 endpointDistance = (fTransportEndPosition << 318 fTransportEndMomentumDir = aFieldTrack.GetMomentumDir() ; >> 319 >> 320 fTransportEndKineticEnergy = aFieldTrack.GetKineticEnergy() ; >> 321 >> 322 // if( fFieldPropagator->GetCurrentFieldManager()->DoesFieldChangeEnergy() ) >> 323 // { >> 324 // // If the field can change energy, then the time must be integrated >> 325 // // - so this should have been updated >> 326 >> 327 fCandidateEndGlobalTime = aFieldTrack.GetLabTimeOfFlight(); >> 328 fEndGlobalTimeComputed = true; >> 329 >> 330 // // was ( fCandidateEndGlobalTime != track.GetGlobalTime() ); >> 331 // // a cleaner way is to have FieldTrack knowing whether time is updated. >> 332 // } >> 333 // else >> 334 // { >> 335 // // The energy should be unchanged by field transport, >> 336 // // - so the time changed will be calculated elsewhere >> 337 // // >> 338 // fEndGlobalTimeComputed = false; >> 339 // >> 340 // // Check that the integration preserved the energy >> 341 // // - and if not correct this! >> 342 // G4double startEnergy= track.GetKineticEnergy(); >> 343 // G4double endEnergy= fTransportEndKineticEnergy; >> 344 // >> 345 // static G4int no_inexact_steps=0, no_large_ediff; >> 346 // G4double absEdiff = std::fabs(startEnergy- endEnergy); >> 347 // if( absEdiff > perMillion * endEnergy ) >> 348 // { >> 349 // no_inexact_steps++; >> 350 // // Possible statistics keeping here ... >> 351 // } >> 352 // if( fVerboseLevel > 1 ) >> 353 // { >> 354 // if( std::fabs(startEnergy- endEnergy) > perThousand * endEnergy ) >> 355 // { >> 356 // static G4int no_warnings= 0, warnModulo=1, moduloFactor= 10; >> 357 // no_large_ediff ++; >> 358 // if( (no_large_ediff% warnModulo) == 0 ) >> 359 // { >> 360 // no_warnings++; >> 361 // G4cout << "WARNING - G4MonopoleTransportation::AlongStepGetPIL() " >> 362 // << " Energy change in Step is above 1^-3 relative value. " << G4endl >> 363 // << " Relative change in 'tracking' step = " >> 364 // << std::setw(15) << (endEnergy-startEnergy)/startEnergy << G4endl >> 365 // << " Starting E= " << std::setw(12) << startEnergy / MeV << " MeV " << G4endl >> 366 // << " Ending E= " << std::setw(12) << endEnergy / MeV << " MeV " << G4endl; >> 367 // G4cout << " Energy has been corrected -- however, review" >> 368 // << " field propagation parameters for accuracy." << G4endl; >> 369 // if( (fVerboseLevel > 2 ) || (no_warnings<4) || (no_large_ediff == warnModulo * moduloFactor) ){ >> 370 // G4cout << " These include EpsilonStepMax(/Min) in G4FieldManager " >> 371 // << " which determine fractional error per step for integrated quantities. " << G4endl >> 372 // << " Note also the influence of the permitted number of integration steps." >> 373 // << G4endl; >> 374 // } >> 375 // G4cerr << "ERROR - G4MonopoleTransportation::AlongStepGetPIL()" << G4endl >> 376 // << " Bad 'endpoint'. Energy change detected" >> 377 // << " and corrected. " >> 378 // << " Has occurred already " >> 379 // << no_large_ediff << " times." << G4endl; >> 380 // if( no_large_ediff == warnModulo * moduloFactor ) >> 381 // { >> 382 // warnModulo *= moduloFactor; >> 383 // } >> 384 // } >> 385 // } >> 386 // } // end of if (fVerboseLevel) >> 387 // >> 388 // // Correct the energy for fields that conserve it >> 389 // // This - hides the integration error >> 390 // // - but gives a better physical answer >> 391 // fTransportEndKineticEnergy= track.GetKineticEnergy(); >> 392 // } >> 393 >> 394 >> 395 fTransportEndSpin = aFieldTrack.GetSpin(); >> 396 fParticleIsLooping = fFieldPropagator->IsParticleLooping() ; >> 397 endpointDistance = (fTransportEndPosition - startPosition).mag() ; 330 } 398 } 331 399 332 // If we are asked to go a step length of 0, 400 // If we are asked to go a step length of 0, and we are on a boundary 333 // then a boundary will also limit the step 401 // then a boundary will also limit the step -> we must flag this. 334 // 402 // 335 if (currentMinimumStep == 0.0) { << 403 if( currentMinimumStep == 0.0 ) 336 if (currentSafety == 0.0) fGeometryLimited << 404 { >> 405 if( currentSafety == 0.0 ) fGeometryLimitedStep = true ; 337 } 406 } 338 407 339 // Update the safety starting from the end-p 408 // Update the safety starting from the end-point, 340 // if it will become negative at the end-poi 409 // if it will become negative at the end-point. 341 // 410 // 342 if (currentSafety < endpointDistance) { << 411 if( currentSafety < endpointDistance ) 343 // if( particleMagneticCharge == 0.0 ) << 412 { 344 // G4cout << " Avoiding call to ComputeS << 413 // if( particleMagneticCharge == 0.0 ) 345 << 414 // G4cout << " Avoiding call to ComputeSafety : charge = 0.0 " << G4endl; 346 if (particleMagneticCharge != 0.0) { << 415 347 G4double endSafety = fLinearNavigator->C << 416 if( particleMagneticCharge != 0.0 ) { 348 currentSafety = endSafety; << 417 349 fPreviousSftOrigin = fTransportEndPositi << 418 G4double endSafety = 350 fPreviousSafety = currentSafety; << 419 fLinearNavigator->ComputeSafety( fTransportEndPosition) ; 351 fpSafetyHelper->SetCurrentSafety(current << 420 currentSafety = endSafety ; 352 << 421 fPreviousSftOrigin = fTransportEndPosition ; 353 // Because the Stepping Manager assumes << 422 fPreviousSafety = currentSafety ; 354 // add the StepLength << 423 fpSafetyHelper->SetCurrentSafety( currentSafety, fTransportEndPosition); 355 // << 424 356 currentSafety += endpointDistance; << 425 // Because the Stepping Manager assumes it is from the start point, 357 << 426 // add the StepLength 358 #ifdef G4DEBUG_TRANSPORT << 427 // 359 G4cout.precision(12); << 428 currentSafety += endpointDistance ; 360 G4cout << "***G4MonopoleTransportation:: << 429 361 G4cout << " Called Navigator->ComputeSa << 430 #ifdef G4DEBUG_TRANSPORT 362 << " and it returned safety= " << 431 G4cout.precision(12) ; 363 G4cout << " Adding endpoint distance " << 432 G4cout << "***G4MonopoleTransportation::AlongStepGPIL ** " << G4endl ; 364 << " to obtain pseudo-safety= " << 433 G4cout << " Called Navigator->ComputeSafety at " << fTransportEndPosition >> 434 << " and it returned safety= " << endSafety << G4endl ; >> 435 G4cout << " Adding endpoint distance " << endpointDistance >> 436 << " to obtain pseudo-safety= " << currentSafety << G4endl ; 365 #endif 437 #endif 366 } << 438 } 367 } << 439 } 368 440 369 fParticleChange.ProposeTrueStepLength(geomet << 441 fParticleChange.ProposeTrueStepLength(geometryStepLength) ; 370 442 371 fMagSetup->SetStepperAndChordFinder(0); << 443 magSetup->SetStepperAndChordFinder(0); 372 // change back to usual equation 444 // change back to usual equation 373 445 374 return geometryStepLength; << 446 return geometryStepLength ; 375 } 447 } 376 448 377 //....oooOO0OOooo........oooOO0OOooo........oo << 449 ////////////////////////////////////////////////////////////////////////// 378 // 450 // 379 // Initialize ParticleChange (by setting al 451 // Initialize ParticleChange (by setting all its members equal 380 // to correspond 452 // to corresponding members in G4Track) 381 453 382 G4VParticleChange* G4MonopoleTransportation::A << 454 G4VParticleChange* G4MonopoleTransportation::AlongStepDoIt( const G4Track& track, 383 << 455 const G4Step& stepData ) 384 { 456 { 385 ++noCalls; << 457 static G4int noCalls=0; >> 458 static const G4ParticleDefinition* fOpticalPhoton = >> 459 G4ParticleTable::GetParticleTable()->FindParticle("opticalphoton"); 386 460 387 if (fGeometryLimitedStep) { << 461 noCalls++; 388 stepData.GetPostStepPoint()->SetStepStatus << 389 } << 390 462 391 fParticleChange.Initialize(track); << 463 fParticleChange.Initialize(track) ; 392 464 393 // Code for specific process << 465 // Code for specific process 394 // 466 // 395 fParticleChange.ProposePosition(fTransportEn << 467 fParticleChange.ProposePosition(fTransportEndPosition) ; 396 fParticleChange.ProposeMomentumDirection(fTr << 468 fParticleChange.ProposeMomentumDirection(fTransportEndMomentumDir) ; 397 fParticleChange.ProposeEnergy(fTransportEndK << 469 fParticleChange.ProposeEnergy(fTransportEndKineticEnergy) ; 398 fParticleChange.SetMomentumChanged(fMomentum << 470 fParticleChange.SetMomentumChanged(fMomentumChanged) ; 399 471 400 fParticleChange.ProposePolarization(fTranspo 472 fParticleChange.ProposePolarization(fTransportEndSpin); 401 << 473 402 G4double deltaTime = 0.0; << 474 G4double deltaTime = 0.0 ; 403 475 404 // Calculate Lab Time of Flight (ONLY if fi 476 // Calculate Lab Time of Flight (ONLY if field Equations used it!) >> 477 // G4double endTime = fCandidateEndGlobalTime; >> 478 // G4double delta_time = endTime - startTime; 405 479 406 G4double startTime = track.GetGlobalTime(); << 480 G4double startTime = track.GetGlobalTime() ; 407 << 481 408 if (!fEndGlobalTimeComputed) { << 482 if (!fEndGlobalTimeComputed) 409 // The time was not integrated .. make the << 483 { 410 // << 484 // The time was not integrated .. make the best estimate possible 411 G4double finalVelocity = track.GetVelocity << 485 // 412 G4double initialVelocity = stepData.GetPre << 486 G4double finalVelocity = track.GetVelocity() ; 413 G4double stepLength = track.GetStepLength( << 487 G4double initialVelocity = stepData.GetPreStepPoint()->GetVelocity() ; 414 << 488 G4double stepLength = track.GetStepLength() ; 415 deltaTime = 0.0; // in case initialVeloci << 489 416 if (finalVelocity > 0.0) { << 490 deltaTime= 0.0; // in case initialVelocity = 0 417 G4double meanInverseVelocity; << 491 const G4DynamicParticle* fpDynamicParticle = track.GetDynamicParticle(); 418 // deltaTime = stepLength/finalVelocity << 492 if (fpDynamicParticle->GetDefinition()== fOpticalPhoton) 419 meanInverseVelocity = 0.5 * (1.0 / initi << 493 { 420 deltaTime = stepLength * meanInverseVelo << 494 // A photon is in the medium of the final point 421 } << 495 // during the step, so it has the final velocity. 422 else if (initialVelocity > 0.0) { << 496 deltaTime = stepLength/finalVelocity ; 423 deltaTime = stepLength / initialVelocity << 497 } 424 } << 498 else if (finalVelocity > 0.0) 425 fCandidateEndGlobalTime = startTime + delt << 499 { >> 500 G4double meanInverseVelocity ; >> 501 // deltaTime = stepLength/finalVelocity ; >> 502 meanInverseVelocity = 0.5 >> 503 * ( 1.0 / initialVelocity + 1.0 / finalVelocity ) ; >> 504 deltaTime = stepLength * meanInverseVelocity ; >> 505 } >> 506 else if( initialVelocity > 0.0 ) >> 507 { >> 508 deltaTime = stepLength/initialVelocity ; >> 509 } >> 510 fCandidateEndGlobalTime = startTime + deltaTime ; 426 } 511 } 427 else { << 512 else 428 deltaTime = fCandidateEndGlobalTime - star << 513 { >> 514 deltaTime = fCandidateEndGlobalTime - startTime ; 429 } 515 } 430 516 431 fParticleChange.ProposeGlobalTime(fCandidate << 517 fParticleChange.ProposeGlobalTime( fCandidateEndGlobalTime ) ; 432 518 433 // Now Correct by Lorentz factor to get "pro 519 // Now Correct by Lorentz factor to get "proper" deltaTime >> 520 >> 521 G4double restMass = track.GetDynamicParticle()->GetMass() ; >> 522 G4double deltaProperTime = deltaTime*( restMass/track.GetTotalEnergy() ) ; 434 523 435 G4double restMass = track.GetDynamicParticle << 524 fParticleChange.ProposeProperTime(track.GetProperTime() + deltaProperTime) ; 436 G4double deltaProperTime = deltaTime * (rest << 525 //fParticleChange. ProposeTrueStepLength( track.GetStepLength() ) ; 437 << 438 fParticleChange.ProposeProperTime(track.GetP << 439 526 440 // If the particle is caught looping or is s 527 // If the particle is caught looping or is stuck (in very difficult 441 // boundaries) in a magnetic field (doing ma << 528 // boundaries) in a magnetic field (doing many steps) 442 // THEN this kills it ... 529 // THEN this kills it ... 443 // 530 // 444 if (fParticleIsLooping) { << 531 if ( fParticleIsLooping ) 445 G4double endEnergy = fTransportEndKineticE << 532 { >> 533 G4double endEnergy= fTransportEndKineticEnergy; 446 534 447 if ((endEnergy < fThreshold_Important_Ener << 535 if( (endEnergy < fThreshold_Important_Energy) 448 // Kill the looping particle << 536 || (fNoLooperTrials >= fThresholdTrials ) ){ 449 // << 537 // Kill the looping particle 450 fParticleChange.ProposeTrackStatus(fStop << 538 // 451 << 539 fParticleChange.ProposeTrackStatus( fStopAndKill ) ; 452 // 'Bare' statistics << 540 453 fSumEnergyKilled += endEnergy; << 541 // 'Bare' statistics 454 if (endEnergy > fMaxEnergyKilled) { << 542 fSumEnergyKilled += endEnergy; 455 fMaxEnergyKilled = endEnergy; << 543 if( endEnergy > fMaxEnergyKilled) { fMaxEnergyKilled= endEnergy; } 456 } << 457 544 458 #ifdef G4VERBOSE 545 #ifdef G4VERBOSE 459 if ((verboseLevel > 1) || (endEnergy > f << 546 if( (fVerboseLevel > 1) || 460 G4cout << " G4MonopoleTransportation i << 547 ( endEnergy > fThreshold_Warning_Energy ) ) { 461 << "that is looping or stuck " << 548 G4cout << " G4MonopoleTransportation is killing track that is looping or stuck " 462 << track.GetKineticEnergy() / M << 549 << G4endl 463 G4cout << " Number of trials = " << << 550 << " This track has " << track.GetKineticEnergy() / MeV 464 << " No of calls to AlongStep << 551 << " MeV energy." << G4endl; 465 } << 552 G4cout << " Number of trials = " << fNoLooperTrials >> 553 << " No of calls to AlongStepDoIt = " << noCalls >> 554 << G4endl; >> 555 } 466 #endif 556 #endif 467 fNoLooperTrials = 0; << 557 fNoLooperTrials=0; 468 } << 469 else { << 470 fNoLooperTrials++; << 471 #ifdef G4VERBOSE << 472 if ((verboseLevel > 2)) { << 473 G4cout << " G4MonopoleTransportation << 474 << "Particle looping - " << 475 << " Number of trials = " << << 476 << G4endl; << 477 } 558 } >> 559 else{ >> 560 fNoLooperTrials ++; >> 561 #ifdef G4VERBOSE >> 562 if( (fVerboseLevel > 2) ){ >> 563 G4cout << " G4MonopoleTransportation::AlongStepDoIt(): Particle looping - " >> 564 << " Number of trials = " << fNoLooperTrials >> 565 << " No of calls to = " << noCalls >> 566 << G4endl; >> 567 } 478 #endif 568 #endif 479 } << 569 } 480 } << 570 }else{ 481 else { << 571 fNoLooperTrials=0; 482 fNoLooperTrials = 0; << 483 } 572 } 484 573 485 // Another (sometimes better way) is to use 574 // Another (sometimes better way) is to use a user-limit maximum Step size 486 // to alleviate this problem .. << 575 // to alleviate this problem .. 487 576 488 // Introduce smooth curved trajectories to p 577 // Introduce smooth curved trajectories to particle-change 489 // 578 // 490 fParticleChange.SetPointerToVectorOfAuxiliar << 579 fParticleChange.SetPointerToVectorOfAuxiliaryPoints 491 fFieldPropagator->GimmeTrajectoryVectorAnd << 580 (fFieldPropagator->GimmeTrajectoryVectorAndForgetIt() ); 492 581 493 return &fParticleChange; << 582 return &fParticleChange ; 494 } 583 } 495 584 496 //....oooOO0OOooo........oooOO0OOooo........oo << 585 ////////////////////////////////////////////////////////////////////////// 497 // 586 // 498 // This ensures that the PostStep action is a 587 // This ensures that the PostStep action is always called, 499 // so that it can do the relocation if it is 588 // so that it can do the relocation if it is needed. 500 // << 589 // 501 590 502 G4double << 591 G4double G4MonopoleTransportation:: 503 G4MonopoleTransportation::PostStepGetPhysicalI << 592 PostStepGetPhysicalInteractionLength( const G4Track&, 504 << 593 G4double, // previousStepSize 505 << 594 G4ForceCondition* pForceCond ) 506 { << 595 { 507 *pForceCond = Forced; << 596 *pForceCond = Forced ; 508 return DBL_MAX; // was kInfinity ; but conv << 597 return DBL_MAX ; // was kInfinity ; but convention now is DBL_MAX 509 } 598 } 510 599 511 //....oooOO0OOooo........oooOO0OOooo........oo << 600 ///////////////////////////////////////////////////////////////////////////// >> 601 // 512 602 513 G4VParticleChange* G4MonopoleTransportation::P << 603 G4VParticleChange* G4MonopoleTransportation::PostStepDoIt( const G4Track& track, >> 604 const G4Step& ) 514 { 605 { 515 G4TouchableHandle retCurrentTouchable; // T << 606 G4TouchableHandle retCurrentTouchable ; // The one to return 516 607 517 // Initialize ParticleChange (by setting al 608 // Initialize ParticleChange (by setting all its members equal 518 // to correspond 609 // to corresponding members in G4Track) 519 // fParticleChange.Initialize(track) ; // T 610 // fParticleChange.Initialize(track) ; // To initialise TouchableChange 520 611 521 fParticleChange.ProposeTrackStatus(track.Get << 612 fParticleChange.ProposeTrackStatus(track.GetTrackStatus()) ; 522 613 523 // If the Step was determined by the volume 614 // If the Step was determined by the volume boundary, 524 // logically relocate the particle 615 // logically relocate the particle 525 << 616 526 if (fGeometryLimitedStep) { << 617 if(fGeometryLimitedStep) 527 // fCurrentTouchable will now become the p << 618 { >> 619 // fCurrentTouchable will now become the previous touchable, 528 // and what was the previous will be freed 620 // and what was the previous will be freed. 529 // (Needed because the preStepPoint can po 621 // (Needed because the preStepPoint can point to the previous touchable) 530 622 531 fLinearNavigator->SetGeometricallyLimitedS << 623 fLinearNavigator->SetGeometricallyLimitedStep() ; 532 fLinearNavigator->LocateGlobalPointAndUpda << 624 fLinearNavigator-> 533 track.GetPosition(), track.GetMomentumDi << 625 LocateGlobalPointAndUpdateTouchableHandle( track.GetPosition(), 534 // Check whether the particle is out of th << 626 track.GetMomentumDirection(), >> 627 fCurrentTouchableHandle, >> 628 true ) ; >> 629 // Check whether the particle is out of the world volume 535 // If so it has exited and must be killed. 630 // If so it has exited and must be killed. 536 // 631 // 537 if (fCurrentTouchableHandle->GetVolume() = << 632 if( fCurrentTouchableHandle->GetVolume() == 0 ) 538 fParticleChange.ProposeTrackStatus(fStop << 633 { >> 634 fParticleChange.ProposeTrackStatus( fStopAndKill ) ; 539 } 635 } 540 retCurrentTouchable = fCurrentTouchableHan << 636 retCurrentTouchable = fCurrentTouchableHandle ; 541 fParticleChange.SetTouchableHandle(fCurren << 637 fParticleChange.SetTouchableHandle( fCurrentTouchableHandle ) ; 542 } 638 } 543 else // fGeometryLimitedStep is false << 639 else // fGeometryLimitedStep is false 544 { << 640 { 545 // This serves only to move the Navigator' 641 // This serves only to move the Navigator's location 546 // 642 // 547 fLinearNavigator->LocateGlobalPointWithinV << 643 fLinearNavigator->LocateGlobalPointWithinVolume( track.GetPosition() ) ; 548 644 549 // The value of the track's current Toucha << 645 // The value of the track's current Touchable is retained. 550 // (and it must be correct because we must 646 // (and it must be correct because we must use it below to 551 // overwrite the (unset) one in particle c 647 // overwrite the (unset) one in particle change) 552 // It must be fCurrentTouchable too ?? 648 // It must be fCurrentTouchable too ?? 553 // 649 // 554 fParticleChange.SetTouchableHandle(track.G << 650 fParticleChange.SetTouchableHandle( track.GetTouchableHandle() ) ; 555 retCurrentTouchable = track.GetTouchableHa << 651 retCurrentTouchable = track.GetTouchableHandle() ; 556 } // endif ( fGeometryLimitedStep ) << 652 } // endif ( fGeometryLimitedStep ) 557 << 653 558 const G4VPhysicalVolume* pNewVol = retCurren << 654 const G4VPhysicalVolume* pNewVol = retCurrentTouchable->GetVolume() ; 559 const G4Material* pNewMaterial = 0; << 655 const G4Material* pNewMaterial = 0 ; 560 G4VSensitiveDetector* pNewSensitiveDetector << 656 const G4VSensitiveDetector* pNewSensitiveDetector = 0 ; 561 if (pNewVol != 0) { << 657 562 pNewMaterial = pNewVol->GetLogicalVolume() << 658 if( pNewVol != 0 ) 563 pNewSensitiveDetector = pNewVol->GetLogica << 659 { >> 660 pNewMaterial= pNewVol->GetLogicalVolume()->GetMaterial(); >> 661 pNewSensitiveDetector= pNewVol->GetLogicalVolume()->GetSensitiveDetector(); 564 } 662 } 565 663 566 // ( <const_cast> pNewMaterial ) ; 664 // ( <const_cast> pNewMaterial ) ; >> 665 // ( <const_cast> pNewSensitiveDetector) ; 567 666 568 fParticleChange.SetMaterialInTouchable((G4Ma << 667 fParticleChange.SetMaterialInTouchable( (G4Material *) pNewMaterial ) ; 569 fParticleChange.SetSensitiveDetectorInToucha << 668 fParticleChange.SetSensitiveDetectorInTouchable( (G4VSensitiveDetector *) pNewSensitiveDetector ) ; 570 669 571 const G4MaterialCutsCouple* pNewMaterialCuts 670 const G4MaterialCutsCouple* pNewMaterialCutsCouple = 0; 572 if (pNewVol != 0) { << 671 if( pNewVol != 0 ) 573 pNewMaterialCutsCouple = pNewVol->GetLogic << 672 { >> 673 pNewMaterialCutsCouple=pNewVol->GetLogicalVolume()->GetMaterialCutsCouple(); 574 } 674 } 575 675 576 if (pNewVol != 0 && pNewMaterialCutsCouple ! << 676 if( pNewVol!=0 && pNewMaterialCutsCouple!=0 && pNewMaterialCutsCouple->GetMaterial()!=pNewMaterial ) 577 && pNewMaterialCutsCouple->GetMaterial() << 578 { 677 { 579 // for parametrized volume 678 // for parametrized volume 580 // 679 // 581 pNewMaterialCutsCouple = G4ProductionCutsT << 680 pNewMaterialCutsCouple = 582 pNewMaterial, pNewMaterialCutsCouple->Ge << 681 G4ProductionCutsTable::GetProductionCutsTable() >> 682 ->GetMaterialCutsCouple(pNewMaterial, >> 683 pNewMaterialCutsCouple->GetProductionCuts()); 583 } 684 } 584 fParticleChange.SetMaterialCutsCoupleInTouch << 685 fParticleChange.SetMaterialCutsCoupleInTouchable( pNewMaterialCutsCouple ); 585 686 586 // temporarily until Get/Set Material of Par << 687 // temporarily until Get/Set Material of ParticleChange, 587 // and StepPoint can be made const. << 688 // and StepPoint can be made const. 588 // Set the touchable in ParticleChange 689 // Set the touchable in ParticleChange 589 // this must always be done because the part 690 // this must always be done because the particle change always 590 // uses this value to overwrite the current 691 // uses this value to overwrite the current touchable pointer. 591 // 692 // 592 fParticleChange.SetTouchableHandle(retCurren << 693 fParticleChange.SetTouchableHandle(retCurrentTouchable) ; 593 << 694 594 return &fParticleChange; << 695 return &fParticleChange ; 595 } 696 } 596 697 597 //....oooOO0OOooo........oooOO0OOooo........oo << 698 // New method takes over the responsibility to reset the state of G4MonopoleTransportation 598 << 699 // object at the start of a new track or the resumption of a suspended track. 599 // New method takes over the responsibility to << 600 // of G4MonopoleTransportation object at the s << 601 // or the resumption of a suspended track. << 602 700 603 void G4MonopoleTransportation::StartTracking(G << 701 void >> 702 G4MonopoleTransportation::StartTracking(G4Track* aTrack) 604 { 703 { 605 G4VProcess::StartTracking(aTrack); 704 G4VProcess::StartTracking(aTrack); 606 705 607 // The actions here are those that were take << 706 // The actions here are those that were taken in AlongStepGPIL 608 // when track.GetCurrentStepNumber()==1 << 707 // when track.GetCurrentStepNumber()==1 609 708 610 // reset safety value and center 709 // reset safety value and center 611 // 710 // 612 fPreviousSafety = 0.0; << 711 fPreviousSafety = 0.0 ; 613 fPreviousSftOrigin = G4ThreeVector(0., 0., 0 << 712 fPreviousSftOrigin = G4ThreeVector(0.,0.,0.) ; 614 << 713 615 // reset looping counter -- for motion in fi 714 // reset looping counter -- for motion in field 616 fNoLooperTrials = 0; << 715 fNoLooperTrials= 0; 617 // Must clear this state .. else it depends 716 // Must clear this state .. else it depends on last track's value 618 // --> a better solution would set this fro << 717 // --> a better solution would set this from state of suspended track TODO ? 619 // Was if( aTrack->GetCurrentStepNumber()==1 718 // Was if( aTrack->GetCurrentStepNumber()==1 ) { .. } 620 719 621 // ChordFinder reset internal state 720 // ChordFinder reset internal state 622 // 721 // 623 if (DoesGlobalFieldExist()) { << 722 if( DoesGlobalFieldExist() ) { 624 fFieldPropagator->ClearPropagatorState(); << 723 fFieldPropagator->ClearPropagatorState(); 625 // Resets all state of field propagator cl << 724 // Resets all state of field propagator class (ONLY) 626 // including safety values << 725 // including safety values (in case of overlaps and to wipe for first track). 627 // in case of overlaps and to wipe for fir << 628 726 629 // G4ChordFinder* chordF= fFieldPropagator << 727 // G4ChordFinder* chordF= fFieldPropagator->GetChordFinder(); 630 // if( chordF ) chordF->ResetStepEstimate( << 728 // if( chordF ) chordF->ResetStepEstimate(); 631 } 729 } 632 730 633 // Make sure to clear the chord finders of a 731 // Make sure to clear the chord finders of all fields (ie managers) 634 G4FieldManagerStore* fieldMgrStore = G4Field << 732 static G4FieldManagerStore* fieldMgrStore= G4FieldManagerStore::GetInstance(); 635 fieldMgrStore->ClearAllChordFindersState(); << 733 fieldMgrStore->ClearAllChordFindersState(); 636 734 637 // Update the current touchable handle (fro 735 // Update the current touchable handle (from the track's) 638 // 736 // 639 fCurrentTouchableHandle = aTrack->GetTouchab 737 fCurrentTouchableHandle = aTrack->GetTouchableHandle(); 640 } 738 } 641 739 642 //....oooOO0OOooo........oooOO0OOooo........oo << 643 740