source: trunk/source/processes/electromagnetic/lowenergy/src/G4VeLowEnergyLoss.cc@ 1058

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

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