source: trunk/source/processes/transportation/src/G4CoupledTransportation.cc@ 1089

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

update to geant4.9.2
File size: 35.6 KB
Line 
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: G4CoupledTransportation.cc,v 1.27 2008/11/26 13:01:28 japost Exp $
28// --> Merged with 1.60.4.2.2.3 2007/05/09 09:30:28 japost
29// GEANT4 tag $Name: geant4-09-02 $
30// ------------------------------------------------------------
31// GEANT 4 class implementation
32// =======================================================================
33// Modified:
34// 13 May 2006, J. Apostolakis: Revised for parallel navigation (PathFinder)
35// 19 Jan 2006, P.MoraDeFreitas: Fix for suspended tracks (StartTracking)
36// 11 Aug 2004, M.Asai: Add G4VSensitiveDetector* for updating stepPoint.
37// 21 June 2003, J.Apostolakis: Calling field manager with
38// track, to enable it to configure its accuracy
39// 13 May 2003, J.Apostolakis: Zero field areas now taken into
40// account correclty in all cases (thanks to W Pokorski).
41// 29 June 2001, J.Apostolakis, D.Cote-Ahern, P.Gumplinger:
42// correction for spin tracking
43// 20 Febr 2001, J.Apostolakis: update for new FieldTrack
44// 22 Sept 2000, V.Grichine: update of Kinetic Energy
45// Created: 19 March 1997, J. Apostolakis
46// =======================================================================
47
48#include "G4CoupledTransportation.hh"
49#include "G4ProductionCutsTable.hh"
50#include "G4ParticleTable.hh"
51#include "G4ChordFinder.hh"
52#include "G4FieldManagerStore.hh"
53class G4VSensitiveDetector;
54
55//////////////////////////////////////////////////////////////////////////
56//
57// Constructor
58
59G4CoupledTransportation::G4CoupledTransportation( G4int verboseLevel )
60 : G4VProcess( G4String("CoupledTransportation"), fTransportation ),
61 fParticleIsLooping( false ),
62 fPreviousSftOrigin (0.,0.,0.),
63 fPreviousMassSafety ( 0.0 ),
64 fPreviousFullSafety ( 0.0 ),
65 fThreshold_Warning_Energy( 100 * MeV ),
66 fThreshold_Important_Energy( 250 * MeV ),
67 fThresholdTrials( 10 ),
68 fUnimportant_Energy( 1 * MeV ),
69 fNoLooperTrials(0),
70 fSumEnergyKilled( 0.0 ), fMaxEnergyKilled( 0.0 ),
71 fVerboseLevel(verboseLevel) // ( verboseLevel ? verboseLevel : 2 ) // or 4
72{
73 G4TransportationManager* transportMgr ;
74
75 transportMgr = G4TransportationManager::GetTransportationManager() ;
76
77 fMassNavigator = transportMgr->GetNavigatorForTracking() ;
78 fFieldPropagator = transportMgr->GetPropagatorInField() ;
79 // fGlobalFieldMgr = transportMgr->GetFieldManager() ;
80 fNavigatorId= transportMgr->ActivateNavigator( fMassNavigator );
81 if( fVerboseLevel > 0 ){
82 G4cout << " G4CoupledTransportation constructor: ----- " << G4endl;
83 G4cout << " Verbose level is " << fVerboseLevel << G4endl;
84 G4cout << " Navigator Id obtained in G4CoupledTransportation constructor "
85 << fNavigatorId << G4endl;
86 }
87 fPathFinder= G4PathFinder::GetInstance();
88 fpSafetyHelper = transportMgr->GetSafetyHelper(); // New
89
90 // Following assignment is to fix small memory leak from simple use of 'new'
91 static G4TouchableHandle nullTouchableHandle; // Points to (G4VTouchable*) 0
92 fCurrentTouchableHandle = nullTouchableHandle;
93 // fCurrentTouchableHandle = G4TouchableHandle( 0 ); // new G4TouchableHistory();
94
95 fEndGlobalTimeComputed = false;
96 fCandidateEndGlobalTime = 0;
97}
98
99//////////////////////////////////////////////////////////////////////////
100
101G4CoupledTransportation::~G4CoupledTransportation()
102{
103 // fCurrentTouchableHandle is a data member - no deletion required
104
105 if( (fVerboseLevel > 0) || (fSumEnergyKilled > 0.0 ) ){
106 G4cout << " G4CoupledTransportation: Statistics for looping particles " << G4endl;
107 G4cout << " Sum of energy of loopers killed: " << fSumEnergyKilled << G4endl;
108 G4cout << " Max energy of loopers killed: " << fMaxEnergyKilled << G4endl;
109 }
110}
111
112//////////////////////////////////////////////////////////////////////////
113//
114// Responsibilities:
115// Find whether the geometry limits the Step, and to what length
116// Calculate the new value of the safety and return it.
117// Store the final time, position and momentum.
118
119G4double G4CoupledTransportation::
120AlongStepGetPhysicalInteractionLength( const G4Track& track,
121 G4double, // previousStepSize
122 G4double currentMinimumStep,
123 G4double& proposedSafetyForStart,
124 G4GPILSelection* selection )
125{
126 G4double geometryStepLength;
127 G4double startMassSafety= 0.0; // estimated safety for start point (mass geometry)
128 G4double startFullSafety= 0.0; // estimated safety for start point (all geometries)
129 G4double safetyProposal= -1.0; // local copy of proposal
130
131 G4ThreeVector EndUnitMomentum ;
132 G4double lengthAlongCurve=0.0 ;
133
134 fParticleIsLooping = false ;
135
136 // Initial actions moved to StartTrack()
137 // --------------------------------------
138 // Note: in case another process changes touchable handle
139 // it will be necessary to add here (for all steps)
140 // fCurrentTouchableHandle = aTrack->GetTouchableHandle();
141
142 // GPILSelection is set to defaule value of CandidateForSelection
143 // It is a return value
144 //
145 *selection = CandidateForSelection ;
146
147 // Get initial Energy/Momentum of the track
148 //
149 const G4DynamicParticle* pParticle = track.GetDynamicParticle() ;
150 const G4ParticleDefinition* pParticleDef = pParticle->GetDefinition() ;
151 G4ThreeVector startMomentumDir = pParticle->GetMomentumDirection() ;
152 G4ThreeVector startPosition = track.GetPosition() ;
153 G4VPhysicalVolume* currentVolume= track.GetVolume();
154
155#ifdef G4DEBUG_TRANSPORT
156 if( fVerboseLevel > 1 ) {
157 G4cout << "G4CoupledTransportation::AlongStepGPIL> called in volume "
158 << currentVolume->GetName() << G4endl;
159 }
160#endif
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 startMassSafety = 0.0;
170 startFullSafety= 0.0;
171
172 // Recall that FullSafety <= MassSafety
173 // Original: if( MagSqShift < sqr(fPreviousMassSafety) ) {
174 if( MagSqShift < sqr(fPreviousFullSafety) ) { // Revision proposed by Alex H, 2 Oct 07
175 G4double mag_shift= std::sqrt(MagSqShift);
176 startMassSafety = std::max( (fPreviousMassSafety - mag_shift), 0.0);
177 startFullSafety = std::max( (fPreviousFullSafety - mag_shift), 0.0);
178 // Need to be consistent between full safety with Mass safety
179 // in order reproduce results in simple case --> use same calculation method
180
181 // Only compute full safety if massSafety > 0. Else it remains 0
182 // startFullSafety = fPathFinder->ComputeSafety( startPosition );
183 }
184
185 // Is the particle charged ?
186 //
187 G4double particleCharge = pParticle->GetCharge() ;
188
189 fMassGeometryLimitedStep = false ; // Set default - alex
190 fAnyGeometryLimitedStep = false;
191
192 // fEndGlobalTimeComputed = false ;
193
194 // There is no need to locate the current volume. It is Done elsewhere:
195 // On track construction
196 // By the tracking, after all AlongStepDoIts, in "Relocation"
197 // Check whether the particle has an (EM) field force exerted upon it
198 //
199 G4FieldManager* fieldMgr=0;
200 G4bool fieldExertsForce = false;
201 if( (particleCharge != 0.0) ) // || (magneticMoment != 0.0 ) )
202 {
203 fieldMgr= fFieldPropagator->FindAndSetFieldManager( currentVolume );
204 if (fieldMgr != 0) {
205 // Message the field Manager, to configure it for this track
206 fieldMgr->ConfigureForTrack( &track );
207 fieldExertsForce = (fieldMgr->GetDetectorField() != 0);
208 }
209 // the PathFinder will recognise whether the field exerts force
210 }
211 G4double momentumMagnitude = pParticle->GetTotalMomentum() ;
212 G4double restMass = pParticleDef->GetPDGMass() ;
213
214 fFieldPropagator->SetChargeMomentumMass( particleCharge, // in e+ units
215 momentumMagnitude, // in Mev/c
216 restMass ) ;
217 // This should be obsolete - the call to SetChargeAndMoments below should do the work
218
219 G4ThreeVector spin = track.GetPolarization() ;
220 G4FieldTrack theFieldTrack = G4FieldTrack( startPosition,
221 track.GetMomentumDirection(),
222 0.0,
223 track.GetKineticEnergy(),
224 restMass,
225 0.0, // UNUSED: track.GetVelocity(),
226 track.GetGlobalTime(), // Lab.
227 track.GetProperTime(), // Part.
228 &spin ) ;
229 theFieldTrack.SetChargeAndMoments( particleCharge ); // EM moments -- future extension
230
231 G4int stepNo= track.GetCurrentStepNumber();
232
233 ELimited limitedStep;
234 G4FieldTrack endTrackState('a'); // Default values
235
236 fMassGeometryLimitedStep = false ; // default
237 fAnyGeometryLimitedStep = false ;
238 if( currentMinimumStep > 0 ) {
239 G4double newMassSafety= 0.0; // temp. for recalculation
240
241 // Do the Transport in the field (non recti-linear)
242 //
243 lengthAlongCurve = fPathFinder->ComputeStep( theFieldTrack,
244 currentMinimumStep,
245 fNavigatorId,
246 stepNo,
247 newMassSafety,
248 limitedStep,
249 endTrackState,
250 currentVolume ) ;
251 // G4cout << " PathFinder ComputeStep returns " << lengthAlongCurve << G4endl;
252
253 G4double newFullSafety= fPathFinder->GetCurrentSafety();
254 // this was estimated already in step above
255 // G4double newFullStep= fPathFinder->GetMinimumStep();
256
257 if( limitedStep == kUnique || limitedStep == kSharedTransport ) {
258 fMassGeometryLimitedStep = true ;
259 }
260
261 fAnyGeometryLimitedStep = (fPathFinder->GetNumberGeometriesLimitingStep() != 0);
262
263 // #ifdef G4DEBUG
264 if( fMassGeometryLimitedStep && !fAnyGeometryLimitedStep ){
265 G4cerr << " Error in determining geometries limiting the step" << G4endl;
266 G4cerr << " Limiting: mass=" << fMassGeometryLimitedStep
267 << " any= " << fAnyGeometryLimitedStep << G4endl;
268 G4Exception("G4CoupledTransportation::AlongStepGetPhysicalInteractionLength()",
269 "PathFinderConfused",
270 FatalException,
271 "Incompatible conditions - was limited by a geometry?");
272 }
273 // #endif
274
275 // Other potential
276 // fAnyGeometryLimitedStep = newFullStep < currentMinimumStep;
277 // ^^^ Not good enough;
278 // Must compare with maximum requested step size
279 // (eg in case another process requested bigger, got this!)
280
281 geometryStepLength = std::min( lengthAlongCurve, currentMinimumStep);
282
283 // Momentum: Magnitude and direction can be changed too now ...
284 //
285 fMomentumChanged = true ;
286 fTransportEndMomentumDir = endTrackState.GetMomentumDir() ;
287
288 // Remember last safety origin & value.
289 fPreviousSftOrigin = startPosition ;
290 fPreviousMassSafety = newMassSafety ;
291 fPreviousFullSafety = newFullSafety ;
292 // fpSafetyHelper->SetCurrentSafety( newFullSafety, startPosition);
293
294#ifdef G4DEBUG_TRANSPORT
295 if( fVerboseLevel > 1 ){
296 G4cout << "G4Transport:CompStep> "
297 << " called the pathfinder for a new step at " << startPosition
298 << " and obtained step = " << lengthAlongCurve << G4endl;
299 G4cout << " New safety (preStep) = " << newMassSafety
300 << " versus precalculated = " << startMassSafety << G4endl;
301 }
302#endif
303
304 // Store as best estimate value
305 startMassSafety = newMassSafety ;
306 startFullSafety = newFullSafety ;
307
308 // Get the End-Position and End-Momentum (Dir-ection)
309 fTransportEndPosition = endTrackState.GetPosition() ;
310 fTransportEndKineticEnergy = endTrackState.GetKineticEnergy() ;
311 } else {
312 geometryStepLength = lengthAlongCurve= 0.0 ;
313 fMomentumChanged = false ;
314 // fMassGeometryLimitedStep = false ; // --- ???
315 // fAnyGeometryLimitedStep = true;
316 fTransportEndMomentumDir = track.GetMomentumDirection();
317 fTransportEndKineticEnergy = track.GetKineticEnergy();
318
319 fTransportEndPosition = startPosition;
320 // If the step length requested is 0, and we are on a boundary
321 // then a boundary will also limit the step.
322 if( startMassSafety == 0.0 ) {
323 fMassGeometryLimitedStep = true ;
324 fAnyGeometryLimitedStep = true;
325 }
326 // TODO: Add explicit logical status for being at a boundary
327 }
328 // G4FieldTrack aTrackState(endTrackState);
329
330 if( !fieldExertsForce )
331 {
332 fParticleIsLooping = false ;
333 fMomentumChanged = false ;
334 fEndGlobalTimeComputed = false ;
335 // G4cout << " global time is false " << G4endl;
336 }
337 else
338 {
339
340#ifdef G4DEBUG_TRANSPORT
341 if( fVerboseLevel > 1 ){
342 G4cout << " G4CT::CS End Position = " << fTransportEndPosition << G4endl;
343 G4cout << " G4CT::CS End Direction = " << fTransportEndMomentumDir << G4endl;
344 }
345#endif
346 // G4cout << " G4CoupledTransportation Before if change energy statement: " << fFieldPropagator->GetCurrentFieldManager()->DoesFieldChangeEnergy() << G4endl;
347 if( fFieldPropagator->GetCurrentFieldManager()->DoesFieldChangeEnergy() )
348 {
349 // If the field can change energy, then the time must be integrated
350 // - so this should have been updated
351 //
352 fCandidateEndGlobalTime = endTrackState.GetLabTimeOfFlight();
353 fEndGlobalTimeComputed = true;
354 // G4cout << " setting global time to true - why? " << G4endl;
355
356 // was ( fCandidateEndGlobalTime != track.GetGlobalTime() );
357 // a cleaner way is to have FieldTrack knowing whether time is updated.
358 }
359 else
360 {
361 // The energy should be unchanged by field transport,
362 // - so the time changed will be calculated elsewhere
363 //
364 fEndGlobalTimeComputed = false;
365
366 // Check that the integration preserved the energy
367 // - and if not correct this!
368 G4double startEnergy= track.GetKineticEnergy();
369 G4double endEnergy= fTransportEndKineticEnergy;
370
371 static G4int no_inexact_steps=0; // , no_large_ediff;
372 G4double absEdiff = std::fabs(startEnergy- endEnergy);
373 if( absEdiff > perMillion * endEnergy ) {
374 no_inexact_steps++;
375 // Possible statistics keeping here ...
376 }
377#ifdef G4VERBOSE
378 if( (fVerboseLevel > 1) && ( absEdiff > perThousand * endEnergy) ){
379 ReportInexactEnergy(startEnergy, endEnergy);
380 } // end of if (fVerboseLevel)
381#endif
382
383 // Correct the energy for fields that conserve it
384 // This - hides the integration error
385 // - but gives a better physical answer
386 fTransportEndKineticEnergy= track.GetKineticEnergy();
387 }
388 }
389
390 endpointDistance = (fTransportEndPosition - startPosition).mag() ;
391 fParticleIsLooping = fFieldPropagator->IsParticleLooping() ;
392
393 fTransportEndSpin = endTrackState.GetSpin();
394
395 // Calculate the safety
396 safetyProposal= startFullSafety; // used to be startMassSafety
397 // Changed to accomodate processes that cannot update the safety -- JA 22 Nov 06
398
399 // Update safety for the end-point, if becomes negative at the end-point.
400
401 if( (startFullSafety < endpointDistance )
402 && ( particleCharge != 0.0 ) ) // Only needed to prepare for Mult Scat.
403 // && !fAnyGeometryLimitedStep ) // To-Try: No safety update if at a boundary
404 {
405 G4double endFullSafety =
406 fPathFinder->ComputeSafety( fTransportEndPosition);
407 // Expected mission -- only mass geometry's safety
408 // fMassNavigator->ComputeSafety( fTransportEndPosition) ;
409 // Yet discrete processes only have poststep -- and this cannot
410 // currently revise the safety
411 // ==> so we use the all-geometry safety as a precaution
412
413 fpSafetyHelper->SetCurrentSafety( endFullSafety, fTransportEndPosition);
414 // Pushing safety to Helper avoids recalculation at this point
415
416 G4ThreeVector centerPt= G4ThreeVector(0.0, 0.0, 0.0); // Used for return value
417 G4double endMassSafety= fPathFinder->ObtainSafety( fNavigatorId, centerPt);
418 // Retrieves the mass value from PathFinder (it calculated it)
419
420 fPreviousMassSafety = endMassSafety ;
421 fPreviousFullSafety = endFullSafety;
422 fPreviousSftOrigin = fTransportEndPosition ;
423
424 // The convention (Stepping Manager's) is safety from the start point
425 //
426 safetyProposal = endFullSafety + endpointDistance;
427 // --> was endMassSafety
428 // Changed to accomodate processes that cannot update the safety -- JA 22 Nov 06
429
430 // #define G4DEBUG_TRANSPORT 1
431
432#ifdef G4DEBUG_TRANSPORT
433 int prec= G4cout.precision(12) ;
434 G4cout << "***Transportation::AlongStepGPIL ** " << G4endl ;
435 G4cout << " Revised Safety at endpoint " << fTransportEndPosition
436 << " give safety values: Mass= " << endMassSafety
437 << " All= " << endFullSafety << G4endl ;
438 G4cout << " Adding endpoint distance " << endpointDistance
439 << " to obtain pseudo-safety= " << safetyProposal << G4endl ;
440 G4cout.precision(prec);
441 }
442 else{
443 int prec= G4cout.precision(12) ;
444 G4cout << "***Transportation::AlongStepGPIL ** " << G4endl ;
445 G4cout << " Quick Safety estimate at endpoint " << fTransportEndPosition
446 << " gives safety endpoint value = " << startFullSafety - endpointDistance
447 << " using start-point value " << startFullSafety
448 << " and endpointDistance " << endpointDistance << G4endl;
449 G4cout.precision(prec);
450#endif
451 }
452
453 proposedSafetyForStart= safetyProposal;
454 fParticleChange.ProposeTrueStepLength(geometryStepLength) ;
455
456 return geometryStepLength ;
457}
458
459//////////////////////////////////////////////////////////////////////////
460
461G4VParticleChange* G4CoupledTransportation::AlongStepDoIt( const G4Track& track,
462 const G4Step& stepData )
463{
464 static G4int noCalls=0;
465 static const G4ParticleDefinition* fOpticalPhoton =
466 G4ParticleTable::GetParticleTable()->FindParticle("opticalphoton");
467
468 noCalls++;
469
470 fParticleChange.Initialize(track) ;
471 // sets all its members to the value of corresponding members in G4Track
472
473 // Code specific for Transport
474 //
475 fParticleChange.ProposePosition(fTransportEndPosition) ;
476 // G4cout << " G4CoupledTransportation::AlongStepDoIt"
477 // << " proposes position = " << fTransportEndPosition
478 // << " and end momentum direction = " << fTransportEndMomentumDir << G4endl;
479 fParticleChange.ProposeMomentumDirection(fTransportEndMomentumDir) ;
480 fParticleChange.ProposeEnergy(fTransportEndKineticEnergy) ;
481 fParticleChange.SetMomentumChanged(fMomentumChanged) ;
482
483 fParticleChange.ProposePolarization(fTransportEndSpin);
484
485 G4double deltaTime = 0.0 ;
486
487 // Calculate Lab Time of Flight (ONLY if field Equations used it!)
488 // G4double endTime = fCandidateEndGlobalTime;
489 // G4double delta_time = endTime - startTime;
490
491 G4double startTime = track.GetGlobalTime() ;
492
493 if (!fEndGlobalTimeComputed)
494 {
495 G4double finalInverseVel= DBL_MAX, initialInverseVel=DBL_MAX;
496
497 // The time was not integrated .. make the best estimate possible
498 //
499 G4double finalVelocity = track.GetVelocity() ;
500 if( finalVelocity > 0.0 ) { finalInverseVel= 1.0 / finalVelocity; }
501 G4double initialVelocity = stepData.GetPreStepPoint()->GetVelocity() ;
502 if( initialVelocity > 0.0 ) { initialInverseVel= 1.0 / initialVelocity; }
503 G4double stepLength = track.GetStepLength() ;
504
505 const G4DynamicParticle* fpDynamicParticle = track.GetDynamicParticle();
506 if (fpDynamicParticle->GetDefinition()== fOpticalPhoton)
507 {
508 // A photon is in the medium of the final point
509 // during the step, so Peter says it has the final velocity.
510 deltaTime = stepLength * finalInverseVel ;
511 // G4cout << " dt = s / finalV " << " s = " << stepLength
512 // << " 1 / finalV= " << finalInverseVel << G4endl;
513 }
514 else if (finalVelocity > 0.0)
515 {
516 // deltaTime = stepLength/finalVelocity ;
517 G4double meanInverseVelocity = 0.5 * ( initialInverseVel + finalInverseVel );
518 deltaTime = stepLength * meanInverseVelocity ;
519 // G4cout << " dt = s * mean(1/v) , with " << " s = " << stepLength
520 // << " mean(1/v)= " << meanInverseVelocity << G4endl;
521 }
522 else
523 {
524 deltaTime = stepLength * initialInverseVel ;
525 // G4cout << " dt = s / initV " << " s = " << stepLength
526 // << " 1 / initV= " << initialInverseVel << G4endl;
527 } // Could do with better estimate for final step (finalVelocity = 0) ?
528
529 fCandidateEndGlobalTime = startTime + deltaTime ;
530
531 // G4cout << " Calculated global time from start = " << startTime << " and "
532 // << " delta time = " << deltaTime << G4endl;
533 }
534 else
535 {
536 deltaTime = fCandidateEndGlobalTime - startTime ;
537 // G4cout << " Calculated global time from candidate end time = "
538 // << fCandidateEndGlobalTime << " and start time = " << startTime << G4endl;
539 }
540
541 // G4cout << " G4CoupledTransportation::AlongStepDoIt "
542 // << " flag whether computed time = " << fEndGlobalTimeComputed << " and "
543 // << " is proposes end time " << fCandidateEndGlobalTime << G4endl;
544 fParticleChange.ProposeGlobalTime( fCandidateEndGlobalTime ) ;
545
546 // Now Correct by Lorentz factor to get "proper" deltaTime
547
548 G4double restMass = track.GetDynamicParticle()->GetMass() ;
549 G4double deltaProperTime = deltaTime*( restMass/track.GetTotalEnergy() ) ;
550
551 fParticleChange.ProposeProperTime(track.GetProperTime() + deltaProperTime) ;
552 //fParticleChange. ProposeTrueStepLength( track.GetStepLength() ) ;
553
554 // If the particle is caught looping or is stuck (in very difficult
555 // boundaries) in a magnetic field (doing many steps)
556 // THEN this kills it ...
557 //
558 if ( fParticleIsLooping )
559 {
560 G4double endEnergy= fTransportEndKineticEnergy;
561
562 if( (endEnergy < fThreshold_Important_Energy)
563 || (fNoLooperTrials >= fThresholdTrials ) ){
564 // Kill the looping particle
565 //
566 fParticleChange.ProposeTrackStatus( fStopAndKill ) ;
567
568 // 'Bare' statistics
569 fSumEnergyKilled += endEnergy;
570 if( endEnergy > fMaxEnergyKilled) { fMaxEnergyKilled= endEnergy; }
571
572#ifdef G4VERBOSE
573 if( (fVerboseLevel > 1) ||
574 ( endEnergy > fThreshold_Warning_Energy ) ) {
575 G4cout << " G4CoupledTransportation is killing track that is looping or stuck " << G4endl
576 << " This track has " << track.GetKineticEnergy() / MeV
577 << " MeV energy." << G4endl;
578 }
579 if( fVerboseLevel > 0 ) {
580 G4cout << " Steps by this track: " << track.GetCurrentStepNumber() << G4endl;
581 }
582#endif
583 fNoLooperTrials=0;
584 } else{
585 fNoLooperTrials ++;
586#ifdef G4VERBOSE
587 if( (fVerboseLevel > 2) ){
588 G4cout << " ** G4CoupledTransportation::AlongStepDoIt(): Particle looping - " << G4endl
589 << " Number of consecutive problem step (this track) = " << fNoLooperTrials << G4endl
590 << " Steps by this track: " << track.GetCurrentStepNumber() << G4endl
591 << " Total no of calls to this method (all tracks) = " << noCalls << G4endl;
592 }
593#endif
594 }
595 }else{
596 fNoLooperTrials=0;
597 }
598
599 // Another (sometimes better way) is to use a user-limit maximum Step size
600 // to alleviate this problem ..
601
602 // Add smooth curved trajectories to particle-change
603 //
604 // fParticleChange.SetPointerToVectorOfAuxiliaryPoints
605 // (fFieldPropagator->GimmeTrajectoryVectorAndForgetIt() );
606
607 return &fParticleChange ;
608}
609
610//////////////////////////////////////////////////////////////////////////
611//
612// This ensures that the PostStep action is always called,
613// so that it can do the relocation if it is needed.
614//
615
616G4double G4CoupledTransportation::
617PostStepGetPhysicalInteractionLength( const G4Track&,
618 G4double, // previousStepSize
619 G4ForceCondition* pForceCond )
620{
621 // Must act as PostStep action -- to relocate particle
622 *pForceCond = Forced ;
623 return DBL_MAX ;
624}
625
626static
627void ReportMove( G4ThreeVector OldVector, G4ThreeVector NewVector, G4String& Quantity )
628{
629 G4ThreeVector moveVec = ( NewVector - OldVector );
630
631 G4cerr << G4endl
632 << "**************************************************************" << G4endl;
633 G4cerr << "Endpoint has moved between value expected from TransportEndPosition "
634 << " and value from Track in PostStepDoIt. " << G4endl
635 << "Change of " << Quantity << " is " << moveVec.mag() / mm << " mm long, "
636 << " and its vector is " << (1.0/mm) * moveVec << " mm " << G4endl
637 << "Endpoint of ComputeStep was " << OldVector
638 << " and current position to locate is " << NewVector << G4endl;
639}
640
641/////////////////////////////////////////////////////////////////////////////
642
643G4VParticleChange* G4CoupledTransportation::PostStepDoIt( const G4Track& track,
644 const G4Step& )
645{
646 G4TouchableHandle retCurrentTouchable ; // The one to return
647
648 // Initialize ParticleChange (by setting all its members equal
649 // to corresponding members in G4Track)
650 // fParticleChange.Initialize(track) ; // To initialise TouchableChange
651
652 fParticleChange.ProposeTrackStatus(track.GetTrackStatus()) ;
653
654 // Check that the end position and direction are preserved
655 // since call to AlongStepDoIt
656 if( (fTransportEndPosition - track.GetPosition()).mag2() >= 1.0e-16 ){
657 static G4String EndLabelString("End of Step Position");
658 ReportMove( track.GetPosition(), fTransportEndPosition, EndLabelString );
659 G4cerr << " Problem in G4CoupledTransportation::PostStepDoIt " << G4endl;
660 }
661
662 // If the Step was determined by the volume boundary, relocate the particle
663 // The pathFinder will know that the geometry limited the step (!?)
664
665 if( fVerboseLevel > 0 ){
666 G4cout << " Calling PathFinder::Locate() from "
667 << " G4CoupledTransportation::PostStepDoIt " << G4endl;
668 G4cout << " fAnyGeometryLimitedStep is " << fAnyGeometryLimitedStep << G4endl;
669
670 }
671 if(fAnyGeometryLimitedStep)
672 {
673 fPathFinder->Locate( track.GetPosition(),
674 track.GetMomentumDirection(),
675 true);
676
677 // fCurrentTouchable will now become the previous touchable,
678 // and what was the previous will be freed.
679 // (Needed because the preStepPoint can point to the previous touchable)
680 if( fVerboseLevel > 0 )
681 G4cout << "G4CoupledTransportation::PostStepDoIt --- fNavigatorId = "
682 << fNavigatorId << G4endl;
683
684 fCurrentTouchableHandle=
685 fPathFinder->CreateTouchableHandle( fNavigatorId );
686
687 // Check whether the particle is out of the world volume
688 // If so it has exited and must be killed.
689 //
690#ifdef G4DEBUG_TRANSPORT
691 if( fVerboseLevel > 1 ){
692 G4VPhysicalVolume* vol= fCurrentTouchableHandle->GetVolume();
693 G4cout << "CHECK !!!!!!!!!!! fCurrentTouchableHandle->GetVolume() = " << vol;
694 if( vol ) { G4cout << "Name=" << vol->GetName(); }
695 G4cout << G4endl;
696 }
697#endif
698 if( fCurrentTouchableHandle->GetVolume() == 0 )
699 {
700 fParticleChange.ProposeTrackStatus( fStopAndKill ) ;
701 }
702 retCurrentTouchable = fCurrentTouchableHandle ;
703 // fParticleChange.SetTouchableHandle( fCurrentTouchableHandle ) ;
704
705 // Notify particle change that this is last step in volume
706 fParticleChange.ProposeLastStepInVolume(true);
707 // Double check that a boundary limited the step, and
708 // if( fLinearNavigator->Get
709
710 }
711 else // fAnyGeometryLimitedStep is false
712 {
713#ifdef G4DEBUG_TRANSPORT
714 if( fVerboseLevel > 1 ){
715 G4cout << "G4CoupledTransportation::PostStepDoIt -- "
716 << " fAnyGeometryLimitedStep = " << fAnyGeometryLimitedStep
717 << " must be false " << G4endl;
718 }
719#endif
720 // This serves only to move each of the Navigator's location
721 //
722 // fLinearNavigator->LocateGlobalPointWithinVolume( track.GetPosition() ) ;
723
724 // G4cout << "G4CoupledTransportation calling PathFinder::ReLocate() " << G4endl;
725 fPathFinder->ReLocate( track.GetPosition() );
726 // track.GetMomentumDirection() );
727
728 // Keep the value of the track's current Touchable is retained,
729 // and use it to overwrite the (unset) one in particle change.
730 // Expect this must be fCurrentTouchable too
731 // - could it be different, eg at the start of a step ?
732 //
733 retCurrentTouchable = track.GetTouchableHandle() ;
734 // fParticleChange.SetTouchableHandle( track.GetTouchableHandle() ) ;
735
736 // Have not reached a boundary
737 fParticleChange.ProposeLastStepInVolume(false);
738 } // endif ( fAnyGeometryLimitedStep )
739
740 const G4VPhysicalVolume* pNewVol = retCurrentTouchable->GetVolume() ;
741 const G4Material* pNewMaterial = 0 ;
742 const G4VSensitiveDetector* pNewSensitiveDetector = 0 ;
743
744 if( pNewVol != 0 )
745 {
746 pNewMaterial= pNewVol->GetLogicalVolume()->GetMaterial();
747 pNewSensitiveDetector= pNewVol->GetLogicalVolume()->GetSensitiveDetector();
748 }
749
750 // ( const_cast<G4Material *> pNewMaterial ) ;
751 // ( const_cast<G4VSensitiveDetetor *> pNewSensitiveDetector) ;
752
753 fParticleChange.SetMaterialInTouchable( (G4Material *) pNewMaterial ) ;
754 fParticleChange.SetSensitiveDetectorInTouchable( (G4VSensitiveDetector *) pNewSensitiveDetector ) ;
755 // "temporarily" until Get/Set Material of ParticleChange,
756 // and StepPoint can be made const.
757
758 const G4MaterialCutsCouple* pNewMaterialCutsCouple = 0;
759 if( pNewVol != 0 )
760 {
761 pNewMaterialCutsCouple=pNewVol->GetLogicalVolume()->GetMaterialCutsCouple();
762 if( pNewMaterialCutsCouple!=0
763 && pNewMaterialCutsCouple->GetMaterial()!=pNewMaterial )
764 {
765 // for parametrized volume
766 //
767 pNewMaterialCutsCouple =
768 G4ProductionCutsTable::GetProductionCutsTable()
769 ->GetMaterialCutsCouple(pNewMaterial,
770 pNewMaterialCutsCouple->GetProductionCuts());
771 }
772 }
773 fParticleChange.SetMaterialCutsCoupleInTouchable( pNewMaterialCutsCouple );
774
775 // Must always set the touchable in ParticleChange, whether relocated or not
776 fParticleChange.SetTouchableHandle(retCurrentTouchable) ;
777
778 return &fParticleChange ;
779}
780
781// New method takes over the responsibility to reset the state of
782// G4CoupledTransportation object:
783// - at the start of a new track, and
784// - on the resumption of a suspended track.
785
786void
787G4CoupledTransportation::StartTracking(G4Track* aTrack)
788{
789
790 static G4TransportationManager* transportMgr=
791 G4TransportationManager::GetTransportationManager();
792
793 // G4VProcess::StartTracking(aTrack);
794
795 // The 'initialising' actions
796 // once taken in AlongStepGPIL -- if ( track.GetCurrentStepNumber()==1 )
797
798 // fStartedNewTrack= true;
799
800 fMassNavigator = transportMgr->GetNavigatorForTracking() ;
801 fNavigatorId= transportMgr->ActivateNavigator( fMassNavigator ); // Confirm it!
802
803 // if( fVerboseLevel > 1 ){
804 // G4cout << " Navigator Id obtained in StartTracking " << fNavigatorId << G4endl;
805 // }
806 G4ThreeVector position = aTrack->GetPosition();
807 G4ThreeVector direction = aTrack->GetMomentumDirection();
808
809 // if( fVerboseLevel > 1 ){
810 // G4cout << " Calling PathFinder::PrepareNewTrack from "
811 // << " G4CoupledTransportation::StartTracking -- which calls Locate()" << G4endl;
812 // }
813 fPathFinder->PrepareNewTrack( position, direction);
814 // This implies a call to fPathFinder->Locate( position, direction );
815
816 // Global field, if any, must exist before tracking is started
817 fGlobalFieldExists= DoesGlobalFieldExist();
818 // reset safety value and center
819 //
820 fPreviousMassSafety = 0.0 ;
821 fPreviousFullSafety = 0.0 ;
822 fPreviousSftOrigin = G4ThreeVector(0.,0.,0.) ;
823
824 // reset looping counter -- for motion in field
825 fNoLooperTrials= 0;
826 // Must clear this state .. else it depends on last track's value
827 // --> a better solution would set this from state of suspended track TODO ?
828 // Was if( aTrack->GetCurrentStepNumber()==1 ) { .. }
829
830 // ChordFinder reset internal state
831 //
832 if( fGlobalFieldExists ) {
833 fFieldPropagator->ClearPropagatorState();
834 // Resets safety values, in case of overlaps.
835
836 G4ChordFinder* chordF= fFieldPropagator->GetChordFinder();
837 if( chordF ) chordF->ResetStepEstimate();
838 }
839 // Clear the chord finders of all fields (ie managers) derived objects
840 static G4FieldManagerStore* fieldMgrStore= G4FieldManagerStore::GetInstance();
841 fieldMgrStore->ClearAllChordFindersState();
842
843#ifdef G4DEBUG_TRANSPORT
844 if( fVerboseLevel > 1 ){
845 G4cout << " Returning touchable handle " << fCurrentTouchableHandle << G4endl;
846 }
847#endif
848
849 // Update the current touchable handle (from the track's)
850 //
851 fCurrentTouchableHandle = aTrack->GetTouchableHandle();
852}
853
854void
855G4CoupledTransportation::
856ReportInexactEnergy(G4double startEnergy, G4double endEnergy)
857{
858 static G4int no_warnings= 0, warnModulo=1, moduloFactor= 10, no_large_ediff= 0;
859
860 if( std::fabs(startEnergy- endEnergy) > perThousand * endEnergy )
861 {
862 no_large_ediff ++;
863 if( (no_large_ediff% warnModulo) == 0 )
864 {
865 no_warnings++;
866 G4cout << "WARNING - G4CoupledTransportation::AlongStepGetPIL() "
867 << " Energy change in Step is above 1^-3 relative value. " << G4endl
868 << " Relative change in 'tracking' step = "
869 << std::setw(15) << (endEnergy-startEnergy)/startEnergy << G4endl
870 << " Starting E= " << std::setw(12) << startEnergy / MeV << " MeV " << G4endl
871 << " Ending E= " << std::setw(12) << endEnergy / MeV << " MeV " << G4endl;
872 G4cout << " Energy has been corrected -- however, review"
873 << " field propagation parameters for accuracy." << G4endl;
874 if( (fVerboseLevel > 2 ) || (no_warnings<4) || (no_large_ediff == warnModulo * moduloFactor) ){
875 G4cout << " These include EpsilonStepMax(/Min) in G4FieldManager "
876 << " which determine fractional error per step for integrated quantities. " << G4endl
877 << " Note also the influence of the permitted number of integration steps."
878 << G4endl;
879 }
880 G4cerr << "ERROR - G4CoupledTransportation::AlongStepGetPIL()" << G4endl
881 << " Bad 'endpoint'. Energy change detected"
882 << " and corrected. "
883 << " Has occurred already "
884 << no_large_ediff << " times." << G4endl;
885 if( no_large_ediff == warnModulo * moduloFactor )
886 {
887 warnModulo *= moduloFactor;
888 }
889 }
890 }
891}
Note: See TracBrowser for help on using the repository browser.