Geant4  9.6.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4NeutronHPVector.hh
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 // 070606 fix with Valgrind by T. Koi
27 // 080409 Fix div0 error with G4FPE by T. Koi
28 // 080811 Comment out unused method SetBlocked and SetBuffered
29 // Add required cleaning up in CleanUp by T. Koi
30 //
31 #ifndef G4NeutronHPVector_h
32 #define G4NeutronHPVector_h 1
33 
34 #include "G4NeutronHPDataPoint.hh"
35 #include "G4PhysicsVector.hh"
36 #include "G4NeutronHPInterpolator.hh"
37 #include "Randomize.hh"
38 #include "G4ios.hh"
39 #include <fstream>
40 #include "G4InterpolationManager.hh"
41 #include "G4NeutronHPInterpolator.hh"
42 #include "G4NeutronHPHash.hh"
43 #include <cmath>
44 #include <vector>
45 
46 #if defined WIN32-VC
47  #include <float.h>
48 #endif
49 
50 class G4NeutronHPVector
51 {
52  friend G4NeutronHPVector & operator + (G4NeutronHPVector & left,
53  G4NeutronHPVector & right);
54 
55  public:
56 
57  G4NeutronHPVector();
58 
59  G4NeutronHPVector(G4int n);
60 
61  ~G4NeutronHPVector();
62 
63  G4NeutronHPVector & operator = (const G4NeutronHPVector & right);
64 
65  inline void SetVerbose(G4int ff)
66  {
67  Verbose = ff;
68  }
69 
70  inline void Times(G4double factor)
71  {
72  G4int i;
73  for(i=0; i<nEntries; i++)
74  {
75  theData[i].SetY(theData[i].GetY()*factor);
76  }
77  if(theIntegral!=0)
78  {
79  theIntegral[i] *= factor;
80  }
81  }
82 
83  inline void SetPoint(G4int i, const G4NeutronHPDataPoint & it)
84  {
85  G4double x = it.GetX();
86  G4double y = it.GetY();
87  SetData(i, x, y);
88  }
89 
90  inline void SetData(G4int i, G4double x, G4double y)
91  {
92 // G4cout <<"G4NeutronHPVector::SetData called"<<nPoints<<" "<<nEntries<<G4endl;
93  Check(i);
94  if(y>maxValue) maxValue=y;
95  theData[i].SetData(x, y);
96  }
97  inline void SetX(G4int i, G4double e)
98  {
99  Check(i);
100  theData[i].SetX(e);
101  }
102  inline void SetEnergy(G4int i, G4double e)
103  {
104  Check(i);
105  theData[i].SetX(e);
106  }
107  inline void SetY(G4int i, G4double x)
108  {
109  Check(i);
110  if(x>maxValue) maxValue=x;
111  theData[i].SetY(x);
112  }
113  inline void SetXsec(G4int i, G4double x)
114  {
115  Check(i);
116  if(x>maxValue) maxValue=x;
117  theData[i].SetY(x);
118  }
119  inline G4double GetEnergy(G4int i) const { return theData[i].GetX(); }
120  inline G4double GetXsec(G4int i) { return theData[i].GetY(); }
121  inline G4double GetX(G4int i) const
122  {
123  if (i<0) i=0;
124  if(i>=GetVectorLength()) i=GetVectorLength()-1;
125  return theData[i].GetX();
126  }
127  inline const G4NeutronHPDataPoint & GetPoint(G4int i) const { return theData[i]; }
128 
129  void Hash()
130  {
131  G4int i;
132  G4double x, y;
133  for(i=0 ; i<nEntries; i++)
134  {
135  if(0 == (i+1)%10)
136  {
137  x = GetX(i);
138  y = GetY(i);
139  theHash.SetData(i, x, y);
140  }
141  }
142  }
143 
144  void ReHash()
145  {
146  theHash.Clear();
147  Hash();
148  }
149 
150  G4double GetXsec(G4double e);
151  G4double GetXsec(G4double e, G4int min)
152  {
153  G4int i;
154  for(i=min ; i<nEntries; i++)
155  {
156  if(theData[i].GetX()>e) break;
157  }
158  G4int low = i-1;
159  G4int high = i;
160  if(i==0)
161  {
162  low = 0;
163  high = 1;
164  }
165  else if(i==nEntries)
166  {
167  low = nEntries-2;
168  high = nEntries-1;
169  }
170  G4double y;
171  if(e<theData[nEntries-1].GetX())
172  {
173  // Protect against doubled-up x values
174  if( (theData[high].GetX()-theData[low].GetX())/theData[high].GetX() < 0.000001)
175  {
176  y = theData[low].GetY();
177  }
178  else
179  {
180  y = theInt.Interpolate(theManager.GetScheme(high), e,
181  theData[low].GetX(), theData[high].GetX(),
182  theData[low].GetY(), theData[high].GetY());
183  }
184  }
185  else
186  {
187  y=theData[nEntries-1].GetY();
188  }
189  return y;
190  }
191 
192  inline G4double GetY(G4double x) {return GetXsec(x);}
193  inline G4int GetVectorLength() const {return nEntries;}
194 
195  inline G4double GetY(G4int i)
196  {
197  if (i<0) i=0;
198  if(i>=GetVectorLength()) i=GetVectorLength()-1;
199  return theData[i].GetY();
200  }
201 
202  inline G4double GetY(G4int i) const
203  {
204  if (i<0) i=0;
205  if(i>=GetVectorLength()) i=GetVectorLength()-1;
206  return theData[i].GetY();
207  }
208  void Dump();
209 
210  inline void InitInterpolation(std::ifstream & aDataFile)
211  {
212  theManager.Init(aDataFile);
213  }
214 
215  void Init(std::ifstream & aDataFile, G4int total, G4double ux=1., G4double uy=1.)
216  {
217  G4double x,y;
218  for (G4int i=0;i<total;i++)
219  {
220  aDataFile >> x >> y;
221  x*=ux;
222  y*=uy;
223  SetData(i,x,y);
224  if(0 == nEntries%10)
225  {
226  theHash.SetData(nEntries-1, x, y);
227  }
228  }
229  }
230 
231  void Init(std::ifstream & aDataFile,G4double ux=1., G4double uy=1.)
232  {
233  G4int total;
234  aDataFile >> total;
235  if(theData!=0) delete [] theData;
236  theData = new G4NeutronHPDataPoint[total];
237  nPoints=total;
238  nEntries=0;
239  theManager.Init(aDataFile);
240  Init(aDataFile, total, ux, uy);
241  }
242 
243  void ThinOut(G4double precision);
244 
245  inline void SetLabel(G4double aLabel)
246  {
247  label = aLabel;
248  }
249 
250  inline G4double GetLabel()
251  {
252  return label;
253  }
254 
255  inline void CleanUp()
256  {
257  nEntries=0;
258  theManager.CleanUp();
259  maxValue = -DBL_MAX;
260  theHash.Clear();
261 //080811 TK DB
262  delete[] theIntegral;
263  theIntegral = NULL;
264  }
265 
266  // merges the vectors active and passive into *this
267  inline void Merge(G4NeutronHPVector * active, G4NeutronHPVector * passive)
268  {
269  CleanUp();
270  G4int s_tmp = 0, n=0, m_tmp=0;
271  G4NeutronHPVector * tmp;
272  G4int a = s_tmp, p = n, t;
273  while (a<active->GetVectorLength()&&p<passive->GetVectorLength())
274  {
275  if(active->GetEnergy(a) <= passive->GetEnergy(p))
276  {
277  G4double xa = active->GetEnergy(a);
278  G4double yy = active->GetXsec(a);
279  SetData(m_tmp, xa, yy);
280  theManager.AppendScheme(m_tmp, active->GetScheme(a));
281  m_tmp++;
282  a++;
283  G4double xp = passive->GetEnergy(p);
284 
285 //080409 TKDB
286  //if( std::abs(std::abs(xp-xa)/xa)<0.001 ) p++;
287  if ( !( xa == 0 ) && std::abs(std::abs(xp-xa)/xa)<0.001 ) p++;
288  } else {
289  tmp = active;
290  t=a;
291  active = passive;
292  a=p;
293  passive = tmp;
294  p=t;
295  }
296  }
297  while (a!=active->GetVectorLength())
298  {
299  SetData(m_tmp, active->GetEnergy(a), active->GetXsec(a));
300  theManager.AppendScheme(m_tmp++, active->GetScheme(a));
301  a++;
302  }
303  while (p!=passive->GetVectorLength())
304  {
305  if(std::abs(GetEnergy(m_tmp-1)-passive->GetEnergy(p))/passive->GetEnergy(p)>0.001)
306  //if(std::abs(GetEnergy(m)-passive->GetEnergy(p))/passive->GetEnergy(p)>0.001)
307  {
308  SetData(m_tmp, passive->GetEnergy(p), passive->GetXsec(p));
309  theManager.AppendScheme(m_tmp++, active->GetScheme(p));
310  }
311  p++;
312  }
313  }
314 
315  void Merge(G4InterpolationScheme aScheme, G4double aValue,
316  G4NeutronHPVector * active, G4NeutronHPVector * passive);
317 
318  G4double SampleLin() // Samples X according to distribution Y, linear int
319  {
320  G4double result;
321  if(theIntegral==0) IntegrateAndNormalise();
322  if(GetVectorLength()==1)
323  {
324  result = theData[0].GetX();
325  }
326  else
327  {
328  G4int i;
329  G4double rand = G4UniformRand();
330 
331  // this was replaced
332 // for(i=1;i<GetVectorLength();i++)
333 // {
334 // if(rand<theIntegral[i]/theIntegral[GetVectorLength()-1]) break;
335 // }
336 
337 // by this (begin)
338  for(i=GetVectorLength()-1; i>=0 ;i--)
339  {
340  if(rand>theIntegral[i]/theIntegral[GetVectorLength()-1]) break;
341  }
342  if(i!=GetVectorLength()-1) i++;
343 // until this (end)
344 
345  G4double x1, x2, y1, y2;
346  y1 = theData[i-1].GetX();
347  x1 = theIntegral[i-1];
348  y2 = theData[i].GetX();
349  x2 = theIntegral[i];
350  if(std::abs((y2-y1)/y2)<0.0000001) // not really necessary, since the case is excluded by construction
351  {
352  y1 = theData[i-2].GetX();
353  x1 = theIntegral[i-2];
354  }
355  result = theLin.Lin(rand, x1, x2, y1, y2);
356  }
357  return result;
358  }
359 
360  G4double Sample(); // Samples X according to distribution Y
361 
362  G4double * Debug()
363  {
364  return theIntegral;
365  }
366 
367  inline void IntegrateAndNormalise()
368  {
369  G4int i;
370  if(theIntegral!=0) return;
371  theIntegral = new G4double[nEntries];
372  if(nEntries == 1)
373  {
374  theIntegral[0] = 1;
375  return;
376  }
377  theIntegral[0] = 0;
378  G4double sum = 0;
379  G4double x1 = 0;
380  G4double x0 = 0;
381  for(i=1;i<GetVectorLength();i++)
382  {
383  x1 = theData[i].GetX();
384  x0 = theData[i-1].GetX();
385  if (std::abs(x1-x0) > std::abs(x1*0.0000001) )
386  {
387  //********************************************************************
388  //EMendoza -> the interpolation scheme is not always lin-lin
389  /*
390  sum+= 0.5*(theData[i].GetY()+theData[i-1].GetY())*
391  (x1-x0);
392  */
393  //********************************************************************
394  G4InterpolationScheme aScheme = theManager.GetScheme(i);
395  G4double y0 = theData[i-1].GetY();
396  G4double y1 = theData[i].GetY();
397  G4double integ=theInt.GetBinIntegral(aScheme,x0,x1,y0,y1);
398 #if defined WIN32-VC
399  if(!_finite(integ)){integ=0;}
400 #elif defined __IBMCPP__
401  if(isinf(integ)||isnan(integ)){integ=0;}
402 #else
403  if(std::isinf(integ)||std::isnan(integ)){integ=0;}
404 #endif
405  sum+=integ;
406  //********************************************************************
407  }
408  theIntegral[i] = sum;
409  }
410  G4double total = theIntegral[GetVectorLength()-1];
411  for(i=1;i<GetVectorLength();i++)
412  {
413  theIntegral[i]/=total;
414  }
415  }
416 
417  inline void Integrate()
418  {
419  G4int i;
420  if(nEntries == 1)
421  {
422  totalIntegral = 0;
423  return;
424  }
425  G4double sum = 0;
426  for(i=1;i<GetVectorLength();i++)
427  {
428  if(std::abs((theData[i].GetX()-theData[i-1].GetX())/theData[i].GetX())>0.0000001)
429  {
430  G4double x1 = theData[i-1].GetX();
431  G4double x2 = theData[i].GetX();
432  G4double y1 = theData[i-1].GetY();
433  G4double y2 = theData[i].GetY();
434  G4InterpolationScheme aScheme = theManager.GetScheme(i);
435  if(aScheme==LINLIN||aScheme==CLINLIN||aScheme==ULINLIN)
436  {
437  sum+= 0.5*(y2+y1)*(x2-x1);
438  }
439  else if(aScheme==LINLOG||aScheme==CLINLOG||aScheme==ULINLOG)
440  {
441  G4double a = y1;
442  G4double b = (y2-y1)/(std::log(x2)-std::log(x1));
443  sum+= (a-b)*(x2-x1) + b*(x2*std::log(x2)-x1*std::log(x1));
444  }
445  else if(aScheme==LOGLIN||aScheme==CLOGLIN||aScheme==ULOGLIN)
446  {
447  G4double a = std::log(y1);
448  G4double b = (std::log(y2)-std::log(y1))/(x2-x1);
449  sum += (std::exp(a)/b)*(std::exp(b*x2)-std::exp(b*x1));
450  }
451  else if(aScheme==HISTO||aScheme==CHISTO||aScheme==UHISTO)
452  {
453  sum+= y1*(x2-x1);
454  }
455  else if(aScheme==LOGLOG||aScheme==CLOGLOG||aScheme==ULOGLOG)
456  {
457  G4double a = std::log(y1);
458  G4double b = (std::log(y2)-std::log(y1))/(std::log(x2)-std::log(x1));
459  sum += (std::exp(a)/(b+1))*(std::pow(x2,b+1)-std::pow(x1,b+1));
460  }
461  else
462  {
463  throw G4HadronicException(__FILE__, __LINE__, "Unknown interpolation scheme in G4NeutronHPVector::Integrate");
464  }
465 
466  }
467  }
468  totalIntegral = sum;
469  }
470 
471  inline G4double GetIntegral() // linear interpolation; use with care
472  {
473  if(totalIntegral<-0.5) Integrate();
474  return totalIntegral;
475  }
476 
477  inline void SetInterpolationManager(const G4InterpolationManager & aManager)
478  {
479  theManager = aManager;
480  }
481 
482  inline const G4InterpolationManager & GetInterpolationManager() const
483  {
484  return theManager;
485  }
486 
487  inline void SetInterpolationManager(G4InterpolationManager & aMan)
488  {
489  theManager = aMan;
490  }
491 
492  inline void SetScheme(G4int aPoint, const G4InterpolationScheme & aScheme)
493  {
494  theManager.AppendScheme(aPoint, aScheme);
495  }
496 
497  inline G4InterpolationScheme GetScheme(G4int anIndex)
498  {
499  return theManager.GetScheme(anIndex);
500  }
501 
502  G4double GetMeanX()
503  {
504  G4double result;
505  G4double running = 0;
506  G4double weighted = 0;
507  for(G4int i=1; i<nEntries; i++)
508  {
509  running += theInt.GetBinIntegral(theManager.GetScheme(i-1),
510  theData[i-1].GetX(), theData[i].GetX(),
511  theData[i-1].GetY(), theData[i].GetY());
512  weighted += theInt.GetWeightedBinIntegral(theManager.GetScheme(i-1),
513  theData[i-1].GetX(), theData[i].GetX(),
514  theData[i-1].GetY(), theData[i].GetY());
515  }
516  result = weighted / running;
517  return result;
518  }
519 
520 /*
521  void Block(G4double aX)
522  {
523  theBlocked.push_back(aX);
524  }
525 
526  void Buffer(G4double aX)
527  {
528  theBuffered.push_back(aX);
529  }
530 */
531 
532  std::vector<G4double> GetBlocked() {return theBlocked;}
533  std::vector<G4double> GetBuffered() {return theBuffered;}
534 
535 // void SetBlocked(const std::vector<G4double> &aBlocked) {theBlocked = aBlocked;}
536 // void SetBuffered(const std::vector<G4double> &aBuffer) {theBuffered = aBuffer;}
537 
538  G4double Get15percentBorder();
539  G4double Get50percentBorder();
540 
541  private:
542 
543  void Check(G4int i);
544 
545  G4bool IsBlocked(G4double aX);
546 
547  private:
548 
549  G4NeutronHPInterpolator theLin;
550 
551  private:
552 
553  G4double totalIntegral;
554 
555  G4NeutronHPDataPoint * theData; // the data
556  G4InterpolationManager theManager; // knows how to interpolate the data.
557  G4double * theIntegral;
558  G4int nEntries;
559  G4int nPoints;
560  G4double label;
561 
562  G4NeutronHPInterpolator theInt;
563  G4int Verbose;
564  // debug only
565  G4int isFreed;
566 
567  G4NeutronHPHash theHash;
568  G4double maxValue;
569 
570  std::vector<G4double> theBlocked;
571  std::vector<G4double> theBuffered;
572  G4double the15percentBorderCash;
573  G4double the50percentBorderCash;
574 
575 };
576 
577 #endif
Generated on Sat May 25 2013 14:34:03 for Geant4 by   doxygen 1.8.4