| Geant4 Cross Reference |
1 // 1 //
2 // ******************************************* 2 // ********************************************************************
3 // * License and Disclaimer 3 // * License and Disclaimer *
4 // * 4 // * *
5 // * The Geant4 software is copyright of th 5 // * The Geant4 software is copyright of the Copyright Holders of *
6 // * the Geant4 Collaboration. It is provided 6 // * the Geant4 Collaboration. It is provided under the terms and *
7 // * conditions of the Geant4 Software License 7 // * conditions of the Geant4 Software License, included in the file *
8 // * LICENSE and available at http://cern.ch/ 8 // * LICENSE and available at http://cern.ch/geant4/license . These *
9 // * include a list of copyright holders. 9 // * include a list of copyright holders. *
10 // * 10 // * *
11 // * Neither the authors of this software syst 11 // * Neither the authors of this software system, nor their employing *
12 // * institutes,nor the agencies providing fin 12 // * institutes,nor the agencies providing financial support for this *
13 // * work make any representation or warran 13 // * work make any representation or warranty, express or implied, *
14 // * regarding this software system or assum 14 // * regarding this software system or assume any liability for its *
15 // * use. Please see the license in the file 15 // * use. Please see the license in the file LICENSE and URL above *
16 // * for the full disclaimer and the limitatio 16 // * for the full disclaimer and the limitation of liability. *
17 // * 17 // * *
18 // * This code implementation is the result 18 // * This code implementation is the result of the scientific and *
19 // * technical work of the GEANT4 collaboratio 19 // * technical work of the GEANT4 collaboration. *
20 // * By using, copying, modifying or distri 20 // * By using, copying, modifying or distributing the software (or *
21 // * any work based on the software) you ag 21 // * any work based on the software) you agree to acknowledge its *
22 // * use in resulting scientific publicati 22 // * use in resulting scientific publications, and indicate your *
23 // * acceptance of all terms of the Geant4 Sof 23 // * acceptance of all terms of the Geant4 Software license. *
24 // ******************************************* 24 // ********************************************************************
25 // 25 //
26 // 26 //
>> 27 // $Id$
>> 28 // GEANT4 tag $Name: $
27 // 29 //
28 // Author: Maria Grazia Pia (Maria.Grazia.Pia@ 30 // Author: Maria Grazia Pia (Maria.Grazia.Pia@cern.ch)
29 // 31 //
30 // History: 32 // History:
31 // ----------- 33 // -----------
32 // 31 Jul 2001 MGP Created 34 // 31 Jul 2001 MGP Created
33 // 35 //
34 // ------------------------------------------- 36 // -------------------------------------------------------------------
35 37
36 #include "G4RDDopplerProfile.hh" 38 #include "G4RDDopplerProfile.hh"
37 #include "G4DataVector.hh" 39 #include "G4DataVector.hh"
38 #include "G4SystemOfUnits.hh" 40 #include "G4SystemOfUnits.hh"
39 #include "G4RDVEMDataSet.hh" 41 #include "G4RDVEMDataSet.hh"
40 #include "G4RDEMDataSet.hh" 42 #include "G4RDEMDataSet.hh"
41 #include "G4RDCompositeEMDataSet.hh" 43 #include "G4RDCompositeEMDataSet.hh"
42 #include "G4RDVDataSetAlgorithm.hh" 44 #include "G4RDVDataSetAlgorithm.hh"
43 #include "G4RDLogLogInterpolation.hh" 45 #include "G4RDLogLogInterpolation.hh"
44 46
45 #include <fstream> 47 #include <fstream>
46 #include <sstream> 48 #include <sstream>
47 #include "Randomize.hh" 49 #include "Randomize.hh"
48 50
49 // The following deprecated header is included 51 // The following deprecated header is included because <functional> seems not to be found on MGP's laptop
50 //#include "function.h" 52 //#include "function.h"
51 53
52 // Constructor 54 // Constructor
53 55
54 G4RDDopplerProfile::G4RDDopplerProfile(G4int m 56 G4RDDopplerProfile::G4RDDopplerProfile(G4int minZ, G4int maxZ)
55 : zMin(minZ), zMax(maxZ) 57 : zMin(minZ), zMax(maxZ)
56 { 58 {
57 nBiggs = 31; 59 nBiggs = 31;
58 60
59 LoadBiggsP("/doppler/p-biggs"); 61 LoadBiggsP("/doppler/p-biggs");
60 62
61 for (G4int Z=zMin; Z<zMax+1; Z++) 63 for (G4int Z=zMin; Z<zMax+1; Z++)
62 { 64 {
63 LoadProfile("/doppler/profile",Z); 65 LoadProfile("/doppler/profile",Z);
64 } 66 }
65 } 67 }
66 68
67 // Destructor 69 // Destructor
68 G4RDDopplerProfile::~G4RDDopplerProfile() 70 G4RDDopplerProfile::~G4RDDopplerProfile()
69 { 71 {
70 std::map<G4int,G4RDVEMDataSet*,std::less<G4i 72 std::map<G4int,G4RDVEMDataSet*,std::less<G4int> >::iterator pos;
71 for (pos = profileMap.begin(); pos != profil 73 for (pos = profileMap.begin(); pos != profileMap.end(); ++pos)
72 { 74 {
73 G4RDVEMDataSet* dataSet = (*pos).second; 75 G4RDVEMDataSet* dataSet = (*pos).second;
74 delete dataSet; 76 delete dataSet;
75 dataSet = 0; 77 dataSet = 0;
76 } 78 }
77 } 79 }
78 80
79 81
80 size_t G4RDDopplerProfile::NumberOfProfiles(G4 82 size_t G4RDDopplerProfile::NumberOfProfiles(G4int Z) const
81 { 83 {
82 G4int n = 0; 84 G4int n = 0;
83 if (Z>= zMin && Z <= zMax) n = nShells[Z-1]; 85 if (Z>= zMin && Z <= zMax) n = nShells[Z-1];
84 return n; 86 return n;
85 } 87 }
86 88
87 89
88 const G4RDVEMDataSet* G4RDDopplerProfile::Prof 90 const G4RDVEMDataSet* G4RDDopplerProfile::Profiles(G4int Z) const
89 { 91 {
90 std::map<G4int,G4RDVEMDataSet*,std::less<G4i 92 std::map<G4int,G4RDVEMDataSet*,std::less<G4int> >::const_iterator pos;
91 if (Z < zMin || Z > zMax) 93 if (Z < zMin || Z > zMax)
92 G4Exception("G4RDDopplerProfile::Profiles( 94 G4Exception("G4RDDopplerProfile::Profiles()", "OutOfRange",
93 FatalException, "Z outside bou 95 FatalException, "Z outside boundaries!");
94 pos = profileMap.find(Z); 96 pos = profileMap.find(Z);
95 G4RDVEMDataSet* dataSet = (*pos).second; 97 G4RDVEMDataSet* dataSet = (*pos).second;
96 return dataSet; 98 return dataSet;
97 } 99 }
98 100
99 101
100 const G4RDVEMDataSet* G4RDDopplerProfile::Prof 102 const G4RDVEMDataSet* G4RDDopplerProfile::Profile(G4int Z, G4int shellIndex) const
101 { 103 {
102 const G4RDVEMDataSet* profis = Profiles(Z); 104 const G4RDVEMDataSet* profis = Profiles(Z);
103 const G4RDVEMDataSet* profi = profis->GetCom 105 const G4RDVEMDataSet* profi = profis->GetComponent(shellIndex);
104 return profi; 106 return profi;
105 } 107 }
106 108
107 109
108 void G4RDDopplerProfile::PrintData() const 110 void G4RDDopplerProfile::PrintData() const
109 { 111 {
110 for (G4int Z=zMin; Z<zMax; Z++) 112 for (G4int Z=zMin; Z<zMax; Z++)
111 { 113 {
112 const G4RDVEMDataSet* profis = Profiles( 114 const G4RDVEMDataSet* profis = Profiles(Z);
113 profis->PrintData(); 115 profis->PrintData();
114 } 116 }
115 } 117 }
116 118
117 119
118 void G4RDDopplerProfile::LoadBiggsP(const G4St 120 void G4RDDopplerProfile::LoadBiggsP(const G4String& fileName)
119 { 121 {
120 std::ostringstream ost; 122 std::ostringstream ost;
121 ost << fileName << ".dat"; 123 ost << fileName << ".dat";
122 G4String name(ost.str()); 124 G4String name(ost.str());
123 125
124 const char* path = G4FindDataDir("G4LEDATA") << 126 char* path = getenv("G4LEDATA");
125 if (!path) 127 if (!path)
126 { 128 {
127 G4String excep("G4LEDATA environment var 129 G4String excep("G4LEDATA environment variable not set!");
128 G4Exception("G4RDDopplerProfile::LoadBig 130 G4Exception("G4RDDopplerProfile::LoadBiggsP()",
129 "InvalidSetup", FatalExcepti 131 "InvalidSetup", FatalException, excep);
130 } 132 }
131 133
132 G4String pathString(path); 134 G4String pathString(path);
133 G4String dirFile = pathString + name; 135 G4String dirFile = pathString + name;
134 std::ifstream file(dirFile); 136 std::ifstream file(dirFile);
135 std::filebuf* lsdp = file.rdbuf(); 137 std::filebuf* lsdp = file.rdbuf();
136 138
137 if (! (lsdp->is_open()) ) 139 if (! (lsdp->is_open()) )
138 { 140 {
139 G4String s1("Data file: "); 141 G4String s1("Data file: ");
140 G4String s2(" not found"); 142 G4String s2(" not found");
141 G4String excep = s1 + dirFile + s2; 143 G4String excep = s1 + dirFile + s2;
142 G4Exception("G4RDDopplerProfile::LoadBig 144 G4Exception("G4RDDopplerProfile::LoadBiggsP()",
143 "DataNotFound", FatalExcepti 145 "DataNotFound", FatalException, excep);
144 } 146 }
145 147
146 G4double p; 148 G4double p;
147 while(!file.eof()) 149 while(!file.eof())
148 { 150 {
149 file >> p; 151 file >> p;
150 biggsP.push_back(p); 152 biggsP.push_back(p);
151 } 153 }
152 154
153 // Make sure that the number of data loaded 155 // Make sure that the number of data loaded corresponds to the number in Biggs' paper
154 if (biggsP.size() != nBiggs) 156 if (biggsP.size() != nBiggs)
155 G4Exception("G4RDDopplerProfile::LoadBiggs 157 G4Exception("G4RDDopplerProfile::LoadBiggsP()", "InvalidCondition",
156 FatalException, "Number of mom 158 FatalException, "Number of momenta read in is not 31!");
157 } 159 }
158 160
159 161
160 void G4RDDopplerProfile::LoadProfile(const G4S 162 void G4RDDopplerProfile::LoadProfile(const G4String& fileName,G4int Z)
161 { 163 {
162 std::ostringstream ost; 164 std::ostringstream ost;
163 ost << fileName << "-" << Z << ".dat"; 165 ost << fileName << "-" << Z << ".dat";
164 G4String name(ost.str()); 166 G4String name(ost.str());
165 167
166 const char* path = G4FindDataDir("G4LEDATA") << 168 char* path = getenv("G4LEDATA");
167 if (!path) 169 if (!path)
168 { 170 {
169 G4String excep("G4LEDATA environment var 171 G4String excep("G4LEDATA environment variable not set!");
170 G4Exception("G4RDDopplerProfile::LoadPro 172 G4Exception("G4RDDopplerProfile::LoadProfile()",
171 "InvalidSetup", FatalExcepti 173 "InvalidSetup", FatalException, excep);
172 } 174 }
173 175
174 G4String pathString(path); 176 G4String pathString(path);
175 G4String dirFile = pathString + name; 177 G4String dirFile = pathString + name;
176 std::ifstream file(dirFile); 178 std::ifstream file(dirFile);
177 std::filebuf* lsdp = file.rdbuf(); 179 std::filebuf* lsdp = file.rdbuf();
178 180
179 if (! (lsdp->is_open()) ) 181 if (! (lsdp->is_open()) )
180 { 182 {
181 G4String s1("Data file: "); 183 G4String s1("Data file: ");
182 G4String s2(" not found"); 184 G4String s2(" not found");
183 G4String excep = s1 + dirFile + s2; 185 G4String excep = s1 + dirFile + s2;
184 G4Exception("G4RDDopplerProfile::LoadPro 186 G4Exception("G4RDDopplerProfile::LoadProfile()",
185 "DataNotFound", FatalExcepti 187 "DataNotFound", FatalException, excep);
186 } 188 }
187 189
188 G4double p; 190 G4double p;
189 G4int nShell = 0; 191 G4int nShell = 0;
190 192
191 // Create CompositeDataSet for the current Z 193 // Create CompositeDataSet for the current Z
192 G4RDVDataSetAlgorithm* interpolation = new G 194 G4RDVDataSetAlgorithm* interpolation = new G4RDLogLogInterpolation;
193 G4RDVEMDataSet* dataSetForZ = new G4RDCompos 195 G4RDVEMDataSet* dataSetForZ = new G4RDCompositeEMDataSet(interpolation,1.,1.,1,1);
194 196
195 while (!file.eof()) 197 while (!file.eof())
196 { 198 {
197 nShell++; 199 nShell++;
198 G4DataVector* profi = new G4DataVector; 200 G4DataVector* profi = new G4DataVector;
199 G4DataVector* biggs = new G4DataVector; 201 G4DataVector* biggs = new G4DataVector;
200 202
201 // Read in profile data for the current 203 // Read in profile data for the current shell
202 for (size_t i=0; i<nBiggs; i++) 204 for (size_t i=0; i<nBiggs; i++)
203 { 205 {
204 file >> p; 206 file >> p;
205 profi->push_back(p); 207 profi->push_back(p);
206 biggs->push_back(biggsP[i]); 208 biggs->push_back(biggsP[i]);
207 // if (i == 16) G4cout << "profile = " 209 // if (i == 16) G4cout << "profile = " << p << G4endl;
208 } 210 }
209 211
210 // Create G4RDEMDataSet for the current 212 // Create G4RDEMDataSet for the current shell
211 G4RDVDataSetAlgorithm* algo = interpolat 213 G4RDVDataSetAlgorithm* algo = interpolation->Clone();
212 G4RDVEMDataSet* dataSet = new G4RDEMData 214 G4RDVEMDataSet* dataSet = new G4RDEMDataSet(Z, biggs, profi, algo, 1., 1., true);
213 215
214 // Add current shell profile component t 216 // Add current shell profile component to G4RDCompositeEMDataSet for the current Z
215 dataSetForZ->AddComponent(dataSet); 217 dataSetForZ->AddComponent(dataSet);
216 } 218 }
217 219
218 // Fill in number of shells for the current 220 // Fill in number of shells for the current Z
219 nShells.push_back(nShell); 221 nShells.push_back(nShell);
220 222
221 profileMap[Z] = dataSetForZ; 223 profileMap[Z] = dataSetForZ;
222 } 224 }
223 225
224 226
225 G4double G4RDDopplerProfile::RandomSelectMome 227 G4double G4RDDopplerProfile::RandomSelectMomentum(G4int Z, G4int shellIndex) const
226 { 228 {
227 G4double value = 0.; 229 G4double value = 0.;
228 const G4RDVEMDataSet* profis = Profiles(Z); 230 const G4RDVEMDataSet* profis = Profiles(Z);
229 value = profis->RandomSelect(shellIndex); 231 value = profis->RandomSelect(shellIndex);
230 return value; 232 return value;
231 } 233 }
232 234