Geant4_10
G4RDDopplerProfile.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$
28 // GEANT4 tag $Name: $
29 //
30 // Author: Maria Grazia Pia (Maria.Grazia.Pia@cern.ch)
31 //
32 // History:
33 // -----------
34 // 31 Jul 2001 MGP Created
35 //
36 // -------------------------------------------------------------------
37 
38 #include "G4RDDopplerProfile.hh"
39 #include "G4DataVector.hh"
40 #include "G4SystemOfUnits.hh"
41 #include "G4RDVEMDataSet.hh"
42 #include "G4RDEMDataSet.hh"
44 #include "G4RDVDataSetAlgorithm.hh"
46 
47 #include <fstream>
48 #include <sstream>
49 #include "Randomize.hh"
50 
51 // The following deprecated header is included because <functional> seems not to be found on MGP's laptop
52 //#include "function.h"
53 
54 // Constructor
55 
57  : zMin(minZ), zMax(maxZ)
58 {
59  nBiggs = 31;
60 
61  LoadBiggsP("/doppler/p-biggs");
62 
63  for (G4int Z=zMin; Z<zMax+1; Z++)
64  {
65  LoadProfile("/doppler/profile",Z);
66  }
67 }
68 
69 // Destructor
71 {
72  std::map<G4int,G4RDVEMDataSet*,std::less<G4int> >::iterator pos;
73  for (pos = profileMap.begin(); pos != profileMap.end(); ++pos)
74  {
75  G4RDVEMDataSet* dataSet = (*pos).second;
76  delete dataSet;
77  dataSet = 0;
78  }
79 }
80 
81 
83 {
84  G4int n = 0;
85  if (Z>= zMin && Z <= zMax) n = nShells[Z-1];
86  return n;
87 }
88 
89 
91 {
92  std::map<G4int,G4RDVEMDataSet*,std::less<G4int> >::const_iterator pos;
93  if (Z < zMin || Z > zMax)
94  G4Exception("G4RDDopplerProfile::Profiles()", "OutOfRange",
95  FatalException, "Z outside boundaries!");
96  pos = profileMap.find(Z);
97  G4RDVEMDataSet* dataSet = (*pos).second;
98  return dataSet;
99 }
100 
101 
103 {
104  const G4RDVEMDataSet* profis = Profiles(Z);
105  const G4RDVEMDataSet* profi = profis->GetComponent(shellIndex);
106  return profi;
107 }
108 
109 
111 {
112  for (G4int Z=zMin; Z<zMax; Z++)
113  {
114  const G4RDVEMDataSet* profis = Profiles(Z);
115  profis->PrintData();
116  }
117 }
118 
119 
120 void G4RDDopplerProfile::LoadBiggsP(const G4String& fileName)
121 {
122  std::ostringstream ost;
123  ost << fileName << ".dat";
124  G4String name(ost.str());
125 
126  char* path = getenv("G4LEDATA");
127  if (!path)
128  {
129  G4String excep("G4LEDATA environment variable not set!");
130  G4Exception("G4RDDopplerProfile::LoadBiggsP()",
131  "InvalidSetup", FatalException, excep);
132  }
133 
134  G4String pathString(path);
135  G4String dirFile = pathString + name;
136  std::ifstream file(dirFile);
137  std::filebuf* lsdp = file.rdbuf();
138 
139  if (! (lsdp->is_open()) )
140  {
141  G4String s1("Data file: ");
142  G4String s2(" not found");
143  G4String excep = s1 + dirFile + s2;
144  G4Exception("G4RDDopplerProfile::LoadBiggsP()",
145  "DataNotFound", FatalException, excep);
146  }
147 
148  G4double p;
149  while(!file.eof())
150  {
151  file >> p;
152  biggsP.push_back(p);
153  }
154 
155  // Make sure that the number of data loaded corresponds to the number in Biggs' paper
156  if (biggsP.size() != nBiggs)
157  G4Exception("G4RDDopplerProfile::LoadBiggsP()", "InvalidCondition",
158  FatalException, "Number of momenta read in is not 31!");
159 }
160 
161 
162 void G4RDDopplerProfile::LoadProfile(const G4String& fileName,G4int Z)
163 {
164  std::ostringstream ost;
165  ost << fileName << "-" << Z << ".dat";
166  G4String name(ost.str());
167 
168  char* path = getenv("G4LEDATA");
169  if (!path)
170  {
171  G4String excep("G4LEDATA environment variable not set!");
172  G4Exception("G4RDDopplerProfile::LoadProfile()",
173  "InvalidSetup", FatalException, excep);
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("G4RDDopplerProfile::LoadProfile()",
187  "DataNotFound", FatalException, excep);
188  }
189 
190  G4double p;
191  G4int nShell = 0;
192 
193  // Create CompositeDataSet for the current Z
194  G4RDVDataSetAlgorithm* interpolation = new G4RDLogLogInterpolation;
195  G4RDVEMDataSet* dataSetForZ = new G4RDCompositeEMDataSet(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 G4RDEMDataSet for the current shell
213  G4RDVDataSetAlgorithm* algo = interpolation->Clone();
214  G4RDVEMDataSet* dataSet = new G4RDEMDataSet(Z, biggs, profi, algo, 1., 1., true);
215 
216  // Add current shell profile component to G4RDCompositeEMDataSet 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 G4RDVEMDataSet* profis = Profiles(Z);
231  value = profis->RandomSelect(shellIndex);
232  return value;
233 }
G4RDDopplerProfile(G4int minZ=1, G4int maxZ=100)
virtual G4RDVDataSetAlgorithm * Clone() const =0
const G4RDVEMDataSet * Profile(G4int Z, G4int ShellIndex) const
const char * p
Definition: xmltok.h:285
const XML_Char * name
Definition: expat.h:151
virtual const G4RDVEMDataSet * GetComponent(G4int componentId) const =0
TFile * file
int G4int
Definition: G4Types.hh:78
size_t NumberOfProfiles(G4int Z) const
Char_t n[5]
const G4RDVEMDataSet * Profiles(G4int Z) const
Float_t Z
Definition: plot.C:39
virtual void PrintData(void) const =0
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4double RandomSelectMomentum(G4int Z, G4int shellIndex) const
const XML_Char int const XML_Char * value
Definition: expat.h:331
double G4double
Definition: G4Types.hh:76
virtual G4double RandomSelect(G4int componentId=0) const =0
virtual void AddComponent(G4RDVEMDataSet *dataSet)=0