Belle II Software development
G4LongLivedNeutralTransportation.cc
1/**************************************************************************
2 * basf2 (Belle II Analysis Software Framework) *
3 * Author: The Belle II Collaboration *
4 * *
5 * See git log for contributors and copyright holders. *
6 * This file is licensed under LGPL-3.0, see LICENSE.md. *
7 **************************************************************************/
8
9// modified from GEANT4 G4Transportation class
10
11#include <framework/logging/Logger.h>
12
13#include "G4TransportationProcessType.hh"
14
15#include "G4PhysicalConstants.hh"
16#include "G4SystemOfUnits.hh"
17#include "G4ProductionCutsTable.hh"
18#include "G4PrimaryParticle.hh"
19#include "G4FieldManager.hh"
20#include "G4Navigator.hh"
21#include "G4Track.hh"
22#include "G4Step.hh"
23
24#include "G4FieldManagerStore.hh"
25#include <simulation/longlivedneutral/G4LongLivedNeutralTransportation.h>
26
27class G4VSensitiveDetector;
28
29using namespace Belle2;
30
32
34//
35// Constructor
36
38 : G4VProcess(G4String("LongLivedNeutralTransportation"), fTransportation),
39 fFieldExertedForce(false),
40 fPreviousSftOrigin(0., 0., 0.),
41 fPreviousSafety(0.0),
43 fShortStepOptimisation(false) // Old default: true (=fast short steps)
44{
45 SetProcessSubType(static_cast<G4int>(TRANSPORTATION));
46 pParticleChange = &fParticleChange; // Required to conform to G4VProcess
47 SetVerboseLevel(verbosity);
48
49 G4TransportationManager* transportMgr ;
50
51 transportMgr = G4TransportationManager::GetTransportationManager() ;
52
53 fLinearNavigator = transportMgr->GetNavigatorForTracking() ;
54
55 fFieldPropagator = transportMgr->GetPropagatorInField() ;
56
57
58
59
61 // Use the old defaults: Warning = 100 MeV, Important = 250 MeV, No Trials = 10;
62
63 // Cannot determine whether a field exists here, as it would
64 // depend on the relative order of creating the detector's
65 // field and this process. That order is not guaranteed.
67 // This value must be updated using DoesAnyFieldExist() at least at the
68 // start of each Run -- for now this is at the Start of every Track. TODO
69
70 static G4ThreadLocal G4TouchableHandle* pNullTouchableHandle = 0;
71 pNullTouchableHandle = new G4TouchableHandle;
72 fCurrentTouchableHandle = *pNullTouchableHandle;
73 // Points to (G4VTouchable*) 0
74
75
76#ifdef G4VERBOSE
77 if (verboseLevel > 0) {
78 G4cout << " G4LongLivedNeutralTransportation constructor> set fShortStepOptimisation to ";
79 if (fShortStepOptimisation) { G4cout << "true" << G4endl; }
80 else { G4cout << "false" << G4endl; }
81 }
82#endif
83}
84
86
94
96
97void
99{
100 outStr << " G4LongLivedNeutralTransportation: Statistics for looping particles " << G4endl;
101 outStr << " No looping tracks found or killed. " << G4endl;
102
103}
104
106//
107// Responsibilities:
108// Find whether the geometry limits the Step, and to what length
109// Calculate the new value of the safety and return it.
110// Store the final time, position and momentum.
111
113AlongStepGetPhysicalInteractionLength(const G4Track& track,
114 G4double, // previousStepSize
115 G4double currentMinimumStep,
116 G4double& currentSafety,
117 G4GPILSelection* selection)
118{
119 G4double geometryStepLength = -1.0, newSafety = -1.0;
120
121 // Initial actions moved to StartTrack()
122 // --------------------------------------
123 // Note: in case another process changes touchable handle
124 // it will be necessary to add here (for all steps)
125 // fCurrentTouchableHandle = aTrack->GetTouchableHandle();
126
127 // GPILSelection is set to default value of CandidateForSelection
128 // It is a return value
129 //
130 *selection = CandidateForSelection ;
131
133 fLastStepInVolume = false;
134 fNewTrack = false;
135
136 fParticleChange.ProposeFirstStepInVolume(fFirstStepInVolume);
137
138 // Get initial Energy/Momentum of the track
139 //
140 const G4DynamicParticle* pParticle = track.GetDynamicParticle() ;
141 const G4PrimaryParticle* primaryParticle = pParticle->GetPrimaryParticle();
142 G4ThreeVector startMomentumDir = primaryParticle->GetMomentumDirection() ;
143 G4ThreeVector startPosition = track.GetPosition() ;
144
145 // The Step Point safety can be limited by other geometries and/or the
146 // assumptions of any process - it's not always the geometrical safety.
147 // We calculate the starting point's isotropic safety here.
148 //
149 G4ThreeVector OriginShift = startPosition - fPreviousSftOrigin ;
150 G4double MagSqShift = OriginShift.mag2() ;
151 if (MagSqShift >= sqr(fPreviousSafety)) {
152 currentSafety = 0.0 ;
153 } else {
154 currentSafety = fPreviousSafety - std::sqrt(MagSqShift) ;
155 }
156
157 // Is the particle charged or has it a magnetic moment?
158 //
159 G4double particleCharge = 0 ;
160
161 fGeometryLimitedStep = false ;
162
163 // There is no need to locate the current volume. It is Done elsewhere:
164 // On track construction
165 // By the tracking, after all AlongStepDoIts, in "Relocation"
166
167 G4double linearStepLength ;
168 if (fShortStepOptimisation && (currentMinimumStep <= currentSafety)) {
169 // The Step is guaranteed to be taken
170 //
171 geometryStepLength = currentMinimumStep ;
172 fGeometryLimitedStep = false ;
173 } else {
174 // Find whether the straight path intersects a volume
175 //
176 linearStepLength = fLinearNavigator->ComputeStep(startPosition,
177 startMomentumDir,
178 currentMinimumStep,
179 newSafety) ;
180 // Remember last safety origin & value.
181 //
182 fPreviousSftOrigin = startPosition ;
183 fPreviousSafety = newSafety ;
184
185
186 currentSafety = newSafety ;
187
188 fGeometryLimitedStep = (linearStepLength <= currentMinimumStep);
190 // The geometry limits the Step size (an intersection was found.)
191 geometryStepLength = linearStepLength ;
192 } else {
193 // The full Step is taken.
194 geometryStepLength = currentMinimumStep ;
195 }
196
197 fEndPointDistance = geometryStepLength ;
198
199 // Calculate final position
200 //
201 fTransportEndPosition = startPosition + geometryStepLength * startMomentumDir;
202
203 // Momentum direction, energy and polarisation are unchanged by transport
204 //
205 fTransportEndMomentumDir = startMomentumDir ;
206 fTransportEndKineticEnergy = primaryParticle->GetKineticEnergy() ;
207 fTransportEndSpin = track.GetPolarization();
208 fMomentumChanged = false ;
209 fEndGlobalTimeComputed = false ;
210 }
211
212
213 // If we are asked to go a step length of 0, and we are on a boundary
214 // then a boundary will also limit the step -> we must flag this.
215 //
216 if (currentMinimumStep == 0.0) {
217 if (currentSafety == 0.0) { fGeometryLimitedStep = true; }
218 }
219
220 // Update the safety starting from the end-point,
221 // if it will become negative at the end-point.
222 //
223 if (currentSafety < fEndPointDistance) {
224 if (particleCharge != 0.0) {
225 G4double endSafety =
227 currentSafety = endSafety ;
229 fPreviousSafety = currentSafety ;
230
231
232 // Because the Stepping Manager assumes it is from the start point,
233 // add the StepLength
234 //
235 currentSafety += fEndPointDistance ;
236
237#ifdef G4DEBUG_TRANSPORT
238 G4cout.precision(12) ;
239 G4cout << "***G4LongLivedNeutralTransportation::AlongStepGPIL ** " << G4endl ;
240 G4cout << " Called Navigator->ComputeSafety at " << fTransportEndPosition
241 << " and it returned safety= " << endSafety << G4endl ;
242 G4cout << " Adding endpoint distance " << fEndPointDistance
243 << " to obtain pseudo-safety= " << currentSafety << G4endl ;
244 } else {
245 G4cout << "***G4LongLivedNeutralTransportation::AlongStepGPIL ** " << G4endl ;
246 G4cout << " Avoiding call to ComputeSafety : " << G4endl;
247 G4cout << " charge = " << particleCharge << G4endl;
248 G4cout << " mag moment = " << magneticMoment << G4endl;
249#endif
250 }
251 }
252
253 fParticleChange.ProposeTrueStepLength(geometryStepLength) ;
254
255 return geometryStepLength ;
256}
257
259//
260// Initialize ParticleChange (by setting all its members equal
261// to corresponding members in G4Track)
262
263G4VParticleChange* G4LongLivedNeutralTransportation::AlongStepDoIt(const G4Track& track,
264 const G4Step& stepData)
265{
266 // static G4ThreadLocal G4long noCallsASDI = 0;
267 // noCallsASDI++;
268
269 fParticleChange.Initialize(track) ;
270
271 // Code for specific process
272 //
273 fParticleChange.ProposePosition(fTransportEndPosition) ;
274 fParticleChange.ProposeMomentumDirection(fTransportEndMomentumDir) ;
276 fParticleChange.SetMomentumChanged(fMomentumChanged) ;
277
278 fParticleChange.ProposePolarization(fTransportEndSpin);
279
280 G4double deltaTime = 0.0 ;
281
282 // Calculate Lab Time of Flight (ONLY if field Equations used it!)
283 // G4double endTime = fCandidateEndGlobalTime;
284 // G4double delta_time = endTime - startTime;
285
286 G4double startTime = track.GetGlobalTime() ;
287
288 const G4DynamicParticle* pParticle = track.GetDynamicParticle() ;
289 const G4PrimaryParticle* primaryParticle = pParticle->GetPrimaryParticle();
290
291
293 // The time was not integrated .. make the best estimate possible
294 //
295
296 // use c*beta from primary particle, llp will not be slowed down anyways
297 // default was: stepData.GetPreStepPoint()->GetVelocity();, but assumes wrong values from particle definition
298 B2DEBUG(0, "Velocity calculated from Step data: " << stepData.GetPreStepPoint()->GetVelocity());
299
300 G4double initialVelocity = c_light * primaryParticle->GetTotalMomentum() / primaryParticle->GetTotalEnergy();
301 G4double stepLength = track.GetStepLength();
302
303 deltaTime = 0.0; // in case initialVelocity = 0
304 if (initialVelocity > 0.0) { deltaTime = stepLength / initialVelocity; }
305
306 fCandidateEndGlobalTime = startTime + deltaTime ;
307 fParticleChange.ProposeLocalTime(track.GetLocalTime() + deltaTime) ;
308 } else {
309 deltaTime = fCandidateEndGlobalTime - startTime ;
310 fParticleChange.ProposeGlobalTime(fCandidateEndGlobalTime) ;
311 }
312
313 // Now Correct by Lorentz factor to get delta "proper" Time
314 G4double restMass = primaryParticle->GetMass() ;
315 G4double deltaProperTime = deltaTime * (restMass / primaryParticle->GetTotalEnergy()) ;
316
317 fParticleChange.ProposeProperTime(track.GetProperTime() + deltaProperTime) ;
318 //fParticleChange.ProposeTrueStepLength( track.GetStepLength() ) ;
319
320 // Another (sometimes better way) is to use a user-limit maximum Step size
321 // to alleviate this problem ..
322
323 // Introduce smooth curved trajectories to particle-change
324 //
325 fParticleChange.SetPointerToVectorOfAuxiliaryPoints
326 (fFieldPropagator->GimmeTrajectoryVectorAndForgetIt());
327
328 return &fParticleChange ;
329}
330
332//
333// This ensures that the PostStep action is always called,
334// so that it can do the relocation if it is needed.
335//
336
339 G4double, // previousStepSize
340 G4ForceCondition* pForceCond)
341{
342 fFieldExertedForce = false; // Not known
343 *pForceCond = Forced ;
344 return DBL_MAX ; // was kInfinity ; but convention now is DBL_MAX
345}
346
348//
349
350G4VParticleChange* G4LongLivedNeutralTransportation::PostStepDoIt(const G4Track& track,
351 const G4Step&)
352{
353 G4TouchableHandle retCurrentTouchable ; // The one to return
354 G4bool isLastStep = false;
355
356 // Initialize ParticleChange (by setting all its members equal
357 // to corresponding members in G4Track)
358 // fParticleChange.Initialize(track) ; // To initialise TouchableChange
359
360 fParticleChange.ProposeTrackStatus(track.GetTrackStatus()) ;
361
362 // If the Step was determined by the volume boundary,
363 // logically relocate the particle
364
366 // fCurrentTouchable will now become the previous touchable,
367 // and what was the previous will be freed.
368 // (Needed because the preStepPoint can point to the previous touchable)
369
370 fLinearNavigator->SetGeometricallyLimitedStep() ;
372 LocateGlobalPointAndUpdateTouchableHandle(track.GetPosition(),
373 track.GetMomentumDirection(),
375 true) ;
376 // Check whether the particle is out of the world volume
377 // If so it has exited and must be killed.
378 //
379 if (fCurrentTouchableHandle->GetVolume() == 0) {
380 fParticleChange.ProposeTrackStatus(fStopAndKill) ;
381 }
382 retCurrentTouchable = fCurrentTouchableHandle ;
383 fParticleChange.SetTouchableHandle(fCurrentTouchableHandle) ;
384
385 // Update the Step flag which identifies the Last Step in a volume
387 isLastStep = fLinearNavigator->ExitedMotherVolume()
388 or fLinearNavigator->EnteredDaughterVolume() ;
389 else
390 isLastStep = fFieldPropagator->IsLastStepInVolume();
391 } else { // fGeometryLimitedStep is false
392 // This serves only to move the Navigator's location
393 //
394 fLinearNavigator->LocateGlobalPointWithinVolume(track.GetPosition()) ;
395
396 // The value of the track's current Touchable is retained.
397 // (and it must be correct because we must use it below to
398 // overwrite the (unset) one in particle change)
399 // It must be fCurrentTouchable too ??
400 //
401 fParticleChange.SetTouchableHandle(track.GetTouchableHandle()) ;
402 retCurrentTouchable = track.GetTouchableHandle() ;
403
404 isLastStep = false;
405 } // endif ( fGeometryLimitedStep )
406 fLastStepInVolume = isLastStep;
407
408 fParticleChange.ProposeFirstStepInVolume(fFirstStepInVolume);
409 fParticleChange.ProposeLastStepInVolume(isLastStep);
410
411 const G4VPhysicalVolume* pNewVol = retCurrentTouchable->GetVolume() ;
412 G4Material* pNewMaterial = nullptr ;
413 G4VSensitiveDetector* pNewSensitiveDetector = nullptr ;
414
415 if (pNewVol != 0) {
416 pNewMaterial = pNewVol->GetLogicalVolume()->GetMaterial();
417 pNewSensitiveDetector = pNewVol->GetLogicalVolume()->GetSensitiveDetector();
418 }
419
420 fParticleChange.SetMaterialInTouchable(pNewMaterial) ;
421 fParticleChange.SetSensitiveDetectorInTouchable(pNewSensitiveDetector) ;
422
423 const G4MaterialCutsCouple* pNewMaterialCutsCouple = 0;
424 if (pNewVol != 0) {
425 pNewMaterialCutsCouple = pNewVol->GetLogicalVolume()->GetMaterialCutsCouple();
426 }
427
428 if (pNewVol != 0 && pNewMaterialCutsCouple != 0
429 && pNewMaterialCutsCouple->GetMaterial() != pNewMaterial) {
430 // for parametrized volume
431 //
432 pNewMaterialCutsCouple =
433 G4ProductionCutsTable::GetProductionCutsTable()
434 ->GetMaterialCutsCouple(pNewMaterial,
435 pNewMaterialCutsCouple->GetProductionCuts());
436 }
437 fParticleChange.SetMaterialCutsCoupleInTouchable(pNewMaterialCutsCouple);
438
439 // temporarily until Get/Set Material of ParticleChange,
440 // and StepPoint can be made const.
441 // Set the touchable in ParticleChange
442 // this must always be done because the particle change always
443 // uses this value to overwrite the current touchable pointer.
444 //
445 fParticleChange.SetTouchableHandle(retCurrentTouchable) ;
446
447 return &fParticleChange ;
448}
449
451// New method takes over the responsibility to reset the state of
452// G4LongLivedNeutralTransportation object at the start of a new track or the resumption
453// of a suspended track.
454//
455
456void
458{
459 G4VProcess::StartTracking(aTrack);
460 fNewTrack = true;
461 fFirstStepInVolume = true;
462 fLastStepInVolume = false;
463
464 // The actions here are those that were taken in AlongStepGPIL
465 // when track.GetCurrentStepNumber()==1
466
467 // Whether field exists should be determined at run level -- TODO
469
470 // reset safety value and center
471 //
472 fPreviousSafety = 0.0 ;
473 fPreviousSftOrigin = G4ThreeVector(0., 0., 0.) ;
474
475 // reset looping counter -- for motion in field
476 fNoLooperTrials = 0;
477 // Must clear this state .. else it depends on last track's value
478 // --> a better solution would set this from state of suspended track TODO ?
479 // Was if( aTrack->GetCurrentStepNumber()==1 ) { .. }
480
481 // ChordFinder reset internal state
482 //
484 fFieldPropagator->ClearPropagatorState();
485 // Resets all state of field propagator class (ONLY) including safety
486 // values (in case of overlaps and to wipe for first track).
487 }
488
489 // Make sure to clear the chord finders of all fields (i.e. managers)
490 //
491 G4FieldManagerStore* fieldMgrStore = G4FieldManagerStore::GetInstance();
492 fieldMgrStore->ClearAllChordFindersState();
493
494 // Update the current touchable handle (from the track's)
495 //
496 fCurrentTouchableHandle = aTrack->GetTouchableHandle();
497
498 // Inform field propagator of new track
499 //
500 if (fFieldPropagator) fFieldPropagator->PrepareNewTrack();
501}
502
504//
505
506
508//
509// Suppress (or not) warnings about 'looping' particles
510
512{
513 fSilenceLooperWarnings = val; // Flag to *Suppress* all 'looper' warnings
514 // G4CoupledTransportation::fSilenceLooperWarnings= val;
515}
516
518//
523
524
526//
528{
529 // Restores the old high values -- potentially appropriate for energy-frontier
530 // HEP experiments.
531 // Caution: All tracks with E < 100 MeV that are found to loop are
532 SetThresholdWarningEnergy(100.0 * CLHEP::MeV); // Warn above this energy
533 SetThresholdImportantEnergy(250.0 * CLHEP::MeV); // Extra trial above this En
534
535 G4int maxTrials = 10;
536 SetThresholdTrials(maxTrials);
537
538
539
540}
541
543void G4LongLivedNeutralTransportation::SetLowLooperThresholds() // Values for low-E applications
544{
545 // These values were the default in Geant4 10.5 - beta
546 SetThresholdWarningEnergy(1.0 * CLHEP::keV); // Warn above this En
547 SetThresholdImportantEnergy(1.0 * CLHEP::MeV); // Extra trials above it
548
549 G4int maxTrials = 30; // A new value - was 10
550 SetThresholdTrials(maxTrials);
551
552}
553
555//
556
558//
559
561//
563
564// StreamInfo(std::ostream& out, const G4ParticleDefinition& part, G4bool rst) const
565
566{
567 G4String indent = " "; // : "");
568 G4int oldPrec = outStream.precision(6);
569 // outStream << std::setprecision(6);
570 outStream << G4endl << indent << GetProcessName() << ": ";
571
572 outStream << " Parameters for looping particles: " << G4endl
573 << " warning-E = " << fThreshold_Warning_Energy / CLHEP::MeV << " MeV " << G4endl
574 << " important E = " << fThreshold_Important_Energy / CLHEP::MeV << " MeV " << G4endl
575 << " thresholdTrials " << fThresholdTrials << G4endl;
576 outStream.precision(oldPrec);
577}
G4TouchableHandle fCurrentTouchableHandle
Current touchable handle.
void PrintStatistics(std::ostream &outStr) const
returns current logging info of the algorithm
G4VParticleChange * PostStepDoIt(const G4Track &track, const G4Step &stepData)
G4VProcess::PostStepDoIt() implementation,.
G4double fCandidateEndGlobalTime
The particle's state after this Step, Store for DoIt.
G4double AlongStepGetPhysicalInteractionLength(const G4Track &track, G4double previousStepSize, G4double currentMinimumStep, G4double &currentSafety, G4GPILSelection *selection)
G4VProcess::AlongStepGetPhysicalInteractionLength() implementation,.
G4int fThresholdTrials
Number of trials an important looper survives.
void SetThresholdWarningEnergy(G4double newEnWarn)
Set fThreshold_Warning_Energy.
void SetThresholdImportantEnergy(G4double newEnImp)
Set fThreshold_Important_Energy.
virtual void ProcessDescription(std::ostream &outStream) const
G4LongLivedNeutralTransportation::ProcessDescription()
G4bool fGeometryLimitedStep
Flag to determine whether a boundary was reached.
G4ThreeVector fTransportEndMomentumDir
The particle's state after this Step, Store for DoIt.
void SetThresholdTrials(G4int newMaxTrials)
Set fThresholdTrials.
G4VParticleChange * AlongStepDoIt(const G4Track &track, const G4Step &stepData)
G4VProcess::AlongStepDoIt() implementation,.
void SetHighLooperThresholds()
Get/Set parameters for killing loopers: Above 'important' energy a 'looping' particle in field will N...
G4double fSumEnergyKilled
Sum of abandoned looping tracks energies.
G4ThreeVector fPreviousSftOrigin
Remember last safety origin.
static void SetSilenceLooperWarnings(G4bool val)
Do not warn about 'looping' particles.
G4double PostStepGetPhysicalInteractionLength(const G4Track &track, G4double previousStepSize, G4ForceCondition *pForceCond)
G4VProcess::PostStepGetPhysicalInteractionLength() implementation.
G4ThreeVector fTransportEndSpin
The particle's state after this Step, Store for DoIt.
G4Navigator * fLinearNavigator
Propagator used to transport the particle.
static G4bool fSilenceLooperWarnings
Flag to Supress all 'looper' warnings.
G4double fThreshold_Important_Energy
Give a few trial above this E for looping particle.
G4ThreeVector fTransportEndPosition
The particle's state after this Step, Store for DoIt.
G4bool fEndGlobalTimeComputed
The particle's state after this Step, Store for DoIt.
G4bool fShortStepOptimisation
Whether to avoid calling G4Navigator for short step ( < safety) If using it, the safety estimate for ...
G4int fNoLooperTrials
Counter for steps in which particle reports 'looping', if it is above 'Important' Energy.
G4bool fMomentumChanged
The particle's state after this Step, Store for DoIt.
G4LongLivedNeutralTransportation(G4int verbosityLevel=1)
Constructor.
void StartTracking(G4Track *aTrack)
Reset state for new (potentially resumed) track.
G4ParticleChangeForTransport fParticleChange
New ParticleChange.
void SetLowLooperThresholds()
Shortcut method - old values (meant for HEP)
G4double fThreshold_Warning_Energy
Warn above this energy about looping particle.
G4bool DoesGlobalFieldExist()
Checks whether a field exists for the "global" field manager.
static G4bool GetSilenceLooperWarnings()
Do not throw exception about 'looping' particles.
G4PropagatorInField * fFieldPropagator
Propagator used to transport the particle.
G4double fTransportEndKineticEnergy
The particle's state after this Step, Store for DoIt.
Abstract base class for different kinds of events.