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