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

Last change on this file was 1347, checked in by garnier, 15 years ago

geant4 tag 9.4
  • Property svn:executable set to *
File size: 9.8 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// $Id: G4CrossSectionIonisationRuddPartial.cc,v 1.5 2009/06/10 13:32:36 mantero Exp $
27// GEANT4 tag $Name: geant4-09-04-ref-00 $
28
29#include "G4CrossSectionIonisationRuddPartial.hh"
30
31//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
32
33G4CrossSectionIonisationRuddPartial::G4CrossSectionIonisationRuddPartial()
34{
35 lowEnergyLimitDefault = 100 * eV;
36 highEnergyLimitDefault = 100 * MeV;
37
38 G4String fileProton("dna/sigma_ionisation_p_rudd");
39 G4String fileHydrogen("dna/sigma_ionisation_h_rudd");
40 G4String fileAlphaPlusPlus("dna/sigma_ionisation_alphaplusplus_rudd");
41 G4String fileAlphaPlus("dna/sigma_ionisation_alphaplus_rudd");
42 G4String fileHelium("dna/sigma_ionisation_he_rudd");
43
44 G4DNAGenericIonsManager *instance;
45 instance = G4DNAGenericIonsManager::Instance();
46 G4ParticleDefinition* protonDef = G4Proton::ProtonDefinition();
47 G4ParticleDefinition* hydrogenDef = instance->GetIon("hydrogen");
48 G4ParticleDefinition* alphaPlusPlusDef = instance->GetIon("alpha++");
49 G4ParticleDefinition* alphaPlusDef = instance->GetIon("alpha+");
50 G4ParticleDefinition* heliumDef = instance->GetIon("helium");
51
52 G4String proton;
53 G4String hydrogen;
54 G4String alphaPlusPlus;
55 G4String alphaPlus;
56 G4String helium;
57
58 G4double scaleFactor = 1 * m*m;
59
60 if (protonDef != 0)
61 {
62 proton = protonDef->GetParticleName();
63 tableFile[proton] = fileProton;
64
65 lowEnergyLimit[proton] = 100. * eV;
66 highEnergyLimit[proton] = 500. * keV;
67
68 G4DNACrossSectionDataSet* tableProton = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, eV,scaleFactor );
69 tableProton->LoadData(fileProton);
70
71 tableData[proton] = tableProton;
72 }
73 else
74 {
75 G4Exception("G4CrossSectionIonisationRudd Constructor: proton is not defined");
76 }
77
78 if (hydrogenDef != 0)
79 {
80 hydrogen = hydrogenDef->GetParticleName();
81 tableFile[hydrogen] = fileHydrogen;
82
83 lowEnergyLimit[hydrogen] = 100. * eV;
84 highEnergyLimit[hydrogen] = 100. * MeV;
85
86 G4DNACrossSectionDataSet* tableHydrogen = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, eV,scaleFactor );
87 tableHydrogen->LoadData(fileHydrogen);
88
89 tableData[hydrogen] = tableHydrogen;
90 }
91 else
92 {
93 G4Exception("G4CrossSectionIonisationRudd Constructor: hydrogen is not defined");
94 }
95
96 if (alphaPlusPlusDef != 0)
97 {
98 alphaPlusPlus = alphaPlusPlusDef->GetParticleName();
99 tableFile[alphaPlusPlus] = fileAlphaPlusPlus;
100
101 lowEnergyLimit[alphaPlusPlus] = 1. * keV;
102 highEnergyLimit[alphaPlusPlus] = 10. * MeV;
103
104 G4DNACrossSectionDataSet* tableAlphaPlusPlus = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, eV,scaleFactor );
105 tableAlphaPlusPlus->LoadData(fileAlphaPlusPlus);
106
107 tableData[alphaPlusPlus] = tableAlphaPlusPlus;
108 }
109 else
110 {
111 G4Exception("G4CrossSectionIonisationRudd Constructor: alphaPlusPlus is not defined");
112 }
113
114 if (alphaPlusDef != 0)
115 {
116 alphaPlus = alphaPlusDef->GetParticleName();
117 tableFile[alphaPlus] = fileAlphaPlus;
118
119 lowEnergyLimit[alphaPlus] = 1. * keV;
120 highEnergyLimit[alphaPlus] = 10. * MeV;
121
122 G4DNACrossSectionDataSet* tableAlphaPlus = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, eV,scaleFactor );
123 tableAlphaPlus->LoadData(fileAlphaPlus);
124
125 tableData[alphaPlus] = tableAlphaPlus;
126 }
127 else
128 {
129 G4Exception("G4CrossSectionIonisationRudd Constructor: alphaPlus is not defined");
130 }
131
132 if (heliumDef != 0)
133 {
134 helium = heliumDef->GetParticleName();
135 tableFile[helium] = fileHelium;
136
137 lowEnergyLimit[helium] = 1. * keV;
138 highEnergyLimit[helium] = 10. * MeV;
139
140 G4DNACrossSectionDataSet* tableHelium = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, eV,scaleFactor );
141 tableHelium->LoadData(fileHelium);
142
143 tableData[helium] = tableHelium;
144 }
145 else
146 {
147 G4Exception("G4CrossSectionIonisationRudd Constructor: helium is not defined");
148 }
149}
150
151//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
152
153G4CrossSectionIonisationRuddPartial::~G4CrossSectionIonisationRuddPartial()
154{
155 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
156 for (pos = tableData.begin(); pos != tableData.end(); ++pos)
157 {
158 G4DNACrossSectionDataSet* table = pos->second;
159 delete table;
160 }
161}
162
163//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
164
165G4int G4CrossSectionIonisationRuddPartial::RandomSelect(G4double k, const G4String& particle )
166{
167
168 // BEGIN PART 1/2 OF ELECTRON CORRECTION
169
170 // add ONE or TWO electron-water excitation for alpha+ and helium
171
172 G4DNAGenericIonsManager *instance;
173 instance = G4DNAGenericIonsManager::Instance();
174 G4double kElectron(0);
175 G4double electronComponent(0);
176 G4DNACrossSectionDataSet * electronDataset = new G4DNACrossSectionDataSet (new G4LogLogInterpolation, eV, (1./3.343e22)*m*m);
177
178 if ( particle == instance->GetIon("alpha+")->GetParticleName()
179 ||
180 particle == instance->GetIon("helium")->GetParticleName()
181 )
182 {
183 electronDataset->LoadData("dna/sigma_ionisation_e_born");
184
185 kElectron = k * 0.511/3728;
186
187 electronComponent = electronDataset->FindValue(kElectron);
188
189 }
190
191 delete electronDataset;
192
193 // END PART 1/2 OF ELECTRON CORRECTION
194
195 G4int level = 0;
196
197 // Retrieve data table corresponding to the current particle type
198
199 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
200 pos = tableData.find(particle);
201
202 if (pos != tableData.end())
203 {
204 G4DNACrossSectionDataSet* table = pos->second;
205
206 if (table != 0)
207 {
208 G4double* valuesBuffer = new G4double[table->NumberOfComponents()];
209
210 const size_t n(table->NumberOfComponents());
211 size_t i(n);
212 G4double value = 0.;
213
214 while (i>0)
215 {
216 i--;
217 valuesBuffer[i] = table->GetComponent(i)->FindValue(k);
218
219 // BEGIN PART 2/2 OF ELECTRON CORRECTION
220
221 if (particle == instance->GetIon("alpha+")->GetParticleName())
222 {valuesBuffer[i]=table->GetComponent(i)->FindValue(k) + electronComponent; }
223
224 if (particle == instance->GetIon("helium")->GetParticleName())
225 {valuesBuffer[i]=table->GetComponent(i)->FindValue(k) + 2*electronComponent; }
226
227 // BEGIN PART 2/2 OF ELECTRON CORRECTION
228
229 value += valuesBuffer[i];
230 }
231
232 value *= G4UniformRand();
233
234 i = n;
235
236 while (i > 0)
237 {
238 i--;
239
240 if (valuesBuffer[i] > value)
241 {
242 delete[] valuesBuffer;
243 return i;
244 }
245 value -= valuesBuffer[i];
246 }
247
248 if (valuesBuffer) delete[] valuesBuffer;
249
250 }
251 }
252 else
253 {
254 G4Exception("G4CrossSectionIonisationRuddPartial: attempting to calculate cross section for wrong particle");
255 }
256
257 return level;
258}
259
260//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
261
262G4double G4CrossSectionIonisationRuddPartial::CrossSection(const G4Track& track )
263{
264 G4double sigma = 0.;
265
266 const G4DynamicParticle* particle = track.GetDynamicParticle();
267 G4double k = particle->GetKineticEnergy();
268
269 G4double lowLim = lowEnergyLimitDefault;
270 G4double highLim = highEnergyLimitDefault;
271
272 const G4String& particleName = particle->GetDefinition()->GetParticleName();
273
274 std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
275 pos1 = lowEnergyLimit.find(particleName);
276
277 if (pos1 != lowEnergyLimit.end())
278 {
279 lowLim = pos1->second;
280 }
281
282 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
283 pos2 = highEnergyLimit.find(particleName);
284
285 if (pos2 != highEnergyLimit.end())
286 {
287 highLim = pos2->second;
288 }
289
290 if (k >= lowLim && k <= highLim)
291 {
292 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
293 pos = tableData.find(particleName);
294
295 if (pos != tableData.end())
296 {
297 G4DNACrossSectionDataSet* table = pos->second;
298 if (table != 0)
299 {
300 sigma = table->FindValue(k);
301 }
302 }
303 else
304 {
305 G4Exception("G4CrossSectionIonisationRuddPartial: attempting to calculate cross section for wrong particle");
306 }
307 }
308
309 return sigma;
310}
311
312//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
313
314G4double G4CrossSectionIonisationRuddPartial::Sum(G4double /* energy */, const G4String& /* particle */)
315{
316 return 0;
317}
Note: See TracBrowser for help on using the repository browser.