| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 //
27
28 // -------------------------------------------
29 // G4RDHadronIonisation
30 //
31 //
32 // Author: Maria Grazia Pia (MariaGrazia.Pia@g
33 //
34 // 08 Sep 2008 - MGP - Created (initially base
35 // Added PIXE capabilities
36 // Partial clean-up of the
37 // Calculation of Microsco
38 // M.G. Pia et al., PIXE Simulation With Geant
39 // IEEE Trans. Nucl. Sci., vol. 56, no. 6, pp.
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"
57 #include "G4AtomicTransitionManager.hh"
58 #include "G4PixeCrossSectionHandler.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"
70 #include "G4hNuclearStoppingModel.hh"
71 #include "G4hBetheBlochModel.hh"
72 #include "G4hParametrisedLossModel.hh"
73 #include "G4QAOLowEnergyLoss.hh"
74 #include "G4hIonEffChargeSquare.hh"
75 #include "G4IonChuFluctuationModel.hh"
76 #include "G4IonYangFluctuationModel.hh"
77
78 #include "G4MaterialCutsCouple.hh"
79 #include "G4Track.hh"
80 #include "G4Step.hh"
81
82 G4hImpactIonisation::G4hImpactIonisation(const
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
105 void G4hImpactIonisation::InitializeMe()
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 ;
114 antiprotonHighEnergy = 2.*MeV ;
115 minGammaEnergy = 100 * eV;
116 minElectronEnergy = 250.* eV;
117 verboseLevel = 0;
118
119 // Min and max energy of incident particle f
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 // Calculated on the fly, but should be init
130 fdEdx = 0.0;
131 fRangeNow = 0.0;
132 charge = 0.0;
133 chargeSquare = 0.0;
134 initialMass = 0.0;
135 fBarkas = 0.0;
136 }
137
138
139 G4hImpactIonisation::~G4hImpactIonisation()
140 {
141 if (theMeanFreePathTable)
142 {
143 theMeanFreePathTable->clearAndDestroy();
144 delete theMeanFreePathTable;
145 }
146
147 if (betheBlochModel) delete betheBlochModel;
148 if (protonModel) delete protonModel;
149 if (antiprotonModel) delete antiprotonModel;
150 if (theNuclearStoppingModel) delete theNucle
151 if (theIonEffChargeModel) delete theIonEffCh
152 if (theIonChuFluctuationModel) delete theIon
153 if (theIonYangFluctuationModel) delete theIo
154
155 delete pixeCrossSectionHandler;
156
157 // ---- MGP ---- The following is to be chec
158 // if (shellVacancy) delete shellVacancy;
159
160 cutForDelta.clear();
161 }
162
163 // -------------------------------------------
164 void G4hImpactIonisation::SetElectronicStoppin
165 const G4String& dedxTable)
166 // This method defines the ionisation parametr
167 {
168 if (particle->GetPDGCharge() > 0 )
169 {
170 // Positive charge
171 SetProtonElectronicStoppingPowerModel(de
172 }
173 else
174 {
175 // Antiprotons
176 SetAntiProtonElectronicStoppingPowerMode
177 }
178 }
179
180
181 // -------------------------------------------
182 void G4hImpactIonisation::InitializeParametris
183
184 {
185 // Define models for parametrisation of elec
186 betheBlochModel = new G4hBetheBlochModel("Be
187 protonModel = new G4hParametrisedLossModel(p
188 protonHighEnergy = std::min(protonHighEnergy
189 antiprotonModel = new G4QAOLowEnergyLoss(ant
190 theNuclearStoppingModel = new G4hNuclearStop
191 theIonEffChargeModel = new G4hIonEffChargeSq
192 theIonChuFluctuationModel = new G4IonChuFluc
193 theIonYangFluctuationModel = new G4IonYangFl
194 }
195
196
197 // -------------------------------------------
198 void G4hImpactIonisation::BuildPhysicsTable(co
199
200 // just call BuildLossTable+BuildLambdaTable
201 {
202
203 // Verbose print-out
204 if(verboseLevel > 0)
205 {
206 G4cout << "G4hImpactIonisation::BuildPhy
207 << particleDef.GetParticleName()
208 << " mass(MeV)= " << particleDef.GetPDG
209 << " charge= " << particleDef.GetPDGCha
210 << " type= " << particleDef.GetParticle
211 << G4endl;
212
213 if(verboseLevel > 1)
214 {
215 G4ProcessVector* pv = particleDef.GetProce
216
217 G4cout << " 0: " << (*pv)[0]->GetProcessNa
218 << " 1: " << (*pv)[1]->GetProcessName() <
219 // << " 2: " << (*pv)[2]->GetProc
220 << G4endl;
221 G4cout << "ionModel= " << theIonEffChargeM
222 << " MFPtable= " << theMeanFreePathTable
223 << " iniMass= " << initialMass
224 << G4endl;
225 }
226 }
227 // End of verbose print-out
228
229 if (particleDef.GetParticleType() == "nucleu
230 particleDef.GetParticleName() != "Generi
231 particleDef.GetParticleSubType() == "gen
232 {
233
234 G4EnergyLossTables::Register(&particleDe
235 theDEDXpTable,
236 theRangepTable,
237 theInverseRangepTable,
238 theLabTimepTable,
239 theProperTimepTable,
240 LowestKineticEnergy, HighestKinetic
241 proton_mass_c2/particleDef.GetPDGMa
242 TotBin);
243
244 return;
245 }
246
247 if( !CutsWhereModified() && theLossTable) re
248
249 InitializeParametrisation() ;
250 G4Proton* proton = G4Proton::Proton();
251 G4AntiProton* antiproton = G4AntiProton::Ant
252
253 charge = particleDef.GetPDGCharge() / eplus;
254 chargeSquare = charge*charge ;
255
256 const G4ProductionCutsTable* theCoupleTable=
257 G4ProductionCutsTable::GetProductionCutsTa
258 G4int numOfCouples = (G4int)theCoupleTable->
259
260 cutForDelta.clear();
261 cutForGamma.clear();
262
263 for (G4int j=0; j<numOfCouples; ++j) {
264
265 // get material parameters needed for the
266 const G4MaterialCutsCouple* couple = theCo
267 const G4Material* material= couple->GetMat
268
269 // the cut cannot be below lowest limit
270 G4double tCut = (*(theCoupleTable->GetEner
271 if(tCut > HighestKineticEnergy) tCut = Hig
272
273 G4double excEnergy = material->GetIonisati
274
275 tCut = std::max(tCut,excEnergy);
276 cutForDelta.push_back(tCut);
277
278 // the cut cannot be below lowest limit
279 tCut = (*(theCoupleTable->GetEnergyCutsVec
280 if(tCut > HighestKineticEnergy) tCut = Hig
281 tCut = std::max(tCut,minGammaEnergy);
282 cutForGamma.push_back(tCut);
283 }
284
285 if(verboseLevel > 0) {
286 G4cout << "Cuts are defined " << G4endl;
287 }
288
289 if(0.0 < charge)
290 {
291 {
292 BuildLossTable(*proton) ;
293
294 // The following vector has a fixed dim
295 // It happended in the past that caused
296 // G4cout << "[NOTE]: __LINE__=" << _
297
298 RecorderOfpProcess[CounterOfpProcess]
299 CounterOfpProcess++;
300 }
301 } else {
302 {
303 BuildLossTable(*antiproton) ;
304
305 // The following vector has a fixed
306 // It happended in the past that ca
307 // G4cout << "[NOTE]: __LINE__="
308
309 RecorderOfpbarProcess[CounterOfpbarProce
310 CounterOfpbarProcess++;
311 }
312 }
313
314 if(verboseLevel > 0) {
315 G4cout << "G4hImpactIonisation::BuildPhysi
316 << "Loss table is built "
317 // << theLossTable
318 << G4endl;
319 }
320
321 BuildLambdaTable(particleDef) ;
322 // BuildDataForFluorescence(particleDef);
323
324 if(verboseLevel > 1) {
325 G4cout << (*theMeanFreePathTable) << G4end
326 }
327
328 if(verboseLevel > 0) {
329 G4cout << "G4hImpactIonisation::BuildPhysi
330 << "DEDX table will be built "
331 // << theDEDXpTable << " " << theDED
332 // << " " << theRangepTable << " " <
333 << G4endl;
334 }
335
336 BuildDEDXTable(particleDef) ;
337
338 if(verboseLevel > 1) {
339 G4cout << (*theDEDXpTable) << G4endl;
340 }
341
342 if((&particleDef == proton) || (&particleDe
343
344 if(verboseLevel > 0) {
345 G4cout << "G4hImpactIonisation::BuildPhysi
346 << particleDef.GetParticleName() <<
347 }
348 }
349
350
351
352
353
354 // -------------------------------------------
355 void G4hImpactIonisation::BuildLossTable(const
356 {
357 // Initialisation
358 G4double lowEdgeEnergy , ionloss, ionlossBB,
359 //G4double lowEnergy, highEnergy;
360 G4double highEnergy;
361 G4Proton* proton = G4Proton::Proton();
362
363 if (particleDef == *proton)
364 {
365 //lowEnergy = protonLowEnergy ;
366 highEnergy = protonHighEnergy ;
367 charge = 1. ;
368 }
369 else
370 {
371 //lowEnergy = antiprotonLowEnergy ;
372 highEnergy = antiprotonHighEnergy ;
373 charge = -1. ;
374 }
375 chargeSquare = 1. ;
376
377 const G4ProductionCutsTable* theCoupleTable=
378 G4ProductionCutsTable::GetProductionCutsTa
379 G4int numOfCouples = (G4int)theCoupleTable->
380
381 if ( theLossTable)
382 {
383 theLossTable->clearAndDestroy();
384 delete theLossTable;
385 }
386
387 theLossTable = new G4PhysicsTable(numOfCoupl
388
389 // loop for materials
390 for (G4int j=0; j<numOfCouples; ++j) {
391
392 // create physics vector and fill it
393 G4PhysicsLogVector* aVector = new G4Physic
394
395
396
397 // get material parameters needed for the
398 const G4MaterialCutsCouple* couple = theCo
399 const G4Material* material= couple->GetMat
400
401 if ( charge > 0.0 ) {
402 ionloss = ProtonParametrisedDEDX(couple,
403 } else {
404 ionloss = AntiProtonParametrisedDEDX(cou
405 }
406
407 ionlossBB = betheBlochModel->TheValue(&par
408 ionlossBB -= DeltaRaysEnergy(couple,highEn
409
410
411 paramB = ionloss/ionlossBB - 1.0 ;
412
413 // now comes the loop for the kinetic ener
414 for (G4int i = 0 ; i < TotBin ; i++) {
415 lowEdgeEnergy = aVector->GetLowEdgeEnerg
416
417 // low energy part for this material, pa
418 if ( lowEdgeEnergy < highEnergy ) {
419
420 if ( charge > 0.0 ) {
421 ionloss = ProtonParametrisedDEDX(cou
422 } else {
423 ionloss = AntiProtonParametrisedDEDX
424 }
425
426 } else {
427
428 // high energy part for this material,
429 ionloss = betheBlochModel->TheValue(pr
430 lowEdgeEnergy) ;
431
432 ionloss -= DeltaRaysEnergy(couple,lowE
433
434 ionloss *= (1.0 + paramB*highEnergy/lowEdgeE
435 }
436
437 // now put the loss into the vector
438 if(verboseLevel > 1) {
439 G4cout << "E(MeV)= " << lowEdgeEnergy/
440 << " dE/dx(MeV/mm)= " << ionlo
441 << " in " << material->GetName(
442 }
443 aVector->PutValue(i,ionloss) ;
444 }
445 // Insert vector for this material into th
446 theLossTable->insert(aVector) ;
447 }
448 }
449
450
451
452 // -------------------------------------------
453 void G4hImpactIonisation::BuildLambdaTable(con
454
455 {
456 // Build mean free path tables for the delta
457 // tables are built for MATERIALS
458
459 if(verboseLevel > 1) {
460 G4cout << "G4hImpactIonisation::BuildLambd
461 << particleDef.GetParticleName() <<
462 }
463
464
465 G4double lowEdgeEnergy, value;
466 charge = particleDef.GetPDGCharge()/eplus ;
467 chargeSquare = charge*charge ;
468 initialMass = particleDef.GetPDGMass();
469
470 const G4ProductionCutsTable* theCoupleTable=
471 G4ProductionCutsTable::GetProductionCutsTa
472 G4int numOfCouples = (G4int)theCoupleTable->
473
474
475 if (theMeanFreePathTable) {
476 theMeanFreePathTable->clearAndDestroy();
477 delete theMeanFreePathTable;
478 }
479
480 theMeanFreePathTable = new G4PhysicsTable(nu
481
482 // loop for materials
483
484 for (G4int j=0 ; j < numOfCouples; ++j) {
485
486 //create physics vector then fill it ....
487 G4PhysicsLogVector* aVector = new G4Physic
488
489
490
491 // compute the (macroscopic) cross section
492 const G4MaterialCutsCouple* couple = theCo
493 const G4Material* material= couple->GetMat
494
495 const G4ElementVector* theElementVector =
496 const G4double* theAtomicNumDensityVector
497 const G4int numberOfElements = (G4int)mate
498
499 // get the electron kinetic energy cut for
500 // it will be used in ComputeMicroscopicC
501 // ( it is the SAME for ALL the ELEMENTS i
502 // -------------------------------------
503
504 G4double deltaCut = cutForDelta[j];
505
506 for ( G4int i = 0 ; i < TotBin ; i++ ) {
507 lowEdgeEnergy = aVector->GetLowEdgeEnerg
508 G4double sigma = 0.0 ;
509 G4int Z;
510
511 for (G4int iel=0; iel<numberOfElements;
512 {
513 Z = (G4int) (*theElementVector)[iel]->GetZ
514 // ---- MGP --- Corrected duplicated cross
515 // G4hLowEnergyIonisation original
516 G4double microCross = MicroscopicCrossSect
517 lowEdgeEnergy,
518 Z,
519 deltaCut ) ;
520 //totalCrossSectionMap [Z] = microCross;
521 sigma += theAtomicNumDensityVector[iel] *
522 }
523
524 // mean free path = 1./macroscopic cross
525
526 value = sigma<=0 ? DBL_MAX : 1./sigma ;
527
528 aVector->PutValue(i, value) ;
529 }
530
531 theMeanFreePathTable->insert(aVector);
532 }
533
534 }
535
536
537 // -------------------------------------------
538 G4double G4hImpactIonisation::MicroscopicCross
539 G4double kineticEnergy,
540 G4double atomicNumber,
541 G4double deltaCutInEnergy) c
542 {
543 //******************************************
544 // cross section formula is OK for spin=0,
545 // ***************************************
546
547 // Calculates the microscopic cross sectio
548 // Formula documented in Geant4 Phys. Ref.
549 // ( it is called for elements, AtomicNumb
550
551 G4double totalCrossSection = 0.;
552
553 // Particle mass and energy
554 G4double particleMass = initialMass;
555 G4double energy = kineticEnergy + particleMa
556
557 // Some kinematics
558 G4double gamma = energy / particleMass;
559 G4double beta2 = 1. - 1. / (gamma * gamma);
560 G4double var = electron_mass_c2 / particleMa
561 G4double tMax = 2. * electron_mass_c2 * (g
562
563 // Calculate the total cross section
564
565 if ( tMax > deltaCutInEnergy )
566 {
567 var = deltaCutInEnergy / tMax;
568 totalCrossSection = (1. - var * (1. - be
569
570 G4double spin = particleDef.GetPDGSpin()
571
572 // +term for spin=1/2 particle
573 if (spin == 0.5)
574 {
575 totalCrossSection += 0.5 * (tMax - deltaC
576 }
577 // +term for spin=1 particle
578 else if (spin > 0.9 )
579 {
580 totalCrossSection += -std::log(var) /
581 (3. * deltaCutInEnergy) + (tMax - deltaC
582 beta2 / (tMax * deltaCutIn
583 }
584 totalCrossSection *= twopi_mc2_rcl2 * at
585 }
586
587 //std::cout << "Microscopic = " << totalCros
588 // << ", e = " << kineticEnergy/MeV <<std
589
590 return totalCrossSection ;
591 }
592
593
594
595 // -------------------------------------------
596 G4double G4hImpactIonisation::GetMeanFreePath(
597 G4double, // previousStepSize
598 enum G4ForceCondition* conditi
599 {
600 const G4DynamicParticle* dynamicParticle = t
601 const G4MaterialCutsCouple* couple = track.G
602 const G4Material* material = couple->GetMate
603
604 G4double meanFreePath = DBL_MAX;
605 // ---- MGP ---- What is the meaning of the
606 G4bool isOutRange = false;
607
608 *condition = NotForced ;
609
610 G4double kineticEnergy = (dynamicParticle->G
611 charge = dynamicParticle->GetCharge()/eplus;
612 chargeSquare = theIonEffChargeModel->TheValu
613
614 if (kineticEnergy < LowestKineticEnergy)
615 {
616 meanFreePath = DBL_MAX;
617 }
618 else
619 {
620 if (kineticEnergy > HighestKineticEnergy
621 meanFreePath = (((*theMeanFreePathTable)
622 GetValue(kineticEnergy,isOutRange))/
623 }
624
625 return meanFreePath ;
626 }
627
628
629 // -------------------------------------------
630 G4double G4hImpactIonisation::GetConstraints(c
631 const G4MaterialCutsCouple* cou
632 {
633 // returns the Step limit
634 // dEdx is calculated as well as the range
635 // based on Effective Charge Approach
636
637 const G4Material* material = couple->GetMate
638 G4Proton* proton = G4Proton::Proton();
639 G4AntiProton* antiproton = G4AntiProton::Ant
640
641 G4double stepLimit = 0.;
642 G4double dx, highEnergy;
643
644 G4double massRatio = proton_mass_c2/(particl
645 G4double kineticEnergy = particle->GetKineti
646
647 // Scale the kinetic energy
648
649 G4double tScaled = kineticEnergy*massRatio ;
650 fBarkas = 0.;
651
652 if (charge > 0.)
653 {
654 highEnergy = protonHighEnergy ;
655 fRangeNow = G4EnergyLossTables::GetRange
656 dx = G4EnergyLossTables::GetRange(proton
657 fdEdx = G4EnergyLossTables::GetDEDX(prot
658 * chargeSquare ;
659
660 // Correction for positive ions
661 if (theBarkas && tScaled > highEnergy)
662 {
663 fBarkas = BarkasTerm(material,tScaled)*std
664 + BlochTerm(material,tScaled,chargeSquar
665 }
666 // Antiprotons and negative hadrons
667 }
668 else
669 {
670 highEnergy = antiprotonHighEnergy ;
671 fRangeNow = G4EnergyLossTables::GetRange
672 dx = G4EnergyLossTables::GetRange(antipr
673 fdEdx = G4EnergyLossTables::GetDEDX(anti
674
675 if (theBarkas && tScaled > highEnergy)
676 {
677 fBarkas = -BarkasTerm(material,tScaled)*st
678 + BlochTerm(material,tScaled,chargeSquar
679 }
680 }
681 /*
682 const G4Material* mat = couple->GetMateria
683 G4double fac = gram/(MeV*cm2*mat->GetDensi
684 G4cout << particle->GetDefinition()->GetPa
685 << " in " << mat->GetName()
686 << " E(MeV)= " << kineticEnergy/MeV
687 << " dedx(MeV*cm^2/g)= " << fdEdx*fac
688 << " barcas(MeV*cm^2/gram)= " << fBarkas*f
689 << " Q^2= " << chargeSquare
690 << G4endl;
691 */
692 // scaling back
693 fRangeNow /= (chargeSquare*massRatio) ;
694 dx /= (chargeSquare*massRatio) ;
695
696 stepLimit = fRangeNow ;
697 G4double r = std::min(finalRange, couple->Ge
698 ->GetProductionCut(idxG4ElectronCut));
699
700 if (fRangeNow > r)
701 {
702 stepLimit = dRoverRange*fRangeNow + r*(1
703 if (stepLimit > fRangeNow) stepLimit = f
704 }
705 // compute the (random) Step limit in stan
706 if(tScaled > highEnergy )
707 {
708 // add Barkas correction directly to ded
709 fdEdx += fBarkas;
710
711 if(stepLimit > fRangeNow - dx*0.9) stepL
712
713 // Step limit in low energy range
714 }
715 else
716 {
717 G4double x = dx*paramStepLimit;
718 if (stepLimit > x) stepLimit = x;
719 }
720 return stepLimit;
721 }
722
723
724 // -------------------------------------------
725 G4VParticleChange* G4hImpactIonisation::AlongS
726 const G4Step& step)
727 {
728 // compute the energy loss after a step
729 G4Proton* proton = G4Proton::Proton();
730 G4AntiProton* antiproton = G4AntiProton::Ant
731 G4double finalT = 0.;
732
733 aParticleChange.Initialize(track) ;
734
735 const G4MaterialCutsCouple* couple = track.G
736 const G4Material* material = couple->GetMate
737
738 // get the actual (true) Step length from st
739 const G4double stepLength = step.GetStepLeng
740
741 const G4DynamicParticle* particle = track.Ge
742
743 G4double kineticEnergy = particle->GetKineti
744 G4double massRatio = proton_mass_c2/(particl
745 G4double tScaled = kineticEnergy * massRatio
746 G4double eLoss = 0.0 ;
747 G4double nLoss = 0.0 ;
748
749
750 // very small particle energy
751 if (kineticEnergy < MinKineticEnergy)
752 {
753 eLoss = kineticEnergy ;
754 // particle energy outside tabulated ene
755 }
756
757 else if( kineticEnergy > HighestKineticEnerg
758 {
759 eLoss = stepLength * fdEdx ;
760 // big step
761 }
762 else if (stepLength >= fRangeNow )
763 {
764 eLoss = kineticEnergy ;
765
766 // tabulated range
767 }
768 else
769 {
770 // step longer than linear step limit
771 if(stepLength > linLossLimit * fRangeNow
772 {
773 G4double rScaled = fRangeNow * massRatio *
774 G4double sScaled = stepLength * massRatio
775
776 if(charge > 0.0)
777 {
778 eLoss = G4EnergyLossTables::GetPrecise
779 G4EnergyLossTables::GetPreciseEnergyFromRa
780
781 }
782 else
783 {
784 // Antiproton
785 eLoss = G4EnergyLossTables::GetPrecise
786 G4EnergyLossTables::GetPreciseEnergyFromRa
787 }
788 eLoss /= massRatio ;
789
790 // Barkas correction at big step
791 eLoss += fBarkas * stepLength;
792
793 // step shorter than linear step limit
794 }
795 else
796 {
797 eLoss = stepLength *fdEdx ;
798 }
799 if (nStopping && tScaled < protonHighEne
800 {
801 nLoss = (theNuclearStoppingModel->TheValue
802 }
803 }
804
805 if (eLoss < 0.0) eLoss = 0.0;
806
807 finalT = kineticEnergy - eLoss - nLoss;
808
809 if ( EnlossFlucFlag && 0.0 < eLoss && finalT
810 {
811
812 // now the electron loss with fluctuati
813 eLoss = ElectronicLossFluctuation(partic
814 if (eLoss < 0.0) eLoss = 0.0;
815 finalT = kineticEnergy - eLoss - nLoss;
816 }
817
818 // stop particle if the kinetic energy <= M
819 if (finalT*massRatio <= MinKineticEnergy )
820 {
821
822 finalT = 0.0;
823 if (!particle->GetDefinition()->GetProce
824 aParticleChange.ProposeTrackStatus(fSt
825 else
826 aParticleChange.ProposeTrackStatus(fSt
827 }
828
829 aParticleChange.ProposeEnergy( finalT );
830 eLoss = kineticEnergy-finalT;
831
832 aParticleChange.ProposeLocalEnergyDeposit(eL
833 return &aParticleChange ;
834 }
835
836
837
838 // -------------------------------------------
839 G4double G4hImpactIonisation::ProtonParametris
840 G4double kineticEnergy) const
841 {
842 const G4Material* material = couple->GetMate
843 G4Proton* proton = G4Proton::Proton();
844 G4double eLoss = 0.;
845
846 // Free Electron Gas Model
847 if(kineticEnergy < protonLowEnergy) {
848 eLoss = (protonModel->TheValue(proton, mat
849 * std::sqrt(kineticEnergy/protonLowEnerg
850
851 // Parametrisation
852 } else {
853 eLoss = protonModel->TheValue(proton, mate
854 }
855
856 // Delta rays energy
857 eLoss -= DeltaRaysEnergy(couple,kineticEnerg
858
859 if(verboseLevel > 2) {
860 G4cout << "p E(MeV)= " << kineticEnergy/Me
861 << " dE/dx(MeV/mm)= " << eLoss*mm/M
862 << " for " << material->GetName()
863 << " model: " << protonModel << G4e
864 }
865
866 if(eLoss < 0.0) eLoss = 0.0 ;
867
868 return eLoss ;
869 }
870
871
872
873 // -------------------------------------------
874 G4double G4hImpactIonisation::AntiProtonParame
875 G4double kineticEnergy) const
876 {
877 const G4Material* material = couple->GetMate
878 G4AntiProton* antiproton = G4AntiProton::Ant
879 G4double eLoss = 0.0 ;
880
881 // Antiproton model is used
882 if(antiprotonModel->IsInCharge(antiproton,ma
883 if(kineticEnergy < antiprotonLowEnergy) {
884 eLoss = antiprotonModel->TheValue(antipr
885 * std::sqrt(kineticEnergy/antiprotonLowEnerg
886
887 // Parametrisation
888 } else {
889 eLoss = antiprotonModel->TheValue(antipr
890 kineticEnergy);
891 }
892
893 // The proton model is used + Barkas corre
894 } else {
895 if(kineticEnergy < protonLowEnergy) {
896 eLoss = protonModel->TheValue(G4Proton::
897 * std::sqrt(kineticEnergy/protonLowEnergy) ;
898
899 // Parametrisation
900 } else {
901 eLoss = protonModel->TheValue(G4Proton::
902 kineticEnergy);
903 }
904 //if(theBarkas) eLoss -= 2.0*BarkasTerm(ma
905 }
906
907 // Delta rays energy
908 eLoss -= DeltaRaysEnergy(couple,kineticEnerg
909
910 if(verboseLevel > 2) {
911 G4cout << "pbar E(MeV)= " << kineticEnergy
912 << " dE/dx(MeV/mm)= " << eLoss*mm/M
913 << " for " << material->GetName()
914 << " model: " << protonModel << G4e
915 }
916
917 if(eLoss < 0.0) eLoss = 0.0 ;
918
919 return eLoss ;
920 }
921
922
923 // -------------------------------------------
924 G4double G4hImpactIonisation::DeltaRaysEnergy(
925 G4double kineticEnergy,
926 G4double particleMass) const
927 {
928 G4double dLoss = 0.;
929
930 G4double deltaCutNow = cutForDelta[(couple->
931 const G4Material* material = couple->GetMate
932 G4double electronDensity = material->GetElec
933 G4double excitationEnergy = material->GetIon
934
935 G4double tau = kineticEnergy / particleMass
936 G4double rateMass = electron_mass_c2/particl
937
938 // some local variables
939
940 G4double gamma = tau + 1.0 ;
941 G4double bg2 = tau*(tau+2.0) ;
942 G4double beta2 = bg2/(gamma*gamma) ;
943 G4double tMax = 2.*electron_mass_c2*bg2/(1.0
944
945 // Validity range for delta electron cross s
946 G4double deltaCut = std::max(deltaCutNow, ex
947
948 if ( deltaCut < tMax)
949 {
950 G4double x = deltaCut / tMax ;
951 dLoss = ( beta2 * (x-1.) - std::log(x) )
952 }
953 return dLoss ;
954 }
955
956
957 // -------------------------------------------
958
959 G4VParticleChange* G4hImpactIonisation::PostSt
960 const G4Step& step)
961 {
962 // Units are expressed in GEANT4 internal un
963
964 // std::cout << "----- Calling PostStepDoIt
965
966 aParticleChange.Initialize(track) ;
967 const G4MaterialCutsCouple* couple = track.G
968 const G4DynamicParticle* aParticle = track.G
969
970 // Some kinematics
971
972 G4ParticleDefinition* definition = track.Get
973 G4double mass = definition->GetPDGMass();
974 G4double kineticEnergy = aParticle->GetKinet
975 G4double totalEnergy = kineticEnergy + mass
976 G4double pSquare = kineticEnergy *( totalEne
977 G4double eSquare = totalEnergy * totalEnergy
978 G4double betaSquare = pSquare / eSquare;
979 G4ThreeVector particleDirection = aParticle-
980
981 G4double gamma = kineticEnergy / mass + 1.;
982 G4double r = electron_mass_c2 / mass;
983 G4double tMax = 2. * electron_mass_c2 *(gamm
984
985 // Validity range for delta electron cross s
986 G4double deltaCut = cutForDelta[couple->GetI
987
988 // This should not be a case
989 if (deltaCut >= tMax)
990 return G4VContinuousDiscreteProcess::PostS
991
992 G4double xc = deltaCut / tMax;
993 G4double rate = tMax / totalEnergy;
994 rate = rate*rate ;
995 G4double spin = aParticle->GetDefinition()->
996
997 // Sampling follows ...
998 G4double x = 0.;
999 G4double gRej = 0.;
1000
1001 do {
1002 x = xc / (1. - (1. - xc) * G4UniformRand(
1003
1004 if (0.0 == spin)
1005 {
1006 gRej = 1.0 - betaSquare * x ;
1007 }
1008 else if (0.5 == spin)
1009 {
1010 gRej = (1. - betaSquare * x + 0.5 * x*x * r
1011 }
1012 else
1013 {
1014 gRej = (1. - betaSquare * x ) * (1. + x/(3.
1015 x*x * rate * (1. + 0.5*x/xc) / 3.0 /
1016 (1. + 1./(3.*xc) + rate *(1.+ 0.5/xc) / 3
1017 }
1018
1019 } while( G4UniformRand() > gRej );
1020
1021 G4double deltaKineticEnergy = x * tMax;
1022 G4double deltaTotalMomentum = std::sqrt(del
1023 (deltaKineticEnergy + 2. * electr
1024 G4double totalMomentum = std::sqrt(pSquare)
1025 G4double cosTheta = deltaKineticEnergy * (t
1026
1027 // protection against cosTheta > 1 or < -1
1028 if ( cosTheta < -1. ) cosTheta = -1.;
1029 if ( cosTheta > 1. ) cosTheta = 1.;
1030
1031 // direction of the delta electron
1032 G4double phi = twopi * G4UniformRand();
1033 G4double sinTheta = std::sqrt(1. - cosTheta
1034 G4double dirX = sinTheta * std::cos(phi);
1035 G4double dirY = sinTheta * std::sin(phi);
1036 G4double dirZ = cosTheta;
1037
1038 G4ThreeVector deltaDirection(dirX,dirY,dirZ
1039 deltaDirection.rotateUz(particleDirection);
1040
1041 // create G4DynamicParticle object for delt
1042 G4DynamicParticle* deltaRay = new G4Dynamic
1043 deltaRay->SetKineticEnergy( deltaKineticEne
1044 deltaRay->SetMomentumDirection(deltaDirecti
1045 deltaDirection.y(),
1046 deltaDirection.z());
1047 deltaRay->SetDefinition(G4Electron::Electro
1048
1049 // fill aParticleChange
1050 G4double finalKineticEnergy = kineticEnergy
1051 std::size_t totalNumber = 1;
1052
1053 // Atomic relaxation
1054
1055 // ---- MGP ---- Temporary limitation: curr
1056
1057 std::size_t nSecondaries = 0;
1058 std::vector<G4DynamicParticle*>* secondaryV
1059
1060 if (definition == G4Proton::ProtonDefinitio
1061 {
1062 const G4Material* material = couple->Ge
1063
1064 // Lazy initialization of pixeCrossSect
1065 if (pixeCrossSectionHandler == 0)
1066 {
1067 // Instantiate pixeCrossSectionHandler wi
1068 // Ownership of interpolation is transfer
1069 G4IInterpolator* interpolation = new G4Lo
1070 pixeCrossSectionHandler =
1071 new G4PixeCrossSectionHandler(interpola
1072 G4String fileName("proton");
1073 pixeCrossSectionHandler->LoadShellData(fi
1074 // pixeCrossSectionHandler->PrintData(
1075 }
1076
1077 // Select an atom in the current materi
1078 G4int Z = pixeCrossSectionHandler->Sele
1079 // std::cout << "G4hImpactIonisati
1080
1081 // G4double microscopicCross = Mic
1082 // kineticEnergy
1083 // Z, deltaCut);
1084 // G4double crossFromShells = pixeCrossS
1085
1086 //std::cout << "G4hImpactIonisation: Z=
1087 // << Z
1088 // << ", energy = "
1089 // << kineticEnergy/MeV
1090 // <<" MeV, microscopic = "
1091 // << microscopicCross/barn
1092 // << " barn, from shells = "
1093 // << crossFromShells/barn
1094 // << " barn"
1095 // << std::endl;
1096
1097 // Select a shell in the target atom ba
1098 G4int shellIndex = pixeCrossSectionHand
1099
1100 G4AtomicTransitionManager* transitionMa
1101 const G4AtomicShell* atomicShell = tran
1102 G4double bindingEnergy = atomicShell->B
1103
1104 // if (verboseLevel > 1)
1105 // {
1106 // G4cout << "G4hImpactIonisation::Po
1107 // << Z
1108 // << ", shell = "
1109 // << shellIndex
1110 // << ", bindingE (keV) = "
1111 // << bindingEnergy/keV
1112 // << G4endl;
1113 // }
1114
1115 // Generate PIXE if binding energy larg
1116
1117 G4ParticleDefinition* type = 0;
1118
1119 if (finalKineticEnergy >= bindingEnergy
1120 // && (bindingEnergy >= minGammaEnergy |
1121 {
1122 // Vacancy in subshell shellIndex; shellI
1123 G4int shellId = atomicShell->ShellId();
1124 // Atomic relaxation: generate secondarie
1125 secondaryVector = atomicDeexcitation.Gene
1126
1127 // ---- Debug ----
1128 //std::cout << "ShellId = "
1129 // <<shellId << " ---- Atomic relaxati
1130 // << secondaryVector->size()
1131 // << std::endl;
1132
1133 // ---- End debug ---
1134
1135 if (secondaryVector != 0)
1136 {
1137 nSecondaries = secondaryVector->size(
1138 for (std::size_t i = 0; i<nSecondarie
1139 {
1140 G4DynamicParticle* aSecondary = (*secon
1141 if (aSecondary)
1142 {
1143 G4double e = aSecondary->GetKinetic
1144 type = aSecondary->GetDefinition();
1145
1146 // ---- Debug ----
1147 //if (type == G4Gamma::GammaDefinit
1148 // {
1149 // std::cout << "Z = " << Z
1150 // << ", shell: " << shellId
1151 // << ", PIXE photon energy
1152 // << std::endl;
1153 // }
1154 // ---- End debug ---
1155
1156 if (e < finalKineticEnergy &&
1157 ((type == G4Gamma::Gamma() && e > min
1158 (type == G4Electron::Electron() && e
1159 {
1160 // Subtract the energy of the emitted
1161 finalKineticEnergy -= e;
1162 totalNumber++;
1163 // ---- Debug ----
1164 //if (type == G4Gamma::GammaDefinitio
1165 // {
1166 // std::cout << "Z = " << Z
1167 // << ", shell: " << shellId
1168 // << ", PIXE photon energy (k
1169 // << std::endl;
1170 // }
1171 // ---- End debug ---
1172 }
1173 else
1174 {
1175 // The atomic relaxation product has
1176 // ---- Debug ----
1177 // if (type == G4Gamma::GammaDefiniti
1178 // {
1179 // std::cout << "Z = " << Z
1180 //
1181 // << ", PIXE photon energy =
1182 // << " keV below threshold
1183 // << std::endl;
1184 // }
1185 // ---- End debug ---
1186
1187 delete aSecondary;
1188 (*secondaryVector)[i] = 0;
1189 }
1190 }
1191 }
1192 }
1193 }
1194 }
1195
1196
1197 // Save delta-electrons
1198
1199 G4double eDeposit = 0.;
1200
1201 if (finalKineticEnergy > MinKineticEnergy)
1202 {
1203 G4double finalPx = totalMomentum*partic
1204 G4double finalPy = totalMomentum*partic
1205 G4double finalPz = totalMomentum*partic
1206 G4double finalMomentum = std::sqrt(fina
1207 finalPx /= finalMomentum;
1208 finalPy /= finalMomentum;
1209 finalPz /= finalMomentum;
1210
1211 aParticleChange.ProposeMomentumDirectio
1212 }
1213 else
1214 {
1215 eDeposit = finalKineticEnergy;
1216 finalKineticEnergy = 0.;
1217 aParticleChange.ProposeMomentumDirectio
1218 particleDirection.y(),
1219 particleDirection.z());
1220 if(!aParticle->GetDefinition()->GetProc
1221 GetAtRestProcessVector()->size())
1222 aParticleChange.ProposeTrackStatus(fS
1223 else
1224 aParticleChange.ProposeTrackStatus(fS
1225 }
1226
1227 aParticleChange.ProposeEnergy(finalKineticE
1228 aParticleChange.ProposeLocalEnergyDeposit (
1229 aParticleChange.SetNumberOfSecondaries((G4i
1230 aParticleChange.AddSecondary(deltaRay);
1231
1232 // ---- Debug ----
1233 // std::cout << "RDHadronIonisation - fina
1234 // << finalKineticEnergy/MeV
1235 // << ", delta KineticEnergy (keV) = "
1236 // << deltaKineticEnergy/keV
1237 // << ", energy deposit (MeV) = "
1238 // << eDeposit/MeV
1239 // << std::endl;
1240 // ---- End debug ---
1241
1242 // Save Fluorescence and Auger
1243
1244 if (secondaryVector != 0)
1245 {
1246 for (std::size_t l = 0; l < nSecondarie
1247 {
1248 G4DynamicParticle* secondary = (*secondar
1249 if (secondary) aParticleChange.AddSeconda
1250
1251 // ---- Debug ----
1252 //if (secondary != 0)
1253 // {
1254 // if (secondary->GetDefinition() == G4
1255 // {
1256 // G4double eX = secondary->GetKinetic
1257 // std::cout << " PIXE photon of energ
1258 // << " keV added to ParticleChang
1259 // << std::endl;
1260 // }
1261 //}
1262 // ---- End debug ---
1263
1264 }
1265 delete secondaryVector;
1266 }
1267
1268 return G4VContinuousDiscreteProcess::PostSt
1269 }
1270
1271 // ------------------------------------------
1272
1273 G4double G4hImpactIonisation::ComputeDEDX(con
1274 const G4MaterialCutsCouple* coupl
1275 G4double kineticEnergy)
1276 {
1277 const G4Material* material = couple->GetMat
1278 G4Proton* proton = G4Proton::Proton();
1279 G4AntiProton* antiproton = G4AntiProton::An
1280 G4double dedx = 0.;
1281
1282 G4double tScaled = kineticEnergy * proton_m
1283 charge = aParticle->GetPDGCharge() ;
1284
1285 if( charge > 0.)
1286 {
1287 if (tScaled > protonHighEnergy)
1288 {
1289 dedx = G4EnergyLossTables::GetDEDX(proton
1290 }
1291 else
1292 {
1293 dedx = ProtonParametrisedDEDX(couple,tSca
1294 }
1295 }
1296 else
1297 {
1298 if (tScaled > antiprotonHighEnergy)
1299 {
1300 dedx = G4EnergyLossTables::GetDEDX(antipr
1301 }
1302 else
1303 {
1304 dedx = AntiProtonParametrisedDEDX(couple,
1305 }
1306 }
1307 dedx *= theIonEffChargeModel->TheValue(aPar
1308
1309 return dedx ;
1310 }
1311
1312
1313 G4double G4hImpactIonisation::BarkasTerm(cons
1314 G4double kineticEnergy) const
1315 //Function to compute the Barkas term for pro
1316 //
1317 //Ref. Z_1^3 effect in the stopping power of
1318 // J.C Ashley and R.H.Ritchie
1319 // Physical review B Vol.5 No.7 1 April 1
1320 //
1321 {
1322 static G4ThreadLocal G4double FTable[47][2]
1323 { 0.02, 21.5},
1324 { 0.03, 20.0},
1325 { 0.04, 18.0},
1326 { 0.05, 15.6},
1327 { 0.06, 15.0},
1328 { 0.07, 14.0},
1329 { 0.08, 13.5},
1330 { 0.09, 13.},
1331 { 0.1, 12.2},
1332 { 0.2, 9.25},
1333 { 0.3, 7.0},
1334 { 0.4, 6.0},
1335 { 0.5, 4.5},
1336 { 0.6, 3.5},
1337 { 0.7, 3.0},
1338 { 0.8, 2.5},
1339 { 0.9, 2.0},
1340 { 1.0, 1.7},
1341 { 1.2, 1.2},
1342 { 1.3, 1.0},
1343 { 1.4, 0.86},
1344 { 1.5, 0.7},
1345 { 1.6, 0.61},
1346 { 1.7, 0.52},
1347 { 1.8, 0.5},
1348 { 1.9, 0.43},
1349 { 2.0, 0.42},
1350 { 2.1, 0.3},
1351 { 2.4, 0.2},
1352 { 3.0, 0.13},
1353 { 3.08, 0.1},
1354 { 3.1, 0.09},
1355 { 3.3, 0.08},
1356 { 3.5, 0.07},
1357 { 3.8, 0.06},
1358 { 4.0, 0.051},
1359 { 4.1, 0.04},
1360 { 4.8, 0.03},
1361 { 5.0, 0.024},
1362 { 5.1, 0.02},
1363 { 6.0, 0.013},
1364 { 6.5, 0.01},
1365 { 7.0, 0.009},
1366 { 7.1, 0.008},
1367 { 8.0, 0.006},
1368 { 9.0, 0.0032},
1369 { 10.0, 0.0025} };
1370
1371 // Information on particle and material
1372 G4double kinE = kineticEnergy ;
1373 if(0.5*MeV > kinE) kinE = 0.5*MeV ;
1374 G4double gamma = 1.0 + kinE / proton_mass_c
1375 G4double beta2 = 1.0 - 1.0/(gamma*gamma) ;
1376 if(0.0 >= beta2) return 0.0;
1377
1378 G4double BTerm = 0.0;
1379 //G4double AMaterial = 0.0;
1380 G4double ZMaterial = 0.0;
1381 const G4ElementVector* theElementVector = m
1382 G4int numberOfElements = (G4int)material->G
1383
1384 for (G4int i = 0; i<numberOfElements; i++)
1385
1386 //AMaterial = (*theElementVector)[i]->Get
1387 ZMaterial = (*theElementVector)[i]->GetZ(
1388
1389 G4double X = 137.0 * 137.0 * beta2 / ZMat
1390
1391 // Variables to compute L_1
1392 G4double Eta0Chi = 0.8;
1393 G4double EtaChi = Eta0Chi * ( 1.0 + 6.02*
1394 G4double W = ( EtaChi * std::pow( ZMateri
1395 G4double FunctionOfW = FTable[46][1]*FTab
1396
1397 for(G4int j=0; j<47; j++) {
1398
1399 if( W < FTable[j][0] ) {
1400
1401 if(0 == j) {
1402 FunctionOfW = FTable[0][1] ;
1403
1404 } else {
1405 FunctionOfW = (FTable[j][1] - FTabl
1406 / (FTable[j][0] - FTable[j-1][0])
1407 + FTable[j-1][1] ;
1408 }
1409
1410 break;
1411 }
1412
1413 }
1414
1415 BTerm += FunctionOfW /( std::sqrt(ZMateri
1416 }
1417
1418 BTerm *= twopi_mc2_rcl2 * (material->GetEle
1419
1420 return BTerm;
1421 }
1422
1423
1424
1425 G4double G4hImpactIonisation::BlochTerm(const
1426 G4double kineticEnergy,
1427 G4double cSquare) const
1428 //Function to compute the Bloch term for prot
1429 //
1430 //Ref. Z_1^3 effect in the stopping power of
1431 // J.C Ashley and R.H.Ritchie
1432 // Physical review B Vol.5 No.7 1 April 1
1433 //
1434 {
1435 G4double eLoss = 0.0 ;
1436 G4double gamma = 1.0 + kineticEnergy / prot
1437 G4double beta2 = 1.0 - 1.0/(gamma*gamma) ;
1438 G4double y = cSquare / (137.0*137.0*beta2)
1439
1440 if(y < 0.05) {
1441 eLoss = 1.202 ;
1442
1443 } else {
1444 eLoss = 1.0 / (1.0 + y) ;
1445 G4double de = eLoss ;
1446
1447 for(G4int i=2; de>eLoss*0.01; i++) {
1448 de = 1.0/( i * (i*i + y)) ;
1449 eLoss += de ;
1450 }
1451 }
1452 eLoss *= -1.0 * y * cSquare * twopi_mc2_rcl
1453 (material->GetElectronDensity()) / beta2
1454
1455 return eLoss;
1456 }
1457
1458
1459
1460 G4double G4hImpactIonisation::ElectronicLossF
1461 const G4DynamicParticle* partic
1462 const G4MaterialCutsCouple* cou
1463 G4double meanLoss,
1464 G4double step) const
1465 // calculate actual loss from the mean loss
1466 // The model used to get the fluctuation is
1467 // as in Glandz in Geant3.
1468 {
1469 // data members to speed up the fluctuation
1470 // G4int imat ;
1471 // G4double f1Fluct,f2Fluct,e1Fluct,e2Fluc
1472 // G4double e1LogFluct,e2LogFluct,ipotLogF
1473
1474 static const G4double minLoss = 1.*eV ;
1475 static const G4double kappa = 10. ;
1476 static const G4double theBohrBeta2 = 50.0 *
1477
1478 const G4Material* material = couple->GetMat
1479 G4int imaterial = couple->GetIndex() ;
1480 G4double ipotFluct = material->GetIonisat
1481 G4double electronDensity = material->GetEle
1482 G4double zeff = electronDensity/(material->
1483
1484 // get particle data
1485 G4double tkin = particle->GetKineticEnerg
1486 G4double particleMass = particle->GetMass()
1487 G4double deltaCutInKineticEnergyNow = cutFo
1488
1489 // shortcut for very very small loss
1490 if(meanLoss < minLoss) return meanLoss ;
1491
1492 // Validity range for delta electron cross
1493 G4double threshold = std::max(deltaCutInKin
1494 G4double loss, siga;
1495
1496 G4double rmass = electron_mass_c2/particleM
1497 G4double tau = tkin/particleMass;
1498 G4double tau1 = tau+1.0;
1499 G4double tau2 = tau*(tau+2.);
1500 G4double tMax = 2.*electron_mass_c2*tau2
1501
1502
1503 if(tMax > threshold) tMax = threshold;
1504 G4double beta2 = tau2/(tau1*tau1);
1505
1506 // Gaussian fluctuation
1507 if(meanLoss > kappa*tMax || tMax < kappa*ip
1508 {
1509 siga = tMax * (1.0-0.5*beta2) * step *
1510 * electronDensity / beta2 ;
1511
1512 // High velocity or negatively charged
1513 if( beta2 > 3.0*theBohrBeta2*zeff || ch
1514 siga = std::sqrt( siga * chargeSquare ) ;
1515
1516 // Low velocity - additional ion charge flu
1517 // Q.Yang et al., NIM B61(1991)149-155.
1518 } else {
1519 G4double chu = theIonChuFluctuationModel->T
1520 G4double yang = theIonYangFluctuationModel-
1521 siga = std::sqrt( siga * (chargeSquare * ch
1522 }
1523
1524 do {
1525 loss = G4RandGauss::shoot(meanLoss,si
1526 } while (loss < 0. || loss > 2.0*meanLo
1527 return loss;
1528 }
1529
1530 // Non Gaussian fluctuation
1531 static const G4double probLim = 0.01 ;
1532 static const G4double sumaLim = -std::log(p
1533 static const G4double alim = 10.;
1534
1535 G4double suma,w1,w2,C,e0,lossc,w;
1536 G4double a1,a2,a3;
1537 G4int p1,p2,p3;
1538 G4int nb;
1539 G4double corrfac, na,alfa,rfac,namean,sa,al
1540 G4double dp3;
1541
1542 G4double f1Fluct = material->GetIonisat
1543 G4double f2Fluct = material->GetIonisat
1544 G4double e1Fluct = material->GetIonisat
1545 G4double e2Fluct = material->GetIonisat
1546 G4double e1LogFluct = material->GetIonisat
1547 G4double e2LogFluct = material->GetIonisat
1548 G4double rateFluct = material->GetIonisat
1549 G4double ipotLogFluct= material->GetIonisat
1550
1551 w1 = tMax/ipotFluct;
1552 w2 = std::log(2.*electron_mass_c2*tau2);
1553
1554 C = meanLoss*(1.-rateFluct)/(w2-ipotLogFluc
1555
1556 a1 = C*f1Fluct*(w2-e1LogFluct-beta2)/e1Fluc
1557 a2 = C*f2Fluct*(w2-e2LogFluct-beta2)/e2Fluc
1558 a3 = rateFluct*meanLoss*(tMax-ipotFluct)/(i
1559 if(a1 < 0.0) a1 = 0.0;
1560 if(a2 < 0.0) a2 = 0.0;
1561 if(a3 < 0.0) a3 = 0.0;
1562
1563 suma = a1+a2+a3;
1564
1565 loss = 0.;
1566
1567
1568 if(suma < sumaLim) // very smal
1569 {
1570 e0 = material->GetIonisation()->GetEner
1571
1572 if(tMax == ipotFluct)
1573 {
1574 a3 = meanLoss/e0;
1575
1576 if(a3>alim)
1577 {
1578 siga=std::sqrt(a3) ;
1579 p3 = std::max(0,G4int(G4RandGauss::sh
1580 }
1581 else
1582 p3 = (G4int)G4Poisson(a3);
1583
1584 loss = p3*e0 ;
1585
1586 if(p3 > 0)
1587 loss += (1.-2.*G4UniformRand())*e0 ;
1588
1589 }
1590 else
1591 {
1592 tMax = tMax-ipotFluct+e0 ;
1593 a3 = meanLoss*(tMax-e0)/(tMax*e0*std::log
1594
1595 if(a3>alim)
1596 {
1597 siga=std::sqrt(a3) ;
1598 p3 = std::max(0,int(G4RandGauss::shoo
1599 }
1600 else
1601 p3 = (G4int)G4Poisson(a3);
1602
1603 if(p3 > 0)
1604 {
1605 w = (tMax-e0)/tMax ;
1606 if(p3 > nmaxCont2)
1607 {
1608 dp3 = G4float(p3) ;
1609 corrfac = dp3/G4float(nmaxCont2) ;
1610 p3 = G4int(nmaxCont2) ;
1611 }
1612 else
1613 corrfac = 1. ;
1614
1615 for(G4int i=0; i<p3; i++) loss += 1./
1616 loss *= e0*corrfac ;
1617 }
1618 }
1619 }
1620
1621 else // not so
1622 {
1623 // excitation type 1
1624 if(a1>alim)
1625 {
1626 siga=std::sqrt(a1) ;
1627 p1 = std::max(0,G4int(G4RandGauss::shoot(
1628 }
1629 else
1630 p1 = (G4int)G4Poisson(a1);
1631
1632 // excitation type 2
1633 if(a2>alim)
1634 {
1635 siga=std::sqrt(a2) ;
1636 p2 = std::max(0,G4int(G4RandGauss::shoot(
1637 }
1638 else
1639 p2 = (G4int)G4Poisson(a2);
1640
1641 loss = p1*e1Fluct+p2*e2Fluct;
1642
1643 // smearing to avoid unphysical peaks
1644 if(p2 > 0)
1645 loss += (1.-2.*G4UniformRand())*e2Flu
1646 else if (loss>0.)
1647 loss += (1.-2.*G4UniformRand())*e1Flu
1648
1649 // ionisation .........................
1650 if(a3 > 0.)
1651 {
1652 if(a3>alim)
1653 {
1654 siga=std::sqrt(a3) ;
1655 p3 = std::max(0,G4int(G4RandGauss::sh
1656 }
1657 else
1658 p3 = (G4int)G4Poisson(a3);
1659
1660 lossc = 0.;
1661 if(p3 > 0)
1662 {
1663 na = 0.;
1664 alfa = 1.;
1665 if (p3 > nmaxCont2)
1666 {
1667 dp3 = G4float(p3);
1668 rfac = dp3/(G4float(nmaxCont2)+dp
1669 namean = G4float(p3)*rfac;
1670 sa = G4float(nmaxCont1)*rfac;
1671 na = G4RandGauss::shoot(namean,
1672 if (na > 0.)
1673 {
1674 alfa = w1*G4float(nmaxCont2+p3)/
1675 (w1*G4float(nmaxCont2)+G4float(p3));
1676 alfa1 = alfa*std::log(alfa)/(alfa-
1677 ea = na*ipotFluct*alfa1;
1678 sea = ipotFluct*std::sqrt(na*(al
1679 lossc += G4RandGauss::shoot(ea,sea)
1680 }
1681 }
1682
1683 nb = G4int(G4float(p3)-na);
1684 if (nb > 0)
1685 {
1686 w2 = alfa*ipotFluct;
1687 w = (tMax-w2)/tMax;
1688 for (G4int k=0; k<nb; k++) lossc += w2/
1689 }
1690 }
1691 loss += lossc;
1692 }
1693 }
1694
1695 return loss ;
1696 }
1697
1698
1699
1700 void G4hImpactIonisation::SetCutForSecondaryP
1701 {
1702 minGammaEnergy = cut;
1703 }
1704
1705
1706
1707 void G4hImpactIonisation::SetCutForAugerElect
1708 {
1709 minElectronEnergy = cut;
1710 }
1711
1712
1713
1714 void G4hImpactIonisation::ActivateAugerElectr
1715 {
1716 atomicDeexcitation.ActivateAugerElectronPro
1717 }
1718
1719
1720
1721 void G4hImpactIonisation::PrintInfoDefinition
1722 {
1723 G4String comments = " Knock-on electron cr
1724 comments += "\n Good description abo
1725 comments += " Delta ray energy sampl
1726
1727 G4cout << G4endl << GetProcessName() << ":
1728 << "\n PhysicsTables from " <
1729 << " to " << HighestKineticEnergy /
1730 << " in " << TotBin << " bins."
1731 << "\n Electronic stopping power mo
1732 << protonTable
1733 << "\n from " << protonLowEne
1734 << " to " << protonHighEnergy / MeV
1735 G4cout << "\n Parametrisation model
1736 << antiprotonTable
1737 << "\n from " << antiprotonLo
1738 << " to " << antiprotonHighEnergy /
1739 if(theBarkas){
1740 G4cout << " Parametrization of the
1741 << G4endl ;
1742 }
1743 if(nStopping) {
1744 G4cout << " Nuclear stopping power
1745 << G4endl ;
1746 }
1747
1748 G4bool printHead = true;
1749
1750 const G4ProductionCutsTable* theCoupleTable
1751 G4ProductionCutsTable::GetProductionCutsT
1752 G4int numOfCouples = (G4int)theCoupleTable-
1753
1754 // loop for materials
1755
1756 for (G4int j=0 ; j < numOfCouples; ++j) {
1757
1758 const G4MaterialCutsCouple* couple = theC
1759 const G4Material* material= couple->GetMa
1760 G4double deltaCutNow = cutForDelta[(coupl
1761 G4double excitationEnergy = material->Get
1762
1763 if(excitationEnergy > deltaCutNow) {
1764 if(printHead) {
1765 printHead = false ;
1766
1767 G4cout << " material
1768 G4cout << G4endl;
1769 }
1770
1771 G4cout << std::setw(20) << material->Ge
1772 << std::setw(15) << excitationEnergy/k
1773 }
1774 }
1775 }
1776