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