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

Last change on this file since 1347 was 1347, checked in by garnier, 13 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.