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