| 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 // Created on 2016/01/18 26 // Created on 2016/01/18
27 // 27 //
28 // Authors: D. Sakata, W.G. Shin, S. Incerti 28 // Authors: D. Sakata, W.G. Shin, S. Incerti
29 // 29 //
30 // Based on a recent release of the ELSEPA cod 30 // Based on a recent release of the ELSEPA code
31 // developed and provided kindly by F. Salvat 31 // developed and provided kindly by F. Salvat et al.
32 // See 32 // See
33 // Computer Physics Communications, 165(2), 15 33 // Computer Physics Communications, 165(2), 157-190. (2005)
34 // http://dx.doi.org/10.1016/j.cpc.2004.09.006 34 // http://dx.doi.org/10.1016/j.cpc.2004.09.006
35 // 35 //
36 36
37 #include "G4DNAELSEPAElasticModel.hh" 37 #include "G4DNAELSEPAElasticModel.hh"
38 #include "G4PhysicalConstants.hh" 38 #include "G4PhysicalConstants.hh"
39 #include "G4SystemOfUnits.hh" 39 #include "G4SystemOfUnits.hh"
40 #include "G4DNAMolecularMaterial.hh" 40 #include "G4DNAMolecularMaterial.hh"
41 #include "G4EmParameters.hh" <<
42 41
43 //....oooOO0OOooo........oooOO0OOooo........oo 42 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
44 43
45 using namespace std; 44 using namespace std;
46 45
47 //....oooOO0OOooo........oooOO0OOooo........oo 46 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
48 47
49 G4DNAELSEPAElasticModel::G4DNAELSEPAElasticMod 48 G4DNAELSEPAElasticModel::G4DNAELSEPAElasticModel(const G4ParticleDefinition*,
50 const G4String& nam) : 49 const G4String& nam) :
51 G4VEmModel(nam) << 50 G4VEmModel(nam), isInitialised(false)
52 { 51 {
53 verboseLevel = 0; 52 verboseLevel = 0;
54 53
55 G4ProductionCutsTable* theCoupleTable = 54 G4ProductionCutsTable* theCoupleTable =
56 G4ProductionCutsTable::GetProductionCutsTabl 55 G4ProductionCutsTable::GetProductionCutsTable();
57 auto numOfCouples = (G4int)theCoupleTable-> << 56 G4int numOfCouples = (G4int)theCoupleTable->GetTableSize();
58 <<
59 fpBaseWater = G4Material::GetMaterial("G4_WA <<
60 57
61 for(G4int i=0; i<numOfCouples; ++i) 58 for(G4int i=0; i<numOfCouples; ++i)
62 { 59 {
63 const G4MaterialCutsCouple* couple = 60 const G4MaterialCutsCouple* couple =
64 theCoupleTable->GetMaterialCutsCouple 61 theCoupleTable->GetMaterialCutsCouple(i);
65 << 62 const G4Material* material = couple->GetMaterial();
66 const G4Material* material = couple->GetMa << 63 G4int nelm = (G4int)material->GetNumberOfElements();
67 if(!material) material = couple->GetMateri << 64 const G4ElementVector* theElementVector = material->GetElementVector();
68 <<
69 auto nelm = (G4int)material->GetNumberOfE <<
70 65
71 if(nelm==1) 66 if(nelm==1)
72 {// Protection: only for single element 67 {// Protection: only for single element
73 G4int Z = 79; 68 G4int Z = 79;
74 const G4ElementVector* theElementVector <<
75 Z = G4lrint((*theElementVector)[0]->Get 69 Z = G4lrint((*theElementVector)[0]->GetZ());
76 // Protection: only for GOLD 70 // Protection: only for GOLD
77 if (Z==79){ 71 if (Z==79){
78 fkillBelowEnergy_Au = 10. * eV; // Ki 72 fkillBelowEnergy_Au = 10. * eV; // Kills e- tracking
79 flowEnergyLimit = 0 * eV; // Must 73 flowEnergyLimit = 0 * eV; // Must stay at zero for killing
80 fhighEnergyLimit = 1 * GeV; // Defau 74 fhighEnergyLimit = 1 * GeV; // Default
81 SetLowEnergyLimit (flowEnergyLimit); 75 SetLowEnergyLimit (flowEnergyLimit);
82 SetHighEnergyLimit(fhighEnergyLimit); 76 SetHighEnergyLimit(fhighEnergyLimit);
83 }else{ 77 }else{
84 //continue; 78 //continue;
85 } 79 }
86 }else{// Protection: H2O only is available 80 }else{// Protection: H2O only is available
87 if(material==fpBaseWater){ << 81 if(material->GetName()=="G4_WATER"){
88 flowEnergyLimit = 10. * eV; << 82 flowEnergyLimit = 10. * eV;
89 fhighEnergyLimit = 1 * MeV; 83 fhighEnergyLimit = 1 * MeV;
90 SetLowEnergyLimit (flowEnergyLimit); 84 SetLowEnergyLimit (flowEnergyLimit);
91 SetHighEnergyLimit(fhighEnergyLimit); 85 SetHighEnergyLimit(fhighEnergyLimit);
92 }else{ 86 }else{
93 //continue; 87 //continue;
94 } 88 }
95 } 89 }
96 90
97 if (verboseLevel > 0) 91 if (verboseLevel > 0)
98 { 92 {
99 G4cout << "ELSEPA Elastic model is const 93 G4cout << "ELSEPA Elastic model is constructed for "
100 << material->GetName() << G4endl 94 << material->GetName() << G4endl
101 << "Energy range: " 95 << "Energy range: "
102 << flowEnergyLimit / eV << " eV - " 96 << flowEnergyLimit / eV << " eV - "
103 << fhighEnergyLimit / MeV << " MeV" 97 << fhighEnergyLimit / MeV << " MeV"
104 << G4endl; 98 << G4endl;
105 } 99 }
106 } 100 }
107 101
108 fParticleChangeForGamma = nullptr; << 102
109 fpMolDensity = nullptr; << 103 fParticleChangeForGamma = 0;
>> 104 fpMolDensity = 0;
110 105
111 fpData_Au=nullptr; 106 fpData_Au=nullptr;
112 fpData_H2O=nullptr; 107 fpData_H2O=nullptr;
113 } 108 }
114 109
115 //....oooOO0OOooo........oooOO0OOooo........oo 110 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
116 111
117 G4DNAELSEPAElasticModel::~G4DNAELSEPAElasticMo 112 G4DNAELSEPAElasticModel::~G4DNAELSEPAElasticModel()
118 { 113 {
119 delete fpData_Au; << 114 //std::map<G4int,G4DNACrossSectionDataSet*,
120 delete fpData_H2O; << 115 // std::less<G4String>>::iterator posZ;
>> 116 //for (posZ = tableZData.begin(); posZ != tableZData.end(); ++posZ)
>> 117 //{
>> 118 // G4DNACrossSectionDataSet* table = posZ->second;
>> 119 // delete table;
>> 120 //}
>> 121 //for (posZ = tableZData_Au.begin(); posZ != tableZData_Au.end(); ++posZ)
>> 122 //{
>> 123 // G4DNACrossSectionDataSet* table = posZ->second;
>> 124 // delete table;
>> 125 //}
>> 126 //for (posZ = tableZData_H2O.begin(); posZ != tableZData_H2O.end(); ++posZ)
>> 127 //{
>> 128 // G4DNACrossSectionDataSet* table = posZ->second;
>> 129 // delete table;
>> 130 //}
>> 131
>> 132 if(fpData_Au) delete fpData_Au;
>> 133 if(fpData_H2O) delete fpData_H2O;
>> 134
>> 135 //eEdummyVecZ.clear();
>> 136 //eCumZ.clear();
>> 137 //fAngleDataZ.clear();
121 138
122 eEdummyVec_Au.clear(); 139 eEdummyVec_Au.clear();
123 eEdummyVec_H2O.clear(); 140 eEdummyVec_H2O.clear();
124 eCum_Au.clear(); 141 eCum_Au.clear();
125 eCum_H2O.clear(); 142 eCum_H2O.clear();
126 fAngleData_Au.clear(); 143 fAngleData_Au.clear();
127 fAngleData_H2O.clear(); 144 fAngleData_H2O.clear();
128 } 145 }
129 146
130 //....oooOO0OOooo........oooOO0OOooo........oo 147 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
131 148
132 void G4DNAELSEPAElasticModel::Initialise(const 149 void G4DNAELSEPAElasticModel::Initialise(const G4ParticleDefinition* particle,
133 const G4DataVector& ) 150 const G4DataVector& )
134 { 151 {
135 if (verboseLevel > 3) 152 if (verboseLevel > 3)
136 G4cout << "Calling G4DNAELSEPAElasticModel:: 153 G4cout << "Calling G4DNAELSEPAElasticModel::Initialise()" << G4endl;
137 154
138 if (isInitialised) {return;} 155 if (isInitialised) {return;}
139 156
140 if(particle->GetParticleName() != "e-") 157 if(particle->GetParticleName() != "e-")
141 { 158 {
142 G4Exception("G4DNAELSEPAElasticModel::Init 159 G4Exception("G4DNAELSEPAElasticModel::Initialise","em0001",
143 FatalException,"Model not applicable to 160 FatalException,"Model not applicable to particle type.");
144 return; 161 return;
145 } 162 }
146 163
147 G4ProductionCutsTable* theCoupleTable = 164 G4ProductionCutsTable* theCoupleTable =
148 G4ProductionCutsTable::GetProductionCutsTabl 165 G4ProductionCutsTable::GetProductionCutsTable();
149 auto numOfCouples = (G4int)theCoupleTable-> << 166 G4int numOfCouples = (G4int)theCoupleTable->GetTableSize();
150 167
151 // UNIT OF TCS 168 // UNIT OF TCS
152 G4double scaleFactor = 1.*cm*cm; 169 G4double scaleFactor = 1.*cm*cm;
153 170
>> 171 //tableZData.clear();
>> 172 //tableZData_Au.clear();
>> 173 //tableZData_H2O.clear();
>> 174
154 fpData_Au=nullptr; 175 fpData_Au=nullptr;
155 fpData_H2O=nullptr; 176 fpData_H2O=nullptr;
156 fpBaseWater = G4Material::GetMaterial("G4_WA <<
157 177
158 for(G4int i=0; i<numOfCouples; ++i) 178 for(G4int i=0; i<numOfCouples; ++i)
159 { 179 {
160 const G4MaterialCutsCouple* couple = 180 const G4MaterialCutsCouple* couple =
161 theCoupleTable->GetMaterialCutsCouple 181 theCoupleTable->GetMaterialCutsCouple(i);
162 const G4Material* material = couple->GetMa << 182 const G4Material* material = couple->GetMaterial();
163 if(!material) material = couple->GetMateri << 183 const G4ElementVector* theElementVector = material->GetElementVector();
164 184
165 auto nelm = (G4int)material->GetNumberOfE << 185 G4int nelm = (G4int)material->GetNumberOfElements();
166 if (nelm==1){// Protection: only for singl 186 if (nelm==1){// Protection: only for single element
167 const G4ElementVector* theElementVector <<
168 G4int Z = G4lrint((*theElementVector)[0 187 G4int Z = G4lrint((*theElementVector)[0]->GetZ());
169 if (Z!=79)// Protection: only for GOLD 188 if (Z!=79)// Protection: only for GOLD
170 { 189 {
171 continue; 190 continue;
172 } 191 }
173 192
174 if (Z>0) 193 if (Z>0)
175 { 194 {
176 G4String fileZElectron("dna/sigma_elas 195 G4String fileZElectron("dna/sigma_elastic_e_elsepa_Z");
177 std::ostringstream oss; 196 std::ostringstream oss;
178 oss.str(""); 197 oss.str("");
179 oss.clear(stringstream::goodbit); 198 oss.clear(stringstream::goodbit);
180 oss << Z; 199 oss << Z;
181 fileZElectron += oss.str()+"_muffintin 200 fileZElectron += oss.str()+"_muffintin";
>> 201
>> 202 //G4DNACrossSectionDataSet* tableZE =
>> 203 // new G4DNACrossSectionDataSet
>> 204 // (new G4LogLogInterpolation, eV,scaleFactor );
>> 205 //tableZE->LoadData(fileZElectron);
>> 206 ////tableZData_Au[0] = tableZE;
>> 207 //tableZData[Z] = tableZE;
182 208
183 fpData_Au = new G4DNACrossSectionDataS 209 fpData_Au = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
184 210 eV,
185 211 scaleFactor );
186 fpData_Au->LoadData(fileZElectron); 212 fpData_Au->LoadData(fileZElectron);
187 213
188 std::ostringstream eFullFileNameZ; 214 std::ostringstream eFullFileNameZ;
189 const char *path = G4EmParameters::Ins << 215 const char *path = G4FindDataDir("G4LEDATA");
190 << 216 if (!path)
191 if (path == nullptr) <<
192 { 217 {
193 G4Exception("G4DNAELSEPAElasticModel 218 G4Exception("G4DNAELSEPAElasticModel::Initialise","em0002",
194 FatalException,"G4LEDATA environme 219 FatalException,"G4LEDATA environment variable not set.");
195 return; 220 return;
196 } 221 }
197 222
198 eFullFileNameZ.str(""); 223 eFullFileNameZ.str("");
199 eFullFileNameZ.clear(stringstream::goo 224 eFullFileNameZ.clear(stringstream::goodbit);
200 225
201 eFullFileNameZ 226 eFullFileNameZ
202 << path 227 << path
203 << "/dna/sigmadiff_cumulated_elastic 228 << "/dna/sigmadiff_cumulated_elastic_e_elsepa_Z"
204 << Z << "_muffintin.dat"; 229 << Z << "_muffintin.dat";
205 230
206 std::ifstream eDiffCrossSectionZ(eFull 231 std::ifstream eDiffCrossSectionZ(eFullFileNameZ.str().c_str());
207 232
208 if (!eDiffCrossSectionZ) 233 if (!eDiffCrossSectionZ)
209 { 234 {
210 G4Exception("G4DNAELSEPAElasticModel 235 G4Exception("G4DNAELSEPAElasticModel::Initialise","em0003",
211 FatalException,"Missing data file 236 FatalException,"Missing data file for cumulated DCS");
212 return; 237 return;
213 } 238 }
>> 239
>> 240 //eEdummyVecZ.clear();
>> 241 //eCumZ.clear();
>> 242 //fAngleDataZ.clear();
214 243
215 eEdummyVec_Au.clear(); 244 eEdummyVec_Au.clear();
216 eCum_Au.clear(); 245 eCum_Au.clear();
217 fAngleData_Au.clear(); 246 fAngleData_Au.clear();
218 247
>> 248 //eEdummyVecZ[Z].push_back(0.);
219 eEdummyVec_Au.push_back(0.); 249 eEdummyVec_Au.push_back(0.);
220 do 250 do
221 { 251 {
222 G4double eDummy; 252 G4double eDummy;
223 G4double cumDummy; 253 G4double cumDummy;
224 eDiffCrossSectionZ>>eDummy>>cumDummy 254 eDiffCrossSectionZ>>eDummy>>cumDummy;
>> 255 //if (eDummy != eEdummyVecZ[Z].back())
225 if (eDummy != eEdummyVec_Au.back()) 256 if (eDummy != eEdummyVec_Au.back())
226 { 257 {
>> 258
>> 259 //eEdummyVecZ[Z].push_back(eDummy);
227 eEdummyVec_Au.push_back(eDummy); 260 eEdummyVec_Au.push_back(eDummy);
>> 261 //eCumZ[Z][eDummy].push_back(0.);
228 eCum_Au[eDummy].push_back(0.); 262 eCum_Au[eDummy].push_back(0.);
229 } 263 }
>> 264 //eDiffCrossSectionZ>>fAngleDataZ[Z][eDummy][cumDummy];
230 eDiffCrossSectionZ>>fAngleData_Au[eD 265 eDiffCrossSectionZ>>fAngleData_Au[eDummy][cumDummy];
>> 266 //if (cumDummy != eCumZ[Z][eDummy].back())
231 if (cumDummy != eCum_Au[eDummy].back 267 if (cumDummy != eCum_Au[eDummy].back())
232 { 268 {
>> 269 //eCumZ[Z][eDummy].push_back(cumDummy);
233 eCum_Au[eDummy].push_back(cumDummy 270 eCum_Au[eDummy].push_back(cumDummy);
234 } 271 }
235 }while(!eDiffCrossSectionZ.eof()); 272 }while(!eDiffCrossSectionZ.eof());
236 } 273 }
237 274
238 }else{// Protection: H2O only is available 275 }else{// Protection: H2O only is available
239 if(material == fpBaseWater && !fpData_H2 << 276 if(material->GetName()=="G4_WATER"){
240 if (LowEnergyLimit() < 10*eV) 277 if (LowEnergyLimit() < 10*eV)
241 { 278 {
242 G4cout<<"G4DNAELSEPAElasticModel: lo 279 G4cout<<"G4DNAELSEPAElasticModel: low energy limit increased from "
243 << LowEnergyLimit()/eV << " eV 280 << LowEnergyLimit()/eV << " eV to " << 10 << " eV"
244 << G4endl; 281 << G4endl;
245 SetLowEnergyLimit(10.*eV); 282 SetLowEnergyLimit(10.*eV);
246 } 283 }
247 284
248 if (HighEnergyLimit() > 1.*MeV) 285 if (HighEnergyLimit() > 1.*MeV)
249 { 286 {
250 G4cout<<"G4DNAELSEPAElasticModel: hi 287 G4cout<<"G4DNAELSEPAElasticModel: high energy limit decreased from "
251 << HighEnergyLimit()/MeV << " 288 << HighEnergyLimit()/MeV << " MeV to " << 1. << " MeV"
252 << G4endl; 289 << G4endl;
253 SetHighEnergyLimit(1.*MeV); 290 SetHighEnergyLimit(1.*MeV);
254 } 291 }
255 292
256 G4String fileZElectron("dna/sigma_elas 293 G4String fileZElectron("dna/sigma_elastic_e_elsepa_muffin");
257 294
>> 295 //G4DNACrossSectionDataSet* tableZE =
>> 296 // new G4DNACrossSectionDataSet(
>> 297 // new G4LogLogInterpolation, eV,scaleFactor );
>> 298 //tableZE->LoadData(fileZElectron);
>> 299 ////tableZData_H2O[0] = tableZE;
>> 300 //tableZData[0] = tableZE;
>> 301
258 fpData_H2O = new G4DNACrossSectionData 302 fpData_H2O = new G4DNACrossSectionDataSet(new G4LogLogInterpolation,
259 303 eV,
260 304 scaleFactor );
261 fpData_H2O->LoadData(fileZElectron); 305 fpData_H2O->LoadData(fileZElectron);
262 306
263 std::ostringstream eFullFileNameZ; 307 std::ostringstream eFullFileNameZ;
264 308
265 const char *path = G4EmParameters::Ins << 309 const char *path = G4FindDataDir("G4LEDATA");
266 if (path == nullptr) << 310 if (!path)
267 { 311 {
268 G4Exception("G4DNAELSEPAElasticModel 312 G4Exception("G4DNAELSEPAElasticModel::Initialise","em0004",
269 FatalException,"G4LEDATA environme 313 FatalException,"G4LEDATA environment variable not set.");
270 return; 314 return;
271 } 315 }
272 316
273 eFullFileNameZ.str(""); 317 eFullFileNameZ.str("");
274 eFullFileNameZ.clear(stringstream::goo 318 eFullFileNameZ.clear(stringstream::goodbit);
275 319
276 eFullFileNameZ 320 eFullFileNameZ
277 << path 321 << path
278 << "/dna/sigmadiff_cumulated_elasti 322 << "/dna/sigmadiff_cumulated_elastic_e_elsepa_muffin.dat";
279 323
280 std::ifstream eDiffCrossSectionZ(eFull 324 std::ifstream eDiffCrossSectionZ(eFullFileNameZ.str().c_str());
281 325
282 if (!eDiffCrossSectionZ) 326 if (!eDiffCrossSectionZ)
283 G4Exception("G4DNAELSEPAElasticModel: 327 G4Exception("G4DNAELSEPAElasticModel::Initialise","em0005",
284 FatalException, 328 FatalException,
285 "Missing data file for cumulated DCS" 329 "Missing data file for cumulated DCS");
286 330
>> 331 //eEdummyVecZ.clear();
>> 332 //eCumZ.clear();
>> 333 //fAngleDataZ.clear();
>> 334
287 eEdummyVec_H2O.clear(); 335 eEdummyVec_H2O.clear();
288 eCum_H2O.clear(); 336 eCum_H2O.clear();
289 fAngleData_H2O.clear(); 337 fAngleData_H2O.clear();
290 338
>> 339 //eEdummyVecZ[0].push_back(0.);
291 eEdummyVec_H2O.push_back(0.); 340 eEdummyVec_H2O.push_back(0.);
292 341
293 do 342 do
294 { 343 {
295 G4double eDummy; 344 G4double eDummy;
296 G4double cumDummy; 345 G4double cumDummy;
297 eDiffCrossSectionZ>>eDummy>>cumDummy 346 eDiffCrossSectionZ>>eDummy>>cumDummy;
>> 347 //if (eDummy != eEdummyVecZ[0].back())
298 if (eDummy != eEdummyVec_H2O.back()) 348 if (eDummy != eEdummyVec_H2O.back())
299 { 349 {
>> 350 //eEdummyVecZ[0].push_back(eDummy);
300 eEdummyVec_H2O.push_back(eDummy); 351 eEdummyVec_H2O.push_back(eDummy);
>> 352 //eCumZ[0][eDummy].push_back(0.);
301 eCum_H2O[eDummy].push_back(0.); 353 eCum_H2O[eDummy].push_back(0.);
302 } 354 }
>> 355 //eDiffCrossSectionZ>>fAngleDataZ[0][eDummy][cumDummy];
303 eDiffCrossSectionZ>>fAngleData_H2O[e 356 eDiffCrossSectionZ>>fAngleData_H2O[eDummy][cumDummy];
>> 357 //if (cumDummy != eCumZ[0][eDummy].back()){
304 if (cumDummy != eCum_H2O[eDummy].bac 358 if (cumDummy != eCum_H2O[eDummy].back()){
>> 359 //eCumZ[0][eDummy].push_back(cumDummy);
305 eCum_H2O[eDummy].push_back(cumDumm 360 eCum_H2O[eDummy].push_back(cumDummy);
306 } 361 }
307 }while(!eDiffCrossSectionZ.eof()); 362 }while(!eDiffCrossSectionZ.eof());
308 } 363 }
309 } 364 }
310 if (verboseLevel > 2) 365 if (verboseLevel > 2)
311 G4cout << "Loaded cross section files of E 366 G4cout << "Loaded cross section files of ELSEPA Elastic model for"
312 << material->GetName() << G4endl; 367 << material->GetName() << G4endl;
313 368
314 if( verboseLevel>0 ) 369 if( verboseLevel>0 )
315 { 370 {
316 G4cout << "ELSEPA elastic model is initi 371 G4cout << "ELSEPA elastic model is initialized " << G4endl
317 << "Energy range: " 372 << "Energy range: "
318 << LowEnergyLimit() / eV << " eV - " 373 << LowEnergyLimit() / eV << " eV - "
319 << HighEnergyLimit()/ MeV << " MeV" 374 << HighEnergyLimit()/ MeV << " MeV"
320 << G4endl; 375 << G4endl;
321 } 376 }
322 } // Loop on couples 377 } // Loop on couples
323 378
324 379
325 fParticleChangeForGamma = GetParticleChangeF 380 fParticleChangeForGamma = GetParticleChangeForGamma();
326 << 381 fpMolDensity = 0;
327 fpMolDensity = <<
328 G4DNAMolecularMaterial::Instance()-> <<
329 GetNumMolPerVolTableFor(G4Material::GetMater <<
330 382
331 isInitialised = true; 383 isInitialised = true;
332 } 384 }
333 385
334 //....oooOO0OOooo........oooOO0OOooo........oo 386 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
335 387
336 G4double G4DNAELSEPAElasticModel::CrossSection 388 G4double G4DNAELSEPAElasticModel::CrossSectionPerVolume
337 (const G4Material* material, 389 (const G4Material* material,
338 const G4ParticleDefinition* particle, 390 const G4ParticleDefinition* particle,
339 G4double ekin, 391 G4double ekin,
340 G4double, 392 G4double,
341 G4double) 393 G4double)
342 { 394 {
343 395
344 if (verboseLevel > 3) 396 if (verboseLevel > 3)
345 { 397 {
346 G4cout << 398 G4cout <<
347 "Calling CrossSectionPerVolume() of G4DNAE 399 "Calling CrossSectionPerVolume() of G4DNAELSEPAElasticModel"
348 << G4endl; 400 << G4endl;
349 } 401 }
350 402
351 G4double atomicNDensity=0.0; 403 G4double atomicNDensity=0.0;
352 G4double sigma=0; 404 G4double sigma=0;
353 405
>> 406 const G4ElementVector* theElementVector = material->GetElementVector();
354 std::size_t nelm = material->GetNumberOfElem 407 std::size_t nelm = material->GetNumberOfElements();
355 if (nelm==1) // Protection: only for single 408 if (nelm==1) // Protection: only for single element
356 { 409 {
357 // Protection: only for GOLD 410 // Protection: only for GOLD
358 if (material->GetZ()!=79) return 0.0; 411 if (material->GetZ()!=79) return 0.0;
359 << 412
360 const G4ElementVector* theElementVector = <<
361 G4int Z = G4lrint((*theElementVector)[0]-> 413 G4int Z = G4lrint((*theElementVector)[0]->GetZ());
362 414
363 const G4String& particleName = particle->G 415 const G4String& particleName = particle->GetParticleName();
364 atomicNDensity = material->GetAtomicNumDen 416 atomicNDensity = material->GetAtomicNumDensityVector()[0];
365 if(atomicNDensity!= 0.0) 417 if(atomicNDensity!= 0.0)
366 { 418 {
367 if (ekin < fhighEnergyLimit) 419 if (ekin < fhighEnergyLimit)
368 { 420 {
369 if (ekin < fkillBelowEnergy_Au) return 421 if (ekin < fkillBelowEnergy_Au) return DBL_MAX;
370 422
>> 423 //std::map< G4int,G4DNACrossSectionDataSet*,
>> 424 // std::less<G4String> >::iterator pos;
>> 425 ////pos = tableZData_Au.find(0);
>> 426 //pos = tableZData.find(Z);
>> 427 //
>> 428 ////if (pos != tableZData_Au.end())
>> 429 //if (pos != tableZData.end())
>> 430 //{
>> 431 // G4DNACrossSectionDataSet* table = pos->second;
>> 432 // if (table != 0)
>> 433 // {
>> 434 // // XS takes its 10 eV value below 10 eV for GOLD
>> 435 // if (ekin < 10*eV) sigma = table->FindValue(10*eV);
>> 436 // else sigma = table->FindValue(ekin);
>> 437 // }
>> 438 //}
>> 439 //else
>> 440 //{
>> 441 // G4Exception("G4DNAELSEPAElasticModel::ComputeCrossSectionPerVolume",
>> 442 // "em0006",FatalException,"Model not applicable to particle type.");
>> 443 //}
>> 444
371 if (ekin < 10*eV) sigma = fpData_Au->F 445 if (ekin < 10*eV) sigma = fpData_Au->FindValue(10*eV);
372 else sigma = fpData_Au->F 446 else sigma = fpData_Au->FindValue(ekin);
373 } 447 }
374 } 448 }
375 if (verboseLevel > 2) 449 if (verboseLevel > 2)
376 { 450 {
377 G4cout << "_____________________________ 451 G4cout << "__________________________________" << G4endl;
378 G4cout << "=== G4DNAELSEPAElasticModel - 452 G4cout << "=== G4DNAELSEPAElasticModel - XS INFO START" << G4endl;
379 G4cout << "=== Material is made of one e 453 G4cout << "=== Material is made of one element with Z =" << Z << G4endl;
380 G4cout << "=== Kinetic energy(eV)=" << e 454 G4cout << "=== Kinetic energy(eV)=" << ekin/eV << " particle : "
381 << particleName << G4endl; 455 << particleName << G4endl;
382 G4cout << "=== Cross section per atom fo 456 G4cout << "=== Cross section per atom for Z="<<Z<<" is (cm^2)"
383 << sigma/cm/cm << G4endl; 457 << sigma/cm/cm << G4endl;
384 G4cout << "=== Cross section per atom fo 458 G4cout << "=== Cross section per atom for Z="<<Z<<" is (cm^-1)="
385 << sigma*atomicNDensity/(1./cm) < 459 << sigma*atomicNDensity/(1./cm) << G4endl;
386 G4cout << "=== G4DNAELSEPAElasticModel - 460 G4cout << "=== G4DNAELSEPAElasticModel - XS INFO END" << G4endl;
387 } 461 }
388 } 462 }
389 else 463 else
390 { 464 {
>> 465 fpMolDensity =
>> 466 G4DNAMolecularMaterial::Instance()->
>> 467 GetNumMolPerVolTableFor(G4Material::GetMaterial("G4_WATER"));
391 atomicNDensity = (*fpMolDensity)[material- 468 atomicNDensity = (*fpMolDensity)[material->GetIndex()];
392 if(atomicNDensity!= 0.0) 469 if(atomicNDensity!= 0.0)
393 { 470 {
394 if (ekin < HighEnergyLimit() && ekin >= 471 if (ekin < HighEnergyLimit() && ekin >= LowEnergyLimit())
395 { 472 {
>> 473 //std::map< G4int,G4DNACrossSectionDataSet*,
>> 474 //std::less<G4String> >::iterator pos;
>> 475 ////pos = tableZData_H2O.find(0); // the data is stored as Z=0
>> 476 //pos = tableZData.find(0); // the data is stored as Z=0
>> 477 ////SI : XS must not be zero
>> 478 //// otherwise sampling of secondaries method ignored
>> 479 ////if (pos != tableZData_H2O.end())
>> 480 //if (pos != tableZData.end())
>> 481 //{
>> 482 // G4DNACrossSectionDataSet* table = pos->second;
>> 483 // if (table != 0)
>> 484 // {
>> 485 // sigma = table->FindValue(ekin);
>> 486 // }
>> 487 //}
>> 488
396 sigma = fpData_H2O->FindValue(ekin); 489 sigma = fpData_H2O->FindValue(ekin);
397 } 490 }
398 } 491 }
399 if (verboseLevel > 2) 492 if (verboseLevel > 2)
400 { 493 {
401 G4cout << "_____________________________ 494 G4cout << "__________________________________" << G4endl;
402 G4cout << "=== G4DNAELSEPAElasticModel - 495 G4cout << "=== G4DNAELSEPAElasticModel - XS INFO START" << G4endl;
403 G4cout << "=== Kinetic energy(eV)=" << e 496 G4cout << "=== Kinetic energy(eV)=" << ekin/eV
404 << " particle : " << particle->Ge 497 << " particle : " << particle->GetParticleName() << G4endl;
405 G4cout << "=== Cross section per water m 498 G4cout << "=== Cross section per water molecule (cm^2)="
406 << sigma/cm/cm << G4endl; 499 << sigma/cm/cm << G4endl;
407 G4cout << "=== Cross section per water m 500 G4cout << "=== Cross section per water molecule (cm^-1)="
408 << sigma*atomicNDensity/(1./cm) < 501 << sigma*atomicNDensity/(1./cm) << G4endl;
409 G4cout << "=== G4DNAELSEPAElasticModel - 502 G4cout << "=== G4DNAELSEPAElasticModel - XS INFO END" << G4endl;
410 } 503 }
411 } 504 }
412 505
413 return sigma*atomicNDensity; 506 return sigma*atomicNDensity;
414 } 507 }
415 508
416 //....oooOO0OOooo........oooOO0OOooo........oo 509 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
417 510
418 void G4DNAELSEPAElasticModel::SampleSecondarie 511 void G4DNAELSEPAElasticModel::SampleSecondaries(
419 std::vector<G4DynamicParticle*>*, 512 std::vector<G4DynamicParticle*>*,
420 const G4MaterialCutsCouple* couple, 513 const G4MaterialCutsCouple* couple,
421 const G4DynamicParticle* aDynamicElectro 514 const G4DynamicParticle* aDynamicElectron,
422 G4double, 515 G4double,
423 G4double) 516 G4double)
424 { 517 {
425 518
426 if (verboseLevel > 3){ 519 if (verboseLevel > 3){
427 G4cout << 520 G4cout <<
428 "Calling SampleSecondaries() of G4DNAELSEP 521 "Calling SampleSecondaries() of G4DNAELSEPAElasticModel"
429 << G4endl; 522 << G4endl;
430 } 523 }
431 524
432 G4double electronEnergy0 = aDynamicElectron- 525 G4double electronEnergy0 = aDynamicElectron->GetKineticEnergy();
433 526
434 const G4Material* material = couple->GetMate << 527 const G4Material* material = couple->GetMaterial();
435 if(!material) material = couple->GetMaterial << 528 const G4ElementVector* theElementVector = material->GetElementVector();
436 <<
437 std::size_t nelm = material->GetNumberOfElem 529 std::size_t nelm = material->GetNumberOfElements();
438 if (nelm==1) // Protection: only for single 530 if (nelm==1) // Protection: only for single element
439 { 531 {
440 const G4ElementVector* theElementVector = <<
441 G4int Z = G4lrint((*theElementVector)[0]- 532 G4int Z = G4lrint((*theElementVector)[0]->GetZ());
442 if (Z!=79) return; 533 if (Z!=79) return;
443 if (electronEnergy0 < fkillBelowEnergy_Au) 534 if (electronEnergy0 < fkillBelowEnergy_Au)
444 { 535 {
445 fParticleChangeForGamma->SetProposedKine 536 fParticleChangeForGamma->SetProposedKineticEnergy(0.);
446 fParticleChangeForGamma->ProposeMomentum 537 fParticleChangeForGamma->ProposeMomentumDirection(G4ThreeVector(0,0,0));
447 fParticleChangeForGamma->ProposeTrackSta 538 fParticleChangeForGamma->ProposeTrackStatus(fStopAndKill);
448 fParticleChangeForGamma->ProposeLocalEne 539 fParticleChangeForGamma->ProposeLocalEnergyDeposit(electronEnergy0);
449 return; 540 return;
450 } 541 }
451 542
452 if(electronEnergy0>= fkillBelowEnergy_Au & 543 if(electronEnergy0>= fkillBelowEnergy_Au && electronEnergy0 < fhighEnergyLimit)
453 { 544 {
454 G4double cosTheta = 0; 545 G4double cosTheta = 0;
455 if (electronEnergy0>=10*eV) 546 if (electronEnergy0>=10*eV)
456 { 547 {
457 cosTheta = RandomizeCosTheta(Z,electro 548 cosTheta = RandomizeCosTheta(Z,electronEnergy0);
458 } 549 }
459 else 550 else
460 { 551 {
461 cosTheta = RandomizeCosTheta(Z,10*eV); 552 cosTheta = RandomizeCosTheta(Z,10*eV);
462 } 553 }
463 554
464 G4double phi = 2. * CLHEP::pi * G4Unifor 555 G4double phi = 2. * CLHEP::pi * G4UniformRand();
465 556
466 G4ThreeVector zVers = aDynamicElectron-> 557 G4ThreeVector zVers = aDynamicElectron->GetMomentumDirection();
467 G4ThreeVector xVers = zVers.orthogonal() 558 G4ThreeVector xVers = zVers.orthogonal();
468 G4ThreeVector yVers = zVers.cross(xVers) 559 G4ThreeVector yVers = zVers.cross(xVers);
469 560
470 G4double xDir = std::sqrt(1. - cosTheta* 561 G4double xDir = std::sqrt(1. - cosTheta*cosTheta);
471 G4double yDir = xDir; 562 G4double yDir = xDir;
472 xDir *= std::cos(phi); 563 xDir *= std::cos(phi);
473 yDir *= std::sin(phi); 564 yDir *= std::sin(phi);
474 565
475 G4ThreeVector zPrimeVers((xDir*xVers + y 566 G4ThreeVector zPrimeVers((xDir*xVers + yDir*yVers + cosTheta*zVers));
476 fParticleChangeForGamma->ProposeMomentum 567 fParticleChangeForGamma->ProposeMomentumDirection(zPrimeVers.unit());
477 fParticleChangeForGamma->SetProposedKine 568 fParticleChangeForGamma->SetProposedKineticEnergy(electronEnergy0);
478 569
479 } 570 }
480 } 571 }
481 else 572 else
482 { 573 {
483 if(material == fpBaseWater) << 574 if(material->GetName()=="G4_WATER")
484 { 575 {
485 //The data for water is stored as Z=0 576 //The data for water is stored as Z=0
486 G4double cosTheta = RandomizeCosTheta(0, 577 G4double cosTheta = RandomizeCosTheta(0,electronEnergy0);
487 578
488 G4double phi = 2. * pi * G4UniformRand() 579 G4double phi = 2. * pi * G4UniformRand();
489 580
490 G4ThreeVector zVers = aDynamicElectron-> 581 G4ThreeVector zVers = aDynamicElectron->GetMomentumDirection();
491 G4ThreeVector xVers = zVers.orthogonal() 582 G4ThreeVector xVers = zVers.orthogonal();
492 G4ThreeVector yVers = zVers.cross(xVers) 583 G4ThreeVector yVers = zVers.cross(xVers);
493 584
494 G4double xDir = std::sqrt(1. - cosTheta* 585 G4double xDir = std::sqrt(1. - cosTheta*cosTheta);
495 G4double yDir = xDir; 586 G4double yDir = xDir;
496 xDir *= std::cos(phi); 587 xDir *= std::cos(phi);
497 yDir *= std::sin(phi); 588 yDir *= std::sin(phi);
498 589
499 G4ThreeVector zPrimeVers((xDir*xVers + y 590 G4ThreeVector zPrimeVers((xDir*xVers + yDir*yVers + cosTheta*zVers));
500 fParticleChangeForGamma->ProposeMomentum 591 fParticleChangeForGamma->ProposeMomentumDirection(zPrimeVers.unit());
501 fParticleChangeForGamma->SetProposedKine 592 fParticleChangeForGamma->SetProposedKineticEnergy(electronEnergy0);
502 } 593 }
503 } 594 }
504 } 595 }
505 596
506 //....oooOO0OOooo........oooOO0OOooo........oo 597 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
507 598
508 G4double G4DNAELSEPAElasticModel::Theta(G4int 599 G4double G4DNAELSEPAElasticModel::Theta(G4int Z,
509 G4ParticleDefinition * particleDefini 600 G4ParticleDefinition * particleDefinition,
510 G4double k, 601 G4double k,
511 G4double integrDiff) 602 G4double integrDiff)
512 { 603 {
513 604
514 G4double theta = 0.; 605 G4double theta = 0.;
515 G4double valueE1 = 0.; 606 G4double valueE1 = 0.;
516 G4double valueE2 = 0.; 607 G4double valueE2 = 0.;
517 G4double valuecum21 = 0.; 608 G4double valuecum21 = 0.;
518 G4double valuecum22 = 0.; 609 G4double valuecum22 = 0.;
519 G4double valuecum12 = 0.; 610 G4double valuecum12 = 0.;
520 G4double valuecum11 = 0.; 611 G4double valuecum11 = 0.;
521 G4double a11 = 0.; 612 G4double a11 = 0.;
522 G4double a12 = 0.; 613 G4double a12 = 0.;
523 G4double a21 = 0.; 614 G4double a21 = 0.;
524 G4double a22 = 0.; 615 G4double a22 = 0.;
525 616
526 if (particleDefinition == G4Electron::Electro 617 if (particleDefinition == G4Electron::ElectronDefinition())
527 { 618 {
528 << 619 //std::vector<G4double>::iterator e2
>> 620 // = std::upper_bound(eEdummyVecZ[Z].begin(),
>> 621 // eEdummyVecZ[Z].end(), k);
529 std::vector<G4double>::iterator e2; 622 std::vector<G4double>::iterator e2;
530 if(Z==0){ 623 if(Z==0){
531 e2 = std::upper_bound(eEdummyVec_H2O.begin 624 e2 = std::upper_bound(eEdummyVec_H2O.begin(),
532 eEdummyVec_H2O.end() 625 eEdummyVec_H2O.end(), k);
533 }else if (Z==79){ 626 }else if (Z==79){
534 e2 = std::upper_bound(eEdummyVec_Au.begin( 627 e2 = std::upper_bound(eEdummyVec_Au.begin(),
535 eEdummyVec_Au.end(), 628 eEdummyVec_Au.end(), k);
536 } 629 }
537 630
538 auto e1 = e2 - 1; << 631 std::vector<G4double>::iterator e1 = e2 - 1;
539 632
>> 633 //std::vector<G4double>::iterator cum12
>> 634 // = std::upper_bound(eCumZ[Z][(*e1)].begin(),
>> 635 // eCumZ[Z][(*e1)].end(),integrDiff);
540 std::vector<G4double>::iterator cum12; 636 std::vector<G4double>::iterator cum12;
541 if(Z==0){ 637 if(Z==0){
542 cum12 = std::upper_bound(eCum_H2O[(*e1)] 638 cum12 = std::upper_bound(eCum_H2O[(*e1)].begin(),
543 eCum_H2O[(*e1)] 639 eCum_H2O[(*e1)].end(),integrDiff);
544 }else if (Z==79){ 640 }else if (Z==79){
545 cum12 = std::upper_bound(eCum_Au[(*e1)]. 641 cum12 = std::upper_bound(eCum_Au[(*e1)].begin(),
546 eCum_Au[(*e1)]. 642 eCum_Au[(*e1)].end(),integrDiff);
547 } 643 }
548 644
549 auto cum11 = cum12 - 1; << 645 std::vector<G4double>::iterator cum11 = cum12 - 1;
550 646
551 //std::vector<G4double>::iterator cum22 647 //std::vector<G4double>::iterator cum22
552 // = std::upper_bound(eCumZ[Z][(*e 648 // = std::upper_bound(eCumZ[Z][(*e2)].begin(),
553 // eCumZ[Z][(*e2)].end(),integrDif 649 // eCumZ[Z][(*e2)].end(),integrDiff);
554 std::vector<G4double>::iterator cum22; 650 std::vector<G4double>::iterator cum22;
555 if(Z==0){ 651 if(Z==0){
556 cum22 = std::upper_bound(eCum_H2O[(*e2)]. 652 cum22 = std::upper_bound(eCum_H2O[(*e2)].begin(),
557 eCum_H2O[(*e2)]. 653 eCum_H2O[(*e2)].end(),integrDiff);
558 }else if(Z==79){ 654 }else if(Z==79){
559 cum22 = std::upper_bound(eCum_Au[(*e2)].b 655 cum22 = std::upper_bound(eCum_Au[(*e2)].begin(),
560 eCum_Au[(*e2)].e 656 eCum_Au[(*e2)].end(),integrDiff);
561 } 657 }
562 658
563 auto cum21 = cum22 - 1; << 659 std::vector<G4double>::iterator cum21 = cum22 - 1;
564 660
565 valueE1 = *e1; 661 valueE1 = *e1;
566 valueE2 = *e2; 662 valueE2 = *e2;
567 valuecum11 = *cum11; 663 valuecum11 = *cum11;
568 valuecum12 = *cum12; 664 valuecum12 = *cum12;
569 valuecum21 = *cum21; 665 valuecum21 = *cum21;
570 valuecum22 = *cum22; 666 valuecum22 = *cum22;
571 667
>> 668
>> 669 //a11 = fAngleDataZ[Z][valueE1][valuecum11];
>> 670 //a12 = fAngleDataZ[Z][valueE1][valuecum12];
>> 671 //a21 = fAngleDataZ[Z][valueE2][valuecum21];
>> 672 //a22 = fAngleDataZ[Z][valueE2][valuecum22];
572 if(Z==0){ 673 if(Z==0){
573 a11 = fAngleData_H2O[valueE1][valuecum11]; 674 a11 = fAngleData_H2O[valueE1][valuecum11];
574 a12 = fAngleData_H2O[valueE1][valuecum12]; 675 a12 = fAngleData_H2O[valueE1][valuecum12];
575 a21 = fAngleData_H2O[valueE2][valuecum21]; 676 a21 = fAngleData_H2O[valueE2][valuecum21];
576 a22 = fAngleData_H2O[valueE2][valuecum22]; 677 a22 = fAngleData_H2O[valueE2][valuecum22];
577 }else if (Z==79){ 678 }else if (Z==79){
578 a11 = fAngleData_Au[valueE1][valuecum11]; 679 a11 = fAngleData_Au[valueE1][valuecum11];
579 a12 = fAngleData_Au[valueE1][valuecum12]; 680 a12 = fAngleData_Au[valueE1][valuecum12];
580 a21 = fAngleData_Au[valueE2][valuecum21]; 681 a21 = fAngleData_Au[valueE2][valuecum21];
581 a22 = fAngleData_Au[valueE2][valuecum22]; 682 a22 = fAngleData_Au[valueE2][valuecum22];
582 } 683 }
583 684
584 } 685 }
585 686
586 if (a11 == 0 && a12 == 0 && a21 == 0 && a22 = 687 if (a11 == 0 && a12 == 0 && a21 == 0 && a22 == 0) return (0.);
587 688
588 theta = QuadInterpolator(valuecum11, valuecum 689 theta = QuadInterpolator(valuecum11, valuecum12, valuecum21, valuecum22,
589 a11, a12,a21, a22, valueE1, valueE2, 690 a11, a12,a21, a22, valueE1, valueE2, k, integrDiff);
590 return theta; 691 return theta;
591 } 692 }
592 693
593 //....oooOO0OOooo........oooOO0OOooo........oo 694 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
594 // 695 //
595 G4double G4DNAELSEPAElasticModel::LogLinInterp 696 G4double G4DNAELSEPAElasticModel::LogLinInterpolate(G4double e1,
596 G4double e2, 697 G4double e2,
597 G4double e, 698 G4double e,
598 G4double xs1, 699 G4double xs1,
599 G4double xs2) 700 G4double xs2)
600 { 701 {
601 G4double value=0.; 702 G4double value=0.;
602 if(e1!=0){ 703 if(e1!=0){
603 G4double a = std::log10(e) - std::log10(e1 704 G4double a = std::log10(e) - std::log10(e1);
604 G4double b = std::log10(e2) - std::log10(e) 705 G4double b = std::log10(e2) - std::log10(e);
605 value = xs1 + a/(a+b)*(xs2-xs1); 706 value = xs1 + a/(a+b)*(xs2-xs1);
606 } 707 }
607 else{ 708 else{
608 G4double d1 = xs1; 709 G4double d1 = xs1;
609 G4double d2 = xs2; 710 G4double d2 = xs2;
610 value = (d1 + (d2 - d1) * (e - e1) / (e2 - 711 value = (d1 + (d2 - d1) * (e - e1) / (e2 - e1));
611 } 712 }
612 713
613 return value; 714 return value;
614 } 715 }
615 716
616 //....oooOO0OOooo........oooOO0OOooo........oo 717 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
617 718
618 G4double G4DNAELSEPAElasticModel::LinLogInterp 719 G4double G4DNAELSEPAElasticModel::LinLogInterpolate(G4double e1,
619 G4double e2, 720 G4double e2,
620 G4double e, 721 G4double e,
621 G4double xs1, 722 G4double xs1,
622 G4double xs2) 723 G4double xs2)
623 { 724 {
624 G4double d1 = std::log10(xs1); 725 G4double d1 = std::log10(xs1);
625 G4double d2 = std::log10(xs2); 726 G4double d2 = std::log10(xs2);
626 G4double value = std::pow(10,(d1 + (d2 - d1) 727 G4double value = std::pow(10,(d1 + (d2 - d1) * (e - e1) / (e2 - e1)));
627 return value; 728 return value;
628 } 729 }
629 730
630 //....oooOO0OOooo........oooOO0OOooo........oo 731 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
631 732
632 G4double G4DNAELSEPAElasticModel::LinLinInterp 733 G4double G4DNAELSEPAElasticModel::LinLinInterpolate(G4double e1,
633 G4double e2, 734 G4double e2,
634 G4double e, 735 G4double e,
635 G4double xs1, 736 G4double xs1,
636 G4double xs2) 737 G4double xs2)
637 { 738 {
638 G4double d1 = xs1; 739 G4double d1 = xs1;
639 G4double d2 = xs2; 740 G4double d2 = xs2;
640 G4double value = (d1 + (d2 - d1) * (e - e1) / 741 G4double value = (d1 + (d2 - d1) * (e - e1) / (e2 - e1));
641 return value; 742 return value;
642 } 743 }
643 744
644 //....oooOO0OOooo........oooOO0OOooo........oo 745 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
645 746
646 G4double G4DNAELSEPAElasticModel::LogLogInterp 747 G4double G4DNAELSEPAElasticModel::LogLogInterpolate(G4double e1,
647 G4double e2, 748 G4double e2,
648 G4double e, 749 G4double e,
649 G4double xs1, 750 G4double xs1,
650 G4double xs2) 751 G4double xs2)
651 { 752 {
652 G4double a = (std::log10(xs2) - std::log10(xs 753 G4double a = (std::log10(xs2) - std::log10(xs1))
653 / (std::log10(e2) - std::log10(e1 754 / (std::log10(e2) - std::log10(e1));
654 G4double b = std::log10(xs2) - a * std::log10 755 G4double b = std::log10(xs2) - a * std::log10(e2);
655 G4double sigma = a * std::log10(e) + b; 756 G4double sigma = a * std::log10(e) + b;
656 G4double value = (std::pow(10., sigma)); 757 G4double value = (std::pow(10., sigma));
657 return value; 758 return value;
658 } 759 }
659 760
660 //....oooOO0OOooo........oooOO0OOooo........oo 761 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
661 762
662 G4double G4DNAELSEPAElasticModel::QuadInterpol 763 G4double G4DNAELSEPAElasticModel::QuadInterpolator(
663 G4double cum11, 764 G4double cum11,
664 G4double cum12, 765 G4double cum12,
665 G4double cum21, 766 G4double cum21,
666 G4double cum22, 767 G4double cum22,
667 G4double a11, 768 G4double a11,
668 G4double a12, 769 G4double a12,
669 G4double a21, 770 G4double a21,
670 G4double a22, 771 G4double a22,
671 G4double e1, 772 G4double e1,
672 G4double e2, 773 G4double e2,
673 G4double t, 774 G4double t,
674 G4double cum) 775 G4double cum)
675 { 776 {
676 G4double value=0; 777 G4double value=0;
677 G4double interpolatedvalue1=0; 778 G4double interpolatedvalue1=0;
678 G4double interpolatedvalue2=0; 779 G4double interpolatedvalue2=0;
679 780
680 if(cum11!=0){ 781 if(cum11!=0){
681 interpolatedvalue1 = LinLogInterpolate(cu 782 interpolatedvalue1 = LinLogInterpolate(cum11, cum12, cum, a11, a12);
682 } 783 }
683 else{ 784 else{
684 interpolatedvalue1 = LinLinInterpolate(cu 785 interpolatedvalue1 = LinLinInterpolate(cum11, cum12, cum, a11, a12);
685 } 786 }
686 if(cum21!=0){ 787 if(cum21!=0){
687 interpolatedvalue2 = LinLogInterpolate(cu 788 interpolatedvalue2 = LinLogInterpolate(cum21, cum22, cum, a21, a22);
688 } 789 }
689 else{ 790 else{
690 interpolatedvalue2 = LinLinInterpolate(cu 791 interpolatedvalue2 = LinLinInterpolate(cum21, cum22, cum, a21, a22);
691 } 792 }
692 793
693 value = LogLinInterpolate(e1,e2,t,interpola 794 value = LogLinInterpolate(e1,e2,t,interpolatedvalue1,interpolatedvalue2);
694 795
695 return value; 796 return value;
696 } 797 }
697 798
698 //....oooOO0OOooo........oooOO0OOooo........oo 799 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
699 800
700 G4double G4DNAELSEPAElasticModel::RandomizeCos 801 G4double G4DNAELSEPAElasticModel::RandomizeCosTheta(G4int Z, G4double k)
701 { 802 {
702 803
703 G4double integrdiff = 0.; 804 G4double integrdiff = 0.;
704 G4double uniformRand = G4UniformRand(); 805 G4double uniformRand = G4UniformRand();
705 integrdiff = uniformRand; 806 integrdiff = uniformRand;
706 807
707 G4double theta = 0.; 808 G4double theta = 0.;
708 G4double cosTheta = 0.; 809 G4double cosTheta = 0.;
709 theta = Theta(Z, G4Electron::ElectronDefinit 810 theta = Theta(Z, G4Electron::ElectronDefinition(), k / eV, integrdiff);
710 811
711 cosTheta = std::cos(theta * CLHEP::pi / 180. 812 cosTheta = std::cos(theta * CLHEP::pi / 180.);
712 813
713 return cosTheta; 814 return cosTheta;
714 } 815 }
715 816
716 //....oooOO0OOooo........oooOO0OOooo........oo 817 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
717 818
718 void G4DNAELSEPAElasticModel::SetKillBelowThre 819 void G4DNAELSEPAElasticModel::SetKillBelowThreshold(G4double threshold)
719 { 820 {
720 fkillBelowEnergy_Au = threshold; 821 fkillBelowEnergy_Au = threshold;
721 822
722 if (threshold < 10 * eV) 823 if (threshold < 10 * eV)
723 { 824 {
724 G4cout<< "*** WARNING : the G4DNAELSEPAEla 825 G4cout<< "*** WARNING : the G4DNAELSEPAElasticModel model is not "
725 "defined below 10 eV !" << G4endl; 826 "defined below 10 eV !" << G4endl;
726 } 827 }
727 } 828 }
728 829