source: trunk/source/error_propagation/src/G4ErrorPropagator.cc @ 1058

Last change on this file since 1058 was 1058, checked in by garnier, 15 years ago

file release beta
File size: 18.5 KB
Line 
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,v 1.7 2007/11/14 17:01:14 gcosmo Exp $
27// GEANT4 tag $Name: geant4-09-02-ref-02 $
28//
29// ------------------------------------------------------------
30//      GEANT 4 class implementation file
31// ------------------------------------------------------------
32//
33
34#include "G4ErrorPropagator.hh"
35#include "G4ErrorPropagatorData.hh"
36#include "G4ErrorFreeTrajState.hh"
37#include "G4ErrorSurfaceTrajState.hh"
38#include "G4ErrorGeomVolumeTarget.hh"
39#include "G4ErrorSurfaceTarget.hh"
40
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//---------------------------------------------------------------------------
57G4ErrorPropagator::G4ErrorPropagator()
58{
59  verbose =  G4ErrorPropagatorData::verbose();
60#ifdef G4EVERBOSE
61   if(verbose >= 5) { G4cout << "G4ErrorPropagator " << this << G4endl; }
62#endif
63
64  theG4Track = 0;
65  fpSteppingManager = G4EventManager::GetEventManager()
66                    ->GetTrackingManager()->GetSteppingManager();
67  thePropIsInitialized = false;
68}
69
70
71//-----------------------------------------------------------------------
72G4int G4ErrorPropagator::Propagate( G4ErrorTrajState* currentTS,
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 =
80    G4ErrorPropagatorData::GetErrorPropagatorData();
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  return ierr;
173}
174
175
176//-----------------------------------------------------------------------
177G4int G4ErrorPropagator::PropagateOneStep( G4ErrorTrajState* currentTS )
178{
179  G4ErrorPropagatorData* g4edata =
180    G4ErrorPropagatorData::GetErrorPropagatorData();
181
182  if ( (g4edata->GetState() == G4ErrorState_PreInit)
183    || (G4StateManager::GetStateManager()->GetCurrentState()
184       != G4State_GeomClosed) )
185  {
186    G4cout << "ERROR - G4ErrorPropagator::PropagateOneStep()" << G4endl
187           << "        Called before initialization is done for this track."
188           << G4endl
189           << "        Please call G4ErrorPropagatorManager::InitGeant4e()."
190           << G4endl;
191    G4Exception("G4ErrorPropagator::PropagateOneStep()",
192                "InvalidCall", FatalException,
193                "Called before initialization is done for this track!");
194  }
195
196  // to start ierror is set to 0 (= OK)
197  //
198  G4int ierr = 0;
199
200  //--- Do not propagate zero or too low energy particles
201  //
202  if( currentTS->GetMomentum().mag() < 1.E-9*MeV )
203  {
204    G4cerr << "ERROR - G4ErrorPropagator::PropagateOneStep()" << G4endl
205           << "        Energy too low to be propagated: "
206           << G4BestUnit(currentTS->GetMomentum().mag(),"Energy") << G4endl;
207    return -3;   
208  }
209
210#ifdef G4EVERBOSE
211  if( verbose >= 1 )
212  {
213    G4cout << " =====> starting GEANT4E tracking for "
214           << currentTS->GetParticleType()
215           << "  Forwards= " << g4edata->GetMode() << G4endl;
216  }
217  if( verbose >= 3 )
218  {
219    G4cout << " G4ErrorPropagator::Propagate initialTS ";
220    G4cout << *currentTS << G4endl;
221  }
222#endif
223
224  //----- If it is the first step, create a track
225  //
226  if( theStepN == 0 ) { theG4Track = InitG4Track( *currentTS ); }
227    // set to 0 by the initialization in G4ErrorPropagatorManager
228  theStepN++;
229
230  //----- Create a G4ErrorFreeTrajState
231  //
232  G4ErrorFreeTrajState* currentTS_FREE = InitFreeTrajState( currentTS );
233
234  //----- Track the particle one step
235  //
236  ierr = MakeOneStep( currentTS_FREE );
237
238  //----- Get the state on target
239  //
240  GetFinalTrajState( currentTS, currentTS_FREE, g4edata->GetTarget() );
241
242  return ierr;
243}
244
245
246//---------------------------------------------------------------------------
247G4Track* G4ErrorPropagator::InitG4Track( G4ErrorTrajState& initialTS )
248{
249  if( verbose >= 5 ) { G4cout << "InitG4Track " << G4endl; }
250
251  //----- Create Particle
252  //
253  const G4String partType = initialTS.GetParticleType();
254  G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
255  G4ParticleDefinition* particle = particleTable->FindParticle(partType); 
256  if( particle == 0)
257  {
258    G4cerr << "ERROR - G4ErrorPropagator::InitG4Track()" << G4endl
259           << "        Particle type not defined " + partType << G4endl;
260    G4Exception( "G4ErrorPropagator::InitG4Track()", "InvalidSetup",
261                 FatalException, "Particle type not defined !" );
262  }
263 
264  G4DynamicParticle* DP = 
265    new G4DynamicParticle(particle,initialTS.GetMomentum());
266
267  DP->SetPolarization(0.,0.,0.);
268
269  // Set Charge
270  //
271  if( particle->GetPDGCharge() < 0 )
272  {
273    DP->SetCharge(-eplus);
274  }
275  else
276  {
277    DP->SetCharge(eplus);
278  }
279
280  //----- Create Track
281  //
282  theG4Track = new G4Track(DP, 0., initialTS.GetPosition() );
283  theG4Track->SetParentID(0);
284
285#ifdef G4EVERBOSE
286  if(verbose >= 3)
287  {
288    G4cout << " G4ErrorPropagator new track of energy: "
289           << theG4Track->GetKineticEnergy() << G4endl;
290  }
291#endif
292 
293  //---- Reproduce G4TrackingManager::ProcessOneTrack initialization
294  InvokePreUserTrackingAction( theG4Track ); 
295
296  if( fpSteppingManager == 0 )
297  {
298    G4Exception("G4ErrorPropagator::InitG4Track()", "InvalidSetup",
299                FatalException, "G4SteppingManager not initialized yet!");
300  }
301  else
302  {
303    fpSteppingManager->SetInitialStep(theG4Track);
304  }
305
306  // Give SteppingManger the maximum number of processes
307  //
308  fpSteppingManager->GetProcessNumber();
309
310  // Give track the pointer to the Step
311  //
312  theG4Track->SetStep(fpSteppingManager->GetStep());
313
314  // Inform beginning of tracking to physics processes
315  //
316  theG4Track->GetDefinition()->GetProcessManager()->StartTracking(theG4Track);
317
318  initialTS.SetG4Track( theG4Track );
319
320  return theG4Track;
321}
322
323
324//-----------------------------------------------------------------------
325G4int G4ErrorPropagator::MakeSteps( G4ErrorFreeTrajState* currentTS_FREE )
326{
327  G4int ierr = 0;
328
329  //----- Track the particle Step-by-Step while it is alive
330  //
331  theStepLength = 0.;
332 
333  while( (theG4Track->GetTrackStatus() == fAlive) ||
334         (theG4Track->GetTrackStatus() == fStopButAlive) )
335  {
336    ierr = MakeOneStep( currentTS_FREE );
337    if( ierr != 0 ) { break; }
338
339    //----- Check if last step for error propagation
340    //
341    if( CheckIfLastStep( theG4Track ) )
342    {
343      break;
344    }
345  }
346  return ierr;
347}
348
349
350//-----------------------------------------------------------------------
351G4int G4ErrorPropagator::MakeOneStep( G4ErrorFreeTrajState* currentTS_FREE )
352{
353  G4ErrorPropagatorData* g4edata =
354    G4ErrorPropagatorData::GetErrorPropagatorData();
355  G4int ierr = 0;
356
357  //---------- Track one step
358#ifdef G4EVERBOSE
359  if(verbose >= 2 )
360  {
361    G4cout << G4endl
362           << "@@@@@@@@@@@@@@@@@@@@@@@@@ NEW STEP " << G4endl;
363  }
364#endif
365 
366  theG4Track->IncrementCurrentStepNumber();
367
368  fpSteppingManager->Stepping();
369 
370  //---------- Check if Target has been reached (and then set G4ErrorState)
371
372  // G4ErrorPropagationNavigator limits the step if target is closer than
373  // boundary (but the winner process is always "Transportation": then
374  // error propagator will stop the track
375
376  if( theG4Track->GetStep()->GetPostStepPoint()
377      ->GetProcessDefinedStep()->GetProcessName() == "Transportation" )
378  {
379    if( g4edata->GetState()
380       == G4ErrorState(G4ErrorState_TargetCloserThanBoundary) )
381    {  // target or step length reached
382     
383#ifdef G4EVERBOSE
384      if(verbose >= 5 )
385      {
386        G4cout << " transportation determined by geant4e " << G4endl;
387      }
388#endif     
389      g4edata->SetState( G4ErrorState_StoppedAtTarget );
390    }
391    else if( g4edata->GetTarget()->GetType() == G4ErrorTarget_GeomVolume )
392    {
393      G4ErrorGeomVolumeTarget* target =
394        (G4ErrorGeomVolumeTarget*)(g4edata->GetTarget());
395      if( static_cast<G4ErrorGeomVolumeTarget*>( target )
396          ->TargetReached( theG4Track->GetStep() ) )
397      {
398        g4edata->SetState( G4ErrorState_StoppedAtTarget ); 
399      } 
400    }
401  }
402  else if( theG4Track->GetStep()->GetPostStepPoint()->GetProcessDefinedStep()
403           ->GetProcessName() == "G4ErrorTrackLengthTarget" )
404  {
405    g4edata->SetState( G4ErrorState_StoppedAtTarget );
406  }
407
408  //---------- Propagate error 
409
410#ifdef G4EVERBOSE
411  if(verbose >= 2 )
412  {
413    G4cout << " propagating error " << *currentTS_FREE << G4endl;
414  }
415#endif
416  const G4Track* cTrack = const_cast<G4Track*>(theG4Track);
417  ierr = currentTS_FREE->PropagateError( cTrack );
418 
419#ifdef G4EVERBOSE
420  if(verbose >= 3 )
421  {
422    G4cout << " PropagateError returns " << ierr << G4endl;
423  }
424#endif
425
426  currentTS_FREE->Update( cTrack );
427 
428  theStepLength += theG4Track->GetStepLength();
429   
430  if(ierr != 0 )
431  {
432    G4cerr << "ERROR - G4ErrorPropagator:MakeOneStep()" << G4endl
433           << "        Error returned: " << ierr << G4endl
434           << "        Geant4 tracking will be stopped !" << G4endl;
435  }
436
437  return ierr; 
438}
439
440
441//-----------------------------------------------------------------------
442G4ErrorFreeTrajState*
443G4ErrorPropagator::InitFreeTrajState( G4ErrorTrajState* currentTS )
444{
445  G4ErrorFreeTrajState* currentTS_FREE = 0;
446
447  //----- Transform the TrajState to Free coordinates if it is OnSurface
448  //
449  if( currentTS->GetTSType() == G4eTS_OS )
450  {
451    G4ErrorSurfaceTrajState* tssd =
452      static_cast<G4ErrorSurfaceTrajState*>(currentTS);
453    currentTS_FREE = new G4ErrorFreeTrajState( *tssd );
454  }
455  else if( currentTS->GetTSType() == G4eTS_FREE )
456  {
457    currentTS_FREE = static_cast<G4ErrorFreeTrajState*>(currentTS);
458  }
459  else
460  {
461    G4cerr << "ERROR - G4ErrorPropagator::InitFreeTrajState()" << G4endl
462           << "WRONG TrajState " + currentTS->GetTSType() << G4endl;
463    G4Exception("G4ErrorPropagator::InitFreeTrajState()", "InvalidState",
464                FatalException, "WRONG trajectory state !");
465  }
466  return currentTS_FREE;
467}
468
469
470//-----------------------------------------------------------------------
471void G4ErrorPropagator::GetFinalTrajState( G4ErrorTrajState* currentTS,
472                                           G4ErrorFreeTrajState* currentTS_FREE,
473                                           const G4ErrorTarget* target ) 
474{
475  G4ErrorPropagatorData* g4edata =
476    G4ErrorPropagatorData::GetErrorPropagatorData();
477
478#ifdef G4EVERBOSE
479  if(verbose >= 1 )
480  {
481    G4cout << " G4ErrorPropagator::Propagate: final state "
482           << G4int(g4edata->GetState()) << " TSType "
483           << currentTS->GetTSType() << G4endl;
484  }
485#endif
486
487  if( (currentTS->GetTSType() == G4eTS_FREE) || 
488      (g4edata->GetState() != G4ErrorState_StoppedAtTarget) )
489  {
490    currentTS = currentTS_FREE;
491  }
492  else if( currentTS->GetTSType() == G4eTS_OS )
493  {
494    if( target->GetType() == G4ErrorTarget_TrkL )
495    {
496      G4Exception("G4ErrorPropagator:GetFinalTrajState()",
497                  "InvalidSetup", FatalException,
498                  "Using a G4ErrorSurfaceTrajState with wrong target");
499    }
500    const G4ErrorTanPlaneTarget* targetWTP =
501      static_cast<const G4ErrorTanPlaneTarget*>(target);
502    *currentTS = G4ErrorSurfaceTrajState(
503                 *(static_cast<G4ErrorFreeTrajState*>(currentTS_FREE)),
504                 targetWTP->GetTangentPlane( currentTS_FREE->GetPosition() ) );
505#ifdef G4EVERBOSE
506    if(verbose >= 1 )
507    {
508      G4cout << currentTS << " returning tssd " << *currentTS << G4endl;
509    }
510#endif
511    delete currentTS_FREE;
512  }
513}
514
515
516//-------------------------------------------------------------------------
517G4bool G4ErrorPropagator::CheckIfLastStep( G4Track* aTrack )
518{
519  G4bool exception = 0;
520  G4bool lastG4eStep = false;
521  G4ErrorPropagatorData* g4edata =
522    G4ErrorPropagatorData::GetErrorPropagatorData();
523
524#ifdef G4EVERBOSE
525  if( verbose >= 4 )
526  {
527    G4cout << " G4ErrorPropagator::CheckIfLastStep G4ErrorState= "
528           << G4int(g4edata->GetState()) << G4endl;
529  }
530#endif
531 
532  //----- Check if this is the last step: track has reached the target
533  //      or the end of world
534  //
535  if(g4edata->GetState() == G4ErrorState(G4ErrorState_StoppedAtTarget) )
536  {
537    lastG4eStep = true;   
538#ifdef G4EVERBOSE
539    if(verbose >= 5 )
540    {
541      G4cout << " G4ErrorPropagator::CheckIfLastStep " << lastG4eStep
542             << " " << G4int(g4edata->GetState()) << G4endl;
543    }
544#endif
545  }
546  else if( aTrack->GetNextVolume() == 0 )
547  {
548    //----- If particle is out of world, without finding the G4ErrorTarget,
549    //      give a n error/warning
550    //
551    lastG4eStep = true;
552    if( exception )
553    {
554      G4cerr << "ERROR - G4ErrorPropagator::CheckIfLastStep()" << G4endl
555             << "        Track extrapolated until end of World" << G4endl
556             << "        without finding the defined target " << G4endl;
557      G4Exception("G4ErrorPropagator::CheckIfLastStep()",
558                  "InvalidSetup", FatalException,
559                  "Track extrapolated without finding the defined target.");
560    }
561    else
562    {
563      if( verbose >= 1 )
564      {
565        G4cerr << "ERROR - G4ErrorPropagator::CheckIfLastStep()" << G4endl
566               << "        Track extrapolated until end of World" << G4endl
567               << "        without finding the defined target " << G4endl;
568      }
569    }
570  }  //----- not last step from G4e, but track is stopped (energy exhausted)
571  else if( aTrack->GetTrackStatus() == fStopAndKill )
572  { 
573    if( exception )
574    {
575      G4cerr << "ERROR - G4ErrorPropagator::CheckIfLastStep()" << G4endl
576             << "        Track extrapolated until energy is exhausted" << G4endl
577             << "        without finding the defined target !" << G4endl;
578      G4Exception("G4ErrorPropagator::CheckIfLastStep()",
579                  "InvalidSetup", FatalException,
580                  "Track extrapolated without finding the defined target.");
581    }
582    else
583    {
584      if( verbose >= 1 )
585      {
586        G4cerr << "ERROR - G4ErrorPropagator::CheckIfLastStep()" << G4endl
587             << "        Track extrapolated until energy is exhausted" << G4endl
588             << "        without finding the defined target !" << G4endl;
589      }
590      lastG4eStep = 1;
591    }
592  }
593
594#ifdef G4EVERBOSE
595  if( verbose >= 5 )
596  {
597    G4cout << " return CheckIfLastStep " << lastG4eStep << G4endl;
598  }
599#endif
600
601  return  lastG4eStep;
602}
603
604
605//---------------------------------------------------------------------------
606void 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//---------------------------------------------------------------------------
619void G4ErrorPropagator::InvokePostUserTrackingAction( G4Track* fpTrack )
620{
621  const G4UserTrackingAction* fpUserTrackingAction =
622    G4EventManager::GetEventManager()->GetUserTrackingAction();
623  if( fpUserTrackingAction != 0 )
624  {
625    const_cast<G4UserTrackingAction*>(fpUserTrackingAction)
626      ->PostUserTrackingAction((fpTrack) );
627  }
628}
Note: See TracBrowser for help on using the repository browser.