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

Last change on this file since 1066 was 1058, checked in by garnier, 17 years ago

file release beta
File size: 8.2 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: G4DNABornExcitationModel.cc,v 1.3 2009/02/16 11:00:11 sincerti Exp $
27// GEANT4 tag $Name: geant4-09-03-beta-cand-01 $
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
43 lowEnergyLimit = 500 * keV;
44 highEnergyLimit = 10 * MeV;
45 SetLowEnergyLimit(lowEnergyLimit);
46 SetHighEnergyLimit(highEnergyLimit);
47
48 verboseLevel= 0;
49 // Verbosity scale:
50 // 0 = nothing
51 // 1 = warning for energy non-conservation
52 // 2 = details of energy budget
53 // 3 = calculation of cross sections, file openings, sampling of atoms
54 // 4 = entering in methods
55
56 //
57
58 table = 0;
59
60 //
61
62 G4cout << "Born excitation model is constructed " << G4endl
63 << "Energy range: "
64 << lowEnergyLimit / keV << " keV - "
65 << highEnergyLimit / MeV << " MeV"
66 << G4endl;
67
68}
69
70//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
71
72G4DNABornExcitationModel::~G4DNABornExcitationModel()
73{
74 // Cross section
75 delete table;
76}
77
78//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
79
80void G4DNABornExcitationModel::Initialise(const G4ParticleDefinition* /*particle*/,
81 const G4DataVector& /*cuts*/)
82{
83
84 if (verboseLevel > 3)
85 G4cout << "Calling G4DNABornExcitationModel::Initialise()" << G4endl;
86
87 // Energy limits
88
89 if (LowEnergyLimit() < lowEnergyLimit)
90 {
91 G4cout << "G4DNABornExcitationModel: low energy limit increased from " <<
92 LowEnergyLimit()/keV << " keV to " << lowEnergyLimit/keV << " keV" << G4endl;
93 SetLowEnergyLimit(lowEnergyLimit);
94 }
95
96 if (HighEnergyLimit() > highEnergyLimit)
97 {
98 G4cout << "G4DNABornExcitationModel: high energy limit decreased from " <<
99 HighEnergyLimit()/MeV << " MeV to " << highEnergyLimit/MeV << " MeV" << G4endl;
100 SetHighEnergyLimit(highEnergyLimit);
101 }
102
103 //
104
105 if (table == 0)
106 {
107 table = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, eV,(1e-22/3.343)*m*m );
108 table->LoadData("dna/sigma_excitation_p_born");
109 }
110
111 G4cout << "Born excitation model is initialized " << G4endl
112 << "Energy range: "
113 << LowEnergyLimit() / keV << " keV - "
114 << HighEnergyLimit() / MeV << " MeV " << G4endl;
115
116 if(!isInitialised)
117 {
118 isInitialised = true;
119
120 if(pParticleChange)
121 fParticleChangeForGamma = reinterpret_cast<G4ParticleChangeForGamma*>(pParticleChange);
122 else
123 fParticleChangeForGamma = new G4ParticleChangeForGamma();
124 }
125
126 // InitialiseElementSelectors(particle,cuts);
127
128 // Test if water material
129
130 flagMaterialIsWater= false;
131 densityWater = 0;
132
133 const G4ProductionCutsTable* theCoupleTable = G4ProductionCutsTable::GetProductionCutsTable();
134
135 if(theCoupleTable)
136 {
137 G4int numOfCouples = theCoupleTable->GetTableSize();
138
139 if(numOfCouples>0)
140 {
141 for (G4int i=0; i<numOfCouples; i++)
142 {
143 const G4MaterialCutsCouple* couple = theCoupleTable->GetMaterialCutsCouple(i);
144 const G4Material* material = couple->GetMaterial();
145
146 size_t j = material->GetNumberOfElements();
147 while (j>0)
148 {
149 j--;
150 const G4Element* element(material->GetElement(j));
151 if (element->GetZ() == 8.)
152 {
153 G4double density = material->GetAtomicNumDensityVector()[j];
154 if (density > 0.)
155 {
156 flagMaterialIsWater = true;
157 densityWater = density;
158
159 if (verboseLevel > 3)
160 G4cout << "Water material is found with density(cm^-3)=" << density/(cm*cm*cm) << G4endl;
161 }
162 }
163 }
164
165 }
166 } // if(numOfCouples>0)
167
168 } // if (theCoupleTable)
169
170}
171
172//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
173
174G4double G4DNABornExcitationModel::CrossSectionPerVolume(const G4Material*,
175 const G4ParticleDefinition* particleDefinition,
176 G4double k,
177 G4double,
178 G4double)
179{
180 if (verboseLevel > 3)
181 G4cout << "Calling CrossSectionPerVolume() of G4DNABornExcitationModel" << G4endl;
182
183 // Calculate total cross section for model
184
185 G4double crossSection=0;
186
187 if (flagMaterialIsWater)
188 {
189 if (particleDefinition == G4Proton::ProtonDefinition())
190 {
191 if (k >= lowEnergyLimit && k < highEnergyLimit)
192 {
193 crossSection = table->FindValue(k);
194 }
195
196 if (verboseLevel > 3)
197 {
198 G4cout << "---> Kinetic energy(keV)=" << k/keV << G4endl;
199 G4cout << " - Cross section per water molecule (cm^2)=" << crossSection/cm/cm << G4endl;
200 G4cout << " - Cross section per water molecule (cm^-1)=" << crossSection*densityWater/(1./cm) << G4endl;
201 }
202 }
203 } // if (flagMaterialIsWater)
204
205 return crossSection*densityWater;
206
207}
208
209//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
210
211void G4DNABornExcitationModel::SampleSecondaries(std::vector<G4DynamicParticle*>* /*fvect*/,
212 const G4MaterialCutsCouple* /*couple*/,
213 const G4DynamicParticle* aDynamicParticle,
214 G4double,
215 G4double)
216{
217
218 if (verboseLevel > 3)
219 G4cout << "Calling SampleSecondaries() of G4DNABornExcitationModel" << G4endl;
220
221 G4double k = aDynamicParticle->GetKineticEnergy();
222
223 G4int level = RandomSelect(k);
224 G4double excitationEnergy = waterStructure.ExcitationEnergy(level);
225 G4double newEnergy = k - excitationEnergy;
226
227 if (newEnergy > 0)
228 {
229 fParticleChangeForGamma->ProposeMomentumDirection(aDynamicParticle->GetMomentumDirection());
230 fParticleChangeForGamma->SetProposedKineticEnergy(newEnergy);
231 fParticleChangeForGamma->ProposeLocalEnergyDeposit(excitationEnergy);
232 }
233
234}
235
236//....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
237
238G4int G4DNABornExcitationModel::RandomSelect(G4double k)
239{
240 G4int level = 0;
241
242 G4double* valuesBuffer = new G4double[table->NumberOfComponents()];
243
244 const size_t n(table->NumberOfComponents());
245 size_t i(n);
246 G4double value = 0.;
247
248 while (i>0)
249 {
250 i--;
251 valuesBuffer[i] = table->GetComponent(i)->FindValue(k);
252 value += valuesBuffer[i];
253 }
254
255 value *= G4UniformRand();
256
257 i = n;
258
259 while (i > 0)
260 {
261 i--;
262
263 if (valuesBuffer[i] > value)
264 {
265 delete[] valuesBuffer;
266 return i;
267 }
268 value -= valuesBuffer[i];
269 }
270
271 if (valuesBuffer) delete[] valuesBuffer;
272
273 return level;
274}
275
276
Note: See TracBrowser for help on using the repository browser.