1 | |
2 | |
3 | |
4 | |
5 | |
6 | |
7 | |
8 | |
9 | |
10 | |
11 | |
12 | |
13 | |
14 | |
15 | |
16 | |
17 | |
18 | |
19 | |
20 | |
21 | |
22 | |
23 | |
24 | |
25 | |
26 | |
27 | |
28 | |
29 | |
30 | |
31 | |
32 | |
33 | |
34 | |
35 | |
36 | |
37 | |
38 | |
39 | |
40 | |
41 | |
42 | |
43 | |
44 | |
45 | |
46 | |
47 | |
48 | #include "G4CoupledTransportation.hh" |
49 | |
50 | #include "G4PhysicalConstants.hh" |
51 | #include "G4SystemOfUnits.hh" |
52 | #include "G4TransportationProcessType.hh" |
53 | #include "G4ProductionCutsTable.hh" |
54 | #include "G4ParticleTable.hh" |
55 | #include "G4ChordFinder.hh" |
56 | #include "G4Field.hh" |
57 | #include "G4FieldTrack.hh" |
58 | #include "G4FieldManagerStore.hh" |
59 | |
60 | class G4VSensitiveDetector; |
61 | |
62 | G4bool G4CoupledTransportation::fSignifyStepInAnyVolume= false; |
63 | |
64 | |
65 | G4bool G4CoupledTransportation::fUseMagneticMoment=false; |
66 | |
67 | |
68 | |
69 | |
70 | G4CoupledTransportation::G4CoupledTransportation( G4int verbosity ) |
71 | : G4VProcess( G4String("CoupledTransportation"), fTransportation ), |
72 | fTransportEndPosition(0.0, 0.0, 0.0), |
73 | fTransportEndMomentumDir(0.0, 0.0, 0.0), |
74 | fTransportEndKineticEnergy(0.0), |
75 | fTransportEndSpin(0.0, 0.0, 0.0), |
76 | fMomentumChanged(false), |
77 | fEndGlobalTimeComputed(false), |
78 | fCandidateEndGlobalTime(0.0), |
79 | fParticleIsLooping( false ), |
80 | fNewTrack( true ), |
81 | fPreviousSftOrigin( 0.,0.,0. ), |
82 | fPreviousMassSafety( 0.0 ), |
83 | fPreviousFullSafety( 0.0 ), |
84 | fMassGeometryLimitedStep( false ), |
85 | fAnyGeometryLimitedStep( false ), |
86 | fEndpointDistance( -1.0 ), |
87 | fThreshold_Warning_Energy( 100 * MeV ), |
88 | fThreshold_Important_Energy( 250 * MeV ), |
89 | fThresholdTrials( 10 ), |
90 | fNoLooperTrials( 0 ), |
91 | fSumEnergyKilled( 0.0 ), fMaxEnergyKilled( 0.0 ), |
92 | fFirstStepInMassVolume( true ), |
93 | fFirstStepInAnyVolume( true ) |
94 | { |
95 | |
96 | SetProcessSubType(static_cast<G4int>(COUPLED_TRANSPORTATION)); |
97 | SetVerboseLevel(verbosity); |
98 | |
99 | G4TransportationManager* transportMgr ; |
100 | |
101 | transportMgr = G4TransportationManager::GetTransportationManager() ; |
102 | |
103 | fMassNavigator = transportMgr->GetNavigatorForTracking() ; |
104 | fFieldPropagator = transportMgr->GetPropagatorInField() ; |
105 | |
106 | fNavigatorId= transportMgr->ActivateNavigator( fMassNavigator ); |
107 | if( verboseLevel > 0 ) |
108 | { |
109 | G4cout(*G4cout_p) << " G4CoupledTransportation constructor: ----- " << G4endlstd::endl; |
110 | G4cout(*G4cout_p) << " Verbose level is " << verboseLevel << G4endlstd::endl; |
111 | G4cout(*G4cout_p) << " Navigator Id obtained in G4CoupledTransportation constructor " |
112 | << fNavigatorId << G4endlstd::endl; |
113 | G4cout(*G4cout_p) << " Reports First/Last in " |
114 | << (fSignifyStepInAnyVolume ? " any " : " mass " ) << " geometry " << G4endlstd::endl; |
115 | } |
116 | fPathFinder= G4PathFinder::GetInstance(); |
117 | fpSafetyHelper = transportMgr->GetSafetyHelper(); |
118 | |
119 | |
120 | static G4ThreadLocalthread_local G4TouchableHandle* pNullTouchableHandle = 0; |
121 | if ( !pNullTouchableHandle) { pNullTouchableHandle = new G4TouchableHandle; } |
122 | fCurrentTouchableHandle = *pNullTouchableHandle; |
123 | |
124 | |
125 | G4FieldManager *globalFieldMgr= transportMgr->GetFieldManager(); |
126 | fGlobalFieldExists= globalFieldMgr ? globalFieldMgr->GetDetectorField() : 0 ; |
127 | } |
128 | |
129 | |
130 | |
131 | G4CoupledTransportation::~G4CoupledTransportation() |
132 | { |
133 | |
134 | |
135 | if( (verboseLevel > 0) || (fSumEnergyKilled > 0.0 ) ) |
136 | { |
137 | G4cout(*G4cout_p) << " G4CoupledTransportation: Statistics for looping particles " << G4endlstd::endl; |
138 | G4cout(*G4cout_p) << " Sum of energy of loopers killed: " << fSumEnergyKilled << G4endlstd::endl; |
139 | G4cout(*G4cout_p) << " Max energy of loopers killed: " << fMaxEnergyKilled << G4endlstd::endl; |
140 | } |
141 | } |
142 | |
143 | |
144 | |
145 | |
146 | |
147 | |
148 | |
149 | |
150 | G4double G4CoupledTransportation:: |
151 | AlongStepGetPhysicalInteractionLength( const G4Track& track, |
152 | G4double, |
153 | G4double currentMinimumStep, |
154 | G4double& proposedSafetyForStart, |
155 | G4GPILSelection* selection ) |
156 | { |
157 | G4double geometryStepLength; |
158 | G4double startMassSafety= 0.0; |
159 | G4double startFullSafety= 0.0; |
160 | G4double safetyProposal= -1.0; |
161 | |
162 | G4ThreeVector EndUnitMomentum ; |
163 | G4double lengthAlongCurve=0.0 ; |
164 | |
165 | fParticleIsLooping = false ; |
166 | |
167 | |
168 | |
169 | |
170 | |
171 | |
172 | |
173 | |
174 | |
175 | |
176 | *selection = CandidateForSelection ; |
177 | |
178 | fFirstStepInMassVolume = fNewTrack | fMassGeometryLimitedStep ; |
179 | fFirstStepInAnyVolume = fNewTrack | fAnyGeometryLimitedStep ; |
180 | |
181 | #ifdef G4DEBUG_TRANSPORT |
182 | G4cout(*G4cout_p) << " CoupledTransport::AlongStep GPIL: " |
183 | << " 1st-step: any= " <<fFirstStepInAnyVolume << " ( geom= " << fAnyGeometryLimitedStep << " ) " |
184 | << " mass= " << fFirstStepInMassVolume << " ( geom= " << fMassGeometryLimitedStep << " ) " |
185 | << " newTrack= " << fNewTrack << G4endlstd::endl; |
186 | |
187 | #endif |
188 | |
189 | |
190 | fNewTrack = false; |
191 | |
192 | |
193 | |
194 | const G4DynamicParticle* pParticle = track.GetDynamicParticle() ; |
195 | const G4ParticleDefinition* pParticleDef = pParticle->GetDefinition() ; |
196 | G4ThreeVector startMomentumDir = pParticle->GetMomentumDirection() ; |
197 | G4ThreeVector startPosition = track.GetPosition() ; |
198 | G4VPhysicalVolume* currentVolume= track.GetVolume(); |
199 | |
200 | #ifdef G4DEBUG_TRANSPORT |
201 | if( verboseLevel > 1 ) |
202 | { |
203 | G4cout(*G4cout_p) << "G4CoupledTransportation::AlongStepGPIL> called in volume " |
204 | << currentVolume->GetName() << G4endlstd::endl; |
205 | } |
206 | #endif |
207 | |
208 | |
209 | |
210 | |
211 | |
212 | |
213 | G4ThreeVector OriginShift = startPosition - fPreviousSftOrigin ; |
214 | G4double MagSqShift = OriginShift.mag2() ; |
215 | startMassSafety = 0.0; |
216 | startFullSafety= 0.0; |
217 | |
218 | |
219 | |
220 | if( MagSqShift < sqr(fPreviousFullSafety) ) |
221 | { |
222 | G4double mag_shift= std::sqrt(MagSqShift); |
223 | startMassSafety = std::max( (fPreviousMassSafety - mag_shift), 0.0); |
224 | startFullSafety = std::max( (fPreviousFullSafety - mag_shift), 0.0); |
225 | |
226 | |
227 | |
228 | |
229 | |
230 | } |
231 | |
232 | |
233 | |
234 | G4double particleCharge = pParticle->GetCharge() ; |
235 | G4double magneticMoment = pParticle->GetMagneticMoment() ; |
236 | G4double restMass = pParticle->GetMass() ; |
237 | |
238 | fMassGeometryLimitedStep = false ; |
239 | fAnyGeometryLimitedStep = false; |
240 | |
241 | |
242 | |
243 | |
244 | |
245 | |
246 | |
247 | |
248 | |
249 | G4FieldManager* fieldMgr=0; |
250 | G4bool fieldExertsForce = false ; |
251 | |
252 | G4bool gravityOn = false; |
253 | const G4Field* ptrField= 0; |
254 | |
255 | fieldMgr = fFieldPropagator->FindAndSetFieldManager( track.GetVolume() ); |
256 | if( fieldMgr != 0 ) |
257 | { |
258 | |
259 | fieldMgr->ConfigureForTrack( &track ); |
260 | |
261 | |
262 | |
263 | |
264 | ptrField= fieldMgr->GetDetectorField(); |
265 | |
266 | if( ptrField != 0) |
267 | { |
268 | gravityOn= ptrField->IsGravityActive(); |
269 | if( (particleCharge != 0.0) |
270 | || (fUseMagneticMoment && (magneticMoment != 0.0) ) |
271 | || (gravityOn && (restMass != 0.0)) |
272 | ) |
273 | { |
274 | fieldExertsForce = true; |
275 | } |
276 | } |
277 | } |
278 | G4double momentumMagnitude = pParticle->GetTotalMomentum() ; |
279 | |
280 | if( fieldExertsForce ) |
281 | { |
282 | auto equationOfMotion= fFieldPropagator->GetCurrentEquationOfMotion(); |
283 | |
284 | G4ChargeState chargeState(particleCharge, |
285 | magneticMoment, |
286 | pParticleDef->GetPDGSpin() ); |
287 | |
288 | |
289 | |
290 | if( equationOfMotion ) |
291 | { |
292 | equationOfMotion->SetChargeMomentumMass( chargeState, |
293 | momentumMagnitude, |
294 | restMass ); |
295 | } |
296 | #ifdef G4DEBUG_TRANSPORT |
297 | else |
298 | { |
299 | G4cerr(*G4cerr_p) << " ERROR in G4CoupledTransportation> " |
300 | << "Cannot find valid Equation of motion: " << G4endlstd::endl; |
301 | << " Unable to pass Charge, Momentum and Mass " << G4endlstd::endl; |
302 | } |
303 | #endif |
304 | } |
305 | |
306 | G4ThreeVector polarizationVec = track.GetPolarization() ; |
307 | G4FieldTrack aFieldTrack = G4FieldTrack( startPosition, |
308 | track.GetGlobalTime(), |
309 | |
310 | track.GetMomentumDirection(), |
311 | track.GetKineticEnergy(), |
312 | restMass, |
313 | particleCharge, |
314 | polarizationVec, |
315 | pParticleDef->GetPDGMagneticMoment(), |
316 | 0.0, |
317 | pParticleDef->GetPDGSpin() |
318 | ) ; |
319 | G4int stepNo= track.GetCurrentStepNumber(); |
320 | |
321 | ELimited limitedStep; |
322 | G4FieldTrack endTrackState('a'); |
323 | |
324 | fMassGeometryLimitedStep = false ; |
325 | fAnyGeometryLimitedStep = false ; |
326 | if( currentMinimumStep > 0 ) |
327 | { |
328 | G4double newMassSafety= 0.0; |
329 | |
330 | |
331 | |
332 | lengthAlongCurve = fPathFinder->ComputeStep( aFieldTrack, |
333 | currentMinimumStep, |
334 | fNavigatorId, |
335 | stepNo, |
336 | newMassSafety, |
337 | limitedStep, |
338 | endTrackState, |
339 | currentVolume ) ; |
340 | |
341 | |
342 | G4double newFullSafety= fPathFinder->GetCurrentSafety(); |
343 | |
344 | |
345 | |
346 | if( limitedStep == kUnique || limitedStep == kSharedTransport ) |
347 | { |
348 | fMassGeometryLimitedStep = true ; |
349 | } |
350 | |
351 | fAnyGeometryLimitedStep = (fPathFinder->GetNumberGeometriesLimitingStep() != 0) ; |
352 | |
353 | #ifdef G4DEBUG_TRANSPORT |
354 | if( fMassGeometryLimitedStep && !fAnyGeometryLimitedStep ) |
355 | { |
356 | std::ostringstream message; |
357 | message << " ERROR in determining geometries limiting the step" << G4endlstd::endl; |
358 | message << " Limiting: mass=" << fMassGeometryLimitedStep |
359 | << " any= " << fAnyGeometryLimitedStep << G4endlstd::endl; |
360 | message << "Incompatible conditions - by which geometries was it limited ?"<<G4endlstd::endl; |
361 | G4Exception("G4CoupledTransportation::AlongStepGetPhysicalInteractionLength()", |
362 | "PathFinderConfused", FatalException, message); |
363 | } |
364 | #endif |
365 | |
366 | geometryStepLength = std::min( lengthAlongCurve, currentMinimumStep); |
367 | |
368 | |
369 | |
370 | fMomentumChanged = true ; |
371 | fTransportEndMomentumDir = endTrackState.GetMomentumDir() ; |
372 | |
373 | |
374 | fPreviousSftOrigin = startPosition ; |
375 | fPreviousMassSafety = newMassSafety ; |
376 | fPreviousFullSafety = newFullSafety ; |
377 | |
378 | |
379 | #ifdef G4DEBUG_TRANSPORT |
380 | if( verboseLevel > 1 ) |
381 | { |
382 | G4cout(*G4cout_p) << "G4Transport:CompStep> " |
383 | << " called the pathfinder for a new step at " << startPosition |
384 | << " and obtained step = " << lengthAlongCurve << G4endlstd::endl; |
385 | G4cout(*G4cout_p) << " New safety (preStep) = " << newMassSafety |
386 | << " versus precalculated = " << startMassSafety << G4endlstd::endl; |
387 | } |
388 | #endif |
389 | |
390 | |
391 | startMassSafety = newMassSafety ; |
| Value stored to 'startMassSafety' is never read |
392 | startFullSafety = newFullSafety ; |
393 | |
394 | |
395 | fTransportEndPosition = endTrackState.GetPosition() ; |
396 | fTransportEndKineticEnergy = endTrackState.GetKineticEnergy() ; |
397 | } |
398 | else |
399 | { |
400 | geometryStepLength = lengthAlongCurve= 0.0 ; |
401 | fMomentumChanged = false ; |
402 | |
403 | |
404 | fTransportEndMomentumDir = track.GetMomentumDirection(); |
405 | fTransportEndKineticEnergy = track.GetKineticEnergy(); |
406 | |
407 | fTransportEndPosition = startPosition; |
408 | |
409 | endTrackState= aFieldTrack; |
410 | |
411 | |
412 | |
413 | if( startMassSafety == 0.0 ) |
414 | { |
415 | fMassGeometryLimitedStep = true ; |
416 | fAnyGeometryLimitedStep = true; |
417 | } |
418 | |
419 | } |
420 | |
421 | |
422 | if( !fieldExertsForce ) |
423 | { |
424 | fParticleIsLooping = false ; |
425 | fMomentumChanged = false ; |
426 | fEndGlobalTimeComputed = false ; |
427 | |
428 | } |
429 | else |
430 | { |
431 | |
432 | #ifdef G4DEBUG_TRANSPORT |
433 | if( verboseLevel > 1 ) |
434 | { |
435 | G4cout(*G4cout_p) << " G4CT::CS End Position = " << fTransportEndPosition << G4endlstd::endl; |
436 | G4cout(*G4cout_p) << " G4CT::CS End Direction = " << fTransportEndMomentumDir << G4endlstd::endl; |
437 | } |
438 | #endif |
439 | if( fFieldPropagator->GetCurrentFieldManager()->DoesFieldChangeEnergy() ) |
440 | { |
441 | |
442 | |
443 | |
444 | fCandidateEndGlobalTime = endTrackState.GetLabTimeOfFlight(); |
445 | fEndGlobalTimeComputed = true; |
446 | |
447 | |
448 | |
449 | } |
450 | else |
451 | { |
452 | |
453 | |
454 | |
455 | fEndGlobalTimeComputed = false; |
456 | |
457 | |
458 | |
459 | G4double startEnergy= track.GetKineticEnergy(); |
460 | G4double endEnergy= fTransportEndKineticEnergy; |
461 | |
462 | static G4ThreadLocalthread_local G4int no_inexact_steps=0; |
463 | G4double absEdiff = std::fabs(startEnergy- endEnergy); |
464 | if( absEdiff > perMillion * endEnergy ) |
465 | { |
466 | no_inexact_steps++; |
467 | |
468 | } |
469 | #ifdef G4VERBOSE1 |
470 | if( (verboseLevel > 1) && ( absEdiff > perThousand * endEnergy) ) |
471 | { |
472 | ReportInexactEnergy(startEnergy, endEnergy); |
473 | } |
474 | #endif |
475 | |
476 | |
477 | |
478 | fTransportEndKineticEnergy= track.GetKineticEnergy(); |
479 | } |
480 | } |
481 | |
482 | fEndpointDistance = (fTransportEndPosition - startPosition).mag() ; |
483 | fParticleIsLooping = fFieldPropagator->IsParticleLooping() ; |
484 | |
485 | fTransportEndSpin = endTrackState.GetSpin(); |
486 | |
487 | |
488 | safetyProposal= startFullSafety; |
489 | |
490 | |
491 | |
492 | |
493 | if( (startFullSafety < fEndpointDistance ) |
494 | && ( particleCharge != 0.0 ) ) |
495 | |
496 | { |
497 | G4double endFullSafety = |
498 | fPathFinder->ComputeSafety( fTransportEndPosition); |
499 | |
500 | |
501 | |
502 | |
503 | |
504 | |
505 | fpSafetyHelper->SetCurrentSafety( endFullSafety, fTransportEndPosition); |
506 | |
507 | |
508 | G4ThreeVector centerPt= G4ThreeVector(0.0, 0.0, 0.0); |
509 | G4double endMassSafety= fPathFinder->ObtainSafety( fNavigatorId, centerPt); |
510 | |
511 | |
512 | fPreviousMassSafety = endMassSafety ; |
513 | fPreviousFullSafety = endFullSafety; |
514 | fPreviousSftOrigin = fTransportEndPosition ; |
515 | |
516 | |
517 | |
518 | safetyProposal = endFullSafety + fEndpointDistance; |
519 | |
520 | |
521 | |
522 | |
523 | |
524 | #ifdef G4DEBUG_TRANSPORT |
525 | G4int prec= G4cout(*G4cout_p).precision(12) ; |
526 | G4cout(*G4cout_p) << "***Transportation::AlongStepGPIL ** " << G4endlstd::endl ; |
527 | G4cout(*G4cout_p) << " Revised Safety at endpoint " << fTransportEndPosition |
528 | << " give safety values: Mass= " << endMassSafety |
529 | << " All= " << endFullSafety << G4endlstd::endl ; |
530 | G4cout(*G4cout_p) << " Adding endpoint distance " << fEndpointDistance |
531 | << " to obtain pseudo-safety= " << safetyProposal << G4endlstd::endl ; |
532 | G4cout(*G4cout_p).precision(prec); |
533 | } |
534 | else |
535 | { |
536 | G4int prec= G4cout(*G4cout_p).precision(12) ; |
537 | G4cout(*G4cout_p) << "***Transportation::AlongStepGPIL ** " << G4endlstd::endl ; |
538 | G4cout(*G4cout_p) << " Quick Safety estimate at endpoint " << fTransportEndPosition |
539 | << " gives safety endpoint value = " << startFullSafety - fEndpointDistance |
540 | << " using start-point value " << startFullSafety |
541 | << " and endpointDistance " << fEndpointDistance << G4endlstd::endl; |
542 | G4cout(*G4cout_p).precision(prec); |
543 | #endif |
544 | } |
545 | |
546 | proposedSafetyForStart= safetyProposal; |
547 | fParticleChange.ProposeTrueStepLength(geometryStepLength) ; |
548 | |
549 | return geometryStepLength ; |
550 | } |
551 | |
552 | |
553 | |
554 | G4VParticleChange* |
555 | G4CoupledTransportation::AlongStepDoIt( const G4Track& track, |
556 | const G4Step& stepData ) |
557 | { |
558 | static G4ThreadLocalthread_local G4int noCalls=0; |
559 | noCalls++; |
560 | |
561 | fParticleChange.Initialize(track) ; |
562 | |
563 | |
564 | |
565 | |
566 | fParticleChange.ProposePosition(fTransportEndPosition) ; |
567 | |
568 | |
569 | |
570 | fParticleChange.ProposeMomentumDirection(fTransportEndMomentumDir) ; |
571 | fParticleChange.ProposeEnergy(fTransportEndKineticEnergy) ; |
572 | fParticleChange.SetMomentumChanged(fMomentumChanged) ; |
573 | |
574 | fParticleChange.ProposePolarization(fTransportEndSpin); |
575 | |
576 | G4double deltaTime = 0.0 ; |
577 | |
578 | |
579 | |
580 | |
581 | |
582 | G4double startTime = track.GetGlobalTime() ; |
583 | |
584 | if (!fEndGlobalTimeComputed) |
585 | { |
586 | G4double finalInverseVel= DBL_MAXstd::numeric_limits<double>::max(), initialInverseVel=DBL_MAXstd::numeric_limits<double>::max(); |
587 | |
588 | |
589 | |
590 | G4double finalVelocity = track.GetVelocity() ; |
591 | if( finalVelocity > 0.0 ) { finalInverseVel= 1.0 / finalVelocity; } |
592 | G4double initialVelocity = stepData.GetPreStepPoint()->GetVelocity() ; |
593 | if( initialVelocity > 0.0 ) { initialInverseVel= 1.0 / initialVelocity; } |
594 | G4double stepLength = track.GetStepLength() ; |
595 | |
596 | if (finalVelocity > 0.0) |
597 | { |
598 | |
599 | G4double meanInverseVelocity = 0.5 * ( initialInverseVel + finalInverseVel ); |
600 | deltaTime = stepLength * meanInverseVelocity ; |
601 | |
602 | |
603 | } |
604 | else |
605 | { |
606 | deltaTime = stepLength * initialInverseVel ; |
607 | |
608 | |
609 | } |
610 | |
611 | fCandidateEndGlobalTime = startTime + deltaTime ; |
612 | fParticleChange.ProposeLocalTime( track.GetLocalTime() + deltaTime) ; |
613 | |
614 | |
615 | |
616 | } |
617 | else |
618 | { |
619 | deltaTime = fCandidateEndGlobalTime - startTime ; |
620 | fParticleChange.ProposeGlobalTime( fCandidateEndGlobalTime ) ; |
621 | |
622 | |
623 | } |
624 | |
625 | |
626 | |
627 | |
628 | |
629 | |
630 | |
631 | G4double restMass = track.GetDynamicParticle()->GetMass() ; |
632 | G4double deltaProperTime = deltaTime*( restMass/track.GetTotalEnergy() ) ; |
633 | |
634 | fParticleChange.ProposeProperTime(track.GetProperTime() + deltaProperTime) ; |
635 | |
636 | |
637 | |
638 | |
639 | |
640 | |
641 | if ( fParticleIsLooping ) |
642 | { |
643 | G4double endEnergy= fTransportEndKineticEnergy; |
644 | |
645 | if( (endEnergy < fThreshold_Important_Energy) |
646 | || (fNoLooperTrials >= fThresholdTrials ) ) |
647 | { |
648 | |
649 | |
650 | fParticleChange.ProposeTrackStatus( fStopAndKill ) ; |
651 | |
652 | |
653 | fSumEnergyKilled += endEnergy; |
654 | if( endEnergy > fMaxEnergyKilled) { fMaxEnergyKilled= endEnergy; } |
655 | |
656 | #ifdef G4VERBOSE1 |
657 | if((verboseLevel > 1) && ( endEnergy > fThreshold_Warning_Energy )) |
658 | { |
659 | G4cout(*G4cout_p) << " G4CoupledTransportation is killing track that is looping or stuck " << G4endlstd::endl |
660 | << " This track has " << track.GetKineticEnergy() / MeV |
661 | << " MeV energy." << G4endlstd::endl; |
662 | } |
663 | if( verboseLevel > 0 ) |
664 | { |
665 | G4cout(*G4cout_p) << " Steps by this track: " << track.GetCurrentStepNumber() << G4endlstd::endl; |
666 | } |
667 | #endif |
668 | fNoLooperTrials=0; |
669 | } |
670 | else |
671 | { |
672 | fNoLooperTrials ++; |
673 | #ifdef G4VERBOSE1 |
674 | if( (verboseLevel > 2) ) |
675 | { |
676 | G4cout(*G4cout_p) << " ** G4CoupledTransportation::AlongStepDoIt(): Particle looping - " << G4endlstd::endl |
677 | << " Number of consecutive problem step (this track) = " << fNoLooperTrials << G4endlstd::endl |
678 | << " Steps by this track: " << track.GetCurrentStepNumber() << G4endlstd::endl |
679 | << " Total no of calls to this method (all tracks) = " << noCalls << G4endlstd::endl; |
680 | } |
681 | #endif |
682 | } |
683 | } |
684 | else |
685 | { |
686 | fNoLooperTrials=0; |
687 | } |
688 | |
689 | |
690 | |
691 | |
692 | |
693 | |
694 | |
695 | |
696 | |
697 | return &fParticleChange ; |
698 | } |
699 | |
700 | |
701 | |
702 | |
703 | |
704 | |
705 | |
706 | G4double G4CoupledTransportation:: |
707 | PostStepGetPhysicalInteractionLength( const G4Track&, |
708 | G4double, |
709 | G4ForceCondition* pForceCond ) |
710 | { |
711 | |
712 | *pForceCond = Forced ; |
713 | return DBL_MAXstd::numeric_limits<double>::max() ; |
714 | } |
715 | |
716 | void G4CoupledTransportation:: |
717 | ReportMove( G4ThreeVector OldVector, G4ThreeVector NewVector, const G4String& Quantity ) |
718 | { |
719 | G4ThreeVector moveVec = ( NewVector - OldVector ); |
720 | |
721 | G4cerr(*G4cerr_p) << G4endlstd::endl |
722 | << "**************************************************************" << G4endlstd::endl; |
723 | G4cerr(*G4cerr_p) << "Endpoint has moved between value expected from TransportEndPosition " |
724 | << " and value from Track in PostStepDoIt. " << G4endlstd::endl |
725 | << "Change of " << Quantity << " is " << moveVec.mag() / mm << " mm long, " |
726 | << " and its vector is " << (1.0/mm) * moveVec << " mm " << G4endlstd::endl |
727 | << "Endpoint of ComputeStep was " << OldVector |
728 | << " and current position to locate is " << NewVector << G4endlstd::endl; |
729 | } |
730 | |
731 | |
732 | |
733 | G4VParticleChange* G4CoupledTransportation::PostStepDoIt( const G4Track& track, |
734 | const G4Step& ) |
735 | { |
736 | G4TouchableHandle retCurrentTouchable ; |
737 | |
738 | |
739 | |
740 | |
741 | |
742 | fParticleChange.ProposeTrackStatus(track.GetTrackStatus()) ; |
743 | |
744 | |
745 | |
746 | if( fSignifyStepInAnyVolume ){ |
747 | fParticleChange.ProposeFirstStepInVolume( fFirstStepInAnyVolume ); |
748 | |
749 | }else{ |
750 | fParticleChange.ProposeFirstStepInVolume( fFirstStepInMassVolume ); |
751 | |
752 | } |
753 | |
754 | |
755 | |
756 | |
757 | #ifdef G4DEBUG_TRANSPORT |
758 | if( ( verboseLevel > 0 ) |
759 | && ((fTransportEndPosition - track.GetPosition()).mag2() >= 1.0e-16) ) |
760 | { |
761 | ReportMove( track.GetPosition(), fTransportEndPosition, "End of Step Position" ); |
762 | G4cerr(*G4cerr_p) << " Problem in G4CoupledTransportation::PostStepDoIt " << G4endlstd::endl; |
763 | } |
764 | |
765 | |
766 | |
767 | |
768 | if( verboseLevel > 0 ) |
769 | { |
770 | G4cout(*G4cout_p) << " Calling PathFinder::Locate() from " |
771 | << " G4CoupledTransportation::PostStepDoIt " << G4endlstd::endl; |
772 | G4cout(*G4cout_p) << " fAnyGeometryLimitedStep is " << fAnyGeometryLimitedStep << G4endlstd::endl; |
773 | |
774 | } |
775 | #endif |
776 | |
777 | if(fAnyGeometryLimitedStep) |
778 | { |
779 | fPathFinder->Locate( track.GetPosition(), |
780 | track.GetMomentumDirection(), |
781 | true); |
782 | |
783 | |
784 | |
785 | |
786 | |
787 | fCurrentTouchableHandle= |
788 | fPathFinder->CreateTouchableHandle( fNavigatorId ); |
789 | |
790 | #ifdef G4DEBUG_TRANSPORT |
791 | if( verboseLevel > 0 ) |
792 | { |
793 | G4cout(*G4cout_p) << "G4CoupledTransportation::PostStepDoIt --- fNavigatorId = " |
794 | << fNavigatorId << G4endlstd::endl; |
795 | } |
796 | if( verboseLevel > 1 ) |
797 | { |
798 | G4VPhysicalVolume* vol= fCurrentTouchableHandle->GetVolume(); |
799 | G4cout(*G4cout_p) << "CHECK !!!!!!!!!!! fCurrentTouchableHandle->GetVolume() = " << vol; |
800 | if( vol ) { G4cout(*G4cout_p) << "Name=" << vol->GetName(); } |
801 | G4cout(*G4cout_p) << G4endlstd::endl; |
802 | } |
803 | #endif |
804 | |
805 | |
806 | |
807 | |
808 | if( fCurrentTouchableHandle->GetVolume() == 0 ) |
809 | { |
810 | fParticleChange.ProposeTrackStatus( fStopAndKill ) ; |
811 | } |
812 | retCurrentTouchable = fCurrentTouchableHandle ; |
813 | |
814 | } |
815 | else |
816 | { |
817 | #ifdef G4DEBUG_TRANSPORT |
818 | if( verboseLevel > 1 ) |
819 | { |
820 | G4cout(*G4cout_p) << "G4CoupledTransportation::PostStepDoIt -- " |
821 | << " fAnyGeometryLimitedStep = " << fAnyGeometryLimitedStep |
822 | << " must be false " << G4endlstd::endl; |
823 | } |
824 | #endif |
825 | |
826 | |
827 | |
828 | |
829 | |
830 | fPathFinder->ReLocate( track.GetPosition() ); |
831 | |
832 | |
833 | |
834 | |
835 | |
836 | |
837 | |
838 | retCurrentTouchable = track.GetTouchableHandle() ; |
839 | |
840 | } |
841 | |
842 | #ifdef G4DEBUG_NAVIGATION |
843 | G4cout(*G4cout_p) << " CoupledTransport::AlongStep GPIL: " |
844 | << " last-step: any= " << fAnyGeometryLimitedStep << " . ..... x . " |
845 | << " mass= " << fMassGeometryLimitedStep |
846 | << G4endlstd::endl; |
847 | #endif |
848 | |
849 | if( fSignifyStepInAnyVolume ) |
850 | fParticleChange.ProposeLastStepInVolume(fAnyGeometryLimitedStep); |
851 | else |
852 | fParticleChange.ProposeLastStepInVolume(fMassGeometryLimitedStep); |
853 | |
854 | const G4VPhysicalVolume* pNewVol = retCurrentTouchable->GetVolume() ; |
855 | const G4Material* pNewMaterial = 0 ; |
856 | const G4VSensitiveDetector* pNewSensitiveDetector = 0 ; |
857 | |
858 | if( pNewVol != 0 ) |
859 | { |
860 | pNewMaterial= pNewVol->GetLogicalVolume()->GetMaterial(); |
861 | pNewSensitiveDetector= pNewVol->GetLogicalVolume()->GetSensitiveDetector(); |
862 | } |
863 | |
864 | |
865 | |
866 | |
867 | fParticleChange.SetMaterialInTouchable( (G4Material *) pNewMaterial ) ; |
868 | fParticleChange.SetSensitiveDetectorInTouchable( (G4VSensitiveDetector *) pNewSensitiveDetector ) ; |
869 | |
870 | |
871 | |
872 | const G4MaterialCutsCouple* pNewMaterialCutsCouple = 0; |
873 | if( pNewVol != 0 ) |
874 | { |
875 | pNewMaterialCutsCouple=pNewVol->GetLogicalVolume()->GetMaterialCutsCouple(); |
876 | if( pNewMaterialCutsCouple!=0 |
877 | && pNewMaterialCutsCouple->GetMaterial()!=pNewMaterial ) |
878 | { |
879 | |
880 | |
881 | pNewMaterialCutsCouple = |
882 | G4ProductionCutsTable::GetProductionCutsTable() |
883 | ->GetMaterialCutsCouple(pNewMaterial, |
884 | pNewMaterialCutsCouple->GetProductionCuts()); |
885 | } |
886 | } |
887 | fParticleChange.SetMaterialCutsCoupleInTouchable( pNewMaterialCutsCouple ); |
888 | |
889 | |
890 | fParticleChange.SetTouchableHandle(retCurrentTouchable) ; |
891 | |
892 | return &fParticleChange ; |
893 | } |
894 | |
895 | |
896 | |
897 | |
898 | |
899 | |
900 | void |
901 | G4CoupledTransportation::StartTracking(G4Track* aTrack) |
902 | { |
903 | |
904 | G4TransportationManager* transportMgr = |
905 | G4TransportationManager::GetTransportationManager(); |
906 | |
907 | |
908 | fNewTrack= true; |
909 | |
910 | |
911 | |
912 | |
913 | |
914 | |
915 | fMassNavigator = transportMgr->GetNavigatorForTracking() ; |
916 | fNavigatorId= transportMgr->ActivateNavigator( fMassNavigator ); |
917 | |
918 | |
919 | |
920 | |
921 | G4ThreeVector position = aTrack->GetPosition(); |
922 | G4ThreeVector direction = aTrack->GetMomentumDirection(); |
923 | |
924 | |
925 | |
926 | |
927 | |
928 | fPathFinder->PrepareNewTrack( position, direction); |
929 | |
930 | |
931 | |
932 | fGlobalFieldExists= DoesGlobalFieldExist(); |
933 | |
934 | |
935 | fPreviousMassSafety = 0.0 ; |
936 | fPreviousFullSafety = 0.0 ; |
937 | fPreviousSftOrigin = G4ThreeVector(0.,0.,0.) ; |
938 | |
939 | |
940 | fNoLooperTrials= 0; |
941 | |
942 | |
943 | |
944 | |
945 | |
946 | |
947 | if( fFieldPropagator ) |
948 | { |
949 | fFieldPropagator->ClearPropagatorState(); |
950 | |
951 | |
952 | G4ChordFinder* chordF= fFieldPropagator->GetChordFinder(); |
953 | if( chordF ) { chordF->ResetStepEstimate(); } |
954 | } |
955 | |
956 | |
957 | |
958 | G4FieldManagerStore* fieldMgrStore = G4FieldManagerStore::GetInstance(); |
959 | fieldMgrStore->ClearAllChordFindersState(); |
960 | |
961 | #ifdef G4DEBUG_TRANSPORT |
962 | if( verboseLevel > 1 ) |
963 | { |
964 | G4cout(*G4cout_p) << " Returning touchable handle " << fCurrentTouchableHandle << G4endlstd::endl; |
965 | } |
966 | #endif |
967 | |
968 | |
969 | |
970 | fCurrentTouchableHandle = aTrack->GetTouchableHandle(); |
971 | } |
972 | |
973 | void |
974 | G4CoupledTransportation::EndTracking() |
975 | { |
976 | G4TransportationManager::GetTransportationManager()->InactivateAll(); |
977 | fPathFinder->EndTrack(); |
978 | } |
979 | |
980 | void |
981 | G4CoupledTransportation:: |
982 | ReportInexactEnergy(G4double startEnergy, G4double endEnergy) |
983 | { |
984 | static G4ThreadLocalthread_local G4int no_warnings= 0, warnModulo=1, moduloFactor= 10, no_large_ediff= 0; |
985 | |
986 | if( std::fabs(startEnergy- endEnergy) > perThousand * endEnergy ) |
987 | { |
988 | no_large_ediff ++; |
989 | if( (no_large_ediff% warnModulo) == 0 ) |
990 | { |
991 | no_warnings++; |
992 | G4cout(*G4cout_p) << "WARNING - G4CoupledTransportation::AlongStepGetPIL() " |
993 | << " Energy change in Step is above 1^-3 relative value. " << G4endlstd::endl |
994 | << " Relative change in 'tracking' step = " |
995 | << std::setw(15) << (endEnergy-startEnergy)/startEnergy << G4endlstd::endl |
996 | << " Starting E= " << std::setw(12) << startEnergy / MeV << " MeV " << G4endlstd::endl |
997 | << " Ending E= " << std::setw(12) << endEnergy / MeV << " MeV " << G4endlstd::endl; |
998 | G4cout(*G4cout_p) << " Energy has been corrected -- however, review" |
999 | << " field propagation parameters for accuracy." << G4endlstd::endl; |
1000 | if( (verboseLevel > 2 ) || (no_warnings<4) || (no_large_ediff == warnModulo * moduloFactor) ) |
1001 | { |
1002 | G4cout(*G4cout_p) << " These include EpsilonStepMax(/Min) in G4FieldManager " |
1003 | << " which determine fractional error per step for integrated quantities. " << G4endlstd::endl |
1004 | << " Note also the influence of the permitted number of integration steps." |
1005 | << G4endlstd::endl; |
1006 | } |
1007 | G4cerr(*G4cerr_p) << "ERROR - G4CoupledTransportation::AlongStepGetPIL()" << G4endlstd::endl |
1008 | << " Bad 'endpoint'. Energy change detected" |
1009 | << " and corrected. " |
1010 | << " Has occurred already " |
1011 | << no_large_ediff << " times." << G4endlstd::endl; |
1012 | if( no_large_ediff == warnModulo * moduloFactor ) |
1013 | { |
1014 | warnModulo *= moduloFactor; |
1015 | } |
1016 | } |
1017 | } |
1018 | } |
1019 | |
1020 | #include "G4Transportation.hh" |
1021 | G4bool G4CoupledTransportation::EnableUseMagneticMoment(G4bool useMoment) |
1022 | { |
1023 | G4bool lastValue= fUseMagneticMoment; |
1024 | fUseMagneticMoment= useMoment; |
1025 | G4Transportation::fUseMagneticMoment= useMoment; |
1026 | return lastValue; |
1027 | } |