source: trunk/source/processes/electromagnetic/utils/src/G4VEnergyLoss.cc

Last change on this file was 1337, checked in by garnier, 14 years ago

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