| Geant4 Cross Reference |
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 // Based on the work described in
27 // Rad Res 163, 98-111 (2005)
28 // D. Emfietzoglou, H. Nikjoo
29 //
30 // Authors of the class (2014):
31 // I. Kyriakou (kyriak@cc.uoi.gr)
32 // D. Emfietzoglou (demfietz@cc.uoi.gr)
33 // S. Incerti (incerti@cenbg.in2p3.fr)
34 //
35
36 #include "G4DNAEmfietzoglouExcitationModel.hh"
37 #include "G4SystemOfUnits.hh"
38 #include "G4DNAChemistryManager.hh"
39 #include "G4DNAMolecularMaterial.hh"
40
41 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
42
43 using namespace std;
44
45 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
46
47 G4DNAEmfietzoglouExcitationModel::G4DNAEmfietzoglouExcitationModel(const G4ParticleDefinition*,
48 const G4String& nam)
49 :G4VEmModel(nam)
50 {
51 fpMolWaterDensity = nullptr;
52
53 verboseLevel= 0;
54 // Verbosity scale:
55 // 0 = nothing
56 // 1 = warning for energy non-conservation
57 // 2 = details of energy budget
58 // 3 = calculation of cross sections, file openings, sampling of atoms
59 // 4 = entering in methods
60
61 if( verboseLevel>0 )
62 {
63 G4cout << "Emfietzoglou excitation model is constructed " << G4endl;
64 }
65 fParticleChangeForGamma = nullptr;
66
67 SetLowEnergyLimit(8.*eV);
68 SetHighEnergyLimit(10.*keV);
69
70 // Selection of stationary mode
71 statCode = false;
72 }
73
74 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
75
76 G4DNAEmfietzoglouExcitationModel::~G4DNAEmfietzoglouExcitationModel()
77 {
78 // Cross section
79
80 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
81 for (pos = tableData.begin(); pos != tableData.end(); ++pos)
82 {
83 G4DNACrossSectionDataSet* table = pos->second;
84 delete table;
85 }
86
87 }
88
89 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
90
91 void G4DNAEmfietzoglouExcitationModel::Initialise(const G4ParticleDefinition* particle,
92 const G4DataVector& /*cuts*/)
93 {
94
95 if (verboseLevel > 3)
96 G4cout << "Calling G4DNAEmfietzoglouExcitationModel::Initialise()" << G4endl;
97
98 G4String fileElectron("dna/sigma_excitation_e_emfietzoglou");
99
100 G4ParticleDefinition* electronDef = G4Electron::ElectronDefinition();
101
102 G4String electron;
103
104 G4double scaleFactor = (1.e-22 / 3.343) * m*m;
105
106 // *** ELECTRON
107
108 electron = electronDef->GetParticleName();
109
110 tableFile[electron] = fileElectron;
111
112 // Cross section
113
114 auto tableE = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, eV,scaleFactor );
115 tableE->LoadData(fileElectron);
116
117 tableData[electron] = tableE;
118
119 //
120
121 if( verboseLevel>0 )
122 {
123 G4cout << "Emfietzoglou excitation model is initialized " << G4endl
124 << "Energy range: "
125 << LowEnergyLimit() / eV << " eV - "
126 << HighEnergyLimit() / keV << " keV for "
127 << particle->GetParticleName()
128 << G4endl;
129 }
130
131 // Initialize water density pointer
132 fpMolWaterDensity = G4DNAMolecularMaterial::Instance()->GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER"));
133
134 if (isInitialised) return;
135 fParticleChangeForGamma = GetParticleChangeForGamma();
136 isInitialised = true;
137 }
138
139 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
140
141 G4double G4DNAEmfietzoglouExcitationModel::CrossSectionPerVolume(const G4Material* material,
142 const G4ParticleDefinition* particleDefinition,
143 G4double ekin,
144 G4double,
145 G4double)
146 {
147 if (verboseLevel > 3)
148 G4cout << "Calling CrossSectionPerVolume() of G4DNAEmfietzoglouExcitationModel" << G4endl;
149
150 if (particleDefinition != G4Electron::ElectronDefinition()) return 0;
151
152 // Calculate total cross section for model
153
154 G4double sigma=0;
155
156 G4double waterDensity = (*fpMolWaterDensity)[material->GetIndex()];
157
158 const G4String& particleName = particleDefinition->GetParticleName();
159
160 if (ekin >= LowEnergyLimit() && ekin <= HighEnergyLimit())
161 {
162 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
163 pos = tableData.find(particleName);
164
165 if (pos != tableData.end())
166 {
167 G4DNACrossSectionDataSet* table = pos->second;
168 if (table != nullptr) sigma = table->FindValue(ekin);
169 }
170 else
171 {
172 G4Exception("G4DNAEmfietzoglouExcitationModel::CrossSectionPerVolume","em0002",
173 FatalException,"Model not applicable to particle type.");
174 }
175 }
176
177 if (verboseLevel > 2)
178 {
179 G4cout << "__________________________________" << G4endl;
180 G4cout << "G4DNAEmfietzoglouExcitationModel - XS INFO START" << G4endl;
181 G4cout << "Kinetic energy(eV)=" << ekin/eV << " particle : " << particleName << G4endl;
182 G4cout << "Cross section per water molecule (cm^2)=" << sigma/cm/cm << G4endl;
183 G4cout << "Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./cm) << G4endl;
184 //G4cout << " Cross section per water molecule (cm^-1)=" <<
185 ///sigma*material->GetAtomicNumDensityVector()[1]/(1./cm) << G4endl;
186 G4cout << "G4DNAEmfietzoglouExcitationModel - XS INFO END" << G4endl;
187 }
188
189 return sigma*waterDensity;
190 }
191
192 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
193
194 void G4DNAEmfietzoglouExcitationModel::SampleSecondaries(std::vector<G4DynamicParticle*>* /*fvect*/,
195 const G4MaterialCutsCouple* /*couple*/,
196 const G4DynamicParticle* aDynamicParticle,
197 G4double,
198 G4double)
199 {
200
201 if (verboseLevel > 3)
202 G4cout << "Calling SampleSecondaries() of G4DNAEmfietzoglouExcitationModel" << G4endl;
203
204 G4double k = aDynamicParticle->GetKineticEnergy();
205
206 const G4String& particleName = aDynamicParticle->GetDefinition()->GetParticleName();
207
208 G4int level = RandomSelect(k,particleName);
209 G4double excitationEnergy = waterStructure.ExcitationEnergy(level);
210 G4double newEnergy = k - excitationEnergy;
211
212 if (newEnergy > 0)
213 {
214 fParticleChangeForGamma->ProposeMomentumDirection(aDynamicParticle->GetMomentumDirection());
215
216 if (!statCode) fParticleChangeForGamma->SetProposedKineticEnergy(newEnergy);
217 else fParticleChangeForGamma->SetProposedKineticEnergy(k);
218
219 fParticleChangeForGamma->ProposeLocalEnergyDeposit(excitationEnergy);
220 }
221
222 const G4Track * theIncomingTrack = fParticleChangeForGamma->GetCurrentTrack();
223 G4DNAChemistryManager::Instance()->CreateWaterMolecule(eExcitedMolecule,
224 level,
225 theIncomingTrack);
226 }
227
228 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
229
230 G4int G4DNAEmfietzoglouExcitationModel::RandomSelect(G4double k, const G4String& particle)
231 {
232
233 G4int level = 0;
234
235 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
236 pos = tableData.find(particle);
237
238 if (pos != tableData.end())
239 {
240 G4DNACrossSectionDataSet* table = pos->second;
241
242 if (table != nullptr)
243 {
244 auto valuesBuffer = new G4double[table->NumberOfComponents()];
245 const auto n = (G4int)table->NumberOfComponents();
246 G4int i(n);
247 G4double value = 0.;
248
249 //Check reading of initial xs file
250 //G4cout << table->GetComponent(0)->FindValue(k)/ ((1.e-22 / 3.343) * m*m) << G4endl;
251 //G4cout << table->GetComponent(1)->FindValue(k)/ ((1.e-22 / 3.343) * m*m) << G4endl;
252 //G4cout << table->GetComponent(2)->FindValue(k)/ ((1.e-22 / 3.343) * m*m) << G4endl;
253 //G4cout << table->GetComponent(3)->FindValue(k)/ ((1.e-22 / 3.343) * m*m) << G4endl;
254 //G4cout << table->GetComponent(4)->FindValue(k)/ ((1.e-22 / 3.343) * m*m) << G4endl;
255 //G4cout << table->GetComponent(5)->FindValue(k)/ ((1.e-22 / 3.343) * m*m) << G4endl;
256 //G4cout << table->GetComponent(6)->FindValue(k)/ ((1.e-22 / 3.343) * m*m) << G4endl;
257 //abort();
258
259 while (i>0)
260 {
261 i--;
262 valuesBuffer[i] = table->GetComponent(i)->FindValue(k);
263 value += valuesBuffer[i];
264 }
265
266 value *= G4UniformRand();
267
268 i = n;
269
270 while (i > 0)
271 {
272 i--;
273
274 if (valuesBuffer[i] > value)
275 {
276 delete[] valuesBuffer;
277 return i;
278 }
279 value -= valuesBuffer[i];
280 }
281
282 delete[] valuesBuffer;
283
284 }
285 }
286 else
287 {
288 G4Exception("G4DNAEmfietzoglouExcitationModel::RandomSelect","em0002",
289 FatalException,"Model not applicable to particle type.");
290 }
291 return level;
292 }
293