| 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 /// \brief This class is a slightly modified v
28 /// initially written by John Apostolak
29 /// But it should use the exact same al
30 //
31 // Contact : Mathieu Karamitros (kara (AT) cen
32 //
33 // History :
34 // -----------
35 // ===========================================
36 // Modified:
37 // 28 Oct 2011, P.Gumpl./J.Ap: Detect gravi
38 // 20 Nov 2008, J.Apostolakis: Push safety
39 // 9 Nov 2007, J.Apostolakis: Flag for sho
40 // 19 Jan 2006, P.MoraDeFreitas: Fix for su
41 // 11 Aug 2004, M.Asai: Add G4VSensitiveDet
42 // 21 June 2003, J.Apostolakis: Calling fiel
43 // track, to enable it to
44 // 13 May 2003, J.Apostolakis: Zero field a
45 // account correclty in al
46 // 29 June 2001, J.Apostolakis, D.Cote-Ahern
47 // correction for spin tra
48 // 20 Febr 2001, J.Apostolakis: update for
49 // 22 Sept 2000, V.Grichine: update of K
50 // ------------------------------------------
51 // 10 Oct 2011, M.Karamitros: G4ITTranspo
52 // Created: 19 March 1997, J. Apostolakis
53 // ===========================================
54 //
55 // -------------------------------------------
56
57 #include "G4ITTransportation.hh"
58
59 #include "G4ChordFinder.hh"
60 #include "G4FieldManager.hh"
61 #include "G4FieldManagerStore.hh"
62 #include "G4IT.hh"
63 #include "G4ITNavigator.hh"
64 #include "G4ITSafetyHelper.hh"
65 #include "G4ITTransportationManager.hh"
66 #include "G4LowEnergyEmProcessSubType.hh"
67 #include "G4ParticleTable.hh"
68 #include "G4ProductionCutsTable.hh"
69 #include "G4PropagatorInField.hh"
70 #include "G4ReferenceCast.hh"
71 #include "G4SystemOfUnits.hh"
72 #include "G4TrackingInformation.hh"
73 #include "G4TransportationManager.hh"
74 #include "G4UnitsTable.hh"
75
76 #include <memory>
77
78 class G4VSensitiveDetector;
79
80 #ifndef PrepareState
81 # define PrepareState()
82 G4ITTransportationState* __state = this->G
83 #endif
84
85 #ifndef State
86 #define State(theXInfo) (__state->theXInfo)
87 #endif
88
89 //#define DEBUG_MEM
90
91 #ifdef DEBUG_MEM
92 #include "G4MemStat.hh"
93 using namespace G4MemStat;
94 using G4MemStat::MemStat;
95 #endif
96
97 //#define G4DEBUG_TRANSPORT 1
98
99 G4ITTransportation::G4ITTransportation(const G
100 G4VITProcess(aName, fTransportation),
101 fThreshold_Warning_Energy(100 * MeV),
102 fThreshold_Important_Energy(250 * MeV),
103
104 fUnimportant_Energy(1 * MeV), // Old defa
105 fVerboseLevel(verbose)
106 {
107 pParticleChange = &fParticleChange;
108 G4TransportationManager* transportMgr;
109 transportMgr = G4TransportationManager::GetT
110 G4ITTransportationManager* ITtransportMgr;
111 ITtransportMgr = G4ITTransportationManager::
112 fLinearNavigator = ITtransportMgr->GetNaviga
113 fFieldPropagator = transportMgr->GetPropagat
114 fpSafetyHelper = ITtransportMgr->GetSafetyHe
115
116 // Cannot determine whether a field exists h
117 // depend on the relative order of creating
118 // field and this process. That order is no
119 // Instead later the method DoesGlobalFieldE
120
121 enableAtRestDoIt = false;
122 enableAlongStepDoIt = true;
123 enablePostStepDoIt = true;
124 SetProcessSubType(fLowEnergyTransportation);
125 SetInstantiateProcessState(true);
126 G4VITProcess::SetInstantiateProcessState(fal
127 fInstantiateProcessState = true;
128
129 G4VITProcess::fpState = std::make_shared<G4I
130 /*
131 if(fTransportationState == 0)
132 {
133 G4cout << "KILL in G4ITTransportation::G4IT
134 abort();
135 }
136 */
137 }
138
139 G4ITTransportation::G4ITTransportation(const G
140 G4VITProcess(right)
141 {
142 // Copy attributes
143 fVerboseLevel = right.fVerboseLevel;
144 fThreshold_Warning_Energy = right.fThreshold
145 fThreshold_Important_Energy = right.fThresho
146 fThresholdTrials = right.fThresholdTrials;
147 fUnimportant_Energy = right.fUnimportant_Ene
148 fSumEnergyKilled = right.fSumEnergyKilled;
149 fMaxEnergyKilled = right.fMaxEnergyKilled;
150 fShortStepOptimisation = right.fShortStepOpt
151
152 // Setup Navigators
153 G4TransportationManager* transportMgr;
154 transportMgr = G4TransportationManager::GetT
155 G4ITTransportationManager* ITtransportMgr;
156 ITtransportMgr = G4ITTransportationManager::
157 fLinearNavigator = ITtransportMgr->GetNaviga
158 fFieldPropagator = transportMgr->GetPropagat
159 fpSafetyHelper = ITtransportMgr->GetSafetyHe
160
161 // Cannot determine whether a field exists h
162 // depend on the relative order of creating
163 // field and this process. That order is no
164 // Instead later the method DoesGlobalFieldE
165
166 enableAtRestDoIt = false;
167 enableAlongStepDoIt = true;
168 enablePostStepDoIt = true;
169
170 pParticleChange = &fParticleChange;
171 SetInstantiateProcessState(true);
172 G4VITProcess::SetInstantiateProcessState(fal
173 fInstantiateProcessState = right.fInstantiat
174 }
175
176 G4ITTransportation& G4ITTransportation::operat
177 {
178 // if (this == &right) return *this;
179 return *this;
180 }
181
182 //////////////////////////////////////////////
183 /// Process State
184 //////////////////////////////////////////////
185 G4ITTransportation::G4ITTransportationState::G
186 fCurrentTouchableHandle(nullptr)
187 {
188 fTransportEndPosition = G4ThreeVector(0, 0,
189 fTransportEndMomentumDir = G4ThreeVector(0,
190 fTransportEndKineticEnergy = -1;
191 fTransportEndSpin = G4ThreeVector(0, 0, 0);
192 fMomentumChanged = false;
193 fEnergyChanged = false;
194 fEndGlobalTimeComputed = false;
195 fCandidateEndGlobalTime = -1;
196 fParticleIsLooping = false;
197
198 static G4ThreadLocal G4TouchableHandle *null
199 if (nullTouchableHandle == nullptr) nullTouc
200 // Points to (G4VTouchable*) 0
201
202 fCurrentTouchableHandle = *nullTouchableHand
203 fGeometryLimitedStep = false;
204 fPreviousSftOrigin = G4ThreeVector(0, 0, 0);
205 fPreviousSafety = 0.0;
206 fNoLooperTrials = 0;
207 fEndPointDistance = -1;
208 }
209
210 G4ITTransportation::G4ITTransportationState::~
211 {
212 ;
213 }
214
215 G4ITTransportation::~G4ITTransportation()
216 {
217 #ifdef G4VERBOSE
218 if ((fVerboseLevel > 0) && (fSumEnergyKilled
219 {
220 G4cout << " G4ITTransportation: Statistics
221 << G4endl;
222 G4cout << " Sum of energy of loopers kil
223 << fSumEnergyKilled << G4endl;
224 G4cout << " Max energy of loopers killed
225 << fMaxEnergyKilled << G4endl;
226 }
227 #endif
228 }
229
230 void G4ITTransportation::BuildPhysicsTable(con
231 {
232 fpSafetyHelper->InitialiseHelper();
233 }
234
235 G4bool G4ITTransportation::DoesGlobalFieldExis
236 {
237 G4TransportationManager* transportMgr;
238 transportMgr = G4TransportationManager::GetT
239
240 // fFieldExists= transportMgr->GetFieldManag
241 // return fFieldExists;
242 return transportMgr->GetFieldManager()->Does
243 }
244
245 //////////////////////////////////////////////
246 //
247 // Responsibilities:
248 // Find whether the geometry limits the Ste
249 // Calculate the new value of the safety an
250 // Store the final time, position and momen
251 G4double
252 G4ITTransportation::
253 AlongStepGetPhysicalInteractionLength(const G4
254 G4double
255 G4double
256 G4double
257 G4GPILSe
258 {
259 PrepareState();
260 G4double geometryStepLength(-1.0), newSafety
261
262 State(fParticleIsLooping) = false;
263 State(fEndGlobalTimeComputed) = false;
264 State(fGeometryLimitedStep) = false;
265
266 // Initial actions moved to StartTrack()
267 // --------------------------------------
268 // Note: in case another process changes tou
269 // it will be necessary to add here (for
270 // State(fCurrentTouchableHandle) = track.Ge
271
272 // GPILSelection is set to defaule value of
273 // It is a return value
274 //
275 *selection = CandidateForSelection;
276
277 // Get initial Energy/Momentum of the track
278 //
279 const G4DynamicParticle* pParticle = track.G
280 // const G4ParticleDefinition* pParticleDef
281 G4ThreeVector startMomentumDir = pParticle->
282 G4ThreeVector startPosition = track.GetPosit
283
284 // G4double theTime = track.GetGlob
285
286 // The Step Point safety can be limited by o
287 // assumptions of any process - it's not alw
288 // We calculate the starting point's isotrop
289 //
290 G4ThreeVector OriginShift = startPosition -
291 G4double MagSqShift = OriginShift.mag2();
292 if (MagSqShift >= sqr(State(fPreviousSafety)
293 {
294 currentSafety = 0.0;
295 }
296 else
297 {
298 currentSafety = State(fPreviousSafety) - s
299 }
300
301 // Is the particle charged ?
302 //
303 G4double particleCharge = pParticle->GetChar
304
305 // There is no need to locate the current vo
306 // On track construction
307 // By the tracking, after all AlongStepDoI
308
309 // Check whether the particle have an (EM) f
310 //
311 G4FieldManager* fieldMgr = nullptr;
312 G4bool fieldExertsForce = false;
313 if ((particleCharge != 0.0))
314 {
315 fieldMgr = fFieldPropagator->FindAndSetFie
316 if (fieldMgr != nullptr)
317 {
318 // Message the field Manager, to configu
319 fieldMgr->ConfigureForTrack(&track);
320 // Moved here, in order to allow a trans
321 // from a zero-field status (with fie
322 // to a finite field status
323
324 // If the field manager has no field, th
325 fieldExertsForce = (fieldMgr->GetDetecto
326 }
327 }
328
329 // G4cout << " G4Transport: field exerts fo
330 // << " fieldMgr= " << fieldMgr << G4endl
331
332 // Choose the calculation of the transportat
333 //
334 if (!fieldExertsForce)
335 {
336 G4double linearStepLength;
337 if (fShortStepOptimisation && (currentMini
338 {
339 // The Step is guaranteed to be taken
340 //
341 geometryStepLength = currentMinimumStep;
342 State(fGeometryLimitedStep) = false;
343 }
344 else
345 {
346 // Find whether the straight path inter
347 //
348 // fLinearNavigator->SetNavigatorState(G
349 linearStepLength = fLinearNavigator->Com
350
351
352
353
354 // G4cout << "linearStepLength : " << G4
355 // << " | currentMinimumStep: " << cu
356 // << " | trackID: " << track.GetTrac
357
358 // Remember last safety origin & value.
359 //
360 State(fPreviousSftOrigin) = startPositio
361 State(fPreviousSafety) = newSafety;
362
363 G4TrackStateManager& trackStateMan = Get
364 ->GetTrackStateManager();
365 fpSafetyHelper->LoadTrackState(trackStat
366 // fpSafetyHelper->SetTrackState(state);
367 fpSafetyHelper->SetCurrentSafety(newSafe
368 State(f
369 fpSafetyHelper->ResetTrackState();
370
371 // The safety at the initial point has b
372 //
373 currentSafety = newSafety;
374
375 State(fGeometryLimitedStep) = (linearSte
376 if (State(fGeometryLimitedStep))
377 {
378 // The geometry limits the Step size (
379 geometryStepLength = linearStepLength;
380 }
381 else
382 {
383 // The full Step is taken.
384 geometryStepLength = currentMinimumSte
385 }
386 }
387 State(fEndPointDistance) = geometryStepLen
388
389 // Calculate final position
390 //
391 State(fTransportEndPosition) = startPositi
392 + geometryStepLength * startMomentumDi
393
394 // Momentum direction, energy and polarisa
395 //
396 State(fTransportEndMomentumDir) = startMom
397 State(fTransportEndKineticEnergy) = track.
398 State(fTransportEndSpin) = track.GetPolari
399 State(fParticleIsLooping) = false;
400 State(fMomentumChanged) = false;
401 State(fEndGlobalTimeComputed) = true;
402 State(theInteractionTimeLeft) = State(fEnd
403 / track.GetVelocity();
404 State(fCandidateEndGlobalTime) = State(the
405 + track.GetGlobalTime();
406 /*
407 G4cout << "track.GetVelocity() : "
408 << track.GetVelocity() << G4endl;
409 G4cout << "State(endpointDistance) : "
410 << G4BestUnit(State(endpointDistanc
411 G4cout << "State(theInteractionTimeLeft) :
412 << G4BestUnit(State(theInteractionT
413 G4cout << "track.GetGlobalTime() : "
414 << G4BestUnit(track.GetGlobalTime()
415 */
416 }
417 else // A field exerts force
418 {
419
420 G4ExceptionDescription exceptionDescriptio
421 exceptionDescription
422 << "ITTransportation does not support
423 exceptionDescription
424 << " If you are dealing with a tradiat
425 exceptionDescription << "please use G4Tran
426
427 G4Exception("G4ITTransportation::AlongStep
428 "NoExternalFieldSupport", Fata
429 /*
430 G4double momentumMagnitude = pParti
431 // G4ThreeVector EndUnitMomentum
432 G4double lengthAlongCurve ;
433 G4double restMass = pParticleDef->G
434
435 fFieldPropagator->SetChargeMomentumMass(
436 momentumMagnitude, // in Mev/c
437 restMass ) ;
438
439 G4ThreeVector spin = track.GetPola
440 G4FieldTrack aFieldTrack = G4FieldTrack(
441 track.GetMomentumDirection(),
442 0.0,
443 track.GetKineticEnergy(),
444 restMass,
445 track.GetVelocity(),
446 track.GetGlobalTime(), // Lab.
447 track.GetProperTime(), // Part.
448 &spin ) ;
449 if( currentMinimumStep > 0 )
450 {
451 // Do the Transport in the field (non rec
452 //
453 lengthAlongCurve = fFieldPropagator->Comp
454 currentMinimumStep,
455 currentSafety,
456 track.GetVolume() ) ;
457 State(fGeometryLimitedStep)= lengthAlongC
458 if( State(fGeometryLimitedStep) )
459 {
460 geometryStepLength = lengthAlongCurve ;
461 }
462 else
463 {
464 geometryStepLength = currentMinimumStep
465 }
466
467 // Remember last safety origin & value.
468 //
469 State(fPreviousSftOrigin) = startPosition
470 State(fPreviousSafety) = currentSafety
471 fpSafetyHelper->SetCurrentSafety( newSafe
472 }
473 else
474 {
475 geometryStepLength = lengthAlongCurve=
476 State(fGeometryLimitedStep) = false ;
477 }
478
479 // Get the End-Position and End-Momentum
480 //
481 State(fTransportEndPosition) = aFieldTrac
482
483 // Momentum: Magnitude and direction can
484 //
485 State(fMomentumChanged) = true ;
486 State(fTransportEndMomentumDir) = aFieldT
487
488 State(fTransportEndKineticEnergy) = aFie
489
490 if( fFieldPropagator->GetCurrentFieldMana
491 {
492 // If the field can change energy, then t
493 // - so this should have been updated
494 //
495 State(fCandidateEndGlobalTime) = aField
496 State(fEndGlobalTimeComputed) = true;
497
498 State(theInteractionTimeLeft) = State(fCa
499 track
500
501 // was ( State(fCandidateEndGlobalTime) !
502 // a cleaner way is to have FieldTrack kn
503 }
504 else
505 {
506 // The energy should be unchanged by fiel
507 // - so the time changed will be calcu
508 //
509 State(fEndGlobalTimeComputed) = false;
510
511 // Check that the integration preserved t
512 // - and if not correct this!
513 G4double startEnergy= track.GetKineticEn
514 G4double endEnergy= State(fTransportEndK
515
516 static G4int no_inexact_steps=0, no_large
517 G4double absEdiff = std::fabs(startEnergy
518 if( absEdiff > perMillion * endEnergy )
519 {
520 no_inexact_steps++;
521 // Possible statistics keeping here ...
522 }
523 #ifdef G4VERBOSE
524 if( fVerboseLevel > 1 )
525 {
526 if( std::fabs(startEnergy- endEnergy) > p
527 {
528 static G4int no_warnings= 0, warnModulo=1
529 no_large_ediff ++;
530 if( (no_large_ediff% warnModulo) == 0 )
531 {
532 no_warnings++;
533 G4cout << "WARNING - G4Transportation::Al
534 << " Energy change in Step is above 1^-
535 << " Relative change in 'tracking' step
536 << std::setw(15) << (endEnergy-startEnerg
537 << " Starting E= " << std::setw(12) <
538 << G4endl
539 << " Ending E= " << std::setw(12) <
540 << G4endl;
541 G4cout << " Energy has been corrected --
542 << " field propagation parameters for acc
543 if( (fVerboseLevel > 2 ) || (no_warnings<
544 (no_large_ediff == warnModulo * modul
545 {
546 G4cout << " These include EpsilonStepMax(
547 << " which determine fractional error per
548 << G4endl
549 << " Note also the influence of the permi
550 << G4endl;
551 }
552 G4cerr << "ERROR - G4Transportation::Alon
553 << " Bad 'endpoint'. Energy change
554 << " and corrected. "
555 << " Has occurred already "
556 << no_large_ediff << " times." << G4endl;
557 if( no_large_ediff == warnModulo * modulo
558 {
559 warnModulo *= moduloFactor;
560 }
561 }
562 }
563 } // end of if (fVerboseLevel)
564 #endif
565 // Correct the energy for fields that con
566 // This - hides the integration error
567 // - but gives a better physical an
568 State(fTransportEndKineticEnergy)= track.
569 }
570
571 State(fTransportEndSpin) = aFieldTrack.Ge
572 State(fParticleIsLooping) = fFieldPropaga
573 State(endpointDistance) = (State(fTrans
574 startPositio
575 // State(theInteractionTimeLeft) =
576 track.GetVelocity()/Sta
577 */
578 }
579
580 // If we are asked to go a step length of 0,
581 // then a boundary will also limit the step
582 //
583 if (currentMinimumStep == 0.0)
584 {
585 if (currentSafety == 0.0)
586 {
587 State(fGeometryLimitedStep) = true;
588 // G4cout << "!!!! Safety is NULL,
589 // G4cout << " Track position : " <
590 // << G4endl;
591 }
592 }
593
594 // Update the safety starting from the end-p
595 // if it will become negative at the end-poi
596 //
597 if (currentSafety < State(fEndPointDistance)
598 {
599 // if( particleCharge == 0.0 )
600 // G4cout << " Avoiding call to Compu
601
602 if (particleCharge != 0.0)
603 {
604
605 G4double endSafety = fLinearNavigator->C
606 State(fTransportEndPosition));
607 currentSafety = endSafety;
608 State(fPreviousSftOrigin) = State(fTrans
609 State(fPreviousSafety) = currentSafety;
610
611 /*
612 G4VTrackStateHandle state =
613 GetIT(track)->GetTrackingInfo()->GetTra
614 */
615 G4TrackStateManager& trackStateMan = Get
616 ->GetTrackStateManager();
617 fpSafetyHelper->LoadTrackState(trackStat
618 // fpSafetyHelper->SetTrackState(state);
619 fpSafetyHelper->SetCurrentSafety(current
620 State(f
621 fpSafetyHelper->ResetTrackState();
622
623 // Because the Stepping Manager assumes
624 // add the StepLength
625 //
626 currentSafety += State(fEndPointDistance
627
628 #ifdef G4DEBUG_TRANSPORT
629 G4cout.precision(12);
630 G4cout << "***G4Transportation::AlongSte
631 G4cout << " Called Navigator->ComputeSa
632 << State(fTransportEndPosition)
633 << " and it returned safety= " <<
634 G4cout << " Adding endpoint distance "
635 << " to obtain pseudo-safety= "
636 #endif
637 }
638 }
639
640 // fParticleChange.ProposeTrueStepLength(geo
641
642 // G4cout << "G4ITTransportation::AlongStepGe
643 // << G4BestUnit(geometryStepLength,"L
644
645 return geometryStepLength;
646 }
647
648 void G4ITTransportation::ComputeStep(const G4T
649 const G4S
650 const dou
651 double& o
652 {
653 PrepareState();
654 const G4DynamicParticle* pParticle = track.G
655 G4ThreeVector startMomentumDir = pParticle->
656 G4ThreeVector startPosition = track.GetPosit
657
658 track.CalculateVelocity();
659 G4double initialVelocity = track.GetVelocity
660
661 State(fGeometryLimitedStep) = false;
662
663 /////////////////////////
664 // !!! CASE NO FIELD !!!
665 /////////////////////////
666 State(fCandidateEndGlobalTime) = timeStep +
667 State(fEndGlobalTimeComputed) = true;
668
669 // Choose the calculation of the transportat
670 //
671 if (!State(fMomentumChanged))
672 {
673 // G4cout << "Momentum has not chan
674 fParticleChange.ProposeVelocity(initialVel
675 oPhysicalStep = initialVelocity * timeStep
676
677 // Calculate final position
678 //
679 State(fTransportEndPosition) = startPositi
680 + oPhysicalStep * startMomentumDir;
681 }
682 }
683
684 //////////////////////////////////////////////
685 //
686 // Initialize ParticleChange (by setting al
687 // to correspond
688 #include "G4ParticleTable.hh"
689 G4VParticleChange* G4ITTransportation::AlongSt
690
691 {
692
693 #if defined (DEBUG_MEM)
694 MemStat mem_first, mem_second, mem_diff;
695 #endif
696
697 #if defined (DEBUG_MEM)
698 mem_first = MemoryUsage();
699 #endif
700
701 PrepareState();
702
703 // G4cout << "G4ITTransportation::AlongStepD
704 // set pdefOpticalPhoton
705 // Andrea Dotti: the following statement sho
706 // G4-MT transformation tools get confused i
707 // If needed contact: adotti@slac.stanford.e
708 static G4ThreadLocal G4ParticleDefinition* p
709 if (pdefOpticalPhoton == nullptr) pdefOptica
710 G4ParticleTable::GetParticleTable()->Fin
711
712 static G4ThreadLocal G4int noCalls = 0;
713 noCalls++;
714
715 fParticleChange.Initialize(track);
716
717 // Code for specific process
718 //
719 fParticleChange.ProposePosition(State(fTrans
720 fParticleChange.ProposeMomentumDirection(Sta
721 fParticleChange.ProposeEnergy(State(fTranspo
722 fParticleChange.SetMomentumChanged(State(fMo
723
724 fParticleChange.ProposePolarization(State(fT
725
726 G4double deltaTime = 0.0;
727
728 // Calculate Lab Time of Flight (ONLY if fi
729 // G4double endTime = State(fCandidateEndG
730 // G4double delta_time = endTime - startTime
731
732 G4double startTime = track.GetGlobalTime();
733 ///_________________________________________
734 /// !!!!!!!
735 /// A REVOIR !!!!
736 if (State(fEndGlobalTimeComputed) == 0)
737 {
738 // The time was not integrated .. make the
739 //
740 G4double initialVelocity = stepData.GetPre
741 G4double stepLength = track.GetStepLength(
742
743 deltaTime = 0.0; // in case initialVelocit
744 if (track.GetParticleDefinition() == pdefO
745 {
746 // For only Optical Photon, final veloci
747 double finalVelocity = track.CalculateVe
748 fParticleChange.ProposeVelocity(finalVel
749 deltaTime = stepLength / finalVelocity;
750 }
751 else if (initialVelocity > 0.0)
752 {
753 deltaTime = stepLength / initialVelocity
754 }
755
756 State(fCandidateEndGlobalTime) = startTime
757 }
758 else
759 {
760 deltaTime = State(fCandidateEndGlobalTime)
761 }
762
763 fParticleChange.ProposeGlobalTime(State(fCan
764 fParticleChange.ProposeLocalTime(track.GetLo
765 /*
766 // Now Correct by Lorentz factor to get del
767
768 G4double restMass = track.GetDynamic
769 G4double deltaProperTime = deltaTime*( rest
770
771 fParticleChange.ProposeProperTime(track.Get
772 */
773
774 fParticleChange.ProposeTrueStepLength(track.
775
776 ///_________________________________________
777 ///
778
779 // If the particle is caught looping or is s
780 // boundaries) in a magnetic field (doing ma
781 // THEN this kills it ...
782 //
783 if (State(fParticleIsLooping))
784 {
785 G4double endEnergy = State(fTransportEndKi
786
787 if ((endEnergy < fThreshold_Important_Ener
788 >= fThresholdTrials))
789 {
790 // Kill the looping particle
791 //
792 // G4cout << "G4ITTransportation will ki
793 fParticleChange.ProposeTrackStatus(fStop
794
795 // 'Bare' statistics
796 fSumEnergyKilled += endEnergy;
797 if (endEnergy > fMaxEnergyKilled)
798 {
799 fMaxEnergyKilled = endEnergy;
800 }
801
802 #ifdef G4VERBOSE
803 if ((fVerboseLevel > 1) || (endEnergy >
804 {
805 G4cout
806 << " G4ITTransportation is killing
807 << G4endl<< " This track has " <
808 << " MeV energy." << G4endl;
809 G4cout << " Number of trials = " <<
810 << " No of calls to AlongStepDoIt =
811 << G4endl;
812 }
813 #endif
814 State(fNoLooperTrials) = 0;
815 }
816 else
817 {
818 State(fNoLooperTrials)++;
819 #ifdef G4VERBOSE
820 if ((fVerboseLevel > 2))
821 {
822 G4cout << " G4ITTransportation::Alon
823 << " Number of trials = " <<
824 << " No of calls to = " << n
825 }
826 #endif
827 }
828 }
829 else
830 {
831 State(fNoLooperTrials)=0;
832 }
833
834 // Another (sometimes better way) is to use
835 // to alleviate this problem ..
836
837 // Introduce smooth curved trajectories to p
838 //
839 fParticleChange.SetPointerToVectorOfAuxiliar
840 fFieldPropagator->GimmeTrajectoryVectorA
841
842 #if defined (DEBUG_MEM)
843 mem_second = MemoryUsage();
844 mem_diff = mem_second-mem_first;
845 G4cout << "\t || MEM || End of G4ITTransport
846 << mem_diff << G4endl;
847 #endif
848
849 return &fParticleChange;
850 }
851
852 //////////////////////////////////////////////
853 //
854 // This ensures that the PostStep action is a
855 // so that it can do the relocation if it is
856 //
857
858 G4double
859 G4ITTransportation::
860 PostStepGetPhysicalInteractionLength(const G4T
861 G4double,
862 G4ForceCo
863 {
864 *pForceCond = Forced;
865 return DBL_MAX; // was kInfinity ; but conve
866 }
867
868 //////////////////////////////////////////////
869 //
870
871 G4VParticleChange* G4ITTransportation::PostSte
872
873 {
874 // G4cout << "G4ITTransportation::PostSte
875
876 PrepareState();
877 G4TouchableHandle retCurrentTouchable; // Th
878 G4bool isLastStep = false;
879
880 // Initialize ParticleChange (by setting al
881 // to correspond
882 fParticleChange.Initialize(track); // To ini
883
884 fParticleChange.ProposeTrackStatus(track.Get
885
886 // If the Step was determined by the volume
887 // logically relocate the particle
888
889 if (State(fGeometryLimitedStep))
890 {
891
892 if(fVerboseLevel != 0)
893 {
894 G4cout << "Step is limited by geometry "
895 << "track ID : " << track.GetTrac
896 }
897
898 // fCurrentTouchable will now become the p
899 // and what was the previous will be freed
900 // (Needed because the preStepPoint can po
901
902 if ( State(fCurrentTouchableHandle)->GetVo
903 {
904 G4ExceptionDescription exceptionDescript
905 exceptionDescription << "No current touc
906 G4Exception(" G4ITTransportation::PostSt
907 FatalErrorInArgument, except
908 }
909
910 fLinearNavigator->SetGeometricallyLimitedS
911 fLinearNavigator->LocateGlobalPointAndUpda
912 track.GetPosition(), track.GetMomentum
913 State(fCurrentTouchableHandle), true);
914 // Check whether the particle is out of th
915 // If so it has exited and must be killed.
916 //
917 if ( State(fCurrentTouchableHandle)->GetVo
918 {
919 // abort();
920 #ifdef G4VERBOSE
921 if (fVerboseLevel > 0)
922 {
923 G4cout << "Track position : " << track
924 << " [nm]" << " Track ID : " <<
925 G4cout << "G4ITTransportation will kil
926 "State(fCurrentTouchableHandle)->G
927 }
928 #endif
929 fParticleChange.ProposeTrackStatus( fSto
930 }
931
932 retCurrentTouchable = State(fCurrentToucha
933
934 // G4cout << "Current volume : " << tra
935 // << " Next volume : "
936 // << (State(fCurrentTouchableHa
937 // State(fCurrentTouchableHandle)->GetVolum
938 // << " Position : " << track.Ge
939 // << " track ID : " << track.Ge
940 // << G4endl;
941
942 fParticleChange.SetTouchableHandle(State(f
943
944 // Update the Step flag which identifies t
945 isLastStep = fLinearNavigator->ExitedMothe
946 || fLinearNavigator->EnteredDau
947
948 #ifdef G4DEBUG_TRANSPORT
949 // Checking first implementation of flagg
950 G4bool exiting = fLinearNavigator->ExitedM
951 G4bool entering = fLinearNavigator->Entere
952
953 if( ! (exiting || entering) )
954 {
955 G4cout << " Transport> : Proposed isLas
956 << " Exiting " << fLinearNavigator->Exit
957 << " Entering " << fLinearNavigator->Ent
958 << " Track position : " << track.GetPosi
959 << G4endl;
960 G4cout << " Track position : " << track.
961 << G4endl;
962 }
963 #endif
964 }
965 else // fGeometryLimitedStep is false
966 {
967 // This serves only to move the Navigator'
968 //
969 // abort();
970 fLinearNavigator->LocateGlobalPointWithinV
971
972 // The value of the track's current Toucha
973 // (and it must be correct because we must
974 // overwrite the (unset) one in particle c
975 // It must be fCurrentTouchable too ??
976 //
977 fParticleChange.SetTouchableHandle(track.G
978 retCurrentTouchable = track.GetTouchableHa
979
980 isLastStep = false;
981 #ifdef G4DEBUG_TRANSPORT
982 // Checking first implementation of flagg
983 G4cout << " Transport> Proposed isLastStep
984 << " Geometry did not limit step. Position
985 << track.GetPosition()/ nanometer << G4end
986 #endif
987 } // endif ( fGeometryLimitedStep )
988
989 fParticleChange.ProposeLastStepInVolume(isLa
990
991 const G4VPhysicalVolume* pNewVol = retCurren
992 const G4Material* pNewMaterial = nullptr;
993 G4VSensitiveDetector* pNewSensitiveDetector
994
995 if (pNewVol != nullptr)
996 {
997 pNewMaterial = pNewVol->GetLogicalVolume()
998 pNewSensitiveDetector = pNewVol->GetLogica
999 }
1000
1001 // ( <const_cast> pNewMaterial ) ;
1002
1003 fParticleChange.SetMaterialInTouchable((G4M
1004 fParticleChange.SetSensitiveDetectorInTouch
1005
1006 const G4MaterialCutsCouple* pNewMaterialCut
1007 if (pNewVol != nullptr)
1008 {
1009 pNewMaterialCutsCouple =
1010 pNewVol->GetLogicalVolume()->GetMater
1011 }
1012
1013 if (pNewVol != nullptr && pNewMaterialCutsC
1014 && pNewMaterialCutsCouple->GetMaterial(
1015 {
1016 // for parametrized volume
1017 //
1018 pNewMaterialCutsCouple = G4ProductionCuts
1019 ->GetMaterialCutsCouple(pNewMaterial,
1020 pNewMaterialC
1021 }
1022 fParticleChange.SetMaterialCutsCoupleInTouc
1023
1024 // temporarily until Get/Set Material of Pa
1025 // and StepPoint can be made const.
1026 // Set the touchable in ParticleChange
1027 // this must always be done because the par
1028 // uses this value to overwrite the current
1029 //
1030 fParticleChange.SetTouchableHandle(retCurre
1031
1032 return &fParticleChange;
1033 }
1034
1035 // New method takes over the responsibility t
1036 // G4Transportation object at the start of a
1037 // a suspended track.
1038
1039 void G4ITTransportation::StartTracking(G4Trac
1040 {
1041 G4VProcess::StartTracking(track);
1042 if (fInstantiateProcessState)
1043 {
1044 // G4VITProcess::fpState = new G4ITTra
1045 G4VITProcess::fpState = std::make_shared<
1046 // Will set in the same time fTransportat
1047 }
1048
1049 fpSafetyHelper->NewTrackState();
1050 fpSafetyHelper->SaveTrackState(
1051 GetIT(track)->GetTrackingInfo()->GetTra
1052
1053 // The actions here are those that were tak
1054 // when track.GetCurrentStepNumber()==1
1055
1056 // reset safety value and center
1057 //
1058 // State(fPreviousSafety) = 0.0 ;
1059 // State(fPreviousSftOrigin) = G4ThreeVe
1060
1061 // reset looping counter -- for motion in f
1062 // State(fNoLooperTrials)= 0;
1063 // Must clear this state .. else it depends
1064 // --> a better solution would set this fr
1065 // Was if( aTrack->GetCurrentStepNumber()==
1066
1067 // ChordFinder reset internal state
1068 //
1069 if (DoesGlobalFieldExist())
1070 {
1071 fFieldPropagator->ClearPropagatorState();
1072 // Resets all state of field propagator c
1073 // including safety values (in case of ov
1074
1075 // G4ChordFinder* chordF= fFieldPropagato
1076 // if( chordF ) chordF->ResetStepEstimate
1077 }
1078
1079 // Make sure to clear the chord finders of
1080 static G4ThreadLocal G4FieldManagerStore* f
1081 if (fieldMgrStore == nullptr) fieldMgrStore
1082 fieldMgrStore->ClearAllChordFindersState();
1083
1084 // Update the current touchable handle (fr
1085 //
1086 PrepareState();
1087 State(fCurrentTouchableHandle) = track->Get
1088
1089 G4VITProcess::StartTracking(track);
1090 }
1091
1092 #undef State
1093 #undef PrepareState
1094