Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/processes/electromagnetic/lowenergy/src/G4EMDataSet.cc

Version: [ ReleaseNotes ] [ 1.0 ] [ 1.1 ] [ 2.0 ] [ 3.0 ] [ 3.1 ] [ 3.2 ] [ 4.0 ] [ 4.0.p1 ] [ 4.0.p2 ] [ 4.1 ] [ 4.1.p1 ] [ 5.0 ] [ 5.0.p1 ] [ 5.1 ] [ 5.1.p1 ] [ 5.2 ] [ 5.2.p1 ] [ 5.2.p2 ] [ 6.0 ] [ 6.0.p1 ] [ 6.1 ] [ 6.2 ] [ 6.2.p1 ] [ 6.2.p2 ] [ 7.0 ] [ 7.0.p1 ] [ 7.1 ] [ 7.1.p1 ] [ 8.0 ] [ 8.0.p1 ] [ 8.1 ] [ 8.1.p1 ] [ 8.1.p2 ] [ 8.2 ] [ 8.2.p1 ] [ 8.3 ] [ 8.3.p1 ] [ 8.3.p2 ] [ 9.0 ] [ 9.0.p1 ] [ 9.0.p2 ] [ 9.1 ] [ 9.1.p1 ] [ 9.1.p2 ] [ 9.1.p3 ] [ 9.2 ] [ 9.2.p1 ] [ 9.2.p2 ] [ 9.2.p3 ] [ 9.2.p4 ] [ 9.3 ] [ 9.3.p1 ] [ 9.3.p2 ] [ 9.4 ] [ 9.4.p1 ] [ 9.4.p2 ] [ 9.4.p3 ] [ 9.4.p4 ] [ 9.5 ] [ 9.5.p1 ] [ 9.5.p2 ] [ 9.6 ] [ 9.6.p1 ] [ 9.6.p2 ] [ 9.6.p3 ] [ 9.6.p4 ] [ 10.0 ] [ 10.0.p1 ] [ 10.0.p2 ] [ 10.0.p3 ] [ 10.0.p4 ] [ 10.1 ] [ 10.1.p1 ] [ 10.1.p2 ] [ 10.1.p3 ] [ 10.2 ] [ 10.2.p1 ] [ 10.2.p2 ] [ 10.2.p3 ] [ 10.3 ] [ 10.3.p1 ] [ 10.3.p2 ] [ 10.3.p3 ] [ 10.4 ] [ 10.4.p1 ] [ 10.4.p2 ] [ 10.4.p3 ] [ 10.5 ] [ 10.5.p1 ] [ 10.6 ] [ 10.6.p1 ] [ 10.6.p2 ] [ 10.6.p3 ] [ 10.7 ] [ 10.7.p1 ] [ 10.7.p2 ] [ 10.7.p3 ] [ 10.7.p4 ] [ 11.0 ] [ 11.0.p1 ] [ 11.0.p2 ] [ 11.0.p3, ] [ 11.0.p4 ] [ 11.1 ] [ 11.1.1 ] [ 11.1.2 ] [ 11.1.3 ] [ 11.2 ] [ 11.2.1 ] [ 11.2.2 ] [ 11.3.0 ]

Diff markup

