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