source: trunk/source/processes/electromagnetic/lowenergy/src/G4CrossSectionIonisationRuddPartial.cc@ 830

Last change on this file since 830 was 819, checked in by garnier, 17 years ago

import all except CVS
  • Property svn:executable set to *
File size: 12.6 KB
RevLine 
[819]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// $Id: G4CrossSectionIonisationRuddPartial.cc,v 1.3 2007/11/09 20:11:04 pia Exp $
28// GEANT4 tag $Name: $
29//
30// Contact Author: Maria Grazia Pia (Maria.Grazia.Pia@cern.ch)
31//
32// Reference: TNS Geant4-DNA paper
33// Reference for implementation model: NIM. 155, pp. 145-156, 1978
34
35// History:
36// -----------
37// Date Name Modification
38// 28 Apr 2007 M.G. Pia Created in compliance with design described in TNS paper
39//
40// -------------------------------------------------------------------
41
42// Class description:
43// Geant4-DNA Cross total cross section for electron elastic scattering in water
44// Reference: TNS Geant4-DNA paper
45// S. Chauvie et al., Geant4 physics processes for microdosimetry simulation:
46// design foundation and implementation of the first set of models,
47// IEEE Trans. Nucl. Sci., vol. 54, no. 6, Dec. 2007.
48// Further documentation available from http://www.ge.infn.it/geant4/dna
49
50// -------------------------------------------------------------------
51
52
53#include "G4CrossSectionIonisationRuddPartial.hh"
54#include "G4ParticleDefinition.hh"
55#include "G4Electron.hh"
56#include "G4Proton.hh"
57#include "G4Track.hh"
58#include "G4LogLogInterpolation.hh"
59#include "G4SystemOfUnits.hh"
60#include "G4DNAGenericIonsManager.hh"
61
62#include "Randomize.hh"
63
64#include <utility>
65
66G4CrossSectionIonisationRuddPartial::G4CrossSectionIonisationRuddPartial()
67{
68 name = "IonisationRudd";
69
70 // Default energy limits (defined for protection against anomalous behaviour only)
71 name = "IonisationRuddPartial";
72 lowEnergyLimitDefault = 100 * eV;
73 highEnergyLimitDefault = 100 * MeV;
74
75 G4String fileProton("dna/sigma_ionisation_p_rudd");
76 G4String fileHydrogen("dna/sigma_ionisation_h_rudd");
77 G4String fileAlphaPlusPlus("dna/sigma_ionisation_alphaplusplus_rudd");
78 G4String fileAlphaPlus("dna/sigma_ionisation_alphaplus_rudd");
79 G4String fileHelium("dna/sigma_ionisation_he_rudd");
80
81 G4DNAGenericIonsManager *instance;
82 instance = G4DNAGenericIonsManager::Instance();
83 G4ParticleDefinition* protonDef = G4Proton::ProtonDefinition();
84 G4ParticleDefinition* hydrogenDef = instance->GetIon("hydrogen");
85 G4ParticleDefinition* alphaPlusPlusDef = instance->GetIon("alpha++");
86 G4ParticleDefinition* alphaPlusDef = instance->GetIon("alpha+");
87 G4ParticleDefinition* heliumDef = instance->GetIon("helium");
88
89 G4String proton;
90 G4String hydrogen;
91 G4String alphaPlusPlus;
92 G4String alphaPlus;
93 G4String helium;
94
95 // Factor to scale microscopic/macroscopic cross section data in water
96 // ---- MGP ---- Hardcoded (taken from prototype code); to be replaced with proper calculation
97 G4double scaleFactor = 1 * m*m;
98
99 // Data members for protons
100
101 if (protonDef != 0)
102 {
103 proton = protonDef->GetParticleName();
104 tableFile[proton] = fileProton;
105
106 // Energy limits
107 lowEnergyLimit[proton] = 100. * eV;
108 highEnergyLimit[proton] = 500. * keV;
109
110 // Create data set with proton cross section data and load values stored in file
111 G4DNACrossSectionDataSet* tableProton = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, eV,scaleFactor );
112 tableProton->LoadData(fileProton);
113
114 // Insert key-table pair in map
115 tableData[proton] = tableProton;
116 }
117 else
118 {
119 G4Exception("G4CrossSectionIonisationRudd Constructor: proton is not defined");
120 }
121
122 // Data members for hydrogen
123
124 if (hydrogenDef != 0)
125 {
126 hydrogen = hydrogenDef->GetParticleName();
127 tableFile[hydrogen] = fileHydrogen;
128
129 // Energy limits
130 lowEnergyLimit[hydrogen] = 100. * eV;
131 highEnergyLimit[hydrogen] = 100. * MeV;
132
133 // Create data set with hydrogen cross section data and load values stored in file
134 G4DNACrossSectionDataSet* tableHydrogen = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, eV,scaleFactor );
135 tableHydrogen->LoadData(fileHydrogen);
136
137 // Insert key-table pair in map
138 tableData[hydrogen] = tableHydrogen;
139 }
140 else
141 {
142 G4Exception("G4CrossSectionIonisationRudd Constructor: hydrogen is not defined");
143 }
144
145 // Data members for alphaPlusPlus
146
147 if (alphaPlusPlusDef != 0)
148 {
149 alphaPlusPlus = alphaPlusPlusDef->GetParticleName();
150 tableFile[alphaPlusPlus] = fileAlphaPlusPlus;
151
152 // Energy limits
153 lowEnergyLimit[alphaPlusPlus] = 1. * keV;
154 highEnergyLimit[alphaPlusPlus] = 10. * MeV;
155
156 // Create data set with hydrogen cross section data and load values stored in file
157 G4DNACrossSectionDataSet* tableAlphaPlusPlus = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, eV,scaleFactor );
158 tableAlphaPlusPlus->LoadData(fileAlphaPlusPlus);
159
160 // Insert key-table pair in map
161 tableData[alphaPlusPlus] = tableAlphaPlusPlus;
162 }
163 else
164 {
165 G4Exception("G4CrossSectionIonisationRudd Constructor: alphaPlusPlus is not defined");
166 }
167
168 // Data members for alphaPlus
169
170 if (alphaPlusDef != 0)
171 {
172 alphaPlus = alphaPlusDef->GetParticleName();
173 tableFile[alphaPlus] = fileAlphaPlus;
174
175 // Energy limits
176 lowEnergyLimit[alphaPlus] = 1. * keV;
177 highEnergyLimit[alphaPlus] = 10. * MeV;
178
179 // Create data set with hydrogen cross section data and load values stored in file
180 G4DNACrossSectionDataSet* tableAlphaPlus = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, eV,scaleFactor );
181 tableAlphaPlus->LoadData(fileAlphaPlus);
182
183 // Insert key-table pair in map
184 tableData[alphaPlus] = tableAlphaPlus;
185 }
186 else
187 {
188 G4Exception("G4CrossSectionIonisationRudd Constructor: alphaPlus is not defined");
189 }
190
191 // Data members for helium
192
193 if (heliumDef != 0)
194 {
195 helium = heliumDef->GetParticleName();
196 tableFile[helium] = fileHelium;
197
198 // Energy limits
199 lowEnergyLimit[helium] = 1. * keV;
200 highEnergyLimit[helium] = 10. * MeV;
201
202 // Create data set with hydrogen cross section data and load values stored in file
203 G4DNACrossSectionDataSet* tableHelium = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, eV,scaleFactor );
204 tableHelium->LoadData(fileHelium);
205
206 // Insert key-table pair in map
207 tableData[helium] = tableHelium;
208 }
209 else
210 {
211 G4Exception("G4CrossSectionIonisationRudd Constructor: helium is not defined");
212 }
213}
214
215
216G4CrossSectionIonisationRuddPartial::~G4CrossSectionIonisationRuddPartial()
217{
218 // Destroy the content of the map
219 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
220 for (pos = tableData.begin(); pos != tableData.end(); ++pos)
221 {
222 G4DNACrossSectionDataSet* table = pos->second;
223 delete table;
224 }
225}
226
227
228G4int G4CrossSectionIonisationRuddPartial::RandomSelect(G4double k, const G4String& particle )
229{
230
231 // BEGIN PART 1/2 OF ELECTRON CORRECTION
232
233 // add ONE or TWO electron-water excitation for alpha+ and helium
234
235 G4DNAGenericIonsManager *instance;
236 instance = G4DNAGenericIonsManager::Instance();
237 G4double kElectron(0);
238 G4double electronComponent(0);
239 G4DNACrossSectionDataSet * electronDataset = new G4DNACrossSectionDataSet (new G4LogLogInterpolation, eV, (1./3.343e22)*m*m);
240
241 if ( particle == instance->GetIon("alpha+")->GetParticleName()
242 ||
243 particle == instance->GetIon("helium")->GetParticleName()
244 )
245 {
246 electronDataset->LoadData("dna/sigma_ionisation_e_born");
247
248 kElectron = k * 0.511/3728;
249
250 electronComponent = electronDataset->FindValue(kElectron);
251
252 }
253
254 delete electronDataset;
255
256 // END PART 1/2 OF ELECTRON CORRECTION
257
258 G4int level = 0;
259
260 // Retrieve data table corresponding to the current particle type
261
262 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
263 pos = tableData.find(particle);
264
265 if (pos != tableData.end())
266 {
267 G4DNACrossSectionDataSet* table = pos->second;
268
269 if (table != 0)
270 {
271 // C-style arrays are used in G4DNACrossSectionDataSet: this design feature was
272 // introduced without authorization and should be replaced by the use of STL containers
273
274 G4double* valuesBuffer = new G4double[table->NumberOfComponents()];
275
276 const size_t n(table->NumberOfComponents());
277 size_t i(n);
278 G4double value = 0.;
279
280 while (i>0)
281 {
282 i--;
283 valuesBuffer[i] = table->GetComponent(i)->FindValue(k);
284
285 // BEGIN PART 2/2 OF ELECTRON CORRECTION
286
287 if (particle == instance->GetIon("alpha+")->GetParticleName())
288 {valuesBuffer[i]=table->GetComponent(i)->FindValue(k) + electronComponent; }
289
290 if (particle == instance->GetIon("helium")->GetParticleName())
291 {valuesBuffer[i]=table->GetComponent(i)->FindValue(k) + 2*electronComponent; }
292
293 // BEGIN PART 2/2 OF ELECTRON CORRECTION
294
295 value += valuesBuffer[i];
296 }
297
298 value *= G4UniformRand();
299
300 i = n;
301
302 while (i > 0)
303 {
304 i--;
305
306 if (valuesBuffer[i] > value)
307 {
308 delete[] valuesBuffer;
309 return i;
310 }
311 value -= valuesBuffer[i];
312 }
313
314 // It should never end up here
315
316 // ---- MGP ---- Is the following line really necessary?
317 if (valuesBuffer) delete[] valuesBuffer;
318
319 }
320 }
321 else
322 {
323 G4Exception("G4CrossSectionIonisationRuddPartial: attempting to calculate cross section for wrong particle");
324 }
325
326 return level;
327}
328
329
330G4double G4CrossSectionIonisationRuddPartial::CrossSection(const G4Track& track )
331{
332 G4double sigma = 0.;
333
334 const G4DynamicParticle* particle = track.GetDynamicParticle();
335 G4double k = particle->GetKineticEnergy();
336
337 // Cross section = 0 outside the energy validity limits set in the constructor
338 // ---- MGP ---- Better handling of these limits to be set in a following design iteration
339
340 G4double lowLim = lowEnergyLimitDefault;
341 G4double highLim = highEnergyLimitDefault;
342
343 const G4String& particleName = particle->GetDefinition()->GetParticleName();
344
345 // Retrieve energy limits for the current particle type
346
347 std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
348 pos1 = lowEnergyLimit.find(particleName);
349
350 // Lower limit
351 if (pos1 != lowEnergyLimit.end())
352 {
353 lowLim = pos1->second;
354 }
355
356 // Upper limit
357 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
358 pos2 = highEnergyLimit.find(particleName);
359
360 if (pos2 != highEnergyLimit.end())
361 {
362 highLim = pos2->second;
363 }
364
365 // Verify that the current track is within the energy limits of validity of the cross section model
366 if (k >= lowLim && k <= highLim)
367 {
368 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
369 pos = tableData.find(particleName);
370
371 if (pos != tableData.end())
372 {
373 G4DNACrossSectionDataSet* table = pos->second;
374 if (table != 0)
375 {
376 // ---- MGP ---- Temporary
377 // table->PrintData();
378
379 // Cross section
380 sigma = table->FindValue(k);
381 }
382 }
383 else
384 {
385 // The track corresponds to a not pertinent particle
386 G4Exception("G4CrossSectionIonisationRuddPartial: attempting to calculate cross section for wrong particle");
387 }
388 }
389
390 return sigma;
391}
392
393
394G4double G4CrossSectionIonisationRuddPartial::Sum(G4double /* energy */, const G4String& /* particle */)
395{
396
397 return 0;
398}
Note: See TracBrowser for help on using the repository browser.