Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/utils/src/G4VMultipleScattering.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 /processes/electromagnetic/utils/src/G4VMultipleScattering.cc (Version 11.3.0) and /processes/electromagnetic/utils/src/G4VMultipleScattering.cc (Version 5.2.p2)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                    <<   3 // * DISCLAIMER                                                       *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th <<   5 // * The following disclaimer summarizes all the specific disclaimers *
  6 // * the Geant4 Collaboration.  It is provided <<   6 // * of contributors to this software. The specific disclaimers,which *
  7 // * conditions of the Geant4 Software License <<   7 // * govern, are listed with their locations in:                      *
  8 // * LICENSE and available at  http://cern.ch/ <<   8 // *   http://cern.ch/geant4/license                                  *
  9 // * include a list of copyright holders.      << 
 10 // *                                                9 // *                                                                  *
 11 // * Neither the authors of this software syst     10 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     11 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     12 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     13 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file  <<  14 // * use.                                                             *
 16 // * for the full disclaimer and the limitatio << 
 17 // *                                               15 // *                                                                  *
 18 // * This  code  implementation is the result  <<  16 // * This  code  implementation is the  intellectual property  of the *
 19 // * technical work of the GEANT4 collaboratio <<  17 // * GEANT4 collaboration.                                            *
 20 // * By using,  copying,  modifying or  distri <<  18 // * By copying,  distributing  or modifying the Program (or any work *
 21 // * any work based  on the software)  you  ag <<  19 // * based  on  the Program)  you indicate  your  acceptance of  this *
 22 // * use  in  resulting  scientific  publicati <<  20 // * statement, and all its terms.                                    *
 23 // * acceptance of all terms of the Geant4 Sof << 
 24 // *******************************************     21 // ********************************************************************
 25 //                                                 22 //
 26 //                                                 23 //
 27 // -------------------------------------------     24 // -------------------------------------------------------------------
 28 //                                                 25 //
 29 // GEANT4 Class file                               26 // GEANT4 Class file
 30 //                                                 27 //
 31 //                                                 28 //
 32 // File name:     G4VMultipleScattering            29 // File name:     G4VMultipleScattering
 33 //                                                 30 //
 34 // Author:        Vladimir Ivanchenko on base      31 // Author:        Vladimir Ivanchenko on base of Laszlo Urban code
 35 //                                                 32 //
 36 // Creation date: 25.03.2003                       33 // Creation date: 25.03.2003
 37 //                                                 34 //
 38 // Modifications:                                  35 // Modifications:
                                                   >>  36 // 
                                                   >>  37 // 13.04.03 Change printout (V.Ivanchenko)
                                                   >>  38 // 04-06-03 Fix compilation warnings (V.Ivanchenko)
 39 //                                                 39 //
 40 // 16-07-03 Use G4VMscModel interface (V.Ivanc << 
 41 // 03-11-03 Fix initialisation problem in Retr << 
 42 // 04-11-03 Update PrintInfoDefinition (V.Ivan << 
 43 // 01-03-04 SampleCosineTheta signature change << 
 44 // 22-04-04 SampleCosineTheta signature change << 
 45 // 27-08-04 Add InitialiseForRun method (V.Iva << 
 46 // 08-11-04 Migration to new interface of Stor << 
 47 // 11-03-05 Shift verbose level by 1 (V.Ivantc << 
 48 // 15-04-05 optimize internal interface (V.Iva << 
 49 // 15-04-05 remove boundary flag (V.Ivanchenko << 
 50 // 27-10-05 introduce virtual function MscStep << 
 51 // 12-04-07 Add verbosity at destruction (V.Iv << 
 52 // 27-10-07 Virtual functions moved to source  << 
 53 // 11-03-08 Set skin value does not effect ste << 
 54 // 24-06-09 Removed hidden bin in G4PhysicsVec << 
 55 // 04-06-13 Adoptation to MT mode (V.Ivanchenk << 
 56 //                                                 40 //
                                                   >>  41 // Class Description:
                                                   >>  42 //
                                                   >>  43 // It is the generic process of multiple scattering it includes common
                                                   >>  44 // part of calculations for all charged particles
 57                                                    45 
 58 // -------------------------------------------     46 // -------------------------------------------------------------------
 59 //                                                 47 //
 60 //....oooOO0OOooo........oooOO0OOooo........oo     48 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 61 //....oooOO0OOooo........oooOO0OOooo........oo     49 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 62                                                    50 
 63 #include "G4VMultipleScattering.hh"                51 #include "G4VMultipleScattering.hh"
 64 #include "G4PhysicalConstants.hh"              << 
 65 #include "G4SystemOfUnits.hh"                  << 
 66 #include "G4LossTableManager.hh"                   52 #include "G4LossTableManager.hh"
 67 #include "G4MaterialCutsCouple.hh"             << 
 68 #include "G4Step.hh"                               53 #include "G4Step.hh"
 69 #include "G4ParticleDefinition.hh"                 54 #include "G4ParticleDefinition.hh"
                                                   >>  55 #include "G4VEmModel.hh"
 70 #include "G4VEmFluctuationModel.hh"                56 #include "G4VEmFluctuationModel.hh"
                                                   >>  57 #include "G4DataVector.hh"
                                                   >>  58 #include "G4PhysicsTable.hh"
                                                   >>  59 #include "G4PhysicsVector.hh"
                                                   >>  60 #include "G4PhysicsLogVector.hh"
 71 #include "G4UnitsTable.hh"                         61 #include "G4UnitsTable.hh"
 72 #include "G4ProductionCutsTable.hh"                62 #include "G4ProductionCutsTable.hh"
 73 #include "G4Electron.hh"                       <<  63 #include "G4Region.hh"
 74 #include "G4GenericIon.hh"                     <<  64 #include "G4RegionStore.hh"
                                                   >>  65 #include "G4Navigator.hh"
 75 #include "G4TransportationManager.hh"              66 #include "G4TransportationManager.hh"
 76 #include "G4SafetyHelper.hh"                   << 
 77 #include "G4ParticleTable.hh"                  << 
 78 #include "G4ProcessVector.hh"                  << 
 79 #include "G4ProcessManager.hh"                 << 
 80 #include "G4LossTableBuilder.hh"               << 
 81 #include "G4EmTableUtil.hh"                    << 
 82 #include <iostream>                            << 
 83                                                    67 
 84 //....oooOO0OOooo........oooOO0OOooo........oo     68 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 85                                                    69 
 86 G4VMultipleScattering::G4VMultipleScattering(c <<  70 G4VMultipleScattering::G4VMultipleScattering(const G4String& name, G4ProcessType type):
 87   : G4VContinuousDiscreteProcess("msc", fElect <<  71                  G4VContinuousDiscreteProcess(name, type),
 88   fNewPosition(0.,0.,0.),                      <<  72   navigator(0),
 89   fNewDirection(0.,0.,1.)                      <<  73   theLambdaTable(0),
                                                   >>  74   currentCouple(0),
                                                   >>  75   nBins(110),
                                                   >>  76   minKinEnergy(1.0*eV),
                                                   >>  77   maxKinEnergy(100.0*GeV),
                                                   >>  78   boundary(false),
                                                   >>  79   latDisplasment(true),
                                                   >>  80   buildLambdaTable(true)
 90 {                                                  81 {
 91   theParameters = G4EmParameters::Instance();  << 
 92   SetVerboseLevel(1);                          << 
 93   SetProcessSubType(fMultipleScattering);      << 
 94                                                << 
 95   lowestKinEnergy = 10*CLHEP::eV;              << 
 96                                                << 
 97   geomMin = 0.05*CLHEP::nm;                    << 
 98   minDisplacement2 = geomMin*geomMin;          << 
 99                                                << 
100   pParticleChange = &fParticleChange;          << 
101                                                << 
102   modelManager = new G4EmModelManager();           82   modelManager = new G4EmModelManager();
103   emManager = G4LossTableManager::Instance();  <<  83   (G4LossTableManager::Instance())->Register(this);
104   mscModels.reserve(2);                        <<  84 
105   emManager->Register(this);                   << 
106 }                                                  85 }
107                                                    86 
108 //....oooOO0OOooo........oooOO0OOooo........oo     87 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
109                                                    88 
110 G4VMultipleScattering::~G4VMultipleScattering(     89 G4VMultipleScattering::~G4VMultipleScattering()
111 {                                                  90 {
                                                   >>  91   (G4LossTableManager::Instance())->DeRegister(this);
112   delete modelManager;                             92   delete modelManager;
113   emManager->DeRegister(this);                 <<  93   if (theLambdaTable) {
                                                   >>  94     theLambdaTable->clearAndDestroy();
                                                   >>  95     delete theLambdaTable;
                                                   >>  96   }
                                                   >>  97   (G4LossTableManager::Instance())->DeRegister(this);
114 }                                                  98 }
115                                                    99 
116 //....oooOO0OOooo........oooOO0OOooo........oo    100 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
117                                                   101 
118 void G4VMultipleScattering::AddEmModel(G4int o << 102 void G4VMultipleScattering::BuildPhysicsTable(const G4ParticleDefinition& part)
119                                        const G << 
120 {                                                 103 {
121   if(nullptr == ptr) { return; }               << 104   currentCouple = 0;
122   G4VEmFluctuationModel* fm = nullptr;         << 105   if(0 < verboseLevel) {
123   modelManager->AddEmModel(order, ptr, fm, reg << 106     G4cout << "G4VMultipleScattering::BuildPhysicsTable() for "
124   ptr->SetParticleChange(pParticleChange);     << 107            << GetProcessName()
125 }                                              << 108            << " and particle " << part.GetParticleName()
                                                   >> 109            << G4endl;
                                                   >> 110   }
                                                   >> 111 
                                                   >> 112   G4bool cutsWasModified = false;
                                                   >> 113   const G4ProductionCutsTable* theCoupleTable=
                                                   >> 114         G4ProductionCutsTable::GetProductionCutsTable();
                                                   >> 115   size_t numOfCouples = theCoupleTable->GetTableSize();
                                                   >> 116 
                                                   >> 117   for (size_t j=0; j<numOfCouples; j++){
                                                   >> 118     if (theCoupleTable->GetMaterialCutsCouple(j)->IsRecalcNeeded()) {
                                                   >> 119       cutsWasModified = true;
                                                   >> 120       break;
                                                   >> 121     }
                                                   >> 122   }
                                                   >> 123   if( !cutsWasModified ) return;
                                                   >> 124   InitialiseProcess(part);
126                                                   125 
127 //....oooOO0OOooo........oooOO0OOooo........oo << 126   if(latDisplasment) navigator = G4TransportationManager::GetTransportationManager()
                                                   >> 127              ->GetNavigatorForTracking();
                                                   >> 128   modelManager->Clear();
                                                   >> 129   const G4DataVector* theCuts = modelManager->Initialise(&part, 0, 10.0, verboseLevel);
128                                                   130 
129 void G4VMultipleScattering::SetEmModel(G4VMscM << 131   if (buildLambdaTable) {
130 {                                              << 
131   if(nullptr == ptr) { return; }               << 
132   if(!mscModels.empty()) {                     << 
133     for(auto & msc : mscModels) { if(msc == pt << 
134   }                                            << 
135   mscModels.push_back(ptr);                    << 
136 }                                              << 
137                                                   132 
138 //....oooOO0OOooo........oooOO0OOooo........oo << 133     theLambdaTable = new G4PhysicsTable(numOfCouples);
139                                                   134 
140 void                                           << 135     for (size_t i=0; i<numOfCouples; i++) {
141 G4VMultipleScattering::PreparePhysicsTable(con << 
142 {                                              << 
143   G4bool master = emManager->IsMaster();       << 
144   if (nullptr == firstParticle) { firstParticl << 
145                                                   136 
146   emManager->PreparePhysicsTable(&part, this); << 137       // create physics vector and fill it
147   currParticle = nullptr;                      << 138       const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(i);
                                                   >> 139       G4PhysicsVector* aVector = PhysicsVector(couple);
                                                   >> 140       modelManager->FillLambdaVector(aVector, couple, false);
148                                                   141 
149   if(firstParticle == &part) {                 << 142       // Insert vector for this material into the table
150     baseMat = emManager->GetTableBuilder()->Ge << 143       theLambdaTable->insert(aVector) ;
151     G4EmTableUtil::PrepareMscProcess(this, par << 
152              stepLimit, facrange,              << 
153              latDisplacement, master,          << 
154              isIon, baseMat);                  << 
155                                                << 
156     numberOfModels = modelManager->NumberOfMod << 
157     currentModel = GetModelByIndex(0);         << 
158                                                << 
159     if (nullptr == safetyHelper) {             << 
160       safetyHelper = G4TransportationManager:: << 
161   ->GetSafetyHelper();                         << 
162       safetyHelper->InitialiseHelper();        << 
163     }                                             144     }
164   }                                            << 
165 }                                              << 
166                                                   145 
167 //....oooOO0OOooo........oooOO0OOooo........oo << 146     if(0 < verboseLevel) {
168                                                << 147       G4cout << "Lambda table is built for "
169 void G4VMultipleScattering::BuildPhysicsTable( << 148              << part.GetParticleName()
170 {                                              << 149              << G4endl;
171   G4bool master = emManager->IsMaster();       << 150     }
                                                   >> 151     if(2 < verboseLevel) G4cout << *theLambdaTable << G4endl;
                                                   >> 152     if(5 < verboseLevel) G4cout << theCuts << G4endl;
172                                                   153 
173   if(firstParticle == &part) {                 << 
174     emManager->BuildPhysicsTable(&part);       << 
175   }                                            << 
176   const G4VMultipleScattering* ptr = this;     << 
177   if(!master) {                                << 
178     ptr = static_cast<const G4VMultipleScatter << 
179   }                                               154   }
180                                                   155 
181   G4EmTableUtil::BuildMscProcess(this, ptr, pa << 156   G4String num = part.GetParticleName();
182          numberOfModels, master);              << 157   if (verboseLevel>0 || num == "e-" || num == "mu+" || num == "proton")
                                                   >> 158            PrintInfoDefinition();
                                                   >> 159 
                                                   >> 160   if(0 < verboseLevel) {
                                                   >> 161     G4cout << "G4VMultipleScattering::BuildPhysicsTable() done for "
                                                   >> 162            << GetProcessName()
                                                   >> 163            << " and particle " << part.GetParticleName()
                                                   >> 164            << G4endl;
                                                   >> 165   }
183 }                                                 166 }
184                                                   167 
185 //....oooOO0OOooo........oooOO0OOooo........oo    168 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
186                                                   169 
187 void G4VMultipleScattering::StreamInfo(std::os << 170 void G4VMultipleScattering::AddEmModel(G4int order, G4VEmModel* p,
188                   const G4ParticleDefinition&  << 171                                  const G4Region* region)
189 {                                                 172 {
190   G4String indent = (rst ? "  " : "");         << 173   G4VEmFluctuationModel* fm = 0;
191   outFile << G4endl << indent << GetProcessNam << 174   modelManager->AddEmModel(order, p, fm, region);
192   if (!rst) outFile << " for " << part.GetPart << 
193   outFile  << "  SubType= " << GetProcessSubTy << 
194   modelManager->DumpModelList(outFile, verbose << 
195 }                                                 175 }
196                                                   176 
197 //....oooOO0OOooo........oooOO0OOooo........oo    177 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
198                                                   178 
199 void G4VMultipleScattering::StartTracking(G4Tr << 179 G4VParticleChange* G4VMultipleScattering::PostStepDoIt(const G4Track& track,
                                                   >> 180                                                        const G4Step& step)
200 {                                                 181 {
201   G4VEnergyLossProcess* eloss = nullptr;       << 182   fParticleChange.Initialize(track);
202   if(track->GetParticleDefinition() != currPar << 183   G4double kineticEnergy = track.GetKineticEnergy();
203     currParticle = track->GetParticleDefinitio << 184   G4double truestep = step.GetStepLength();
204     fIonisation = emManager->GetEnergyLossProc << 
205     eloss = fIonisation;                       << 
206   }                                            << 
207   for(G4int i=0; i<numberOfModels; ++i) {      << 
208     G4VMscModel* msc = GetModelByIndex(i);     << 
209     msc->StartTracking(track);                 << 
210     if(nullptr != eloss) {                     << 
211       msc->SetIonisation(eloss, currParticle); << 
212     }                                          << 
213   }                                            << 
214 }                                              << 
215                                                   185 
216 //....oooOO0OOooo........oooOO0OOooo........oo << 186   if (kineticEnergy > 0.0) {
                                                   >> 187     G4double cth  = currentModel->SampleCosineTheta(truestep);
                                                   >> 188     G4double sth  = sqrt(1.-cth*cth);
                                                   >> 189     G4double phi  = twopi*G4UniformRand();
                                                   >> 190     G4double dirx = sth*cos(phi);
                                                   >> 191     G4double diry = sth*sin(phi);
                                                   >> 192 
                                                   >> 193     G4ThreeVector oldDirection = track.GetMomentumDirection();
                                                   >> 194     G4ThreeVector newDirection(dirx,diry,cth);
                                                   >> 195     newDirection.rotateUz(oldDirection);
                                                   >> 196     fParticleChange.SetMomentumChange(newDirection);
217                                                   197 
218 G4double G4VMultipleScattering::AlongStepGetPh << 
219                              const G4Track& tr << 
220                              G4double,         << 
221                              G4double currentM << 
222                              G4double&,        << 
223                              G4GPILSelection*  << 
224 {                                              << 
225   // get Step limit proposed by the process    << 
226   *selection = NotCandidateForSelection;       << 
227   physStepLimit = gPathLength = tPathLength =  << 
228                                                   198 
229   G4double ekin = track.GetKineticEnergy();    << 
230   /*                                              199   /*
231   G4cout << "MSC::AlongStepGPIL: Ekin= " << ek << 200   if(0 < verboseLevel) {
232          << "  " << currParticle->GetParticleN << 201     const G4ParticleDefinition* pd = dynParticle->GetDefinition();
233          << " currMod " << currentModel        << 202     G4cout << "G4VMultipleScattering::PostStepDoIt: Sample secondary; E= " << finalT/MeV
234          << G4endl;                            << 203            << " MeV; model= (" << currentModel->LowEnergyLimit(pd)
235   */                                           << 204            << ", " <<  currentModel->HighEnergyLimit(pd) << ")"
236   // isIon flag is used only to select a model << 205            << G4endl;
237   if(isIon) {                                  << 
238     ekin *= proton_mass_c2/track.GetParticleDe << 
239   }                                            << 
240   const G4MaterialCutsCouple* couple = track.G << 
241                                                << 
242   // select new model, static cast is possible << 
243   if(1 < numberOfModels) {                     << 
244     currentModel =                             << 
245       static_cast<G4VMscModel*>(SelectModel(ek << 
246   }                                            << 
247   currentModel->SetCurrentCouple(couple);      << 
248   // msc is active is model is active, energy  << 
249   // and step size is above the limit;         << 
250   // if it is active msc may limit the step    << 
251   if(currentModel->IsActive(ekin) && tPathLeng << 
252      && ekin >= lowestKinEnergy) {             << 
253     isActive = true;                           << 
254     tPathLength =                              << 
255       currentModel->ComputeTruePathLengthLimit << 
256     if (tPathLength < physStepLimit) {         << 
257       *selection = CandidateForSelection;      << 
258     }                                          << 
259   } else {                                     << 
260     isActive = false;                          << 
261     gPathLength = DBL_MAX;                     << 
262   }                                               206   }
263                                                << 
264   //if(currParticle->GetPDGMass() > GeV)       << 
265   /*                                           << 
266   G4cout << "MSC::AlongStepGPIL: Ekin= " << ek << 
267          << "  " << currParticle->GetParticleN << 
268          << " gPathLength= " << gPathLength    << 
269          << " tPathLength= " << tPathLength    << 
270          << " currentMinimalStep= " << current << 
271          << " isActive " << isActive << G4endl << 
272   */                                              207   */
273   return gPathLength;                          << 
274 }                                              << 
275                                                   208 
276 //....oooOO0OOooo........oooOO0OOooo........oo << 209     //    G4cout << "PostStep: sth= " << sth << " trueLength= " << truestep << " tLast= " << truePathLength << G4endl;
277                                                   210 
278 G4double                                       << 211     if (latDisplasment) {
279 G4VMultipleScattering::PostStepGetPhysicalInte << 
280               const G4Track&, G4double, G4Forc << 
281 {                                              << 
282   *condition = NotForced;                      << 
283   return DBL_MAX;                              << 
284 }                                              << 
285                                                   212 
286 //....oooOO0OOooo........oooOO0OOooo........oo << 213       G4double safety = step.GetPostStepPoint()->GetSafety();
                                                   >> 214       if ( safety > 0.0) {
                                                   >> 215         G4double r = currentModel->SampleDisplacement();
                                                   >> 216         if (r > safety) r = safety;
287                                                   217 
288 G4VParticleChange*                             << 218   //    G4cout << "r= " << r << " safety= " << safety << G4endl;
289 G4VMultipleScattering::AlongStepDoIt(const G4T << 
290 {                                              << 
291   fParticleChange.InitialiseMSC(track, step);  << 
292   fNewPosition = fParticleChange.GetProposedPo << 
293   fPositionChanged = false;                    << 
294                                                   219 
295   G4double geomLength = step.GetStepLength();  << 220         // sample direction of lateral displacement
                                                   >> 221         G4double phi  = twopi*G4UniformRand();
                                                   >> 222         G4double dirx = cos(phi);
                                                   >> 223         G4double diry = sin(phi);
296                                                   224 
297   // very small step - no msc                  << 225         G4ThreeVector latDirection(dirx,diry,0.0);
298   if(!isActive) {                              << 226         latDirection.rotateUz(oldDirection);
299     tPathLength = geomLength;                  << 
300                                                   227 
301     // sample msc                              << 228         // compute new endpoint of the Step
302   } else {                                     << 229         G4ThreeVector newPosition = (step.GetPostStepPoint())->GetPosition()
303     G4double range =                           << 230               + r*latDirection;
304       currentModel->GetRange(currParticle,trac << 231 
305                              track.GetMaterial << 232         navigator->LocateGlobalPointWithinVolume(newPosition);
306                                                << 233 
307     tPathLength = currentModel->ComputeTrueSte << 234         fParticleChange.SetPositionChange(newPosition);
308                                                << 
309     /*                                         << 
310     if(currParticle->GetPDGMass() > 0.9*GeV)   << 
311     G4cout << "G4VMsc::AlongStepDoIt: GeomLeng << 
312            << geomLength                       << 
313            << " tPathLength= " << tPathLength  << 
314            << " physStepLimit= " << physStepLi << 
315            << " dr= " << range - tPathLength   << 
316            << " ekin= " << track.GetKineticEne << 
317     */                                         << 
318     // protection against wrong t->g->t conver << 
319     tPathLength = std::min(tPathLength, physSt << 
320                                                << 
321     // do not sample scattering at the last or << 
322     if(tPathLength < range && tPathLength > ge << 
323                                                << 
324       static const G4double minSafety = 1.20*C << 
325       static const G4double sFact = 0.99;      << 
326                                                << 
327       G4ThreeVector displacement = currentMode << 
328         step.GetPostStepPoint()->GetMomentumDi << 
329                                                << 
330       G4double r2 = displacement.mag2();       << 
331       //G4cout << "    R= " << sqrt(r2) << " R << 
332       //     << " flag= " << fDispBeyondSafety << 
333       if(r2 > minDisplacement2) {              << 
334                                                << 
335         fPositionChanged = true;               << 
336         G4double dispR = std::sqrt(r2);        << 
337         G4double postSafety =                  << 
338           sFact*safetyHelper->ComputeSafety(fN << 
339         //G4cout<<"    R= "<< dispR<<" postSaf << 
340                                                << 
341         // far away from geometry boundary     << 
342         if(postSafety > 0.0 && dispR <= postSa << 
343           fNewPosition += displacement;        << 
344                                                << 
345           //near the boundary                  << 
346         } else {                               << 
347           // displaced point is definitely wit << 
348           //G4cout<<"    R= "<<dispR<<" postSa << 
349           if(dispR < postSafety) {             << 
350             fNewPosition += displacement;      << 
351                                                << 
352             // reduced displacement            << 
353           } else if(postSafety > geomMin) {    << 
354             fNewPosition += displacement*(post << 
355                                                << 
356             // very small postSafety           << 
357           } else {                             << 
358             fPositionChanged = false;          << 
359           }                                    << 
360         }                                      << 
361         if(fPositionChanged) {                 << 
362           safetyHelper->ReLocateWithinVolume(f << 
363           fParticleChange.ProposePosition(fNew << 
364         }                                      << 
365       }                                           235       }
366     }                                             236     }
367   }                                               237   }
368   fParticleChange.ProposeTrueStepLength(tPathL << 
369   return &fParticleChange;                     << 
370 }                                              << 
371                                                   238 
372 //....oooOO0OOooo........oooOO0OOooo........oo << 239   return &fParticleChange;
373                                                << 
374 G4double G4VMultipleScattering::GetContinuousS << 
375                                        const G << 
376                                        G4doubl << 
377                                        G4doubl << 
378                                        G4doubl << 
379 {                                              << 
380   G4GPILSelection selection = NotCandidateForS << 
381   G4double x = AlongStepGetPhysicalInteraction << 
382                                                << 
383                                                << 
384                                                << 
385   return x;                                    << 
386 }                                                 240 }
387                                                   241 
388 //....oooOO0OOooo........oooOO0OOooo........oo    242 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
389                                                   243 
390 G4double G4VMultipleScattering::ContinuousStep << 244 void G4VMultipleScattering::PrintInfoDefinition()
391                                        const G << 
392                                        G4doubl << 
393                                        G4doubl << 
394                                        G4doubl << 
395 {                                                 245 {
396   return GetContinuousStepLimit(track,previous << 246   G4cout << G4endl << GetProcessName() << ":  Model variant of multiple scattering " << G4endl;
397                                 currentSafety) << 247   if (theLambdaTable) {
                                                   >> 248     G4cout << "      Lambda tables from "
                                                   >> 249            << G4BestUnit(MinKinEnergy(),"Energy")
                                                   >> 250      << " to "
                                                   >> 251            << G4BestUnit(MaxKinEnergy(),"Energy")
                                                   >> 252            << " in " << nBins << " bins."
                                                   >> 253            << G4endl;
                                                   >> 254   }
                                                   >> 255   if (1 < verboseLevel) {
                                                   >> 256       G4cout << "LambdaTable address= " << theLambdaTable << G4endl;
                                                   >> 257       if(theLambdaTable) G4cout << (*theLambdaTable) << G4endl;
                                                   >> 258   }
398 }                                                 259 }
399                                                   260 
400 //....oooOO0OOooo........oooOO0OOooo........oo    261 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
401                                                   262 
402 G4double G4VMultipleScattering::GetMeanFreePat << 263 G4PhysicsVector* G4VMultipleScattering::PhysicsVector(const G4MaterialCutsCouple* couple)
403               const G4Track&, G4double, G4Forc << 
404 {                                                 264 {
405   *condition = Forced;                         << 265   G4int nbins = 3;
406   return DBL_MAX;                              << 266   //G4int nbins =  nDEDXBins;
                                                   >> 267   if( couple->IsUsed() ) nbins = nBins;
                                                   >> 268   //  G4double xmax = maxKinEnergy*exp( log(maxKinEnergy/minKinEnergy) / ((G4double)(nbins-1)) );
                                                   >> 269   G4PhysicsVector* v = new G4PhysicsLogVector(minKinEnergy, maxKinEnergy, nbins);
                                                   >> 270   return v;
407 }                                                 271 }
408                                                   272 
409 //....oooOO0OOooo........oooOO0OOooo........oo << 273 G4bool G4VMultipleScattering::StorePhysicsTable(G4ParticleDefinition* part,
410                                                << 274              const G4String& directory,
411 G4bool                                         << 275                    G4bool ascii)
412 G4VMultipleScattering::StorePhysicsTable(const << 
413                                          const << 
414                                          G4boo << 
415 {                                                 276 {
416   G4bool yes = true;                           << 277   G4bool res = true;
417   if(part != firstParticle || !emManager->IsMa << 278   if ( theLambdaTable ) {
                                                   >> 279     const G4String name = GetPhysicsTableFileName(part,directory,"Lambda",ascii);
                                                   >> 280     G4bool yes = theLambdaTable->StorePhysicsTable(name,ascii);
                                                   >> 281     if( !yes ) res = false;
                                                   >> 282   }
418                                                   283 
419   return G4EmTableUtil::StoreMscTable(this, pa << 284   if ( res ) {
420               numberOfModels, verboseLevel,    << 285     G4cout << "Physics table are stored for " << part->GetParticleName()
421                                       ascii);  << 286            << " and process " << GetProcessName()
                                                   >> 287      << " in the directory <" << directory
                                                   >> 288      << "> " << G4endl;
                                                   >> 289   } else {
                                                   >> 290     G4cout << "Fail to store Physics Table for " << part->GetParticleName()
                                                   >> 291            << " and process " << GetProcessName()
                                                   >> 292      << " in the directory <" << directory
                                                   >> 293      << "> " << G4endl;
                                                   >> 294   }
                                                   >> 295   return res;
422 }                                                 296 }
423                                                   297 
424 //....oooOO0OOooo........oooOO0OOooo........oo    298 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
425                                                   299 
426 G4bool                                         << 300 G4bool G4VMultipleScattering::RetrievePhysicsTable(G4ParticleDefinition* part,
427 G4VMultipleScattering::RetrievePhysicsTable(co << 301                   const G4String& directory,
428                                             co << 302                         G4bool ascii)
429                                             G4 << 303 {
430 {                                              << 304   if(0 < verboseLevel) {
431   return true;                                 << 305     G4cout << "G4VMultipleScattering::RetrievePhysicsTable() for "
432 }                                              << 306            << part->GetParticleName() << " and process "
433                                                << 307      << GetProcessName() << G4endl;
434 //....oooOO0OOooo........oooOO0OOooo........oo << 308   }
435                                                << 309   G4bool res = true;
436 void G4VMultipleScattering::ProcessDescription << 310   if(!buildLambdaTable) return res;
437 {                                              << 311 
438   if(nullptr != firstParticle) {               << 312   G4String num = part->GetParticleName();
439     StreamInfo(outFile, *firstParticle, true); << 313   const G4ProductionCutsTable* theCoupleTable=
                                                   >> 314           G4ProductionCutsTable::GetProductionCutsTable();
                                                   >> 315   size_t numOfCouples = theCoupleTable->GetTableSize();
                                                   >> 316 
                                                   >> 317   G4String filename = GetPhysicsTableFileName(part,directory,"Lambda",ascii);
                                                   >> 318   theLambdaTable = new G4PhysicsTable(numOfCouples);
                                                   >> 319   res = theLambdaTable->RetrievePhysicsTable(filename,ascii);
                                                   >> 320   if ( res ) {
                                                   >> 321     if (0 < verboseLevel) {
                                                   >> 322         G4cout << "Lambda table for " << num << " is retrieved from <"
                                                   >> 323                << filename << ">"
                                                   >> 324                << G4endl;
                                                   >> 325     }
                                                   >> 326   } else {
                                                   >> 327     theLambdaTable->clearAndDestroy();
                                                   >> 328     theLambdaTable = 0;
                                                   >> 329     if (0 < verboseLevel) {
                                                   >> 330         G4cout << "Lambda table for " << num << " in file <"
                                                   >> 331                << filename << "> is not exist"
                                                   >> 332                << G4endl;
                                                   >> 333     }
440   }                                               334   }
                                                   >> 335 
                                                   >> 336   if (verboseLevel>0 || num == "e-" || num == "mu+" || num == "proton")
                                                   >> 337            PrintInfoDefinition();
                                                   >> 338   return res;
441 }                                                 339 }
442                                                   340 
443 //....oooOO0OOooo........oooOO0OOooo........oo    341 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
444                                                   342 
445                                                   343