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