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