source: trunk/source/processes/electromagnetic/pii/src/G4hImpactIonisation.cc@ 1357

Last change on this file since 1357 was 1350, checked in by garnier, 15 years ago

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