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