| 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 //
27 // -------------------------------------------
28 //
29 // GEANT4 Class implementation file
30 //
31 // File name: G4GoudsmitSaundersonMscModel
32 //
33 // Author: Mihaly Novak / (Omrane Kadri
34 //
35 // Creation date: 20.02.2009
36 //
37 // Modifications:
38 // 04.03.2009 V.Ivanchenko cleanup and format
39 // 12.05.2010 O.Kadri: adding Qn1 and Qn12 as
40 // 18.05.2015 M. Novak provide PRELIMINARY ver
41 // This class has been revised and
42 // A new version of Kawrakow-Bielaj
43 // based on the screened Rutherford
44 // electrons/positrons has been int
45 // angular distributions over a 2D
46 // and the CDFs are now stored in a
47 // together with the corresponding
48 // These angular distributions are
49 // G4GoudsmitSaundersonTable class
50 // it was no, single, few or multip
51 // angular deflection (i.e. cos(the
52 // Two screening options are provid
53 // - if fIsUsePWATotalXsecData=TRU
54 // was called before initialisat
55 // determined such that the firs
56 // computed according to the scr
57 // scattering will reproduce the
58 // and first transport mean free
59 // - if fIsUsePWATotalXsecData=FAL
60 // SetOptionPWAScreening(FALSE)
61 // screening parameter value A i
62 // formula (by using material de
63 // precomputed for each material
64 // G4GoudsmitSaundersonTable) [3
65 // Elastic and first trasport mean
66 // The new version is self-consiste
67 // robust and accurate compared to
68 // Spin effects as well as a more a
69 // computations of Lewis moments wi
70 // 02.09.2015 M. Novak: first version of new s
71 // fUseSafetyPlus corresponds to Ur
72 // fUseDistanceToBoundary correspon
73 // fUseSafety corresponds to EGSnr
74 // Range factor can be significantl
75 // 23.08.2017 M. Novak: added corrections to a
76 // It can be activated by setting t
77 // before initialization using the
78 // The fMottCorrection member is re
79 // correction (rejection) functions
80 // Goudsmit-Saunderson agnular dist
81 // effects and screening correction
82 // GS angular distributions is: DCS
83 // # DCS_{SR} is the relativisti
84 // solution of the Klein-Gordo
85 // scattering of spinless e- o
86 // note: the default (without
87 // are based on this DCS_{SR}
88 // # DCS_{R} is the Rutherford D
89 // screening
90 // # DCS_{Mott} is the Mott DCS
91 // Coulomb potential i.e. scat
92 // point-like unscreened Coulo
93 // # moreover, the screening par
94 // the DCS_{cor} with this cor
95 // transport cross sections ob
96 // (i.e. from elsepa [4])
97 // Unlike the default GS, the Mott-
98 // (different for e- and e+ <= the
99 // (Z and material) dependent.
100 // 27.10.2017 M. Novak:
101 // - Mott-correction flag is set no
102 // - new form of PWA correction to
103 // - changed step limit flag conven
104 // # fUseSafety corresponds to U
105 // # fUseDistanceToBoundary corr
106 // # fUseSafetyPlus corresponds
107 // 02.02.2018 M. Novak: implemented CrossSecti
108 //
109 // Class description:
110 // Kawrakow-Bielajew Goudsmit-Saunderson MSC
111 // for elastic scattering of e-/e+. Option,
112 // also available now (SetOptionMottCorrecti
113 // algorithm (UseSafety) is available beyond
114 // and true to geomerty and geometry to true
115 // from the Urban model[5]. The most accurat
116 // UseSafetyPlus MSC step limit with Mott-co
117 // are expected to be set through the G4EmPa
118 // # G4EmParameters::Instance()->SetMscStep
119 // # G4EmParameters::Instance()->SetUseMott
120 //
121 //
122 // References:
123 // [1] A.F.Bielajew, NIMB 111 (1996) 195-208
124 // [2] I.Kawrakow, A.F.Bielajew, NIMB 134(19
125 // [3] I.Kawrakow, E.Mainegra-Hing, D.W.O.Ro
126 // Report PIRS-701 (2013)
127 // [4] F.Salvat, A.Jablonski, C.J. Powell, C
128 // [5] L.Urban, Preprint CERN-OPEN-2006-077
129 //
130 // -------------------------------------------
131
132
133 #include "G4GoudsmitSaundersonMscModel.hh"
134
135 #include "G4GoudsmitSaundersonTable.hh"
136 #include "G4GSPWACorrections.hh"
137
138 #include "G4PhysicalConstants.hh"
139 #include "G4SystemOfUnits.hh"
140
141 #include "G4ParticleChangeForMSC.hh"
142 #include "G4DynamicParticle.hh"
143 #include "G4Electron.hh"
144 #include "G4Positron.hh"
145
146 #include "G4LossTableManager.hh"
147 #include "G4EmParameters.hh"
148 #include "G4Track.hh"
149 #include "G4PhysicsTable.hh"
150 #include "Randomize.hh"
151 #include "G4Log.hh"
152 #include "G4Exp.hh"
153 #include "G4Pow.hh"
154 #include <fstream>
155
156
157 // set accurate energy loss and dispalcement s
158 G4bool G4GoudsmitSaundersonMscModel::gIsUseAcc
159 // set the usual optimization to be always act
160 G4bool G4GoudsmitSaundersonMscModel::gIsOptimi
161
162
163 G4GoudsmitSaundersonMscModel::G4GoudsmitSaunde
164 : G4VMscModel(nam) {
165 charge = 0;
166 currentMaterialIndex = -1;
167 //
168 fr = 0.1;
169 rangeinit = 1.e+21;
170 geombig = 1.e+50*mm;
171 geomlimit = geombig;
172 tgeom = geombig;
173 tlimit = 1.e+10*mm;
174 presafety = 0.*mm;
175 //
176 particle = nullptr;
177 theManager = G4LossTableManager:
178 firstStep = true;
179 currentKinEnergy = 0.0;
180 currentRange = 0.0;
181 //
182 tlimitminfix2 = 1.*nm;
183 tausmall = 1.e-16;
184 mass = electron_mass_c2;
185 taulim = 1.e-6;
186 //
187 currentCouple = nullptr;
188 fParticleChange = nullptr;
189 //
190 fZeff = 1.;
191 //
192 par1 = 0.;
193 par2 = 0.;
194 par3 = 0.;
195 //
196 // Moliere screeing parameter will be used a
197 // appalied to the integrated quantities (sc
198 // and second moments) derived from the corr
199 // this PWA correction is ignored if Mott-co
200 // Mott-correction contains all these correc
201 fIsUsePWACorrection = true;
202 //
203 fIsUseMottCorrection = false;
204 //
205 fLambda0 = 0.0; // elastic mea
206 fLambda1 = 0.0; // first trans
207 fScrA = 0.0; // screening p
208 fG1 = 0.0; // first trans
209 //
210 fMCtoScrA = 1.0;
211 fMCtoQ1 = 1.0;
212 fMCtoG2PerG1 = 1.0;
213 //
214 fTheTrueStepLenght = 0.;
215 fTheTransportDistance = 0.;
216 fTheZPathLenght = 0.;
217 //
218 fTheDisplacementVector.set(0.,0.,0.);
219 fTheNewDirection.set(0.,0.,1.);
220 //
221 fIsEverythingWasDone = false;
222 fIsMultipleSacettring = false;
223 fIsSingleScattering = false;
224 fIsEndedUpOnBoundary = false;
225 fIsNoScatteringInMSC = false;
226 fIsNoDisplace = false;
227 fIsInsideSkin = false;
228 fIsWasOnBoundary = false;
229 fIsFirstRealStep = false;
230 rndmEngineMod = G4Random::getTheEng
231 //
232 fGSTable = nullptr;
233 fPWACorrection = nullptr;
234 }
235
236
237 G4GoudsmitSaundersonMscModel::~G4GoudsmitSaund
238 if (IsMaster()) {
239 if (fGSTable) {
240 delete fGSTable;
241 fGSTable = nullptr;
242 }
243 if (fPWACorrection) {
244 delete fPWACorrection;
245 fPWACorrection = nullptr;
246 }
247 }
248 }
249
250
251 void G4GoudsmitSaundersonMscModel::Initialise(
252 SetParticle(p);
253 InitialiseParameters(p);
254 // -create GoudsmitSaundersonTable and init
255 // Mott-correction was required
256 if (IsMaster()) {
257 // get the Mott-correction flag from EmPar
258 if (G4EmParameters::Instance()->UseMottCor
259 fIsUseMottCorrection = true;
260 }
261 // Mott-correction includes other way of P
262 // when Mott-correction is activated by th
263 if (fIsUseMottCorrection) {
264 fIsUsePWACorrection = false;
265 }
266 // clear GS-table
267 if (fGSTable) {
268 delete fGSTable;
269 fGSTable = nullptr;
270 }
271 // clear PWA corrections table if any
272 if (fPWACorrection) {
273 delete fPWACorrection;
274 fPWACorrection = nullptr;
275 }
276 // create GS-table
277 G4bool isElectron = true;
278 if (p->GetPDGCharge()>0.) {
279 isElectron = false;
280 }
281 fGSTable = new G4GoudsmitSaundersonTable(i
282 // G4GSTable will be initialised:
283 // - Screened-Rutherford DCS based GS angu
284 // - Mott-correction will be initialised i
285 fGSTable->SetOptionMottCorrection(fIsUseMo
286 // - set PWA correction (correction to int
287 fGSTable->SetOptionPWACorrection(fIsUsePWA
288 // init
289 fGSTable->Initialise(LowEnergyLimit(),High
290 // create PWA corrections table if it was
291 if (fIsUsePWACorrection) {
292 fPWACorrection = new G4GSPWACorrections(
293 fPWACorrection->Initialise();
294 }
295 }
296 fParticleChange = GetParticleChangeForMSC(p)
297 }
298
299
300 void G4GoudsmitSaundersonMscModel::InitialiseL
301 fGSTable = static_cast<G4Goud
302 fIsUseMottCorrection = static_cast<G4Goud
303 fIsUsePWACorrection = static_cast<G4Goud
304 fPWACorrection = static_cast<G4Goud
305 }
306
307
308 // computes macroscopic first transport cross
309 G4double G4GoudsmitSaundersonMscModel::CrossSe
310 const
311 G4dou
312 G4dou
313 G4dou
314 G4double xsecTr1 = 0.; // cross section per
315 //
316 fLambda0 = 0.0; // elastic mean free path
317 fLambda1 = 0.0; // first transport mean free
318 fScrA = 0.0; // screening parameter
319 fG1 = 0.0; // first transport coef.
320 // use Moliere's screening (with Mott-corret
321 G4double efEnergy = std::max(kineticEnergy,
322 // total mometum square
323 G4double pt2 = efEnergy*(efEnergy+2.0*el
324 // beta square
325 G4double beta2 = pt2/(pt2+electron_mass_c2
326 // current material index
327 G4int matindx = (G4int)mat->GetIndex();
328 // Moliere's b_c
329 G4double bc = fGSTable->GetMoliereBc(ma
330 // get the Mott-correcton factors if Mott-co
331 fMCtoScrA = 1.0;
332 fMCtoQ1 = 1.0;
333 fMCtoG2PerG1 = 1.0;
334 G4double scpCor = 1.0;
335 if (fIsUseMottCorrection) {
336 fGSTable->GetMottCorrectionFactors(G4Log(e
337 // ! no scattering power correction since
338 // scpCor = fGSTable->ComputeScatteringPow
339 } else if (fIsUsePWACorrection) {
340 fPWACorrection->GetPWACorrectionFactors(G4
341 // scpCor = fGSTable->ComputeScatteringPow
342 }
343 // screening parameter:
344 // - if Mott-corretioncorrection: the Screen
345 // screening parameter gives back the (els
346 // - if PWA correction: he Screened-Rutherfo
347 // gives back the (elsepa) PWA first trans
348 fScrA = fGSTable->GetMoliereXc2(matindx)/
349 // elastic mean free path in Geant4 internal
350 // (if Mott-corretion: the corrected screeni
351 // corrected with the screening parameter co
352 fLambda0 = beta2*(1.+fScrA)*fMCtoScrA/bc/scp
353 // first transport coefficient (if Mott-corr
354 // consistent with the one used during the p
355 fG1 = 2.0*fScrA*((1.0+fScrA)*G4Log(1.0/
356 // first transport mean free path
357 fLambda1 = fLambda0/fG1;
358 xsecTr1 = 1./fLambda1;
359 return xsecTr1;
360 }
361
362
363 // gives back the first transport mean free pa
364 G4double
365 G4GoudsmitSaundersonMscModel::GetTransportMean
366
367 // kinetic energy is assumed to be in Geant4
368 G4double efEnergy = kineticEnergy;
369 //
370 const G4Material* mat = currentCouple->GetM
371 //
372 fLambda0 = 0.0; // elastic mean free path
373 fLambda1 = 0.0; // first transport mean free
374 fScrA = 0.0; // screening parameter
375 fG1 = 0.0; // first transport coef.
376
377 // use Moliere's screening (with Mott-corret
378 if (efEnergy<10.*CLHEP::eV) efEnergy = 10.*
379 // total mometum square
380 G4double pt2 = efEnergy*(efEnergy+2.0*el
381 // beta square
382 G4double beta2 = pt2/(pt2+electron_mass_c2
383 // current material index
384 G4int matindx = (G4int)mat->GetIndex();
385 // Moliere's b_c
386 G4double bc = fGSTable->GetMoliereBc(ma
387 // get the Mott-correcton factors if Mott-co
388 fMCtoScrA = 1.0;
389 fMCtoQ1 = 1.0;
390 fMCtoG2PerG1 = 1.0;
391 G4double scpCor = 1.0;
392 if (fIsUseMottCorrection) {
393 fGSTable->GetMottCorrectionFactors(G4Log(e
394 scpCor = fGSTable->ComputeScatteringPowerC
395 } else if (fIsUsePWACorrection) {
396 fPWACorrection->GetPWACorrectionFactors(G4
397 // scpCor = fGSTable->ComputeScatteringPow
398 }
399 // screening parameter:
400 // - if Mott-corretioncorrection: the Screen
401 // screening parameter gives back the (els
402 // - if PWA correction: he Screened-Rutherfo
403 // gives back the (elsepa) PWA first trans
404 fScrA = fGSTable->GetMoliereXc2(matindx)/
405 // elastic mean free path in Geant4 internal
406 // (if Mott-corretion: the corrected screeni
407 // corrected with the screening parameter co
408 fLambda0 = beta2*(1.+fScrA)*fMCtoScrA/bc/scp
409 // first transport coefficient (if Mott-corr
410 // consistent with the one used during the p
411 fG1 = 2.0*fScrA*((1.0+fScrA)*G4Log(1.0/
412 // first transport mean free path
413 fLambda1 = fLambda0/fG1;
414
415 return fLambda1;
416 }
417
418
419 G4double
420 G4GoudsmitSaundersonMscModel::GetTransportMean
421
422 // kinetic energy is assumed to be in Geant4
423 G4double efEnergy = kineticEnergy;
424 //
425 const G4Material* mat = currentCouple->GetM
426 //
427 G4double lambda0 = 0.0; // elastc mean free
428 G4double lambda1 = 0.0; // first transport m
429 G4double scrA = 0.0; // screening paramet
430 G4double g1 = 0.0; // first transport m
431
432 // use Moliere's screening (with Mott-corret
433 if (efEnergy<10.*CLHEP::eV) efEnergy = 10.*
434 // total mometum square in Geant4 internal e
435 G4double pt2 = efEnergy*(efEnergy+2.0*el
436 G4double beta2 = pt2/(pt2+electron_mass_c2
437 G4int matindx = (G4int)mat->GetIndex();
438 G4double bc = fGSTable->GetMoliereBc(ma
439 // get the Mott-correcton factors if Mott-co
440 G4double mctoScrA = 1.0;
441 G4double mctoQ1 = 1.0;
442 G4double mctoG2PerG1 = 1.0;
443 G4double scpCor = 1.0;
444 if (fIsUseMottCorrection) {
445 fGSTable->GetMottCorrectionFactors(G4Log(e
446 scpCor = fGSTable->ComputeScatteringPowerC
447 } else if (fIsUsePWACorrection) {
448 fPWACorrection->GetPWACorrectionFactors(G4
449 // scpCor = fGSTable->ComputeScatteringPow
450 }
451 scrA = fGSTable->GetMoliereXc2(matindx)/(
452 // total elastic mean free path in Geant4 in
453 lambda0 = beta2*(1.+scrA)*mctoScrA/bc/scpCor
454 g1 = 2.0*scrA*((1.0+scrA)*G4Log(1.0/scr
455 lambda1 = lambda0/g1;
456
457 return lambda1;
458 }
459
460
461 void G4GoudsmitSaundersonMscModel::StartTracki
462 SetParticle(track->GetDynamicParticle()->Get
463 firstStep = true;
464 tlimit = tgeom = rangeinit = geombig;
465 rangeinit = 1.e+21;
466 }
467
468
469 G4double G4GoudsmitSaundersonMscModel::Compute
470
471 G4double skindepth = 0.;
472 //
473 const G4DynamicParticle* dp = track.GetDynam
474 G4StepPoint* sp = track.GetStep(
475 G4StepStatus stepStatus = sp->GetStepSta
476 currentCouple = track.GetMater
477 SetCurrentCouple(currentCouple);
478 currentMaterialIndex = (G4int)current
479 currentKinEnergy = dp->GetKinetic
480 currentRange = GetRange(parti
481 dp->G
482 // elastic and first transport mfp, screenin
483 // (Mott-correction will be used if it was r
484 fLambda1 = GetTransportMeanFreePath(particle
485 // Set initial values:
486 // : lengths are initialised to currentMini
487 // step length from all other physics
488 fTheTrueStepLenght = currentMinimalStep;
489 fTheTransportDistance = currentMinimalStep;
490 fTheZPathLenght = currentMinimalStep;
491 fTheDisplacementVector.set(0.,0.,0.);
492 fTheNewDirection.set(0.,0.,1.);
493
494 // Can everything be done in the step limit
495 fIsEverythingWasDone = false;
496 // Multiple scattering needs to be sample ?
497 fIsMultipleSacettring = false;
498 // Single scattering needs to be sample ?
499 fIsSingleScattering = false;
500 // Was zero deflection in multiple scatterin
501 fIsNoScatteringInMSC = false;
502 // Do not care about displacement in MSC sam
503 // ( used only in the case of gIsOptimizatio
504 fIsNoDisplace = false;
505 // get pre-step point safety
506 presafety = sp->GetSafety();
507 //
508 fZeff = currentCouple->GetMaterial()->GetIon
509 // distance will take into account max-fluct
510 G4double distance = currentRange;
511 distance *= (1.20-fZeff*(1.62e-2-9.22e-5*fZe
512 //
513 // Possible optimization : if the distance i
514 // particle will never leave this volume ->
515 // as the effect of multiple elastic scatter
516 // Important : this optimization can cause p
517 // in a bigger volume since MSC won't be don
518 // distance < safety so don't use optimized-
519 if (gIsOptimizationOn && (distance<presafety
520 // Indicate that we need to do MSC after
521 fIsMultipleSacettring = true;
522 fIsNoDisplace = true;
523 } else if (steppingAlgorithm==fUseDistanceTo
524 //Compute geomlimit (and presafety) :
525 // - geomlimit will be:
526 // == the straight line distance to the
527 // longer than that
528 // == a big value [geombig = 1.e50*mm]
529 // the straight line distance to the
530 // - presafety will be updated as well
531 // So the particle can travell 'gemlimit'
532 // line!) in its current direction:
533 // (1) before reaching a boundary (geomli
534 // (2) before reaching its current range
535 geomlimit = ComputeGeomLimit(track, presaf
536 // Record that the particle is on a bounda
537 if ( (stepStatus==fGeomBoundary) || (stepS
538 fIsWasOnBoundary = true;
539 }
540 // Set skin depth = skin x elastic_mean_fr
541 skindepth = skin*fLambda0;
542 // Init the flag that indicates that the p
543 // distance from a boundary
544 fIsInsideSkin = false;
545 // Check if we can try Single Scattering b
546 // distance from/to a boundary OR the curr
547 // shorter than skindepth. NOTICE: the lat
548 // because the MSC angular sampling is fin
549 // faster to try single scattering in case
550 if ((stepStatus==fGeomBoundary) || (presaf
551 // check if we are within skindepth dist
552 if ((stepStatus == fGeomBoundary) || (pr
553 fIsInsideSkin = true;
554 fIsWasOnBoundary = true;
555 }
556 //Try single scattering:
557 // - sample distance to next single scat
558 // - compare to current minimum length
559 // == if sslimit is the shorter:
560 // - set the step length to ssl
561 // - indicate that single scatt
562 // == else : nothing to do
563 //- in both cases, the step length was v
564 // true path length are the same
565 G4double sslimit = -1.*fLambda0*G4Log(G4
566 // compare to current minimum step lengt
567 if (sslimit<fTheTrueStepLenght) {
568 fTheTrueStepLenght = sslimit;
569 fIsSingleScattering = true;
570 }
571 // short step -> true step length equal
572 fTheZPathLenght = fTheTrueStepLeng
573 // Set taht everything is done in step-l
574 // We will check if we need to perform t
575 // sampling i.e. if single elastic scatt
576 fIsEverythingWasDone = true;
577 } else {
578 // After checking we know that we cannot
579 // need to make an MSC step
580 // Indicate that we need to make and MSC
581 // do it now i.e. if presafety>final_tru
582 // fIsEverythingWasDone = false which in
583 // MSC after transportation.
584 fIsMultipleSacettring = true;
585 // Init the first-real-step falg: it wil
586 // non-single scattering step in this vo
587 fIsFirstRealStep = false;
588 // If previously the partcile was on bou
589 // well. When it is not within skin anym
590 // so we make the first real MSC step wi
591 if (fIsWasOnBoundary && !fIsInsideSkin)
592 // reset the 'was on boundary' indicat
593 fIsWasOnBoundary = false;
594 fIsFirstRealStep = true;
595 }
596 // If this is the first-real msc step (t
597 // skin) or this is the first step with
598 // primary):
599 // - set the initial range that will b
600 // (only in this volume, because aft
601 // first-real MSC step we will reset
602 // - don't let the partcile to cross th
603 if (firstStep || fIsFirstRealStep || ran
604 rangeinit = currentRange;
605 // If geomlimit < geombig than the par
606 // along its initial direction before
607 // Otherwise we can be sure that the p
608 // before reaching the boundary along
609 // geometrical limit appalied. [Howeve
610 // first or the first-real MSC step. A
611 // MSC step the direction will change
612 // But we will try to end up within sk
613 // the actual value of geomlimit(See l
614 // boundary).]
615 if (geomlimit<geombig) {
616 // transfrom straight line distance
617 // length based on the mean values (
618 // first-transport mean free path i.
619 if ((1.-geomlimit/fLambda1)> 0.) {
620 geomlimit = -fLambda1*G4Log(1.-geo
621 }
622 // the 2-different case that could l
623 if (firstStep) {
624 tgeom = 2.*geomlimit/facgeom;
625 } else {
626 tgeom = geomlimit/facgeom;
627 }
628 } else {
629 tgeom = geombig;
630 }
631 }
632 // True step length limit from range fac
633 // range is used that was set at the fir
634 // in this volume with this particle.
635 tlimit = facrange*rangeinit;
636 // Take the minimum of the true step len
637 // geometrical constraint or range-facto
638 tlimit = std::min(tlimit,tgeom);
639 // Step reduction close to boundary: we
640 // from the boundary ( Notice: in case o
641 // because geomlimit is the straigth lin
642 // the currect direction (if geomlimit<g
643 // change the initial direction. So te p
644 // before in a different direction. Howe
645 // path length to this (straight line) l
646 // transport distance (straight line) wi
647 // geomlimit-0.999*skindepth after the c
648 if (geomlimit<geombig) {
649 tlimit = std::min(tlimit, geomlimit-0.
650 }
651 // randomize 1st step or 1st 'normal' st
652 if (firstStep || fIsFirstRealStep) {
653 fTheTrueStepLenght = std::min(fTheTrue
654 } else {
655 fTheTrueStepLenght = std::min(fTheTrue
656 }
657 }
658 } else if (steppingAlgorithm==fUseSafetyPlus
659 presafety = ComputeSafety(sp->GetPosition
660 geomlimit = presafety;
661 // Set skin depth = skin x elastic_mean_fr
662 skindepth = skin*fLambda0;
663 // Check if we can try Single Scattering b
664 // distance from/to a boundary OR the curr
665 // shorter than skindepth. NOTICE: the lat
666 // because the MSC angular sampling is fin
667 // faster to try single scattering in case
668 if ((stepStatus==fGeomBoundary) || (presaf
669 //Try single scattering:
670 // - sample distance to next single scat
671 // - compare to current minimum length
672 // == if sslimit is the shorter:
673 // - set the step length to ssl
674 // - indicate that single scatt
675 // == else : nothing to do
676 //- in both cases, the step length was v
677 // true path length are the same
678 G4double sslimit = -1.*fLambda0*G4Log(G4
679 // compare to current minimum step lengt
680 if (sslimit<fTheTrueStepLenght) {
681 fTheTrueStepLenght = sslimit;
682 fIsSingleScattering = true;
683 }
684 // short step -> true step length equal
685 fTheZPathLenght = fTheTrueStepLeng
686 // Set taht everything is done in step-l
687 // We will check if we need to perform t
688 // sampling i.e. if single elastic scatt
689 fIsEverythingWasDone = true;
690 } else {
691 // After checking we know that we cannot
692 // need to make an MSC step
693 // Indicate that we need to make and MSC
694 fIsMultipleSacettring = true;
695 fIsEverythingWasDone = true;
696 // limit from range factor
697 fTheTrueStepLenght = std::min(fTheTru
698 // never let the particle go further tha
699 // if we are here we are out of the skin
700 if (fTheTrueStepLenght>presafety) {
701 fTheTrueStepLenght = std::min(fTheTrue
702 }
703 // make sure that we are still within th
704 // i.e. true step length is not longer t
705 // We schould take into account energy l
706 // step length as well. So let it 0.5 x
707 fTheTrueStepLenght = std::min(fTheTrueSt
708 }
709 } else {
710 // This is the default stepping algorithm:
711 // accurate that corresponds to fUseSafety
712 // model can handle any short steps so we
713 //
714 // NO single scattering in case of skin or
715 // model will be single or even no scatter
716 // compared to the elastic mean free path.
717 //
718 // indicate that MSC needs to be done (alw
719 fIsMultipleSacettring = true;
720 if (stepStatus!=fGeomBoundary) {
721 presafety = ComputeSafety(sp->GetPositio
722 }
723 // Far from boundary-> in optimized mode d
724 if ((distance<presafety) && (gIsOptimizati
725 fIsNoDisplace = true;
726 } else {
727 // Urban like
728 if (firstStep || (stepStatus==fGeomBound
729 rangeinit = currentRange;
730 fr = facrange;
731 // We don't use this: we won't converge to the
732 // decreasing range-factor.
733 // rangeinit = std::max(rangeinit
734 // if(fLambda1 > lambdalimit) {
735 // fr *= (0.75+0.25*fLambda1/la
736 // }
737
738 }
739 //step limit
740 tlimit = std::max(fr*rangeinit, facsafet
741 // first step randomization
742 if (firstStep || stepStatus==fGeomBounda
743 fTheTrueStepLenght = std::min(fTheTrue
744 } else {
745 fTheTrueStepLenght = std::min(fTheTrue
746 }
747 }
748 }
749 //
750 // unset first-step
751 firstStep =false;
752 // performe single scattering, multiple scat
753 if (fIsEverythingWasDone) {
754 if (fIsSingleScattering) {
755 // sample single scattering
756 //G4double ekin = 0.5*(currentKinEnerg
757 G4double lekin = G4Log(currentKinEnergy
758 G4double pt2 = currentKinEnergy*(curr
759 G4double beta2 = pt2/(pt2+CLHEP::electr
760 G4double cost = fGSTable->SingleScatte
761 // protection
762 if (cost<-1.) cost = -1.;
763 if (cost> 1.) cost = 1.;
764 // compute sint
765 G4double dum = 1.-cost;
766 G4double sint = std::sqrt(dum*(2.-dum)
767 G4double phi = CLHEP::twopi*G4Uniform
768 G4double sinPhi = std::sin(phi);
769 G4double cosPhi = std::cos(phi);
770 fTheNewDirection.set(sint*cosPhi,sint*si
771 } else if (fIsMultipleSacettring) {
772 // sample multiple scattering
773 SampleMSC(); // fTheZPathLenght, fTheDis
774 } // and if single scattering but it was l
775 } //else { do nothing here but after transpo
776 //
777 return ConvertTrueToGeom(fTheTrueStepLenght,
778 }
779
780
781 G4double G4GoudsmitSaundersonMscModel::Compute
782 // convert true ->geom
783 // It is called from the step limitation Com
784 // !fIsEverythingWasDone but protect:
785 par1 = -1.;
786 par2 = par3 = 0.;
787 // if fIsEverythingWasDone = TRUE => fTheZP
788 // so return with the already known value
789 // Otherwise:
790 if (!fIsEverythingWasDone) {
791 // this correction needed to run MSC with
792 // and makes no harm for a normal run
793 fTheTrueStepLenght = std::min(fTheTrueStep
794 // do the true -> geom transformation
795 fTheZPathLenght = fTheTrueStepLenght;
796 // z = t for very small true-path-length
797 if (fTheTrueStepLenght<tlimitminfix2) {
798 return fTheZPathLenght;
799 }
800 G4double tau = fTheTrueStepLenght/fLambda1
801 if (tau<=tausmall) {
802 fTheZPathLenght = std::min(fTheTrueStepL
803 } else if (fTheTrueStepLenght<currentRang
804 if (tau<taulim) fTheZPathLenght = fTheTr
805 else fTheZPathLenght = fLambd
806 } else if (currentKinEnergy<mass || fTheTr
807 par1 = 1./currentRange ; // alpha =1
808 par2 = 1./(par1*fLambda1) ; // 1/(alpha
809 par3 = 1.+par2 ; // 1+1/
810 if (fTheTrueStepLenght<currentRange) {
811 fTheZPathLenght = 1./(par1*par3) * (1.
812 } else {
813 fTheZPathLenght = 1./(par1*par3);
814 }
815 } else {
816 G4double rfin = std::max(currentRange
817 G4double T1 = GetEnergy(particle,rf
818 G4double lambda1 = GetTransportMeanFreeP
819 //
820 par1 = (fLambda1-lambda1)/(fLambda1*fThe
821 par2 = 1./(par1*fLambda1);
822 par3 = 1.+par2 ;
823 G4Pow *g4calc = G4Pow::GetInstance();
824 fTheZPathLenght = 1./(par1*par3) * (1.-g
825 }
826 }
827 fTheZPathLenght = std::min(fTheZPathLenght,
828 //
829 return fTheZPathLenght;
830 }
831
832
833 G4double G4GoudsmitSaundersonMscModel::Compute
834 // init
835 fIsEndedUpOnBoundary = false;
836 // step defined other than transportation
837 if (geomStepLength==fTheZPathLenght) {
838 return fTheTrueStepLenght;
839 }
840 // else ::
841 // - set the flag that transportation was th
842 // - convert geom -> true by using the mean
843 fIsEndedUpOnBoundary = true; // OR LAST STEP
844 fTheZPathLenght = geomStepLength;
845 // was a short single scattering step
846 if (fIsEverythingWasDone && !fIsMultipleSace
847 fTheTrueStepLenght = geomStepLength;
848 return fTheTrueStepLenght;
849 }
850 // t = z for very small step
851 if (geomStepLength<tlimitminfix2) {
852 fTheTrueStepLenght = geomStepLength;
853 // recalculation
854 } else {
855 G4double tlength = geomStepLength;
856 if (geomStepLength>fLambda1*tausmall) {
857 if (par1< 0.) {
858 tlength = -fLambda1*G4Log(1.-geomStepL
859 } else {
860 if (par1*par3*geomStepLength<1.) {
861 G4Pow *g4calc = G4Pow::GetInstance()
862 tlength = (1.-g4calc->powA( 1.-par1*
863 } else {
864 tlength = currentRange;
865 }
866 }
867 if (tlength<geomStepLength || tlength>fT
868 tlength = geomStepLength;
869 }
870 }
871 fTheTrueStepLenght = tlength;
872 }
873 //
874 return fTheTrueStepLenght;
875 }
876
877 G4ThreeVector&
878 G4GoudsmitSaundersonMscModel::SampleScattering
879 if (steppingAlgorithm==fUseDistanceToBoundar
880 // single scattering was and scattering ha
881 fTheNewDirection.rotateUz(oldDirection);
882 fParticleChange->ProposeMomentumDirection(
883 return fTheDisplacementVector;
884 } else if (steppingAlgorithm==fUseSafetyPlus
885 if (fIsEndedUpOnBoundary) { // do nothing
886 return fTheDisplacementVector;
887 } else if (fIsEverythingWasDone) { // evry
888 // check single scattering and see if it
889 if (fIsSingleScattering) {
890 fTheNewDirection.rotateUz(oldDirection
891 fParticleChange->ProposeMomentumDirect
892 return fTheDisplacementVector;
893 }
894 // check if multiple scattering happened
895 if (fIsMultipleSacettring && !fIsNoScatt
896 fTheNewDirection.rotateUz(oldDirect
897 fTheDisplacementVector.rotateUz(old
898 fParticleChange->ProposeMomentumDir
899 }
900 // The only thing that could happen if w
901 // is that single scattering was tried
902 // So no displacement and no scattering
903 return fTheDisplacementVector;
904 }
905 //
906 // The only thing that could still happen
907 // optimization branch: so sample MSC angl
908 }
909 //else MSC needs to be done here
910 SampleMSC();
911 if (!fIsNoScatteringInMSC) {
912 fTheNewDirection.rotateUz(oldDirection);
913 fParticleChange->ProposeMomentumDirection(
914 if (!fIsNoDisplace) {
915 fTheDisplacementVector.rotateUz(oldDirec
916 }
917 }
918 //
919 return fTheDisplacementVector;
920 }
921
922
923 void G4GoudsmitSaundersonMscModel::SampleMSC()
924 fIsNoScatteringInMSC = false;
925 // kinetic energy is assumed to be in Geant4
926 G4double kineticEnergy = currentKinEnergy;
927 //
928 // Energy loss correction: 2 version
929 G4double eloss = 0.0;
930 // if (fTheTrueStepLenght > currentRange*dtrl
931 eloss = kineticEnergy - GetEnergy(particle,c
932 // } else {
933 // eloss = fTheTrueStepLenght*GetDEDX(parti
934 // }
935
936 G4double tau = 0.;// = kineticEnergy/ele
937 G4double tau2 = 0.;// = tau*tau;
938 G4double eps0 = 0.;// = eloss/kineticEnerg
939 G4double epsm = 0.;// = eloss/kineticEnerg
940
941 // - init.
942 G4double efEnergy = kineticEnergy;
943 G4double efStep = fTheTrueStepLenght;
944
945 G4double kineticEnergy0 = kineticEnergy;
946 if (gIsUseAccurate) { // - use accurate e
947 kineticEnergy -= 0.5*eloss; // mean ener
948 // other parameters for energy loss correc
949 tau = kineticEnergy/electron_m
950 tau2 = tau*tau;
951 eps0 = eloss/kineticEnergy0; //
952 epsm = eloss/kineticEnergy; //
953
954 efEnergy = kineticEnergy * (1.-epsm
955 G4double dum = 0.166666*(4.+tau*(6.+tau
956 efStep = fTheTrueStepLenght*(1.-d
957 } else { // - t
958 kineticEnergy -= 0.5*eloss; // mean ener
959 efEnergy = kineticEnergy;
960 G4double factor = 1./(1.+0.9784671*kinetic
961 eps0 = eloss/kineticEnergy0;
962 epsm = eps0/(1.-0.5*eps0);
963 G4double temp = 0.3*(1 -factor*(1.-0.333
964 efStep = fTheTrueStepLenght*(1.+t
965 }
966 //
967 // compute elastic mfp, first transport mfp,
968 // if it was requested by the user)
969 fLambda1 = GetTransportMeanFreePath(particle
970 // s/lambda_el
971 G4double lambdan=0.;
972 if (fLambda0>0.0) {
973 lambdan=efStep/fLambda0;
974 }
975 if (lambdan<=1.0e-12) {
976 if (fIsEverythingWasDone) {
977 fTheZPathLenght = fTheTrueStepLenght;
978 }
979 fIsNoScatteringInMSC = true;
980 return;
981 }
982 // first moment: 2.* lambdan *scrA*((1.+scrA
983 G4double Qn1 = lambdan *fG1;
984 // sample scattering angles
985 // new direction, relative to the orriginal
986 G4double cosTheta1 = 1.0, sinTheta1 = 0.0, c
987 G4double cosPhi1 = 1.0, sinPhi1 = 0.0, c
988 G4double uss = 0.0, vss = 0.0, w
989 G4double x_coord = 0.0, y_coord = 0.0, z
990 G4double u2 = 0.0, v2 = 0.0;
991 // if we are above the upper grid limit with
992 // => izotropic distribution: lambG1_max =7.
993 if (0.5*Qn1 > 7.0){
994 cosTheta1 = 1.-2.*G4UniformRand();
995 sinTheta1 = std::sqrt((1.-cosTheta1)*(1.+c
996 cosTheta2 = 1.-2.*G4UniformRand();
997 sinTheta2 = std::sqrt((1.-cosTheta2)*(1.+c
998 } else {
999 // sample 2 scattering cost1, sint1, cost
1000 G4double lekin = G4Log(efEnergy);
1001 G4double pt2 = efEnergy*(efEnergy+2.0
1002 G4double beta2 = pt2/(pt2+CLHEP::electr
1003 // backup GS angular dtr pointer (kineti
1004 // if the first was an msc sampling (the
1005 G4GoudsmitSaundersonTable::GSMSCAngularD
1006 G4int mcEkinIdx = -1;
1007 G4int mcDeltIdx = -1;
1008 G4double transfPar = 0.;
1009 G4bool isMsc = fGSTable->Sampling(0.5*la
1010 curren
1011 true);
1012 fGSTable->Sampling(0.5*lambdan, 0.5*Qn1,
1013 currentMaterialIndex,
1014 if (cosTheta1+cosTheta2==2.) { // no sca
1015 if (fIsEverythingWasDone)
1016 fTheZPathLenght = fTheTrueStepLeng
1017 fIsNoScatteringInMSC = true;
1018 return;
1019 }
1020 }
1021 // sample 2 azimuthal angles
1022 G4double phi1 = CLHEP::twopi*G4UniformRand(
1023 sinPhi1 = std::sin(phi1);
1024 cosPhi1 = std::cos(phi1);
1025 G4double phi2 = CLHEP::twopi*G4UniformRand(
1026 sinPhi2 = std::sin(phi2);
1027 cosPhi2 = std::cos(phi2);
1028
1029 // compute final direction realtive to z-di
1030 u2 = sinTheta2*cosPhi2;
1031 v2 = sinTheta2*sinPhi2;
1032 G4double u2p = cosTheta1*u2 + sinTheta1*cos
1033 uss = u2p*cosPhi1 - v2*sinPhi1;
1034 vss = u2p*sinPhi1 + v2*cosPhi1;
1035 wss = cosTheta1*cosTheta2 - sinTheta1*u2;
1036
1037 // set new direction (is scattering frame)
1038 fTheNewDirection.set(uss,vss,wss);
1039
1040 // set the fTheZPathLenght if we don't samp
1041 // we should do everything at the step-limi
1042 if(fIsNoDisplace && fIsEverythingWasDone)
1043 fTheZPathLenght = fTheTrueStepLenght;
1044
1045 // in optimized-mode if the current-safety
1046 if(fIsNoDisplace)
1047 return;
1048
1049 ///////////////////////////////////////////
1050 // Compute final position
1051 Qn1 *= fMCtoQ1;
1052 if (gIsUseAccurate) {
1053 // correction parameter
1054 G4double par =1.;
1055 if(Qn1<0.7) par = 1.;
1056 else if (Qn1<7.0) par = -0.031376*Qn1+1.
1057 else par = 0.79;
1058
1059 // Moments with energy loss correction
1060 // --first the uncorrected (for energy l
1061 // gamma = G_2/G_1 based on G2 computed
1062 G4double loga = G4Log(1.0+1.0/fScrA);
1063 G4double gamma = 6.0*fScrA*(1.0 + fScrA
1064 gamma *= fMCtoG2PerG1;
1065 // sample eta from p(eta)=2*eta i.e. P(e
1066 G4double eta = std::sqrt(G4UniformRan
1067 G4double eta1 = 0.5*(1 - eta); // use
1068 // 0.5 +sqrt(6)/6 = 0.9082483;
1069 // 1/(4*sqrt(6)) = 0.1020621;
1070 // (4-sqrt(6)/(24*sqrt(6))) = 0.02637471
1071 // delta = 0.9082483-(0.1020621-0.026374
1072 G4double delta = 0.9082483-(0.1020621-0
1073
1074 // compute alpha1 and alpha2 for energy
1075 G4double temp1 = 2.0 + tau;
1076 G4double temp = (2.0+tau*temp1)/((tau+1
1077 //Take logarithmic dependence
1078 temp = temp - (tau+1.0)/((tau+2.0)*(loga
1079 temp = temp * epsm;
1080 temp1 = 1.0 - temp;
1081 delta = delta + 0.40824829*(eps0*(tau+1.
1082 (loga*(1.0+fScrA)-1.0)*(loga*(1.
1083 G4double b = eta*delta;
1084 G4double c = eta*(1.0-delta);
1085
1086 //Calculate transport direction cosines:
1087 // ut,vt,wt is the final position divide
1088 G4double w1v2 = cosTheta1*v2;
1089 G4double ut = b*sinTheta1*cosPhi1 + c*
1090 G4double vt = b*sinTheta1*sinPhi1 + c*
1091 G4double wt = eta1*(1+temp) + b*
1092
1093 // long step correction
1094 ut *=par;
1095 vt *=par;
1096 wt *=par;
1097
1098 // final position relative to the pre-st
1099 // ut = x_f/s so needs to multiply by s
1100 x_coord = ut*fTheTrueStepLenght;
1101 y_coord = vt*fTheTrueStepLenght;
1102 z_coord = wt*fTheTrueStepLenght;
1103
1104 if(fIsEverythingWasDone){
1105 // We sample in the step limit so set
1106 // and lateral displacement (x_coord,y
1107 //Calculate transport distance
1108 G4double transportDistance = std::sqr
1109 // protection
1110 if(transportDistance>fTheTrueStepLengh
1111 transportDistance = fTheTrueStepLen
1112 fTheZPathLenght = transportDistance;
1113 }
1114 // else:: we sample in the DoIt so
1115 // the fTheZPathLenght was already
1116 fTheDisplacementVector.set(x_coord,y_coo
1117 } else {
1118 // compute zz = <z>/tPathLength
1119 // s -> true-path-length
1120 // z -> geom-path-length:: when PRESTA i
1121 // r -> lateral displacement = s/2 sin(t
1122 G4double zz = 0.0;
1123 if(fIsEverythingWasDone){
1124 // We sample in the step limit so set
1125 // and lateral displacement (x_coord,
1126 if(Qn1<0.1) { // use 3-order Taylor a
1127 zz = 1.0 - Qn1*(0.5 - Qn1*(0.166666
1128 } else {
1129 zz = (1.-G4Exp(-Qn1))/Qn1;
1130 }
1131 } else {
1132 // we sample in the DoIt so
1133 // the fTheZPathLenght was already se
1134 zz = fTheZPathLenght/fTheTrueStepLeng
1135 }
1136
1137 G4double rr = (1.-zz*zz)/(1.-wss*wss); /
1138 if(rr >= 0.25) rr = 0.25; //
1139 G4double rperp = fTheTrueStepLenght*std:
1140 x_coord = rperp*uss;
1141 y_coord = rperp*vss;
1142 z_coord = zz*fTheTrueStepLenght;
1143
1144 if(fIsEverythingWasDone){
1145 G4double transportDistance = std::sqrt
1146 fTheZPathLenght = transportDistance;
1147 }
1148
1149 fTheDisplacementVector.set(x_coord,y_coo
1150 }
1151 }
1152