Geant4  10.01
G4hImpactIonisation.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 
28 // ------------------------------------------------------------
29 // G4RDHadronIonisation
30 //
31 // $Id: G4hImpactIonisation.cc 70904 2013-06-07 10:34:25Z gcosmo $
32 //
33 // Author: Maria Grazia Pia (MariaGrazia.Pia@ge.infn.it)
34 //
35 // 08 Sep 2008 - MGP - Created (initially based on G4hLowEnergyIonisation)
36 // Added PIXE capabilities
37 // Partial clean-up of the implementation (more needed)
38 // Calculation of MicroscopicCrossSection delegated to specialised cla// Documentation available in:
39 // M.G. Pia et al., PIXE Simulation With Geant4,
40 // IEEE Trans. Nucl. Sci., vol. 56, no. 6, pp. 3614-3649, Dec. 2009.
41 
42 //
43 // ------------------------------------------------------------
44 
45 #include "G4hImpactIonisation.hh"
46 #include "globals.hh"
47 #include "G4ios.hh"
48 #include "Randomize.hh"
49 #include "G4PhysicalConstants.hh"
50 #include "G4SystemOfUnits.hh"
51 #include "G4Poisson.hh"
52 #include "G4UnitsTable.hh"
53 #include "G4EnergyLossTables.hh"
54 #include "G4Material.hh"
55 #include "G4DynamicParticle.hh"
56 #include "G4ParticleDefinition.hh"
57 #include "G4AtomicDeexcitation.hh"
60 #include "G4IInterpolator.hh"
61 #include "G4LogLogInterpolator.hh"
62 #include "G4Gamma.hh"
63 #include "G4Electron.hh"
64 #include "G4Proton.hh"
65 #include "G4ProcessManager.hh"
66 #include "G4ProductionCutsTable.hh"
67 #include "G4PhysicsLogVector.hh"
68 #include "G4PhysicsLinearVector.hh"
69 
70 #include "G4VLowEnergyModel.hh"
72 #include "G4hBetheBlochModel.hh"
74 #include "G4QAOLowEnergyLoss.hh"
75 #include "G4hIonEffChargeSquare.hh"
78 
79 #include "G4MaterialCutsCouple.hh"
80 #include "G4Track.hh"
81 #include "G4Step.hh"
82 
84  : G4hRDEnergyLoss(processName),
85  betheBlochModel(0),
86  protonModel(0),
87  antiprotonModel(0),
88  theIonEffChargeModel(0),
89  theNuclearStoppingModel(0),
90  theIonChuFluctuationModel(0),
91  theIonYangFluctuationModel(0),
92  protonTable("ICRU_R49p"),
93  antiprotonTable("ICRU_R49p"),
94  theNuclearTable("ICRU_R49"),
95  nStopping(true),
96  theBarkas(true),
97  theMeanFreePathTable(0),
98  paramStepLimit (0.005),
99  pixeCrossSectionHandler(0)
100 {
101  InitializeMe();
102 }
103 
104 
105 
107 {
108  LowestKineticEnergy = 10.0*eV ;
109  HighestKineticEnergy = 100.0*GeV ;
110  MinKineticEnergy = 10.0*eV ;
111  TotBin = 360 ;
112  protonLowEnergy = 1.*keV ;
113  protonHighEnergy = 100.*MeV ;
114  antiprotonLowEnergy = 25.*keV ;
116  minGammaEnergy = 100 * eV;
117  minElectronEnergy = 250.* eV;
118  verboseLevel = 0;
119 
120  // Min and max energy of incident particle for the calculation of shell cross sections
121  // for PIXE generation
122  eMinPixe = 1.* keV;
123  eMaxPixe = 200. * MeV;
124 
125  G4String defaultPixeModel("ecpssr");
126  modelK = defaultPixeModel;
127  modelL = defaultPixeModel;
128  modelM = defaultPixeModel;
129 }
130 
131 
133 {
135  {
137  delete theMeanFreePathTable;
138  }
139 
140  if (betheBlochModel) delete betheBlochModel;
141  if (protonModel) delete protonModel;
142  if (antiprotonModel) delete antiprotonModel;
147 
149 
150  // ---- MGP ---- The following is to be checked
151  // if (shellVacancy) delete shellVacancy;
152 
153  cutForDelta.clear();
154 }
155 
156 // --------------------------------------------------------------------
158  const G4String& dedxTable)
159 // This method defines the ionisation parametrisation method via its name
160 {
161  if (particle->GetPDGCharge() > 0 )
162  {
163  // Positive charge
165  }
166  else
167  {
168  // Antiprotons
170  }
171 }
172 
173 
174 // --------------------------------------------------------------------
176 
177 {
178  // Define models for parametrisation of electronic energy losses
179  betheBlochModel = new G4hBetheBlochModel("Bethe-Bloch") ;
184  theIonEffChargeModel = new G4hIonEffChargeSquare("Ziegler1988") ;
187 }
188 
189 
190 // --------------------------------------------------------------------
192 
193 // just call BuildLossTable+BuildLambdaTable
194 {
195 
196  // Verbose print-out
197  if(verboseLevel > 0)
198  {
199  G4cout << "G4hImpactIonisation::BuildPhysicsTable for "
200  << particleDef.GetParticleName()
201  << " mass(MeV)= " << particleDef.GetPDGMass()/MeV
202  << " charge= " << particleDef.GetPDGCharge()/eplus
203  << " type= " << particleDef.GetParticleType()
204  << G4endl;
205 
206  if(verboseLevel > 1)
207  {
208  G4ProcessVector* pv = particleDef.GetProcessManager()->GetProcessList();
209 
210  G4cout << " 0: " << (*pv)[0]->GetProcessName() << " " << (*pv)[0]
211  << " 1: " << (*pv)[1]->GetProcessName() << " " << (*pv)[1]
212  // << " 2: " << (*pv)[2]->GetProcessName() << " " << (*pv)[2]
213  << G4endl;
214  G4cout << "ionModel= " << theIonEffChargeModel
215  << " MFPtable= " << theMeanFreePathTable
216  << " iniMass= " << initialMass
217  << G4endl;
218  }
219  }
220  // End of verbose print-out
221 
222  if (particleDef.GetParticleType() == "nucleus" &&
223  particleDef.GetParticleName() != "GenericIon" &&
224  particleDef.GetParticleSubType() == "generic")
225  {
226 
227  G4EnergyLossTables::Register(&particleDef,
234  proton_mass_c2/particleDef.GetPDGMass(),
235  TotBin);
236 
237  return;
238  }
239 
240  if( !CutsWhereModified() && theLossTable) return;
241 
244  G4AntiProton* antiproton = G4AntiProton::AntiProton();
245 
246  charge = particleDef.GetPDGCharge() / eplus;
248 
249  const G4ProductionCutsTable* theCoupleTable=
251  size_t numOfCouples = theCoupleTable->GetTableSize();
252 
253  cutForDelta.clear();
254  cutForGamma.clear();
255 
256  for (size_t j=0; j<numOfCouples; j++) {
257 
258  // get material parameters needed for the energy loss calculation
259  const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(j);
260  const G4Material* material= couple->GetMaterial();
261 
262  // the cut cannot be below lowest limit
263  G4double tCut = (*(theCoupleTable->GetEnergyCutsVector(1)))[j];
264  if(tCut > HighestKineticEnergy) tCut = HighestKineticEnergy;
265 
266  G4double excEnergy = material->GetIonisation()->GetMeanExcitationEnergy();
267 
268  tCut = std::max(tCut,excEnergy);
269  cutForDelta.push_back(tCut);
270 
271  // the cut cannot be below lowest limit
272  tCut = (*(theCoupleTable->GetEnergyCutsVector(0)))[j];
273  if(tCut > HighestKineticEnergy) tCut = HighestKineticEnergy;
274  tCut = std::max(tCut,minGammaEnergy);
275  cutForGamma.push_back(tCut);
276  }
277 
278  if(verboseLevel > 0) {
279  G4cout << "Cuts are defined " << G4endl;
280  }
281 
282  if(0.0 < charge)
283  {
284  {
285  BuildLossTable(*proton) ;
286 
287  // The following vector has a fixed dimension (see src/G4hImpactLoss.cc for more details)
288  // It happended in the past that caused memory corruption errors. The problem is still pending, even if temporary solved
289  // G4cout << "[NOTE]: __LINE__=" << __LINE__ << ", particleDef=" << particleDef.GetParticleName() << ", proton=" << proton << ", theLossTable=" << theLossTable << ", CounterOfpProcess=" << CounterOfpProcess << G4endl;
290 
293  }
294  } else {
295  {
296  BuildLossTable(*antiproton) ;
297 
298  // The following vector has a fixed dimension (see src/G4hImpactLoss.cc for more details)
299  // It happended in the past that caused memory corruption errors. The problem is still pending, even if temporary solved
300  // G4cout << "[NOTE]: __LINE__=" << __LINE__ << ", particleDef=" << particleDef.GetParticleName() << ", antiproton=" << antiproton << ", theLossTable=" << theLossTable << ", CounterOfpbarProcess=" << CounterOfpbarProcess << G4endl;
301 
304  }
305  }
306 
307  if(verboseLevel > 0) {
308  G4cout << "G4hImpactIonisation::BuildPhysicsTable: "
309  << "Loss table is built "
310  // << theLossTable
311  << G4endl;
312  }
313 
314  BuildLambdaTable(particleDef) ;
315  // BuildDataForFluorescence(particleDef);
316 
317  if(verboseLevel > 1) {
318  G4cout << (*theMeanFreePathTable) << G4endl;
319  }
320 
321  if(verboseLevel > 0) {
322  G4cout << "G4hImpactIonisation::BuildPhysicsTable: "
323  << "DEDX table will be built "
324  // << theDEDXpTable << " " << theDEDXpbarTable
325  // << " " << theRangepTable << " " << theRangepbarTable
326  << G4endl;
327  }
328 
329  BuildDEDXTable(particleDef) ;
330 
331  if(verboseLevel > 1) {
332  G4cout << (*theDEDXpTable) << G4endl;
333  }
334 
335  if((&particleDef == proton) || (&particleDef == antiproton)) PrintInfoDefinition() ;
336 
337  if(verboseLevel > 0) {
338  G4cout << "G4hImpactIonisation::BuildPhysicsTable: end for "
339  << particleDef.GetParticleName() << G4endl;
340  }
341 }
342 
343 
344 
345 
346 
347 // --------------------------------------------------------------------
349 {
350  // Initialisation
351  G4double lowEdgeEnergy , ionloss, ionlossBB, paramB ;
352  //G4double lowEnergy, highEnergy;
353  G4double highEnergy;
355 
356  if (particleDef == *proton)
357  {
358  //lowEnergy = protonLowEnergy ;
359  highEnergy = protonHighEnergy ;
360  charge = 1. ;
361  }
362  else
363  {
364  //lowEnergy = antiprotonLowEnergy ;
365  highEnergy = antiprotonHighEnergy ;
366  charge = -1. ;
367  }
368  chargeSquare = 1. ;
369 
370  const G4ProductionCutsTable* theCoupleTable=
372  size_t numOfCouples = theCoupleTable->GetTableSize();
373 
374  if ( theLossTable)
375  {
377  delete theLossTable;
378  }
379 
380  theLossTable = new G4PhysicsTable(numOfCouples);
381 
382  // loop for materials
383  for (size_t j=0; j<numOfCouples; j++) {
384 
385  // create physics vector and fill it
388  TotBin);
389 
390  // get material parameters needed for the energy loss calculation
391  const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(j);
392  const G4Material* material= couple->GetMaterial();
393 
394  if ( charge > 0.0 ) {
395  ionloss = ProtonParametrisedDEDX(couple,highEnergy) ;
396  } else {
397  ionloss = AntiProtonParametrisedDEDX(couple,highEnergy) ;
398  }
399 
400  ionlossBB = betheBlochModel->TheValue(&particleDef,material,highEnergy) ;
401  ionlossBB -= DeltaRaysEnergy(couple,highEnergy,proton_mass_c2) ;
402 
403 
404  paramB = ionloss/ionlossBB - 1.0 ;
405 
406  // now comes the loop for the kinetic energy values
407  for (G4int i = 0 ; i < TotBin ; i++) {
408  lowEdgeEnergy = aVector->GetLowEdgeEnergy(i) ;
409 
410  // low energy part for this material, parametrised energy loss formulae
411  if ( lowEdgeEnergy < highEnergy ) {
412 
413  if ( charge > 0.0 ) {
414  ionloss = ProtonParametrisedDEDX(couple,lowEdgeEnergy) ;
415  } else {
416  ionloss = AntiProtonParametrisedDEDX(couple,lowEdgeEnergy) ;
417  }
418 
419  } else {
420 
421  // high energy part for this material, Bethe-Bloch formula
422  ionloss = betheBlochModel->TheValue(proton,material,
423  lowEdgeEnergy) ;
424 
425  ionloss -= DeltaRaysEnergy(couple,lowEdgeEnergy,proton_mass_c2) ;
426 
427  ionloss *= (1.0 + paramB*highEnergy/lowEdgeEnergy) ;
428  }
429 
430  // now put the loss into the vector
431  if(verboseLevel > 1) {
432  G4cout << "E(MeV)= " << lowEdgeEnergy/MeV
433  << " dE/dx(MeV/mm)= " << ionloss*mm/MeV
434  << " in " << material->GetName() << G4endl;
435  }
436  aVector->PutValue(i,ionloss) ;
437  }
438  // Insert vector for this material into the table
439  theLossTable->insert(aVector) ;
440  }
441 }
442 
443 
444 
445 // --------------------------------------------------------------------
447 
448 {
449  // Build mean free path tables for the delta ray production process
450  // tables are built for MATERIALS
451 
452  if(verboseLevel > 1) {
453  G4cout << "G4hImpactIonisation::BuildLambdaTable for "
454  << particleDef.GetParticleName() << " is started" << G4endl;
455  }
456 
457 
458  G4double lowEdgeEnergy, value;
459  charge = particleDef.GetPDGCharge()/eplus ;
461  initialMass = particleDef.GetPDGMass();
462 
463  const G4ProductionCutsTable* theCoupleTable=
465  size_t numOfCouples = theCoupleTable->GetTableSize();
466 
467 
468  if (theMeanFreePathTable) {
470  delete theMeanFreePathTable;
471  }
472 
473  theMeanFreePathTable = new G4PhysicsTable(numOfCouples);
474 
475  // loop for materials
476 
477  for (size_t j=0 ; j < numOfCouples; j++) {
478 
479  //create physics vector then fill it ....
482  TotBin);
483 
484  // compute the (macroscopic) cross section first
485  const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(j);
486  const G4Material* material= couple->GetMaterial();
487 
488  const G4ElementVector* theElementVector = material->GetElementVector() ;
489  const G4double* theAtomicNumDensityVector = material->GetAtomicNumDensityVector();
490  const G4int numberOfElements = material->GetNumberOfElements() ;
491 
492  // get the electron kinetic energy cut for the actual material,
493  // it will be used in ComputeMicroscopicCrossSection
494  // ( it is the SAME for ALL the ELEMENTS in THIS MATERIAL )
495  // ------------------------------------------------------
496 
497  G4double deltaCut = cutForDelta[j];
498 
499  for ( G4int i = 0 ; i < TotBin ; i++ ) {
500  lowEdgeEnergy = aVector->GetLowEdgeEnergy(i) ;
501  G4double sigma = 0.0 ;
502  G4int Z;
503 
504  for (G4int iel=0; iel<numberOfElements; iel++ )
505  {
506  Z = (G4int) (*theElementVector)[iel]->GetZ();
507  // ---- MGP --- Corrected duplicated cross section calculation here from
508  // G4hLowEnergyIonisation original
509  G4double microCross = MicroscopicCrossSection( particleDef,
510  lowEdgeEnergy,
511  Z,
512  deltaCut ) ;
513  //totalCrossSectionMap [Z] = microCross;
514  sigma += theAtomicNumDensityVector[iel] * microCross ;
515  }
516 
517  // mean free path = 1./macroscopic cross section
518 
519  value = sigma<=0 ? DBL_MAX : 1./sigma ;
520 
521  aVector->PutValue(i, value) ;
522  }
523 
524  theMeanFreePathTable->insert(aVector);
525  }
526 
527 }
528 
529 
530 // --------------------------------------------------------------------
532  G4double kineticEnergy,
533  G4double atomicNumber,
534  G4double deltaCutInEnergy) const
535 {
536  //******************************************************************
537  // cross section formula is OK for spin=0, 1/2, 1 only !
538  // *****************************************************************
539 
540  // Calculates the microscopic cross section in GEANT4 internal units
541  // Formula documented in Geant4 Phys. Ref. Manual
542  // ( it is called for elements, AtomicNumber = z )
543 
544  G4double totalCrossSection = 0.;
545 
546  // Particle mass and energy
547  G4double particleMass = initialMass;
548  G4double energy = kineticEnergy + particleMass;
549 
550  // Some kinematics
551  G4double gamma = energy / particleMass;
552  G4double beta2 = 1. - 1. / (gamma * gamma);
553  G4double var = electron_mass_c2 / particleMass;
554  G4double tMax = 2. * electron_mass_c2 * (gamma*gamma - 1.) / (1. + 2.* gamma*var + var*var);
555 
556  // Calculate the total cross section
557 
558  if ( tMax > deltaCutInEnergy )
559  {
560  var = deltaCutInEnergy / tMax;
561  totalCrossSection = (1. - var * (1. - beta2 * std::log(var))) / deltaCutInEnergy ;
562 
563  G4double spin = particleDef.GetPDGSpin() ;
564 
565  // +term for spin=1/2 particle
566  if (spin == 0.5)
567  {
568  totalCrossSection += 0.5 * (tMax - deltaCutInEnergy) / (energy*energy);
569  }
570  // +term for spin=1 particle
571  else if (spin > 0.9 )
572  {
573  totalCrossSection += -std::log(var) /
574  (3. * deltaCutInEnergy) + (tMax - deltaCutInEnergy) * ( (5. + 1. /var)*0.25 / (energy*energy) -
575  beta2 / (tMax * deltaCutInEnergy) ) / 3. ;
576  }
577  totalCrossSection *= twopi_mc2_rcl2 * atomicNumber / beta2 ;
578  }
579 
580  //std::cout << "Microscopic = " << totalCrossSection/barn
581  // << ", e = " << kineticEnergy/MeV <<std:: endl;
582 
583  return totalCrossSection ;
584 }
585 
586 
587 
588 // --------------------------------------------------------------------
590  G4double, // previousStepSize
592 {
593  const G4DynamicParticle* dynamicParticle = track.GetDynamicParticle();
594  const G4MaterialCutsCouple* couple = track.GetMaterialCutsCouple();
595  const G4Material* material = couple->GetMaterial();
596 
597  G4double meanFreePath = DBL_MAX;
598  // ---- MGP ---- What is the meaning of the local variable isOutOfRange?
599  G4bool isOutRange = false;
600 
601  *condition = NotForced ;
602 
603  G4double kineticEnergy = (dynamicParticle->GetKineticEnergy())*initialMass/(dynamicParticle->GetMass());
604  charge = dynamicParticle->GetCharge()/eplus;
605  chargeSquare = theIonEffChargeModel->TheValue(dynamicParticle, material);
606 
607  if (kineticEnergy < LowestKineticEnergy)
608  {
609  meanFreePath = DBL_MAX;
610  }
611  else
612  {
613  if (kineticEnergy > HighestKineticEnergy) kineticEnergy = HighestKineticEnergy;
614  meanFreePath = (((*theMeanFreePathTable)(couple->GetIndex()))->
615  GetValue(kineticEnergy,isOutRange))/chargeSquare;
616  }
617 
618  return meanFreePath ;
619 }
620 
621 
622 // --------------------------------------------------------------------
624  const G4MaterialCutsCouple* couple)
625 {
626  // returns the Step limit
627  // dEdx is calculated as well as the range
628  // based on Effective Charge Approach
629 
630  const G4Material* material = couple->GetMaterial();
632  G4AntiProton* antiproton = G4AntiProton::AntiProton();
633 
634  G4double stepLimit = 0.;
635  G4double dx, highEnergy;
636 
637  G4double massRatio = proton_mass_c2/(particle->GetMass()) ;
638  G4double kineticEnergy = particle->GetKineticEnergy() ;
639 
640  // Scale the kinetic energy
641 
642  G4double tScaled = kineticEnergy*massRatio ;
643  fBarkas = 0.;
644 
645  if (charge > 0.)
646  {
647  highEnergy = protonHighEnergy ;
648  fRangeNow = G4EnergyLossTables::GetRange(proton, tScaled, couple);
649  dx = G4EnergyLossTables::GetRange(proton, highEnergy, couple);
650  fdEdx = G4EnergyLossTables::GetDEDX(proton, tScaled, couple)
651  * chargeSquare ;
652 
653  // Correction for positive ions
654  if (theBarkas && tScaled > highEnergy)
655  {
656  fBarkas = BarkasTerm(material,tScaled)*std::sqrt(chargeSquare)*chargeSquare
657  + BlochTerm(material,tScaled,chargeSquare);
658  }
659  // Antiprotons and negative hadrons
660  }
661  else
662  {
663  highEnergy = antiprotonHighEnergy ;
664  fRangeNow = G4EnergyLossTables::GetRange(antiproton, tScaled, couple);
665  dx = G4EnergyLossTables::GetRange(antiproton, highEnergy, couple);
666  fdEdx = G4EnergyLossTables::GetDEDX(antiproton, tScaled, couple) * chargeSquare ;
667 
668  if (theBarkas && tScaled > highEnergy)
669  {
670  fBarkas = -BarkasTerm(material,tScaled)*std::sqrt(chargeSquare)*chargeSquare
671  + BlochTerm(material,tScaled,chargeSquare);
672  }
673  }
674  /*
675  const G4Material* mat = couple->GetMaterial();
676  G4double fac = gram/(MeV*cm2*mat->GetDensity());
677  G4cout << particle->GetDefinition()->GetParticleName()
678  << " in " << mat->GetName()
679  << " E(MeV)= " << kineticEnergy/MeV
680  << " dedx(MeV*cm^2/g)= " << fdEdx*fac
681  << " barcas(MeV*cm^2/gram)= " << fBarkas*fac
682  << " Q^2= " << chargeSquare
683  << G4endl;
684  */
685  // scaling back
686  fRangeNow /= (chargeSquare*massRatio) ;
687  dx /= (chargeSquare*massRatio) ;
688 
689  stepLimit = fRangeNow ;
692 
693  if (fRangeNow > r)
694  {
695  stepLimit = dRoverRange*fRangeNow + r*(1.0 - dRoverRange)*(2.0 - r/fRangeNow);
696  if (stepLimit > fRangeNow) stepLimit = fRangeNow;
697  }
698  // compute the (random) Step limit in standard energy range
699  if(tScaled > highEnergy )
700  {
701  // add Barkas correction directly to dedx
702  fdEdx += fBarkas;
703 
704  if(stepLimit > fRangeNow - dx*0.9) stepLimit = fRangeNow - dx*0.9 ;
705 
706  // Step limit in low energy range
707  }
708  else
709  {
710  G4double x = dx*paramStepLimit;
711  if (stepLimit > x) stepLimit = x;
712  }
713  return stepLimit;
714 }
715 
716 
717 // --------------------------------------------------------------------
719  const G4Step& step)
720 {
721  // compute the energy loss after a step
723  G4AntiProton* antiproton = G4AntiProton::AntiProton();
724  G4double finalT = 0.;
725 
726  aParticleChange.Initialize(track) ;
727 
728  const G4MaterialCutsCouple* couple = track.GetMaterialCutsCouple();
729  const G4Material* material = couple->GetMaterial();
730 
731  // get the actual (true) Step length from step
732  const G4double stepLength = step.GetStepLength() ;
733 
734  const G4DynamicParticle* particle = track.GetDynamicParticle() ;
735 
736  G4double kineticEnergy = particle->GetKineticEnergy() ;
737  G4double massRatio = proton_mass_c2/(particle->GetMass()) ;
738  G4double tScaled = kineticEnergy * massRatio ;
739  G4double eLoss = 0.0 ;
740  G4double nLoss = 0.0 ;
741 
742 
743  // very small particle energy
744  if (kineticEnergy < MinKineticEnergy)
745  {
746  eLoss = kineticEnergy ;
747  // particle energy outside tabulated energy range
748  }
749 
750  else if( kineticEnergy > HighestKineticEnergy)
751  {
752  eLoss = stepLength * fdEdx ;
753  // big step
754  }
755  else if (stepLength >= fRangeNow )
756  {
757  eLoss = kineticEnergy ;
758 
759  // tabulated range
760  }
761  else
762  {
763  // step longer than linear step limit
764  if(stepLength > linLossLimit * fRangeNow)
765  {
766  G4double rScaled = fRangeNow * massRatio * chargeSquare ;
767  G4double sScaled = stepLength * massRatio * chargeSquare ;
768 
769  if(charge > 0.0)
770  {
771  eLoss = G4EnergyLossTables::GetPreciseEnergyFromRange(proton,rScaled, couple) -
772  G4EnergyLossTables::GetPreciseEnergyFromRange(proton,rScaled-sScaled,couple) ;
773 
774  }
775  else
776  {
777  // Antiproton
778  eLoss = G4EnergyLossTables::GetPreciseEnergyFromRange(antiproton,rScaled,couple) -
779  G4EnergyLossTables::GetPreciseEnergyFromRange(antiproton,rScaled-sScaled,couple) ;
780  }
781  eLoss /= massRatio ;
782 
783  // Barkas correction at big step
784  eLoss += fBarkas * stepLength;
785 
786  // step shorter than linear step limit
787  }
788  else
789  {
790  eLoss = stepLength *fdEdx ;
791  }
792  if (nStopping && tScaled < protonHighEnergy)
793  {
794  nLoss = (theNuclearStoppingModel->TheValue(particle, material)) * stepLength;
795  }
796  }
797 
798  if (eLoss < 0.0) eLoss = 0.0;
799 
800  finalT = kineticEnergy - eLoss - nLoss;
801 
802  if ( EnlossFlucFlag && 0.0 < eLoss && finalT > MinKineticEnergy)
803  {
804 
805  // now the electron loss with fluctuation
806  eLoss = ElectronicLossFluctuation(particle, couple, eLoss, stepLength) ;
807  if (eLoss < 0.0) eLoss = 0.0;
808  finalT = kineticEnergy - eLoss - nLoss;
809  }
810 
811  // stop particle if the kinetic energy <= MinKineticEnergy
812  if (finalT*massRatio <= MinKineticEnergy )
813  {
814 
815  finalT = 0.0;
818  else
820  }
821 
822  aParticleChange.ProposeEnergy( finalT );
823  eLoss = kineticEnergy-finalT;
824 
826  return &aParticleChange ;
827 }
828 
829 
830 
831 // --------------------------------------------------------------------
833  G4double kineticEnergy) const
834 {
835  const G4Material* material = couple->GetMaterial();
837  G4double eLoss = 0.;
838 
839  // Free Electron Gas Model
840  if(kineticEnergy < protonLowEnergy) {
841  eLoss = (protonModel->TheValue(proton, material, protonLowEnergy))
842  * std::sqrt(kineticEnergy/protonLowEnergy) ;
843 
844  // Parametrisation
845  } else {
846  eLoss = protonModel->TheValue(proton, material, kineticEnergy) ;
847  }
848 
849  // Delta rays energy
850  eLoss -= DeltaRaysEnergy(couple,kineticEnergy,proton_mass_c2) ;
851 
852  if(verboseLevel > 2) {
853  G4cout << "p E(MeV)= " << kineticEnergy/MeV
854  << " dE/dx(MeV/mm)= " << eLoss*mm/MeV
855  << " for " << material->GetName()
856  << " model: " << protonModel << G4endl;
857  }
858 
859  if(eLoss < 0.0) eLoss = 0.0 ;
860 
861  return eLoss ;
862 }
863 
864 
865 
866 // --------------------------------------------------------------------
868  G4double kineticEnergy) const
869 {
870  const G4Material* material = couple->GetMaterial();
871  G4AntiProton* antiproton = G4AntiProton::AntiProton();
872  G4double eLoss = 0.0 ;
873 
874  // Antiproton model is used
875  if(antiprotonModel->IsInCharge(antiproton,material)) {
876  if(kineticEnergy < antiprotonLowEnergy) {
877  eLoss = antiprotonModel->TheValue(antiproton,material,antiprotonLowEnergy)
878  * std::sqrt(kineticEnergy/antiprotonLowEnergy) ;
879 
880  // Parametrisation
881  } else {
882  eLoss = antiprotonModel->TheValue(antiproton,material,
883  kineticEnergy);
884  }
885 
886  // The proton model is used + Barkas correction
887  } else {
888  if(kineticEnergy < protonLowEnergy) {
890  * std::sqrt(kineticEnergy/protonLowEnergy) ;
891 
892  // Parametrisation
893  } else {
894  eLoss = protonModel->TheValue(G4Proton::Proton(),material,
895  kineticEnergy);
896  }
897  //if(theBarkas) eLoss -= 2.0*BarkasTerm(material, kineticEnergy);
898  }
899 
900  // Delta rays energy
901  eLoss -= DeltaRaysEnergy(couple,kineticEnergy,proton_mass_c2) ;
902 
903  if(verboseLevel > 2) {
904  G4cout << "pbar E(MeV)= " << kineticEnergy/MeV
905  << " dE/dx(MeV/mm)= " << eLoss*mm/MeV
906  << " for " << material->GetName()
907  << " model: " << protonModel << G4endl;
908  }
909 
910  if(eLoss < 0.0) eLoss = 0.0 ;
911 
912  return eLoss ;
913 }
914 
915 
916 // --------------------------------------------------------------------
918  G4double kineticEnergy,
919  G4double particleMass) const
920 {
921  G4double dLoss = 0.;
922 
923  G4double deltaCutNow = cutForDelta[(couple->GetIndex())] ;
924  const G4Material* material = couple->GetMaterial();
925  G4double electronDensity = material->GetElectronDensity();
926  G4double excitationEnergy = material->GetIonisation()->GetMeanExcitationEnergy();
927 
928  G4double tau = kineticEnergy / particleMass ;
929  G4double rateMass = electron_mass_c2/particleMass ;
930 
931  // some local variables
932 
933  G4double gamma = tau + 1.0 ;
934  G4double bg2 = tau*(tau+2.0) ;
935  G4double beta2 = bg2/(gamma*gamma) ;
936  G4double tMax = 2.*electron_mass_c2*bg2/(1.0+2.0*gamma*rateMass+rateMass*rateMass) ;
937 
938  // Validity range for delta electron cross section
939  G4double deltaCut = std::max(deltaCutNow, excitationEnergy);
940 
941  if ( deltaCut < tMax)
942  {
943  G4double x = deltaCut / tMax ;
944  dLoss = ( beta2 * (x-1.) - std::log(x) ) * twopi_mc2_rcl2 * electronDensity / beta2 ;
945  }
946  return dLoss ;
947 }
948 
949 
950 // -------------------------------------------------------------------------
951 
953  const G4Step& step)
954 {
955  // Units are expressed in GEANT4 internal units.
956 
957  // std::cout << "----- Calling PostStepDoIt ----- " << std::endl;
958 
959  aParticleChange.Initialize(track) ;
960  const G4MaterialCutsCouple* couple = track.GetMaterialCutsCouple();
961  const G4DynamicParticle* aParticle = track.GetDynamicParticle() ;
962 
963  // Some kinematics
964 
965  G4ParticleDefinition* definition = track.GetDefinition();
966  G4double mass = definition->GetPDGMass();
967  G4double kineticEnergy = aParticle->GetKineticEnergy();
968  G4double totalEnergy = kineticEnergy + mass ;
969  G4double pSquare = kineticEnergy *( totalEnergy + mass) ;
970  G4double eSquare = totalEnergy * totalEnergy;
971  G4double betaSquare = pSquare / eSquare;
972  G4ThreeVector particleDirection = aParticle->GetMomentumDirection() ;
973 
974  G4double gamma = kineticEnergy / mass + 1.;
975  G4double r = electron_mass_c2 / mass;
976  G4double tMax = 2. * electron_mass_c2 *(gamma*gamma - 1.) / (1. + 2.*gamma*r + r*r);
977 
978  // Validity range for delta electron cross section
979  G4double deltaCut = cutForDelta[couple->GetIndex()];
980 
981  // This should not be a case
982  if (deltaCut >= tMax)
984 
985  G4double xc = deltaCut / tMax;
986  G4double rate = tMax / totalEnergy;
987  rate = rate*rate ;
988  G4double spin = aParticle->GetDefinition()->GetPDGSpin() ;
989 
990  // Sampling follows ...
991  G4double x = 0.;
992  G4double gRej = 0.;
993 
994  do {
995  x = xc / (1. - (1. - xc) * G4UniformRand());
996 
997  if (0.0 == spin)
998  {
999  gRej = 1.0 - betaSquare * x ;
1000  }
1001  else if (0.5 == spin)
1002  {
1003  gRej = (1. - betaSquare * x + 0.5 * x*x * rate) / (1. + 0.5 * rate) ;
1004  }
1005  else
1006  {
1007  gRej = (1. - betaSquare * x ) * (1. + x/(3.*xc)) +
1008  x*x * rate * (1. + 0.5*x/xc) / 3.0 /
1009  (1. + 1./(3.*xc) + rate *(1.+ 0.5/xc) / 3.);
1010  }
1011 
1012  } while( G4UniformRand() > gRej );
1013 
1014  G4double deltaKineticEnergy = x * tMax;
1015  G4double deltaTotalMomentum = std::sqrt(deltaKineticEnergy *
1016  (deltaKineticEnergy + 2. * electron_mass_c2 ));
1017  G4double totalMomentum = std::sqrt(pSquare) ;
1018  G4double cosTheta = deltaKineticEnergy * (totalEnergy + electron_mass_c2) / (deltaTotalMomentum*totalMomentum);
1019 
1020  // protection against cosTheta > 1 or < -1
1021  if ( cosTheta < -1. ) cosTheta = -1.;
1022  if ( cosTheta > 1. ) cosTheta = 1.;
1023 
1024  // direction of the delta electron
1025  G4double phi = twopi * G4UniformRand();
1026  G4double sinTheta = std::sqrt(1. - cosTheta*cosTheta);
1027  G4double dirX = sinTheta * std::cos(phi);
1028  G4double dirY = sinTheta * std::sin(phi);
1029  G4double dirZ = cosTheta;
1030 
1031  G4ThreeVector deltaDirection(dirX,dirY,dirZ);
1032  deltaDirection.rotateUz(particleDirection);
1033 
1034  // create G4DynamicParticle object for delta ray
1035  G4DynamicParticle* deltaRay = new G4DynamicParticle;
1036  deltaRay->SetKineticEnergy( deltaKineticEnergy );
1037  deltaRay->SetMomentumDirection(deltaDirection.x(),
1038  deltaDirection.y(),
1039  deltaDirection.z());
1040  deltaRay->SetDefinition(G4Electron::Electron());
1041 
1042  // fill aParticleChange
1043  G4double finalKineticEnergy = kineticEnergy - deltaKineticEnergy;
1044  size_t totalNumber = 1;
1045 
1046  // Atomic relaxation
1047 
1048  // ---- MGP ---- Temporary limitation: currently PIXE only for incident protons
1049 
1050  size_t nSecondaries = 0;
1051  std::vector<G4DynamicParticle*>* secondaryVector = 0;
1052 
1053  if (definition == G4Proton::ProtonDefinition())
1054  {
1055  const G4Material* material = couple->GetMaterial();
1056 
1057  // Lazy initialization of pixeCrossSectionHandler
1058  if (pixeCrossSectionHandler == 0)
1059  {
1060  // Instantiate pixeCrossSectionHandler with selected shell cross section models
1061  // Ownership of interpolation is transferred to pixeCrossSectionHandler
1062  G4IInterpolator* interpolation = new G4LogLogInterpolator();
1065  G4String fileName("proton");
1067  // pixeCrossSectionHandler->PrintData();
1068  }
1069 
1070  // Select an atom in the current material based on the total shell cross sections
1071  G4int Z = pixeCrossSectionHandler->SelectRandomAtom(material,kineticEnergy);
1072  // std::cout << "G4hImpactIonisation::PostStepDoIt - Z = " << Z << std::endl;
1073 
1074  // G4double microscopicCross = MicroscopicCrossSection(*definition,
1075  // kineticEnergy,
1076  // Z, deltaCut);
1077  // G4double crossFromShells = pixeCrossSectionHandler->FindValue(Z,kineticEnergy);
1078 
1079  //std::cout << "G4hImpactIonisation: Z= "
1080  // << Z
1081  // << ", energy = "
1082  // << kineticEnergy/MeV
1083  // <<" MeV, microscopic = "
1084  // << microscopicCross/barn
1085  // << " barn, from shells = "
1086  // << crossFromShells/barn
1087  // << " barn"
1088  // << std::endl;
1089 
1090  // Select a shell in the target atom based on the individual shell cross sections
1091  G4int shellIndex = pixeCrossSectionHandler->SelectRandomShell(Z,kineticEnergy);
1092 
1094  const G4AtomicShell* atomicShell = transitionManager->Shell(Z,shellIndex);
1095  G4double bindingEnergy = atomicShell->BindingEnergy();
1096 
1097  // if (verboseLevel > 1)
1098  // {
1099  // G4cout << "G4hImpactIonisation::PostStepDoIt - Z = "
1100  // << Z
1101  // << ", shell = "
1102  // << shellIndex
1103  // << ", bindingE (keV) = "
1104  // << bindingEnergy/keV
1105  // << G4endl;
1106  // }
1107 
1108  // Generate PIXE if binding energy larger than cut for photons or electrons
1109 
1110  G4ParticleDefinition* type = 0;
1111 
1112  if (finalKineticEnergy >= bindingEnergy)
1113  // && (bindingEnergy >= minGammaEnergy || bindingEnergy >= minElectronEnergy) )
1114  {
1115  // Vacancy in subshell shellIndex; shellId is the subshell identifier in EADL jargon
1116  G4int shellId = atomicShell->ShellId();
1117  // Atomic relaxation: generate secondaries
1118  secondaryVector = atomicDeexcitation.GenerateParticles(Z, shellId);
1119 
1120  // ---- Debug ----
1121  //std::cout << "ShellId = "
1122  // <<shellId << " ---- Atomic relaxation secondaries: ---- "
1123  // << secondaryVector->size()
1124  // << std::endl;
1125 
1126  // ---- End debug ---
1127 
1128  if (secondaryVector != 0)
1129  {
1130  nSecondaries = secondaryVector->size();
1131  for (size_t i = 0; i<nSecondaries; i++)
1132  {
1133  G4DynamicParticle* aSecondary = (*secondaryVector)[i];
1134  if (aSecondary)
1135  {
1136  G4double e = aSecondary->GetKineticEnergy();
1137  type = aSecondary->GetDefinition();
1138 
1139  // ---- Debug ----
1140  //if (type == G4Gamma::GammaDefinition())
1141  // {
1142  // std::cout << "Z = " << Z
1143  // << ", shell: " << shellId
1144  // << ", PIXE photon energy (keV) = " << e/keV
1145  // << std::endl;
1146  // }
1147  // ---- End debug ---
1148 
1149  if (e < finalKineticEnergy &&
1150  ((type == G4Gamma::Gamma() && e > minGammaEnergy ) ||
1151  (type == G4Electron::Electron() && e > minElectronEnergy )))
1152  {
1153  // Subtract the energy of the emitted secondary from the primary
1154  finalKineticEnergy -= e;
1155  totalNumber++;
1156  // ---- Debug ----
1157  //if (type == G4Gamma::GammaDefinition())
1158  // {
1159  // std::cout << "Z = " << Z
1160  // << ", shell: " << shellId
1161  // << ", PIXE photon energy (keV) = " << e/keV
1162  // << std::endl;
1163  // }
1164  // ---- End debug ---
1165  }
1166  else
1167  {
1168  // The atomic relaxation product has energy below the cut
1169  // ---- Debug ----
1170  // if (type == G4Gamma::GammaDefinition())
1171  // {
1172  // std::cout << "Z = " << Z
1173  //
1174  // << ", PIXE photon energy = " << e/keV
1175  // << " keV below threshold " << minGammaEnergy/keV << " keV"
1176  // << std::endl;
1177  // }
1178  // ---- End debug ---
1179 
1180  delete aSecondary;
1181  (*secondaryVector)[i] = 0;
1182  }
1183  }
1184  }
1185  }
1186  }
1187  }
1188 
1189 
1190  // Save delta-electrons
1191 
1192  G4double eDeposit = 0.;
1193 
1194  if (finalKineticEnergy > MinKineticEnergy)
1195  {
1196  G4double finalPx = totalMomentum*particleDirection.x() - deltaTotalMomentum*deltaDirection.x();
1197  G4double finalPy = totalMomentum*particleDirection.y() - deltaTotalMomentum*deltaDirection.y();
1198  G4double finalPz = totalMomentum*particleDirection.z() - deltaTotalMomentum*deltaDirection.z();
1199  G4double finalMomentum = std::sqrt(finalPx*finalPx + finalPy*finalPy + finalPz*finalPz) ;
1200  finalPx /= finalMomentum;
1201  finalPy /= finalMomentum;
1202  finalPz /= finalMomentum;
1203 
1204  aParticleChange.ProposeMomentumDirection( finalPx,finalPy,finalPz );
1205  }
1206  else
1207  {
1208  eDeposit = finalKineticEnergy;
1209  finalKineticEnergy = 0.;
1210  aParticleChange.ProposeMomentumDirection(particleDirection.x(),
1211  particleDirection.y(),
1212  particleDirection.z());
1213  if(!aParticle->GetDefinition()->GetProcessManager()->
1214  GetAtRestProcessVector()->size())
1216  else
1218  }
1219 
1220  aParticleChange.ProposeEnergy(finalKineticEnergy);
1223  aParticleChange.AddSecondary(deltaRay);
1224 
1225  // ---- Debug ----
1226  // std::cout << "RDHadronIonisation - finalKineticEnergy (MeV) = "
1227  // << finalKineticEnergy/MeV
1228  // << ", delta KineticEnergy (keV) = "
1229  // << deltaKineticEnergy/keV
1230  // << ", energy deposit (MeV) = "
1231  // << eDeposit/MeV
1232  // << std::endl;
1233  // ---- End debug ---
1234 
1235  // Save Fluorescence and Auger
1236 
1237  if (secondaryVector != 0)
1238  {
1239  for (size_t l = 0; l < nSecondaries; l++)
1240  {
1241  G4DynamicParticle* secondary = (*secondaryVector)[l];
1242  if (secondary) aParticleChange.AddSecondary(secondary);
1243 
1244  // ---- Debug ----
1245  //if (secondary != 0)
1246  // {
1247  // if (secondary->GetDefinition() == G4Gamma::GammaDefinition())
1248  // {
1249  // G4double eX = secondary->GetKineticEnergy();
1250  // std::cout << " PIXE photon of energy " << eX/keV
1251  // << " keV added to ParticleChange; total number of secondaries is " << totalNumber
1252  // << std::endl;
1253  // }
1254  //}
1255  // ---- End debug ---
1256 
1257  }
1258  delete secondaryVector;
1259  }
1260 
1261  return G4VContinuousDiscreteProcess::PostStepDoIt(track,step);
1262 }
1263 
1264 // -------------------------------------------------------------------------
1265 
1267  const G4MaterialCutsCouple* couple,
1268  G4double kineticEnergy)
1269 {
1270  const G4Material* material = couple->GetMaterial();
1272  G4AntiProton* antiproton = G4AntiProton::AntiProton();
1273  G4double dedx = 0.;
1274 
1275  G4double tScaled = kineticEnergy * proton_mass_c2 / (aParticle->GetPDGMass()) ;
1276  charge = aParticle->GetPDGCharge() ;
1277 
1278  if( charge > 0.)
1279  {
1280  if (tScaled > protonHighEnergy)
1281  {
1282  dedx = G4EnergyLossTables::GetDEDX(proton,tScaled,couple) ;
1283  }
1284  else
1285  {
1286  dedx = ProtonParametrisedDEDX(couple,tScaled) ;
1287  }
1288  }
1289  else
1290  {
1291  if (tScaled > antiprotonHighEnergy)
1292  {
1293  dedx = G4EnergyLossTables::GetDEDX(antiproton,tScaled,couple);
1294  }
1295  else
1296  {
1297  dedx = AntiProtonParametrisedDEDX(couple,tScaled) ;
1298  }
1299  }
1300  dedx *= theIonEffChargeModel->TheValue(aParticle, material, kineticEnergy) ;
1301 
1302  return dedx ;
1303 }
1304 
1305 
1307  G4double kineticEnergy) const
1308 //Function to compute the Barkas term for protons:
1309 //
1310 //Ref. Z_1^3 effect in the stopping power of matter for charged particles
1311 // J.C Ashley and R.H.Ritchie
1312 // Physical review B Vol.5 No.7 1 April 1972 pagg. 2393-2397
1313 //
1314 {
1315  static G4ThreadLocal G4double FTable[47][2] = {
1316  { 0.02, 21.5},
1317  { 0.03, 20.0},
1318  { 0.04, 18.0},
1319  { 0.05, 15.6},
1320  { 0.06, 15.0},
1321  { 0.07, 14.0},
1322  { 0.08, 13.5},
1323  { 0.09, 13.},
1324  { 0.1, 12.2},
1325  { 0.2, 9.25},
1326  { 0.3, 7.0},
1327  { 0.4, 6.0},
1328  { 0.5, 4.5},
1329  { 0.6, 3.5},
1330  { 0.7, 3.0},
1331  { 0.8, 2.5},
1332  { 0.9, 2.0},
1333  { 1.0, 1.7},
1334  { 1.2, 1.2},
1335  { 1.3, 1.0},
1336  { 1.4, 0.86},
1337  { 1.5, 0.7},
1338  { 1.6, 0.61},
1339  { 1.7, 0.52},
1340  { 1.8, 0.5},
1341  { 1.9, 0.43},
1342  { 2.0, 0.42},
1343  { 2.1, 0.3},
1344  { 2.4, 0.2},
1345  { 3.0, 0.13},
1346  { 3.08, 0.1},
1347  { 3.1, 0.09},
1348  { 3.3, 0.08},
1349  { 3.5, 0.07},
1350  { 3.8, 0.06},
1351  { 4.0, 0.051},
1352  { 4.1, 0.04},
1353  { 4.8, 0.03},
1354  { 5.0, 0.024},
1355  { 5.1, 0.02},
1356  { 6.0, 0.013},
1357  { 6.5, 0.01},
1358  { 7.0, 0.009},
1359  { 7.1, 0.008},
1360  { 8.0, 0.006},
1361  { 9.0, 0.0032},
1362  { 10.0, 0.0025} };
1363 
1364  // Information on particle and material
1365  G4double kinE = kineticEnergy ;
1366  if(0.5*MeV > kinE) kinE = 0.5*MeV ;
1367  G4double gamma = 1.0 + kinE / proton_mass_c2 ;
1368  G4double beta2 = 1.0 - 1.0/(gamma*gamma) ;
1369  if(0.0 >= beta2) return 0.0;
1370 
1371  G4double BTerm = 0.0;
1372  //G4double AMaterial = 0.0;
1373  G4double ZMaterial = 0.0;
1374  const G4ElementVector* theElementVector = material->GetElementVector();
1375  G4int numberOfElements = material->GetNumberOfElements();
1376 
1377  for (G4int i = 0; i<numberOfElements; i++) {
1378 
1379  //AMaterial = (*theElementVector)[i]->GetA()*mole/g;
1380  ZMaterial = (*theElementVector)[i]->GetZ();
1381 
1382  G4double X = 137.0 * 137.0 * beta2 / ZMaterial;
1383 
1384  // Variables to compute L_1
1385  G4double Eta0Chi = 0.8;
1386  G4double EtaChi = Eta0Chi * ( 1.0 + 6.02*std::pow( ZMaterial,-1.19 ) );
1387  G4double W = ( EtaChi * std::pow( ZMaterial,1.0/6.0 ) ) / std::sqrt(X);
1388  G4double FunctionOfW = FTable[46][1]*FTable[46][0]/W ;
1389 
1390  for(G4int j=0; j<47; j++) {
1391 
1392  if( W < FTable[j][0] ) {
1393 
1394  if(0 == j) {
1395  FunctionOfW = FTable[0][1] ;
1396 
1397  } else {
1398  FunctionOfW = (FTable[j][1] - FTable[j-1][1]) * (W - FTable[j-1][0])
1399  / (FTable[j][0] - FTable[j-1][0])
1400  + FTable[j-1][1] ;
1401  }
1402 
1403  break;
1404  }
1405 
1406  }
1407 
1408  BTerm += FunctionOfW /( std::sqrt(ZMaterial * X) * X);
1409  }
1410 
1411  BTerm *= twopi_mc2_rcl2 * (material->GetElectronDensity()) / beta2 ;
1412 
1413  return BTerm;
1414 }
1415 
1416 
1417 
1419  G4double kineticEnergy,
1420  G4double cSquare) const
1421 //Function to compute the Bloch term for protons:
1422 //
1423 //Ref. Z_1^3 effect in the stopping power of matter for charged particles
1424 // J.C Ashley and R.H.Ritchie
1425 // Physical review B Vol.5 No.7 1 April 1972 pagg. 2393-2397
1426 //
1427 {
1428  G4double eLoss = 0.0 ;
1429  G4double gamma = 1.0 + kineticEnergy / proton_mass_c2 ;
1430  G4double beta2 = 1.0 - 1.0/(gamma*gamma) ;
1431  G4double y = cSquare / (137.0*137.0*beta2) ;
1432 
1433  if(y < 0.05) {
1434  eLoss = 1.202 ;
1435 
1436  } else {
1437  eLoss = 1.0 / (1.0 + y) ;
1438  G4double de = eLoss ;
1439 
1440  for(G4int i=2; de>eLoss*0.01; i++) {
1441  de = 1.0/( i * (i*i + y)) ;
1442  eLoss += de ;
1443  }
1444  }
1445  eLoss *= -1.0 * y * cSquare * twopi_mc2_rcl2 *
1446  (material->GetElectronDensity()) / beta2 ;
1447 
1448  return eLoss;
1449 }
1450 
1451 
1452 
1454  const G4DynamicParticle* particle,
1455  const G4MaterialCutsCouple* couple,
1456  G4double meanLoss,
1457  G4double step) const
1458 // calculate actual loss from the mean loss
1459 // The model used to get the fluctuation is essentially the same
1460 // as in Glandz in Geant3.
1461 {
1462  // data members to speed up the fluctuation calculation
1463  // G4int imat ;
1464  // G4double f1Fluct,f2Fluct,e1Fluct,e2Fluct,rateFluct,ipotFluct;
1465  // G4double e1LogFluct,e2LogFluct,ipotLogFluct;
1466 
1467  static const G4double minLoss = 1.*eV ;
1468  static const G4double kappa = 10. ;
1469  static const G4double theBohrBeta2 = 50.0 * keV/proton_mass_c2 ;
1470 
1471  const G4Material* material = couple->GetMaterial();
1472  G4int imaterial = couple->GetIndex() ;
1473  G4double ipotFluct = material->GetIonisation()->GetMeanExcitationEnergy() ;
1474  G4double electronDensity = material->GetElectronDensity() ;
1475  G4double zeff = electronDensity/(material->GetTotNbOfAtomsPerVolume()) ;
1476 
1477  // get particle data
1478  G4double tkin = particle->GetKineticEnergy();
1479  G4double particleMass = particle->GetMass() ;
1480  G4double deltaCutInKineticEnergyNow = cutForDelta[imaterial];
1481 
1482  // shortcut for very very small loss
1483  if(meanLoss < minLoss) return meanLoss ;
1484 
1485  // Validity range for delta electron cross section
1486  G4double threshold = std::max(deltaCutInKineticEnergyNow,ipotFluct);
1487  G4double loss, siga;
1488 
1489  G4double rmass = electron_mass_c2/particleMass;
1490  G4double tau = tkin/particleMass;
1491  G4double tau1 = tau+1.0;
1492  G4double tau2 = tau*(tau+2.);
1493  G4double tMax = 2.*electron_mass_c2*tau2/(1.+2.*tau1*rmass+rmass*rmass);
1494 
1495 
1496  if(tMax > threshold) tMax = threshold;
1497  G4double beta2 = tau2/(tau1*tau1);
1498 
1499  // Gaussian fluctuation
1500  if(meanLoss > kappa*tMax || tMax < kappa*ipotFluct )
1501  {
1502  siga = tMax * (1.0-0.5*beta2) * step * twopi_mc2_rcl2
1503  * electronDensity / beta2 ;
1504 
1505  // High velocity or negatively charged particle
1506  if( beta2 > 3.0*theBohrBeta2*zeff || charge < 0.0) {
1507  siga = std::sqrt( siga * chargeSquare ) ;
1508 
1509  // Low velocity - additional ion charge fluctuations according to
1510  // Q.Yang et al., NIM B61(1991)149-155.
1511  } else {
1512  G4double chu = theIonChuFluctuationModel->TheValue(particle, material);
1513  G4double yang = theIonYangFluctuationModel->TheValue(particle, material);
1514  siga = std::sqrt( siga * (chargeSquare * chu + yang)) ;
1515  }
1516 
1517  do {
1518  loss = G4RandGauss::shoot(meanLoss,siga);
1519  } while (loss < 0. || loss > 2.0*meanLoss);
1520  return loss;
1521  }
1522 
1523  // Non Gaussian fluctuation
1524  static const G4double probLim = 0.01 ;
1525  static const G4double sumaLim = -std::log(probLim) ;
1526  static const G4double alim = 10.;
1527 
1528  G4double suma,w1,w2,C,e0,lossc,w;
1529  G4double a1,a2,a3;
1530  G4int p1,p2,p3;
1531  G4int nb;
1532  G4double corrfac, na,alfa,rfac,namean,sa,alfa1,ea,sea;
1533  G4double dp3;
1534 
1535  G4double f1Fluct = material->GetIonisation()->GetF1fluct();
1536  G4double f2Fluct = material->GetIonisation()->GetF2fluct();
1537  G4double e1Fluct = material->GetIonisation()->GetEnergy1fluct();
1538  G4double e2Fluct = material->GetIonisation()->GetEnergy2fluct();
1539  G4double e1LogFluct = material->GetIonisation()->GetLogEnergy1fluct();
1540  G4double e2LogFluct = material->GetIonisation()->GetLogEnergy2fluct();
1541  G4double rateFluct = material->GetIonisation()->GetRateionexcfluct();
1542  G4double ipotLogFluct= material->GetIonisation()->GetLogMeanExcEnergy();
1543 
1544  w1 = tMax/ipotFluct;
1545  w2 = std::log(2.*electron_mass_c2*tau2);
1546 
1547  C = meanLoss*(1.-rateFluct)/(w2-ipotLogFluct-beta2);
1548 
1549  a1 = C*f1Fluct*(w2-e1LogFluct-beta2)/e1Fluct;
1550  a2 = C*f2Fluct*(w2-e2LogFluct-beta2)/e2Fluct;
1551  a3 = rateFluct*meanLoss*(tMax-ipotFluct)/(ipotFluct*tMax*std::log(w1));
1552  if(a1 < 0.0) a1 = 0.0;
1553  if(a2 < 0.0) a2 = 0.0;
1554  if(a3 < 0.0) a3 = 0.0;
1555 
1556  suma = a1+a2+a3;
1557 
1558  loss = 0.;
1559 
1560 
1561  if(suma < sumaLim) // very small Step
1562  {
1563  e0 = material->GetIonisation()->GetEnergy0fluct();
1564 
1565  if(tMax == ipotFluct)
1566  {
1567  a3 = meanLoss/e0;
1568 
1569  if(a3>alim)
1570  {
1571  siga=std::sqrt(a3) ;
1572  p3 = std::max(0,G4int(G4RandGauss::shoot(a3,siga)+0.5));
1573  }
1574  else
1575  p3 = G4Poisson(a3);
1576 
1577  loss = p3*e0 ;
1578 
1579  if(p3 > 0)
1580  loss += (1.-2.*G4UniformRand())*e0 ;
1581 
1582  }
1583  else
1584  {
1585  tMax = tMax-ipotFluct+e0 ;
1586  a3 = meanLoss*(tMax-e0)/(tMax*e0*std::log(tMax/e0));
1587 
1588  if(a3>alim)
1589  {
1590  siga=std::sqrt(a3) ;
1591  p3 = std::max(0,int(G4RandGauss::shoot(a3,siga)+0.5));
1592  }
1593  else
1594  p3 = G4Poisson(a3);
1595 
1596  if(p3 > 0)
1597  {
1598  w = (tMax-e0)/tMax ;
1599  if(p3 > nmaxCont2)
1600  {
1601  dp3 = G4float(p3) ;
1602  corrfac = dp3/G4float(nmaxCont2) ;
1603  p3 = nmaxCont2 ;
1604  }
1605  else
1606  corrfac = 1. ;
1607 
1608  for(G4int i=0; i<p3; i++) loss += 1./(1.-w*G4UniformRand()) ;
1609  loss *= e0*corrfac ;
1610  }
1611  }
1612  }
1613 
1614  else // not so small Step
1615  {
1616  // excitation type 1
1617  if(a1>alim)
1618  {
1619  siga=std::sqrt(a1) ;
1620  p1 = std::max(0,G4int(G4RandGauss::shoot(a1,siga)+0.5));
1621  }
1622  else
1623  p1 = G4Poisson(a1);
1624 
1625  // excitation type 2
1626  if(a2>alim)
1627  {
1628  siga=std::sqrt(a2) ;
1629  p2 = std::max(0,G4int(G4RandGauss::shoot(a2,siga)+0.5));
1630  }
1631  else
1632  p2 = G4Poisson(a2);
1633 
1634  loss = p1*e1Fluct+p2*e2Fluct;
1635 
1636  // smearing to avoid unphysical peaks
1637  if(p2 > 0)
1638  loss += (1.-2.*G4UniformRand())*e2Fluct;
1639  else if (loss>0.)
1640  loss += (1.-2.*G4UniformRand())*e1Fluct;
1641 
1642  // ionisation .......................................
1643  if(a3 > 0.)
1644  {
1645  if(a3>alim)
1646  {
1647  siga=std::sqrt(a3) ;
1648  p3 = std::max(0,G4int(G4RandGauss::shoot(a3,siga)+0.5));
1649  }
1650  else
1651  p3 = G4Poisson(a3);
1652 
1653  lossc = 0.;
1654  if(p3 > 0)
1655  {
1656  na = 0.;
1657  alfa = 1.;
1658  if (p3 > nmaxCont2)
1659  {
1660  dp3 = G4float(p3);
1661  rfac = dp3/(G4float(nmaxCont2)+dp3);
1662  namean = G4float(p3)*rfac;
1663  sa = G4float(nmaxCont1)*rfac;
1664  na = G4RandGauss::shoot(namean,sa);
1665  if (na > 0.)
1666  {
1667  alfa = w1*G4float(nmaxCont2+p3)/
1668  (w1*G4float(nmaxCont2)+G4float(p3));
1669  alfa1 = alfa*std::log(alfa)/(alfa-1.);
1670  ea = na*ipotFluct*alfa1;
1671  sea = ipotFluct*std::sqrt(na*(alfa-alfa1*alfa1));
1672  lossc += G4RandGauss::shoot(ea,sea);
1673  }
1674  }
1675 
1676  nb = G4int(G4float(p3)-na);
1677  if (nb > 0)
1678  {
1679  w2 = alfa*ipotFluct;
1680  w = (tMax-w2)/tMax;
1681  for (G4int k=0; k<nb; k++) lossc += w2/(1.-w*G4UniformRand());
1682  }
1683  }
1684  loss += lossc;
1685  }
1686  }
1687 
1688  return loss ;
1689 }
1690 
1691 
1692 
1694 {
1695  minGammaEnergy = cut;
1696 }
1697 
1698 
1699 
1701 {
1702  minElectronEnergy = cut;
1703 }
1704 
1705 
1706 
1708 {
1710 }
1711 
1712 
1713 
1715 {
1716  G4String comments = " Knock-on electron cross sections . ";
1717  comments += "\n Good description above the mean excitation energy.\n";
1718  comments += " Delta ray energy sampled from differential Xsection.";
1719 
1720  G4cout << G4endl << GetProcessName() << ": " << comments
1721  << "\n PhysicsTables from " << LowestKineticEnergy / eV << " eV "
1722  << " to " << HighestKineticEnergy / TeV << " TeV "
1723  << " in " << TotBin << " bins."
1724  << "\n Electronic stopping power model is "
1725  << protonTable
1726  << "\n from " << protonLowEnergy / keV << " keV "
1727  << " to " << protonHighEnergy / MeV << " MeV " << "." << G4endl ;
1728  G4cout << "\n Parametrisation model for antiprotons is "
1729  << antiprotonTable
1730  << "\n from " << antiprotonLowEnergy / keV << " keV "
1731  << " to " << antiprotonHighEnergy / MeV << " MeV " << "." << G4endl ;
1732  if(theBarkas){
1733  G4cout << " Parametrization of the Barkas effect is switched on."
1734  << G4endl ;
1735  }
1736  if(nStopping) {
1737  G4cout << " Nuclear stopping power model is " << theNuclearTable
1738  << G4endl ;
1739  }
1740 
1741  G4bool printHead = true;
1742 
1743  const G4ProductionCutsTable* theCoupleTable=
1745  size_t numOfCouples = theCoupleTable->GetTableSize();
1746 
1747  // loop for materials
1748 
1749  for (size_t j=0 ; j < numOfCouples; j++) {
1750 
1751  const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(j);
1752  const G4Material* material= couple->GetMaterial();
1753  G4double deltaCutNow = cutForDelta[(couple->GetIndex())] ;
1754  G4double excitationEnergy = material->GetIonisation()->GetMeanExcitationEnergy();
1755 
1756  if(excitationEnergy > deltaCutNow) {
1757  if(printHead) {
1758  printHead = false ;
1759 
1760  G4cout << " material min.delta energy(keV) " << G4endl;
1761  G4cout << G4endl;
1762  }
1763 
1764  G4cout << std::setw(20) << material->GetName()
1765  << std::setw(15) << excitationEnergy/keV << G4endl;
1766  }
1767  }
1768 }
G4double condition(const G4ErrorSymMatrix &m)
G4ParticleDefinition * GetDefinition() const
void BuildPhysicsTable(const G4ParticleDefinition &aParticleType)
void BuildLossTable(const G4ParticleDefinition &aParticleType)
ThreeVector shoot(const G4int Ap, const G4int Af)
G4IonisParamMat * GetIonisation() const
Definition: G4Material.hh:224
void SetProtonElectronicStoppingPowerModel(const G4String &dedxTable)
static G4ThreadLocal G4PhysicsTable * theProperTimepTable
static const double MeV
Definition: G4SIunits.hh:193
static G4ThreadLocal G4double LowestKineticEnergy
G4long G4Poisson(G4double mean)
Definition: G4Poisson.hh:51
static G4ThreadLocal G4PhysicsTable * theLabTimepTable
G4int ShellId() const
const std::vector< G4double > * GetEnergyCutsVector(size_t pcIdx) const
G4int verboseLevel
Definition: G4VProcess.hh:368
std::vector< G4Element * > G4ElementVector
void insert(G4PhysicsVector *)
G4double GetKineticEnergy() const
CLHEP::Hep3Vector G4ThreeVector
void ActivateAugerElectronProduction(G4bool val)
G4double GetStepLength() const
const G4DynamicParticle * GetDynamicParticle() const
G4double GetEnergy2fluct() const
static G4Material * GetMaterial(const G4String &name, G4bool warning=true)
Definition: G4Material.cc:603
G4double DeltaRaysEnergy(const G4MaterialCutsCouple *couple, G4double kineticEnergy, G4double particleMass) const
G4double GetConstraints(const G4DynamicParticle *particle, const G4MaterialCutsCouple *couple)
size_t GetIndex() const
Definition: G4Material.hh:260
G4double GetProductionCut(G4int index) const
static const G4double a1
const G4String & GetName() const
Definition: G4Material.hh:176
float G4float
Definition: G4Types.hh:77
G4VLowEnergyModel * theNuclearStoppingModel
static G4Proton * ProtonDefinition()
Definition: G4Proton.cc:88
const G4MaterialCutsCouple * GetMaterialCutsCouple() const
G4double GetLogEnergy2fluct() const
void SetMomentumDirection(const G4ThreeVector &aDirection)
G4double AntiProtonParametrisedDEDX(const G4MaterialCutsCouple *couple, G4double kineticEnergy) const
G4ParticleDefinition * GetDefinition() const
static G4ThreadLocal G4PhysicsTable ** RecorderOfpbarProcess
G4double BindingEnergy() const
G4double GetLogMeanExcEnergy() const
G4hImpactIonisation(const G4String &processName="hImpactIoni")
const G4String & GetParticleSubType() const
static G4ThreadLocal G4double finalRange
G4double GetLowEdgeEnergy(size_t binNumber) const
G4VLowEnergyModel * theIonEffChargeModel
#define G4ThreadLocal
Definition: tls.hh:84
G4ProcessManager * GetProcessManager() const
const G4ElementVector * GetElementVector() const
Definition: G4Material.hh:188
int G4int
Definition: G4Types.hh:78
const G4String & GetParticleName() const
const G4double paramStepLimit
G4VLowEnergyModel * theIonChuFluctuationModel
void ProposeLocalEnergyDeposit(G4double anEnergyPart)
G4double MicroscopicCrossSection(const G4ParticleDefinition &aParticleType, G4double kineticEnergy, G4double atomicNumber, G4double deltaCutInEnergy) const
virtual G4bool IsInCharge(const G4DynamicParticle *particle, const G4Material *material) const =0
G4VLowEnergyModel * theIonYangFluctuationModel
static G4ThreadLocal G4int TotBin
G4double GetEnergy0fluct() const
G4double MinKineticEnergy
static G4ThreadLocal G4PhysicsTable ** RecorderOfpProcess
G4ProcessVector * GetAtRestProcessVector(G4ProcessVectorTypeIndex typ=typeGPIL) const
void SetCutForSecondaryPhotons(G4double cut)
#define G4UniformRand()
Definition: Randomize.hh:95
G4GLOB_DLL std::ostream G4cout
static G4ThreadLocal G4PhysicsTable * theRangepTable
void BuildLambdaTable(const G4ParticleDefinition &aParticleType)
G4double GetElectronDensity() const
Definition: G4Material.hh:215
static G4AntiProton * AntiProton()
Definition: G4AntiProton.cc:93
G4double GetMass() const
G4PixeCrossSectionHandler * pixeCrossSectionHandler
G4AtomicDeexcitation atomicDeexcitation
void SetCutForAugerElectrons(G4double cut)
G4double BarkasTerm(const G4Material *material, G4double kineticEnergy) const
bool G4bool
Definition: G4Types.hh:79
G4double ElectronicLossFluctuation(const G4DynamicParticle *particle, const G4MaterialCutsCouple *material, G4double meanLoss, G4double step) const
const G4ThreeVector & GetMomentumDirection() const
virtual G4double HighEnergyLimit(const G4ParticleDefinition *aParticle, const G4Material *material) const =0
G4double ProtonParametrisedDEDX(const G4MaterialCutsCouple *couple, G4double kineticEnergy) const
G4double GetCharge() const
void PutValue(size_t index, G4double theValue)
void SetElectronicStoppingPowerModel(const G4ParticleDefinition *aParticle, const G4String &dedxTable)
static G4Proton * Proton()
Definition: G4Proton.cc:93
G4int SelectRandomAtom(const G4Material *material, G4double e) const
static const double GeV
Definition: G4SIunits.hh:196
const G4String & GetParticleType() const
static G4ThreadLocal G4double dRoverRange
Definition: G4Step.hh:76
void ActivateAugerElectronProduction(G4bool val)
void LoadShellData(const G4String &dataFile)
static G4ThreadLocal G4int CounterOfpProcess
static G4Gamma * Gamma()
Definition: G4Gamma.cc:86
G4PhysicsTable * theMeanFreePathTable
const G4String & GetProcessName() const
Definition: G4VProcess.hh:408
const G4double * GetAtomicNumDensityVector() const
Definition: G4Material.hh:214
void SetKineticEnergy(G4double aEnergy)
static const G4double a3
static void BuildDEDXTable(const G4ParticleDefinition &aParticleType)
static G4ThreadLocal G4double HighestKineticEnergy
G4int size() const
G4double GetTotNbOfAtomsPerVolume() const
Definition: G4Material.hh:207
virtual void Initialize(const G4Track &)
static G4ProductionCutsTable * GetProductionCutsTable()
static const double eV
Definition: G4SIunits.hh:194
static G4ThreadLocal G4int CounterOfpbarProcess
G4int SelectRandomShell(G4int Z, G4double e) const
G4double GetLogEnergy1fluct() const
G4VLowEnergyModel * protonModel
G4double GetRateionexcfluct() const
G4double GetPDGMass() const
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
T max(const T t1, const T t2)
brief Return the largest of the two arguments
G4double energy(const ThreeVector &p, const G4double m)
void SetNumberOfSecondaries(G4int totSecondaries)
G4bool CutsWhereModified()
static G4ThreadLocal G4bool EnlossFlucFlag
void ProposeEnergy(G4double finalEnergy)
G4ParticleChange aParticleChange
Definition: G4VProcess.hh:289
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
void AddSecondary(G4Track *aSecondary)
virtual G4double TheValue(const G4DynamicParticle *particle, const G4Material *material)=0
G4VLowEnergyModel * antiprotonModel
G4double GetPDGSpin() const
G4double GetMeanExcitationEnergy() const
G4PhysicsTable * theLossTable
static G4Electron * Electron()
Definition: G4Electron.cc:94
#define G4endl
Definition: G4ios.hh:61
static const double TeV
Definition: G4SIunits.hh:197
G4double GetF2fluct() const
void SetAntiProtonElectronicStoppingPowerModel(const G4String &dedxTable)
size_t GetNumberOfElements() const
Definition: G4Material.hh:184
std::vector< G4DynamicParticle * > * GenerateParticles(G4int Z, G4int shellId)
static const double keV
Definition: G4SIunits.hh:195
virtual G4VParticleChange * PostStepDoIt(const G4Track &, const G4Step &)
static G4double GetDEDX(const G4ParticleDefinition *aParticle, G4double KineticEnergy, const G4Material *aMaterial)
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
double G4double
Definition: G4Types.hh:76
void ProposeTrackStatus(G4TrackStatus status)
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)
void SetDefinition(const G4ParticleDefinition *aParticleDefinition)
G4double GetMeanFreePath(const G4Track &track, G4double previousStepSize, enum G4ForceCondition *condition)
static const double eplus
Definition: G4SIunits.hh:178
G4ForceCondition
G4double BlochTerm(const G4Material *material, G4double kineticEnergy, G4double cSquare) const
G4double GetPDGCharge() const
static G4double GetPreciseEnergyFromRange(const G4ParticleDefinition *aParticle, G4double range, const G4Material *aMaterial)
G4double bindingEnergy(G4int A, G4int Z)
static G4ThreadLocal G4PhysicsTable * theInverseRangepTable
static G4AtomicTransitionManager * Instance()
G4ProductionCuts * GetProductionCuts() const
G4double ComputeDEDX(const G4ParticleDefinition *aParticle, const G4MaterialCutsCouple *couple, G4double kineticEnergy)
#define DBL_MAX
Definition: templates.hh:83
G4VParticleChange * PostStepDoIt(const G4Track &track, const G4Step &Step)
static const double mm
Definition: G4SIunits.hh:102
G4AtomicShell * Shell(G4int Z, size_t shellIndex) const
static G4double GetRange(const G4ParticleDefinition *aParticle, G4double KineticEnergy, const G4Material *aMaterial)
G4VParticleChange * AlongStepDoIt(const G4Track &trackData, const G4Step &stepData)
G4double GetF1fluct() const
static const G4double a2
void clearAndDestroy()
static G4ThreadLocal G4PhysicsTable * theDEDXpTable
G4ProcessVector * GetProcessList() const
G4VLowEnergyModel * betheBlochModel
const G4Material * GetMaterial() const
G4double GetEnergy1fluct() const
void PrintInfoDefinition() const