Geant4 Cross Reference

Cross-Referencing   Geant4
Geant4/examples/advanced/eRosita/physics/src/G4RDEMDataSet.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 /examples/advanced/eRosita/physics/src/G4RDEMDataSet.cc (Version 11.3.0) and /examples/advanced/eRosita/physics/src/G4RDEMDataSet.cc (Version 1.0)


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