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