| 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