ECCE @ EIC Software
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4LivermoreIonisationModel.cc
Go to the documentation of this file. Or view the newest version in sPHENIX GitHub for file G4LivermoreIonisationModel.cc
1 //
2 // ********************************************************************
3 // * License and Disclaimer *
4 // * *
5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. *
10 // * *
11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitation of liability. *
17 // * *
18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Software license. *
24 // ********************************************************************
25 //
26 //
27 // Author: Luciano Pandola
28 // on base of G4LowEnergyIonisation developed by A.Forti and V.Ivanchenko
29 //
30 // History:
31 // --------
32 // 12 Jan 2009 L Pandola Migration from process to model
33 // 03 Mar 2009 L Pandola Bug fix (release memory in the destructor)
34 // 15 Apr 2009 V Ivanchenko Cleanup initialisation and generation of secondaries:
35 // - apply internal high-energy limit only in constructor
36 // - do not apply low-energy limit (default is 0)
37 // - simplify sampling of deexcitation by using cut in energy
38 // - set activation of Auger "false"
39 // - remove initialisation of element selectors
40 // 19 May 2009 L Pandola Explicitely set to zero pointers deleted in
41 // Initialise(), since they might be checked later on
42 // 23 Oct 2009 L Pandola
43 // - atomic deexcitation managed via G4VEmModel::DeexcitationFlag() is
44 // set as "true" (default would be false)
45 // 12 Oct 2010 L Pandola
46 // - add debugging information about energy in
47 // SampleDeexcitationAlongStep()
48 // - generate fluorescence SampleDeexcitationAlongStep() only above
49 // the cuts.
50 // 01 Jun 2011 V Ivanchenko general cleanup - all old deexcitation code removed
51 //
52 
53 #include "G4LivermoreIonisationModel.hh"
54 #include "G4PhysicalConstants.hh"
55 #include "G4SystemOfUnits.hh"
56 #include "G4ParticleDefinition.hh"
57 #include "G4MaterialCutsCouple.hh"
58 #include "G4ProductionCutsTable.hh"
59 #include "G4DynamicParticle.hh"
60 #include "G4Element.hh"
61 #include "G4ParticleChangeForLoss.hh"
62 #include "G4Electron.hh"
63 #include "G4CrossSectionHandler.hh"
64 #include "G4VEMDataSet.hh"
65 #include "G4eIonisationCrossSectionHandler.hh"
66 #include "G4eIonisationSpectrum.hh"
67 #include "G4VEnergySpectrum.hh"
68 #include "G4SemiLogInterpolation.hh"
69 #include "G4AtomicTransitionManager.hh"
70 #include "G4AtomicShell.hh"
71 #include "G4DeltaAngle.hh"
72 
73 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
74 
75 
76 G4LivermoreIonisationModel::G4LivermoreIonisationModel(const G4ParticleDefinition*,
77  const G4String& nam) :
78  G4VEmModel(nam), fParticleChange(0),
79  isInitialised(false),
80  crossSectionHandler(0), energySpectrum(0)
81 {
82  fIntrinsicLowEnergyLimit = 12.*eV;
83  fIntrinsicHighEnergyLimit = 100.0*GeV;
84 
85  verboseLevel = 0;
86  SetAngularDistribution(new G4DeltaAngle());
87 
88  transitionManager = G4AtomicTransitionManager::Instance();
89 }
90 
91 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
92 
93 G4LivermoreIonisationModel::~G4LivermoreIonisationModel()
94 {
95  delete energySpectrum;
96  delete crossSectionHandler;
97 }
98 
99 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
100 
101 void G4LivermoreIonisationModel::Initialise(const G4ParticleDefinition* particle,
102  const G4DataVector& cuts)
103 {
104  //Check that the Livermore Ionisation is NOT attached to e+
105  if (particle != G4Electron::Electron())
106  {
107  G4Exception("G4LivermoreIonisationModel::Initialise",
108  "em0002",FatalException,
109  "Livermore Ionisation Model is applicable only to electrons");
110  }
111 
112  transitionManager->Initialise();
113 
114  //Read energy spectrum
115  if (energySpectrum)
116  {
117  delete energySpectrum;
118  energySpectrum = 0;
119  }
120  energySpectrum = new G4eIonisationSpectrum();
121  if (verboseLevel > 3)
122  G4cout << "G4VEnergySpectrum is initialized" << G4endl;
123 
124  //Initialize cross section handler
125  if (crossSectionHandler)
126  {
127  delete crossSectionHandler;
128  crossSectionHandler = 0;
129  }
130 
131  const size_t nbins = 20;
132  G4double emin = LowEnergyLimit();
133  G4double emax = HighEnergyLimit();
134  G4int ndec = G4int(std::log10(emax/emin) + 0.5);
135  if(ndec <= 0) { ndec = 1; }
136 
137  G4VDataSetAlgorithm* interpolation = new G4SemiLogInterpolation();
138  crossSectionHandler =
139  new G4eIonisationCrossSectionHandler(energySpectrum,interpolation,
140  emin,emax,nbins*ndec);
141  crossSectionHandler->Clear();
142  crossSectionHandler->LoadShellData("ioni/ion-ss-cs-");
143  //This is used to retrieve cross section values later on
144  G4VEMDataSet* emdata =
145  crossSectionHandler->BuildMeanFreePathForMaterials(&cuts);
146  //The method BuildMeanFreePathForMaterials() is required here only to force
147  //the building of an internal table: the output pointer can be deleted
148  delete emdata;
149 
150  if (verboseLevel > 0)
151  {
152  G4cout << "Livermore Ionisation model is initialized " << G4endl
153  << "Energy range: "
154  << LowEnergyLimit() / keV << " keV - "
155  << HighEnergyLimit() / GeV << " GeV"
156  << G4endl;
157  }
158 
159  if (verboseLevel > 3)
160  {
161  G4cout << "Cross section data: " << G4endl;
162  crossSectionHandler->PrintData();
163  G4cout << "Parameters: " << G4endl;
164  energySpectrum->PrintData();
165  }
166 
167  if(isInitialised) { return; }
168  fParticleChange = GetParticleChangeForLoss();
169  isInitialised = true;
170 }
171 
172 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
173 
174 G4double
175 G4LivermoreIonisationModel::ComputeCrossSectionPerAtom(
176  const G4ParticleDefinition*,
177  G4double energy,
178  G4double Z, G4double,
179  G4double cutEnergy,
180  G4double)
181 {
182  G4int iZ = (G4int) Z;
183  if (!crossSectionHandler)
184  {
185  G4Exception("G4LivermoreIonisationModel::ComputeCrossSectionPerAtom",
186  "em1007",FatalException,
187  "The cross section handler is not correctly initialized");
188  return 0;
189  }
190 
191  //The cut is already included in the crossSectionHandler
192  G4double cs =
193  crossSectionHandler->GetCrossSectionAboveThresholdForElement(energy,
194  cutEnergy,
195  iZ);
196 
197  if (verboseLevel > 1)
198  {
199  G4cout << "G4LivermoreIonisationModel " << G4endl;
200  G4cout << "Cross section for delta emission > "
201  << cutEnergy/keV << " keV at "
202  << energy/keV << " keV and Z = " << iZ << " --> "
203  << cs/barn << " barn" << G4endl;
204  }
205  return cs;
206 }
207 
208 
209 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
210 
211 G4double
212 G4LivermoreIonisationModel::ComputeDEDXPerVolume(const G4Material* material,
213  const G4ParticleDefinition*,
214  G4double kineticEnergy,
215  G4double cutEnergy)
216 {
217  G4double sPower = 0.0;
218 
219  const G4ElementVector* theElementVector = material->GetElementVector();
220  size_t NumberOfElements = material->GetNumberOfElements() ;
221  const G4double* theAtomicNumDensityVector =
222  material->GetAtomicNumDensityVector();
223 
224  // loop for elements in the material
225  for (size_t iel=0; iel<NumberOfElements; iel++ )
226  {
227  G4int iZ = (G4int)((*theElementVector)[iel]->GetZ());
228  G4int nShells = transitionManager->NumberOfShells(iZ);
229  for (G4int n=0; n<nShells; n++)
230  {
231  G4double e = energySpectrum->AverageEnergy(iZ, 0.0,cutEnergy,
232  kineticEnergy, n);
233  G4double cs= crossSectionHandler->FindValue(iZ,kineticEnergy, n);
234  sPower += e * cs * theAtomicNumDensityVector[iel];
235  }
236  G4double esp = energySpectrum->Excitation(iZ,kineticEnergy);
237  sPower += esp * theAtomicNumDensityVector[iel];
238  }
239 
240  if (verboseLevel > 2)
241  {
242  G4cout << "G4LivermoreIonisationModel " << G4endl;
243  G4cout << "Stopping power < " << cutEnergy/keV
244  << " keV at " << kineticEnergy/keV << " keV = "
245  << sPower/(keV/mm) << " keV/mm" << G4endl;
246  }
247 
248  return sPower;
249 }
250 
251 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
252 
253 void G4LivermoreIonisationModel::SampleSecondaries(
254  std::vector<G4DynamicParticle*>* fvect,
255  const G4MaterialCutsCouple* couple,
256  const G4DynamicParticle* aDynamicParticle,
257  G4double cutE,
258  G4double maxE)
259 {
260 
261  G4double kineticEnergy = aDynamicParticle->GetKineticEnergy();
262 
263  if (kineticEnergy <= fIntrinsicLowEnergyLimit)
264  {
265  fParticleChange->SetProposedKineticEnergy(0.);
266  fParticleChange->ProposeLocalEnergyDeposit(kineticEnergy);
267  return;
268  }
269 
270  // Select atom and shell
271  G4int Z = crossSectionHandler->SelectRandomAtom(couple, kineticEnergy);
272  G4int shellIndex = crossSectionHandler->SelectRandomShell(Z, kineticEnergy);
273  const G4AtomicShell* shell = transitionManager->Shell(Z,shellIndex);
274  G4double bindingEnergy = shell->BindingEnergy();
275 
276  // Sample delta energy using energy interval for delta-electrons
277  G4double energyMax =
278  std::min(maxE,energySpectrum->MaxEnergyOfSecondaries(kineticEnergy));
279  G4double energyDelta = energySpectrum->SampleEnergy(Z, cutE, energyMax,
280  kineticEnergy, shellIndex);
281 
282  if (energyDelta == 0.) //nothing happens
283  { return; }
284 
285  const G4ParticleDefinition* electron = G4Electron::Electron();
286  G4DynamicParticle* delta = new G4DynamicParticle(electron,
287  GetAngularDistribution()->SampleDirectionForShell(aDynamicParticle, energyDelta,
288  Z, shellIndex,
289  couple->GetMaterial()),
290  energyDelta);
291 
292  fvect->push_back(delta);
293 
294  // Change kinematics of primary particle
295  G4ThreeVector direction = aDynamicParticle->GetMomentumDirection();
296  G4double totalMomentum = std::sqrt(kineticEnergy*(kineticEnergy + 2*electron_mass_c2));
297 
298  G4ThreeVector finalP = totalMomentum*direction - delta->GetMomentum();
299  finalP = finalP.unit();
300 
301  //This is the amount of energy available for fluorescence
302  G4double theEnergyDeposit = bindingEnergy;
303 
304  // fill ParticleChange
305  // changed energy and momentum of the actual particle
306  G4double finalKinEnergy = kineticEnergy - energyDelta - theEnergyDeposit;
307  if(finalKinEnergy < 0.0)
308  {
309  theEnergyDeposit += finalKinEnergy;
310  finalKinEnergy = 0.0;
311  }
312  else
313  {
314  fParticleChange->ProposeMomentumDirection(finalP);
315  }
316  fParticleChange->SetProposedKineticEnergy(finalKinEnergy);
317 
318  if (theEnergyDeposit < 0)
319  {
320  G4cout << "G4LivermoreIonisationModel: Negative energy deposit: "
321  << theEnergyDeposit/eV << " eV" << G4endl;
322  theEnergyDeposit = 0.0;
323  }
324 
325  //Assign local energy deposit
326  fParticleChange->ProposeLocalEnergyDeposit(theEnergyDeposit);
327 
328  if (verboseLevel > 1)
329  {
330  G4cout << "-----------------------------------------------------------" << G4endl;
331  G4cout << "Energy balance from G4LivermoreIonisation" << G4endl;
332  G4cout << "Incoming primary energy: " << kineticEnergy/keV << " keV" << G4endl;
333  G4cout << "-----------------------------------------------------------" << G4endl;
334  G4cout << "Outgoing primary energy: " << finalKinEnergy/keV << " keV" << G4endl;
335  G4cout << "Delta ray " << energyDelta/keV << " keV" << G4endl;
336  G4cout << "Fluorescence: " << (bindingEnergy-theEnergyDeposit)/keV << " keV" << G4endl;
337  G4cout << "Local energy deposit " << theEnergyDeposit/keV << " keV" << G4endl;
338  G4cout << "Total final state: " << (finalKinEnergy+energyDelta+bindingEnergy)
339  << " keV" << G4endl;
340  G4cout << "-----------------------------------------------------------" << G4endl;
341  }
342  return;
343 }
344 
345 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
346