G4ErrorPropagator.cc

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 // $Id: G4ErrorPropagator.cc 69834 2013-05-16 08:03:15Z gcosmo $
00027 //
00028 // ------------------------------------------------------------
00029 //      GEANT 4 class implementation file 
00030 // ------------------------------------------------------------
00031 //
00032 
00033 #include "G4ErrorPropagator.hh"
00034 #include "G4ErrorPropagatorData.hh"
00035 #include "G4ErrorFreeTrajState.hh"
00036 #include "G4ErrorSurfaceTrajState.hh"
00037 #include "G4ErrorGeomVolumeTarget.hh"
00038 #include "G4ErrorSurfaceTarget.hh"
00039 
00040 #include "G4SystemOfUnits.hh"
00041 #include "G4DynamicParticle.hh"
00042 #include "G4Track.hh"
00043 #include "G4SteppingManager.hh"
00044 #include "G4EventManager.hh"
00045 #include "G4TrackingManager.hh"
00046 #include "G4ParticleTable.hh"
00047 #include "G4StateManager.hh"
00048 
00049 #include "G4VPhysicalVolume.hh"
00050 #include "G4PhysicalVolumeStore.hh"
00051 #include "G4UnitsTable.hh"
00052 
00053 #include <vector>
00054 
00055 
00056 //---------------------------------------------------------------------------
00057 G4ErrorPropagator::G4ErrorPropagator()
00058   : theStepLength(0.), theInitialTrajState(0), theStepN(0), theG4Track(0)
00059 {
00060   verbose =  G4ErrorPropagatorData::verbose();
00061 #ifdef G4EVERBOSE
00062    if(verbose >= 5) { G4cout << "G4ErrorPropagator " << this << G4endl; }
00063 #endif
00064 
00065   fpSteppingManager = G4EventManager::GetEventManager()
00066                     ->GetTrackingManager()->GetSteppingManager();
00067   thePropIsInitialized = false;
00068 }
00069 
00070 
00071 //-----------------------------------------------------------------------
00072 G4int G4ErrorPropagator::Propagate( G4ErrorTrajState* currentTS,
00073                               const G4ErrorTarget* target, G4ErrorMode mode )
00074 {
00075   // to start ierror is set to 1 (= OK)
00076   //
00077   G4int ierr = 1;
00078 
00079   G4ErrorPropagatorData* g4edata =
00080     G4ErrorPropagatorData::GetErrorPropagatorData();
00081 
00082   //----- Do not propagate zero or too low energy particles
00083   //
00084   if( currentTS->GetMomentum().mag() < 1.E-9*MeV )
00085   {
00086     G4cerr << "ERROR - G4ErrorPropagator::Propagate()" << G4endl
00087            << "        Energy too low to be propagated: "
00088            << G4BestUnit(currentTS->GetMomentum().mag(),"Energy") << G4endl;
00089     return -3; 
00090   }
00091 
00092   g4edata->SetMode( mode );
00093 
00094 #ifdef G4EVERBOSE
00095   if( verbose >= 1 )
00096   {
00097      G4cout << " =====> starting GEANT4E tracking for "
00098             << currentTS->GetParticleType()
00099             << "  Forwards= " << g4edata->GetMode() << G4endl;
00100   }
00101   if(verbose >= 1 )
00102   {
00103      G4cout << G4endl << "@@@@@@@@@@@@@@@@@@@@@@@@@ NEW STEP " << G4endl;
00104   }
00105 
00106   if( verbose >= 3 )
00107   {
00108     G4cout << " G4ErrorPropagator::Propagate initialTS ";
00109     G4cout << *currentTS << G4endl;
00110     target->Dump(G4String(" to target "));
00111   }
00112 #endif
00113 
00114   g4edata->SetTarget( target );
00115 
00116   //----- Create a track
00117   //
00118   if( theG4Track != 0 ) { delete theG4Track; }
00119   theG4Track = InitG4Track( *currentTS );
00120 
00121   //----- Create a G4ErrorFreeTrajState
00122   //
00123   G4ErrorFreeTrajState* currentTS_FREE = InitFreeTrajState( currentTS );
00124 
00125   //----- Track the particle
00126   //
00127   ierr = MakeSteps( currentTS_FREE );
00128 
00129   //------ Tracking ended, check if target has been reached
00130   //       if target not found
00131   //
00132   if( g4edata->GetState() != G4ErrorState_StoppedAtTarget )
00133   {
00134     if( theG4Track->GetKineticEnergy() > 0. )
00135     {
00136       ierr = -ierr - 10;
00137     }
00138     else
00139     {
00140       ierr = -ierr - 20;
00141     }
00142     *currentTS = *currentTS_FREE;
00143     if(verbose >= 0 )
00144     {
00145        G4cerr << "ERROR - G4ErrorPropagator::Propagate()" << G4endl
00146               << "        Particle does not reach target: " << *currentTS
00147               << G4endl;
00148     }
00149   }
00150   else
00151   {
00152     GetFinalTrajState( currentTS, currentTS_FREE, target );
00153   }
00154 
00155 #ifdef G4EVERBOSE
00156   if( verbose >= 1 )
00157   {
00158     G4cout << " G4ErrorPropagator: propagation ended " << G4endl;
00159   }
00160   if( verbose >= 2 )
00161   {
00162     G4cout << " Current TrajState " << currentTS << G4endl;
00163   }
00164 #endif
00165  
00166   // Inform end of tracking to physics processes
00167   //
00168   theG4Track->GetDefinition()->GetProcessManager()->EndTracking();
00169 
00170   InvokePostUserTrackingAction( theG4Track );
00171 
00172   // delete currentTS_FREE;
00173 
00174   return ierr;
00175 }
00176 
00177 
00178 //-----------------------------------------------------------------------
00179 G4int G4ErrorPropagator::PropagateOneStep( G4ErrorTrajState* currentTS )
00180 {
00181   G4ErrorPropagatorData* g4edata =
00182     G4ErrorPropagatorData::GetErrorPropagatorData();
00183 
00184   if ( (g4edata->GetState() == G4ErrorState_PreInit)
00185     || (G4StateManager::GetStateManager()->GetCurrentState()
00186        != G4State_GeomClosed) )
00187   {
00188     std::ostringstream message;
00189     message << "Called before initialization is done for this track!";
00190     G4Exception("G4ErrorPropagator::PropagateOneStep()",
00191                 "InvalidCall", FatalException, message,
00192                 "Please call G4ErrorPropagatorManager::InitGeant4e().");
00193   }
00194 
00195   // to start ierror is set to 0 (= OK)
00196   //
00197   G4int ierr = 0;
00198 
00199   //--- Do not propagate zero or too low energy particles
00200   //
00201   if( currentTS->GetMomentum().mag() < 1.E-9*MeV )
00202   {
00203     G4cerr << "ERROR - G4ErrorPropagator::PropagateOneStep()" << G4endl
00204            << "        Energy too low to be propagated: "
00205            << G4BestUnit(currentTS->GetMomentum().mag(),"Energy") << G4endl;
00206     return -3;   
00207   }
00208 
00209 #ifdef G4EVERBOSE
00210   if( verbose >= 1 )
00211   {
00212     G4cout << " =====> starting GEANT4E tracking for "
00213            << currentTS->GetParticleType()
00214            << "  Forwards= " << g4edata->GetMode() << G4endl;
00215   }
00216   if( verbose >= 3 )
00217   {
00218     G4cout << " G4ErrorPropagator::Propagate initialTS ";
00219     G4cout << *currentTS << G4endl;
00220   }
00221 #endif
00222 
00223   //----- If it is the first step, create a track
00224   //
00225   if( theStepN == 0 )
00226   {
00227     if( theG4Track != 0 ) { delete theG4Track; }
00228     theG4Track = InitG4Track( *currentTS );
00229   }
00230     // set to 0 by the initialization in G4ErrorPropagatorManager
00231   theStepN++;
00232 
00233   //----- Create a G4ErrorFreeTrajState
00234   //
00235   G4ErrorFreeTrajState* currentTS_FREE = InitFreeTrajState( currentTS );
00236 
00237   //----- Track the particle one step
00238   //
00239   ierr = MakeOneStep( currentTS_FREE );
00240 
00241   //----- Get the state on target
00242   //
00243   GetFinalTrajState( currentTS, currentTS_FREE, g4edata->GetTarget() );
00244 
00245   return ierr;
00246 }
00247 
00248 
00249 //---------------------------------------------------------------------------
00250 G4Track* G4ErrorPropagator::InitG4Track( G4ErrorTrajState& initialTS )
00251 {
00252   if( verbose >= 5 ) { G4cout << "InitG4Track " << G4endl; }
00253 
00254   //----- Create Particle
00255   //
00256   const G4String partType = initialTS.GetParticleType();
00257   G4ParticleTable* particleTable = G4ParticleTable::GetParticleTable();
00258   G4ParticleDefinition* particle = particleTable->FindParticle(partType); 
00259   if( particle == 0)
00260   {
00261     std::ostringstream message;
00262     message << "Particle type not defined: " << partType;
00263     G4Exception( "G4ErrorPropagator::InitG4Track()", "InvalidSetup",
00264                  FatalException, message );
00265   }
00266  
00267   G4DynamicParticle* DP = 
00268     new G4DynamicParticle(particle,initialTS.GetMomentum());
00269 
00270   DP->SetPolarization(0.,0.,0.);
00271 
00272   // Set Charge
00273   //
00274   if( particle->GetPDGCharge() < 0 )
00275   {
00276     DP->SetCharge(-eplus);
00277   }
00278   else
00279   {
00280     DP->SetCharge(eplus);
00281   }
00282 
00283   //----- Create Track 
00284   //
00285   theG4Track = new G4Track(DP, 0., initialTS.GetPosition() );
00286   theG4Track->SetParentID(0);
00287 
00288 #ifdef G4EVERBOSE
00289   if(verbose >= 3)
00290   {
00291     G4cout << " G4ErrorPropagator new track of energy: "
00292            << theG4Track->GetKineticEnergy() << G4endl;
00293   }
00294 #endif
00295   
00296   //---- Reproduce G4TrackingManager::ProcessOneTrack initialization
00297   InvokePreUserTrackingAction( theG4Track );  
00298 
00299   if( fpSteppingManager == 0 )
00300   {
00301     G4Exception("G4ErrorPropagator::InitG4Track()", "InvalidSetup",
00302                 FatalException, "G4SteppingManager not initialized yet!");
00303   }
00304   else
00305   {
00306     fpSteppingManager->SetInitialStep(theG4Track);
00307   }
00308 
00309   // Give SteppingManger the maximum number of processes
00310   //
00311   fpSteppingManager->GetProcessNumber();
00312 
00313   // Give track the pointer to the Step
00314   //
00315   theG4Track->SetStep(fpSteppingManager->GetStep());
00316 
00317   // Inform beginning of tracking to physics processes
00318   //
00319   theG4Track->GetDefinition()->GetProcessManager()->StartTracking(theG4Track);
00320 
00321   initialTS.SetG4Track( theG4Track );
00322 
00323   return theG4Track;
00324 }
00325 
00326 
00327 //-----------------------------------------------------------------------
00328 G4int G4ErrorPropagator::MakeSteps( G4ErrorFreeTrajState* currentTS_FREE )
00329 {
00330   G4int ierr = 0;
00331 
00332   //----- Track the particle Step-by-Step while it is alive
00333   //
00334   theStepLength = 0.;
00335   
00336   while( (theG4Track->GetTrackStatus() == fAlive) ||
00337          (theG4Track->GetTrackStatus() == fStopButAlive) )
00338   {
00339     ierr = MakeOneStep( currentTS_FREE );
00340     if( ierr != 0 ) { break; }
00341 
00342     //----- Check if last step for error propagation
00343     //
00344     if( CheckIfLastStep( theG4Track ) )
00345     {
00346       break;
00347     }
00348   }
00349   return ierr;
00350 }
00351 
00352 
00353 //-----------------------------------------------------------------------
00354 G4int G4ErrorPropagator::MakeOneStep( G4ErrorFreeTrajState* currentTS_FREE )
00355 {
00356   G4ErrorPropagatorData* g4edata =
00357     G4ErrorPropagatorData::GetErrorPropagatorData();
00358   G4int ierr = 0;
00359 
00360   //---------- Track one step
00361 #ifdef G4EVERBOSE
00362   if(verbose >= 2 )
00363   {
00364     G4cout << G4endl
00365            << "@@@@@@@@@@@@@@@@@@@@@@@@@ NEW STEP " << G4endl;
00366   }
00367 #endif
00368   
00369   theG4Track->IncrementCurrentStepNumber();
00370 
00371   fpSteppingManager->Stepping();
00372   
00373   //---------- Check if Target has been reached (and then set G4ErrorState)
00374 
00375   // G4ErrorPropagationNavigator limits the step if target is closer than
00376   // boundary (but the winner process is always "Transportation": then
00377   // error propagator will stop the track
00378 
00379   if( theG4Track->GetStep()->GetPostStepPoint()
00380       ->GetProcessDefinedStep()->GetProcessName() == "Transportation" )
00381   {
00382     if( g4edata->GetState()
00383        == G4ErrorState(G4ErrorState_TargetCloserThanBoundary) )
00384     {  // target or step length reached
00385       
00386 #ifdef G4EVERBOSE
00387       if(verbose >= 5 )
00388       {
00389         G4cout << " transportation determined by geant4e " << G4endl;
00390       }
00391 #endif      
00392       g4edata->SetState( G4ErrorState_StoppedAtTarget );
00393     }
00394     else if( g4edata->GetTarget()->GetType() == G4ErrorTarget_GeomVolume )
00395     {
00396       G4ErrorGeomVolumeTarget* target =
00397         (G4ErrorGeomVolumeTarget*)(g4edata->GetTarget());
00398       if( static_cast<G4ErrorGeomVolumeTarget*>( target )
00399           ->TargetReached( theG4Track->GetStep() ) )
00400       {
00401         g4edata->SetState( G4ErrorState_StoppedAtTarget ); 
00402       } 
00403     }
00404   }
00405   else if( theG4Track->GetStep()->GetPostStepPoint()->GetProcessDefinedStep()
00406            ->GetProcessName() == "G4ErrorTrackLengthTarget" )
00407   {
00408     g4edata->SetState( G4ErrorState_StoppedAtTarget );
00409   }
00410 
00411   //---------- Propagate error  
00412 
00413 #ifdef G4EVERBOSE
00414   if(verbose >= 2 )
00415   {
00416     G4cout << " propagating error " << *currentTS_FREE << G4endl;
00417   }
00418 #endif
00419   const G4Track* cTrack = const_cast<G4Track*>(theG4Track);
00420   ierr = currentTS_FREE->PropagateError( cTrack );
00421   
00422 #ifdef G4EVERBOSE
00423   if(verbose >= 3 )
00424   {
00425     G4cout << " PropagateError returns " << ierr << G4endl;
00426   }
00427 #endif
00428 
00429   currentTS_FREE->Update( cTrack );
00430   
00431   theStepLength += theG4Track->GetStepLength();
00432    
00433   if(ierr != 0 )
00434   {
00435     G4cerr << "ERROR - G4ErrorPropagator:MakeOneStep()" << G4endl
00436            << "        Error returned: " << ierr << G4endl
00437            << "        Geant4 tracking will be stopped !" << G4endl;
00438   }
00439 
00440   return ierr; 
00441 }
00442 
00443 
00444 //-----------------------------------------------------------------------
00445 G4ErrorFreeTrajState*
00446 G4ErrorPropagator::InitFreeTrajState( G4ErrorTrajState* currentTS )
00447 {
00448   G4ErrorFreeTrajState* currentTS_FREE = 0;
00449 
00450   //----- Transform the TrajState to Free coordinates if it is OnSurface
00451   //
00452   if( currentTS->GetTSType() == G4eTS_OS )
00453   {
00454     G4ErrorSurfaceTrajState* tssd =
00455       static_cast<G4ErrorSurfaceTrajState*>(currentTS);
00456     currentTS_FREE = new G4ErrorFreeTrajState( *tssd );
00457   }
00458   else if( currentTS->GetTSType() == G4eTS_FREE )
00459   {
00460     currentTS_FREE = static_cast<G4ErrorFreeTrajState*>(currentTS);
00461   }
00462   else
00463   {
00464     std::ostringstream message;
00465     message << "Wrong trajectory state: " << currentTS->GetTSType();
00466     G4Exception("G4ErrorPropagator::InitFreeTrajState()", "InvalidState",
00467                 FatalException, message);
00468   }
00469   return currentTS_FREE;
00470 }
00471 
00472 
00473 //-----------------------------------------------------------------------
00474 void G4ErrorPropagator::GetFinalTrajState( G4ErrorTrajState* currentTS,
00475                                            G4ErrorFreeTrajState* currentTS_FREE,
00476                                            const G4ErrorTarget* target ) 
00477 {
00478   G4ErrorPropagatorData* g4edata =
00479     G4ErrorPropagatorData::GetErrorPropagatorData();
00480 
00481 #ifdef G4EVERBOSE
00482   if(verbose >= 1 )
00483   {
00484     G4cout << " G4ErrorPropagator::Propagate: final state "
00485            << G4int(g4edata->GetState()) << " TSType "
00486            << currentTS->GetTSType() << G4endl;
00487   }
00488 #endif
00489 
00490   if( (currentTS->GetTSType() == G4eTS_FREE) || 
00491       (g4edata->GetState() != G4ErrorState_StoppedAtTarget) )
00492   {
00493     currentTS = currentTS_FREE;
00494   }
00495   else if( currentTS->GetTSType() == G4eTS_OS )
00496   {
00497     if( target->GetType() == G4ErrorTarget_TrkL )
00498     {
00499       G4Exception("G4ErrorPropagator:GetFinalTrajState()",
00500                   "InvalidSetup", FatalException,
00501                   "Using a G4ErrorSurfaceTrajState with wrong target");
00502     }
00503     const G4ErrorTanPlaneTarget* targetWTP =
00504       static_cast<const G4ErrorTanPlaneTarget*>(target);
00505     *currentTS = G4ErrorSurfaceTrajState(
00506                  *(static_cast<G4ErrorFreeTrajState*>(currentTS_FREE)),
00507                  targetWTP->GetTangentPlane( currentTS_FREE->GetPosition() ) );
00508 #ifdef G4EVERBOSE
00509     if(verbose >= 1 )
00510     {
00511       G4cout << currentTS << " returning tssd " << *currentTS << G4endl;
00512     }
00513 #endif
00514     delete currentTS_FREE;
00515   }
00516 }
00517 
00518 
00519 //-------------------------------------------------------------------------
00520 G4bool G4ErrorPropagator::CheckIfLastStep( G4Track* aTrack )
00521 {
00522   G4bool exception = 0;
00523   G4bool lastG4eStep = false;
00524   G4ErrorPropagatorData* g4edata =
00525     G4ErrorPropagatorData::GetErrorPropagatorData();
00526 
00527 #ifdef G4EVERBOSE
00528   if( verbose >= 4 )
00529   {
00530     G4cout << " G4ErrorPropagator::CheckIfLastStep G4ErrorState= "
00531            << G4int(g4edata->GetState()) << G4endl;
00532   }
00533 #endif
00534   
00535   //----- Check if this is the last step: track has reached the target
00536   //      or the end of world
00537   //
00538   if(g4edata->GetState() == G4ErrorState(G4ErrorState_StoppedAtTarget) )
00539   {
00540     lastG4eStep = true;    
00541 #ifdef G4EVERBOSE
00542     if(verbose >= 5 )
00543     {
00544       G4cout << " G4ErrorPropagator::CheckIfLastStep " << lastG4eStep
00545              << " " << G4int(g4edata->GetState()) << G4endl;
00546     }
00547 #endif
00548   }
00549   else if( aTrack->GetNextVolume() == 0 )
00550   {
00551     //----- If particle is out of world, without finding the G4ErrorTarget,
00552     //      give a n error/warning
00553     //
00554     lastG4eStep = true;
00555     if( exception )
00556     {
00557       std::ostringstream message;
00558       message << "Track extrapolated until end of World" << G4endl
00559               << "without finding the defined target!";
00560       G4Exception("G4ErrorPropagator::CheckIfLastStep()",
00561                   "InvalidSetup", FatalException, message);
00562     }
00563     else
00564     {
00565       if( verbose >= 1 )
00566       {
00567         G4cerr << "ERROR - G4ErrorPropagator::CheckIfLastStep()" << G4endl
00568                << "        Track extrapolated until end of World" << G4endl
00569                << "        without finding the defined target " << G4endl;
00570       }
00571     }
00572   }  //----- not last step from G4e, but track is stopped (energy exhausted)
00573   else if( aTrack->GetTrackStatus() == fStopAndKill )
00574   { 
00575     if( exception )
00576     {
00577       std::ostringstream message;
00578       message << "Track extrapolated until energy is exhausted" << G4endl
00579               << "without finding the defined target!";
00580       G4Exception("G4ErrorPropagator::CheckIfLastStep()",
00581                   "InvalidSetup", FatalException, message);
00582     }
00583     else
00584     {
00585       if( verbose >= 1 )
00586       {
00587         G4cerr << "ERROR - G4ErrorPropagator::CheckIfLastStep()" << G4endl
00588              << "        Track extrapolated until energy is exhausted" << G4endl
00589              << "        without finding the defined target !" << G4endl;
00590       }
00591       lastG4eStep = 1;
00592     }
00593   }
00594 
00595 #ifdef G4EVERBOSE
00596   if( verbose >= 5 )
00597   {
00598     G4cout << " return CheckIfLastStep " << lastG4eStep << G4endl;
00599   }
00600 #endif
00601 
00602   return  lastG4eStep;
00603 }
00604 
00605 
00606 //---------------------------------------------------------------------------
00607 void G4ErrorPropagator::InvokePreUserTrackingAction( G4Track* fpTrack )
00608 {
00609   const G4UserTrackingAction* fpUserTrackingAction =
00610     G4EventManager::GetEventManager()->GetUserTrackingAction();
00611   if( fpUserTrackingAction != 0 )
00612   {
00613     const_cast<G4UserTrackingAction*>(fpUserTrackingAction)
00614       ->PreUserTrackingAction((fpTrack) );
00615   }
00616 }
00617 
00618 
00619 //---------------------------------------------------------------------------
00620 void G4ErrorPropagator::InvokePostUserTrackingAction( G4Track* fpTrack )
00621 {
00622   const G4UserTrackingAction* fpUserTrackingAction =
00623     G4EventManager::GetEventManager()->GetUserTrackingAction();
00624   if( fpUserTrackingAction != 0 )
00625   {
00626     const_cast<G4UserTrackingAction*>(fpUserTrackingAction)
00627       ->PostUserTrackingAction((fpTrack) );
00628   }
00629 }

Generated on Mon May 27 17:48:11 2013 for Geant4 by  doxygen 1.4.7