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

Last change on this file since 1296 was 1228, checked in by garnier, 16 years ago

update geant4.9.3 tag
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-03 $
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.