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