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