source: trunk/source/processes/transportation/src/G4Transportation.cc@ 880

Last change on this file since 880 was 819, checked in by garnier, 17 years ago

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