Geant4 Cross Reference |
1 // 1 2 // ******************************************* 3 // * License and Disclaimer 4 // * 5 // * The Geant4 software is copyright of th 6 // * the Geant4 Collaboration. It is provided 7 // * conditions of the Geant4 Software License 8 // * LICENSE and available at http://cern.ch/ 9 // * include a list of copyright holders. 10 // * 11 // * Neither the authors of this software syst 12 // * institutes,nor the agencies providing fin 13 // * work make any representation or warran 14 // * regarding this software system or assum 15 // * use. Please see the license in the file 16 // * for the full disclaimer and the limitatio 17 // * 18 // * This code implementation is the result 19 // * technical work of the GEANT4 collaboratio 20 // * By using, copying, modifying or distri 21 // * any work based on the software) you ag 22 // * use in resulting scientific publicati 23 // * acceptance of all terms of the Geant4 Sof 24 // ******************************************* 25 // 26 // G4TransportationWithMsc 27 // 28 // Class Description: 29 // 30 // It is a generic process of transportation w 31 // in the step limitation and propagation. 32 // 33 // Original author: Jonas Hahnfeld, 2022 34 35 // ------------------------------------------- 36 // 37 //....oooOO0OOooo........oooOO0OOooo........oo 38 //....oooOO0OOooo........oooOO0OOooo........oo 39 40 #include "G4TransportationWithMsc.hh" 41 42 #include "G4DynamicParticle.hh" 43 #include "G4Electron.hh" 44 #include "G4EmConfigurator.hh" 45 #include "G4EmDataHandler.hh" 46 #include "G4LossTableBuilder.hh" 47 #include "G4LossTableManager.hh" 48 #include "G4ParticleChangeForGamma.hh" 49 #include "G4ParticleChangeForMSC.hh" 50 #include "G4PhysicalConstants.hh" 51 #include "G4PhysicsTableHelper.hh" 52 #include "G4PhysicsVector.hh" 53 #include "G4ProductionCutsTable.hh" 54 #include "G4Step.hh" 55 #include "G4StepPoint.hh" 56 #include "G4StepStatus.hh" 57 #include "G4Track.hh" 58 #include "G4TransportationProcessType.hh" 59 #include "G4VMscModel.hh" 60 61 //....oooOO0OOooo........oooOO0OOooo........oo 62 63 static constexpr G4double kLowestKinEnergy = 1 64 static constexpr G4double kGeomMin = 0.05 * CL 65 static constexpr G4double kMinDisplacement2 = 66 67 //....oooOO0OOooo........oooOO0OOooo........oo 68 69 G4TransportationWithMsc::G4TransportationWithM 70 : G4Transportation(verbosity, "Transportatio 71 { 72 SetProcessSubType(static_cast<G4int>(TRANSPO 73 SetVerboseLevel(1); 74 75 fEmManager = G4LossTableManager::Instance(); 76 fModelManager = new G4EmModelManager; 77 78 if (type == ScatteringType::MultipleScatteri 79 fParticleChangeForMSC = new G4ParticleChan 80 } 81 else if (type == ScatteringType::SingleScatt 82 fParticleChangeForSS = new G4ParticleChang 83 fSecondariesSS = new std::vector<G4Dynamic 84 } 85 86 G4ThreeVector zero; 87 fSubStepDynamicParticle = new G4DynamicParti 88 fSubStepTrack = new G4Track(fSubStepDynamicP 89 fSubStep = new G4Step; 90 fSubStepTrack->SetStep(fSubStep); 91 } 92 93 //....oooOO0OOooo........oooOO0OOooo........oo 94 95 G4TransportationWithMsc::~G4TransportationWith 96 { 97 delete fModelManager; 98 delete fParticleChangeForMSC; 99 delete fEmData; 100 delete fParticleChangeForSS; 101 delete fSecondariesSS; 102 103 // fSubStepDynamicParticle is owned and also 104 delete fSubStepTrack; 105 delete fSubStep; 106 } 107 108 //....oooOO0OOooo........oooOO0OOooo........oo 109 110 void G4TransportationWithMsc::AddMscModel(G4VM 111 cons 112 { 113 if (fType != ScatteringType::MultipleScatter 114 G4Exception("G4TransportationWithMsc::AddM 115 "not allowed unless type == Mu 116 } 117 118 fModelManager->AddEmModel(order, mscModel, n 119 mscModel->SetParticleChange(fParticleChangeF 120 } 121 122 //....oooOO0OOooo........oooOO0OOooo........oo 123 124 void G4TransportationWithMsc::AddSSModel(G4VEm 125 { 126 if (fType != ScatteringType::SingleScatterin 127 G4Exception("G4TransportationWithMsc::AddS 128 "not allowed unless type == Si 129 } 130 131 fModelManager->AddEmModel(order, model, null 132 model->SetPolarAngleLimit(0.0); 133 model->SetParticleChange(fParticleChangeForS 134 } 135 136 //....oooOO0OOooo........oooOO0OOooo........oo 137 138 void G4TransportationWithMsc::PreparePhysicsTa 139 { 140 if (nullptr == fFirstParticle) { 141 fFirstParticle = ∂ 142 G4VMultipleScattering* ptr = nullptr; 143 auto emConfigurator = fEmManager->EmConfig 144 emConfigurator->PrepareModels(&part, ptr, 145 } 146 147 if (fFirstParticle == &part) { 148 G4bool master = fEmManager->IsMaster(); 149 G4LossTableBuilder* bld = fEmManager->GetT 150 G4bool baseMat = bld->GetBaseMaterialFlag( 151 const auto* theParameters = G4EmParameters 152 153 if (master) { 154 SetVerboseLevel(theParameters->Verbose() 155 } 156 else { 157 SetVerboseLevel(theParameters->WorkerVer 158 } 159 160 const G4int numberOfModels = fModelManager 161 if (fType == ScatteringType::MultipleScatt 162 for (G4int i = 0; i < numberOfModels; ++ 163 auto msc = static_cast<G4VMscModel*>(f 164 msc->SetPolarAngleLimit(theParameters- 165 G4double emax = std::min(msc->HighEner 166 msc->SetHighEnergyLimit(emax); 167 msc->SetUseBaseMaterials(baseMat); 168 } 169 } 170 else if (fType == ScatteringType::SingleSc 171 if (master) { 172 if (fEmData == nullptr) { 173 fEmData = new G4EmDataHandler(2); 174 } 175 176 fLambdaTable = fEmData->MakeTable(0); 177 bld->InitialiseBaseMaterials(fLambdaTa 178 } 179 } 180 181 fCuts = fModelManager->Initialise(fFirstPa 182 } 183 } 184 185 //....oooOO0OOooo........oooOO0OOooo........oo 186 187 void G4TransportationWithMsc::BuildPhysicsTabl 188 { 189 if (fFirstParticle == &part) { 190 fEmManager->BuildPhysicsTable(fFirstPartic 191 192 if (fEmManager->IsMaster()) { 193 if (fType == ScatteringType::SingleScatt 194 const auto* theParameters = G4EmParame 195 G4LossTableBuilder* bld = fEmManager-> 196 const G4ProductionCutsTable* theCouple 197 G4ProductionCutsTable::GetProduction 198 std::size_t numOfCouples = theCoupleTa 199 200 G4double emin = theParameters->MinKinE 201 G4double emax = theParameters->MaxKinE 202 203 G4double scale = emax / emin; 204 G4int nbin = theParameters->NumberOfBi 205 scale = nbin / G4Log(scale); 206 207 G4int bin = G4lrint(scale * G4Log(emax 208 bin = std::max(bin, 5); 209 210 for (std::size_t i = 0; i < numOfCoupl 211 if (!bld->GetFlag(i)) continue; 212 213 // Create physics vector and fill it 214 const G4MaterialCutsCouple* couple = 215 216 auto* aVector = new G4PhysicsLogVect 217 fModelManager->FillLambdaVector(aVec 218 aVector->FillSecondDerivatives(); 219 G4PhysicsTableHelper::SetPhysicsVect 220 } 221 } 222 } 223 else { 224 const auto masterProcess = static_cast<c 225 226 // Initialisation of models. 227 const G4int numberOfModels = fModelManag 228 if (fType == ScatteringType::MultipleSca 229 for (G4int i = 0; i < numberOfModels; 230 auto msc = static_cast<G4VMscModel*> 231 auto msc0 = static_cast<G4VMscModel* 232 msc->SetCrossSectionTable(msc0->GetC 233 msc->InitialiseLocal(fFirstParticle, 234 } 235 } 236 else if (fType == ScatteringType::Single 237 this->fLambdaTable = masterProcess->fL 238 } 239 } 240 } 241 242 if (!G4EmParameters::Instance()->IsPrintLock 243 G4cout << G4endl; 244 G4cout << GetProcessName() << ": for " << 245 if (fMultipleSteps) { 246 G4cout << " (multipleSteps: 1)"; 247 } 248 G4cout << G4endl; 249 fModelManager->DumpModelList(G4cout, verbo 250 } 251 } 252 253 //....oooOO0OOooo........oooOO0OOooo........oo 254 255 void G4TransportationWithMsc::StartTracking(G4 256 { 257 auto* currParticle = track->GetParticleDefin 258 fIonisation = fEmManager->GetEnergyLossProce 259 260 fSubStepDynamicParticle->SetDefinition(currP 261 262 const G4int numberOfModels = fModelManager-> 263 if (fType == ScatteringType::MultipleScatter 264 for (G4int i = 0; i < numberOfModels; ++i) 265 auto msc = static_cast<G4VMscModel*>(fMo 266 msc->StartTracking(track); 267 msc->SetIonisation(fIonisation, currPart 268 } 269 } 270 271 // Ensure that field propagation state is al 272 G4Transportation::StartTracking(track); 273 } 274 275 //....oooOO0OOooo........oooOO0OOooo........oo 276 277 G4double G4TransportationWithMsc::AlongStepGet 278 279 280 281 282 { 283 *selection = NotCandidateForSelection; 284 285 const G4double physStepLimit = currentMinimu 286 287 switch (fType) { 288 case ScatteringType::MultipleScattering: { 289 // Select the MSC model for the current 290 G4VMscModel* mscModel = nullptr; 291 const G4double ekin = track.GetKineticEn 292 const auto* couple = track.GetMaterialCu 293 const auto* particleDefinition = track.G 294 if (physStepLimit > kGeomMin) { 295 G4double ekinForSelection = ekin; 296 G4double pdgMass = particleDefinition- 297 if (pdgMass > CLHEP::GeV) { 298 ekinForSelection *= proton_mass_c2 / 299 } 300 301 if (ekinForSelection >= kLowestKinEner 302 mscModel = static_cast<G4VMscModel*> 303 fModelManager->SelectModel(ekinFor 304 if (mscModel == nullptr) { 305 G4Exception("G4TransportationWithM 306 "no MSC model found"); 307 } 308 if (!mscModel->IsActive(ekinForSelec 309 mscModel = nullptr; 310 } 311 } 312 } 313 314 // Call the MSC model to potentially lim 315 // geometric path length. 316 if (mscModel != nullptr) { 317 mscModel->SetCurrentCouple(couple); 318 319 // Use the provided track for the firs 320 const G4Track* currentTrackPtr = &trac 321 322 G4double currentSafety = proposedSafet 323 G4double currentEnergy = ekin; 324 325 G4double stepLimitLeft = physStepLimit 326 G4double totalGeometryStepLength = 0, 327 G4bool firstStep = true, continueStepp 328 329 do { 330 G4double gPathLength = stepLimitLeft 331 G4double tPathLength = 332 mscModel->ComputeTruePathLengthLim 333 G4bool mscLimitsStep = (tPathLength 334 if (!fMultipleSteps && mscLimitsStep 335 // MSC limits the step. 336 *selection = CandidateForSelection 337 } 338 339 if (!firstStep) { 340 // Move the navigator to where the 341 fLinearNavigator->LocateGlobalPoin 342 } 343 344 G4GPILSelection transportSelection; 345 G4double geometryStepLength = G4Tran 346 *currentTrackPtr, previousStepSize 347 if (geometryStepLength < gPathLength 348 // Transportation limits the step, 349 *selection = CandidateForSelection 350 continueStepping = false; 351 } 352 if (fTransportEndKineticEnergy != cu 353 // Field propagation changed the e 354 // estimate the continuous energy 355 continueStepping = false; 356 } 357 358 if (firstStep) { 359 proposedSafety = currentSafety; 360 } 361 totalGeometryStepLength += geometryS 362 363 // Sample MSC direction change and d 364 const G4double range = mscModel->Get 365 366 tPathLength = mscModel->ComputeTrueS 367 368 // Protect against wrong t->g->t con 369 tPathLength = std::min(tPathLength, 370 371 totalTruePathLength += tPathLength; 372 if (*selection != CandidateForSelect 373 // If neither MSC nor transportati 374 // distance we want - make sure we 375 continueStepping = false; 376 } 377 else if (tPathLength >= range) { 378 // The particle will stop, exit th 379 continueStepping = false; 380 } 381 else { 382 stepLimitLeft -= tPathLength; 383 } 384 385 // Do not sample scattering at the l 386 if (tPathLength < range && tPathLeng 387 static constexpr G4double minSafet 388 static constexpr G4double sFact = 389 390 // The call to SampleScattering() 391 // direction into fParticleChangeF 392 // 1) Make sure the momentum direc 393 fParticleChangeForMSC->ProposeMome 394 // 2) Call SampleScattering(), whi 395 const G4ThreeVector displacement = 396 mscModel->SampleScattering(fTran 397 // 3) Get the changed direction. 398 fTransportEndMomentumDir = *fParti 399 400 const G4double r2 = displacement.m 401 if (r2 > kMinDisplacement2) { 402 G4bool positionChanged = true; 403 G4double dispR = std::sqrt(r2); 404 G4double postSafety = 405 sFact * fpSafetyHelper->Comput 406 407 // Far away from geometry bounda 408 if (postSafety > 0.0 && dispR <= 409 fTransportEndPosition += displ 410 411 // Near the boundary 412 } 413 else { 414 // displaced point is definite 415 if (dispR < postSafety) { 416 fTransportEndPosition += dis 417 418 // reduced displacement 419 } 420 else if (postSafety > kGeomMin 421 fTransportEndPosition += dis 422 423 // very small postSafety 424 } 425 else { 426 positionChanged = false; 427 } 428 } 429 if (positionChanged) { 430 fpSafetyHelper->ReLocateWithin 431 } 432 } 433 } 434 435 if (continueStepping) { 436 // Update safety according to the 437 if (currentSafety < fEndPointDista 438 currentSafety = 0; 439 } 440 else { 441 currentSafety -= fEndPointDistan 442 } 443 444 // Update the kinetic energy accor 445 currentEnergy = mscModel->GetEnerg 446 447 // From now on, use the track that 448 currentTrackPtr = fSubStepTrack; 449 450 fSubStepDynamicParticle->SetKineti 451 fSubStepDynamicParticle->SetMoment 452 fSubStepTrack->SetPosition(fTransp 453 454 G4StepPoint& subPreStepPoint = *fS 455 subPreStepPoint.SetMaterialCutsCou 456 subPreStepPoint.SetPosition(fTrans 457 subPreStepPoint.SetSafety(currentS 458 subPreStepPoint.SetStepStatus(fAlo 459 } 460 firstStep = false; 461 } while (continueStepping); 462 463 // Note: currentEnergy is only updated 464 // In case field propagation changed t 465 // immediately set to false and curren 466 // initial kinetic energy stored in ek 467 if (currentEnergy != ekin) { 468 // If field propagation didn't chang 469 // did multiple steps, reset the ene 470 // propose to not subtract the energ 471 fTransportEndKineticEnergy = ekin; 472 // Also ask for the range again with 473 // correctly cached in the G4VEnergy 474 // FIXME: Asking for a range should 475 (void)mscModel->GetRange(particleDef 476 } 477 478 fParticleChange.ProposeTrueStepLength( 479 480 // Inform G4Transportation that the mo 481 // to scattering. We do this unconditi 482 // where the last step is done without 483 // the flag, for example when running 484 fMomentumChanged = true; 485 486 return totalGeometryStepLength; 487 } 488 break; 489 } 490 491 case ScatteringType::SingleScattering: { 492 // Select the model for the current kine 493 const G4double ekin = track.GetKineticEn 494 const auto* couple = track.GetMaterialCu 495 const auto* particleDefinition = track.G 496 497 G4double ekinForSelection = ekin; 498 G4double pdgMass = particleDefinition->G 499 if (pdgMass > CLHEP::GeV) { 500 ekinForSelection *= proton_mass_c2 / p 501 } 502 503 G4VEmModel* currentModel = fModelManager 504 if (currentModel == nullptr) { 505 G4Exception("G4TransportationWithMsc:: 506 "no scattering model found 507 } 508 if (!currentModel->IsActive(ekinForSelec 509 currentModel = nullptr; 510 } 511 512 if (currentModel != nullptr) { 513 currentModel->SetCurrentCouple(couple) 514 G4int coupleIndex = couple->GetIndex() 515 516 // Compute mean free path. 517 G4double logEkin = track.GetDynamicPar 518 G4double lambda = ((*fLambdaTable)[cou 519 if (lambda > 0.0) { 520 // Assume that the mean free path an 521 // step, which is a valid approximat 522 G4double meanFreePath = 1.0 / lambda 523 G4double dedx = fIonisation->GetDEDX 524 525 G4double currentSafety = proposedSaf 526 G4double currentEnergy = ekin; 527 528 // Use the provided track for the fi 529 const G4Track* currentTrackPtr = &tr 530 531 G4double stepLimitLeft = physStepLim 532 G4double totalStepLength = 0; 533 G4bool firstStep = true, continueSte 534 535 do { 536 G4double interactionLength = meanF 537 538 G4bool ssLimitsStep = (interaction 539 G4double gPathLength = stepLimitLe 540 if (ssLimitsStep) { 541 if (!fMultipleSteps) { 542 // Scattering limits the step. 543 *selection = CandidateForSelec 544 } 545 gPathLength = interactionLength; 546 } 547 548 if (!firstStep) { 549 // Move the navigator to where t 550 fLinearNavigator->LocateGlobalPo 551 } 552 553 G4GPILSelection transportSelection 554 G4double geometryStepLength = G4Tr 555 *currentTrackPtr, previousStepSi 556 if (geometryStepLength < gPathLeng 557 // Transportation limits the ste 558 *selection = CandidateForSelecti 559 ssLimitsStep = false; 560 continueStepping = false; 561 } 562 if (fTransportEndKineticEnergy != 563 // Field propagation changed the 564 // estimate the continuous energ 565 continueStepping = false; 566 } 567 568 if (firstStep) { 569 proposedSafety = currentSafety; 570 } 571 totalStepLength += geometryStepLen 572 573 if (*selection != CandidateForSele 574 // If neither scattering nor tra 575 // got the distance we want - ma 576 continueStepping = false; 577 } 578 else { 579 stepLimitLeft -= geometryStepLen 580 } 581 582 // Update the kinetic energy accor 583 G4double energyAfterLinearLoss = 584 fTransportEndKineticEnergy - geo 585 586 if (ssLimitsStep) { 587 fSubStepDynamicParticle->SetKine 588 589 // The call to SampleSecondaries 590 // direction into fParticleChang 591 // 1) Set the momentum direction 592 fSubStepDynamicParticle->SetMome 593 // 2) Call SampleSecondaries(), 594 currentModel->SampleSecondaries( 595 596 // 3) Get the changed direction. 597 fTransportEndMomentumDir = fPart 598 599 // Check that the model neither 600 // a local energy deposit becaus 601 // to handle these cases. 602 if (fSecondariesSS->size() > 0) 603 G4Exception("G4TransportationW 604 "scattering model 605 } 606 if (fParticleChangeForSS->GetLoc 607 G4Exception("G4TransportationW 608 "scattering model 609 } 610 } 611 612 if (continueStepping) { 613 // Update safety according to th 614 if (currentSafety < fEndPointDis 615 currentSafety = 0; 616 } 617 else { 618 currentSafety -= fEndPointDist 619 } 620 621 // Update the energy taking cont 622 currentEnergy = energyAfterLinea 623 624 // From now on, use the track th 625 currentTrackPtr = fSubStepTrack; 626 627 fSubStepDynamicParticle->SetKine 628 fSubStepDynamicParticle->SetMome 629 fSubStepTrack->SetPosition(fTran 630 631 G4StepPoint& subPreStepPoint = * 632 subPreStepPoint.SetMaterialCutsC 633 subPreStepPoint.SetPosition(fTra 634 subPreStepPoint.SetSafety(curren 635 subPreStepPoint.SetStepStatus(fA 636 } 637 firstStep = false; 638 } while (continueStepping); 639 640 // Note: currentEnergy is only updat 641 // In case field propagation changed 642 // immediately set to false and curr 643 // initial kinetic energy stored in 644 if (currentEnergy != ekin) { 645 // If field propagation didn't cha 646 // did multiple steps, reset the e 647 // propose to not subtract the ene 648 fTransportEndKineticEnergy = ekin; 649 } 650 651 fParticleChange.ProposeTrueStepLengt 652 653 // Inform G4Transportation that the 654 // to scattering, even if there is n 655 fMomentumChanged = true; 656 657 return totalStepLength; 658 } 659 } 660 break; 661 } 662 } 663 664 // If we get here, no scattering has happene 665 return G4Transportation::AlongStepGetPhysica 666 track, previousStepSize, currentMinimumSte 667 } 668 669 //....oooOO0OOooo........oooOO0OOooo........oo 670