G4CoupledTransportation.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: G4CoupledTransportation.cc 69887 2013-05-17 08:17:02Z gcosmo $
00028 //
00029 // ------------------------------------------------------------
00030 //  GEANT 4 class implementation
00031 // =======================================================================
00032 // Modified:
00033 //            13 May  2006, J. Apostolakis: Revised for parallel navigation (PathFinder)
00034 //            19 Jan  2006, P.MoraDeFreitas: Fix for suspended tracks (StartTracking)
00035 //            11 Aug  2004, M.Asai: Add G4VSensitiveDetector* for updating stepPoint.
00036 //            21 June 2003, J.Apostolakis: Calling field manager with 
00037 //                            track, to enable it to configure its accuracy
00038 //            13 May  2003, J.Apostolakis: Zero field areas now taken into
00039 //                            account correclty in all cases (thanks to W Pokorski).
00040 //            29 June 2001, J.Apostolakis, D.Cote-Ahern, P.Gumplinger: 
00041 //                          correction for spin tracking   
00042 //            20 Febr 2001, J.Apostolakis:  update for new FieldTrack
00043 //            22 Sept 2000, V.Grichine:     update of Kinetic Energy
00044 // Created:  19 March 1997, J. Apostolakis
00045 // =======================================================================
00046 
00047 #include "G4CoupledTransportation.hh"
00048 
00049 #include "G4PhysicalConstants.hh"
00050 #include "G4SystemOfUnits.hh"
00051 #include "G4TransportationProcessType.hh"
00052 #include "G4ProductionCutsTable.hh"
00053 #include "G4ParticleTable.hh"
00054 #include "G4ChordFinder.hh"
00055 #include "G4Field.hh"
00056 #include "G4FieldManagerStore.hh"
00057 
00058 class G4VSensitiveDetector;
00059 
00061 //
00062 // Constructor
00063 
00064 G4CoupledTransportation::G4CoupledTransportation( G4int verbosity )
00065   : G4VProcess( G4String("CoupledTransportation"), fTransportation ),
00066     fParticleIsLooping( false ),
00067     fPreviousSftOrigin( 0.,0.,0. ),
00068     fPreviousMassSafety( 0.0 ),
00069     fPreviousFullSafety( 0.0 ),
00070     fThreshold_Warning_Energy( 100 * MeV ),  
00071     fThreshold_Important_Energy( 250 * MeV ), 
00072     fThresholdTrials( 10 ), 
00073     fNoLooperTrials( 0 ),
00074     fSumEnergyKilled( 0.0 ), fMaxEnergyKilled( 0.0 ), 
00075     fUseMagneticMoment( false ),  
00076     fVerboseLevel( verbosity )
00077 {
00078   // set Process Sub Type
00079   SetProcessSubType(static_cast<G4int>(COUPLED_TRANSPORTATION));
00080 
00081   G4TransportationManager* transportMgr ; 
00082 
00083   transportMgr = G4TransportationManager::GetTransportationManager() ; 
00084 
00085   fMassNavigator = transportMgr->GetNavigatorForTracking() ; 
00086   fFieldPropagator = transportMgr->GetPropagatorInField() ;
00087   // fGlobalFieldMgr = transportMgr->GetFieldManager() ;
00088   fNavigatorId= transportMgr->ActivateNavigator( fMassNavigator ); 
00089   if( fVerboseLevel > 0 )
00090   {
00091     G4cout << " G4CoupledTransportation constructor: ----- " << G4endl;
00092     G4cout << " Verbose level is " << fVerboseLevel << G4endl;
00093     G4cout << " Navigator Id obtained in G4CoupledTransportation constructor " 
00094            << fNavigatorId << G4endl;
00095   }
00096   fPathFinder=  G4PathFinder::GetInstance(); 
00097   fpSafetyHelper = transportMgr->GetSafetyHelper();  // New 
00098 
00099   // Following assignment is to fix small memory leak from simple use of 'new'
00100   static G4TouchableHandle nullTouchableHandle;  // Points to (G4VTouchable*) 0
00101   fCurrentTouchableHandle = nullTouchableHandle; 
00102 
00103   fEndGlobalTimeComputed  = false;
00104   fCandidateEndGlobalTime = 0;
00105 }
00106 
00108 
00109 G4CoupledTransportation::~G4CoupledTransportation()
00110 {
00111   // fCurrentTouchableHandle is a data member - no deletion required
00112 
00113   if( (fVerboseLevel > 0) || (fSumEnergyKilled > 0.0 ) )
00114   { 
00115     G4cout << " G4CoupledTransportation: Statistics for looping particles " << G4endl;
00116     G4cout << "   Sum of energy of loopers killed: " <<  fSumEnergyKilled << G4endl;
00117     G4cout << "   Max energy of loopers killed: " <<  fMaxEnergyKilled << G4endl;
00118   } 
00119 }
00120 
00122 //
00123 // Responsibilities:
00124 //    Find whether the geometry limits the Step, and to what length
00125 //    Calculate the new value of the safety and return it.
00126 //    Store the final time, position and momentum.
00127 
00128 G4double G4CoupledTransportation::
00129 AlongStepGetPhysicalInteractionLength( const G4Track&  track,
00130                                              G4double, //  previousStepSize
00131                                              G4double  currentMinimumStep,
00132                                              G4double& proposedSafetyForStart,
00133                                              G4GPILSelection* selection )
00134 {
00135   G4double geometryStepLength; 
00136   G4double startMassSafety= 0.0;   //  estimated safety for start point (mass geometry)
00137   G4double startFullSafety= 0.0;   //  estimated safety for start point (all geometries)
00138   G4double safetyProposal= -1.0;   //  local copy of proposal 
00139 
00140   G4ThreeVector  EndUnitMomentum ;
00141   G4double       lengthAlongCurve=0.0 ;
00142  
00143   fParticleIsLooping = false ;
00144 
00145   // Initial actions moved to  StartTrack()   
00146   // --------------------------------------
00147   // Note: in case another process changes touchable handle
00148   //    it will be necessary to add here (for all steps)   
00149   // fCurrentTouchableHandle = aTrack->GetTouchableHandle();
00150 
00151   // GPILSelection is set to defaule value of CandidateForSelection
00152   // It is a return value
00153   //
00154   *selection = CandidateForSelection ;
00155 
00156   // Get initial Energy/Momentum of the track
00157   //
00158   const G4DynamicParticle*    pParticle  = track.GetDynamicParticle() ;
00159   const G4ParticleDefinition* pParticleDef   = pParticle->GetDefinition() ;
00160   G4ThreeVector startMomentumDir       = pParticle->GetMomentumDirection() ;
00161   G4ThreeVector startPosition          = track.GetPosition() ;
00162   G4VPhysicalVolume* currentVolume= track.GetVolume(); 
00163 
00164 #ifdef G4DEBUG_TRANSPORT
00165   if( fVerboseLevel > 1 )
00166   {
00167     G4cout << "G4CoupledTransportation::AlongStepGPIL> called in volume " 
00168            << currentVolume->GetName() << G4endl; 
00169   }
00170 #endif
00171   // G4double   theTime        = track.GetGlobalTime() ;
00172 
00173   // The Step Point safety can be limited by other geometries and/or the 
00174   // assumptions of any process - it's not always the geometrical safety.
00175   // We calculate the starting point's isotropic safety here.
00176   //
00177   G4ThreeVector OriginShift = startPosition - fPreviousSftOrigin ;
00178   G4double      MagSqShift  = OriginShift.mag2() ;
00179   startMassSafety = 0.0; 
00180   startFullSafety= 0.0; 
00181 
00182   //  Recall that FullSafety <= MassSafety 
00183   // Original: if( MagSqShift < sqr(fPreviousMassSafety) ) {
00184   if( MagSqShift < sqr(fPreviousFullSafety) )   // Revision proposed by Alex H, 2 Oct 07
00185   {
00186      G4double mag_shift= std::sqrt(MagSqShift); 
00187      startMassSafety = std::max( (fPreviousMassSafety - mag_shift), 0.0); 
00188      startFullSafety = std::max( (fPreviousFullSafety - mag_shift), 0.0);
00189        // Need to be consistent between full safety with Mass safety
00190        //   in order reproduce results in simple case  --> use same calculation method
00191 
00192      // Only compute full safety if massSafety > 0.  Else it remains 0
00193      //   startFullSafety = fPathFinder->ComputeSafety( startPosition ); 
00194   }
00195 
00196   // Is the particle charged or has it a magnetic moment?
00197   //
00198   G4double particleCharge = pParticle->GetCharge() ;
00199   G4double magneticMoment = pParticle->GetMagneticMoment() ;
00200   G4double       restMass = pParticleDef->GetPDGMass() ; 
00201 
00202   fMassGeometryLimitedStep = false ; //  Set default - alex
00203   fAnyGeometryLimitedStep = false; 
00204 
00205   // fEndGlobalTimeComputed = false ;
00206 
00207   // There is no need to locate the current volume. It is Done elsewhere:
00208   //   On track construction 
00209   //   By the tracking, after all AlongStepDoIts, in "Relocation"
00210 
00211   // Check if the particle has a force, EM or gravitational, exerted on it
00212   //
00213   G4FieldManager* fieldMgr=0;
00214   G4bool          fieldExertsForce = false ;
00215 
00216   G4bool gravityOn = false;
00217   const G4Field* ptrField= 0;
00218 
00219   fieldMgr = fFieldPropagator->FindAndSetFieldManager( track.GetVolume() );
00220   if( fieldMgr != 0 )
00221   {
00222      // Message the field Manager, to configure it for this track
00223      fieldMgr->ConfigureForTrack( &track );
00224      // Here it can transition from a null field-ptr to a finite field 
00225 
00226      // If the field manager has no field ptr, the field is zero 
00227      //     by definition ( = there is no field ! )
00228      ptrField= fieldMgr->GetDetectorField();
00229  
00230      if( ptrField != 0)
00231      { 
00232         gravityOn= ptrField->IsGravityActive();
00233         if(  (particleCharge != 0.0) 
00234              || (fUseMagneticMoment && (magneticMoment != 0.0) )
00235              || (gravityOn && (restMass != 0.0))
00236           )
00237         {
00238            fieldExertsForce = true;
00239         }
00240      }
00241   }
00242   G4double momentumMagnitude = pParticle->GetTotalMomentum() ;
00243  
00244   fFieldPropagator->SetChargeMomentumMass( particleCharge,    // in e+ units
00245                                            momentumMagnitude, // in Mev/c 
00246                                            restMass           ) ;  
00247   // This should be obsolete - the call to SetChargeAndMoments below should do the work
00248 
00249   G4ThreeVector spin        = track.GetPolarization() ;
00250   G4FieldTrack  theFieldTrack = G4FieldTrack( startPosition, 
00251                                             track.GetMomentumDirection(),
00252                                             0.0, 
00253                                             track.GetKineticEnergy(),
00254                                             restMass,
00255                                             0.0,    // UNUSED: track.GetVelocity(),
00256                                             track.GetGlobalTime(), // Lab.
00257                                             track.GetProperTime(), // Part.
00258                                             &spin                  ) ;
00259   theFieldTrack.SetChargeAndMoments( particleCharge ); // EM moments -- future extension 
00260 
00261   G4int stepNo= track.GetCurrentStepNumber(); 
00262 
00263   ELimited limitedStep; 
00264   G4FieldTrack endTrackState('a');  //  Default values
00265 
00266   fMassGeometryLimitedStep = false ;    //  default 
00267   fAnyGeometryLimitedStep  = false ;
00268   if( currentMinimumStep > 0 )
00269   {
00270       G4double newMassSafety= 0.0;     //  temp. for recalculation
00271 
00272       // Do the Transport in the field (non recti-linear)
00273       //
00274       lengthAlongCurve = fPathFinder->ComputeStep( theFieldTrack,
00275                                                    currentMinimumStep, 
00276                                                    fNavigatorId,
00277                                                    stepNo,
00278                                                    newMassSafety,
00279                                                    limitedStep,
00280                                                    endTrackState,
00281                                                    currentVolume ) ;
00282       // G4cout << " PathFinder ComputeStep returns " << lengthAlongCurve << G4endl; 
00283 
00284       G4double newFullSafety= fPathFinder->GetCurrentSafety();  
00285                // this was estimated already in step above
00286       // G4double newFullStep= fPathFinder->GetMinimumStep(); 
00287 
00288       if( limitedStep == kUnique || limitedStep == kSharedTransport )
00289       {
00290          fMassGeometryLimitedStep = true ;
00291       }
00292 
00293       fAnyGeometryLimitedStep = (fPathFinder->GetNumberGeometriesLimitingStep() != 0); 
00294 
00295 //#ifdef G4DEBUG_TRANSPORT
00296       if( fMassGeometryLimitedStep && !fAnyGeometryLimitedStep )
00297       {
00298          G4cerr << " Error in determining geometries limiting the step" << G4endl;
00299          G4cerr << "  Limiting:  mass=" << fMassGeometryLimitedStep
00300                 << " any= " << fAnyGeometryLimitedStep << G4endl;
00301          G4Exception("G4CoupledTransportation::AlongStepGetPhysicalInteractionLength()", 
00302                      "PathFinderConfused", FatalException, 
00303                      "Incompatible conditions - was limited by a geometry?");
00304       }
00305 //#endif
00306 
00307       // Other potential 
00308       // fAnyGeometryLimitedStep = newFullStep < currentMinimumStep; 
00309       //                                      ^^^ Not good enough; 
00310       //          Must compare with maximum requested step size
00311       //           (eg in case another process requested bigger, got this!)
00312 
00313       geometryStepLength = std::min( lengthAlongCurve, currentMinimumStep); 
00314 
00315       // Momentum:  Magnitude and direction can be changed too now ...
00316       //
00317       fMomentumChanged         = true ; 
00318       fTransportEndMomentumDir = endTrackState.GetMomentumDir() ;
00319 
00320       // Remember last safety origin & value.
00321       fPreviousSftOrigin  = startPosition ;
00322       fPreviousMassSafety = newMassSafety ;         
00323       fPreviousFullSafety = newFullSafety ; 
00324       // fpSafetyHelper->SetCurrentSafety( newFullSafety, startPosition);
00325 
00326 #ifdef G4DEBUG_TRANSPORT
00327       if( fVerboseLevel > 1 )
00328       {
00329         G4cout << "G4Transport:CompStep> " 
00330                << " called the pathfinder for a new step at " << startPosition
00331                << " and obtained step = " << lengthAlongCurve << G4endl;
00332         G4cout << "  New safety (preStep) = " << newMassSafety 
00333                << " versus precalculated = "  << startMassSafety << G4endl; 
00334       }
00335 #endif
00336 
00337       // Store as best estimate value
00338       startMassSafety    = newMassSafety ; 
00339       startFullSafety    = newFullSafety ; 
00340 
00341       // Get the End-Position and End-Momentum (Dir-ection)
00342       fTransportEndPosition = endTrackState.GetPosition() ;
00343       fTransportEndKineticEnergy  = endTrackState.GetKineticEnergy() ; 
00344   }
00345   else
00346   {
00347       geometryStepLength   = lengthAlongCurve= 0.0 ;
00348       fMomentumChanged         = false ; 
00349       // fMassGeometryLimitedStep = false ;   //  --- ???
00350       // fAnyGeometryLimitedStep = true;
00351       fTransportEndMomentumDir = track.GetMomentumDirection();
00352       fTransportEndKineticEnergy  = track.GetKineticEnergy();
00353 
00354       fTransportEndPosition = startPosition;
00355       // If the step length requested is 0, and we are on a boundary
00356       //   then a boundary will also limit the step.
00357       if( startMassSafety == 0.0 )
00358       {
00359          fMassGeometryLimitedStep = true ;
00360          fAnyGeometryLimitedStep = true;
00361       }
00362       //   TODO:  Add explicit logical status for being at a boundary
00363   }
00364   // G4FieldTrack aTrackState(endTrackState);  
00365 
00366   if( !fieldExertsForce ) 
00367   { 
00368       fParticleIsLooping         = false ; 
00369       fMomentumChanged           = false ; 
00370       fEndGlobalTimeComputed     = false ; 
00371       // G4cout << " global time is false " << G4endl; 
00372   } 
00373   else 
00374   { 
00375   
00376 #ifdef G4DEBUG_TRANSPORT
00377       if( fVerboseLevel > 1 )
00378       {
00379         G4cout << " G4CT::CS End Position = "  << fTransportEndPosition << G4endl; 
00380         G4cout << " G4CT::CS End Direction = " << fTransportEndMomentumDir << G4endl; 
00381       }
00382 #endif
00383       if( fFieldPropagator->GetCurrentFieldManager()->DoesFieldChangeEnergy() )
00384       {
00385           // If the field can change energy, then the time must be integrated
00386           //    - so this should have been updated
00387           //
00388           fCandidateEndGlobalTime   = endTrackState.GetLabTimeOfFlight();
00389           fEndGlobalTimeComputed    = true;
00390   
00391           // was ( fCandidateEndGlobalTime != track.GetGlobalTime() );
00392           // a cleaner way is to have FieldTrack knowing whether time is updated.
00393       }
00394       else
00395       {
00396           // The energy should be unchanged by field transport,
00397           //    - so the time changed will be calculated elsewhere
00398           //
00399           fEndGlobalTimeComputed = false;
00400   
00401           // Check that the integration preserved the energy 
00402           //     -  and if not correct this!
00403           G4double  startEnergy= track.GetKineticEnergy();
00404           G4double  endEnergy= fTransportEndKineticEnergy; 
00405       
00406           static G4int no_inexact_steps=0; // , no_large_ediff;
00407           G4double absEdiff = std::fabs(startEnergy- endEnergy);
00408           if( absEdiff > perMillion * endEnergy )
00409           {
00410             no_inexact_steps++;
00411             // Possible statistics keeping here ...
00412           }
00413 #ifdef G4VERBOSE
00414           if( (fVerboseLevel > 1) && ( absEdiff > perThousand * endEnergy) )
00415           {
00416             ReportInexactEnergy(startEnergy, endEnergy); 
00417           }  // end of if (fVerboseLevel)
00418 #endif
00419           // Correct the energy for fields that conserve it
00420           //  This - hides the integration error
00421           //       - but gives a better physical answer
00422           fTransportEndKineticEnergy= track.GetKineticEnergy(); 
00423       }
00424   }
00425 
00426   endpointDistance   = (fTransportEndPosition - startPosition).mag() ;
00427   fParticleIsLooping = fFieldPropagator->IsParticleLooping() ;
00428 
00429   fTransportEndSpin = endTrackState.GetSpin();
00430 
00431   // Calculate the safety 
00432   safetyProposal= startFullSafety;   // used to be startMassSafety
00433      // Changed to accomodate processes that cannot update the safety -- JA 22 Nov 06
00434 
00435   // Update safety for the end-point, if becomes negative at the end-point.
00436 
00437   if(   (startFullSafety < endpointDistance ) 
00438         && ( particleCharge != 0.0 ) )        //  Only needed to prepare for Mult Scat.
00439    //   && !fAnyGeometryLimitedStep )          // To-Try:  No safety update if at a boundary
00440   {
00441       G4double endFullSafety =
00442         fPathFinder->ComputeSafety( fTransportEndPosition); 
00443         // Expected mission -- only mass geometry's safety
00444         //   fMassNavigator->ComputeSafety( fTransportEndPosition) ;
00445         // Yet discrete processes only have poststep -- and this cannot 
00446         //   currently revise the safety  
00447         //   ==> so we use the all-geometry safety as a precaution
00448 
00449       fpSafetyHelper->SetCurrentSafety( endFullSafety, fTransportEndPosition);
00450         // Pushing safety to Helper avoids recalculation at this point
00451 
00452       G4ThreeVector centerPt= G4ThreeVector(0.0, 0.0, 0.0);  // Used for return value
00453       G4double endMassSafety= fPathFinder->ObtainSafety( fNavigatorId, centerPt); 
00454         //  Retrieves the mass value from PathFinder (it calculated it)
00455 
00456       fPreviousMassSafety = endMassSafety ; 
00457       fPreviousFullSafety = endFullSafety; 
00458       fPreviousSftOrigin = fTransportEndPosition ;
00459 
00460       // The convention (Stepping Manager's) is safety from the start point
00461       //
00462       safetyProposal = endFullSafety + endpointDistance;
00463           //  --> was endMassSafety
00464       // Changed to accomodate processes that cannot update the safety -- JA 22 Nov 06
00465 
00466       // #define G4DEBUG_TRANSPORT 1
00467 
00468 #ifdef G4DEBUG_TRANSPORT 
00469       G4int prec= G4cout.precision(12) ;
00470       G4cout << "***Transportation::AlongStepGPIL ** " << G4endl  ;
00471       G4cout << "  Revised Safety at endpoint "  << fTransportEndPosition
00472              << "   give safety values: Mass= " << endMassSafety 
00473              << "  All= " << endFullSafety << G4endl ; 
00474       G4cout << "  Adding endpoint distance " << endpointDistance 
00475              << "   to obtain pseudo-safety= " << safetyProposal << G4endl ; 
00476       G4cout.precision(prec); 
00477   }  
00478   else
00479   {
00480       G4int prec= G4cout.precision(12) ;
00481       G4cout << "***Transportation::AlongStepGPIL ** " << G4endl  ;
00482       G4cout << "  Quick Safety estimate at endpoint "  << fTransportEndPosition
00483              << "   gives safety endpoint value = " << startFullSafety - endpointDistance
00484              << "  using start-point value " << startFullSafety 
00485              << "  and endpointDistance " << endpointDistance << G4endl; 
00486       G4cout.precision(prec); 
00487 #endif
00488   }          
00489 
00490   proposedSafetyForStart= safetyProposal; 
00491   fParticleChange.ProposeTrueStepLength(geometryStepLength) ;
00492 
00493   return geometryStepLength ;
00494 }
00495 
00497 
00498 G4VParticleChange*
00499 G4CoupledTransportation::AlongStepDoIt( const G4Track& track,
00500                                         const G4Step&  stepData )
00501 {
00502   static G4int noCalls=0;
00503   noCalls++;
00504 
00505   fParticleChange.Initialize(track) ;
00506       // sets all its members to the value of corresponding members in G4Track
00507 
00508   //  Code specific for Transport
00509   //
00510   fParticleChange.ProposePosition(fTransportEndPosition) ;
00511   // G4cout << " G4CoupledTransportation::AlongStepDoIt" 
00512   //     << " proposes position = " << fTransportEndPosition  
00513   //     << " and end momentum direction  = " << fTransportEndMomentumDir <<  G4endl;
00514   fParticleChange.ProposeMomentumDirection(fTransportEndMomentumDir) ;
00515   fParticleChange.ProposeEnergy(fTransportEndKineticEnergy) ;
00516   fParticleChange.SetMomentumChanged(fMomentumChanged) ;
00517 
00518   fParticleChange.ProposePolarization(fTransportEndSpin);
00519   
00520   G4double deltaTime = 0.0 ;
00521 
00522   // Calculate  Lab Time of Flight (ONLY if field Equations used it!)
00523      // G4double endTime   = fCandidateEndGlobalTime;
00524      // G4double delta_time = endTime - startTime;
00525 
00526   G4double startTime = track.GetGlobalTime() ;
00527   
00528   if (!fEndGlobalTimeComputed)
00529   {
00530      G4double finalInverseVel= DBL_MAX, initialInverseVel=DBL_MAX; 
00531 
00532      // The time was not integrated .. make the best estimate possible
00533      //
00534      G4double finalVelocity   = track.GetVelocity() ;
00535      if( finalVelocity > 0.0 ) { finalInverseVel= 1.0 / finalVelocity; }
00536      G4double initialVelocity = stepData.GetPreStepPoint()->GetVelocity() ;
00537      if( initialVelocity > 0.0 ) { initialInverseVel= 1.0 / initialVelocity; }
00538      G4double stepLength      = track.GetStepLength() ;
00539 
00540      if (finalVelocity > 0.0)
00541      {
00542         // deltaTime = stepLength/finalVelocity ;
00543         G4double meanInverseVelocity = 0.5 * ( initialInverseVel + finalInverseVel );
00544         deltaTime = stepLength * meanInverseVelocity ;
00545         // G4cout << " dt = s * mean(1/v) , with " << "  s = " << stepLength
00546         //     << "  mean(1/v)= "  << meanInverseVelocity << G4endl;
00547      }
00548      else
00549      {
00550         deltaTime = stepLength * initialInverseVel ;
00551         // G4cout << " dt = s / initV "  << "  s = "   << stepLength
00552         //        << " 1 / initV= " << initialInverseVel << G4endl; 
00553      }  //  Could do with better estimate for final step (finalVelocity = 0) ?
00554 
00555      fCandidateEndGlobalTime   = startTime + deltaTime ;
00556      fParticleChange.ProposeLocalTime(  track.GetLocalTime() + deltaTime) ;
00557 
00558      // G4cout << " Calculated global time from start = " << startTime << " and "
00559      //        << " delta time = " << deltaTime << G4endl;
00560   }
00561   else
00562   {
00563      deltaTime = fCandidateEndGlobalTime - startTime ;
00564      fParticleChange.ProposeGlobalTime( fCandidateEndGlobalTime ) ;
00565      // G4cout << " Calculated global time from candidate end time = "
00566      //    << fCandidateEndGlobalTime << " and start time = " << startTime << G4endl;
00567   }
00568 
00569   // G4cout << " G4CoupledTransportation::AlongStepDoIt  "
00570   // << " flag whether computed time = " << fEndGlobalTimeComputed  << " and " 
00571   // << " is proposes end time " << fCandidateEndGlobalTime << G4endl; 
00572 
00573   // Now Correct by Lorentz factor to get "proper" deltaTime
00574   
00575   G4double  restMass       = track.GetDynamicParticle()->GetMass() ;
00576   G4double deltaProperTime = deltaTime*( restMass/track.GetTotalEnergy() ) ;
00577 
00578   fParticleChange.ProposeProperTime(track.GetProperTime() + deltaProperTime) ;
00579   //fParticleChange. ProposeTrueStepLength( track.GetStepLength() ) ;
00580 
00581   // If the particle is caught looping or is stuck (in very difficult
00582   // boundaries) in a magnetic field (doing many steps) 
00583   //   THEN this kills it ...
00584   //
00585   if ( fParticleIsLooping )
00586   {
00587      G4double endEnergy= fTransportEndKineticEnergy;
00588 
00589      if( (endEnergy < fThreshold_Important_Energy) 
00590           || (fNoLooperTrials >= fThresholdTrials ) )
00591      {
00592         // Kill the looping particle 
00593         //
00594         fParticleChange.ProposeTrackStatus( fStopAndKill )  ;
00595 
00596         // 'Bare' statistics
00597         fSumEnergyKilled += endEnergy; 
00598         if( endEnergy > fMaxEnergyKilled) { fMaxEnergyKilled= endEnergy; }
00599 
00600 #ifdef G4VERBOSE
00601         if((fVerboseLevel > 1) && ( endEnergy > fThreshold_Warning_Energy ))
00602         { 
00603           G4cout << " G4CoupledTransportation is killing track that is looping or stuck " << G4endl
00604                  << "   This track has " << track.GetKineticEnergy() / MeV
00605                  << " MeV energy." << G4endl;
00606         }
00607         if( fVerboseLevel > 0 )
00608         { 
00609           G4cout << "   Steps by this track: " << track.GetCurrentStepNumber() << G4endl;
00610         }
00611 #endif
00612         fNoLooperTrials=0; 
00613       }
00614       else
00615       { 
00616         fNoLooperTrials ++; 
00617 #ifdef G4VERBOSE
00618         if( (fVerboseLevel > 2) )
00619         {
00620           G4cout << "  ** G4CoupledTransportation::AlongStepDoIt(): Particle looping -  " << G4endl
00621                  << "   Number of consecutive problem step (this track) = " << fNoLooperTrials << G4endl
00622                  << "   Steps by this track: " << track.GetCurrentStepNumber() << G4endl
00623                  << "   Total no of calls to this method (all tracks) = " << noCalls << G4endl;
00624         }
00625 #endif
00626       }
00627   }
00628   else
00629   { 
00630       fNoLooperTrials=0; 
00631   }
00632 
00633   // Another (sometimes better way) is to use a user-limit maximum Step size
00634   // to alleviate this problem .. 
00635 
00636   // Add smooth curved trajectories to particle-change
00637   //
00638   // fParticleChange.SetPointerToVectorOfAuxiliaryPoints
00639   //   (fFieldPropagator->GimmeTrajectoryVectorAndForgetIt() );
00640 
00641   return &fParticleChange ;
00642 }
00643 
00645 //
00646 //  This ensures that the PostStep action is always called,
00647 //  so that it can do the relocation if it is needed.
00648 // 
00649 
00650 G4double G4CoupledTransportation::
00651 PostStepGetPhysicalInteractionLength( const G4Track&,
00652                                             G4double, // previousStepSize
00653                                             G4ForceCondition* pForceCond )
00654 { 
00655   // Must act as PostStep action -- to relocate particle
00656   *pForceCond = Forced ;    
00657   return DBL_MAX ;
00658 }
00659 
00660 void G4CoupledTransportation::
00661 ReportMove( G4ThreeVector OldVector, G4ThreeVector NewVector, const G4String& Quantity )
00662 {
00663     G4ThreeVector moveVec = ( NewVector - OldVector );
00664 
00665     G4cerr << G4endl
00666            << "**************************************************************" << G4endl;
00667     G4cerr << "Endpoint has moved between value expected from TransportEndPosition "
00668            << " and value from Track in PostStepDoIt. " << G4endl
00669            << "Change of " << Quantity << " is " << moveVec.mag() / mm << " mm long, "
00670            << " and its vector is " << (1.0/mm) * moveVec << " mm " << G4endl
00671            << "Endpoint of ComputeStep was " << OldVector
00672            << " and current position to locate is " << NewVector << G4endl;
00673 }
00674 
00676 
00677 G4VParticleChange* G4CoupledTransportation::PostStepDoIt( const G4Track& track,
00678                                                           const G4Step& )
00679 {
00680   G4TouchableHandle retCurrentTouchable ;   // The one to return
00681 
00682   // Initialize ParticleChange  (by setting all its members equal
00683   //                             to corresponding members in G4Track)
00684   // fParticleChange.Initialize(track) ;  // To initialise TouchableChange
00685 
00686   fParticleChange.ProposeTrackStatus(track.GetTrackStatus()) ;
00687 
00688   // Check that the end position and direction are preserved 
00689   // since call to AlongStepDoIt
00690 
00691 #ifdef G4DEBUG_TRANSPORT
00692   if( ( fVerboseLevel > 0 ) && ((fTransportEndPosition - track.GetPosition()).mag2() >= 1.0e-16) )
00693   {
00694      ReportMove( track.GetPosition(), fTransportEndPosition, "End of Step Position" ); 
00695      G4cerr << " Problem in G4CoupledTransportation::PostStepDoIt " << G4endl; 
00696   }
00697 
00698   // If the Step was determined by the volume boundary, relocate the particle
00699   // The pathFinder will know that the geometry limited the step (!?)
00700 
00701   if( fVerboseLevel > 0 )
00702   {
00703      G4cout << " Calling PathFinder::Locate() from " 
00704             << " G4CoupledTransportation::PostStepDoIt " << G4endl;
00705      G4cout << "  fAnyGeometryLimitedStep is " << fAnyGeometryLimitedStep << G4endl;
00706 
00707   }
00708 #endif
00709 
00710   if(fAnyGeometryLimitedStep)
00711   {  
00712     fPathFinder->Locate( track.GetPosition(), 
00713                        track.GetMomentumDirection(),
00714                        true); 
00715 
00716     // fCurrentTouchable will now become the previous touchable, 
00717     // and what was the previous will be freed.
00718     // (Needed because the preStepPoint can point to the previous touchable)
00719 
00720     fCurrentTouchableHandle= 
00721       fPathFinder->CreateTouchableHandle( fNavigatorId );
00722 
00723 #ifdef G4DEBUG_TRANSPORT
00724     if( fVerboseLevel > 0 )
00725     {
00726       G4cout << "G4CoupledTransportation::PostStepDoIt --- fNavigatorId = " 
00727              << fNavigatorId << G4endl;
00728     }
00729     if( fVerboseLevel > 1 )
00730     {
00731        G4VPhysicalVolume* vol= fCurrentTouchableHandle->GetVolume(); 
00732        G4cout << "CHECK !!!!!!!!!!! fCurrentTouchableHandle->GetVolume() = " << vol;
00733        if( vol ) { G4cout << "Name=" << vol->GetName(); }
00734        G4cout << G4endl;
00735     }
00736 #endif
00737 
00738     // Check whether the particle is out of the world volume 
00739     // If so it has exited and must be killed.
00740     //
00741     if( fCurrentTouchableHandle->GetVolume() == 0 )
00742     {
00743        fParticleChange.ProposeTrackStatus( fStopAndKill ) ;
00744     }
00745     retCurrentTouchable = fCurrentTouchableHandle ;
00746     // fParticleChange.SetTouchableHandle( fCurrentTouchableHandle ) ;
00747 
00748     // Notify particle change that this is last step in volume
00749     fParticleChange.ProposeLastStepInVolume(true);
00750   }
00751   else                 // fAnyGeometryLimitedStep  is false
00752   { 
00753 #ifdef G4DEBUG_TRANSPORT
00754     if( fVerboseLevel > 1 )
00755     {
00756        G4cout << "G4CoupledTransportation::PostStepDoIt -- "
00757               << " fAnyGeometryLimitedStep  = " << fAnyGeometryLimitedStep  
00758               << " must be false " << G4endl;
00759     }
00760 #endif
00761     // This serves only to move each of the Navigator's location
00762     //
00763     // fLinearNavigator->LocateGlobalPointWithinVolume( track.GetPosition() ) ;
00764 
00765     // G4cout << "G4CoupledTransportation calling PathFinder::ReLocate() " << G4endl;
00766     fPathFinder->ReLocate( track.GetPosition() );
00767                            // track.GetMomentumDirection() ); 
00768 
00769     // Keep the value of the track's current Touchable is retained,
00770     //  and use it to overwrite the (unset) one in particle change.
00771     // Expect this must be fCurrentTouchable too
00772     //   - could it be different, eg at the start of a step ?
00773     //
00774     retCurrentTouchable = track.GetTouchableHandle() ;
00775     // fParticleChange.SetTouchableHandle( track.GetTouchableHandle() ) ;
00776 
00777     // Have not reached a boundary
00778     fParticleChange.ProposeLastStepInVolume(false);
00779   }         // endif ( fAnyGeometryLimitedStep ) 
00780 
00781   const G4VPhysicalVolume* pNewVol = retCurrentTouchable->GetVolume() ;
00782   const G4Material* pNewMaterial   = 0 ;
00783   const G4VSensitiveDetector* pNewSensitiveDetector   = 0 ;
00784                                                                                        
00785   if( pNewVol != 0 )
00786   {
00787     pNewMaterial= pNewVol->GetLogicalVolume()->GetMaterial();
00788     pNewSensitiveDetector= pNewVol->GetLogicalVolume()->GetSensitiveDetector();
00789   }
00790 
00791   // ( const_cast<G4Material *> pNewMaterial ) ;
00792   // ( const_cast<G4VSensitiveDetetor *> pNewSensitiveDetector) ;
00793 
00794   fParticleChange.SetMaterialInTouchable( (G4Material *) pNewMaterial ) ;
00795   fParticleChange.SetSensitiveDetectorInTouchable( (G4VSensitiveDetector *) pNewSensitiveDetector ) ;
00796              // "temporarily" until Get/Set Material of ParticleChange, 
00797              // and StepPoint can be made const. 
00798 
00799   const G4MaterialCutsCouple* pNewMaterialCutsCouple = 0;
00800   if( pNewVol != 0 )
00801   {
00802     pNewMaterialCutsCouple=pNewVol->GetLogicalVolume()->GetMaterialCutsCouple();
00803     if( pNewMaterialCutsCouple!=0 
00804         && pNewMaterialCutsCouple->GetMaterial()!=pNewMaterial )
00805       {
00806         // for parametrized volume
00807         //
00808         pNewMaterialCutsCouple =
00809           G4ProductionCutsTable::GetProductionCutsTable()
00810                        ->GetMaterialCutsCouple(pNewMaterial,
00811                                                pNewMaterialCutsCouple->GetProductionCuts());
00812       }
00813   }
00814   fParticleChange.SetMaterialCutsCoupleInTouchable( pNewMaterialCutsCouple );
00815 
00816   // Must always set the touchable in ParticleChange, whether relocated or not
00817   fParticleChange.SetTouchableHandle(retCurrentTouchable) ;
00818 
00819   return &fParticleChange ;
00820 }
00821 
00822 // New method takes over the responsibility to reset the state of 
00823 //   G4CoupledTransportation object:
00824 //      - at the start of a new track,  and
00825 //      - on the resumption of a suspended track. 
00826 
00827 void 
00828 G4CoupledTransportation::StartTracking(G4Track* aTrack)
00829 {
00830 
00831   G4TransportationManager* transportMgr =
00832     G4TransportationManager::GetTransportationManager();
00833 
00834   // G4VProcess::StartTracking(aTrack);
00835 
00836   //  The 'initialising' actions
00837   //     once taken in AlongStepGPIL -- if ( track.GetCurrentStepNumber()==1 )
00838 
00839   // fStartedNewTrack= true; 
00840 
00841   fMassNavigator = transportMgr->GetNavigatorForTracking() ; 
00842   fNavigatorId= transportMgr->ActivateNavigator( fMassNavigator );  // Confirm it!
00843 
00844   // if( fVerboseLevel > 1 ){
00845   //  G4cout << " Navigator Id obtained in StartTracking " << fNavigatorId << G4endl;
00846   // }
00847   G4ThreeVector position = aTrack->GetPosition(); 
00848   G4ThreeVector direction = aTrack->GetMomentumDirection();
00849 
00850   // if( fVerboseLevel > 1 ){
00851   //   G4cout << " Calling PathFinder::PrepareNewTrack from    " 
00852   //   << " G4CoupledTransportation::StartTracking -- which calls Locate()" << G4endl;
00853   // }
00854   fPathFinder->PrepareNewTrack( position, direction); 
00855   // This implies a call to fPathFinder->Locate( position, direction ); 
00856 
00857   // Global field, if any, must exist before tracking is started
00858   fGlobalFieldExists= DoesGlobalFieldExist(); 
00859   // reset safety value and center
00860   //
00861   fPreviousMassSafety  = 0.0 ; 
00862   fPreviousFullSafety  = 0.0 ; 
00863   fPreviousSftOrigin = G4ThreeVector(0.,0.,0.) ;
00864   
00865   // reset looping counter -- for motion in field  
00866   fNoLooperTrials= 0; 
00867   // Must clear this state .. else it depends on last track's value
00868   //  --> a better solution would set this from state of suspended track TODO ? 
00869   // Was if( aTrack->GetCurrentStepNumber()==1 ) { .. }
00870 
00871   // ChordFinder reset internal state
00872   //
00873   if( fGlobalFieldExists )
00874   {
00875      fFieldPropagator->ClearPropagatorState();   
00876        // Resets safety values, in case of overlaps.  
00877 
00878      G4ChordFinder* chordF= fFieldPropagator->GetChordFinder();
00879      if( chordF )  { chordF->ResetStepEstimate(); }
00880   }
00881 
00882   // Clear the chord finders of all fields (ie managers) derived objects
00883   //
00884   G4FieldManagerStore* fieldMgrStore = G4FieldManagerStore::GetInstance();
00885   fieldMgrStore->ClearAllChordFindersState(); 
00886 
00887 #ifdef G4DEBUG_TRANSPORT
00888   if( fVerboseLevel > 1 )
00889   {
00890     G4cout << " Returning touchable handle " << fCurrentTouchableHandle << G4endl;
00891   }
00892 #endif
00893 
00894   // Update the current touchable handle  (from the track's)
00895   //
00896   fCurrentTouchableHandle = aTrack->GetTouchableHandle();  
00897 }
00898 
00899 void 
00900 G4CoupledTransportation::EndTracking()
00901 {
00902   G4TransportationManager::GetTransportationManager()->InactivateAll();
00903 }
00904 
00905 void
00906 G4CoupledTransportation::
00907 ReportInexactEnergy(G4double startEnergy, G4double endEnergy)
00908 {
00909   static G4int no_warnings= 0, warnModulo=1,  moduloFactor= 10, no_large_ediff= 0; 
00910 
00911   if( std::fabs(startEnergy- endEnergy) > perThousand * endEnergy )
00912   {
00913     no_large_ediff ++;
00914     if( (no_large_ediff% warnModulo) == 0 )
00915     {
00916       no_warnings++;
00917       G4cout << "WARNING - G4CoupledTransportation::AlongStepGetPIL() " 
00918              << "   Energy change in Step is above 1^-3 relative value. " << G4endl
00919              << "   Relative change in 'tracking' step = " 
00920              << std::setw(15) << (endEnergy-startEnergy)/startEnergy << G4endl
00921              << "     Starting E= " << std::setw(12) << startEnergy / MeV << " MeV " << G4endl
00922              << "     Ending   E= " << std::setw(12) << endEnergy   / MeV << " MeV " << G4endl;       
00923       G4cout << " Energy has been corrected -- however, review"
00924              << " field propagation parameters for accuracy."  << G4endl;
00925       if( (fVerboseLevel > 2 ) || (no_warnings<4) || (no_large_ediff == warnModulo * moduloFactor) )
00926       {
00927         G4cout << " These include EpsilonStepMax(/Min) in G4FieldManager "
00928                << " which determine fractional error per step for integrated quantities. " << G4endl
00929                << " Note also the influence of the permitted number of integration steps."
00930                << G4endl;
00931       }
00932       G4cerr << "ERROR - G4CoupledTransportation::AlongStepGetPIL()" << G4endl
00933              << "        Bad 'endpoint'. Energy change detected"
00934              << " and corrected. " 
00935              << " Has occurred already "
00936              << no_large_ediff << " times." << G4endl;
00937       if( no_large_ediff == warnModulo * moduloFactor )
00938       {
00939         warnModulo *= moduloFactor;
00940       }
00941     }
00942   }
00943 }

Generated on Mon May 27 17:47:58 2013 for Geant4 by  doxygen 1.4.7