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

Last change on this file since 1201 was 1196, checked in by garnier, 16 years ago

update CVS release candidate geant4.9.3.01
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-03-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 //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.