| Geant4 Cross Reference |
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 // $Id: G4Transportation.cc,v 1.60 2006/11/22 18:47:02 japost Exp $
28 // GEANT4 tag $Name: geant4-08-02 $
29 //
30 // ------------------------------------------------------------
31 // GEANT 4 include file implementation
32 //
33 // ------------------------------------------------------------
34 //
35 // This class is a process responsible for the transportation of
36 // a particle, ie the geometrical propagation that encounters the
37 // geometrical sub-volumes of the detectors.
38 //
39 // It is also tasked with part of updating the "safety".
40 //
41 // =======================================================================
42 // Modified:
43 // 19 Jan 2006, P.MoraDeFreitas: Fix for suspended tracks (StartTracking)
44 // 11 Aug 2004, M.Asai: Add G4VSensitiveDetector* for updating stepPoint.
45 // 21 June 2003, J.Apostolakis: Calling field manager with
46 // track, to enable it to configure its accuracy
47 // 13 May 2003, J.Apostolakis: Zero field areas now taken into
48 // account correclty in all cases (thanks to W Pokorski).
49 // 29 June 2001, J.Apostolakis, D.Cote-Ahern, P.Gumplinger:
50 // correction for spin tracking
51 // 20 Febr 2001, J.Apostolakis: update for new FieldTrack
52 // 22 Sept 2000, V.Grichine: update of Kinetic Energy
53 // 9 June 1999, J.Apostolakis & S.Giani: protect full relocation
54 // in DEBUG for track starting on surface that
55 // goes step < tolerance.
56 // Created: 19 March 1997, J. Apostolakis
57 // =======================================================================
58
59 #include "G4Transportation.hh"
60 #include "G4ProductionCutsTable.hh"
61 #include "G4ParticleTable.hh"
62 #include "G4ChordFinder.hh"
63 class G4VSensitiveDetector;
64
65 //////////////////////////////////////////////////////////////////////////
66 //
67 // Constructor
68
69 G4Transportation::G4Transportation( G4int verboseLevel )
70 : G4VProcess( G4String("Transportation"), fTransportation ),
71 fParticleIsLooping( false ),
72 fPreviousSftOrigin (0.,0.,0.),
73 fPreviousSafety ( 0.0 ),
74 fThreshold_Warning_Energy( 100 * MeV ),
75 fThreshold_Important_Energy( 250 * MeV ),
76 fThresholdTrials( 10 ),
77 fUnimportant_Energy( 1 * MeV ),
78 fNoLooperTrials(0),
79 fSumEnergyKilled( 0.0 ), fMaxEnergyKilled( 0.0 ),
80 fVerboseLevel( verboseLevel )
81 {
82 G4TransportationManager* transportMgr ;
83
84 transportMgr = G4TransportationManager::GetTransportationManager() ;
85
86 fLinearNavigator = transportMgr->GetNavigatorForTracking() ;
87
88 // fGlobalFieldMgr = transportMgr->GetFieldManager() ;
89
90 fFieldPropagator = transportMgr->GetPropagatorInField() ;
91
92 // Cannot determine whether a field exists here,
93 // because it would only work if the field manager has informed
94 // about the detector's field before this transportation process
95 // is constructed.
96 // Instead later the method DoesGlobalFieldExist() is called
97
98 // Small memory leak from this new
99 // fCurrentTouchableHandle = new G4TouchableHistory();
100 // Fixes:
101 // 1) static G4TouchableHistory sfNullTouchableHistory = new G4TouchableHistory();
102 // fCurrentTouchableHandle = new G4TouchableHandle( sfNullTouchableHistory );
103 // 2) set it to (G4TouchableHistory*) 0
104 // fCurrentTouchableHandle = new G4TouchableHandle( (G4TouchableHistory*) 0 );
105 // 3) Below:
106 static G4TouchableHandle nullTouchableHandle; // Points to (G4VTouchable*) 0
107 fCurrentTouchableHandle = nullTouchableHandle;
108
109 fEndGlobalTimeComputed = false;
110 fCandidateEndGlobalTime = 0;
111 }
112
113 //////////////////////////////////////////////////////////////////////////
114
115 G4Transportation::~G4Transportation()
116 {
117
118 // --- Alternative code to delete 'junk data' in touchable handle
119 // ** Corresponds to the case that the constructor has
120 // fCurrentTouchableHandle = new G4TouchableHistory();
121 // ** Likely incorrect - as at the end of the simulation
122 // the handle no longer holds the original 'null' history
123 // G4VTouchable* pTouchable= fCurrentTouchableHandle(); // Incorrect
124 // delete pTouchable; // Incorrect
125
126 if( (fVerboseLevel > 0) && (fSumEnergyKilled > 0.0 ) ){
127 G4cout << " G4Transportation: Statistics for looping particles " << G4endl;
128 G4cout << " Sum of energy of loopers killed: " << fSumEnergyKilled << G4endl;
129 G4cout << " Max energy of loopers killed: " << fMaxEnergyKilled << G4endl;
130 }
131 }
132
133 //////////////////////////////////////////////////////////////////////////
134 //
135 // Responsibilities:
136 // Find whether the geometry limits the Step, and to what length
137 // Calculate the new value of the safety and return it.
138 // Store the final time, position and momentum.
139
140 G4double G4Transportation::
141 AlongStepGetPhysicalInteractionLength( const G4Track& track,
142 G4double, // previousStepSize
143 G4double currentMinimumStep,
144 G4double& currentSafety,
145 G4GPILSelection* selection )
146 {
147 G4double geometryStepLength, newSafety ;
148 fParticleIsLooping = false ;
149
150 // Initial actions moved to StartTrack()
151 // --------------------------------------
152 // Note: in case another process changes touchable handle
153 // it will be necessary to add here (for all steps)
154 // fCurrentTouchableHandle = aTrack->GetTouchableHandle();
155
156 // GPILSelection is set to defaule value of CandidateForSelection
157 // It is a return value
158 //
159 *selection = CandidateForSelection ;
160
161 // Get initial Energy/Momentum of the track
162 //
163 const G4DynamicParticle* pParticle = track.GetDynamicParticle() ;
164 const G4ParticleDefinition* pParticleDef = pParticle->GetDefinition() ;
165 G4ThreeVector startMomentumDir = pParticle->GetMomentumDirection() ;
166 G4ThreeVector startPosition = track.GetPosition() ;
167
168 // G4double theTime = track.GetGlobalTime() ;
169
170 // The Step Point safety can be limited by other geometries and/or the
171 // assumptions of any process - it's not always the geometrical safety.
172 // We calculate the starting point's isotropic safety here.
173 //
174 G4ThreeVector OriginShift = startPosition - fPreviousSftOrigin ;
175 G4double MagSqShift = OriginShift.mag2() ;
176 if( MagSqShift >= sqr(fPreviousSafety) )
177 {
178 currentSafety = 0.0 ;
179 }
180 else
181 {
182 currentSafety = fPreviousSafety - std::sqrt(MagSqShift) ;
183 }
184
185 // Is the particle charged ?
186 //
187 G4double particleCharge = pParticle->GetCharge() ;
188
189 fGeometryLimitedStep = false ;
190 // fEndGlobalTimeComputed = false ;
191
192 // There is no need to locate the current volume. It is Done elsewhere:
193 // On track construction
194 // By the tracking, after all AlongStepDoIts, in "Relocation"
195
196 // Check whether the particle have an (EM) field force exerting upon it
197 //
198 G4FieldManager* fieldMgr=0;
199 G4bool fieldExertsForce = false ;
200 if( (particleCharge != 0.0) )
201 {
202 fieldMgr= fFieldPropagator->FindAndSetFieldManager( track.GetVolume() );
203 if (fieldMgr != 0) {
204 // If the field manager has no field, there is no field !
205 fieldExertsForce = (fieldMgr->GetDetectorField() != 0);
206 }
207 }
208
209 // G4cout << " G4Transport: field exerts force= " << fieldExertsForce
210 // << " fieldMgr= " << fieldMgr << G4endl;
211
212 // Choose the calculation of the transportation: Field or not
213 //
214 if( !fieldExertsForce )
215 {
216 G4double linearStepLength ;
217 if( currentMinimumStep <= currentSafety )
218 {
219 // The Step is guaranteed to be taken
220 //
221 geometryStepLength = currentMinimumStep ;
222 fGeometryLimitedStep = false ;
223 }
224 else
225 {
226 // Find whether the straight path intersects a volume
227 //
228 linearStepLength = fLinearNavigator->ComputeStep( startPosition,
229 startMomentumDir,
230 currentMinimumStep,
231 newSafety) ;
232 // Remember last safety origin & value.
233 //
234 fPreviousSftOrigin = startPosition ;
235 fPreviousSafety = newSafety ;
236
237 // The safety at the initial point has been re-calculated:
238 //
239 currentSafety = newSafety ;
240
241 if( linearStepLength <= currentMinimumStep)
242 {
243 // The geometry limits the Step size (an intersection was found.)
244 //
245 geometryStepLength = linearStepLength ;
246 fGeometryLimitedStep = true ;
247 }
248 else
249 {
250 // The full Step is taken.
251 //
252 geometryStepLength = currentMinimumStep ;
253 fGeometryLimitedStep = false ;
254 }
255 }
256 endpointDistance = geometryStepLength ;
257
258 // Calculate final position
259 //
260 fTransportEndPosition = startPosition+geometryStepLength*startMomentumDir ;
261
262 // Momentum direction, energy and polarisation are unchanged by transport
263 //
264 fTransportEndMomentumDir = startMomentumDir ;
265 fTransportEndKineticEnergy = track.GetKineticEnergy() ;
266 fTransportEndSpin = track.GetPolarization();
267 fParticleIsLooping = false ;
268 fMomentumChanged = false ;
269 fEndGlobalTimeComputed = false ;
270 }
271 else // A field exerts force
272 {
273 G4double momentumMagnitude = pParticle->GetTotalMomentum() ;
274 G4ThreeVector EndUnitMomentum ;
275 G4double lengthAlongCurve ;
276 G4double restMass = pParticleDef->GetPDGMass() ;
277
278 fFieldPropagator->SetChargeMomentumMass( particleCharge, // in e+ units
279 momentumMagnitude, // in Mev/c
280 restMass ) ;
281
282 // Message the field Manager, to configure it for this track
283 fieldMgr->ConfigureForTrack( &track );
284
285 G4ThreeVector spin = track.GetPolarization() ;
286 G4FieldTrack aFieldTrack = G4FieldTrack( startPosition,
287 track.GetMomentumDirection(),
288 0.0,
289 track.GetKineticEnergy(),
290 restMass,
291 track.GetVelocity(),
292 track.GetGlobalTime(), // Lab.
293 track.GetProperTime(), // Part.
294 &spin ) ;
295 if( currentMinimumStep > 0 )
296 {
297 // Do the Transport in the field (non recti-linear)
298 //
299 lengthAlongCurve = fFieldPropagator->ComputeStep( aFieldTrack,
300 currentMinimumStep,
301 currentSafety,
302 track.GetVolume() ) ;
303 if( lengthAlongCurve < currentMinimumStep)
304 {
305 geometryStepLength = lengthAlongCurve ;
306 fGeometryLimitedStep = true ;
307 }
308 else
309 {
310 geometryStepLength = currentMinimumStep ;
311 fGeometryLimitedStep = false ;
312 }
313 }
314 else
315 {
316 geometryStepLength = lengthAlongCurve= 0.0 ;
317 fGeometryLimitedStep = false ;
318 }
319
320 // Remember last safety origin & value.
321 //
322 fPreviousSftOrigin = startPosition ;
323 fPreviousSafety = currentSafety ;
324
325 // Get the End-Position and End-Momentum (Dir-ection)
326 //
327 fTransportEndPosition = aFieldTrack.GetPosition() ;
328
329 // Momentum: Magnitude and direction can be changed too now ...
330 //
331 fMomentumChanged = true ;
332 fTransportEndMomentumDir = aFieldTrack.GetMomentumDir() ;
333
334 fTransportEndKineticEnergy = aFieldTrack.GetKineticEnergy() ;
335
336 if( fFieldPropagator->GetCurrentFieldManager()->DoesFieldChangeEnergy() )
337 {
338 // If the field can change energy, then the time must be integrated
339 // - so this should have been updated
340 //
341 fCandidateEndGlobalTime = aFieldTrack.GetLabTimeOfFlight();
342 fEndGlobalTimeComputed = true;
343
344 // was ( fCandidateEndGlobalTime != track.GetGlobalTime() );
345 // a cleaner way is to have FieldTrack knowing whether time is updated.
346 }
347 else
348 {
349 // The energy should be unchanged by field transport,
350 // - so the time changed will be calculated elsewhere
351 //
352 fEndGlobalTimeComputed = false;
353
354 // Check that the integration preserved the energy
355 // - and if not correct this!
356 G4double startEnergy= track.GetKineticEnergy();
357 G4double endEnergy= fTransportEndKineticEnergy;
358
359 static G4int no_inexact_steps=0, no_large_ediff;
360 G4double absEdiff = std::fabs(startEnergy- endEnergy);
361 if( absEdiff > perMillion * endEnergy )
362 {
363 no_inexact_steps++;
364 // Possible statistics keeping here ...
365 }
366 if( fVerboseLevel > 1 )
367 {
368 if( std::fabs(startEnergy- endEnergy) > perThousand * endEnergy )
369 {
370 static G4int no_warnings= 0, warnModulo=1, moduloFactor= 10;
371 no_large_ediff ++;
372 if( (no_large_ediff% warnModulo) == 0 )
373 {
374 no_warnings++;
375 G4cout << "WARNING - G4Transportation::AlongStepGetPIL() "
376 << " Energy change in Step is above 1^-3 relative value. " << G4endl
377 << " Relative change in 'tracking' step = "
378 << std::setw(15) << (endEnergy-startEnergy)/startEnergy << G4endl
379 << " Starting E= " << std::setw(12) << startEnergy / MeV << " MeV " << G4endl
380 << " Ending E= " << std::setw(12) << endEnergy / MeV << " MeV " << G4endl;
381 G4cout << " Energy has been corrected -- however, review"
382 << " field propagation parameters for accuracy." << G4endl;
383 if( (fVerboseLevel > 2 ) || (no_warnings<4) || (no_large_ediff == warnModulo * moduloFactor) ){
384 G4cout << " These include EpsilonStepMax(/Min) in G4FieldManager "
385 << " which determine fractional error per step for integrated quantities. " << G4endl
386 << " Note also the influence of the permitted number of integration steps."
387 << G4endl;
388 }
389 G4cerr << "ERROR - G4Transportation::AlongStepGetPIL()" << G4endl
390 << " Bad 'endpoint'. Energy change detected"
391 << " and corrected. "
392 << " Has occurred already "
393 << no_large_ediff << " times." << G4endl;
394 if( no_large_ediff == warnModulo * moduloFactor )
395 {
396 warnModulo *= moduloFactor;
397 }
398 }
399 }
400 } // end of if (fVerboseLevel)
401
402 // Correct the energy for fields that conserve it
403 // This - hides the integration error
404 // - but gives a better physical answer
405 fTransportEndKineticEnergy= track.GetKineticEnergy();
406 }
407
408 fTransportEndSpin = aFieldTrack.GetSpin();
409 fParticleIsLooping = fFieldPropagator->IsParticleLooping() ;
410 endpointDistance = (fTransportEndPosition - startPosition).mag() ;
411 }
412
413 // If we are asked to go a step length of 0, and we are on a boundary
414 // then a boundary will also limit the step -> we must flag this.
415 //
416 if( currentMinimumStep == 0.0 )
417 {
418 if( currentSafety == 0.0 ) fGeometryLimitedStep = true ;
419 }
420
421 // Update the safety starting from the end-point,
422 // if it will become negative at the end-point.
423 //
424 if( currentSafety < endpointDistance )
425 {
426 G4double endSafety =
427 fLinearNavigator->ComputeSafety( fTransportEndPosition) ;
428 currentSafety = endSafety ;
429 fPreviousSftOrigin = fTransportEndPosition ;
430 fPreviousSafety = currentSafety ;
431
432 // Because the Stepping Manager assumes it is from the start point,
433 // add the StepLength
434 //
435 currentSafety += endpointDistance ;
436
437 #ifdef G4DEBUG_TRANSPORT
438 G4cout.precision(16) ;
439 G4cout << "***Transportation::AlongStepGPIL ** " << G4endl ;
440 G4cout << " Called Navigator->ComputeSafety at "
441 << fTransportEndPosition
442 << " and it returned safety= " << endSafety << G4endl ;
443 G4cout << " Adding endpoint distance " << endpointDistance
444 << " we obtain pseudo-safety= " << currentSafety << G4endl ;
445 #endif
446 }
447
448 fParticleChange.ProposeTrueStepLength(geometryStepLength) ;
449
450 return geometryStepLength ;
451 }
452
453 //////////////////////////////////////////////////////////////////////////
454 //
455 // Initialize ParticleChange (by setting all its members equal
456 // to corresponding members in G4Track)
457
458 G4VParticleChange* G4Transportation::AlongStepDoIt( const G4Track& track,
459 const G4Step& stepData )
460 {
461 static G4int noCalls=0;
462 static const G4ParticleDefinition* fOpticalPhoton =
463 G4ParticleTable::GetParticleTable()->FindParticle("opticalphoton");
464
465 noCalls++;
466
467 fParticleChange.Initialize(track) ;
468
469 // Code for specific process
470 //
471 fParticleChange.ProposePosition(fTransportEndPosition) ;
472 fParticleChange.ProposeMomentumDirection(fTransportEndMomentumDir) ;
473 fParticleChange.ProposeEnergy(fTransportEndKineticEnergy) ;
474 fParticleChange.SetMomentumChanged(fMomentumChanged) ;
475
476 fParticleChange.ProposePolarization(fTransportEndSpin);
477
478 G4double deltaTime = 0.0 ;
479
480 // Calculate Lab Time of Flight (ONLY if field Equations used it!)
481 // G4double endTime = fCandidateEndGlobalTime;
482 // G4double delta_time = endTime - startTime;
483
484 G4double startTime = track.GetGlobalTime() ;
485
486 if (!fEndGlobalTimeComputed)
487 {
488 // The time was not integrated .. make the best estimate possible
489 //
490 G4double finalVelocity = track.GetVelocity() ;
491 G4double initialVelocity = stepData.GetPreStepPoint()->GetVelocity() ;
492 G4double stepLength = track.GetStepLength() ;
493
494 const G4DynamicParticle* fpDynamicParticle = track.GetDynamicParticle();
495 if (fpDynamicParticle->GetDefinition()== fOpticalPhoton)
496 {
497 // A photon is in the medium of the final point
498 // during the step, so it has the final velocity.
499 deltaTime = stepLength/finalVelocity ;
500 }
501 else if (finalVelocity > 0.0)
502 {
503 G4double meanInverseVelocity ;
504 // deltaTime = stepLength/finalVelocity ;
505 meanInverseVelocity = 0.5
506 * ( 1.0 / initialVelocity + 1.0 / finalVelocity ) ;
507 deltaTime = stepLength * meanInverseVelocity ;
508 }
509 else
510 {
511 deltaTime = stepLength/initialVelocity ;
512 }
513 fCandidateEndGlobalTime = startTime + deltaTime ;
514 }
515 else
516 {
517 deltaTime = fCandidateEndGlobalTime - startTime ;
518 }
519
520 fParticleChange.ProposeGlobalTime( fCandidateEndGlobalTime ) ;
521
522 // Now Correct by Lorentz factor to get "proper" deltaTime
523
524 G4double restMass = track.GetDynamicParticle()->GetMass() ;
525 G4double deltaProperTime = deltaTime*( restMass/track.GetTotalEnergy() ) ;
526
527 fParticleChange.ProposeProperTime(track.GetProperTime() + deltaProperTime) ;
528 //fParticleChange. ProposeTrueStepLength( track.GetStepLength() ) ;
529
530 // If the particle is caught looping or is stuck (in very difficult
531 // boundaries) in a magnetic field (doing many steps)
532 // THEN this kills it ...
533 //
534 if ( fParticleIsLooping )
535 {
536 G4double endEnergy= fTransportEndKineticEnergy;
537
538 if( (endEnergy < fThreshold_Important_Energy)
539 || (fNoLooperTrials >= fThresholdTrials ) ){
540 // Kill the looping particle
541 //
542 fParticleChange.ProposeTrackStatus( fStopAndKill ) ;
543
544 // 'Bare' statistics
545 fSumEnergyKilled += endEnergy;
546 if( endEnergy > fMaxEnergyKilled) { fMaxEnergyKilled= endEnergy; }
547
548 #ifdef G4VERBOSE
549 if( (fVerboseLevel > 1) ||
550 ( endEnergy > fThreshold_Warning_Energy ) ) {
551 G4cout << " G4Transportation is killing track that is looping or stuck "
552 << G4endl
553 << " This track has " << track.GetKineticEnergy() / MeV
554 << " MeV energy." << G4endl;
555 G4cout << " Number of trials = " << fNoLooperTrials
556 << " No of calls to AlongStepDoIt = " << noCalls
557 << G4endl;
558 }
559 #endif
560 fNoLooperTrials=0;
561 }
562 else{
563 fNoLooperTrials ++;
564 #ifdef G4VERBOSE
565 if( (fVerboseLevel > 2) ){
566 G4cout << " Transportation::AlongStepDoIt(): Particle looping - "
567 << " Number of trials = " << fNoLooperTrials
568 << " No of calls to = " << noCalls
569 << G4endl;
570 }
571 #endif
572 }
573 }else{
574 fNoLooperTrials=0;
575 }
576
577 // Another (sometimes better way) is to use a user-limit maximum Step size
578 // to alleviate this problem ..
579
580 // Introduce smooth curved trajectories to particle-change
581 //
582 fParticleChange.SetPointerToVectorOfAuxiliaryPoints
583 (fFieldPropagator->GimmeTrajectoryVectorAndForgetIt() );
584
585 return &fParticleChange ;
586 }
587
588 //////////////////////////////////////////////////////////////////////////
589 //
590 // This ensures that the PostStep action is always called,
591 // so that it can do the relocation if it is needed.
592 //
593
594 G4double G4Transportation::
595 PostStepGetPhysicalInteractionLength( const G4Track&,
596 G4double, // previousStepSize
597 G4ForceCondition* pForceCond )
598 {
599 *pForceCond = Forced ;
600 return DBL_MAX ; // was kInfinity ; but convention now is DBL_MAX
601 }
602
603 /////////////////////////////////////////////////////////////////////////////
604 //
605
606 G4VParticleChange* G4Transportation::PostStepDoIt( const G4Track& track,
607 const G4Step& )
608 {
609 G4TouchableHandle retCurrentTouchable ; // The one to return
610
611 // Initialize ParticleChange (by setting all its members equal
612 // to corresponding members in G4Track)
613 // fParticleChange.Initialize(track) ; // To initialise TouchableChange
614
615 fParticleChange.ProposeTrackStatus(track.GetTrackStatus()) ;
616
617 // If the Step was determined by the volume boundary,
618 // logically relocate the particle
619
620 if(fGeometryLimitedStep)
621 {
622 // fCurrentTouchable will now become the previous touchable,
623 // and what was the previous will be freed.
624 // (Needed because the preStepPoint can point to the previous touchable)
625
626 fLinearNavigator->SetGeometricallyLimitedStep() ;
627 fLinearNavigator->
628 LocateGlobalPointAndUpdateTouchableHandle( track.GetPosition(),
629 track.GetMomentumDirection(),
630 fCurrentTouchableHandle,
631 true ) ;
632 // Check whether the particle is out of the world volume
633 // If so it has exited and must be killed.
634 //
635 if( fCurrentTouchableHandle->GetVolume() == 0 )
636 {
637 fParticleChange.ProposeTrackStatus( fStopAndKill ) ;
638 }
639 retCurrentTouchable = fCurrentTouchableHandle ;
640 fParticleChange.SetTouchableHandle( fCurrentTouchableHandle ) ;
641 }
642 else // fGeometryLimitedStep is false
643 {
644 // This serves only to move the Navigator's location
645 //
646 fLinearNavigator->LocateGlobalPointWithinVolume( track.GetPosition() ) ;
647
648 // The value of the track's current Touchable is retained.
649 // (and it must be correct because we must use it below to
650 // overwrite the (unset) one in particle change)
651 // It must be fCurrentTouchable too ??
652 //
653 fParticleChange.SetTouchableHandle( track.GetTouchableHandle() ) ;
654 retCurrentTouchable = track.GetTouchableHandle() ;
655 } // endif ( fGeometryLimitedStep )
656
657 const G4VPhysicalVolume* pNewVol = retCurrentTouchable->GetVolume() ;
658 const G4Material* pNewMaterial = 0 ;
659 const G4VSensitiveDetector* pNewSensitiveDetector = 0 ;
660
661 if( pNewVol != 0 )
662 {
663 pNewMaterial= pNewVol->GetLogicalVolume()->GetMaterial();
664 pNewSensitiveDetector= pNewVol->GetLogicalVolume()->GetSensitiveDetector();
665 }
666
667 // ( <const_cast> pNewMaterial ) ;
668 // ( <const_cast> pNewSensitiveDetector) ;
669
670 fParticleChange.SetMaterialInTouchable( (G4Material *) pNewMaterial ) ;
671 fParticleChange.SetSensitiveDetectorInTouchable( (G4VSensitiveDetector *) pNewSensitiveDetector ) ;
672
673 const G4MaterialCutsCouple* pNewMaterialCutsCouple = 0;
674 if( pNewVol != 0 )
675 {
676 pNewMaterialCutsCouple=pNewVol->GetLogicalVolume()->GetMaterialCutsCouple();
677 }
678
679 if( pNewVol!=0 && pNewMaterialCutsCouple!=0 && pNewMaterialCutsCouple->GetMaterial()!=pNewMaterial )
680 {
681 // for parametrized volume
682 //
683 pNewMaterialCutsCouple =
684 G4ProductionCutsTable::GetProductionCutsTable()
685 ->GetMaterialCutsCouple(pNewMaterial,
686 pNewMaterialCutsCouple->GetProductionCuts());
687 }
688 fParticleChange.SetMaterialCutsCoupleInTouchable( pNewMaterialCutsCouple );
689
690 // temporarily until Get/Set Material of ParticleChange,
691 // and StepPoint can be made const.
692 // Set the touchable in ParticleChange
693 // this must always be done because the particle change always
694 // uses this value to overwrite the current touchable pointer.
695 //
696 fParticleChange.SetTouchableHandle(retCurrentTouchable) ;
697
698 return &fParticleChange ;
699 }
700
701 // New method takes over the responsibility to reset the state of G4Transportation
702 // object at the start of a new track or the resumption of a suspended track.
703
704 void
705 G4Transportation::StartTracking(G4Track* aTrack)
706 {
707 G4VProcess::StartTracking(aTrack);
708
709 // The actions here are those that were taken in AlongStepGPIL
710 // when track.GetCurrentStepNumber()==1
711
712 // reset safety value and center
713 //
714 fPreviousSafety = 0.0 ;
715 fPreviousSftOrigin = G4ThreeVector(0.,0.,0.) ;
716
717 // reset looping counter -- for motion in field
718 if( aTrack->GetCurrentStepNumber()==1 ) {
719 fNoLooperTrials= 0;
720 }
721
722 // ChordFinder reset internal state
723 //
724 if( DoesGlobalFieldExist() ) {
725 G4ChordFinder* chordF= fFieldPropagator->GetChordFinder();
726 if( chordF ) chordF->ResetStepEstimate();
727 }
728
729 // Update the current touchable handle (from the track's)
730 //
731 fCurrentTouchableHandle = aTrack->GetTouchableHandle();
732 }
733
734