source: trunk/source/processes/electromagnetic/lowenergy/src/G4hLowEnergyLoss.cc @ 1066

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

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