Geant4 Cross Reference |
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 // 26 // 27 // ------------------------------------------- 27 // ------------------------------------------------------------------- 28 // 28 // 29 // GEANT4 Class file 29 // GEANT4 Class file 30 // 30 // 31 // 31 // 32 // File name: G4VMultipleScattering 32 // File name: G4VMultipleScattering 33 // 33 // 34 // Author: Vladimir Ivanchenko on base 34 // Author: Vladimir Ivanchenko on base of Laszlo Urban code 35 // 35 // 36 // Creation date: 25.03.2003 36 // Creation date: 25.03.2003 37 // 37 // 38 // Modifications: 38 // Modifications: 39 // 39 // >> 40 // 13.04.03 Change printout (V.Ivanchenko) >> 41 // 04-06-03 Fix compilation warnings (V.Ivanchenko) 40 // 16-07-03 Use G4VMscModel interface (V.Ivanc 42 // 16-07-03 Use G4VMscModel interface (V.Ivanchenko) 41 // 03-11-03 Fix initialisation problem in Retr 43 // 03-11-03 Fix initialisation problem in RetrievePhysicsTable (V.Ivanchenko) 42 // 04-11-03 Update PrintInfoDefinition (V.Ivan 44 // 04-11-03 Update PrintInfoDefinition (V.Ivanchenko) 43 // 01-03-04 SampleCosineTheta signature change 45 // 01-03-04 SampleCosineTheta signature changed 44 // 22-04-04 SampleCosineTheta signature change 46 // 22-04-04 SampleCosineTheta signature changed back to original 45 // 27-08-04 Add InitialiseForRun method (V.Iva 47 // 27-08-04 Add InitialiseForRun method (V.Ivanchneko) 46 // 08-11-04 Migration to new interface of Stor 48 // 08-11-04 Migration to new interface of Store/Retrieve tables (V.Ivantchenko) 47 // 11-03-05 Shift verbose level by 1 (V.Ivantc 49 // 11-03-05 Shift verbose level by 1 (V.Ivantchenko) 48 // 15-04-05 optimize internal interface (V.Iva 50 // 15-04-05 optimize internal interface (V.Ivanchenko) 49 // 15-04-05 remove boundary flag (V.Ivanchenko 51 // 15-04-05 remove boundary flag (V.Ivanchenko) 50 // 27-10-05 introduce virtual function MscStep 52 // 27-10-05 introduce virtual function MscStepLimitation() (V.Ivanchenko) 51 // 12-04-07 Add verbosity at destruction (V.Iv 53 // 12-04-07 Add verbosity at destruction (V.Ivanchenko) 52 // 27-10-07 Virtual functions moved to source 54 // 27-10-07 Virtual functions moved to source (V.Ivanchenko) 53 // 11-03-08 Set skin value does not effect ste 55 // 11-03-08 Set skin value does not effect step limit type (V.Ivanchenko) 54 // 24-06-09 Removed hidden bin in G4PhysicsVec 56 // 24-06-09 Removed hidden bin in G4PhysicsVector (V.Ivanchenko) 55 // 04-06-13 Adoptation to MT mode (V.Ivanchenk 57 // 04-06-13 Adoptation to MT mode (V.Ivanchenko) 56 // 58 // >> 59 // Class Description: >> 60 // >> 61 // It is the generic process of multiple scattering it includes common >> 62 // part of calculations for all charged particles 57 63 58 // ------------------------------------------- 64 // ------------------------------------------------------------------- 59 // 65 // 60 //....oooOO0OOooo........oooOO0OOooo........oo 66 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 61 //....oooOO0OOooo........oooOO0OOooo........oo 67 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 62 68 63 #include "G4VMultipleScattering.hh" 69 #include "G4VMultipleScattering.hh" 64 #include "G4PhysicalConstants.hh" 70 #include "G4PhysicalConstants.hh" 65 #include "G4SystemOfUnits.hh" 71 #include "G4SystemOfUnits.hh" 66 #include "G4LossTableManager.hh" 72 #include "G4LossTableManager.hh" 67 #include "G4MaterialCutsCouple.hh" 73 #include "G4MaterialCutsCouple.hh" 68 #include "G4Step.hh" 74 #include "G4Step.hh" 69 #include "G4ParticleDefinition.hh" 75 #include "G4ParticleDefinition.hh" 70 #include "G4VEmFluctuationModel.hh" 76 #include "G4VEmFluctuationModel.hh" 71 #include "G4UnitsTable.hh" 77 #include "G4UnitsTable.hh" 72 #include "G4ProductionCutsTable.hh" 78 #include "G4ProductionCutsTable.hh" 73 #include "G4Electron.hh" 79 #include "G4Electron.hh" 74 #include "G4GenericIon.hh" 80 #include "G4GenericIon.hh" 75 #include "G4TransportationManager.hh" 81 #include "G4TransportationManager.hh" 76 #include "G4SafetyHelper.hh" 82 #include "G4SafetyHelper.hh" 77 #include "G4ParticleTable.hh" 83 #include "G4ParticleTable.hh" 78 #include "G4ProcessVector.hh" 84 #include "G4ProcessVector.hh" 79 #include "G4ProcessManager.hh" 85 #include "G4ProcessManager.hh" 80 #include "G4LossTableBuilder.hh" << 81 #include "G4EmTableUtil.hh" << 82 #include <iostream> 86 #include <iostream> 83 87 84 //....oooOO0OOooo........oooOO0OOooo........oo 88 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 85 89 86 G4VMultipleScattering::G4VMultipleScattering(c << 90 G4VMultipleScattering::G4VMultipleScattering(const G4String& name, G4ProcessType) 87 : G4VContinuousDiscreteProcess("msc", fElect 91 : G4VContinuousDiscreteProcess("msc", fElectromagnetic), >> 92 numberOfModels(0), >> 93 firstParticle(nullptr), >> 94 currParticle(nullptr), >> 95 stepLimit(fUseSafety), >> 96 facrange(0.04), >> 97 latDisplacement(true), >> 98 isIon(false), 88 fNewPosition(0.,0.,0.), 99 fNewPosition(0.,0.,0.), 89 fNewDirection(0.,0.,1.) 100 fNewDirection(0.,0.,1.) 90 { 101 { 91 theParameters = G4EmParameters::Instance(); 102 theParameters = G4EmParameters::Instance(); 92 SetVerboseLevel(1); 103 SetVerboseLevel(1); 93 SetProcessSubType(fMultipleScattering); 104 SetProcessSubType(fMultipleScattering); >> 105 if("ionmsc" == name) { firstParticle = G4GenericIon::GenericIon(); } 94 106 95 lowestKinEnergy = 10*CLHEP::eV; 107 lowestKinEnergy = 10*CLHEP::eV; 96 108 97 geomMin = 0.05*CLHEP::nm; << 109 physStepLimit = gPathLength = tPathLength = 0.0; >> 110 fIonisation = nullptr; >> 111 >> 112 geomMin = 0.05*CLHEP::nm; 98 minDisplacement2 = geomMin*geomMin; 113 minDisplacement2 = geomMin*geomMin; 99 114 100 pParticleChange = &fParticleChange; 115 pParticleChange = &fParticleChange; >> 116 safetyHelper = nullptr; >> 117 fPositionChanged = false; >> 118 isActive = false; 101 119 >> 120 currentModel = nullptr; 102 modelManager = new G4EmModelManager(); 121 modelManager = new G4EmModelManager(); 103 emManager = G4LossTableManager::Instance(); 122 emManager = G4LossTableManager::Instance(); 104 mscModels.reserve(2); 123 mscModels.reserve(2); 105 emManager->Register(this); 124 emManager->Register(this); 106 } 125 } 107 126 108 //....oooOO0OOooo........oooOO0OOooo........oo 127 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 109 128 110 G4VMultipleScattering::~G4VMultipleScattering( 129 G4VMultipleScattering::~G4VMultipleScattering() 111 { 130 { >> 131 /* >> 132 if(1 < verboseLevel) { >> 133 G4cout << "G4VMultipleScattering destruct " << GetProcessName() >> 134 << G4endl; >> 135 } >> 136 */ 112 delete modelManager; 137 delete modelManager; 113 emManager->DeRegister(this); 138 emManager->DeRegister(this); 114 } 139 } 115 140 116 //....oooOO0OOooo........oooOO0OOooo........oo 141 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 117 142 118 void G4VMultipleScattering::AddEmModel(G4int o << 143 void G4VMultipleScattering::AddEmModel(G4int order, G4VEmModel* p, 119 const G 144 const G4Region* region) 120 { 145 { 121 if(nullptr == ptr) { return; } << 146 if(!p) { return; } 122 G4VEmFluctuationModel* fm = nullptr; 147 G4VEmFluctuationModel* fm = nullptr; 123 modelManager->AddEmModel(order, ptr, fm, reg << 148 modelManager->AddEmModel(order, p, fm, region); 124 ptr->SetParticleChange(pParticleChange); << 149 p->SetParticleChange(pParticleChange); 125 } 150 } 126 151 127 //....oooOO0OOooo........oooOO0OOooo........oo 152 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 128 153 129 void G4VMultipleScattering::SetEmModel(G4VMscM << 154 void G4VMultipleScattering::SetEmModel(G4VMscModel* p, size_t) 130 { 155 { 131 if(nullptr == ptr) { return; } << 156 for(auto & msc : mscModels) { if(msc == p) { return; } } 132 if(!mscModels.empty()) { << 157 mscModels.push_back(p); 133 for(auto & msc : mscModels) { if(msc == pt << 158 } 134 } << 159 135 mscModels.push_back(ptr); << 160 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... >> 161 >> 162 G4VMscModel* G4VMultipleScattering::EmModel(size_t index) const >> 163 { >> 164 return (index < mscModels.size()) ? mscModels[index] : nullptr; 136 } 165 } 137 166 138 //....oooOO0OOooo........oooOO0OOooo........oo 167 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 139 168 140 void 169 void 141 G4VMultipleScattering::PreparePhysicsTable(con 170 G4VMultipleScattering::PreparePhysicsTable(const G4ParticleDefinition& part) 142 { 171 { >> 172 if(1 < verboseLevel) { >> 173 G4cout << "### G4VMultipleScattering::PrepearPhysicsTable() for " >> 174 << GetProcessName() >> 175 << " and particle " << part.GetParticleName() >> 176 << G4endl; >> 177 } 143 G4bool master = emManager->IsMaster(); 178 G4bool master = emManager->IsMaster(); 144 if (nullptr == firstParticle) { firstParticl << 145 179 146 emManager->PreparePhysicsTable(&part, this); << 180 if(!firstParticle) { firstParticle = ∂ } >> 181 if(part.GetParticleType() == "nucleus") { >> 182 stepLimit = fMinimal; >> 183 latDisplacement = false; >> 184 facrange = 0.2; >> 185 G4String pname = part.GetParticleName(); >> 186 if(pname != "deuteron" && pname != "triton" && >> 187 pname != "alpha+" && pname != "helium" && >> 188 pname != "alpha" && pname != "He3" && >> 189 pname != "hydrogen") { >> 190 >> 191 const G4ParticleDefinition* theGenericIon = >> 192 G4ParticleTable::GetParticleTable()->FindParticle("GenericIon"); >> 193 if(&part == theGenericIon) { isIon = true; } >> 194 >> 195 if(theGenericIon && firstParticle != theGenericIon) { >> 196 G4ProcessManager* pm = theGenericIon->GetProcessManager(); >> 197 G4ProcessVector* v = pm->GetAlongStepProcessVector(); >> 198 size_t n = v->size(); >> 199 for(size_t j=0; j<n; ++j) { >> 200 if((*v)[j] == this) { >> 201 firstParticle = theGenericIon; >> 202 isIon = true; >> 203 break; >> 204 } >> 205 } >> 206 } >> 207 } >> 208 } >> 209 >> 210 emManager->PreparePhysicsTable(&part, this, master); 147 currParticle = nullptr; 211 currParticle = nullptr; 148 212 >> 213 if(1 < verboseLevel) { >> 214 G4cout << "### G4VMultipleScattering::PrepearPhysicsTable() for " >> 215 << GetProcessName() >> 216 << " and particle " << part.GetParticleName() >> 217 << " local particle " << firstParticle->GetParticleName() >> 218 << " isIon: " << isIon << " isMaster: " << master >> 219 << G4endl; >> 220 } >> 221 149 if(firstParticle == &part) { 222 if(firstParticle == &part) { 150 baseMat = emManager->GetTableBuilder()->Ge << 151 G4EmTableUtil::PrepareMscProcess(this, par << 152 stepLimit, facrange, << 153 latDisplacement, master, << 154 isIon, baseMat); << 155 223 >> 224 // initialise process >> 225 InitialiseProcess(firstParticle); >> 226 >> 227 // heavy particles and not ions >> 228 if(!isIon) { >> 229 if(part.GetPDGMass() > MeV) { >> 230 stepLimit = theParameters->MscMuHadStepLimitType(); >> 231 facrange = theParameters->MscMuHadRangeFactor(); >> 232 latDisplacement = theParameters->MuHadLateralDisplacement(); >> 233 } else { >> 234 stepLimit = theParameters->MscStepLimitType(); >> 235 facrange = theParameters->MscRangeFactor(); >> 236 latDisplacement = theParameters->LateralDisplacement(); >> 237 } >> 238 } >> 239 if(master) { SetVerboseLevel(theParameters->Verbose()); } >> 240 else { SetVerboseLevel(theParameters->WorkerVerbose()); } >> 241 >> 242 // initialisation of models 156 numberOfModels = modelManager->NumberOfMod 243 numberOfModels = modelManager->NumberOfModels(); 157 currentModel = GetModelByIndex(0); << 244 /* >> 245 G4cout << "### G4VMultipleScattering::PreparePhysicsTable() for " >> 246 << GetProcessName() >> 247 << " and particle " << part.GetParticleName() >> 248 << " Nmod= " << numberOfModels << " " << this >> 249 << G4endl; >> 250 */ >> 251 for(G4int i=0; i<numberOfModels; ++i) { >> 252 G4VMscModel* msc = static_cast<G4VMscModel*>(GetModelByIndex(i)); >> 253 if(!msc) { continue; } >> 254 msc->SetIonisation(nullptr, firstParticle); >> 255 msc->SetMasterThread(master); >> 256 currentModel = msc; >> 257 msc->SetPolarAngleLimit(theParameters->MscThetaLimit()); >> 258 G4double emax = >> 259 std::min(msc->HighEnergyLimit(),theParameters->MaxKinEnergy()); >> 260 msc->SetHighEnergyLimit(emax); >> 261 } 158 262 159 if (nullptr == safetyHelper) { << 263 modelManager->Initialise(firstParticle, G4Electron::Electron(), >> 264 10.0, verboseLevel); >> 265 >> 266 if(!safetyHelper) { 160 safetyHelper = G4TransportationManager:: 267 safetyHelper = G4TransportationManager::GetTransportationManager() 161 ->GetSafetyHelper(); << 268 ->GetSafetyHelper(); 162 safetyHelper->InitialiseHelper(); 269 safetyHelper->InitialiseHelper(); 163 } 270 } 164 } 271 } 165 } 272 } 166 273 167 //....oooOO0OOooo........oooOO0OOooo........oo 274 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 168 275 169 void G4VMultipleScattering::BuildPhysicsTable( 276 void G4VMultipleScattering::BuildPhysicsTable(const G4ParticleDefinition& part) 170 { 277 { >> 278 G4String num = part.GetParticleName(); 171 G4bool master = emManager->IsMaster(); 279 G4bool master = emManager->IsMaster(); 172 << 280 if(1 < verboseLevel) { 173 if(firstParticle == &part) { << 281 G4cout << "### G4VMultipleScattering::BuildPhysicsTable() for " 174 emManager->BuildPhysicsTable(&part); << 282 << GetProcessName() >> 283 << " and particle " << num << " isIon: " << isIon >> 284 << " IsMaster: " << master << G4endl; 175 } 285 } 176 const G4VMultipleScattering* ptr = this; << 286 const G4VMultipleScattering* masterProcess = 177 if(!master) { << 287 static_cast<const G4VMultipleScattering*>(GetMasterProcess()); 178 ptr = static_cast<const G4VMultipleScatter << 288 >> 289 if(firstParticle == &part) { >> 290 /* >> 291 G4cout << "### G4VMultipleScattering::BuildPhysicsTable() for " >> 292 << GetProcessName() >> 293 << " and particle " << num >> 294 << " IsMaster= " << G4LossTableManager::Instance()->IsMaster() >> 295 << " " << this >> 296 << G4endl; >> 297 */ >> 298 emManager->BuildPhysicsTable(firstParticle); >> 299 >> 300 if(!master) { >> 301 // initialisation of models >> 302 /* >> 303 G4cout << "### G4VMultipleScattering::BuildPhysicsTable() for " >> 304 << GetProcessName() >> 305 << " and particle " << num >> 306 << " Nmod= " << numberOfModels << " " << this >> 307 << G4endl; >> 308 */ >> 309 for(G4int i=0; i<numberOfModels; ++i) { >> 310 G4VMscModel* msc = static_cast<G4VMscModel*>(GetModelByIndex(i)); >> 311 if(!msc) { continue; } >> 312 G4VMscModel* msc0= >> 313 static_cast<G4VMscModel*>(masterProcess->GetModelByIndex(i)); >> 314 msc->SetCrossSectionTable(msc0->GetCrossSectionTable(), false); >> 315 msc->InitialiseLocal(firstParticle, msc0); >> 316 } >> 317 } 179 } 318 } 180 319 181 G4EmTableUtil::BuildMscProcess(this, ptr, pa << 320 // explicitly defined printout by particle name 182 numberOfModels, master); << 321 if(1 < verboseLevel || >> 322 (0 < verboseLevel && (num == "e-" || >> 323 num == "e+" || num == "mu+" || >> 324 num == "mu-" || num == "proton"|| >> 325 num == "pi+" || num == "pi-" || >> 326 num == "kaon+" || num == "kaon-" || >> 327 num == "alpha" || num == "anti_proton" || >> 328 num == "GenericIon" || num == "alpha+" || >> 329 num == "alpha++" ))) >> 330 { >> 331 StreamInfo(G4cout, part); >> 332 } >> 333 >> 334 if(1 < verboseLevel) { >> 335 G4cout << "### G4VMultipleScattering::BuildPhysicsTable() done for " >> 336 << GetProcessName() >> 337 << " and particle " << num >> 338 << G4endl; >> 339 } 183 } 340 } 184 341 185 //....oooOO0OOooo........oooOO0OOooo........oo 342 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 186 343 187 void G4VMultipleScattering::StreamInfo(std::os 344 void G4VMultipleScattering::StreamInfo(std::ostream& outFile, 188 const G4ParticleDefinition& 345 const G4ParticleDefinition& part, G4bool rst) const 189 { 346 { 190 G4String indent = (rst ? " " : ""); 347 G4String indent = (rst ? " " : ""); 191 outFile << G4endl << indent << GetProcessNam 348 outFile << G4endl << indent << GetProcessName() << ": "; 192 if (!rst) outFile << " for " << part.GetPart 349 if (!rst) outFile << " for " << part.GetParticleName(); 193 outFile << " SubType= " << GetProcessSubTy 350 outFile << " SubType= " << GetProcessSubType() << G4endl; >> 351 //StreamProcessInfo(outFile); 194 modelManager->DumpModelList(outFile, verbose 352 modelManager->DumpModelList(outFile, verboseLevel); 195 } 353 } 196 354 197 //....oooOO0OOooo........oooOO0OOooo........oo 355 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 198 356 199 void G4VMultipleScattering::StartTracking(G4Tr 357 void G4VMultipleScattering::StartTracking(G4Track* track) 200 { 358 { 201 G4VEnergyLossProcess* eloss = nullptr; 359 G4VEnergyLossProcess* eloss = nullptr; 202 if(track->GetParticleDefinition() != currPar 360 if(track->GetParticleDefinition() != currParticle) { 203 currParticle = track->GetParticleDefinitio 361 currParticle = track->GetParticleDefinition(); 204 fIonisation = emManager->GetEnergyLossProc 362 fIonisation = emManager->GetEnergyLossProcess(currParticle); 205 eloss = fIonisation; 363 eloss = fIonisation; 206 } 364 } >> 365 /* >> 366 G4cout << "G4VMultipleScattering::StartTracking Nmod= " << numberOfModels >> 367 << " " << currParticle->GetParticleName() >> 368 << " E(MeV)= " << track->GetKineticEnergy() >> 369 << " Ion= " << eloss << " " << fIonisation << " IsMaster= " >> 370 << G4LossTableManager::Instance()->IsMaster() >> 371 << G4endl; >> 372 */ 207 for(G4int i=0; i<numberOfModels; ++i) { 373 for(G4int i=0; i<numberOfModels; ++i) { 208 G4VMscModel* msc = GetModelByIndex(i); << 374 /* >> 375 G4cout << "Next model " << i << " " << msc >> 376 << " Emin= " << msc->LowEnergyLimit() >> 377 << " Emax= " << msc->HighEnergyLimit() >> 378 << " Eact= " << msc->LowEnergyActivationLimit() << G4endl; >> 379 */ >> 380 G4VEmModel* msc = GetModelByIndex(i); 209 msc->StartTracking(track); 381 msc->StartTracking(track); 210 if(nullptr != eloss) { << 382 if(eloss) { 211 msc->SetIonisation(eloss, currParticle); << 383 G4VMscModel* mscmod = static_cast<G4VMscModel*>(msc); >> 384 if(mscmod) { mscmod->SetIonisation(fIonisation, currParticle); } 212 } 385 } 213 } 386 } 214 } << 387 } 215 388 216 //....oooOO0OOooo........oooOO0OOooo........oo 389 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 217 390 218 G4double G4VMultipleScattering::AlongStepGetPh 391 G4double G4VMultipleScattering::AlongStepGetPhysicalInteractionLength( 219 const G4Track& tr 392 const G4Track& track, 220 G4double, 393 G4double, 221 G4double currentM 394 G4double currentMinimalStep, 222 G4double&, 395 G4double&, 223 G4GPILSelection* 396 G4GPILSelection* selection) 224 { 397 { 225 // get Step limit proposed by the process 398 // get Step limit proposed by the process 226 *selection = NotCandidateForSelection; 399 *selection = NotCandidateForSelection; 227 physStepLimit = gPathLength = tPathLength = 400 physStepLimit = gPathLength = tPathLength = currentMinimalStep; 228 401 229 G4double ekin = track.GetKineticEnergy(); 402 G4double ekin = track.GetKineticEnergy(); 230 /* 403 /* 231 G4cout << "MSC::AlongStepGPIL: Ekin= " << ek 404 G4cout << "MSC::AlongStepGPIL: Ekin= " << ekin 232 << " " << currParticle->GetParticleN 405 << " " << currParticle->GetParticleName() 233 << " currMod " << currentModel 406 << " currMod " << currentModel 234 << G4endl; 407 << G4endl; 235 */ 408 */ 236 // isIon flag is used only to select a model 409 // isIon flag is used only to select a model 237 if(isIon) { 410 if(isIon) { 238 ekin *= proton_mass_c2/track.GetParticleDe 411 ekin *= proton_mass_c2/track.GetParticleDefinition()->GetPDGMass(); 239 } 412 } 240 const G4MaterialCutsCouple* couple = track.G << 241 413 242 // select new model, static cast is possible << 414 // select new model 243 if(1 < numberOfModels) { 415 if(1 < numberOfModels) { 244 currentModel = << 416 currentModel = static_cast<G4VMscModel*>( 245 static_cast<G4VMscModel*>(SelectModel(ek << 417 SelectModel(ekin,track.GetMaterialCutsCouple()->GetIndex())); 246 } 418 } 247 currentModel->SetCurrentCouple(couple); << 248 // msc is active is model is active, energy 419 // msc is active is model is active, energy above the limit, 249 // and step size is above the limit; 420 // and step size is above the limit; 250 // if it is active msc may limit the step 421 // if it is active msc may limit the step 251 if(currentModel->IsActive(ekin) && tPathLeng 422 if(currentModel->IsActive(ekin) && tPathLength > geomMin 252 && ekin >= lowestKinEnergy) { 423 && ekin >= lowestKinEnergy) { 253 isActive = true; 424 isActive = true; 254 tPathLength = 425 tPathLength = 255 currentModel->ComputeTruePathLengthLimit 426 currentModel->ComputeTruePathLengthLimit(track, gPathLength); 256 if (tPathLength < physStepLimit) { 427 if (tPathLength < physStepLimit) { 257 *selection = CandidateForSelection; 428 *selection = CandidateForSelection; 258 } 429 } 259 } else { << 430 } else { isActive = false; } 260 isActive = false; << 261 gPathLength = DBL_MAX; << 262 } << 263 431 264 //if(currParticle->GetPDGMass() > GeV) 432 //if(currParticle->GetPDGMass() > GeV) 265 /* 433 /* 266 G4cout << "MSC::AlongStepGPIL: Ekin= " << ek 434 G4cout << "MSC::AlongStepGPIL: Ekin= " << ekin 267 << " " << currParticle->GetParticleN 435 << " " << currParticle->GetParticleName() 268 << " gPathLength= " << gPathLength 436 << " gPathLength= " << gPathLength 269 << " tPathLength= " << tPathLength 437 << " tPathLength= " << tPathLength 270 << " currentMinimalStep= " << current 438 << " currentMinimalStep= " << currentMinimalStep 271 << " isActive " << isActive << G4endl 439 << " isActive " << isActive << G4endl; 272 */ 440 */ 273 return gPathLength; 441 return gPathLength; 274 } 442 } 275 443 276 //....oooOO0OOooo........oooOO0OOooo........oo 444 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 277 445 278 G4double 446 G4double 279 G4VMultipleScattering::PostStepGetPhysicalInte 447 G4VMultipleScattering::PostStepGetPhysicalInteractionLength( 280 const G4Track&, G4double, G4Forc 448 const G4Track&, G4double, G4ForceCondition* condition) 281 { 449 { 282 *condition = NotForced; 450 *condition = NotForced; 283 return DBL_MAX; 451 return DBL_MAX; 284 } 452 } 285 453 286 //....oooOO0OOooo........oooOO0OOooo........oo 454 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 287 455 288 G4VParticleChange* 456 G4VParticleChange* 289 G4VMultipleScattering::AlongStepDoIt(const G4T 457 G4VMultipleScattering::AlongStepDoIt(const G4Track& track, const G4Step& step) 290 { 458 { 291 fParticleChange.InitialiseMSC(track, step); << 459 fParticleChange.ProposeMomentumDirection( 292 fNewPosition = fParticleChange.GetProposedPo << 460 step.GetPostStepPoint()->GetMomentumDirection()); >> 461 fNewPosition = step.GetPostStepPoint()->GetPosition(); >> 462 fParticleChange.ProposePosition(fNewPosition); 293 fPositionChanged = false; 463 fPositionChanged = false; 294 464 295 G4double geomLength = step.GetStepLength(); 465 G4double geomLength = step.GetStepLength(); 296 466 297 // very small step - no msc 467 // very small step - no msc 298 if(!isActive) { 468 if(!isActive) { 299 tPathLength = geomLength; 469 tPathLength = geomLength; 300 470 301 // sample msc 471 // sample msc 302 } else { 472 } else { 303 G4double range = 473 G4double range = 304 currentModel->GetRange(currParticle,trac 474 currentModel->GetRange(currParticle,track.GetKineticEnergy(), 305 track.GetMaterial 475 track.GetMaterialCutsCouple()); 306 476 307 tPathLength = currentModel->ComputeTrueSte 477 tPathLength = currentModel->ComputeTrueStepLength(geomLength); 308 478 309 /* 479 /* 310 if(currParticle->GetPDGMass() > 0.9*GeV) 480 if(currParticle->GetPDGMass() > 0.9*GeV) 311 G4cout << "G4VMsc::AlongStepDoIt: GeomLeng 481 G4cout << "G4VMsc::AlongStepDoIt: GeomLength= " 312 << geomLength 482 << geomLength 313 << " tPathLength= " << tPathLength 483 << " tPathLength= " << tPathLength 314 << " physStepLimit= " << physStepLi 484 << " physStepLimit= " << physStepLimit 315 << " dr= " << range - tPathLength 485 << " dr= " << range - tPathLength 316 << " ekin= " << track.GetKineticEne 486 << " ekin= " << track.GetKineticEnergy() << G4endl; 317 */ 487 */ 318 // protection against wrong t->g->t conver 488 // protection against wrong t->g->t conversion 319 tPathLength = std::min(tPathLength, physSt 489 tPathLength = std::min(tPathLength, physStepLimit); 320 490 321 // do not sample scattering at the last or 491 // do not sample scattering at the last or at a small step 322 if(tPathLength < range && tPathLength > ge 492 if(tPathLength < range && tPathLength > geomMin) { 323 493 324 static const G4double minSafety = 1.20*C 494 static const G4double minSafety = 1.20*CLHEP::nm; 325 static const G4double sFact = 0.99; 495 static const G4double sFact = 0.99; 326 496 327 G4ThreeVector displacement = currentMode 497 G4ThreeVector displacement = currentModel->SampleScattering( 328 step.GetPostStepPoint()->GetMomentumDi 498 step.GetPostStepPoint()->GetMomentumDirection(),minSafety); 329 499 330 G4double r2 = displacement.mag2(); 500 G4double r2 = displacement.mag2(); 331 //G4cout << " R= " << sqrt(r2) << " R 501 //G4cout << " R= " << sqrt(r2) << " Rmin= " << sqrt(minDisplacement2) 332 // << " flag= " << fDispBeyondSafety 502 // << " flag= " << fDispBeyondSafety << G4endl; 333 if(r2 > minDisplacement2) { 503 if(r2 > minDisplacement2) { 334 504 335 fPositionChanged = true; 505 fPositionChanged = true; 336 G4double dispR = std::sqrt(r2); 506 G4double dispR = std::sqrt(r2); 337 G4double postSafety = 507 G4double postSafety = 338 sFact*safetyHelper->ComputeSafety(fN 508 sFact*safetyHelper->ComputeSafety(fNewPosition, dispR); 339 //G4cout<<" R= "<< dispR<<" postSaf 509 //G4cout<<" R= "<< dispR<<" postSafety= "<<postSafety<<G4endl; 340 510 341 // far away from geometry boundary 511 // far away from geometry boundary 342 if(postSafety > 0.0 && dispR <= postSa 512 if(postSafety > 0.0 && dispR <= postSafety) { 343 fNewPosition += displacement; 513 fNewPosition += displacement; 344 514 345 //near the boundary 515 //near the boundary 346 } else { 516 } else { 347 // displaced point is definitely wit 517 // displaced point is definitely within the volume 348 //G4cout<<" R= "<<dispR<<" postSa 518 //G4cout<<" R= "<<dispR<<" postSafety= "<<postSafety<<G4endl; 349 if(dispR < postSafety) { 519 if(dispR < postSafety) { 350 fNewPosition += displacement; 520 fNewPosition += displacement; 351 521 352 // reduced displacement 522 // reduced displacement 353 } else if(postSafety > geomMin) { 523 } else if(postSafety > geomMin) { 354 fNewPosition += displacement*(post 524 fNewPosition += displacement*(postSafety/dispR); 355 525 356 // very small postSafety 526 // very small postSafety 357 } else { 527 } else { 358 fPositionChanged = false; 528 fPositionChanged = false; 359 } 529 } 360 } 530 } 361 if(fPositionChanged) { 531 if(fPositionChanged) { 362 safetyHelper->ReLocateWithinVolume(f 532 safetyHelper->ReLocateWithinVolume(fNewPosition); 363 fParticleChange.ProposePosition(fNew 533 fParticleChange.ProposePosition(fNewPosition); 364 } 534 } 365 } 535 } 366 } 536 } 367 } 537 } 368 fParticleChange.ProposeTrueStepLength(tPathL 538 fParticleChange.ProposeTrueStepLength(tPathLength); 369 return &fParticleChange; 539 return &fParticleChange; 370 } 540 } 371 541 372 //....oooOO0OOooo........oooOO0OOooo........oo 542 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 373 543 >> 544 G4VParticleChange* >> 545 G4VMultipleScattering::PostStepDoIt(const G4Track& track, const G4Step&) >> 546 { >> 547 fParticleChange.Initialize(track); >> 548 return &fParticleChange; >> 549 } >> 550 >> 551 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... >> 552 374 G4double G4VMultipleScattering::GetContinuousS 553 G4double G4VMultipleScattering::GetContinuousStepLimit( 375 const G 554 const G4Track& track, 376 G4doubl 555 G4double previousStepSize, 377 G4doubl 556 G4double currentMinimalStep, 378 G4doubl 557 G4double& currentSafety) 379 { 558 { 380 G4GPILSelection selection = NotCandidateForS 559 G4GPILSelection selection = NotCandidateForSelection; 381 G4double x = AlongStepGetPhysicalInteraction 560 G4double x = AlongStepGetPhysicalInteractionLength(track,previousStepSize, 382 561 currentMinimalStep, 383 562 currentSafety, 384 563 &selection); 385 return x; 564 return x; 386 } 565 } 387 566 388 //....oooOO0OOooo........oooOO0OOooo........oo 567 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 389 568 390 G4double G4VMultipleScattering::ContinuousStep 569 G4double G4VMultipleScattering::ContinuousStepLimit( 391 const G 570 const G4Track& track, 392 G4doubl 571 G4double previousStepSize, 393 G4doubl 572 G4double currentMinimalStep, 394 G4doubl 573 G4double& currentSafety) 395 { 574 { 396 return GetContinuousStepLimit(track,previous 575 return GetContinuousStepLimit(track,previousStepSize,currentMinimalStep, 397 currentSafety) 576 currentSafety); 398 } 577 } 399 578 400 //....oooOO0OOooo........oooOO0OOooo........oo 579 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 401 580 402 G4double G4VMultipleScattering::GetMeanFreePat 581 G4double G4VMultipleScattering::GetMeanFreePath( 403 const G4Track&, G4double, G4Forc 582 const G4Track&, G4double, G4ForceCondition* condition) 404 { 583 { 405 *condition = Forced; 584 *condition = Forced; 406 return DBL_MAX; 585 return DBL_MAX; 407 } 586 } 408 587 409 //....oooOO0OOooo........oooOO0OOooo........oo 588 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 410 589 411 G4bool 590 G4bool 412 G4VMultipleScattering::StorePhysicsTable(const 591 G4VMultipleScattering::StorePhysicsTable(const G4ParticleDefinition* part, 413 const 592 const G4String& directory, 414 G4boo 593 G4bool ascii) 415 { 594 { 416 G4bool yes = true; 595 G4bool yes = true; 417 if(part != firstParticle || !emManager->IsMa << 596 if(part != firstParticle) { return yes; } 418 << 597 const G4VMultipleScattering* masterProcess = 419 return G4EmTableUtil::StoreMscTable(this, pa << 598 static_cast<const G4VMultipleScattering*>(GetMasterProcess()); 420 numberOfModels, verboseLevel, << 599 if(masterProcess && masterProcess != this) { return yes; } 421 ascii); << 600 >> 601 G4int nmod = modelManager->NumberOfModels(); >> 602 static const G4String ss[4] = {"1","2","3","4"}; >> 603 for(G4int i=0; i<nmod; ++i) { >> 604 G4VEmModel* msc = modelManager->GetModel(i); >> 605 yes = true; >> 606 G4PhysicsTable* table = msc->GetCrossSectionTable(); >> 607 if (table) { >> 608 G4int j = std::min(i,3); >> 609 G4String name = >> 610 GetPhysicsTableFileName(part,directory,"LambdaMod"+ss[j],ascii); >> 611 yes = table->StorePhysicsTable(name,ascii); >> 612 >> 613 if ( yes ) { >> 614 if ( verboseLevel>0 ) { >> 615 G4cout << "Physics table are stored for " >> 616 << part->GetParticleName() >> 617 << " and process " << GetProcessName() >> 618 << " with a name <" << name << "> " << G4endl; >> 619 } >> 620 } else { >> 621 G4cout << "Fail to store Physics Table for " >> 622 << part->GetParticleName() >> 623 << " and process " << GetProcessName() >> 624 << " in the directory <" << directory >> 625 << "> " << G4endl; >> 626 } >> 627 } >> 628 } >> 629 return yes; 422 } 630 } 423 631 424 //....oooOO0OOooo........oooOO0OOooo........oo 632 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 425 633 426 G4bool 634 G4bool 427 G4VMultipleScattering::RetrievePhysicsTable(co 635 G4VMultipleScattering::RetrievePhysicsTable(const G4ParticleDefinition*, 428 co 636 const G4String&, 429 G4 637 G4bool) 430 { 638 { 431 return true; 639 return true; 432 } 640 } 433 641 434 //....oooOO0OOooo........oooOO0OOooo........oo 642 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 435 643 >> 644 void G4VMultipleScattering::SetIonisation(G4VEnergyLossProcess* p) >> 645 { >> 646 for(G4int i=0; i<numberOfModels; ++i) { >> 647 G4VMscModel* msc = static_cast<G4VMscModel*>(GetModelByIndex(i, true)); >> 648 if(msc) { msc->SetIonisation(p, firstParticle); } >> 649 } >> 650 } >> 651 >> 652 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... >> 653 436 void G4VMultipleScattering::ProcessDescription 654 void G4VMultipleScattering::ProcessDescription(std::ostream& outFile) const 437 { 655 { 438 if(nullptr != firstParticle) { << 656 if(firstParticle) { 439 StreamInfo(outFile, *firstParticle, true); 657 StreamInfo(outFile, *firstParticle, true); 440 } 658 } 441 } 659 } 442 660 443 //....oooOO0OOooo........oooOO0OOooo........oo 661 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 444 662 445 663