| Geant4 Cross Reference |
1 // 1
2 // *******************************************
3 // * License and Disclaimer
4 // *
5 // * The Geant4 software is copyright of th
6 // * the Geant4 Collaboration. It is provided
7 // * conditions of the Geant4 Software License
8 // * LICENSE and available at http://cern.ch/
9 // * include a list of copyright holders.
10 // *
11 // * Neither the authors of this software syst
12 // * institutes,nor the agencies providing fin
13 // * work make any representation or warran
14 // * regarding this software system or assum
15 // * use. Please see the license in the file
16 // * for the full disclaimer and the limitatio
17 // *
18 // * This code implementation is the result
19 // * technical work of the GEANT4 collaboratio
20 // * By using, copying, modifying or distri
21 // * any work based on the software) you ag
22 // * use in resulting scientific publicati
23 // * acceptance of all terms of the Geant4 Sof
24 // *******************************************
25 //
26 //
27
28 #include "G4DNABornExcitationModel2.hh"
29 #include "G4SystemOfUnits.hh"
30 #include "G4DNAChemistryManager.hh"
31 #include "G4DNAMolecularMaterial.hh"
32 #include "G4PhysicsTable.hh"
33 #include "G4PhysicsVector.hh"
34 #include "G4UnitsTable.hh"
35 #include <map>
36
37 //....oooOO0OOooo........oooOO0OOooo........oo
38
39 using namespace std;
40
41 //....oooOO0OOooo........oooOO0OOooo........oo
42
43 G4DNABornExcitationModel2::G4DNABornExcitation
44
45 G4VEmModel(nam)
46 {
47 fpMolWaterDensity = nullptr;
48 fHighEnergy = 0;
49 fLowEnergy = 0;
50 fParticleDefinition = nullptr;
51
52 verboseLevel = 0;
53 // Verbosity scale:
54 // 0 = nothing
55 // 1 = warning for energy non-conservation
56 // 2 = details of energy budget
57 // 3 = calculation of cross sections, file o
58 // 4 = entering in methods
59
60 if (verboseLevel > 0)
61 {
62 G4cout << "Born excitation model is constr
63 }
64 fParticleChangeForGamma = nullptr;
65 fLastBinCallForFinalXS = 0;
66 fTotalXS = nullptr;
67 fTableData = nullptr;
68
69 // Selection of stationary mode
70
71 statCode = false;
72 }
73
74 //....oooOO0OOooo........oooOO0OOooo........oo
75
76 G4DNABornExcitationModel2::~G4DNABornExcitatio
77 {
78 // Cross section
79
80 delete fTableData;
81 }
82
83 //....oooOO0OOooo........oooOO0OOooo........oo
84
85 void G4DNABornExcitationModel2::Initialise(con
86 con
87 {
88
89 if (verboseLevel > 3)
90 {
91 G4cout << "Calling G4DNABornExcitationMode
92 }
93
94 if(fParticleDefinition != nullptr && fPartic
95 {
96 G4Exception("G4DNABornExcitationModel2::In
97 FatalException,"Model already initiali
98 }
99
100 fParticleDefinition = particle;
101
102 std::ostringstream fullFileName;
103 const char* path = G4FindDataDir("G4LEDATA")
104
105 if(G4String(path).empty())
106 {
107 G4Exception("G4DNABornExcitationModel2::In
108 FatalException, "G4LEDATA not defined
109 }
110
111 fullFileName << path;
112
113 if(particle->GetParticleName() == "e-")
114 {
115 fullFileName << "/dna/bornExcitation-e.dat
116 fLowEnergy = 9*eV;
117 fHighEnergy = 1*MeV;
118 }
119 else if(particle->GetParticleName() == "prot
120 {
121 fullFileName << "/dna/bornExcitation-p.dat
122 fLowEnergy = 500. * keV;
123 fHighEnergy = 100. * MeV;
124 }
125
126 SetLowEnergyLimit(fLowEnergy);
127 SetHighEnergyLimit(fHighEnergy);
128
129 //G4double scaleFactor = (1.e-22 / 3.343) *
130
131 fTableData = new G4PhysicsTable();
132 fTableData->RetrievePhysicsTable(fullFileNam
133 /*
134 for(std::size_t level = 0; level<fTableData-
135 {
136 //(*fTableData)(level)->ScaleVector(1,scal
137 }
138 */
139 std::size_t finalBin_i = 2000;
140 G4double E_min = fLowEnergy;
141 G4double E_max = fHighEnergy;
142 fTotalXS = new G4PhysicsLogVector(E_min, E_m
143 G4double energy;
144 G4double finalXS;
145
146 for(std::size_t energy_i = 0; energy_i < fin
147 {
148 energy = fTotalXS->Energy(energy_i);
149 finalXS = 0;
150
151 for(std::size_t level = 0; level<fTableDat
152 {
153 finalXS += (*fTableData)(level)->Value(e
154 }
155 fTotalXS->PutValue(energy_i, finalXS);
156 //G4cout << "energy = " << energy << " " <
157 // << " " << energy_i << " " << fina
158 }
159
160 // for(energy = LowEnergyLimit() ; energy
161 // {
162 // G4cout << "energy = " << energy << "
163 // }
164
165 if( verboseLevel>0 )
166 {
167 G4cout << "Born excitation model is initia
168 << "Energy range: "
169 << LowEnergyLimit() / eV << " eV - "
170 << HighEnergyLimit() / keV << " keV for "
171 << particle->GetParticleName()
172 << G4endl;
173 }
174
175 // Initialize water density pointer
176 fpMolWaterDensity = G4DNAMolecularMaterial::
177
178 if (isInitialised)
179 { return;}
180 fParticleChangeForGamma = GetParticleChangeF
181 isInitialised = true;
182 }
183
184 //....oooOO0OOooo........oooOO0OOooo........oo
185
186 G4double G4DNABornExcitationModel2::CrossSecti
187
188
189
190
191 {
192 if (verboseLevel > 3)
193 {
194 G4cout << "Calling CrossSectionPerVolume()
195 << G4endl;
196 }
197
198 if(particleDefinition != fParticleDefinition
199
200 // Calculate total cross section for model
201
202 G4double sigma=0;
203
204 G4double waterDensity = (*fpMolWaterDensity)
205
206 if (ekin >= fLowEnergy && ekin <= fHighEnerg
207 {
208 sigma = fTotalXS->Value(ekin, fLastBinCall
209
210 // for(std::size_t i = 0; i < 5; ++i)
211 // sigma += (*fTableData)[i]->Value(ekin);
212
213 if(sigma == 0)
214 {
215 G4cerr << "PROBLEM SIGMA = 0 at " << G4B
216 }
217 }
218
219 if (verboseLevel > 2)
220 {
221 G4cout << "_______________________________
222 G4cout << "G4DNABornExcitationModel2 - XS
223 G4cout << "Kinetic energy(eV)=" << ekin/eV
224 G4cout << "Cross section per water molecul
225 G4cout << "Cross section per water molecul
226 G4cout << "G4DNABornExcitationModel2 - XS
227 }
228
229 return sigma*waterDensity;
230 }
231
232 //....oooOO0OOooo........oooOO0OOooo........oo
233
234 void G4DNABornExcitationModel2::SampleSecondar
235
236
237
238
239 {
240
241 if (verboseLevel > 3)
242 {
243 G4cout << "Calling SampleSecondaries() of
244 << G4endl;
245 }
246
247 G4double k = aDynamicParticle->GetKineticEne
248
249 G4int level = RandomSelect(k);
250 G4double excitationEnergy = waterStructure.E
251 G4double newEnergy = k - excitationEnergy;
252
253 if (newEnergy > 0)
254 {
255 fParticleChangeForGamma->ProposeMomentumDi
256
257 if (!statCode) fParticleChangeForGamma->Se
258 else fParticleChangeForGamma->SetProposedK
259
260 fParticleChangeForGamma->ProposeLocalEnerg
261 }
262
263 const G4Track * theIncomingTrack = fParticle
264 G4DNAChemistryManager::Instance()->CreateWat
265 level,
266 theIncomingTrack);
267 }
268
269 //....oooOO0OOooo........oooOO0OOooo........oo
270
271 G4double G4DNABornExcitationModel2::GetPartial
272
273
274
275 {
276 if (fParticleDefinition != particle)
277 {
278 G4Exception("G4DNABornExcitationModel2::Ge
279 "bornParticleType",
280 FatalException,
281 "Model initialized for another
282 }
283
284 return (*fTableData)(level)->Value(kineticEn
285 }
286
287 //....oooOO0OOooo........oooOO0OOooo........oo
288
289 G4int G4DNABornExcitationModel2::RandomSelect(
290 {
291 const std::size_t n(fTableData->size());
292 std::size_t i(n);
293
294 G4double value = fTotalXS->Value(k, fLastBin
295
296 value *= G4UniformRand();
297 i = n;
298
299 G4double partialXS;
300
301 while (i > 0)
302 {
303 i--;
304
305 partialXS = (*fTableData)(i)->Value(k);
306 if (partialXS > value)
307 {
308 return (G4int)i;
309 }
310 value -= partialXS;
311 }
312
313 return 0;
314 }
315
316