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 9.5)


  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 // $Id: G4VMultipleScattering.cc,v 1.86 2010-10-26 11:30:46 vnivanch Exp $
                                                   >>  27 // GEANT4 tag $Name: not supported by cvs2svn $
 26 //                                                 28 //
 27 // -------------------------------------------     29 // -------------------------------------------------------------------
 28 //                                                 30 //
 29 // GEANT4 Class file                               31 // GEANT4 Class file
 30 //                                                 32 //
 31 //                                                 33 //
 32 // File name:     G4VMultipleScattering            34 // File name:     G4VMultipleScattering
 33 //                                                 35 //
 34 // Author:        Vladimir Ivanchenko on base      36 // Author:        Vladimir Ivanchenko on base of Laszlo Urban code
 35 //                                                 37 //
 36 // Creation date: 25.03.2003                       38 // Creation date: 25.03.2003
 37 //                                                 39 //
 38 // Modifications:                                  40 // Modifications:
 39 //                                                 41 //
                                                   >>  42 // 13.04.03 Change printout (V.Ivanchenko)
                                                   >>  43 // 04-06-03 Fix compilation warnings (V.Ivanchenko)
 40 // 16-07-03 Use G4VMscModel interface (V.Ivanc     44 // 16-07-03 Use G4VMscModel interface (V.Ivanchenko)
 41 // 03-11-03 Fix initialisation problem in Retr     45 // 03-11-03 Fix initialisation problem in RetrievePhysicsTable (V.Ivanchenko)
 42 // 04-11-03 Update PrintInfoDefinition (V.Ivan     46 // 04-11-03 Update PrintInfoDefinition (V.Ivanchenko)
 43 // 01-03-04 SampleCosineTheta signature change     47 // 01-03-04 SampleCosineTheta signature changed
 44 // 22-04-04 SampleCosineTheta signature change     48 // 22-04-04 SampleCosineTheta signature changed back to original
 45 // 27-08-04 Add InitialiseForRun method (V.Iva     49 // 27-08-04 Add InitialiseForRun method (V.Ivanchneko)
 46 // 08-11-04 Migration to new interface of Stor     50 // 08-11-04 Migration to new interface of Store/Retrieve tables (V.Ivantchenko)
 47 // 11-03-05 Shift verbose level by 1 (V.Ivantc     51 // 11-03-05 Shift verbose level by 1 (V.Ivantchenko)
 48 // 15-04-05 optimize internal interface (V.Iva     52 // 15-04-05 optimize internal interface (V.Ivanchenko)
 49 // 15-04-05 remove boundary flag (V.Ivanchenko     53 // 15-04-05 remove boundary flag (V.Ivanchenko)
 50 // 27-10-05 introduce virtual function MscStep     54 // 27-10-05 introduce virtual function MscStepLimitation() (V.Ivanchenko)
 51 // 12-04-07 Add verbosity at destruction (V.Iv     55 // 12-04-07 Add verbosity at destruction (V.Ivanchenko)
 52 // 27-10-07 Virtual functions moved to source      56 // 27-10-07 Virtual functions moved to source (V.Ivanchenko)
 53 // 11-03-08 Set skin value does not effect ste     57 // 11-03-08 Set skin value does not effect step limit type (V.Ivanchenko)
 54 // 24-06-09 Removed hidden bin in G4PhysicsVec     58 // 24-06-09 Removed hidden bin in G4PhysicsVector (V.Ivanchenko)
 55 // 04-06-13 Adoptation to MT mode (V.Ivanchenk << 
 56 //                                                 59 //
                                                   >>  60 // Class Description:
                                                   >>  61 //
                                                   >>  62 // It is the generic process of multiple scattering it includes common
                                                   >>  63 // part of calculations for all charged particles
 57                                                    64 
 58 // -------------------------------------------     65 // -------------------------------------------------------------------
 59 //                                                 66 //
 60 //....oooOO0OOooo........oooOO0OOooo........oo     67 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 61 //....oooOO0OOooo........oooOO0OOooo........oo     68 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 62                                                    69 
 63 #include "G4VMultipleScattering.hh"                70 #include "G4VMultipleScattering.hh"
 64 #include "G4PhysicalConstants.hh"              << 
 65 #include "G4SystemOfUnits.hh"                  << 
 66 #include "G4LossTableManager.hh"                   71 #include "G4LossTableManager.hh"
 67 #include "G4MaterialCutsCouple.hh"             <<  72 #include "G4LossTableBuilder.hh"
 68 #include "G4Step.hh"                               73 #include "G4Step.hh"
 69 #include "G4ParticleDefinition.hh"                 74 #include "G4ParticleDefinition.hh"
 70 #include "G4VEmFluctuationModel.hh"                75 #include "G4VEmFluctuationModel.hh"
                                                   >>  76 #include "G4DataVector.hh"
                                                   >>  77 #include "G4PhysicsTable.hh"
                                                   >>  78 #include "G4PhysicsVector.hh"
                                                   >>  79 #include "G4PhysicsLogVector.hh"
 71 #include "G4UnitsTable.hh"                         80 #include "G4UnitsTable.hh"
 72 #include "G4ProductionCutsTable.hh"                81 #include "G4ProductionCutsTable.hh"
 73 #include "G4Electron.hh"                       <<  82 #include "G4Region.hh"
                                                   >>  83 #include "G4RegionStore.hh"
                                                   >>  84 #include "G4PhysicsTableHelper.hh"
 74 #include "G4GenericIon.hh"                         85 #include "G4GenericIon.hh"
 75 #include "G4TransportationManager.hh"          <<  86 #include "G4Electron.hh"
 76 #include "G4SafetyHelper.hh"                   <<  87 #include "G4EmConfigurator.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                                                    88 
 84 //....oooOO0OOooo........oooOO0OOooo........oo     89 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
 85                                                    90 
 86 G4VMultipleScattering::G4VMultipleScattering(c <<  91 G4VMultipleScattering::G4VMultipleScattering(const G4String& name, 
 87   : G4VContinuousDiscreteProcess("msc", fElect <<  92                G4ProcessType type):
 88   fNewPosition(0.,0.,0.),                      <<  93   G4VContinuousDiscreteProcess(name, type),
 89   fNewDirection(0.,0.,1.)                      <<  94   buildLambdaTable(true),
                                                   >>  95   theLambdaTable(0),
                                                   >>  96   firstParticle(0),
                                                   >>  97   theDensityFactor(0),
                                                   >>  98   theDensityIdx(0),
                                                   >>  99   stepLimit(fUseSafety),
                                                   >> 100   skin(1.0),
                                                   >> 101   facrange(0.04),
                                                   >> 102   facgeom(2.5),
                                                   >> 103   latDisplasment(true),
                                                   >> 104   isIon(false),
                                                   >> 105   currentParticle(0),
                                                   >> 106   currentCouple(0)
 90 {                                                 107 {
 91   theParameters = G4EmParameters::Instance();  << 
 92   SetVerboseLevel(1);                             108   SetVerboseLevel(1);
 93   SetProcessSubType(fMultipleScattering);         109   SetProcessSubType(fMultipleScattering);
 94                                                   110 
 95   lowestKinEnergy = 10*CLHEP::eV;              << 111   // Size of tables assuming spline
                                                   >> 112   minKinEnergy = 0.1*keV;
                                                   >> 113   maxKinEnergy = 10.0*TeV;
                                                   >> 114   nBins        = 77;
 96                                                   115 
 97   geomMin = 0.05*CLHEP::nm;                    << 116   // default limit on polar angle
 98   minDisplacement2 = geomMin*geomMin;          << 117   polarAngleLimit = 0.0;
 99                                                   118 
100   pParticleChange = &fParticleChange;             119   pParticleChange = &fParticleChange;
101                                                   120 
102   modelManager = new G4EmModelManager();          121   modelManager = new G4EmModelManager();
103   emManager = G4LossTableManager::Instance();  << 122   (G4LossTableManager::Instance())->Register(this);
104   mscModels.reserve(2);                        << 123 
105   emManager->Register(this);                   << 
106 }                                                 124 }
107                                                   125 
108 //....oooOO0OOooo........oooOO0OOooo........oo    126 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
109                                                   127 
110 G4VMultipleScattering::~G4VMultipleScattering(    128 G4VMultipleScattering::~G4VMultipleScattering()
111 {                                                 129 {
                                                   >> 130   if(1 < verboseLevel) {
                                                   >> 131     G4cout << "G4VMultipleScattering destruct " << GetProcessName() 
                                                   >> 132      << G4endl;
                                                   >> 133   }
112   delete modelManager;                            134   delete modelManager;
113   emManager->DeRegister(this);                 << 135   if (theLambdaTable) {
                                                   >> 136     theLambdaTable->clearAndDestroy();
                                                   >> 137     delete theLambdaTable;
                                                   >> 138   }
                                                   >> 139   (G4LossTableManager::Instance())->DeRegister(this);
114 }                                                 140 }
115                                                   141 
116 //....oooOO0OOooo........oooOO0OOooo........oo    142 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
117                                                   143 
118 void G4VMultipleScattering::AddEmModel(G4int o << 144 void G4VMultipleScattering::AddEmModel(G4int order, G4VEmModel* p,
119                                        const G << 145                const G4Region* region)
120 {                                                 146 {
121   if(nullptr == ptr) { return; }               << 147   G4VEmFluctuationModel* fm = 0;
122   G4VEmFluctuationModel* fm = nullptr;         << 148   modelManager->AddEmModel(order, p, fm, region);
123   modelManager->AddEmModel(order, ptr, fm, reg << 149   if(p) { p->SetParticleChange(pParticleChange); }
124   ptr->SetParticleChange(pParticleChange);     << 
125 }                                                 150 }
126                                                << 
127 //....oooOO0OOooo........oooOO0OOooo........oo    151 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
128                                                   152 
129 void G4VMultipleScattering::SetEmModel(G4VMscM << 153 void G4VMultipleScattering::SetModel(G4VMscModel* p, G4int index)
130 {                                                 154 {
131   if(nullptr == ptr) { return; }               << 155   G4int n = mscModels.size();
132   if(!mscModels.empty()) {                     << 156   if(index >= n) { for(G4int i=n; i<=index; ++i) {mscModels.push_back(0);} }
133     for(auto & msc : mscModels) { if(msc == pt << 157   mscModels[index] = p;
134   }                                            << 
135   mscModels.push_back(ptr);                    << 
136 }                                                 158 }
137                                                   159 
138 //....oooOO0OOooo........oooOO0OOooo........oo    160 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
139                                                   161 
140 void                                           << 162 G4VMscModel* G4VMultipleScattering::Model(G4int index)
141 G4VMultipleScattering::PreparePhysicsTable(con << 
142 {                                                 163 {
143   G4bool master = emManager->IsMaster();       << 164   G4VMscModel* p = 0;
144   if (nullptr == firstParticle) { firstParticl << 165   if(index >= 0 && index <  G4int(mscModels.size())) { p = mscModels[index]; }
                                                   >> 166   return p;
                                                   >> 167 }
145                                                   168 
146   emManager->PreparePhysicsTable(&part, this); << 169 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
147   currParticle = nullptr;                      << 
148                                                   170 
149   if(firstParticle == &part) {                 << 171 G4VEmModel* 
150     baseMat = emManager->GetTableBuilder()->Ge << 172 G4VMultipleScattering::GetModelByIndex(G4int idx, G4bool ver) const
151     G4EmTableUtil::PrepareMscProcess(this, par << 173 {
152              stepLimit, facrange,              << 174   return modelManager->GetModel(idx, ver);
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     }                                          << 
164   }                                            << 
165 }                                                 175 }
166                                                   176 
167 //....oooOO0OOooo........oooOO0OOooo........oo    177 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
168                                                   178 
169 void G4VMultipleScattering::BuildPhysicsTable( << 179 void G4VMultipleScattering::PreparePhysicsTable(const G4ParticleDefinition& part)
170 {                                                 180 {
171   G4bool master = emManager->IsMaster();       << 181   if (!firstParticle) {
172                                                << 182     currentCouple = 0;
173   if(firstParticle == &part) {                 << 183     if(part.GetParticleType() == "nucleus" && 
174     emManager->BuildPhysicsTable(&part);       << 184        part.GetParticleSubType() == "generic") {
                                                   >> 185       firstParticle = G4GenericIon::GenericIon();
                                                   >> 186       isIon = true;
                                                   >> 187     } else {
                                                   >> 188       firstParticle = &part;
                                                   >> 189       if(part.GetParticleType() == "nucleus" || 
                                                   >> 190    part.GetPDGMass() > GeV) {isIon = true;} 
                                                   >> 191     } 
                                                   >> 192     // limitations for ions
                                                   >> 193     if(isIon) {
                                                   >> 194       SetStepLimitType(fMinimal);
                                                   >> 195       SetLateralDisplasmentFlag(false);
                                                   >> 196       SetBuildLambdaTable(false);
                                                   >> 197     }
                                                   >> 198     currentParticle = &part;
175   }                                               199   }
176   const G4VMultipleScattering* ptr = this;     << 200 
177   if(!master) {                                << 201   G4LossTableManager* man = G4LossTableManager::Instance();
178     ptr = static_cast<const G4VMultipleScatter << 202   man->PreparePhysicsTable(&part, this);
                                                   >> 203 
                                                   >> 204   if(1 < verboseLevel) {
                                                   >> 205     G4cout << "### G4VMultipleScattering::PrepearPhysicsTable() for "
                                                   >> 206            << GetProcessName()
                                                   >> 207            << " and particle " << part.GetParticleName()
                                                   >> 208            << " local particle " << firstParticle->GetParticleName()
                                                   >> 209            << G4endl;
179   }                                               210   }
180                                                   211 
181   G4EmTableUtil::BuildMscProcess(this, ptr, pa << 212   if(firstParticle == &part) {
182          numberOfModels, master);              << 213 
                                                   >> 214     InitialiseProcess(firstParticle);
                                                   >> 215 
                                                   >> 216 
                                                   >> 217     // initialisation of models
                                                   >> 218     G4int nmod = modelManager->NumberOfModels();
                                                   >> 219     for(G4int i=0; i<nmod; ++i) {
                                                   >> 220       G4VMscModel* msc = static_cast<G4VMscModel*>(modelManager->GetModel(i));
                                                   >> 221       if(isIon) {
                                                   >> 222   msc->SetStepLimitType(fMinimal);
                                                   >> 223   msc->SetLateralDisplasmentFlag(false);
                                                   >> 224   msc->SetRangeFactor(0.2);
                                                   >> 225       } else {
                                                   >> 226   msc->SetStepLimitType(StepLimitType());
                                                   >> 227   msc->SetLateralDisplasmentFlag(LateralDisplasmentFlag());
                                                   >> 228   msc->SetSkin(Skin());
                                                   >> 229   msc->SetRangeFactor(RangeFactor());
                                                   >> 230   msc->SetGeomFactor(GeomFactor());
                                                   >> 231       }
                                                   >> 232       msc->SetPolarAngleLimit(polarAngleLimit);
                                                   >> 233       if(msc->HighEnergyLimit() > maxKinEnergy) {
                                                   >> 234   msc->SetHighEnergyLimit(maxKinEnergy);
                                                   >> 235       }
                                                   >> 236     }
                                                   >> 237 
                                                   >> 238     modelManager->Initialise(firstParticle, G4Electron::Electron(), 
                                                   >> 239            10.0, verboseLevel);
                                                   >> 240 
                                                   >> 241     // prepare tables
                                                   >> 242     G4LossTableBuilder* bld = man->GetTableBuilder();
                                                   >> 243     if(buildLambdaTable) {
                                                   >> 244       theLambdaTable = G4PhysicsTableHelper::PreparePhysicsTable(theLambdaTable);
                                                   >> 245       bld->InitialiseBaseMaterials(theLambdaTable);
                                                   >> 246     }
                                                   >> 247     theDensityFactor = bld->GetDensityFactors();
                                                   >> 248     theDensityIdx = bld->GetCoupleIndexes();
                                                   >> 249   }
183 }                                                 250 }
184                                                   251 
185 //....oooOO0OOooo........oooOO0OOooo........oo    252 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
186                                                   253 
187 void G4VMultipleScattering::StreamInfo(std::os << 254 void G4VMultipleScattering::BuildPhysicsTable(const G4ParticleDefinition& part)
188                   const G4ParticleDefinition&  << 
189 {                                                 255 {
190   G4String indent = (rst ? "  " : "");         << 256   G4String num = part.GetParticleName();
191   outFile << G4endl << indent << GetProcessNam << 257   if(1 < verboseLevel) {
192   if (!rst) outFile << " for " << part.GetPart << 258     G4cout << "### G4VMultipleScattering::BuildPhysicsTable() for "
193   outFile  << "  SubType= " << GetProcessSubTy << 259            << GetProcessName()
194   modelManager->DumpModelList(outFile, verbose << 260            << " and particle " << num
                                                   >> 261            << G4endl;
                                                   >> 262   }
                                                   >> 263 
                                                   >> 264   (G4LossTableManager::Instance())->BuildPhysicsTable(firstParticle);
                                                   >> 265 
                                                   >> 266   if (buildLambdaTable && firstParticle == &part) {
                                                   >> 267 
                                                   >> 268     const G4ProductionCutsTable* theCoupleTable=
                                                   >> 269           G4ProductionCutsTable::GetProductionCutsTable();
                                                   >> 270     size_t numOfCouples = theCoupleTable->GetTableSize();
                                                   >> 271 
                                                   >> 272     //G4LossTableBuilder* bld = (G4LossTableManager::Instance())->GetTableBuilder();
                                                   >> 273     G4bool splineFlag = (G4LossTableManager::Instance())->SplineFlag();
                                                   >> 274 
                                                   >> 275     G4PhysicsLogVector* aVector = 0;
                                                   >> 276     G4PhysicsLogVector* bVector = 0;
                                                   >> 277 
                                                   >> 278     for (size_t i=0; i<numOfCouples; ++i) {
                                                   >> 279 
                                                   >> 280       if (theLambdaTable->GetFlag(i)/*bld->GetFlag(i)*/) {
                                                   >> 281         // create physics vector and fill it
                                                   >> 282         const G4MaterialCutsCouple* couple = 
                                                   >> 283     theCoupleTable->GetMaterialCutsCouple(i);
                                                   >> 284   if(!bVector) {
                                                   >> 285     aVector = static_cast<G4PhysicsLogVector*>(PhysicsVector(couple));
                                                   >> 286     bVector = aVector;
                                                   >> 287   } else {
                                                   >> 288     aVector = new G4PhysicsLogVector(*bVector);
                                                   >> 289   }       
                                                   >> 290   aVector->SetSpline(splineFlag);
                                                   >> 291         modelManager->FillLambdaVector(aVector, couple, false);
                                                   >> 292   if(splineFlag) { aVector->FillSecondDerivatives(); }
                                                   >> 293         G4PhysicsTableHelper::SetPhysicsVector(theLambdaTable, i, aVector);
                                                   >> 294       }
                                                   >> 295     }
                                                   >> 296 
                                                   >> 297     if(1 < verboseLevel) {
                                                   >> 298       G4cout << "Lambda table is built for "
                                                   >> 299              << num
                                                   >> 300              << G4endl;
                                                   >> 301     }
                                                   >> 302   }
                                                   >> 303   if(verboseLevel>0 && ( num == "e-" || num == "mu+" ||  
                                                   >> 304                          num == "proton" || num == "pi-" || 
                                                   >> 305        num == "kaon+" || num == "GenericIon")) {
                                                   >> 306     PrintInfoDefinition();
                                                   >> 307     if(2 < verboseLevel && theLambdaTable) 
                                                   >> 308       { G4cout << *theLambdaTable << G4endl; }
                                                   >> 309   }
                                                   >> 310 
                                                   >> 311   if(1 < verboseLevel) {
                                                   >> 312     G4cout << "### G4VMultipleScattering::BuildPhysicsTable() done for "
                                                   >> 313            << GetProcessName()
                                                   >> 314            << " and particle " << num
                                                   >> 315            << G4endl;
                                                   >> 316   }
195 }                                                 317 }
196                                                   318 
197 //....oooOO0OOooo........oooOO0OOooo........oo    319 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
198                                                   320 
199 void G4VMultipleScattering::StartTracking(G4Tr << 321 void G4VMultipleScattering::PrintInfoDefinition()
200 {                                              << 322 {
201   G4VEnergyLossProcess* eloss = nullptr;       << 323   if (0 < verboseLevel) {
202   if(track->GetParticleDefinition() != currPar << 324     G4cout << G4endl << GetProcessName() 
203     currParticle = track->GetParticleDefinitio << 325      << ":   for " << firstParticle->GetParticleName()
204     fIonisation = emManager->GetEnergyLossProc << 326      << "    SubType= " << GetProcessSubType() 
205     eloss = fIonisation;                       << 327      << G4endl;
206   }                                            << 328     if (theLambdaTable) {
207   for(G4int i=0; i<numberOfModels; ++i) {      << 329       G4cout << "      Lambda tables from "
208     G4VMscModel* msc = GetModelByIndex(i);     << 330        << G4BestUnit(MinKinEnergy(),"Energy")
209     msc->StartTracking(track);                 << 331        << " to "
210     if(nullptr != eloss) {                     << 332        << G4BestUnit(MaxKinEnergy(),"Energy")
211       msc->SetIonisation(eloss, currParticle); << 333        << " in " << nBins << " bins, spline: " 
                                                   >> 334        << (G4LossTableManager::Instance())->SplineFlag()
                                                   >> 335        << G4endl;
                                                   >> 336     }
                                                   >> 337     PrintInfo();
                                                   >> 338     modelManager->DumpModelList(verboseLevel);
                                                   >> 339     if (2 < verboseLevel) {
                                                   >> 340       G4cout << "LambdaTable address= " << theLambdaTable << G4endl;
                                                   >> 341       if(theLambdaTable) { G4cout << (*theLambdaTable) << G4endl; }
212     }                                             342     }
213   }                                               343   }
214 }                                                 344 }
215                                                   345 
216 //....oooOO0OOooo........oooOO0OOooo........oo    346 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
217                                                   347 
218 G4double G4VMultipleScattering::AlongStepGetPh    348 G4double G4VMultipleScattering::AlongStepGetPhysicalInteractionLength(
219                              const G4Track& tr    349                              const G4Track& track,
220                              G4double,            350                              G4double,
221                              G4double currentM    351                              G4double currentMinimalStep,
222                              G4double&,           352                              G4double&,
223                              G4GPILSelection*     353                              G4GPILSelection* selection)
224 {                                                 354 {
225   // get Step limit proposed by the process       355   // get Step limit proposed by the process
226   *selection = NotCandidateForSelection;          356   *selection = NotCandidateForSelection;
227   physStepLimit = gPathLength = tPathLength =  << 357   G4double x = currentMinimalStep;
228                                                << 358   DefineMaterial(track.GetMaterialCutsCouple());
229   G4double ekin = track.GetKineticEnergy();       359   G4double ekin = track.GetKineticEnergy();
230   /*                                           << 360   if(isIon) { ekin *= proton_mass_c2/track.GetParticleDefinition()->GetPDGMass(); }
231   G4cout << "MSC::AlongStepGPIL: Ekin= " << ek << 361   currentModel = static_cast<G4VMscModel*>(SelectModel(ekin));
232          << "  " << currParticle->GetParticleN << 362 
233          << " currMod " << currentModel        << 363   // define ionisation process
234          << G4endl;                            << 364   if(!currentModel->GetIonisation()) {
235   */                                           << 365     currentModel->SetIonisation(G4LossTableManager::Instance()
236   // isIon flag is used only to select a model << 366         ->GetEnergyLossProcess(track.GetParticleDefinition()));
237   if(isIon) {                                  << 367   }  
238     ekin *= proton_mass_c2/track.GetParticleDe << 368 
239   }                                            << 369   if(x > 0.0 && ekin > 0.0 && currentModel->IsActive(ekin)) {
240   const G4MaterialCutsCouple* couple = track.G << 370     G4double tPathLength = 
241                                                << 371       currentModel->ComputeTruePathLengthLimit(track, theLambdaTable, x);
242   // select new model, static cast is possible << 372     if (tPathLength < x) { *selection = CandidateForSelection; }
243   if(1 < numberOfModels) {                     << 373     x = currentModel->ComputeGeomPathLength(tPathLength);
244     currentModel =                             << 374     //  G4cout << "tPathLength= " << tPathLength
245       static_cast<G4VMscModel*>(SelectModel(ek << 375     //         << " stepLimit= " << x
246   }                                            << 376     //        << " currentMinimalStep= " << currentMinimalStep<< G4endl;
247   currentModel->SetCurrentCouple(couple);      << 377   }
248   // msc is active is model is active, energy  << 378   return x;
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   }                                            << 
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   */                                           << 
273   return gPathLength;                          << 
274 }                                                 379 }
275                                                   380 
276 //....oooOO0OOooo........oooOO0OOooo........oo    381 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
277                                                   382 
278 G4double                                          383 G4double 
279 G4VMultipleScattering::PostStepGetPhysicalInte    384 G4VMultipleScattering::PostStepGetPhysicalInteractionLength(
280               const G4Track&, G4double, G4Forc    385               const G4Track&, G4double, G4ForceCondition* condition)
281 {                                                 386 {
282   *condition = NotForced;                      << 387   *condition = Forced;
283   return DBL_MAX;                                 388   return DBL_MAX;
284 }                                                 389 }
285                                                   390 
286 //....oooOO0OOooo........oooOO0OOooo........oo    391 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
287                                                   392 
288 G4VParticleChange*                                393 G4VParticleChange* 
289 G4VMultipleScattering::AlongStepDoIt(const G4T    394 G4VMultipleScattering::AlongStepDoIt(const G4Track& track, const G4Step& step)
290 {                                                 395 {
291   fParticleChange.InitialiseMSC(track, step);  << 396   if(currentModel->IsActive(track.GetKineticEnergy())) {
292   fNewPosition = fParticleChange.GetProposedPo << 397     fParticleChange.ProposeTrueStepLength(currentModel->ComputeTrueStepLength(step.GetStepLength()));
293   fPositionChanged = false;                    << 398   } else {
294                                                << 399     fParticleChange.ProposeTrueStepLength(step.GetStepLength());
295   G4double geomLength = step.GetStepLength();  << 400   }
                                                   >> 401   return &fParticleChange;
                                                   >> 402 }
296                                                   403 
297   // very small step - no msc                  << 404 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
298   if(!isActive) {                              << 
299     tPathLength = geomLength;                  << 
300                                                   405 
301     // sample msc                              << 406 G4VParticleChange* 
302   } else {                                     << 407 G4VMultipleScattering::PostStepDoIt(const G4Track& track, const G4Step& step)
303     G4double range =                           << 408 {
304       currentModel->GetRange(currParticle,trac << 409   fParticleChange.Initialize(track);
305                              track.GetMaterial << 410   if(currentModel->IsActive(track.GetKineticEnergy())) {
306                                                << 411     currentModel->SampleScattering(track.GetDynamicParticle(),
307     tPathLength = currentModel->ComputeTrueSte << 412            step.GetPostStepPoint()->GetSafety());
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       }                                        << 
366     }                                          << 
367   }                                               413   }
368   fParticleChange.ProposeTrueStepLength(tPathL << 
369   return &fParticleChange;                        414   return &fParticleChange;
370 }                                                 415 }
371                                                   416 
372 //....oooOO0OOooo........oooOO0OOooo........oo    417 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
373                                                   418 
374 G4double G4VMultipleScattering::GetContinuousS    419 G4double G4VMultipleScattering::GetContinuousStepLimit(
375                                        const G    420                                        const G4Track& track,
376                                        G4doubl    421                                        G4double previousStepSize,
377                                        G4doubl    422                                        G4double currentMinimalStep,
378                                        G4doubl    423                                        G4double& currentSafety)
379 {                                                 424 {
380   G4GPILSelection selection = NotCandidateForS    425   G4GPILSelection selection = NotCandidateForSelection;
381   G4double x = AlongStepGetPhysicalInteraction    426   G4double x = AlongStepGetPhysicalInteractionLength(track,previousStepSize,
382                                                << 427                  currentMinimalStep,
383                                                << 428                  currentSafety, &selection);
384                                                << 
385   return x;                                       429   return x;
386 }                                                 430 }
387                                                   431 
388 //....oooOO0OOooo........oooOO0OOooo........oo    432 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
389                                                   433 
390 G4double G4VMultipleScattering::ContinuousStep    434 G4double G4VMultipleScattering::ContinuousStepLimit(
391                                        const G    435                                        const G4Track& track,
392                                        G4doubl    436                                        G4double previousStepSize,
393                                        G4doubl    437                                        G4double currentMinimalStep,
394                                        G4doubl    438                                        G4double& currentSafety)
395 {                                                 439 {
396   return GetContinuousStepLimit(track,previous    440   return GetContinuousStepLimit(track,previousStepSize,currentMinimalStep,
397                                 currentSafety) << 441         currentSafety);
398 }                                                 442 }
399                                                   443 
400 //....oooOO0OOooo........oooOO0OOooo........oo    444 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
401                                                   445 
402 G4double G4VMultipleScattering::GetMeanFreePat    446 G4double G4VMultipleScattering::GetMeanFreePath(
403               const G4Track&, G4double, G4Forc    447               const G4Track&, G4double, G4ForceCondition* condition)
404 {                                                 448 {
405   *condition = Forced;                            449   *condition = Forced;
406   return DBL_MAX;                                 450   return DBL_MAX;
407 }                                                 451 }
408                                                   452 
                                                   >> 453 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
                                                   >> 454 
                                                   >> 455 G4PhysicsVector* 
                                                   >> 456 G4VMultipleScattering::PhysicsVector(const G4MaterialCutsCouple* couple)
                                                   >> 457 {
                                                   >> 458   G4int nbins = 3;
                                                   >> 459   if( couple->IsUsed() ) { nbins = nBins; }
                                                   >> 460   G4PhysicsVector* v = new G4PhysicsLogVector(minKinEnergy, maxKinEnergy, nbins);
                                                   >> 461   v->SetSpline((G4LossTableManager::Instance())->SplineFlag());
                                                   >> 462   return v;
                                                   >> 463 }
                                                   >> 464 
409 //....oooOO0OOooo........oooOO0OOooo........oo    465 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
410                                                   466 
411 G4bool                                         << 467 G4bool G4VMultipleScattering::StorePhysicsTable(const G4ParticleDefinition* part,
412 G4VMultipleScattering::StorePhysicsTable(const << 468                         const G4String& directory,
413                                          const << 469                               G4bool ascii)
414                                          G4boo << 
415 {                                                 470 {
416   G4bool yes = true;                              471   G4bool yes = true;
417   if(part != firstParticle || !emManager->IsMa << 472   if ( theLambdaTable && part == firstParticle) {
418                                                << 473     const G4String name = GetPhysicsTableFileName(part,directory,"Lambda",ascii);
419   return G4EmTableUtil::StoreMscTable(this, pa << 474     G4bool yes = theLambdaTable->StorePhysicsTable(name,ascii);
420               numberOfModels, verboseLevel,    << 475 
421                                       ascii);  << 476     if ( yes ) {
                                                   >> 477       if ( verboseLevel>0 ) {
                                                   >> 478         G4cout << "Physics table are stored for " << part->GetParticleName()
                                                   >> 479                << " and process " << GetProcessName()
                                                   >> 480          << " in the directory <" << directory
                                                   >> 481          << "> " << G4endl;
                                                   >> 482       }
                                                   >> 483     } else {
                                                   >> 484       G4cout << "Fail to store Physics Table for " << part->GetParticleName()
                                                   >> 485              << " and process " << GetProcessName()
                                                   >> 486        << " in the directory <" << directory
                                                   >> 487        << "> " << G4endl;
                                                   >> 488     }
                                                   >> 489   }
                                                   >> 490   return yes;
422 }                                                 491 }
423                                                   492 
424 //....oooOO0OOooo........oooOO0OOooo........oo    493 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
425                                                   494 
426 G4bool                                            495 G4bool 
427 G4VMultipleScattering::RetrievePhysicsTable(co << 496 G4VMultipleScattering::RetrievePhysicsTable(const G4ParticleDefinition* part,
428                                             co << 497               const G4String& directory,
429                                             G4 << 498               G4bool ascii)
430 {                                              << 499 {
431   return true;                                 << 500   if(0 < verboseLevel) {
432 }                                              << 501     G4cout << "G4VMultipleScattering::RetrievePhysicsTable() for "
                                                   >> 502            << part->GetParticleName() << " and process "
                                                   >> 503      << GetProcessName() << G4endl;
                                                   >> 504   }
                                                   >> 505   G4bool yes = true;
433                                                   506 
434 //....oooOO0OOooo........oooOO0OOooo........oo << 507   if(!buildLambdaTable || firstParticle != part) return yes;
435                                                   508 
436 void G4VMultipleScattering::ProcessDescription << 509   const G4String particleName = part->GetParticleName();
437 {                                              << 510 
438   if(nullptr != firstParticle) {               << 511   G4String filename = GetPhysicsTableFileName(part,directory,"Lambda",ascii);
439     StreamInfo(outFile, *firstParticle, true); << 512   yes = 
                                                   >> 513     G4PhysicsTableHelper::RetrievePhysicsTable(theLambdaTable,filename,ascii);
                                                   >> 514   if ( yes ) {
                                                   >> 515     if (0 < verboseLevel) {
                                                   >> 516         G4cout << "Lambda table for " << part->GetParticleName() 
                                                   >> 517          << " is retrieved from <"
                                                   >> 518                << filename << ">"
                                                   >> 519                << G4endl;
                                                   >> 520     }
                                                   >> 521     if((G4LossTableManager::Instance())->SplineFlag()) {
                                                   >> 522       size_t n = theLambdaTable->length();
                                                   >> 523       for(size_t i=0; i<n; ++i) {
                                                   >> 524         if((* theLambdaTable)[i]) {
                                                   >> 525     (* theLambdaTable)[i]->SetSpline(true);
                                                   >> 526   }
                                                   >> 527       }
                                                   >> 528     }
                                                   >> 529   } else {
                                                   >> 530     if (1 < verboseLevel) {
                                                   >> 531         G4cout << "Lambda table for " << part->GetParticleName() 
                                                   >> 532          << " in file <"
                                                   >> 533                << filename << "> is not exist"
                                                   >> 534                << G4endl;
                                                   >> 535     }
440   }                                               536   }
                                                   >> 537   return yes;
441 }                                                 538 }
442                                                   539 
443 //....oooOO0OOooo........oooOO0OOooo........oo    540 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
444                                                   541 
445                                                   542