source: trunk/source/processes/electromagnetic/lowenergy/src/G4DNABornExcitationModel.cc@ 1346

Last change on this file since 1346 was 1340, checked in by garnier, 15 years ago

update ti head
File size: 9.6 KB
RevLine 
[1058]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//
[1340]26// $Id: G4DNABornExcitationModel.cc,v 1.10 2010/08/24 13:51:06 sincerti Exp $
27// GEANT4 tag $Name: emlowen-V09-03-54 $
[1058]28//
29
30#include "G4DNABornExcitationModel.hh"
31
32//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
33
34using namespace std;
35
36//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
37
38G4DNABornExcitationModel::G4DNABornExcitationModel(const G4ParticleDefinition*,
39 const G4String& nam)
40:G4VEmModel(nam),isInitialised(false)
41{
42 verboseLevel= 0;
43 // Verbosity scale:
44 // 0 = nothing
45 // 1 = warning for energy non-conservation
46 // 2 = details of energy budget
47 // 3 = calculation of cross sections, file openings, sampling of atoms
48 // 4 = entering in methods
49
[1192]50 if( verboseLevel>0 )
51 {
[1315]52 G4cout << "Born excitation model is constructed " << G4endl;
[1192]53 }
54
[1058]55}
56
57//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
58
59G4DNABornExcitationModel::~G4DNABornExcitationModel()
60{
61 // Cross section
[1315]62
63 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
64 for (pos = tableData.begin(); pos != tableData.end(); ++pos)
65 {
66 G4DNACrossSectionDataSet* table = pos->second;
67 delete table;
68 }
69
[1058]70}
71
72//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
73
[1315]74void G4DNABornExcitationModel::Initialise(const G4ParticleDefinition* particle,
[1058]75 const G4DataVector& /*cuts*/)
76{
77
78 if (verboseLevel > 3)
79 G4cout << "Calling G4DNABornExcitationModel::Initialise()" << G4endl;
80
[1315]81 G4String fileElectron("dna/sigma_excitation_e_born");
82 G4String fileProton("dna/sigma_excitation_p_born");
83
84 G4ParticleDefinition* electronDef = G4Electron::ElectronDefinition();
85 G4ParticleDefinition* protonDef = G4Proton::ProtonDefinition();
86
87 G4String electron;
88 G4String proton;
[1058]89
[1315]90 G4double scaleFactor = (1.e-22 / 3.343) * m*m;
91
92 if (electronDef != 0)
[1058]93 {
[1315]94 electron = electronDef->GetParticleName();
95
96 tableFile[electron] = fileElectron;
97
98 lowEnergyLimit[electron] = 9. * eV;
99 highEnergyLimit[electron] = 1. * MeV;
100
101 // Cross section
102
103 G4DNACrossSectionDataSet* tableE = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, eV,scaleFactor );
104 tableE->LoadData(fileElectron);
105
106 tableData[electron] = tableE;
107
[1058]108 }
[1315]109 else
110 {
111 G4Exception("G4DNABornExcitationModel::Initialise(): electron is not defined");
112 }
[1058]113
[1315]114 if (protonDef != 0)
[1058]115 {
[1315]116 proton = protonDef->GetParticleName();
117
118 tableFile[proton] = fileProton;
119
120 lowEnergyLimit[proton] = 500. * keV;
121 highEnergyLimit[proton] = 100. * MeV;
122
123 // Cross section
124
125 G4DNACrossSectionDataSet* tableP = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, eV,scaleFactor );
126 tableP->LoadData(fileProton);
127
128 tableData[proton] = tableP;
129
[1058]130 }
[1315]131 else
132 {
133 G4Exception("G4DNABornExcitationModel::Initialise(): proton is not defined");
134 }
[1058]135
[1315]136 if (particle==electronDef)
[1058]137 {
[1315]138 SetLowEnergyLimit(lowEnergyLimit[electron]);
139 SetHighEnergyLimit(highEnergyLimit[electron]);
[1058]140 }
141
[1315]142 if (particle==protonDef)
143 {
144 SetLowEnergyLimit(lowEnergyLimit[proton]);
145 SetHighEnergyLimit(highEnergyLimit[proton]);
146 }
147
[1192]148 if( verboseLevel>0 )
149 {
150 G4cout << "Born excitation model is initialized " << G4endl
151 << "Energy range: "
[1315]152 << LowEnergyLimit() / eV << " eV - "
153 << HighEnergyLimit() / keV << " keV for "
154 << particle->GetParticleName()
155 << G4endl;
[1192]156 }
157
[1315]158
[1058]159 if(!isInitialised)
160 {
161 isInitialised = true;
162
163 if(pParticleChange)
164 fParticleChangeForGamma = reinterpret_cast<G4ParticleChangeForGamma*>(pParticleChange);
165 else
166 fParticleChangeForGamma = new G4ParticleChangeForGamma();
167 }
168
169 // InitialiseElementSelectors(particle,cuts);
170
171}
172
173//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
174
[1315]175G4double G4DNABornExcitationModel::CrossSectionPerVolume(const G4Material* material,
[1058]176 const G4ParticleDefinition* particleDefinition,
[1315]177 G4double ekin,
[1058]178 G4double,
179 G4double)
180{
181 if (verboseLevel > 3)
182 G4cout << "Calling CrossSectionPerVolume() of G4DNABornExcitationModel" << G4endl;
183
[1315]184 if (
185 particleDefinition != G4Proton::ProtonDefinition()
186 &&
187 particleDefinition != G4Electron::ElectronDefinition()
188 )
189
190 return 0;
191
[1058]192 // Calculate total cross section for model
193
[1315]194 G4double lowLim = 0;
195 G4double highLim = 0;
196 G4double sigma=0;
197
198 if (material->GetName() == "G4_WATER")
[1058]199 {
[1315]200 const G4String& particleName = particleDefinition->GetParticleName();
201
202 std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
203 pos1 = lowEnergyLimit.find(particleName);
204 if (pos1 != lowEnergyLimit.end())
[1058]205 {
[1315]206 lowLim = pos1->second;
207 }
208
209 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
210 pos2 = highEnergyLimit.find(particleName);
211 if (pos2 != highEnergyLimit.end())
212 {
213 highLim = pos2->second;
214 }
215
216 if (ekin >= lowLim && ekin < highLim)
217 {
218 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
219 pos = tableData.find(particleName);
220
221 if (pos != tableData.end())
[1058]222 {
[1315]223 G4DNACrossSectionDataSet* table = pos->second;
224 if (table != 0)
225 {
226 sigma = table->FindValue(ekin);
227 }
[1058]228 }
[1315]229 else
[1058]230 {
[1315]231 G4Exception("G4DNABornExcitationModel::CrossSectionPerVolume: attempting to calculate cross section for wrong particle");
232 }
[1058]233 }
234
[1315]235 if (verboseLevel > 3)
236 {
237 G4cout << "---> Kinetic energy(eV)=" << ekin/eV << G4endl;
238 G4cout << " - Cross section per water molecule (cm^2)=" << sigma/cm/cm << G4endl;
239 G4cout << " - Cross section per water molecule (cm^-1)=" << sigma*material->GetAtomicNumDensityVector()[1]/(1./cm) << G4endl;
240 }
241
242 } // if (waterMaterial)
243
244 return sigma*material->GetAtomicNumDensityVector()[1];
245
[1058]246
247}
248
249//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
250
251void G4DNABornExcitationModel::SampleSecondaries(std::vector<G4DynamicParticle*>* /*fvect*/,
252 const G4MaterialCutsCouple* /*couple*/,
253 const G4DynamicParticle* aDynamicParticle,
254 G4double,
255 G4double)
256{
257
258 if (verboseLevel > 3)
259 G4cout << "Calling SampleSecondaries() of G4DNABornExcitationModel" << G4endl;
260
261 G4double k = aDynamicParticle->GetKineticEnergy();
262
[1315]263 const G4String& particleName = aDynamicParticle->GetDefinition()->GetParticleName();
264
265 G4int level = RandomSelect(k,particleName);
[1058]266 G4double excitationEnergy = waterStructure.ExcitationEnergy(level);
267 G4double newEnergy = k - excitationEnergy;
268
269 if (newEnergy > 0)
270 {
271 fParticleChangeForGamma->ProposeMomentumDirection(aDynamicParticle->GetMomentumDirection());
272 fParticleChangeForGamma->SetProposedKineticEnergy(newEnergy);
273 fParticleChangeForGamma->ProposeLocalEnergyDeposit(excitationEnergy);
274 }
275
276}
277
278//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
279
[1315]280G4int G4DNABornExcitationModel::RandomSelect(G4double k, const G4String& particle)
[1058]281{
282 G4int level = 0;
283
[1315]284 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
285 pos = tableData.find(particle);
[1058]286
[1315]287 if (pos != tableData.end())
[1058]288 {
[1315]289 G4DNACrossSectionDataSet* table = pos->second;
290
291 if (table != 0)
[1058]292 {
[1315]293 G4double* valuesBuffer = new G4double[table->NumberOfComponents()];
294 const size_t n(table->NumberOfComponents());
295 size_t i(n);
296 G4double value = 0.;
297
298 while (i>0)
299 {
300 i--;
301 valuesBuffer[i] = table->GetComponent(i)->FindValue(k);
302 value += valuesBuffer[i];
303 }
304
305 value *= G4UniformRand();
306
307 i = n;
308
309 while (i > 0)
310 {
311 i--;
312
313 if (valuesBuffer[i] > value)
314 {
315 delete[] valuesBuffer;
316 return i;
317 }
318 value -= valuesBuffer[i];
319 }
320
321 if (valuesBuffer) delete[] valuesBuffer;
322
[1058]323 }
324 }
[1315]325 else
326 {
327 G4Exception("G4DNABornExcitationModel::RandomSelect attempting to calculate cross section for wrong particle");
328 }
[1058]329 return level;
330}
331
332
Note: See TracBrowser for help on using the repository browser.