Back to home page

EIC code displayed by LXR

 
 

    


File indexing completed on 2025-02-23 09:21:13

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 /// \file exoticphysics/monopole/src/G4MonopoleTransportation.cc
0027 /// \brief Implementation of the G4MonopoleTransportation class
0028 //
0029 //
0030 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0031 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0032 //
0033 // This class is a process responsible for the transportation of
0034 // magnetic monopoles, ie the geometrical propagation that encounters the
0035 // geometrical sub-volumes of the detectors.
0036 //
0037 // For monopoles, uses a different equation of motion and ignores energy
0038 // conservation.
0039 //
0040 
0041 // =======================================================================
0042 // Created:  3 May 2010, J. Apostolakis, B. Bozsogi
0043 // =======================================================================
0044 
0045 #include "G4MonopoleTransportation.hh"
0046 
0047 #include "DetectorConstruction.hh"
0048 
0049 #include "G4ChordFinder.hh"
0050 #include "G4FieldManagerStore.hh"
0051 #include "G4Monopole.hh"
0052 #include "G4ParticleTable.hh"
0053 #include "G4ProductionCutsTable.hh"
0054 #include "G4RunManager.hh"
0055 #include "G4SafetyHelper.hh"
0056 #include "G4SystemOfUnits.hh"
0057 #include "G4TransportationProcessType.hh"
0058 
0059 class G4VSensitiveDetector;
0060 
0061 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0062 
0063 G4MonopoleTransportation::G4MonopoleTransportation(const G4Monopole* mpl, G4int verb)
0064   : G4VProcess(G4String("MonopoleTransportation"), fTransportation),
0065     fParticleDef(mpl),
0066     fMagSetup(0),
0067     fLinearNavigator(0),
0068     fFieldPropagator(0),
0069     fParticleIsLooping(false),
0070     fPreviousSftOrigin(0., 0., 0.),
0071     fPreviousSafety(0.0),
0072     fThreshold_Warning_Energy(100 * MeV),
0073     fThreshold_Important_Energy(250 * MeV),
0074     fThresholdTrials(10),
0075     // fUnimportant_Energy( 1 * MeV ),
0076     fNoLooperTrials(0),
0077     fSumEnergyKilled(0.0),
0078     fMaxEnergyKilled(0.0),
0079     fShortStepOptimisation(false),  // Old default: true (=fast short steps)
0080     fpSafetyHelper(0),
0081     noCalls(0)
0082 {
0083   verboseLevel = verb;
0084 
0085   // set Process Sub Type
0086   SetProcessSubType(TRANSPORTATION);
0087 
0088   // Do not finalize the G4MonopoleTransportation class
0089   if (G4Threading::IsMasterThread() && G4Threading::IsMultithreadedApplication()) {
0090     return;
0091   }
0092 
0093   const DetectorConstruction* detector = static_cast<const DetectorConstruction*>(
0094     G4RunManager::GetRunManager()->GetUserDetectorConstruction());
0095 
0096   fMagSetup = detector->GetMonopoleFieldSetup();
0097 
0098   G4TransportationManager* transportMgr = G4TransportationManager::GetTransportationManager();
0099 
0100   fLinearNavigator = transportMgr->GetNavigatorForTracking();
0101 
0102   fFieldPropagator = transportMgr->GetPropagatorInField();
0103   fpSafetyHelper = transportMgr->GetSafetyHelper();
0104 
0105   // New
0106 
0107   // Cannot determine whether a field exists here,
0108   //  because it would only work if the field manager has informed
0109   //  about the detector's field before this transportation process
0110   //  is constructed.
0111   // Instead later the method DoesGlobalFieldExist() is called
0112   fCurrentTouchableHandle = nullptr;
0113 
0114   fEndGlobalTimeComputed = false;
0115   fCandidateEndGlobalTime = 0;
0116 }
0117 
0118 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0119 
0120 G4MonopoleTransportation::~G4MonopoleTransportation()
0121 {
0122   if ((verboseLevel > 0) && (fSumEnergyKilled > 0.0)) {
0123     G4cout << " G4MonopoleTransportation: Statistics for looping particles " << G4endl;
0124     G4cout << "   Sum of energy of loopers killed: " << fSumEnergyKilled << G4endl;
0125     G4cout << "   Max energy of loopers killed: " << fMaxEnergyKilled << G4endl;
0126   }
0127 }
0128 
0129 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0130 //
0131 // Responsibilities:
0132 //    Find whether the geometry limits the Step, and to what length
0133 //    Calculate the new value of the safety and return it.
0134 //    Store the final time, position and momentum.
0135 
0136 G4double G4MonopoleTransportation::AlongStepGetPhysicalInteractionLength(
0137   const G4Track& track,
0138   G4double,  //  previousStepSize
0139   G4double currentMinimumStep, G4double& currentSafety, G4GPILSelection* selection)
0140 {
0141   fMagSetup->SetStepperAndChordFinder(1);
0142   // change to monopole equation
0143 
0144   G4double geometryStepLength, newSafety;
0145   fParticleIsLooping = false;
0146 
0147   // Initial actions moved to  StartTrack()
0148   // --------------------------------------
0149   // Note: in case another process changes touchable handle
0150   //    it will be necessary to add here (for all steps)
0151   // fCurrentTouchableHandle = aTrack->GetTouchableHandle();
0152 
0153   // GPILSelection is set to defaule value of CandidateForSelection
0154   // It is a return value
0155   //
0156   *selection = CandidateForSelection;
0157 
0158   // Get initial Energy/Momentum of the track
0159   //
0160   const G4DynamicParticle* pParticle = track.GetDynamicParticle();
0161   G4ThreeVector startMomentumDir = pParticle->GetMomentumDirection();
0162   G4ThreeVector startPosition = track.GetPosition();
0163 
0164   // G4double   theTime        = track.GetGlobalTime() ;
0165 
0166   // The Step Point safety can be limited by other geometries and/or the
0167   // assumptions of any process - it's not always the geometrical safety.
0168   // We calculate the starting point's isotropic safety here.
0169   //
0170   G4ThreeVector OriginShift = startPosition - fPreviousSftOrigin;
0171   G4double MagSqShift = OriginShift.mag2();
0172   if (MagSqShift >= sqr(fPreviousSafety)) {
0173     currentSafety = 0.0;
0174   }
0175   else {
0176     currentSafety = fPreviousSafety - std::sqrt(MagSqShift);
0177   }
0178 
0179   // Is the monopole charged ?
0180   //
0181   G4double particleMagneticCharge = fParticleDef->MagneticCharge();
0182   G4double particleElectricCharge = pParticle->GetCharge();
0183 
0184   fGeometryLimitedStep = false;
0185   // fEndGlobalTimeComputed = false ;
0186 
0187   // There is no need to locate the current volume. It is Done elsewhere:
0188   //   On track construction
0189   //   By the tracking, after all AlongStepDoIts, in "Relocation"
0190 
0191   // Check whether the particle have an (EM) field force exerting upon it
0192   //
0193   G4FieldManager* fieldMgr = 0;
0194   G4bool fieldExertsForce = false;
0195 
0196   if ((particleMagneticCharge != 0.0)) {
0197     fieldMgr = fFieldPropagator->FindAndSetFieldManager(track.GetVolume());
0198     if (fieldMgr != 0) {
0199       // Message the field Manager, to configure it for this track
0200       fieldMgr->ConfigureForTrack(&track);
0201       // Moved here, in order to allow a transition
0202       //   from a zero-field  status (with fieldMgr->(field)0
0203       //   to a finite field  status
0204 
0205       // If the field manager has no field, there is no field !
0206       fieldExertsForce = (fieldMgr->GetDetectorField() != 0);
0207     }
0208   }
0209 
0210   // G4cout << " G4Transport:  field exerts force= " << fieldExertsForce
0211   //          << "  fieldMgr= " << fieldMgr << G4endl;
0212 
0213   // Choose the calculation of the transportation: Field or not
0214   //
0215   if (!fieldExertsForce) {
0216     G4double linearStepLength;
0217     if (fShortStepOptimisation && (currentMinimumStep <= currentSafety)) {
0218       // The Step is guaranteed to be taken
0219       //
0220       geometryStepLength = currentMinimumStep;
0221       fGeometryLimitedStep = false;
0222     }
0223     else {
0224       //  Find whether the straight path intersects a volume
0225       //
0226       linearStepLength = fLinearNavigator->ComputeStep(startPosition, startMomentumDir,
0227                                                        currentMinimumStep, newSafety);
0228       // Remember last safety origin & value.
0229       //
0230       fPreviousSftOrigin = startPosition;
0231       fPreviousSafety = newSafety;
0232       // fpSafetyHelper->SetCurrentSafety( newSafety, startPosition);
0233 
0234       // The safety at the initial point has been re-calculated:
0235       //
0236       currentSafety = newSafety;
0237 
0238       fGeometryLimitedStep = (linearStepLength <= currentMinimumStep);
0239       if (fGeometryLimitedStep) {
0240         // The geometry limits the Step size (an intersection was found.)
0241         geometryStepLength = linearStepLength;
0242       }
0243       else {
0244         // The full Step is taken.
0245         geometryStepLength = currentMinimumStep;
0246       }
0247     }
0248     endpointDistance = geometryStepLength;
0249 
0250     // Calculate final position
0251     //
0252     fTransportEndPosition = startPosition + geometryStepLength * startMomentumDir;
0253 
0254     // Momentum direction, energy and polarisation are unchanged by transport
0255     //
0256     fTransportEndMomentumDir = startMomentumDir;
0257     fTransportEndKineticEnergy = track.GetKineticEnergy();
0258     fTransportEndSpin = track.GetPolarization();
0259     fParticleIsLooping = false;
0260     fMomentumChanged = false;
0261     fEndGlobalTimeComputed = false;
0262   }
0263   else  //  A field exerts force
0264   {
0265     G4double momentumMagnitude = pParticle->GetTotalMomentum();
0266     G4ThreeVector EndUnitMomentum;
0267     G4double lengthAlongCurve;
0268     G4double restMass = fParticleDef->GetPDGMass();
0269 
0270     G4ChargeState chargeState(particleElectricCharge,  // The charge can change
0271                               fParticleDef->GetPDGSpin(),
0272                               0,  //  Magnetic moment:
0273                               0,  //  Electric Dipole moment
0274                               particleMagneticCharge);  // in Mev/c
0275 
0276     G4EquationOfMotion* equationOfMotion =
0277       fFieldPropagator->GetChordFinder()->GetIntegrationDriver()->GetEquationOfMotion();
0278 
0279     equationOfMotion->SetChargeMomentumMass(chargeState, momentumMagnitude, restMass);
0280     // SetChargeMomentumMass now passes both the electric and magnetic
0281     // charge in chargeState
0282 
0283     G4ThreeVector spin = track.GetPolarization();
0284     G4FieldTrack aFieldTrack = G4FieldTrack(startPosition, track.GetMomentumDirection(), 0.0,
0285                                             track.GetKineticEnergy(), restMass, track.GetVelocity(),
0286                                             track.GetGlobalTime(),  // Lab.
0287                                             track.GetProperTime(),  // Part.
0288                                             &spin);
0289     if (currentMinimumStep > 0) {
0290       // Do the Transport in the field (non recti-linear)
0291       //
0292       lengthAlongCurve = fFieldPropagator->ComputeStep(aFieldTrack, currentMinimumStep,
0293                                                        currentSafety, track.GetVolume());
0294       fGeometryLimitedStep = lengthAlongCurve < currentMinimumStep;
0295       if (fGeometryLimitedStep) {
0296         geometryStepLength = lengthAlongCurve;
0297       }
0298       else {
0299         geometryStepLength = currentMinimumStep;
0300       }
0301     }
0302     else {
0303       geometryStepLength = lengthAlongCurve = 0.0;
0304       fGeometryLimitedStep = false;
0305     }
0306 
0307     // Remember last safety origin & value.
0308     //
0309     fPreviousSftOrigin = startPosition;
0310     fPreviousSafety = currentSafety;
0311     // fpSafetyHelper->SetCurrentSafety( newSafety, startPosition);
0312 
0313     // Get the End-Position and End-Momentum (Dir-ection)
0314     //
0315     fTransportEndPosition = aFieldTrack.GetPosition();
0316 
0317     // Momentum:  Magnitude and direction can be changed too now ...
0318     //
0319     fMomentumChanged = true;
0320     fTransportEndMomentumDir = aFieldTrack.GetMomentumDir();
0321 
0322     fTransportEndKineticEnergy = aFieldTrack.GetKineticEnergy();
0323 
0324     fCandidateEndGlobalTime = aFieldTrack.GetLabTimeOfFlight();
0325     fEndGlobalTimeComputed = true;
0326 
0327     fTransportEndSpin = aFieldTrack.GetSpin();
0328     fParticleIsLooping = fFieldPropagator->IsParticleLooping();
0329     endpointDistance = (fTransportEndPosition - startPosition).mag();
0330   }
0331 
0332   // If we are asked to go a step length of 0, and we are on a boundary
0333   // then a boundary will also limit the step -> we must flag this.
0334   //
0335   if (currentMinimumStep == 0.0) {
0336     if (currentSafety == 0.0) fGeometryLimitedStep = true;
0337   }
0338 
0339   // Update the safety starting from the end-point,
0340   // if it will become negative at the end-point.
0341   //
0342   if (currentSafety < endpointDistance) {
0343     // if( particleMagneticCharge == 0.0 )
0344     // G4cout  << "  Avoiding call to ComputeSafety : charge = 0.0 " << G4endl;
0345 
0346     if (particleMagneticCharge != 0.0) {
0347       G4double endSafety = fLinearNavigator->ComputeSafety(fTransportEndPosition);
0348       currentSafety = endSafety;
0349       fPreviousSftOrigin = fTransportEndPosition;
0350       fPreviousSafety = currentSafety;
0351       fpSafetyHelper->SetCurrentSafety(currentSafety, fTransportEndPosition);
0352 
0353       // Because the Stepping Manager assumes it is from the start point,
0354       //  add the StepLength
0355       //
0356       currentSafety += endpointDistance;
0357 
0358 #ifdef G4DEBUG_TRANSPORT
0359       G4cout.precision(12);
0360       G4cout << "***G4MonopoleTransportation::AlongStepGPIL ** " << G4endl;
0361       G4cout << "  Called Navigator->ComputeSafety at " << fTransportEndPosition
0362              << "    and it returned safety= " << endSafety << G4endl;
0363       G4cout << "  Adding endpoint distance " << endpointDistance
0364              << "   to obtain pseudo-safety= " << currentSafety << G4endl;
0365 #endif
0366     }
0367   }
0368 
0369   fParticleChange.ProposeTrueStepLength(geometryStepLength);
0370 
0371   fMagSetup->SetStepperAndChordFinder(0);
0372   // change back to usual equation
0373 
0374   return geometryStepLength;
0375 }
0376 
0377 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0378 //
0379 //   Initialize ParticleChange  (by setting all its members equal
0380 //                               to corresponding members in G4Track)
0381 
0382 G4VParticleChange* G4MonopoleTransportation::AlongStepDoIt(const G4Track& track,
0383                                                            const G4Step& stepData)
0384 {
0385   ++noCalls;
0386 
0387   if (fGeometryLimitedStep) {
0388     stepData.GetPostStepPoint()->SetStepStatus(fGeomBoundary);
0389   }
0390 
0391   fParticleChange.Initialize(track);
0392 
0393   //  Code for specific process
0394   //
0395   fParticleChange.ProposePosition(fTransportEndPosition);
0396   fParticleChange.ProposeMomentumDirection(fTransportEndMomentumDir);
0397   fParticleChange.ProposeEnergy(fTransportEndKineticEnergy);
0398   fParticleChange.SetMomentumChanged(fMomentumChanged);
0399 
0400   fParticleChange.ProposePolarization(fTransportEndSpin);
0401 
0402   G4double deltaTime = 0.0;
0403 
0404   // Calculate  Lab Time of Flight (ONLY if field Equations used it!)
0405 
0406   G4double startTime = track.GetGlobalTime();
0407 
0408   if (!fEndGlobalTimeComputed) {
0409     // The time was not integrated .. make the best estimate possible
0410     //
0411     G4double finalVelocity = track.GetVelocity();
0412     G4double initialVelocity = stepData.GetPreStepPoint()->GetVelocity();
0413     G4double stepLength = track.GetStepLength();
0414 
0415     deltaTime = 0.0;  // in case initialVelocity = 0
0416     if (finalVelocity > 0.0) {
0417       G4double meanInverseVelocity;
0418       // deltaTime = stepLength/finalVelocity ;
0419       meanInverseVelocity = 0.5 * (1.0 / initialVelocity + 1.0 / finalVelocity);
0420       deltaTime = stepLength * meanInverseVelocity;
0421     }
0422     else if (initialVelocity > 0.0) {
0423       deltaTime = stepLength / initialVelocity;
0424     }
0425     fCandidateEndGlobalTime = startTime + deltaTime;
0426   }
0427   else {
0428     deltaTime = fCandidateEndGlobalTime - startTime;
0429   }
0430 
0431   fParticleChange.ProposeGlobalTime(fCandidateEndGlobalTime);
0432 
0433   // Now Correct by Lorentz factor to get "proper" deltaTime
0434 
0435   G4double restMass = track.GetDynamicParticle()->GetMass();
0436   G4double deltaProperTime = deltaTime * (restMass / track.GetTotalEnergy());
0437 
0438   fParticleChange.ProposeProperTime(track.GetProperTime() + deltaProperTime);
0439 
0440   // If the particle is caught looping or is stuck (in very difficult
0441   // boundaries) in a magnetic field (doing many steps)
0442   //   THEN this kills it ...
0443   //
0444   if (fParticleIsLooping) {
0445     G4double endEnergy = fTransportEndKineticEnergy;
0446 
0447     if ((endEnergy < fThreshold_Important_Energy) || (fNoLooperTrials >= fThresholdTrials)) {
0448       // Kill the looping particle
0449       //
0450       fParticleChange.ProposeTrackStatus(fStopAndKill);
0451 
0452       // 'Bare' statistics
0453       fSumEnergyKilled += endEnergy;
0454       if (endEnergy > fMaxEnergyKilled) {
0455         fMaxEnergyKilled = endEnergy;
0456       }
0457 
0458 #ifdef G4VERBOSE
0459       if ((verboseLevel > 1) || (endEnergy > fThreshold_Warning_Energy)) {
0460         G4cout << " G4MonopoleTransportation is killing track "
0461                << "that is looping or stuck " << G4endl << "   This track has "
0462                << track.GetKineticEnergy() / MeV << " MeV energy." << G4endl;
0463         G4cout << "   Number of trials = " << fNoLooperTrials
0464                << "   No of calls to AlongStepDoIt = " << noCalls << G4endl;
0465       }
0466 #endif
0467       fNoLooperTrials = 0;
0468     }
0469     else {
0470       fNoLooperTrials++;
0471 #ifdef G4VERBOSE
0472       if ((verboseLevel > 2)) {
0473         G4cout << "   G4MonopoleTransportation::AlongStepDoIt(): "
0474                << "Particle looping - "
0475                << "   Number of trials = " << fNoLooperTrials << "   No of calls to  = " << noCalls
0476                << G4endl;
0477       }
0478 #endif
0479     }
0480   }
0481   else {
0482     fNoLooperTrials = 0;
0483   }
0484 
0485   // Another (sometimes better way) is to use a user-limit maximum Step size
0486   // to alleviate this problem ..
0487 
0488   // Introduce smooth curved trajectories to particle-change
0489   //
0490   fParticleChange.SetPointerToVectorOfAuxiliaryPoints(
0491     fFieldPropagator->GimmeTrajectoryVectorAndForgetIt());
0492 
0493   return &fParticleChange;
0494 }
0495 
0496 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0497 //
0498 //  This ensures that the PostStep action is always called,
0499 //  so that it can do the relocation if it is needed.
0500 //
0501 
0502 G4double
0503 G4MonopoleTransportation::PostStepGetPhysicalInteractionLength(const G4Track&,
0504                                                                G4double,  // previousStepSize
0505                                                                G4ForceCondition* pForceCond)
0506 {
0507   *pForceCond = Forced;
0508   return DBL_MAX;  // was kInfinity ; but convention now is DBL_MAX
0509 }
0510 
0511 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0512 
0513 G4VParticleChange* G4MonopoleTransportation::PostStepDoIt(const G4Track& track, const G4Step&)
0514 {
0515   G4TouchableHandle retCurrentTouchable;  // The one to return
0516 
0517   // Initialize ParticleChange  (by setting all its members equal
0518   //                             to corresponding members in G4Track)
0519   // fParticleChange.Initialize(track) ;  // To initialise TouchableChange
0520 
0521   fParticleChange.ProposeTrackStatus(track.GetTrackStatus());
0522 
0523   // If the Step was determined by the volume boundary,
0524   // logically relocate the particle
0525 
0526   if (fGeometryLimitedStep) {
0527     // fCurrentTouchable will now become the previous touchable,
0528     // and what was the previous will be freed.
0529     // (Needed because the preStepPoint can point to the previous touchable)
0530 
0531     fLinearNavigator->SetGeometricallyLimitedStep();
0532     fLinearNavigator->LocateGlobalPointAndUpdateTouchableHandle(
0533       track.GetPosition(), track.GetMomentumDirection(), fCurrentTouchableHandle, true);
0534     // Check whether the particle is out of the world volume
0535     // If so it has exited and must be killed.
0536     //
0537     if (fCurrentTouchableHandle->GetVolume() == 0) {
0538       fParticleChange.ProposeTrackStatus(fStopAndKill);
0539     }
0540     retCurrentTouchable = fCurrentTouchableHandle;
0541     fParticleChange.SetTouchableHandle(fCurrentTouchableHandle);
0542   }
0543   else  // fGeometryLimitedStep  is false
0544   {
0545     // This serves only to move the Navigator's location
0546     //
0547     fLinearNavigator->LocateGlobalPointWithinVolume(track.GetPosition());
0548 
0549     // The value of the track's current Touchable is retained.
0550     // (and it must be correct because we must use it below to
0551     // overwrite the (unset) one in particle change)
0552     //  It must be fCurrentTouchable too ??
0553     //
0554     fParticleChange.SetTouchableHandle(track.GetTouchableHandle());
0555     retCurrentTouchable = track.GetTouchableHandle();
0556   }  // endif ( fGeometryLimitedStep )
0557 
0558   const G4VPhysicalVolume* pNewVol = retCurrentTouchable->GetVolume();
0559   const G4Material* pNewMaterial = 0;
0560   G4VSensitiveDetector* pNewSensitiveDetector = 0;
0561   if (pNewVol != 0) {
0562     pNewMaterial = pNewVol->GetLogicalVolume()->GetMaterial();
0563     pNewSensitiveDetector = pNewVol->GetLogicalVolume()->GetSensitiveDetector();
0564   }
0565 
0566   // ( <const_cast> pNewMaterial ) ;
0567 
0568   fParticleChange.SetMaterialInTouchable((G4Material*)pNewMaterial);
0569   fParticleChange.SetSensitiveDetectorInTouchable(pNewSensitiveDetector);
0570 
0571   const G4MaterialCutsCouple* pNewMaterialCutsCouple = 0;
0572   if (pNewVol != 0) {
0573     pNewMaterialCutsCouple = pNewVol->GetLogicalVolume()->GetMaterialCutsCouple();
0574   }
0575 
0576   if (pNewVol != 0 && pNewMaterialCutsCouple != 0
0577       && pNewMaterialCutsCouple->GetMaterial() != pNewMaterial)
0578   {
0579     // for parametrized volume
0580     //
0581     pNewMaterialCutsCouple = G4ProductionCutsTable::GetProductionCutsTable()->GetMaterialCutsCouple(
0582       pNewMaterial, pNewMaterialCutsCouple->GetProductionCuts());
0583   }
0584   fParticleChange.SetMaterialCutsCoupleInTouchable(pNewMaterialCutsCouple);
0585 
0586   // temporarily until Get/Set Material of ParticleChange,
0587   // and StepPoint can be made const.
0588   // Set the touchable in ParticleChange
0589   // this must always be done because the particle change always
0590   // uses this value to overwrite the current touchable pointer.
0591   //
0592   fParticleChange.SetTouchableHandle(retCurrentTouchable);
0593 
0594   return &fParticleChange;
0595 }
0596 
0597 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
0598 
0599 // New method takes over the responsibility to reset the state
0600 // of G4MonopoleTransportation object at the start of a new track
0601 // or the resumption of a suspended track.
0602 
0603 void G4MonopoleTransportation::StartTracking(G4Track* aTrack)
0604 {
0605   G4VProcess::StartTracking(aTrack);
0606 
0607   // The actions here are those that were taken in AlongStepGPIL
0608   //   when track.GetCurrentStepNumber()==1
0609 
0610   // reset safety value and center
0611   //
0612   fPreviousSafety = 0.0;
0613   fPreviousSftOrigin = G4ThreeVector(0., 0., 0.);
0614 
0615   // reset looping counter -- for motion in field
0616   fNoLooperTrials = 0;
0617   // Must clear this state .. else it depends on last track's value
0618   //  --> a better solution would set this from state of suspended track TODO ?
0619   // Was if( aTrack->GetCurrentStepNumber()==1 ) { .. }
0620 
0621   // ChordFinder reset internal state
0622   //
0623   if (DoesGlobalFieldExist()) {
0624     fFieldPropagator->ClearPropagatorState();
0625     // Resets all state of field propagator class (ONLY)
0626     // including safety values
0627     // in case of overlaps and to wipe for first track).
0628 
0629     // G4ChordFinder* chordF= fFieldPropagator->GetChordFinder();
0630     // if( chordF ) chordF->ResetStepEstimate();
0631   }
0632 
0633   // Make sure to clear the chord finders of all fields (ie managers)
0634   G4FieldManagerStore* fieldMgrStore = G4FieldManagerStore::GetInstance();
0635   fieldMgrStore->ClearAllChordFindersState();
0636 
0637   // Update the current touchable handle  (from the track's)
0638   //
0639   fCurrentTouchableHandle = aTrack->GetTouchableHandle();
0640 }
0641 
0642 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......