physics/src/G4RDDopplerProfile.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 // 31 Jul 2001   MGP        Created
0033 //
0034 // -------------------------------------------------------------------
0035
0036 #include "G4RDDopplerProfile.hh"
0037 #include "G4DataVector.hh"
0038 #include "G4SystemOfUnits.hh"
0039 #include "G4RDVEMDataSet.hh"
0040 #include "G4RDEMDataSet.hh"
0041 #include "G4RDCompositeEMDataSet.hh"
0042 #include "G4RDVDataSetAlgorithm.hh"
0043 #include "G4RDLogLogInterpolation.hh"
0044
0045 #include <fstream>
0046 #include <sstream>
0047 #include "Randomize.hh"
0048
0049 // The following deprecated header is included because <functional> seems not to be found on MGP's laptop
0050 //#include "function.h"
0051
0052 // Constructor
0053
0054 G4RDDopplerProfile::G4RDDopplerProfile(G4int minZ, G4int maxZ)
0055   : zMin(minZ), zMax(maxZ)
0056 {
0057   nBiggs = 31;
0058
0059   LoadBiggsP("/doppler/p-biggs");
0060
0061   for (G4int Z=zMin; Z<zMax+1; Z++)
0062     {
0063       LoadProfile("/doppler/profile",Z);
0064     }
0065 }
0066
0067 // Destructor
0068 G4RDDopplerProfile::~G4RDDopplerProfile()
0069 {
0070   std::map<G4int,G4RDVEMDataSet*,std::less<G4int> >::iterator pos;
0071   for (pos = profileMap.begin(); pos != profileMap.end(); ++pos)
0072     {
0073       G4RDVEMDataSet* dataSet = (*pos).second;
0074       delete dataSet;
0075       dataSet = 0;
0076     }
0077 }
0078
0079
0080 size_t G4RDDopplerProfile::NumberOfProfiles(G4int Z) const
0081 {
0082   G4int n = 0;
0083   if (Z>= zMin && Z <= zMax) n = nShells[Z-1];
0084   return n;
0085 }
0086
0087
0088 const G4RDVEMDataSet* G4RDDopplerProfile::Profiles(G4int Z) const
0089 {
0090   std::map<G4int,G4RDVEMDataSet*,std::less<G4int> >::const_iterator pos;
0091   if (Z < zMin || Z > zMax)
0092     G4Exception("G4RDDopplerProfile::Profiles()", "OutOfRange",
0093                 FatalException, "Z outside boundaries!");
0094   pos = profileMap.find(Z);
0095   G4RDVEMDataSet* dataSet = (*pos).second;
0096   return dataSet;
0097 }
0098
0099
0100 const G4RDVEMDataSet* G4RDDopplerProfile::Profile(G4int Z, G4int shellIndex) const
0101 {
0102   const G4RDVEMDataSet* profis = Profiles(Z);
0103   const G4RDVEMDataSet* profi = profis->GetComponent(shellIndex);
0104   return profi;
0105 }
0106
0107
0108 void G4RDDopplerProfile::PrintData() const
0109 {
0110   for (G4int Z=zMin; Z<zMax; Z++)
0111     {
0112       const G4RDVEMDataSet* profis = Profiles(Z);
0113       profis->PrintData();
0114     }
0115 }
0116
0117
0118 void G4RDDopplerProfile::LoadBiggsP(const G4String& fileName)
0119 {
0120   std::ostringstream ost;
0121   ost << fileName << ".dat";
0122   G4String name(ost.str());
0123
0124   const char* path = G4FindDataDir("G4LEDATA");
0125   if (!path)
0126     {
0127       G4String excep("G4LEDATA environment variable not set!");
0128       G4Exception("G4RDDopplerProfile::LoadBiggsP()",
0129                   "InvalidSetup", FatalException, excep);
0130     }
0131
0132   G4String pathString(path);
0133   G4String dirFile = pathString + name;
0134   std::ifstream file(dirFile);
0135   std::filebuf* lsdp = file.rdbuf();
0136
0137   if (! (lsdp->is_open()) )
0138     {
0139       G4String s1("Data file: ");
0140       G4String s2(" not found");
0141       G4String excep = s1 + dirFile + s2;
0142       G4Exception("G4RDDopplerProfile::LoadBiggsP()",
0143                   "DataNotFound", FatalException, excep);
0144     }
0145
0146   G4double p;
0147   while(!file.eof())
0148     {
0149       file >> p;
0150       biggsP.push_back(p);
0151     }
0152
0153   // Make sure that the number of data loaded corresponds to the number in Biggs' paper
0154   if (biggsP.size() != nBiggs)
0155     G4Exception("G4RDDopplerProfile::LoadBiggsP()", "InvalidCondition",
0156                 FatalException, "Number of momenta read in is not 31!");
0157 }
0158
0159
0160 void G4RDDopplerProfile::LoadProfile(const G4String& fileName,G4int Z)
0161 {
0162   std::ostringstream ost;
0163   ost << fileName << "-" << Z << ".dat";
0164   G4String name(ost.str());
0165
0166   const char* path = G4FindDataDir("G4LEDATA");
0167   if (!path)
0168     {
0169       G4String excep("G4LEDATA environment variable not set!");
0170       G4Exception("G4RDDopplerProfile::LoadProfile()",
0171                   "InvalidSetup", FatalException, excep);
0172     }
0173
0174   G4String pathString(path);
0175   G4String dirFile = pathString + name;
0176   std::ifstream file(dirFile);
0177   std::filebuf* lsdp = file.rdbuf();
0178
0179   if (! (lsdp->is_open()) )
0180     {
0181       G4String s1("Data file: ");
0182       G4String s2(" not found");
0183       G4String excep = s1 + dirFile + s2;
0184       G4Exception("G4RDDopplerProfile::LoadProfile()",
0185                   "DataNotFound", FatalException, excep);
0186     }
0187
0188   G4double p;
0189   G4int nShell = 0;
0190
0191   // Create CompositeDataSet for the current Z
0192   G4RDVDataSetAlgorithm* interpolation = new G4RDLogLogInterpolation;
0193   G4RDVEMDataSet* dataSetForZ = new G4RDCompositeEMDataSet(interpolation,1.,1.,1,1);
0194
0195   while (!file.eof())
0196     {
0197       nShell++;
0198       G4DataVector* profi = new G4DataVector;
0199       G4DataVector* biggs = new G4DataVector;
0200
0201       // Read in profile data for the current shell
0202       for (size_t i=0; i<nBiggs; i++)
0203     {
0204       file >> p;
0205       profi->push_back(p);
0206           biggs->push_back(biggsP[i]);
0207       //      if (i == 16) G4cout << "profile = " << p << G4endl;
0208     }
0209
0210       // Create G4RDEMDataSet for the current shell
0211       G4RDVDataSetAlgorithm* algo = interpolation->Clone();
0212       G4RDVEMDataSet* dataSet = new G4RDEMDataSet(Z, biggs, profi, algo, 1., 1., true);
0213
0214       // Add current shell profile component to G4RDCompositeEMDataSet for the current Z
0215       dataSetForZ->AddComponent(dataSet);
0216     }
0217
0218   // Fill in number of shells for the current Z
0219   nShells.push_back(nShell);
0220
0221   profileMap[Z] = dataSetForZ;
0222 }
0223
0224
0225  G4double G4RDDopplerProfile::RandomSelectMomentum(G4int Z, G4int shellIndex) const
0226 {
0227   G4double value = 0.;
0228   const G4RDVEMDataSet* profis = Profiles(Z);
0229   value = profis->RandomSelect(shellIndex);
0230   return value;
0231 }