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

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

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