| 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 // 27 //
28 // Author: Maria Grazia Pia (Maria.Grazia.Pia@ 28 // Author: Maria Grazia Pia (Maria.Grazia.Pia@cern.ch)
29 // 29 //
30 // History: 30 // History:
31 // ----------- 31 // -----------
32 // 31 Jul 2001 MGP Created 32 // 31 Jul 2001 MGP Created
33 // 33 //
34 // ------------------------------------------- 34 // -------------------------------------------------------------------
35 35
36 #include "G4DopplerProfile.hh" 36 #include "G4DopplerProfile.hh"
37 #include "G4DataVector.hh" 37 #include "G4DataVector.hh"
38 #include "G4SystemOfUnits.hh" 38 #include "G4SystemOfUnits.hh"
39 #include "G4VEMDataSet.hh" 39 #include "G4VEMDataSet.hh"
40 #include "G4EMDataSet.hh" 40 #include "G4EMDataSet.hh"
41 #include "G4CompositeEMDataSet.hh" 41 #include "G4CompositeEMDataSet.hh"
42 #include "G4VDataSetAlgorithm.hh" 42 #include "G4VDataSetAlgorithm.hh"
43 #include "G4LogLogInterpolation.hh" 43 #include "G4LogLogInterpolation.hh"
44 44
45 #include <fstream> 45 #include <fstream>
46 #include <sstream> 46 #include <sstream>
47 #include "Randomize.hh" 47 #include "Randomize.hh"
48 48
49 //....oooOO0OOooo........oooOO0OOooo........oo 49 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
50 50
51 G4DopplerProfile::G4DopplerProfile(G4int minZ, 51 G4DopplerProfile::G4DopplerProfile(G4int minZ, G4int maxZ)
52 : zMin(minZ), zMax(maxZ) 52 : zMin(minZ), zMax(maxZ)
53 { 53 {
54 LoadBiggsP("/doppler/p-biggs"); 54 LoadBiggsP("/doppler/p-biggs");
55 55
56 for (G4int Z=zMin; Z<zMax+1; Z++) 56 for (G4int Z=zMin; Z<zMax+1; Z++)
57 { 57 {
58 LoadProfile("/doppler/profile",Z); 58 LoadProfile("/doppler/profile",Z);
59 } 59 }
60 } 60 }
61 61
62 //....oooOO0OOooo........oooOO0OOooo........oo 62 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
63 63
64 // Destructor 64 // Destructor
65 G4DopplerProfile::~G4DopplerProfile() 65 G4DopplerProfile::~G4DopplerProfile()
66 { 66 {
67 for (auto& pos : profileMap) 67 for (auto& pos : profileMap)
68 { 68 {
69 G4VEMDataSet* dataSet = pos.second; 69 G4VEMDataSet* dataSet = pos.second;
70 delete dataSet; 70 delete dataSet;
71 dataSet = nullptr; 71 dataSet = nullptr;
72 } 72 }
73 } 73 }
74 74
75 //....oooOO0OOooo........oooOO0OOooo........oo 75 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
76 76
77 size_t G4DopplerProfile::NumberOfProfiles(G4in 77 size_t G4DopplerProfile::NumberOfProfiles(G4int Z) const
78 { 78 {
79 G4int n = 0; 79 G4int n = 0;
80 if (Z>= zMin && Z <= zMax) n = nShells[Z-1]; 80 if (Z>= zMin && Z <= zMax) n = nShells[Z-1];
81 return n; 81 return n;
82 } 82 }
83 83
84 //....oooOO0OOooo........oooOO0OOooo........oo 84 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
85 85
86 const G4VEMDataSet* G4DopplerProfile::Profiles 86 const G4VEMDataSet* G4DopplerProfile::Profiles(G4int Z) const
87 { 87 {
88 if (Z < zMin || Z > zMax) 88 if (Z < zMin || Z > zMax)
89 G4Exception("G4DopplerProfile::Profiles", 89 G4Exception("G4DopplerProfile::Profiles",
90 "em1005",FatalException,"Z outside bou 90 "em1005",FatalException,"Z outside boundaries");
91 auto pos = profileMap.find(Z); 91 auto pos = profileMap.find(Z);
92 G4VEMDataSet* dataSet = (*pos).second; 92 G4VEMDataSet* dataSet = (*pos).second;
93 return dataSet; 93 return dataSet;
94 } 94 }
95 95
96 //....oooOO0OOooo........oooOO0OOooo........oo 96 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
97 97
98 const G4VEMDataSet* G4DopplerProfile::Profile( 98 const G4VEMDataSet* G4DopplerProfile::Profile(G4int Z, G4int shellIndex) const
99 { 99 {
100 const G4VEMDataSet* profis = Profiles(Z); 100 const G4VEMDataSet* profis = Profiles(Z);
101 const G4VEMDataSet* profi = profis->GetCompo 101 const G4VEMDataSet* profi = profis->GetComponent(shellIndex);
102 return profi; 102 return profi;
103 } 103 }
104 104
105 //....oooOO0OOooo........oooOO0OOooo........oo 105 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
106 106
107 void G4DopplerProfile::PrintData() const 107 void G4DopplerProfile::PrintData() const
108 { 108 {
109 for (G4int Z=zMin; Z<zMax; Z++) 109 for (G4int Z=zMin; Z<zMax; Z++)
110 { 110 {
111 const G4VEMDataSet* profis = Profiles(Z) 111 const G4VEMDataSet* profis = Profiles(Z);
112 profis->PrintData(); 112 profis->PrintData();
113 } 113 }
114 } 114 }
115 115
116 //....oooOO0OOooo........oooOO0OOooo........oo 116 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
117 117
118 void G4DopplerProfile::LoadBiggsP(const G4Stri 118 void G4DopplerProfile::LoadBiggsP(const G4String& fileName)
119 { 119 {
120 std::ostringstream ost; 120 std::ostringstream ost;
121 ost << fileName << ".dat"; 121 ost << fileName << ".dat";
122 G4String name(ost.str()); 122 G4String name(ost.str());
123 123
124 const char* path = G4FindDataDir("G4LEDATA") << 124 char* path = std::getenv("G4LEDATA");
125 if (!path) 125 if (!path)
126 { 126 {
127 G4Exception("G4DopplerProfile::LoadBiggs 127 G4Exception("G4DopplerProfile::LoadBiggsP",
128 "em0006",FatalException,"G4LEDATA envi 128 "em0006",FatalException,"G4LEDATA environment variable not set");
129 return; 129 return;
130 } 130 }
131 131
132 G4String pathString(path); 132 G4String pathString(path);
133 G4String dirFile = pathString + name; 133 G4String dirFile = pathString + name;
134 std::ifstream file(dirFile); 134 std::ifstream file(dirFile);
135 std::filebuf* lsdp = file.rdbuf(); 135 std::filebuf* lsdp = file.rdbuf();
136 136
137 if (! (lsdp->is_open()) ) 137 if (! (lsdp->is_open()) )
138 { 138 {
139 G4String s1("data file: "); 139 G4String s1("data file: ");
140 G4String s2(" not found"); 140 G4String s2(" not found");
141 G4String excep = s1 + dirFile + s2; 141 G4String excep = s1 + dirFile + s2;
142 G4Exception("G4DopplerProfile::LoadBiggs 142 G4Exception("G4DopplerProfile::LoadBiggsP",
143 "em0003",FatalException,excep); 143 "em0003",FatalException,excep);
144 } 144 }
145 145
146 G4double p; 146 G4double p;
147 while(!file.eof()) 147 while(!file.eof())
148 { 148 {
149 file >> p; 149 file >> p;
150 biggsP.push_back(p); 150 biggsP.push_back(p);
151 } 151 }
152 152
153 // Make sure that the number of data loaded 153 // Make sure that the number of data loaded corresponds to the number in Biggs' paper
154 if (biggsP.size() != nBiggs) 154 if (biggsP.size() != nBiggs)
155 G4Exception("G4DopplerProfile::LoadBiggsP" 155 G4Exception("G4DopplerProfile::LoadBiggsP",
156 "em1006",FatalException,"Number of mom 156 "em1006",FatalException,"Number of momenta read in is not 31");
157 } 157 }
158 158
159 //....oooOO0OOooo........oooOO0OOooo........oo 159 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
160 160
161 void G4DopplerProfile::LoadProfile(const G4Str 161 void G4DopplerProfile::LoadProfile(const G4String& fileName,G4int Z)
162 { 162 {
163 std::ostringstream ost; 163 std::ostringstream ost;
164 ost << fileName << "-" << Z << ".dat"; 164 ost << fileName << "-" << Z << ".dat";
165 G4String name(ost.str()); 165 G4String name(ost.str());
166 166
167 const char* path = G4FindDataDir("G4LEDATA") << 167 char* path = std::getenv("G4LEDATA");
168 if (!path) 168 if (!path)
169 { 169 {
170 G4String excep("G4LEDATA environment var 170 G4String excep("G4LEDATA environment variable not set");
171 G4Exception("G4DopplerProfile::LoadProfi 171 G4Exception("G4DopplerProfile::LoadProfile",
172 "em0006",FatalException,excep); 172 "em0006",FatalException,excep);
173 return; 173 return;
174 } 174 }
175 175
176 G4String pathString(path); 176 G4String pathString(path);
177 G4String dirFile = pathString + name; 177 G4String dirFile = pathString + name;
178 std::ifstream file(dirFile); 178 std::ifstream file(dirFile);
179 std::filebuf* lsdp = file.rdbuf(); 179 std::filebuf* lsdp = file.rdbuf();
180 180
181 if (! (lsdp->is_open()) ) 181 if (! (lsdp->is_open()) )
182 { 182 {
183 G4String s1("data file: "); 183 G4String s1("data file: ");
184 G4String s2(" not found"); 184 G4String s2(" not found");
185 G4String excep = s1 + dirFile + s2; 185 G4String excep = s1 + dirFile + s2;
186 G4Exception("G4DopplerProfile::LoadProfi 186 G4Exception("G4DopplerProfile::LoadProfile",
187 "em0003",FatalException,excep); 187 "em0003",FatalException,excep);
188 } 188 }
189 189
190 G4double p; 190 G4double p;
191 G4int nShell = 0; 191 G4int nShell = 0;
192 192
193 // Create CompositeDataSet for the current Z 193 // Create CompositeDataSet for the current Z
194 G4VDataSetAlgorithm* interpolation = new G4L 194 G4VDataSetAlgorithm* interpolation = new G4LogLogInterpolation;
195 G4VEMDataSet* dataSetForZ = new G4CompositeE 195 G4VEMDataSet* dataSetForZ = new G4CompositeEMDataSet(interpolation,1.,1.,1,1);
196 196
197 while (!file.eof()) 197 while (!file.eof())
198 { 198 {
199 nShell++; 199 nShell++;
200 G4DataVector* profi = new G4DataVector; 200 G4DataVector* profi = new G4DataVector;
201 G4DataVector* biggs = new G4DataVector; 201 G4DataVector* biggs = new G4DataVector;
202 202
203 // Read in profile data for the current 203 // Read in profile data for the current shell
204 for (size_t i=0; i<nBiggs; i++) 204 for (size_t i=0; i<nBiggs; i++)
205 { 205 {
206 file >> p; 206 file >> p;
207 profi->push_back(p); 207 profi->push_back(p);
208 biggs->push_back(biggsP[i]); 208 biggs->push_back(biggsP[i]);
209 } 209 }
210 210
211 // Create G4EMDataSet for the current sh 211 // Create G4EMDataSet for the current shell
212 G4VDataSetAlgorithm* algo = interpolatio 212 G4VDataSetAlgorithm* algo = interpolation->Clone();
213 G4VEMDataSet* dataSet = new G4EMDataSet( 213 G4VEMDataSet* dataSet = new G4EMDataSet(Z, biggs, profi, algo, 1., 1., true);
214 214
215 // Add current shell profile component t 215 // Add current shell profile component to G4CompositeEMDataSet for the current Z
216 dataSetForZ->AddComponent(dataSet); 216 dataSetForZ->AddComponent(dataSet);
217 } 217 }
218 218
219 // Fill in number of shells for the current 219 // Fill in number of shells for the current Z
220 nShells.push_back(nShell); 220 nShells.push_back(nShell);
221 221
222 profileMap[Z] = dataSetForZ; 222 profileMap[Z] = dataSetForZ;
223 } 223 }
224 224
225 //....oooOO0OOooo........oooOO0OOooo........oo 225 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo....
226 226
227 G4double G4DopplerProfile::RandomSelectMoment 227 G4double G4DopplerProfile::RandomSelectMomentum(G4int Z, G4int shellIndex) const
228 { 228 {
229 G4double value = 0.; 229 G4double value = 0.;
230 const G4VEMDataSet* profis = Profiles(Z); 230 const G4VEMDataSet* profis = Profiles(Z);
231 value = profis->RandomSelect(shellIndex); 231 value = profis->RandomSelect(shellIndex);
232 return value; 232 return value;
233 } 233 }
234 234