Geant4  10.02.p02
G4DiscreteScatteringModel.cc
Go to the documentation of this file.
3 #include "G4AtomicShells.hh"
4 #include "G4Track.hh"
5 #include "G4Material.hh"
6 #include "G4UnitsTable.hh"
9 #include "G4ElementData.hh"
10 #include "G4Exp.hh"
11 #include "G4Log.hh"
12 
13 #include <iostream>
14 #include <string>
15 
16 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
17 
18 using namespace std;
19 using namespace CLHEP;
20 
23 
24 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
25 
27  : G4VEmModel("DiscrScat"), fParticleChange(nullptr),fAnalogModel("pwe"),
28  fNumAngles(iNumAngles), fLowEnergyLimit(2*keV)
29 {
30  SetHighEnergyLimit(100.*MeV);
32  if(IsMaster() && fCdf == nullptr) {
33  fCdf = new G4ElementData();
34  fTcs = new G4ElementData();
35  }
36 }
37 
38 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
39 
41 {
42  if(IsMaster()) {
43  delete fCdf;
44  delete fTcs;
45  fCdf = nullptr;
46  fTcs = nullptr;
47  }
48 }
49 
50 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
51 
53  const G4DataVector&)
54 {
55  if(nullptr == fParticleChange) {
57  }
58  if(!IsMaster()) { return; }
59 
60  G4cout << "G4DiscreteScatteringModel::Initialise start"<<G4endl;
61 
62  const G4int maxZ = 100;
63 
64  char *path = getenv("G4GBFPDATA");
65  if (!path)
66  {
67  G4Exception("G4DiscreteScatteringModel::Initialise","em0006",
68  FatalException,"G4GBFPDATA environment variable not set.");
69  return;
70  }
71 
72  std::ostringstream eFullFileName;
73  eFullFileName << path;
74 
75  G4ProductionCutsTable* theCoupleTable =
77 
78  G4int numOfCouples = theCoupleTable->GetTableSize();
79 
80  for(G4int i=0; i<numOfCouples; ++i)
81  {
82  const G4Material* material =
83  theCoupleTable->GetMaterialCutsCouple(i)->GetMaterial();
84  const G4ElementVector* theElementVector = material->GetElementVector();
85  G4int nelm = material->GetNumberOfElements();
86 
87  for (G4int j=0; j<nelm; ++j)
88  {
89  G4int Z = G4lrint((*theElementVector)[j]->GetZ());
90  if(Z < 1) { Z = 1; }
91  else if(Z > maxZ) { Z = maxZ; }
92  if(!fTcs->GetElementData(Z)) { ReadData(Z, path); }
93  }
94  }
95  G4cout << "G4DiscreteScatteringModel::Initialise completed"<<G4endl;
96 }
97 
98 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
99 
101 {
102 
103  // Set strings with paths to tcs data file. Currently not a variable but this
104  // this will change
105  G4String fullFileName(argFileName);
106 
107  stringstream ss1;//create a stringstream
108  ss1 << Z;//add number to the stream
109  stringstream ss2;//create a stringstream
110  ss2 << fNumAngles;
111  G4String tcsPath = fullFileName+"/log_gbfp_"+fAnalogModel
112  +"_tcs_"+ss1.str()+"_"+ss2.str()+".dat";
113 
114  // Do some error checking and exiting if data does not exist
115  std::ifstream in(tcsPath, std::ifstream::binary|std::ifstream::in);
116  if (!in.is_open())
117  {
118  G4String message("Data file \"");
119  message+=tcsPath;
120  message+="\" not found";
121  G4Exception("G4DiscreteScatteringModel::LoadData","em0003",
122  FatalException,message);
123  return;
124  }
125  //G4cout<<"Reading in data file "<< tcsPath<<G4endl;
126 
127  // Create a temporary G4PhysicsVector object pointer
128  G4PhysicsVector* tempData = new G4PhysicsVector(false);
129 
130  // Use retrieve to read in the total cross section (tcs) data
131  tempData->Retrieve(in,true);
132 
133  // Convert tcs from cm^2 to mm^2
134  //tempData->ScaleVector(1.0, 100.0);
135 
136  // store pass this data to the tcs object and initialise for current element
137  fTcs->InitialiseForElement(Z,tempData);
138  in.close();
139 
140  // Set strings with paths to cdf data files. Currently not a variable
141  // but this this will change
142  G4String cdfPath = fullFileName+"/gbfp_"+fAnalogModel+"_cdf_"
143  +ss1.str()+"_"+ss2.str()+".dat";
144 
145  // Do some error checking and exiting if data does not exist
146  std::ifstream in2(cdfPath, std::ifstream::binary|std::ifstream::in);
147  if (!in2.is_open())
148  {
149  G4String message("Data file \"");
150  message+=cdfPath;
151  message+="\" not found";
152  G4Exception("G4DiscreteScatteringModel::LoadData","em0003",
153  FatalException,message);
154  return;
155  }
156 
157  //G4cout<<"Reading in data file "<< cdfPath<<G4endl;
158 
159  // The cumulative distribution functions (CDF) for energy E_j and X_i
160  // on (0,1) is C(E_j,X_i) and read-in/stored at this time.
161  // For the purposes of this model, each energy grid point where the CDF
162  // is evaluated is considered a component. The number of energy grid
163  // points is consistent with the fTcs data, so the following int is set
164  // by calling G4PhysicsVector::GetVectorLength().
165  G4int numEnergies = fTcs->GetElementData(Z)->GetVectorLength();
166 
167  // The ElementData object pointer is then initialized by
168  fCdf->InitialiseForComponent(Z, numEnergies);
169 
170  // Now the data files are read in for all energies. At each energy,
171  // there are fNumAngles angles and fNumAngles CDF values.
172 
173  std::vector<G4PhysicsVector*> tempDataCDF;
174  // Loop through each energy
175  for (int j=0; j<numEnergies; j++)
176  {
177  // Push back a new G4PhysicsVector for the jth energy
178  tempDataCDF.push_back(new G4PhysicsVector(false));
179 
180  // For use with David's PhysicsVector class
181  //tempDataCDF.push_back(new G4PhysicsVector(false,false,true));
182 
183  //tempDataCDF.push_back(new G4PhysicsVector(false,false,false));
184  //tempDataCDF.push_back(new G4PhysicsVector());
185  }
186 
187  // Open a temporary stream. The data for the jth energy group is copied
188  // to the temp file and then the temp file is sent to
189  // G4PhysicsVector::Retrieve(). Once the data is stored in the
190  // G4PhysicsVector, tempDataCDF, it is then passed to the ElementData, cdf,
191  // which is the container for all of the data for each element and energy.
192 
193  //static G4Mutex m = G4MUTEX_INITIALIZER;
194  //G4AutoLock l(&m);
195 
196  // Open the stream and call file "tempDataFile.dat"
197  std::ofstream file("tempDataFile.dat",
198  std::fstream::out | std::fstream::trunc);
199 
200  // Write to file the lower/upper bounds and the number of grid points
201  // for each column of data. The data in tempfile is 2xfNumAngles,
202  // hence fNumAngles fNumAngles.
203  // file<<lowerBound<<" "<<upperBound<<" "<<fNumAngles<<" "
204  // <<fNumAngles<<G4endl;
205  file<<"-1. 1. "<<fNumAngles<<" "<<fNumAngles<<G4endl;
206 
207  // Start while loop over the entire data file opened above
208  // e.g. pwe_cdf_79.dat. This data contains all of the data for each energy
209  G4int cntr = 0;
210  G4int j = 0;
211  G4double temp1, temp2;
212  while(in2>>temp1>>temp2)
213  {
214  // Write data to temporary file
215  file<<setprecision(16)<<temp1<<" "<<temp2<<G4endl;
216  cntr++;
217 
218  // When the first fNumAngles data points are copied to tempDataFile.dat,
219  // store data in tempDataCDF[j]. Then increment the energy index, j,
220  // close and clear the streams, and then reopen the streams such that
221  // the next fNumAngles data points are copied to tempDataFile.dat.
222  if (cntr==fNumAngles)
223  {
224  cntr=0;
225  std::ifstream inTemp("tempDataFile.dat",
226  std::ifstream::binary|std::ifstream::in);
227  tempDataCDF[j]->Retrieve(inTemp,true);
228  fCdf->AddComponent(Z,j,tempDataCDF[j]);
229  j++;
230  inTemp.close(), inTemp.clear(), file.close(), file.clear();
231  file.open("tempDataFile.dat",std::fstream::out | std::fstream::trunc);
232  file<<"-1. 1. "<<fNumAngles<<" "<<fNumAngles<<G4endl;
233  }
234  }
235 
236  in2.close();
237 
238  return;
239 }
240 
241 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
242 
244  std::vector<G4DynamicParticle*>*,
245  const G4MaterialCutsCouple* couple,
246  const G4DynamicParticle* p,
247  G4double cutEnergy,
248  G4double)
249 {
250  G4double E = p->GetKineticEnergy();
251 
252  if(E < fLowEnergyLimit) return;
253 
254  // Select random atom
256  E,cutEnergy,E)->GetZ());
257 
258  // Determine the energy bin
259  G4double logE = G4Log(E);
260  G4ThreeVector dir = p->GetMomentumDirection(); //old direction
261  G4int j = fTcs->GetElementData(Z)->FindBin(logE,0);
262 
263  //-------------------------------------------------------------------------
264  // it would be nice to have the following block of code in G4PhysicsVector.
265  // it could be a simple function.
266 
267  // v.GetVectorLength() - number of points
268  // v.Energy(size_t idx) - value x(i)
269  // v[i] - value y(i)
270 
271  // This is a monte carlo interpolation scheme
273  G4double e2 = fTcs->GetElementData(Z)->Energy(j+1);
274 
275  // This is a monte carlo interpolation scheme
276  G4double pie1 = (logE-e1)/(e2-e1);
277  G4double r = G4UniformRand();
278  if (r<pie1){ ++j; }
279 
280  //-------------------------------------------------------------------------
281  // it would be nice to have the following block of code in G4PhysicsVector.
282  //---------------
283  // Given the energy grid value associated with the DCS,
284  // sample a deflection cosine
285  r = G4UniformRand();
286  G4int k = -1;
287  // First test if the angle is the most probable angle
288  if (r>(*fCdf->GetComponentDataByIndex(Z,j)).Energy(fNumAngles-2))
289  { k = fNumAngles-1; }
290  // or the least probable... not sure why I do this (maybe because
291  // it is a simple check)
292  else if (r<=(*fCdf->GetComponentDataByIndex(Z,j)).Energy(0)) { k = 0; }
293  // if neither then loop through remaining angles, break when locating angle
294  else {
295  for (G4int i=fNumAngles-2; i>0; --i) {
296  if ( (r>(*fCdf->GetComponentDataByIndex(Z,j)).Energy(i-1))
297  && (r<= (*fCdf->GetComponentDataByIndex(Z,j)).Energy(i)) )
298  { k=i; break;}
299  }
300  }
301  // Throw an error if an angle was not sampled, data is probably no good
302  if(k<0)
303  {
304  G4cout << "G4DiscreteScatteringModel::SampleSecondaries():"
305  << " CDF was not inverted properly "<<k<<G4endl;
306  for (G4int i=0;i<fNumAngles;++i) {
307  G4cout<<i<<" "<<(*fCdf->GetComponentDataByIndex(Z,j)).Energy(i)<<" "<<r<<G4endl;
308  }
309  }
310  //-------------------------------------------------------------------------
311 
312  // Otherwise, go get the angle and pass it too local method GetNewDirection.
313  // Then do transformation and update fParticleChange.
314 
315  //G4cout<<"------------------"<<G4endl;
316  //G4cout<<"Sample Secondaries"<<G4endl;
317  //G4cout<<G4Exp(e1)<<" "<<E<<" "<<G4Exp(e2)<<" "<<pie1<<" "
318  // <<(*fCdf->GetComponentDataByIndex(Z,j)).Energy(k)<<G4endl;
319  //G4cout<<"------------------"<<G4endl;
320  //G4cout<<" "<<G4endl;
321 
322  G4ThreeVector newDirection =
324  newDirection.rotateUz(dir);
326 }
327 
328 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
329 
333 {
334  // super simple, very compact look-up for total cross section
335  if (E==0.) {return 1e30;}
336  //G4cout<<"------------------"<<G4endl;
337  //G4cout<<"ComputeCrossSectionPerAtom"<<G4endl;
338  //G4cout<< E<<" "<<G4Log(E)<<" "<< G4Exp(fTcs->GetValueForElement(Z,log(E)))
339  //<<" "<<fTcs->GetValueForElement(Z,log(E))<<G4endl;
340  //G4cout<<"------------------"<<G4endl;
341  //G4cout<<" "<<G4endl;
342 
343  //G4cout<< E<<" "<<G4Exp(fTcs->GetValueForElement(Z,log(E)))<<G4endl;
344 
345  return G4Exp(fTcs->GetValueForElement(Z,G4Log(E)));
346 }
347 
348 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
349 
351 {
352  G4ThreeVector dir(0.0,0.0,1.0);
353 
354  G4double sint = sin(acos(z1));
355  G4double cost = sqrt(1.0 - sint*sint);
356  G4double phi = twopi* G4UniformRand();
357  G4double dirx = sint*cos(phi);
358  G4double diry = sint*sin(phi);
359  G4double dirz = cost;
360 
361  dir.set(dirx,diry,dirz);
362  return dir;
363 }
364 
365 //....oooOO0OOooo........oooOO0OOooo........oooOO0OOooo........oooOO0OOooo......
366 
367 
static const double MeV
Definition: G4SIunits.hh:211
G4PhysicsVector * GetComponentDataByIndex(G4int Z, size_t idx)
void ReadData(G4int, const G4String &argFileName)
void AddComponent(G4int Z, G4int id, G4PhysicsVector *v)
std::vector< G4Element * > G4ElementVector
G4double GetKineticEnergy() const
CLHEP::Hep3Vector G4ThreeVector
G4ThreeVector GetNewDirection(G4double z1)
G4PhysicsVector * GetElementData(G4int Z)
static const G4double e2
G4bool IsMaster() const
Definition: G4VEmModel.hh:718
G4ParticleDefinition * GetDefinition() const
G4DiscreteScatteringModel(G4int iNumAngles=1)
size_t GetVectorLength() const
const G4ElementVector * GetElementVector() const
Definition: G4Material.hh:190
int G4int
Definition: G4Types.hh:78
void ProposeMomentumDirection(G4double Px, G4double Py, G4double Pz)
void SetHighEnergyLimit(G4double)
Definition: G4VEmModel.hh:725
virtual G4double ComputeCrossSectionPerAtom(const G4ParticleDefinition *, G4double E, G4double Z, G4double A=0., G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
virtual void Initialise(const G4ParticleDefinition *, const G4DataVector &)
void InitialiseForElement(G4int Z, G4PhysicsVector *v)
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
const G4ThreeVector & GetMomentumDirection() const
virtual void SampleSecondaries(std::vector< G4DynamicParticle * > *, const G4MaterialCutsCouple *, const G4DynamicParticle *, G4double tmin=0.0, G4double maxEnergy=DBL_MAX)
static const double twopi
Definition: G4SIunits.hh:75
size_t FindBin(G4double energy, size_t idx) const
G4double Energy(size_t index) const
static const G4double e1
void G4Exception(const char *originOfException, const char *exceptionCode, G4ExceptionSeverity severity, const char *comments)
Definition: G4Exception.cc:41
G4double G4Log(G4double x)
Definition: G4Log.hh:230
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
static G4ProductionCutsTable * GetProductionCutsTable()
const G4MaterialCutsCouple * GetMaterialCutsCouple(G4int i) const
int G4lrint(double ad)
Definition: templates.hh:163
G4double GetValueForElement(G4int Z, G4double kinEnergy)
G4ParticleChangeForGamma * fParticleChange
virtual G4bool Retrieve(std::ifstream &fIn, G4bool ascii=false)
#define G4endl
Definition: G4ios.hh:61
size_t GetNumberOfElements() const
Definition: G4Material.hh:186
static const double keV
Definition: G4SIunits.hh:213
double G4double
Definition: G4Types.hh:76
void SetLowEnergyLimit(G4double)
Definition: G4VEmModel.hh:732
void InitialiseForComponent(G4int Z, G4int nComponents=0)
const G4Element * SelectRandomAtom(const G4MaterialCutsCouple *, const G4ParticleDefinition *, G4double kineticEnergy, G4double cutEnergy=0.0, G4double maxEnergy=DBL_MAX)
Definition: G4VEmModel.hh:544
G4ParticleChangeForGamma * GetParticleChangeForGamma()
Definition: G4VEmModel.cc:134
const G4Material * GetMaterial() const