| 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 // 26 //
>> 27 // $Id: G4RDVeLowEnergyLoss.cc,v 1.25 2006/06/29 19:41:50 gunter Exp $
>> 28 // GEANT4 tag $Name: geant4-09-01-ref-00 $
27 // 29 //
28 // 30 //
29 // ------------------------------------------- 31 // --------------------------------------------------------------
30 // GEANT 4 class implementation file 32 // GEANT 4 class implementation file
31 // 33 //
32 // History: first implementation, based on ob 34 // History: first implementation, based on object model of
33 // 2nd December 1995, G.Cosmo 35 // 2nd December 1995, G.Cosmo
34 // ------------------------------------------- 36 // --------------------------------------------------------------
35 // 37 //
36 // Modifications: 38 // Modifications:
37 // 20/09/00 update fluctuations V.Ivanchenko 39 // 20/09/00 update fluctuations V.Ivanchenko
38 // 22/11/00 minor fix in fluctuations V.Ivanch 40 // 22/11/00 minor fix in fluctuations V.Ivanchenko
39 // 10/05/01 V.Ivanchenko Clean up againist Li 41 // 10/05/01 V.Ivanchenko Clean up againist Linux compilation with -Wall
40 // 22/05/01 V.Ivanchenko Update range calcula 42 // 22/05/01 V.Ivanchenko Update range calculation
41 // 23/11/01 V.Ivanchenko Move static member-f 43 // 23/11/01 V.Ivanchenko Move static member-functions from header to source
42 // 22/01/03 V.Ivanchenko Cut per region 44 // 22/01/03 V.Ivanchenko Cut per region
43 // 11/02/03 V.Ivanchenko Add limits to fluctu 45 // 11/02/03 V.Ivanchenko Add limits to fluctuations
44 // 24/04/03 V.Ivanchenko Fix the problem of t 46 // 24/04/03 V.Ivanchenko Fix the problem of table size
45 // 47 //
46 // ------------------------------------------- 48 // --------------------------------------------------------------
47 49
48 #include "G4RDVeLowEnergyLoss.hh" 50 #include "G4RDVeLowEnergyLoss.hh"
49 #include "G4PhysicalConstants.hh" <<
50 #include "G4SystemOfUnits.hh" <<
51 #include "G4ProductionCutsTable.hh" 51 #include "G4ProductionCutsTable.hh"
52 52
53 G4double G4RDVeLowEnergyLoss::ParticleMass 53 G4double G4RDVeLowEnergyLoss::ParticleMass ;
54 G4double G4RDVeLowEnergyLoss::taulow 54 G4double G4RDVeLowEnergyLoss::taulow ;
55 G4double G4RDVeLowEnergyLoss::tauhigh 55 G4double G4RDVeLowEnergyLoss::tauhigh ;
56 G4double G4RDVeLowEnergyLoss::ltaulow 56 G4double G4RDVeLowEnergyLoss::ltaulow ;
57 G4double G4RDVeLowEnergyLoss::ltauhigh 57 G4double G4RDVeLowEnergyLoss::ltauhigh ;
58 58
59 59
60 G4bool G4RDVeLowEnergyLoss::rndmStepFlag 60 G4bool G4RDVeLowEnergyLoss::rndmStepFlag = false;
61 G4bool G4RDVeLowEnergyLoss::EnlossFlucFl 61 G4bool G4RDVeLowEnergyLoss::EnlossFlucFlag = true;
62 G4double G4RDVeLowEnergyLoss::dRoverRange 62 G4double G4RDVeLowEnergyLoss::dRoverRange = 20*perCent;
63 G4double G4RDVeLowEnergyLoss::finalRange 63 G4double G4RDVeLowEnergyLoss::finalRange = 200*micrometer;
64 G4double G4RDVeLowEnergyLoss::c1lim = dRov 64 G4double G4RDVeLowEnergyLoss::c1lim = dRoverRange ;
65 G4double G4RDVeLowEnergyLoss::c2lim = 2.*( 65 G4double G4RDVeLowEnergyLoss::c2lim = 2.*(1.-dRoverRange)*finalRange ;
66 G4double G4RDVeLowEnergyLoss::c3lim = -(1. 66 G4double G4RDVeLowEnergyLoss::c3lim = -(1.-dRoverRange)*finalRange*finalRange;
67 67
68 68
69 // 69 //
70 70
71 G4RDVeLowEnergyLoss::G4RDVeLowEnergyLoss() 71 G4RDVeLowEnergyLoss::G4RDVeLowEnergyLoss()
72 :G4VContinuousDiscreteProce 72 :G4VContinuousDiscreteProcess("No Name Loss Process"),
73 lastMaterial(0), 73 lastMaterial(0),
74 nmaxCont1(4), 74 nmaxCont1(4),
75 nmaxCont2(16) 75 nmaxCont2(16)
76 { 76 {
77 G4Exception("G4RDVeLowEnergyLoss::G4RDVeLowE 77 G4Exception("G4RDVeLowEnergyLoss::G4RDVeLowEnergyLoss()", "InvalidCall",
78 FatalException, "Default constru 78 FatalException, "Default constructor is called!");
79 } 79 }
80 80
81 // 81 //
82 82
83 G4RDVeLowEnergyLoss::G4RDVeLowEnergyLoss(const 83 G4RDVeLowEnergyLoss::G4RDVeLowEnergyLoss(const G4String& aName, G4ProcessType aType)
84 : G4VContinuousDiscreteProce 84 : G4VContinuousDiscreteProcess(aName, aType),
85 lastMaterial(0), 85 lastMaterial(0),
86 nmaxCont1(4), 86 nmaxCont1(4),
87 nmaxCont2(16) 87 nmaxCont2(16)
88 { 88 {
89 } 89 }
90 90
91 // 91 //
92 92
93 G4RDVeLowEnergyLoss::~G4RDVeLowEnergyLoss() 93 G4RDVeLowEnergyLoss::~G4RDVeLowEnergyLoss()
94 { 94 {
95 } 95 }
96 96
97 // 97 //
98 98
99 G4RDVeLowEnergyLoss::G4RDVeLowEnergyLoss(G4RDV 99 G4RDVeLowEnergyLoss::G4RDVeLowEnergyLoss(G4RDVeLowEnergyLoss& right)
100 : G4VContinuousDiscreteProce 100 : G4VContinuousDiscreteProcess(right),
101 lastMaterial(0), 101 lastMaterial(0),
102 nmaxCont1(4), 102 nmaxCont1(4),
103 nmaxCont2(16) 103 nmaxCont2(16)
104 { 104 {
105 } 105 }
106 106
107 void G4RDVeLowEnergyLoss::SetRndmStep(G4bool v 107 void G4RDVeLowEnergyLoss::SetRndmStep(G4bool value)
108 { 108 {
109 rndmStepFlag = value; 109 rndmStepFlag = value;
110 } 110 }
111 111
112 void G4RDVeLowEnergyLoss::SetEnlossFluc(G4bool 112 void G4RDVeLowEnergyLoss::SetEnlossFluc(G4bool value)
113 { 113 {
114 EnlossFlucFlag = value; 114 EnlossFlucFlag = value;
115 } 115 }
116 116
117 void G4RDVeLowEnergyLoss::SetStepFunction (G4d 117 void G4RDVeLowEnergyLoss::SetStepFunction (G4double c1, G4double c2)
118 { 118 {
119 dRoverRange = c1; 119 dRoverRange = c1;
120 finalRange = c2; 120 finalRange = c2;
121 c1lim=dRoverRange; 121 c1lim=dRoverRange;
122 c2lim=2.*(1-dRoverRange)*finalRange; 122 c2lim=2.*(1-dRoverRange)*finalRange;
123 c3lim=-(1.-dRoverRange)*finalRange*finalRan 123 c3lim=-(1.-dRoverRange)*finalRange*finalRange;
124 } 124 }
125 125
126 G4PhysicsTable* G4RDVeLowEnergyLoss::BuildRang 126 G4PhysicsTable* G4RDVeLowEnergyLoss::BuildRangeTable(G4PhysicsTable* theDEDXTable,
127 G4PhysicsTable* theRangeTable, 127 G4PhysicsTable* theRangeTable,
128 G4double lowestKineticEnergy, 128 G4double lowestKineticEnergy,
129 G4double highestKineticEnergy, 129 G4double highestKineticEnergy,
130 G4int TotBin) 130 G4int TotBin)
131 // Build range table from the energy loss tabl 131 // Build range table from the energy loss table
132 { 132 {
133 133
134 G4int numOfCouples = theDEDXTable->length() 134 G4int numOfCouples = theDEDXTable->length();
135 135
136 if(theRangeTable) 136 if(theRangeTable)
137 { theRangeTable->clearAndDestroy(); 137 { theRangeTable->clearAndDestroy();
138 delete theRangeTable; } 138 delete theRangeTable; }
139 theRangeTable = new G4PhysicsTable(numOfCou 139 theRangeTable = new G4PhysicsTable(numOfCouples);
140 140
141 // loop for materials 141 // loop for materials
142 142
143 for (G4int J=0; J<numOfCouples; J++) 143 for (G4int J=0; J<numOfCouples; J++)
144 { 144 {
145 G4PhysicsLogVector* aVector; 145 G4PhysicsLogVector* aVector;
146 aVector = new G4PhysicsLogVector(lowestKi 146 aVector = new G4PhysicsLogVector(lowestKineticEnergy,
147 highestKineticEn 147 highestKineticEnergy,TotBin);
148 BuildRangeVector(theDEDXTable,lowestKinet 148 BuildRangeVector(theDEDXTable,lowestKineticEnergy,highestKineticEnergy,
149 TotBin,J,aVector); 149 TotBin,J,aVector);
150 theRangeTable->insert(aVector); 150 theRangeTable->insert(aVector);
151 } 151 }
152 return theRangeTable ; 152 return theRangeTable ;
153 } 153 }
154 154
155 // 155 //
156 156
157 void G4RDVeLowEnergyLoss::BuildRangeVector(G4P 157 void G4RDVeLowEnergyLoss::BuildRangeVector(G4PhysicsTable* theDEDXTable,
158 G4dou 158 G4double lowestKineticEnergy,
159 G4dou 159 G4double,
160 G4int 160 G4int TotBin,
161 G4int 161 G4int materialIndex,
162 G4Phy 162 G4PhysicsLogVector* rangeVector)
163 163
164 // create range vector for a material 164 // create range vector for a material
165 { 165 {
166 G4bool isOut; 166 G4bool isOut;
167 G4PhysicsVector* physicsVector= (*theDEDXTab 167 G4PhysicsVector* physicsVector= (*theDEDXTable)[materialIndex];
168 G4double energy1 = lowestKineticEnergy; 168 G4double energy1 = lowestKineticEnergy;
169 G4double dedx = physicsVector->GetValue(e 169 G4double dedx = physicsVector->GetValue(energy1,isOut);
170 G4double range = 0.5*energy1/dedx; 170 G4double range = 0.5*energy1/dedx;
171 rangeVector->PutValue(0,range); 171 rangeVector->PutValue(0,range);
172 G4int n = 100; 172 G4int n = 100;
173 G4double del = 1.0/(G4double)n ; 173 G4double del = 1.0/(G4double)n ;
174 174
175 for (G4int j=1; j<TotBin; j++) { 175 for (G4int j=1; j<TotBin; j++) {
176 176
177 G4double energy2 = rangeVector->GetLowEdge 177 G4double energy2 = rangeVector->GetLowEdgeEnergy(j);
178 G4double de = (energy2 - energy1) * del ; 178 G4double de = (energy2 - energy1) * del ;
179 G4double dedx1 = dedx ; 179 G4double dedx1 = dedx ;
180 180
181 for (G4int i=1; i<n; i++) { 181 for (G4int i=1; i<n; i++) {
182 G4double energy = energy1 + i*de ; 182 G4double energy = energy1 + i*de ;
183 G4double dedx2 = physicsVector->GetValu 183 G4double dedx2 = physicsVector->GetValue(energy,isOut);
184 range += 0.5*de*(1.0/dedx1 + 1.0/dedx2) 184 range += 0.5*de*(1.0/dedx1 + 1.0/dedx2);
185 dedx1 = dedx2; 185 dedx1 = dedx2;
186 } 186 }
187 rangeVector->PutValue(j,range); 187 rangeVector->PutValue(j,range);
188 dedx = dedx1 ; 188 dedx = dedx1 ;
189 energy1 = energy2 ; 189 energy1 = energy2 ;
190 } 190 }
191 } 191 }
192 192
193 // 193 //
194 194
195 G4double G4RDVeLowEnergyLoss::RangeIntLin(G4Ph 195 G4double G4RDVeLowEnergyLoss::RangeIntLin(G4PhysicsVector* physicsVector,
196 G4int nbin) 196 G4int nbin)
197 // num. integration, linear binning 197 // num. integration, linear binning
198 { 198 {
199 G4double dtau,Value,taui,ti,lossi,ci; 199 G4double dtau,Value,taui,ti,lossi,ci;
200 G4bool isOut; 200 G4bool isOut;
201 dtau = (tauhigh-taulow)/nbin; 201 dtau = (tauhigh-taulow)/nbin;
202 Value = 0.; 202 Value = 0.;
203 203
204 for (G4int i=0; i<=nbin; i++) 204 for (G4int i=0; i<=nbin; i++)
205 { 205 {
206 taui = taulow + dtau*i ; 206 taui = taulow + dtau*i ;
207 ti = ParticleMass*taui; 207 ti = ParticleMass*taui;
208 lossi = physicsVector->GetValue(ti,isOut); 208 lossi = physicsVector->GetValue(ti,isOut);
209 if(i==0) 209 if(i==0)
210 ci=0.5; 210 ci=0.5;
211 else 211 else
212 { 212 {
213 if(i<nbin) 213 if(i<nbin)
214 ci=1.; 214 ci=1.;
215 else 215 else
216 ci=0.5; 216 ci=0.5;
217 } 217 }
218 Value += ci/lossi; 218 Value += ci/lossi;
219 } 219 }
220 Value *= ParticleMass*dtau; 220 Value *= ParticleMass*dtau;
221 return Value; 221 return Value;
222 } 222 }
223 223
224 // 224 //
225 225
226 G4double G4RDVeLowEnergyLoss::RangeIntLog(G4Ph 226 G4double G4RDVeLowEnergyLoss::RangeIntLog(G4PhysicsVector* physicsVector,
227 G4int nbin) 227 G4int nbin)
228 // num. integration, logarithmic binning 228 // num. integration, logarithmic binning
229 { 229 {
230 G4double ltt,dltau,Value,ui,taui,ti,lossi,ci 230 G4double ltt,dltau,Value,ui,taui,ti,lossi,ci;
231 G4bool isOut; 231 G4bool isOut;
232 ltt = ltauhigh-ltaulow; 232 ltt = ltauhigh-ltaulow;
233 dltau = ltt/nbin; 233 dltau = ltt/nbin;
234 Value = 0.; 234 Value = 0.;
235 235
236 for (G4int i=0; i<=nbin; i++) 236 for (G4int i=0; i<=nbin; i++)
237 { 237 {
238 ui = ltaulow+dltau*i; 238 ui = ltaulow+dltau*i;
239 taui = std::exp(ui); 239 taui = std::exp(ui);
240 ti = ParticleMass*taui; 240 ti = ParticleMass*taui;
241 lossi = physicsVector->GetValue(ti,isOut); 241 lossi = physicsVector->GetValue(ti,isOut);
242 if(i==0) 242 if(i==0)
243 ci=0.5; 243 ci=0.5;
244 else 244 else
245 { 245 {
246 if(i<nbin) 246 if(i<nbin)
247 ci=1.; 247 ci=1.;
248 else 248 else
249 ci=0.5; 249 ci=0.5;
250 } 250 }
251 Value += ci*taui/lossi; 251 Value += ci*taui/lossi;
252 } 252 }
253 Value *= ParticleMass*dltau; 253 Value *= ParticleMass*dltau;
254 return Value; 254 return Value;
255 } 255 }
256 256
257 257
258 // 258 //
259 259
260 G4PhysicsTable* G4RDVeLowEnergyLoss::BuildLabT 260 G4PhysicsTable* G4RDVeLowEnergyLoss::BuildLabTimeTable(G4PhysicsTable* theDEDXTable,
261 G4PhysicsTable* theLabTimeTab 261 G4PhysicsTable* theLabTimeTable,
262 G4double lowestKineticEnergy, 262 G4double lowestKineticEnergy,
263 G4double highestKineticEnergy 263 G4double highestKineticEnergy,G4int TotBin)
264 264
265 { 265 {
266 266
267 G4int numOfCouples = G4ProductionCutsTable:: 267 G4int numOfCouples = G4ProductionCutsTable::GetProductionCutsTable()->GetTableSize();
268 268
269 if(theLabTimeTable) 269 if(theLabTimeTable)
270 { theLabTimeTable->clearAndDestroy(); 270 { theLabTimeTable->clearAndDestroy();
271 delete theLabTimeTable; } 271 delete theLabTimeTable; }
272 theLabTimeTable = new G4PhysicsTable(numOfCo 272 theLabTimeTable = new G4PhysicsTable(numOfCouples);
273 273
274 274
275 for (G4int J=0; J<numOfCouples; J++) 275 for (G4int J=0; J<numOfCouples; J++)
276 { 276 {
277 G4PhysicsLogVector* aVector; 277 G4PhysicsLogVector* aVector;
278 278
279 aVector = new G4PhysicsLogVector(lowestKin 279 aVector = new G4PhysicsLogVector(lowestKineticEnergy,
280 highestKineticEner 280 highestKineticEnergy,TotBin);
281 281
282 BuildLabTimeVector(theDEDXTable, 282 BuildLabTimeVector(theDEDXTable,
283 lowestKineticEnergy,highestKinet 283 lowestKineticEnergy,highestKineticEnergy,TotBin,J,aVector);
284 theLabTimeTable->insert(aVector); 284 theLabTimeTable->insert(aVector);
285 285
286 286
287 } 287 }
288 return theLabTimeTable ; 288 return theLabTimeTable ;
289 } 289 }
290 290
291 // 291 //
292 292
293 G4PhysicsTable* G4RDVeLowEnergyLoss::BuildProp 293 G4PhysicsTable* G4RDVeLowEnergyLoss::BuildProperTimeTable(G4PhysicsTable* theDEDXTable,
294 G4PhysicsTable* theProperTimeTab 294 G4PhysicsTable* theProperTimeTable,
295 G4double lowestKineticEnergy, 295 G4double lowestKineticEnergy,
296 G4double highestKineticEnergy,G4 296 G4double highestKineticEnergy,G4int TotBin)
297 297
298 { 298 {
299 299
300 G4int numOfCouples = G4ProductionCutsTable:: 300 G4int numOfCouples = G4ProductionCutsTable::GetProductionCutsTable()->GetTableSize();
301 301
302 if(theProperTimeTable) 302 if(theProperTimeTable)
303 { theProperTimeTable->clearAndDestroy(); 303 { theProperTimeTable->clearAndDestroy();
304 delete theProperTimeTable; } 304 delete theProperTimeTable; }
305 theProperTimeTable = new G4PhysicsTable(numO 305 theProperTimeTable = new G4PhysicsTable(numOfCouples);
306 306
307 307
308 for (G4int J=0; J<numOfCouples; J++) 308 for (G4int J=0; J<numOfCouples; J++)
309 { 309 {
310 G4PhysicsLogVector* aVector; 310 G4PhysicsLogVector* aVector;
311 311
312 aVector = new G4PhysicsLogVector(lowestKin 312 aVector = new G4PhysicsLogVector(lowestKineticEnergy,
313 highestKineticEner 313 highestKineticEnergy,TotBin);
314 314
315 BuildProperTimeVector(theDEDXTable, 315 BuildProperTimeVector(theDEDXTable,
316 lowestKineticEnergy,highestKinet 316 lowestKineticEnergy,highestKineticEnergy,TotBin,J,aVector);
317 theProperTimeTable->insert(aVector); 317 theProperTimeTable->insert(aVector);
318 318
319 319
320 } 320 }
321 return theProperTimeTable ; 321 return theProperTimeTable ;
322 } 322 }
323 323
324 // 324 //
325 325
326 void G4RDVeLowEnergyLoss::BuildLabTimeVector(G 326 void G4RDVeLowEnergyLoss::BuildLabTimeVector(G4PhysicsTable* theDEDXTable,
327 G4double, // lowestKineticEnergy, 327 G4double, // lowestKineticEnergy,
328 G4double, // highestKineticEnergy 328 G4double, // highestKineticEnergy,
329 G4i 329 G4int TotBin,
330 G4int materialIndex, G4PhysicsLog 330 G4int materialIndex, G4PhysicsLogVector* timeVector)
331 // create lab time vector for a material 331 // create lab time vector for a material
332 { 332 {
333 333
334 G4int nbin=100; 334 G4int nbin=100;
335 G4bool isOut; 335 G4bool isOut;
336 G4double tlim=5.*keV,parlowen=0.4,ppar=0.5-p 336 G4double tlim=5.*keV,parlowen=0.4,ppar=0.5-parlowen ;
337 G4double losslim,clim,taulim,timelim, 337 G4double losslim,clim,taulim,timelim,
338 LowEdgeEnergy,tau,Value ; 338 LowEdgeEnergy,tau,Value ;
339 339
340 G4PhysicsVector* physicsVector= (*theDEDXTab 340 G4PhysicsVector* physicsVector= (*theDEDXTable)[materialIndex];
341 341
342 // low energy part first... 342 // low energy part first...
343 losslim = physicsVector->GetValue(tlim,isOut 343 losslim = physicsVector->GetValue(tlim,isOut);
344 taulim=tlim/ParticleMass ; 344 taulim=tlim/ParticleMass ;
345 clim=std::sqrt(ParticleMass*tlim/2.)/(c_ligh 345 clim=std::sqrt(ParticleMass*tlim/2.)/(c_light*losslim*ppar) ;
346 346
347 G4int i=-1; 347 G4int i=-1;
348 G4double oldValue = 0. ; 348 G4double oldValue = 0. ;
349 G4double tauold ; 349 G4double tauold ;
350 do 350 do
351 { 351 {
352 i += 1 ; 352 i += 1 ;
353 LowEdgeEnergy = timeVector->GetLowEdgeEner 353 LowEdgeEnergy = timeVector->GetLowEdgeEnergy(i);
354 tau = LowEdgeEnergy/ParticleMass ; 354 tau = LowEdgeEnergy/ParticleMass ;
355 if ( tau <= taulim ) 355 if ( tau <= taulim )
356 { 356 {
357 Value = clim*std::exp(ppar*std::log(tau/ 357 Value = clim*std::exp(ppar*std::log(tau/taulim)) ;
358 } 358 }
359 else 359 else
360 { 360 {
361 timelim=clim ; 361 timelim=clim ;
362 ltaulow = std::log(taulim); 362 ltaulow = std::log(taulim);
363 ltauhigh = std::log(tau); 363 ltauhigh = std::log(tau);
364 Value = timelim+LabTimeIntLog(physicsVec 364 Value = timelim+LabTimeIntLog(physicsVector,nbin);
365 } 365 }
366 timeVector->PutValue(i,Value); 366 timeVector->PutValue(i,Value);
367 oldValue = Value ; 367 oldValue = Value ;
368 tauold = tau ; 368 tauold = tau ;
369 } while (tau<=taulim) ; 369 } while (tau<=taulim) ;
370 i += 1 ; 370 i += 1 ;
371 for (G4int j=i; j<TotBin; j++) 371 for (G4int j=i; j<TotBin; j++)
372 { 372 {
373 LowEdgeEnergy = timeVector->GetLowEdgeEner 373 LowEdgeEnergy = timeVector->GetLowEdgeEnergy(j);
374 tau = LowEdgeEnergy/ParticleMass ; 374 tau = LowEdgeEnergy/ParticleMass ;
375 ltaulow = std::log(tauold); 375 ltaulow = std::log(tauold);
376 ltauhigh = std::log(tau); 376 ltauhigh = std::log(tau);
377 Value = oldValue+LabTimeIntLog(physicsVect 377 Value = oldValue+LabTimeIntLog(physicsVector,nbin);
378 timeVector->PutValue(j,Value); 378 timeVector->PutValue(j,Value);
379 oldValue = Value ; 379 oldValue = Value ;
380 tauold = tau ; 380 tauold = tau ;
381 } 381 }
382 } 382 }
383 383
384 // 384 //
385 385
386 void G4RDVeLowEnergyLoss::BuildProperTimeVecto 386 void G4RDVeLowEnergyLoss::BuildProperTimeVector(G4PhysicsTable* theDEDXTable,
387 G4double, // lowestKineticEner 387 G4double, // lowestKineticEnergy,
388 G4double, // highestKineticEne 388 G4double, // highestKineticEnergy,
389 389 G4int TotBin,
390 G4int materialIndex, G4Physics 390 G4int materialIndex, G4PhysicsLogVector* timeVector)
391 // create proper time vector for a material 391 // create proper time vector for a material
392 { 392 {
393 G4int nbin=100; 393 G4int nbin=100;
394 G4bool isOut; 394 G4bool isOut;
395 G4double tlim=5.*keV,parlowen=0.4,ppar=0.5-p 395 G4double tlim=5.*keV,parlowen=0.4,ppar=0.5-parlowen ;
396 G4double losslim,clim,taulim,timelim, 396 G4double losslim,clim,taulim,timelim,
397 LowEdgeEnergy,tau,Value ; 397 LowEdgeEnergy,tau,Value ;
398 398
399 G4PhysicsVector* physicsVector= (*theDEDXTab 399 G4PhysicsVector* physicsVector= (*theDEDXTable)[materialIndex];
400 //const G4MaterialTable* theMaterialTable = 400 //const G4MaterialTable* theMaterialTable = G4Material::GetMaterialTable();
401 401
402 // low energy part first... 402 // low energy part first...
403 losslim = physicsVector->GetValue(tlim,isOut 403 losslim = physicsVector->GetValue(tlim,isOut);
404 taulim=tlim/ParticleMass ; 404 taulim=tlim/ParticleMass ;
405 clim=std::sqrt(ParticleMass*tlim/2.)/(c_ligh 405 clim=std::sqrt(ParticleMass*tlim/2.)/(c_light*losslim*ppar) ;
406 406
407 G4int i=-1; 407 G4int i=-1;
408 G4double oldValue = 0. ; 408 G4double oldValue = 0. ;
409 G4double tauold ; 409 G4double tauold ;
410 do 410 do
411 { 411 {
412 i += 1 ; 412 i += 1 ;
413 LowEdgeEnergy = timeVector->GetLowEdgeEner 413 LowEdgeEnergy = timeVector->GetLowEdgeEnergy(i);
414 tau = LowEdgeEnergy/ParticleMass ; 414 tau = LowEdgeEnergy/ParticleMass ;
415 if ( tau <= taulim ) 415 if ( tau <= taulim )
416 { 416 {
417 Value = clim*std::exp(ppar*std::log(tau/ 417 Value = clim*std::exp(ppar*std::log(tau/taulim)) ;
418 } 418 }
419 else 419 else
420 { 420 {
421 timelim=clim ; 421 timelim=clim ;
422 ltaulow = std::log(taulim); 422 ltaulow = std::log(taulim);
423 ltauhigh = std::log(tau); 423 ltauhigh = std::log(tau);
424 Value = timelim+ProperTimeIntLog(physics 424 Value = timelim+ProperTimeIntLog(physicsVector,nbin);
425 } 425 }
426 timeVector->PutValue(i,Value); 426 timeVector->PutValue(i,Value);
427 oldValue = Value ; 427 oldValue = Value ;
428 tauold = tau ; 428 tauold = tau ;
429 } while (tau<=taulim) ; 429 } while (tau<=taulim) ;
430 i += 1 ; 430 i += 1 ;
431 for (G4int j=i; j<TotBin; j++) 431 for (G4int j=i; j<TotBin; j++)
432 { 432 {
433 LowEdgeEnergy = timeVector->GetLowEdgeEner 433 LowEdgeEnergy = timeVector->GetLowEdgeEnergy(j);
434 tau = LowEdgeEnergy/ParticleMass ; 434 tau = LowEdgeEnergy/ParticleMass ;
435 ltaulow = std::log(tauold); 435 ltaulow = std::log(tauold);
436 ltauhigh = std::log(tau); 436 ltauhigh = std::log(tau);
437 Value = oldValue+ProperTimeIntLog(physicsV 437 Value = oldValue+ProperTimeIntLog(physicsVector,nbin);
438 timeVector->PutValue(j,Value); 438 timeVector->PutValue(j,Value);
439 oldValue = Value ; 439 oldValue = Value ;
440 tauold = tau ; 440 tauold = tau ;
441 } 441 }
442 } 442 }
443 443
444 // 444 //
445 445
446 G4double G4RDVeLowEnergyLoss::LabTimeIntLog(G4 446 G4double G4RDVeLowEnergyLoss::LabTimeIntLog(G4PhysicsVector* physicsVector,
447 G4int nbin) 447 G4int nbin)
448 // num. integration, logarithmic binning 448 // num. integration, logarithmic binning
449 { 449 {
450 G4double ltt,dltau,Value,ui,taui,ti,lossi,ci 450 G4double ltt,dltau,Value,ui,taui,ti,lossi,ci;
451 G4bool isOut; 451 G4bool isOut;
452 ltt = ltauhigh-ltaulow; 452 ltt = ltauhigh-ltaulow;
453 dltau = ltt/nbin; 453 dltau = ltt/nbin;
454 Value = 0.; 454 Value = 0.;
455 455
456 for (G4int i=0; i<=nbin; i++) 456 for (G4int i=0; i<=nbin; i++)
457 { 457 {
458 ui = ltaulow+dltau*i; 458 ui = ltaulow+dltau*i;
459 taui = std::exp(ui); 459 taui = std::exp(ui);
460 ti = ParticleMass*taui; 460 ti = ParticleMass*taui;
461 lossi = physicsVector->GetValue(ti,isOut); 461 lossi = physicsVector->GetValue(ti,isOut);
462 if(i==0) 462 if(i==0)
463 ci=0.5; 463 ci=0.5;
464 else 464 else
465 { 465 {
466 if(i<nbin) 466 if(i<nbin)
467 ci=1.; 467 ci=1.;
468 else 468 else
469 ci=0.5; 469 ci=0.5;
470 } 470 }
471 Value += ci*taui*(ti+ParticleMass)/(std::s 471 Value += ci*taui*(ti+ParticleMass)/(std::sqrt(ti*(ti+2.*ParticleMass))*lossi);
472 } 472 }
473 Value *= ParticleMass*dltau/c_light; 473 Value *= ParticleMass*dltau/c_light;
474 return Value; 474 return Value;
475 } 475 }
476 476
477 // 477 //
478 478
479 G4double G4RDVeLowEnergyLoss::ProperTimeIntLog 479 G4double G4RDVeLowEnergyLoss::ProperTimeIntLog(G4PhysicsVector* physicsVector,
480 G4int nbin) 480 G4int nbin)
481 // num. integration, logarithmic binning 481 // num. integration, logarithmic binning
482 { 482 {
483 G4double ltt,dltau,Value,ui,taui,ti,lossi,ci 483 G4double ltt,dltau,Value,ui,taui,ti,lossi,ci;
484 G4bool isOut; 484 G4bool isOut;
485 ltt = ltauhigh-ltaulow; 485 ltt = ltauhigh-ltaulow;
486 dltau = ltt/nbin; 486 dltau = ltt/nbin;
487 Value = 0.; 487 Value = 0.;
488 488
489 for (G4int i=0; i<=nbin; i++) 489 for (G4int i=0; i<=nbin; i++)
490 { 490 {
491 ui = ltaulow+dltau*i; 491 ui = ltaulow+dltau*i;
492 taui = std::exp(ui); 492 taui = std::exp(ui);
493 ti = ParticleMass*taui; 493 ti = ParticleMass*taui;
494 lossi = physicsVector->GetValue(ti,isOut); 494 lossi = physicsVector->GetValue(ti,isOut);
495 if(i==0) 495 if(i==0)
496 ci=0.5; 496 ci=0.5;
497 else 497 else
498 { 498 {
499 if(i<nbin) 499 if(i<nbin)
500 ci=1.; 500 ci=1.;
501 else 501 else
502 ci=0.5; 502 ci=0.5;
503 } 503 }
504 Value += ci*taui*ParticleMass/(std::sqrt(t 504 Value += ci*taui*ParticleMass/(std::sqrt(ti*(ti+2.*ParticleMass))*lossi);
505 } 505 }
506 Value *= ParticleMass*dltau/c_light; 506 Value *= ParticleMass*dltau/c_light;
507 return Value; 507 return Value;
508 } 508 }
509 509
510 // 510 //
511 511
512 G4PhysicsTable* G4RDVeLowEnergyLoss::BuildInve 512 G4PhysicsTable* G4RDVeLowEnergyLoss::BuildInverseRangeTable(G4PhysicsTable* theRangeTable,
513 G4PhysicsTable*, 513 G4PhysicsTable*,
514 G4PhysicsTable*, 514 G4PhysicsTable*,
515 G4PhysicsTable*, 515 G4PhysicsTable*,
516 G4PhysicsTable* theInverseRang 516 G4PhysicsTable* theInverseRangeTable,
517 G4double, // lowestKineticEner 517 G4double, // lowestKineticEnergy,
518 G4double, // highestKineticEne 518 G4double, // highestKineticEnergy
519 G4int ) // nbins 519 G4int ) // nbins
520 // Build inverse table of the range table 520 // Build inverse table of the range table
521 { 521 {
522 G4bool b; 522 G4bool b;
523 523
524 G4int numOfCouples = G4ProductionCutsTable:: 524 G4int numOfCouples = G4ProductionCutsTable::GetProductionCutsTable()->GetTableSize();
525 525
526 if(theInverseRangeTable) 526 if(theInverseRangeTable)
527 { theInverseRangeTable->clearAndDestroy(); 527 { theInverseRangeTable->clearAndDestroy();
528 delete theInverseRangeTable; } 528 delete theInverseRangeTable; }
529 theInverseRangeTable = new G4PhysicsTable( 529 theInverseRangeTable = new G4PhysicsTable(numOfCouples);
530 530
531 // loop for materials 531 // loop for materials
532 for (G4int i=0; i<numOfCouples; i++) 532 for (G4int i=0; i<numOfCouples; i++)
533 { 533 {
534 534
535 G4PhysicsVector* pv = (*theRangeTable)[i]; 535 G4PhysicsVector* pv = (*theRangeTable)[i];
536 size_t nbins = pv->GetVectorLength(); 536 size_t nbins = pv->GetVectorLength();
537 G4double elow = pv->GetLowEdgeEnergy(0); 537 G4double elow = pv->GetLowEdgeEnergy(0);
538 G4double ehigh = pv->GetLowEdgeEnergy(nbin 538 G4double ehigh = pv->GetLowEdgeEnergy(nbins-1);
539 G4double rlow = pv->GetValue(elow, b); 539 G4double rlow = pv->GetValue(elow, b);
540 G4double rhigh = pv->GetValue(ehigh, b); 540 G4double rhigh = pv->GetValue(ehigh, b);
541 541
542 rhigh *= std::exp(std::log(rhigh/rlow)/((G 542 rhigh *= std::exp(std::log(rhigh/rlow)/((G4double)(nbins-1)));
543 543
544 G4PhysicsLogVector* v = new G4PhysicsLogVe 544 G4PhysicsLogVector* v = new G4PhysicsLogVector(rlow, rhigh, nbins);
545 545
546 v->PutValue(0,elow); 546 v->PutValue(0,elow);
547 G4double energy1 = elow; 547 G4double energy1 = elow;
548 G4double range1 = rlow; 548 G4double range1 = rlow;
549 G4double energy2 = elow; 549 G4double energy2 = elow;
550 G4double range2 = rlow; 550 G4double range2 = rlow;
551 size_t ilow = 0; 551 size_t ilow = 0;
552 size_t ihigh; 552 size_t ihigh;
553 553
554 for (size_t j=1; j<nbins; j++) { 554 for (size_t j=1; j<nbins; j++) {
555 555
556 G4double range = v->GetLowEdgeEnergy(j); 556 G4double range = v->GetLowEdgeEnergy(j);
557 557
558 for (ihigh=ilow+1; ihigh<nbins; ihigh++) 558 for (ihigh=ilow+1; ihigh<nbins; ihigh++) {
559 energy2 = pv->GetLowEdgeEnergy(ihigh); 559 energy2 = pv->GetLowEdgeEnergy(ihigh);
560 range2 = pv->GetValue(energy2, b); 560 range2 = pv->GetValue(energy2, b);
561 if(range2 >= range || ihigh == nbins-1 561 if(range2 >= range || ihigh == nbins-1) {
562 ilow = ihigh - 1; 562 ilow = ihigh - 1;
563 energy1 = pv->GetLowEdgeEnergy(ilow) 563 energy1 = pv->GetLowEdgeEnergy(ilow);
564 range1 = pv->GetValue(energy1, b); 564 range1 = pv->GetValue(energy1, b);
565 break; 565 break;
566 } 566 }
567 } 567 }
568 568
569 G4double e = std::log(energy1) + std::lo 569 G4double e = std::log(energy1) + std::log(energy2/energy1)*std::log(range/range1)/std::log(range2/range1);
570 570
571 v->PutValue(j,std::exp(e)); 571 v->PutValue(j,std::exp(e));
572 } 572 }
573 theInverseRangeTable->insert(v); 573 theInverseRangeTable->insert(v);
574 574
575 } 575 }
576 return theInverseRangeTable ; 576 return theInverseRangeTable ;
577 } 577 }
578 578
579 // 579 //
580 580
581 void G4RDVeLowEnergyLoss::InvertRangeVector(G4 581 void G4RDVeLowEnergyLoss::InvertRangeVector(G4PhysicsTable* theRangeTable,
582 G4PhysicsTable* theRangeCoeffATabl 582 G4PhysicsTable* theRangeCoeffATable,
583 G4PhysicsTable* theRangeCoeffBTabl 583 G4PhysicsTable* theRangeCoeffBTable,
584 G4PhysicsTable* theRangeCoeffCTabl 584 G4PhysicsTable* theRangeCoeffCTable,
585 G4double lowestKineticEnergy, 585 G4double lowestKineticEnergy,
586 G4double highestKineticEnergy, G4i 586 G4double highestKineticEnergy, G4int TotBin,
587 G4int materialIndex, G4PhysicsLog 587 G4int materialIndex, G4PhysicsLogVector* aVector)
588 // invert range vector for a material 588 // invert range vector for a material
589 { 589 {
590 G4double LowEdgeRange,A,B,C,discr,KineticEne 590 G4double LowEdgeRange,A,B,C,discr,KineticEnergy ;
591 G4double RTable = std::exp(std::log(highestK 591 G4double RTable = std::exp(std::log(highestKineticEnergy/lowestKineticEnergy)/TotBin) ;
592 G4double Tbin = lowestKineticEnergy/RTable ; 592 G4double Tbin = lowestKineticEnergy/RTable ;
593 G4double rangebin = 0.0 ; 593 G4double rangebin = 0.0 ;
594 G4int binnumber = -1 ; 594 G4int binnumber = -1 ;
595 G4bool isOut ; 595 G4bool isOut ;
596 596
597 //loop for range values 597 //loop for range values
598 for( G4int i=0; i<TotBin; i++) 598 for( G4int i=0; i<TotBin; i++)
599 { 599 {
600 LowEdgeRange = aVector->GetLowEdgeEnergy(i 600 LowEdgeRange = aVector->GetLowEdgeEnergy(i) ; //i.e. GetLowEdgeValue(i)
601 if( rangebin < LowEdgeRange ) 601 if( rangebin < LowEdgeRange )
602 { 602 {
603 do 603 do
604 { 604 {
605 binnumber += 1 ; 605 binnumber += 1 ;
606 Tbin *= RTable ; 606 Tbin *= RTable ;
607 rangebin = (*theRangeTable)(materialIn 607 rangebin = (*theRangeTable)(materialIndex)->GetValue(Tbin,isOut) ;
608 } 608 }
609 while ((rangebin < LowEdgeRange) && (bin 609 while ((rangebin < LowEdgeRange) && (binnumber < TotBin )) ;
610 } 610 }
611 611
612 if(binnumber == 0) 612 if(binnumber == 0)
613 KineticEnergy = lowestKineticEnergy ; 613 KineticEnergy = lowestKineticEnergy ;
614 else if(binnumber == TotBin-1) 614 else if(binnumber == TotBin-1)
615 KineticEnergy = highestKineticEnergy ; 615 KineticEnergy = highestKineticEnergy ;
616 else 616 else
617 { 617 {
618 A = (*(*theRangeCoeffATable)(materialInd 618 A = (*(*theRangeCoeffATable)(materialIndex))(binnumber-1) ;
619 B = (*(*theRangeCoeffBTable)(materialInd 619 B = (*(*theRangeCoeffBTable)(materialIndex))(binnumber-1) ;
620 C = (*(*theRangeCoeffCTable)(materialInd 620 C = (*(*theRangeCoeffCTable)(materialIndex))(binnumber-1) ;
621 if(A==0.) 621 if(A==0.)
622 KineticEnergy = (LowEdgeRange -C )/B 622 KineticEnergy = (LowEdgeRange -C )/B ;
623 else 623 else
624 { 624 {
625 discr = B*B - 4.*A*(C-LowEdgeRange); 625 discr = B*B - 4.*A*(C-LowEdgeRange);
626 discr = discr>0. ? std::sqrt(discr) : 626 discr = discr>0. ? std::sqrt(discr) : 0.;
627 KineticEnergy = 0.5*(discr-B)/A ; 627 KineticEnergy = 0.5*(discr-B)/A ;
628 } 628 }
629 } 629 }
630 630
631 aVector->PutValue(i,KineticEnergy) ; 631 aVector->PutValue(i,KineticEnergy) ;
632 } 632 }
633 } 633 }
634 634
635 // 635 //
636 636
637 G4PhysicsTable* G4RDVeLowEnergyLoss::BuildRang 637 G4PhysicsTable* G4RDVeLowEnergyLoss::BuildRangeCoeffATable(G4PhysicsTable* theRangeTable,
638 G4PhysicsTable* theRangeCoeffAT 638 G4PhysicsTable* theRangeCoeffATable,
639 G4double lowestKineticEnergy, 639 G4double lowestKineticEnergy,
640 G4double highestKineticEnergy, 640 G4double highestKineticEnergy, G4int TotBin)
641 // Build tables of coefficients for the energy 641 // Build tables of coefficients for the energy loss calculation
642 // create table for coefficients "A" 642 // create table for coefficients "A"
643 { 643 {
644 644
645 G4int numOfCouples = G4ProductionCutsTable:: 645 G4int numOfCouples = G4ProductionCutsTable::GetProductionCutsTable()->GetTableSize();
646 646
647 if(theRangeCoeffATable) 647 if(theRangeCoeffATable)
648 { theRangeCoeffATable->clearAndDestroy(); 648 { theRangeCoeffATable->clearAndDestroy();
649 delete theRangeCoeffATable; } 649 delete theRangeCoeffATable; }
650 theRangeCoeffATable = new G4PhysicsTable(num 650 theRangeCoeffATable = new G4PhysicsTable(numOfCouples);
651 651
652 G4double RTable = std::exp(std::log(highestK 652 G4double RTable = std::exp(std::log(highestKineticEnergy/lowestKineticEnergy)/TotBin) ;
653 G4double R2 = RTable*RTable ; 653 G4double R2 = RTable*RTable ;
654 G4double R1 = RTable+1.; 654 G4double R1 = RTable+1.;
655 G4double w = R1*(RTable-1.)*(RTable-1.); 655 G4double w = R1*(RTable-1.)*(RTable-1.);
656 G4double w1 = RTable/w , w2 = -RTable*R1/w , 656 G4double w1 = RTable/w , w2 = -RTable*R1/w , w3 = R2/w ;
657 G4double Ti , Tim , Tip , Ri , Rim , Rip , V 657 G4double Ti , Tim , Tip , Ri , Rim , Rip , Value ;
658 G4bool isOut; 658 G4bool isOut;
659 659
660 // loop for materials 660 // loop for materials
661 for (G4int J=0; J<numOfCouples; J++) 661 for (G4int J=0; J<numOfCouples; J++)
662 { 662 {
663 G4int binmax=TotBin ; 663 G4int binmax=TotBin ;
664 G4PhysicsLinearVector* aVector = 664 G4PhysicsLinearVector* aVector =
665 new G4PhysicsLinear 665 new G4PhysicsLinearVector(0.,binmax, TotBin);
666 Ti = lowestKineticEnergy ; 666 Ti = lowestKineticEnergy ;
667 G4PhysicsVector* rangeVector= (*theRangeTa 667 G4PhysicsVector* rangeVector= (*theRangeTable)[J];
668 668
669 for ( G4int i=0; i<TotBin; i++) 669 for ( G4int i=0; i<TotBin; i++)
670 { 670 {
671 Ri = rangeVector->GetValue(Ti,isOut) ; 671 Ri = rangeVector->GetValue(Ti,isOut) ;
672 if ( i==0 ) 672 if ( i==0 )
673 Rim = 0. ; 673 Rim = 0. ;
674 else 674 else
675 { 675 {
676 Tim = Ti/RTable ; 676 Tim = Ti/RTable ;
677 Rim = rangeVector->GetValue(Tim,isOut) 677 Rim = rangeVector->GetValue(Tim,isOut);
678 } 678 }
679 if ( i==(TotBin-1)) 679 if ( i==(TotBin-1))
680 Rip = Ri ; 680 Rip = Ri ;
681 else 681 else
682 { 682 {
683 Tip = Ti*RTable ; 683 Tip = Ti*RTable ;
684 Rip = rangeVector->GetValue(Tip,isOut) 684 Rip = rangeVector->GetValue(Tip,isOut);
685 } 685 }
686 Value = (w1*Rip + w2*Ri + w3*Rim)/(Ti*Ti 686 Value = (w1*Rip + w2*Ri + w3*Rim)/(Ti*Ti) ;
687 687
688 aVector->PutValue(i,Value); 688 aVector->PutValue(i,Value);
689 Ti = RTable*Ti ; 689 Ti = RTable*Ti ;
690 } 690 }
691 691
692 theRangeCoeffATable->insert(aVector); 692 theRangeCoeffATable->insert(aVector);
693 } 693 }
694 return theRangeCoeffATable ; 694 return theRangeCoeffATable ;
695 } 695 }
696 696
697 // 697 //
698 698
699 G4PhysicsTable* G4RDVeLowEnergyLoss::BuildRang 699 G4PhysicsTable* G4RDVeLowEnergyLoss::BuildRangeCoeffBTable(G4PhysicsTable* theRangeTable,
700 G4PhysicsTable* theRangeCoeffBT 700 G4PhysicsTable* theRangeCoeffBTable,
701 G4double lowestKineticEnergy, 701 G4double lowestKineticEnergy,
702 G4double highestKineticEnergy, 702 G4double highestKineticEnergy, G4int TotBin)
703 // Build tables of coefficients for the energy 703 // Build tables of coefficients for the energy loss calculation
704 // create table for coefficients "B" 704 // create table for coefficients "B"
705 { 705 {
706 706
707 G4int numOfCouples = G4ProductionCutsTable:: 707 G4int numOfCouples = G4ProductionCutsTable::GetProductionCutsTable()->GetTableSize();
708 708
709 if(theRangeCoeffBTable) 709 if(theRangeCoeffBTable)
710 { theRangeCoeffBTable->clearAndDestroy(); 710 { theRangeCoeffBTable->clearAndDestroy();
711 delete theRangeCoeffBTable; } 711 delete theRangeCoeffBTable; }
712 theRangeCoeffBTable = new G4PhysicsTable(num 712 theRangeCoeffBTable = new G4PhysicsTable(numOfCouples);
713 713
714 G4double RTable = std::exp(std::log(highestK 714 G4double RTable = std::exp(std::log(highestKineticEnergy/lowestKineticEnergy)/TotBin) ;
715 G4double R2 = RTable*RTable ; 715 G4double R2 = RTable*RTable ;
716 G4double R1 = RTable+1.; 716 G4double R1 = RTable+1.;
717 G4double w = R1*(RTable-1.)*(RTable-1.); 717 G4double w = R1*(RTable-1.)*(RTable-1.);
718 G4double w1 = -R1/w , w2 = R1*(R2+1.)/w , w3 718 G4double w1 = -R1/w , w2 = R1*(R2+1.)/w , w3 = -R2*R1/w ;
719 G4double Ti , Tim , Tip , Ri , Rim , Rip , V 719 G4double Ti , Tim , Tip , Ri , Rim , Rip , Value ;
720 G4bool isOut; 720 G4bool isOut;
721 721
722 // loop for materials 722 // loop for materials
723 for (G4int J=0; J<numOfCouples; J++) 723 for (G4int J=0; J<numOfCouples; J++)
724 { 724 {
725 G4int binmax=TotBin ; 725 G4int binmax=TotBin ;
726 G4PhysicsLinearVector* aVector = 726 G4PhysicsLinearVector* aVector =
727 new G4PhysicsLinearVec 727 new G4PhysicsLinearVector(0.,binmax, TotBin);
728 Ti = lowestKineticEnergy ; 728 Ti = lowestKineticEnergy ;
729 G4PhysicsVector* rangeVector= (*theRangeTa 729 G4PhysicsVector* rangeVector= (*theRangeTable)[J];
730 730
731 for ( G4int i=0; i<TotBin; i++) 731 for ( G4int i=0; i<TotBin; i++)
732 { 732 {
733 Ri = rangeVector->GetValue(Ti,isOut) ; 733 Ri = rangeVector->GetValue(Ti,isOut) ;
734 if ( i==0 ) 734 if ( i==0 )
735 Rim = 0. ; 735 Rim = 0. ;
736 else 736 else
737 { 737 {
738 Tim = Ti/RTable ; 738 Tim = Ti/RTable ;
739 Rim = rangeVector->GetValue(Tim,isOut) 739 Rim = rangeVector->GetValue(Tim,isOut);
740 } 740 }
741 if ( i==(TotBin-1)) 741 if ( i==(TotBin-1))
742 Rip = Ri ; 742 Rip = Ri ;
743 else 743 else
744 { 744 {
745 Tip = Ti*RTable ; 745 Tip = Ti*RTable ;
746 Rip = rangeVector->GetValue(Tip,isOut) 746 Rip = rangeVector->GetValue(Tip,isOut);
747 } 747 }
748 Value = (w1*Rip + w2*Ri + w3*Rim)/Ti; 748 Value = (w1*Rip + w2*Ri + w3*Rim)/Ti;
749 749
750 aVector->PutValue(i,Value); 750 aVector->PutValue(i,Value);
751 Ti = RTable*Ti ; 751 Ti = RTable*Ti ;
752 } 752 }
753 theRangeCoeffBTable->insert(aVector); 753 theRangeCoeffBTable->insert(aVector);
754 } 754 }
755 return theRangeCoeffBTable ; 755 return theRangeCoeffBTable ;
756 } 756 }
757 757
758 // 758 //
759 759
760 G4PhysicsTable* G4RDVeLowEnergyLoss::BuildRang 760 G4PhysicsTable* G4RDVeLowEnergyLoss::BuildRangeCoeffCTable(G4PhysicsTable* theRangeTable,
761 G4PhysicsTable* theRangeCoeffCT 761 G4PhysicsTable* theRangeCoeffCTable,
762 G4double lowestKineticEnergy, 762 G4double lowestKineticEnergy,
763 G4double highestKineticEnergy, 763 G4double highestKineticEnergy, G4int TotBin)
764 // Build tables of coefficients for the energy 764 // Build tables of coefficients for the energy loss calculation
765 // create table for coefficients "C" 765 // create table for coefficients "C"
766 { 766 {
767 767
768 G4int numOfCouples = G4ProductionCutsTable:: 768 G4int numOfCouples = G4ProductionCutsTable::GetProductionCutsTable()->GetTableSize();
769 769
770 if(theRangeCoeffCTable) 770 if(theRangeCoeffCTable)
771 { theRangeCoeffCTable->clearAndDestroy(); 771 { theRangeCoeffCTable->clearAndDestroy();
772 delete theRangeCoeffCTable; } 772 delete theRangeCoeffCTable; }
773 theRangeCoeffCTable = new G4PhysicsTable(num 773 theRangeCoeffCTable = new G4PhysicsTable(numOfCouples);
774 774
775 G4double RTable = std::exp(std::log(highestK 775 G4double RTable = std::exp(std::log(highestKineticEnergy/lowestKineticEnergy)/TotBin) ;
776 G4double R2 = RTable*RTable ; 776 G4double R2 = RTable*RTable ;
777 G4double R1 = RTable+1.; 777 G4double R1 = RTable+1.;
778 G4double w = R1*(RTable-1.)*(RTable-1.); 778 G4double w = R1*(RTable-1.)*(RTable-1.);
779 G4double w1 = 1./w , w2 = -RTable*R1/w , w3 779 G4double w1 = 1./w , w2 = -RTable*R1/w , w3 = RTable*R2/w ;
780 G4double Ti , Tim , Tip , Ri , Rim , Rip , V 780 G4double Ti , Tim , Tip , Ri , Rim , Rip , Value ;
781 G4bool isOut; 781 G4bool isOut;
782 782
783 // loop for materials 783 // loop for materials
784 for (G4int J=0; J<numOfCouples; J++) 784 for (G4int J=0; J<numOfCouples; J++)
785 { 785 {
786 G4int binmax=TotBin ; 786 G4int binmax=TotBin ;
787 G4PhysicsLinearVector* aVector = 787 G4PhysicsLinearVector* aVector =
788 new G4PhysicsLinearVecto 788 new G4PhysicsLinearVector(0.,binmax, TotBin);
789 Ti = lowestKineticEnergy ; 789 Ti = lowestKineticEnergy ;
790 G4PhysicsVector* rangeVector= (*theRangeTa 790 G4PhysicsVector* rangeVector= (*theRangeTable)[J];
791 791
792 for ( G4int i=0; i<TotBin; i++) 792 for ( G4int i=0; i<TotBin; i++)
793 { 793 {
794 Ri = rangeVector->GetValue(Ti,isOut) ; 794 Ri = rangeVector->GetValue(Ti,isOut) ;
795 if ( i==0 ) 795 if ( i==0 )
796 Rim = 0. ; 796 Rim = 0. ;
797 else 797 else
798 { 798 {
799 Tim = Ti/RTable ; 799 Tim = Ti/RTable ;
800 Rim = rangeVector->GetValue(Tim,isOut) 800 Rim = rangeVector->GetValue(Tim,isOut);
801 } 801 }
802 if ( i==(TotBin-1)) 802 if ( i==(TotBin-1))
803 Rip = Ri ; 803 Rip = Ri ;
804 else 804 else
805 { 805 {
806 Tip = Ti*RTable ; 806 Tip = Ti*RTable ;
807 Rip = rangeVector->GetValue(Tip,isOut) 807 Rip = rangeVector->GetValue(Tip,isOut);
808 } 808 }
809 Value = w1*Rip + w2*Ri + w3*Rim ; 809 Value = w1*Rip + w2*Ri + w3*Rim ;
810 810
811 aVector->PutValue(i,Value); 811 aVector->PutValue(i,Value);
812 Ti = RTable*Ti ; 812 Ti = RTable*Ti ;
813 } 813 }
814 theRangeCoeffCTable->insert(aVector); 814 theRangeCoeffCTable->insert(aVector);
815 } 815 }
816 return theRangeCoeffCTable ; 816 return theRangeCoeffCTable ;
817 } 817 }
818 818
819 // 819 //
820 820
821 G4double G4RDVeLowEnergyLoss::GetLossWithFluct 821 G4double G4RDVeLowEnergyLoss::GetLossWithFluct(const G4DynamicParticle* aParticle,
822 c 822 const G4MaterialCutsCouple* couple,
823 G4double MeanLoss, 823 G4double MeanLoss,
824 G4double step) 824 G4double step)
825 // calculate actual loss from the mean loss 825 // calculate actual loss from the mean loss
826 // The model used to get the fluctuation is e 826 // The model used to get the fluctuation is essentially the same as in Glandz in Geant3.
827 { 827 {
828 static const G4double minLoss = 1.*eV ; 828 static const G4double minLoss = 1.*eV ;
829 static const G4double probLim = 0.01 ; 829 static const G4double probLim = 0.01 ;
830 static const G4double sumaLim = -std::log(p 830 static const G4double sumaLim = -std::log(probLim) ;
831 static const G4double alim=10.; 831 static const G4double alim=10.;
832 static const G4double kappa = 10. ; 832 static const G4double kappa = 10. ;
833 static const G4double factor = twopi_mc2_rc 833 static const G4double factor = twopi_mc2_rcl2 ;
834 const G4Material* aMaterial = couple->GetMat 834 const G4Material* aMaterial = couple->GetMaterial();
835 835
836 // check if the material has changed ( cache 836 // check if the material has changed ( cache mechanism)
837 837
838 if (aMaterial != lastMaterial) 838 if (aMaterial != lastMaterial)
839 { 839 {
840 lastMaterial = aMaterial; 840 lastMaterial = aMaterial;
841 imat = couple->GetIndex(); 841 imat = couple->GetIndex();
842 f1Fluct = aMaterial->GetIonisation( 842 f1Fluct = aMaterial->GetIonisation()->GetF1fluct();
843 f2Fluct = aMaterial->GetIonisation( 843 f2Fluct = aMaterial->GetIonisation()->GetF2fluct();
844 e1Fluct = aMaterial->GetIonisation( 844 e1Fluct = aMaterial->GetIonisation()->GetEnergy1fluct();
845 e2Fluct = aMaterial->GetIonisation( 845 e2Fluct = aMaterial->GetIonisation()->GetEnergy2fluct();
846 e1LogFluct = aMaterial->GetIonisation( 846 e1LogFluct = aMaterial->GetIonisation()->GetLogEnergy1fluct();
847 e2LogFluct = aMaterial->GetIonisation( 847 e2LogFluct = aMaterial->GetIonisation()->GetLogEnergy2fluct();
848 rateFluct = aMaterial->GetIonisation( 848 rateFluct = aMaterial->GetIonisation()->GetRateionexcfluct();
849 ipotFluct = aMaterial->GetIonisation( 849 ipotFluct = aMaterial->GetIonisation()->GetMeanExcitationEnergy();
850 ipotLogFluct = aMaterial->GetIonisation( 850 ipotLogFluct = aMaterial->GetIonisation()->GetLogMeanExcEnergy();
851 } 851 }
852 G4double threshold,w1,w2,C, 852 G4double threshold,w1,w2,C,
853 beta2,suma,e0,loss,lossc,w; 853 beta2,suma,e0,loss,lossc,w;
854 G4double a1,a2,a3; 854 G4double a1,a2,a3;
855 G4int p1,p2,p3; 855 G4int p1,p2,p3;
856 G4int nb; 856 G4int nb;
857 G4double Corrfac, na,alfa,rfac,namean,sa,alf 857 G4double Corrfac, na,alfa,rfac,namean,sa,alfa1,ea,sea;
858 // G4double dp1; 858 // G4double dp1;
859 G4double dp3; 859 G4double dp3;
860 G4double siga ; 860 G4double siga ;
861 861
862 // shortcut for very very small loss 862 // shortcut for very very small loss
863 if(MeanLoss < minLoss) return MeanLoss ; 863 if(MeanLoss < minLoss) return MeanLoss ;
864 864
865 // get particle data 865 // get particle data
866 G4double Tkin = aParticle->GetKineticEnerg 866 G4double Tkin = aParticle->GetKineticEnergy();
867 867
868 // G4cout << "MGP -- Fluc Tkin " << Tkin/ke 868 // G4cout << "MGP -- Fluc Tkin " << Tkin/keV << " keV " << " MeanLoss = " << MeanLoss/keV << G4endl;
869 869
870 threshold = (*((G4ProductionCutsTable::GetPr 870 threshold = (*((G4ProductionCutsTable::GetProductionCutsTable())
871 ->GetEnergyCutsVector(1)))[ima 871 ->GetEnergyCutsVector(1)))[imat];
872 G4double rmass = electron_mass_c2/ParticleMa 872 G4double rmass = electron_mass_c2/ParticleMass;
873 G4double tau = Tkin/ParticleMass, tau1 = t 873 G4double tau = Tkin/ParticleMass, tau1 = tau+1., tau2 = tau*(tau+2.);
874 G4double Tm = 2.*electron_mass_c2*tau2/(1 874 G4double Tm = 2.*electron_mass_c2*tau2/(1.+2.*tau1*rmass+rmass*rmass);
875 875
876 // G4cout << "MGP Particle mass " << Particl 876 // G4cout << "MGP Particle mass " << ParticleMass/MeV << " Tm " << Tm << G4endl;
877 877
878 if(Tm > threshold) Tm = threshold; 878 if(Tm > threshold) Tm = threshold;
879 beta2 = tau2/(tau1*tau1); 879 beta2 = tau2/(tau1*tau1);
880 880
881 // Gaussian fluctuation ? 881 // Gaussian fluctuation ?
882 if(MeanLoss >= kappa*Tm || MeanLoss <= kappa 882 if(MeanLoss >= kappa*Tm || MeanLoss <= kappa*ipotFluct)
883 { 883 {
884 G4double electronDensity = aMaterial->GetE 884 G4double electronDensity = aMaterial->GetElectronDensity() ;
885 siga = std::sqrt(Tm*(1.0-0.5*beta2)*step* 885 siga = std::sqrt(Tm*(1.0-0.5*beta2)*step*
886 factor*electronDensity/beta2) 886 factor*electronDensity/beta2) ;
887 do { 887 do {
888 loss = G4RandGauss::shoot(MeanLoss,siga) 888 loss = G4RandGauss::shoot(MeanLoss,siga) ;
889 } while (loss < 0. || loss > 2.0*MeanLoss) 889 } while (loss < 0. || loss > 2.0*MeanLoss);
890 return loss ; 890 return loss ;
891 } 891 }
892 892
893 w1 = Tm/ipotFluct; 893 w1 = Tm/ipotFluct;
894 w2 = std::log(2.*electron_mass_c2*tau2); 894 w2 = std::log(2.*electron_mass_c2*tau2);
895 895
896 C = MeanLoss*(1.-rateFluct)/(w2-ipotLogFluct 896 C = MeanLoss*(1.-rateFluct)/(w2-ipotLogFluct-beta2);
897 897
898 a1 = C*f1Fluct*(w2-e1LogFluct-beta2)/e1Fluct 898 a1 = C*f1Fluct*(w2-e1LogFluct-beta2)/e1Fluct;
899 a2 = C*f2Fluct*(w2-e2LogFluct-beta2)/e2Fluct 899 a2 = C*f2Fluct*(w2-e2LogFluct-beta2)/e2Fluct;
900 a3 = rateFluct*MeanLoss*(Tm-ipotFluct)/(ipot 900 a3 = rateFluct*MeanLoss*(Tm-ipotFluct)/(ipotFluct*Tm*std::log(w1));
901 901
902 suma = a1+a2+a3; 902 suma = a1+a2+a3;
903 903
904 loss = 0. ; 904 loss = 0. ;
905 905
906 if(suma < sumaLim) // very small 906 if(suma < sumaLim) // very small Step
907 { 907 {
908 e0 = aMaterial->GetIonisation()->GetEner 908 e0 = aMaterial->GetIonisation()->GetEnergy0fluct();
909 // G4cout << "MGP e0 = " << e0/keV << G4 909 // G4cout << "MGP e0 = " << e0/keV << G4endl;
910 910
911 if(Tm == ipotFluct) 911 if(Tm == ipotFluct)
912 { 912 {
913 a3 = MeanLoss/e0; 913 a3 = MeanLoss/e0;
914 914
915 if(a3>alim) 915 if(a3>alim)
916 { 916 {
917 siga=std::sqrt(a3) ; 917 siga=std::sqrt(a3) ;
918 p3 = std::max(0,G4int(G4RandGauss::sho 918 p3 = std::max(0,G4int(G4RandGauss::shoot(a3,siga)+0.5));
919 } 919 }
920 else p3 = G4Poisson(a3); 920 else p3 = G4Poisson(a3);
921 921
922 loss = p3*e0 ; 922 loss = p3*e0 ;
923 923
924 if(p3 > 0) loss += (1.-2.*G4UniformRand()) 924 if(p3 > 0) loss += (1.-2.*G4UniformRand())*e0 ;
925 // G4cout << "MGP very small step " << los 925 // G4cout << "MGP very small step " << loss/keV << G4endl;
926 } 926 }
927 else 927 else
928 { 928 {
929 // G4cout << "MGP old Tm = " << Tm << " 929 // G4cout << "MGP old Tm = " << Tm << " " << ipotFluct << " " << e0 << G4endl;
930 Tm = Tm-ipotFluct+e0 ; 930 Tm = Tm-ipotFluct+e0 ;
931 931
932 // MGP ---- workaround to avoid log argume 932 // MGP ---- workaround to avoid log argument<0, TO BE CHECKED
933 if (Tm <= 0.) 933 if (Tm <= 0.)
934 { 934 {
935 loss = MeanLoss; 935 loss = MeanLoss;
936 p3 = 0; 936 p3 = 0;
937 // G4cout << "MGP correction loss = Me 937 // G4cout << "MGP correction loss = MeanLoss " << loss/keV << G4endl;
938 } 938 }
939 else 939 else
940 { 940 {
941 a3 = MeanLoss*(Tm-e0)/(Tm*e0*std::log( 941 a3 = MeanLoss*(Tm-e0)/(Tm*e0*std::log(Tm/e0));
942 942
943 // G4cout << "MGP new Tm = " << Tm << 943 // G4cout << "MGP new Tm = " << Tm << " " << ipotFluct << " " << e0 << " a3= " << a3 << G4endl;
944 944
945 if(a3>alim) 945 if(a3>alim)
946 { 946 {
947 siga=std::sqrt(a3) ; 947 siga=std::sqrt(a3) ;
948 p3 = std::max(0,G4int(G4RandGauss::shoot 948 p3 = std::max(0,G4int(G4RandGauss::shoot(a3,siga)+0.5));
949 } 949 }
950 else 950 else
951 p3 = G4Poisson(a3); 951 p3 = G4Poisson(a3);
952 //G4cout << "MGP p3 " << p3 << G4endl; 952 //G4cout << "MGP p3 " << p3 << G4endl;
953 953
954 } 954 }
955 955
956 if(p3 > 0) 956 if(p3 > 0)
957 { 957 {
958 w = (Tm-e0)/Tm ; 958 w = (Tm-e0)/Tm ;
959 if(p3 > nmaxCont2) 959 if(p3 > nmaxCont2)
960 { 960 {
961 // G4cout << "MGP dp3 " << dp3 << " p3 " 961 // G4cout << "MGP dp3 " << dp3 << " p3 " << p3 << " " << nmaxCont2 << G4endl;
962 dp3 = G4double(p3) ; 962 dp3 = G4double(p3) ;
963 Corrfac = dp3/G4double(nmaxCont2) ; 963 Corrfac = dp3/G4double(nmaxCont2) ;
964 p3 = nmaxCont2 ; 964 p3 = nmaxCont2 ;
965 } 965 }
966 else 966 else
967 Corrfac = 1. ; 967 Corrfac = 1. ;
968 968
969 for(G4int i=0; i<p3; i++) loss += 1./( 969 for(G4int i=0; i<p3; i++) loss += 1./(1.-w*G4UniformRand()) ;
970 loss *= e0*Corrfac ; 970 loss *= e0*Corrfac ;
971 // G4cout << "MGP Corrfac = " << Corrf 971 // G4cout << "MGP Corrfac = " << Corrfac << " e0 = " << e0/keV << " loss = " << loss/keV << G4endl;
972 } 972 }
973 } 973 }
974 } 974 }
975 975
976 else // not so 976 else // not so small Step
977 { 977 {
978 // excitation type 1 978 // excitation type 1
979 if(a1>alim) 979 if(a1>alim)
980 { 980 {
981 siga=std::sqrt(a1) ; 981 siga=std::sqrt(a1) ;
982 p1 = std::max(0,int(G4RandGauss::shoot 982 p1 = std::max(0,int(G4RandGauss::shoot(a1,siga)+0.5));
983 } 983 }
984 else 984 else
985 p1 = G4Poisson(a1); 985 p1 = G4Poisson(a1);
986 986
987 // excitation type 2 987 // excitation type 2
988 if(a2>alim) 988 if(a2>alim)
989 { 989 {
990 siga=std::sqrt(a2) ; 990 siga=std::sqrt(a2) ;
991 p2 = std::max(0,int(G4RandGauss::shoot 991 p2 = std::max(0,int(G4RandGauss::shoot(a2,siga)+0.5));
992 } 992 }
993 else 993 else
994 p2 = G4Poisson(a2); 994 p2 = G4Poisson(a2);
995 995
996 loss = p1*e1Fluct+p2*e2Fluct; 996 loss = p1*e1Fluct+p2*e2Fluct;
997 997
998 // smearing to avoid unphysical peaks 998 // smearing to avoid unphysical peaks
999 if(p2 > 0) 999 if(p2 > 0)
1000 loss += (1.-2.*G4UniformRand())*e2Flu 1000 loss += (1.-2.*G4UniformRand())*e2Fluct;
1001 else if (loss>0.) 1001 else if (loss>0.)
1002 loss += (1.-2.*G4UniformRand())*e1Flu 1002 loss += (1.-2.*G4UniformRand())*e1Fluct;
1003 1003
1004 // ionisation ......................... 1004 // ionisation .......................................
1005 if(a3 > 0.) 1005 if(a3 > 0.)
1006 { 1006 {
1007 if(a3>alim) 1007 if(a3>alim)
1008 { 1008 {
1009 siga=std::sqrt(a3) ; 1009 siga=std::sqrt(a3) ;
1010 p3 = std::max(0,int(G4RandGauss::shoo 1010 p3 = std::max(0,int(G4RandGauss::shoot(a3,siga)+0.5));
1011 } 1011 }
1012 else 1012 else
1013 p3 = G4Poisson(a3); 1013 p3 = G4Poisson(a3);
1014 1014
1015 lossc = 0.; 1015 lossc = 0.;
1016 if(p3 > 0) 1016 if(p3 > 0)
1017 { 1017 {
1018 na = 0.; 1018 na = 0.;
1019 alfa = 1.; 1019 alfa = 1.;
1020 if (p3 > nmaxCont2) 1020 if (p3 > nmaxCont2)
1021 { 1021 {
1022 dp3 = G4double(p3); 1022 dp3 = G4double(p3);
1023 rfac = dp3/(G4double(nmaxCont2)+d 1023 rfac = dp3/(G4double(nmaxCont2)+dp3);
1024 namean = G4double(p3)*rfac; 1024 namean = G4double(p3)*rfac;
1025 sa = G4double(nmaxCont1)*rfac; 1025 sa = G4double(nmaxCont1)*rfac;
1026 na = G4RandGauss::shoot(namean, 1026 na = G4RandGauss::shoot(namean,sa);
1027 if (na > 0.) 1027 if (na > 0.)
1028 { 1028 {
1029 alfa = w1*G4double(nmaxCont2+p3)/ 1029 alfa = w1*G4double(nmaxCont2+p3)/(w1*G4double(nmaxCont2)+G4double(p3));
1030 alfa1 = alfa*std::log(alfa)/(alfa- 1030 alfa1 = alfa*std::log(alfa)/(alfa-1.);
1031 ea = na*ipotFluct*alfa1; 1031 ea = na*ipotFluct*alfa1;
1032 sea = ipotFluct*std::sqrt(na*(al 1032 sea = ipotFluct*std::sqrt(na*(alfa-alfa1*alfa1));
1033 lossc += G4RandGauss::shoot(ea,sea) 1033 lossc += G4RandGauss::shoot(ea,sea);
1034 } 1034 }
1035 } 1035 }
1036 1036
1037 nb = G4int(G4double(p3)-na); 1037 nb = G4int(G4double(p3)-na);
1038 if (nb > 0) 1038 if (nb > 0)
1039 { 1039 {
1040 w2 = alfa*ipotFluct; 1040 w2 = alfa*ipotFluct;
1041 w = (Tm-w2)/Tm; 1041 w = (Tm-w2)/Tm;
1042 for (G4int k=0; k<nb; k++) lossc += w2/ 1042 for (G4int k=0; k<nb; k++) lossc += w2/(1.-w*G4UniformRand());
1043 } 1043 }
1044 } 1044 }
1045 1045
1046 loss += lossc; 1046 loss += lossc;
1047 } 1047 }
1048 } 1048 }
1049 1049
1050 return loss ; 1050 return loss ;
1051 } 1051 }
1052 1052
1053 // 1053 //
1054 1054
1055 1055