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