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

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

update to last version 4.9.4
File size: 9.7 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: G4DNASancheExcitationModel.cc,v 1.4 2010/11/11 22:32:22 sincerti Exp $
27// GEANT4 tag $Name: geant4-09-04-ref-00 $
28//
29
30// Created by Z. Francis
31
32#include "G4DNASancheExcitationModel.hh"
33
34//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
35
36using namespace std;
37
38//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
39
40G4DNASancheExcitationModel::G4DNASancheExcitationModel(const G4ParticleDefinition*,
41 const G4String& nam)
42:G4VEmModel(nam),isInitialised(false)
43{
44
45 lowEnergyLimit = 2 * eV;
46 highEnergyLimit = 100 * eV;
47 SetLowEnergyLimit(lowEnergyLimit);
48 SetHighEnergyLimit(highEnergyLimit);
49 nLevels = 9;
50
51 verboseLevel= 0;
52 // Verbosity scale:
53 // 0 = nothing
54 // 1 = warning for energy non-conservation
55 // 2 = details of energy budget
56 // 3 = calculation of cross sections, file openings, sampling of atoms
57 // 4 = entering in methods
58
59 if (verboseLevel > 0)
60 {
61 G4cout << "Sanche Excitation model is constructed " << G4endl
62 << "Energy range: "
63 << lowEnergyLimit / eV << " eV - "
64 << highEnergyLimit / eV << " eV"
65 << G4endl;
66 }
67
68}
69
70//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
71
72G4DNASancheExcitationModel::~G4DNASancheExcitationModel()
73{}
74
75//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
76
77void G4DNASancheExcitationModel::Initialise(const G4ParticleDefinition* /*particle*/,
78 const G4DataVector& /*cuts*/)
79{
80
81
82 if (verboseLevel > 3)
83 G4cout << "Calling G4DNASancheExcitationModel::Initialise()" << G4endl;
84
85 // Energy limits
86
87 if (LowEnergyLimit() < lowEnergyLimit)
88 {
89 G4cout << "G4DNASancheExcitationModel: low energy limit increased from " <<
90 LowEnergyLimit()/eV << " eV to " << lowEnergyLimit/eV << " eV" << G4endl;
91 SetLowEnergyLimit(lowEnergyLimit);
92 }
93
94 if (HighEnergyLimit() > highEnergyLimit)
95 {
96 G4cout << "G4DNASancheExcitationModel: high energy limit decreased from " <<
97 HighEnergyLimit()/eV << " eV to " << highEnergyLimit/eV << " eV" << G4endl;
98 SetHighEnergyLimit(highEnergyLimit);
99 }
100
101 //
102
103 if (verboseLevel > 0)
104
105 G4cout << "Sanche Excitation model is initialized " << G4endl
106 << "Energy range: "
107 << LowEnergyLimit() / eV << " eV - "
108 << HighEnergyLimit() / eV << " eV"
109 << G4endl;
110
111 if(!isInitialised)
112 {
113 isInitialised = true;
114
115 if(pParticleChange)
116 fParticleChangeForGamma = reinterpret_cast<G4ParticleChangeForGamma*>(pParticleChange);
117 else
118 fParticleChangeForGamma = new G4ParticleChangeForGamma();
119 }
120
121 // InitialiseElementSelectors(particle,cuts);
122
123 char *path = getenv("G4LEDATA");
124 std::ostringstream eFullFileName;
125 eFullFileName << path << "/dna/sigma_excitationvib_e_sanche.dat";
126 std::ifstream input(eFullFileName.str().c_str());
127
128 if (!input)
129 {
130 G4Exception("G4DNASancheExcitationModel:::ERROR OPENING XS DATA FILE");
131 }
132
133 while(!input.eof())
134 {
135 double t;
136 input>>t;
137 tdummyVec.push_back(t);
138 input>>map1[t][0]>>map1[t][1]>>map1[t][2]>>map1[t][3]>>map1[t][4]>>map1[t][5]>>map1[t][6]>>map1[t][7]>>map1[t][8];
139 //G4cout<<t<<" "<<map1[t][0]<<map1[t][1]<<map1[t][2]<<map1[t][3]<<map1[t][4]<<map1[t][5]<<map1[t][6]<<map1[t][7]<<map1[t][8]<<G4endl;
140 }
141
142}
143
144//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
145
146G4double G4DNASancheExcitationModel::CrossSectionPerVolume(const G4Material* material,
147 const G4ParticleDefinition* particleDefinition,
148 G4double ekin,
149 G4double,
150 G4double)
151{
152 if (verboseLevel > 3)
153 G4cout << "Calling CrossSectionPerVolume() of G4DNASancheExcitationModel" << G4endl;
154
155 // Calculate total cross section for model
156
157 G4double sigma=0;
158
159 if (material->GetName() == "G4_WATER")
160 {
161
162 if (particleDefinition == G4Electron::ElectronDefinition())
163 {
164 if (ekin >= lowEnergyLimit && ekin < highEnergyLimit)
165 {
166 sigma = Sum(ekin);
167 }
168 }
169
170 if (verboseLevel > 3)
171 {
172 G4cout << "---> Kinetic energy(eV)=" << ekin/eV << G4endl;
173 G4cout << " - Cross section per water molecule (cm^2)=" << sigma/cm/cm << G4endl;
174 G4cout << " - Cross section per water molecule (cm^-1)=" << sigma*material->GetAtomicNumDensityVector()[1]/(1./cm) << G4endl;
175 }
176
177 } // if water
178
179
180 return sigma*2*material->GetAtomicNumDensityVector()[1];
181 // see papers for factor 2 description
182
183}
184
185//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
186
187void G4DNASancheExcitationModel::SampleSecondaries(std::vector<G4DynamicParticle*>*,
188 const G4MaterialCutsCouple*,
189 const G4DynamicParticle* aDynamicElectron,
190 G4double,
191 G4double)
192{
193
194 if (verboseLevel > 3)
195 G4cout << "Calling SampleSecondaries() of G4DNASancheExcitationModel" << G4endl;
196
197 G4double electronEnergy0 = aDynamicElectron->GetKineticEnergy();
198 G4int level = RandomSelect(electronEnergy0);
199 G4double excitationEnergy = VibrationEnergy(level); // levels go from 0 to 8
200 G4double newEnergy = electronEnergy0 - excitationEnergy;
201
202/*
203 if (electronEnergy0 < highEnergyLimit)
204 {
205 if (newEnergy >= lowEnergyLimit)
206 {
207 fParticleChangeForGamma->ProposeMomentumDirection(aDynamicElectron->GetMomentumDirection());
208 fParticleChangeForGamma->SetProposedKineticEnergy(newEnergy);
209 fParticleChangeForGamma->ProposeLocalEnergyDeposit(excitationEnergy);
210 }
211
212 else
213 {
214 fParticleChangeForGamma->ProposeTrackStatus(fStopAndKill);
215 fParticleChangeForGamma->ProposeLocalEnergyDeposit(electronEnergy0);
216 }
217 }
218*/
219
220 if (electronEnergy0 < highEnergyLimit && newEnergy>0.)
221 {
222 fParticleChangeForGamma->ProposeMomentumDirection(aDynamicElectron->GetMomentumDirection());
223 fParticleChangeForGamma->SetProposedKineticEnergy(newEnergy);
224 fParticleChangeForGamma->ProposeLocalEnergyDeposit(excitationEnergy);
225 }
226
227//
228}
229
230//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
231
232G4double G4DNASancheExcitationModel::PartialCrossSection(G4double t, G4int level)
233{
234 std::vector<double>::iterator t2 = std::upper_bound(tdummyVec.begin(),tdummyVec.end(), t/eV);
235 std::vector<double>::iterator t1 = t2-1;
236
237 double sigma = LinInterpolate((*t1), (*t2), t/eV, map1[*t1][level], map1[*t2][level]);
238 sigma*=1e-16*cm*cm;
239 if(sigma==0.)sigma=1e-30;
240 return (sigma);
241}
242
243//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
244
245G4double G4DNASancheExcitationModel::VibrationEnergy(G4int level)
246{
247 G4double energies[9] = {0.01, 0.024, 0.061, 0.092, 0.204, 0.417, 0.460, 0.500, 0.835};
248 return(energies[level]*eV);
249}
250
251//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
252
253G4int G4DNASancheExcitationModel::RandomSelect(G4double k)
254{
255
256 // Level Selection Counting can be done here !
257
258 G4int i = nLevels;
259 G4double value = 0.;
260 std::deque<double> values;
261
262 while (i > 0)
263 {
264 i--;
265 G4double partial = PartialCrossSection(k,i);
266 values.push_front(partial);
267 value += partial;
268 }
269
270 value *= G4UniformRand();
271
272 i = nLevels;
273
274 while (i > 0)
275 {
276 i--;
277 if (values[i] > value)
278 {
279 //outcount<<i<<" "<<VibrationEnergy(i)<<G4endl;
280 return i;
281 }
282 value -= values[i];
283 }
284
285 //outcount<<0<<" "<<VibrationEnergy(0)<<G4endl;
286
287 return 0;
288}
289
290//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
291
292G4double G4DNASancheExcitationModel::Sum(G4double k)
293{
294 G4double totalCrossSection = 0.;
295
296 for (G4int i=0; i<nLevels; i++)
297 {
298 totalCrossSection += PartialCrossSection(k,i);
299 }
300 return totalCrossSection;
301}
302
303//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
304
305G4double G4DNASancheExcitationModel::LinInterpolate(G4double e1,
306 G4double e2,
307 G4double e,
308 G4double xs1,
309 G4double xs2)
310{
311 G4double a = (xs2 - xs1) / (e2 - e1);
312 G4double b = xs2 - a*e2;
313 G4double value = a*e + b;
314 // G4cout<<"interP >> "<<e1<<" "<<e2<<" "<<e<<" "<<xs1<<" "<<xs2<<" "<<a<<" "<<b<<" "<<value<<G4endl;
315
316 return value;
317}
318
Note: See TracBrowser for help on using the repository browser.