Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/examples/extended/exoticphysics/monopole/src/G4MonopoleTransportation.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /examples/extended/exoticphysics/monopole/src/G4MonopoleTransportation.cc (Version 11.3.0) and /examples/extended/exoticphysics/monopole/src/G4MonopoleTransportation.cc (Version 10.1.p3)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                         3 // * License and Disclaimer                                           *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th      5 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
  6 // * the Geant4 Collaboration.  It is provided      6 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
  7 // * conditions of the Geant4 Software License      7 // * conditions of the Geant4 Software License,  included in the file *
  8 // * LICENSE and available at  http://cern.ch/      8 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
  9 // * include a list of copyright holders.           9 // * include a list of copyright holders.                             *
 10 // *                                               10 // *                                                                  *
 11 // * Neither the authors of this software syst     11 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     12 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     13 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     14 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file      15 // * use.  Please see the license in the file  LICENSE  and URL above *
 16 // * for the full disclaimer and the limitatio     16 // * for the full disclaimer and the limitation of liability.         *
 17 // *                                               17 // *                                                                  *
 18 // * This  code  implementation is the result      18 // * This  code  implementation is the result of  the  scientific and *
 19 // * technical work of the GEANT4 collaboratio     19 // * technical work of the GEANT4 collaboration.                      *
 20 // * By using,  copying,  modifying or  distri     20 // * By using,  copying,  modifying or  distributing the software (or *
 21 // * any work based  on the software)  you  ag     21 // * any work based  on the software)  you  agree  to acknowledge its *
 22 // * use  in  resulting  scientific  publicati     22 // * use  in  resulting  scientific  publications,  and indicate your *
 23 // * acceptance of all terms of the Geant4 Sof     23 // * acceptance of all terms of the Geant4 Software license.          *
 24 // *******************************************     24 // ********************************************************************
 25 //                                                 25 //
 26 /// \file exoticphysics/monopole/src/G4Monopol     26 /// \file exoticphysics/monopole/src/G4MonopoleTransportation.cc
 27 /// \brief Implementation of the G4MonopoleTra     27 /// \brief Implementation of the G4MonopoleTransportation class
 28 //                                                 28 //
                                                   >>  29 // $Id: G4MonopoleTransportation.cc 84606 2014-10-17 07:50:04Z gcosmo $
 29 //                                                 30 //
 30 //....oooOO0OOooo........oooOO0OOooo........oo     31 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 31 //....oooOO0OOooo........oooOO0OOooo........oo     32 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 32 //                                                 33 //
 33 // This class is a process responsible for the <<  34 // This class is a process responsible for the transportation of 
 34 // magnetic monopoles, ie the geometrical prop <<  35 // magnetic monopoles, ie the geometrical propagation that encounters the 
 35 // geometrical sub-volumes of the detectors.       36 // geometrical sub-volumes of the detectors.
 36 //                                                 37 //
 37 // For monopoles, uses a different equation of     38 // For monopoles, uses a different equation of motion and ignores energy
 38 // conservation.                               <<  39 // conservation. 
 39 //                                                 40 //
 40                                                    41 
 41 // ===========================================     42 // =======================================================================
 42 // Created:  3 May 2010, J. Apostolakis, B. Bo     43 // Created:  3 May 2010, J. Apostolakis, B. Bozsogi
 43 // ===========================================     44 // =======================================================================
 44                                                    45 
 45 #include "G4MonopoleTransportation.hh"             46 #include "G4MonopoleTransportation.hh"
 46                                                <<  47 #include "G4ProductionCutsTable.hh"
 47 #include "DetectorConstruction.hh"             <<  48 #include "G4ParticleTable.hh"
 48                                                << 
 49 #include "G4ChordFinder.hh"                        49 #include "G4ChordFinder.hh"
                                                   >>  50 #include "G4SafetyHelper.hh"
 50 #include "G4FieldManagerStore.hh"                  51 #include "G4FieldManagerStore.hh"
 51 #include "G4Monopole.hh"                           52 #include "G4Monopole.hh"
 52 #include "G4ParticleTable.hh"                  << 
 53 #include "G4ProductionCutsTable.hh"            << 
 54 #include "G4RunManager.hh"                     << 
 55 #include "G4SafetyHelper.hh"                   << 
 56 #include "G4SystemOfUnits.hh"                  << 
 57 #include "G4TransportationProcessType.hh"          53 #include "G4TransportationProcessType.hh"
                                                   >>  54 #include "G4SystemOfUnits.hh"
 58                                                    55 
 59 class G4VSensitiveDetector;                        56 class G4VSensitiveDetector;
 60                                                    57 
 61 //....oooOO0OOooo........oooOO0OOooo........oo     58 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
 62                                                    59 
 63 G4MonopoleTransportation::G4MonopoleTransporta <<  60 G4MonopoleTransportation::G4MonopoleTransportation( const G4Monopole* mpl,
 64   : G4VProcess(G4String("MonopoleTransportatio <<  61                                                     G4int verb)
                                                   >>  62   : G4VProcess( G4String("MonopoleTransportation"), fTransportation ),
 65     fParticleDef(mpl),                             63     fParticleDef(mpl),
 66     fMagSetup(0),                                  64     fMagSetup(0),
 67     fLinearNavigator(0),                           65     fLinearNavigator(0),
 68     fFieldPropagator(0),                           66     fFieldPropagator(0),
 69     fParticleIsLooping(false),                 <<  67     fParticleIsLooping( false ),
 70     fPreviousSftOrigin(0., 0., 0.),            <<  68     fPreviousSftOrigin (0.,0.,0.),
 71     fPreviousSafety(0.0),                      <<  69     fPreviousSafety    ( 0.0 ),
 72     fThreshold_Warning_Energy(100 * MeV),      <<  70     fThreshold_Warning_Energy( 100 * MeV ),  
 73     fThreshold_Important_Energy(250 * MeV),    <<  71     fThreshold_Important_Energy( 250 * MeV ), 
 74     fThresholdTrials(10),                      <<  72     fThresholdTrials( 10 ), 
 75     // fUnimportant_Energy( 1 * MeV ),         <<  73     // fUnimportant_Energy( 1 * MeV ), 
 76     fNoLooperTrials(0),                            74     fNoLooperTrials(0),
 77     fSumEnergyKilled(0.0),                     <<  75     fSumEnergyKilled( 0.0 ), fMaxEnergyKilled( 0.0 ), 
 78     fMaxEnergyKilled(0.0),                     <<  76     fShortStepOptimisation(false),    // Old default: true (=fast short steps)
 79     fShortStepOptimisation(false),  // Old def <<  77     fpSafetyHelper(0)
 80     fpSafetyHelper(0),                         << 
 81     noCalls(0)                                 << 
 82 {                                                  78 {
 83   verboseLevel = verb;                             79   verboseLevel = verb;
 84                                                    80 
 85   // set Process Sub Type                          81   // set Process Sub Type
 86   SetProcessSubType(TRANSPORTATION);           <<  82   SetProcessSubType(TRANSPORTATION);  
 87                                                    83 
 88   // Do not finalize the G4MonopoleTransportat <<  84   fMagSetup = G4MonopoleFieldSetup::GetMonopoleFieldSetup();
 89   if (G4Threading::IsMasterThread() && G4Threa <<  85   
 90     return;                                    <<  86   G4TransportationManager* transportMgr ; 
 91   }                                            << 
 92                                                << 
 93   const DetectorConstruction* detector = stati << 
 94     G4RunManager::GetRunManager()->GetUserDete << 
 95                                                    87 
 96   fMagSetup = detector->GetMonopoleFieldSetup( <<  88   transportMgr = G4TransportationManager::GetTransportationManager() ; 
 97                                                    89 
 98   G4TransportationManager* transportMgr = G4Tr <<  90   fLinearNavigator = transportMgr->GetNavigatorForTracking() ; 
 99                                                    91 
100   fLinearNavigator = transportMgr->GetNavigato <<  92   // fGlobalFieldMgr = transportMgr->GetFieldManager() ;
101                                                    93 
102   fFieldPropagator = transportMgr->GetPropagat <<  94   fFieldPropagator = transportMgr->GetPropagatorInField() ;
103   fpSafetyHelper = transportMgr->GetSafetyHelp << 
104                                                    95 
105   // New                                       <<  96   fpSafetyHelper =   transportMgr->GetSafetyHelper();  // New 
106                                                    97 
107   // Cannot determine whether a field exists h     98   // Cannot determine whether a field exists here,
108   //  because it would only work if the field  <<  99   //  because it would only work if the field manager has informed 
109   //  about the detector's field before this t << 100   //  about the detector's field before this transportation process 
110   //  is constructed.                             101   //  is constructed.
111   // Instead later the method DoesGlobalFieldE    102   // Instead later the method DoesGlobalFieldExist() is called
112   fCurrentTouchableHandle = nullptr;           << 
113                                                   103 
114   fEndGlobalTimeComputed = false;              << 104   static G4TouchableHandle nullTouchableHandle;  // Points to (G4VTouchable*) 0
                                                   >> 105   fCurrentTouchableHandle = nullTouchableHandle; 
                                                   >> 106 
                                                   >> 107   fEndGlobalTimeComputed  = false;
115   fCandidateEndGlobalTime = 0;                    108   fCandidateEndGlobalTime = 0;
116 }                                                 109 }
117                                                   110 
118 //....oooOO0OOooo........oooOO0OOooo........oo    111 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
119                                                   112 
120 G4MonopoleTransportation::~G4MonopoleTransport    113 G4MonopoleTransportation::~G4MonopoleTransportation()
121 {                                                 114 {
122   if ((verboseLevel > 0) && (fSumEnergyKilled  << 115   if( (verboseLevel > 0) && (fSumEnergyKilled > 0.0 ) ){ 
123     G4cout << " G4MonopoleTransportation: Stat << 116     G4cout << " G4MonopoleTransportation: Statistics for looping particles " 
124     G4cout << "   Sum of energy of loopers kil << 117            << G4endl;
125     G4cout << "   Max energy of loopers killed << 118     G4cout << "   Sum of energy of loopers killed: " <<  fSumEnergyKilled << G4endl;
126   }                                            << 119     G4cout << "   Max energy of loopers killed: " <<  fMaxEnergyKilled << G4endl;
                                                   >> 120   } 
127 }                                                 121 }
128                                                   122 
129 //....oooOO0OOooo........oooOO0OOooo........oo    123 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
130 //                                                124 //
131 // Responsibilities:                              125 // Responsibilities:
132 //    Find whether the geometry limits the Ste    126 //    Find whether the geometry limits the Step, and to what length
133 //    Calculate the new value of the safety an    127 //    Calculate the new value of the safety and return it.
134 //    Store the final time, position and momen    128 //    Store the final time, position and momentum.
135                                                   129 
136 G4double G4MonopoleTransportation::AlongStepGe << 130 G4double G4MonopoleTransportation::
137   const G4Track& track,                        << 131 AlongStepGetPhysicalInteractionLength( const G4Track&  track,
138   G4double,  //  previousStepSize              << 132                                              G4double, //  previousStepSize
139   G4double currentMinimumStep, G4double& curre << 133                                              G4double  currentMinimumStep,
140 {                                              << 134                                              G4double& currentSafety,
141   fMagSetup->SetStepperAndChordFinder(1);      << 135                                              G4GPILSelection* selection )
                                                   >> 136 {  
                                                   >> 137   fMagSetup->SetStepperAndChordFinder(1); 
142   // change to monopole equation                  138   // change to monopole equation
143                                                   139 
144   G4double geometryStepLength, newSafety;      << 140   G4double geometryStepLength, newSafety ; 
145   fParticleIsLooping = false;                  << 141   fParticleIsLooping = false ;
146                                                   142 
147   // Initial actions moved to  StartTrack()    << 143   // Initial actions moved to  StartTrack()   
148   // --------------------------------------       144   // --------------------------------------
149   // Note: in case another process changes tou    145   // Note: in case another process changes touchable handle
150   //    it will be necessary to add here (for  << 146   //    it will be necessary to add here (for all steps)   
151   // fCurrentTouchableHandle = aTrack->GetTouc    147   // fCurrentTouchableHandle = aTrack->GetTouchableHandle();
152                                                   148 
153   // GPILSelection is set to defaule value of     149   // GPILSelection is set to defaule value of CandidateForSelection
154   // It is a return value                         150   // It is a return value
155   //                                              151   //
156   *selection = CandidateForSelection;          << 152   *selection = CandidateForSelection ;
157                                                   153 
158   // Get initial Energy/Momentum of the track     154   // Get initial Energy/Momentum of the track
159   //                                              155   //
160   const G4DynamicParticle* pParticle = track.G << 156   const G4DynamicParticle* pParticle      = track.GetDynamicParticle() ;
161   G4ThreeVector startMomentumDir = pParticle-> << 157   G4ThreeVector startMomentumDir          = pParticle->GetMomentumDirection() ;
162   G4ThreeVector startPosition = track.GetPosit << 158   G4ThreeVector startPosition             = track.GetPosition() ;
163                                                   159 
164   // G4double   theTime        = track.GetGlob    160   // G4double   theTime        = track.GetGlobalTime() ;
165                                                   161 
166   // The Step Point safety can be limited by o << 162   // The Step Point safety can be limited by other geometries and/or the 
167   // assumptions of any process - it's not alw    163   // assumptions of any process - it's not always the geometrical safety.
168   // We calculate the starting point's isotrop    164   // We calculate the starting point's isotropic safety here.
169   //                                              165   //
170   G4ThreeVector OriginShift = startPosition -  << 166   G4ThreeVector OriginShift = startPosition - fPreviousSftOrigin ;
171   G4double MagSqShift = OriginShift.mag2();    << 167   G4double      MagSqShift  = OriginShift.mag2() ;
172   if (MagSqShift >= sqr(fPreviousSafety)) {    << 168   if( MagSqShift >= sqr(fPreviousSafety) )
173     currentSafety = 0.0;                       << 169   {
                                                   >> 170      currentSafety = 0.0 ;
174   }                                               171   }
175   else {                                       << 172   else
176     currentSafety = fPreviousSafety - std::sqr << 173   {
                                                   >> 174      currentSafety = fPreviousSafety - std::sqrt(MagSqShift) ;
177   }                                               175   }
178                                                   176 
179   // Is the monopole charged ?                    177   // Is the monopole charged ?
180   //                                              178   //
181   G4double particleMagneticCharge = fParticleD << 179   G4double particleMagneticCharge = fParticleDef->MagneticCharge() ; 
182   G4double particleElectricCharge = pParticle-    180   G4double particleElectricCharge = pParticle->GetCharge();
183                                                   181 
184   fGeometryLimitedStep = false;                << 182   fGeometryLimitedStep = false ;
185   // fEndGlobalTimeComputed = false ;             183   // fEndGlobalTimeComputed = false ;
186                                                   184 
187   // There is no need to locate the current vo    185   // There is no need to locate the current volume. It is Done elsewhere:
188   //   On track construction                   << 186   //   On track construction 
189   //   By the tracking, after all AlongStepDoI    187   //   By the tracking, after all AlongStepDoIts, in "Relocation"
190                                                   188 
191   // Check whether the particle have an (EM) f    189   // Check whether the particle have an (EM) field force exerting upon it
192   //                                              190   //
193   G4FieldManager* fieldMgr = 0;                << 191   G4FieldManager* fieldMgr=0;
194   G4bool fieldExertsForce = false;             << 192   G4bool          fieldExertsForce = false ;
195                                                << 193     
196   if ((particleMagneticCharge != 0.0)) {       << 194   if( (particleMagneticCharge != 0.0) )
197     fieldMgr = fFieldPropagator->FindAndSetFie << 195   {      
198     if (fieldMgr != 0) {                       << 196      fieldMgr= fFieldPropagator->FindAndSetFieldManager( track.GetVolume() ); 
199       // Message the field Manager, to configu << 197      if (fieldMgr != 0) {
200       fieldMgr->ConfigureForTrack(&track);     << 198        // Message the field Manager, to configure it for this track
201       // Moved here, in order to allow a trans << 199        fieldMgr->ConfigureForTrack( &track );
202       //   from a zero-field  status (with fie << 200        // Moved here, in order to allow a transition
203       //   to a finite field  status           << 201        //   from a zero-field  status (with fieldMgr->(field)0
204                                                << 202        //   to a finite field  status
205       // If the field manager has no field, th << 203 
206       fieldExertsForce = (fieldMgr->GetDetecto << 204        // If the field manager has no field, there is no field !
207     }                                          << 205        fieldExertsForce = (fieldMgr->GetDetectorField() != 0);
                                                   >> 206      }      
208   }                                               207   }
209                                                   208 
210   // G4cout << " G4Transport:  field exerts fo    209   // G4cout << " G4Transport:  field exerts force= " << fieldExertsForce
211   //          << "  fieldMgr= " << fieldMgr <<    210   //          << "  fieldMgr= " << fieldMgr << G4endl;
212                                                   211 
213   // Choose the calculation of the transportat << 212   // Choose the calculation of the transportation: Field or not 
214   //                                              213   //
215   if (!fieldExertsForce) {                     << 214   if( !fieldExertsForce ) 
216     G4double linearStepLength;                 << 215   {
217     if (fShortStepOptimisation && (currentMini << 216      G4double linearStepLength ;
218       // The Step is guaranteed to be taken    << 217      if( fShortStepOptimisation && (currentMinimumStep <= currentSafety) )
219       //                                       << 218      {
220       geometryStepLength = currentMinimumStep; << 219        // The Step is guaranteed to be taken
221       fGeometryLimitedStep = false;            << 220        //
222     }                                          << 221        geometryStepLength   = currentMinimumStep ;
223     else {                                     << 222        fGeometryLimitedStep = false ;
224       //  Find whether the straight path inter << 223      }
225       //                                       << 224      else
226       linearStepLength = fLinearNavigator->Com << 225      {
227                                                << 226        //  Find whether the straight path intersects a volume
228       // Remember last safety origin & value.  << 227        //
229       //                                       << 228        linearStepLength = fLinearNavigator->ComputeStep( startPosition, 
230       fPreviousSftOrigin = startPosition;      << 229                                                          startMomentumDir,
231       fPreviousSafety = newSafety;             << 230                                                          currentMinimumStep, 
232       // fpSafetyHelper->SetCurrentSafety( new << 231                                                          newSafety) ;
233                                                << 232        // Remember last safety origin & value.
234       // The safety at the initial point has b << 233        //
235       //                                       << 234        fPreviousSftOrigin = startPosition ;
236       currentSafety = newSafety;               << 235        fPreviousSafety    = newSafety ; 
237                                                << 236        // fpSafetyHelper->SetCurrentSafety( newSafety, startPosition);
238       fGeometryLimitedStep = (linearStepLength << 237 
239       if (fGeometryLimitedStep) {              << 238        // The safety at the initial point has been re-calculated:
240         // The geometry limits the Step size ( << 239        //
241         geometryStepLength = linearStepLength; << 240        currentSafety = newSafety ;
242       }                                        << 241           
243       else {                                   << 242        fGeometryLimitedStep= (linearStepLength <= currentMinimumStep); 
244         // The full Step is taken.             << 243        if( fGeometryLimitedStep )
245         geometryStepLength = currentMinimumSte << 244        {
246       }                                        << 245          // The geometry limits the Step size (an intersection was found.)
247     }                                          << 246          geometryStepLength   = linearStepLength ;
248     endpointDistance = geometryStepLength;     << 247        } 
249                                                << 248        else
250     // Calculate final position                << 249        {
251     //                                         << 250          // The full Step is taken.
252     fTransportEndPosition = startPosition + ge << 251          geometryStepLength   = currentMinimumStep ;
253                                                << 252        }
254     // Momentum direction, energy and polarisa << 253      }
255     //                                         << 254      endpointDistance = geometryStepLength ;
256     fTransportEndMomentumDir = startMomentumDi << 255 
257     fTransportEndKineticEnergy = track.GetKine << 256      // Calculate final position
258     fTransportEndSpin = track.GetPolarization( << 257      //
259     fParticleIsLooping = false;                << 258      fTransportEndPosition = startPosition+geometryStepLength*startMomentumDir ;
260     fMomentumChanged = false;                  << 259 
261     fEndGlobalTimeComputed = false;            << 260      // Momentum direction, energy and polarisation are unchanged by transport
262   }                                            << 261      //
263   else  //  A field exerts force               << 262      fTransportEndMomentumDir   = startMomentumDir ; 
264   {                                            << 263      fTransportEndKineticEnergy = track.GetKineticEnergy() ;
265     G4double momentumMagnitude = pParticle->Ge << 264      fTransportEndSpin          = track.GetPolarization();
266     G4ThreeVector EndUnitMomentum;             << 265      fParticleIsLooping         = false ;
267     G4double lengthAlongCurve;                 << 266      fMomentumChanged           = false ; 
268     G4double restMass = fParticleDef->GetPDGMa << 267      fEndGlobalTimeComputed     = false ;
269                                                << 268   }
270     G4ChargeState chargeState(particleElectric << 269   else   //  A field exerts force
271                               fParticleDef->Ge << 270   {
272                               0,  //  Magnetic << 271      G4double       momentumMagnitude = pParticle->GetTotalMomentum() ;
273                               0,  //  Electric << 272      G4ThreeVector  EndUnitMomentum ;
274                               particleMagnetic << 273      G4double       lengthAlongCurve ;
275                                                << 274      G4double       restMass = fParticleDef->GetPDGMass() ;
276     G4EquationOfMotion* equationOfMotion =     << 275 
277       fFieldPropagator->GetChordFinder()->GetI << 276      G4ChargeState chargeState(particleElectricCharge,  // The charge can change (dynamic)
278                                                << 277                                fParticleDef->GetPDGSpin(),
279     equationOfMotion->SetChargeMomentumMass(ch << 278                                0,   //  Magnetic moment:  pParticleDef->GetMagneticMoment(),
280     // SetChargeMomentumMass now passes both t << 279                                0,   //  Electric Dipole moment - not in Particle Definition 
281     // charge in chargeState                   << 280                                particleMagneticCharge );   // in Mev/c 
282                                                << 281 
283     G4ThreeVector spin = track.GetPolarization << 282      G4EquationOfMotion* equationOfMotion = 
284     G4FieldTrack aFieldTrack = G4FieldTrack(st << 283      (fFieldPropagator->GetChordFinder()->GetIntegrationDriver()->GetStepper())
285                                             tr << 284      ->GetEquationOfMotion();
286                                             tr << 285 
287                                             tr << 286      equationOfMotion
288                                             &s << 287        ->SetChargeMomentumMass( chargeState,       //  Was particleMagneticCharge - in Mev/c
289     if (currentMinimumStep > 0) {              << 288                                 momentumMagnitude, //  Was particleElectricCharge 
290       // Do the Transport in the field (non re << 289                                 restMass ) ;  
291       //                                       << 290      // SetChargeMomentumMass now passes both the electric and magnetic charge - in chargeState
292       lengthAlongCurve = fFieldPropagator->Com << 291 
293                                                << 292      G4ThreeVector spin        = track.GetPolarization() ;
294       fGeometryLimitedStep = lengthAlongCurve  << 293      G4FieldTrack  aFieldTrack = G4FieldTrack( startPosition, 
295       if (fGeometryLimitedStep) {              << 294                                                track.GetMomentumDirection(),
296         geometryStepLength = lengthAlongCurve; << 295                                                0.0, 
297       }                                        << 296                                                track.GetKineticEnergy(),
298       else {                                   << 297                                                restMass,
299         geometryStepLength = currentMinimumSte << 298                                                track.GetVelocity(),
300       }                                        << 299                                                track.GetGlobalTime(), // Lab.
301     }                                          << 300                                                track.GetProperTime(), // Part.
302     else {                                     << 301                                                &spin                  ) ;
303       geometryStepLength = lengthAlongCurve =  << 302      if( currentMinimumStep > 0 ) 
304       fGeometryLimitedStep = false;            << 303      {
305     }                                          << 304         // Do the Transport in the field (non recti-linear)
306                                                << 305         //
307     // Remember last safety origin & value.    << 306         lengthAlongCurve = fFieldPropagator->ComputeStep( aFieldTrack,
308     //                                         << 307                                                           currentMinimumStep, 
309     fPreviousSftOrigin = startPosition;        << 308                                                           currentSafety,
310     fPreviousSafety = currentSafety;           << 309                                                           track.GetVolume() ) ;
311     // fpSafetyHelper->SetCurrentSafety( newSa << 310         fGeometryLimitedStep= lengthAlongCurve < currentMinimumStep; 
312                                                << 311         if( fGeometryLimitedStep ) {
313     // Get the End-Position and End-Momentum ( << 312            geometryStepLength   = lengthAlongCurve ;
314     //                                         << 313         } else {
315     fTransportEndPosition = aFieldTrack.GetPos << 314            geometryStepLength   = currentMinimumStep ;
316                                                << 315         }
317     // Momentum:  Magnitude and direction can  << 316      }
318     //                                         << 317      else
319     fMomentumChanged = true;                   << 318      {
320     fTransportEndMomentumDir = aFieldTrack.Get << 319         geometryStepLength   = lengthAlongCurve= 0.0 ;
321                                                << 320         fGeometryLimitedStep = false ;
322     fTransportEndKineticEnergy = aFieldTrack.G << 321      }
323                                                << 322 
324     fCandidateEndGlobalTime = aFieldTrack.GetL << 323      // Remember last safety origin & value.
325     fEndGlobalTimeComputed = true;             << 324      //
326                                                << 325      fPreviousSftOrigin = startPosition ;
327     fTransportEndSpin = aFieldTrack.GetSpin(); << 326      fPreviousSafety    = currentSafety ;         
328     fParticleIsLooping = fFieldPropagator->IsP << 327      // fpSafetyHelper->SetCurrentSafety( newSafety, startPosition);
329     endpointDistance = (fTransportEndPosition  << 328        
                                                   >> 329      // Get the End-Position and End-Momentum (Dir-ection)
                                                   >> 330      //
                                                   >> 331      fTransportEndPosition = aFieldTrack.GetPosition() ;
                                                   >> 332 
                                                   >> 333      // Momentum:  Magnitude and direction can be changed too now ...
                                                   >> 334      //
                                                   >> 335      fMomentumChanged         = true ; 
                                                   >> 336      fTransportEndMomentumDir = aFieldTrack.GetMomentumDir() ;
                                                   >> 337 
                                                   >> 338      fTransportEndKineticEnergy  = aFieldTrack.GetKineticEnergy() ; 
                                                   >> 339 
                                                   >> 340      fCandidateEndGlobalTime   = aFieldTrack.GetLabTimeOfFlight();
                                                   >> 341      fEndGlobalTimeComputed    = true;
                                                   >> 342 
                                                   >> 343      fTransportEndSpin = aFieldTrack.GetSpin();
                                                   >> 344      fParticleIsLooping = fFieldPropagator->IsParticleLooping() ;
                                                   >> 345      endpointDistance   = (fTransportEndPosition - startPosition).mag() ;
330   }                                               346   }
331                                                   347 
332   // If we are asked to go a step length of 0,    348   // If we are asked to go a step length of 0, and we are on a boundary
333   // then a boundary will also limit the step     349   // then a boundary will also limit the step -> we must flag this.
334   //                                              350   //
335   if (currentMinimumStep == 0.0) {             << 351   if( currentMinimumStep == 0.0 ) 
336     if (currentSafety == 0.0) fGeometryLimited << 352   {
                                                   >> 353       if( currentSafety == 0.0 )  fGeometryLimitedStep = true ;
337   }                                               354   }
338                                                   355 
339   // Update the safety starting from the end-p    356   // Update the safety starting from the end-point,
340   // if it will become negative at the end-poi    357   // if it will become negative at the end-point.
341   //                                              358   //
342   if (currentSafety < endpointDistance) {      << 359   if( currentSafety < endpointDistance ) 
343     // if( particleMagneticCharge == 0.0 )     << 360   {
344     // G4cout  << "  Avoiding call to ComputeS << 361       // if( particleMagneticCharge == 0.0 ) 
345                                                << 362       //    G4cout  << "  Avoiding call to ComputeSafety : charge = 0.0 " << G4endl;
346     if (particleMagneticCharge != 0.0) {       << 363  
347       G4double endSafety = fLinearNavigator->C << 364       if( particleMagneticCharge != 0.0 ) {
348       currentSafety = endSafety;               << 365 
349       fPreviousSftOrigin = fTransportEndPositi << 366          G4double endSafety =
350       fPreviousSafety = currentSafety;         << 367                fLinearNavigator->ComputeSafety( fTransportEndPosition) ;
351       fpSafetyHelper->SetCurrentSafety(current << 368          currentSafety      = endSafety ;
352                                                << 369          fPreviousSftOrigin = fTransportEndPosition ;
353       // Because the Stepping Manager assumes  << 370          fPreviousSafety    = currentSafety ; 
354       //  add the StepLength                   << 371          fpSafetyHelper->SetCurrentSafety( currentSafety, fTransportEndPosition);
355       //                                       << 372 
356       currentSafety += endpointDistance;       << 373          // Because the Stepping Manager assumes it is from the start point, 
357                                                << 374          //  add the StepLength
358 #ifdef G4DEBUG_TRANSPORT                       << 375          //
359       G4cout.precision(12);                    << 376          currentSafety     += endpointDistance ;
360       G4cout << "***G4MonopoleTransportation:: << 377 
361       G4cout << "  Called Navigator->ComputeSa << 378 #ifdef G4DEBUG_TRANSPORT 
362              << "    and it returned safety= " << 379          G4cout.precision(12) ;
363       G4cout << "  Adding endpoint distance "  << 380          G4cout << "***G4MonopoleTransportation::AlongStepGPIL ** " << G4endl  ;
364              << "   to obtain pseudo-safety= " << 381          G4cout << "  Called Navigator->ComputeSafety at " << fTransportEndPosition
                                                   >> 382                 << "    and it returned safety= " << endSafety << G4endl ; 
                                                   >> 383          G4cout << "  Adding endpoint distance " << endpointDistance 
                                                   >> 384                 << "   to obtain pseudo-safety= " << currentSafety << G4endl ; 
365 #endif                                            385 #endif
366     }                                          << 386       }
367   }                                            << 387   }            
368                                                   388 
369   fParticleChange.ProposeTrueStepLength(geomet << 389   fParticleChange.ProposeTrueStepLength(geometryStepLength) ;
370                                                   390 
371   fMagSetup->SetStepperAndChordFinder(0);         391   fMagSetup->SetStepperAndChordFinder(0);
372   // change back to usual equation                392   // change back to usual equation
373                                                   393 
374   return geometryStepLength;                   << 394   return geometryStepLength ;
375 }                                                 395 }
376                                                   396 
377 //....oooOO0OOooo........oooOO0OOooo........oo    397 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
378 //                                                398 //
379 //   Initialize ParticleChange  (by setting al    399 //   Initialize ParticleChange  (by setting all its members equal
380 //                               to correspond    400 //                               to corresponding members in G4Track)
381                                                   401 
382 G4VParticleChange* G4MonopoleTransportation::A << 402 G4VParticleChange* G4MonopoleTransportation::AlongStepDoIt( const G4Track& track,
383                                                << 403                                                     const G4Step&  stepData )
384 {                                                 404 {
385   ++noCalls;                                   << 405   static G4int noCalls=0;
                                                   >> 406   static const G4ParticleDefinition* fOpticalPhoton =
                                                   >> 407            G4ParticleTable::GetParticleTable()->FindParticle("opticalphoton");
386                                                   408 
387   if (fGeometryLimitedStep) {                  << 409   noCalls++;
388     stepData.GetPostStepPoint()->SetStepStatus << 
389   }                                            << 
390                                                   410 
391   fParticleChange.Initialize(track);           << 411   fParticleChange.Initialize(track) ;
392                                                   412 
393   //  Code for specific process                << 413   //  Code for specific process 
394   //                                              414   //
395   fParticleChange.ProposePosition(fTransportEn << 415   fParticleChange.ProposePosition(fTransportEndPosition) ;
396   fParticleChange.ProposeMomentumDirection(fTr << 416   fParticleChange.ProposeMomentumDirection(fTransportEndMomentumDir) ;
397   fParticleChange.ProposeEnergy(fTransportEndK << 417   fParticleChange.ProposeEnergy(fTransportEndKineticEnergy) ;
398   fParticleChange.SetMomentumChanged(fMomentum << 418   fParticleChange.SetMomentumChanged(fMomentumChanged) ;
399                                                   419 
400   fParticleChange.ProposePolarization(fTranspo    420   fParticleChange.ProposePolarization(fTransportEndSpin);
401                                                << 421   
402   G4double deltaTime = 0.0;                    << 422   G4double deltaTime = 0.0 ;
403                                                   423 
404   // Calculate  Lab Time of Flight (ONLY if fi    424   // Calculate  Lab Time of Flight (ONLY if field Equations used it!)
                                                   >> 425      // G4double endTime   = fCandidateEndGlobalTime;
                                                   >> 426      // G4double delta_time = endTime - startTime;
405                                                   427 
406   G4double startTime = track.GetGlobalTime();  << 428   G4double startTime = track.GetGlobalTime() ;
407                                                << 429   
408   if (!fEndGlobalTimeComputed) {               << 430   if (!fEndGlobalTimeComputed)
409     // The time was not integrated .. make the << 431   {
410     //                                         << 432      // The time was not integrated .. make the best estimate possible
411     G4double finalVelocity = track.GetVelocity << 433      //
412     G4double initialVelocity = stepData.GetPre << 434      G4double finalVelocity   = track.GetVelocity() ;
413     G4double stepLength = track.GetStepLength( << 435      G4double initialVelocity = stepData.GetPreStepPoint()->GetVelocity() ;
414                                                << 436      G4double stepLength      = track.GetStepLength() ;
415     deltaTime = 0.0;  // in case initialVeloci << 437 
416     if (finalVelocity > 0.0) {                 << 438      deltaTime= 0.0;  // in case initialVelocity = 0 
417       G4double meanInverseVelocity;            << 439      const G4DynamicParticle* fpDynamicParticle = track.GetDynamicParticle();
418       // deltaTime = stepLength/finalVelocity  << 440      if (fpDynamicParticle->GetDefinition()== fOpticalPhoton)
419       meanInverseVelocity = 0.5 * (1.0 / initi << 441      {
420       deltaTime = stepLength * meanInverseVelo << 442         //  A photon is in the medium of the final point
421     }                                          << 443         //  during the step, so it has the final velocity.
422     else if (initialVelocity > 0.0) {          << 444         deltaTime = stepLength/finalVelocity ;
423       deltaTime = stepLength / initialVelocity << 445      }
424     }                                          << 446      else if (finalVelocity > 0.0)
425     fCandidateEndGlobalTime = startTime + delt << 447      {
                                                   >> 448         G4double meanInverseVelocity ;
                                                   >> 449         // deltaTime = stepLength/finalVelocity ;
                                                   >> 450         meanInverseVelocity = 0.5
                                                   >> 451                             * ( 1.0 / initialVelocity + 1.0 / finalVelocity ) ;
                                                   >> 452         deltaTime = stepLength * meanInverseVelocity ;
                                                   >> 453      }
                                                   >> 454      else if( initialVelocity > 0.0 )
                                                   >> 455      {
                                                   >> 456         deltaTime = stepLength/initialVelocity ;
                                                   >> 457      }
                                                   >> 458      fCandidateEndGlobalTime   = startTime + deltaTime ;
426   }                                               459   }
427   else {                                       << 460   else
428     deltaTime = fCandidateEndGlobalTime - star << 461   {
                                                   >> 462      deltaTime = fCandidateEndGlobalTime - startTime ;
429   }                                               463   }
430                                                   464 
431   fParticleChange.ProposeGlobalTime(fCandidate << 465   fParticleChange.ProposeGlobalTime( fCandidateEndGlobalTime ) ;
432                                                   466 
433   // Now Correct by Lorentz factor to get "pro    467   // Now Correct by Lorentz factor to get "proper" deltaTime
                                                   >> 468   
                                                   >> 469   G4double  restMass       = track.GetDynamicParticle()->GetMass() ;
                                                   >> 470   G4double deltaProperTime = deltaTime*( restMass/track.GetTotalEnergy() ) ;
434                                                   471 
435   G4double restMass = track.GetDynamicParticle << 472   fParticleChange.ProposeProperTime(track.GetProperTime() + deltaProperTime) ;
436   G4double deltaProperTime = deltaTime * (rest << 473   //fParticleChange. ProposeTrueStepLength( track.GetStepLength() ) ;
437                                                << 
438   fParticleChange.ProposeProperTime(track.GetP << 
439                                                   474 
440   // If the particle is caught looping or is s    475   // If the particle is caught looping or is stuck (in very difficult
441   // boundaries) in a magnetic field (doing ma << 476   // boundaries) in a magnetic field (doing many steps) 
442   //   THEN this kills it ...                     477   //   THEN this kills it ...
443   //                                              478   //
444   if (fParticleIsLooping) {                    << 479   if ( fParticleIsLooping )
445     G4double endEnergy = fTransportEndKineticE << 480   {
                                                   >> 481       G4double endEnergy= fTransportEndKineticEnergy;
446                                                   482 
447     if ((endEnergy < fThreshold_Important_Ener << 483       if( (endEnergy < fThreshold_Important_Energy) 
448       // Kill the looping particle             << 484           || (fNoLooperTrials >= fThresholdTrials ) ){
449       //                                       << 485         // Kill the looping particle 
450       fParticleChange.ProposeTrackStatus(fStop << 486         //
451                                                << 487         fParticleChange.ProposeTrackStatus( fStopAndKill )  ;
452       // 'Bare' statistics                     << 488 
453       fSumEnergyKilled += endEnergy;           << 489         // 'Bare' statistics
454       if (endEnergy > fMaxEnergyKilled) {      << 490         fSumEnergyKilled += endEnergy; 
455         fMaxEnergyKilled = endEnergy;          << 491         if( endEnergy > fMaxEnergyKilled) { fMaxEnergyKilled= endEnergy; }
456       }                                        << 
457                                                   492 
458 #ifdef G4VERBOSE                                  493 #ifdef G4VERBOSE
459       if ((verboseLevel > 1) || (endEnergy > f << 494         if( (verboseLevel > 1) || 
460         G4cout << " G4MonopoleTransportation i << 495             ( endEnergy > fThreshold_Warning_Energy )  ) { 
461                << "that is looping or stuck "  << 496           G4cout << " G4MonopoleTransportation is killing track that is looping or stuck "
462                << track.GetKineticEnergy() / M << 497                  << G4endl
463         G4cout << "   Number of trials = " <<  << 498                  << "   This track has " << track.GetKineticEnergy() / MeV
464                << "   No of calls to AlongStep << 499                  << " MeV energy." << G4endl;
465       }                                        << 500           G4cout << "   Number of trials = " << fNoLooperTrials 
                                                   >> 501                  << "   No of calls to AlongStepDoIt = " << noCalls 
                                                   >> 502                  << G4endl;
                                                   >> 503         }
466 #endif                                            504 #endif
467       fNoLooperTrials = 0;                     << 505         fNoLooperTrials=0; 
468     }                                          << 
469     else {                                     << 
470       fNoLooperTrials++;                       << 
471 #ifdef G4VERBOSE                               << 
472       if ((verboseLevel > 2)) {                << 
473         G4cout << "   G4MonopoleTransportation << 
474                << "Particle looping - "        << 
475                << "   Number of trials = " <<  << 
476                << G4endl;                      << 
477       }                                           506       }
                                                   >> 507       else{
                                                   >> 508         fNoLooperTrials ++; 
                                                   >> 509 #ifdef G4VERBOSE
                                                   >> 510         if( (verboseLevel > 2) ){
                                                   >> 511           G4cout << "   G4MonopoleTransportation::AlongStepDoIt(): Particle looping -  "
                                                   >> 512                  << "   Number of trials = " << fNoLooperTrials 
                                                   >> 513                  << "   No of calls to  = " << noCalls 
                                                   >> 514                  << G4endl;
                                                   >> 515         }
478 #endif                                            516 #endif
479     }                                          << 517       }
480   }                                            << 518   }else{
481   else {                                       << 519       fNoLooperTrials=0; 
482     fNoLooperTrials = 0;                       << 
483   }                                               520   }
484                                                   521 
485   // Another (sometimes better way) is to use     522   // Another (sometimes better way) is to use a user-limit maximum Step size
486   // to alleviate this problem ..              << 523   // to alleviate this problem .. 
487                                                   524 
488   // Introduce smooth curved trajectories to p    525   // Introduce smooth curved trajectories to particle-change
489   //                                              526   //
490   fParticleChange.SetPointerToVectorOfAuxiliar << 527   fParticleChange.SetPointerToVectorOfAuxiliaryPoints
491     fFieldPropagator->GimmeTrajectoryVectorAnd << 528     (fFieldPropagator->GimmeTrajectoryVectorAndForgetIt() );
492                                                   529 
493   return &fParticleChange;                     << 530   return &fParticleChange ;
494 }                                                 531 }
495                                                   532 
496 //....oooOO0OOooo........oooOO0OOooo........oo    533 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
497 //                                                534 //
498 //  This ensures that the PostStep action is a    535 //  This ensures that the PostStep action is always called,
499 //  so that it can do the relocation if it is     536 //  so that it can do the relocation if it is needed.
500 //                                             << 537 // 
501                                                   538 
502 G4double                                       << 539 G4double G4MonopoleTransportation::
503 G4MonopoleTransportation::PostStepGetPhysicalI << 540 PostStepGetPhysicalInteractionLength( const G4Track&,
504                                                << 541                                             G4double, // previousStepSize
505                                                << 542                                             G4ForceCondition* pForceCond )
506 {                                              << 543 {  
507   *pForceCond = Forced;                        << 544   *pForceCond = Forced ; 
508   return DBL_MAX;  // was kInfinity ; but conv << 545   return DBL_MAX ;  // was kInfinity ; but convention now is DBL_MAX
509 }                                                 546 }
510                                                   547 
511 //....oooOO0OOooo........oooOO0OOooo........oo    548 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
512                                                   549 
513 G4VParticleChange* G4MonopoleTransportation::P << 550 G4VParticleChange* G4MonopoleTransportation::PostStepDoIt( const G4Track& track,
                                                   >> 551                                                    const G4Step& )
514 {                                                 552 {
515   G4TouchableHandle retCurrentTouchable;  // T << 553   G4TouchableHandle retCurrentTouchable ;   // The one to return
516                                                   554 
517   // Initialize ParticleChange  (by setting al    555   // Initialize ParticleChange  (by setting all its members equal
518   //                             to correspond    556   //                             to corresponding members in G4Track)
519   // fParticleChange.Initialize(track) ;  // T    557   // fParticleChange.Initialize(track) ;  // To initialise TouchableChange
520                                                   558 
521   fParticleChange.ProposeTrackStatus(track.Get << 559   fParticleChange.ProposeTrackStatus(track.GetTrackStatus()) ;
522                                                   560 
523   // If the Step was determined by the volume     561   // If the Step was determined by the volume boundary,
524   // logically relocate the particle              562   // logically relocate the particle
525                                                << 563   
526   if (fGeometryLimitedStep) {                  << 564   if(fGeometryLimitedStep)
527     // fCurrentTouchable will now become the p << 565   {  
                                                   >> 566     // fCurrentTouchable will now become the previous touchable, 
528     // and what was the previous will be freed    567     // and what was the previous will be freed.
529     // (Needed because the preStepPoint can po    568     // (Needed because the preStepPoint can point to the previous touchable)
530                                                   569 
531     fLinearNavigator->SetGeometricallyLimitedS << 570     fLinearNavigator->SetGeometricallyLimitedStep() ;
532     fLinearNavigator->LocateGlobalPointAndUpda << 571     fLinearNavigator->
533       track.GetPosition(), track.GetMomentumDi << 572     LocateGlobalPointAndUpdateTouchableHandle( track.GetPosition(),
534     // Check whether the particle is out of th << 573                                                track.GetMomentumDirection(),
                                                   >> 574                                                fCurrentTouchableHandle,
                                                   >> 575                                                true                      ) ;
                                                   >> 576     // Check whether the particle is out of the world volume 
535     // If so it has exited and must be killed.    577     // If so it has exited and must be killed.
536     //                                            578     //
537     if (fCurrentTouchableHandle->GetVolume() = << 579     if( fCurrentTouchableHandle->GetVolume() == 0 )
538       fParticleChange.ProposeTrackStatus(fStop << 580     {
                                                   >> 581        fParticleChange.ProposeTrackStatus( fStopAndKill ) ;
539     }                                             582     }
540     retCurrentTouchable = fCurrentTouchableHan << 583     retCurrentTouchable = fCurrentTouchableHandle ;
541     fParticleChange.SetTouchableHandle(fCurren << 584     fParticleChange.SetTouchableHandle( fCurrentTouchableHandle ) ;
542   }                                               585   }
543   else  // fGeometryLimitedStep  is false      << 586   else                 // fGeometryLimitedStep  is false
544   {                                            << 587   {                    
545     // This serves only to move the Navigator'    588     // This serves only to move the Navigator's location
546     //                                            589     //
547     fLinearNavigator->LocateGlobalPointWithinV << 590     fLinearNavigator->LocateGlobalPointWithinVolume( track.GetPosition() ) ;
548                                                   591 
549     // The value of the track's current Toucha << 592     // The value of the track's current Touchable is retained. 
550     // (and it must be correct because we must    593     // (and it must be correct because we must use it below to
551     // overwrite the (unset) one in particle c    594     // overwrite the (unset) one in particle change)
552     //  It must be fCurrentTouchable too ??       595     //  It must be fCurrentTouchable too ??
553     //                                            596     //
554     fParticleChange.SetTouchableHandle(track.G << 597     fParticleChange.SetTouchableHandle( track.GetTouchableHandle() ) ;
555     retCurrentTouchable = track.GetTouchableHa << 598     retCurrentTouchable = track.GetTouchableHandle() ;
556   }  // endif ( fGeometryLimitedStep )         << 599   }         // endif ( fGeometryLimitedStep ) 
557                                                << 600 
558   const G4VPhysicalVolume* pNewVol = retCurren << 601   const G4VPhysicalVolume* pNewVol = retCurrentTouchable->GetVolume() ;
559   const G4Material* pNewMaterial = 0;          << 602   const G4Material* pNewMaterial   = 0 ;
560   G4VSensitiveDetector* pNewSensitiveDetector  << 603   const G4VSensitiveDetector* pNewSensitiveDetector   = 0 ;
561   if (pNewVol != 0) {                          << 604                                                                                        
562     pNewMaterial = pNewVol->GetLogicalVolume() << 605   if( pNewVol != 0 )
563     pNewSensitiveDetector = pNewVol->GetLogica << 606   {
                                                   >> 607     pNewMaterial= pNewVol->GetLogicalVolume()->GetMaterial();
                                                   >> 608     pNewSensitiveDetector= pNewVol->GetLogicalVolume()->GetSensitiveDetector();
564   }                                               609   }
565                                                   610 
566   // ( <const_cast> pNewMaterial ) ;              611   // ( <const_cast> pNewMaterial ) ;
                                                   >> 612   // ( <const_cast> pNewSensitiveDetector) ;
567                                                   613 
568   fParticleChange.SetMaterialInTouchable((G4Ma << 614   fParticleChange.SetMaterialInTouchable( 
569   fParticleChange.SetSensitiveDetectorInToucha << 615                      (G4Material *) pNewMaterial ) ;
                                                   >> 616   fParticleChange.SetSensitiveDetectorInTouchable( 
                                                   >> 617                      (G4VSensitiveDetector *) pNewSensitiveDetector ) ;
570                                                   618 
571   const G4MaterialCutsCouple* pNewMaterialCuts    619   const G4MaterialCutsCouple* pNewMaterialCutsCouple = 0;
572   if (pNewVol != 0) {                          << 620   if( pNewVol != 0 )
573     pNewMaterialCutsCouple = pNewVol->GetLogic << 621   {
                                                   >> 622     pNewMaterialCutsCouple=pNewVol->GetLogicalVolume()->GetMaterialCutsCouple();
574   }                                               623   }
575                                                   624 
576   if (pNewVol != 0 && pNewMaterialCutsCouple ! << 625   if( pNewVol!=0 && pNewMaterialCutsCouple!=0 && 
577       && pNewMaterialCutsCouple->GetMaterial() << 626       pNewMaterialCutsCouple->GetMaterial()!=pNewMaterial )
578   {                                               627   {
579     // for parametrized volume                    628     // for parametrized volume
580     //                                            629     //
581     pNewMaterialCutsCouple = G4ProductionCutsT << 630     pNewMaterialCutsCouple =
582       pNewMaterial, pNewMaterialCutsCouple->Ge << 631       G4ProductionCutsTable::GetProductionCutsTable()
                                                   >> 632                              ->GetMaterialCutsCouple(pNewMaterial,
                                                   >> 633                                pNewMaterialCutsCouple->GetProductionCuts());
583   }                                               634   }
584   fParticleChange.SetMaterialCutsCoupleInTouch << 635   fParticleChange.SetMaterialCutsCoupleInTouchable( pNewMaterialCutsCouple );
585                                                   636 
586   // temporarily until Get/Set Material of Par << 637   // temporarily until Get/Set Material of ParticleChange, 
587   // and StepPoint can be made const.          << 638   // and StepPoint can be made const. 
588   // Set the touchable in ParticleChange          639   // Set the touchable in ParticleChange
589   // this must always be done because the part    640   // this must always be done because the particle change always
590   // uses this value to overwrite the current     641   // uses this value to overwrite the current touchable pointer.
591   //                                              642   //
592   fParticleChange.SetTouchableHandle(retCurren << 643   fParticleChange.SetTouchableHandle(retCurrentTouchable) ;
593                                                << 644   
594   return &fParticleChange;                     << 645   return &fParticleChange ;
595 }                                                 646 }
596                                                   647 
597 //....oooOO0OOooo........oooOO0OOooo........oo    648 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
598                                                   649 
599 // New method takes over the responsibility to << 650 // New method takes over the responsibility to reset the state 
600 // of G4MonopoleTransportation object at the s << 651 // of G4MonopoleTransportation object at the start of a new track 
601 // or the resumption of a suspended track.     << 652 // or the resumption of a suspended track. 
602                                                   653 
603 void G4MonopoleTransportation::StartTracking(G << 654 void 
                                                   >> 655 G4MonopoleTransportation::StartTracking(G4Track* aTrack)
604 {                                                 656 {
605   G4VProcess::StartTracking(aTrack);              657   G4VProcess::StartTracking(aTrack);
606                                                   658 
607   // The actions here are those that were take << 659 // The actions here are those that were taken in AlongStepGPIL
608   //   when track.GetCurrentStepNumber()==1    << 660 //   when track.GetCurrentStepNumber()==1
609                                                   661 
610   // reset safety value and center                662   // reset safety value and center
611   //                                              663   //
612   fPreviousSafety = 0.0;                       << 664   fPreviousSafety    = 0.0 ; 
613   fPreviousSftOrigin = G4ThreeVector(0., 0., 0 << 665   fPreviousSftOrigin = G4ThreeVector(0.,0.,0.) ;
614                                                << 666   
615   // reset looping counter -- for motion in fi    667   // reset looping counter -- for motion in field
616   fNoLooperTrials = 0;                         << 668   fNoLooperTrials= 0; 
617   // Must clear this state .. else it depends     669   // Must clear this state .. else it depends on last track's value
618   //  --> a better solution would set this fro << 670   //  --> a better solution would set this from state of suspended track TODO ? 
619   // Was if( aTrack->GetCurrentStepNumber()==1    671   // Was if( aTrack->GetCurrentStepNumber()==1 ) { .. }
620                                                   672 
621   // ChordFinder reset internal state             673   // ChordFinder reset internal state
622   //                                              674   //
623   if (DoesGlobalFieldExist()) {                << 675   if( DoesGlobalFieldExist() ) {
624     fFieldPropagator->ClearPropagatorState();  << 676      fFieldPropagator->ClearPropagatorState();   
625     // Resets all state of field propagator cl << 677        // Resets all state of field propagator class (ONLY)
626     // including safety values                 << 678        //  including safety values (in case of overlaps and to wipe for first track).
627     // in case of overlaps and to wipe for fir << 
628                                                   679 
629     // G4ChordFinder* chordF= fFieldPropagator << 680      // G4ChordFinder* chordF= fFieldPropagator->GetChordFinder();
630     // if( chordF ) chordF->ResetStepEstimate( << 681      // if( chordF ) chordF->ResetStepEstimate();
631   }                                               682   }
632                                                   683 
633   // Make sure to clear the chord finders of a    684   // Make sure to clear the chord finders of all fields (ie managers)
634   G4FieldManagerStore* fieldMgrStore = G4Field << 685   static G4FieldManagerStore* fieldMgrStore= G4FieldManagerStore::GetInstance();
635   fieldMgrStore->ClearAllChordFindersState();  << 686   fieldMgrStore->ClearAllChordFindersState(); 
636                                                   687 
637   // Update the current touchable handle  (fro    688   // Update the current touchable handle  (from the track's)
638   //                                              689   //
639   fCurrentTouchableHandle = aTrack->GetTouchab    690   fCurrentTouchableHandle = aTrack->GetTouchableHandle();
640 }                                                 691 }
641                                                   692 
642 //....oooOO0OOooo........oooOO0OOooo........oo    693 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
                                                   >> 694 
643                                                   695