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