ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4hRDEnergyLoss.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4hRDEnergyLoss.cc
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.
87 
90 
93 
96 
107 
114 
120 
124 
125 //const G4Proton* G4hRDEnergyLoss::theProton=G4Proton::Proton() ;
126 //const G4AntiProton* G4hRDEnergyLoss::theAntiProton=G4AntiProton::AntiProton() ;
127 
131 
137 
139 G4ThreadLocal G4double G4hRDEnergyLoss::finalRange = 0.2; // 200.*micrometer
140 
143  // 2.*(1.-(0.20))*(200.*micrometer)
145  // -(1.-(0.20))*(200.*micrometer)*(200.*micrometer)
146 
148 
151 
156 
157 //--------------------------------------------------------------------------------
158 
159 // constructor and destructor
160 
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 {
175 }
176 
177 //--------------------------------------------------------------------------------
178 
180 {
181  if(theLossTable)
182  {
184  delete theLossTable;
185  }
186 }
187 
188 //--------------------------------------------------------------------------------
189 
191 {
192  return NumberOfProcesses;
193 }
194 
195 //--------------------------------------------------------------------------------
196 
198 {
199  NumberOfProcesses=number;
200 }
201 
202 //--------------------------------------------------------------------------------
203 
205 {
207 }
208 
209 //--------------------------------------------------------------------------------
210 
212 {
214 }
215 
216 //--------------------------------------------------------------------------------
217 
219 {
220  dRoverRange = value;
221 }
222 
223 //--------------------------------------------------------------------------------
224 
226 {
228 }
229 
230 //--------------------------------------------------------------------------------
231 
233 {
235 }
236 
237 //--------------------------------------------------------------------------------
238 
240 {
241  dRoverRange = c1;
242  finalRange = c2;
246 }
247 
248 //--------------------------------------------------------------------------------
249 
251 {
252  // calculate data members TotBin,LOGRTable,RTable first
253 
257 
258  const G4ProductionCutsTable* theCoupleTable=
260  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;}
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  {
280 
282  {
283  if(theDEDXpTable)
285  delete theDEDXpTable; }
286  theDEDXpTable = new G4PhysicsTable(numOfCouples);
288  }
289  }
290  else
291  {
294 
296  {
297  if(theDEDXpbarTable)
299  delete theDEDXpbarTable; }
300  theDEDXpbarTable = new G4PhysicsTable(numOfCouples);
302  }
303  }
304 
306  {
307  // loop for materials
308  G4double LowEdgeEnergy , Value ;
309  G4bool isOutRange ;
310  G4PhysicsTable* pointer ;
311 
312  for (size_t J=0; J<numOfCouples; J++)
313  {
314  // create physics vector and fill it
315  G4PhysicsLogVector* aVector =
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.)
342  else
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)?
366  (Charge>0)?
368  (Charge>0)?
370  (Charge>0)?
372  (Charge>0)?
375  proton_mass_c2/aParticleType.GetPDGMass(),
376  TotBin);
377 
378 }
379 
380 //--------------------------------------------------------------------------------
381 
383 {
384  // Build range table from the energy loss table
385 
386  Mass = aParticleType.GetPDGMass();
387 
388  const G4ProductionCutsTable* theCoupleTable=
390  size_t numOfCouples = theCoupleTable->GetTableSize();
391 
392  if( Charge >0.)
393  {
394  if(theRangepTable)
396  delete theRangepTable; }
397  theRangepTable = new G4PhysicsTable(numOfCouples);
399  }
400  else
401  {
404  delete theRangepbarTable; }
405  theRangepbarTable = new G4PhysicsTable(numOfCouples);
407  }
408 
409  // loop for materials
410 
411  for (size_t J=0; J<numOfCouples; J++)
412  {
413  G4PhysicsLogVector* aVector;
416 
417  BuildRangeVector(J, aVector);
418  theRangeTable->insert(aVector);
419  }
420 }
421 
422 //--------------------------------------------------------------------------------
423 
425 {
426  const G4ProductionCutsTable* theCoupleTable=
428  size_t numOfCouples = theCoupleTable->GetTableSize();
429 
430  if(&aParticleType == G4Proton::Proton())
431  {
432  if(theLabTimepTable)
434  delete theLabTimepTable; }
435  theLabTimepTable = new G4PhysicsTable(numOfCouples);
437 
440  delete theProperTimepTable; }
441  theProperTimepTable = new G4PhysicsTable(numOfCouples);
443  }
444 
445  if(&aParticleType == G4AntiProton::AntiProton())
446  {
449  delete theLabTimepbarTable; }
450  theLabTimepbarTable = new G4PhysicsTable(numOfCouples);
452 
455  delete theProperTimepbarTable; }
456  theProperTimepbarTable = new G4PhysicsTable(numOfCouples);
458  }
459 
460  for (size_t J=0; J<numOfCouples; J++)
461  {
462  G4PhysicsLogVector* aVector;
463  G4PhysicsLogVector* bVector;
464 
467 
468  BuildLabTimeVector(J, aVector);
469  theLabTimeTable->insert(aVector);
470 
473 
474  BuildProperTimeVector(J, bVector);
475  theProperTimeTable->insert(bVector);
476  }
477 }
478 
479 //--------------------------------------------------------------------------------
480 
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 
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 
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 
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 
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 
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 
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 
772 {
773  // Build tables of coefficients for the energy loss calculation
774  // create table for coefficients "A"
775 
777 
778  if(Charge>0.)
779  {
782  delete thepRangeCoeffATable; }
783  thepRangeCoeffATable = new G4PhysicsTable(numOfCouples);
786  }
787  else
788  {
791  delete thepbarRangeCoeffATable; }
792  thepbarRangeCoeffATable = new G4PhysicsTable(numOfCouples);
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 
855 {
856  // Build tables of coefficients for the energy loss calculation
857  // create table for coefficients "B"
858 
859  G4int numOfCouples =
861 
862  if(Charge>0.)
863  {
866  delete thepRangeCoeffBTable; }
867  thepRangeCoeffBTable = new G4PhysicsTable(numOfCouples);
870  }
871  else
872  {
875  delete thepbarRangeCoeffBTable; }
876  thepbarRangeCoeffBTable = new G4PhysicsTable(numOfCouples);
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 
931 {
932  // Build tables of coefficients for the energy loss calculation
933  // create table for coefficients "C"
934 
935  G4int numOfCouples =
937 
938  if(Charge>0.)
939  {
942  delete thepRangeCoeffCTable; }
943  thepRangeCoeffCTable = new G4PhysicsTable(numOfCouples);
946  }
947  else
948  {
951  delete thepbarRangeCoeffCTable; }
952  thepbarRangeCoeffCTable = new G4PhysicsTable(numOfCouples);
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::
1003 {
1004  // Build inverse table of the range table
1005 
1006  G4bool b;
1007 
1008  const G4ProductionCutsTable* theCoupleTable=
1010  size_t numOfCouples = theCoupleTable->GetTableSize();
1011 
1012  if(&aParticleType == G4Proton::Proton())
1013  {
1016  delete theInverseRangepTable; }
1017  theInverseRangepTable = new G4PhysicsTable(numOfCouples);
1024  }
1025 
1026  if(&aParticleType == G4AntiProton::AntiProton())
1027  {
1030  delete theInverseRangepbarTable; }
1031  theInverseRangepbarTable = new G4PhysicsTable(numOfCouples);
1038  }
1039 
1040  // loop for materials
1041  for (size_t i=0; i<numOfCouples; i++)
1042  {
1043 
1044  G4PhysicsVector* pv = (*theRangeTable)[i];
1045  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  size_t ilow = 0;
1073  size_t ihigh;
1074 
1075  for (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  }
1095 
1096  }
1097 }
1098 
1099 //--------------------------------------------------------------------------------
1100 
1102  G4PhysicsLogVector* aVector)
1103 {
1104  // invert range vector for a material
1105 
1106  G4double LowEdgeRange,A,B,C,discr,KineticEnergy ;
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 
1155 {
1156  G4bool wasModified = false;
1157  const G4ProductionCutsTable* theCoupleTable=
1159  size_t numOfCouples = theCoupleTable->GetTableSize();
1160 
1161  for (size_t 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 //}