Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-01-31 09:22:07 

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 //
0027 //
0028 // Author: Maria Grazia Pia (Maria.Grazia.Pia@cern.ch)
0029 //
0030 // History:
0031 // -----------
0032 // 31 Jul 2001   MGP        Created
0033 //
0034 // -------------------------------------------------------------------
0035 
0036 #include "G4RDEMDataSet.hh"
0037 #include "G4RDVDataSetAlgorithm.hh"
0038 #include <fstream>
0039 #include <sstream>
0040 #include "G4Integrator.hh"
0041 #include "Randomize.hh"
0042 #include "G4RDLinInterpolation.hh"
0043 
0044 
0045 G4RDEMDataSet::G4RDEMDataSet(G4int Z, 
0046              G4RDVDataSetAlgorithm* algo, 
0047              G4double xUnit, 
0048              G4double yUnit,
0049              G4bool random): 
0050   z(Z),
0051   energies(0),
0052   data(0),
0053   algorithm(algo),
0054   unitEnergies(xUnit),
0055   unitData(yUnit),
0056   pdf(0),
0057   randomSet(random)
0058 {
0059   if (algorithm == 0)
0060     G4Exception("G4RDEMDataSet::G4RDEMDataSet()",
0061                 "InvalidSetup", FatalException, "Interpolation == 0!");
0062   if (randomSet) BuildPdf();
0063 }
0064 
0065 G4RDEMDataSet::G4RDEMDataSet(G4int argZ, 
0066              G4DataVector* dataX, 
0067              G4DataVector* dataY, 
0068              G4RDVDataSetAlgorithm* algo, 
0069              G4double xUnit, 
0070              G4double yUnit,
0071              G4bool random):
0072   z(argZ),
0073   energies(dataX),
0074   data(dataY),
0075   algorithm(algo),
0076   unitEnergies(xUnit),
0077   unitData(yUnit),
0078   pdf(0),
0079   randomSet(random)
0080 {
0081   if (algorithm == 0)
0082     G4Exception("G4RDEMDataSet::G4RDEMDataSet()",
0083                 "InvalidSetup", FatalException, "Interpolation == 0!");
0084 
0085   if ((energies == 0) ^ (data == 0))
0086     G4Exception("G4RDEMDataSet::G4RDEMDataSet()", "InvalidSetup",
0087            FatalException, "Different size for energies and data (zero case)!");
0088 
0089   if (energies == 0) return;
0090   
0091   if (energies->size() != data->size()) 
0092     G4Exception("G4RDEMDataSet::G4RDEMDataSet()", "InvalidSetup",
0093            FatalException, "Different size for energies and data!");
0094 
0095   if (randomSet) BuildPdf();
0096 }
0097 
0098 G4RDEMDataSet::~G4RDEMDataSet()
0099 { 
0100   delete algorithm;
0101   if (energies) delete energies;
0102   if (data) delete data;
0103   if (pdf) delete pdf;
0104 }
0105 
0106 G4double G4RDEMDataSet::FindValue(G4double energy, G4int /* componentId */) const
0107 {
0108   if (!energies)
0109     G4Exception("G4RDEMDataSet::FindValue()", "InvalidSetup",
0110                 FatalException, "Energies == 0!");
0111   if (energies->empty()) return 0;
0112   if (energy <= (*energies)[0]) return (*data)[0];
0113 
0114   size_t i = energies->size()-1;
0115   if (energy >= (*energies)[i]) return (*data)[i];
0116 
0117   return algorithm->Calculate(energy, FindLowerBound(energy), *energies, *data);
0118 }
0119 
0120 
0121 void G4RDEMDataSet::PrintData(void) const
0122 {
0123   if (!energies)
0124     {
0125       G4cout << "Data not available." << G4endl;
0126     }
0127   else
0128     {
0129       size_t size = energies->size();
0130       for (size_t i(0); i<size; i++)
0131     {
0132       G4cout << "Point: " << ((*energies)[i]/unitEnergies)
0133          << " - Data value: " << ((*data)[i]/unitData);
0134       if (pdf != 0) G4cout << " - PDF : " << (*pdf)[i];
0135       G4cout << G4endl; 
0136     }
0137     }
0138 }
0139 
0140 
0141 void G4RDEMDataSet::SetEnergiesData(G4DataVector* dataX, 
0142                   G4DataVector* dataY, 
0143                   G4int /* componentId */)
0144 {
0145   if (energies) delete energies;
0146   energies = dataX;
0147 
0148   if (data) delete data;
0149   data = dataY;
0150  
0151   if ((energies == 0) ^ (data==0)) 
0152     G4Exception("G4RDEMDataSet::SetEnergiesData()", "InvalidSetup",
0153        FatalException, "Different size for energies and data (zero case)!");
0154 
0155   if (energies == 0) return;
0156   
0157   if (energies->size() != data->size()) 
0158     G4Exception("G4RDEMDataSet::SetEnergiesData()", "InvalidSetup",
0159        FatalException, "Different size for energies and data!");
0160 }
0161 
0162 G4bool G4RDEMDataSet::LoadData(const G4String& fileName)
0163 {
0164   // The file is organized into two columns:
0165   // 1st column is the energy
0166   // 2nd column is the corresponding value
0167   // The file terminates with the pattern: -1   -1
0168   //                                       -2   -2
0169  
0170   G4String fullFileName(FullFileName(fileName));
0171   std::ifstream in(fullFileName);
0172 
0173   if (!in.is_open())
0174     {
0175       G4String message("Data file \"");
0176       message += fullFileName;
0177       message += "\" not found";
0178       G4Exception("G4RDEMDataSet::LoadData()", "DataNotFound",
0179                   FatalException, message);
0180     }
0181 
0182   G4DataVector* argEnergies=new G4DataVector;
0183   G4DataVector* argData=new G4DataVector;
0184 
0185   G4double a;
0186   bool energyColumn(true);
0187 
0188   do
0189     {
0190       in >> a;
0191   
0192       if (a!=-1 && a!=-2)
0193     {
0194       if (energyColumn)
0195         argEnergies->push_back(a*unitEnergies);
0196       else
0197         argData->push_back(a*unitData);
0198       energyColumn=(!energyColumn);
0199     }
0200     }
0201   while (a != -2);
0202  
0203   SetEnergiesData(argEnergies, argData, 0);
0204   if (randomSet) BuildPdf();
0205  
0206   return true;
0207 }
0208 
0209 G4bool G4RDEMDataSet::SaveData(const G4String& name) const
0210 {
0211   // The file is organized into two columns:
0212   // 1st column is the energy
0213   // 2nd column is the corresponding value
0214   // The file terminates with the pattern: -1   -1
0215   //                                       -2   -2
0216  
0217   G4String fullFileName(FullFileName(name));
0218   std::ofstream out(fullFileName);
0219 
0220   if (!out.is_open())
0221     {
0222       G4String message("Cannot open \"");
0223       message+=fullFileName;
0224       message+="\"";
0225       G4Exception("G4RDEMDataSet::SaveData()", "CannotOpenFile",
0226                   FatalException, message);
0227     }
0228  
0229   out.precision(10);
0230   out.width(15);
0231   out.setf(std::ofstream::left);
0232   
0233   if (energies!=0 && data!=0)
0234     {
0235       G4DataVector::const_iterator i(energies->begin());
0236       G4DataVector::const_iterator endI(energies->end());
0237       G4DataVector::const_iterator j(data->begin());
0238   
0239       while (i!=endI)
0240     {
0241       out.precision(10);
0242       out.width(15);
0243       out.setf(std::ofstream::left);
0244       out << ((*i)/unitEnergies) << ' ';
0245 
0246       out.precision(10);
0247       out.width(15);
0248       out.setf(std::ofstream::left);
0249       out << ((*j)/unitData) << std::endl;
0250 
0251       i++;
0252       j++;
0253     }
0254     }
0255  
0256   out.precision(10);
0257   out.width(15);
0258   out.setf(std::ofstream::left);
0259   out << -1.f << ' ';
0260 
0261   out.precision(10);
0262   out.width(15);
0263   out.setf(std::ofstream::left);
0264   out << -1.f << std::endl;
0265 
0266   out.precision(10);
0267   out.width(15);
0268   out.setf(std::ofstream::left);
0269   out << -2.f << ' ';
0270 
0271   out.precision(10);
0272   out.width(15);
0273   out.setf(std::ofstream::left);
0274   out << -2.f << std::endl;
0275 
0276   return true;
0277 }
0278 
0279 size_t G4RDEMDataSet::FindLowerBound(G4double x) const
0280 {
0281   size_t lowerBound = 0;
0282   size_t upperBound(energies->size() - 1);
0283   
0284   while (lowerBound <= upperBound) 
0285     {
0286       size_t midBin((lowerBound + upperBound) / 2);
0287 
0288       if (x < (*energies)[midBin]) upperBound = midBin - 1;
0289       else lowerBound = midBin + 1;
0290     }
0291   
0292   return upperBound;
0293 }
0294 
0295 
0296 size_t G4RDEMDataSet::FindLowerBound(G4double x, G4DataVector* values) const
0297 {
0298   size_t lowerBound = 0;;
0299   size_t upperBound(values->size() - 1);
0300   
0301   while (lowerBound <= upperBound) 
0302     {
0303       size_t midBin((lowerBound + upperBound) / 2);
0304 
0305       if (x < (*values)[midBin]) upperBound = midBin - 1;
0306       else lowerBound = midBin + 1;
0307     }
0308   
0309   return upperBound;
0310 }
0311 
0312 
0313 G4String G4RDEMDataSet::FullFileName(const G4String& name) const
0314 {
0315   const char* path = G4FindDataDir("G4LEDATA");
0316   if (!path)
0317     G4Exception("G4RDEMDataSet::FullFileName()", "InvalidSetup",
0318                 FatalException, "G4LEDATA environment variable not set!");
0319   
0320   std::ostringstream fullFileName;
0321   fullFileName << path << '/' << name << z << ".dat";
0322                       
0323   return G4String(fullFileName.str().c_str());
0324 }
0325 
0326 
0327 void G4RDEMDataSet::BuildPdf() 
0328 {
0329   pdf = new G4DataVector;
0330   G4Integrator <G4RDEMDataSet, G4double(G4RDEMDataSet::*)(G4double)> integrator;
0331 
0332   G4int nData = data->size();
0333   pdf->push_back(0.);
0334 
0335   // Integrate the data distribution 
0336   G4int i;
0337   G4double totalSum = 0.;
0338   for (i=1; i<nData; i++)
0339     {
0340       G4double xLow = (*energies)[i-1];
0341       G4double xHigh = (*energies)[i];
0342       G4double sum = integrator.Legendre96(this, &G4RDEMDataSet::IntegrationFunction, xLow, xHigh);
0343       totalSum = totalSum + sum;
0344       pdf->push_back(totalSum);
0345     }
0346 
0347   // Normalize to the last bin
0348   G4double tot = 0.;
0349   if (totalSum > 0.) tot = 1. / totalSum;
0350   for (i=1;  i<nData; i++)
0351     {
0352       (*pdf)[i] = (*pdf)[i] * tot;
0353     }
0354 }
0355 
0356 
0357 G4double G4RDEMDataSet::RandomSelect(G4int /* componentId */) const
0358 {
0359   // Random select a X value according to the cumulative probability distribution
0360   // derived from the data
0361 
0362   if (!pdf) G4Exception("G4RDEMDataSet::RandomSelect()", "InvalidSetup",
0363            FatalException, "PDF has not been created for this data set");
0364 
0365   G4double value = 0.;
0366   G4double x = G4UniformRand();
0367 
0368   // Locate the random value in the X vector based on the PDF
0369   size_t bin = FindLowerBound(x,pdf);
0370 
0371   // Interpolate the PDF to calculate the X value: 
0372   // linear interpolation in the first bin (to avoid problem with 0),
0373   // interpolation with associated data set algorithm in other bins
0374 
0375   G4RDLinInterpolation linearAlgo;
0376   if (bin == 0) value = linearAlgo.Calculate(x, bin, *pdf, *energies);
0377   else value = algorithm->Calculate(x, bin, *pdf, *energies);
0378 
0379   //  G4cout << x << " random bin "<< bin << " - " << value << G4endl;
0380   return value;
0381 }
0382 
0383 G4double G4RDEMDataSet::IntegrationFunction(G4double x)
0384 {
0385   // This function is needed by G4Integrator to calculate the integral of the data distribution
0386 
0387   G4double y = 0;
0388 
0389  // Locate the random value in the X vector based on the PDF
0390   size_t bin = FindLowerBound(x);
0391 
0392   // Interpolate to calculate the X value: 
0393   // linear interpolation in the first bin (to avoid problem with 0),
0394   // interpolation with associated algorithm in other bins
0395 
0396   G4RDLinInterpolation linearAlgo;
0397   
0398   if (bin == 0) y = linearAlgo.Calculate(x, bin, *energies, *data);
0399   else y = algorithm->Calculate(x, bin, *energies, *data);
0400  
0401   return y;
0402 }
0403
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!