| Geant4 Cross Reference |
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"
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 G4String& processName)
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 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 // Calculated on the fly, but should be initialized to sensible values
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 theNuclearStoppingModel;
151 if (theIonEffChargeModel) delete theIonEffChargeModel;
152 if (theIonChuFluctuationModel) delete theIonChuFluctuationModel;
153 if (theIonYangFluctuationModel) delete theIonYangFluctuationModel;
154
155 delete pixeCrossSectionHandler;
156
157 // ---- MGP ---- The following is to be checked
158 // if (shellVacancy) delete shellVacancy;
159
160 cutForDelta.clear();
161 }
162
163 // --------------------------------------------------------------------
164 void G4hImpactIonisation::SetElectronicStoppingPowerModel(const G4ParticleDefinition* particle,
165 const G4String& dedxTable)
166 // This method defines the ionisation parametrisation method via its name
167 {
168 if (particle->GetPDGCharge() > 0 )
169 {
170 // Positive charge
171 SetProtonElectronicStoppingPowerModel(dedxTable) ;
172 }
173 else
174 {
175 // Antiprotons
176 SetAntiProtonElectronicStoppingPowerModel(dedxTable) ;
177 }
178 }
179
180
181 // --------------------------------------------------------------------
182 void G4hImpactIonisation::InitializeParametrisation()
183
184 {
185 // Define models for parametrisation of electronic energy losses
186 betheBlochModel = new G4hBetheBlochModel("Bethe-Bloch") ;
187 protonModel = new G4hParametrisedLossModel(protonTable) ;
188 protonHighEnergy = std::min(protonHighEnergy,protonModel->HighEnergyLimit(0, 0));
189 antiprotonModel = new G4QAOLowEnergyLoss(antiprotonTable) ;
190 theNuclearStoppingModel = new G4hNuclearStoppingModel(theNuclearTable) ;
191 theIonEffChargeModel = new G4hIonEffChargeSquare("Ziegler1988") ;
192 theIonChuFluctuationModel = new G4IonChuFluctuationModel("Chu") ;
193 theIonYangFluctuationModel = new G4IonYangFluctuationModel("Yang") ;
194 }
195
196
197 // --------------------------------------------------------------------
198 void G4hImpactIonisation::BuildPhysicsTable(const G4ParticleDefinition& particleDef)
199
200 // just call BuildLossTable+BuildLambdaTable
201 {
202
203 // Verbose print-out
204 if(verboseLevel > 0)
205 {
206 G4cout << "G4hImpactIonisation::BuildPhysicsTable for "
207 << particleDef.GetParticleName()
208 << " mass(MeV)= " << particleDef.GetPDGMass()/MeV
209 << " charge= " << particleDef.GetPDGCharge()/eplus
210 << " type= " << particleDef.GetParticleType()
211 << G4endl;
212
213 if(verboseLevel > 1)
214 {
215 G4ProcessVector* pv = particleDef.GetProcessManager()->GetProcessList();
216
217 G4cout << " 0: " << (*pv)[0]->GetProcessName() << " " << (*pv)[0]
218 << " 1: " << (*pv)[1]->GetProcessName() << " " << (*pv)[1]
219 // << " 2: " << (*pv)[2]->GetProcessName() << " " << (*pv)[2]
220 << G4endl;
221 G4cout << "ionModel= " << theIonEffChargeModel
222 << " MFPtable= " << theMeanFreePathTable
223 << " iniMass= " << initialMass
224 << G4endl;
225 }
226 }
227 // End of verbose print-out
228
229 if (particleDef.GetParticleType() == "nucleus" &&
230 particleDef.GetParticleName() != "GenericIon" &&
231 particleDef.GetParticleSubType() == "generic")
232 {
233
234 G4EnergyLossTables::Register(&particleDef,
235 theDEDXpTable,
236 theRangepTable,
237 theInverseRangepTable,
238 theLabTimepTable,
239 theProperTimepTable,
240 LowestKineticEnergy, HighestKineticEnergy,
241 proton_mass_c2/particleDef.GetPDGMass(),
242 TotBin);
243
244 return;
245 }
246
247 if( !CutsWhereModified() && theLossTable) return;
248
249 InitializeParametrisation() ;
250 G4Proton* proton = G4Proton::Proton();
251 G4AntiProton* antiproton = G4AntiProton::AntiProton();
252
253 charge = particleDef.GetPDGCharge() / eplus;
254 chargeSquare = charge*charge ;
255
256 const G4ProductionCutsTable* theCoupleTable=
257 G4ProductionCutsTable::GetProductionCutsTable();
258 G4int numOfCouples = (G4int)theCoupleTable->GetTableSize();
259
260 cutForDelta.clear();
261 cutForGamma.clear();
262
263 for (G4int j=0; j<numOfCouples; ++j) {
264
265 // get material parameters needed for the energy loss calculation
266 const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(j);
267 const G4Material* material= couple->GetMaterial();
268
269 // the cut cannot be below lowest limit
270 G4double tCut = (*(theCoupleTable->GetEnergyCutsVector(1)))[j];
271 if(tCut > HighestKineticEnergy) tCut = HighestKineticEnergy;
272
273 G4double excEnergy = material->GetIonisation()->GetMeanExcitationEnergy();
274
275 tCut = std::max(tCut,excEnergy);
276 cutForDelta.push_back(tCut);
277
278 // the cut cannot be below lowest limit
279 tCut = (*(theCoupleTable->GetEnergyCutsVector(0)))[j];
280 if(tCut > HighestKineticEnergy) tCut = HighestKineticEnergy;
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 dimension (see src/G4hImpactLoss.cc for more details)
295 // It happended in the past that caused memory corruption errors. The problem is still pending, even if temporary solved
296 // G4cout << "[NOTE]: __LINE__=" << __LINE__ << ", particleDef=" << particleDef.GetParticleName() << ", proton=" << proton << ", theLossTable=" << theLossTable << ", CounterOfpProcess=" << CounterOfpProcess << G4endl;
297
298 RecorderOfpProcess[CounterOfpProcess] = theLossTable ;
299 CounterOfpProcess++;
300 }
301 } else {
302 {
303 BuildLossTable(*antiproton) ;
304
305 // The following vector has a fixed dimension (see src/G4hImpactLoss.cc for more details)
306 // It happended in the past that caused memory corruption errors. The problem is still pending, even if temporary solved
307 // G4cout << "[NOTE]: __LINE__=" << __LINE__ << ", particleDef=" << particleDef.GetParticleName() << ", antiproton=" << antiproton << ", theLossTable=" << theLossTable << ", CounterOfpbarProcess=" << CounterOfpbarProcess << G4endl;
308
309 RecorderOfpbarProcess[CounterOfpbarProcess] = theLossTable ;
310 CounterOfpbarProcess++;
311 }
312 }
313
314 if(verboseLevel > 0) {
315 G4cout << "G4hImpactIonisation::BuildPhysicsTable: "
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) << G4endl;
326 }
327
328 if(verboseLevel > 0) {
329 G4cout << "G4hImpactIonisation::BuildPhysicsTable: "
330 << "DEDX table will be built "
331 // << theDEDXpTable << " " << theDEDXpbarTable
332 // << " " << theRangepTable << " " << theRangepbarTable
333 << G4endl;
334 }
335
336 BuildDEDXTable(particleDef) ;
337
338 if(verboseLevel > 1) {
339 G4cout << (*theDEDXpTable) << G4endl;
340 }
341
342 if((&particleDef == proton) || (&particleDef == antiproton)) PrintInfoDefinition() ;
343
344 if(verboseLevel > 0) {
345 G4cout << "G4hImpactIonisation::BuildPhysicsTable: end for "
346 << particleDef.GetParticleName() << G4endl;
347 }
348 }
349
350
351
352
353
354 // --------------------------------------------------------------------
355 void G4hImpactIonisation::BuildLossTable(const G4ParticleDefinition& particleDef)
356 {
357 // Initialisation
358 G4double lowEdgeEnergy , ionloss, ionlossBB, paramB ;
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::GetProductionCutsTable();
379 G4int numOfCouples = (G4int)theCoupleTable->GetTableSize();
380
381 if ( theLossTable)
382 {
383 theLossTable->clearAndDestroy();
384 delete theLossTable;
385 }
386
387 theLossTable = new G4PhysicsTable(numOfCouples);
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 G4PhysicsLogVector(LowestKineticEnergy,
394 HighestKineticEnergy,
395 TotBin);
396
397 // get material parameters needed for the energy loss calculation
398 const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(j);
399 const G4Material* material= couple->GetMaterial();
400
401 if ( charge > 0.0 ) {
402 ionloss = ProtonParametrisedDEDX(couple,highEnergy) ;
403 } else {
404 ionloss = AntiProtonParametrisedDEDX(couple,highEnergy) ;
405 }
406
407 ionlossBB = betheBlochModel->TheValue(&particleDef,material,highEnergy) ;
408 ionlossBB -= DeltaRaysEnergy(couple,highEnergy,proton_mass_c2) ;
409
410
411 paramB = ionloss/ionlossBB - 1.0 ;
412
413 // now comes the loop for the kinetic energy values
414 for (G4int i = 0 ; i < TotBin ; i++) {
415 lowEdgeEnergy = aVector->GetLowEdgeEnergy(i) ;
416
417 // low energy part for this material, parametrised energy loss formulae
418 if ( lowEdgeEnergy < highEnergy ) {
419
420 if ( charge > 0.0 ) {
421 ionloss = ProtonParametrisedDEDX(couple,lowEdgeEnergy) ;
422 } else {
423 ionloss = AntiProtonParametrisedDEDX(couple,lowEdgeEnergy) ;
424 }
425
426 } else {
427
428 // high energy part for this material, Bethe-Bloch formula
429 ionloss = betheBlochModel->TheValue(proton,material,
430 lowEdgeEnergy) ;
431
432 ionloss -= DeltaRaysEnergy(couple,lowEdgeEnergy,proton_mass_c2) ;
433
434 ionloss *= (1.0 + paramB*highEnergy/lowEdgeEnergy) ;
435 }
436
437 // now put the loss into the vector
438 if(verboseLevel > 1) {
439 G4cout << "E(MeV)= " << lowEdgeEnergy/MeV
440 << " dE/dx(MeV/mm)= " << ionloss*mm/MeV
441 << " in " << material->GetName() << G4endl;
442 }
443 aVector->PutValue(i,ionloss) ;
444 }
445 // Insert vector for this material into the table
446 theLossTable->insert(aVector) ;
447 }
448 }
449
450
451
452 // --------------------------------------------------------------------
453 void G4hImpactIonisation::BuildLambdaTable(const G4ParticleDefinition& particleDef)
454
455 {
456 // Build mean free path tables for the delta ray production process
457 // tables are built for MATERIALS
458
459 if(verboseLevel > 1) {
460 G4cout << "G4hImpactIonisation::BuildLambdaTable for "
461 << particleDef.GetParticleName() << " is started" << G4endl;
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::GetProductionCutsTable();
472 G4int numOfCouples = (G4int)theCoupleTable->GetTableSize();
473
474
475 if (theMeanFreePathTable) {
476 theMeanFreePathTable->clearAndDestroy();
477 delete theMeanFreePathTable;
478 }
479
480 theMeanFreePathTable = new G4PhysicsTable(numOfCouples);
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 G4PhysicsLogVector(LowestKineticEnergy,
488 HighestKineticEnergy,
489 TotBin);
490
491 // compute the (macroscopic) cross section first
492 const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(j);
493 const G4Material* material= couple->GetMaterial();
494
495 const G4ElementVector* theElementVector = material->GetElementVector() ;
496 const G4double* theAtomicNumDensityVector = material->GetAtomicNumDensityVector();
497 const G4int numberOfElements = (G4int)material->GetNumberOfElements() ;
498
499 // get the electron kinetic energy cut for the actual material,
500 // it will be used in ComputeMicroscopicCrossSection
501 // ( it is the SAME for ALL the ELEMENTS in THIS MATERIAL )
502 // ------------------------------------------------------
503
504 G4double deltaCut = cutForDelta[j];
505
506 for ( G4int i = 0 ; i < TotBin ; i++ ) {
507 lowEdgeEnergy = aVector->GetLowEdgeEnergy(i) ;
508 G4double sigma = 0.0 ;
509 G4int Z;
510
511 for (G4int iel=0; iel<numberOfElements; iel++ )
512 {
513 Z = (G4int) (*theElementVector)[iel]->GetZ();
514 // ---- MGP --- Corrected duplicated cross section calculation here from
515 // G4hLowEnergyIonisation original
516 G4double microCross = MicroscopicCrossSection( particleDef,
517 lowEdgeEnergy,
518 Z,
519 deltaCut ) ;
520 //totalCrossSectionMap [Z] = microCross;
521 sigma += theAtomicNumDensityVector[iel] * microCross ;
522 }
523
524 // mean free path = 1./macroscopic cross section
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::MicroscopicCrossSection(const G4ParticleDefinition& particleDef,
539 G4double kineticEnergy,
540 G4double atomicNumber,
541 G4double deltaCutInEnergy) const
542 {
543 //******************************************************************
544 // cross section formula is OK for spin=0, 1/2, 1 only !
545 // *****************************************************************
546
547 // Calculates the microscopic cross section in GEANT4 internal units
548 // Formula documented in Geant4 Phys. Ref. Manual
549 // ( it is called for elements, AtomicNumber = z )
550
551 G4double totalCrossSection = 0.;
552
553 // Particle mass and energy
554 G4double particleMass = initialMass;
555 G4double energy = kineticEnergy + particleMass;
556
557 // Some kinematics
558 G4double gamma = energy / particleMass;
559 G4double beta2 = 1. - 1. / (gamma * gamma);
560 G4double var = electron_mass_c2 / particleMass;
561 G4double tMax = 2. * electron_mass_c2 * (gamma*gamma - 1.) / (1. + 2.* gamma*var + var*var);
562
563 // Calculate the total cross section
564
565 if ( tMax > deltaCutInEnergy )
566 {
567 var = deltaCutInEnergy / tMax;
568 totalCrossSection = (1. - var * (1. - beta2 * std::log(var))) / deltaCutInEnergy ;
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 - deltaCutInEnergy) / (energy*energy);
576 }
577 // +term for spin=1 particle
578 else if (spin > 0.9 )
579 {
580 totalCrossSection += -std::log(var) /
581 (3. * deltaCutInEnergy) + (tMax - deltaCutInEnergy) * ( (5. + 1. /var)*0.25 / (energy*energy) -
582 beta2 / (tMax * deltaCutInEnergy) ) / 3. ;
583 }
584 totalCrossSection *= twopi_mc2_rcl2 * atomicNumber / beta2 ;
585 }
586
587 //std::cout << "Microscopic = " << totalCrossSection/barn
588 // << ", e = " << kineticEnergy/MeV <<std:: endl;
589
590 return totalCrossSection ;
591 }
592
593
594
595 // --------------------------------------------------------------------
596 G4double G4hImpactIonisation::GetMeanFreePath(const G4Track& track,
597 G4double, // previousStepSize
598 enum G4ForceCondition* condition)
599 {
600 const G4DynamicParticle* dynamicParticle = track.GetDynamicParticle();
601 const G4MaterialCutsCouple* couple = track.GetMaterialCutsCouple();
602 const G4Material* material = couple->GetMaterial();
603
604 G4double meanFreePath = DBL_MAX;
605 // ---- MGP ---- What is the meaning of the local variable isOutOfRange?
606 G4bool isOutRange = false;
607
608 *condition = NotForced ;
609
610 G4double kineticEnergy = (dynamicParticle->GetKineticEnergy())*initialMass/(dynamicParticle->GetMass());
611 charge = dynamicParticle->GetCharge()/eplus;
612 chargeSquare = theIonEffChargeModel->TheValue(dynamicParticle, material);
613
614 if (kineticEnergy < LowestKineticEnergy)
615 {
616 meanFreePath = DBL_MAX;
617 }
618 else
619 {
620 if (kineticEnergy > HighestKineticEnergy) kineticEnergy = HighestKineticEnergy;
621 meanFreePath = (((*theMeanFreePathTable)(couple->GetIndex()))->
622 GetValue(kineticEnergy,isOutRange))/chargeSquare;
623 }
624
625 return meanFreePath ;
626 }
627
628
629 // --------------------------------------------------------------------
630 G4double G4hImpactIonisation::GetConstraints(const G4DynamicParticle* particle,
631 const G4MaterialCutsCouple* couple)
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->GetMaterial();
638 G4Proton* proton = G4Proton::Proton();
639 G4AntiProton* antiproton = G4AntiProton::AntiProton();
640
641 G4double stepLimit = 0.;
642 G4double dx, highEnergy;
643
644 G4double massRatio = proton_mass_c2/(particle->GetMass()) ;
645 G4double kineticEnergy = particle->GetKineticEnergy() ;
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(proton, tScaled, couple);
656 dx = G4EnergyLossTables::GetRange(proton, highEnergy, couple);
657 fdEdx = G4EnergyLossTables::GetDEDX(proton, tScaled, couple)
658 * chargeSquare ;
659
660 // Correction for positive ions
661 if (theBarkas && tScaled > highEnergy)
662 {
663 fBarkas = BarkasTerm(material,tScaled)*std::sqrt(chargeSquare)*chargeSquare
664 + BlochTerm(material,tScaled,chargeSquare);
665 }
666 // Antiprotons and negative hadrons
667 }
668 else
669 {
670 highEnergy = antiprotonHighEnergy ;
671 fRangeNow = G4EnergyLossTables::GetRange(antiproton, tScaled, couple);
672 dx = G4EnergyLossTables::GetRange(antiproton, highEnergy, couple);
673 fdEdx = G4EnergyLossTables::GetDEDX(antiproton, tScaled, couple) * chargeSquare ;
674
675 if (theBarkas && tScaled > highEnergy)
676 {
677 fBarkas = -BarkasTerm(material,tScaled)*std::sqrt(chargeSquare)*chargeSquare
678 + BlochTerm(material,tScaled,chargeSquare);
679 }
680 }
681 /*
682 const G4Material* mat = couple->GetMaterial();
683 G4double fac = gram/(MeV*cm2*mat->GetDensity());
684 G4cout << particle->GetDefinition()->GetParticleName()
685 << " in " << mat->GetName()
686 << " E(MeV)= " << kineticEnergy/MeV
687 << " dedx(MeV*cm^2/g)= " << fdEdx*fac
688 << " barcas(MeV*cm^2/gram)= " << fBarkas*fac
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->GetProductionCuts()
698 ->GetProductionCut(idxG4ElectronCut));
699
700 if (fRangeNow > r)
701 {
702 stepLimit = dRoverRange*fRangeNow + r*(1.0 - dRoverRange)*(2.0 - r/fRangeNow);
703 if (stepLimit > fRangeNow) stepLimit = fRangeNow;
704 }
705 // compute the (random) Step limit in standard energy range
706 if(tScaled > highEnergy )
707 {
708 // add Barkas correction directly to dedx
709 fdEdx += fBarkas;
710
711 if(stepLimit > fRangeNow - dx*0.9) stepLimit = fRangeNow - dx*0.9 ;
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::AlongStepDoIt(const G4Track& track,
726 const G4Step& step)
727 {
728 // compute the energy loss after a step
729 G4Proton* proton = G4Proton::Proton();
730 G4AntiProton* antiproton = G4AntiProton::AntiProton();
731 G4double finalT = 0.;
732
733 aParticleChange.Initialize(track) ;
734
735 const G4MaterialCutsCouple* couple = track.GetMaterialCutsCouple();
736 const G4Material* material = couple->GetMaterial();
737
738 // get the actual (true) Step length from step
739 const G4double stepLength = step.GetStepLength() ;
740
741 const G4DynamicParticle* particle = track.GetDynamicParticle() ;
742
743 G4double kineticEnergy = particle->GetKineticEnergy() ;
744 G4double massRatio = proton_mass_c2/(particle->GetMass()) ;
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 energy range
755 }
756
757 else if( kineticEnergy > HighestKineticEnergy)
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 * chargeSquare ;
774 G4double sScaled = stepLength * massRatio * chargeSquare ;
775
776 if(charge > 0.0)
777 {
778 eLoss = G4EnergyLossTables::GetPreciseEnergyFromRange(proton,rScaled, couple) -
779 G4EnergyLossTables::GetPreciseEnergyFromRange(proton,rScaled-sScaled,couple) ;
780
781 }
782 else
783 {
784 // Antiproton
785 eLoss = G4EnergyLossTables::GetPreciseEnergyFromRange(antiproton,rScaled,couple) -
786 G4EnergyLossTables::GetPreciseEnergyFromRange(antiproton,rScaled-sScaled,couple) ;
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 < protonHighEnergy)
800 {
801 nLoss = (theNuclearStoppingModel->TheValue(particle, material)) * stepLength;
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 > MinKineticEnergy)
810 {
811
812 // now the electron loss with fluctuation
813 eLoss = ElectronicLossFluctuation(particle, couple, eLoss, stepLength) ;
814 if (eLoss < 0.0) eLoss = 0.0;
815 finalT = kineticEnergy - eLoss - nLoss;
816 }
817
818 // stop particle if the kinetic energy <= MinKineticEnergy
819 if (finalT*massRatio <= MinKineticEnergy )
820 {
821
822 finalT = 0.0;
823 if (!particle->GetDefinition()->GetProcessManager()->GetAtRestProcessVector()->size())
824 aParticleChange.ProposeTrackStatus(fStopAndKill);
825 else
826 aParticleChange.ProposeTrackStatus(fStopButAlive);
827 }
828
829 aParticleChange.ProposeEnergy( finalT );
830 eLoss = kineticEnergy-finalT;
831
832 aParticleChange.ProposeLocalEnergyDeposit(eLoss);
833 return &aParticleChange ;
834 }
835
836
837
838 // --------------------------------------------------------------------
839 G4double G4hImpactIonisation::ProtonParametrisedDEDX(const G4MaterialCutsCouple* couple,
840 G4double kineticEnergy) const
841 {
842 const G4Material* material = couple->GetMaterial();
843 G4Proton* proton = G4Proton::Proton();
844 G4double eLoss = 0.;
845
846 // Free Electron Gas Model
847 if(kineticEnergy < protonLowEnergy) {
848 eLoss = (protonModel->TheValue(proton, material, protonLowEnergy))
849 * std::sqrt(kineticEnergy/protonLowEnergy) ;
850
851 // Parametrisation
852 } else {
853 eLoss = protonModel->TheValue(proton, material, kineticEnergy) ;
854 }
855
856 // Delta rays energy
857 eLoss -= DeltaRaysEnergy(couple,kineticEnergy,proton_mass_c2) ;
858
859 if(verboseLevel > 2) {
860 G4cout << "p E(MeV)= " << kineticEnergy/MeV
861 << " dE/dx(MeV/mm)= " << eLoss*mm/MeV
862 << " for " << material->GetName()
863 << " model: " << protonModel << G4endl;
864 }
865
866 if(eLoss < 0.0) eLoss = 0.0 ;
867
868 return eLoss ;
869 }
870
871
872
873 // --------------------------------------------------------------------
874 G4double G4hImpactIonisation::AntiProtonParametrisedDEDX(const G4MaterialCutsCouple* couple,
875 G4double kineticEnergy) const
876 {
877 const G4Material* material = couple->GetMaterial();
878 G4AntiProton* antiproton = G4AntiProton::AntiProton();
879 G4double eLoss = 0.0 ;
880
881 // Antiproton model is used
882 if(antiprotonModel->IsInCharge(antiproton,material)) {
883 if(kineticEnergy < antiprotonLowEnergy) {
884 eLoss = antiprotonModel->TheValue(antiproton,material,antiprotonLowEnergy)
885 * std::sqrt(kineticEnergy/antiprotonLowEnergy) ;
886
887 // Parametrisation
888 } else {
889 eLoss = antiprotonModel->TheValue(antiproton,material,
890 kineticEnergy);
891 }
892
893 // The proton model is used + Barkas correction
894 } else {
895 if(kineticEnergy < protonLowEnergy) {
896 eLoss = protonModel->TheValue(G4Proton::Proton(),material,protonLowEnergy)
897 * std::sqrt(kineticEnergy/protonLowEnergy) ;
898
899 // Parametrisation
900 } else {
901 eLoss = protonModel->TheValue(G4Proton::Proton(),material,
902 kineticEnergy);
903 }
904 //if(theBarkas) eLoss -= 2.0*BarkasTerm(material, kineticEnergy);
905 }
906
907 // Delta rays energy
908 eLoss -= DeltaRaysEnergy(couple,kineticEnergy,proton_mass_c2) ;
909
910 if(verboseLevel > 2) {
911 G4cout << "pbar E(MeV)= " << kineticEnergy/MeV
912 << " dE/dx(MeV/mm)= " << eLoss*mm/MeV
913 << " for " << material->GetName()
914 << " model: " << protonModel << G4endl;
915 }
916
917 if(eLoss < 0.0) eLoss = 0.0 ;
918
919 return eLoss ;
920 }
921
922
923 // --------------------------------------------------------------------
924 G4double G4hImpactIonisation::DeltaRaysEnergy(const G4MaterialCutsCouple* couple,
925 G4double kineticEnergy,
926 G4double particleMass) const
927 {
928 G4double dLoss = 0.;
929
930 G4double deltaCutNow = cutForDelta[(couple->GetIndex())] ;
931 const G4Material* material = couple->GetMaterial();
932 G4double electronDensity = material->GetElectronDensity();
933 G4double excitationEnergy = material->GetIonisation()->GetMeanExcitationEnergy();
934
935 G4double tau = kineticEnergy / particleMass ;
936 G4double rateMass = electron_mass_c2/particleMass ;
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+2.0*gamma*rateMass+rateMass*rateMass) ;
944
945 // Validity range for delta electron cross section
946 G4double deltaCut = std::max(deltaCutNow, excitationEnergy);
947
948 if ( deltaCut < tMax)
949 {
950 G4double x = deltaCut / tMax ;
951 dLoss = ( beta2 * (x-1.) - std::log(x) ) * twopi_mc2_rcl2 * electronDensity / beta2 ;
952 }
953 return dLoss ;
954 }
955
956
957 // -------------------------------------------------------------------------
958
959 G4VParticleChange* G4hImpactIonisation::PostStepDoIt(const G4Track& track,
960 const G4Step& step)
961 {
962 // Units are expressed in GEANT4 internal units.
963
964 // std::cout << "----- Calling PostStepDoIt ----- " << std::endl;
965
966 aParticleChange.Initialize(track) ;
967 const G4MaterialCutsCouple* couple = track.GetMaterialCutsCouple();
968 const G4DynamicParticle* aParticle = track.GetDynamicParticle() ;
969
970 // Some kinematics
971
972 G4ParticleDefinition* definition = track.GetDefinition();
973 G4double mass = definition->GetPDGMass();
974 G4double kineticEnergy = aParticle->GetKineticEnergy();
975 G4double totalEnergy = kineticEnergy + mass ;
976 G4double pSquare = kineticEnergy *( totalEnergy + mass) ;
977 G4double eSquare = totalEnergy * totalEnergy;
978 G4double betaSquare = pSquare / eSquare;
979 G4ThreeVector particleDirection = aParticle->GetMomentumDirection() ;
980
981 G4double gamma = kineticEnergy / mass + 1.;
982 G4double r = electron_mass_c2 / mass;
983 G4double tMax = 2. * electron_mass_c2 *(gamma*gamma - 1.) / (1. + 2.*gamma*r + r*r);
984
985 // Validity range for delta electron cross section
986 G4double deltaCut = cutForDelta[couple->GetIndex()];
987
988 // This should not be a case
989 if (deltaCut >= tMax)
990 return G4VContinuousDiscreteProcess::PostStepDoIt(track,step);
991
992 G4double xc = deltaCut / tMax;
993 G4double rate = tMax / totalEnergy;
994 rate = rate*rate ;
995 G4double spin = aParticle->GetDefinition()->GetPDGSpin() ;
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 * rate) / (1. + 0.5 * rate) ;
1011 }
1012 else
1013 {
1014 gRej = (1. - betaSquare * x ) * (1. + x/(3.*xc)) +
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(deltaKineticEnergy *
1023 (deltaKineticEnergy + 2. * electron_mass_c2 ));
1024 G4double totalMomentum = std::sqrt(pSquare) ;
1025 G4double cosTheta = deltaKineticEnergy * (totalEnergy + electron_mass_c2) / (deltaTotalMomentum*totalMomentum);
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*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 delta ray
1042 G4DynamicParticle* deltaRay = new G4DynamicParticle;
1043 deltaRay->SetKineticEnergy( deltaKineticEnergy );
1044 deltaRay->SetMomentumDirection(deltaDirection.x(),
1045 deltaDirection.y(),
1046 deltaDirection.z());
1047 deltaRay->SetDefinition(G4Electron::Electron());
1048
1049 // fill aParticleChange
1050 G4double finalKineticEnergy = kineticEnergy - deltaKineticEnergy;
1051 std::size_t totalNumber = 1;
1052
1053 // Atomic relaxation
1054
1055 // ---- MGP ---- Temporary limitation: currently PIXE only for incident protons
1056
1057 std::size_t nSecondaries = 0;
1058 std::vector<G4DynamicParticle*>* secondaryVector = 0;
1059
1060 if (definition == G4Proton::ProtonDefinition())
1061 {
1062 const G4Material* material = couple->GetMaterial();
1063
1064 // Lazy initialization of pixeCrossSectionHandler
1065 if (pixeCrossSectionHandler == 0)
1066 {
1067 // Instantiate pixeCrossSectionHandler with selected shell cross section models
1068 // Ownership of interpolation is transferred to pixeCrossSectionHandler
1069 G4IInterpolator* interpolation = new G4LogLogInterpolator();
1070 pixeCrossSectionHandler =
1071 new G4PixeCrossSectionHandler(interpolation,modelK,modelL,modelM,eMinPixe,eMaxPixe);
1072 G4String fileName("proton");
1073 pixeCrossSectionHandler->LoadShellData(fileName);
1074 // pixeCrossSectionHandler->PrintData();
1075 }
1076
1077 // Select an atom in the current material based on the total shell cross sections
1078 G4int Z = pixeCrossSectionHandler->SelectRandomAtom(material,kineticEnergy);
1079 // std::cout << "G4hImpactIonisation::PostStepDoIt - Z = " << Z << std::endl;
1080
1081 // G4double microscopicCross = MicroscopicCrossSection(*definition,
1082 // kineticEnergy,
1083 // Z, deltaCut);
1084 // G4double crossFromShells = pixeCrossSectionHandler->FindValue(Z,kineticEnergy);
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 based on the individual shell cross sections
1098 G4int shellIndex = pixeCrossSectionHandler->SelectRandomShell(Z,kineticEnergy);
1099
1100 G4AtomicTransitionManager* transitionManager = G4AtomicTransitionManager::Instance();
1101 const G4AtomicShell* atomicShell = transitionManager->Shell(Z,shellIndex);
1102 G4double bindingEnergy = atomicShell->BindingEnergy();
1103
1104 // if (verboseLevel > 1)
1105 // {
1106 // G4cout << "G4hImpactIonisation::PostStepDoIt - Z = "
1107 // << Z
1108 // << ", shell = "
1109 // << shellIndex
1110 // << ", bindingE (keV) = "
1111 // << bindingEnergy/keV
1112 // << G4endl;
1113 // }
1114
1115 // Generate PIXE if binding energy larger than cut for photons or electrons
1116
1117 G4ParticleDefinition* type = 0;
1118
1119 if (finalKineticEnergy >= bindingEnergy)
1120 // && (bindingEnergy >= minGammaEnergy || bindingEnergy >= minElectronEnergy) )
1121 {
1122 // Vacancy in subshell shellIndex; shellId is the subshell identifier in EADL jargon
1123 G4int shellId = atomicShell->ShellId();
1124 // Atomic relaxation: generate secondaries
1125 secondaryVector = atomicDeexcitation.GenerateParticles(Z, shellId);
1126
1127 // ---- Debug ----
1128 //std::cout << "ShellId = "
1129 // <<shellId << " ---- Atomic relaxation secondaries: ---- "
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<nSecondaries; i++)
1139 {
1140 G4DynamicParticle* aSecondary = (*secondaryVector)[i];
1141 if (aSecondary)
1142 {
1143 G4double e = aSecondary->GetKineticEnergy();
1144 type = aSecondary->GetDefinition();
1145
1146 // ---- Debug ----
1147 //if (type == G4Gamma::GammaDefinition())
1148 // {
1149 // std::cout << "Z = " << Z
1150 // << ", shell: " << shellId
1151 // << ", PIXE photon energy (keV) = " << e/keV
1152 // << std::endl;
1153 // }
1154 // ---- End debug ---
1155
1156 if (e < finalKineticEnergy &&
1157 ((type == G4Gamma::Gamma() && e > minGammaEnergy ) ||
1158 (type == G4Electron::Electron() && e > minElectronEnergy )))
1159 {
1160 // Subtract the energy of the emitted secondary from the primary
1161 finalKineticEnergy -= e;
1162 totalNumber++;
1163 // ---- Debug ----
1164 //if (type == G4Gamma::GammaDefinition())
1165 // {
1166 // std::cout << "Z = " << Z
1167 // << ", shell: " << shellId
1168 // << ", PIXE photon energy (keV) = " << e/keV
1169 // << std::endl;
1170 // }
1171 // ---- End debug ---
1172 }
1173 else
1174 {
1175 // The atomic relaxation product has energy below the cut
1176 // ---- Debug ----
1177 // if (type == G4Gamma::GammaDefinition())
1178 // {
1179 // std::cout << "Z = " << Z
1180 //
1181 // << ", PIXE photon energy = " << e/keV
1182 // << " keV below threshold " << minGammaEnergy/keV << " keV"
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*particleDirection.x() - deltaTotalMomentum*deltaDirection.x();
1204 G4double finalPy = totalMomentum*particleDirection.y() - deltaTotalMomentum*deltaDirection.y();
1205 G4double finalPz = totalMomentum*particleDirection.z() - deltaTotalMomentum*deltaDirection.z();
1206 G4double finalMomentum = std::sqrt(finalPx*finalPx + finalPy*finalPy + finalPz*finalPz) ;
1207 finalPx /= finalMomentum;
1208 finalPy /= finalMomentum;
1209 finalPz /= finalMomentum;
1210
1211 aParticleChange.ProposeMomentumDirection( finalPx,finalPy,finalPz );
1212 }
1213 else
1214 {
1215 eDeposit = finalKineticEnergy;
1216 finalKineticEnergy = 0.;
1217 aParticleChange.ProposeMomentumDirection(particleDirection.x(),
1218 particleDirection.y(),
1219 particleDirection.z());
1220 if(!aParticle->GetDefinition()->GetProcessManager()->
1221 GetAtRestProcessVector()->size())
1222 aParticleChange.ProposeTrackStatus(fStopAndKill);
1223 else
1224 aParticleChange.ProposeTrackStatus(fStopButAlive);
1225 }
1226
1227 aParticleChange.ProposeEnergy(finalKineticEnergy);
1228 aParticleChange.ProposeLocalEnergyDeposit (eDeposit);
1229 aParticleChange.SetNumberOfSecondaries((G4int)totalNumber);
1230 aParticleChange.AddSecondary(deltaRay);
1231
1232 // ---- Debug ----
1233 // std::cout << "RDHadronIonisation - finalKineticEnergy (MeV) = "
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 < nSecondaries; l++)
1247 {
1248 G4DynamicParticle* secondary = (*secondaryVector)[l];
1249 if (secondary) aParticleChange.AddSecondary(secondary);
1250
1251 // ---- Debug ----
1252 //if (secondary != 0)
1253 // {
1254 // if (secondary->GetDefinition() == G4Gamma::GammaDefinition())
1255 // {
1256 // G4double eX = secondary->GetKineticEnergy();
1257 // std::cout << " PIXE photon of energy " << eX/keV
1258 // << " keV added to ParticleChange; total number of secondaries is " << totalNumber
1259 // << std::endl;
1260 // }
1261 //}
1262 // ---- End debug ---
1263
1264 }
1265 delete secondaryVector;
1266 }
1267
1268 return G4VContinuousDiscreteProcess::PostStepDoIt(track,step);
1269 }
1270
1271 // -------------------------------------------------------------------------
1272
1273 G4double G4hImpactIonisation::ComputeDEDX(const G4ParticleDefinition* aParticle,
1274 const G4MaterialCutsCouple* couple,
1275 G4double kineticEnergy)
1276 {
1277 const G4Material* material = couple->GetMaterial();
1278 G4Proton* proton = G4Proton::Proton();
1279 G4AntiProton* antiproton = G4AntiProton::AntiProton();
1280 G4double dedx = 0.;
1281
1282 G4double tScaled = kineticEnergy * proton_mass_c2 / (aParticle->GetPDGMass()) ;
1283 charge = aParticle->GetPDGCharge() ;
1284
1285 if( charge > 0.)
1286 {
1287 if (tScaled > protonHighEnergy)
1288 {
1289 dedx = G4EnergyLossTables::GetDEDX(proton,tScaled,couple) ;
1290 }
1291 else
1292 {
1293 dedx = ProtonParametrisedDEDX(couple,tScaled) ;
1294 }
1295 }
1296 else
1297 {
1298 if (tScaled > antiprotonHighEnergy)
1299 {
1300 dedx = G4EnergyLossTables::GetDEDX(antiproton,tScaled,couple);
1301 }
1302 else
1303 {
1304 dedx = AntiProtonParametrisedDEDX(couple,tScaled) ;
1305 }
1306 }
1307 dedx *= theIonEffChargeModel->TheValue(aParticle, material, kineticEnergy) ;
1308
1309 return dedx ;
1310 }
1311
1312
1313 G4double G4hImpactIonisation::BarkasTerm(const G4Material* material,
1314 G4double kineticEnergy) const
1315 //Function to compute the Barkas term for protons:
1316 //
1317 //Ref. Z_1^3 effect in the stopping power of matter for charged particles
1318 // J.C Ashley and R.H.Ritchie
1319 // Physical review B Vol.5 No.7 1 April 1972 pagg. 2393-2397
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_c2 ;
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 = material->GetElementVector();
1382 G4int numberOfElements = (G4int)material->GetNumberOfElements();
1383
1384 for (G4int i = 0; i<numberOfElements; i++) {
1385
1386 //AMaterial = (*theElementVector)[i]->GetA()*mole/g;
1387 ZMaterial = (*theElementVector)[i]->GetZ();
1388
1389 G4double X = 137.0 * 137.0 * beta2 / ZMaterial;
1390
1391 // Variables to compute L_1
1392 G4double Eta0Chi = 0.8;
1393 G4double EtaChi = Eta0Chi * ( 1.0 + 6.02*std::pow( ZMaterial,-1.19 ) );
1394 G4double W = ( EtaChi * std::pow( ZMaterial,1.0/6.0 ) ) / std::sqrt(X);
1395 G4double FunctionOfW = FTable[46][1]*FTable[46][0]/W ;
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] - FTable[j-1][1]) * (W - FTable[j-1][0])
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(ZMaterial * X) * X);
1416 }
1417
1418 BTerm *= twopi_mc2_rcl2 * (material->GetElectronDensity()) / beta2 ;
1419
1420 return BTerm;
1421 }
1422
1423
1424
1425 G4double G4hImpactIonisation::BlochTerm(const G4Material* material,
1426 G4double kineticEnergy,
1427 G4double cSquare) const
1428 //Function to compute the Bloch term for protons:
1429 //
1430 //Ref. Z_1^3 effect in the stopping power of matter for charged particles
1431 // J.C Ashley and R.H.Ritchie
1432 // Physical review B Vol.5 No.7 1 April 1972 pagg. 2393-2397
1433 //
1434 {
1435 G4double eLoss = 0.0 ;
1436 G4double gamma = 1.0 + kineticEnergy / proton_mass_c2 ;
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_rcl2 *
1453 (material->GetElectronDensity()) / beta2 ;
1454
1455 return eLoss;
1456 }
1457
1458
1459
1460 G4double G4hImpactIonisation::ElectronicLossFluctuation(
1461 const G4DynamicParticle* particle,
1462 const G4MaterialCutsCouple* couple,
1463 G4double meanLoss,
1464 G4double step) const
1465 // calculate actual loss from the mean loss
1466 // The model used to get the fluctuation is essentially the same
1467 // as in Glandz in Geant3.
1468 {
1469 // data members to speed up the fluctuation calculation
1470 // G4int imat ;
1471 // G4double f1Fluct,f2Fluct,e1Fluct,e2Fluct,rateFluct,ipotFluct;
1472 // G4double e1LogFluct,e2LogFluct,ipotLogFluct;
1473
1474 static const G4double minLoss = 1.*eV ;
1475 static const G4double kappa = 10. ;
1476 static const G4double theBohrBeta2 = 50.0 * keV/proton_mass_c2 ;
1477
1478 const G4Material* material = couple->GetMaterial();
1479 G4int imaterial = couple->GetIndex() ;
1480 G4double ipotFluct = material->GetIonisation()->GetMeanExcitationEnergy() ;
1481 G4double electronDensity = material->GetElectronDensity() ;
1482 G4double zeff = electronDensity/(material->GetTotNbOfAtomsPerVolume()) ;
1483
1484 // get particle data
1485 G4double tkin = particle->GetKineticEnergy();
1486 G4double particleMass = particle->GetMass() ;
1487 G4double deltaCutInKineticEnergyNow = cutForDelta[imaterial];
1488
1489 // shortcut for very very small loss
1490 if(meanLoss < minLoss) return meanLoss ;
1491
1492 // Validity range for delta electron cross section
1493 G4double threshold = std::max(deltaCutInKineticEnergyNow,ipotFluct);
1494 G4double loss, siga;
1495
1496 G4double rmass = electron_mass_c2/particleMass;
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/(1.+2.*tau1*rmass+rmass*rmass);
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*ipotFluct )
1508 {
1509 siga = tMax * (1.0-0.5*beta2) * step * twopi_mc2_rcl2
1510 * electronDensity / beta2 ;
1511
1512 // High velocity or negatively charged particle
1513 if( beta2 > 3.0*theBohrBeta2*zeff || charge < 0.0) {
1514 siga = std::sqrt( siga * chargeSquare ) ;
1515
1516 // Low velocity - additional ion charge fluctuations according to
1517 // Q.Yang et al., NIM B61(1991)149-155.
1518 } else {
1519 G4double chu = theIonChuFluctuationModel->TheValue(particle, material);
1520 G4double yang = theIonYangFluctuationModel->TheValue(particle, material);
1521 siga = std::sqrt( siga * (chargeSquare * chu + yang)) ;
1522 }
1523
1524 do {
1525 loss = G4RandGauss::shoot(meanLoss,siga);
1526 } while (loss < 0. || loss > 2.0*meanLoss);
1527 return loss;
1528 }
1529
1530 // Non Gaussian fluctuation
1531 static const G4double probLim = 0.01 ;
1532 static const G4double sumaLim = -std::log(probLim) ;
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,alfa1,ea,sea;
1540 G4double dp3;
1541
1542 G4double f1Fluct = material->GetIonisation()->GetF1fluct();
1543 G4double f2Fluct = material->GetIonisation()->GetF2fluct();
1544 G4double e1Fluct = material->GetIonisation()->GetEnergy1fluct();
1545 G4double e2Fluct = material->GetIonisation()->GetEnergy2fluct();
1546 G4double e1LogFluct = material->GetIonisation()->GetLogEnergy1fluct();
1547 G4double e2LogFluct = material->GetIonisation()->GetLogEnergy2fluct();
1548 G4double rateFluct = material->GetIonisation()->GetRateionexcfluct();
1549 G4double ipotLogFluct= material->GetIonisation()->GetLogMeanExcEnergy();
1550
1551 w1 = tMax/ipotFluct;
1552 w2 = std::log(2.*electron_mass_c2*tau2);
1553
1554 C = meanLoss*(1.-rateFluct)/(w2-ipotLogFluct-beta2);
1555
1556 a1 = C*f1Fluct*(w2-e1LogFluct-beta2)/e1Fluct;
1557 a2 = C*f2Fluct*(w2-e2LogFluct-beta2)/e2Fluct;
1558 a3 = rateFluct*meanLoss*(tMax-ipotFluct)/(ipotFluct*tMax*std::log(w1));
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 small Step
1569 {
1570 e0 = material->GetIonisation()->GetEnergy0fluct();
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::shoot(a3,siga)+0.5));
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(tMax/e0));
1594
1595 if(a3>alim)
1596 {
1597 siga=std::sqrt(a3) ;
1598 p3 = std::max(0,int(G4RandGauss::shoot(a3,siga)+0.5));
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./(1.-w*G4UniformRand()) ;
1616 loss *= e0*corrfac ;
1617 }
1618 }
1619 }
1620
1621 else // not so small Step
1622 {
1623 // excitation type 1
1624 if(a1>alim)
1625 {
1626 siga=std::sqrt(a1) ;
1627 p1 = std::max(0,G4int(G4RandGauss::shoot(a1,siga)+0.5));
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(a2,siga)+0.5));
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())*e2Fluct;
1646 else if (loss>0.)
1647 loss += (1.-2.*G4UniformRand())*e1Fluct;
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::shoot(a3,siga)+0.5));
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)+dp3);
1669 namean = G4float(p3)*rfac;
1670 sa = G4float(nmaxCont1)*rfac;
1671 na = G4RandGauss::shoot(namean,sa);
1672 if (na > 0.)
1673 {
1674 alfa = w1*G4float(nmaxCont2+p3)/
1675 (w1*G4float(nmaxCont2)+G4float(p3));
1676 alfa1 = alfa*std::log(alfa)/(alfa-1.);
1677 ea = na*ipotFluct*alfa1;
1678 sea = ipotFluct*std::sqrt(na*(alfa-alfa1*alfa1));
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/(1.-w*G4UniformRand());
1689 }
1690 }
1691 loss += lossc;
1692 }
1693 }
1694
1695 return loss ;
1696 }
1697
1698
1699
1700 void G4hImpactIonisation::SetCutForSecondaryPhotons(G4double cut)
1701 {
1702 minGammaEnergy = cut;
1703 }
1704
1705
1706
1707 void G4hImpactIonisation::SetCutForAugerElectrons(G4double cut)
1708 {
1709 minElectronEnergy = cut;
1710 }
1711
1712
1713
1714 void G4hImpactIonisation::ActivateAugerElectronProduction(G4bool val)
1715 {
1716 atomicDeexcitation.ActivateAugerElectronProduction(val);
1717 }
1718
1719
1720
1721 void G4hImpactIonisation::PrintInfoDefinition() const
1722 {
1723 G4String comments = " Knock-on electron cross sections . ";
1724 comments += "\n Good description above the mean excitation energy.\n";
1725 comments += " Delta ray energy sampled from differential Xsection.";
1726
1727 G4cout << G4endl << GetProcessName() << ": " << comments
1728 << "\n PhysicsTables from " << LowestKineticEnergy / eV << " eV "
1729 << " to " << HighestKineticEnergy / TeV << " TeV "
1730 << " in " << TotBin << " bins."
1731 << "\n Electronic stopping power model is "
1732 << protonTable
1733 << "\n from " << protonLowEnergy / keV << " keV "
1734 << " to " << protonHighEnergy / MeV << " MeV " << "." << G4endl ;
1735 G4cout << "\n Parametrisation model for antiprotons is "
1736 << antiprotonTable
1737 << "\n from " << antiprotonLowEnergy / keV << " keV "
1738 << " to " << antiprotonHighEnergy / MeV << " MeV " << "." << G4endl ;
1739 if(theBarkas){
1740 G4cout << " Parametrization of the Barkas effect is switched on."
1741 << G4endl ;
1742 }
1743 if(nStopping) {
1744 G4cout << " Nuclear stopping power model is " << theNuclearTable
1745 << G4endl ;
1746 }
1747
1748 G4bool printHead = true;
1749
1750 const G4ProductionCutsTable* theCoupleTable=
1751 G4ProductionCutsTable::GetProductionCutsTable();
1752 G4int numOfCouples = (G4int)theCoupleTable->GetTableSize();
1753
1754 // loop for materials
1755
1756 for (G4int j=0 ; j < numOfCouples; ++j) {
1757
1758 const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(j);
1759 const G4Material* material= couple->GetMaterial();
1760 G4double deltaCutNow = cutForDelta[(couple->GetIndex())] ;
1761 G4double excitationEnergy = material->GetIonisation()->GetMeanExcitationEnergy();
1762
1763 if(excitationEnergy > deltaCutNow) {
1764 if(printHead) {
1765 printHead = false ;
1766
1767 G4cout << " material min.delta energy(keV) " << G4endl;
1768 G4cout << G4endl;
1769 }
1770
1771 G4cout << std::setw(20) << material->GetName()
1772 << std::setw(15) << excitationEnergy/keV << G4endl;
1773 }
1774 }
1775 }
1776