ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4MonopoleTransportation.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4MonopoleTransportation.cc
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
28 //
29 //
30 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
31 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
32 //
33 // This class is a process responsible for the transportation of
34 // magnetic monopoles, ie the geometrical propagation that encounters the
35 // geometrical sub-volumes of the detectors.
36 //
37 // For monopoles, uses a different equation of motion and ignores energy
38 // conservation.
39 //
40 
41 // =======================================================================
42 // Created: 3 May 2010, J. Apostolakis, B. Bozsogi
43 // =======================================================================
44 
46 #include "G4ProductionCutsTable.hh"
47 #include "G4ParticleTable.hh"
48 #include "G4ChordFinder.hh"
49 #include "G4SafetyHelper.hh"
50 #include "G4FieldManagerStore.hh"
51 #include "G4Monopole.hh"
53 #include "G4SystemOfUnits.hh"
54 
55 #include "G4RunManager.hh"
56 #include "DetectorConstruction.hh"
57 
59 
60 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
61 
63  G4int verb)
64  : G4VProcess( G4String("MonopoleTransportation"), fTransportation ),
65  fParticleDef(mpl),
66  fMagSetup(0),
67  fLinearNavigator(0),
68  fFieldPropagator(0),
69  fParticleIsLooping( false ),
70  fPreviousSftOrigin (0.,0.,0.),
71  fPreviousSafety ( 0.0 ),
72  fThreshold_Warning_Energy( 100 * MeV ),
73  fThreshold_Important_Energy( 250 * MeV ),
74  fThresholdTrials( 10 ),
75  // fUnimportant_Energy( 1 * MeV ),
76  fNoLooperTrials(0),
77  fSumEnergyKilled( 0.0 ), fMaxEnergyKilled( 0.0 ),
78  fShortStepOptimisation(false), // Old default: true (=fast short steps)
79  fpSafetyHelper(0),
80  noCalls(0)
81 {
82  verboseLevel = verb;
83 
84  // set Process Sub Type
86 
87 
88  // Do not finalize the G4MonopoleTransportation class
90  {
91  return;
92  }
93 
94  const DetectorConstruction* detector =
95  static_cast<const DetectorConstruction*>
97 
98  fMagSetup = detector->GetMonopoleFieldSetup();
99 
100  G4TransportationManager* transportMgr =
102 
103  fLinearNavigator = transportMgr->GetNavigatorForTracking() ;
104 
105  fFieldPropagator = transportMgr->GetPropagatorInField() ;
106  fpSafetyHelper = transportMgr->GetSafetyHelper();
107 
108 // New
109 
110  // Cannot determine whether a field exists here,
111  // because it would only work if the field manager has informed
112  // about the detector's field before this transportation process
113  // is constructed.
114  // Instead later the method DoesGlobalFieldExist() is called
115  fCurrentTouchableHandle = nullptr;
116 
117  fEndGlobalTimeComputed = false;
119 }
120 
121 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
122 
124 {
125  if( (verboseLevel > 0) && (fSumEnergyKilled > 0.0 ) ){
126  G4cout << " G4MonopoleTransportation: Statistics for looping particles "
127  << G4endl;
128  G4cout << " Sum of energy of loopers killed: "
129  << fSumEnergyKilled << G4endl;
130  G4cout << " Max energy of loopers killed: "
131  << fMaxEnergyKilled << G4endl;
132  }
133 }
134 
135 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
136 //
137 // Responsibilities:
138 // Find whether the geometry limits the Step, and to what length
139 // Calculate the new value of the safety and return it.
140 // Store the final time, position and momentum.
141 
144  G4double, // previousStepSize
145  G4double currentMinimumStep,
146  G4double& currentSafety,
147  G4GPILSelection* selection )
148 {
150  // change to monopole equation
151 
152  G4double geometryStepLength, newSafety ;
154 
155  // Initial actions moved to StartTrack()
156  // --------------------------------------
157  // Note: in case another process changes touchable handle
158  // it will be necessary to add here (for all steps)
159  // fCurrentTouchableHandle = aTrack->GetTouchableHandle();
160 
161  // GPILSelection is set to defaule value of CandidateForSelection
162  // It is a return value
163  //
164  *selection = CandidateForSelection ;
165 
166  // Get initial Energy/Momentum of the track
167  //
168  const G4DynamicParticle* pParticle = track.GetDynamicParticle() ;
169  G4ThreeVector startMomentumDir = pParticle->GetMomentumDirection() ;
170  G4ThreeVector startPosition = track.GetPosition() ;
171 
172  // G4double theTime = track.GetGlobalTime() ;
173 
174  // The Step Point safety can be limited by other geometries and/or the
175  // assumptions of any process - it's not always the geometrical safety.
176  // We calculate the starting point's isotropic safety here.
177  //
178  G4ThreeVector OriginShift = startPosition - fPreviousSftOrigin ;
179  G4double MagSqShift = OriginShift.mag2() ;
180  if( MagSqShift >= sqr(fPreviousSafety) )
181  {
182  currentSafety = 0.0 ;
183  }
184  else
185  {
186  currentSafety = fPreviousSafety - std::sqrt(MagSqShift) ;
187  }
188 
189  // Is the monopole charged ?
190  //
191  G4double particleMagneticCharge = fParticleDef->MagneticCharge() ;
192  G4double particleElectricCharge = pParticle->GetCharge();
193 
195  // fEndGlobalTimeComputed = false ;
196 
197  // There is no need to locate the current volume. It is Done elsewhere:
198  // On track construction
199  // By the tracking, after all AlongStepDoIts, in "Relocation"
200 
201  // Check whether the particle have an (EM) field force exerting upon it
202  //
203  G4FieldManager* fieldMgr=0;
204  G4bool fieldExertsForce = false ;
205 
206  if( (particleMagneticCharge != 0.0) )
207  {
208  fieldMgr= fFieldPropagator->FindAndSetFieldManager( track.GetVolume() );
209  if (fieldMgr != 0) {
210  // Message the field Manager, to configure it for this track
211  fieldMgr->ConfigureForTrack( &track );
212  // Moved here, in order to allow a transition
213  // from a zero-field status (with fieldMgr->(field)0
214  // to a finite field status
215 
216  // If the field manager has no field, there is no field !
217  fieldExertsForce = (fieldMgr->GetDetectorField() != 0);
218  }
219  }
220 
221  // G4cout << " G4Transport: field exerts force= " << fieldExertsForce
222  // << " fieldMgr= " << fieldMgr << G4endl;
223 
224  // Choose the calculation of the transportation: Field or not
225  //
226  if( !fieldExertsForce )
227  {
228  G4double linearStepLength ;
229  if( fShortStepOptimisation && (currentMinimumStep <= currentSafety) )
230  {
231  // The Step is guaranteed to be taken
232  //
233  geometryStepLength = currentMinimumStep ;
235  }
236  else
237  {
238  // Find whether the straight path intersects a volume
239  //
240  linearStepLength = fLinearNavigator->ComputeStep( startPosition,
241  startMomentumDir,
242  currentMinimumStep,
243  newSafety) ;
244  // Remember last safety origin & value.
245  //
246  fPreviousSftOrigin = startPosition ;
247  fPreviousSafety = newSafety ;
248  // fpSafetyHelper->SetCurrentSafety( newSafety, startPosition);
249 
250  // The safety at the initial point has been re-calculated:
251  //
252  currentSafety = newSafety ;
253 
254  fGeometryLimitedStep= (linearStepLength <= currentMinimumStep);
256  {
257  // The geometry limits the Step size (an intersection was found.)
258  geometryStepLength = linearStepLength ;
259  }
260  else
261  {
262  // The full Step is taken.
263  geometryStepLength = currentMinimumStep ;
264  }
265  }
266  endpointDistance = geometryStepLength ;
267 
268  // Calculate final position
269  //
270  fTransportEndPosition = startPosition+geometryStepLength*startMomentumDir ;
271 
272  // Momentum direction, energy and polarisation are unchanged by transport
273  //
274  fTransportEndMomentumDir = startMomentumDir ;
280  }
281  else // A field exerts force
282  {
283  G4double momentumMagnitude = pParticle->GetTotalMomentum() ;
284  G4ThreeVector EndUnitMomentum ;
285  G4double lengthAlongCurve ;
286  G4double restMass = fParticleDef->GetPDGMass() ;
287 
288  G4ChargeState chargeState(particleElectricCharge, // The charge can change
290  0, // Magnetic moment:
291  0, // Electric Dipole moment
292  particleMagneticCharge ); // in Mev/c
293 
294  G4EquationOfMotion* equationOfMotion = fFieldPropagator->
295  GetChordFinder()->GetIntegrationDriver()->GetEquationOfMotion();
296 
297  equationOfMotion->SetChargeMomentumMass( chargeState,
298  momentumMagnitude,
299  restMass );
300  // SetChargeMomentumMass now passes both the electric and magnetic
301  // charge in chargeState
302 
303  G4ThreeVector spin = track.GetPolarization() ;
304  G4FieldTrack aFieldTrack = G4FieldTrack( startPosition,
305  track.GetMomentumDirection(),
306  0.0,
307  track.GetKineticEnergy(),
308  restMass,
309  track.GetVelocity(),
310  track.GetGlobalTime(), // Lab.
311  track.GetProperTime(), // Part.
312  &spin ) ;
313  if( currentMinimumStep > 0 )
314  {
315  // Do the Transport in the field (non recti-linear)
316  //
317  lengthAlongCurve = fFieldPropagator->ComputeStep( aFieldTrack,
318  currentMinimumStep,
319  currentSafety,
320  track.GetVolume() ) ;
321  fGeometryLimitedStep= lengthAlongCurve < currentMinimumStep;
322  if( fGeometryLimitedStep ) {
323  geometryStepLength = lengthAlongCurve ;
324  } else {
325  geometryStepLength = currentMinimumStep ;
326  }
327  }
328  else
329  {
330  geometryStepLength = lengthAlongCurve= 0.0 ;
332  }
333 
334  // Remember last safety origin & value.
335  //
336  fPreviousSftOrigin = startPosition ;
337  fPreviousSafety = currentSafety ;
338  // fpSafetyHelper->SetCurrentSafety( newSafety, startPosition);
339 
340  // Get the End-Position and End-Momentum (Dir-ection)
341  //
342  fTransportEndPosition = aFieldTrack.GetPosition() ;
343 
344  // Momentum: Magnitude and direction can be changed too now ...
345  //
347  fTransportEndMomentumDir = aFieldTrack.GetMomentumDir() ;
348 
350 
352  fEndGlobalTimeComputed = true;
353 
354  fTransportEndSpin = aFieldTrack.GetSpin();
356  endpointDistance = (fTransportEndPosition - startPosition).mag() ;
357  }
358 
359  // If we are asked to go a step length of 0, and we are on a boundary
360  // then a boundary will also limit the step -> we must flag this.
361  //
362  if( currentMinimumStep == 0.0 )
363  {
364  if( currentSafety == 0.0 ) fGeometryLimitedStep = true ;
365  }
366 
367  // Update the safety starting from the end-point,
368  // if it will become negative at the end-point.
369  //
370  if( currentSafety < endpointDistance )
371  {
372  // if( particleMagneticCharge == 0.0 )
373  //G4cout << " Avoiding call to ComputeSafety : charge = 0.0 " << G4endl;
374 
375  if( particleMagneticCharge != 0.0 ) {
376 
377  G4double endSafety =
379  currentSafety = endSafety ;
380  fPreviousSftOrigin = fTransportEndPosition ;
381  fPreviousSafety = currentSafety ;
383 
384  // Because the Stepping Manager assumes it is from the start point,
385  // add the StepLength
386  //
387  currentSafety += endpointDistance ;
388 
389 #ifdef G4DEBUG_TRANSPORT
390  G4cout.precision(12) ;
391  G4cout << "***G4MonopoleTransportation::AlongStepGPIL ** " << G4endl;
392  G4cout << " Called Navigator->ComputeSafety at "
394  << " and it returned safety= " << endSafety << G4endl;
395  G4cout << " Adding endpoint distance " << endpointDistance
396  << " to obtain pseudo-safety= " << currentSafety << G4endl;
397 #endif
398  }
399  }
400 
401  fParticleChange.ProposeTrueStepLength(geometryStepLength) ;
402 
404  // change back to usual equation
405 
406  return geometryStepLength ;
407 }
408 
409 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
410 //
411 // Initialize ParticleChange (by setting all its members equal
412 // to corresponding members in G4Track)
413 
415  const G4Track& track, const G4Step& stepData )
416 {
417  ++noCalls;
418 
419  fParticleChange.Initialize(track) ;
420 
421  // Code for specific process
422  //
427 
429 
430  G4double deltaTime = 0.0 ;
431 
432  // Calculate Lab Time of Flight (ONLY if field Equations used it!)
433 
434  G4double startTime = track.GetGlobalTime() ;
435 
437  {
438  // The time was not integrated .. make the best estimate possible
439  //
440  G4double finalVelocity = track.GetVelocity() ;
441  G4double initialVelocity = stepData.GetPreStepPoint()->GetVelocity() ;
442  G4double stepLength = track.GetStepLength() ;
443 
444  deltaTime= 0.0; // in case initialVelocity = 0
445  if (finalVelocity > 0.0)
446  {
447  G4double meanInverseVelocity ;
448  // deltaTime = stepLength/finalVelocity ;
449  meanInverseVelocity = 0.5
450  * ( 1.0 / initialVelocity + 1.0 / finalVelocity );
451  deltaTime = stepLength * meanInverseVelocity ;
452  }
453  else if( initialVelocity > 0.0 )
454  {
455  deltaTime = stepLength/initialVelocity ;
456  }
457  fCandidateEndGlobalTime = startTime + deltaTime ;
458  }
459  else
460  {
461  deltaTime = fCandidateEndGlobalTime - startTime ;
462  }
463 
465 
466  // Now Correct by Lorentz factor to get "proper" deltaTime
467 
468  G4double restMass = track.GetDynamicParticle()->GetMass();
469  G4double deltaProperTime = deltaTime*( restMass/track.GetTotalEnergy() );
470 
471  fParticleChange.ProposeProperTime(track.GetProperTime() + deltaProperTime);
472 
473  // If the particle is caught looping or is stuck (in very difficult
474  // boundaries) in a magnetic field (doing many steps)
475  // THEN this kills it ...
476  //
477  if ( fParticleIsLooping )
478  {
480 
481  if( (endEnergy < fThreshold_Important_Energy)
483  // Kill the looping particle
484  //
486 
487  // 'Bare' statistics
488  fSumEnergyKilled += endEnergy;
489  if( endEnergy > fMaxEnergyKilled) { fMaxEnergyKilled= endEnergy; }
490 
491 #ifdef G4VERBOSE
492  if( (verboseLevel > 1) ||
493  ( endEnergy > fThreshold_Warning_Energy ) ) {
494  G4cout << " G4MonopoleTransportation is killing track "
495  << "that is looping or stuck " << G4endl
496  << " This track has " << track.GetKineticEnergy() / MeV
497  << " MeV energy." << G4endl;
498  G4cout << " Number of trials = " << fNoLooperTrials
499  << " No of calls to AlongStepDoIt = " << noCalls
500  << G4endl;
501  }
502 #endif
503  fNoLooperTrials=0;
504  }
505  else{
506  fNoLooperTrials ++;
507 #ifdef G4VERBOSE
508  if( (verboseLevel > 2) ){
509  G4cout << " G4MonopoleTransportation::AlongStepDoIt(): "
510  << "Particle looping - "
511  << " Number of trials = " << fNoLooperTrials
512  << " No of calls to = " << noCalls
513  << G4endl;
514  }
515 #endif
516  }
517  }else{
518  fNoLooperTrials=0;
519  }
520 
521  // Another (sometimes better way) is to use a user-limit maximum Step size
522  // to alleviate this problem ..
523 
524  // Introduce smooth curved trajectories to particle-change
525  //
528 
529  return &fParticleChange ;
530 }
531 
532 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
533 //
534 // This ensures that the PostStep action is always called,
535 // so that it can do the relocation if it is needed.
536 //
537 
540  G4double, // previousStepSize
541  G4ForceCondition* pForceCond )
542 {
543  *pForceCond = Forced ;
544  return DBL_MAX ; // was kInfinity ; but convention now is DBL_MAX
545 }
546 
547 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
548 
550  const G4Step& )
551 {
552  G4TouchableHandle retCurrentTouchable ; // The one to return
553 
554  // Initialize ParticleChange (by setting all its members equal
555  // to corresponding members in G4Track)
556  // fParticleChange.Initialize(track) ; // To initialise TouchableChange
557 
559 
560  // If the Step was determined by the volume boundary,
561  // logically relocate the particle
562 
564  {
565  // fCurrentTouchable will now become the previous touchable,
566  // and what was the previous will be freed.
567  // (Needed because the preStepPoint can point to the previous touchable)
568 
571  LocateGlobalPointAndUpdateTouchableHandle( track.GetPosition(),
572  track.GetMomentumDirection(),
574  true ) ;
575  // Check whether the particle is out of the world volume
576  // If so it has exited and must be killed.
577  //
578  if( fCurrentTouchableHandle->GetVolume() == 0 )
579  {
581  }
582  retCurrentTouchable = fCurrentTouchableHandle ;
584  }
585  else // fGeometryLimitedStep is false
586  {
587  // This serves only to move the Navigator's location
588  //
590 
591  // The value of the track's current Touchable is retained.
592  // (and it must be correct because we must use it below to
593  // overwrite the (unset) one in particle change)
594  // It must be fCurrentTouchable too ??
595  //
597  retCurrentTouchable = track.GetTouchableHandle() ;
598  } // endif ( fGeometryLimitedStep )
599 
600  const G4VPhysicalVolume* pNewVol = retCurrentTouchable->GetVolume() ;
601  const G4Material* pNewMaterial = 0 ;
602  const G4VSensitiveDetector* pNewSensitiveDetector = 0 ;
603  if( pNewVol != 0 )
604  {
605  pNewMaterial= pNewVol->GetLogicalVolume()->GetMaterial();
606  pNewSensitiveDetector= pNewVol->GetLogicalVolume()->GetSensitiveDetector();
607  }
608 
609  // ( <const_cast> pNewMaterial ) ;
610  // ( <const_cast> pNewSensitiveDetector) ;
611 
613  (G4Material *) pNewMaterial ) ;
615  (G4VSensitiveDetector *) pNewSensitiveDetector ) ;
616 
617  const G4MaterialCutsCouple* pNewMaterialCutsCouple = 0;
618  if( pNewVol != 0 )
619  {
620  pNewMaterialCutsCouple=pNewVol->GetLogicalVolume()->GetMaterialCutsCouple();
621  }
622 
623  if( pNewVol!=0 && pNewMaterialCutsCouple!=0 &&
624  pNewMaterialCutsCouple->GetMaterial()!=pNewMaterial )
625  {
626  // for parametrized volume
627  //
628  pNewMaterialCutsCouple =
630  ->GetMaterialCutsCouple(pNewMaterial,
631  pNewMaterialCutsCouple->GetProductionCuts());
632  }
633  fParticleChange.SetMaterialCutsCoupleInTouchable( pNewMaterialCutsCouple );
634 
635  // temporarily until Get/Set Material of ParticleChange,
636  // and StepPoint can be made const.
637  // Set the touchable in ParticleChange
638  // this must always be done because the particle change always
639  // uses this value to overwrite the current touchable pointer.
640  //
641  fParticleChange.SetTouchableHandle(retCurrentTouchable) ;
642 
643  return &fParticleChange ;
644 }
645 
646 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
647 
648 // New method takes over the responsibility to reset the state
649 // of G4MonopoleTransportation object at the start of a new track
650 // or the resumption of a suspended track.
651 
652 void
654 {
656 
657 // The actions here are those that were taken in AlongStepGPIL
658 // when track.GetCurrentStepNumber()==1
659 
660  // reset safety value and center
661  //
662  fPreviousSafety = 0.0 ;
663  fPreviousSftOrigin = G4ThreeVector(0.,0.,0.) ;
664 
665  // reset looping counter -- for motion in field
666  fNoLooperTrials= 0;
667  // Must clear this state .. else it depends on last track's value
668  // --> a better solution would set this from state of suspended track TODO ?
669  // Was if( aTrack->GetCurrentStepNumber()==1 ) { .. }
670 
671  // ChordFinder reset internal state
672  //
673  if( DoesGlobalFieldExist() ) {
675  // Resets all state of field propagator class (ONLY)
676  // including safety values
677  // in case of overlaps and to wipe for first track).
678 
679  // G4ChordFinder* chordF= fFieldPropagator->GetChordFinder();
680  // if( chordF ) chordF->ResetStepEstimate();
681  }
682 
683  // Make sure to clear the chord finders of all fields (ie managers)
685  fieldMgrStore->ClearAllChordFindersState();
686 
687  // Update the current touchable handle (from the track's)
688  //
690 }
691 
692 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
693