00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031
00032
00033
00034
00035
00036
00037 #include "G4DopplerProfile.hh"
00038 #include "G4DataVector.hh"
00039 #include "G4SystemOfUnits.hh"
00040 #include "G4VEMDataSet.hh"
00041 #include "G4EMDataSet.hh"
00042 #include "G4CompositeEMDataSet.hh"
00043 #include "G4VDataSetAlgorithm.hh"
00044 #include "G4LogLogInterpolation.hh"
00045
00046 #include <fstream>
00047 #include <sstream>
00048 #include "Randomize.hh"
00049
00050
00051
00052
00053
00054
00055 G4DopplerProfile::G4DopplerProfile(G4int minZ, G4int maxZ)
00056 : zMin(minZ), zMax(maxZ)
00057 {
00058 nBiggs = 31;
00059
00060 LoadBiggsP("/doppler/p-biggs");
00061
00062 for (G4int Z=zMin; Z<zMax+1; Z++)
00063 {
00064 LoadProfile("/doppler/profile",Z);
00065 }
00066 }
00067
00068
00069 G4DopplerProfile::~G4DopplerProfile()
00070 {
00071 std::map<G4int,G4VEMDataSet*,std::less<G4int> >::iterator pos;
00072 for (pos = profileMap.begin(); pos != profileMap.end(); ++pos)
00073 {
00074 G4VEMDataSet* dataSet = (*pos).second;
00075 delete dataSet;
00076 dataSet = 0;
00077 }
00078 }
00079
00080
00081 size_t G4DopplerProfile::NumberOfProfiles(G4int Z) const
00082 {
00083 G4int n = 0;
00084 if (Z>= zMin && Z <= zMax) n = nShells[Z-1];
00085 return n;
00086 }
00087
00088
00089 const G4VEMDataSet* G4DopplerProfile::Profiles(G4int Z) const
00090 {
00091 std::map<G4int,G4VEMDataSet*,std::less<G4int> >::const_iterator pos;
00092 if (Z < zMin || Z > zMax)
00093 G4Exception("G4DopplerProfile::Profiles",
00094 "em1005",FatalException,"Z outside boundaries");
00095 pos = profileMap.find(Z);
00096 G4VEMDataSet* dataSet = (*pos).second;
00097 return dataSet;
00098 }
00099
00100
00101 const G4VEMDataSet* G4DopplerProfile::Profile(G4int Z, G4int shellIndex) const
00102 {
00103 const G4VEMDataSet* profis = Profiles(Z);
00104 const G4VEMDataSet* profi = profis->GetComponent(shellIndex);
00105 return profi;
00106 }
00107
00108
00109 void G4DopplerProfile::PrintData() const
00110 {
00111 for (G4int Z=zMin; Z<zMax; Z++)
00112 {
00113 const G4VEMDataSet* profis = Profiles(Z);
00114 profis->PrintData();
00115 }
00116 }
00117
00118
00119 void G4DopplerProfile::LoadBiggsP(const G4String& fileName)
00120 {
00121 std::ostringstream ost;
00122 ost << fileName << ".dat";
00123 G4String name(ost.str());
00124
00125 char* path = getenv("G4LEDATA");
00126 if (!path)
00127 {
00128 G4Exception("G4DopplerProfile::LoadBiggsP",
00129 "em0006",FatalException,"G4LEDATA environment variable not set");
00130 return;
00131 }
00132
00133 G4String pathString(path);
00134 G4String dirFile = pathString + name;
00135 std::ifstream file(dirFile);
00136 std::filebuf* lsdp = file.rdbuf();
00137
00138 if (! (lsdp->is_open()) )
00139 {
00140 G4String s1("data file: ");
00141 G4String s2(" not found");
00142 G4String excep = s1 + dirFile + s2;
00143 G4Exception("G4DopplerProfile::LoadBiggsP",
00144 "em0003",FatalException,excep);
00145 }
00146
00147 G4double p;
00148 while(!file.eof())
00149 {
00150 file >> p;
00151 biggsP.push_back(p);
00152 }
00153
00154
00155 if (biggsP.size() != nBiggs)
00156 G4Exception("G4DopplerProfile::LoadBiggsP",
00157 "em1006",FatalException,"Number of momenta read in is not 31");
00158 }
00159
00160
00161 void G4DopplerProfile::LoadProfile(const G4String& fileName,G4int Z)
00162 {
00163 std::ostringstream ost;
00164 ost << fileName << "-" << Z << ".dat";
00165 G4String name(ost.str());
00166
00167 char* path = getenv("G4LEDATA");
00168 if (!path)
00169 {
00170 G4String excep("G4LEDATA environment variable not set");
00171 G4Exception("G4DopplerProfile::LoadProfile",
00172 "em0006",FatalException,excep);
00173 return;
00174 }
00175
00176 G4String pathString(path);
00177 G4String dirFile = pathString + name;
00178 std::ifstream file(dirFile);
00179 std::filebuf* lsdp = file.rdbuf();
00180
00181 if (! (lsdp->is_open()) )
00182 {
00183 G4String s1("data file: ");
00184 G4String s2(" not found");
00185 G4String excep = s1 + dirFile + s2;
00186 G4Exception("G4DopplerProfile::LoadProfile",
00187 "em0003",FatalException,excep);
00188 }
00189
00190 G4double p;
00191 G4int nShell = 0;
00192
00193
00194 G4VDataSetAlgorithm* interpolation = new G4LogLogInterpolation;
00195 G4VEMDataSet* dataSetForZ = new G4CompositeEMDataSet(interpolation,1.,1.,1,1);
00196
00197 while (!file.eof())
00198 {
00199 nShell++;
00200 G4DataVector* profi = new G4DataVector;
00201 G4DataVector* biggs = new G4DataVector;
00202
00203
00204 for (size_t i=0; i<nBiggs; i++)
00205 {
00206 file >> p;
00207 profi->push_back(p);
00208 biggs->push_back(biggsP[i]);
00209
00210 }
00211
00212
00213 G4VDataSetAlgorithm* algo = interpolation->Clone();
00214 G4VEMDataSet* dataSet = new G4EMDataSet(Z, biggs, profi, algo, 1., 1., true);
00215
00216
00217 dataSetForZ->AddComponent(dataSet);
00218 }
00219
00220
00221 nShells.push_back(nShell);
00222
00223 profileMap[Z] = dataSetForZ;
00224 }
00225
00226
00227 G4double G4DopplerProfile::RandomSelectMomentum(G4int Z, G4int shellIndex) const
00228 {
00229 G4double value = 0.;
00230 const G4VEMDataSet* profis = Profiles(Z);
00231 value = profis->RandomSelect(shellIndex);
00232 return value;
00233 }