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