Differences between /processes/electromagnetic/lowenergy/src/G4EMDataSet.cc (Version 11.3.0) and /processes/electromagnetic/lowenergy/src/G4EMDataSet.cc (Version 6.2.p1)


  1 //                                                  1 //
  2 // *******************************************      2 // ********************************************************************
  3 // * License and Disclaimer                    <<   3 // * DISCLAIMER                                                       *
  4 // *                                                4 // *                                                                  *
  5 // * The  Geant4 software  is  copyright of th <<   5 // * The following disclaimer summarizes all the specific disclaimers *
  6 // * the Geant4 Collaboration.  It is provided <<   6 // * of contributors to this software. The specific disclaimers,which *
  7 // * conditions of the Geant4 Software License <<   7 // * govern, are listed with their locations in:                      *
  8 // * LICENSE and available at  http://cern.ch/ <<   8 // *   http://cern.ch/geant4/license                                  *
  9 // * include a list of copyright holders.      << 
 10 // *                                                9 // *                                                                  *
 11 // * Neither the authors of this software syst     10 // * Neither the authors of this software system, nor their employing *
 12 // * institutes,nor the agencies providing fin     11 // * institutes,nor the agencies providing financial support for this *
 13 // * work  make  any representation or  warran     12 // * work  make  any representation or  warranty, express or implied, *
 14 // * regarding  this  software system or assum     13 // * regarding  this  software system or assume any liability for its *
 15 // * use.  Please see the license in the file  <<  14 // * use.                                                             *
 16 // * for the full disclaimer and the limitatio << 
 17 // *                                               15 // *                                                                  *
 18 // * This  code  implementation is the result  <<  16 // * This  code  implementation is the  intellectual property  of the *
 19 // * technical work of the GEANT4 collaboratio <<  17 // * GEANT4 collaboration.                                            *
 20 // * By using,  copying,  modifying or  distri <<  18 // * By copying,  distributing  or modifying the Program (or any work *
 21 // * any work based  on the software)  you  ag <<  19 // * based  on  the Program)  you indicate  your  acceptance of  this *
 22 // * use  in  resulting  scientific  publicati <<  20 // * statement, and all its terms.                                    *
 23 // * acceptance of all terms of the Geant4 Sof << 
 24 // *******************************************     21 // ********************************************************************
 25 //                                                 22 //
 26 // History:                                    << 
 27 // -----------                                 << 
 28 // 31 Jul 2001   MGP        Created            << 
 29 //                                             << 
 30 // 15 Jul 2009   Nicolas A. Karakatsanis       << 
 31 //                                             << 
 32 //                           - New Constructor << 
 33 //                             to enhance CPU  << 
 34 //                                                 23 //
 35 //                           - Destructor was  <<  24 // $Id: G4EMDataSet.cc,v 1.9 2003/06/16 17:00:07 gunter Exp $
                                                   >>  25 // GEANT4 tag $Name: geant4-05-02-patch-01 $
 36 //                                                 26 //
 37 //                           - LoadNonLogData  <<  27 // Author: Maria Grazia Pia (Maria.Grazia.Pia@cern.ch)
 38 //                             dataset. It is  << 
 39 //                                                 28 //
 40 //                           - LoadData method <<  29 // History:
 41 //                             It retrieves th <<  30 // -----------
 42 //                             respective log  <<  31 // 31 Jul 2001   MGP        Created
 43 //                                             << 
 44 //                           - SetLogEnergiesD << 
 45 //                             The EM data set << 
 46 //                             These initializ << 
 47 //                             operations      << 
 48 //                                             << 
 49 // 26 Dec 2010 V.Ivanchenko Fixed Coverity war << 
 50 //                                                 32 //
 51 // -------------------------------------------     33 // -------------------------------------------------------------------
 52                                                    34 
 53 #include "G4EMDataSet.hh"                          35 #include "G4EMDataSet.hh"
 54 #include "G4VDataSetAlgorithm.hh"                  36 #include "G4VDataSetAlgorithm.hh"
 55 #include <fstream>                                 37 #include <fstream>
 56 #include <sstream>                             <<  38 #include <strstream>
 57 #include "G4Integrator.hh"                     << 
 58 #include "Randomize.hh"                        << 
 59 #include "G4LinInterpolation.hh"               << 
 60                                                << 
 61 G4EMDataSet::G4EMDataSet(G4int Z,              << 
 62        G4VDataSetAlgorithm* algo,              << 
 63        G4double xUnit,                         << 
 64        G4double yUnit,                         << 
 65        G4bool random):                         << 
 66   energies(nullptr),                           << 
 67   data(nullptr),                               << 
 68   log_energies(nullptr),                       << 
 69   log_data(nullptr),                           << 
 70   algorithm(algo),                             << 
 71   pdf(nullptr),                                << 
 72   unitEnergies(xUnit),                         << 
 73   unitData(yUnit),                             << 
 74   z(Z),                                        << 
 75   randomSet(random)                            << 
 76 {                                              << 
 77   if (algorithm == nullptr) {                  << 
 78     G4Exception("G4EMDataSet::G4EMDataSet",    << 
 79         "em1012",FatalException,"interpolation << 
 80   } else if (randomSet) { BuildPdf(); }        << 
 81 }                                              << 
 82                                                    39 
 83 //....oooOO0OOooo........oooOO0OOooo........oo <<  40 // Constructor
 84                                                    41 
 85 G4EMDataSet::G4EMDataSet(G4int argZ,           <<  42 G4EMDataSet::G4EMDataSet(G4int Z,
 86        G4DataVector* dataX,                    <<  43        G4DataVector* points, 
 87        G4DataVector* dataY,                    <<  44        G4DataVector* values,
 88        G4VDataSetAlgorithm* algo,              <<  45        const G4VDataSetAlgorithm* interpolation,
 89        G4double xUnit,                         <<  46        G4double unitE, G4double unitData)
 90        G4double yUnit,                         <<  47   :z(Z), energies(points), data(values), algorithm(interpolation)
 91        G4bool random):                         <<  48 {
 92   energies(dataX),                             <<  49   numberOfBins = energies->size();
 93   data(dataY),                                 <<  50   unit1 = unitE;
 94   log_energies(nullptr),                       <<  51   unit2 = unitData;
 95   log_data(nullptr),                           <<  52   if (interpolation == 0) 
 96   algorithm(algo),                             <<  53     G4Exception("G4EMDataSet::G4EMDataSet - interpolation algorithm = 0");
 97   pdf(nullptr),                                <<  54 }
 98   unitEnergies(xUnit),                         <<  55 
 99   unitData(yUnit),                             <<  56 G4EMDataSet:: G4EMDataSet(G4int Z, 
100   z(argZ),                                     <<  57         const G4String& dataFile,
101   randomSet(random)                            <<  58         const G4VDataSetAlgorithm* interpolation,
                                                   >>  59         G4double unitE, G4double unitData)
                                                   >>  60       :z(Z), algorithm(interpolation)
102 {                                                  61 {
103   if (!algorithm || !energies || !data) {      <<  62   energies = new G4DataVector;
104     G4Exception("G4EMDataSet::G4EMDataSet",    <<  63   data = new G4DataVector;
105         "em1012",FatalException,"interpolation <<  64   unit1 = unitE;
106   } else {                                     <<  65   unit2 = unitData;  
107     if (energies->size() != data->size()) {    <<  66   LoadData(dataFile);
108       G4Exception("G4EMDataSet::G4EMDataSet",  <<  67   numberOfBins = energies->size();
109       "em1012",FatalException,"different size  <<  68   if (interpolation == 0) 
110     } else if (randomSet) {                    <<  69     G4Exception("G4EMDataSet::G4EMDataSet - interpolation algorithm = 0");
111       BuildPdf();                              << 
112     }                                          << 
113   }                                            << 
114 }                                              << 
115                                                << 
116 //....oooOO0OOooo........oooOO0OOooo........oo << 
117                                                    70 
118 G4EMDataSet::G4EMDataSet(G4int argZ,           << 
119        G4DataVector* dataX,                    << 
120        G4DataVector* dataY,                    << 
121                          G4DataVector* dataLog << 
122                          G4DataVector* dataLog << 
123        G4VDataSetAlgorithm* algo,              << 
124        G4double xUnit,                         << 
125        G4double yUnit,                         << 
126        G4bool random):                         << 
127   energies(dataX),                             << 
128   data(dataY),                                 << 
129   log_energies(dataLogX),                      << 
130   log_data(dataLogY),                          << 
131   algorithm(algo),                             << 
132   pdf(nullptr),                                << 
133   unitEnergies(xUnit),                         << 
134   unitData(yUnit),                             << 
135   z(argZ),                                     << 
136   randomSet(random)                            << 
137 {                                              << 
138   if (!algorithm || !energies || !data || !log << 
139     G4Exception("G4EMDataSet::G4EMDataSet",    << 
140         "em1012",FatalException,"interpolation << 
141   } else {                                     << 
142     if ((energies->size() != data->size()) ||  << 
143   (energies->size() != log_energies->size()) | << 
144   (energies->size() != log_data->size())) {    << 
145       G4Exception("G4EMDataSet::G4EMDataSet",  << 
146       "em1012",FatalException,"different size  << 
147     } else if (randomSet) {                    << 
148       BuildPdf();                              << 
149     }                                          << 
150   }                                            << 
151 }                                                  71 }
152                                                    72 
153 //....oooOO0OOooo........oooOO0OOooo........oo <<  73 // Destructor
154                                                    74 
155 G4EMDataSet::~G4EMDataSet()                        75 G4EMDataSet::~G4EMDataSet()
156 {                                                  76 { 
157   delete algorithm;                                77   delete algorithm;
158   delete energies;                             <<  78   delete energies;
159   delete data;                                 <<  79   delete data;
160   delete pdf;                                  << 
161   delete log_energies;                         << 
162   delete log_data;                             << 
163 }                                              << 
164                                                << 
165 //....oooOO0OOooo........oooOO0OOooo........oo << 
166                                                << 
167 G4double G4EMDataSet::FindValue(G4double energ << 
168 {                                              << 
169   if (energy <= (*energies)[0]) {              << 
170     return (*data)[0];                         << 
171   }                                            << 
172   std::size_t i = energies->size()-1;          << 
173   if (energy >= (*energies)[i]) { return (*dat << 
174                                                << 
175   //Nicolas A. Karakatsanis: Check if the loga << 
176   //                         which Interpolati << 
177   if (log_energies != nullptr)                 << 
178    {                                           << 
179      return algorithm->Calculate(energy, (G4in << 
180    }                                           << 
181                                                << 
182   return algorithm->Calculate(energy, (G4int)F << 
183 }                                              << 
184                                                << 
185 //....oooOO0OOooo........oooOO0OOooo........oo << 
186                                                << 
187 void G4EMDataSet::PrintData(void) const        << 
188 {                                              << 
189   std::size_t size = energies->size();         << 
190   for (std::size_t i(0); i<size; ++i)          << 
191     {                                          << 
192       G4cout << "Point: " << ((*energies)[i]/u << 
193        << " - Data value: " << ((*data)[i]/uni << 
194       if (pdf != 0) G4cout << " - PDF : " << ( << 
195       G4cout << G4endl;                        << 
196     }                                          << 
197 }                                              << 
198                                                << 
199 //....oooOO0OOooo........oooOO0OOooo........oo << 
200                                                << 
201 void G4EMDataSet::SetEnergiesData(G4DataVector << 
202           G4DataVector* dataY,                 << 
203           G4int /* componentId */)             << 
204 {                                              << 
205   if(!dataX || !dataY) {                       << 
206     G4Exception("G4EMDataSet::SetEnergiesData" << 
207         "em1012",FatalException,"new interpola << 
208   } else {                                     << 
209     if (dataX->size() != dataY->size()) {      << 
210       G4Exception("G4EMDataSet::SetEnergiesDat << 
211       "em1012",FatalException,"different size  << 
212     } else {                                   << 
213                                                << 
214       delete energies;                         << 
215       energies = dataX;                        << 
216                                                << 
217       delete data;                             << 
218       data = dataY;                            << 
219       //G4cout << "Size of energies: " << ener << 
220       //<< "Size of data: " << data->size() << << 
221     }                                          << 
222   }                                            << 
223 }                                              << 
224                                                << 
225 //....oooOO0OOooo........oooOO0OOooo........oo << 
226                                                << 
227 void G4EMDataSet::SetLogEnergiesData(G4DataVec << 
228                                      G4DataVec << 
229                                      G4DataVec << 
230              G4DataVector* data_logY,          << 
231                                      G4int /*  << 
232 {                                              << 
233   if(!dataX || !dataY || !data_logX || !data_l << 
234     G4Exception("G4EMDataSet::SetEnergiesData" << 
235         "em1012",FatalException,"new interpola << 
236   } else {                                     << 
237     if (dataX->size() != dataY->size() ||      << 
238   dataX->size() != data_logX->size() ||        << 
239   dataX->size() != data_logY->size()) {        << 
240       G4Exception("G4EMDataSet::SetEnergiesDat << 
241       "em1012",FatalException,"different size  << 
242     } else {                                   << 
243                                                << 
244       delete energies;                         << 
245       energies = dataX;                        << 
246                                                << 
247       delete data;                             << 
248       data = dataY;                            << 
249                                                << 
250       delete log_energies;                     << 
251       log_energies = data_logX;                << 
252                                                << 
253       delete log_data;                         << 
254       log_data = data_logY;                    << 
255       //G4cout << "Size of energies: " << ener << 
256       //<< "Size of data: " << data->size() << << 
257     }                                          << 
258   }                                            << 
259 }                                                  80 }
260                                                    81 
261 //....oooOO0OOooo........oooOO0OOooo........oo << 
262                                                    82 
263 G4bool G4EMDataSet::LoadData(const G4String& f <<  83 G4double G4EMDataSet::FindValue(G4double e, G4int) const
264 {                                                  84 {
265  // The file is organized into four columns:   <<  85   G4double value = 0.;
266  // 1st column contains the values of energy   <<  86   if ( !(energies->empty()) )
267  // 2nd column contains the corresponding data << 
268  // The file terminates with the pattern: -1   << 
269  //                                       -2   << 
270                                                << 
271   G4String fullFileName(FullFileName(fileName) << 
272   std::ifstream in(fullFileName);              << 
273                                                << 
274   if (!in.is_open())                           << 
275     {                                          << 
276       G4String message("data file \"");        << 
277       message += fullFileName;                 << 
278       message += "\" not found";               << 
279       G4Exception("G4EMDataSet::LoadData",     << 
280         "em1012",FatalException,message);      << 
281       return false;                            << 
282     }                                          << 
283                                                << 
284   delete energies;                             << 
285   delete data;                                 << 
286   delete log_energies;                         << 
287   delete log_data;                             << 
288   energies = new G4DataVector;                 << 
289   data = new G4DataVector;                     << 
290   log_energies = new G4DataVector;             << 
291   log_data = new G4DataVector;                 << 
292                                                << 
293   G4double a, b;                               << 
294   do                                           << 
295     {                                              87     {
296       in >> a >> b;                            <<  88       G4double e0 = (*energies)[0];
297                                                <<  89       // Protections
298       if (a != -1 && a != -2)                  <<  90       size_t bin = FindBinLocation(e);
                                                   >>  91       if (bin == numberOfBins)
299   {                                                92   {
300     if (a==0.) { a=1e-300; }                   <<  93     //      G4cout << "WARNING - G4EMDataSet::FindValue: energy outside upper boundary"
301     if (b==0.) { b=1e-300; }                   <<  94     //     << G4endl;
302           a *= unitEnergies;                   <<  95     value = (*data)[bin];
303           b *= unitData;                       << 
304     energies->push_back(a);                    << 
305     log_energies->push_back(std::log10(a));    << 
306     data->push_back(b);                        << 
307     log_data->push_back(std::log10(b));        << 
308   }                                                96   }
309     }                                          <<  97       else if (e <= e0)
310   while (a != -2);                             << 
311                                                << 
312   if (randomSet) { BuildPdf(); }               << 
313                                                << 
314   return true;                                 << 
315 }                                              << 
316                                                << 
317 //....oooOO0OOooo........oooOO0OOooo........oo << 
318                                                << 
319 G4bool G4EMDataSet::LoadNonLogData(const G4Str << 
320 {                                              << 
321  // The file is organized into four columns:   << 
322  // 1st column contains the values of energy   << 
323  // 2nd column contains the corresponding data << 
324  // The file terminates with the pattern: -1   << 
325  //                                       -2   << 
326                                                << 
327   G4String fullFileName(FullFileName(fileName) << 
328   std::ifstream in(fullFileName);              << 
329   if (!in.is_open())                           << 
330     {                                          << 
331       G4String message("data file \"");        << 
332       message += fullFileName;                 << 
333       message += "\" not found";               << 
334       G4Exception("G4EMDataSet::LoadNonLogData << 
335         "em1012",FatalException,message);      << 
336     }                                          << 
337                                                << 
338   G4DataVector* argEnergies=new G4DataVector;  << 
339   G4DataVector* argData=new G4DataVector;      << 
340                                                << 
341   G4double a;                                  << 
342   G4int k = 0;                                 << 
343   G4int nColumns = 2;                          << 
344                                                << 
345   do                                           << 
346     {                                          << 
347       in >> a;                                 << 
348                                                << 
349       if (a != -1 && a != -2)                  << 
350   {                                                98   {
351     if (k%nColumns == 0)                       <<  99     //     G4cout << "WARNING - G4EMDataSet::FindValue: energy outside lower boundary"
352             {                                  << 100     //     << G4endl;
353        argEnergies->push_back(a*unitEnergies); << 101     value = (*data)[0];
354             }                                  << 
355     else if (k%nColumns == 1)                  << 
356             {                                  << 
357        argData->push_back(a*unitData);         << 
358             }                                  << 
359           k++;                                 << 
360   }                                               102   }
361     }                                          << 103       else
362   while (a != -2);                             << 
363                                                << 
364   SetEnergiesData(argEnergies, argData, 0);    << 
365                                                << 
366   if (randomSet) BuildPdf();                   << 
367                                                << 
368   return true;                                 << 
369 }                                              << 
370                                                << 
371 //....oooOO0OOooo........oooOO0OOooo........oo << 
372                                                << 
373 G4bool G4EMDataSet::SaveData(const G4String& n << 
374 {                                              << 
375   // The file is organized into two columns:   << 
376   // 1st column is the energy                  << 
377   // 2nd column is the corresponding value     << 
378   // The file terminates with the pattern: -1  << 
379   //                                       -2  << 
380                                                << 
381   G4String fullFileName(FullFileName(name));   << 
382   std::ofstream out(fullFileName);             << 
383                                                << 
384   if (!out.is_open())                          << 
385     {                                          << 
386       G4String message("cannot open \"");      << 
387       message+=fullFileName;                   << 
388       message+="\"";                           << 
389       G4Exception("G4EMDataSet::SaveData",     << 
390         "em1012",FatalException,message);      << 
391     }                                          << 
392                                                << 
393   out.precision(10);                           << 
394   out.width(15);                               << 
395   out.setf(std::ofstream::left);               << 
396                                                << 
397   if (energies!=0 && data!=0)                  << 
398     {                                          << 
399       G4DataVector::const_iterator i(energies- << 
400       G4DataVector::const_iterator endI(energi << 
401       G4DataVector::const_iterator j(data->beg << 
402                                                << 
403       while (i!=endI)                          << 
404   {                                               104   {
405     out.precision(10);                         << 105     if (algorithm == 0) G4Exception("G4EMDataSet::FindValue - interpolation algorithm = 0");
406     out.width(15);                             << 106     value = algorithm->Calculate(e,bin,*energies,*data);
407     out.setf(std::ofstream::left);             << 
408     out << ((*i)/unitEnergies) << ' ';         << 
409                                                << 
410     out.precision(10);                         << 
411     out.width(15);                             << 
412     out.setf(std::ofstream::left);             << 
413     out << ((*j)/unitData) << std::endl;       << 
414                                                << 
415     i++;                                       << 
416     j++;                                       << 
417   }                                               107   }
418     }                                             108     }
419                                                << 109   return value;
420   out.precision(10);                           << 
421   out.width(15);                               << 
422   out.setf(std::ofstream::left);               << 
423   out << -1.f << ' ';                          << 
424                                                << 
425   out.precision(10);                           << 
426   out.width(15);                               << 
427   out.setf(std::ofstream::left);               << 
428   out << -1.f << std::endl;                    << 
429                                                << 
430   out.precision(10);                           << 
431   out.width(15);                               << 
432   out.setf(std::ofstream::left);               << 
433   out << -2.f << ' ';                          << 
434                                                << 
435   out.precision(10);                           << 
436   out.width(15);                               << 
437   out.setf(std::ofstream::left);               << 
438   out << -2.f << std::endl;                    << 
439                                                << 
440   return true;                                 << 
441 }                                                 110 }
442                                                   111 
443 //....oooOO0OOooo........oooOO0OOooo........oo << 112 G4int G4EMDataSet::FindBinLocation(G4double energy) const
444                                                << 
445 std::size_t G4EMDataSet::FindLowerBound(G4doub << 
446 {                                                 113 {
447   std::size_t lowerBound = 0;                  << 114   // Protection against call outside allowed range
448   std::size_t upperBound(energies->size() - 1) << 115   G4double e0 = (*energies)[0];
                                                   >> 116   if (energy < e0)
                                                   >> 117     {
                                                   >> 118       //  G4cout << z
                                                   >> 119       //     << " - WARNING - G4EMDataSet::FindBinLocation called with argument " 
                                                   >> 120       //     << energy 
                                                   >> 121       //    << " outside lower limit " 
                                                   >> 122       //   << e0
                                                   >> 123       //   << "; replaced with lower limit" 
                                                   >> 124       //  << G4endl;
                                                   >> 125       energy = e0;
                                                   >> 126     }
                                                   >> 127 
                                                   >> 128   size_t lowerBound = 0;
                                                   >> 129   size_t upperBound = numberOfBins - 1;
449                                                   130   
                                                   >> 131   // Binary search
450   while (lowerBound <= upperBound)                132   while (lowerBound <= upperBound) 
451     {                                             133     {
452       std::size_t midBin((lowerBound + upperBo << 134       size_t midBin = (lowerBound + upperBound)/2;
453                                                << 135       if ( energy < (*energies)[midBin] ) upperBound = midBin-1;
454       if (x < (*energies)[midBin]) upperBound  << 136       else lowerBound = midBin+1;
455       else lowerBound = midBin + 1;            << 137   }
456     }                                          << 
457                                                   138   
458   return upperBound;                              139   return upperBound;
459 }                                                 140 }
460                                                   141 
461 //....oooOO0OOooo........oooOO0OOooo........oo << 142 void G4EMDataSet::LoadData(const G4String& fileName)
462                                                << 
463 std::size_t G4EMDataSet::FindLowerBound(G4doub << 
464 {                                                 143 {
465   std::size_t lowerBound = 0;;                 << 144   // Build the complete string identifying the file with the data set
466   std::size_t upperBound(values->size() - 1);  << 
467                                                   145   
468   while (lowerBound <= upperBound)             << 146   char nameChar[100] = {""};
469     {                                          << 147   std::ostrstream ost(nameChar, 100, std::ios::out);
470       std::size_t midBin((lowerBound + upperBo << 148   
471                                                << 149   ost << fileName << z << ".dat";
472       if (x < (*values)[midBin]) upperBound =  << 150   
473       else lowerBound = midBin + 1;            << 151   G4String name(nameChar);
                                                   >> 152   
                                                   >> 153   char* path = getenv("G4LEDATA");
                                                   >> 154   if (!path)
                                                   >> 155     { 
                                                   >> 156       G4String excep("G4EMDataSet - G4LEDATA environment variable not set");
                                                   >> 157       G4Exception(excep);
474     }                                             158     }
475                                                   159   
476   return upperBound;                           << 160   G4String pathString(path);
477 }                                              << 161   G4String separator("/");
478                                                   162 
479 //....oooOO0OOooo........oooOO0OOooo........oo << 163   G4String dirFile = pathString + separator + name;
                                                   >> 164   std::ifstream file(dirFile);
                                                   >> 165   std::filebuf* lsdp = file.rdbuf();
                                                   >> 166   
                                                   >> 167   if (! (lsdp->is_open()) )
                                                   >> 168   {
                                                   >> 169     G4String s1("G4EMDataSet - data file: ");
                                                   >> 170           G4String s2(" not found");
                                                   >> 171     G4String excep = s1 + dirFile + s2;
                                                   >> 172     G4Exception(excep);
                                                   >> 173   }
                                                   >> 174   G4double a = 0;
                                                   >> 175   G4int k = 1;
480                                                   176 
481 G4String G4EMDataSet::FullFileName(const G4Str << 177   do
482 {                                              << 178     {
483   const char* path = G4FindDataDir("G4LEDATA") << 179       file >> a;
484   if (!path) {                                 << 180       G4int nColumns = 2;
485      G4Exception("G4EMDataSet::FullFileName",  << 181       // The file is organized into two columns:
486         "em0006",FatalException,"G4LEDATA envi << 182       // 1st column is the energy
487     return "";                                 << 183       // 2nd column is the corresponding value
488   }                                            << 184       // The file terminates with the pattern: -1   -1
489   std::ostringstream fullFileName;             << 185       //                                       -2   -2
490   fullFileName << path << '/' << name << z <<  << 186       if (a == -1 || a == -2)
491                                                << 187   {
492   return G4String(fullFileName.str().c_str()); << 188   }
                                                   >> 189       else
                                                   >> 190   {
                                                   >> 191     if (k%nColumns != 0)
                                                   >> 192       { 
                                                   >> 193         G4double e = a * unit1;
                                                   >> 194         energies->push_back(e);
                                                   >> 195         k++;
                                                   >> 196       }
                                                   >> 197     else if (k%nColumns == 0)
                                                   >> 198       {
                                                   >> 199         G4double value = a * unit2;
                                                   >> 200         data->push_back(value);
                                                   >> 201         k = 1;
                                                   >> 202       }
                                                   >> 203   }
                                                   >> 204       
                                                   >> 205     } while (a != -2); // end of file
                                                   >> 206   
                                                   >> 207   file.close();
493 }                                                 208 }
494                                                   209 
495 //....oooOO0OOooo........oooOO0OOooo........oo << 210 void G4EMDataSet::PrintData() const
496                                                << 
497 void G4EMDataSet::BuildPdf()                   << 
498 {                                                 211 {
499   pdf = new G4DataVector;                      << 212   size_t size = numberOfBins;
500   G4Integrator <G4EMDataSet, G4double(G4EMData << 213   for (size_t i=0; i<size; i++)
501                                                << 
502   std::size_t nData = data->size();            << 
503   pdf->push_back(0.);                          << 
504                                                << 
505   // Integrate the data distribution           << 
506   std::size_t i;                               << 
507   G4double totalSum = 0.;                      << 
508   for (i=1; i<nData; ++i)                      << 
509     {                                             214     {
510       G4double xLow = (*energies)[i-1];        << 215       G4double e = (*energies)[i]  / unit1;
511       G4double xHigh = (*energies)[i];         << 216       G4double sigma = (*data)[i] / unit2 ;
512       G4double sum = integrator.Legendre96(thi << 217       G4cout << "Point: "
513       totalSum = totalSum + sum;               << 218        << e
514       pdf->push_back(totalSum);                << 219        << " - Data value : "
515     }                                          << 220        << sigma 
516                                                << 221        << G4endl; 
517   // Normalize to the last bin                 << 
518   G4double tot = 0.;                           << 
519   if (totalSum > 0.) tot = 1. / totalSum;      << 
520   for (i=1;  i<nData; ++i)                     << 
521     {                                          << 
522       (*pdf)[i] = (*pdf)[i] * tot;             << 
523     }                                             222     }
524 }                                                 223 }
525                                                   224 
526 //....oooOO0OOooo........oooOO0OOooo........oo << 
527                                                   225 
528 G4double G4EMDataSet::RandomSelect(G4int /* co << 226 const G4VEMDataSet* G4EMDataSet::GetComponent(G4int) const 
529 {                                              << 227 { return 0; }
530   G4double value = 0.;                         << 
531   // Random select a X value according to the  << 
532   // derived from the data                     << 
533                                                << 
534   if (!pdf) {                                  << 
535     G4Exception("G4EMDataSet::RandomSelect",   << 
536         "em1012",FatalException,"PDF has not b << 
537     return value;                              << 
538   }                                            << 
539                                                << 
540   G4double x = G4UniformRand();                << 
541                                                << 
542   // Locate the random value in the X vector b << 
543   G4int bin = (G4int)FindLowerBound(x,pdf);    << 
544                                                << 
545   // Interpolate the PDF to calculate the X va << 
546   // linear interpolation in the first bin (to << 
547   // interpolation with associated data set al << 
548   G4LinInterpolation linearAlgo;               << 
549   if (bin == 0) value = linearAlgo.Calculate(x << 
550   else value = algorithm->Calculate(x, bin, *p << 
551                                                << 
552   return value;                                << 
553 }                                              << 
554                                                << 
555 //....oooOO0OOooo........oooOO0OOooo........oo << 
556                                                << 
557 G4double G4EMDataSet::IntegrationFunction(G4do << 
558 {                                              << 
559   // This function is needed by G4Integrator t << 
560   G4double y = 0;                              << 
561                                                   228 
562   // Locate the random value in the X vector b << 229 void G4EMDataSet::AddComponent(G4VEMDataSet*) 
563   G4int bin = (G4int)FindLowerBound(x);        << 230 { }
564                                                << 
565   // Interpolate to calculate the X value:     << 
566   // linear interpolation in the first bin (to << 
567   // interpolation with associated algorithm i << 
568                                                << 
569   G4LinInterpolation linearAlgo;               << 
570                                                << 
571   if (bin == 0) y = linearAlgo.Calculate(x, bi << 
572   else y = algorithm->Calculate(x, bin, *energ << 
573                                                << 
574   return y;                                    << 
575 }                                              << 
576                                                   231 
                                                   >> 232 size_t G4EMDataSet::NumberOfComponents() const 
                                                   >> 233 { return 0; }
577                                                   234