physics/src/G4RDShellEMDataSet.cc

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 //  1 Aug 2001   MGP        Created
0033 // 09.10.01   V.Ivanchenko Add case z=0
0034 //  9 Mar 2008   MGP        Cleaned up unreadable code modified by former developer
0035 //                          (Further clean-up needed)
0036 //
0037 // -------------------------------------------------------------------
0038
0039 #include "G4RDShellEMDataSet.hh"
0040 #include "G4RDEMDataSet.hh"
0041 #include "G4RDVDataSetAlgorithm.hh"
0042 #include <fstream>
0043 #include <sstream>
0044
0045
0046 G4RDShellEMDataSet::G4RDShellEMDataSet(G4int zeta, G4RDVDataSetAlgorithm* algo,
0047                    G4double eUnit,
0048                    G4double dataUnit)
0049   :
0050   z(zeta),
0051   algorithm(algo),
0052   unitEnergies(eUnit),
0053   unitData(dataUnit)
0054 {
0055   if (algorithm == 0)
0056     G4Exception("G4RDShellEMDataSet::G4RDShellEMDataSet()", "InvalidSetup",
0057                 FatalException, "Interpolation == 0!");
0058 }
0059
0060
0061 G4RDShellEMDataSet::~G4RDShellEMDataSet()
0062 {
0063   CleanUpComponents();
0064   if (algorithm) delete algorithm;
0065 }
0066
0067
0068 G4double G4RDShellEMDataSet::FindValue(G4double energy, G4int /* componentId */) const
0069 {
0070   // Returns the sum over the shells corresponding to e
0071   G4double value = 0.;
0072
0073   std::vector<G4RDVEMDataSet *>::const_iterator i(components.begin());
0074   std::vector<G4RDVEMDataSet *>::const_iterator end(components.end());
0075
0076   while (i != end)
0077     {
0078       value += (*i)->FindValue(energy);
0079       i++;
0080     }
0081
0082   return value;
0083 }
0084
0085
0086 void G4RDShellEMDataSet::PrintData(void) const
0087 {
0088   const size_t n = NumberOfComponents();
0089
0090   G4cout << "The data set has " << n << " components" << G4endl;
0091   G4cout << G4endl;
0092
0093   size_t i = 0;
0094
0095   while (i < n)
0096     {
0097       G4cout << "--- Component " << i << " ---" << G4endl;
0098       GetComponent(i)->PrintData();
0099       i++;
0100     }
0101 }
0102
0103
0104 void G4RDShellEMDataSet::SetEnergiesData(G4DataVector* energies,
0105                        G4DataVector* data,
0106                        G4int componentId)
0107 {
0108   G4RDVEMDataSet* component = components[componentId];
0109
0110   if (component)
0111     {
0112       component->SetEnergiesData(energies, data, 0);
0113       return;
0114     }
0115
0116   std::ostringstream message;
0117   message << "Component " << componentId << " not found";
0118
0119   G4Exception("G4RDShellEMDataSet::SetEnergiesData()", "DataNotFound",
0120               FatalException, message.str().c_str());
0121 }
0122
0123
0124 G4bool G4RDShellEMDataSet::LoadData(const G4String& file)
0125 {
0126   CleanUpComponents();
0127
0128   G4String fullFileName = FullFileName(file);
0129   std::ifstream in(fullFileName);
0130
0131   if (!in.is_open())
0132     {
0133       G4String message("Data file \"");
0134       message += fullFileName;
0135       message += "\" not found";
0136       G4Exception("G4RDShellEMDataSet::LoadData()", "DataNotFound",
0137                   FatalException, message);
0138     }
0139
0140   G4DataVector* energies = 0;
0141   G4DataVector* data = 0;
0142
0143   G4double a = 0.;
0144   G4int shellIndex = 0;
0145   bool energyColumn = true;
0146
0147   do
0148     {
0149       in >> a;
0150
0151       if (a == -1)
0152     {
0153       if (energyColumn && energies!=0)
0154         {
0155           AddComponent(new G4RDEMDataSet(shellIndex, energies, data, algorithm->Clone(), unitEnergies, unitData));
0156           energies = 0;
0157           data = 0;
0158         }
0159
0160       energyColumn = (!energyColumn);
0161     }
0162       else if (a != -2)
0163     {
0164       if (energies == 0)
0165         {
0166           energies = new G4DataVector;
0167           data = new G4DataVector;
0168         }
0169
0170       if (energyColumn)
0171         energies->push_back(a * unitEnergies);
0172       else
0173         data->push_back(a * unitData);
0174
0175       energyColumn = (!energyColumn);
0176     }
0177     }
0178   while (a != -2);
0179
0180   return true;
0181 }
0182
0183
0184 G4bool G4RDShellEMDataSet::SaveData(const G4String& file) const
0185 {
0186   G4String fullFileName = FullFileName(file);
0187   std::ofstream out(fullFileName);
0188
0189   if (!out.is_open())
0190     {
0191       G4String message("Cannot open \"");
0192       message += fullFileName;
0193       message += "\"";
0194       G4Exception("G4RDEMDataSet::SaveData()", "CannotOpenFile",
0195                   FatalException, message);
0196     }
0197
0198   const size_t n = NumberOfComponents();
0199   size_t k = 0;
0200
0201   while (k < n)
0202     {
0203       const G4RDVEMDataSet* component = GetComponent(k);
0204
0205       if (component)
0206     {
0207       const G4DataVector& energies = component->GetEnergies(0);
0208       const G4DataVector& data = component->GetData(0);
0209
0210       G4DataVector::const_iterator i = energies.begin();
0211       G4DataVector::const_iterator endI = energies.end();
0212       G4DataVector::const_iterator j = data.begin();
0213
0214       while (i != endI)
0215         {
0216           out.precision(10);
0217           out.width(15);
0218           out.setf(std::ofstream::left);
0219           out << ((*i)/unitEnergies) << ' ';
0220
0221           out.precision(10);
0222           out.width(15);
0223           out.setf(std::ofstream::left);
0224           out << ((*j)/unitData) << std::endl;
0225           i++;
0226           j++;
0227         }
0228     }
0229
0230       out.precision(10);
0231       out.width(15);
0232       out.setf(std::ofstream::left);
0233       out << -1.f << ' ';
0234
0235       out.precision(10);
0236       out.width(15);
0237       out.setf(std::ofstream::left);
0238       out << -1.f << std::endl;
0239
0240       k++;
0241     }
0242
0243   out.precision(10);
0244   out.width(15);
0245   out.setf(std::ofstream::left);
0246   out << -2.f << ' ';
0247
0248   out.precision(10);
0249   out.width(15);
0250   out.setf(std::ofstream::left);
0251   out << -2.f << std::endl;
0252
0253   return true;
0254 }
0255
0256
0257 void G4RDShellEMDataSet::CleanUpComponents(void)
0258 {
0259   while (!components.empty())
0260     {
0261       if (components.back()) delete components.back();
0262       components.pop_back();
0263     }
0264 }
0265
0266
0267 G4String G4RDShellEMDataSet::FullFileName(const G4String& fileName) const
0268 {
0269   const char* path = G4FindDataDir("G4LEDATA");
0270   if (!path)
0271     G4Exception("G4RDShellEMDataSet::FullFileName()", "InvalidSetup",
0272                 FatalException, "G4LEDATA environment variable not set!");
0273
0274   std::ostringstream fullFileName;
0275
0276   fullFileName << path << '/' << fileName << z << ".dat";
0277
0278   return G4String(fullFileName.str().c_str());
0279 }