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

Last change on this file since 1347 was 1347, checked in by garnier, 13 years ago

geant4 tag 9.4
File size: 34.8 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.30 2009/07/23 09:15:37 vnivanch 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 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  //G4cout << "BuildDEDXTable for " << aParticleType.GetParticleName() << G4endl;
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  //G4cout << "BuildDEDXTable done " << G4endl;
347
348  G4EnergyLossTables::Register(&aParticleType,
349    (Charge>0)?
350      theDEDXpTable: theDEDXpbarTable,
351    (Charge>0)?
352      theRangepTable: theRangepbarTable,
353    (Charge>0)?
354      theInverseRangepTable: theInverseRangepbarTable,
355    (Charge>0)?
356      theLabTimepTable: theLabTimepbarTable,
357    (Charge>0)?
358      theProperTimepTable: theProperTimepbarTable,
359    LowestKineticEnergy, HighestKineticEnergy,
360    proton_mass_c2/aParticleType.GetPDGMass(),
361    TotBin);
362
363}
364
365//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
366
367void G4hLowEnergyLoss::BuildRangeTable(
368                             const G4ParticleDefinition& aParticleType)
369// Build range table from the energy loss table
370{
371  //G4cout << "BuildRangeTable for " << aParticleType.GetParticleName() << G4endl;
372   Mass = aParticleType.GetPDGMass();
373
374   const G4ProductionCutsTable* theCoupleTable=
375         G4ProductionCutsTable::GetProductionCutsTable();
376   size_t numOfCouples = theCoupleTable->GetTableSize();
377
378   if( Charge >0.)
379   {
380     if(theRangepTable)
381     { theRangepTable->clearAndDestroy();
382       delete theRangepTable; }
383     theRangepTable = new G4PhysicsTable(numOfCouples);
384     theRangeTable = theRangepTable ;
385   }
386   else
387   {
388     if(theRangepbarTable)
389     { theRangepbarTable->clearAndDestroy();
390       delete theRangepbarTable; }
391     theRangepbarTable = new G4PhysicsTable(numOfCouples);
392     theRangeTable = theRangepbarTable ;
393   }
394
395   // loop for materials
396
397   for (size_t J=0;  J<numOfCouples; J++)
398   {
399     G4PhysicsLogVector* aVector;
400     aVector = new G4PhysicsLogVector(LowestKineticEnergy,
401                              HighestKineticEnergy,TotBin);
402
403     BuildRangeVector(J, aVector);
404     theRangeTable->insert(aVector);
405   }
406}
407
408//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
409
410void G4hLowEnergyLoss::BuildTimeTables(
411                             const G4ParticleDefinition& aParticleType)
412{
413  //G4cout << "BuildTimeTable for " << aParticleType.GetParticleName() << G4endl;
414
415  const G4ProductionCutsTable* theCoupleTable=
416          G4ProductionCutsTable::GetProductionCutsTable();
417  size_t numOfCouples = theCoupleTable->GetTableSize();
418
419  if(&aParticleType == G4Proton::Proton())
420  {
421    if(theLabTimepTable)
422    { theLabTimepTable->clearAndDestroy();
423      delete theLabTimepTable; }
424    theLabTimepTable = new G4PhysicsTable(numOfCouples);
425    theLabTimeTable = theLabTimepTable ;
426
427    if(theProperTimepTable)
428    { theProperTimepTable->clearAndDestroy();
429      delete theProperTimepTable; }
430    theProperTimepTable = new G4PhysicsTable(numOfCouples);
431    theProperTimeTable = theProperTimepTable ;
432  }
433
434  if(&aParticleType == G4AntiProton::AntiProton())
435  {
436    if(theLabTimepbarTable)
437    { theLabTimepbarTable->clearAndDestroy();
438      delete theLabTimepbarTable; }
439    theLabTimepbarTable = new G4PhysicsTable(numOfCouples);
440    theLabTimeTable = theLabTimepbarTable ;
441
442    if(theProperTimepbarTable)
443    { theProperTimepbarTable->clearAndDestroy();
444      delete theProperTimepbarTable; }
445    theProperTimepbarTable = new G4PhysicsTable(numOfCouples);
446    theProperTimeTable = theProperTimepbarTable ;
447  }
448  //G4cout << "numOfCouples= " << numOfCouples << G4endl;
449  for (size_t J=0;  J<numOfCouples; J++)
450  {
451    G4PhysicsLogVector* aVector;
452    G4PhysicsLogVector* bVector;
453
454    aVector = new G4PhysicsLogVector(LowestKineticEnergy,
455                            HighestKineticEnergy,TotBin);
456
457    BuildLabTimeVector(J, aVector);
458    //G4cout << "LabTime OK " << J << G4endl;
459    theLabTimeTable->insert(aVector);
460
461    bVector = new G4PhysicsLogVector(LowestKineticEnergy,
462                            HighestKineticEnergy,TotBin);
463
464    BuildProperTimeVector(J, bVector);
465    //G4cout << "PropTime OK " << J << G4endl;
466    theProperTimeTable->insert(bVector);
467  }
468}
469
470//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
471
472void G4hLowEnergyLoss::BuildRangeVector(G4int materialIndex,
473                                        G4PhysicsLogVector* rangeVector)
474//  create range vector for a material
475{
476  G4bool isOut;
477  G4PhysicsVector* physicsVector= (*theDEDXTable)[materialIndex];
478  G4double energy1 = rangeVector->GetLowEdgeEnergy(0);
479  G4double dedx    = physicsVector->GetValue(energy1,isOut);
480  G4double range   = 0.5*energy1/dedx;
481  rangeVector->PutValue(0,range);
482  G4int n = 100;
483  G4double del = 1.0/(G4double)n ;
484
485  for (G4int j=1; j<=TotBin; j++) {
486
487    G4double energy2 = rangeVector->GetLowEdgeEnergy(j);
488    G4double de = (energy2 - energy1) * del ;
489    G4double dedx1 = dedx ;
490
491    for (G4int i=1; i<n; i++) {
492      G4double energy = energy1 + i*de ;
493      G4double dedx2  = physicsVector->GetValue(energy,isOut);
494      range  += 0.5*de*(1.0/dedx1 + 1.0/dedx2);
495      dedx1   = dedx2;
496    }
497    rangeVector->PutValue(j,range);
498    dedx = dedx1 ;
499    energy1 = energy2 ;
500  }
501}   
502
503//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
504
505void G4hLowEnergyLoss::BuildLabTimeVector(G4int materialIndex,
506                                          G4PhysicsLogVector* timeVector)
507//  create lab time vector for a material
508{
509
510  G4int nbin=100;
511  G4bool isOut;
512  G4double tlim=5.*keV,parlowen=0.4,ppar=0.5-parlowen ;
513  G4double losslim,clim,taulim,timelim,ltaulim,ltaumax,
514           LowEdgeEnergy,tau,Value ;
515
516  G4PhysicsVector* physicsVector= (*theDEDXTable)[materialIndex];
517
518  // low energy part first...
519  losslim = physicsVector->GetValue(tlim,isOut);
520  taulim=tlim/ParticleMass ;
521  clim=std::sqrt(ParticleMass*tlim/2.)/(c_light*losslim*ppar) ;
522  ltaulim = std::log(taulim);
523  ltaumax = std::log(HighestKineticEnergy/ParticleMass) ;
524
525  G4int i=-1;
526  G4double oldValue = 0. ;
527  G4double tauold ;
528  do 
529  {
530    i += 1 ;
531    LowEdgeEnergy = timeVector->GetLowEdgeEnergy(i);
532    tau = LowEdgeEnergy/ParticleMass ;
533    if ( tau <= taulim )
534    {
535      Value = clim*std::exp(ppar*std::log(tau/taulim)) ;
536    }
537    else
538    {
539      timelim=clim ;
540      ltaulow = std::log(taulim);
541      ltauhigh = std::log(tau);
542      Value = timelim+LabTimeIntLog(physicsVector,nbin);
543    }
544    timeVector->PutValue(i,Value);
545    oldValue = Value ;
546    tauold = tau ;
547  } while (tau<=taulim) ;
548
549  i += 1 ;
550  //G4cout << "do is OK i= " << i << G4endl;
551  for (G4int j=i; j<=TotBin; j++)
552  {
553    LowEdgeEnergy = timeVector->GetLowEdgeEnergy(j);
554    tau = LowEdgeEnergy/ParticleMass ;
555    //G4cout << "j= " << j << " tauold= " << tauold << " tau= " << tau << G4endl;
556    ltaulow = std::log(tauold);
557    ltauhigh = std::log(tau);
558    Value = oldValue+LabTimeIntLog(physicsVector,nbin);
559    timeVector->PutValue(j,Value);
560    oldValue = Value ;
561    tauold = tau ;
562  }
563  // G4cout << "LabTime OK for  " << materialIndex << G4endl;
564}
565
566//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
567
568void G4hLowEnergyLoss::BuildProperTimeVector(G4int materialIndex,
569                                             G4PhysicsLogVector* timeVector)
570//  create proper time vector for a material
571{ 
572  G4int nbin=100;
573  G4bool isOut;
574  G4double tlim=5.*keV,parlowen=0.4,ppar=0.5-parlowen ;
575  G4double losslim,clim,taulim,timelim,ltaulim,ltaumax,
576           LowEdgeEnergy,tau,Value ;
577
578  G4PhysicsVector* physicsVector= (*theDEDXTable)[materialIndex];
579 
580  // low energy part first...
581  losslim = physicsVector->GetValue(tlim,isOut);
582  taulim=tlim/ParticleMass ;
583  clim=std::sqrt(ParticleMass*tlim/2.)/(c_light*losslim*ppar) ;
584  ltaulim = std::log(taulim);
585  ltaumax = std::log(HighestKineticEnergy/ParticleMass) ;
586
587  G4int i=-1;
588  G4double oldValue = 0. ;
589  G4double tauold ;
590  do
591  {
592    i += 1 ;
593    LowEdgeEnergy = timeVector->GetLowEdgeEnergy(i);
594    tau = LowEdgeEnergy/ParticleMass ;
595    if ( tau <= taulim )
596    {
597      Value = clim*std::exp(ppar*std::log(tau/taulim)) ;
598    }
599    else
600    {
601      timelim=clim ;
602      ltaulow = std::log(taulim);
603      ltauhigh = std::log(tau);
604      Value = timelim+ProperTimeIntLog(physicsVector,nbin);
605    }
606    timeVector->PutValue(i,Value);
607    oldValue = Value ;
608    tauold = tau ;
609  } while (tau<=taulim) ;
610
611  i += 1 ;
612  for (G4int j=i; j<=TotBin; j++)
613  {
614    LowEdgeEnergy = timeVector->GetLowEdgeEnergy(j);
615    tau = LowEdgeEnergy/ParticleMass ;
616    ltaulow = std::log(tauold);
617    ltauhigh = std::log(tau);
618    Value = oldValue+ProperTimeIntLog(physicsVector,nbin);
619    timeVector->PutValue(j,Value);
620    oldValue = Value ;
621    tauold = tau ;
622  }
623}
624
625//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
626
627G4double G4hLowEnergyLoss::RangeIntLin(G4PhysicsVector* physicsVector,
628                                       G4int nbin)
629//  num. integration, linear binning
630{
631  G4double dtau,Value,taui,ti,lossi,ci;
632  G4bool isOut;
633  dtau = (tauhigh-taulow)/nbin;
634  Value = 0.;
635
636  for (G4int i=0; i<nbin; i++)
637  {
638    taui = taulow + dtau*i ;
639    ti = Mass*taui;
640    lossi = physicsVector->GetValue(ti,isOut);
641    if(i==0)
642      ci=0.5;
643    else
644    {
645      if(i<nbin-1)
646        ci=1.;
647      else
648        ci=0.5;
649    }
650    Value += ci/lossi;
651  }
652  Value *= Mass*dtau;
653  return Value;
654}
655
656//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
657
658G4double G4hLowEnergyLoss::RangeIntLog(G4PhysicsVector* physicsVector,
659                                       G4int nbin)
660//  num. integration, logarithmic binning
661{
662  G4double ltt,dltau,Value,ui,taui,ti,lossi,ci;
663  G4bool isOut;
664  ltt = ltauhigh-ltaulow;
665  dltau = ltt/nbin;
666  Value = 0.;
667
668  for (G4int i=0; i<nbin; i++)
669  {
670    ui = ltaulow+dltau*i;
671    taui = std::exp(ui);
672    ti = Mass*taui;
673    lossi = physicsVector->GetValue(ti,isOut);
674    if(i==0)
675      ci=0.5;
676    else
677    {
678      if(i<nbin-1)
679        ci=1.;
680      else
681        ci=0.5;
682    }
683    Value += ci*taui/lossi;
684  }
685  Value *= Mass*dltau;
686  return Value;
687}
688
689//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
690
691G4double G4hLowEnergyLoss::LabTimeIntLog(G4PhysicsVector* physicsVector,
692                                         G4int nbin)
693//  num. integration, logarithmic binning
694{
695  G4double ltt,dltau,Value,ui,taui,ti,lossi,ci;
696  G4bool isOut;
697  ltt = ltauhigh-ltaulow;
698  dltau = ltt/nbin;
699  Value = 0.;
700
701  for (G4int i=0; i<nbin; i++)
702  {
703    ui = ltaulow+dltau*i;
704    taui = std::exp(ui);
705    ti = ParticleMass*taui;
706    lossi = physicsVector->GetValue(ti,isOut);
707    if(i==0)
708      ci=0.5;
709    else
710    {
711      if(i<nbin-1)
712        ci=1.;
713      else
714        ci=0.5;
715    }
716    Value += ci*taui*(ti+ParticleMass)/(std::sqrt(ti*(ti+2.*ParticleMass))*lossi);
717  }
718  Value *= ParticleMass*dltau/c_light;
719  return Value;
720}
721
722//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
723
724G4double G4hLowEnergyLoss::ProperTimeIntLog(G4PhysicsVector* physicsVector,
725                                            G4int nbin)
726//  num. integration, logarithmic binning
727{
728  G4double ltt,dltau,Value,ui,taui,ti,lossi,ci;
729  G4bool isOut;
730  ltt = ltauhigh-ltaulow;
731  dltau = ltt/nbin;
732  Value = 0.;
733
734  for (G4int i=0; i<nbin; i++)
735  {
736    ui = ltaulow+dltau*i;
737    taui = std::exp(ui);
738    ti = ParticleMass*taui;
739    lossi = physicsVector->GetValue(ti,isOut);
740    if(i==0)
741      ci=0.5;
742    else
743    {
744      if(i<nbin-1)
745        ci=1.;
746      else
747        ci=0.5;
748    }
749    Value += ci*taui*ParticleMass/(std::sqrt(ti*(ti+2.*ParticleMass))*lossi);
750  }
751  Value *= ParticleMass*dltau/c_light;
752  return Value;
753}
754
755//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
756
757void G4hLowEnergyLoss::BuildRangeCoeffATable(
758                                             const G4ParticleDefinition& )
759// Build tables of coefficients for the energy loss calculation
760//  create table for coefficients "A"
761{
762  //G4cout << "BuildRangeCoeffATable for " << G4endl;
763
764  G4int numOfCouples = G4ProductionCutsTable::GetProductionCutsTable()->GetTableSize();
765
766  if(Charge>0.)
767  {
768    if(thepRangeCoeffATable)
769    { thepRangeCoeffATable->clearAndDestroy();
770      delete thepRangeCoeffATable; }
771    thepRangeCoeffATable = new G4PhysicsTable(numOfCouples);
772    theRangeCoeffATable = thepRangeCoeffATable ;
773    theRangeTable = theRangepTable ;
774  }
775  else
776  {
777    if(thepbarRangeCoeffATable)
778    { thepbarRangeCoeffATable->clearAndDestroy();
779      delete thepbarRangeCoeffATable; }
780    thepbarRangeCoeffATable = new G4PhysicsTable(numOfCouples);
781    theRangeCoeffATable = thepbarRangeCoeffATable ;
782    theRangeTable = theRangepbarTable ;
783  }
784  G4double R2 = RTable*RTable ;
785  G4double R1 = RTable+1.;
786  G4double w = R1*(RTable-1.)*(RTable-1.);
787  G4double w1 = RTable/w , w2 = -RTable*R1/w , w3 = R2/w ;
788  G4double Ti , Tim , Tip , Ri , Rim , Rip , Value ;
789  G4bool isOut;
790
791  //  loop for materials
792  for (G4int J=0; J<numOfCouples; J++)
793  {
794    G4int binmax=TotBin ;
795    G4PhysicsLinearVector* aVector = 
796                           new G4PhysicsLinearVector(0.,binmax, TotBin);
797    Ti = LowestKineticEnergy ;
798    G4PhysicsVector* rangeVector= (*theRangeTable)[J];
799
800    for ( G4int i=0; i<=TotBin; i++)
801    {
802      Ri = rangeVector->GetValue(Ti,isOut) ;
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)
821        Rip = Ri ;
822      else
823      {
824        Tip = Ti*RTable ;
825        Rip = rangeVector->GetValue(Tip,isOut);
826      }
827      if (Ti!=0) Value = (w1*Rip + w2*Ri + w3*Rim)/(Ti*Ti); else Value=0; 
828
829      aVector->PutValue(i,Value);
830      Ti = RTable*Ti ;
831    }
832 
833    theRangeCoeffATable->insert(aVector);
834  } 
835}
836
837//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
838
839void G4hLowEnergyLoss::BuildRangeCoeffBTable(
840                                             const G4ParticleDefinition& )
841// Build tables of coefficients for the energy loss calculation
842//  create table for coefficients "B"
843{
844  //G4cout << "BuildRangeCoeffBTable for " << G4endl;
845
846  G4int numOfCouples = G4ProductionCutsTable::GetProductionCutsTable()->GetTableSize();
847
848  if(Charge>0.)
849  {
850    if(thepRangeCoeffBTable)
851    { thepRangeCoeffBTable->clearAndDestroy();
852      delete thepRangeCoeffBTable; }
853    thepRangeCoeffBTable = new G4PhysicsTable(numOfCouples);
854    theRangeCoeffBTable = thepRangeCoeffBTable ;
855    theRangeTable = theRangepTable ;
856  }
857  else
858  {
859    if(thepbarRangeCoeffBTable)
860    { thepbarRangeCoeffBTable->clearAndDestroy();
861      delete thepbarRangeCoeffBTable; }
862    thepbarRangeCoeffBTable = new G4PhysicsTable(numOfCouples);
863    theRangeCoeffBTable = thepbarRangeCoeffBTable ;
864    theRangeTable = theRangepbarTable ;
865  }
866
867  G4double R2 = RTable*RTable ;
868  G4double R1 = RTable+1.;
869  G4double w = R1*(RTable-1.)*(RTable-1.);
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    G4PhysicsVector* rangeVector= (*theRangeTable)[J];
882   
883    for ( G4int i=0; i<=TotBin; i++)
884    {
885      Ri = rangeVector->GetValue(Ti,isOut) ;
886      if ( i==0 )
887         Rim = 0. ;
888      else
889      {
890        if (RTable!=0) Tim = Ti/RTable ; else Tim =0;
891        Rim = rangeVector->GetValue(Tim,isOut);
892      }
893      if ( i==TotBin)
894        Rip = Ri ;
895      else
896      {
897        Tip = Ti*RTable ;
898        Rip = rangeVector->GetValue(Tip,isOut);
899      }
900      if (Ti!=0) Value = (w1*Rip + w2*Ri + w3*Rim)/Ti; else Value=0;
901
902      aVector->PutValue(i,Value);
903      Ti = RTable*Ti ;
904    }
905    theRangeCoeffBTable->insert(aVector);
906  } 
907}
908
909//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
910
911void G4hLowEnergyLoss::BuildRangeCoeffCTable(
912                                             const G4ParticleDefinition& )
913// Build tables of coefficients for the energy loss calculation
914//  create table for coefficients "C"
915{
916  //G4cout << "BuildRangeCoeffCTable for " << G4endl;
917
918  G4int numOfCouples = G4ProductionCutsTable::GetProductionCutsTable()->GetTableSize();
919
920  if(Charge>0.)
921  {
922    if(thepRangeCoeffCTable)
923    { thepRangeCoeffCTable->clearAndDestroy();
924      delete thepRangeCoeffCTable; }
925    thepRangeCoeffCTable = new G4PhysicsTable(numOfCouples);
926    theRangeCoeffCTable = thepRangeCoeffCTable ;
927    theRangeTable = theRangepTable ;
928  }
929  else
930  {
931    if(thepbarRangeCoeffCTable)
932    { thepbarRangeCoeffCTable->clearAndDestroy();
933      delete thepbarRangeCoeffCTable; }
934    thepbarRangeCoeffCTable = new G4PhysicsTable(numOfCouples);
935    theRangeCoeffCTable = thepbarRangeCoeffCTable ;
936    theRangeTable = theRangepbarTable ;
937  }
938 
939  G4double R2 = RTable*RTable ;
940  G4double R1 = RTable+1.;
941  G4double w = R1*(RTable-1.)*(RTable-1.);
942  G4double w1 = 1./w , w2 = -RTable*R1/w , w3 = RTable*R2/w ;
943  G4double Ti , Tim , Tip , Ri , Rim , Rip , Value ;
944  G4bool isOut;
945
946  //  loop for materials
947  for (G4int J=0; J<numOfCouples; J++)
948  {
949    G4int binmax=TotBin ;
950    G4PhysicsLinearVector* aVector =
951                      new G4PhysicsLinearVector(0.,binmax, TotBin);
952    Ti = LowestKineticEnergy ;
953    G4PhysicsVector* rangeVector= (*theRangeTable)[J];
954   
955    for ( G4int i=0; i<=TotBin; i++)
956    {
957      Ri = rangeVector->GetValue(Ti,isOut) ;
958      if ( i==0 )
959        Rim = 0. ;
960      else
961      {
962        if (RTable!=0) Tim = Ti/RTable ; else Tim=0;
963        Rim = rangeVector->GetValue(Tim,isOut);
964      }
965      if ( i==TotBin)
966        Rip = Ri ;
967      else
968      {
969        Tip = Ti*RTable ;
970        Rip = rangeVector->GetValue(Tip,isOut);
971      }
972      Value = w1*Rip + w2*Ri + w3*Rim ;
973
974      aVector->PutValue(i,Value);
975      Ti = RTable*Ti ;
976    }
977    theRangeCoeffCTable->insert(aVector);
978  } 
979}
980
981//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
982
983void G4hLowEnergyLoss::BuildInverseRangeTable(
984                             const G4ParticleDefinition& aParticleType)
985// Build inverse table of the range table
986{
987  //G4cout << "BuildInverseRangeTable for " << aParticleType.GetParticleName() << G4endl;
988  G4bool b;
989
990  const G4ProductionCutsTable* theCoupleTable=
991        G4ProductionCutsTable::GetProductionCutsTable();
992  size_t numOfCouples = theCoupleTable->GetTableSize();
993
994  if(&aParticleType == G4Proton::Proton())
995  {
996    if(theInverseRangepTable)
997    { theInverseRangepTable->clearAndDestroy();
998      delete theInverseRangepTable; }
999    theInverseRangepTable = new G4PhysicsTable(numOfCouples);
1000    theInverseRangeTable = theInverseRangepTable ;
1001    theRangeTable = theRangepTable ;
1002    theDEDXTable =  theDEDXpTable ;
1003    theRangeCoeffATable = thepRangeCoeffATable ;
1004    theRangeCoeffBTable = thepRangeCoeffBTable ;
1005    theRangeCoeffCTable = thepRangeCoeffCTable ;
1006  }
1007
1008  if(&aParticleType == G4AntiProton::AntiProton())
1009  {
1010    if(theInverseRangepbarTable)
1011    { theInverseRangepbarTable->clearAndDestroy();
1012      delete theInverseRangepbarTable; }
1013    theInverseRangepbarTable = new G4PhysicsTable(numOfCouples);
1014    theInverseRangeTable = theInverseRangepbarTable ;
1015    theRangeTable = theRangepbarTable ;
1016    theDEDXTable =  theDEDXpbarTable ;
1017    theRangeCoeffATable = thepbarRangeCoeffATable ;
1018    theRangeCoeffBTable = thepbarRangeCoeffBTable ;
1019    theRangeCoeffCTable = thepbarRangeCoeffCTable ;
1020  }
1021
1022  // loop for materials
1023  for (size_t i=0;  i<numOfCouples; i++)
1024  {
1025
1026    G4PhysicsVector* pv = (*theRangeTable)[i];
1027    size_t nbins   = pv->GetVectorLength();
1028    G4double elow  = pv->GetLowEdgeEnergy(0);
1029    G4double ehigh = pv->GetLowEdgeEnergy(nbins-1);
1030    G4double rlow  = pv->GetValue(elow, b);
1031    G4double rhigh = pv->GetValue(ehigh, b);
1032    //G4cout << "elow= " << elow << " ehigh= " << ehigh << " rlow= " << rlow << " rhigh= " << rhigh << G4endl;
1033    //    rhigh *= std::exp(std::log(rhigh/rlow)/((G4double)(nbins-1)));
1034
1035    G4PhysicsLogVector* v = new G4PhysicsLogVector(rlow, rhigh, nbins-1);
1036
1037    v->PutValue(0,elow);
1038    G4double energy1 = elow;
1039    G4double range1  = rlow;
1040    G4double energy2 = elow;
1041    G4double range2  = rlow;
1042    size_t ilow      = 0;
1043    size_t ihigh;
1044
1045    for (size_t j=1; j<nbins; j++) {
1046
1047      G4double range = v->GetLowEdgeEnergy(j);
1048
1049      for (ihigh=ilow+1; ihigh<nbins; ihigh++) {
1050        energy2 = pv->GetLowEdgeEnergy(ihigh);
1051        range2  = pv->GetValue(energy2, b);
1052        if(range2 >= range || ihigh == nbins-1) {
1053          ilow = ihigh - 1;
1054          energy1 = pv->GetLowEdgeEnergy(ilow);
1055          range1  = pv->GetValue(energy1, b);
1056          break;
1057        }
1058      }
1059
1060      G4double e = std::log(energy1) + std::log(energy2/energy1)*std::log(range/range1)/std::log(range2/range1);
1061
1062      v->PutValue(j,std::exp(e));
1063    }
1064    theInverseRangeTable->insert(v);
1065
1066  }
1067}
1068
1069//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1070
1071void G4hLowEnergyLoss::InvertRangeVector(G4int materialIndex,
1072                                         G4PhysicsLogVector* aVector)
1073//  invert range vector for a material
1074{
1075  G4double LowEdgeRange,A,B,C,discr,KineticEnergy ;
1076  G4double Tbin = 0;
1077  if (RTable !=0.) Tbin = LowestKineticEnergy/RTable ;
1078  G4double rangebin = 0.0 ;
1079  G4int binnumber = -1 ;
1080  G4bool isOut ;
1081
1082
1083  //loop for range values
1084  for( G4int i=0; i<=TotBin; i++)
1085  {
1086    LowEdgeRange = aVector->GetLowEdgeEnergy(i) ;  //i.e. GetLowEdgeValue(i)
1087
1088    if( rangebin < LowEdgeRange )
1089    {
1090      do
1091      {
1092        binnumber += 1 ;
1093        Tbin *= RTable ;
1094        rangebin = (*theRangeTable)(materialIndex)->GetValue(Tbin,isOut) ;
1095      }
1096      while ((rangebin < LowEdgeRange) && (binnumber < TotBin )) ;
1097    }
1098
1099    if(binnumber == 0)
1100      KineticEnergy = LowestKineticEnergy ;
1101    else if(binnumber == TotBin)
1102      KineticEnergy = HighestKineticEnergy ;
1103    else
1104    {
1105      A = (*(*theRangeCoeffATable)(materialIndex))(binnumber-1) ;
1106      B = (*(*theRangeCoeffBTable)(materialIndex))(binnumber-1) ;
1107      C = (*(*theRangeCoeffCTable)(materialIndex))(binnumber-1) ;
1108      if(A==0.)
1109         KineticEnergy = (LowEdgeRange -C )/B ;
1110      else
1111      {
1112         discr = B*B - 4.*A*(C-LowEdgeRange);
1113         discr = discr>0. ? std::sqrt(discr) : 0.;
1114         KineticEnergy = 0.5*(discr-B)/A ;
1115      }
1116    }
1117
1118    aVector->PutValue(i,KineticEnergy) ;
1119  }
1120}
1121
1122//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1123
1124G4bool G4hLowEnergyLoss::CutsWhereModified()
1125{
1126  G4bool wasModified = false;
1127  const G4ProductionCutsTable* theCoupleTable=
1128        G4ProductionCutsTable::GetProductionCutsTable();
1129  size_t numOfCouples = theCoupleTable->GetTableSize();
1130
1131  for (size_t j=0; j<numOfCouples; j++){
1132    if (theCoupleTable->GetMaterialCutsCouple(j)->IsRecalcNeeded()) {
1133      wasModified = true;
1134      break;
1135    }
1136  }
1137  return wasModified;
1138}
1139
1140//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
1141
1142
Note: See TracBrowser for help on using the repository browser.