| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // Based on the work described in 26 // Based on the work described in
27 // Rad Res 163, 98-111 (2005) 27 // Rad Res 163, 98-111 (2005)
28 // D. Emfietzoglou, H. Nikjoo 28 // D. Emfietzoglou, H. Nikjoo
29 // << 29 //
30 // Authors of the class (2014): 30 // Authors of the class (2014):
31 // I. Kyriakou (kyriak@cc.uoi.gr) 31 // I. Kyriakou (kyriak@cc.uoi.gr)
32 // D. Emfietzoglou (demfietz@cc.uoi.gr) 32 // D. Emfietzoglou (demfietz@cc.uoi.gr)
33 // S. Incerti (incerti@cenbg.in2p3.fr) 33 // S. Incerti (incerti@cenbg.in2p3.fr)
34 // 34 //
35 35
36 #include "G4DNAEmfietzoglouExcitationModel.hh" 36 #include "G4DNAEmfietzoglouExcitationModel.hh"
37 #include "G4SystemOfUnits.hh" 37 #include "G4SystemOfUnits.hh"
38 #include "G4DNAChemistryManager.hh" 38 #include "G4DNAChemistryManager.hh"
39 #include "G4DNAMolecularMaterial.hh" 39 #include "G4DNAMolecularMaterial.hh"
40 40
41 //....oooOO0OOooo........oooOO0OOooo........oo 41 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
42 42
43 using namespace std; 43 using namespace std;
44 44
45 //....oooOO0OOooo........oooOO0OOooo........oo 45 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
46 46
47 G4DNAEmfietzoglouExcitationModel::G4DNAEmfietz 47 G4DNAEmfietzoglouExcitationModel::G4DNAEmfietzoglouExcitationModel(const G4ParticleDefinition*,
48 48 const G4String& nam)
49 :G4VEmModel(nam) << 49 :G4VEmModel(nam),isInitialised(false)
50 { 50 {
51 fpMolWaterDensity = nullptr; << 51
>> 52 // nistwater = G4NistManager::Instance()->FindOrBuildMaterial("G4_WATER");
>> 53 fpMolWaterDensity = 0;
52 54
53 verboseLevel= 0; 55 verboseLevel= 0;
54 // Verbosity scale: 56 // Verbosity scale:
55 // 0 = nothing 57 // 0 = nothing
56 // 1 = warning for energy non-conservation 58 // 1 = warning for energy non-conservation
57 // 2 = details of energy budget 59 // 2 = details of energy budget
58 // 3 = calculation of cross sections, file 60 // 3 = calculation of cross sections, file openings, sampling of atoms
59 // 4 = entering in methods 61 // 4 = entering in methods
60 62
61 if( verboseLevel>0 ) 63 if( verboseLevel>0 )
62 { 64 {
63 G4cout << "Emfietzoglou excitation model << 65 G4cout << "Emfietzoglou excitation model is constructed " << G4endl;
64 } 66 }
65 fParticleChangeForGamma = nullptr; << 67 fParticleChangeForGamma = 0;
66 << 68 SetLowEnergyLimit(8. * eV);
67 SetLowEnergyLimit(8.*eV); << 69 SetHighEnergyLimit(10. * keV);
68 SetHighEnergyLimit(10.*keV); <<
69 70
70 // Selection of stationary mode 71 // Selection of stationary mode
>> 72
71 statCode = false; 73 statCode = false;
72 } 74 }
73 75
74 //....oooOO0OOooo........oooOO0OOooo........oo 76 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
75 77
76 G4DNAEmfietzoglouExcitationModel::~G4DNAEmfiet 78 G4DNAEmfietzoglouExcitationModel::~G4DNAEmfietzoglouExcitationModel()
77 { << 79 {
78 // Cross section 80 // Cross section
79 81
80 std::map< G4String,G4DNACrossSectionDataSe 82 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
81 for (pos = tableData.begin(); pos != table 83 for (pos = tableData.begin(); pos != tableData.end(); ++pos)
82 { 84 {
83 G4DNACrossSectionDataSet* table = pos- 85 G4DNACrossSectionDataSet* table = pos->second;
84 delete table; 86 delete table;
85 } 87 }
86 88
87 } 89 }
88 90
89 //....oooOO0OOooo........oooOO0OOooo........oo 91 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
90 92
91 void G4DNAEmfietzoglouExcitationModel::Initial 93 void G4DNAEmfietzoglouExcitationModel::Initialise(const G4ParticleDefinition* particle,
92 cons 94 const G4DataVector& /*cuts*/)
93 { 95 {
94 96
95 if (verboseLevel > 3) 97 if (verboseLevel > 3)
96 G4cout << "Calling G4DNAEmfietzoglouEx 98 G4cout << "Calling G4DNAEmfietzoglouExcitationModel::Initialise()" << G4endl;
97 99
98 G4String fileElectron("dna/sigma_excitatio 100 G4String fileElectron("dna/sigma_excitation_e_emfietzoglou");
99 101
100 G4ParticleDefinition* electronDef = G4Elec 102 G4ParticleDefinition* electronDef = G4Electron::ElectronDefinition();
101 103
102 G4String electron; 104 G4String electron;
103 105
104 G4double scaleFactor = (1.e-22 / 3.343) * 106 G4double scaleFactor = (1.e-22 / 3.343) * m*m;
105 107
106 // *** ELECTRON 108 // *** ELECTRON
107 109
108 electron = electronDef->GetParticleName(); 110 electron = electronDef->GetParticleName();
109 111
110 tableFile[electron] = fileElectron; 112 tableFile[electron] = fileElectron;
111 113
112 // Cross section << 114 lowEnergyLimit[electron] = 8. * eV;
>> 115 highEnergyLimit[electron] = 10. * keV;
113 116
114 auto tableE = new G4DNACrossSectionDataSe << 117 // Cross section
>> 118
>> 119 G4DNACrossSectionDataSet* tableE = new G4DNACrossSectionDataSet(new G4LogLogInterpolation, eV,scaleFactor );
115 tableE->LoadData(fileElectron); 120 tableE->LoadData(fileElectron);
116 121
117 tableData[electron] = tableE; 122 tableData[electron] = tableE;
118 << 123
119 // 124 //
120 125
>> 126 if (particle==electronDef)
>> 127 {
>> 128 SetLowEnergyLimit(lowEnergyLimit[electron]);
>> 129 SetHighEnergyLimit(highEnergyLimit[electron]);
>> 130 }
>> 131
121 if( verboseLevel>0 ) 132 if( verboseLevel>0 )
122 { 133 {
123 G4cout << "Emfietzoglou excitation model << 134 G4cout << "Emfietzoglou excitation model is initialized " << G4endl
124 << "Energy range: " << 135 << "Energy range: "
125 << LowEnergyLimit() / eV << " eV << 136 << LowEnergyLimit() / eV << " eV - "
126 << HighEnergyLimit() / keV << " k << 137 << HighEnergyLimit() / keV << " keV for "
127 << particle->GetParticleName() << 138 << particle->GetParticleName()
128 << G4endl; << 139 << G4endl;
129 } 140 }
130 141
131 // Initialize water density pointer 142 // Initialize water density pointer
132 fpMolWaterDensity = G4DNAMolecularMaterial 143 fpMolWaterDensity = G4DNAMolecularMaterial::Instance()->GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER"));
133 144
134 if (isInitialised) return; << 145 if (isInitialised) { return; }
135 fParticleChangeForGamma = GetParticleChang 146 fParticleChangeForGamma = GetParticleChangeForGamma();
136 isInitialised = true; 147 isInitialised = true;
137 } 148 }
138 149
139 //....oooOO0OOooo........oooOO0OOooo........oo 150 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
140 151
141 G4double G4DNAEmfietzoglouExcitationModel::Cro 152 G4double G4DNAEmfietzoglouExcitationModel::CrossSectionPerVolume(const G4Material* material,
142 153 const G4ParticleDefinition* particleDefinition,
143 154 G4double ekin,
144 155 G4double,
145 156 G4double)
146 { 157 {
147 if (verboseLevel > 3) 158 if (verboseLevel > 3)
148 G4cout << "Calling CrossSectionPerVolu 159 G4cout << "Calling CrossSectionPerVolume() of G4DNAEmfietzoglouExcitationModel" << G4endl;
149 160
150 if (particleDefinition != G4Electron::Elec 161 if (particleDefinition != G4Electron::ElectronDefinition()) return 0;
151 162
152 // Calculate total cross section for model 163 // Calculate total cross section for model
153 164
>> 165 G4double lowLim = 0;
>> 166 G4double highLim = 0;
154 G4double sigma=0; 167 G4double sigma=0;
155 168
156 G4double waterDensity = (*fpMolWaterDensit 169 G4double waterDensity = (*fpMolWaterDensity)[material->GetIndex()];
157 170
158 const G4String& particleName = particleDef << 171 if(waterDensity!= 0.0)
159 <<
160 if (ekin >= LowEnergyLimit() && ekin <= Hi <<
161 { 172 {
162 std::map< G4String,G4DNACrossSectionData << 173 const G4String& particleName = particleDefinition->GetParticleName();
163 pos = tableData.find(particleName); <<
164 174
165 if (pos != tableData.end()) << 175 std::map< G4String,G4double,std::less<G4String> >::iterator pos1;
166 { << 176 pos1 = lowEnergyLimit.find(particleName);
167 G4DNACrossSectionDataSet* table = pos- << 177 if (pos1 != lowEnergyLimit.end())
168 if (table != nullptr) sigma = table->F << 178 {
169 } << 179 lowLim = pos1->second;
170 else << 180 }
171 { << 181
172 G4Exception("G4DNAEmfietzoglouExcitati << 182 std::map< G4String,G4double,std::less<G4String> >::iterator pos2;
>> 183 pos2 = highEnergyLimit.find(particleName);
>> 184 if (pos2 != highEnergyLimit.end())
>> 185 {
>> 186 highLim = pos2->second;
>> 187 }
>> 188
>> 189 if (ekin >= lowLim && ekin < highLim)
>> 190 {
>> 191 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
>> 192 pos = tableData.find(particleName);
>> 193
>> 194 if (pos != tableData.end())
>> 195 {
>> 196 G4DNACrossSectionDataSet* table = pos->second;
>> 197 if (table != 0)
>> 198 {
>> 199 sigma = table->FindValue(ekin);
>> 200 }
>> 201 }
>> 202 else
>> 203 {
>> 204 G4Exception("G4DNAEmfietzoglouExcitationModel::CrossSectionPerVolume","em0002",
173 FatalException,"Mo 205 FatalException,"Model not applicable to particle type.");
174 } << 206 }
175 } << 207 }
176 208
177 if (verboseLevel > 2) << 209 if (verboseLevel > 2)
178 { << 210 {
179 G4cout << "_____________________________ << 211 G4cout << "__________________________________" << G4endl;
180 G4cout << "G4DNAEmfietzoglouExcitationMo << 212 G4cout << "G4DNAEmfietzoglouExcitationModel - XS INFO START" << G4endl;
181 G4cout << "Kinetic energy(eV)=" << ekin/ << 213 G4cout << "Kinetic energy(eV)=" << ekin/eV << " particle : " << particleName << G4endl;
182 G4cout << "Cross section per water molec << 214 G4cout << "Cross section per water molecule (cm^2)=" << sigma/cm/cm << G4endl;
183 G4cout << "Cross section per water molec << 215 G4cout << "Cross section per water molecule (cm^-1)=" << sigma*waterDensity/(1./cm) << G4endl;
184 //G4cout << " Cross section per water << 216 // G4cout << " Cross section per water molecule (cm^-1)=" << sigma*material->GetAtomicNumDensityVector()[1]/(1./cm) << G4endl;
185 ///sigma*material->GetAtomicNumDensityVe << 217 G4cout << "G4DNAEmfietzoglouExcitationModel - XS INFO END" << G4endl;
186 G4cout << "G4DNAEmfietzoglouExcitationMo << 218 }
187 } << 219
>> 220 } // if (waterMaterial)
188 221
189 return sigma*waterDensity; 222 return sigma*waterDensity;
190 } 223 }
191 224
192 //....oooOO0OOooo........oooOO0OOooo........oo 225 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
193 226
194 void G4DNAEmfietzoglouExcitationModel::SampleS 227 void G4DNAEmfietzoglouExcitationModel::SampleSecondaries(std::vector<G4DynamicParticle*>* /*fvect*/,
195 228 const G4MaterialCutsCouple* /*couple*/,
196 229 const G4DynamicParticle* aDynamicParticle,
197 230 G4double,
198 231 G4double)
199 { 232 {
200 233
201 if (verboseLevel > 3) 234 if (verboseLevel > 3)
202 G4cout << "Calling SampleSecondaries() 235 G4cout << "Calling SampleSecondaries() of G4DNAEmfietzoglouExcitationModel" << G4endl;
203 236
204 G4double k = aDynamicParticle->GetKineticE 237 G4double k = aDynamicParticle->GetKineticEnergy();
205 238
206 const G4String& particleName = aDynamicPar 239 const G4String& particleName = aDynamicParticle->GetDefinition()->GetParticleName();
207 240
208 G4int level = RandomSelect(k,particleName) 241 G4int level = RandomSelect(k,particleName);
209 G4double excitationEnergy = waterStructure 242 G4double excitationEnergy = waterStructure.ExcitationEnergy(level);
210 G4double newEnergy = k - excitationEnergy; 243 G4double newEnergy = k - excitationEnergy;
211 244
212 if (newEnergy > 0) 245 if (newEnergy > 0)
213 { 246 {
214 fParticleChangeForGamma->ProposeMoment 247 fParticleChangeForGamma->ProposeMomentumDirection(aDynamicParticle->GetMomentumDirection());
215 248
216 if (!statCode) fParticleChangeForGamma 249 if (!statCode) fParticleChangeForGamma->SetProposedKineticEnergy(newEnergy);
217 else fParticleChangeForGamma->SetPropo 250 else fParticleChangeForGamma->SetProposedKineticEnergy(k);
218 251
219 fParticleChangeForGamma->ProposeLocalE 252 fParticleChangeForGamma->ProposeLocalEnergyDeposit(excitationEnergy);
220 } 253 }
221 254
222 const G4Track * theIncomingTrack = fPartic 255 const G4Track * theIncomingTrack = fParticleChangeForGamma->GetCurrentTrack();
223 G4DNAChemistryManager::Instance()->CreateW 256 G4DNAChemistryManager::Instance()->CreateWaterMolecule(eExcitedMolecule,
224 257 level,
225 258 theIncomingTrack);
226 } 259 }
227 260
228 //....oooOO0OOooo........oooOO0OOooo........oo 261 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
229 262
230 G4int G4DNAEmfietzoglouExcitationModel::Random 263 G4int G4DNAEmfietzoglouExcitationModel::RandomSelect(G4double k, const G4String& particle)
231 { << 264 {
232 << 265
233 G4int level = 0; << 266 G4int level = 0;
234 267
235 std::map< G4String,G4DNACrossSectionDataSe 268 std::map< G4String,G4DNACrossSectionDataSet*,std::less<G4String> >::iterator pos;
236 pos = tableData.find(particle); 269 pos = tableData.find(particle);
237 270
238 if (pos != tableData.end()) 271 if (pos != tableData.end())
239 { 272 {
240 G4DNACrossSectionDataSet* table = pos- 273 G4DNACrossSectionDataSet* table = pos->second;
241 274
242 if (table != nullptr) << 275 if (table != 0)
243 { 276 {
244 auto valuesBuffer = new G4double[ << 277 G4double* valuesBuffer = new G4double[table->NumberOfComponents()];
245 const auto n = (G4int)table->Numb << 278 const size_t n(table->NumberOfComponents());
246 G4int i(n); << 279 size_t i(n);
247 G4double value = 0.; 280 G4double value = 0.;
248 281
249 //Check reading of initial xs file 282 //Check reading of initial xs file
250 //G4cout << table->GetComponent(0)-> << 283 //G4cout << table->GetComponent(0)->FindValue(k)/ ((1.e-22 / 3.343) * m*m) << G4endl;
251 //G4cout << table->GetComponent(1) 284 //G4cout << table->GetComponent(1)->FindValue(k)/ ((1.e-22 / 3.343) * m*m) << G4endl;
252 //G4cout << table->GetComponent(2) 285 //G4cout << table->GetComponent(2)->FindValue(k)/ ((1.e-22 / 3.343) * m*m) << G4endl;
253 //G4cout << table->GetComponent(3) 286 //G4cout << table->GetComponent(3)->FindValue(k)/ ((1.e-22 / 3.343) * m*m) << G4endl;
254 //G4cout << table->GetComponent(4) 287 //G4cout << table->GetComponent(4)->FindValue(k)/ ((1.e-22 / 3.343) * m*m) << G4endl;
255 //G4cout << table->GetComponent(5) 288 //G4cout << table->GetComponent(5)->FindValue(k)/ ((1.e-22 / 3.343) * m*m) << G4endl;
256 //G4cout << table->GetComponent(6) 289 //G4cout << table->GetComponent(6)->FindValue(k)/ ((1.e-22 / 3.343) * m*m) << G4endl;
257 //abort(); 290 //abort();
258 << 291
259 while (i>0) 292 while (i>0)
260 { 293 {
261 i--; 294 i--;
262 valuesBuffer[i] = table->GetCo 295 valuesBuffer[i] = table->GetComponent(i)->FindValue(k);
263 value += valuesBuffer[i]; 296 value += valuesBuffer[i];
264 } 297 }
265 298
266 value *= G4UniformRand(); 299 value *= G4UniformRand();
267 300
268 i = n; 301 i = n;
269 302
270 while (i > 0) 303 while (i > 0)
271 { 304 {
272 i--; 305 i--;
273 306
274 if (valuesBuffer[i] > value) 307 if (valuesBuffer[i] > value)
275 { 308 {
276 delete[] valuesBuffer; 309 delete[] valuesBuffer;
277 return i; 310 return i;
278 } 311 }
279 value -= valuesBuffer[i]; 312 value -= valuesBuffer[i];
280 } 313 }
281 314
282 delete[] valuesBuffer; << 315 if (valuesBuffer) delete[] valuesBuffer;
283 316
284 } 317 }
285 } 318 }
286 else 319 else
287 { 320 {
288 G4Exception("G4DNAEmfietzoglouExcitati 321 G4Exception("G4DNAEmfietzoglouExcitationModel::RandomSelect","em0002",
289 FatalException,"Model not 322 FatalException,"Model not applicable to particle type.");
290 } 323 }
291 return level; 324 return level;
292 } 325 }
293 326