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

Last change on this file since 830 was 819, checked in by garnier, 17 years ago

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