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

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

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