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), >> 99 fDispBeyondSafety(false), 88 fNewPosition(0.,0.,0.), 100 fNewPosition(0.,0.,0.), 89 fNewDirection(0.,0.,1.) 101 fNewDirection(0.,0.,1.) 90 { 102 { 91 theParameters = G4EmParameters::Instance(); 103 theParameters = G4EmParameters::Instance(); 92 SetVerboseLevel(1); 104 SetVerboseLevel(1); 93 SetProcessSubType(fMultipleScattering); 105 SetProcessSubType(fMultipleScattering); >> 106 if("ionmsc" == name) { firstParticle = G4GenericIon::GenericIon(); } 94 107 95 lowestKinEnergy = 10*CLHEP::eV; 108 lowestKinEnergy = 10*CLHEP::eV; 96 109 97 geomMin = 0.05*CLHEP::nm; << 110 physStepLimit = gPathLength = tPathLength = 0.0; >> 111 fIonisation = nullptr; >> 112 >> 113 geomMin = 0.05*CLHEP::nm; 98 minDisplacement2 = geomMin*geomMin; 114 minDisplacement2 = geomMin*geomMin; 99 115 100 pParticleChange = &fParticleChange; 116 pParticleChange = &fParticleChange; >> 117 safetyHelper = nullptr; >> 118 fPositionChanged = false; >> 119 isActive = false; 101 120 >> 121 currentModel = nullptr; 102 modelManager = new G4EmModelManager(); 122 modelManager = new G4EmModelManager(); 103 emManager = G4LossTableManager::Instance(); 123 emManager = G4LossTableManager::Instance(); 104 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; } << 122 G4VEmFluctuationModel* fm = nullptr; 146 G4VEmFluctuationModel* fm = nullptr; 123 modelManager->AddEmModel(order, ptr, fm, reg << 147 modelManager->AddEmModel(order, p, fm, region); 124 ptr->SetParticleChange(pParticleChange); << 148 if(p) { p->SetParticleChange(pParticleChange); } 125 } 149 } 126 150 127 //....oooOO0OOooo........oooOO0OOooo........oo 151 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 128 152 129 void G4VMultipleScattering::SetEmModel(G4VMscM << 153 void G4VMultipleScattering::SetEmModel(G4VMscModel* p, size_t) 130 { 154 { 131 if(nullptr == ptr) { return; } << 155 for(auto & msc : mscModels) { if(msc == p) { return; } } 132 if(!mscModels.empty()) { << 156 mscModels.push_back(p); 133 for(auto & msc : mscModels) { if(msc == pt << 157 } 134 } << 158 135 mscModels.push_back(ptr); << 159 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... >> 160 >> 161 G4VMscModel* G4VMultipleScattering::EmModel(size_t index) const >> 162 { >> 163 return (index < mscModels.size()) ? mscModels[index] : nullptr; 136 } 164 } 137 165 138 //....oooOO0OOooo........oooOO0OOooo........oo 166 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 139 167 140 void 168 void 141 G4VMultipleScattering::PreparePhysicsTable(con 169 G4VMultipleScattering::PreparePhysicsTable(const G4ParticleDefinition& part) 142 { 170 { >> 171 if(1 < verboseLevel) { >> 172 G4cout << "### G4VMultipleScattering::PrepearPhysicsTable() for " >> 173 << GetProcessName() >> 174 << " and particle " << part.GetParticleName() >> 175 << G4endl; >> 176 } 143 G4bool master = emManager->IsMaster(); 177 G4bool master = emManager->IsMaster(); 144 if (nullptr == firstParticle) { firstParticl << 145 178 146 emManager->PreparePhysicsTable(&part, this); << 179 if(!firstParticle) { firstParticle = ∂ } >> 180 if(part.GetParticleType() == "nucleus") { >> 181 stepLimit = fMinimal; >> 182 latDisplacement = false; >> 183 facrange = 0.2; >> 184 G4String pname = part.GetParticleName(); >> 185 if(pname != "deuteron" && pname != "triton" && >> 186 pname != "alpha+" && pname != "helium" && >> 187 pname != "alpha" && pname != "He3" && >> 188 pname != "hydrogen") { >> 189 >> 190 const G4ParticleDefinition* theGenericIon = >> 191 G4ParticleTable::GetParticleTable()->FindParticle("GenericIon"); >> 192 if(&part == theGenericIon) { isIon = true; } >> 193 >> 194 if(theGenericIon && firstParticle != theGenericIon) { >> 195 G4ProcessManager* pm = theGenericIon->GetProcessManager(); >> 196 G4ProcessVector* v = pm->GetAlongStepProcessVector(); >> 197 size_t n = v->size(); >> 198 for(size_t j=0; j<n; ++j) { >> 199 if((*v)[j] == this) { >> 200 firstParticle = theGenericIon; >> 201 isIon = true; >> 202 break; >> 203 } >> 204 } >> 205 } >> 206 } >> 207 } >> 208 >> 209 emManager->PreparePhysicsTable(&part, this, master); 147 currParticle = nullptr; 210 currParticle = nullptr; 148 211 >> 212 if(1 < verboseLevel) { >> 213 G4cout << "### G4VMultipleScattering::PrepearPhysicsTable() for " >> 214 << GetProcessName() >> 215 << " and particle " << part.GetParticleName() >> 216 << " local particle " << firstParticle->GetParticleName() >> 217 << " isIon: " << isIon << " isMaster: " << master >> 218 << G4endl; >> 219 } >> 220 149 if(firstParticle == &part) { 221 if(firstParticle == &part) { 150 baseMat = emManager->GetTableBuilder()->Ge << 151 G4EmTableUtil::PrepareMscProcess(this, par << 152 stepLimit, facrange, << 153 latDisplacement, master, << 154 isIon, baseMat); << 155 222 >> 223 // initialise process >> 224 InitialiseProcess(firstParticle); >> 225 >> 226 // heavy particles and not ions >> 227 if(!isIon) { >> 228 if(part.GetPDGMass() > MeV) { >> 229 stepLimit = theParameters->MscMuHadStepLimitType(); >> 230 facrange = theParameters->MscMuHadRangeFactor(); >> 231 latDisplacement = theParameters->MuHadLateralDisplacement(); >> 232 } else { >> 233 stepLimit = theParameters->MscStepLimitType(); >> 234 facrange = theParameters->MscRangeFactor(); >> 235 latDisplacement = theParameters->LateralDisplacement(); >> 236 } >> 237 if(latDisplacement) { >> 238 fDispBeyondSafety = theParameters->LatDisplacementBeyondSafety(); >> 239 } >> 240 } >> 241 if(master) { SetVerboseLevel(theParameters->Verbose()); } >> 242 else { SetVerboseLevel(theParameters->WorkerVerbose()); } >> 243 >> 244 // initialisation of models 156 numberOfModels = modelManager->NumberOfMod 245 numberOfModels = modelManager->NumberOfModels(); 157 currentModel = GetModelByIndex(0); << 246 /* >> 247 G4cout << "### G4VMultipleScattering::PreparePhysicsTable() for " >> 248 << GetProcessName() >> 249 << " and particle " << part.GetParticleName() >> 250 << " Nmod= " << numberOfModels << " " << this >> 251 << G4endl; >> 252 */ >> 253 for(G4int i=0; i<numberOfModels; ++i) { >> 254 G4VMscModel* msc = static_cast<G4VMscModel*>(GetModelByIndex(i)); >> 255 if(!msc) { continue; } >> 256 msc->SetIonisation(nullptr, firstParticle); >> 257 msc->SetMasterThread(master); >> 258 currentModel = msc; >> 259 msc->SetStepLimitType(stepLimit); >> 260 msc->SetLateralDisplasmentFlag(latDisplacement); >> 261 msc->SetSkin(theParameters->MscSkin()); >> 262 msc->SetRangeFactor(facrange); >> 263 msc->SetGeomFactor(theParameters->MscGeomFactor()); >> 264 msc->SetPolarAngleLimit(theParameters->MscThetaLimit()); >> 265 G4double emax = >> 266 std::min(msc->HighEnergyLimit(),theParameters->MaxKinEnergy()); >> 267 msc->SetHighEnergyLimit(emax); >> 268 } 158 269 159 if (nullptr == safetyHelper) { << 270 modelManager->Initialise(firstParticle, G4Electron::Electron(), >> 271 10.0, verboseLevel); >> 272 >> 273 if(!safetyHelper) { 160 safetyHelper = G4TransportationManager:: 274 safetyHelper = G4TransportationManager::GetTransportationManager() 161 ->GetSafetyHelper(); << 275 ->GetSafetyHelper(); 162 safetyHelper->InitialiseHelper(); 276 safetyHelper->InitialiseHelper(); 163 } 277 } 164 } 278 } 165 } 279 } 166 280 167 //....oooOO0OOooo........oooOO0OOooo........oo 281 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 168 282 169 void G4VMultipleScattering::BuildPhysicsTable( 283 void G4VMultipleScattering::BuildPhysicsTable(const G4ParticleDefinition& part) 170 { 284 { >> 285 G4String num = part.GetParticleName(); 171 G4bool master = emManager->IsMaster(); 286 G4bool master = emManager->IsMaster(); 172 << 287 if(1 < verboseLevel) { 173 if(firstParticle == &part) { << 288 G4cout << "### G4VMultipleScattering::BuildPhysicsTable() for " 174 emManager->BuildPhysicsTable(&part); << 289 << GetProcessName() >> 290 << " and particle " << num << " isIon: " << isIon >> 291 << " IsMaster: " << master << G4endl; 175 } 292 } 176 const G4VMultipleScattering* ptr = this; << 293 const G4VMultipleScattering* masterProcess = 177 if(!master) { << 294 static_cast<const G4VMultipleScattering*>(GetMasterProcess()); 178 ptr = static_cast<const G4VMultipleScatter << 295 >> 296 if(firstParticle == &part) { >> 297 /* >> 298 G4cout << "### G4VMultipleScattering::BuildPhysicsTable() for " >> 299 << GetProcessName() >> 300 << " and particle " << num >> 301 << " IsMaster= " << G4LossTableManager::Instance()->IsMaster() >> 302 << " " << this >> 303 << G4endl; >> 304 */ >> 305 emManager->BuildPhysicsTable(firstParticle); >> 306 >> 307 if(!master) { >> 308 // initialisation of models >> 309 /* >> 310 G4cout << "### G4VMultipleScattering::BuildPhysicsTable() for " >> 311 << GetProcessName() >> 312 << " and particle " << num >> 313 << " Nmod= " << numberOfModels << " " << this >> 314 << G4endl; >> 315 */ >> 316 for(G4int i=0; i<numberOfModels; ++i) { >> 317 G4VMscModel* msc = static_cast<G4VMscModel*>(GetModelByIndex(i)); >> 318 if(!msc) { continue; } >> 319 G4VMscModel* msc0= >> 320 static_cast<G4VMscModel*>(masterProcess->GetModelByIndex(i)); >> 321 msc->SetCrossSectionTable(msc0->GetCrossSectionTable(), false); >> 322 msc->InitialiseLocal(firstParticle, msc0); >> 323 } >> 324 } 179 } 325 } 180 326 181 G4EmTableUtil::BuildMscProcess(this, ptr, pa << 327 // explicitly defined printout by particle name 182 numberOfModels, master); << 328 if(1 < verboseLevel || >> 329 (0 < verboseLevel && (num == "e-" || >> 330 num == "e+" || num == "mu+" || >> 331 num == "mu-" || num == "proton"|| >> 332 num == "pi+" || num == "pi-" || >> 333 num == "kaon+" || num == "kaon-" || >> 334 num == "alpha" || num == "anti_proton" || >> 335 num == "GenericIon" || num == "alpha+" || >> 336 num == "alpha++" ))) >> 337 { >> 338 StreamInfo(G4cout, part); >> 339 } >> 340 >> 341 if(1 < verboseLevel) { >> 342 G4cout << "### G4VMultipleScattering::BuildPhysicsTable() done for " >> 343 << GetProcessName() >> 344 << " and particle " << num >> 345 << G4endl; >> 346 } 183 } 347 } 184 348 185 //....oooOO0OOooo........oooOO0OOooo........oo 349 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 186 350 187 void G4VMultipleScattering::StreamInfo(std::os 351 void G4VMultipleScattering::StreamInfo(std::ostream& outFile, 188 const G4ParticleDefinition& 352 const G4ParticleDefinition& part, G4bool rst) const 189 { 353 { 190 G4String indent = (rst ? " " : ""); 354 G4String indent = (rst ? " " : ""); 191 outFile << G4endl << indent << GetProcessNam 355 outFile << G4endl << indent << GetProcessName() << ": "; 192 if (!rst) outFile << " for " << part.GetPart << 356 if (!rst) outFile << " for " << part.GetParticleName(); 193 outFile << " SubType= " << GetProcessSubTy << 357 outFile << " SubType= " << GetProcessSubType() << G4endl; >> 358 StreamProcessInfo(outFile); 194 modelManager->DumpModelList(outFile, verbose 359 modelManager->DumpModelList(outFile, verboseLevel); 195 } 360 } 196 361 197 //....oooOO0OOooo........oooOO0OOooo........oo 362 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 198 363 199 void G4VMultipleScattering::StartTracking(G4Tr 364 void G4VMultipleScattering::StartTracking(G4Track* track) 200 { 365 { 201 G4VEnergyLossProcess* eloss = nullptr; 366 G4VEnergyLossProcess* eloss = nullptr; 202 if(track->GetParticleDefinition() != currPar 367 if(track->GetParticleDefinition() != currParticle) { 203 currParticle = track->GetParticleDefinitio 368 currParticle = track->GetParticleDefinition(); 204 fIonisation = emManager->GetEnergyLossProc 369 fIonisation = emManager->GetEnergyLossProcess(currParticle); 205 eloss = fIonisation; 370 eloss = fIonisation; 206 } 371 } >> 372 /* >> 373 G4cout << "G4VMultipleScattering::StartTracking Nmod= " << numberOfModels >> 374 << " " << currParticle->GetParticleName() >> 375 << " E(MeV)= " << track->GetKineticEnergy() >> 376 << " Ion= " << eloss << " " << fIonisation << " IsMaster= " >> 377 << G4LossTableManager::Instance()->IsMaster() >> 378 << G4endl; >> 379 */ 207 for(G4int i=0; i<numberOfModels; ++i) { 380 for(G4int i=0; i<numberOfModels; ++i) { 208 G4VMscModel* msc = GetModelByIndex(i); << 381 /* >> 382 G4cout << "Next model " << i << " " << msc >> 383 << " Emin= " << msc->LowEnergyLimit() >> 384 << " Emax= " << msc->HighEnergyLimit() >> 385 << " Eact= " << msc->LowEnergyActivationLimit() << G4endl; >> 386 */ >> 387 G4VEmModel* msc = GetModelByIndex(i); 209 msc->StartTracking(track); 388 msc->StartTracking(track); 210 if(nullptr != eloss) { << 389 if(eloss) { 211 msc->SetIonisation(eloss, currParticle); << 390 G4VMscModel* mscmod = static_cast<G4VMscModel*>(msc); >> 391 if(mscmod) { mscmod->SetIonisation(fIonisation, currParticle); } 212 } 392 } 213 } 393 } 214 } << 394 } 215 395 216 //....oooOO0OOooo........oooOO0OOooo........oo 396 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 217 397 218 G4double G4VMultipleScattering::AlongStepGetPh 398 G4double G4VMultipleScattering::AlongStepGetPhysicalInteractionLength( 219 const G4Track& tr 399 const G4Track& track, 220 G4double, 400 G4double, 221 G4double currentM 401 G4double currentMinimalStep, 222 G4double&, 402 G4double&, 223 G4GPILSelection* 403 G4GPILSelection* selection) 224 { 404 { 225 // get Step limit proposed by the process 405 // get Step limit proposed by the process 226 *selection = NotCandidateForSelection; 406 *selection = NotCandidateForSelection; 227 physStepLimit = gPathLength = tPathLength = 407 physStepLimit = gPathLength = tPathLength = currentMinimalStep; 228 408 229 G4double ekin = track.GetKineticEnergy(); 409 G4double ekin = track.GetKineticEnergy(); 230 /* 410 /* 231 G4cout << "MSC::AlongStepGPIL: Ekin= " << ek 411 G4cout << "MSC::AlongStepGPIL: Ekin= " << ekin 232 << " " << currParticle->GetParticleN 412 << " " << currParticle->GetParticleName() 233 << " currMod " << currentModel 413 << " currMod " << currentModel 234 << G4endl; 414 << G4endl; 235 */ 415 */ 236 // isIon flag is used only to select a model 416 // isIon flag is used only to select a model 237 if(isIon) { 417 if(isIon) { 238 ekin *= proton_mass_c2/track.GetParticleDe 418 ekin *= proton_mass_c2/track.GetParticleDefinition()->GetPDGMass(); 239 } 419 } 240 const G4MaterialCutsCouple* couple = track.G << 241 420 242 // select new model, static cast is possible << 421 // select new model 243 if(1 < numberOfModels) { 422 if(1 < numberOfModels) { 244 currentModel = << 423 currentModel = static_cast<G4VMscModel*>( 245 static_cast<G4VMscModel*>(SelectModel(ek << 424 SelectModel(ekin,track.GetMaterialCutsCouple()->GetIndex())); 246 } 425 } 247 currentModel->SetCurrentCouple(couple); << 248 // msc is active is model is active, energy 426 // msc is active is model is active, energy above the limit, 249 // and step size is above the limit; 427 // and step size is above the limit; 250 // if it is active msc may limit the step 428 // if it is active msc may limit the step 251 if(currentModel->IsActive(ekin) && tPathLeng 429 if(currentModel->IsActive(ekin) && tPathLength > geomMin 252 && ekin >= lowestKinEnergy) { 430 && ekin >= lowestKinEnergy) { 253 isActive = true; 431 isActive = true; 254 tPathLength = 432 tPathLength = 255 currentModel->ComputeTruePathLengthLimit 433 currentModel->ComputeTruePathLengthLimit(track, gPathLength); 256 if (tPathLength < physStepLimit) { 434 if (tPathLength < physStepLimit) { 257 *selection = CandidateForSelection; 435 *selection = CandidateForSelection; 258 } 436 } 259 } else { << 437 } else { isActive = false; } 260 isActive = false; << 261 gPathLength = DBL_MAX; << 262 } << 263 438 264 //if(currParticle->GetPDGMass() > GeV) 439 //if(currParticle->GetPDGMass() > GeV) 265 /* 440 /* 266 G4cout << "MSC::AlongStepGPIL: Ekin= " << ek 441 G4cout << "MSC::AlongStepGPIL: Ekin= " << ekin 267 << " " << currParticle->GetParticleN 442 << " " << currParticle->GetParticleName() 268 << " gPathLength= " << gPathLength 443 << " gPathLength= " << gPathLength 269 << " tPathLength= " << tPathLength 444 << " tPathLength= " << tPathLength 270 << " currentMinimalStep= " << current 445 << " currentMinimalStep= " << currentMinimalStep 271 << " isActive " << isActive << G4endl 446 << " isActive " << isActive << G4endl; 272 */ 447 */ 273 return gPathLength; 448 return gPathLength; 274 } 449 } 275 450 276 //....oooOO0OOooo........oooOO0OOooo........oo 451 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 277 452 278 G4double 453 G4double 279 G4VMultipleScattering::PostStepGetPhysicalInte 454 G4VMultipleScattering::PostStepGetPhysicalInteractionLength( 280 const G4Track&, G4double, G4Forc 455 const G4Track&, G4double, G4ForceCondition* condition) 281 { 456 { 282 *condition = NotForced; 457 *condition = NotForced; 283 return DBL_MAX; 458 return DBL_MAX; 284 } 459 } 285 460 286 //....oooOO0OOooo........oooOO0OOooo........oo 461 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 287 462 288 G4VParticleChange* 463 G4VParticleChange* 289 G4VMultipleScattering::AlongStepDoIt(const G4T 464 G4VMultipleScattering::AlongStepDoIt(const G4Track& track, const G4Step& step) 290 { 465 { 291 fParticleChange.InitialiseMSC(track, step); << 466 fParticleChange.ProposeMomentumDirection( 292 fNewPosition = fParticleChange.GetProposedPo << 467 step.GetPostStepPoint()->GetMomentumDirection()); >> 468 fNewPosition = step.GetPostStepPoint()->GetPosition(); >> 469 fParticleChange.ProposePosition(fNewPosition); 293 fPositionChanged = false; 470 fPositionChanged = false; 294 471 295 G4double geomLength = step.GetStepLength(); 472 G4double geomLength = step.GetStepLength(); 296 473 297 // very small step - no msc 474 // very small step - no msc 298 if(!isActive) { 475 if(!isActive) { 299 tPathLength = geomLength; 476 tPathLength = geomLength; 300 477 301 // sample msc 478 // sample msc 302 } else { 479 } else { 303 G4double range = 480 G4double range = 304 currentModel->GetRange(currParticle,trac 481 currentModel->GetRange(currParticle,track.GetKineticEnergy(), 305 track.GetMaterial 482 track.GetMaterialCutsCouple()); 306 483 307 tPathLength = currentModel->ComputeTrueSte 484 tPathLength = currentModel->ComputeTrueStepLength(geomLength); 308 485 309 /* 486 /* 310 if(currParticle->GetPDGMass() > 0.9*GeV) 487 if(currParticle->GetPDGMass() > 0.9*GeV) 311 G4cout << "G4VMsc::AlongStepDoIt: GeomLeng 488 G4cout << "G4VMsc::AlongStepDoIt: GeomLength= " 312 << geomLength 489 << geomLength 313 << " tPathLength= " << tPathLength 490 << " tPathLength= " << tPathLength 314 << " physStepLimit= " << physStepLi 491 << " physStepLimit= " << physStepLimit 315 << " dr= " << range - tPathLength 492 << " dr= " << range - tPathLength 316 << " ekin= " << track.GetKineticEne 493 << " ekin= " << track.GetKineticEnergy() << G4endl; 317 */ 494 */ 318 // protection against wrong t->g->t conver 495 // protection against wrong t->g->t conversion 319 tPathLength = std::min(tPathLength, physSt 496 tPathLength = std::min(tPathLength, physStepLimit); 320 497 321 // do not sample scattering at the last or 498 // do not sample scattering at the last or at a small step 322 if(tPathLength < range && tPathLength > ge 499 if(tPathLength < range && tPathLength > geomMin) { 323 500 324 static const G4double minSafety = 1.20*C 501 static const G4double minSafety = 1.20*CLHEP::nm; 325 static const G4double sFact = 0.99; 502 static const G4double sFact = 0.99; 326 503 327 G4ThreeVector displacement = currentMode 504 G4ThreeVector displacement = currentModel->SampleScattering( 328 step.GetPostStepPoint()->GetMomentumDi 505 step.GetPostStepPoint()->GetMomentumDirection(),minSafety); 329 506 330 G4double r2 = displacement.mag2(); 507 G4double r2 = displacement.mag2(); 331 //G4cout << " R= " << sqrt(r2) << " R 508 //G4cout << " R= " << sqrt(r2) << " Rmin= " << sqrt(minDisplacement2) 332 // << " flag= " << fDispBeyondSafety 509 // << " flag= " << fDispBeyondSafety << G4endl; 333 if(r2 > minDisplacement2) { 510 if(r2 > minDisplacement2) { 334 511 335 fPositionChanged = true; 512 fPositionChanged = true; 336 G4double dispR = std::sqrt(r2); << 513 G4double dispR = std::sqrt(r2); 337 G4double postSafety = 514 G4double postSafety = 338 sFact*safetyHelper->ComputeSafety(fN << 515 sFact*safetyHelper->ComputeSafety(fNewPosition, dispR); 339 //G4cout<<" R= "<< dispR<<" postSaf 516 //G4cout<<" R= "<< dispR<<" postSafety= "<<postSafety<<G4endl; 340 517 341 // far away from geometry boundary 518 // far away from geometry boundary 342 if(postSafety > 0.0 && dispR <= postSa 519 if(postSafety > 0.0 && dispR <= postSafety) { 343 fNewPosition += displacement; 520 fNewPosition += displacement; 344 521 345 //near the boundary << 522 //near the boundary 346 } else { 523 } else { 347 // displaced point is definitely wit 524 // displaced point is definitely within the volume 348 //G4cout<<" R= "<<dispR<<" postSa 525 //G4cout<<" R= "<<dispR<<" postSafety= "<<postSafety<<G4endl; 349 if(dispR < postSafety) { 526 if(dispR < postSafety) { 350 fNewPosition += displacement; 527 fNewPosition += displacement; 351 528 352 // reduced displacement << 529 // optional extra mechanism is applied only if a particle 353 } else if(postSafety > geomMin) { << 530 // is stopped by the boundary 354 fNewPosition += displacement*(post << 531 } else if(fDispBeyondSafety && 0.0 == postSafety) { >> 532 fNewPosition += displacement; >> 533 G4double maxshift = >> 534 std::min(2.0*dispR, geomLength*(physStepLimit/tPathLength - 1.0)); >> 535 G4double dist = 0.0; >> 536 G4double safety = postSafety + dispR; >> 537 fNewDirection = *(fParticleChange.GetMomentumDirection()); >> 538 /* >> 539 G4cout << "##MSC before Recheck maxshift= " << maxshift >> 540 << " postsafety= " << postSafety >> 541 << " Ekin= " << track.GetKineticEnergy() >> 542 << " " << track.GetDefinition()->GetParticleName() >> 543 << G4endl; >> 544 */ >> 545 // check if it is possible to shift to the boundary >> 546 // and the shift is not large >> 547 if(safetyHelper->RecheckDistanceToCurrentBoundary(fNewPosition, >> 548 fNewDirection, maxshift, &dist, &safety) >> 549 && std::abs(dist) < maxshift) { >> 550 /* >> 551 G4cout << "##MSC after Recheck dist= " << dist >> 552 << " postsafety= " << postSafety >> 553 << " t= " << tPathLength >> 554 << " g= " << geomLength >> 555 << " p= " << physStepLimit >> 556 << G4endl; >> 557 */ >> 558 // shift is positive >> 559 if(dist >= 0.0) { >> 560 tPathLength *= (1.0 + dist/geomLength); >> 561 fNewPosition += dist*fNewDirection; >> 562 >> 563 // shift is negative cannot be larger than geomLength >> 564 } else { >> 565 maxshift = std::min(maxshift, geomLength); >> 566 if(0.0 < maxshift + dist) { >> 567 const G4ThreeVector& postpoint = step.GetPostStepPoint()->GetPosition(); >> 568 G4ThreeVector point = fNewPosition + dist*fNewDirection; >> 569 G4double R2 = (postpoint - point).mag2(); >> 570 G4double newdist = dist; >> 571 // check not more than 10 extra boundaries >> 572 for(G4int i=0; i<10; ++i) { >> 573 dist = 0.0; >> 574 if(safetyHelper->RecheckDistanceToCurrentBoundary( >> 575 point, fNewDirection, maxshift, &dist, &safety) >> 576 && std::abs(newdist + dist) < maxshift) { >> 577 point += dist*fNewDirection; >> 578 G4double R2new = (postpoint - point).mag2(); >> 579 //G4cout << "Backward i= " << i << " dist= " << dist >> 580 // << " R2= " << R2new << G4endl; >> 581 if(dist >= 0.0 || R2new > R2) { break; } >> 582 R2 = R2new; >> 583 fNewPosition = point; >> 584 newdist += dist; >> 585 } else { >> 586 break; >> 587 } >> 588 } >> 589 tPathLength *= (1.0 + newdist/geomLength); >> 590 // shift on boundary is not possible for negative disp >> 591 } else { >> 592 fNewPosition += displacement*(postSafety/dispR - 1.0); >> 593 } >> 594 } >> 595 // shift on boundary is not possible for any disp >> 596 } else { >> 597 fNewPosition += displacement*(postSafety/dispR - 1.0); >> 598 } >> 599 // reduced displacement >> 600 } else if(postSafety > geomMin) { >> 601 fNewPosition += displacement*(postSafety/dispR); 355 602 356 // very small postSafety 603 // very small postSafety 357 } else { 604 } else { 358 fPositionChanged = false; 605 fPositionChanged = false; 359 } 606 } 360 } 607 } 361 if(fPositionChanged) { << 608 if(fPositionChanged) { 362 safetyHelper->ReLocateWithinVolume(f << 609 safetyHelper->ReLocateWithinVolume(fNewPosition); 363 fParticleChange.ProposePosition(fNew << 610 fParticleChange.ProposePosition(fNewPosition); 364 } << 611 } 365 } 612 } 366 } 613 } 367 } 614 } 368 fParticleChange.ProposeTrueStepLength(tPathL 615 fParticleChange.ProposeTrueStepLength(tPathLength); 369 return &fParticleChange; 616 return &fParticleChange; 370 } 617 } 371 618 372 //....oooOO0OOooo........oooOO0OOooo........oo 619 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 373 620 >> 621 G4VParticleChange* >> 622 G4VMultipleScattering::PostStepDoIt(const G4Track& track, const G4Step&) >> 623 { >> 624 fParticleChange.Initialize(track); >> 625 return &fParticleChange; >> 626 } >> 627 >> 628 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... >> 629 374 G4double G4VMultipleScattering::GetContinuousS 630 G4double G4VMultipleScattering::GetContinuousStepLimit( 375 const G 631 const G4Track& track, 376 G4doubl 632 G4double previousStepSize, 377 G4doubl 633 G4double currentMinimalStep, 378 G4doubl 634 G4double& currentSafety) 379 { 635 { 380 G4GPILSelection selection = NotCandidateForS 636 G4GPILSelection selection = NotCandidateForSelection; 381 G4double x = AlongStepGetPhysicalInteraction 637 G4double x = AlongStepGetPhysicalInteractionLength(track,previousStepSize, 382 638 currentMinimalStep, 383 639 currentSafety, 384 640 &selection); 385 return x; 641 return x; 386 } 642 } 387 643 388 //....oooOO0OOooo........oooOO0OOooo........oo 644 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 389 645 390 G4double G4VMultipleScattering::ContinuousStep 646 G4double G4VMultipleScattering::ContinuousStepLimit( 391 const G 647 const G4Track& track, 392 G4doubl 648 G4double previousStepSize, 393 G4doubl 649 G4double currentMinimalStep, 394 G4doubl 650 G4double& currentSafety) 395 { 651 { 396 return GetContinuousStepLimit(track,previous 652 return GetContinuousStepLimit(track,previousStepSize,currentMinimalStep, 397 currentSafety) 653 currentSafety); 398 } 654 } 399 655 400 //....oooOO0OOooo........oooOO0OOooo........oo 656 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 401 657 402 G4double G4VMultipleScattering::GetMeanFreePat 658 G4double G4VMultipleScattering::GetMeanFreePath( 403 const G4Track&, G4double, G4Forc 659 const G4Track&, G4double, G4ForceCondition* condition) 404 { 660 { 405 *condition = Forced; 661 *condition = Forced; 406 return DBL_MAX; 662 return DBL_MAX; 407 } 663 } 408 664 409 //....oooOO0OOooo........oooOO0OOooo........oo 665 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 410 666 411 G4bool 667 G4bool 412 G4VMultipleScattering::StorePhysicsTable(const 668 G4VMultipleScattering::StorePhysicsTable(const G4ParticleDefinition* part, 413 const 669 const G4String& directory, 414 G4boo 670 G4bool ascii) 415 { 671 { 416 G4bool yes = true; 672 G4bool yes = true; 417 if(part != firstParticle || !emManager->IsMa << 673 if(part != firstParticle) { return yes; } 418 << 674 const G4VMultipleScattering* masterProcess = 419 return G4EmTableUtil::StoreMscTable(this, pa << 675 static_cast<const G4VMultipleScattering*>(GetMasterProcess()); 420 numberOfModels, verboseLevel, << 676 if(masterProcess && masterProcess != this) { return yes; } 421 ascii); << 677 >> 678 G4int nmod = modelManager->NumberOfModels(); >> 679 static const G4String ss[4] = {"1","2","3","4"}; >> 680 for(G4int i=0; i<nmod; ++i) { >> 681 G4VEmModel* msc = modelManager->GetModel(i); >> 682 yes = true; >> 683 G4PhysicsTable* table = msc->GetCrossSectionTable(); >> 684 if (table) { >> 685 G4int j = std::min(i,3); >> 686 G4String name = >> 687 GetPhysicsTableFileName(part,directory,"LambdaMod"+ss[j],ascii); >> 688 yes = table->StorePhysicsTable(name,ascii); >> 689 >> 690 if ( yes ) { >> 691 if ( verboseLevel>0 ) { >> 692 G4cout << "Physics table are stored for " >> 693 << part->GetParticleName() >> 694 << " and process " << GetProcessName() >> 695 << " with a name <" << name << "> " << G4endl; >> 696 } >> 697 } else { >> 698 G4cout << "Fail to store Physics Table for " >> 699 << part->GetParticleName() >> 700 << " and process " << GetProcessName() >> 701 << " in the directory <" << directory >> 702 << "> " << G4endl; >> 703 } >> 704 } >> 705 } >> 706 return yes; 422 } 707 } 423 708 424 //....oooOO0OOooo........oooOO0OOooo........oo 709 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo...... 425 710 426 G4bool 711 G4bool 427 G4VMultipleScattering::RetrievePhysicsTable(co 712 G4VMultipleScattering::RetrievePhysicsTable(const G4ParticleDefinition*, 428 co 713 const G4String&, 429 G4 714 G4bool) 430 { 715 { 431 return true; 716 return true; 432 } 717 } 433 718 434 //....oooOO0OOooo........oooOO0OOooo........oo 719 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 435 720 >> 721 void G4VMultipleScattering::SetIonisation(G4VEnergyLossProcess* p) >> 722 { >> 723 for(G4int i=0; i<numberOfModels; ++i) { >> 724 G4VMscModel* msc = static_cast<G4VMscModel*>(GetModelByIndex(i, true)); >> 725 if(msc) { msc->SetIonisation(p, firstParticle); } >> 726 } >> 727 } >> 728 >> 729 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... >> 730 436 void G4VMultipleScattering::ProcessDescription 731 void G4VMultipleScattering::ProcessDescription(std::ostream& outFile) const 437 { 732 { 438 if(nullptr != firstParticle) { << 733 if(firstParticle) { 439 StreamInfo(outFile, *firstParticle, true); 734 StreamInfo(outFile, *firstParticle, true); 440 } 735 } 441 } 736 } 442 737 443 //....oooOO0OOooo........oooOO0OOooo........oo 738 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo.... 444 739 445 740