source: trunk/source/processes/electromagnetic/pii/src/G4hRDEnergyLoss.cc

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

update to last version 4.9.4
File size: 34.7 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: G4hRDEnergyLoss.cc,v 1.3 2010/11/25 19:49:43 pia Exp $
28// GEANT4 tag $Name: geant4-09-04-ref-00 $
29//
30// -----------------------------------------------------------
31//      GEANT 4 class implementation file
32//
33//      History: based on object model of
34//      2nd December 1995, G.Cosmo
35//      ---------- G4hEnergyLoss physics process -----------
36//                by Laszlo Urban, 30 May 1997
37//
38// **************************************************************
39// It is the first implementation of the NEW UNIFIED ENERGY LOSS PROCESS.
40// It calculates the energy loss of charged hadrons.
41// **************************************************************
42//
43// 7/10/98    bug fixes + some cleanup , L.Urban
44// 22/10/98   cleanup , L.Urban
45// 07/12/98   works for ions as well+ bug corrected, L.Urban
46// 02/02/99   several bugs fixed, L.Urban
47// 31/03/00   rename to lowenergy as G4hRDEnergyLoss.cc V.Ivanchenko
48// 05/11/00   new method to calculate particle ranges
49// 10/05/01   V.Ivanchenko Clean up againist Linux compilation with -Wall
50// 23/11/01   V.Ivanchenko Move static member-functions from header to source
51// 28/10/02   V.Ivanchenko Optimal binning for dE/dx
52// 21/01/03   V.Ivanchenko Cut per region
53// 23/01/03   V.Ivanchenko Fix in table build
54
55// 31 Jul 2008 MGP     Short term supply of energy loss of hadrons through clone of
56//                     former G4hLowEnergyLoss (with some initial cleaning)
57//                     To be replaced by reworked class to deal with condensed/discrete
58//                     issues properly
59
60// 25 Nov 2010 MGP     Added some protections for FPE (mostly division by 0)
61//                     The whole energy loss domain would profit from a design iteration
62// --------------------------------------------------------------
63
64#include "G4hRDEnergyLoss.hh"
65#include "G4EnergyLossTables.hh"
66#include "G4Poisson.hh"
67#include "G4ProductionCutsTable.hh"
68
69//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
70
71// Initialisation of static members ******************************************
72// contributing processes : ion.loss ->NumberOfProcesses is initialized
73//   to 1 . YOU DO NOT HAVE TO CHANGE this variable for a 'normal' run.
74// You have to change NumberOfProcesses
75// if you invent a new process contributing to the cont. energy loss,
76//   NumberOfProcesses should be 2 in this case,
77//  or for debugging purposes.
78//  The NumberOfProcesses data member can be changed using the (public static)
79//  functions Get/Set/Plus/MinusNumberOfProcesses (see G4hRDEnergyLoss.hh)
80
81
82// The following vectors have a fixed dimension this means that if they are
83// filled in with more than 100 elements will corrupt the memory.
84G4int G4hRDEnergyLoss::NumberOfProcesses = 1 ;
85
86G4int            G4hRDEnergyLoss::CounterOfProcess = 0 ;
87G4PhysicsTable** G4hRDEnergyLoss::RecorderOfProcess =
88  new G4PhysicsTable*[100] ;
89
90G4int            G4hRDEnergyLoss::CounterOfpProcess = 0 ;
91G4PhysicsTable** G4hRDEnergyLoss::RecorderOfpProcess =
92  new G4PhysicsTable*[100] ;
93
94G4int            G4hRDEnergyLoss::CounterOfpbarProcess = 0 ;
95G4PhysicsTable** G4hRDEnergyLoss::RecorderOfpbarProcess =
96  new G4PhysicsTable*[100] ;
97
98G4PhysicsTable* G4hRDEnergyLoss::theDEDXpTable = 0 ;
99G4PhysicsTable* G4hRDEnergyLoss::theDEDXpbarTable = 0 ;
100G4PhysicsTable* G4hRDEnergyLoss::theRangepTable = 0 ;
101G4PhysicsTable* G4hRDEnergyLoss::theRangepbarTable = 0 ;
102G4PhysicsTable* G4hRDEnergyLoss::theInverseRangepTable = 0 ;
103G4PhysicsTable* G4hRDEnergyLoss::theInverseRangepbarTable = 0 ;
104G4PhysicsTable* G4hRDEnergyLoss::theLabTimepTable = 0 ;
105G4PhysicsTable* G4hRDEnergyLoss::theLabTimepbarTable = 0 ;
106G4PhysicsTable* G4hRDEnergyLoss::theProperTimepTable = 0 ;
107G4PhysicsTable* G4hRDEnergyLoss::theProperTimepbarTable = 0 ;
108
109G4PhysicsTable* G4hRDEnergyLoss::thepRangeCoeffATable = 0 ;
110G4PhysicsTable* G4hRDEnergyLoss::thepRangeCoeffBTable = 0 ;
111G4PhysicsTable* G4hRDEnergyLoss::thepRangeCoeffCTable = 0 ;
112G4PhysicsTable* G4hRDEnergyLoss::thepbarRangeCoeffATable = 0 ;
113G4PhysicsTable* G4hRDEnergyLoss::thepbarRangeCoeffBTable = 0 ;
114G4PhysicsTable* G4hRDEnergyLoss::thepbarRangeCoeffCTable = 0 ;
115
116G4PhysicsTable* G4hRDEnergyLoss::theDEDXTable = 0 ;
117G4PhysicsTable* G4hRDEnergyLoss::theRangeTable = 0 ;
118G4PhysicsTable* G4hRDEnergyLoss::theInverseRangeTable = 0 ;
119G4PhysicsTable* G4hRDEnergyLoss::theLabTimeTable = 0 ;
120G4PhysicsTable* G4hRDEnergyLoss::theProperTimeTable = 0 ;
121
122G4PhysicsTable* G4hRDEnergyLoss::theRangeCoeffATable = 0 ;
123G4PhysicsTable* G4hRDEnergyLoss::theRangeCoeffBTable = 0 ;
124G4PhysicsTable* G4hRDEnergyLoss::theRangeCoeffCTable = 0 ;
125
126//const G4Proton* G4hRDEnergyLoss::theProton=G4Proton::Proton() ;
127//const G4AntiProton* G4hRDEnergyLoss::theAntiProton=G4AntiProton::AntiProton() ;
128
129G4double G4hRDEnergyLoss::ParticleMass;
130G4double G4hRDEnergyLoss::ptableElectronCutInRange = 0.0*mm ;
131G4double G4hRDEnergyLoss::pbartableElectronCutInRange = 0.0*mm ;
132
133G4double G4hRDEnergyLoss::Mass,
134  G4hRDEnergyLoss::taulow, 
135  G4hRDEnergyLoss::tauhigh, 
136  G4hRDEnergyLoss::ltaulow, 
137  G4hRDEnergyLoss::ltauhigh; 
138
139G4double G4hRDEnergyLoss::dRoverRange = 0.20 ;
140G4double G4hRDEnergyLoss::finalRange = 200.*micrometer ;
141
142G4double     G4hRDEnergyLoss::c1lim = dRoverRange ;
143G4double     G4hRDEnergyLoss::c2lim = 2.*(1.-dRoverRange)*finalRange ;
144G4double     G4hRDEnergyLoss::c3lim = -(1.-dRoverRange)*finalRange*finalRange;
145
146G4double         G4hRDEnergyLoss::Charge ;   
147
148G4bool   G4hRDEnergyLoss::rndmStepFlag   = false ;
149G4bool   G4hRDEnergyLoss::EnlossFlucFlag = true ;
150
151G4double G4hRDEnergyLoss::LowestKineticEnergy = 10.*eV;
152G4double G4hRDEnergyLoss::HighestKineticEnergy= 100.*GeV;
153G4int    G4hRDEnergyLoss::TotBin = 360;
154G4double G4hRDEnergyLoss::RTable,G4hRDEnergyLoss::LOGRTable;
155
156//--------------------------------------------------------------------------------
157
158// constructor and destructor
159 
160G4hRDEnergyLoss::G4hRDEnergyLoss(const G4String& processName)
161  : G4VContinuousDiscreteProcess (processName),
162    MaxExcitationNumber (1.e6),
163    probLimFluct (0.01),
164    nmaxDirectFluct (100),
165    nmaxCont1(4),
166    nmaxCont2(16),
167    theLossTable(0),
168    linLossLimit(0.05),
169    MinKineticEnergy(0.0) 
170{;}
171
172//--------------------------------------------------------------------------------
173
174G4hRDEnergyLoss::~G4hRDEnergyLoss() 
175{
176  if(theLossTable) {
177    theLossTable->clearAndDestroy();
178    delete theLossTable;
179  }
180}
181
182//--------------------------------------------------------------------------------
183
184G4int G4hRDEnergyLoss::GetNumberOfProcesses()   
185{   
186  return NumberOfProcesses; 
187} 
188
189//--------------------------------------------------------------------------------
190
191void G4hRDEnergyLoss::SetNumberOfProcesses(G4int number)
192{
193  NumberOfProcesses=number; 
194} 
195
196//--------------------------------------------------------------------------------
197
198void G4hRDEnergyLoss::PlusNumberOfProcesses()
199{ 
200  NumberOfProcesses++; 
201} 
202
203//--------------------------------------------------------------------------------
204
205void G4hRDEnergyLoss::MinusNumberOfProcesses()
206{ 
207  NumberOfProcesses--; 
208} 
209
210//--------------------------------------------------------------------------------
211
212void G4hRDEnergyLoss::SetdRoverRange(G4double value) 
213{
214  dRoverRange = value;
215}
216
217//--------------------------------------------------------------------------------
218
219void G4hRDEnergyLoss::SetRndmStep (G4bool   value) 
220{
221  rndmStepFlag = value;
222}
223
224//--------------------------------------------------------------------------------
225
226void G4hRDEnergyLoss::SetEnlossFluc (G4bool value) 
227{
228  EnlossFlucFlag = value;
229}
230
231//--------------------------------------------------------------------------------
232
233void G4hRDEnergyLoss::SetStepFunction (G4double c1, G4double c2)
234{
235  dRoverRange = c1; 
236  finalRange = c2;
237  c1lim=dRoverRange;
238  c2lim=2.*(1-dRoverRange)*finalRange;
239  c3lim=-(1.-dRoverRange)*finalRange*finalRange;
240}
241
242//--------------------------------------------------------------------------------
243 
244void G4hRDEnergyLoss::BuildDEDXTable(
245                                     const G4ParticleDefinition& aParticleType)
246{
247  //  calculate data members TotBin,LOGRTable,RTable first
248
249  const G4ProductionCutsTable* theCoupleTable=
250    G4ProductionCutsTable::GetProductionCutsTable();
251  size_t numOfCouples = theCoupleTable->GetTableSize();
252
253  // create/fill proton or antiproton tables depending on the charge
254  Charge = aParticleType.GetPDGCharge()/eplus;
255  ParticleMass = aParticleType.GetPDGMass() ;
256
257  if (Charge>0.) {theDEDXTable= theDEDXpTable;}
258  else           {theDEDXTable= theDEDXpbarTable;}
259
260  if( ((Charge>0.) && (theDEDXTable==0)) ||
261      ((Charge<0.) && (theDEDXTable==0)) 
262      )
263    {
264
265      // Build energy loss table as a sum of the energy loss due to the
266      //              different processes.
267      if( Charge >0.)
268        {
269          RecorderOfProcess=RecorderOfpProcess;
270          CounterOfProcess=CounterOfpProcess;
271
272          if(CounterOfProcess == NumberOfProcesses)
273            {
274              if(theDEDXpTable)
275                { theDEDXpTable->clearAndDestroy();
276                  delete theDEDXpTable; }
277              theDEDXpTable = new G4PhysicsTable(numOfCouples);
278              theDEDXTable = theDEDXpTable;
279            }
280        }
281      else
282        {
283          RecorderOfProcess=RecorderOfpbarProcess;
284          CounterOfProcess=CounterOfpbarProcess;
285
286          if(CounterOfProcess == NumberOfProcesses)
287            {
288              if(theDEDXpbarTable)
289                { theDEDXpbarTable->clearAndDestroy();
290                  delete theDEDXpbarTable; }
291              theDEDXpbarTable = new G4PhysicsTable(numOfCouples);
292              theDEDXTable = theDEDXpbarTable;
293            }
294        }
295
296      if(CounterOfProcess == NumberOfProcesses)
297        {
298          //  loop for materials
299          G4double LowEdgeEnergy , Value ;
300          G4bool isOutRange ;
301          G4PhysicsTable* pointer ;
302
303          for (size_t J=0; J<numOfCouples; J++)
304            {
305              // create physics vector and fill it
306              G4PhysicsLogVector* aVector = new G4PhysicsLogVector(
307                                                                   LowestKineticEnergy, HighestKineticEnergy, TotBin);
308
309              // loop for the kinetic energy
310              for (G4int i=0; i<TotBin; i++)
311                {
312                  LowEdgeEnergy = aVector->GetLowEdgeEnergy(i) ;
313                  Value = 0. ;
314
315                  // loop for the contributing processes
316                  for (G4int process=0; process < NumberOfProcesses; process++)
317                    {
318                      pointer= RecorderOfProcess[process];
319                      Value += (*pointer)[J]->
320                        GetValue(LowEdgeEnergy,isOutRange) ;
321                    }
322
323                  aVector->PutValue(i,Value) ;
324                }
325
326              theDEDXTable->insert(aVector) ;
327            }
328
329          //  reset counter to zero ..................
330          if( Charge >0.)
331            CounterOfpProcess=0 ;
332          else
333            CounterOfpbarProcess=0 ;
334
335          // Build range table
336          BuildRangeTable( aParticleType);
337
338          // Build lab/proper time tables
339          BuildTimeTables( aParticleType) ;
340
341          // Build coeff tables for the energy loss calculation
342          BuildRangeCoeffATable( aParticleType);
343          BuildRangeCoeffBTable( aParticleType);
344          BuildRangeCoeffCTable( aParticleType);
345
346          // invert the range table
347
348          BuildInverseRangeTable(aParticleType);
349        }
350    }
351  // make the energy loss and the range table available
352
353  G4EnergyLossTables::Register(&aParticleType,
354                               (Charge>0)?
355                               theDEDXpTable: theDEDXpbarTable,
356                               (Charge>0)?
357                               theRangepTable: theRangepbarTable,
358                               (Charge>0)?
359                               theInverseRangepTable: theInverseRangepbarTable,
360                               (Charge>0)?
361                               theLabTimepTable: theLabTimepbarTable,
362                               (Charge>0)?
363                               theProperTimepTable: theProperTimepbarTable,
364                               LowestKineticEnergy, HighestKineticEnergy,
365                               proton_mass_c2/aParticleType.GetPDGMass(),
366                               TotBin);
367
368}
369
370//--------------------------------------------------------------------------------
371
372void G4hRDEnergyLoss::BuildRangeTable(
373                                      const G4ParticleDefinition& aParticleType)
374// Build range table from the energy loss table
375{
376  Mass = aParticleType.GetPDGMass();
377
378  const G4ProductionCutsTable* theCoupleTable=
379    G4ProductionCutsTable::GetProductionCutsTable();
380  size_t numOfCouples = theCoupleTable->GetTableSize();
381
382  if( Charge >0.)
383    {
384      if(theRangepTable)
385        { theRangepTable->clearAndDestroy();
386          delete theRangepTable; }
387      theRangepTable = new G4PhysicsTable(numOfCouples);
388      theRangeTable = theRangepTable ;
389    }
390  else
391    {
392      if(theRangepbarTable)
393        { theRangepbarTable->clearAndDestroy();
394          delete theRangepbarTable; }
395      theRangepbarTable = new G4PhysicsTable(numOfCouples);
396      theRangeTable = theRangepbarTable ;
397    }
398
399  // loop for materials
400
401  for (size_t J=0;  J<numOfCouples; J++)
402    {
403      G4PhysicsLogVector* aVector;
404      aVector = new G4PhysicsLogVector(LowestKineticEnergy,
405                                       HighestKineticEnergy,TotBin);
406
407      BuildRangeVector(J, aVector);
408      theRangeTable->insert(aVector);
409    }
410}
411
412//--------------------------------------------------------------------------------
413
414void G4hRDEnergyLoss::BuildTimeTables(
415                                      const G4ParticleDefinition& aParticleType)
416{
417
418  const G4ProductionCutsTable* theCoupleTable=
419    G4ProductionCutsTable::GetProductionCutsTable();
420  size_t numOfCouples = theCoupleTable->GetTableSize();
421
422  if(&aParticleType == G4Proton::Proton())
423    {
424      if(theLabTimepTable)
425        { theLabTimepTable->clearAndDestroy();
426          delete theLabTimepTable; }
427      theLabTimepTable = new G4PhysicsTable(numOfCouples);
428      theLabTimeTable = theLabTimepTable ;
429
430      if(theProperTimepTable)
431        { theProperTimepTable->clearAndDestroy();
432          delete theProperTimepTable; }
433      theProperTimepTable = new G4PhysicsTable(numOfCouples);
434      theProperTimeTable = theProperTimepTable ;
435    }
436
437  if(&aParticleType == G4AntiProton::AntiProton())
438    {
439      if(theLabTimepbarTable)
440        { theLabTimepbarTable->clearAndDestroy();
441          delete theLabTimepbarTable; }
442      theLabTimepbarTable = new G4PhysicsTable(numOfCouples);
443      theLabTimeTable = theLabTimepbarTable ;
444
445      if(theProperTimepbarTable)
446        { theProperTimepbarTable->clearAndDestroy();
447          delete theProperTimepbarTable; }
448      theProperTimepbarTable = new G4PhysicsTable(numOfCouples);
449      theProperTimeTable = theProperTimepbarTable ;
450    }
451
452  for (size_t J=0;  J<numOfCouples; J++)
453    {
454      G4PhysicsLogVector* aVector;
455      G4PhysicsLogVector* bVector;
456
457      aVector = new G4PhysicsLogVector(LowestKineticEnergy,
458                                       HighestKineticEnergy,TotBin);
459
460      BuildLabTimeVector(J, aVector);
461      theLabTimeTable->insert(aVector);
462
463      bVector = new G4PhysicsLogVector(LowestKineticEnergy,
464                                       HighestKineticEnergy,TotBin);
465
466      BuildProperTimeVector(J, bVector);
467      theProperTimeTable->insert(bVector);
468    }
469}
470
471//--------------------------------------------------------------------------------
472
473void G4hRDEnergyLoss::BuildRangeVector(G4int materialIndex,
474                                       G4PhysicsLogVector* rangeVector)
475//  create range vector for a material
476{
477  G4bool isOut;
478  G4PhysicsVector* physicsVector= (*theDEDXTable)[materialIndex];
479  G4double energy1 = rangeVector->GetLowEdgeEnergy(0);
480  G4double dedx    = physicsVector->GetValue(energy1,isOut);
481  G4double range   = 0.5*energy1/dedx;
482  rangeVector->PutValue(0,range);
483  G4int n = 100;
484  G4double del = 1.0/(G4double)n ;
485
486  for (G4int j=1; j<TotBin; j++) {
487
488    G4double energy2 = rangeVector->GetLowEdgeEnergy(j);
489    G4double de = (energy2 - energy1) * del ;
490    G4double dedx1 = dedx ;
491
492    for (G4int i=1; i<n; i++) {
493      G4double energy = energy1 + i*de ;
494      G4double dedx2  = physicsVector->GetValue(energy,isOut);
495      range  += 0.5*de*(1.0/dedx1 + 1.0/dedx2);
496      dedx1   = dedx2;
497    }
498    rangeVector->PutValue(j,range);
499    dedx = dedx1 ;
500    energy1 = energy2 ;
501  }
502}   
503
504//--------------------------------------------------------------------------------
505
506void G4hRDEnergyLoss::BuildLabTimeVector(G4int materialIndex,
507                                         G4PhysicsLogVector* timeVector)
508//  create lab time vector for a material
509{
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+LabTimeIntLog(physicsVector,nbin);
544        }
545      timeVector->PutValue(i,Value);
546      oldValue = Value ;
547      tauold = tau ;
548    } while (tau<=taulim) ;
549
550  i += 1 ;
551  for (G4int j=i; j<TotBin; j++)
552    {
553      LowEdgeEnergy = timeVector->GetLowEdgeEnergy(j);
554      tau = LowEdgeEnergy/ParticleMass ;
555      ltaulow = std::log(tauold);
556      ltauhigh = std::log(tau);
557      Value = oldValue+LabTimeIntLog(physicsVector,nbin);
558      timeVector->PutValue(j,Value);
559      oldValue = Value ;
560      tauold = tau ;
561    }
562}
563
564//--------------------------------------------------------------------------------
565
566void G4hRDEnergyLoss::BuildProperTimeVector(G4int materialIndex,
567                                            G4PhysicsLogVector* timeVector)
568//  create proper time vector for a material
569{
570  G4int nbin=100;
571  G4bool isOut;
572  G4double tlim=5.*keV,parlowen=0.4,ppar=0.5-parlowen ;
573  G4double losslim,clim,taulim,timelim,ltaulim,ltaumax,
574    LowEdgeEnergy,tau,Value ;
575
576  G4PhysicsVector* physicsVector= (*theDEDXTable)[materialIndex];
577 
578  // low energy part first...
579  losslim = physicsVector->GetValue(tlim,isOut);
580  taulim=tlim/ParticleMass ;
581  clim=std::sqrt(ParticleMass*tlim/2.)/(c_light*losslim*ppar) ;
582  ltaulim = std::log(taulim);
583  ltaumax = std::log(HighestKineticEnergy/ParticleMass) ;
584
585  G4int i=-1;
586  G4double oldValue = 0. ;
587  G4double tauold ;
588  do
589    {
590      i += 1 ;
591      LowEdgeEnergy = timeVector->GetLowEdgeEnergy(i);
592      tau = LowEdgeEnergy/ParticleMass ;
593      if ( tau <= taulim )
594        {
595          Value = clim*std::exp(ppar*std::log(tau/taulim)) ;
596        }
597      else
598        {
599          timelim=clim ;
600          ltaulow = std::log(taulim);
601          ltauhigh = std::log(tau);
602          Value = timelim+ProperTimeIntLog(physicsVector,nbin);
603        }
604      timeVector->PutValue(i,Value);
605      oldValue = Value ;
606      tauold = tau ;
607    } while (tau<=taulim) ;
608
609  i += 1 ;
610  for (G4int j=i; j<TotBin; j++)
611    {
612      LowEdgeEnergy = timeVector->GetLowEdgeEnergy(j);
613      tau = LowEdgeEnergy/ParticleMass ;
614      ltaulow = std::log(tauold);
615      ltauhigh = std::log(tau);
616      Value = oldValue+ProperTimeIntLog(physicsVector,nbin);
617      timeVector->PutValue(j,Value);
618      oldValue = Value ;
619      tauold = tau ;
620    }
621}
622
623//--------------------------------------------------------------------------------
624
625G4double G4hRDEnergyLoss::RangeIntLin(G4PhysicsVector* physicsVector,
626                                      G4int nbin)
627//  num. integration, linear binning
628{
629  G4double dtau,Value,taui,ti,lossi,ci;
630  G4bool isOut;
631  dtau = (tauhigh-taulow)/nbin;
632  Value = 0.;
633
634  for (G4int i=0; i<=nbin; i++)
635    {
636      taui = taulow + dtau*i ;
637      ti = Mass*taui;
638      lossi = physicsVector->GetValue(ti,isOut);
639      if(i==0)
640        ci=0.5;
641      else
642        {
643          if(i<nbin)
644            ci=1.;
645          else
646            ci=0.5;
647        }
648      Value += ci/lossi;
649    }
650  Value *= Mass*dtau;
651  return Value;
652}
653
654//--------------------------------------------------------------------------------
655
656G4double G4hRDEnergyLoss::RangeIntLog(G4PhysicsVector* physicsVector,
657                                      G4int nbin)
658//  num. integration, logarithmic binning
659{
660  G4double ltt,dltau,Value,ui,taui,ti,lossi,ci;
661  G4bool isOut;
662  ltt = ltauhigh-ltaulow;
663  dltau = ltt/nbin;
664  Value = 0.;
665
666  for (G4int i=0; i<=nbin; i++)
667    {
668      ui = ltaulow+dltau*i;
669      taui = std::exp(ui);
670      ti = Mass*taui;
671      lossi = physicsVector->GetValue(ti,isOut);
672      if(i==0)
673        ci=0.5;
674      else
675        {
676          if(i<nbin)
677            ci=1.;
678          else
679            ci=0.5;
680        }
681      Value += ci*taui/lossi;
682    }
683  Value *= Mass*dltau;
684  return Value;
685}
686
687//--------------------------------------------------------------------------------
688
689G4double G4hRDEnergyLoss::LabTimeIntLog(G4PhysicsVector* physicsVector,
690                                        G4int nbin)
691//  num. integration, logarithmic binning
692{
693  G4double ltt,dltau,Value,ui,taui,ti,lossi,ci;
694  G4bool isOut;
695  ltt = ltauhigh-ltaulow;
696  dltau = ltt/nbin;
697  Value = 0.;
698
699  for (G4int i=0; i<=nbin; i++)
700    {
701      ui = ltaulow+dltau*i;
702      taui = std::exp(ui);
703      ti = ParticleMass*taui;
704      lossi = physicsVector->GetValue(ti,isOut);
705      if(i==0)
706        ci=0.5;
707      else
708        {
709          if(i<nbin)
710            ci=1.;
711          else
712            ci=0.5;
713        }
714      Value += ci*taui*(ti+ParticleMass)/(std::sqrt(ti*(ti+2.*ParticleMass))*lossi);
715    }
716  Value *= ParticleMass*dltau/c_light;
717  return Value;
718}
719
720//--------------------------------------------------------------------------------
721
722G4double G4hRDEnergyLoss::ProperTimeIntLog(G4PhysicsVector* physicsVector,
723                                           G4int nbin)
724//  num. integration, logarithmic binning
725{
726  G4double ltt,dltau,Value,ui,taui,ti,lossi,ci;
727  G4bool isOut;
728  ltt = ltauhigh-ltaulow;
729  dltau = ltt/nbin;
730  Value = 0.;
731
732  for (G4int i=0; i<=nbin; i++)
733    {
734      ui = ltaulow+dltau*i;
735      taui = std::exp(ui);
736      ti = ParticleMass*taui;
737      lossi = physicsVector->GetValue(ti,isOut);
738      if(i==0)
739        ci=0.5;
740      else
741        {
742          if(i<nbin)
743            ci=1.;
744          else
745            ci=0.5;
746        }
747      Value += ci*taui*ParticleMass/(std::sqrt(ti*(ti+2.*ParticleMass))*lossi);
748    }
749  Value *= ParticleMass*dltau/c_light;
750  return Value;
751}
752
753//--------------------------------------------------------------------------------
754
755void G4hRDEnergyLoss::BuildRangeCoeffATable(
756                                            const G4ParticleDefinition& )
757// Build tables of coefficients for the energy loss calculation
758//  create table for coefficients "A"
759{
760
761  G4int numOfCouples = G4ProductionCutsTable::GetProductionCutsTable()->GetTableSize();
762
763  if(Charge>0.)
764    {
765      if(thepRangeCoeffATable)
766        { thepRangeCoeffATable->clearAndDestroy();
767          delete thepRangeCoeffATable; }
768      thepRangeCoeffATable = new G4PhysicsTable(numOfCouples);
769      theRangeCoeffATable = thepRangeCoeffATable ;
770      theRangeTable = theRangepTable ;
771    }
772  else
773    {
774      if(thepbarRangeCoeffATable)
775        { thepbarRangeCoeffATable->clearAndDestroy();
776          delete thepbarRangeCoeffATable; }
777      thepbarRangeCoeffATable = new G4PhysicsTable(numOfCouples);
778      theRangeCoeffATable = thepbarRangeCoeffATable ;
779      theRangeTable = theRangepbarTable ;
780    }
781 
782  G4double R2 = RTable*RTable ;
783  G4double R1 = RTable+1.;
784  G4double w = R1*(RTable-1.)*(RTable-1.);
785  G4double w1 = RTable/w , w2 = -RTable*R1/w , w3 = R2/w ;
786  G4double Ti , Tim , Tip , Ri , Rim , Rip , Value ;
787  G4bool isOut;
788
789  //  loop for materials
790  for (G4int J=0; J<numOfCouples; J++)
791    {
792      G4int binmax=TotBin ;
793      G4PhysicsLinearVector* aVector = 
794        new G4PhysicsLinearVector(0.,binmax, TotBin);
795      Ti = LowestKineticEnergy ;
796      if (Ti < DBL_MIN) Ti = 1.e-8;
797      G4PhysicsVector* rangeVector= (*theRangeTable)[J];
798
799      for ( G4int i=0; i<TotBin; i++)
800        {
801          Ri = rangeVector->GetValue(Ti,isOut) ;
802          if (Ti < DBL_MIN) Ti = 1.e-8;
803          if ( i==0 )
804            Rim = 0. ;
805          else
806            {
807              // ---- MGP ---- Modified to avoid a floating point exception
808              // The correction is just a temporary patch, the whole class should be redesigned
809              // Original: Tim = Ti/RTable  results in 0./0.
810              if (RTable != 0.)
811                {
812                  Tim = Ti/RTable ;
813                }
814              else
815                {
816                  Tim = 0.;
817                }
818              Rim = rangeVector->GetValue(Tim,isOut);
819            }
820          if ( i==(TotBin-1))
821            Rip = Ri ;
822          else
823            {
824              Tip = Ti*RTable ;
825              Rip = rangeVector->GetValue(Tip,isOut);
826            }
827          Value = (w1*Rip + w2*Ri + w3*Rim)/(Ti*Ti) ; 
828
829          aVector->PutValue(i,Value);
830          Ti = RTable*Ti ;
831        }
832 
833      theRangeCoeffATable->insert(aVector);
834    } 
835}
836
837//--------------------------------------------------------------------------------
838
839void G4hRDEnergyLoss::BuildRangeCoeffBTable(
840                                            const G4ParticleDefinition& )
841// Build tables of coefficients for the energy loss calculation
842//  create table for coefficients "B"
843{
844
845  G4int numOfCouples = G4ProductionCutsTable::GetProductionCutsTable()->GetTableSize();
846
847  if(Charge>0.)
848    {
849      if(thepRangeCoeffBTable)
850        { thepRangeCoeffBTable->clearAndDestroy();
851          delete thepRangeCoeffBTable; }
852      thepRangeCoeffBTable = new G4PhysicsTable(numOfCouples);
853      theRangeCoeffBTable = thepRangeCoeffBTable ;
854      theRangeTable = theRangepTable ;
855    }
856  else
857    {
858      if(thepbarRangeCoeffBTable)
859        { thepbarRangeCoeffBTable->clearAndDestroy();
860          delete thepbarRangeCoeffBTable; }
861      thepbarRangeCoeffBTable = new G4PhysicsTable(numOfCouples);
862      theRangeCoeffBTable = thepbarRangeCoeffBTable ;
863      theRangeTable = theRangepbarTable ;
864    }
865
866  G4double R2 = RTable*RTable ;
867  G4double R1 = RTable+1.;
868  G4double w = R1*(RTable-1.)*(RTable-1.);
869  if (w < DBL_MIN) w = DBL_MIN;
870  G4double w1 = -R1/w , w2 = R1*(R2+1.)/w , w3 = -R2*R1/w ;
871  G4double Ti , Tim , Tip , Ri , Rim , Rip , Value ;
872  G4bool isOut;
873
874  //  loop for materials
875  for (G4int J=0; J<numOfCouples; J++)
876    {
877      G4int binmax=TotBin ;
878      G4PhysicsLinearVector* aVector =
879        new G4PhysicsLinearVector(0.,binmax, TotBin);
880      Ti = LowestKineticEnergy ;
881      if (Ti < DBL_MIN) Ti = 1.e-8;
882      G4PhysicsVector* rangeVector= (*theRangeTable)[J];
883   
884      for ( G4int i=0; i<TotBin; i++)
885        {
886          Ri = rangeVector->GetValue(Ti,isOut) ;
887          if (Ti < DBL_MIN) Ti = 1.e-8;
888          if ( i==0 )
889            Rim = 0. ;
890          else
891            {
892              if (RTable < DBL_MIN) RTable = DBL_MIN;
893              Tim = Ti/RTable ;
894              Rim = rangeVector->GetValue(Tim,isOut);
895            }
896          if ( i==(TotBin-1))
897            Rip = Ri ;
898          else
899            {
900              Tip = Ti*RTable ;
901              Rip = rangeVector->GetValue(Tip,isOut);
902            }
903          if (Ti < DBL_MIN) Ti = DBL_MIN;
904          Value = (w1*Rip + w2*Ri + w3*Rim)/Ti;
905
906          aVector->PutValue(i,Value);
907          Ti = RTable*Ti ;
908        }
909      theRangeCoeffBTable->insert(aVector);
910    } 
911}
912
913//--------------------------------------------------------------------------------
914
915void G4hRDEnergyLoss::BuildRangeCoeffCTable(
916                                            const G4ParticleDefinition& )
917// Build tables of coefficients for the energy loss calculation
918//  create table for coefficients "C"
919{
920
921  G4int numOfCouples = G4ProductionCutsTable::GetProductionCutsTable()->GetTableSize();
922
923  if(Charge>0.)
924    {
925      if(thepRangeCoeffCTable)
926        { thepRangeCoeffCTable->clearAndDestroy();
927          delete thepRangeCoeffCTable; }
928      thepRangeCoeffCTable = new G4PhysicsTable(numOfCouples);
929      theRangeCoeffCTable = thepRangeCoeffCTable ;
930      theRangeTable = theRangepTable ;
931    }
932  else
933    {
934      if(thepbarRangeCoeffCTable)
935        { thepbarRangeCoeffCTable->clearAndDestroy();
936          delete thepbarRangeCoeffCTable; }
937      thepbarRangeCoeffCTable = new G4PhysicsTable(numOfCouples);
938      theRangeCoeffCTable = thepbarRangeCoeffCTable ;
939      theRangeTable = theRangepbarTable ;
940    }
941 
942  G4double R2 = RTable*RTable ;
943  G4double R1 = RTable+1.;
944  G4double w = R1*(RTable-1.)*(RTable-1.);
945  G4double w1 = 1./w , w2 = -RTable*R1/w , w3 = RTable*R2/w ;
946  G4double Ti , Tim , Tip , Ri , Rim , Rip , Value ;
947  G4bool isOut;
948
949  //  loop for materials
950  for (G4int J=0; J<numOfCouples; J++)
951    {
952      G4int binmax=TotBin ;
953      G4PhysicsLinearVector* aVector =
954        new G4PhysicsLinearVector(0.,binmax, TotBin);
955      Ti = LowestKineticEnergy ;
956      G4PhysicsVector* rangeVector= (*theRangeTable)[J];
957   
958      for ( G4int i=0; i<TotBin; i++)
959        {
960          Ri = rangeVector->GetValue(Ti,isOut) ;
961          if ( i==0 )
962            Rim = 0. ;
963          else
964            {
965              Tim = Ti/RTable ;
966              Rim = rangeVector->GetValue(Tim,isOut);
967            }
968          if ( i==(TotBin-1))
969            Rip = Ri ;
970          else
971            {
972              Tip = Ti*RTable ;
973              Rip = rangeVector->GetValue(Tip,isOut);
974            }
975          Value = w1*Rip + w2*Ri + w3*Rim ;
976
977          aVector->PutValue(i,Value);
978          Ti = RTable*Ti ;
979        }
980      theRangeCoeffCTable->insert(aVector);
981    } 
982}
983
984//--------------------------------------------------------------------------------
985
986void G4hRDEnergyLoss::BuildInverseRangeTable(
987                                             const G4ParticleDefinition& aParticleType)
988// Build inverse table of the range table
989{
990  G4bool b;
991
992  const G4ProductionCutsTable* theCoupleTable=
993    G4ProductionCutsTable::GetProductionCutsTable();
994  size_t numOfCouples = theCoupleTable->GetTableSize();
995
996  if(&aParticleType == G4Proton::Proton())
997    {
998      if(theInverseRangepTable)
999        { theInverseRangepTable->clearAndDestroy();
1000          delete theInverseRangepTable; }
1001      theInverseRangepTable = new G4PhysicsTable(numOfCouples);
1002      theInverseRangeTable = theInverseRangepTable ;
1003      theRangeTable = theRangepTable ;
1004      theDEDXTable =  theDEDXpTable ;
1005      theRangeCoeffATable = thepRangeCoeffATable ;
1006      theRangeCoeffBTable = thepRangeCoeffBTable ;
1007      theRangeCoeffCTable = thepRangeCoeffCTable ;
1008    }
1009
1010  if(&aParticleType == G4AntiProton::AntiProton())
1011    {
1012      if(theInverseRangepbarTable)
1013        { theInverseRangepbarTable->clearAndDestroy();
1014          delete theInverseRangepbarTable; }
1015      theInverseRangepbarTable = new G4PhysicsTable(numOfCouples);
1016      theInverseRangeTable = theInverseRangepbarTable ;
1017      theRangeTable = theRangepbarTable ;
1018      theDEDXTable =  theDEDXpbarTable ;
1019      theRangeCoeffATable = thepbarRangeCoeffATable ;
1020      theRangeCoeffBTable = thepbarRangeCoeffBTable ;
1021      theRangeCoeffCTable = thepbarRangeCoeffCTable ;
1022    }
1023
1024  // loop for materials
1025  for (size_t i=0;  i<numOfCouples; i++)
1026    {
1027
1028      G4PhysicsVector* pv = (*theRangeTable)[i];
1029      size_t nbins   = pv->GetVectorLength();
1030      G4double elow  = pv->GetLowEdgeEnergy(0);
1031      G4double ehigh = pv->GetLowEdgeEnergy(nbins-1);
1032      G4double rlow  = pv->GetValue(elow, b);
1033      G4double rhigh = pv->GetValue(ehigh, b);
1034
1035      if (rlow <DBL_MIN) rlow = 1.e-8;
1036      if (rhigh > 1.e16) rhigh = 1.e16;
1037      if (rhigh < 1.e-8) rhigh =1.e-8;
1038      G4double tmpTrick = rhigh/rlow;
1039
1040      //std::cout << nbins << ", elow " << elow << ", ehigh " << ehigh
1041      //                << ", rlow " << rlow << ", rhigh " << rhigh << ", trick " << tmpTrick << std::endl;
1042
1043      if (tmpTrick <= 0. || tmpTrick < DBL_MIN) tmpTrick = 1.e-8; 
1044      if (tmpTrick > 1.e16) tmpTrick = 1.e16; 
1045     
1046      rhigh *= std::exp(std::log(tmpTrick)/((G4double)(nbins-1)));
1047
1048      G4PhysicsLogVector* v = new G4PhysicsLogVector(rlow, rhigh, nbins);
1049
1050      v->PutValue(0,elow);
1051      G4double energy1 = elow;
1052      G4double range1  = rlow;
1053      G4double energy2 = elow;
1054      G4double range2  = rlow;
1055      size_t ilow      = 0;
1056      size_t ihigh;
1057
1058      for (size_t j=1; j<nbins; j++) {
1059
1060        G4double range = v->GetLowEdgeEnergy(j);
1061
1062        for (ihigh=ilow+1; ihigh<nbins; ihigh++) {
1063          energy2 = pv->GetLowEdgeEnergy(ihigh);
1064          range2  = pv->GetValue(energy2, b);
1065          if(range2 >= range || ihigh == nbins-1) {
1066            ilow = ihigh - 1;
1067            energy1 = pv->GetLowEdgeEnergy(ilow);
1068            range1  = pv->GetValue(energy1, b);
1069            break;
1070          }
1071        }
1072
1073        G4double e = std::log(energy1) + std::log(energy2/energy1)*std::log(range/range1)/std::log(range2/range1);
1074
1075        v->PutValue(j,std::exp(e));
1076      }
1077      theInverseRangeTable->insert(v);
1078
1079    }
1080}
1081
1082//--------------------------------------------------------------------------------
1083
1084void G4hRDEnergyLoss::InvertRangeVector(G4int materialIndex,
1085                                        G4PhysicsLogVector* aVector)
1086//  invert range vector for a material
1087{
1088  G4double LowEdgeRange,A,B,C,discr,KineticEnergy ;
1089  G4double Tbin = LowestKineticEnergy/RTable ;
1090  G4double rangebin = 0.0 ;
1091  G4int binnumber = -1 ;
1092  G4bool isOut ;
1093
1094
1095  //loop for range values
1096  for( G4int i=0; i<TotBin; i++)
1097    {
1098      LowEdgeRange = aVector->GetLowEdgeEnergy(i) ;  //i.e. GetLowEdgeValue(i)
1099
1100      if( rangebin < LowEdgeRange )
1101        {
1102          do
1103            {
1104              binnumber += 1 ;
1105              Tbin *= RTable ;
1106              rangebin = (*theRangeTable)(materialIndex)->GetValue(Tbin,isOut) ;
1107            }
1108          while ((rangebin < LowEdgeRange) && (binnumber < TotBin )) ;
1109        }
1110
1111      if(binnumber == 0)
1112        KineticEnergy = LowestKineticEnergy ;
1113      else if(binnumber == TotBin-1)
1114        KineticEnergy = HighestKineticEnergy ;
1115      else
1116        {
1117          A = (*(*theRangeCoeffATable)(materialIndex))(binnumber-1) ;
1118          B = (*(*theRangeCoeffBTable)(materialIndex))(binnumber-1) ;
1119          C = (*(*theRangeCoeffCTable)(materialIndex))(binnumber-1) ;
1120          if(A==0.)
1121            KineticEnergy = (LowEdgeRange -C )/B ;
1122          else
1123            {
1124              discr = B*B - 4.*A*(C-LowEdgeRange);
1125              discr = discr>0. ? std::sqrt(discr) : 0.;
1126              KineticEnergy = 0.5*(discr-B)/A ;
1127            }
1128        }
1129
1130      aVector->PutValue(i,KineticEnergy) ;
1131    }
1132}
1133
1134//------------------------------------------------------------------------------
1135
1136G4bool G4hRDEnergyLoss::CutsWhereModified()
1137{
1138  G4bool wasModified = false;
1139  const G4ProductionCutsTable* theCoupleTable=
1140    G4ProductionCutsTable::GetProductionCutsTable();
1141  size_t numOfCouples = theCoupleTable->GetTableSize();
1142
1143  for (size_t j=0; j<numOfCouples; j++){
1144    if (theCoupleTable->GetMaterialCutsCouple(j)->IsRecalcNeeded()) {
1145      wasModified = true;
1146      break;
1147    }
1148  }
1149  return wasModified;
1150}
1151
1152//-----------------------------------------------------------------------
1153//--- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- --- ---
1154
1155//G4bool G4hRDEnergyLoss::IsApplicable(const G4ParticleDefinition&
1156//                                                     particle)
1157//{
1158//   return((particle.GetPDGCharge()!= 0.) && (particle.GetLeptonNumber() == 0));
1159//}
1160         
1161//G4double G4hRDEnergyLoss::GetContinuousStepLimit(
1162//                                               const G4Track& track,
1163//                                               G4double,
1164//                                               G4double currentMinimumStep,
1165//                                               G4double&)
1166//{
1167//
1168//  G4double Step =
1169//    GetConstraints(track.GetDynamicParticle(),track.GetMaterial()) ;
1170//
1171//  if((Step>0.0)&&(Step<currentMinimumStep))
1172//     currentMinimumStep = Step ;
1173//
1174//  return Step ;
1175//}
Note: See TracBrowser for help on using the repository browser.