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

Last change on this file since 1058 was 1055, checked in by garnier, 17 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.