source: trunk/source/processes/electromagnetic/lowenergy/src/G4CrossSectionChargeIncreasePartial.cc@ 1057

Last change on this file since 1057 was 1007, checked in by garnier, 17 years ago

update to geant4.9.2
  • Property svn:executable set to *
File size: 6.9 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: G4CrossSectionChargeIncreasePartial.cc,v 1.2 2008/07/14 20:47:34 sincerti Exp $
27// GEANT4 tag $Name: geant4-09-02 $
28
29#include "G4CrossSectionChargeIncreasePartial.hh"
30
31//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
32
33G4CrossSectionChargeIncreasePartial::G4CrossSectionChargeIncreasePartial()
34{
35 //ALPHA+
36
37 f0[0][0]=1.;
38 a0[0][0]=2.25;
39 a1[0][0]=-0.75;
40 b0[0][0]=-32.10;
41 c0[0][0]=0.600;
42 d0[0][0]=2.40;
43 x0[0][0]=4.60;
44
45 x1[0][0]=-1.;
46 b1[0][0]=-1.;
47
48 numberOfPartialCrossSections[0]=1;
49
50 //HELIUM
51
52 f0[0][1]=1.;
53 a0[0][1]=2.25;
54 a1[0][1]=-0.75;
55 b0[0][1]=-30.93;
56 c0[0][1]=0.590;
57 d0[0][1]=2.35;
58 x0[0][1]=4.29;
59
60 f0[1][1]=1.;
61 a0[1][1]=2.25;
62 a1[1][1]=-0.75;
63 b0[1][1]=-32.61;
64 c0[1][1]=0.435;
65 d0[1][1]=2.70;
66 x0[1][1]=4.45;
67
68 x1[0][1]=-1.;
69 b1[0][1]=-1.;
70
71 x1[1][1]=-1.;
72 b1[1][1]=-1.;
73
74 numberOfPartialCrossSections[1]=2;
75}
76
77//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
78
79G4CrossSectionChargeIncreasePartial::~G4CrossSectionChargeIncreasePartial()
80{}
81
82//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
83
84G4double G4CrossSectionChargeIncreasePartial::CrossSection(G4double k, G4int index,
85 const G4ParticleDefinition* particleDefinition)
86{
87
88 G4int particleTypeIndex = 0;
89 G4DNAGenericIonsManager *instance;
90 instance = G4DNAGenericIonsManager::Instance();
91
92 if (particleDefinition == instance->GetIon("alpha+")) particleTypeIndex=0;
93 if (particleDefinition == instance->GetIon("helium")) particleTypeIndex=1;
94
95 //
96 // sigma(T) = f0 10 ^ y(log10(T/eV))
97 //
98 // / a0 x + b0 x < x0
99 // |
100 // y(x) = < a0 x + b0 - c0 (x - x0)^d0 x0 <= x < x1
101 // |
102 // \ a1 x + b1 x >= x1
103 //
104 //
105 // f0, a0, a1, b0, b1, c0, d0, x0, x1 are parameters that change for protons and helium (0, +, ++)
106 //
107 // f0 has been added to the code in order to manage partial (shell-dependent) cross sections (if no shell dependence is present. f0=1. Sum of f0 over the considered shells should give 1)
108 //
109 // From Rad. Phys. and Chem. 59 (2000) 255-275, M. Dingfelder et al.
110 // Inelastic-collision cross sections of liquid water for interactions of energetic proton
111 //
112
113 if (x1[index][particleTypeIndex]<x0[index][particleTypeIndex])
114 {
115 //
116 // if x1 < x0 means that x1 and b1 will be calculated with the following formula (this piece of code is run on all alphas and not on protons)
117 //
118 // x1 = x0 + ((a0 - a1)/(c0 * d0)) ^ (1 / (d0 - 1))
119 //
120 // b1 = (a0 - a1) * x1 + b0 - c0 * (x1 - x0) ^ d0
121 //
122
123 x1[index][particleTypeIndex]=x0[index][particleTypeIndex] + std::pow((a0[index][particleTypeIndex] - a1[index][particleTypeIndex]) / (c0[index][particleTypeIndex] * d0[index][particleTypeIndex]), 1. / (d0[index][particleTypeIndex] - 1.));
124 b1[index][particleTypeIndex]=(a0[index][particleTypeIndex] - a1[index][particleTypeIndex]) * x1[index][particleTypeIndex] + b0[index][particleTypeIndex] - c0[index][particleTypeIndex] * std::pow(x1[index][particleTypeIndex] - x0[index][particleTypeIndex], d0[index][particleTypeIndex]);
125 }
126
127 G4double x(std::log10(k/eV));
128 G4double y;
129
130 if (x<x0[index][particleTypeIndex])
131 y=a0[index][particleTypeIndex] * x + b0[index][particleTypeIndex];
132 else if (x<x1[index][particleTypeIndex])
133 y=a0[index][particleTypeIndex] * x + b0[index][particleTypeIndex] - c0[index][particleTypeIndex] * std::pow(x - x0[index][particleTypeIndex], d0[index][particleTypeIndex]);
134 else
135 y=a1[index][particleTypeIndex] * x + b1[index][particleTypeIndex];
136
137 return f0[index][particleTypeIndex] * std::pow(10., y)*m*m;
138
139}
140
141//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
142
143G4int G4CrossSectionChargeIncreasePartial::RandomSelect(G4double k,
144 const G4ParticleDefinition* particleDefinition)
145{
146
147 G4int particleTypeIndex = 0;
148 G4DNAGenericIonsManager *instance;
149 instance = G4DNAGenericIonsManager::Instance();
150
151 if (particleDefinition == instance->GetIon("hydrogen")) return 0;
152 if (particleDefinition == instance->GetIon("alpha+")) particleTypeIndex=0;
153 if (particleDefinition == instance->GetIon("helium")) particleTypeIndex=1;
154
155 const G4int n = numberOfPartialCrossSections[particleTypeIndex];
156 G4double* values(new G4double[n]);
157 G4double value = 0;
158 G4int i = n;
159
160 while (i>0)
161 {
162 i--;
163 values[i]=CrossSection(k, i, particleDefinition);
164 value+=values[i];
165 }
166
167 value*=G4UniformRand();
168
169 i=n;
170 while (i>0)
171 {
172 i--;
173
174 if (values[i]>value)
175 break;
176
177 value-=values[i];
178 }
179
180 delete[] values;
181
182 return i;
183}
184
185//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
186
187G4double G4CrossSectionChargeIncreasePartial::Sum(G4double k, const G4ParticleDefinition* particleDefinition)
188{
189 G4int particleTypeIndex = 0;
190 G4DNAGenericIonsManager *instance;
191 instance = G4DNAGenericIonsManager::Instance();
192
193 if (particleDefinition == instance->GetIon("alpha+")) particleTypeIndex=0;
194 if (particleDefinition == instance->GetIon("helium")) particleTypeIndex=1;
195
196 G4double totalCrossSection = 0.;
197
198 for (G4int i=0; i<numberOfPartialCrossSections[particleTypeIndex]; i++)
199 {
200 totalCrossSection += CrossSection(k,i,particleDefinition);
201 }
202 return totalCrossSection;
203}
204
Note: See TracBrowser for help on using the repository browser.