G4ITTransportation.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: G4ITTransportation.cc 64374 2012-10-31 16:37:23Z gcosmo $
00027 //
00031 //
00032 // Contact : Mathieu Karamitros (kara (AT) cenbg . in2p3 . fr)
00033 //
00034 // History :
00035 // -----------
00036 // =======================================================================
00037 // Modified:
00038 //   28 Oct  2011, P.Gumpl./J.Ap: Detect gravity field, use magnetic moment
00039 //   20 Nov  2008, J.Apostolakis: Push safety to helper - after ComputeSafety
00040 //    9 Nov  2007, J.Apostolakis: Flag for short steps, push safety to helper
00041 //   19 Jan  2006, P.MoraDeFreitas: Fix for suspended tracks (StartTracking)
00042 //   11 Aug  2004, M.Asai: Add G4VSensitiveDetector* for updating stepPoint.
00043 //   21 June 2003, J.Apostolakis: Calling field manager with
00044 //                     track, to enable it to configure its accuracy
00045 //   13 May  2003, J.Apostolakis: Zero field areas now taken into
00046 //                     account correclty in all cases (thanks to W Pokorski).
00047 //   29 June 2001, J.Apostolakis, D.Cote-Ahern, P.Gumplinger:
00048 //                     correction for spin tracking
00049 //   20 Febr 2001, J.Apostolakis:  update for new FieldTrack
00050 //   22 Sept 2000, V.Grichine:     update of Kinetic Energy
00051 //  ---------------------------------------------------
00052 //   10 Oct  2011, M.Karamitros:   G4ITTransportation created
00053 // Created:  19 March 1997, J. Apostolakis
00054 // =======================================================================
00055 //
00056 // -------------------------------------------------------------------
00057 
00058 #include "G4ITTransportation.hh"
00059 #include "G4SystemOfUnits.hh"
00060 #include "G4TransportationManager.hh"
00061 #include "G4ITTransportationManager.hh"
00062 #include "G4ProductionCutsTable.hh"
00063 #include "G4ParticleTable.hh"
00064 #include "G4ITNavigator.hh"
00065 #include "G4PropagatorInField.hh"
00066 #include "G4FieldManager.hh"
00067 #include "G4ChordFinder.hh"
00068 #include "G4SafetyHelper.hh"
00069 #include "G4FieldManagerStore.hh"
00070 
00071 #include "G4UnitsTable.hh"
00072 
00073 class G4VSensitiveDetector;
00074 
00075 #ifndef State
00076 #define State(theXInfo) (fTransportationState->theXInfo)
00077 #endif
00078 
00079 //#define G4DEBUG_TRANSPORT 1
00080 
00081 G4ITTransportation::G4ITTransportation(const G4String& aName, int verbose) :
00082     G4VITProcess(aName, fTransportation),
00083     InitProcessState(fTransportationState, fpState),
00084     fThreshold_Warning_Energy( 100 * MeV ),
00085     fThreshold_Important_Energy( 250 * MeV ),
00086     fThresholdTrials( 10 ),
00087     fUnimportant_Energy( 1 * MeV ),  // Not used
00088     fSumEnergyKilled( 0.0 ), fMaxEnergyKilled( 0.0 ),
00089     fShortStepOptimisation(false),    // Old default: true (=fast short steps)
00090     fVerboseLevel( verbose )
00091 {
00092     pParticleChange = &fParticleChange;
00093     G4TransportationManager* transportMgr ;
00094     transportMgr = G4TransportationManager::GetTransportationManager() ;
00095     G4ITTransportationManager* ITtransportMgr ;
00096     ITtransportMgr = G4ITTransportationManager::GetTransportationManager() ;
00097     fLinearNavigator = ITtransportMgr->GetNavigatorForTracking() ;
00098     fFieldPropagator = transportMgr->GetPropagatorInField() ;
00099     fpSafetyHelper =   0; // transportMgr->GetSafetyHelper();  // New
00100 
00101     // Cannot determine whether a field exists here, as it would
00102     //  depend on the relative order of creating the detector's
00103     //  field and this process. That order is not guaranted.
00104     // Instead later the method DoesGlobalFieldExist() is called
00105 
00106     enableAtRestDoIt    = false;
00107     enableAlongStepDoIt = true;
00108     enablePostStepDoIt  = true;
00109     SetProcessSubType(60);
00110     SetInstantiateProcessState(true);
00111     G4VITProcess::SetInstantiateProcessState(false);
00112     fInstantiateProcessState = true;
00113 }
00114 
00115 
00116 G4ITTransportation::G4ITTransportation(const G4ITTransportation& right) :
00117     G4VITProcess(right),
00118     InitProcessState(fTransportationState, fpState)
00119 {
00120     // Copy attributes
00121     fVerboseLevel               = right.fVerboseLevel ;
00122     fThreshold_Warning_Energy   = right.fThreshold_Warning_Energy;
00123     fThreshold_Important_Energy = right.fThreshold_Important_Energy;
00124     fThresholdTrials            = right.fThresholdTrials;
00125     fUnimportant_Energy         = right.fUnimportant_Energy;
00126     fSumEnergyKilled            = right.fSumEnergyKilled;
00127     fMaxEnergyKilled            = right.fMaxEnergyKilled;
00128     fShortStepOptimisation      = right.fShortStepOptimisation;
00129 
00130     // Setup Navigators
00131     G4TransportationManager* transportMgr ;
00132     transportMgr = G4TransportationManager::GetTransportationManager() ;
00133     G4ITTransportationManager* ITtransportMgr ;
00134     ITtransportMgr = G4ITTransportationManager::GetTransportationManager() ;
00135     fLinearNavigator = ITtransportMgr->GetNavigatorForTracking() ;
00136     fFieldPropagator = transportMgr->GetPropagatorInField() ;
00137     fpSafetyHelper =   0; //transportMgr->GetSafetyHelper();  // New
00138 
00139     // Cannot determine whether a field exists here, as it would
00140     //  depend on the relative order of creating the detector's
00141     //  field and this process. That order is not guaranted.
00142     // Instead later the method DoesGlobalFieldExist() is called
00143 
00144     enableAtRestDoIt    = false;
00145     enableAlongStepDoIt = true;
00146     enablePostStepDoIt  = true;
00147 
00148     pParticleChange = &fParticleChange;
00149     SetInstantiateProcessState(true);
00150     G4VITProcess::SetInstantiateProcessState(false);
00151     fInstantiateProcessState = right.fInstantiateProcessState;
00152 }
00153 
00154 G4ITTransportation& G4ITTransportation::operator=(const G4ITTransportation& right)
00155 {
00156     if(this == &right) return *this;
00157     return *this;
00158 }
00159 
00163 G4ITTransportation::G4ITTransportationState::G4ITTransportationState() : G4ProcessState(),
00164     fCurrentTouchableHandle(0)
00165 {
00166     fTransportEndPosition = G4ThreeVector(0,0,0);
00167     fTransportEndMomentumDir = G4ThreeVector(0,0,0);
00168     fTransportEndKineticEnergy = -1;
00169     fTransportEndSpin = G4ThreeVector(0,0,0);
00170     fMomentumChanged = false;
00171     fEnergyChanged = false;
00172     fEndGlobalTimeComputed = false;
00173     fCandidateEndGlobalTime = -1;
00174     fParticleIsLooping = false;
00175     static G4TouchableHandle nullTouchableHandle; // Points to (G4VTouchable*) 0
00176     fCurrentTouchableHandle = nullTouchableHandle;
00177     fGeometryLimitedStep = false;
00178     fPreviousSftOrigin = G4ThreeVector(0,0,0);
00179     fPreviousSafety = 0.0;
00180     fNoLooperTrials = false;
00181     endpointDistance= -1;
00182 }
00183 
00184 G4ITTransportation::G4ITTransportationState::~G4ITTransportationState()
00185 {
00186     ;
00187 }
00188 
00189 G4ITTransportation::~G4ITTransportation()
00190 {
00191 #ifdef G4VERBOSE
00192     if( (fVerboseLevel > 0) && (fSumEnergyKilled > 0.0 ) )
00193     {
00194         G4cout << " G4ITTransportation: Statistics for looping particles " << G4endl;
00195         G4cout << "   Sum of energy of loopers killed: " <<  fSumEnergyKilled << G4endl;
00196         G4cout << "   Max energy of loopers killed: " <<  fMaxEnergyKilled << G4endl;
00197     }
00198 #endif
00199 }
00200 
00201 G4bool G4ITTransportation::DoesGlobalFieldExist()
00202 {
00203     G4TransportationManager* transportMgr;
00204     transportMgr= G4TransportationManager::GetTransportationManager();
00205 
00206     // fFieldExists= transportMgr->GetFieldManager()->DoesFieldExist();
00207     // return fFieldExists;
00208     return transportMgr->GetFieldManager()->DoesFieldExist();
00209 }
00210 
00211 
00213 //
00214 // Responsibilities:
00215 //    Find whether the geometry limits the Step, and to what length
00216 //    Calculate the new value of the safety and return it.
00217 //    Store the final time, position and momentum.
00218 G4double G4ITTransportation::AlongStepGetPhysicalInteractionLength(
00219         const G4Track& track,
00220         G4double ,
00221         G4double currentMinimumStep,
00222         G4double& currentSafety,
00223         G4GPILSelection* selection)
00224 {
00225     G4double geometryStepLength(-1.0), newSafety(-1.0) ;
00226 
00227     State(fParticleIsLooping) = false ;
00228     State(fEndGlobalTimeComputed) = false ;
00229     State(fGeometryLimitedStep) = false ;
00230 
00231     // Initial actions moved to  StartTrack()
00232     // --------------------------------------
00233     // Note: in case another process changes touchable handle
00234     //    it will be necessary to add here (for all steps)
00235     State(fCurrentTouchableHandle) = track.GetTouchableHandle();
00236 
00237     // GPILSelection is set to defaule value of CandidateForSelection
00238     // It is a return value
00239     //
00240     *selection = CandidateForSelection ;
00241 
00242     // Get initial Energy/Momentum of the track
00243     //
00244     const G4DynamicParticle*    pParticle  = track.GetDynamicParticle() ;
00245 //    const G4ParticleDefinition* pParticleDef   = pParticle->GetDefinition() ;
00246     G4ThreeVector startMomentumDir       = pParticle->GetMomentumDirection() ;
00247     G4ThreeVector startPosition          = track.GetPosition() ;
00248 
00249     // G4double   theTime        = track.GetGlobalTime() ;
00250 
00251     // The Step Point safety can be limited by other geometries and/or the
00252     // assumptions of any process - it's not always the geometrical safety.
00253     // We calculate the starting point's isotropic safety here.
00254     //
00255     G4ThreeVector OriginShift = startPosition - State(fPreviousSftOrigin) ;
00256     G4double      MagSqShift  = OriginShift.mag2() ;
00257     if( MagSqShift >= sqr(State(fPreviousSafety)) )
00258     {
00259         currentSafety = 0.0 ;
00260     }
00261     else
00262     {
00263         currentSafety = State(fPreviousSafety) - std::sqrt(MagSqShift) ;
00264     }
00265 
00266     // Is the particle charged ?
00267     //
00268     G4double              particleCharge = pParticle->GetCharge() ;
00269 
00270 
00271     // There is no need to locate the current volume. It is Done elsewhere:
00272     //   On track construction
00273     //   By the tracking, after all AlongStepDoIts, in "Relocation"
00274 
00275     // Check whether the particle have an (EM) field force exerting upon it
00276     //
00277     G4FieldManager* fieldMgr=0;
00278     G4bool          fieldExertsForce = false ;
00279     if( (particleCharge != 0.0) )
00280     {
00281         fieldMgr= fFieldPropagator->FindAndSetFieldManager( track.GetVolume() );
00282         if (fieldMgr != 0)
00283         {
00284             // Message the field Manager, to configure it for this track
00285             fieldMgr->ConfigureForTrack( &track );
00286             // Moved here, in order to allow a transition
00287             //   from a zero-field  status (with fieldMgr->(field)0
00288             //   to a finite field  status
00289 
00290             // If the field manager has no field, there is no field !
00291             fieldExertsForce = (fieldMgr->GetDetectorField() != 0);
00292         }
00293     }
00294 
00295     // G4cout << " G4Transport:  field exerts force= " << fieldExertsForce
00296     //   << "  fieldMgr= " << fieldMgr << G4endl;
00297 
00298     // Choose the calculation of the transportation: Field or not
00299     //
00300     if( !fieldExertsForce )
00301     {
00302         G4double linearStepLength ;
00303         if( fShortStepOptimisation && (currentMinimumStep <= currentSafety) )
00304         {
00305             // The Step is guaranteed to be taken
00306             //
00307             geometryStepLength   = currentMinimumStep ;
00308             State(fGeometryLimitedStep) = false ;
00309         }
00310         else
00311         {
00312             //  Find whether the straight path intersects a volume
00313             //
00314             linearStepLength = fLinearNavigator->ComputeStep( startPosition,
00315                                                               startMomentumDir,
00316                                                               currentMinimumStep,
00317                                                               newSafety) ;
00318             // Remember last safety origin & value.
00319             //
00320             State(fPreviousSftOrigin) = startPosition ;
00321             State(fPreviousSafety)    = newSafety ;
00322             // fpSafetyHelper->SetCurrentSafety( newSafety, startPosition);
00323 
00324             // The safety at the initial point has been re-calculated:
00325             //
00326             currentSafety = newSafety ;
00327 
00328             State(fGeometryLimitedStep)= (linearStepLength <= currentMinimumStep);
00329             if( State(fGeometryLimitedStep) )
00330             {
00331                 // The geometry limits the Step size (an intersection was found.)
00332                 geometryStepLength   = linearStepLength ;
00333             }
00334             else
00335             {
00336                 // The full Step is taken.
00337                 geometryStepLength   = currentMinimumStep ;
00338             }
00339         }
00340         State(endpointDistance) = geometryStepLength ;
00341 
00342         // Calculate final position
00343         //
00344         State(fTransportEndPosition) = startPosition+geometryStepLength*startMomentumDir ;
00345 
00346         // Momentum direction, energy and polarisation are unchanged by transport
00347         //
00348         State(fTransportEndMomentumDir)   = startMomentumDir ;
00349         State(fTransportEndKineticEnergy) = track.GetKineticEnergy() ;
00350         State(fTransportEndSpin)          = track.GetPolarization();
00351         State(fParticleIsLooping)         = false ;
00352         State(fMomentumChanged)           = false ;
00353         State(fEndGlobalTimeComputed)     = true ;
00354         State(theInteractionTimeLeft)     = State(endpointDistance)/track.GetVelocity();
00355         State(fCandidateEndGlobalTime)    = State(theInteractionTimeLeft)+track.GetGlobalTime();
00356 
00357         //        G4cout << "track.GetVelocity() : " << track.GetVelocity() << G4endl;
00358         //        G4cout << "State(endpointDistance) : " << G4BestUnit(State(endpointDistance),"Length") << G4endl;
00359         //        G4cout << "State(theInteractionTimeLeft) : " << G4BestUnit(State(theInteractionTimeLeft),"Time") << G4endl;
00360         //        G4cout << "track.GetGlobalTime() : " << G4BestUnit(track.GetGlobalTime(),"Time") << G4endl;
00361 
00362     }
00363     else   //  A field exerts force
00364     {
00365 
00366         G4ExceptionDescription exceptionDescription;
00367         exceptionDescription << "ITTransportation does not support external fields.";
00368         exceptionDescription << " If you are dealing with a tradiational MC simulation, ";
00369         exceptionDescription << "please use G4Transportation.";
00370 
00371         G4Exception("G4ITTransportation::AlongStepGetPhysicalInteractionLength","NoExternalFieldSupport",
00372                     FatalException,exceptionDescription);
00373  /*
00374         G4double       momentumMagnitude = pParticle->GetTotalMomentum() ;
00375         //        G4ThreeVector  EndUnitMomentum ;
00376         G4double       lengthAlongCurve ;
00377         G4double       restMass = pParticleDef->GetPDGMass() ;
00378 
00379         fFieldPropagator->SetChargeMomentumMass( particleCharge,    // in e+ units
00380                                                  momentumMagnitude, // in Mev/c
00381                                                  restMass           ) ;
00382 
00383         G4ThreeVector spin        = track.GetPolarization() ;
00384         G4FieldTrack  aFieldTrack = G4FieldTrack( startPosition,
00385                                                   track.GetMomentumDirection(),
00386                                                   0.0,
00387                                                   track.GetKineticEnergy(),
00388                                                   restMass,
00389                                                   track.GetVelocity(),
00390                                                   track.GetGlobalTime(), // Lab.
00391                                                   track.GetProperTime(), // Part.
00392                                                   &spin                  ) ;
00393         if( currentMinimumStep > 0 )
00394         {
00395             // Do the Transport in the field (non recti-linear)
00396             //
00397             lengthAlongCurve = fFieldPropagator->ComputeStep( aFieldTrack,
00398                                                               currentMinimumStep,
00399                                                               currentSafety,
00400                                                               track.GetVolume() ) ;
00401             State(fGeometryLimitedStep)= lengthAlongCurve < currentMinimumStep;
00402             if( State(fGeometryLimitedStep) )
00403             {
00404                 geometryStepLength   = lengthAlongCurve ;
00405             }
00406             else
00407             {
00408                 geometryStepLength   = currentMinimumStep ;
00409             }
00410 
00411             // Remember last safety origin & value.
00412             //
00413             State(fPreviousSftOrigin) = startPosition ;
00414             State(fPreviousSafety)    = currentSafety ;
00415             fpSafetyHelper->SetCurrentSafety( newSafety, startPosition);
00416         }
00417         else
00418         {
00419             geometryStepLength   = lengthAlongCurve= 0.0 ;
00420             State(fGeometryLimitedStep) = false ;
00421         }
00422 
00423         // Get the End-Position and End-Momentum (Dir-ection)
00424         //
00425         State(fTransportEndPosition) = aFieldTrack.GetPosition() ;
00426 
00427         // Momentum:  Magnitude and direction can be changed too now ...
00428         //
00429         State(fMomentumChanged)         = true ;
00430         State(fTransportEndMomentumDir) = aFieldTrack.GetMomentumDir() ;
00431 
00432         State(fTransportEndKineticEnergy)  = aFieldTrack.GetKineticEnergy() ;
00433 
00434         if( fFieldPropagator->GetCurrentFieldManager()->DoesFieldChangeEnergy() )
00435         {
00436             // If the field can change energy, then the time must be integrated
00437             //    - so this should have been updated
00438             //
00439             State(fCandidateEndGlobalTime)   = aFieldTrack.GetLabTimeOfFlight();
00440             State(fEndGlobalTimeComputed)    = true;
00441 
00442             State(theInteractionTimeLeft) = State(fCandidateEndGlobalTime) - track.GetGlobalTime() ;
00443 
00444             // was ( State(fCandidateEndGlobalTime) != track.GetGlobalTime() );
00445             // a cleaner way is to have FieldTrack knowing whether time is updated.
00446         }
00447         else
00448         {
00449             // The energy should be unchanged by field transport,
00450             //    - so the time changed will be calculated elsewhere
00451             //
00452             State(fEndGlobalTimeComputed) = false;
00453 
00454             // Check that the integration preserved the energy
00455             //     -  and if not correct this!
00456             G4double  startEnergy= track.GetKineticEnergy();
00457             G4double  endEnergy= State(fTransportEndKineticEnergy);
00458 
00459             static G4int no_inexact_steps=0, no_large_ediff;
00460             G4double absEdiff = std::fabs(startEnergy- endEnergy);
00461             if( absEdiff > perMillion * endEnergy )
00462             {
00463                 no_inexact_steps++;
00464                 // Possible statistics keeping here ...
00465             }
00466 #ifdef G4VERBOSE
00467             if( fVerboseLevel > 1 )
00468             {
00469                 if( std::fabs(startEnergy- endEnergy) > perThousand * endEnergy )
00470                 {
00471                     static G4int no_warnings= 0, warnModulo=1,  moduloFactor= 10;
00472                     no_large_ediff ++;
00473                     if( (no_large_ediff% warnModulo) == 0 )
00474                     {
00475                         no_warnings++;
00476                         G4cout << "WARNING - G4Transportation::AlongStepGetPIL() "
00477                                << "   Energy change in Step is above 1^-3 relative value. " << G4endl
00478                                << "   Relative change in 'tracking' step = "
00479                                << std::setw(15) << (endEnergy-startEnergy)/startEnergy << G4endl
00480                                << "     Starting E= " << std::setw(12) << startEnergy / MeV << " MeV " << G4endl
00481                                << "     Ending   E= " << std::setw(12) << endEnergy   / MeV << " MeV " << G4endl;
00482                         G4cout << " Energy has been corrected -- however, review"
00483                                << " field propagation parameters for accuracy."  << G4endl;
00484                         if( (fVerboseLevel > 2 ) || (no_warnings<4) || (no_large_ediff == warnModulo * moduloFactor) )
00485                         {
00486                             G4cout << " These include EpsilonStepMax(/Min) in G4FieldManager "
00487                                    << " which determine fractional error per step for integrated quantities. " << G4endl
00488                                    << " Note also the influence of the permitted number of integration steps."
00489                                    << G4endl;
00490                         }
00491                         G4cerr << "ERROR - G4Transportation::AlongStepGetPIL()" << G4endl
00492                                << "        Bad 'endpoint'. Energy change detected"
00493                                << " and corrected. "
00494                                << " Has occurred already "
00495                                << no_large_ediff << " times." << G4endl;
00496                         if( no_large_ediff == warnModulo * moduloFactor )
00497                         {
00498                             warnModulo *= moduloFactor;
00499                         }
00500                     }
00501                 }
00502             }  // end of if (fVerboseLevel)
00503 #endif
00504             // Correct the energy for fields that conserve it
00505             //  This - hides the integration error
00506             //       - but gives a better physical answer
00507             State(fTransportEndKineticEnergy)= track.GetKineticEnergy();
00508         }
00509 
00510         State(fTransportEndSpin) = aFieldTrack.GetSpin();
00511         State(fParticleIsLooping) = fFieldPropagator->IsParticleLooping() ;
00512         State(endpointDistance)   = (State(fTransportEndPosition) - startPosition).mag() ;
00513         //        State(theInteractionTimeLeft) = track.GetVelocity()/State(endpointDistance) ;
00514 */
00515     }
00516 
00517     // If we are asked to go a step length of 0, and we are on a boundary
00518     // then a boundary will also limit the step -> we must flag this.
00519     //
00520     if( currentMinimumStep == 0.0 )
00521     {
00522         if( currentSafety == 0.0 )
00523         {
00524             State(fGeometryLimitedStep) = true ;
00525 //            G4cout << "!!!! Safety is NULL, on the Boundary !!!!!" << G4endl;
00526 //            G4cout << " Track position : " << track.GetPosition() /nanometer << G4endl;
00527         }
00528     }
00529 
00530     // Update the safety starting from the end-point,
00531     // if it will become negative at the end-point.
00532     //
00533     if( currentSafety < State(endpointDistance) )
00534     {
00535         // if( particleCharge == 0.0 )
00536         //    G4cout  << "  Avoiding call to ComputeSafety : charge = 0.0 " << G4endl;
00537 
00538         if( particleCharge != 0.0 )
00539         {
00540 
00541             G4double endSafety =
00542                     fLinearNavigator->ComputeSafety( State(fTransportEndPosition)) ;
00543             currentSafety      = endSafety ;
00544             State(fPreviousSftOrigin) = State(fTransportEndPosition) ;
00545             State(fPreviousSafety)    = currentSafety ;
00546             fpSafetyHelper->SetCurrentSafety( currentSafety, State(fTransportEndPosition));
00547 
00548             // Because the Stepping Manager assumes it is from the start point,
00549             //  add the StepLength
00550             //
00551             currentSafety     += State(endpointDistance) ;
00552 
00553 #ifdef G4DEBUG_TRANSPORT
00554             G4cout.precision(12) ;
00555             G4cout << "***G4Transportation::AlongStepGPIL ** " << G4endl  ;
00556             G4cout << "  Called Navigator->ComputeSafety at " << State(fTransportEndPosition)
00557                    << "    and it returned safety= " << endSafety << G4endl ;
00558             G4cout << "  Adding endpoint distance " << State(endpointDistance)
00559                    << "   to obtain pseudo-safety= " << currentSafety << G4endl ;
00560 #endif
00561         }
00562     }
00563 
00564     //    fParticleChange.ProposeTrueStepLength(geometryStepLength) ;
00565 
00566     return geometryStepLength ;
00567 }
00568 
00569 
00570 void G4ITTransportation::ComputeStep(const G4Track& track,
00571                                      const G4Step& /*step*/,
00572                                      const double timeStep,
00573                                      double& oPhysicalStep)
00574 {
00575     const G4DynamicParticle*    pParticle  = track.GetDynamicParticle() ;
00576     G4ThreeVector startMomentumDir       = pParticle->GetMomentumDirection() ;
00577     G4ThreeVector startPosition          = track.GetPosition() ;
00578 
00579     track.CalculateVelocity();
00580     G4double initialVelocity   = track.GetVelocity() ;
00581 
00582     State(fGeometryLimitedStep) = false;
00583 
00585     // !!! CASE NO FIELD !!!
00587     State(fCandidateEndGlobalTime) = timeStep + track.GetGlobalTime();
00588     State(fEndGlobalTimeComputed) = true ;
00589 
00590     // Choose the calculation of the transportation: Field or not
00591     //
00592     if( !State(fMomentumChanged) )
00593     {
00594         //        G4cout << "Momentum has not changed" << G4endl;
00595         fParticleChange.ProposeVelocity(initialVelocity);
00596         oPhysicalStep =  initialVelocity*timeStep ;
00597 
00598         // Calculate final position
00599         //
00600         State(fTransportEndPosition) = startPosition + oPhysicalStep*startMomentumDir ;
00601     }
00602 }
00603 
00604 
00606 //
00607 //   Initialize ParticleChange  (by setting all its members equal
00608 //                               to corresponding members in G4Track)
00609 #include "G4ParticleTable.hh"
00610 G4VParticleChange* G4ITTransportation::AlongStepDoIt( const G4Track& track,
00611                                                       const G4Step&  stepData )
00612 {
00613 
00614     //    G4cout << "G4ITTransportation::AlongStepDoIt" << G4endl;
00615     // set  pdefOpticalPhoton
00616     static  G4ParticleDefinition* pdefOpticalPhoton =
00617             G4ParticleTable::GetParticleTable()->FindParticle("opticalphoton");
00618 
00619     static G4int noCalls=0;
00620     noCalls++;
00621 
00622     fParticleChange.Initialize(track) ;
00623 
00624     //  Code for specific process
00625     //
00626     fParticleChange.ProposePosition(State(fTransportEndPosition)) ;
00627     fParticleChange.ProposeMomentumDirection(State(fTransportEndMomentumDir)) ;
00628     fParticleChange.ProposeEnergy(State(fTransportEndKineticEnergy)) ;
00629     fParticleChange.SetMomentumChanged(State(fMomentumChanged)) ;
00630 
00631     fParticleChange.ProposePolarization(State(fTransportEndSpin));
00632 
00633     G4double deltaTime = 0.0 ;
00634 
00635     // Calculate  Lab Time of Flight (ONLY if field Equations used it!)
00636     // G4double endTime   = State(fCandidateEndGlobalTime);
00637     // G4double delta_time = endTime - startTime;
00638 
00639     G4double startTime = track.GetGlobalTime() ;
00643     if (State(fEndGlobalTimeComputed) == false)
00644     {
00645         // The time was not integrated .. make the best estimate possible
00646         //
00647         G4double initialVelocity = stepData.GetPreStepPoint()->GetVelocity() ;
00648         G4double stepLength      = track.GetStepLength() ;
00649 
00650         deltaTime= 0.0;  // in case initialVelocity = 0
00651         if (track.GetParticleDefinition() == pdefOpticalPhoton)
00652         {
00653             // For only Optical Photon, final velocity is used
00654             double finalVelocity = track.CalculateVelocityForOpticalPhoton();
00655             fParticleChange.ProposeVelocity(finalVelocity);
00656             deltaTime = stepLength/finalVelocity ;
00657         }
00658         else if( initialVelocity > 0.0 )
00659         {
00660             deltaTime = stepLength/initialVelocity ;
00661         }
00662 
00663         State(fCandidateEndGlobalTime)   = startTime + deltaTime ;
00664     }
00665     else
00666     {
00667         deltaTime = State(fCandidateEndGlobalTime) - startTime ;
00668     }
00669 
00670     fParticleChange.ProposeGlobalTime( State(fCandidateEndGlobalTime) ) ;
00671     fParticleChange.ProposeLocalTime(  track.GetLocalTime() + deltaTime) ;
00672     /*
00673     // Now Correct by Lorentz factor to get delta "proper" Time
00674 
00675     G4double  restMass       = track.GetDynamicParticle()->GetMass() ;
00676     G4double deltaProperTime = deltaTime*( restMass/track.GetTotalEnergy() ) ;
00677 
00678     fParticleChange.ProposeProperTime(track.GetProperTime() + deltaProperTime) ;
00679 */
00680 
00681     fParticleChange. ProposeTrueStepLength( track.GetStepLength() ) ;
00682 
00685 
00686     // If the particle is caught looping or is stuck (in very difficult
00687     // boundaries) in a magnetic field (doing many steps)
00688     //   THEN this kills it ...
00689     //
00690     if ( State(fParticleIsLooping) )
00691     {
00692         G4double endEnergy= State(fTransportEndKineticEnergy);
00693 
00694         if( (endEnergy < fThreshold_Important_Energy)
00695                 || (State(fNoLooperTrials) >= fThresholdTrials ) )
00696         {
00697             // Kill the looping particle
00698             //
00699             //            G4cout << "G4ITTransportation will killed the molecule"<< G4endl;
00700             fParticleChange.ProposeTrackStatus( fStopAndKill )  ;
00701 
00702             // 'Bare' statistics
00703             fSumEnergyKilled += endEnergy;
00704             if( endEnergy > fMaxEnergyKilled)
00705             {
00706                 fMaxEnergyKilled= endEnergy;
00707             }
00708 
00709 #ifdef G4VERBOSE
00710             if( (fVerboseLevel > 1) ||
00711                     ( endEnergy > fThreshold_Warning_Energy )  )
00712             {
00713                 G4cout << " G4ITTransportation is killing track that is looping or stuck "
00714                        << G4endl
00715                        << "   This track has " << track.GetKineticEnergy() / MeV
00716                        << " MeV energy." << G4endl;
00717                 G4cout << "   Number of trials = " << State(fNoLooperTrials)
00718                        << "   No of calls to AlongStepDoIt = " << noCalls
00719                        << G4endl;
00720             }
00721 #endif
00722             State(fNoLooperTrials)=0;
00723         }
00724         else
00725         {
00726             State(fNoLooperTrials) ++;
00727 #ifdef G4VERBOSE
00728             if( (fVerboseLevel > 2) )
00729             {
00730                 G4cout << "   G4ITTransportation::AlongStepDoIt(): Particle looping -  "
00731                        << "   Number of trials = " << State(fNoLooperTrials)
00732                        << "   No of calls to  = " << noCalls
00733                        << G4endl;
00734             }
00735 #endif
00736         }
00737     }
00738     else
00739     {
00740         State(fNoLooperTrials)=0;
00741     }
00742 
00743     // Another (sometimes better way) is to use a user-limit maximum Step size
00744     // to alleviate this problem ..
00745 
00746     // Introduce smooth curved trajectories to particle-change
00747     //
00748     fParticleChange.SetPointerToVectorOfAuxiliaryPoints
00749             (fFieldPropagator->GimmeTrajectoryVectorAndForgetIt() );
00750 
00751     return &fParticleChange ;
00752 }
00753 
00755 //
00756 //  This ensures that the PostStep action is always called,
00757 //  so that it can do the relocation if it is needed.
00758 //
00759 
00760 G4double G4ITTransportation::
00761 PostStepGetPhysicalInteractionLength(
00762         const G4Track&  , // track
00763         G4double , // previousStepSize
00764         G4ForceCondition* pForceCond
00765         )
00766 {
00767     *pForceCond = Forced ;
00768     return DBL_MAX ;  // was kInfinity ; but convention now is DBL_MAX
00769 }
00770 
00772 //
00773 
00774 G4VParticleChange* G4ITTransportation::PostStepDoIt( const G4Track& track,
00775                                                      const G4Step& )
00776 {
00777     //    G4cout << "G4ITTransportation::PostStepDoIt" << G4endl;
00778     G4TouchableHandle retCurrentTouchable ;   // The one to return
00779     G4bool isLastStep= false;
00780 
00781     // Initialize ParticleChange  (by setting all its members equal
00782     //                             to corresponding members in G4Track)
00783     fParticleChange.Initialize(track) ;  // To initialise TouchableChange
00784 
00785     fParticleChange.ProposeTrackStatus(track.GetTrackStatus()) ;
00786 
00787     // If the Step was determined by the volume boundary,
00788     // logically relocate the particle
00789 
00790     if(State(fGeometryLimitedStep))
00791     {
00792 
00793         //        G4cout << "Step is limited by geometry " <<  "track ID : " << track.GetTrackID() << G4endl;
00794 
00795         // fCurrentTouchable will now become the previous touchable,
00796         // and what was the previous will be freed.
00797         // (Needed because the preStepPoint can point to the previous touchable)
00798 
00799         if( State(fCurrentTouchableHandle)->GetVolume() == 0 )
00800         {
00801             G4ExceptionDescription exceptionDescription ;
00802             exceptionDescription << "No current touchable found " ;
00803             G4Exception(" G4ITTransportation::PostStepDoIt","G4ITTransportation001",
00804                         FatalErrorInArgument,exceptionDescription);
00805         }
00806 
00807         fLinearNavigator->SetGeometricallyLimitedStep() ;
00808         fLinearNavigator->
00809                 LocateGlobalPointAndUpdateTouchableHandle( track.GetPosition(),
00810                                                            track.GetMomentumDirection(),
00811                                                            State(fCurrentTouchableHandle),
00812                                                            true                      ) ;
00813         // Check whether the particle is out of the world volume
00814         // If so it has exited and must be killed.
00815         //
00816         if( State(fCurrentTouchableHandle)->GetVolume() == 0 )
00817         {
00818 #ifdef G4VERBOSE
00819             if(fVerboseLevel > 0)
00820             {
00821                 G4cout << "Track position : " << track.GetPosition() / nanometer << " [nm]"
00822                        << " Track ID : " << track.GetTrackID()<< G4endl;
00823                 G4cout << "G4ITTransportation will killed the track because State(fCurrentTouchableHandle)->GetVolume() == 0"<< G4endl;
00824             }
00825 #endif
00826             fParticleChange.ProposeTrackStatus( fStopAndKill ) ;
00827         }
00828 
00829         retCurrentTouchable = State(fCurrentTouchableHandle) ;
00830 
00831 //        G4cout << "Current volume : " << track.GetVolume()->GetName()
00832 //               << " Next volume : "
00833 //               << (State(fCurrentTouchableHandle)->GetVolume() ? State(fCurrentTouchableHandle)->GetVolume()->GetName():"OutWorld")
00834 //               << " Position : " << track.GetPosition() / nanometer
00835 //               << " track ID : " << track.GetTrackID()
00836 //               << G4endl;
00837 
00838         fParticleChange.SetTouchableHandle( State(fCurrentTouchableHandle) ) ;
00839 
00840         // Update the Step flag which identifies the Last Step in a volume
00841         isLastStep =  fLinearNavigator->ExitedMotherVolume()
00842                 | fLinearNavigator->EnteredDaughterVolume() ;
00843 
00844 #ifdef G4DEBUG_TRANSPORT
00845         //  Checking first implementation of flagging Last Step in Volume
00846         G4bool exiting =  fLinearNavigator->ExitedMotherVolume();
00847         G4bool entering = fLinearNavigator->EnteredDaughterVolume();
00848 
00849         if( ! (exiting || entering) )
00850         {
00851             G4cout << " Transport> :  Proposed isLastStep= " << isLastStep
00852                    << " Exiting "  << fLinearNavigator->ExitedMotherVolume()
00853                    << " Entering " << fLinearNavigator->EnteredDaughterVolume()
00854                    << " Track position : " << track.GetPosition() /nanometer << " [nm]"
00855                    << G4endl;
00856             G4cout << " Track position : " << track.GetPosition() /nanometer << G4endl;
00857         }
00858 #endif
00859     }
00860     else                 // fGeometryLimitedStep  is false
00861     {
00862         // This serves only to move the Navigator's location
00863         //
00864         fLinearNavigator->LocateGlobalPointWithinVolume( track.GetPosition() ) ;
00865 
00866         // The value of the track's current Touchable is retained.
00867         // (and it must be correct because we must use it below to
00868         // overwrite the (unset) one in particle change)
00869         //  It must be fCurrentTouchable too ??
00870         //
00871         fParticleChange.SetTouchableHandle( track.GetTouchableHandle() ) ;
00872         retCurrentTouchable = track.GetTouchableHandle() ;
00873 
00874         isLastStep= false;
00875 #ifdef G4DEBUG_TRANSPORT
00876         //  Checking first implementation of flagging Last Step in Volume
00877         G4cout << " Transport> Proposed isLastStep= " << isLastStep
00878                << " Geometry did not limit step. Position : "
00879                << track.GetPosition()/ nanometer << G4endl;
00880 #endif
00881     }         // endif ( fGeometryLimitedStep )
00882 
00883     fParticleChange.ProposeLastStepInVolume(isLastStep);
00884 
00885     const G4VPhysicalVolume* pNewVol = retCurrentTouchable->GetVolume() ;
00886     const G4Material* pNewMaterial   = 0 ;
00887     const G4VSensitiveDetector* pNewSensitiveDetector   = 0 ;
00888 
00889     if( pNewVol != 0 )
00890     {
00891         pNewMaterial= pNewVol->GetLogicalVolume()->GetMaterial();
00892         pNewSensitiveDetector= pNewVol->GetLogicalVolume()->GetSensitiveDetector();
00893     }
00894 
00895     // ( <const_cast> pNewMaterial ) ;
00896     // ( <const_cast> pNewSensitiveDetector) ;
00897 
00898     fParticleChange.SetMaterialInTouchable( (G4Material *) pNewMaterial ) ;
00899     fParticleChange.SetSensitiveDetectorInTouchable( (G4VSensitiveDetector *) pNewSensitiveDetector ) ;
00900 
00901     const G4MaterialCutsCouple* pNewMaterialCutsCouple = 0;
00902     if( pNewVol != 0 )
00903     {
00904         pNewMaterialCutsCouple=pNewVol->GetLogicalVolume()->GetMaterialCutsCouple();
00905     }
00906 
00907     if( pNewVol!=0 && pNewMaterialCutsCouple!=0 && pNewMaterialCutsCouple->GetMaterial()!=pNewMaterial )
00908     {
00909         // for parametrized volume
00910         //
00911         pNewMaterialCutsCouple =
00912                 G4ProductionCutsTable::GetProductionCutsTable()
00913                 ->GetMaterialCutsCouple(pNewMaterial,
00914                                         pNewMaterialCutsCouple->GetProductionCuts());
00915     }
00916     fParticleChange.SetMaterialCutsCoupleInTouchable( pNewMaterialCutsCouple );
00917 
00918     // temporarily until Get/Set Material of ParticleChange,
00919     // and StepPoint can be made const.
00920     // Set the touchable in ParticleChange
00921     // this must always be done because the particle change always
00922     // uses this value to overwrite the current touchable pointer.
00923     //
00924     fParticleChange.SetTouchableHandle(retCurrentTouchable) ;
00925 
00926     return &fParticleChange ;
00927 }
00928 
00929 // New method takes over the responsibility to reset the state of G4Transportation
00930 //   object at the start of a new track or the resumption of a suspended track.
00931 
00932 void
00933 G4ITTransportation::StartTracking(G4Track* track)
00934 {
00935     G4VProcess::StartTracking(track);
00936     if(fInstantiateProcessState)
00937         G4VITProcess::fpState = new G4ITTransportationState();
00938     // Will set in the same time fTransportationState
00939 
00940     // The actions here are those that were taken in AlongStepGPIL
00941     //   when track.GetCurrentStepNumber()==1
00942 
00943     // reset safety value and center
00944     //
00945     //    State(fPreviousSafety)    = 0.0 ;
00946     //    State(fPreviousSftOrigin) = G4ThreeVector(0.,0.,0.) ;
00947 
00948     // reset looping counter -- for motion in field
00949     //    State(fNoLooperTrials)= 0;
00950     // Must clear this state .. else it depends on last track's value
00951     //  --> a better solution would set this from state of suspended track TODO ?
00952     // Was if( aTrack->GetCurrentStepNumber()==1 ) { .. }
00953 
00954     // ChordFinder reset internal state
00955     //
00956     if( DoesGlobalFieldExist() )
00957     {
00958         fFieldPropagator->ClearPropagatorState();
00959         // Resets all state of field propagator class (ONLY)
00960         //  including safety values (in case of overlaps and to wipe for first track).
00961 
00962         // G4ChordFinder* chordF= fFieldPropagator->GetChordFinder();
00963         // if( chordF ) chordF->ResetStepEstimate();
00964     }
00965 
00966     // Make sure to clear the chord finders of all fields (ie managers)
00967     static G4FieldManagerStore* fieldMgrStore= G4FieldManagerStore::GetInstance();
00968     fieldMgrStore->ClearAllChordFindersState();
00969 
00970     // Update the current touchable handle  (from the track's)
00971     //
00972     State(fCurrentTouchableHandle) = track->GetTouchableHandle();
00973 
00974     G4VITProcess::StartTracking(track);
00975 }
00976 
00977 #undef State

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