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