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

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

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