| 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 /// \file exoticphysics/monopole/src/G4Monopol
27 /// \brief Implementation of the G4MonopoleTra
28 //
29 //
30 //....oooOO0OOooo........oooOO0OOooo........oo
31 //....oooOO0OOooo........oooOO0OOooo........oo
32 //
33 // This class is a process responsible for the
34 // magnetic monopoles, ie the geometrical prop
35 // geometrical sub-volumes of the detectors.
36 //
37 // For monopoles, uses a different equation of
38 // conservation.
39 //
40
41 // ===========================================
42 // Created: 3 May 2010, J. Apostolakis, B. Bo
43 // ===========================================
44
45 #include "G4MonopoleTransportation.hh"
46
47 #include "DetectorConstruction.hh"
48
49 #include "G4ChordFinder.hh"
50 #include "G4FieldManagerStore.hh"
51 #include "G4Monopole.hh"
52 #include "G4ParticleTable.hh"
53 #include "G4ProductionCutsTable.hh"
54 #include "G4RunManager.hh"
55 #include "G4SafetyHelper.hh"
56 #include "G4SystemOfUnits.hh"
57 #include "G4TransportationProcessType.hh"
58
59 class G4VSensitiveDetector;
60
61 //....oooOO0OOooo........oooOO0OOooo........oo
62
63 G4MonopoleTransportation::G4MonopoleTransporta
64 : G4VProcess(G4String("MonopoleTransportatio
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),
78 fMaxEnergyKilled(0.0),
79 fShortStepOptimisation(false), // Old def
80 fpSafetyHelper(0),
81 noCalls(0)
82 {
83 verboseLevel = verb;
84
85 // set Process Sub Type
86 SetProcessSubType(TRANSPORTATION);
87
88 // Do not finalize the G4MonopoleTransportat
89 if (G4Threading::IsMasterThread() && G4Threa
90 return;
91 }
92
93 const DetectorConstruction* detector = stati
94 G4RunManager::GetRunManager()->GetUserDete
95
96 fMagSetup = detector->GetMonopoleFieldSetup(
97
98 G4TransportationManager* transportMgr = G4Tr
99
100 fLinearNavigator = transportMgr->GetNavigato
101
102 fFieldPropagator = transportMgr->GetPropagat
103 fpSafetyHelper = transportMgr->GetSafetyHelp
104
105 // New
106
107 // Cannot determine whether a field exists h
108 // because it would only work if the field
109 // about the detector's field before this t
110 // is constructed.
111 // Instead later the method DoesGlobalFieldE
112 fCurrentTouchableHandle = nullptr;
113
114 fEndGlobalTimeComputed = false;
115 fCandidateEndGlobalTime = 0;
116 }
117
118 //....oooOO0OOooo........oooOO0OOooo........oo
119
120 G4MonopoleTransportation::~G4MonopoleTransport
121 {
122 if ((verboseLevel > 0) && (fSumEnergyKilled
123 G4cout << " G4MonopoleTransportation: Stat
124 G4cout << " Sum of energy of loopers kil
125 G4cout << " Max energy of loopers killed
126 }
127 }
128
129 //....oooOO0OOooo........oooOO0OOooo........oo
130 //
131 // Responsibilities:
132 // Find whether the geometry limits the Ste
133 // Calculate the new value of the safety an
134 // Store the final time, position and momen
135
136 G4double G4MonopoleTransportation::AlongStepGe
137 const G4Track& track,
138 G4double, // previousStepSize
139 G4double currentMinimumStep, G4double& curre
140 {
141 fMagSetup->SetStepperAndChordFinder(1);
142 // change to monopole equation
143
144 G4double geometryStepLength, newSafety;
145 fParticleIsLooping = false;
146
147 // Initial actions moved to StartTrack()
148 // --------------------------------------
149 // Note: in case another process changes tou
150 // it will be necessary to add here (for
151 // fCurrentTouchableHandle = aTrack->GetTouc
152
153 // GPILSelection is set to defaule value of
154 // It is a return value
155 //
156 *selection = CandidateForSelection;
157
158 // Get initial Energy/Momentum of the track
159 //
160 const G4DynamicParticle* pParticle = track.G
161 G4ThreeVector startMomentumDir = pParticle->
162 G4ThreeVector startPosition = track.GetPosit
163
164 // G4double theTime = track.GetGlob
165
166 // The Step Point safety can be limited by o
167 // assumptions of any process - it's not alw
168 // We calculate the starting point's isotrop
169 //
170 G4ThreeVector OriginShift = startPosition -
171 G4double MagSqShift = OriginShift.mag2();
172 if (MagSqShift >= sqr(fPreviousSafety)) {
173 currentSafety = 0.0;
174 }
175 else {
176 currentSafety = fPreviousSafety - std::sqr
177 }
178
179 // Is the monopole charged ?
180 //
181 G4double particleMagneticCharge = fParticleD
182 G4double particleElectricCharge = pParticle-
183
184 fGeometryLimitedStep = false;
185 // fEndGlobalTimeComputed = false ;
186
187 // There is no need to locate the current vo
188 // On track construction
189 // By the tracking, after all AlongStepDoI
190
191 // Check whether the particle have an (EM) f
192 //
193 G4FieldManager* fieldMgr = 0;
194 G4bool fieldExertsForce = false;
195
196 if ((particleMagneticCharge != 0.0)) {
197 fieldMgr = fFieldPropagator->FindAndSetFie
198 if (fieldMgr != 0) {
199 // Message the field Manager, to configu
200 fieldMgr->ConfigureForTrack(&track);
201 // Moved here, in order to allow a trans
202 // from a zero-field status (with fie
203 // to a finite field status
204
205 // If the field manager has no field, th
206 fieldExertsForce = (fieldMgr->GetDetecto
207 }
208 }
209
210 // G4cout << " G4Transport: field exerts fo
211 // << " fieldMgr= " << fieldMgr <<
212
213 // Choose the calculation of the transportat
214 //
215 if (!fieldExertsForce) {
216 G4double linearStepLength;
217 if (fShortStepOptimisation && (currentMini
218 // The Step is guaranteed to be taken
219 //
220 geometryStepLength = currentMinimumStep;
221 fGeometryLimitedStep = false;
222 }
223 else {
224 // Find whether the straight path inter
225 //
226 linearStepLength = fLinearNavigator->Com
227
228 // Remember last safety origin & value.
229 //
230 fPreviousSftOrigin = startPosition;
231 fPreviousSafety = newSafety;
232 // fpSafetyHelper->SetCurrentSafety( new
233
234 // The safety at the initial point has b
235 //
236 currentSafety = newSafety;
237
238 fGeometryLimitedStep = (linearStepLength
239 if (fGeometryLimitedStep) {
240 // The geometry limits the Step size (
241 geometryStepLength = linearStepLength;
242 }
243 else {
244 // The full Step is taken.
245 geometryStepLength = currentMinimumSte
246 }
247 }
248 endpointDistance = geometryStepLength;
249
250 // Calculate final position
251 //
252 fTransportEndPosition = startPosition + ge
253
254 // Momentum direction, energy and polarisa
255 //
256 fTransportEndMomentumDir = startMomentumDi
257 fTransportEndKineticEnergy = track.GetKine
258 fTransportEndSpin = track.GetPolarization(
259 fParticleIsLooping = false;
260 fMomentumChanged = false;
261 fEndGlobalTimeComputed = false;
262 }
263 else // A field exerts force
264 {
265 G4double momentumMagnitude = pParticle->Ge
266 G4ThreeVector EndUnitMomentum;
267 G4double lengthAlongCurve;
268 G4double restMass = fParticleDef->GetPDGMa
269
270 G4ChargeState chargeState(particleElectric
271 fParticleDef->Ge
272 0, // Magnetic
273 0, // Electric
274 particleMagnetic
275
276 G4EquationOfMotion* equationOfMotion =
277 fFieldPropagator->GetChordFinder()->GetI
278
279 equationOfMotion->SetChargeMomentumMass(ch
280 // SetChargeMomentumMass now passes both t
281 // charge in chargeState
282
283 G4ThreeVector spin = track.GetPolarization
284 G4FieldTrack aFieldTrack = G4FieldTrack(st
285 tr
286 tr
287 tr
288 &s
289 if (currentMinimumStep > 0) {
290 // Do the Transport in the field (non re
291 //
292 lengthAlongCurve = fFieldPropagator->Com
293
294 fGeometryLimitedStep = lengthAlongCurve
295 if (fGeometryLimitedStep) {
296 geometryStepLength = lengthAlongCurve;
297 }
298 else {
299 geometryStepLength = currentMinimumSte
300 }
301 }
302 else {
303 geometryStepLength = lengthAlongCurve =
304 fGeometryLimitedStep = false;
305 }
306
307 // Remember last safety origin & value.
308 //
309 fPreviousSftOrigin = startPosition;
310 fPreviousSafety = currentSafety;
311 // fpSafetyHelper->SetCurrentSafety( newSa
312
313 // Get the End-Position and End-Momentum (
314 //
315 fTransportEndPosition = aFieldTrack.GetPos
316
317 // Momentum: Magnitude and direction can
318 //
319 fMomentumChanged = true;
320 fTransportEndMomentumDir = aFieldTrack.Get
321
322 fTransportEndKineticEnergy = aFieldTrack.G
323
324 fCandidateEndGlobalTime = aFieldTrack.GetL
325 fEndGlobalTimeComputed = true;
326
327 fTransportEndSpin = aFieldTrack.GetSpin();
328 fParticleIsLooping = fFieldPropagator->IsP
329 endpointDistance = (fTransportEndPosition
330 }
331
332 // If we are asked to go a step length of 0,
333 // then a boundary will also limit the step
334 //
335 if (currentMinimumStep == 0.0) {
336 if (currentSafety == 0.0) fGeometryLimited
337 }
338
339 // Update the safety starting from the end-p
340 // if it will become negative at the end-poi
341 //
342 if (currentSafety < endpointDistance) {
343 // if( particleMagneticCharge == 0.0 )
344 // G4cout << " Avoiding call to ComputeS
345
346 if (particleMagneticCharge != 0.0) {
347 G4double endSafety = fLinearNavigator->C
348 currentSafety = endSafety;
349 fPreviousSftOrigin = fTransportEndPositi
350 fPreviousSafety = currentSafety;
351 fpSafetyHelper->SetCurrentSafety(current
352
353 // Because the Stepping Manager assumes
354 // add the StepLength
355 //
356 currentSafety += endpointDistance;
357
358 #ifdef G4DEBUG_TRANSPORT
359 G4cout.precision(12);
360 G4cout << "***G4MonopoleTransportation::
361 G4cout << " Called Navigator->ComputeSa
362 << " and it returned safety= "
363 G4cout << " Adding endpoint distance "
364 << " to obtain pseudo-safety= "
365 #endif
366 }
367 }
368
369 fParticleChange.ProposeTrueStepLength(geomet
370
371 fMagSetup->SetStepperAndChordFinder(0);
372 // change back to usual equation
373
374 return geometryStepLength;
375 }
376
377 //....oooOO0OOooo........oooOO0OOooo........oo
378 //
379 // Initialize ParticleChange (by setting al
380 // to correspond
381
382 G4VParticleChange* G4MonopoleTransportation::A
383
384 {
385 ++noCalls;
386
387 if (fGeometryLimitedStep) {
388 stepData.GetPostStepPoint()->SetStepStatus
389 }
390
391 fParticleChange.Initialize(track);
392
393 // Code for specific process
394 //
395 fParticleChange.ProposePosition(fTransportEn
396 fParticleChange.ProposeMomentumDirection(fTr
397 fParticleChange.ProposeEnergy(fTransportEndK
398 fParticleChange.SetMomentumChanged(fMomentum
399
400 fParticleChange.ProposePolarization(fTranspo
401
402 G4double deltaTime = 0.0;
403
404 // Calculate Lab Time of Flight (ONLY if fi
405
406 G4double startTime = track.GetGlobalTime();
407
408 if (!fEndGlobalTimeComputed) {
409 // The time was not integrated .. make the
410 //
411 G4double finalVelocity = track.GetVelocity
412 G4double initialVelocity = stepData.GetPre
413 G4double stepLength = track.GetStepLength(
414
415 deltaTime = 0.0; // in case initialVeloci
416 if (finalVelocity > 0.0) {
417 G4double meanInverseVelocity;
418 // deltaTime = stepLength/finalVelocity
419 meanInverseVelocity = 0.5 * (1.0 / initi
420 deltaTime = stepLength * meanInverseVelo
421 }
422 else if (initialVelocity > 0.0) {
423 deltaTime = stepLength / initialVelocity
424 }
425 fCandidateEndGlobalTime = startTime + delt
426 }
427 else {
428 deltaTime = fCandidateEndGlobalTime - star
429 }
430
431 fParticleChange.ProposeGlobalTime(fCandidate
432
433 // Now Correct by Lorentz factor to get "pro
434
435 G4double restMass = track.GetDynamicParticle
436 G4double deltaProperTime = deltaTime * (rest
437
438 fParticleChange.ProposeProperTime(track.GetP
439
440 // If the particle is caught looping or is s
441 // boundaries) in a magnetic field (doing ma
442 // THEN this kills it ...
443 //
444 if (fParticleIsLooping) {
445 G4double endEnergy = fTransportEndKineticE
446
447 if ((endEnergy < fThreshold_Important_Ener
448 // Kill the looping particle
449 //
450 fParticleChange.ProposeTrackStatus(fStop
451
452 // 'Bare' statistics
453 fSumEnergyKilled += endEnergy;
454 if (endEnergy > fMaxEnergyKilled) {
455 fMaxEnergyKilled = endEnergy;
456 }
457
458 #ifdef G4VERBOSE
459 if ((verboseLevel > 1) || (endEnergy > f
460 G4cout << " G4MonopoleTransportation i
461 << "that is looping or stuck "
462 << track.GetKineticEnergy() / M
463 G4cout << " Number of trials = " <<
464 << " No of calls to AlongStep
465 }
466 #endif
467 fNoLooperTrials = 0;
468 }
469 else {
470 fNoLooperTrials++;
471 #ifdef G4VERBOSE
472 if ((verboseLevel > 2)) {
473 G4cout << " G4MonopoleTransportation
474 << "Particle looping - "
475 << " Number of trials = " <<
476 << G4endl;
477 }
478 #endif
479 }
480 }
481 else {
482 fNoLooperTrials = 0;
483 }
484
485 // Another (sometimes better way) is to use
486 // to alleviate this problem ..
487
488 // Introduce smooth curved trajectories to p
489 //
490 fParticleChange.SetPointerToVectorOfAuxiliar
491 fFieldPropagator->GimmeTrajectoryVectorAnd
492
493 return &fParticleChange;
494 }
495
496 //....oooOO0OOooo........oooOO0OOooo........oo
497 //
498 // This ensures that the PostStep action is a
499 // so that it can do the relocation if it is
500 //
501
502 G4double
503 G4MonopoleTransportation::PostStepGetPhysicalI
504
505
506 {
507 *pForceCond = Forced;
508 return DBL_MAX; // was kInfinity ; but conv
509 }
510
511 //....oooOO0OOooo........oooOO0OOooo........oo
512
513 G4VParticleChange* G4MonopoleTransportation::P
514 {
515 G4TouchableHandle retCurrentTouchable; // T
516
517 // Initialize ParticleChange (by setting al
518 // to correspond
519 // fParticleChange.Initialize(track) ; // T
520
521 fParticleChange.ProposeTrackStatus(track.Get
522
523 // If the Step was determined by the volume
524 // logically relocate the particle
525
526 if (fGeometryLimitedStep) {
527 // fCurrentTouchable will now become the p
528 // and what was the previous will be freed
529 // (Needed because the preStepPoint can po
530
531 fLinearNavigator->SetGeometricallyLimitedS
532 fLinearNavigator->LocateGlobalPointAndUpda
533 track.GetPosition(), track.GetMomentumDi
534 // Check whether the particle is out of th
535 // If so it has exited and must be killed.
536 //
537 if (fCurrentTouchableHandle->GetVolume() =
538 fParticleChange.ProposeTrackStatus(fStop
539 }
540 retCurrentTouchable = fCurrentTouchableHan
541 fParticleChange.SetTouchableHandle(fCurren
542 }
543 else // fGeometryLimitedStep is false
544 {
545 // This serves only to move the Navigator'
546 //
547 fLinearNavigator->LocateGlobalPointWithinV
548
549 // The value of the track's current Toucha
550 // (and it must be correct because we must
551 // overwrite the (unset) one in particle c
552 // It must be fCurrentTouchable too ??
553 //
554 fParticleChange.SetTouchableHandle(track.G
555 retCurrentTouchable = track.GetTouchableHa
556 } // endif ( fGeometryLimitedStep )
557
558 const G4VPhysicalVolume* pNewVol = retCurren
559 const G4Material* pNewMaterial = 0;
560 G4VSensitiveDetector* pNewSensitiveDetector
561 if (pNewVol != 0) {
562 pNewMaterial = pNewVol->GetLogicalVolume()
563 pNewSensitiveDetector = pNewVol->GetLogica
564 }
565
566 // ( <const_cast> pNewMaterial ) ;
567
568 fParticleChange.SetMaterialInTouchable((G4Ma
569 fParticleChange.SetSensitiveDetectorInToucha
570
571 const G4MaterialCutsCouple* pNewMaterialCuts
572 if (pNewVol != 0) {
573 pNewMaterialCutsCouple = pNewVol->GetLogic
574 }
575
576 if (pNewVol != 0 && pNewMaterialCutsCouple !
577 && pNewMaterialCutsCouple->GetMaterial()
578 {
579 // for parametrized volume
580 //
581 pNewMaterialCutsCouple = G4ProductionCutsT
582 pNewMaterial, pNewMaterialCutsCouple->Ge
583 }
584 fParticleChange.SetMaterialCutsCoupleInTouch
585
586 // temporarily until Get/Set Material of Par
587 // and StepPoint can be made const.
588 // Set the touchable in ParticleChange
589 // this must always be done because the part
590 // uses this value to overwrite the current
591 //
592 fParticleChange.SetTouchableHandle(retCurren
593
594 return &fParticleChange;
595 }
596
597 //....oooOO0OOooo........oooOO0OOooo........oo
598
599 // New method takes over the responsibility to
600 // of G4MonopoleTransportation object at the s
601 // or the resumption of a suspended track.
602
603 void G4MonopoleTransportation::StartTracking(G
604 {
605 G4VProcess::StartTracking(aTrack);
606
607 // The actions here are those that were take
608 // when track.GetCurrentStepNumber()==1
609
610 // reset safety value and center
611 //
612 fPreviousSafety = 0.0;
613 fPreviousSftOrigin = G4ThreeVector(0., 0., 0
614
615 // reset looping counter -- for motion in fi
616 fNoLooperTrials = 0;
617 // Must clear this state .. else it depends
618 // --> a better solution would set this fro
619 // Was if( aTrack->GetCurrentStepNumber()==1
620
621 // ChordFinder reset internal state
622 //
623 if (DoesGlobalFieldExist()) {
624 fFieldPropagator->ClearPropagatorState();
625 // Resets all state of field propagator cl
626 // including safety values
627 // in case of overlaps and to wipe for fir
628
629 // G4ChordFinder* chordF= fFieldPropagator
630 // if( chordF ) chordF->ResetStepEstimate(
631 }
632
633 // Make sure to clear the chord finders of a
634 G4FieldManagerStore* fieldMgrStore = G4Field
635 fieldMgrStore->ClearAllChordFindersState();
636
637 // Update the current touchable handle (fro
638 //
639 fCurrentTouchableHandle = aTrack->GetTouchab
640 }
641
642 //....oooOO0OOooo........oooOO0OOooo........oo
643