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

Generated on Mon May 27 17:50:02 2013 for Geant4 by  doxygen 1.4.7