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