| Geant4 Cross Reference |
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 //
26 //
27 // ------------------------------------------------------------
28 // GEANT 4 class implementation file
29 // ------------------------------------------------------------
30 //
31
32 #include "G4ErrorPropagator.hh"
33 #include "G4ErrorPropagatorData.hh"
34 #include "G4ErrorFreeTrajState.hh"
35 #include "G4ErrorSurfaceTrajState.hh"
36 #include "G4ErrorGeomVolumeTarget.hh"
37 #include "G4ErrorSurfaceTarget.hh"
38
39 #include "G4SystemOfUnits.hh"
40 #include "G4DynamicParticle.hh"
41 #include "G4Track.hh"
42 #include "G4SteppingManager.hh"
43 #include "G4EventManager.hh"
44 #include "G4TrackingManager.hh"
45 #include "G4ParticleTable.hh"
46 #include "G4StateManager.hh"
47
48 #include "G4VPhysicalVolume.hh"
49 #include "G4PhysicalVolumeStore.hh"
50 #include "G4UnitsTable.hh"
51
52 #include <vector>
53
54 //---------------------------------------------------------------------------
55 G4ErrorPropagator::G4ErrorPropagator()
56 : theStepLength(0.)
57 , theInitialTrajState(0)
58 , theStepN(0)
59 , theG4Track(0)
60 {
61 verbose = G4ErrorPropagatorData::verbose();
62 #ifdef G4EVERBOSE
63 if(verbose >= 5)
64 {
65 G4cout << "G4ErrorPropagator " << this << G4endl;
66 }
67 #endif
68
69 fpSteppingManager = G4EventManager::GetEventManager()
70 ->GetTrackingManager()
71 ->GetSteppingManager();
72 thePropIsInitialized = false;
73 }
74
75 //-----------------------------------------------------------------------
76 G4int G4ErrorPropagator::Propagate(G4ErrorTrajState* currentTS,
77 const G4ErrorTarget* target,
78 G4ErrorMode mode)
79 {
80 // to start ierror is set to 1 (= OK)
81 //
82 G4int ierr = 1;
83
84 G4ErrorPropagatorData* g4edata =
85 G4ErrorPropagatorData::GetErrorPropagatorData();
86
87 //----- Do not propagate zero or too low energy particles
88 //
89 if(currentTS->GetMomentum().mag() < 1.E-9 * MeV)
90 {
91 std::ostringstream message;
92 message << "Energy too low to be propagated: "
93 << G4BestUnit(currentTS->GetMomentum().mag(), "Energy");
94 G4Exception("G4ErrorPropagator::Propagate()", "GEANT4e-Notification",
95 JustWarning, message);
96 return -3;
97 }
98
99 g4edata->SetMode(mode);
100
101 #ifdef G4EVERBOSE
102 if(verbose >= 1)
103 {
104 G4cout << " =====> starting GEANT4E tracking for "
105 << currentTS->GetParticleType()
106 << " Forwards= " << g4edata->GetMode() << G4endl;
107 }
108 if(verbose >= 1)
109 {
110 G4cout << G4endl << "@@@@@@@@@@@@@@@@@@@@@@@@@ NEW STEP " << G4endl;
111 }
112
113 if(verbose >= 3)
114 {
115 G4cout << " G4ErrorPropagator::Propagate initialTS ";
116 G4cout << *currentTS << G4endl;
117 target->Dump(G4String(" to target "));
118 }
119 #endif
120
121 g4edata->SetTarget(target);
122
123 //----- Create a track
124 //
125 if(theG4Track != 0)
126 {
127 delete theG4Track;
128 }
129 theG4Track = InitG4Track(*currentTS);
130
131 //----- Create a G4ErrorFreeTrajState
132 //
133 G4ErrorFreeTrajState* currentTS_FREE = InitFreeTrajState(currentTS);
134
135 //----- Track the particle
136 //
137 ierr = MakeSteps(currentTS_FREE);
138
139 //------ Tracking ended, check if target has been reached
140 // if target not found
141 //
142 if(g4edata->GetState() != G4ErrorState_StoppedAtTarget)
143 {
144 if(theG4Track->GetKineticEnergy() > 0.)
145 {
146 ierr = -ierr - 10;
147 }
148 else
149 {
150 ierr = -ierr - 20;
151 }
152 *currentTS = *currentTS_FREE;
153 if(verbose >= 0)
154 {
155 std::ostringstream message;
156 message << "Particle does not reach target: " << *currentTS;
157 G4Exception("G4ErrorPropagator::Propagate()", "GEANT4e-Notification",
158 JustWarning, message);
159 }
160 }
161 else
162 {
163 GetFinalTrajState(currentTS, currentTS_FREE, target);
164 }
165
166 #ifdef G4EVERBOSE
167 if(verbose >= 1)
168 {
169 G4cout << " G4ErrorPropagator: propagation ended " << G4endl;
170 }
171 if(verbose >= 2)
172 {
173 G4cout << " Current TrajState " << currentTS << G4endl;
174 }
175 #endif
176
177 // Inform end of tracking to physics processes
178 //
179 theG4Track->GetDefinition()->GetProcessManager()->EndTracking();
180
181 InvokePostUserTrackingAction(theG4Track);
182
183 // delete currentTS_FREE;
184
185 return ierr;
186 }
187
188 //-----------------------------------------------------------------------
189 G4int G4ErrorPropagator::PropagateOneStep(G4ErrorTrajState* currentTS)
190 {
191 G4ErrorPropagatorData* g4edata =
192 G4ErrorPropagatorData::GetErrorPropagatorData();
193
194 if((g4edata->GetState() == G4ErrorState_PreInit) ||
195 (G4StateManager::GetStateManager()->GetCurrentState() !=
196 G4State_GeomClosed))
197 {
198 std::ostringstream message;
199 message << "Called before initialization is done for this track!";
200 G4Exception("G4ErrorPropagator::PropagateOneStep()", "InvalidCall",
201 FatalException, message,
202 "Please call G4ErrorPropagatorManager::InitGeant4e().");
203 }
204
205 // to start ierror is set to 0 (= OK)
206 //
207 G4int ierr = 0;
208
209 //--- Do not propagate zero or too low energy particles
210 //
211 if(currentTS->GetMomentum().mag() < 1.E-9 * MeV)
212 {
213 std::ostringstream message;
214 message << "Energy too low to be propagated: "
215 << G4BestUnit(currentTS->GetMomentum().mag(), "Energy");
216 G4Exception("G4ErrorPropagator::PropagateOneStep()", "GEANT4e-Notification",
217 JustWarning, message);
218 return -3;
219 }
220
221 #ifdef G4EVERBOSE
222 if(verbose >= 1)
223 {
224 G4cout << " =====> starting GEANT4E tracking for "
225 << currentTS->GetParticleType()
226 << " Forwards= " << g4edata->GetMode() << G4endl;
227 }
228 if(verbose >= 3)
229 {
230 G4cout << " G4ErrorPropagator::Propagate initialTS ";
231 G4cout << *currentTS << G4endl;
232 }
233 #endif
234
235 //----- If it is the first step, create a track
236 //
237 if(theStepN == 0)
238 {
239 if(theG4Track != 0)
240 {
241 delete theG4Track;
242 }
243 theG4Track = InitG4Track(*currentTS);
244 }
245 // set to 0 by the initialization in G4ErrorPropagatorManager
246 theStepN++;
247
248 //----- Create a G4ErrorFreeTrajState
249 //
250 G4ErrorFreeTrajState* currentTS_FREE = InitFreeTrajState(currentTS);
251
252 //----- Track the particle one step
253 //
254 ierr = MakeOneStep(currentTS_FREE);
255
256 //----- Get the state on target
257 //
258 GetFinalTrajState(currentTS, currentTS_FREE, g4edata->GetTarget());
259
260 return ierr;
261 }
262
263 //---------------------------------------------------------------------------
264 G4Track* G4ErrorPropagator::InitG4Track(G4ErrorTrajState& initialTS)
265 {
266 if(verbose >= 5)
267 {
268 G4cout << "InitG4Track " << G4endl;
269 }
270
271 //----- Create Particle
272 //
273 const G4String partType = initialTS.GetParticleType();
274 G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
275 G4ParticleDefinition* particle = particleTable->FindParticle(partType);
276 if(particle == 0)
277 {
278 std::ostringstream message;
279 message << "Particle type not defined: " << partType;
280 G4Exception("G4ErrorPropagator::InitG4Track()", "InvalidSetup",
281 FatalException, message);
282 }
283
284 G4DynamicParticle* DP =
285 new G4DynamicParticle(particle, initialTS.GetMomentum());
286
287 DP->SetPolarization(0., 0., 0.);
288
289 // Set Charge
290 //
291 if(particle->GetPDGCharge() < 0)
292 {
293 DP->SetCharge(-eplus);
294 }
295 else
296 {
297 DP->SetCharge(eplus);
298 }
299
300 //----- Create Track
301 //
302 theG4Track = new G4Track(DP, 0., initialTS.GetPosition());
303 theG4Track->SetParentID(0);
304
305 #ifdef G4EVERBOSE
306 if(verbose >= 3)
307 {
308 G4cout << " G4ErrorPropagator new track of energy: "
309 << theG4Track->GetKineticEnergy() << G4endl;
310 }
311 #endif
312
313 //---- Reproduce G4TrackingManager::ProcessOneTrack initialization
314 InvokePreUserTrackingAction(theG4Track);
315
316 if(fpSteppingManager == 0)
317 {
318 G4Exception("G4ErrorPropagator::InitG4Track()", "InvalidSetup",
319 FatalException, "G4SteppingManager not initialized yet!");
320 }
321 else
322 {
323 fpSteppingManager->SetInitialStep(theG4Track);
324 }
325
326 // Give SteppingManger the maximum number of processes
327 //
328 fpSteppingManager->GetProcessNumber();
329
330 // Give track the pointer to the Step
331 //
332 theG4Track->SetStep(fpSteppingManager->GetStep());
333
334 // Inform beginning of tracking to physics processes
335 //
336 theG4Track->GetDefinition()->GetProcessManager()->StartTracking(theG4Track);
337
338 initialTS.SetG4Track(theG4Track);
339
340 return theG4Track;
341 }
342
343 //-----------------------------------------------------------------------
344 G4int G4ErrorPropagator::MakeSteps(G4ErrorFreeTrajState* currentTS_FREE)
345 {
346 G4int ierr = 0;
347
348 //----- Track the particle Step-by-Step while it is alive
349 //
350 theStepLength = 0.;
351
352 while((theG4Track->GetTrackStatus() == fAlive) ||
353 (theG4Track->GetTrackStatus() == fStopButAlive))
354 {
355 ierr = MakeOneStep(currentTS_FREE);
356 if(ierr != 0)
357 {
358 break;
359 }
360
361 //----- Check if last step for error propagation
362 //
363 if(CheckIfLastStep(theG4Track))
364 {
365 break;
366 }
367 } // Loop checking, 06.08.2015, G.Cosmo
368 return ierr;
369 }
370
371 //-----------------------------------------------------------------------
372 G4int G4ErrorPropagator::MakeOneStep(G4ErrorFreeTrajState* currentTS_FREE)
373 {
374 G4ErrorPropagatorData* g4edata =
375 G4ErrorPropagatorData::GetErrorPropagatorData();
376 G4int ierr = 0;
377
378 //---------- Track one step
379 #ifdef G4EVERBOSE
380 if(verbose >= 2)
381 {
382 G4cout << G4endl << "@@@@@@@@@@@@@@@@@@@@@@@@@ NEW STEP " << G4endl;
383 }
384 #endif
385
386 theG4Track->IncrementCurrentStepNumber();
387
388 fpSteppingManager->Stepping();
389
390 //---------- Check if Target has been reached (and then set G4ErrorState)
391
392 // G4ErrorPropagationNavigator limits the step if target is closer than
393 // boundary (but the winner process is always "Transportation": then
394 // error propagator will stop the track
395
396 if(theG4Track->GetStep()
397 ->GetPostStepPoint()
398 ->GetProcessDefinedStep()
399 ->GetProcessName() == "Transportation")
400 {
401 if(g4edata->GetState() ==
402 G4ErrorState(G4ErrorState_TargetCloserThanBoundary))
403 { // target or step length reached
404
405 #ifdef G4EVERBOSE
406 if(verbose >= 5)
407 {
408 G4cout << " transportation determined by geant4e " << G4endl;
409 }
410 #endif
411 g4edata->SetState(G4ErrorState_StoppedAtTarget);
412 }
413 else if(g4edata->GetTarget()->GetType() == G4ErrorTarget_GeomVolume)
414 {
415 G4ErrorGeomVolumeTarget* target =
416 (G4ErrorGeomVolumeTarget*) (g4edata->GetTarget());
417 if(static_cast<G4ErrorGeomVolumeTarget*>(target)->TargetReached(
418 theG4Track->GetStep()))
419 {
420 g4edata->SetState(G4ErrorState_StoppedAtTarget);
421 }
422 }
423 }
424 else if(theG4Track->GetStep()
425 ->GetPostStepPoint()
426 ->GetProcessDefinedStep()
427 ->GetProcessName() == "G4ErrorTrackLengthTarget")
428 {
429 g4edata->SetState(G4ErrorState_StoppedAtTarget);
430 }
431
432 //---------- Propagate error
433
434 #ifdef G4EVERBOSE
435 if(verbose >= 2)
436 {
437 G4cout << " propagating error " << *currentTS_FREE << G4endl;
438 }
439 #endif
440 const G4Track* cTrack = const_cast<G4Track*>(theG4Track);
441 ierr = currentTS_FREE->PropagateError(cTrack);
442
443 #ifdef G4EVERBOSE
444 if(verbose >= 3)
445 {
446 G4cout << " PropagateError returns " << ierr << G4endl;
447 }
448 #endif
449
450 currentTS_FREE->Update(cTrack);
451
452 theStepLength += theG4Track->GetStepLength();
453
454 if(ierr != 0)
455 {
456 std::ostringstream message;
457 message << "Error returned: " << ierr;
458 G4Exception("G4ErrorPropagator::MakeOneStep()", "GEANT4e-Notification",
459 JustWarning, message, "Geant4 tracking will be stopped !");
460 }
461
462 return ierr;
463 }
464
465 //-----------------------------------------------------------------------
466 G4ErrorFreeTrajState* G4ErrorPropagator::InitFreeTrajState(
467 G4ErrorTrajState* currentTS)
468 {
469 G4ErrorFreeTrajState* currentTS_FREE = 0;
470
471 //----- Transform the TrajState to Free coordinates if it is OnSurface
472 //
473 if(currentTS->GetTSType() == G4eTS_OS)
474 {
475 G4ErrorSurfaceTrajState* tssd =
476 static_cast<G4ErrorSurfaceTrajState*>(currentTS);
477 currentTS_FREE = new G4ErrorFreeTrajState(*tssd);
478 }
479 else if(currentTS->GetTSType() == G4eTS_FREE)
480 {
481 currentTS_FREE = static_cast<G4ErrorFreeTrajState*>(currentTS);
482 }
483 else
484 {
485 std::ostringstream message;
486 message << "Wrong trajectory state: " << currentTS->GetTSType();
487 G4Exception("G4ErrorPropagator::InitFreeTrajState()", "InvalidState",
488 FatalException, message);
489 }
490 return currentTS_FREE;
491 }
492
493 //-----------------------------------------------------------------------
494 void G4ErrorPropagator::GetFinalTrajState(G4ErrorTrajState* currentTS,
495 G4ErrorFreeTrajState* currentTS_FREE,
496 const G4ErrorTarget* target)
497 {
498 G4ErrorPropagatorData* g4edata =
499 G4ErrorPropagatorData::GetErrorPropagatorData();
500
501 #ifdef G4EVERBOSE
502 if(verbose >= 1)
503 {
504 G4cout << " G4ErrorPropagator::Propagate: final state "
505 << G4int(g4edata->GetState()) << " TSType " << currentTS->GetTSType()
506 << G4endl;
507 }
508 #endif
509
510 if((currentTS->GetTSType() == G4eTS_FREE) ||
511 (g4edata->GetState() != G4ErrorState_StoppedAtTarget))
512 {
513 currentTS = currentTS_FREE;
514 }
515 else if(currentTS->GetTSType() == G4eTS_OS)
516 {
517 if(target->GetType() == G4ErrorTarget_TrkL)
518 {
519 G4Exception("G4ErrorPropagator:GetFinalTrajState()", "InvalidSetup",
520 FatalException,
521 "Using a G4ErrorSurfaceTrajState with wrong target");
522 }
523 const G4ErrorTanPlaneTarget* targetWTP =
524 static_cast<const G4ErrorTanPlaneTarget*>(target);
525 *currentTS = G4ErrorSurfaceTrajState(
526 *(static_cast<G4ErrorFreeTrajState*>(currentTS_FREE)),
527 targetWTP->GetTangentPlane(currentTS_FREE->GetPosition()));
528 #ifdef G4EVERBOSE
529 if(verbose >= 1)
530 {
531 G4cout << currentTS << " returning tssd " << *currentTS << G4endl;
532 }
533 #endif
534 delete currentTS_FREE;
535 }
536 }
537
538 //-------------------------------------------------------------------------
539 G4bool G4ErrorPropagator::CheckIfLastStep(G4Track* aTrack)
540 {
541 G4bool lastG4eStep = false;
542 G4ErrorPropagatorData* g4edata =
543 G4ErrorPropagatorData::GetErrorPropagatorData();
544
545 #ifdef G4EVERBOSE
546 if(verbose >= 4)
547 {
548 G4cout << " G4ErrorPropagator::CheckIfLastStep G4ErrorState= "
549 << G4int(g4edata->GetState()) << G4endl;
550 }
551 #endif
552
553 //----- Check if this is the last step: track has reached the target
554 // or the end of world
555 //
556 if(g4edata->GetState() == G4ErrorState(G4ErrorState_StoppedAtTarget))
557 {
558 lastG4eStep = true;
559 #ifdef G4EVERBOSE
560 if(verbose >= 5)
561 {
562 G4cout << " G4ErrorPropagator::CheckIfLastStep " << lastG4eStep << " "
563 << G4int(g4edata->GetState()) << G4endl;
564 }
565 #endif
566 }
567 else if(aTrack->GetNextVolume() == 0)
568 {
569 //----- If particle is out of world, without finding the G4ErrorTarget,
570 // give a n error/warning
571 //
572 lastG4eStep = true;
573 if(verbose >= 1)
574 {
575 std::ostringstream message;
576 message << "Track extrapolated until end of World" << G4endl
577 << "without finding the defined target.";
578 G4Exception("G4ErrorPropagator::CheckIfLastStep()",
579 "GEANT4e-Notification", JustWarning, message);
580 }
581 } //----- not last step from G4e, but track is stopped (energy exhausted)
582 else if(aTrack->GetTrackStatus() == fStopAndKill)
583 {
584 if(verbose >= 1)
585 {
586 std::ostringstream message;
587 message << "Track extrapolated until energy is exhausted" << G4endl
588 << "without finding the defined target.";
589 G4Exception("G4ErrorPropagator::CheckIfLastStep()",
590 "GEANT4e-Notification", JustWarning, message);
591 }
592 lastG4eStep = 1;
593 }
594
595 #ifdef G4EVERBOSE
596 if(verbose >= 5)
597 {
598 G4cout << " return CheckIfLastStep " << lastG4eStep << G4endl;
599 }
600 #endif
601
602 return lastG4eStep;
603 }
604
605 //---------------------------------------------------------------------------
606 void G4ErrorPropagator::InvokePreUserTrackingAction(G4Track* fpTrack)
607 {
608 const G4UserTrackingAction* fpUserTrackingAction =
609 G4EventManager::GetEventManager()->GetUserTrackingAction();
610 if(fpUserTrackingAction != 0)
611 {
612 const_cast<G4UserTrackingAction*>(fpUserTrackingAction)
613 ->PreUserTrackingAction((fpTrack));
614 }
615 }
616
617 //---------------------------------------------------------------------------
618 void G4ErrorPropagator::InvokePostUserTrackingAction(G4Track* fpTrack)
619 {
620 const G4UserTrackingAction* fpUserTrackingAction =
621 G4EventManager::GetEventManager()->GetUserTrackingAction();
622 if(fpUserTrackingAction != 0)
623 {
624 const_cast<G4UserTrackingAction*>(fpUserTrackingAction)
625 ->PostUserTrackingAction((fpTrack));
626 }
627 }
628