source: trunk/source/processes/electromagnetic/lowenergy/src/G4hLowEnergyIonisation.cc@ 1346

Last change on this file since 1346 was 1315, checked in by garnier, 15 years ago

update geant4-09-04-beta-cand-01 interfaces-V09-03-09 vis-V09-03-08
File size: 66.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// GEANT 4 class implementation file
29//
30// History: based on object model of
31// 2nd December 1995, G.Cosmo
32// ---------- G4hLowEnergyIonisation physics process -------
33// by Vladimir Ivanchenko, 14 July 1999
34// was made on the base of G4hIonisation class
35// developed by Laszlo Urban
36// ************************************************************
37// It is the extention of the ionisation process for the slow
38// charged hadrons.
39// ************************************************************
40// 28 July 1999 V.Ivanchenko cleen up
41// 17 August 1999 G.Mancinelli added ICRU parametrisations for protons
42// 20 August 1999 G.Mancinelli added ICRU tables for alpha
43// 31 August 1999 V.Ivanchenko update and cleen up
44// 30 Sept. 1999 V.Ivanchenko minor upgrade
45// 12 Dec. 1999 S. Chauvie added Barkas correction
46// 19 Jan. 2000 V.Ivanchenko minor changing in Barkas corrections
47// 02 April 2000 S. Chauvie linearization of Barkas effect
48// 03 April 2000 V.Ivanchenko Nuclear Stopping power for antiprotons
49// 23 May 2000 MG Pia Clean up for QAO model
50// 24 May 2000 MG Pia Code properly indented to improve legibility
51// 17 July 2000 V.Ivanchenko Bug in scaling AlongStepDoIt method
52// 25 July 2000 V.Ivanchenko New design iteration
53// 17 August 2000 V.Ivanchenko Add ion fluctuation models
54// 18 August 2000 V.Ivanchenko Bug fixed in GetConstrain
55// 22 August 2000 V.Ivanchenko Insert paramStepLimit and
56// reorganise access to Barkas and Bloch terms
57// 04 Sept. 2000 V.Ivanchenko rename fluctuations
58// 05 Sept. 2000 V.Ivanchenko clean up
59// 03 Oct. 2000 V.Ivanchenko CodeWizard clean up
60// 03 Nov. 2000 V.Ivanchenko MinKineticEnergy=LowestKineticEnergy=10eV
61// 05 Nov. 2000 MG Pia - Removed const cast previously introduced to get
62// the code compiled (const G4Material* now introduced in
63// electromagnetic/utils utils-V02-00-03 tag)
64// (this is going back and forth, to cope with Michel's
65// utils tag not being accepted yet by system testing)
66// 21 Nov. 2000 V.Ivanchenko Fix a problem in fluctuations
67// 23 Nov. 2000 V.Ivanchenko Ion type fluctuations only for charge>0
68// 10 May 2001 V.Ivanchenko Clean up againist Linux compilation with -Wall
69// 23 May 2001 V.Ivanchenko Minor fix in PostStepDoIt
70// 07 June 2001 V.Ivanchenko Clean up AntiProtonDEDX + add print out
71// 18 June 2001 V.Ivanchenko Cleanup print out
72// 18 Oct. 2001 V.Ivanchenko Add fluorescence
73// 30 Oct. 2001 V.Ivanchenko Add minGammaEnergy and minElectronEnergy
74// 07 Dec 2001 V.Ivanchenko Add SetFluorescence method
75// 15 Feb 2002 V.Ivanchenko Fix problem of Generic Ions
76// 25 Mar 2002 V.Ivanchenko Fix problem of fluorescence below threshold
77// 28 Mar 2002 V.Ivanchenko Set fluorescence off by default
78// 09 Apr 2002 V.Ivanchenko Fix table problem of GenericIons
79// 28 May 2002 V.Ivanchenko Remove flag fStopAndKill
80// 31 May 2002 V.Ivanchenko Add path of Fluo + Auger cuts to
81// AtomicDeexcitation
82// 03 Jun 2002 MGP Restore fStopAndKill
83// 10 Jun 2002 V.Ivanchenko Restore fStopButAlive
84// 12 Jun 2002 V.Ivanchenko Fix in fluctuations - if tmax<2*Ipot Gaussian
85// fluctuations enables
86// 20 Sept 2002 V.Ivanchenko Clean up energy ranges for models
87// 07 Oct 2002 V.Ivanchenko Clean up initialisation of fluorescence
88// 28 Oct 2002 V.Ivanchenko Optimal binning for dE/dx
89// 10 Dec 2002 V.Ivanchenko antiProtonLowEnergy -> 25 keV, QEG model below
90// 21 Jan 2003 V.Ivanchenko Cut per region
91// 10 Mar 2003 V.Ivanchenko Use SubTypes for ions
92// 12 Apr 2003 V.Ivanchenko Cut per region for fluo AlongStep
93// 18 Apr 2003 V.Ivanchenko finalRange redefinition
94// 26 Apr 2003 V.Ivanchenko fix for stepLimit
95// 30 Mar 2004 S.Saliceti add shellCS data member and expFlag variable,
96// atom total cross section for the Empiric Model
97// 28 May 2004 V.Ivanchenko fix for ionisation of antiprotons in complex materials
98// 30 Aug 2004 V.Ivanchenko use energy limit for parameterisation from model
99// 03 Oct 2005 V.Ivanchenko change logic of definition of high energy limit for
100// parametrised proton model: min(user value, model limit)
101// 26 Jan 2005 S. Chauvie added PrintInfoDefinition() for antiproton
102// 30 Sep 2009 A.Mantero Removed dependencies to old shell Ionisation XS models
103// 07 Jun 2010 Code Celaning for June beta Release
104// -----------------------------------------------------------------------
105
106//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
107//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
108
109#include "G4hLowEnergyIonisation.hh"
110#include "globals.hh"
111#include "G4ios.hh"
112#include "Randomize.hh"
113#include "G4Poisson.hh"
114#include "G4UnitsTable.hh"
115#include "G4EnergyLossTables.hh"
116#include "G4Material.hh"
117#include "G4DynamicParticle.hh"
118#include "G4ParticleDefinition.hh"
119#include "G4AtomicDeexcitation.hh"
120#include "G4AtomicTransitionManager.hh"
121#include "G4ShellVacancy.hh"
122#include "G4VhShellCrossSection.hh"
123#include "G4VEMDataSet.hh"
124#include "G4EMDataSet.hh"
125#include "G4CompositeEMDataSet.hh"
126#include "G4Gamma.hh"
127#include "G4LogLogInterpolation.hh"
128#include "G4SemiLogInterpolation.hh"
129#include "G4ProcessManager.hh"
130#include "G4ProductionCutsTable.hh"
131#include "G4teoCrossSection.hh"
132#include "G4empCrossSection.hh"
133//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
134
135G4hLowEnergyIonisation::G4hLowEnergyIonisation(const G4String& processName)
136 : G4hLowEnergyLoss(processName),
137 theBetheBlochModel(0),
138 theProtonModel(0),
139 theAntiProtonModel(0),
140 theIonEffChargeModel(0),
141 theNuclearStoppingModel(0),
142 theIonChuFluctuationModel(0),
143 theIonYangFluctuationModel(0),
144 theProtonTable("ICRU_R49p"),
145 theAntiProtonTable("ICRU_R49p"),
146 theNuclearTable("ICRU_R49"),
147 nStopping(true),
148 theBarkas(true),
149 theMeanFreePathTable(0),
150 paramStepLimit (0.005),
151 shellVacancy(0),
152 shellCS(0),
153 theFluo(false)
154{
155 InitializeMe();
156}
157
158//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
159
160void G4hLowEnergyIonisation::InitializeMe()
161{
162 LowestKineticEnergy = 10.0*eV ;
163 HighestKineticEnergy = 100.0*GeV ;
164 MinKineticEnergy = 10.0*eV ;
165 TotBin = 360 ;
166 protonLowEnergy = 1.*keV ;
167 protonHighEnergy = 100.*MeV ;
168 antiProtonLowEnergy = 25.*keV ;
169 antiProtonHighEnergy = 2.*MeV ;
170 minGammaEnergy = 25.*keV;
171 minElectronEnergy = 25.*keV;
172 verboseLevel = 0;
173
174 shellCS = new G4teoCrossSection("analytical");
175
176}
177
178//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
179
180G4hLowEnergyIonisation::~G4hLowEnergyIonisation()
181{
182 if (theMeanFreePathTable) {
183 theMeanFreePathTable->clearAndDestroy();
184 delete theMeanFreePathTable;
185 }
186 if(theBetheBlochModel)delete theBetheBlochModel;
187 if(theProtonModel)delete theProtonModel;
188 if(theAntiProtonModel)delete theAntiProtonModel;
189 if(theNuclearStoppingModel)delete theNuclearStoppingModel;
190 if(theIonEffChargeModel)delete theIonEffChargeModel;
191 if(theIonChuFluctuationModel)delete theIonChuFluctuationModel;
192 if(theIonYangFluctuationModel)delete theIonYangFluctuationModel;
193 if(shellVacancy) delete shellVacancy;
194 if(shellCS) delete shellCS;
195 cutForDelta.clear();
196 G4int length = zFluoDataVector.size();
197 if(length) {
198 for(G4int i=0; i<length; i++) {
199 delete zFluoDataVector[i];
200 }
201 zFluoDataVector.clear();
202 }
203}
204
205//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
206
207void G4hLowEnergyIonisation::SetElectronicStoppingPowerModel(
208 const G4ParticleDefinition* aParticle,
209 const G4String& dedxTable)
210 // This method defines the ionisation parametrisation method via its name
211{
212 if(0 < aParticle->GetPDGCharge()) {
213 SetProtonElectronicStoppingPowerModel(dedxTable) ;
214 } else {
215 SetAntiProtonElectronicStoppingPowerModel(dedxTable) ;
216 }
217}
218
219//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
220
221void G4hLowEnergyIonisation::InitializeParametrisation()
222
223{
224 // Define models for parametrisation of electronic energy losses
225 theBetheBlochModel = new G4hBetheBlochModel("Bethe-Bloch") ;
226 theProtonModel = new G4hParametrisedLossModel(theProtonTable) ;
227 protonHighEnergy = std::min(protonHighEnergy,theProtonModel->HighEnergyLimit(0, 0));
228 theAntiProtonModel = new G4QAOLowEnergyLoss(theAntiProtonTable) ;
229 theNuclearStoppingModel = new G4hNuclearStoppingModel(theNuclearTable) ;
230 theIonEffChargeModel = new G4hIonEffChargeSquare("Ziegler1988") ;
231 theIonChuFluctuationModel = new G4IonChuFluctuationModel("Chu") ;
232 theIonYangFluctuationModel = new G4IonYangFluctuationModel("Yang") ;
233}
234
235//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
236
237void G4hLowEnergyIonisation::BuildPhysicsTable(
238 const G4ParticleDefinition& aParticleType)
239
240 // just call BuildLossTable+BuildLambdaTable
241{
242 if(verboseLevel > 0) {
243 G4cout << "G4hLowEnergyIonisation::BuildPhysicsTable for "
244 << aParticleType.GetParticleName()
245 << " mass(MeV)= " << aParticleType.GetPDGMass()/MeV
246 << " charge= " << aParticleType.GetPDGCharge()/eplus
247 << " type= " << aParticleType.GetParticleType()
248 << G4endl;
249
250 if(verboseLevel > 1) {
251 G4ProcessVector* pv = aParticleType.GetProcessManager()->GetProcessList();
252
253 G4cout << " 0: " << (*pv)[0]->GetProcessName() << " " << (*pv)[0]
254 << " 1: " << (*pv)[1]->GetProcessName() << " " << (*pv)[1]
255 // << " 2: " << (*pv)[2]->GetProcessName() << " " << (*pv)[2]
256 << G4endl;
257 G4cout << "ionModel= " << theIonEffChargeModel
258 << " MFPtable= " << theMeanFreePathTable
259 << " iniMass= " << initialMass
260 << G4endl;
261 }
262 }
263
264 if(aParticleType.GetParticleType() == "nucleus" &&
265 aParticleType.GetParticleName() != "GenericIon" &&
266 aParticleType.GetParticleSubType() == "generic")
267 {
268
269 G4EnergyLossTables::Register(&aParticleType,
270 theDEDXpTable,
271 theRangepTable,
272 theInverseRangepTable,
273 theLabTimepTable,
274 theProperTimepTable,
275 LowestKineticEnergy, HighestKineticEnergy,
276 proton_mass_c2/aParticleType.GetPDGMass(),
277 TotBin);
278
279 return;
280 }
281
282 if( !CutsWhereModified() && theLossTable) return;
283
284 InitializeParametrisation() ;
285 G4Proton* theProton = G4Proton::Proton();
286 G4AntiProton* theAntiProton = G4AntiProton::AntiProton();
287
288 charge = aParticleType.GetPDGCharge()/eplus;
289 chargeSquare = charge*charge ;
290
291 const G4ProductionCutsTable* theCoupleTable=
292 G4ProductionCutsTable::GetProductionCutsTable();
293 size_t numOfCouples = theCoupleTable->GetTableSize();
294
295 cutForDelta.clear();
296 cutForGamma.clear();
297
298 for (size_t j=0; j<numOfCouples; j++) {
299
300 // get material parameters needed for the energy loss calculation
301 const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(j);
302 const G4Material* material= couple->GetMaterial();
303
304 // the cut cannot be below lowest limit
305 G4double tCut = (*(theCoupleTable->GetEnergyCutsVector(1)))[j];
306 if(tCut > HighestKineticEnergy) tCut = HighestKineticEnergy;
307
308 G4double excEnergy = material->GetIonisation()->GetMeanExcitationEnergy();
309
310 tCut = std::max(tCut,excEnergy);
311 cutForDelta.push_back(tCut);
312
313 // the cut cannot be below lowest limit
314 tCut = (*(theCoupleTable->GetEnergyCutsVector(0)))[j];
315 if(tCut > HighestKineticEnergy) tCut = HighestKineticEnergy;
316 tCut = std::max(tCut,minGammaEnergy);
317 cutForGamma.push_back(tCut);
318 }
319
320 if(verboseLevel > 0) {
321 G4cout << "Cuts are defined " << G4endl;
322 }
323
324 if(0.0 < charge)
325 {
326 {
327 BuildLossTable(*theProton) ;
328
329// The following vector has a fixed dimension (see src/G4hLowEnergyLoss.cc for more details)
330// It happended in the past that caused memory corruption errors. The problem is still pending, even if temporary solved
331// G4cout << "[NOTE]: __LINE__=" << __LINE__ << ", aParticleType=" << aParticleType.GetParticleName() << ", theProton=" << theProton << ", theLossTable=" << theLossTable << ", CounterOfpProcess=" << CounterOfpProcess << G4endl;
332
333 RecorderOfpProcess[CounterOfpProcess] = theLossTable ;
334 CounterOfpProcess++;
335 }
336 } else {
337 {
338 BuildLossTable(*theAntiProton) ;
339
340// The following vector has a fixed dimension (see src/G4hLowEnergyLoss.cc for more details)
341// It happended in the past that caused memory corruption errors. The problem is still pending, even if temporary solved
342// G4cout << "[NOTE]: __LINE__=" << __LINE__ << ", aParticleType=" << aParticleType.GetParticleName() << ", theAntiProton=" << theAntiProton << ", theLossTable=" << theLossTable << ", CounterOfpbarProcess=" << CounterOfpbarProcess << G4endl;
343
344 RecorderOfpbarProcess[CounterOfpbarProcess] = theLossTable ;
345 CounterOfpbarProcess++;
346 }
347 }
348
349 if(verboseLevel > 0) {
350 G4cout << "G4hLowEnergyIonisation::BuildPhysicsTable: "
351 << "Loss table is built "
352// << theLossTable
353 << G4endl;
354 }
355
356 BuildLambdaTable(aParticleType) ;
357 BuildDataForFluorescence(aParticleType);
358
359 if(verboseLevel > 1) {
360 G4cout << (*theMeanFreePathTable) << G4endl;
361 }
362
363 if(verboseLevel > 0) {
364 G4cout << "G4hLowEnergyIonisation::BuildPhysicsTable: "
365 << "DEDX table will be built "
366// << theDEDXpTable << " " << theDEDXpbarTable
367// << " " << theRangepTable << " " << theRangepbarTable
368 << G4endl;
369 }
370
371 BuildDEDXTable(aParticleType) ;
372
373 if(verboseLevel > 1) {
374 G4cout << (*theDEDXpTable) << G4endl;
375 }
376
377 if((&aParticleType == theProton) || (&aParticleType == theAntiProton)) PrintInfoDefinition() ;
378
379 if(verboseLevel > 0) {
380 G4cout << "G4hLowEnergyIonisation::BuildPhysicsTable: end for "
381 << aParticleType.GetParticleName() << G4endl;
382 }
383}
384
385//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
386
387void G4hLowEnergyIonisation::BuildLossTable(
388 const G4ParticleDefinition& aParticleType)
389{
390
391 // Initialisation
392 G4double lowEdgeEnergy , ionloss, ionlossBB, paramB ;
393 G4double lowEnergy, highEnergy;
394 G4Proton* theProton = G4Proton::Proton();
395
396 if(aParticleType == *theProton) {
397 lowEnergy = protonLowEnergy ;
398 highEnergy = protonHighEnergy ;
399 charge = 1.0 ;
400 } else {
401 lowEnergy = antiProtonLowEnergy ;
402 highEnergy = antiProtonHighEnergy ;
403 charge = -1.0 ;
404 }
405 chargeSquare = 1.0 ;
406
407 const G4ProductionCutsTable* theCoupleTable=
408 G4ProductionCutsTable::GetProductionCutsTable();
409 size_t numOfCouples = theCoupleTable->GetTableSize();
410
411 if ( theLossTable) {
412 theLossTable->clearAndDestroy();
413 delete theLossTable;
414 }
415
416 theLossTable = new G4PhysicsTable(numOfCouples);
417
418 // loop for materials
419 for (size_t j=0; j<numOfCouples; j++) {
420
421 // create physics vector and fill it
422 G4PhysicsLogVector* aVector = new G4PhysicsLogVector(LowestKineticEnergy,
423 HighestKineticEnergy,
424 TotBin);
425
426 // get material parameters needed for the energy loss calculation
427 const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(j);
428 const G4Material* material= couple->GetMaterial();
429
430 if ( charge > 0.0 ) {
431 ionloss = ProtonParametrisedDEDX(couple,highEnergy) ;
432 } else {
433 ionloss = AntiProtonParametrisedDEDX(couple,highEnergy) ;
434 }
435
436 ionlossBB = theBetheBlochModel->TheValue(&aParticleType,material,highEnergy) ;
437 ionlossBB -= DeltaRaysEnergy(couple,highEnergy,proton_mass_c2) ;
438
439
440 paramB = ionloss/ionlossBB - 1.0 ;
441
442 // now comes the loop for the kinetic energy values
443 for (G4int i = 0 ; i < TotBin ; i++) {
444 lowEdgeEnergy = aVector->GetLowEdgeEnergy(i) ;
445
446 // low energy part for this material, parametrised energy loss formulae
447 if ( lowEdgeEnergy < highEnergy ) {
448
449 if ( charge > 0.0 ) {
450 ionloss = ProtonParametrisedDEDX(couple,lowEdgeEnergy) ;
451 } else {
452 ionloss = AntiProtonParametrisedDEDX(couple,lowEdgeEnergy) ;
453 }
454
455 } else {
456
457 // high energy part for this material, Bethe-Bloch formula
458 ionloss = theBetheBlochModel->TheValue(theProton,material,
459 lowEdgeEnergy) ;
460
461 ionloss -= DeltaRaysEnergy(couple,lowEdgeEnergy,proton_mass_c2) ;
462
463 ionloss *= (1.0 + paramB*highEnergy/lowEdgeEnergy) ;
464 }
465
466 // now put the loss into the vector
467 if(verboseLevel > 1) {
468 G4cout << "E(MeV)= " << lowEdgeEnergy/MeV
469 << " dE/dx(MeV/mm)= " << ionloss*mm/MeV
470 << " in " << material->GetName() << G4endl;
471 }
472 aVector->PutValue(i,ionloss) ;
473 }
474 // Insert vector for this material into the table
475 theLossTable->insert(aVector) ;
476 }
477}
478
479//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
480
481void G4hLowEnergyIonisation::BuildDataForFluorescence(
482 const G4ParticleDefinition& aParticleType)
483{
484
485 if(verboseLevel > 1) {
486 G4cout << "G4hLowEnergyIonisation::BuildDataForFluorescence for "
487 << aParticleType.GetParticleName() << " is started" << G4endl;
488 }
489
490 // fill data for fluorescence
491
492 deexcitationManager.SetCutForSecondaryPhotons(minGammaEnergy);
493 deexcitationManager.SetCutForAugerElectrons(minElectronEnergy);
494
495 G4double mass = aParticleType.GetPDGMass();
496 const G4ProductionCutsTable* theCoupleTable=
497 G4ProductionCutsTable::GetProductionCutsTable();
498 size_t numOfCouples = theCoupleTable->GetTableSize();
499
500 if (shellVacancy != 0) delete shellVacancy;
501 shellVacancy = new G4ShellVacancy();
502 G4DataVector* ksi = 0;
503 G4DataVector* ksi1 = 0;
504 G4DataVector* energy = 0;
505 G4DataVector* energy1 = 0;
506 size_t binForFluo = TotBin/10;
507 G4int length = zFluoDataVector.size();
508 if(length > 0) {
509 for(G4int i=0; i<length; i++) {
510 G4VEMDataSet* x = zFluoDataVector[i];
511 delete x;
512 }
513 zFluoDataVector.clear();
514 }
515
516 G4PhysicsLogVector* bVector = new G4PhysicsLogVector(LowestKineticEnergy,
517 HighestKineticEnergy,
518 binForFluo);
519 const G4AtomicTransitionManager* transitionManager =
520 G4AtomicTransitionManager::Instance();
521
522 G4double bindingEnergy;
523 // G4double x;
524 // G4double y;
525
526 // loop for materials
527 for (size_t j=0; j<numOfCouples; j++) {
528
529 // get material parameters needed for the energy loss calculation
530 const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(j);
531 const G4Material* material= couple->GetMaterial();
532
533 const G4ElementVector* theElementVector = material->GetElementVector();
534 size_t NumberOfElements = material->GetNumberOfElements() ;
535 const G4double* theAtomicNumDensityVector =
536 material->GetAtomicNumDensityVector();
537 G4VDataSetAlgorithm* interp = new G4SemiLogInterpolation();
538 G4VEMDataSet* xsis = new G4CompositeEMDataSet(interp, 1., 1.);
539 G4VDataSetAlgorithm* interp1 = new G4SemiLogInterpolation();
540 G4VEMDataSet* xsis1 = new G4CompositeEMDataSet(interp1, 1., 1.);
541
542 G4double tCut = cutForDelta[j];
543 G4double elDensity = 1.;
544
545 for (size_t iel=0; iel<NumberOfElements; iel++ ) {
546
547 G4int Z = (G4int)((*theElementVector)[iel]->GetZ());
548 G4int nShells = transitionManager->NumberOfShells(Z);
549 energy = new G4DataVector();
550 ksi = new G4DataVector();
551 energy1= new G4DataVector();
552 ksi1 = new G4DataVector();
553 //if(NumberOfElements > 1)
554 elDensity = theAtomicNumDensityVector[iel]/((G4double)nShells);
555
556 for (size_t j = 0; j<binForFluo; j++) {
557
558 G4double tkin = bVector->GetLowEdgeEnergy(j);
559 G4double gamma = tkin/mass + 1.;
560 G4double beta2 = 1.0 - 1.0/(gamma*gamma);
561 G4double r = electron_mass_c2/mass;
562 G4double tmax = 2.*electron_mass_c2*(gamma*gamma - 1.)/(1. + 2.*gamma*r + r*r);
563 G4double cross = 0.;
564 G4double cross1 = 0.;
565 G4double eAverage= 0.;
566 G4double tmin = std::min(tCut,tmax);
567 G4double rel;
568
569 for (G4int n=0; n<nShells; n++) {
570
571 bindingEnergy = transitionManager->Shell(Z, n)->BindingEnergy();
572 if (tmin > bindingEnergy) {
573 rel = std::log(tmin/bindingEnergy);
574 eAverage += rel - beta2*(tmin - bindingEnergy)/tmax;
575 cross += 1.0/bindingEnergy - 1.0/tmin - beta2*rel/tmax;
576 }
577 if (tmax > tmin) {
578 cross1 += 1.0/tmin - 1.0/tmax - beta2*std::log(tmax/tmin)/tmax;
579 }
580 }
581
582 cross1 *= elDensity;
583 energy1->push_back(tkin);
584 ksi1->push_back(cross1);
585
586 if(eAverage > 0.) cross /= eAverage;
587 else cross = 0.;
588
589 energy->push_back(tkin);
590 ksi->push_back(cross);
591 }
592 G4VDataSetAlgorithm* algo = interp->Clone();
593 G4VEMDataSet* set = new G4EMDataSet(Z,energy,ksi,algo,1.,1.);
594 xsis->AddComponent(set);
595 G4VDataSetAlgorithm* algo1 = interp1->Clone();
596 G4VEMDataSet* set1 = new G4EMDataSet(Z,energy1,ksi1,algo1,1.,1.);
597 xsis1->AddComponent(set1);
598 }
599 if(verboseLevel > 1) {
600 G4cout << "### Shell inverse cross sections for "
601 << material->GetName() << G4endl;
602 xsis->PrintData();
603 G4cout << "### Atom cross sections for "
604 << material->GetName() << G4endl;
605 xsis1->PrintData();
606 }
607 shellVacancy->AddXsiTable(xsis);
608 zFluoDataVector.push_back(xsis1);
609 }
610 delete bVector;
611}
612
613//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
614
615void G4hLowEnergyIonisation::BuildLambdaTable(
616 const G4ParticleDefinition& aParticleType)
617
618{
619 // Build mean free path tables for the delta ray production process
620 // tables are built for MATERIALS
621
622 if(verboseLevel > 1) {
623 G4cout << "G4hLowEnergyIonisation::BuildLambdaTable for "
624 << aParticleType.GetParticleName() << " is started" << G4endl;
625 }
626
627
628 G4double lowEdgeEnergy, value;
629 charge = aParticleType.GetPDGCharge()/eplus ;
630 chargeSquare = charge*charge ;
631 initialMass = aParticleType.GetPDGMass();
632
633 const G4ProductionCutsTable* theCoupleTable=
634 G4ProductionCutsTable::GetProductionCutsTable();
635 size_t numOfCouples = theCoupleTable->GetTableSize();
636
637
638 if (theMeanFreePathTable) {
639 theMeanFreePathTable->clearAndDestroy();
640 delete theMeanFreePathTable;
641 }
642
643 theMeanFreePathTable = new G4PhysicsTable(numOfCouples);
644
645 // loop for materials
646
647 for (size_t J=0 ; J < numOfCouples; J++) {
648
649 //create physics vector then fill it ....
650 G4PhysicsLogVector* aVector = new G4PhysicsLogVector(LowestKineticEnergy,
651 HighestKineticEnergy,
652 TotBin);
653
654 // compute the (macroscopic) cross section first
655 const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(J);
656 const G4Material* material= couple->GetMaterial();
657
658 const G4ElementVector* theElementVector =
659 material->GetElementVector() ;
660 const G4double* theAtomicNumDensityVector =
661 material->GetAtomicNumDensityVector();
662 const G4int NumberOfElements = material->GetNumberOfElements() ;
663
664 // get the electron kinetic energy cut for the actual material,
665 // it will be used in ComputeMicroscopicCrossSection
666 // ( it is the SAME for ALL the ELEMENTS in THIS MATERIAL )
667 // ------------------------------------------------------
668
669 G4double deltaCut = cutForDelta[J];
670
671 for ( G4int i = 0 ; i < TotBin ; i++ ) {
672 lowEdgeEnergy = aVector->GetLowEdgeEnergy(i) ;
673 G4double sigma = 0.0 ;
674 G4int Z;
675
676 for (G4int iel=0; iel<NumberOfElements; iel++ ) {
677 Z = (G4int) (*theElementVector)[iel]->GetZ();
678 totalCrossSectionMap [Z] = ComputeMicroscopicCrossSection(
679 aParticleType,
680 lowEdgeEnergy,
681 Z,
682 deltaCut ) ;
683 sigma += theAtomicNumDensityVector[iel]*ComputeMicroscopicCrossSection(
684 aParticleType,
685 lowEdgeEnergy,
686 Z,
687 deltaCut ) ;
688
689 }
690
691 // mean free path = 1./macroscopic cross section
692
693 value = sigma<=0 ? DBL_MAX : 1./sigma ;
694
695 aVector->PutValue(i, value) ;
696 }
697
698 theMeanFreePathTable->insert(aVector);
699 }
700
701}
702
703
704//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
705
706G4double G4hLowEnergyIonisation::ComputeMicroscopicCrossSection(
707 const G4ParticleDefinition& aParticleType,
708 G4double kineticEnergy,
709 G4double atomicNumber,
710 G4double deltaCutInEnergy) const
711{
712 //******************************************************************
713 // cross section formula is OK for spin=0, 1/2, 1 only !
714 // *****************************************************************
715
716 // calculates the microscopic cross section in GEANT4 internal units
717 // ( it is called for elements , AtomicNumber = z )
718
719 G4double energy, gamma, beta2, tmax, var;
720 G4double totalCrossSection = 0.0 ;
721
722 G4double particleMass = initialMass;
723
724 // get particle data ...................................
725
726 energy = kineticEnergy + particleMass;
727
728 // some kinematics......................
729
730 gamma = energy/particleMass;
731 beta2 = 1.0 - 1.0/(gamma*gamma);
732 var = electron_mass_c2/particleMass;
733 tmax = 2.*electron_mass_c2*(gamma*gamma - 1.)/(1. + 2.*gamma*var + var*var);
734
735 // now you can calculate the total cross section
736
737 if( tmax > deltaCutInEnergy ) {
738
739 var=deltaCutInEnergy/tmax;
740 totalCrossSection = (1.0 - var*(1.0 - beta2*std::log(var))) / deltaCutInEnergy ;
741 G4double spin = aParticleType.GetPDGSpin() ;
742
743 // +term for spin=1/2 particle
744 if( 0.5 == spin )
745 totalCrossSection += 0.5 * (tmax - deltaCutInEnergy) / (energy*energy);
746
747 // +term for spin=1 particle
748 else if( 0.9 < spin )
749 totalCrossSection += -std::log(var)/(3.0*deltaCutInEnergy) +
750 (tmax - deltaCutInEnergy) * ( (5.0+ 1.0/var)*0.25 / (energy*energy) -
751 beta2 / (tmax * deltaCutInEnergy) ) / 3.0 ;
752
753 totalCrossSection *= twopi_mc2_rcl2 * atomicNumber / beta2 ;
754 }
755
756 return totalCrossSection ;
757}
758
759//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
760
761G4double G4hLowEnergyIonisation::GetMeanFreePath(const G4Track& trackData,
762 G4double, // previousStepSize
763 enum G4ForceCondition* condition)
764{
765 const G4DynamicParticle* aParticle = trackData.GetDynamicParticle();
766 const G4MaterialCutsCouple* couple = trackData.GetMaterialCutsCouple();
767 const G4Material* material = couple->GetMaterial();
768 G4double meanFreePath;
769 G4bool isOutRange ;
770
771 *condition = NotForced ;
772
773 G4double kineticEnergy = (aParticle->GetKineticEnergy())*initialMass/(aParticle->GetMass());
774 charge = aParticle->GetCharge()/eplus;
775 chargeSquare = theIonEffChargeModel->TheValue(aParticle, material);
776
777 if(kineticEnergy < LowestKineticEnergy) meanFreePath = DBL_MAX;
778
779 else {
780 if(kineticEnergy > HighestKineticEnergy)
781 kineticEnergy = HighestKineticEnergy;
782 meanFreePath = (((*theMeanFreePathTable)(couple->GetIndex()))->
783 GetValue(kineticEnergy,isOutRange))/chargeSquare;
784 }
785
786 return meanFreePath ;
787}
788
789//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
790
791G4double G4hLowEnergyIonisation::GetConstraints(
792 const G4DynamicParticle* particle,
793 const G4MaterialCutsCouple* couple)
794{
795 // returns the Step limit
796 // dEdx is calculated as well as the range
797 // based on Effective Charge Approach
798
799 const G4Material* material = couple->GetMaterial();
800 G4Proton* theProton = G4Proton::Proton();
801 G4AntiProton* theAntiProton = G4AntiProton::AntiProton();
802
803 G4double stepLimit = 0.0 ;
804 G4double dx, highEnergy;
805
806 G4double massRatio = proton_mass_c2/(particle->GetMass()) ;
807 G4double kineticEnergy = particle->GetKineticEnergy() ;
808
809 // Scale the kinetic energy
810
811 G4double tscaled = kineticEnergy*massRatio ;
812 fBarkas = 0.0;
813
814 if(charge > 0.0) {
815
816 highEnergy = protonHighEnergy ;
817
818 fRangeNow = G4EnergyLossTables::GetRange(theProton, tscaled, couple);
819 dx = G4EnergyLossTables::GetRange(theProton, highEnergy, couple);
820 fdEdx = G4EnergyLossTables::GetDEDX(theProton, tscaled, couple)
821 * chargeSquare ;
822
823 // Correction for positive ions
824 if(theBarkas && tscaled > highEnergy) {
825 fBarkas = BarkasTerm(material,tscaled)*std::sqrt(chargeSquare)*chargeSquare
826 + BlochTerm(material,tscaled,chargeSquare);
827 }
828 // Antiprotons and negative hadrons
829 } else {
830
831 highEnergy = antiProtonHighEnergy ;
832 fRangeNow = G4EnergyLossTables::GetRange(theAntiProton, tscaled, couple);
833 dx = G4EnergyLossTables::GetRange(theAntiProton, highEnergy, couple);
834 fdEdx = G4EnergyLossTables::GetDEDX(theAntiProton, tscaled, couple)
835 * chargeSquare ;
836
837 if(theBarkas && tscaled > highEnergy) {
838 fBarkas = -BarkasTerm(material,tscaled)*std::sqrt(chargeSquare)*chargeSquare
839 + BlochTerm(material,tscaled,chargeSquare);
840 }
841 }
842 /*
843 const G4Material* mat = couple->GetMaterial();
844 G4double fac = gram/(MeV*cm2*mat->GetDensity());
845 G4cout << particle->GetDefinition()->GetParticleName()
846 << " in " << mat->GetName()
847 << " E(MeV)= " << kineticEnergy/MeV
848 << " dedx(MeV*cm^2/g)= " << fdEdx*fac
849 << " barcas(MeV*cm^2/gram)= " << fBarkas*fac
850 << " Q^2= " << chargeSquare
851 << G4endl;
852 */
853 // scaling back
854 fRangeNow /= (chargeSquare*massRatio) ;
855 dx /= (chargeSquare*massRatio) ;
856
857 stepLimit = fRangeNow ;
858 G4double r = std::min(finalRange, couple->GetProductionCuts()
859 ->GetProductionCut(idxG4ElectronCut));
860
861 if (fRangeNow > r) {
862 stepLimit = dRoverRange*fRangeNow + r*(1.0 - dRoverRange)*(2.0 - r/fRangeNow);
863 if (stepLimit > fRangeNow) stepLimit = fRangeNow;
864 }
865 // compute the (random) Step limit in standard energy range
866 if(tscaled > highEnergy ) {
867
868 // add Barkas correction directly to dedx
869 fdEdx += fBarkas;
870
871 if(stepLimit > fRangeNow - dx*0.9) stepLimit = fRangeNow - dx*0.9 ;
872
873 // Step limit in low energy range
874 } else {
875 G4double x = dx*paramStepLimit;
876 if (stepLimit > x) stepLimit = x;
877 }
878 return stepLimit ;
879}
880
881//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
882
883G4VParticleChange* G4hLowEnergyIonisation::AlongStepDoIt(
884 const G4Track& trackData,
885 const G4Step& stepData)
886{
887 // compute the energy loss after a step
888 G4Proton* theProton = G4Proton::Proton();
889 G4AntiProton* theAntiProton = G4AntiProton::AntiProton();
890 G4double finalT = 0.0 ;
891
892 aParticleChange.Initialize(trackData) ;
893
894 const G4MaterialCutsCouple* couple = trackData.GetMaterialCutsCouple();
895 const G4Material* material = couple->GetMaterial();
896
897 // get the actual (true) Step length from stepData
898 const G4double step = stepData.GetStepLength() ;
899
900 const G4DynamicParticle* particle = trackData.GetDynamicParticle() ;
901
902 G4double kineticEnergy = particle->GetKineticEnergy() ;
903 G4double massRatio = proton_mass_c2/(particle->GetMass()) ;
904 G4double tscaled= kineticEnergy*massRatio ;
905 G4double eloss = 0.0 ;
906 G4double nloss = 0.0 ;
907
908
909 // very small particle energy
910 if(kineticEnergy < MinKineticEnergy) {
911
912 eloss = kineticEnergy ;
913
914 // particle energy outside tabulated energy range
915 } else if( kineticEnergy > HighestKineticEnergy) {
916 eloss = step*fdEdx ;
917
918 // big step
919 } else if(step >= fRangeNow ) {
920 eloss = kineticEnergy ;
921
922 // tabulated range
923 } else {
924
925 // step longer than linear step limit
926 if(step > linLossLimit*fRangeNow) {
927
928 G4double rscaled= fRangeNow*massRatio*chargeSquare ;
929 G4double sscaled= step *massRatio*chargeSquare ;
930
931 if(charge > 0.0) {
932 eloss = G4EnergyLossTables::GetPreciseEnergyFromRange(
933 theProton,rscaled, couple) -
934 G4EnergyLossTables::GetPreciseEnergyFromRange(
935 theProton,rscaled-sscaled,couple) ;
936
937 } else {
938 eloss = G4EnergyLossTables::GetPreciseEnergyFromRange(
939 theAntiProton,rscaled,couple) -
940 G4EnergyLossTables::GetPreciseEnergyFromRange(
941 theAntiProton,rscaled-sscaled,couple) ;
942 }
943 eloss /= massRatio ;
944
945 // Barkas correction at big step
946 eloss += fBarkas*step;
947
948 // step shorter than linear step limit
949 } else {
950 eloss = step*fdEdx ;
951 }
952 if(nStopping && tscaled < protonHighEnergy) {
953 nloss = (theNuclearStoppingModel->TheValue(particle, material))*step;
954 }
955 }
956
957 if(eloss < 0.0) eloss = 0.0;
958
959 finalT = kineticEnergy - eloss - nloss;
960
961 if( EnlossFlucFlag && 0.0 < eloss && finalT > MinKineticEnergy) {
962
963 // now the electron loss with fluctuation
964 eloss = ElectronicLossFluctuation(particle, couple, eloss, step) ;
965 if(eloss < 0.0) eloss = 0.0;
966 finalT = kineticEnergy - eloss - nloss;
967 }
968
969 // stop particle if the kinetic energy <= MinKineticEnergy
970 if (finalT*massRatio <= MinKineticEnergy ) {
971
972 finalT = 0.0;
973 if(!particle->GetDefinition()->GetProcessManager()->
974 GetAtRestProcessVector()->size())
975 aParticleChange.ProposeTrackStatus(fStopAndKill);
976 else
977 aParticleChange.ProposeTrackStatus(fStopButAlive);
978 }
979
980 aParticleChange.ProposeEnergy( finalT );
981 eloss = kineticEnergy-finalT;
982
983 // Deexcitation only of ionised atoms
984 G4double hMass = particle->GetMass();
985 std::vector<G4DynamicParticle*>* newpart = 0;
986 G4DynamicParticle* part = 0;
987
988 if(theFluo) newpart = DeexciteAtom(couple, kineticEnergy, hMass, eloss);
989
990 if(newpart != 0) {
991
992 size_t nSecondaries = newpart->size();
993 aParticleChange.SetNumberOfSecondaries(nSecondaries);
994 G4Track* newtrack = 0;
995 const G4StepPoint* preStep = stepData.GetPreStepPoint();
996 const G4StepPoint* postStep = stepData.GetPostStepPoint();
997 G4ThreeVector r = preStep->GetPosition();
998 G4ThreeVector deltaR = postStep->GetPosition();
999 deltaR -= r;
1000 G4double t = preStep->GetGlobalTime();
1001 G4double deltaT = postStep->GetGlobalTime();
1002 deltaT -= t;
1003 G4double time, q, e;
1004 G4ThreeVector position;
1005
1006 for(size_t i=0; i<nSecondaries; i++) {
1007
1008 part = (*newpart)[i];
1009 if(part) {
1010
1011 e = part->GetKineticEnergy();
1012 if(e <= eloss) {
1013
1014 eloss -= e;
1015 q = G4UniformRand();
1016 time = deltaT*q + t;
1017 position = deltaR*q;
1018 position += r;
1019 newtrack = new G4Track(part, time, position);
1020 aParticleChange.AddSecondary(newtrack);
1021
1022 } else {
1023
1024 delete part;
1025
1026 }
1027 }
1028 }
1029 delete newpart;
1030 }
1031
1032 aParticleChange.ProposeLocalEnergyDeposit(eloss);
1033 return &aParticleChange ;
1034}
1035
1036//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1037
1038G4double G4hLowEnergyIonisation::ProtonParametrisedDEDX(
1039 const G4MaterialCutsCouple* couple,
1040 G4double kineticEnergy) const
1041{
1042 const G4Material* material = couple->GetMaterial();
1043 G4Proton* theProton = G4Proton::Proton();
1044 G4double eloss = 0.0;
1045
1046 // Free Electron Gas Model
1047 if(kineticEnergy < protonLowEnergy) {
1048 eloss = (theProtonModel->TheValue(theProton, material, protonLowEnergy))
1049 * std::sqrt(kineticEnergy/protonLowEnergy) ;
1050
1051 // Parametrisation
1052 } else {
1053 eloss = theProtonModel->TheValue(theProton, material, kineticEnergy) ;
1054 }
1055
1056 // Delta rays energy
1057 eloss -= DeltaRaysEnergy(couple,kineticEnergy,proton_mass_c2) ;
1058
1059 if(verboseLevel > 2) {
1060 G4cout << "p E(MeV)= " << kineticEnergy/MeV
1061 << " dE/dx(MeV/mm)= " << eloss*mm/MeV
1062 << " for " << material->GetName()
1063 << " model: " << theProtonModel << G4endl;
1064 }
1065
1066 if(eloss < 0.0) eloss = 0.0 ;
1067
1068 return eloss ;
1069}
1070
1071//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1072
1073G4double G4hLowEnergyIonisation::AntiProtonParametrisedDEDX(
1074 const G4MaterialCutsCouple* couple,
1075 G4double kineticEnergy) const
1076{
1077 const G4Material* material = couple->GetMaterial();
1078 G4AntiProton* theAntiProton = G4AntiProton::AntiProton();
1079 G4double eloss = 0.0 ;
1080
1081 // Antiproton model is used
1082 if(theAntiProtonModel->IsInCharge(theAntiProton,material)) {
1083 if(kineticEnergy < antiProtonLowEnergy) {
1084 eloss = theAntiProtonModel->TheValue(theAntiProton,material,antiProtonLowEnergy)
1085 * std::sqrt(kineticEnergy/antiProtonLowEnergy) ;
1086
1087 // Parametrisation
1088 } else {
1089 eloss = theAntiProtonModel->TheValue(theAntiProton,material,
1090 kineticEnergy);
1091 }
1092
1093 // The proton model is used + Barkas correction
1094 } else {
1095 if(kineticEnergy < protonLowEnergy) {
1096 eloss = theProtonModel->TheValue(G4Proton::Proton(),material,protonLowEnergy)
1097 * std::sqrt(kineticEnergy/protonLowEnergy) ;
1098
1099 // Parametrisation
1100 } else {
1101 eloss = theProtonModel->TheValue(G4Proton::Proton(),material,
1102 kineticEnergy);
1103 }
1104 //if(theBarkas) eloss -= 2.0*BarkasTerm(material, kineticEnergy);
1105 }
1106
1107 // Delta rays energy
1108 eloss -= DeltaRaysEnergy(couple,kineticEnergy,proton_mass_c2) ;
1109
1110 if(verboseLevel > 2) {
1111 G4cout << "pbar E(MeV)= " << kineticEnergy/MeV
1112 << " dE/dx(MeV/mm)= " << eloss*mm/MeV
1113 << " for " << material->GetName()
1114 << " model: " << theProtonModel << G4endl;
1115 }
1116
1117 if(eloss < 0.0) eloss = 0.0 ;
1118
1119 return eloss ;
1120}
1121
1122//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1123
1124G4double G4hLowEnergyIonisation::DeltaRaysEnergy(
1125 const G4MaterialCutsCouple* couple,
1126 G4double kineticEnergy,
1127 G4double particleMass) const
1128{
1129 G4double dloss = 0.0 ;
1130
1131 G4double deltaCutNow = cutForDelta[(couple->GetIndex())] ;
1132 const G4Material* material = couple->GetMaterial();
1133 G4double electronDensity = material->GetElectronDensity();
1134 G4double eexc = material->GetIonisation()->GetMeanExcitationEnergy();
1135
1136 G4double tau = kineticEnergy/particleMass ;
1137 G4double rateMass = electron_mass_c2/particleMass ;
1138
1139 // some local variables
1140
1141 G4double gamma,bg2,beta2,tmax,x ;
1142
1143 gamma = tau + 1.0 ;
1144 bg2 = tau*(tau+2.0) ;
1145 beta2 = bg2/(gamma*gamma) ;
1146 tmax = 2.*electron_mass_c2*bg2/(1.0+2.0*gamma*rateMass+rateMass*rateMass) ;
1147
1148 // Validity range for delta electron cross section
1149 G4double deltaCut = std::max(deltaCutNow, eexc);
1150
1151 if ( deltaCut < tmax) {
1152 x = deltaCut / tmax ;
1153 dloss = ( beta2 * (x - 1.0) - std::log(x) ) * twopi_mc2_rcl2
1154 * electronDensity / beta2 ;
1155 }
1156 return dloss ;
1157}
1158
1159//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1160
1161G4VParticleChange* G4hLowEnergyIonisation::PostStepDoIt(
1162 const G4Track& trackData,
1163 const G4Step& stepData)
1164{
1165 // Units are expressed in GEANT4 internal units.
1166
1167 G4double KineticEnergy,TotalEnergy,TotalMomentum,betasquare,
1168 DeltaKineticEnergy,DeltaTotalMomentum,costheta,sintheta,phi,
1169 dirx,diry,dirz,finalKineticEnergy,finalPx,finalPy,finalPz,
1170 x,xc,grej,Psquare,Esquare,rate,finalMomentum ;
1171
1172 aParticleChange.Initialize(trackData) ;
1173 const G4MaterialCutsCouple* couple = trackData.GetMaterialCutsCouple();
1174
1175 const G4DynamicParticle* aParticle = trackData.GetDynamicParticle() ;
1176
1177 // some kinematics
1178
1179 ParticleMass=aParticle->GetDefinition()->GetPDGMass();
1180 KineticEnergy=aParticle->GetKineticEnergy();
1181 TotalEnergy=KineticEnergy + ParticleMass ;
1182 Psquare=KineticEnergy*(TotalEnergy+ParticleMass) ;
1183 Esquare=TotalEnergy*TotalEnergy;
1184 betasquare=Psquare/Esquare;
1185 G4ThreeVector ParticleDirection = aParticle->GetMomentumDirection() ;
1186
1187 G4double gamma= KineticEnergy/ParticleMass + 1.;
1188 G4double r = electron_mass_c2/ParticleMass;
1189 G4double tmax = 2.*electron_mass_c2*(gamma*gamma - 1.)/(1. + 2.*gamma*r + r*r);
1190
1191 // Validity range for delta electron cross section
1192 G4double DeltaCut = cutForDelta[couple->GetIndex()];
1193
1194 // This should not be a case
1195 if(DeltaCut >= tmax)
1196 return G4VContinuousDiscreteProcess::PostStepDoIt(trackData,stepData);
1197
1198 xc = DeltaCut / tmax;
1199 rate = tmax / TotalEnergy;
1200 rate = rate*rate ;
1201 G4double spin = aParticle->GetDefinition()->GetPDGSpin() ;
1202
1203 // sampling follows ...
1204 do {
1205 x=xc/(1.-(1.-xc)*G4UniformRand());
1206
1207 if(0.0 == spin) {
1208 grej = 1.0 - betasquare * x ;
1209
1210 } else if (0.5 == spin) {
1211 grej = (1.0 - betasquare * x + 0.5*x*x*rate) / (1.0 + 0.5 * rate) ;
1212
1213 } else {
1214 grej = (1.0 - betasquare * x ) * (1.0 + x/ (3.0*xc)) +
1215 x * x * rate * (1.0 + 0.5 * x / xc) / 3.0 /
1216 (1.0 + 1.0/(3.0*xc) + rate *(1.0+ 0.5/xc) /3.0) ;
1217 }
1218
1219 } while( G4UniformRand() > grej );
1220
1221
1222 DeltaKineticEnergy = x * tmax;
1223
1224 DeltaTotalMomentum = std::sqrt(DeltaKineticEnergy * (DeltaKineticEnergy +
1225 2. * electron_mass_c2 )) ;
1226 TotalMomentum = std::sqrt(Psquare) ;
1227 costheta = DeltaKineticEnergy * (TotalEnergy + electron_mass_c2)
1228 /(DeltaTotalMomentum * TotalMomentum) ;
1229
1230 // protection against costheta > 1 or < -1 ---------------
1231 if ( costheta < -1. )
1232 costheta = -1. ;
1233 if ( costheta > +1. )
1234 costheta = +1. ;
1235
1236 // direction of the delta electron ........
1237 phi = twopi * G4UniformRand() ;
1238 sintheta = std::sqrt(1. - costheta*costheta);
1239 dirx = sintheta * std::cos(phi) ;
1240 diry = sintheta * std::sin(phi) ;
1241 dirz = costheta ;
1242
1243 G4ThreeVector DeltaDirection(dirx,diry,dirz) ;
1244 DeltaDirection.rotateUz(ParticleDirection) ;
1245
1246 // create G4DynamicParticle object for delta ray
1247 G4DynamicParticle *theDeltaRay = new G4DynamicParticle;
1248 theDeltaRay->SetKineticEnergy( DeltaKineticEnergy );
1249 theDeltaRay->SetMomentumDirection(DeltaDirection.x(),
1250 DeltaDirection.y(),
1251 DeltaDirection.z());
1252 theDeltaRay->SetDefinition(G4Electron::Electron());
1253
1254 // fill aParticleChange
1255 finalKineticEnergy = KineticEnergy - DeltaKineticEnergy ;
1256
1257 // Generation of Fluorescence and Auger
1258 size_t nSecondaries = 0;
1259 size_t totalNumber = 1;
1260 std::vector<G4DynamicParticle*>* secondaryVector = 0;
1261 G4DynamicParticle* aSecondary = 0;
1262 G4ParticleDefinition* type = 0;
1263
1264 // Select atom and shell
1265 G4int Z = SelectRandomAtom(couple, KineticEnergy);
1266
1267 // G4cout << "Fluorescence is switched :" << theFluo << G4endl;
1268
1269 // Fluorescence data start from element 6
1270 if(theFluo && Z > 5) {
1271
1272
1273
1274 // Atom total cross section
1275 shellCS->SetTotalCS(totalCrossSectionMap[Z]);
1276
1277 G4int shell = shellCS->SelectRandomShell(Z, KineticEnergy,ParticleMass,DeltaKineticEnergy);
1278
1279 if (shell!=-1) {
1280
1281 const G4AtomicShell* atomicShell =
1282 (G4AtomicTransitionManager::Instance())->Shell(Z, shell);
1283 G4double bindingEnergy = atomicShell->BindingEnergy();
1284
1285 if(verboseLevel > 1) {
1286 G4cout << "PostStep Z= " << Z << " shell= " << shell
1287 << " bindingE(keV)= " << bindingEnergy/keV
1288 << " finalE(keV)= " << finalKineticEnergy/keV
1289 << G4endl;
1290 }
1291
1292
1293
1294 if (finalKineticEnergy >= bindingEnergy
1295 && (bindingEnergy >= minGammaEnergy
1296 || bindingEnergy >= minElectronEnergy) ) {
1297
1298 G4int shellId = atomicShell->ShellId();
1299 secondaryVector = deexcitationManager.GenerateParticles(Z, shellId);
1300
1301 if (secondaryVector != 0) {
1302
1303 nSecondaries = secondaryVector->size();
1304 for (size_t i = 0; i<nSecondaries; i++) {
1305
1306 aSecondary = (*secondaryVector)[i];
1307 if (aSecondary) {
1308
1309 G4double e = aSecondary->GetKineticEnergy();
1310 type = aSecondary->GetDefinition();
1311 if (e < finalKineticEnergy &&
1312 ((type == G4Gamma::Gamma() && e > minGammaEnergy ) ||
1313 (type == G4Electron::Electron() && e > minElectronEnergy ))) {
1314
1315 finalKineticEnergy -= e;
1316 totalNumber++;
1317
1318 } else {
1319
1320 delete aSecondary;
1321 (*secondaryVector)[i] = 0;
1322 }
1323 }
1324 }
1325 }
1326 }
1327 }
1328 }
1329
1330 // Save delta-electrons
1331
1332 G4double edep = 0.0;
1333
1334 if (finalKineticEnergy > MinKineticEnergy)
1335 {
1336 finalPx = TotalMomentum*ParticleDirection.x()
1337 - DeltaTotalMomentum*DeltaDirection.x();
1338 finalPy = TotalMomentum*ParticleDirection.y()
1339 - DeltaTotalMomentum*DeltaDirection.y();
1340 finalPz = TotalMomentum*ParticleDirection.z()
1341 - DeltaTotalMomentum*DeltaDirection.z();
1342 finalMomentum =
1343 std::sqrt(finalPx*finalPx+finalPy*finalPy+finalPz*finalPz) ;
1344 finalPx /= finalMomentum ;
1345 finalPy /= finalMomentum ;
1346 finalPz /= finalMomentum ;
1347
1348 aParticleChange.ProposeMomentumDirection( finalPx,finalPy,finalPz );
1349 }
1350 else
1351 {
1352 edep = finalKineticEnergy;
1353 finalKineticEnergy = 0.;
1354 aParticleChange.ProposeMomentumDirection(ParticleDirection.x(),
1355 ParticleDirection.y(),ParticleDirection.z());
1356 if(!aParticle->GetDefinition()->GetProcessManager()->
1357 GetAtRestProcessVector()->size())
1358 aParticleChange.ProposeTrackStatus(fStopAndKill);
1359 else
1360 aParticleChange.ProposeTrackStatus(fStopButAlive);
1361 }
1362
1363 aParticleChange.ProposeEnergy( finalKineticEnergy );
1364 aParticleChange.ProposeLocalEnergyDeposit (edep);
1365 aParticleChange.SetNumberOfSecondaries(totalNumber);
1366 aParticleChange.AddSecondary(theDeltaRay);
1367
1368 // Save Fluorescence and Auger
1369
1370 if (secondaryVector) {
1371
1372 for (size_t l = 0; l < nSecondaries; l++) {
1373
1374 aSecondary = (*secondaryVector)[l];
1375 if(aSecondary) {
1376 aParticleChange.AddSecondary(aSecondary);
1377 }
1378 }
1379 delete secondaryVector;
1380 }
1381
1382 return G4VContinuousDiscreteProcess::PostStepDoIt(trackData,stepData);
1383}
1384
1385//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1386
1387
1388
1389void G4hLowEnergyIonisation::SelectShellIonisationCS(G4String val) {
1390
1391 if (val == "analytical" ) {
1392 if (shellCS) delete shellCS;
1393 shellCS = new G4teoCrossSection(val);
1394 }
1395 else if (val == "empirical") {
1396 if (shellCS) delete shellCS;
1397 shellCS = new G4empCrossSection();
1398 }
1399}
1400
1401
1402
1403//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1404
1405
1406std::vector<G4DynamicParticle*>* G4hLowEnergyIonisation::DeexciteAtom(const G4MaterialCutsCouple* couple,
1407 G4double incidentEnergy,
1408 G4double hMass,
1409 G4double eLoss)
1410{
1411
1412 if (verboseLevel > 1) {
1413 G4cout << "DeexciteAtom: cutForPhotons(keV)= " << minGammaEnergy/keV
1414 << " cutForElectrons(keV)= " << minElectronEnergy/keV
1415 << " eLoss(MeV)= " << eLoss
1416 << G4endl;
1417 }
1418
1419
1420
1421 if(eLoss < minGammaEnergy && eLoss < minElectronEnergy) return 0;
1422
1423 const G4Material* material = couple->GetMaterial();
1424 G4int index = couple->GetIndex();
1425 // G4double eexc = material->GetIonisation()->GetMeanExcitationEnergy();
1426 G4double gamma = incidentEnergy/hMass + 1;
1427 G4double beta2 = 1.0 - 1.0/(gamma*gamma);
1428 G4double r = electron_mass_c2/hMass;
1429 G4double tmax = 2.*electron_mass_c2*(gamma*gamma - 1.)/(1. + 2.*gamma*r + r*r);
1430 G4double tcut = std::min(tmax,cutForDelta[index]);
1431 const G4AtomicTransitionManager* transitionManager =
1432 G4AtomicTransitionManager::Instance();
1433
1434 size_t nElements = material->GetNumberOfElements();
1435 const G4ElementVector* theElementVector = material->GetElementVector();
1436 G4bool stop = true;
1437
1438 for (size_t j=0; j<nElements; j++) {
1439
1440 G4int Z = (G4int)((*theElementVector)[j]->GetZ());
1441 G4double maxE = transitionManager->Shell(Z, 0)->BindingEnergy();
1442
1443 if (Z > 5 && maxE < tcut && (maxE > minGammaEnergy || maxE > minElectronEnergy) ) {
1444 stop = false;
1445 break;
1446 }
1447 }
1448
1449 if(stop) return 0;
1450
1451 // create vector of tracks of secondary particles
1452
1453 std::vector<G4DynamicParticle*>* partVector =
1454 new std::vector<G4DynamicParticle*>;
1455 std::vector<G4DynamicParticle*>* secVector = 0;
1456 G4DynamicParticle* aSecondary = 0;
1457 G4ParticleDefinition* type = 0;
1458 G4double e, tkin, grej;
1459 G4ThreeVector position;
1460 G4int shell, shellId;
1461
1462 // sample secondaries
1463
1464 G4double etot = 0.0;
1465 std::vector<G4int> n = shellVacancy->GenerateNumberOfIonisations(couple,
1466 incidentEnergy, eLoss);
1467
1468 for (size_t i=0; i<nElements; i++) {
1469
1470 size_t nVacancies = n[i];
1471 G4int Z = (G4int)((*theElementVector)[i]->GetZ());
1472 G4double maxE = transitionManager->Shell(Z, 0)->BindingEnergy();
1473
1474 if (nVacancies && Z > 5 && maxE < tcut && (maxE > minGammaEnergy || maxE > minElectronEnergy)) {
1475 for(size_t j=0; j<nVacancies; j++) {
1476
1477 // sampling follows
1478 do {
1479 tkin = tcut/(1.0 + (tcut/maxE - 1.0)*G4UniformRand());
1480 grej = 1.0 - beta2 * tkin/tmax;
1481
1482 } while( G4UniformRand() > grej );
1483
1484 shell = shellCS->SelectRandomShell(Z,incidentEnergy,hMass,tkin);
1485
1486 shellId = transitionManager->Shell(Z, shell)->ShellId();
1487 G4double maxE = transitionManager->Shell(Z, shell)->BindingEnergy();
1488
1489 if (maxE>minGammaEnergy || maxE>minElectronEnergy ) {
1490 secVector = deexcitationManager.GenerateParticles(Z, shellId);
1491 } else {
1492 secVector = 0;
1493 }
1494
1495 if (secVector) {
1496
1497 for (size_t l = 0; l<secVector->size(); l++) {
1498
1499 aSecondary = (*secVector)[l];
1500 if(aSecondary) {
1501
1502 e = aSecondary->GetKineticEnergy();
1503 type = aSecondary->GetDefinition();
1504 if ( etot + e <= eLoss &&
1505 ( (type == G4Gamma::Gamma() && e > minGammaEnergy ) ||
1506 (type == G4Electron::Electron() && e > minElectronEnergy) ) ) {
1507
1508 etot += e;
1509 partVector->push_back(aSecondary);
1510
1511 } else {
1512 delete aSecondary;
1513 }
1514 }
1515 }
1516 delete secVector;
1517 }
1518 }
1519 }
1520 }
1521
1522 if(partVector->empty()) {
1523 delete partVector;
1524 return 0;
1525 }
1526
1527 return partVector;
1528}
1529
1530//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1531
1532G4int G4hLowEnergyIonisation::SelectRandomAtom(const G4MaterialCutsCouple* couple,
1533 G4double kineticEnergy) const
1534{
1535 const G4Material* material = couple->GetMaterial();
1536 G4int nElements = material->GetNumberOfElements();
1537 G4int Z = 0;
1538
1539 if(nElements == 1) {
1540 Z = (G4int)(material->GetZ());
1541 return Z;
1542 }
1543
1544 const G4ElementVector* theElementVector = material->GetElementVector();
1545 std::vector<G4double> p;
1546 G4int index = couple->GetIndex();
1547
1548 G4double norm = 0.0;
1549 for (G4int j=0; j<nElements; j++) {
1550
1551 const G4VEMDataSet* set = (zFluoDataVector[index])->GetComponent(j);
1552 G4double cross = set->FindValue(kineticEnergy);
1553
1554 p.push_back(cross);
1555 norm += cross;
1556 }
1557
1558 if(norm == 0.0) return 0;
1559
1560 G4double q = norm*G4UniformRand();
1561
1562 for (G4int i=0; i<nElements; i++) {
1563
1564 if(p[i] > q) {
1565 Z = (G4int)((*theElementVector)[i]->GetZ());
1566 break;
1567 }
1568 q -= p[i];
1569 }
1570
1571 return Z;
1572}
1573
1574//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1575
1576G4double G4hLowEnergyIonisation::ComputeDEDX(
1577 const G4ParticleDefinition* aParticle,
1578 const G4MaterialCutsCouple* couple,
1579 G4double kineticEnergy)
1580{
1581 const G4Material* material = couple->GetMaterial();
1582 G4Proton* theProton = G4Proton::Proton();
1583 G4AntiProton* theAntiProton = G4AntiProton::AntiProton();
1584 G4double dedx = 0.0 ;
1585
1586 G4double tscaled = kineticEnergy*proton_mass_c2/(aParticle->GetPDGMass()) ;
1587 charge = aParticle->GetPDGCharge() ;
1588
1589 if(charge>0.0) {
1590 if(tscaled > protonHighEnergy) {
1591 dedx=G4EnergyLossTables::GetDEDX(theProton,tscaled,couple) ;
1592
1593 } else {
1594 dedx=ProtonParametrisedDEDX(couple,tscaled) ;
1595 }
1596
1597 } else {
1598 if(tscaled > antiProtonHighEnergy) {
1599 dedx=G4EnergyLossTables::GetDEDX(theAntiProton,tscaled,couple);
1600
1601 } else {
1602 dedx=AntiProtonParametrisedDEDX(couple,tscaled) ;
1603 }
1604 }
1605 dedx *= theIonEffChargeModel->TheValue(aParticle, material, kineticEnergy) ;
1606
1607 return dedx ;
1608}
1609
1610//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1611
1612G4double G4hLowEnergyIonisation::BarkasTerm(const G4Material* material,
1613 G4double kineticEnergy) const
1614//Function to compute the Barkas term for protons:
1615//
1616//Ref. Z_1^3 effect in the stopping power of matter for charged particles
1617// J.C Ashley and R.H.Ritchie
1618// Physical review B Vol.5 No.7 1 April 1972 pagg. 2393-2397
1619//
1620{
1621 static double FTable[47][2] = {
1622 { 0.02, 21.5},
1623 { 0.03, 20.0},
1624 { 0.04, 18.0},
1625 { 0.05, 15.6},
1626 { 0.06, 15.0},
1627 { 0.07, 14.0},
1628 { 0.08, 13.5},
1629 { 0.09, 13.},
1630 { 0.1, 12.2},
1631 { 0.2, 9.25},
1632 { 0.3, 7.0},
1633 { 0.4, 6.0},
1634 { 0.5, 4.5},
1635 { 0.6, 3.5},
1636 { 0.7, 3.0},
1637 { 0.8, 2.5},
1638 { 0.9, 2.0},
1639 { 1.0, 1.7},
1640 { 1.2, 1.2},
1641 { 1.3, 1.0},
1642 { 1.4, 0.86},
1643 { 1.5, 0.7},
1644 { 1.6, 0.61},
1645 { 1.7, 0.52},
1646 { 1.8, 0.5},
1647 { 1.9, 0.43},
1648 { 2.0, 0.42},
1649 { 2.1, 0.3},
1650 { 2.4, 0.2},
1651 { 3.0, 0.13},
1652 { 3.08, 0.1},
1653 { 3.1, 0.09},
1654 { 3.3, 0.08},
1655 { 3.5, 0.07},
1656 { 3.8, 0.06},
1657 { 4.0, 0.051},
1658 { 4.1, 0.04},
1659 { 4.8, 0.03},
1660 { 5.0, 0.024},
1661 { 5.1, 0.02},
1662 { 6.0, 0.013},
1663 { 6.5, 0.01},
1664 { 7.0, 0.009},
1665 { 7.1, 0.008},
1666 { 8.0, 0.006},
1667 { 9.0, 0.0032},
1668 { 10.0, 0.0025} };
1669
1670 // Information on particle and material
1671 G4double kinE = kineticEnergy ;
1672 if(0.5*MeV > kinE) kinE = 0.5*MeV ;
1673 G4double gamma = 1.0 + kinE / proton_mass_c2 ;
1674 G4double beta2 = 1.0 - 1.0/(gamma*gamma) ;
1675 if(0.0 >= beta2) return 0.0;
1676
1677 G4double BarkasTerm = 0.0;
1678 G4double AMaterial = 0.0;
1679 G4double ZMaterial = 0.0;
1680 const G4ElementVector* theElementVector = material->GetElementVector();
1681 G4int numberOfElements = material->GetNumberOfElements();
1682
1683 for (G4int i = 0; i<numberOfElements; i++) {
1684
1685 AMaterial = (*theElementVector)[i]->GetA()*mole/g;
1686 ZMaterial = (*theElementVector)[i]->GetZ();
1687
1688 G4double X = 137.0 * 137.0 * beta2 / ZMaterial;
1689
1690 // Variables to compute L_1
1691 G4double Eta0Chi = 0.8;
1692 G4double EtaChi = Eta0Chi * ( 1.0 + 6.02*std::pow( ZMaterial,-1.19 ) );
1693 G4double W = ( EtaChi * std::pow( ZMaterial,1.0/6.0 ) ) / std::sqrt(X);
1694 G4double FunctionOfW = FTable[46][1]*FTable[46][0]/W ;
1695
1696 for(G4int j=0; j<47; j++) {
1697
1698 if( W < FTable[j][0] ) {
1699
1700 if(0 == j) {
1701 FunctionOfW = FTable[0][1] ;
1702
1703 } else {
1704 FunctionOfW = (FTable[j][1] - FTable[j-1][1]) * (W - FTable[j-1][0])
1705 / (FTable[j][0] - FTable[j-1][0])
1706 + FTable[j-1][1] ;
1707 }
1708
1709 break;
1710 }
1711
1712 }
1713
1714 BarkasTerm += FunctionOfW /( std::sqrt(ZMaterial * X) * X);
1715 }
1716
1717 BarkasTerm *= twopi_mc2_rcl2 * (material->GetElectronDensity()) / beta2 ;
1718
1719 return BarkasTerm;
1720}
1721
1722//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1723
1724G4double G4hLowEnergyIonisation::BlochTerm(const G4Material* material,
1725 G4double kineticEnergy,
1726 G4double cSquare) const
1727//Function to compute the Bloch term for protons:
1728//
1729//Ref. Z_1^3 effect in the stopping power of matter for charged particles
1730// J.C Ashley and R.H.Ritchie
1731// Physical review B Vol.5 No.7 1 April 1972 pagg. 2393-2397
1732//
1733{
1734 G4double eloss = 0.0 ;
1735 G4double gamma = 1.0 + kineticEnergy / proton_mass_c2 ;
1736 G4double beta2 = 1.0 - 1.0/(gamma*gamma) ;
1737 G4double y = cSquare / (137.0*137.0*beta2) ;
1738
1739 if(y < 0.05) {
1740 eloss = 1.202 ;
1741
1742 } else {
1743 eloss = 1.0 / (1.0 + y) ;
1744 G4double de = eloss ;
1745
1746 for(G4int i=2; de>eloss*0.01; i++) {
1747 de = 1.0/( i * (i*i + y)) ;
1748 eloss += de ;
1749 }
1750 }
1751 eloss *= -1.0 * y * cSquare * twopi_mc2_rcl2 *
1752 (material->GetElectronDensity()) / beta2 ;
1753
1754 return eloss;
1755}
1756
1757//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1758
1759G4double G4hLowEnergyIonisation::ElectronicLossFluctuation(
1760 const G4DynamicParticle* particle,
1761 const G4MaterialCutsCouple* couple,
1762 G4double meanLoss,
1763 G4double step) const
1764// calculate actual loss from the mean loss
1765// The model used to get the fluctuation is essentially the same
1766// as in Glandz in Geant3.
1767{
1768 // data members to speed up the fluctuation calculation
1769 // G4int imat ;
1770 // G4double f1Fluct,f2Fluct,e1Fluct,e2Fluct,rateFluct,ipotFluct;
1771 // G4double e1LogFluct,e2LogFluct,ipotLogFluct;
1772
1773 static const G4double minLoss = 1.*eV ;
1774 static const G4double kappa = 10. ;
1775 static const G4double theBohrBeta2 = 50.0 * keV/proton_mass_c2 ;
1776
1777 const G4Material* material = couple->GetMaterial();
1778 G4int imaterial = couple->GetIndex() ;
1779 G4double ipotFluct = material->GetIonisation()->GetMeanExcitationEnergy() ;
1780 G4double electronDensity = material->GetElectronDensity() ;
1781 G4double zeff = electronDensity/(material->GetTotNbOfAtomsPerVolume()) ;
1782
1783 // get particle data
1784 G4double tkin = particle->GetKineticEnergy();
1785 G4double particleMass = particle->GetMass() ;
1786 G4double deltaCutInKineticEnergyNow = cutForDelta[imaterial];
1787
1788 // shortcut for very very small loss
1789 if(meanLoss < minLoss) return meanLoss ;
1790
1791 // Validity range for delta electron cross section
1792 G4double threshold = std::max(deltaCutInKineticEnergyNow,ipotFluct);
1793 G4double loss, siga;
1794
1795 G4double rmass = electron_mass_c2/particleMass;
1796 G4double tau = tkin/particleMass;
1797 G4double tau1 = tau+1.0;
1798 G4double tau2 = tau*(tau+2.);
1799 G4double tmax = 2.*electron_mass_c2*tau2/(1.+2.*tau1*rmass+rmass*rmass);
1800
1801
1802 if(tmax > threshold) tmax = threshold;
1803 G4double beta2 = tau2/(tau1*tau1);
1804
1805 // Gaussian fluctuation
1806 if(meanLoss > kappa*tmax || tmax < kappa*ipotFluct )
1807 {
1808 siga = tmax * (1.0-0.5*beta2) * step * twopi_mc2_rcl2
1809 * electronDensity / beta2 ;
1810
1811 // High velocity or negatively charged particle
1812 if( beta2 > 3.0*theBohrBeta2*zeff || charge < 0.0) {
1813 siga = std::sqrt( siga * chargeSquare ) ;
1814
1815 // Low velocity - additional ion charge fluctuations according to
1816 // Q.Yang et al., NIM B61(1991)149-155.
1817 } else {
1818 G4double chu = theIonChuFluctuationModel->TheValue(particle, material);
1819 G4double yang = theIonYangFluctuationModel->TheValue(particle, material);
1820 siga = std::sqrt( siga * (chargeSquare * chu + yang)) ;
1821 }
1822
1823 do {
1824 loss = G4RandGauss::shoot(meanLoss,siga);
1825 } while (loss < 0. || loss > 2.0*meanLoss);
1826 return loss;
1827 }
1828
1829 // Non Gaussian fluctuation
1830 static const G4double probLim = 0.01 ;
1831 static const G4double sumaLim = -std::log(probLim) ;
1832 static const G4double alim = 10.;
1833
1834 G4double suma,w1,w2,C,e0,lossc,w;
1835 G4double a1,a2,a3;
1836 G4int p1,p2,p3;
1837 G4int nb;
1838 G4double corrfac, na,alfa,rfac,namean,sa,alfa1,ea,sea;
1839 G4double dp3;
1840
1841 G4double f1Fluct = material->GetIonisation()->GetF1fluct();
1842 G4double f2Fluct = material->GetIonisation()->GetF2fluct();
1843 G4double e1Fluct = material->GetIonisation()->GetEnergy1fluct();
1844 G4double e2Fluct = material->GetIonisation()->GetEnergy2fluct();
1845 G4double e1LogFluct = material->GetIonisation()->GetLogEnergy1fluct();
1846 G4double e2LogFluct = material->GetIonisation()->GetLogEnergy2fluct();
1847 G4double rateFluct = material->GetIonisation()->GetRateionexcfluct();
1848 G4double ipotLogFluct= material->GetIonisation()->GetLogMeanExcEnergy();
1849
1850 w1 = tmax/ipotFluct;
1851 w2 = std::log(2.*electron_mass_c2*tau2);
1852
1853 C = meanLoss*(1.-rateFluct)/(w2-ipotLogFluct-beta2);
1854
1855 a1 = C*f1Fluct*(w2-e1LogFluct-beta2)/e1Fluct;
1856 a2 = C*f2Fluct*(w2-e2LogFluct-beta2)/e2Fluct;
1857 a3 = rateFluct*meanLoss*(tmax-ipotFluct)/(ipotFluct*tmax*std::log(w1));
1858 if(a1 < 0.0) a1 = 0.0;
1859 if(a2 < 0.0) a2 = 0.0;
1860 if(a3 < 0.0) a3 = 0.0;
1861
1862 suma = a1+a2+a3;
1863
1864 loss = 0.;
1865
1866
1867 if(suma < sumaLim) // very small Step
1868 {
1869 e0 = material->GetIonisation()->GetEnergy0fluct();
1870
1871 if(tmax == ipotFluct)
1872 {
1873 a3 = meanLoss/e0;
1874
1875 if(a3>alim)
1876 {
1877 siga=std::sqrt(a3) ;
1878 p3 = std::max(0,G4int(G4RandGauss::shoot(a3,siga)+0.5));
1879 }
1880 else
1881 p3 = G4Poisson(a3);
1882
1883 loss = p3*e0 ;
1884
1885 if(p3 > 0)
1886 loss += (1.-2.*G4UniformRand())*e0 ;
1887
1888 }
1889 else
1890 {
1891 tmax = tmax-ipotFluct+e0 ;
1892 a3 = meanLoss*(tmax-e0)/(tmax*e0*std::log(tmax/e0));
1893
1894 if(a3>alim)
1895 {
1896 siga=std::sqrt(a3) ;
1897 p3 = std::max(0,int(G4RandGauss::shoot(a3,siga)+0.5));
1898 }
1899 else
1900 p3 = G4Poisson(a3);
1901
1902 if(p3 > 0)
1903 {
1904 w = (tmax-e0)/tmax ;
1905 if(p3 > nmaxCont2)
1906 {
1907 dp3 = G4float(p3) ;
1908 corrfac = dp3/G4float(nmaxCont2) ;
1909 p3 = nmaxCont2 ;
1910 }
1911 else
1912 corrfac = 1. ;
1913
1914 for(G4int i=0; i<p3; i++) loss += 1./(1.-w*G4UniformRand()) ;
1915 loss *= e0*corrfac ;
1916 }
1917 }
1918 }
1919
1920 else // not so small Step
1921 {
1922 // excitation type 1
1923 if(a1>alim)
1924 {
1925 siga=std::sqrt(a1) ;
1926 p1 = std::max(0,G4int(G4RandGauss::shoot(a1,siga)+0.5));
1927 }
1928 else
1929 p1 = G4Poisson(a1);
1930
1931 // excitation type 2
1932 if(a2>alim)
1933 {
1934 siga=std::sqrt(a2) ;
1935 p2 = std::max(0,G4int(G4RandGauss::shoot(a2,siga)+0.5));
1936 }
1937 else
1938 p2 = G4Poisson(a2);
1939
1940 loss = p1*e1Fluct+p2*e2Fluct;
1941
1942 // smearing to avoid unphysical peaks
1943 if(p2 > 0)
1944 loss += (1.-2.*G4UniformRand())*e2Fluct;
1945 else if (loss>0.)
1946 loss += (1.-2.*G4UniformRand())*e1Fluct;
1947
1948 // ionisation .......................................
1949 if(a3 > 0.)
1950 {
1951 if(a3>alim)
1952 {
1953 siga=std::sqrt(a3) ;
1954 p3 = std::max(0,G4int(G4RandGauss::shoot(a3,siga)+0.5));
1955 }
1956 else
1957 p3 = G4Poisson(a3);
1958
1959 lossc = 0.;
1960 if(p3 > 0)
1961 {
1962 na = 0.;
1963 alfa = 1.;
1964 if (p3 > nmaxCont2)
1965 {
1966 dp3 = G4float(p3);
1967 rfac = dp3/(G4float(nmaxCont2)+dp3);
1968 namean = G4float(p3)*rfac;
1969 sa = G4float(nmaxCont1)*rfac;
1970 na = G4RandGauss::shoot(namean,sa);
1971 if (na > 0.)
1972 {
1973 alfa = w1*G4float(nmaxCont2+p3)/
1974 (w1*G4float(nmaxCont2)+G4float(p3));
1975 alfa1 = alfa*std::log(alfa)/(alfa-1.);
1976 ea = na*ipotFluct*alfa1;
1977 sea = ipotFluct*std::sqrt(na*(alfa-alfa1*alfa1));
1978 lossc += G4RandGauss::shoot(ea,sea);
1979 }
1980 }
1981
1982 nb = G4int(G4float(p3)-na);
1983 if (nb > 0)
1984 {
1985 w2 = alfa*ipotFluct;
1986 w = (tmax-w2)/tmax;
1987 for (G4int k=0; k<nb; k++) lossc += w2/(1.-w*G4UniformRand());
1988 }
1989 }
1990 loss += lossc;
1991 }
1992 }
1993
1994 return loss ;
1995}
1996
1997//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
1998
1999void G4hLowEnergyIonisation::SetCutForSecondaryPhotons(G4double cut)
2000{
2001 minGammaEnergy = cut;
2002}
2003
2004//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
2005
2006void G4hLowEnergyIonisation::SetCutForAugerElectrons(G4double cut)
2007{
2008 minElectronEnergy = cut;
2009}
2010
2011//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
2012
2013void G4hLowEnergyIonisation::ActivateAugerElectronProduction(G4bool val)
2014{
2015 deexcitationManager.ActivateAugerElectronProduction(val);
2016}
2017
2018//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
2019
2020void G4hLowEnergyIonisation::PrintInfoDefinition() const
2021{
2022 G4String comments = " Knock-on electron cross sections . ";
2023 comments += "\n Good description above the mean excitation energy.\n";
2024 comments += " Delta ray energy sampled from differential Xsection.";
2025
2026 G4cout << G4endl << GetProcessName() << ": " << comments
2027 << "\n PhysicsTables from " << LowestKineticEnergy / eV << " eV "
2028 << " to " << HighestKineticEnergy / TeV << " TeV "
2029 << " in " << TotBin << " bins."
2030 << "\n Electronic stopping power model is "
2031 << theProtonTable
2032 << "\n from " << protonLowEnergy / keV << " keV "
2033 << " to " << protonHighEnergy / MeV << " MeV " << "." << G4endl ;
2034 G4cout << "\n Parametrisation model for antiprotons is "
2035 << theAntiProtonTable
2036 << "\n from " << antiProtonLowEnergy / keV << " keV "
2037 << " to " << antiProtonHighEnergy / MeV << " MeV " << "." << G4endl ;
2038 if(theBarkas){
2039 G4cout << " Parametrization of the Barkas effect is switched on."
2040 << G4endl ;
2041 }
2042 if(nStopping) {
2043 G4cout << " Nuclear stopping power model is " << theNuclearTable
2044 << G4endl ;
2045 }
2046
2047 G4bool printHead = true;
2048
2049 const G4ProductionCutsTable* theCoupleTable=
2050 G4ProductionCutsTable::GetProductionCutsTable();
2051 size_t numOfCouples = theCoupleTable->GetTableSize();
2052
2053 // loop for materials
2054
2055 for (size_t j=0 ; j < numOfCouples; j++) {
2056
2057 const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(j);
2058 const G4Material* material= couple->GetMaterial();
2059 G4double deltaCutNow = cutForDelta[(couple->GetIndex())] ;
2060 G4double eexc = material->GetIonisation()->GetMeanExcitationEnergy();
2061
2062 if(eexc > deltaCutNow) {
2063 if(printHead) {
2064 printHead = false ;
2065
2066 G4cout << " material min.delta energy(keV) " << G4endl;
2067 G4cout << G4endl;
2068 }
2069
2070 G4cout << std::setw(20) << material->GetName()
2071 << std::setw(15) << eexc/keV << G4endl;
2072 }
2073 }
2074}
Note: See TracBrowser for help on using the repository browser.