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

Last change on this file since 1229 was 1228, checked in by garnier, 16 years ago

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