| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 //
27 //
28 // -------------------------------------------
29 // GEANT 4 class implementation
30 //
31 // ===========================================
32 // Created: 19 March 1997, J. Apostolakis
33 // ===========================================
34
35 #include "G4CoupledTransportation.hh"
36 #include "G4TransportationProcessType.hh"
37 #include "G4TransportationLogger.hh"
38
39 #include "G4PhysicalConstants.hh"
40 #include "G4SystemOfUnits.hh"
41 #include "G4ProductionCutsTable.hh"
42 #include "G4ParticleTable.hh"
43 #include "G4ChordFinder.hh"
44 #include "G4Field.hh"
45 #include "G4FieldTrack.hh"
46 #include "G4FieldManagerStore.hh"
47 #include "G4PathFinder.hh"
48
49 #include "G4PropagatorInField.hh"
50 #include "G4TransportationManager.hh"
51
52 class G4VSensitiveDetector;
53
54 G4bool G4CoupledTransportation::fSignifyStepIn
55 // This mode must apply to all threads
56
57 //////////////////////////////////////////////
58 //
59 // Constructor
60
61 G4CoupledTransportation::G4CoupledTransportati
62 : G4Transportation( verbosity, "CoupledTrans
63 fPreviousMassSafety( 0.0 ),
64 fPreviousFullSafety( 0.0 ),
65 fMassGeometryLimitedStep( false ),
66 fFirstStepInMassVolume( true )
67 {
68 SetProcessSubType(static_cast<G4int>(COUPLED
69 // SetVerboseLevel is called in the construc
70
71 if( verboseLevel > 0 )
72 {
73 G4cout << " G4CoupledTransportation constr
74 G4cout << " Verbose level is " << verboseL
75 G4cout << " Reports First/Last in "
76 << (fSignifyStepInAnyVolume ? " any
77 << " geometry " << G4endl;
78 }
79 fPathFinder= G4PathFinder::GetInstance();
80 }
81
82 //////////////////////////////////////////////
83
84 G4CoupledTransportation::~G4CoupledTransportat
85 {
86 }
87
88 //////////////////////////////////////////////
89 //
90 // Responsibilities:
91 // Find whether the geometry limits the Ste
92 // Calculate the new value of the safety an
93 // Store the final time, position and momen
94
95 G4double G4CoupledTransportation::
96 AlongStepGetPhysicalInteractionLength( const G
97 G
98 G
99 G
100 G
101 {
102 G4double geometryStepLength;
103 G4double startFullSafety= 0.0; // estimated
104 G4double safetyProposal= -1.0; // local copy
105
106 G4ThreeVector EndUnitMomentum ;
107 G4double lengthAlongCurve = 0.0 ;
108
109 fParticleIsLooping = false ;
110
111 // Initial actions moved to StartTrack()
112 // --------------------------------------
113 // Note: in case another process changes tou
114 // it will be necessary to add here (for
115 // fCurrentTouchableHandle = aTrack->GetTouc
116
117 // GPILSelection is set to defaule value of
118 // It is a return value
119 //
120 *selection = CandidateForSelection ;
121
122 fFirstStepInMassVolume = fNewTrack || fMassG
123 fFirstStepInVolume = fNewTrack || fGeome
124
125 #ifdef G4DEBUG_TRANSPORT
126 G4cout << " CoupledTransport::AlongStep GPI
127 << " 1st-step: any= " <<fFirstStep
128 << fGeometryLimitedStep << " ) "
129 << " mass= " << fFirstStepI
130 << fMassGeometryLimitedStep << " ) "
131 << " newTrack= " << fNewTrack << G4e
132 #endif
133
134 // fLastStepInVolume= false;
135 fNewTrack = false;
136
137 // Get initial Energy/Momentum of the track
138 //
139 const G4DynamicParticle* pParticle = tra
140 const G4ParticleDefinition* pParticleDef =
141 G4ThreeVector startMomentumDir = pPart
142 G4ThreeVector startPosition = track
143 G4VPhysicalVolume* currentVolume= track.GetV
144
145 #ifdef G4DEBUG_TRANSPORT
146 if( verboseLevel > 1 )
147 {
148 G4cout << "G4CoupledTransportation::AlongS
149 << currentVolume->GetName() << G4en
150 }
151 #endif
152 // G4double theTime = track.GetGlob
153
154 // The Step Point safety can be limited by o
155 // assumptions of any process - it's not alw
156 // We calculate the starting point's isotrop
157 //
158 G4ThreeVector OriginShift = startPosition -
159 G4double MagSqShift = OriginShift.mag2
160 startFullSafety= 0.0;
161
162 // Recall that FullSafety <= MassSafety
163 // Original: if( MagSqShift < sqr(fPreviousM
164 if( MagSqShift < sqr(fPreviousFullSafety) )
165 {
166 G4double mag_shift= std::sqrt(MagSqShift)
167 startFullSafety = std::max( (fPreviousFul
168 // Need to be consistent between full s
169 // in order reproduce results in simple
170 // --> use same calculation method
171
172 // Only compute full safety if massSafety
173 // startFullSafety = fPathFinder->Compute
174 }
175
176 // Is the particle charged or has it a magne
177 //
178 G4double particleCharge = pParticle->GetChar
179 G4double magneticMoment = pParticle->GetMagn
180 G4double restMass = pParticle->GetMass
181
182 fMassGeometryLimitedStep = false ; // Set d
183 fGeometryLimitedStep = false;
184
185 // There is no need to locate the current vo
186 // On track construction
187 // By the tracking, after all AlongStepDoI
188
189 // Check if the particle has a force, EM or
190 //
191 G4FieldManager* fieldMgr= nullptr;
192 G4bool fieldExertsForce = false ;
193
194 const G4Field* ptrField= nullptr;
195
196 fieldMgr = fFieldPropagator->FindAndSetField
197 G4bool eligibleEM = (particleCharge != 0.0)
198 || ( fUseMagneticMoment &&
199 G4bool eligibleGrav = fUseGravity && (restM
200
201 if( (fieldMgr!=nullptr) && (eligibleEM||elig
202 {
203 // Message the field Manager, to configur
204 //
205 fieldMgr->ConfigureForTrack( &track );
206
207 // The above call can transition from a n
208 // If the field manager has no field ptr,
209 // by definition ( = there is no field !
210 //
211 ptrField= fieldMgr->GetDetectorField();
212
213 if( ptrField != nullptr)
214 {
215 fieldExertsForce = eligibleEM
216 || ( eligibleGrav && ptrField->I
217 }
218 }
219 G4double momentumMagnitude = pParticle->GetT
220
221 if( fieldExertsForce )
222 {
223 auto equationOfMotion= fFieldPropagator->
224
225 G4ChargeState chargeState(particleCharge,
226 magneticMoment,
227 pParticleDef->G
228 if( equationOfMotion )
229 {
230 equationOfMotion->SetChargeMomentumMas
231
232
233 }
234 #ifdef G4DEBUG_TRANSPORT
235 else
236 {
237 G4cerr << " ERROR in G4CoupledTranspor
238 << "Cannot find valid Equation
239 << " Unable to pass Charge, Mom
240 }
241 #endif
242 }
243
244 G4ThreeVector polarizationVec = track.GetPo
245 G4FieldTrack aFieldTrack = G4FieldTrack(sta
246 tra
247 tra
248 tra
249 res
250 par
251 pol
252 pPa
253 0.0
254 pPa
255 G4int stepNo= track.GetCurrentStepNumber();
256
257 ELimited limitedStep;
258 G4FieldTrack endTrackState('a'); // Defaul
259
260 fMassGeometryLimitedStep = false ; // de
261 fGeometryLimitedStep = false;
262 if( currentMinimumStep > 0 )
263 {
264 G4double newMassSafety= 0.0; // tem
265
266 // Do the Transport in the field (non re
267 //
268 lengthAlongCurve = fPathFinder->ComputeS
269
270
271
272
273
274
275
276
277 G4double newFullSafety= fPathFinder->Get
278 // this was estimated already in step
279
280 if( limitedStep == kUnique || limitedSte
281 {
282 fMassGeometryLimitedStep = true ;
283 }
284
285 fGeometryLimitedStep = (fPathFinder->Get
286
287 #ifdef G4DEBUG_TRANSPORT
288 if( fMassGeometryLimitedStep && !fGeomet
289 {
290 std::ostringstream message;
291 message << " ERROR in determining geom
292 message << " Limiting: mass=" << fMa
293 << " any= " << fGeometryLimite
294 message << "Incompatible conditions -
295 G4Exception("G4CoupledTransportation::
296 "PathFinderConfused", Fata
297 }
298 #endif
299
300 geometryStepLength = std::min( lengthAlo
301
302 // Momentum: Magnitude and direction ca
303 //
304 fMomentumChanged = true ;
305 fTransportEndMomentumDir = endTrackState
306
307 // Remember last safety origin & value.
308 fPreviousSftOrigin = startPosition ;
309 fPreviousMassSafety = newMassSafety ;
310 fPreviousFullSafety = newFullSafety ;
311 // fpSafetyHelper->SetCurrentSafety( new
312
313 #ifdef G4DEBUG_TRANSPORT
314 if( verboseLevel > 1 )
315 {
316 G4cout << "G4Transport:CompStep> "
317 << " called the pathfinder for
318 << " and obtained step = " << l
319 G4cout << " New safety (preStep) = "
320 }
321 #endif
322
323 // Store as best estimate value
324 startFullSafety = newFullSafety ;
325
326 // Get the End-Position and End-Momentum
327 fTransportEndPosition = endTrackState.Ge
328 fTransportEndKineticEnergy = endTrackSt
329 }
330 else
331 {
332 geometryStepLength = lengthAlongCurve=
333 fMomentumChanged = false ;
334 // fMassGeometryLimitedStep = false ;
335 // fGeometryLimitedStep = true;
336 fTransportEndMomentumDir = track.GetMome
337 fTransportEndKineticEnergy = track.GetK
338
339 fTransportEndPosition = startPosition;
340
341 endTrackState= aFieldTrack; // Ensures
342 }
343 // G4FieldTrack aTrackState(endTrackState);
344
345 if( !fieldExertsForce )
346 {
347 fParticleIsLooping = false ;
348 fMomentumChanged = false ;
349 fEndGlobalTimeComputed = false ;
350 }
351 else
352 {
353 fParticleIsLooping = fFieldPropagator->I
354
355 #ifdef G4DEBUG_TRANSPORT
356 if( verboseLevel > 1 )
357 {
358 G4cout << " G4CT::CS End Position = "
359 << fTransportEndPosition << G4e
360 G4cout << " G4CT::CS End Direction = "
361 << fTransportEndMomentumDir <<
362 }
363 #endif
364 if( fFieldPropagator->GetCurrentFieldMan
365 {
366 // If the field can change energy, t
367 // - so this should have been upd
368 //
369 fCandidateEndGlobalTime = endTrack
370 fEndGlobalTimeComputed = true;
371
372 // was ( fCandidateEndGlobalTime !=
373 // a cleaner way is to have FieldTra
374 // is updated
375 }
376 else
377 {
378 // The energy should be unchanged by
379 // - so the time changed will be
380 //
381 fEndGlobalTimeComputed = false;
382
383 #ifdef G4VERBOSE
384 // Check that the integration preser
385 // - and if not correct this!
386 G4double startEnergy= track.GetKine
387 G4double endEnergy= fTransportEndKi
388
389 G4double absEdiff = std::fabs(startE
390 if( (verboseLevel > 1) && ( absEdiff
391 {
392 ReportInexactEnergy(startEnergy, e
393 } // end of if (verboseLevel)
394 #endif
395 // Correct the energy for fields tha
396 // This - hides the integration err
397 // - but gives a better physic
398 fTransportEndKineticEnergy= track.Ge
399 }
400 }
401
402 fEndPointDistance = (fTransportEndPosition
403 fTransportEndSpin = endTrackState.GetSpin();
404
405 // Calculate the safety
406
407 safetyProposal= startFullSafety; // used t
408 // Changed to accomodate processes that ca
409
410 // Update safety for the end-point, if becom
411
412 if( (startFullSafety < fEndPointDistance )
413 && ( particleCharge != 0.0 ) ) // Onl
414 // && !fGeometryLimitedStep ) // To-Try:
415 {
416 G4double endFullSafety =
417 fPathFinder->ComputeSafety( fTransport
418 // Expected mission -- only mass geome
419 // fLinearNavigator->ComputeSafety(
420 // Yet discrete processes only have po
421 // currently revise the safety
422 // ==> so we use the all-geometry sa
423
424 fpSafetyHelper->SetCurrentSafety( endFul
425 // Pushing safety to Helper avoids rec
426
427 G4ThreeVector centerPt= G4ThreeVector(0.
428 G4double endMassSafety= fPathFinder->Obt
429 // Retrieves the mass value from Path
430
431 fPreviousMassSafety = endMassSafety ;
432 fPreviousFullSafety = endFullSafety;
433 fPreviousSftOrigin = fTransportEndPositi
434
435 // The convention (Stepping Manager's) i
436 //
437 safetyProposal = endFullSafety + fEndPoi
438 // --> was endMassSafety
439 // Changed to accomodate processes that
440
441 #ifdef G4DEBUG_TRANSPORT
442 G4int prec= G4cout.precision(12) ;
443 G4cout << "***CoupledTransportation::Alo
444 G4cout << " Revised Safety at endpoint
445 << " give safety values: Mass=
446 << " All= " << endFullSafety <<
447 G4cout << " Adding endpoint distance "
448 << " to obtain pseudo-safety= "
449 G4cout.precision(prec);
450 }
451 else
452 {
453 G4int prec= G4cout.precision(12) ;
454 G4cout << "***CoupledTransportation::Alo
455 G4cout << " Quick Safety estimate at en
456 << fTransportEndPosition
457 << " gives safety endpoint valu
458 << startFullSafety - fEndPointDis
459 << " using start-point value " <
460 << " and endpointDistance " << f
461 G4cout.precision(prec);
462 #endif
463 }
464
465 proposedSafetyForStart= safetyProposal;
466 fParticleChange.ProposeTrueStepLength(geomet
467
468 return geometryStepLength ;
469 }
470
471 //////////////////////////////////////////////
472
473 void G4CoupledTransportation::
474 ReportMove( G4ThreeVector OldVector, G4ThreeVe
475 const G4String& Quantity )
476 {
477 G4ThreeVector moveVec = ( NewVector - OldV
478
479 G4cerr << G4endl
480 << "*******************************
481 << G4endl;
482 G4cerr << "Endpoint has moved between valu
483 << " and value from Track in PostSt
484 << "Change of " << Quantity << " is
485 << " mm long, "
486 << " and its vector is " << (1.0/mm
487 << "Endpoint of ComputeStep was " <
488 << " and current position to locate
489 }
490
491 //////////////////////////////////////////////
492
493 G4VParticleChange* G4CoupledTransportation::Po
494
495 {
496 G4TouchableHandle retCurrentTouchable ; //
497
498 // Initialize ParticleChange (by setting al
499 // to correspond
500 // fParticleChange.Initialize(track) ; // T
501
502 fParticleChange.ProposeTrackStatus(track.Get
503
504 if( fSignifyStepInAnyVolume )
505 {
506 fParticleChange.ProposeFirstStepInVolume(
507 }
508 else
509 {
510 fParticleChange.ProposeFirstStepInVolume(
511 }
512
513 // Check that the end position and direction
514 // since call to AlongStepDoIt
515
516 #ifdef G4DEBUG_TRANSPORT
517 if( ( verboseLevel > 0 )
518 && ((fTransportEndPosition - track.GetPos
519 {
520 ReportMove( track.GetPosition(), fTranspo
521 "End of Step Position" );
522 G4cerr << " Problem in G4CoupledTransport
523 }
524
525 // If the Step was determined by the volume
526 // The pathFinder will know that the geometr
527
528 if( verboseLevel > 0 )
529 {
530 G4cout << " Calling PathFinder::Locate()
531 << " G4CoupledTransportation::Post
532 G4cout << " fGeometryLimitedStep is " <
533 }
534 #endif
535
536 if(fGeometryLimitedStep)
537 {
538 fPathFinder->Locate( track.GetPosition(),
539 track.GetMomentumDire
540 true);
541
542 // fCurrentTouchable will now become the p
543 // and what was the previous will be freed
544 // (Needed because the preStepPoint can po
545
546 fCurrentTouchableHandle=
547 fPathFinder->CreateTouchableHandle( G4Tr
548
549 #ifdef G4DEBUG_TRANSPORT
550 if( verboseLevel > 1 )
551 {
552 G4VPhysicalVolume* vol= fCurrentTouchab
553 G4cout << "CHECK !!!!!!!!!!! fCurrentTo
554 << vol;
555 if( vol ) { G4cout << "Name=" << vol->G
556 G4cout << G4endl;
557 }
558 #endif
559
560 // Check whether the particle is out of th
561 // If so it has exited and must be killed.
562 //
563 if( fCurrentTouchableHandle->GetVolume() =
564 {
565 fParticleChange.ProposeTrackStatus( fSt
566 }
567 retCurrentTouchable = fCurrentTouchableHan
568 // fParticleChange.SetTouchableHandle( fCu
569 }
570 else // fGeometryLimitedStep is false
571 {
572 #ifdef G4DEBUG_TRANSPORT
573 if( verboseLevel > 1 )
574 {
575 G4cout << "G4CoupledTransportation::Post
576 << " fGeometryLimitedStep = " <<
577 << " must be false " << G4endl;
578 }
579 #endif
580 // This serves only to move each of the Na
581 //
582 // fLinearNavigator->LocateGlobalPointWith
583
584 fPathFinder->ReLocate( track.GetPosition()
585 // track.GetMomentu
586
587 // Keep the value of the track's current T
588 // and use it to overwrite the (unset) on
589 // Expect this must be fCurrentTouchable t
590 // - could it be different, eg at the st
591 //
592 retCurrentTouchable = track.GetTouchableHa
593 // fParticleChange.SetTouchableHandle( tra
594 } // endif ( fGeometryLimitedStep )
595
596 #ifdef G4DEBUG_NAVIGATION
597 G4cout << " CoupledTransport::AlongStep GPI
598 << " last-step: any= " << fGeometryL
599 << " mass= " << fMassGeometryLimitedS
600 #endif
601
602 if( fSignifyStepInAnyVolume )
603 fParticleChange.ProposeLastStepInVolume(fG
604 else
605 fParticleChange.ProposeLastStepInVolume(f
606
607 SetTouchableInformation(retCurrentTouchable)
608
609 return &fParticleChange ;
610 }
611
612 //////////////////////////////////////////////
613 // New method takes over the responsibility to
614 // G4CoupledTransportation object:
615 // - at the start of a new track, and
616 // - on the resumption of a suspended tra
617 //
618 void
619 G4CoupledTransportation::StartTracking(G4Track
620 {
621 G4Transportation::StartTracking(aTrack);
622
623 G4ThreeVector position = aTrack->GetPosition
624 G4ThreeVector direction = aTrack->GetMomentu
625
626 fPathFinder->PrepareNewTrack( position, dire
627 // This implies a call to fPathFinder->Locat
628
629 // reset safety value and center
630 //
631 fPreviousMassSafety = 0.0 ;
632 fPreviousFullSafety = 0.0 ;
633 fPreviousSftOrigin = G4ThreeVector(0.,0.,0.)
634 }
635
636 //////////////////////////////////////////////
637
638 void
639 G4CoupledTransportation::EndTracking()
640 {
641 G4TransportationManager::GetTransportationMa
642 fPathFinder->EndTrack();
643 // Resets TransportationManager to use ord
644 }
645
646 //////////////////////////////////////////////
647
648 void
649 G4CoupledTransportation::
650 ReportInexactEnergy(G4double startEnergy, G4do
651 {
652 static G4ThreadLocal G4int no_warnings= 0, w
653 moduloFactor= 10,
654
655 if( std::fabs(startEnergy- endEnergy) > perT
656 {
657 no_large_ediff ++;
658 if( (no_large_ediff% warnModulo) == 0 )
659 {
660 no_warnings++;
661 std::ostringstream message;
662 message << "Energy change in Step is abo
663 << G4endl
664 << " Relative change in 'track
665 << std::setw(15) << (endEnergy-s
666 << G4endl
667 << " Starting E= " << std::set
668 << " MeV " << G4endl
669 << " Ending E= " << std::set
670 << " MeV " << G4endl
671 << "Energy has been corrected --
672 << " field propagation parameter
673 if ( (verboseLevel > 2 ) || (no_warnings
674 || (no_large_ediff == warnModulo * mod
675 {
676 message << "These include EpsilonStepM
677 << G4endl
678 << "which determine fractional
679 << G4endl
680 << "Note also the influence of
681 << G4endl;
682 }
683 message << "Bad 'endpoint'. Energy chang
684 << G4endl
685 << "Has occurred already " << no
686 G4Exception("G4CoupledTransportation::Al
687 "EnergyChange", JustWarning,
688 if( no_large_ediff == warnModulo * modul
689 {
690 warnModulo *= moduloFactor;
691 }
692 }
693 }
694 }
695