source: trunk/source/processes/electromagnetic/utils/src/G4VEnergyLoss.cc@ 1193

Last change on this file since 1193 was 1007, checked in by garnier, 17 years ago

update to geant4.9.2
File size: 35.3 KB
Line 
1//
2// ********************************************************************
3// * License and Disclaimer *
4// * *
5// * The Geant4 software is copyright of the Copyright Holders of *
6// * the Geant4 Collaboration. It is provided under the terms and *
7// * conditions of the Geant4 Software License, included in the file *
8// * LICENSE and available at http://cern.ch/geant4/license . These *
9// * include a list of copyright holders. *
10// * *
11// * Neither the authors of this software system, nor their employing *
12// * institutes,nor the agencies providing financial support for this *
13// * work make any representation or warranty, express or implied, *
14// * regarding this software system or assume any liability for its *
15// * use. Please see the license in the file LICENSE and URL above *
16// * for the full disclaimer and the limitation of liability. *
17// * *
18// * This code implementation is the result of the scientific and *
19// * technical work of the GEANT4 collaboration. *
20// * By using, copying, modifying or distributing the software (or *
21// * any work based on the software) you agree to acknowledge its *
22// * use in resulting scientific publications, and indicate your *
23// * acceptance of all terms of the Geant4 Software license. *
24// ********************************************************************
25//
26//
27// $Id: G4VEnergyLoss.cc,v 1.46 2006/06/29 19:55:21 gunter Exp $
28// GEANT4 tag $Name: geant4-09-02 $
29//
30
31// --------------------------------------------------------------
32//
33// bug fixed in fluct., L.Urban 01/02/01
34// bug fixed in fluct., L.Urban 26/05/00
35// bug fixed in fluct., L.Urban 22/11/00
36// bugfix in fluct.
37// (some variables are doubles instead of ints now),L.Urban 23/03/01
38// 18/05/01 V.Ivanchenko Clean up againist Linux ANSI compilation
39// 17-09-01 migration of Materials to pure STL (mma)
40// 26-10-01 static inline functions moved from .hh file (mma)
41// 08.11.01 some static methods,data members are not static L.Urban
42// 11.02.02 subSecFlag = false --> No sucutoff generation (mma)
43// 14.02.02 initial value of data member finalRange has been changed L.Urban
44// 26.02.02 initial value of data member finalRange = 1 mm (mma)
45// 21.07.02 V.Ivanchenko Fix at low energies - if tmax below ionisation
46// potential then only Gaussian fluctuations are sampled.
47// 15.01.03 Migrade to cut per region (V.Ivanchenko)
48// 05.02.03 Minor fix for several region case (V.Ivanchenko)
49// 25.03.03 add finalRangeRequested (mma)
50//
51//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
52//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
53
54#include "G4VEnergyLoss.hh"
55#include "G4EnergyLossMessenger.hh"
56#include "G4ProductionCutsTable.hh"
57#include "G4LossTableManager.hh"
58
59//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
60
61G4bool G4VEnergyLoss::rndmStepFlag = false;
62G4bool G4VEnergyLoss::EnlossFlucFlag = true;
63
64G4bool G4VEnergyLoss::subSecFlag = false;
65G4double G4VEnergyLoss::MinDeltaCutInRange = 0.100*mm;
66G4double* G4VEnergyLoss::MinDeltaEnergy = 0;
67G4bool* G4VEnergyLoss::LowerLimitForced = 0;
68G4bool G4VEnergyLoss::setMinDeltaCutInRange = false;
69
70G4double G4VEnergyLoss::dRoverRange = 20*perCent;
71G4double G4VEnergyLoss::finalRange = 1*mm;
72G4double G4VEnergyLoss::finalRangeRequested = -1*mm;
73G4double G4VEnergyLoss::c1lim = dRoverRange;
74G4double G4VEnergyLoss::c2lim = 2.*(1.-dRoverRange)*finalRange;
75G4double G4VEnergyLoss::c3lim = -(1.-dRoverRange)*finalRange*finalRange;
76
77G4EnergyLossMessenger* G4VEnergyLoss::ELossMessenger = 0;
78
79//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
80
81G4VEnergyLoss::G4VEnergyLoss()
82 :G4VContinuousDiscreteProcess("No Name Loss Process"),
83 lastMaterial(NULL),
84 nmaxCont1(4),
85 nmaxCont2(16)
86{
87 G4Exception("G4VEnergyLoss:: default constructor is called");
88}
89
90//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
91
92G4VEnergyLoss::G4VEnergyLoss(const G4String& aName , G4ProcessType aType)
93 : G4VContinuousDiscreteProcess(aName, aType),
94 lastMaterial(NULL),
95 nmaxCont1(4),
96 nmaxCont2(16)
97{
98 //create (only once) EnergyLoss messenger
99// if(!ELossMessenger) ELossMessenger = new G4EnergyLossMessenger();
100 if(!ELossMessenger) ELossMessenger = G4LossTableManager::Instance()->GetMessenger();
101}
102
103//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
104
105G4VEnergyLoss::~G4VEnergyLoss()
106{}
107
108//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
109
110G4VEnergyLoss::G4VEnergyLoss(G4VEnergyLoss& right)
111 : G4VContinuousDiscreteProcess(right),
112 lastMaterial(NULL),
113 nmaxCont1(4),
114 nmaxCont2(16)
115{}
116
117//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
118
119void G4VEnergyLoss::SetRndmStep(G4bool value) {rndmStepFlag = value;}
120
121//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
122
123void G4VEnergyLoss::SetEnlossFluc(G4bool value) {EnlossFlucFlag = value;}
124
125//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
126
127void G4VEnergyLoss::SetSubSec(G4bool value) {subSecFlag = value;}
128
129//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
130
131void G4VEnergyLoss::SetMinDeltaCutInRange(G4double value)
132{MinDeltaCutInRange = value; setMinDeltaCutInRange = true;}
133
134//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
135
136void G4VEnergyLoss::SetStepFunction(G4double c1, G4double c2)
137{
138 dRoverRange = c1; finalRangeRequested = c2;
139 c1lim=dRoverRange;
140 c2lim=2.*(1-dRoverRange)*finalRangeRequested;
141 c3lim=-(1.-dRoverRange)*finalRangeRequested*finalRangeRequested;
142}
143
144//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
145
146G4PhysicsTable* G4VEnergyLoss::BuildRangeTable(
147 G4PhysicsTable* theDEDXTable,G4PhysicsTable* theRangeTable,
148 G4double LowestKineticEnergy,G4double HighestKineticEnergy,G4int TotBin)
149// Build range table from the energy loss table
150{
151 size_t numOfMaterials = theDEDXTable->length();
152
153 if(theRangeTable)
154 { theRangeTable->clearAndDestroy();
155 delete theRangeTable; }
156 theRangeTable = new G4PhysicsTable(numOfMaterials);
157
158 // loop for materials
159
160 for (size_t J=0; J<numOfMaterials; J++)
161 {
162 G4PhysicsLogVector* aVector;
163 aVector = new G4PhysicsLogVector(LowestKineticEnergy,
164 HighestKineticEnergy,TotBin);
165
166 BuildRangeVector(theDEDXTable,LowestKineticEnergy,HighestKineticEnergy,
167 TotBin,J,aVector);
168 theRangeTable->insert(aVector);
169 }
170 return theRangeTable ;
171}
172
173//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
174
175void G4VEnergyLoss::BuildRangeVector(G4PhysicsTable* theDEDXTable,
176 G4double,G4double HighestKineticEnergy,G4int TotBin,
177 G4int materialIndex,G4PhysicsLogVector* rangeVector)
178// create range vector for a material
179{
180 G4int nbin=100,i;
181 G4bool isOut;
182
183 const G4double small = 1.e-6 ;
184 const G4double masslimit = 0.52*MeV ;
185
186 G4double tlim=2.*MeV,t1=0.1*MeV,t2=0.025*MeV ;
187 G4double tlime=0.2*keV,factor=2.*electron_mass_c2 ;
188 G4double loss1,loss2,ca,cb,cba ;
189 G4double losslim,clim ;
190 G4double taulim,rangelim,ltaulim,ltaumax,
191 LowEdgeEnergy,tau,Value,tau1,sqtau1 ;
192 G4double oldValue,tauold ;
193
194 G4PhysicsVector* physicsVector= (*theDEDXTable)[materialIndex];
195
196 // cure 'accidental' 0. dE/dx vales first .....
197 G4double lossmin = +1.e10 ;
198 for (G4int i1=0; i1<TotBin; i1++)
199 {
200 LowEdgeEnergy = rangeVector->GetLowEdgeEnergy(i1) ;
201 Value = physicsVector->GetValue(LowEdgeEnergy,isOut);
202 if((Value < lossmin)&&(Value > 0.))
203 lossmin = Value ;
204 }
205 for (G4int i2=0; i2<TotBin; i2++)
206 {
207 LowEdgeEnergy = rangeVector->GetLowEdgeEnergy(i2) ;
208 Value = physicsVector->GetValue(LowEdgeEnergy,isOut);
209 if(Value < lossmin)
210 physicsVector->PutValue(i2,small*lossmin) ;
211 }
212
213 // low energy part first...
214 // heavy particle
215 if(ParticleMass > masslimit)
216 {
217 loss1 = physicsVector->GetValue(t1,isOut);
218 loss2 = physicsVector->GetValue(t2,isOut);
219 tau1 = t1/ParticleMass ;
220 sqtau1 = std::sqrt(tau1) ;
221 ca = (4.*loss2-loss1)/sqtau1 ;
222 cb = (2.*loss1-4.*loss2)/tau1 ;
223 cba = cb/ca ;
224 taulim = tlim/ParticleMass ;
225 ltaulim = std::log(taulim) ;
226 ltaumax = std::log(HighestKineticEnergy/ParticleMass) ;
227
228 i=-1;
229 oldValue = 0. ;
230
231 do
232 {
233 i += 1 ;
234 LowEdgeEnergy = rangeVector->GetLowEdgeEnergy(i);
235 tau = LowEdgeEnergy/ParticleMass;
236 if ( tau <= tau1 )
237 {
238 Value = 2.*ParticleMass*std::log(1.+cba*std::sqrt(tau))/cb ;
239 }
240 else
241 {
242 Value = 2.*ParticleMass*std::log(1.+cba*sqtau1)/cb ;
243 if(tau<=taulim)
244 {
245 taulow = tau1 ;
246 tauhigh = tau ;
247 Value += RangeIntLin(physicsVector,nbin);
248 }
249 else
250 {
251
252 taulow = tau1 ;
253 tauhigh = taulim ;
254 Value += RangeIntLin(physicsVector,nbin) ;
255 ltaulow = ltaulim ;
256 ltauhigh = std::log(tau) ;
257 Value += RangeIntLog(physicsVector,nbin);
258 }
259 }
260
261 rangeVector->PutValue(i,Value);
262 oldValue = Value ;
263 tauold = tau ;
264 } while (tau<=taulim) ;
265 }
266 else
267 // electron/positron
268 {
269 losslim = physicsVector->GetValue(tlime,isOut) ;
270
271 taulim = tlime/electron_mass_c2;
272 clim = losslim;
273 ltaulim = std::log(taulim);
274 ltaumax = std::log(HighestKineticEnergy/electron_mass_c2);
275
276 i=-1;
277 oldValue = 0.;
278
279 do
280 {
281 i += 1 ;
282 LowEdgeEnergy = rangeVector->GetLowEdgeEnergy(i);
283 tau = LowEdgeEnergy/electron_mass_c2;
284 if (tau <= taulim) Value = factor*std::sqrt(tau*taulim)/clim;
285 else {
286 rangelim = factor*taulim/losslim ;
287 ltaulow = std::log(taulim);
288 ltauhigh = std::log(tau);
289 Value = rangelim+RangeIntLog(physicsVector,nbin);
290 }
291 rangeVector->PutValue(i,Value);
292 oldValue = Value;
293 tauold = tau;
294
295 } while (tau<=taulim);
296
297 }
298
299 i += 1 ;
300 for (G4int j=i; j<TotBin; j++)
301 {
302 LowEdgeEnergy = rangeVector->GetLowEdgeEnergy(j);
303 tau = LowEdgeEnergy/ParticleMass;
304 ltaulow = std::log(tauold);
305 ltauhigh = std::log(tau);
306 Value = oldValue+RangeIntLog(physicsVector,nbin);
307 rangeVector->PutValue(j,Value);
308 oldValue = Value ;
309
310 tauold = tau ;
311 }
312}
313
314//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
315
316G4double G4VEnergyLoss::RangeIntLin(G4PhysicsVector* physicsVector,
317 G4int nbin)
318// num. integration, linear binning
319{
320 G4double dtau,Value,taui,ti,lossi,ci;
321 G4bool isOut;
322 dtau = (tauhigh-taulow)/nbin;
323 Value = 0.;
324
325 for (G4int i=0; i<=nbin; i++)
326 {
327 taui = taulow + dtau*i ;
328 ti = ParticleMass*taui;
329 lossi = physicsVector->GetValue(ti,isOut);
330 if(i==0)
331 ci=0.5;
332 else
333 {
334 if(i<nbin)
335 ci=1.;
336 else
337 ci=0.5;
338 }
339 Value += ci/lossi;
340 }
341 Value *= ParticleMass*dtau;
342 return Value;
343}
344
345//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
346
347G4double G4VEnergyLoss::RangeIntLog(G4PhysicsVector* physicsVector,
348 G4int nbin)
349// num. integration, logarithmic binning
350{
351 G4double ltt,dltau,Value,ui,taui,ti,lossi,ci;
352 G4bool isOut;
353 ltt = ltauhigh-ltaulow;
354 dltau = ltt/nbin;
355 Value = 0.;
356
357 for (G4int i=0; i<=nbin; i++)
358 {
359 ui = ltaulow+dltau*i;
360 taui = std::exp(ui);
361 ti = ParticleMass*taui;
362 lossi = physicsVector->GetValue(ti,isOut);
363 if(i==0)
364 ci=0.5;
365 else
366 {
367 if(i<nbin)
368 ci=1.;
369 else
370 ci=0.5;
371 }
372 Value += ci*taui/lossi;
373 }
374 Value *= ParticleMass*dltau;
375 return Value;
376}
377
378//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
379
380G4PhysicsTable* G4VEnergyLoss::BuildLabTimeTable(G4PhysicsTable* theDEDXTable,
381 G4PhysicsTable* theLabTimeTable,
382 G4double LowestKineticEnergy,
383 G4double HighestKineticEnergy,G4int TotBin)
384
385{
386 size_t numOfMaterials = theDEDXTable->length();
387
388 if(theLabTimeTable)
389 { theLabTimeTable->clearAndDestroy();
390 delete theLabTimeTable; }
391 theLabTimeTable = new G4PhysicsTable(numOfMaterials);
392
393
394 for (size_t J=0; J<numOfMaterials; J++)
395 {
396 G4PhysicsLogVector* aVector;
397
398 aVector = new G4PhysicsLogVector(LowestKineticEnergy,
399 HighestKineticEnergy,TotBin);
400
401 BuildLabTimeVector(theDEDXTable,
402 LowestKineticEnergy,HighestKineticEnergy,TotBin,J,aVector);
403 theLabTimeTable->insert(aVector);
404
405
406 }
407 return theLabTimeTable ;
408}
409
410//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
411
412G4PhysicsTable* G4VEnergyLoss::BuildProperTimeTable(G4PhysicsTable* theDEDXTable,
413 G4PhysicsTable* theProperTimeTable,
414 G4double LowestKineticEnergy,
415 G4double HighestKineticEnergy,G4int TotBin)
416
417{
418 size_t numOfMaterials = theDEDXTable->length();
419
420 if(theProperTimeTable)
421 { theProperTimeTable->clearAndDestroy();
422 delete theProperTimeTable; }
423 theProperTimeTable = new G4PhysicsTable(numOfMaterials);
424
425
426 for (size_t J=0; J<numOfMaterials; J++)
427 {
428 G4PhysicsLogVector* aVector;
429
430 aVector = new G4PhysicsLogVector(LowestKineticEnergy,
431 HighestKineticEnergy,TotBin);
432
433 BuildProperTimeVector(theDEDXTable,
434 LowestKineticEnergy,HighestKineticEnergy,TotBin,J,aVector);
435 theProperTimeTable->insert(aVector);
436
437
438 }
439 return theProperTimeTable ;
440}
441
442//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
443
444void G4VEnergyLoss::BuildLabTimeVector(G4PhysicsTable* theDEDXTable,
445 G4double,
446 G4double HighestKineticEnergy,G4int TotBin,
447 G4int materialIndex, G4PhysicsLogVector* timeVector)
448// create lab time vector for a material
449{
450
451 G4int nbin=100;
452 G4bool isOut;
453 G4double tlim=5.*keV,parlowen=0.4,ppar=0.5-parlowen ;
454 G4double losslim,clim,taulim,timelim,ltaulim,ltaumax,
455 LowEdgeEnergy,tau,Value ;
456
457 G4PhysicsVector* physicsVector= (*theDEDXTable)[materialIndex];
458
459 // low energy part first...
460 losslim = physicsVector->GetValue(tlim,isOut);
461 taulim=tlim/ParticleMass ;
462 clim=std::sqrt(ParticleMass*tlim/2.)/(c_light*losslim*ppar) ;
463 ltaulim = std::log(taulim);
464 ltaumax = std::log(HighestKineticEnergy/ParticleMass) ;
465
466 G4int i=-1;
467 G4double oldValue = 0. ;
468 G4double tauold ;
469 do
470 {
471 i += 1 ;
472 LowEdgeEnergy = timeVector->GetLowEdgeEnergy(i);
473 tau = LowEdgeEnergy/ParticleMass ;
474 if ( tau <= taulim )
475 {
476 Value = clim*std::exp(ppar*std::log(tau/taulim)) ;
477 }
478 else
479 {
480 timelim=clim ;
481 ltaulow = std::log(taulim);
482 ltauhigh = std::log(tau);
483 Value = timelim+LabTimeIntLog(physicsVector,nbin);
484 }
485 timeVector->PutValue(i,Value);
486 oldValue = Value ;
487 tauold = tau ;
488 } while (tau<=taulim) ;
489 i += 1 ;
490 for (G4int j=i; j<TotBin; j++)
491 {
492 LowEdgeEnergy = timeVector->GetLowEdgeEnergy(j);
493 tau = LowEdgeEnergy/ParticleMass ;
494 ltaulow = std::log(tauold);
495 ltauhigh = std::log(tau);
496 Value = oldValue+LabTimeIntLog(physicsVector,nbin);
497 timeVector->PutValue(j,Value);
498 oldValue = Value ;
499 tauold = tau ;
500 }
501}
502
503//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
504
505void G4VEnergyLoss::BuildProperTimeVector(G4PhysicsTable* theDEDXTable,
506 G4double,
507 G4double HighestKineticEnergy,G4int TotBin,
508 G4int materialIndex, G4PhysicsLogVector* timeVector)
509// create proper time vector for a material
510{
511 G4int nbin=100;
512 G4bool isOut;
513 G4double tlim=5.*keV,parlowen=0.4,ppar=0.5-parlowen ;
514 G4double losslim,clim,taulim,timelim,ltaulim,ltaumax,
515 LowEdgeEnergy,tau,Value ;
516
517 G4PhysicsVector* physicsVector= (*theDEDXTable)[materialIndex];
518
519 // low energy part first...
520 losslim = physicsVector->GetValue(tlim,isOut);
521 taulim=tlim/ParticleMass ;
522 clim=std::sqrt(ParticleMass*tlim/2.)/(c_light*losslim*ppar) ;
523 ltaulim = std::log(taulim);
524 ltaumax = std::log(HighestKineticEnergy/ParticleMass) ;
525
526 G4int i=-1;
527 G4double oldValue = 0. ;
528 G4double tauold ;
529 do
530 {
531 i += 1 ;
532 LowEdgeEnergy = timeVector->GetLowEdgeEnergy(i);
533 tau = LowEdgeEnergy/ParticleMass ;
534 if ( tau <= taulim )
535 {
536 Value = clim*std::exp(ppar*std::log(tau/taulim)) ;
537 }
538 else
539 {
540 timelim=clim ;
541 ltaulow = std::log(taulim);
542 ltauhigh = std::log(tau);
543 Value = timelim+ProperTimeIntLog(physicsVector,nbin);
544 }
545 timeVector->PutValue(i,Value);
546 oldValue = Value ;
547 tauold = tau ;
548 } while (tau<=taulim) ;
549 i += 1 ;
550 for (G4int j=i; j<TotBin; j++)
551 {
552 LowEdgeEnergy = timeVector->GetLowEdgeEnergy(j);
553 tau = LowEdgeEnergy/ParticleMass ;
554 ltaulow = std::log(tauold);
555 ltauhigh = std::log(tau);
556 Value = oldValue+ProperTimeIntLog(physicsVector,nbin);
557 timeVector->PutValue(j,Value);
558 oldValue = Value ;
559 tauold = tau ;
560 }
561}
562
563//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
564
565G4double G4VEnergyLoss::LabTimeIntLog(G4PhysicsVector* physicsVector,
566 G4int nbin)
567// num. integration, logarithmic binning
568{
569 G4double ltt,dltau,Value,ui,taui,ti,lossi,ci;
570 G4bool isOut;
571 ltt = ltauhigh-ltaulow;
572 dltau = ltt/nbin;
573 Value = 0.;
574
575 for (G4int i=0; i<=nbin; i++)
576 {
577 ui = ltaulow+dltau*i;
578 taui = std::exp(ui);
579 ti = ParticleMass*taui;
580 lossi = physicsVector->GetValue(ti,isOut);
581 if(i==0)
582 ci=0.5;
583 else
584 {
585 if(i<nbin)
586 ci=1.;
587 else
588 ci=0.5;
589 }
590 Value += ci*taui*(ti+ParticleMass)/(std::sqrt(ti*(ti+2.*ParticleMass))*lossi);
591 }
592 Value *= ParticleMass*dltau/c_light;
593 return Value;
594}
595
596//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
597
598G4double G4VEnergyLoss::ProperTimeIntLog(G4PhysicsVector* physicsVector,
599 G4int nbin)
600// num. integration, logarithmic binning
601{
602 G4double ltt,dltau,Value,ui,taui,ti,lossi,ci;
603 G4bool isOut;
604 ltt = ltauhigh-ltaulow;
605 dltau = ltt/nbin;
606 Value = 0.;
607
608 for (G4int i=0; i<=nbin; i++)
609 {
610 ui = ltaulow+dltau*i;
611 taui = std::exp(ui);
612 ti = ParticleMass*taui;
613 lossi = physicsVector->GetValue(ti,isOut);
614 if(i==0)
615 ci=0.5;
616 else
617 {
618 if(i<nbin)
619 ci=1.;
620 else
621 ci=0.5;
622 }
623 Value += ci*taui*ParticleMass/(std::sqrt(ti*(ti+2.*ParticleMass))*lossi);
624 }
625 Value *= ParticleMass*dltau/c_light;
626 return Value;
627}
628
629//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
630
631G4PhysicsTable* G4VEnergyLoss::BuildInverseRangeTable(G4PhysicsTable* theRangeTable,
632 G4PhysicsTable* theRangeCoeffATable,
633 G4PhysicsTable* theRangeCoeffBTable,
634 G4PhysicsTable* theRangeCoeffCTable,
635 G4PhysicsTable* theInverseRangeTable,
636 G4double LowestKineticEnergy,
637 G4double HighestKineticEnergy,G4int TotBin)
638// Build inverse table of the range table
639{
640 G4double SmallestRange,BiggestRange ;
641 G4bool isOut ;
642 size_t numOfMaterials = theRangeTable->length();
643
644 if(theInverseRangeTable)
645 { theInverseRangeTable->clearAndDestroy();
646 delete theInverseRangeTable; }
647 theInverseRangeTable = new G4PhysicsTable(numOfMaterials);
648
649 // loop for materials
650 for (size_t J=0; J<numOfMaterials; J++)
651 {
652 SmallestRange = (*theRangeTable)(J)->
653 GetValue(LowestKineticEnergy,isOut) ;
654 BiggestRange = (*theRangeTable)(J)->
655 GetValue(HighestKineticEnergy,isOut) ;
656 G4PhysicsLogVector* aVector;
657 aVector = new G4PhysicsLogVector(SmallestRange,
658 BiggestRange,TotBin);
659
660 InvertRangeVector(theRangeTable,
661 theRangeCoeffATable,
662 theRangeCoeffBTable,
663 theRangeCoeffCTable,
664 LowestKineticEnergy,HighestKineticEnergy,TotBin,J, aVector);
665
666 theInverseRangeTable->insert(aVector);
667 }
668 return theInverseRangeTable ;
669}
670
671//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
672
673void G4VEnergyLoss::InvertRangeVector(G4PhysicsTable* theRangeTable,
674 G4PhysicsTable* theRangeCoeffATable,
675 G4PhysicsTable* theRangeCoeffBTable,
676 G4PhysicsTable* theRangeCoeffCTable,
677 G4double LowestKineticEnergy,
678 G4double HighestKineticEnergy,G4int TotBin,
679 G4int materialIndex,G4PhysicsLogVector* aVector)
680// invert range vector for a material
681{
682 G4double LowEdgeRange,A,B,C,discr,KineticEnergy ;
683 G4double RTable = std::exp(std::log(HighestKineticEnergy/LowestKineticEnergy)/TotBin) ;
684 G4double Tbin = LowestKineticEnergy/RTable ;
685 G4double rangebin = 0.0 ;
686 G4int binnumber = -1 ;
687 G4bool isOut ;
688
689 //loop for range values
690 for( G4int i=0; i<TotBin; i++)
691 {
692 LowEdgeRange = aVector->GetLowEdgeEnergy(i) ; //i.e. GetLowEdgeValue(i)
693 if( rangebin < LowEdgeRange )
694 {
695 do
696 {
697 binnumber += 1 ;
698 Tbin *= RTable ;
699 rangebin = (*theRangeTable)(materialIndex)->GetValue(Tbin,isOut) ;
700 }
701 while ((rangebin < LowEdgeRange) && (binnumber < TotBin )) ;
702 }
703
704 if(binnumber == 0)
705 KineticEnergy = LowestKineticEnergy ;
706 else if(binnumber == TotBin-1)
707 KineticEnergy = HighestKineticEnergy ;
708 else
709 {
710 A = (*(*theRangeCoeffATable)(materialIndex))(binnumber-1) ;
711 B = (*(*theRangeCoeffBTable)(materialIndex))(binnumber-1) ;
712 C = (*(*theRangeCoeffCTable)(materialIndex))(binnumber-1) ;
713 if(A==0.)
714 KineticEnergy = (LowEdgeRange -C )/B ;
715 else
716 {
717 discr = B*B - 4.*A*(C-LowEdgeRange);
718 discr = discr>0. ? std::sqrt(discr) : 0.;
719 KineticEnergy = 0.5*(discr-B)/A ;
720 }
721 }
722
723 aVector->PutValue(i,KineticEnergy) ;
724 }
725}
726
727//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
728
729G4PhysicsTable* G4VEnergyLoss::BuildRangeCoeffATable(G4PhysicsTable* theRangeTable,
730 G4PhysicsTable* theRangeCoeffATable,
731 G4double LowestKineticEnergy,
732 G4double HighestKineticEnergy,G4int TotBin)
733// Build tables of coefficients for the energy loss calculation
734// create table for coefficients "A"
735{
736 G4int numOfMaterials = theRangeTable->length();
737
738 if(theRangeCoeffATable)
739 { theRangeCoeffATable->clearAndDestroy();
740 delete theRangeCoeffATable; }
741 theRangeCoeffATable = new G4PhysicsTable(numOfMaterials);
742
743 G4double RTable = std::exp(std::log(HighestKineticEnergy/LowestKineticEnergy)/TotBin) ;
744 G4double R2 = RTable*RTable ;
745 G4double R1 = RTable+1.;
746 G4double w = R1*(RTable-1.)*(RTable-1.);
747 G4double w1 = RTable/w , w2 = -RTable*R1/w , w3 = R2/w ;
748 G4double Ti , Tim , Tip , Ri , Rim , Rip , Value ;
749 G4bool isOut;
750
751 // loop for materials
752 for (G4int J=0; J<numOfMaterials; J++)
753 {
754 G4int binmax=TotBin ;
755 G4PhysicsLinearVector* aVector =
756 new G4PhysicsLinearVector(0.,binmax, TotBin);
757 Ti = LowestKineticEnergy ;
758 G4PhysicsVector* rangeVector= (*theRangeTable)[J];
759
760 for ( G4int i=0; i<TotBin; i++)
761 {
762 Ri = rangeVector->GetValue(Ti,isOut) ;
763 if ( i==0 )
764 Rim = 0. ;
765 else
766 {
767 Tim = Ti/RTable ;
768 Rim = rangeVector->GetValue(Tim,isOut);
769 }
770 if ( i==(TotBin-1))
771 Rip = Ri ;
772 else
773 {
774 Tip = Ti*RTable ;
775 Rip = rangeVector->GetValue(Tip,isOut);
776 }
777 Value = (w1*Rip + w2*Ri + w3*Rim)/(Ti*Ti) ;
778
779 aVector->PutValue(i,Value);
780 Ti = RTable*Ti ;
781 }
782
783 theRangeCoeffATable->insert(aVector);
784 }
785 return theRangeCoeffATable ;
786}
787
788//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
789
790G4PhysicsTable* G4VEnergyLoss::BuildRangeCoeffBTable(G4PhysicsTable* theRangeTable,
791 G4PhysicsTable* theRangeCoeffBTable,
792 G4double LowestKineticEnergy,
793 G4double HighestKineticEnergy,G4int TotBin)
794// Build tables of coefficients for the energy loss calculation
795// create table for coefficients "B"
796{
797 G4int numOfMaterials = theRangeTable->length();
798
799 if(theRangeCoeffBTable)
800 { theRangeCoeffBTable->clearAndDestroy();
801 delete theRangeCoeffBTable; }
802 theRangeCoeffBTable = new G4PhysicsTable(numOfMaterials);
803
804 G4double RTable = std::exp(std::log(HighestKineticEnergy/LowestKineticEnergy)/TotBin) ;
805 G4double R2 = RTable*RTable ;
806 G4double R1 = RTable+1.;
807 G4double w = R1*(RTable-1.)*(RTable-1.);
808 G4double w1 = -R1/w , w2 = R1*(R2+1.)/w , w3 = -R2*R1/w ;
809 G4double Ti , Tim , Tip , Ri , Rim , Rip , Value ;
810 G4bool isOut;
811
812 // loop for materials
813 for (G4int J=0; J<numOfMaterials; J++)
814 {
815 G4int binmax=TotBin ;
816 G4PhysicsLinearVector* aVector =
817 new G4PhysicsLinearVector(0.,binmax, TotBin);
818 Ti = LowestKineticEnergy ;
819 G4PhysicsVector* rangeVector= (*theRangeTable)[J];
820
821 for ( G4int i=0; i<TotBin; i++)
822 {
823 Ri = rangeVector->GetValue(Ti,isOut) ;
824 if ( i==0 )
825 Rim = 0. ;
826 else
827 {
828 Tim = Ti/RTable ;
829 Rim = rangeVector->GetValue(Tim,isOut);
830 }
831 if ( i==(TotBin-1))
832 Rip = Ri ;
833 else
834 {
835 Tip = Ti*RTable ;
836 Rip = rangeVector->GetValue(Tip,isOut);
837 }
838 Value = (w1*Rip + w2*Ri + w3*Rim)/Ti;
839
840 aVector->PutValue(i,Value);
841 Ti = RTable*Ti ;
842 }
843 theRangeCoeffBTable->insert(aVector);
844 }
845 return theRangeCoeffBTable ;
846}
847
848//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
849
850G4PhysicsTable* G4VEnergyLoss::BuildRangeCoeffCTable(G4PhysicsTable* theRangeTable,
851 G4PhysicsTable* theRangeCoeffCTable,
852 G4double LowestKineticEnergy,
853 G4double HighestKineticEnergy,G4int TotBin)
854// Build tables of coefficients for the energy loss calculation
855// create table for coefficients "C"
856{
857 G4int numOfMaterials = theRangeTable->length();
858
859 if(theRangeCoeffCTable)
860 { theRangeCoeffCTable->clearAndDestroy();
861 delete theRangeCoeffCTable; }
862 theRangeCoeffCTable = new G4PhysicsTable(numOfMaterials);
863
864 G4double RTable = std::exp(std::log(HighestKineticEnergy/LowestKineticEnergy)/TotBin) ;
865 G4double R2 = RTable*RTable ;
866 G4double R1 = RTable+1.;
867 G4double w = R1*(RTable-1.)*(RTable-1.);
868 G4double w1 = 1./w , w2 = -RTable*R1/w , w3 = RTable*R2/w ;
869 G4double Ti , Tim , Tip , Ri , Rim , Rip , Value ;
870 G4bool isOut;
871
872 // loop for materials
873 for (G4int J=0; J<numOfMaterials; J++)
874 {
875 G4int binmax=TotBin ;
876 G4PhysicsLinearVector* aVector =
877 new G4PhysicsLinearVector(0.,binmax, TotBin);
878 Ti = LowestKineticEnergy ;
879 G4PhysicsVector* rangeVector= (*theRangeTable)[J];
880
881 for ( G4int i=0; i<TotBin; i++)
882 {
883 Ri = rangeVector->GetValue(Ti,isOut) ;
884 if ( i==0 )
885 Rim = 0. ;
886 else
887 {
888 Tim = Ti/RTable ;
889 Rim = rangeVector->GetValue(Tim,isOut);
890 }
891 if ( i==(TotBin-1))
892 Rip = Ri ;
893 else
894 {
895 Tip = Ti*RTable ;
896 Rip = rangeVector->GetValue(Tip,isOut);
897 }
898 Value = w1*Rip + w2*Ri + w3*Rim ;
899
900 aVector->PutValue(i,Value);
901 Ti = RTable*Ti ;
902 }
903 theRangeCoeffCTable->insert(aVector);
904 }
905 return theRangeCoeffCTable ;
906}
907
908//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
909
910G4double G4VEnergyLoss::GetLossWithFluct(const G4DynamicParticle* aParticle,
911 const G4MaterialCutsCouple* couple,
912 G4double ChargeSquare,
913 G4double MeanLoss,
914 G4double step )
915// calculate actual loss from the mean loss
916// The model used to get the fluctuation is essentially the same as in Glandz in Geant3.
917{
918 const G4double minLoss = 1.*eV ;
919 const G4double probLim = 0.01 ;
920 const G4double sumaLim = -std::log(probLim) ;
921 const G4double alim=10.;
922 const G4double kappa = 10. ;
923 const G4double factor = twopi_mc2_rcl2 ;
924 const G4Material* aMaterial = couple->GetMaterial();
925
926 // check if the material has changed ( cache mechanism)
927
928 if (aMaterial != lastMaterial)
929 {
930 lastMaterial = aMaterial;
931 imat = couple->GetIndex();
932 f1Fluct = aMaterial->GetIonisation()->GetF1fluct();
933 f2Fluct = aMaterial->GetIonisation()->GetF2fluct();
934 e1Fluct = aMaterial->GetIonisation()->GetEnergy1fluct();
935 e2Fluct = aMaterial->GetIonisation()->GetEnergy2fluct();
936 e1LogFluct = aMaterial->GetIonisation()->GetLogEnergy1fluct();
937 e2LogFluct = aMaterial->GetIonisation()->GetLogEnergy2fluct();
938 rateFluct = aMaterial->GetIonisation()->GetRateionexcfluct();
939 ipotFluct = aMaterial->GetIonisation()->GetMeanExcitationEnergy();
940 ipotLogFluct = aMaterial->GetIonisation()->GetLogMeanExcEnergy();
941 }
942 G4double threshold,w1,w2,C,
943 beta2,suma,e0,loss,lossc ,w,electronDensity;
944 G4double a1,a2,a3;
945 G4double p1,p2,p3 ;
946 G4int nb;
947 G4double Corrfac, na,alfa,rfac,namean,sa,alfa1,ea,sea;
948 G4double dp3;
949 G4double siga ;
950
951 // shortcut for very very small loss
952 if(MeanLoss < minLoss) return MeanLoss ;
953
954 // get particle data
955 G4double Tkin = aParticle->GetKineticEnergy();
956 ParticleMass = aParticle->GetMass() ;
957
958 threshold = (*((G4ProductionCutsTable::GetProductionCutsTable())
959 ->GetEnergyCutsVector(1)))[imat];
960 G4double rmass = electron_mass_c2/ParticleMass;
961 G4double tau = Tkin/ParticleMass, tau1 = tau+1., tau2 = tau*(tau+2.);
962 G4double Tm = 2.*electron_mass_c2*tau2/(1.+2.*tau1*rmass+rmass*rmass);
963
964 if(Tm > threshold) Tm = threshold;
965
966 beta2 = tau2/(tau1*tau1);
967
968 // Gaussian fluctuation ?
969 if(MeanLoss >= kappa*Tm || Tm < kappa*ipotFluct)
970 {
971 electronDensity = aMaterial->GetElectronDensity() ;
972 siga = std::sqrt(Tm*(1.0-0.5*beta2)*step*
973 factor*electronDensity*ChargeSquare/beta2) ;
974 do {
975 loss = G4RandGauss::shoot(MeanLoss,siga) ;
976 } while (loss < 0. || loss > 2.*MeanLoss);
977 return loss;
978 }
979
980 w1 = Tm/ipotFluct;
981 w2 = std::log(2.*electron_mass_c2*tau2);
982
983 C = MeanLoss*(1.-rateFluct)/(w2-ipotLogFluct-beta2);
984
985 a1 = C*f1Fluct*(w2-e1LogFluct-beta2)/e1Fluct;
986 a2 = C*f2Fluct*(w2-e2LogFluct-beta2)/e2Fluct;
987 a3 = rateFluct*MeanLoss*(Tm-ipotFluct)/(ipotFluct*Tm*std::log(w1));
988 if(a1 < 0.) a1 = 0.;
989 if(a2 < 0.) a2 = 0.;
990 if(a3 < 0.) a3 = 0.;
991
992 suma = a1+a2+a3;
993
994 loss = 0. ;
995
996 if(suma < sumaLim) // very small Step
997 {
998 e0 = aMaterial->GetIonisation()->GetEnergy0fluct();
999
1000 if(Tm == ipotFluct)
1001 {
1002 a3 = MeanLoss/e0;
1003
1004 if(a3>alim)
1005 {
1006 siga=std::sqrt(a3) ;
1007 p3 = std::max(0.,G4RandGauss::shoot(a3,siga)+0.5);
1008 }
1009 else
1010 p3 = G4float(G4Poisson(a3));
1011
1012 loss = p3*e0 ;
1013
1014 if(p3 > 0)
1015 loss += (1.-2.*G4UniformRand())*e0 ;
1016
1017 }
1018 else
1019 {
1020 Tm = Tm-ipotFluct+e0 ;
1021 a3 = MeanLoss*(Tm-e0)/(Tm*e0*std::log(Tm/e0));
1022
1023 if(a3>alim)
1024 {
1025 siga=std::sqrt(a3) ;
1026 p3 = std::max(0.,G4RandGauss::shoot(a3,siga)+0.5);
1027 }
1028 else
1029 p3 = G4float(G4Poisson(a3));
1030
1031 if(p3 > 0)
1032 {
1033 w = (Tm-e0)/Tm ;
1034 if(p3 > G4float(nmaxCont2))
1035 {
1036 dp3 = p3 ;
1037 Corrfac = dp3/G4float(nmaxCont2) ;
1038 p3 = G4float(nmaxCont2) ;
1039 }
1040 else
1041 Corrfac = 1. ;
1042
1043 for(G4int i=0; i<G4int(p3); i++) loss += 1./(1.-w*G4UniformRand()) ;
1044 loss *= e0*Corrfac ;
1045 }
1046 }
1047 }
1048
1049 else // not so small Step
1050 {
1051 // excitation type 1
1052 if(a1>alim)
1053 {
1054 siga=std::sqrt(a1) ;
1055 p1 = std::max(0.,G4RandGauss::shoot(a1,siga)+0.5);
1056 }
1057 else
1058 p1 = G4float(G4Poisson(a1));
1059
1060 // excitation type 2
1061 if(a2>alim)
1062 {
1063 siga=std::sqrt(a2) ;
1064 p2 = std::max(0.,G4RandGauss::shoot(a2,siga)+0.5);
1065 }
1066 else
1067 p2 = G4float(G4Poisson(a2));
1068
1069 loss = p1*e1Fluct+p2*e2Fluct;
1070 // smearing to avoid unphysical peaks
1071 if(p2 > 0)
1072 loss += (1.-2.*G4UniformRand())*e2Fluct;
1073 else if (loss>0.)
1074 loss += (1.-2.*G4UniformRand())*e1Fluct;
1075
1076 // ionisation .......................................
1077 if(a3 > 0.)
1078 {
1079 if(a3>alim)
1080 {
1081 siga=std::sqrt(a3) ;
1082 p3 = std::max(0.,G4RandGauss::shoot(a3,siga)+0.5);
1083 }
1084 else
1085 p3 = G4float(G4Poisson(a3));
1086
1087 lossc = 0.;
1088 if(p3 > 0)
1089 {
1090 na = 0.;
1091 alfa = 1.;
1092 if (p3 > G4float(nmaxCont2))
1093 {
1094 dp3 = p3;
1095 rfac = dp3/(G4float(nmaxCont2)+dp3);
1096 namean = p3*rfac;
1097 sa = G4float(nmaxCont1)*rfac;
1098 na = G4RandGauss::shoot(namean,sa);
1099 if (na > 0.)
1100 {
1101 alfa = w1*(G4float(nmaxCont2)+p3)/(w1*G4float(nmaxCont2)+p3);
1102 alfa1 = alfa*std::log(alfa)/(alfa-1.);
1103 ea = na*ipotFluct*alfa1;
1104 sea = ipotFluct*std::sqrt(na*(alfa-alfa1*alfa1));
1105 lossc += G4RandGauss::shoot(ea,sea);
1106 }
1107 }
1108
1109 nb = G4int(p3-na);
1110 if (nb > 0)
1111 {
1112 w2 = alfa*ipotFluct;
1113 w = (Tm-w2)/Tm;
1114 for (G4int k=0; k<nb; k++) lossc += w2/(1.-w*G4UniformRand());
1115 }
1116 }
1117 loss += lossc;
1118 }
1119 }
1120
1121 if( loss < 0.)
1122 loss = 0.;
1123
1124 return loss ;
1125}
1126
1127//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1128
1129G4bool G4VEnergyLoss::EqualCutVectors( G4double* vec1, G4double* vec2 )
1130{
1131 if ( (vec1==0 ) || (vec2==0) ) return false;
1132
1133 G4bool flag = true;
1134
1135 for (size_t j=0; flag && j<G4Material::GetNumberOfMaterials(); j++){
1136 flag = (vec1[j] == vec2[j]);
1137 }
1138
1139 return flag;
1140}
1141
1142//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1143
1144G4double* G4VEnergyLoss::CopyCutVectors( G4double* dest, G4double* source )
1145{
1146 if ( dest != 0) delete [] dest;
1147 dest = new G4double [G4Material::GetNumberOfMaterials()];
1148 for (size_t j=0; j<G4Material::GetNumberOfMaterials(); j++){
1149 dest[j] = source[j];
1150 }
1151 return dest;
1152}
1153
1154//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1155
1156G4bool G4VEnergyLoss::CutsWhereModified()
1157{
1158 G4bool wasModified = false;
1159 const G4ProductionCutsTable* theCoupleTable=
1160 G4ProductionCutsTable::GetProductionCutsTable();
1161 size_t numOfCouples = theCoupleTable->GetTableSize();
1162
1163 for (size_t j=0; j<numOfCouples; j++){
1164 if (theCoupleTable->GetMaterialCutsCouple(j)->IsRecalcNeeded()) {
1165 wasModified = true;
1166 break;
1167 }
1168 }
1169 return wasModified;
1170}
1171
1172//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
Note: See TracBrowser for help on using the repository browser.