| 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 // $Id: G4GoudsmitSaundersonMscModel.cc 79188 2014-02-20 09:22:48Z gcosmo $
26 // 27 //
27 // ------------------------------------------- << 28 // -------------------------------------------------------------------
28 // 29 //
29 // GEANT4 Class implementation file << 30 // GEANT4 Class file
30 // 31 //
31 // File name: G4GoudsmitSaundersonMscModel 32 // File name: G4GoudsmitSaundersonMscModel
32 // 33 //
33 // Author: Mihaly Novak / (Omrane Kadri << 34 // Author: Omrane Kadri
34 // 35 //
35 // Creation date: 20.02.2009 36 // Creation date: 20.02.2009
36 // 37 //
37 // Modifications: 38 // Modifications:
38 // 04.03.2009 V.Ivanchenko cleanup and format 39 // 04.03.2009 V.Ivanchenko cleanup and format according to Geant4 EM style
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 // 40 //
109 // Class description: << 41 // 15.04.2009 O.Kadri: cleanup: discard no scattering and single scattering theta
110 // Kawrakow-Bielajew Goudsmit-Saunderson MSC << 42 // sampling from SampleCosineTheta() which means the splitting
111 // for elastic scattering of e-/e+. Option, << 43 // step into two sub-steps occur only for msc regime
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 // 44 //
>> 45 // 12.06.2009 O.Kadri: linear log-log extrapolation of lambda0 & lambda1 between 1 GeV - 100 TeV
>> 46 // adding a theta min limit due to screening effect of the atomic nucleus
>> 47 // 26.08.2009 O.Kadri: Cubic Spline interpolation was replaced with polynomial method
>> 48 // within CalculateIntegrals method
>> 49 // 05.10.2009 O.Kadri: tuning small angle theta distributions
>> 50 // assuming the case of lambdan<1 as single scattering regime
>> 51 // tuning theta sampling for theta below the screening angle
>> 52 // 08.02.2010 O.Kadri: bugfix in compound xsection calculation and small angle computation
>> 53 // adding a rejection condition to hard collision angular sampling
>> 54 // ComputeTruePathLengthLimit was taken from G4WentzelVIModel
>> 55 // 26.03.2010 O.Kadri: direct xsection calculation not inverse of the inverse
>> 56 // angular sampling without large angle rejection method
>> 57 // longitudinal displacement is computed exactly from <z>
>> 58 // 12.05.2010 O.Kadri: exchange between target and projectile has as a condition the particle type (e-/e-)
>> 59 // some cleanup to minimize time consuming (adding lamdan12 & Qn12, changing the error to 1.0e-12 for scrA)
121 // 60 //
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 // 61 //
130 // ------------------------------------------- << 62 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
131 << 63 //REFERENCES:
132 << 64 //Ref.1:E. Benedito et al.,"Mixed simulation ... cross-sections", NIMB 174 (2001) pp 91-110;
>> 65 //Ref.2:I. Kawrakow et al.,"On the condensed ... transport",NIMB 142 (1998) pp 253-280;
>> 66 //Ref.3:I. Kawrakow et al.,"On the representation ... calculations",NIMB 134 (1998) pp 325-336;
>> 67 //Ref.4:Bielajew et al.,".....", NIMB 173 (2001) 332-343;
>> 68 //Ref.5:F. Salvat et al.,"ELSEPA--Dirac partial ...molecules", Comp.Phys.Comm.165 (2005) pp 157-190;
>> 69 //Ref.6:G4UrbanMscModel G4 9.2;
>> 70 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
133 #include "G4GoudsmitSaundersonMscModel.hh" 71 #include "G4GoudsmitSaundersonMscModel.hh"
134 <<
135 #include "G4GoudsmitSaundersonTable.hh" 72 #include "G4GoudsmitSaundersonTable.hh"
136 #include "G4GSPWACorrections.hh" <<
137 73
138 #include "G4PhysicalConstants.hh" 74 #include "G4PhysicalConstants.hh"
139 #include "G4SystemOfUnits.hh" 75 #include "G4SystemOfUnits.hh"
140 76
141 #include "G4ParticleChangeForMSC.hh" 77 #include "G4ParticleChangeForMSC.hh"
>> 78 #include "G4MaterialCutsCouple.hh"
142 #include "G4DynamicParticle.hh" 79 #include "G4DynamicParticle.hh"
143 #include "G4Electron.hh" 80 #include "G4Electron.hh"
144 #include "G4Positron.hh" 81 #include "G4Positron.hh"
145 82
146 #include "G4LossTableManager.hh" 83 #include "G4LossTableManager.hh"
147 #include "G4EmParameters.hh" <<
148 #include "G4Track.hh" 84 #include "G4Track.hh"
149 #include "G4PhysicsTable.hh" 85 #include "G4PhysicsTable.hh"
150 #include "Randomize.hh" 86 #include "Randomize.hh"
151 #include "G4Log.hh" 87 #include "G4Log.hh"
152 #include "G4Exp.hh" 88 #include "G4Exp.hh"
153 #include "G4Pow.hh" <<
154 #include <fstream> 89 #include <fstream>
155 90
>> 91 using namespace std;
156 92
157 // set accurate energy loss and dispalcement s << 93 G4double G4GoudsmitSaundersonMscModel::ener[] = {-1.};
158 G4bool G4GoudsmitSaundersonMscModel::gIsUseAcc << 94 G4double G4GoudsmitSaundersonMscModel::TCSE[103][106] ;
159 // set the usual optimization to be always act << 95 G4double G4GoudsmitSaundersonMscModel::FTCSE[103][106] ;
160 G4bool G4GoudsmitSaundersonMscModel::gIsOptimi << 96 G4double G4GoudsmitSaundersonMscModel::TCSP[103][106] ;
161 << 97 G4double G4GoudsmitSaundersonMscModel::FTCSP[103][106] ;
162 98
>> 99 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
163 G4GoudsmitSaundersonMscModel::G4GoudsmitSaunde 100 G4GoudsmitSaundersonMscModel::G4GoudsmitSaundersonMscModel(const G4String& nam)
164 : G4VMscModel(nam) { << 101 : G4VMscModel(nam),lowKEnergy(0.1*keV),highKEnergy(100.*TeV)
165 charge = 0; << 102 {
166 currentMaterialIndex = -1; << 103 currentKinEnergy=currentRange=skindepth=par1=par2=par3
167 // << 104 =zPathLength=truePathLength
168 fr = 0.1; << 105 =tausmall=taulim=tlimit=charge=lambdalimit=tPathLength=lambda0=lambda1
169 rangeinit = 1.e+21; << 106 =lambda11=mass=0.0;
170 geombig = 1.e+50*mm; << 107 currentMaterialIndex = -1;
171 geomlimit = geombig; << 108
172 tgeom = geombig; << 109 fr=0.02,rangeinit=0.,masslimite=0.6*MeV,
173 tlimit = 1.e+10*mm; << 110 particle=0;tausmall=1.e-16;taulim=1.e-6;tlimit=1.e10*mm;
174 presafety = 0.*mm; << 111 tlimitmin=10.e-6*mm;geombig=1.e50*mm;geommin=1.e-3*mm,tgeom=geombig;
175 // << 112 tlimitminfix=1.e-6*mm;stepmin=tlimitminfix;lambdalimit=1.*mm;smallstep=1.e10;
176 particle = nullptr; << 113 theManager=G4LossTableManager::Instance();
177 theManager = G4LossTableManager: << 114 inside=false;insideskin=false;
178 firstStep = true; << 115 samplez=false;
179 currentKinEnergy = 0.0; << 116 firstStep = true;
180 currentRange = 0.0; << 117
181 // << 118 GSTable = new G4GoudsmitSaundersonTable();
182 tlimitminfix2 = 1.*nm; << 119
183 tausmall = 1.e-16; << 120 if(ener[0] < 0.0){
184 mass = electron_mass_c2; << 121 G4cout << "### G4GoudsmitSaundersonMscModel loading ELSEPA data" << G4endl;
185 taulim = 1.e-6; << 122 LoadELSEPAXSections();
186 // << 123 }
187 currentCouple = nullptr; << 124 }
188 fParticleChange = nullptr; << 125 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
189 // << 126 G4GoudsmitSaundersonMscModel::~G4GoudsmitSaundersonMscModel()
190 fZeff = 1.; << 127 {
191 // << 128 delete GSTable;
192 par1 = 0.; << 129 }
193 par2 = 0.; << 130 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
194 par3 = 0.; << 131 void G4GoudsmitSaundersonMscModel::Initialise(const G4ParticleDefinition* p,
195 // << 132 const G4DataVector&)
196 // Moliere screeing parameter will be used a << 133 {
197 // appalied to the integrated quantities (sc << 134 skindepth=skin*stepmin;
198 // and second moments) derived from the corr << 135 SetParticle(p);
199 // this PWA correction is ignored if Mott-co << 136 fParticleChange = GetParticleChangeForMSC(p);
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 } 137 }
235 138
>> 139 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
236 140
237 G4GoudsmitSaundersonMscModel::~G4GoudsmitSaund << 141 G4double
238 if (IsMaster()) { << 142 G4GoudsmitSaundersonMscModel::ComputeCrossSectionPerAtom(const G4ParticleDefinition* p,
239 if (fGSTable) { << 143 G4double kineticEnergy,G4double Z, G4double, G4double, G4double)
240 delete fGSTable; << 144 {
241 fGSTable = nullptr; << 145 G4double kinEnergy = kineticEnergy;
242 } << 146 if(kinEnergy<lowKEnergy) kinEnergy=lowKEnergy;
243 if (fPWACorrection) { << 147 if(kinEnergy>highKEnergy)kinEnergy=highKEnergy;
244 delete fPWACorrection; << 148
245 fPWACorrection = nullptr; << 149 G4double cs(0.0), cs0(0.0);
246 } << 150 CalculateIntegrals(p,Z,kinEnergy,cs0,cs);
>> 151
>> 152 return cs;
>> 153 }
>> 154 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 155
>> 156 G4ThreeVector&
>> 157 G4GoudsmitSaundersonMscModel::SampleScattering(const G4ThreeVector& oldDirection, G4double)
>> 158 {
>> 159 fDisplacement.set(0.0,0.0,0.0);
>> 160 G4double kineticEnergy = currentKinEnergy;
>> 161 //dynParticle->GetKineticEnergy();
>> 162 if((kineticEnergy <= 0.0) || (tPathLength <= tlimitminfix)||
>> 163 (tPathLength/tausmall < lambda1)) { return fDisplacement; }
>> 164
>> 165 ///////////////////////////////////////////
>> 166 // Effective energy
>> 167 G4double eloss = 0.0;
>> 168 if (tPathLength > currentRange*dtrl) {
>> 169 eloss = kineticEnergy -
>> 170 GetEnergy(particle,currentRange-tPathLength,currentCouple);
>> 171 } else {
>> 172 eloss = tPathLength*GetDEDX(particle,kineticEnergy,currentCouple);
247 } 173 }
248 } << 174 /*
>> 175 G4double ttau = kineticEnergy/electron_mass_c2;
>> 176 G4double ttau2 = ttau*ttau;
>> 177 G4double epsilonpp = eloss/kineticEnergy;
>> 178 G4double cst1 = epsilonpp*epsilonpp*(6+10*ttau+5*ttau2)/(24*ttau2+48*ttau+72);
>> 179 kineticEnergy *= (1 - cst1);
>> 180 */
>> 181 kineticEnergy -= 0.5*eloss;
>> 182
>> 183 ///////////////////////////////////////////
>> 184 // additivity rule for mixture and compound xsection's
>> 185 const G4Material* mat = currentCouple->GetMaterial();
>> 186 const G4ElementVector* theElementVector = mat->GetElementVector();
>> 187 const G4double* theAtomNumDensityVector = mat->GetVecNbOfAtomsPerVolume();
>> 188 G4int nelm = mat->GetNumberOfElements();
>> 189 G4double s0(0.0), s1(0.0);
>> 190 lambda0 = 0.0;
>> 191 for(G4int i=0;i<nelm;i++)
>> 192 {
>> 193 CalculateIntegrals(particle,(*theElementVector)[i]->GetZ(),kineticEnergy,s0,s1);
>> 194 lambda0 += (theAtomNumDensityVector[i]*s0);
>> 195 }
>> 196 if(lambda0>0.0) { lambda0 =1./lambda0; }
>> 197
>> 198 // Newton-Raphson root's finding method of scrA from:
>> 199 // Sig1(PWA)/Sig0(PWA)=g1=2*scrA*((1+scrA)*log(1+1/scrA)-1)
>> 200 G4double g1=0.0;
>> 201 if(lambda1>0.0) { g1 = lambda0/lambda1; }
>> 202
>> 203 G4double logx0,x1,delta;
>> 204 G4double x0=g1*0.5;
>> 205 // V.Ivanchenko added limit of the loop
>> 206 for(G4int i=0;i<1000;++i)
>> 207 {
>> 208 logx0 = G4Log(1.+1./x0);
>> 209 x1 = x0-(x0*((1.+x0)*logx0-1.0)-g1*0.5)/( (1.+2.*x0)*logx0-2.0);
>> 210
>> 211 // V.Ivanchenko cut step size of iterative procedure
>> 212 if(x1 < 0.0) { x1 = 0.5*x0; }
>> 213 else if(x1 > 2*x0) { x1 = 2*x0; }
>> 214 else if(x1 < 0.5*x0) { x1 = 0.5*x0; }
>> 215 delta = std::fabs( x1 - x0 );
>> 216 x0 = x1;
>> 217 if(delta < 1.0e-3*x1) { break;}
>> 218 }
>> 219 G4double scrA = x1;
249 220
>> 221 G4double lambdan=0.;
>> 222
>> 223 if(lambda0>0.0) { lambdan=tPathLength/lambda0; }
>> 224 if(lambdan<=1.0e-12) { return fDisplacement; }
>> 225
>> 226 //G4cout << "E(eV)= " << kineticEnergy/eV << " L0= " << lambda0
>> 227 // << " L1= " << lambda1 << G4endl;
>> 228
>> 229 G4double Qn1 = lambdan *g1;//2.* lambdan *scrA*((1.+scrA)*log(1.+1./scrA)-1.);
>> 230 G4double Qn12 = 0.5*Qn1;
>> 231
>> 232 G4double cosTheta1,sinTheta1,cosTheta2,sinTheta2;
>> 233 G4double cosPhi1=1.0,sinPhi1=0.0,cosPhi2=1.0,sinPhi2=0.0;
>> 234 G4double us=0.0,vs=0.0,ws=1.0,wss=0.,x_coord=0.0,y_coord=0.0,z_coord=1.0;
>> 235
>> 236 G4double epsilon1=G4UniformRand();
>> 237 G4double expn = G4Exp(-lambdan);
>> 238
>> 239 if(epsilon1<expn)// no scattering
>> 240 { return fDisplacement; }
>> 241 else if((epsilon1<((1.+lambdan)*expn))||(lambdan<1.))//single or plural scattering (Rutherford DCS's)
>> 242 {
>> 243 G4double xi=G4UniformRand();
>> 244 xi= 2.*scrA*xi/(1.-xi + scrA);
>> 245 if(xi<0.)xi=0.;
>> 246 else if(xi>2.)xi=2.;
>> 247 ws=(1. - xi);
>> 248 wss=std::sqrt(xi*(2.-xi));
>> 249 G4double phi0=CLHEP::twopi*G4UniformRand();
>> 250 us=wss*cos(phi0);
>> 251 vs=wss*sin(phi0);
>> 252 }
>> 253 else // multiple scattering
>> 254 {
>> 255 // Ref.2 subsection 4.4 "The best solution found"
>> 256 // Sample first substep scattering angle
>> 257 SampleCosineTheta(0.5*lambdan,scrA,cosTheta1,sinTheta1);
>> 258 G4double phi1 = CLHEP::twopi*G4UniformRand();
>> 259 cosPhi1 = cos(phi1);
>> 260 sinPhi1 = sin(phi1);
>> 261
>> 262 // Sample second substep scattering angle
>> 263 SampleCosineTheta(0.5*lambdan,scrA,cosTheta2,sinTheta2);
>> 264 G4double phi2 = CLHEP::twopi*G4UniformRand();
>> 265 cosPhi2 = cos(phi2);
>> 266 sinPhi2 = sin(phi2);
>> 267
>> 268 // Overall scattering direction
>> 269 us = sinTheta2*(cosTheta1*cosPhi1*cosPhi2 - sinPhi1*sinPhi2) + cosTheta2*sinTheta1*cosPhi1;
>> 270 vs = sinTheta2*(cosTheta1*sinPhi1*cosPhi2 + cosPhi1*sinPhi2) + cosTheta2*sinTheta1*sinPhi1;
>> 271 ws = cosTheta1*cosTheta2 - sinTheta1*sinTheta2*cosPhi2;
>> 272 G4double sqrtA=sqrt(scrA);
>> 273 if(acos(ws)<sqrtA)//small angle approximation for theta less than screening angle
>> 274 {
>> 275 G4int i=0;
>> 276 do{i++;
>> 277 ws=1.+Qn12*G4Log(G4UniformRand());
>> 278 }while((fabs(ws)>1.)&&(i<20));//i<20 to avoid time consuming during the run
>> 279 if(i>=19)ws=cos(sqrtA);
>> 280 wss=std::sqrt((1.-ws*ws));
>> 281 us=wss*std::cos(phi1);
>> 282 vs=wss*std::sin(phi1);
>> 283 }
>> 284 }
>> 285
>> 286 //G4ThreeVector oldDirection = dynParticle->GetMomentumDirection();
>> 287 G4ThreeVector newDirection(us,vs,ws);
>> 288 newDirection.rotateUz(oldDirection);
>> 289 fParticleChange->ProposeMomentumDirection(newDirection);
>> 290
>> 291 // corresponding to error less than 1% in the exact formula of <z>
>> 292 if(Qn1<0.02) { z_coord = 1.0 - Qn1*(0.5 - Qn1/6.); }
>> 293 else { z_coord = (1.-G4Exp(-Qn1))/Qn1; }
>> 294 G4double rr = zPathLength*std::sqrt((1.- z_coord*z_coord)/(1.-ws*ws));
>> 295 x_coord = rr*us;
>> 296 y_coord = rr*vs;
>> 297
>> 298 // displacement is computed relatively to the end point
>> 299 z_coord -= 1.0;
>> 300 z_coord *= zPathLength;
>> 301 /*
>> 302 G4cout << "G4GS::SampleSecondaries: e(MeV)= " << kineticEnergy
>> 303 << " sinTheta= " << sqrt(1.0 - ws*ws)
>> 304 << " trueStep(mm)= " << tPathLength
>> 305 << " geomStep(mm)= " << zPathLength
>> 306 << G4endl;
>> 307 */
>> 308
>> 309 fDisplacement.set(x_coord,y_coord,z_coord);
>> 310 fDisplacement.rotateUz(oldDirection);
>> 311
>> 312 return fDisplacement;
>> 313 }
>> 314
>> 315 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 316
>> 317 void
>> 318 G4GoudsmitSaundersonMscModel::SampleCosineTheta(G4double lambdan, G4double scrA,
>> 319 G4double &cost, G4double &sint)
>> 320 {
>> 321 G4double r1,tet,xi=0.;
>> 322 G4double Qn1 = 2.* lambdan;
>> 323 if(scrA < 10.) { Qn1 *= scrA*((1.+scrA)*G4Log(1.+1./scrA)-1.); }
>> 324 else { Qn1*= (1.0 - 0.5/scrA - 0.5/(scrA*scrA)) ; }
>> 325 if (Qn1<0.001)
>> 326 {
>> 327 do{
>> 328 r1=G4UniformRand();
>> 329 xi=-0.5*Qn1*G4Log(G4UniformRand());
>> 330 tet=acos(1.-xi);
>> 331 }while(tet*r1*r1>sin(tet));
>> 332 }
>> 333 else if(Qn1>0.5) { xi=2.*G4UniformRand(); }//isotropic distribution
>> 334 else{ xi=2.*(GSTable->SampleTheta(lambdan,scrA,G4UniformRand()));}
>> 335
>> 336
>> 337 if(xi<0.)xi=0.;
>> 338 else if(xi>2.)xi=2.;
>> 339
>> 340 cost=(1. - xi);
>> 341 sint=sqrt(xi*(2.-xi));
>> 342
>> 343 }
>> 344
>> 345 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 346 // Polynomial log-log interpolation of Lambda0 and Lambda1 between 100 eV - 1 GeV
>> 347 // linear log-log extrapolation between 1 GeV - 100 TeV
>> 348
>> 349 void
>> 350 G4GoudsmitSaundersonMscModel::CalculateIntegrals(const G4ParticleDefinition* p,G4double Z,
>> 351 G4double kinEnergy,G4double &Sig0,
>> 352 G4double &Sig1)
>> 353 {
>> 354 G4double x1,x2,y1,y2,acoeff,bcoeff;
>> 355 G4double kineticE = kinEnergy;
>> 356 if(kineticE<lowKEnergy)kineticE=lowKEnergy;
>> 357 if(kineticE>highKEnergy)kineticE=highKEnergy;
>> 358 kineticE /= eV;
>> 359 G4double logE=G4Log(kineticE);
>> 360
>> 361 G4int iZ = G4int(Z);
>> 362 if(iZ > 103) iZ = 103;
>> 363
>> 364 G4int enerInd=0;
>> 365 for(G4int i=0;i<105;i++)
>> 366 {
>> 367 if((logE>=ener[i])&&(logE<ener[i+1])){enerInd=i;break;}
>> 368 }
>> 369
>> 370 if(p==G4Electron::Electron())
>> 371 {
>> 372 if(kineticE<=1.0e+9)//Interpolation of the form y=ax²+b
>> 373 {
>> 374 x1=ener[enerInd];
>> 375 x2=ener[enerInd+1];
>> 376 y1=TCSE[iZ-1][enerInd];
>> 377 y2=TCSE[iZ-1][enerInd+1];
>> 378 acoeff=(y2-y1)/(x2*x2-x1*x1);
>> 379 bcoeff=y2-acoeff*x2*x2;
>> 380 Sig0=acoeff*logE*logE+bcoeff;
>> 381 Sig0 =G4Exp(Sig0);
>> 382 y1=FTCSE[iZ-1][enerInd];
>> 383 y2=FTCSE[iZ-1][enerInd+1];
>> 384 acoeff=(y2-y1)/(x2*x2-x1*x1);
>> 385 bcoeff=y2-acoeff*x2*x2;
>> 386 Sig1=acoeff*logE*logE+bcoeff;
>> 387 Sig1=G4Exp(Sig1);
>> 388 }
>> 389 else //Interpolation of the form y=ax+b
>> 390 {
>> 391 x1=ener[104];
>> 392 x2=ener[105];
>> 393 y1=TCSE[iZ-1][104];
>> 394 y2=TCSE[iZ-1][105];
>> 395 Sig0=(y2-y1)*(logE-x1)/(x2-x1)+y1;
>> 396 Sig0=G4Exp(Sig0);
>> 397 y1=FTCSE[iZ-1][104];
>> 398 y2=FTCSE[iZ-1][105];
>> 399 Sig1=(y2-y1)*(logE-x1)/(x2-x1)+y1;
>> 400 Sig1=G4Exp(Sig1);
>> 401 }
>> 402 }
>> 403 if(p==G4Positron::Positron())
>> 404 {
>> 405 if(kinEnergy<=1.0e+9)
>> 406 {
>> 407 x1=ener[enerInd];
>> 408 x2=ener[enerInd+1];
>> 409 y1=TCSP[iZ-1][enerInd];
>> 410 y2=TCSP[iZ-1][enerInd+1];
>> 411 acoeff=(y2-y1)/(x2*x2-x1*x1);
>> 412 bcoeff=y2-acoeff*x2*x2;
>> 413 Sig0=acoeff*logE*logE+bcoeff;
>> 414 Sig0 =G4Exp(Sig0);
>> 415 y1=FTCSP[iZ-1][enerInd];
>> 416 y2=FTCSP[iZ-1][enerInd+1];
>> 417 acoeff=(y2-y1)/(x2*x2-x1*x1);
>> 418 bcoeff=y2-acoeff*x2*x2;
>> 419 Sig1=acoeff*logE*logE+bcoeff;
>> 420 Sig1=G4Exp(Sig1);
>> 421 }
>> 422 else
>> 423 {
>> 424 x1=ener[104];
>> 425 x2=ener[105];
>> 426 y1=TCSP[iZ-1][104];
>> 427 y2=TCSP[iZ-1][105];
>> 428 Sig0=(y2-y1)*(logE-x1)/(x2-x1)+y1;
>> 429 Sig0 =G4Exp(Sig0);
>> 430 y1=FTCSP[iZ-1][104];
>> 431 y2=FTCSP[iZ-1][105];
>> 432 Sig1=(y2-y1)*(logE-x1)/(x2-x1)+y1;
>> 433 Sig1=G4Exp(Sig1);
>> 434 }
>> 435 }
>> 436
>> 437 Sig0 *= barn;
>> 438 Sig1 *= barn;
250 439
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 } 440 }
298 441
>> 442 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
299 443
300 void G4GoudsmitSaundersonMscModel::InitialiseL << 444 void G4GoudsmitSaundersonMscModel::StartTracking(G4Track* track)
301 fGSTable = static_cast<G4Goud << 445 {
302 fIsUseMottCorrection = static_cast<G4Goud << 446 SetParticle(track->GetDynamicParticle()->GetDefinition());
303 fIsUsePWACorrection = static_cast<G4Goud << 447 firstStep = true;
304 fPWACorrection = static_cast<G4Goud << 448 inside = false;
>> 449 insideskin = false;
>> 450 tlimit = geombig;
305 } 451 }
306 452
>> 453 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 454 //t->g->t step transformations taken from Ref.6
>> 455
>> 456 G4double
>> 457 G4GoudsmitSaundersonMscModel::ComputeTruePathLengthLimit(const G4Track& track,
>> 458 G4double& currentMinimalStep)
>> 459 {
>> 460 tPathLength = currentMinimalStep;
>> 461 const G4DynamicParticle* dp = track.GetDynamicParticle();
>> 462 G4StepPoint* sp = track.GetStep()->GetPreStepPoint();
>> 463 G4StepStatus stepStatus = sp->GetStepStatus();
>> 464 currentCouple = track.GetMaterialCutsCouple();
>> 465 SetCurrentCouple(currentCouple);
>> 466 currentMaterialIndex = currentCouple->GetIndex();
>> 467 currentKinEnergy = dp->GetKineticEnergy();
>> 468 currentRange = GetRange(particle,currentKinEnergy,currentCouple);
>> 469
>> 470 lambda1 = GetTransportMeanFreePath(particle,currentKinEnergy);
>> 471
>> 472 // stop here if small range particle
>> 473 if(inside || tPathLength < tlimitminfix) {
>> 474 return ConvertTrueToGeom(tPathLength, currentMinimalStep);
>> 475 }
>> 476 if(tPathLength > currentRange) tPathLength = currentRange;
>> 477
>> 478 G4double presafety = sp->GetSafety();
>> 479
>> 480 //G4cout << "G4GS::StepLimit tPathLength= "
>> 481 // <<tPathLength<<" safety= " << presafety
>> 482 // << " range= " <<currentRange<< " lambda= "<<lambda1
>> 483 // << " Alg: " << steppingAlgorithm <<G4endl;
>> 484
>> 485 // far from geometry boundary
>> 486 if(currentRange < presafety)
>> 487 {
>> 488 inside = true;
>> 489 return ConvertTrueToGeom(tPathLength, currentMinimalStep);
>> 490 }
>> 491
>> 492 // standard version
>> 493 //
>> 494 if (steppingAlgorithm == fUseDistanceToBoundary)
>> 495 {
>> 496 //compute geomlimit and presafety
>> 497 G4double geomlimit = ComputeGeomLimit(track, presafety, tPathLength);
>> 498
>> 499 // is far from boundary
>> 500 if(currentRange <= presafety)
>> 501 {
>> 502 inside = true;
>> 503 return ConvertTrueToGeom(tPathLength, currentMinimalStep);
>> 504 }
>> 505
>> 506 smallstep += 1.;
>> 507 insideskin = false;
>> 508
>> 509 if(firstStep || stepStatus == fGeomBoundary)
>> 510 {
>> 511 rangeinit = currentRange;
>> 512 if(firstStep) smallstep = 1.e10;
>> 513 else smallstep = 1.;
>> 514
>> 515 //define stepmin here (it depends on lambda!)
>> 516 //rough estimation of lambda_elastic/lambda_transport
>> 517 G4double rat = currentKinEnergy/MeV ;
>> 518 rat = 1.e-3/(rat*(10.+rat)) ;
>> 519 //stepmin ~ lambda_elastic
>> 520 stepmin = rat*lambda1;
>> 521 skindepth = skin*stepmin;
>> 522 //define tlimitmin
>> 523 tlimitmin = 10.*stepmin;
>> 524 if(tlimitmin < tlimitminfix) tlimitmin = tlimitminfix;
>> 525
>> 526 //G4cout << "rangeinit= " << rangeinit << " stepmin= " << stepmin
>> 527 // << " tlimitmin= " << tlimitmin << " geomlimit= " << geomlimit <<G4endl;
>> 528 // constraint from the geometry
>> 529 if((geomlimit < geombig) && (geomlimit > geommin))
>> 530 {
>> 531 if(stepStatus == fGeomBoundary)
>> 532 tgeom = geomlimit/facgeom;
>> 533 else
>> 534 tgeom = 2.*geomlimit/facgeom;
>> 535 }
>> 536 else
>> 537 tgeom = geombig;
307 538
308 // computes macroscopic first transport cross << 539 }
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 540
>> 541 //step limit
>> 542 tlimit = facrange*rangeinit;
>> 543 if(tlimit < facsafety*presafety)
>> 544 tlimit = facsafety*presafety;
>> 545
>> 546 //lower limit for tlimit
>> 547 if(tlimit < tlimitmin) tlimit = tlimitmin;
>> 548
>> 549 if(tlimit > tgeom) tlimit = tgeom;
>> 550
>> 551 //G4cout << "tgeom= " << tgeom << " geomlimit= " << geomlimit
>> 552 // << " tlimit= " << tlimit << " presafety= " << presafety << G4endl;
>> 553
>> 554 // shortcut
>> 555 if((tPathLength < tlimit) && (tPathLength < presafety) &&
>> 556 (smallstep >= skin) && (tPathLength < geomlimit-0.999*skindepth))
>> 557 return ConvertTrueToGeom(tPathLength, currentMinimalStep);
>> 558
>> 559 // step reduction near to boundary
>> 560 if(smallstep < skin)
>> 561 {
>> 562 tlimit = stepmin;
>> 563 insideskin = true;
>> 564 }
>> 565 else if(geomlimit < geombig)
>> 566 {
>> 567 if(geomlimit > skindepth)
>> 568 {
>> 569 if(tlimit > geomlimit-0.999*skindepth)
>> 570 tlimit = geomlimit-0.999*skindepth;
>> 571 }
>> 572 else
>> 573 {
>> 574 insideskin = true;
>> 575 if(tlimit > stepmin) tlimit = stepmin;
>> 576 }
>> 577 }
>> 578
>> 579 if(tlimit < stepmin) tlimit = stepmin;
>> 580
>> 581 if(tPathLength > tlimit) tPathLength = tlimit;
>> 582
>> 583 }
>> 584 // for 'normal' simulation with or without magnetic field
>> 585 // there no small step/single scattering at boundaries
>> 586 else if(steppingAlgorithm == fUseSafety)
>> 587 {
>> 588 // compute presafety again if presafety <= 0 and no boundary
>> 589 // i.e. when it is needed for optimization purposes
>> 590 if((stepStatus != fGeomBoundary) && (presafety < tlimitminfix))
>> 591 presafety = ComputeSafety(sp->GetPosition(),tPathLength);
>> 592
>> 593 // is far from boundary
>> 594 if(currentRange < presafety)
>> 595 {
>> 596 inside = true;
>> 597 return ConvertTrueToGeom(tPathLength, currentMinimalStep);
>> 598 }
362 599
363 // gives back the first transport mean free pa << 600 if(firstStep || stepStatus == fGeomBoundary)
364 G4double << 601 {
365 G4GoudsmitSaundersonMscModel::GetTransportMean << 602 rangeinit = currentRange;
366 << 603 fr = facrange;
367 // kinetic energy is assumed to be in Geant4 << 604 // 9.1 like stepping for e+/e- only (not for muons,hadrons)
368 G4double efEnergy = kineticEnergy; << 605 if(mass < masslimite)
369 // << 606 {
370 const G4Material* mat = currentCouple->GetM << 607 if(lambda1 > currentRange)
371 // << 608 rangeinit = lambda1;
372 fLambda0 = 0.0; // elastic mean free path << 609 if(lambda1 > lambdalimit)
373 fLambda1 = 0.0; // first transport mean free << 610 fr *= 0.75+0.25*lambda1/lambdalimit;
374 fScrA = 0.0; // screening parameter << 611 }
375 fG1 = 0.0; // first transport coef. << 612
376 << 613 //lower limit for tlimit
377 // use Moliere's screening (with Mott-corret << 614 G4double rat = currentKinEnergy/MeV ;
378 if (efEnergy<10.*CLHEP::eV) efEnergy = 10.* << 615 rat = 1.e-3/(rat*(10.+rat)) ;
379 // total mometum square << 616 tlimitmin = 10.*lambda1*rat;
380 G4double pt2 = efEnergy*(efEnergy+2.0*el << 617 if(tlimitmin < tlimitminfix) tlimitmin = tlimitminfix;
381 // beta square << 618 }
382 G4double beta2 = pt2/(pt2+electron_mass_c2 << 619 //step limit
383 // current material index << 620 tlimit = fr*rangeinit;
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 621
415 return fLambda1; << 622 if(tlimit < facsafety*presafety)
416 } << 623 tlimit = facsafety*presafety;
417 624
>> 625 //lower limit for tlimit
>> 626 if(tlimit < tlimitmin) tlimit = tlimitmin;
418 627
419 G4double << 628 if(tPathLength > tlimit) tPathLength = tlimit;
420 G4GoudsmitSaundersonMscModel::GetTransportMean << 629 }
421 << 630
422 // kinetic energy is assumed to be in Geant4 << 631 // version similar to 7.1 (needed for some experiments)
423 G4double efEnergy = kineticEnergy; << 632 else
424 // << 633 {
425 const G4Material* mat = currentCouple->GetM << 634 if (stepStatus == fGeomBoundary)
426 // << 635 {
427 G4double lambda0 = 0.0; // elastc mean free << 636 if (currentRange > lambda1) tlimit = facrange*currentRange;
428 G4double lambda1 = 0.0; // first transport m << 637 else tlimit = facrange*lambda1;
429 G4double scrA = 0.0; // screening paramet << 638
430 G4double g1 = 0.0; // first transport m << 639 if(tlimit < tlimitmin) tlimit = tlimitmin;
431 << 640 if(tPathLength > tlimit) tPathLength = tlimit;
432 // use Moliere's screening (with Mott-corret << 641 }
433 if (efEnergy<10.*CLHEP::eV) efEnergy = 10.* << 642 }
434 // total mometum square in Geant4 internal e << 643 //G4cout << "tPathLength= " << tPathLength
435 G4double pt2 = efEnergy*(efEnergy+2.0*el << 644 // << " currentMinimalStep= " << currentMinimalStep << G4endl;
436 G4double beta2 = pt2/(pt2+electron_mass_c2 << 645 return ConvertTrueToGeom(tPathLength, currentMinimalStep);
437 G4int matindx = (G4int)mat->GetIndex(); << 646 }
438 G4double bc = fGSTable->GetMoliereBc(ma << 647
439 // get the Mott-correcton factors if Mott-co << 648 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
440 G4double mctoScrA = 1.0; << 649 // taken from Ref.6
441 G4double mctoQ1 = 1.0; << 650 G4double G4GoudsmitSaundersonMscModel::ComputeGeomPathLength(G4double)
442 G4double mctoG2PerG1 = 1.0; << 651 {
443 G4double scpCor = 1.0; << 652 firstStep = false;
444 if (fIsUseMottCorrection) { << 653 par1 = -1. ;
445 fGSTable->GetMottCorrectionFactors(G4Log(e << 654 par2 = par3 = 0. ;
446 scpCor = fGSTable->ComputeScatteringPowerC << 655
447 } else if (fIsUsePWACorrection) { << 656 // do the true -> geom transformation
448 fPWACorrection->GetPWACorrectionFactors(G4 << 657 zPathLength = tPathLength;
449 // scpCor = fGSTable->ComputeScatteringPow << 658
450 } << 659 // z = t for very small tPathLength
451 scrA = fGSTable->GetMoliereXc2(matindx)/( << 660 if(tPathLength < tlimitminfix) { return zPathLength; }
452 // total elastic mean free path in Geant4 in << 661
453 lambda0 = beta2*(1.+scrA)*mctoScrA/bc/scpCor << 662 // this correction needed to run MSC with eIoni and eBrem inactivated
454 g1 = 2.0*scrA*((1.0+scrA)*G4Log(1.0/scr << 663 // and makes no harm for a normal run
455 lambda1 = lambda0/g1; << 664 if(tPathLength > currentRange)
>> 665 { tPathLength = currentRange; }
>> 666
>> 667 G4double tau = tPathLength/lambda1 ;
>> 668
>> 669 if ((tau <= tausmall) || insideskin) {
>> 670 zPathLength = tPathLength;
>> 671 if(zPathLength > lambda1) { zPathLength = lambda1; }
>> 672 return zPathLength;
>> 673 }
>> 674
>> 675 G4double zmean = tPathLength;
>> 676 if (tPathLength < currentRange*dtrl) {
>> 677 if(tau < taulim) zmean = tPathLength*(1.-0.5*tau) ;
>> 678 else zmean = lambda1*(1.-G4Exp(-tau));
>> 679 } else if(currentKinEnergy < mass || tPathLength == currentRange) {
>> 680 par1 = 1./currentRange ;
>> 681 par2 = 1./(par1*lambda1) ;
>> 682 par3 = 1.+par2 ;
>> 683 if(tPathLength < currentRange)
>> 684 zmean = (1.-G4Exp(par3*G4Log(1.-tPathLength/currentRange)))/(par1*par3) ;
>> 685 else
>> 686 zmean = 1./(par1*par3) ;
>> 687 } else {
>> 688 G4double T1 = GetEnergy(particle,currentRange-tPathLength,currentCouple);
456 689
457 return lambda1; << 690 lambda11 = GetTransportMeanFreePath(particle,T1);
458 } <<
459 691
>> 692 par1 = (lambda1-lambda11)/(lambda1*tPathLength) ;
>> 693 par2 = 1./(par1*lambda1) ;
>> 694 par3 = 1.+par2 ;
>> 695 zmean = (1.-G4Exp(par3*G4Log(lambda11/lambda1)))/(par1*par3) ;
>> 696 }
460 697
461 void G4GoudsmitSaundersonMscModel::StartTracki << 698 zPathLength = zmean ;
462 SetParticle(track->GetDynamicParticle()->Get << 699 // sample z
463 firstStep = true; << 700 if(samplez) {
464 tlimit = tgeom = rangeinit = geombig; <<
465 rangeinit = 1.e+21; <<
466 } <<
467 701
>> 702 const G4double ztmax = 0.99;
>> 703 G4double zt = zmean/tPathLength ;
468 704
469 G4double G4GoudsmitSaundersonMscModel::Compute << 705 if (tPathLength > stepmin && zt < ztmax) {
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 706
738 } << 707 G4double u,cz1;
739 //step limit << 708 if(zt >= 0.333333333) {
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 709
>> 710 G4double cz = 0.5*(3.*zt-1.)/(1.-zt) ;
>> 711 cz1 = 1.+cz ;
>> 712 G4double u0 = cz/cz1 ;
>> 713 G4double grej ;
>> 714 do {
>> 715 u = G4Exp(G4Log(G4UniformRand())/cz1) ;
>> 716 grej = exp(cz*G4Log(u/u0))*(1.-u)/(1.-u0) ;
>> 717 } while (grej < G4UniformRand()) ;
780 718
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 { 719 } else {
813 fTheZPathLenght = 1./(par1*par3); << 720 cz1 = 1./zt-1.;
>> 721 u = 1.-G4Exp(G4Log(G4UniformRand())/cz1) ;
814 } 722 }
815 } else { << 723 zPathLength = tPathLength*u ;
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 } 724 }
826 } 725 }
827 fTheZPathLenght = std::min(fTheZPathLenght, << 726 if(zPathLength > lambda1) zPathLength = lambda1;
828 // << 727 //G4cout << "zPathLength= " << zPathLength << " lambda1= " << lambda1 << G4endl;
829 return fTheZPathLenght; << 728
>> 729 return zPathLength;
830 } 730 }
831 731
>> 732 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 733 // taken from Ref.6
>> 734 G4double
>> 735 G4GoudsmitSaundersonMscModel::ComputeTrueStepLength(G4double geomStepLength)
>> 736 {
>> 737 // step defined other than transportation
>> 738 if(geomStepLength == zPathLength && tPathLength <= currentRange)
>> 739 return tPathLength;
832 740
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 741 // t = z for very small step
851 if (geomStepLength<tlimitminfix2) { << 742 zPathLength = geomStepLength;
852 fTheTrueStepLenght = geomStepLength; << 743 tPathLength = geomStepLength;
>> 744 if(geomStepLength < tlimitminfix) return tPathLength;
>> 745
853 // recalculation 746 // recalculation
854 } else { << 747 if((geomStepLength > lambda1*tausmall) && !insideskin)
855 G4double tlength = geomStepLength; << 748 {
856 if (geomStepLength>fLambda1*tausmall) { << 749 if(par1 < 0.)
857 if (par1< 0.) { << 750 tPathLength = -lambda1*G4Log(1.-geomStepLength/lambda1) ;
858 tlength = -fLambda1*G4Log(1.-geomStepL << 751 else
>> 752 {
>> 753 if(par1*par3*geomStepLength < 1.)
>> 754 tPathLength = (1.-G4Exp(G4Log(1.-par1*par3*geomStepLength)/par3))/par1 ;
>> 755 else
>> 756 tPathLength = currentRange;
>> 757 }
>> 758 }
>> 759 if(tPathLength < geomStepLength) tPathLength = geomStepLength;
>> 760 //G4cout << "tPathLength= " << tPathLength << " step= " << geomStepLength << G4endl;
>> 761
>> 762 return tPathLength;
>> 763 }
>> 764
>> 765 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
>> 766 //Total & first transport x sections for e-/e+ generated from ELSEPA code
>> 767
>> 768 void G4GoudsmitSaundersonMscModel::LoadELSEPAXSections()
>> 769 {
>> 770 G4String filename = "XSECTIONS.dat";
>> 771
>> 772 char* path = getenv("G4LEDATA");
>> 773 if (!path)
>> 774 {
>> 775 G4Exception("G4GoudsmitSaundersonMscModel::LoadELSEPAXSections()","em0006",
>> 776 FatalException,
>> 777 "Environment variable G4LEDATA not defined");
>> 778 return;
>> 779 }
>> 780
>> 781 G4String pathString(path);
>> 782 G4String dirFile = pathString + "/msc_GS/" + filename;
>> 783 std::ifstream infile(dirFile, std::ios::in);
>> 784 if( !infile.is_open()) {
>> 785 G4ExceptionDescription ed;
>> 786 ed << "Data file <" + dirFile + "> is not opened!";
>> 787 G4Exception("G4GoudsmitSaundersonMscModel::LoadELSEPAXSections()",
>> 788 "em0003",FatalException,ed);
>> 789 return;
>> 790 }
>> 791
>> 792 // Read parameters from tables and take logarithms
>> 793 G4double aRead;
>> 794 for(G4int i=0 ; i<106 ;i++){
>> 795 infile >> aRead;
>> 796 if(infile.fail()) {
>> 797 G4ExceptionDescription ed;
>> 798 ed << "Error reading <" + dirFile + "> loop #1 i= " << i;
>> 799 G4Exception("G4GoudsmitSaundersonMscModel::LoadELSEPAXSections()",
>> 800 "em0003",FatalException,ed);
>> 801 return;
>> 802 } else {
>> 803 if(aRead > 0.0) { aRead = G4Log(aRead); }
>> 804 else { aRead = 0.0; }
>> 805 }
>> 806 ener[i]=aRead;
>> 807 }
>> 808 for(G4int j=0;j<103;j++){
>> 809 for(G4int i=0;i<106;i++){
>> 810 infile >> aRead;
>> 811 if(infile.fail()) {
>> 812 G4ExceptionDescription ed;
>> 813 ed << "Error reading <" + dirFile + "> loop #2 j= " << j
>> 814 << "; i= " << i;
>> 815 G4Exception("G4GoudsmitSaundersonMscModel::LoadELSEPAXSections()",
>> 816 "em0003",FatalException,ed);
>> 817 return;
859 } else { 818 } else {
860 if (par1*par3*geomStepLength<1.) { << 819 if(aRead > 0.0) { aRead = G4Log(aRead); }
861 G4Pow *g4calc = G4Pow::GetInstance() << 820 else { aRead = 0.0; }
862 tlength = (1.-g4calc->powA( 1.-par1* <<
863 } else { <<
864 tlength = currentRange; <<
865 } <<
866 } <<
867 if (tlength<geomStepLength || tlength>fT <<
868 tlength = geomStepLength; <<
869 } 821 }
870 } << 822 TCSE[j][i]=aRead;
871 fTheTrueStepLenght = tlength; << 823 }
872 } 824 }
873 // << 825 for(G4int j=0;j<103;j++){
874 return fTheTrueStepLenght; << 826 for(G4int i=0;i<106;i++){
875 } << 827 infile >> aRead;
876 << 828 if(infile.fail()) {
877 G4ThreeVector& << 829 G4ExceptionDescription ed;
878 G4GoudsmitSaundersonMscModel::SampleScattering << 830 ed << "Error reading <" + dirFile + "> loop #3 j= " << j
879 if (steppingAlgorithm==fUseDistanceToBoundar << 831 << "; i= " << i;
880 // single scattering was and scattering ha << 832 G4Exception("G4GoudsmitSaundersonMscModel::LoadELSEPAXSections()",
881 fTheNewDirection.rotateUz(oldDirection); << 833 "em0003",FatalException,ed);
882 fParticleChange->ProposeMomentumDirection( << 834 return;
883 return fTheDisplacementVector; << 835 } else {
884 } else if (steppingAlgorithm==fUseSafetyPlus << 836 if(aRead > 0.0) { aRead = G4Log(aRead); }
885 if (fIsEndedUpOnBoundary) { // do nothing << 837 else { aRead = 0.0; }
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 } 838 }
894 // check if multiple scattering happened << 839 FTCSE[j][i]=aRead;
895 if (fIsMultipleSacettring && !fIsNoScatt << 840 }
896 fTheNewDirection.rotateUz(oldDirect << 841 }
897 fTheDisplacementVector.rotateUz(old << 842 for(G4int j=0;j<103;j++){
898 fParticleChange->ProposeMomentumDir << 843 for(G4int i=0;i<106;i++){
>> 844 infile >> aRead;
>> 845 if(infile.fail()) {
>> 846 G4ExceptionDescription ed;
>> 847 ed << "Error reading <" + dirFile + "> loop #4 j= " << j
>> 848 << "; i= " << i;
>> 849 G4Exception("G4GoudsmitSaundersonMscModel::LoadELSEPAXSections()",
>> 850 "em0003",FatalException,ed);
>> 851 return;
>> 852 } else {
>> 853 if(aRead > 0.0) { aRead = G4Log(aRead); }
>> 854 else { aRead = 0.0; }
899 } 855 }
900 // The only thing that could happen if w << 856 TCSP[j][i]=aRead;
901 // is that single scattering was tried << 857 }
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 } 858 }
918 // << 859 for(G4int j=0;j<103;j++){
919 return fTheDisplacementVector; << 860 for(G4int i=0;i<106;i++){
920 } << 861 infile >> aRead;
921 << 862 if(infile.fail()) {
922 << 863 G4ExceptionDescription ed;
923 void G4GoudsmitSaundersonMscModel::SampleMSC() << 864 ed << "Error reading <" + dirFile + "> loop #5 j= " << j
924 fIsNoScatteringInMSC = false; << 865 << "; i= " << i;
925 // kinetic energy is assumed to be in Geant4 << 866 G4Exception("G4GoudsmitSaundersonMscModel::LoadELSEPAXSections()",
926 G4double kineticEnergy = currentKinEnergy; << 867 "em0003",FatalException,ed);
927 // << 868 return;
928 // Energy loss correction: 2 version << 869 } else {
929 G4double eloss = 0.0; << 870 if(aRead > 0.0) { aRead = G4Log(aRead); }
930 // if (fTheTrueStepLenght > currentRange*dtrl << 871 else { aRead = 0.0; }
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 } 872 }
979 fIsNoScatteringInMSC = true; << 873 FTCSP[j][i]=aRead;
980 return; << 874 }
981 } 875 }
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 876
1049 /////////////////////////////////////////// << 877 infile.close();
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 } 878 }
>> 879
>> 880 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1152 881