| 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 // 26 //
27 // 27 //
28 // Author: Maria Grazia Pia (Maria.Grazia.Pia@ 28 // Author: Maria Grazia Pia (Maria.Grazia.Pia@cern.ch)
29 // 29 //
30 // History: 30 // History:
31 // ----------- 31 // -----------
32 // 31 Jul 2001 MGP Created 32 // 31 Jul 2001 MGP Created
33 // 33 //
34 // ------------------------------------------- 34 // -------------------------------------------------------------------
35 35
36 #include "G4RDEMDataSet.hh" 36 #include "G4RDEMDataSet.hh"
37 #include "G4RDVDataSetAlgorithm.hh" 37 #include "G4RDVDataSetAlgorithm.hh"
38 #include <fstream> 38 #include <fstream>
39 #include <sstream> 39 #include <sstream>
40 #include "G4Integrator.hh" 40 #include "G4Integrator.hh"
41 #include "Randomize.hh" 41 #include "Randomize.hh"
42 #include "G4RDLinInterpolation.hh" 42 #include "G4RDLinInterpolation.hh"
43 43
44 44
45 G4RDEMDataSet::G4RDEMDataSet(G4int Z, 45 G4RDEMDataSet::G4RDEMDataSet(G4int Z,
46 G4RDVDataSetAlgorithm* algo, 46 G4RDVDataSetAlgorithm* algo,
47 G4double xUnit, 47 G4double xUnit,
48 G4double yUnit, 48 G4double yUnit,
49 G4bool random): 49 G4bool random):
50 z(Z), 50 z(Z),
51 energies(0), 51 energies(0),
52 data(0), 52 data(0),
53 algorithm(algo), 53 algorithm(algo),
54 unitEnergies(xUnit), 54 unitEnergies(xUnit),
55 unitData(yUnit), 55 unitData(yUnit),
56 pdf(0), 56 pdf(0),
57 randomSet(random) 57 randomSet(random)
58 { 58 {
59 if (algorithm == 0) 59 if (algorithm == 0)
60 G4Exception("G4RDEMDataSet::G4RDEMDataSet( 60 G4Exception("G4RDEMDataSet::G4RDEMDataSet()",
61 "InvalidSetup", FatalException 61 "InvalidSetup", FatalException, "Interpolation == 0!");
62 if (randomSet) BuildPdf(); 62 if (randomSet) BuildPdf();
63 } 63 }
64 64
65 G4RDEMDataSet::G4RDEMDataSet(G4int argZ, 65 G4RDEMDataSet::G4RDEMDataSet(G4int argZ,
66 G4DataVector* dataX, 66 G4DataVector* dataX,
67 G4DataVector* dataY, 67 G4DataVector* dataY,
68 G4RDVDataSetAlgorithm* algo, 68 G4RDVDataSetAlgorithm* algo,
69 G4double xUnit, 69 G4double xUnit,
70 G4double yUnit, 70 G4double yUnit,
71 G4bool random): 71 G4bool random):
72 z(argZ), 72 z(argZ),
73 energies(dataX), 73 energies(dataX),
74 data(dataY), 74 data(dataY),
75 algorithm(algo), 75 algorithm(algo),
76 unitEnergies(xUnit), 76 unitEnergies(xUnit),
77 unitData(yUnit), 77 unitData(yUnit),
78 pdf(0), 78 pdf(0),
79 randomSet(random) 79 randomSet(random)
80 { 80 {
81 if (algorithm == 0) 81 if (algorithm == 0)
82 G4Exception("G4RDEMDataSet::G4RDEMDataSet( 82 G4Exception("G4RDEMDataSet::G4RDEMDataSet()",
83 "InvalidSetup", FatalException 83 "InvalidSetup", FatalException, "Interpolation == 0!");
84 84
85 if ((energies == 0) ^ (data == 0)) 85 if ((energies == 0) ^ (data == 0))
86 G4Exception("G4RDEMDataSet::G4RDEMDataSet( 86 G4Exception("G4RDEMDataSet::G4RDEMDataSet()", "InvalidSetup",
87 FatalException, "Different size for 87 FatalException, "Different size for energies and data (zero case)!");
88 88
89 if (energies == 0) return; 89 if (energies == 0) return;
90 90
91 if (energies->size() != data->size()) 91 if (energies->size() != data->size())
92 G4Exception("G4RDEMDataSet::G4RDEMDataSet( 92 G4Exception("G4RDEMDataSet::G4RDEMDataSet()", "InvalidSetup",
93 FatalException, "Different size for 93 FatalException, "Different size for energies and data!");
94 94
95 if (randomSet) BuildPdf(); 95 if (randomSet) BuildPdf();
96 } 96 }
97 97
98 G4RDEMDataSet::~G4RDEMDataSet() 98 G4RDEMDataSet::~G4RDEMDataSet()
99 { 99 {
100 delete algorithm; 100 delete algorithm;
101 if (energies) delete energies; 101 if (energies) delete energies;
102 if (data) delete data; 102 if (data) delete data;
103 if (pdf) delete pdf; 103 if (pdf) delete pdf;
104 } 104 }
105 105
106 G4double G4RDEMDataSet::FindValue(G4double ene 106 G4double G4RDEMDataSet::FindValue(G4double energy, G4int /* componentId */) const
107 { 107 {
108 if (!energies) 108 if (!energies)
109 G4Exception("G4RDEMDataSet::FindValue()", 109 G4Exception("G4RDEMDataSet::FindValue()", "InvalidSetup",
110 FatalException, "Energies == 0 110 FatalException, "Energies == 0!");
111 if (energies->empty()) return 0; 111 if (energies->empty()) return 0;
112 if (energy <= (*energies)[0]) return (*data) 112 if (energy <= (*energies)[0]) return (*data)[0];
113 113
114 size_t i = energies->size()-1; 114 size_t i = energies->size()-1;
115 if (energy >= (*energies)[i]) return (*data) 115 if (energy >= (*energies)[i]) return (*data)[i];
116 116
117 return algorithm->Calculate(energy, FindLowe 117 return algorithm->Calculate(energy, FindLowerBound(energy), *energies, *data);
118 } 118 }
119 119
120 120
121 void G4RDEMDataSet::PrintData(void) const 121 void G4RDEMDataSet::PrintData(void) const
122 { 122 {
123 if (!energies) 123 if (!energies)
124 { 124 {
125 G4cout << "Data not available." << G4end 125 G4cout << "Data not available." << G4endl;
126 } 126 }
127 else 127 else
128 { 128 {
129 size_t size = energies->size(); 129 size_t size = energies->size();
130 for (size_t i(0); i<size; i++) 130 for (size_t i(0); i<size; i++)
131 { 131 {
132 G4cout << "Point: " << ((*energies)[i]/uni 132 G4cout << "Point: " << ((*energies)[i]/unitEnergies)
133 << " - Data value: " << ((*data)[i]/unitD 133 << " - Data value: " << ((*data)[i]/unitData);
134 if (pdf != 0) G4cout << " - PDF : " << (*p 134 if (pdf != 0) G4cout << " - PDF : " << (*pdf)[i];
135 G4cout << G4endl; 135 G4cout << G4endl;
136 } 136 }
137 } 137 }
138 } 138 }
139 139
140 140
141 void G4RDEMDataSet::SetEnergiesData(G4DataVect 141 void G4RDEMDataSet::SetEnergiesData(G4DataVector* dataX,
142 G4DataVector* dataY, 142 G4DataVector* dataY,
143 G4int /* componentId */) 143 G4int /* componentId */)
144 { 144 {
145 if (energies) delete energies; 145 if (energies) delete energies;
146 energies = dataX; 146 energies = dataX;
147 147
148 if (data) delete data; 148 if (data) delete data;
149 data = dataY; 149 data = dataY;
150 150
151 if ((energies == 0) ^ (data==0)) 151 if ((energies == 0) ^ (data==0))
152 G4Exception("G4RDEMDataSet::SetEnergiesDat 152 G4Exception("G4RDEMDataSet::SetEnergiesData()", "InvalidSetup",
153 FatalException, "Different size for ene 153 FatalException, "Different size for energies and data (zero case)!");
154 154
155 if (energies == 0) return; 155 if (energies == 0) return;
156 156
157 if (energies->size() != data->size()) 157 if (energies->size() != data->size())
158 G4Exception("G4RDEMDataSet::SetEnergiesDat 158 G4Exception("G4RDEMDataSet::SetEnergiesData()", "InvalidSetup",
159 FatalException, "Different size for ene 159 FatalException, "Different size for energies and data!");
160 } 160 }
161 161
162 G4bool G4RDEMDataSet::LoadData(const G4String& 162 G4bool G4RDEMDataSet::LoadData(const G4String& fileName)
163 { 163 {
164 // The file is organized into two columns: 164 // The file is organized into two columns:
165 // 1st column is the energy 165 // 1st column is the energy
166 // 2nd column is the corresponding value 166 // 2nd column is the corresponding value
167 // The file terminates with the pattern: -1 167 // The file terminates with the pattern: -1 -1
168 // -2 168 // -2 -2
169 169
170 G4String fullFileName(FullFileName(fileName) 170 G4String fullFileName(FullFileName(fileName));
171 std::ifstream in(fullFileName); 171 std::ifstream in(fullFileName);
172 172
173 if (!in.is_open()) 173 if (!in.is_open())
174 { 174 {
175 G4String message("Data file \""); 175 G4String message("Data file \"");
176 message += fullFileName; 176 message += fullFileName;
177 message += "\" not found"; 177 message += "\" not found";
178 G4Exception("G4RDEMDataSet::LoadData()", 178 G4Exception("G4RDEMDataSet::LoadData()", "DataNotFound",
179 FatalException, message); 179 FatalException, message);
180 } 180 }
181 181
182 G4DataVector* argEnergies=new G4DataVector; 182 G4DataVector* argEnergies=new G4DataVector;
183 G4DataVector* argData=new G4DataVector; 183 G4DataVector* argData=new G4DataVector;
184 184
185 G4double a; 185 G4double a;
186 bool energyColumn(true); 186 bool energyColumn(true);
187 187
188 do 188 do
189 { 189 {
190 in >> a; 190 in >> a;
191 191
192 if (a!=-1 && a!=-2) 192 if (a!=-1 && a!=-2)
193 { 193 {
194 if (energyColumn) 194 if (energyColumn)
195 argEnergies->push_back(a*unitEnergies); 195 argEnergies->push_back(a*unitEnergies);
196 else 196 else
197 argData->push_back(a*unitData); 197 argData->push_back(a*unitData);
198 energyColumn=(!energyColumn); 198 energyColumn=(!energyColumn);
199 } 199 }
200 } 200 }
201 while (a != -2); 201 while (a != -2);
202 202
203 SetEnergiesData(argEnergies, argData, 0); 203 SetEnergiesData(argEnergies, argData, 0);
204 if (randomSet) BuildPdf(); 204 if (randomSet) BuildPdf();
205 205
206 return true; 206 return true;
207 } 207 }
208 208
209 G4bool G4RDEMDataSet::SaveData(const G4String& 209 G4bool G4RDEMDataSet::SaveData(const G4String& name) const
210 { 210 {
211 // The file is organized into two columns: 211 // The file is organized into two columns:
212 // 1st column is the energy 212 // 1st column is the energy
213 // 2nd column is the corresponding value 213 // 2nd column is the corresponding value
214 // The file terminates with the pattern: -1 214 // The file terminates with the pattern: -1 -1
215 // -2 215 // -2 -2
216 216
217 G4String fullFileName(FullFileName(name)); 217 G4String fullFileName(FullFileName(name));
218 std::ofstream out(fullFileName); 218 std::ofstream out(fullFileName);
219 219
220 if (!out.is_open()) 220 if (!out.is_open())
221 { 221 {
222 G4String message("Cannot open \""); 222 G4String message("Cannot open \"");
223 message+=fullFileName; 223 message+=fullFileName;
224 message+="\""; 224 message+="\"";
225 G4Exception("G4RDEMDataSet::SaveData()", 225 G4Exception("G4RDEMDataSet::SaveData()", "CannotOpenFile",
226 FatalException, message); 226 FatalException, message);
227 } 227 }
228 228
229 out.precision(10); 229 out.precision(10);
230 out.width(15); 230 out.width(15);
231 out.setf(std::ofstream::left); 231 out.setf(std::ofstream::left);
232 232
233 if (energies!=0 && data!=0) 233 if (energies!=0 && data!=0)
234 { 234 {
235 G4DataVector::const_iterator i(energies- 235 G4DataVector::const_iterator i(energies->begin());
236 G4DataVector::const_iterator endI(energi 236 G4DataVector::const_iterator endI(energies->end());
237 G4DataVector::const_iterator j(data->beg 237 G4DataVector::const_iterator j(data->begin());
238 238
239 while (i!=endI) 239 while (i!=endI)
240 { 240 {
241 out.precision(10); 241 out.precision(10);
242 out.width(15); 242 out.width(15);
243 out.setf(std::ofstream::left); 243 out.setf(std::ofstream::left);
244 out << ((*i)/unitEnergies) << ' '; 244 out << ((*i)/unitEnergies) << ' ';
245 245
246 out.precision(10); 246 out.precision(10);
247 out.width(15); 247 out.width(15);
248 out.setf(std::ofstream::left); 248 out.setf(std::ofstream::left);
249 out << ((*j)/unitData) << std::endl; 249 out << ((*j)/unitData) << std::endl;
250 250
251 i++; 251 i++;
252 j++; 252 j++;
253 } 253 }
254 } 254 }
255 255
256 out.precision(10); 256 out.precision(10);
257 out.width(15); 257 out.width(15);
258 out.setf(std::ofstream::left); 258 out.setf(std::ofstream::left);
259 out << -1.f << ' '; 259 out << -1.f << ' ';
260 260
261 out.precision(10); 261 out.precision(10);
262 out.width(15); 262 out.width(15);
263 out.setf(std::ofstream::left); 263 out.setf(std::ofstream::left);
264 out << -1.f << std::endl; 264 out << -1.f << std::endl;
265 265
266 out.precision(10); 266 out.precision(10);
267 out.width(15); 267 out.width(15);
268 out.setf(std::ofstream::left); 268 out.setf(std::ofstream::left);
269 out << -2.f << ' '; 269 out << -2.f << ' ';
270 270
271 out.precision(10); 271 out.precision(10);
272 out.width(15); 272 out.width(15);
273 out.setf(std::ofstream::left); 273 out.setf(std::ofstream::left);
274 out << -2.f << std::endl; 274 out << -2.f << std::endl;
275 275
276 return true; 276 return true;
277 } 277 }
278 278
279 size_t G4RDEMDataSet::FindLowerBound(G4double 279 size_t G4RDEMDataSet::FindLowerBound(G4double x) const
280 { 280 {
281 size_t lowerBound = 0; 281 size_t lowerBound = 0;
282 size_t upperBound(energies->size() - 1); 282 size_t upperBound(energies->size() - 1);
283 283
284 while (lowerBound <= upperBound) 284 while (lowerBound <= upperBound)
285 { 285 {
286 size_t midBin((lowerBound + upperBound) 286 size_t midBin((lowerBound + upperBound) / 2);
287 287
288 if (x < (*energies)[midBin]) upperBound 288 if (x < (*energies)[midBin]) upperBound = midBin - 1;
289 else lowerBound = midBin + 1; 289 else lowerBound = midBin + 1;
290 } 290 }
291 291
292 return upperBound; 292 return upperBound;
293 } 293 }
294 294
295 295
296 size_t G4RDEMDataSet::FindLowerBound(G4double 296 size_t G4RDEMDataSet::FindLowerBound(G4double x, G4DataVector* values) const
297 { 297 {
298 size_t lowerBound = 0;; 298 size_t lowerBound = 0;;
299 size_t upperBound(values->size() - 1); 299 size_t upperBound(values->size() - 1);
300 300
301 while (lowerBound <= upperBound) 301 while (lowerBound <= upperBound)
302 { 302 {
303 size_t midBin((lowerBound + upperBound) 303 size_t midBin((lowerBound + upperBound) / 2);
304 304
305 if (x < (*values)[midBin]) upperBound = 305 if (x < (*values)[midBin]) upperBound = midBin - 1;
306 else lowerBound = midBin + 1; 306 else lowerBound = midBin + 1;
307 } 307 }
308 308
309 return upperBound; 309 return upperBound;
310 } 310 }
311 311
312 312
313 G4String G4RDEMDataSet::FullFileName(const G4S 313 G4String G4RDEMDataSet::FullFileName(const G4String& name) const
314 { 314 {
315 const char* path = G4FindDataDir("G4LEDATA") 315 const char* path = G4FindDataDir("G4LEDATA");
316 if (!path) 316 if (!path)
317 G4Exception("G4RDEMDataSet::FullFileName() 317 G4Exception("G4RDEMDataSet::FullFileName()", "InvalidSetup",
318 FatalException, "G4LEDATA envi 318 FatalException, "G4LEDATA environment variable not set!");
319 319
320 std::ostringstream fullFileName; 320 std::ostringstream fullFileName;
321 fullFileName << path << '/' << name << z << 321 fullFileName << path << '/' << name << z << ".dat";
322 322
323 return G4String(fullFileName.str().c_str()); 323 return G4String(fullFileName.str().c_str());
324 } 324 }
325 325
326 326
327 void G4RDEMDataSet::BuildPdf() 327 void G4RDEMDataSet::BuildPdf()
328 { 328 {
329 pdf = new G4DataVector; 329 pdf = new G4DataVector;
330 G4Integrator <G4RDEMDataSet, G4double(G4RDEM 330 G4Integrator <G4RDEMDataSet, G4double(G4RDEMDataSet::*)(G4double)> integrator;
331 331
332 G4int nData = data->size(); 332 G4int nData = data->size();
333 pdf->push_back(0.); 333 pdf->push_back(0.);
334 334
335 // Integrate the data distribution 335 // Integrate the data distribution
336 G4int i; 336 G4int i;
337 G4double totalSum = 0.; 337 G4double totalSum = 0.;
338 for (i=1; i<nData; i++) 338 for (i=1; i<nData; i++)
339 { 339 {
340 G4double xLow = (*energies)[i-1]; 340 G4double xLow = (*energies)[i-1];
341 G4double xHigh = (*energies)[i]; 341 G4double xHigh = (*energies)[i];
342 G4double sum = integrator.Legendre96(thi 342 G4double sum = integrator.Legendre96(this, &G4RDEMDataSet::IntegrationFunction, xLow, xHigh);
343 totalSum = totalSum + sum; 343 totalSum = totalSum + sum;
344 pdf->push_back(totalSum); 344 pdf->push_back(totalSum);
345 } 345 }
346 346
347 // Normalize to the last bin 347 // Normalize to the last bin
348 G4double tot = 0.; 348 G4double tot = 0.;
349 if (totalSum > 0.) tot = 1. / totalSum; 349 if (totalSum > 0.) tot = 1. / totalSum;
350 for (i=1; i<nData; i++) 350 for (i=1; i<nData; i++)
351 { 351 {
352 (*pdf)[i] = (*pdf)[i] * tot; 352 (*pdf)[i] = (*pdf)[i] * tot;
353 } 353 }
354 } 354 }
355 355
356 356
357 G4double G4RDEMDataSet::RandomSelect(G4int /* 357 G4double G4RDEMDataSet::RandomSelect(G4int /* componentId */) const
358 { 358 {
359 // Random select a X value according to the 359 // Random select a X value according to the cumulative probability distribution
360 // derived from the data 360 // derived from the data
361 361
362 if (!pdf) G4Exception("G4RDEMDataSet::Random 362 if (!pdf) G4Exception("G4RDEMDataSet::RandomSelect()", "InvalidSetup",
363 FatalException, "PDF has not been c 363 FatalException, "PDF has not been created for this data set");
364 364
365 G4double value = 0.; 365 G4double value = 0.;
366 G4double x = G4UniformRand(); 366 G4double x = G4UniformRand();
367 367
368 // Locate the random value in the X vector b 368 // Locate the random value in the X vector based on the PDF
369 size_t bin = FindLowerBound(x,pdf); 369 size_t bin = FindLowerBound(x,pdf);
370 370
371 // Interpolate the PDF to calculate the X va 371 // Interpolate the PDF to calculate the X value:
372 // linear interpolation in the first bin (to 372 // linear interpolation in the first bin (to avoid problem with 0),
373 // interpolation with associated data set al 373 // interpolation with associated data set algorithm in other bins
374 374
375 G4RDLinInterpolation linearAlgo; 375 G4RDLinInterpolation linearAlgo;
376 if (bin == 0) value = linearAlgo.Calculate(x 376 if (bin == 0) value = linearAlgo.Calculate(x, bin, *pdf, *energies);
377 else value = algorithm->Calculate(x, bin, *p 377 else value = algorithm->Calculate(x, bin, *pdf, *energies);
378 378
379 // G4cout << x << " random bin "<< bin << " 379 // G4cout << x << " random bin "<< bin << " - " << value << G4endl;
380 return value; 380 return value;
381 } 381 }
382 382
383 G4double G4RDEMDataSet::IntegrationFunction(G4 383 G4double G4RDEMDataSet::IntegrationFunction(G4double x)
384 { 384 {
385 // This function is needed by G4Integrator t 385 // This function is needed by G4Integrator to calculate the integral of the data distribution
386 386
387 G4double y = 0; 387 G4double y = 0;
388 388
389 // Locate the random value in the X vector ba 389 // Locate the random value in the X vector based on the PDF
390 size_t bin = FindLowerBound(x); 390 size_t bin = FindLowerBound(x);
391 391
392 // Interpolate to calculate the X value: 392 // Interpolate to calculate the X value:
393 // linear interpolation in the first bin (to 393 // linear interpolation in the first bin (to avoid problem with 0),
394 // interpolation with associated algorithm i 394 // interpolation with associated algorithm in other bins
395 395
396 G4RDLinInterpolation linearAlgo; 396 G4RDLinInterpolation linearAlgo;
397 397
398 if (bin == 0) y = linearAlgo.Calculate(x, bi 398 if (bin == 0) y = linearAlgo.Calculate(x, bin, *energies, *data);
399 else y = algorithm->Calculate(x, bin, *energ 399 else y = algorithm->Calculate(x, bin, *energies, *data);
400 400
401 return y; 401 return y;
402 } 402 }
403 403
404 404