Geant4  10.03.p01
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
G4ParticleHPVector Class Reference

#include <G4ParticleHPVector.hh>

Public Member Functions

 G4ParticleHPVector ()
 
 G4ParticleHPVector (G4int n)
 
 ~G4ParticleHPVector ()
 
G4ParticleHPVectoroperator= (const G4ParticleHPVector &right)
 
void SetVerbose (G4int ff)
 
void Times (G4double factor)
 
void SetPoint (G4int i, const G4ParticleHPDataPoint &it)
 
void SetData (G4int i, G4double x, G4double y)
 
void SetX (G4int i, G4double e)
 
void SetEnergy (G4int i, G4double e)
 
void SetY (G4int i, G4double x)
 
void SetXsec (G4int i, G4double x)
 
G4double GetEnergy (G4int i) const
 
G4double GetXsec (G4int i)
 
G4double GetX (G4int i) const
 
const G4ParticleHPDataPointGetPoint (G4int i) const
 
void Hash ()
 
void ReHash ()
 
G4double GetXsec (G4double e)
 
G4double GetXsec (G4double e, G4int min)
 
G4double GetY (G4double x)
 
G4int GetVectorLength () const
 
G4double GetY (G4int i)
 
G4double GetY (G4int i) const
 
void Dump ()
 
void InitInterpolation (std::istream &aDataFile)
 
void Init (std::istream &aDataFile, G4int total, G4double ux=1., G4double uy=1.)
 
void Init (std::istream &aDataFile, G4double ux=1., G4double uy=1.)
 
void ThinOut (G4double precision)
 
void SetLabel (G4double aLabel)
 
G4double GetLabel ()
 
void CleanUp ()
 
void Merge (G4ParticleHPVector *active, G4ParticleHPVector *passive)
 
void Merge (G4InterpolationScheme aScheme, G4double aValue, G4ParticleHPVector *active, G4ParticleHPVector *passive)
 
G4double SampleLin ()
 
G4double Sample ()
 
G4doubleDebug ()
 
void IntegrateAndNormalise ()
 
void Integrate ()
 
G4double GetIntegral ()
 
void SetInterpolationManager (const G4InterpolationManager &aManager)
 
const G4InterpolationManagerGetInterpolationManager () const
 
void SetInterpolationManager (G4InterpolationManager &aMan)
 
void SetScheme (G4int aPoint, const G4InterpolationScheme &aScheme)
 
G4InterpolationScheme GetScheme (G4int anIndex)
 
G4double GetMeanX ()
 
std::vector< G4doubleGetBlocked ()
 
std::vector< G4doubleGetBuffered ()
 
G4double Get15percentBorder ()
 
G4double Get50percentBorder ()
 

Friends

G4ParticleHPVectoroperator+ (G4ParticleHPVector &left, G4ParticleHPVector &right)
 

Detailed Description

Definition at line 56 of file G4ParticleHPVector.hh.

Constructor & Destructor Documentation

G4ParticleHPVector::G4ParticleHPVector ( )

Definition at line 83 of file G4ParticleHPVector.cc.

84  {
85  theData = new G4ParticleHPDataPoint[20];
86  nPoints=20;
87  nEntries=0;
88  Verbose=0;
89  theIntegral=0;
90  totalIntegral=-1;
91  isFreed = 0;
92  maxValue = -DBL_MAX;
93  the15percentBorderCash = -DBL_MAX;
94  the50percentBorderCash = -DBL_MAX;
95  label = -DBL_MAX;
96  }
#define DBL_MAX
Definition: templates.hh:83
G4ParticleHPVector::G4ParticleHPVector ( G4int  n)

Definition at line 98 of file G4ParticleHPVector.cc.

99  {
100  nPoints=std::max(n, 20);
101  theData = new G4ParticleHPDataPoint[nPoints];
102  nEntries=0;
103  Verbose=0;
104  theIntegral=0;
105  totalIntegral=-1;
106  isFreed = 0;
107  maxValue = -DBL_MAX;
108  the15percentBorderCash = -DBL_MAX;
109  the50percentBorderCash = -DBL_MAX;
110  label = -DBL_MAX;
111  }
T max(const T t1, const T t2)
brief Return the largest of the two arguments
#define DBL_MAX
Definition: templates.hh:83

Here is the call graph for this function:

G4ParticleHPVector::~G4ParticleHPVector ( )

Definition at line 113 of file G4ParticleHPVector.cc.

114  {
115 // if(Verbose==1)G4cout <<"G4ParticleHPVector::~G4ParticleHPVector"<<G4endl;
116  delete [] theData;
117 // if(Verbose==1)G4cout <<"Vector: delete theData"<<G4endl;
118  delete [] theIntegral;
119 // if(Verbose==1)G4cout <<"Vector: delete theIntegral"<<G4endl;
120  theHash.Clear();
121  isFreed = 1;
122  }

Here is the call graph for this function:

Member Function Documentation

void G4ParticleHPVector::CleanUp ( )
inline

Definition at line 261 of file G4ParticleHPVector.hh.

262  {
263  nEntries=0;
264  theManager.CleanUp();
265  maxValue = -DBL_MAX;
266  theHash.Clear();
267 //080811 TK DB
268  delete[] theIntegral;
269  theIntegral = NULL;
270  }
#define DBL_MAX
Definition: templates.hh:83

Here is the call graph for this function:

Here is the caller graph for this function:

G4double* G4ParticleHPVector::Debug ( )
inline

Definition at line 368 of file G4ParticleHPVector.hh.

369  {
370  return theIntegral;
371  }
void G4ParticleHPVector::Dump ( )

Definition at line 205 of file G4ParticleHPVector.cc.

206  {
207  G4cout << nEntries<<G4endl;
208  for(G4int i=0; i<nEntries; i++)
209  {
210  G4cout << theData[i].GetX()<<" ";
211  G4cout << theData[i].GetY()<<" ";
212 // if (i!=1&&i==5*(i/5)) G4cout << G4endl;
213  G4cout << G4endl;
214  }
215  G4cout << G4endl;
216  }
int G4int
Definition: G4Types.hh:78
G4GLOB_DLL std::ostream G4cout
#define G4endl
Definition: G4ios.hh:61

Here is the call graph for this function:

G4double G4ParticleHPVector::Get15percentBorder ( )

Definition at line 480 of file G4ParticleHPVector.cc.

481  {
482  if(the15percentBorderCash>-DBL_MAX/2.) return the15percentBorderCash;
484  if(GetVectorLength()==1)
485  {
486  result = theData[0].GetX();
487  the15percentBorderCash = result;
488  }
489  else
490  {
491  if(theIntegral==0) { IntegrateAndNormalise(); }
492  G4int i;
493  result = theData[GetVectorLength()-1].GetX();
494  for(i=0;i<GetVectorLength();i++)
495  {
496  if(theIntegral[i]/theIntegral[GetVectorLength()-1]>0.15)
497  {
498  result = theData[std::min(i+1, GetVectorLength()-1)].GetX();
499  the15percentBorderCash = result;
500  break;
501  }
502  }
503  the15percentBorderCash = result;
504  }
505  return result;
506  }
G4double G4ParticleHPJENDLHEData::G4double result
G4int GetVectorLength() const
int G4int
Definition: G4Types.hh:78
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
double G4double
Definition: G4Types.hh:76
#define DBL_MAX
Definition: templates.hh:83

Here is the call graph for this function:

G4double G4ParticleHPVector::Get50percentBorder ( )

Definition at line 508 of file G4ParticleHPVector.cc.

509  {
510  if(the50percentBorderCash>-DBL_MAX/2.) return the50percentBorderCash;
512  if(GetVectorLength()==1)
513  {
514  result = theData[0].GetX();
515  the50percentBorderCash = result;
516  }
517  else
518  {
519  if(theIntegral==0) { IntegrateAndNormalise(); }
520  G4int i;
521  G4double x = 0.5;
522  result = theData[GetVectorLength()-1].GetX();
523  for(i=0;i<GetVectorLength();i++)
524  {
525  if(theIntegral[i]/theIntegral[GetVectorLength()-1]>x)
526  {
527  G4int it;
528  it = i;
529  if(it == GetVectorLength()-1)
530  {
531  result = theData[GetVectorLength()-1].GetX();
532  }
533  else
534  {
535  G4double x1, x2, y1, y2;
536  x1 = theIntegral[i-1]/theIntegral[GetVectorLength()-1];
537  x2 = theIntegral[i]/theIntegral[GetVectorLength()-1];
538  y1 = theData[i-1].GetX();
539  y2 = theData[i].GetX();
540  result = theLin.Lin(x, x1, x2, y1, y2);
541  }
542  the50percentBorderCash = result;
543  break;
544  }
545  }
546  the50percentBorderCash = result;
547  }
548  return result;
549  }
G4double G4ParticleHPJENDLHEData::G4double result
G4int GetVectorLength() const
G4double Lin(G4double x, G4double x1, G4double x2, G4double y1, G4double y2)
int G4int
Definition: G4Types.hh:78
double G4double
Definition: G4Types.hh:76
#define DBL_MAX
Definition: templates.hh:83

Here is the call graph for this function:

std::vector<G4double> G4ParticleHPVector::GetBlocked ( )
inline

Definition at line 538 of file G4ParticleHPVector.hh.

538 {return theBlocked;}
std::vector<G4double> G4ParticleHPVector::GetBuffered ( )
inline

Definition at line 539 of file G4ParticleHPVector.hh.

539 {return theBuffered;}
G4double G4ParticleHPVector::GetEnergy ( G4int  i) const
inline

Definition at line 125 of file G4ParticleHPVector.hh.

125 { return theData[i].GetX(); }

Here is the call graph for this function:

Here is the caller graph for this function:

G4double G4ParticleHPVector::GetIntegral ( )
inline

Definition at line 477 of file G4ParticleHPVector.hh.

478  {
479  if(totalIntegral<-0.5) Integrate();
480  return totalIntegral;
481  }

Here is the call graph for this function:

Here is the caller graph for this function:

const G4InterpolationManager& G4ParticleHPVector::GetInterpolationManager ( ) const
inline

Definition at line 488 of file G4ParticleHPVector.hh.

489  {
490  return theManager;
491  }

Here is the caller graph for this function:

G4double G4ParticleHPVector::GetLabel ( )
inline

Definition at line 256 of file G4ParticleHPVector.hh.

257  {
258  return label;
259  }

Here is the caller graph for this function:

G4double G4ParticleHPVector::GetMeanX ( )
inline

Definition at line 508 of file G4ParticleHPVector.hh.

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());
518  weighted += theInt.GetWeightedBinIntegral(theManager.GetScheme(i-1),
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  }
G4double G4ParticleHPJENDLHEData::G4double result
int G4int
Definition: G4Types.hh:78
G4InterpolationScheme GetScheme(G4int index) const
G4double GetWeightedBinIntegral(const G4InterpolationScheme &aScheme, const G4double x1, const G4double x2, const G4double y1, const G4double y2)
G4double GetBinIntegral(const G4InterpolationScheme &aScheme, const G4double x1, const G4double x2, const G4double y1, const G4double y2)
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

Here is the caller graph for this function:

const G4ParticleHPDataPoint& G4ParticleHPVector::GetPoint ( G4int  i) const
inline

Definition at line 133 of file G4ParticleHPVector.hh.

133 { return theData[i]; }

Here is the caller graph for this function:

G4InterpolationScheme G4ParticleHPVector::GetScheme ( G4int  anIndex)
inline

Definition at line 503 of file G4ParticleHPVector.hh.

504  {
505  return theManager.GetScheme(anIndex);
506  }
G4InterpolationScheme GetScheme(G4int index) const

Here is the call graph for this function:

Here is the caller graph for this function:

G4int G4ParticleHPVector::GetVectorLength ( ) const
inline

Definition at line 199 of file G4ParticleHPVector.hh.

199 {return nEntries;}

Here is the caller graph for this function:

G4double G4ParticleHPVector::GetX ( G4int  i) const
inline

Definition at line 127 of file G4ParticleHPVector.hh.

128  {
129  if (i<0) i=0;
130  if(i>=GetVectorLength()) i=GetVectorLength()-1;
131  return theData[i].GetX();
132  }
G4int GetVectorLength() const

Here is the call graph for this function:

Here is the caller graph for this function:

G4double G4ParticleHPVector::GetXsec ( G4int  i)
inline

Definition at line 126 of file G4ParticleHPVector.hh.

126 { return theData[i].GetY(); }

Here is the call graph for this function:

Here is the caller graph for this function:

G4double G4ParticleHPVector::GetXsec ( G4double  e)

Definition at line 149 of file G4ParticleHPVector.cc.

150  {
151  if(nEntries == 0) return 0;
152  //if(!theHash.Prepared()) Hash();
153  if ( !theHash.Prepared() ) {
154  if ( G4Threading::IsWorkerThread() ) {
155  ;
156  } else {
157  Hash();
158  }
159  }
160  G4int min = theHash.GetMinIndex(e);
161  G4int i;
162  for(i=min ; i<nEntries; i++)
163  {
164  //if(theData[i].GetX()>e) break;
165  if(theData[i].GetX() >= e) break;
166  }
167  G4int low = i-1;
168  G4int high = i;
169  if(i==0)
170  {
171  low = 0;
172  high = 1;
173  }
174  else if(i==nEntries)
175  {
176  low = nEntries-2;
177  high = nEntries-1;
178  }
179  G4double y;
180  if(e<theData[nEntries-1].GetX())
181  {
182  // Protect against doubled-up x values
183  //if( (theData[high].GetX()-theData[low].GetX())/theData[high].GetX() < 0.000001)
184  if ( theData[high].GetX() !=0
185  //080808 TKDB
186  //&&( theData[high].GetX()-theData[low].GetX())/theData[high].GetX() < 0.000001)
187  &&( std::abs( (theData[high].GetX()-theData[low].GetX())/theData[high].GetX() ) < 0.000001 ) )
188  {
189  y = theData[low].GetY();
190  }
191  else
192  {
193  y = theInt.Interpolate(theManager.GetScheme(high), e,
194  theData[low].GetX(), theData[high].GetX(),
195  theData[low].GetY(), theData[high].GetY());
196  }
197  }
198  else
199  {
200  y=theData[nEntries-1].GetY();
201  }
202  return y;
203  }
G4bool Prepared() const
int G4int
Definition: G4Types.hh:78
G4InterpolationScheme GetScheme(G4int index) const
G4double GetX(G4int i) const
G4double Interpolate(G4InterpolationScheme aScheme, G4double x, G4double x1, G4double x2, G4double y1, G4double y2) const
G4bool IsWorkerThread()
Definition: G4Threading.cc:145
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
double G4double
Definition: G4Types.hh:76
G4int GetMinIndex(G4double e) const

Here is the call graph for this function:

G4double G4ParticleHPVector::GetXsec ( G4double  e,
G4int  min 
)
inline

Definition at line 157 of file G4ParticleHPVector.hh.

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  {
186  y = theInt.Interpolate(theManager.GetScheme(high), e,
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  }
int G4int
Definition: G4Types.hh:78
G4InterpolationScheme GetScheme(G4int index) const
G4double GetX(G4int i) const
G4double Interpolate(G4InterpolationScheme aScheme, G4double x, G4double x1, G4double x2, G4double y1, G4double y2) const
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

G4double G4ParticleHPVector::GetY ( G4double  x)
inline

Definition at line 198 of file G4ParticleHPVector.hh.

198 {return GetXsec(x);}
G4double GetXsec(G4int i)

Here is the call graph for this function:

Here is the caller graph for this function:

G4double G4ParticleHPVector::GetY ( G4int  i)
inline

Definition at line 201 of file G4ParticleHPVector.hh.

202  {
203  if (i<0) i=0;
204  if(i>=GetVectorLength()) i=GetVectorLength()-1;
205  return theData[i].GetY();
206  }
G4int GetVectorLength() const

Here is the call graph for this function:

G4double G4ParticleHPVector::GetY ( G4int  i) const
inline

Definition at line 208 of file G4ParticleHPVector.hh.

209  {
210  if (i<0) i=0;
211  if(i>=GetVectorLength()) i=GetVectorLength()-1;
212  return theData[i].GetY();
213  }
G4int GetVectorLength() const

Here is the call graph for this function:

void G4ParticleHPVector::Hash ( )
inline

Definition at line 135 of file G4ParticleHPVector.hh.

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  }
int G4int
Definition: G4Types.hh:78
void SetData(G4int index, G4double x, G4double y)
G4double GetX(G4int i) const
G4double GetY(G4double x)
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

Here is the caller graph for this function:

void G4ParticleHPVector::Init ( std::istream &  aDataFile,
G4int  total,
G4double  ux = 1.,
G4double  uy = 1. 
)
inline

Definition at line 221 of file G4ParticleHPVector.hh.

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  }
void SetData(G4int i, G4double x, G4double y)
int G4int
Definition: G4Types.hh:78
void SetData(G4int index, G4double x, G4double y)
G4double total(Particle const *const p1, Particle const *const p2)
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

Here is the caller graph for this function:

void G4ParticleHPVector::Init ( std::istream &  aDataFile,
G4double  ux = 1.,
G4double  uy = 1. 
)
inline

Definition at line 237 of file G4ParticleHPVector.hh.

238  {
239  G4int total;
240  aDataFile >> total;
241  if(theData!=0) delete [] theData;
242  theData = new G4ParticleHPDataPoint[total];
243  nPoints=total;
244  nEntries=0;
245  theManager.Init(aDataFile);
246  Init(aDataFile, total, ux, uy);
247  }
void Init(G4int aScheme, G4int aRange)
int G4int
Definition: G4Types.hh:78
void Init(std::istream &aDataFile, G4int total, G4double ux=1., G4double uy=1.)
G4double total(Particle const *const p1, Particle const *const p2)

Here is the call graph for this function:

void G4ParticleHPVector::InitInterpolation ( std::istream &  aDataFile)
inline

Definition at line 216 of file G4ParticleHPVector.hh.

217  {
218  theManager.Init(aDataFile);
219  }
void Init(G4int aScheme, G4int aRange)

Here is the call graph for this function:

void G4ParticleHPVector::Integrate ( )
inline

Definition at line 423 of file G4ParticleHPVector.hh.

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();
440  G4InterpolationScheme aScheme = theManager.GetScheme(i);
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  }
static G4Pow * GetInstance()
Definition: G4Pow.cc:55
G4double powA(G4double A, G4double y) const
Definition: G4Pow.hh:259
G4int GetVectorLength() const
int G4int
Definition: G4Types.hh:78
G4InterpolationScheme GetScheme(G4int index) const
G4double GetX(G4int i) const
G4double G4Log(G4double x)
Definition: G4Log.hh:230
G4double G4Exp(G4double initial_x)
Exponential Function double precision.
Definition: G4Exp.hh:183
G4InterpolationScheme
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

Here is the caller graph for this function:

void G4ParticleHPVector::IntegrateAndNormalise ( )
inline

Definition at line 373 of file G4ParticleHPVector.hh.

374  {
375  G4int i;
376  if(theIntegral!=0) return;
377  theIntegral = new G4double[nEntries];
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  //********************************************************************
400  G4InterpolationScheme aScheme = theManager.GetScheme(i);
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  }
416  G4double total = theIntegral[GetVectorLength()-1];
417  for(i=1;i<GetVectorLength();i++)
418  {
419  theIntegral[i]/=total;
420  }
421  }
G4int GetVectorLength() const
int G4int
Definition: G4Types.hh:78
G4InterpolationScheme GetScheme(G4int index) const
G4InterpolationScheme
G4double total(Particle const *const p1, Particle const *const p2)
G4double GetBinIntegral(const G4InterpolationScheme &aScheme, const G4double x1, const G4double x2, const G4double y1, const G4double y2)
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

Here is the caller graph for this function:

void G4ParticleHPVector::Merge ( G4ParticleHPVector active,
G4ParticleHPVector passive 
)
inline

Definition at line 273 of file G4ParticleHPVector.hh.

274  {
275  CleanUp();
276  G4int s_tmp = 0, n=0, m_tmp=0;
277  G4ParticleHPVector * tmp;
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  }
G4double GetEnergy(G4int i) const
G4int GetVectorLength() const
void SetData(G4int i, G4double x, G4double y)
const char * p
Definition: xmltok.h:285
int G4int
Definition: G4Types.hh:78
G4double GetXsec(G4int i)
void AppendScheme(G4int aPoint, const G4InterpolationScheme &aScheme)
double G4double
Definition: G4Types.hh:76
G4InterpolationScheme GetScheme(G4int anIndex)

Here is the call graph for this function:

Here is the caller graph for this function:

void G4ParticleHPVector::Merge ( G4InterpolationScheme  aScheme,
G4double  aValue,
G4ParticleHPVector active,
G4ParticleHPVector passive 
)

Definition at line 233 of file G4ParticleHPVector.cc.

235  {
236  // interpolate between labels according to aScheme, cut at aValue,
237  // continue in unknown areas by substraction of the last difference.
238 
239  CleanUp();
240  G4int s_tmp = 0, n=0, m_tmp=0;
241  G4ParticleHPVector * tmp;
242  G4int a = s_tmp, p = n, t;
243  while ( a<active->GetVectorLength() ) // Loop checking, 11.05.2015, T. Koi
244  {
245  if(active->GetEnergy(a) <= passive->GetEnergy(p))
246  {
247  G4double xa = active->GetEnergy(a);
248  G4double yy = theInt.Interpolate(aScheme, aValue, active->GetLabel(), passive->GetLabel(),
249  active->GetXsec(a), passive->GetXsec(xa));
250  SetData(m_tmp, xa, yy);
251  theManager.AppendScheme(m_tmp, active->GetScheme(a));
252  m_tmp++;
253  a++;
254  G4double xp = passive->GetEnergy(p);
255  //if( std::abs(std::abs(xp-xa)/xa)<0.0000001&&a<active->GetVectorLength() )
256  if ( xa != 0
257  && std::abs(std::abs(xp-xa)/xa) < 0.0000001
258  && a < active->GetVectorLength() )
259  {
260  p++;
261  tmp = active; t=a;
262  active = passive; a=p;
263  passive = tmp; p=t;
264  }
265  } else {
266  tmp = active; t=a;
267  active = passive; a=p;
268  passive = tmp; p=t;
269  }
270  }
271 
272  G4double deltaX = passive->GetXsec(GetEnergy(m_tmp-1)) - GetXsec(m_tmp-1);
273  while (p!=passive->GetVectorLength()&&passive->GetEnergy(p)<=aValue) // Loop checking, 11.05.2015, T. Koi
274  {
275  G4double anX;
276  anX = passive->GetXsec(p)-deltaX;
277  if(anX>0)
278  {
279  //if(std::abs(GetEnergy(m-1)-passive->GetEnergy(p))/passive->GetEnergy(p)>0.0000001)
280  if ( passive->GetEnergy(p) == 0
281  || std::abs(GetEnergy(m_tmp-1)-passive->GetEnergy(p))/passive->GetEnergy(p) > 0.0000001 )
282  {
283  SetData(m_tmp, passive->GetEnergy(p), anX);
284  theManager.AppendScheme(m_tmp++, passive->GetScheme(p));
285  }
286  }
287  p++;
288  }
289  // Rebuild the Hash;
290  if(theHash.Prepared())
291  {
292  ReHash();
293  }
294  }
G4double GetEnergy(G4int i) const
G4int GetVectorLength() const
G4bool Prepared() const
void SetData(G4int i, G4double x, G4double y)
const char * p
Definition: xmltok.h:285
int G4int
Definition: G4Types.hh:78
G4double GetXsec(G4int i)
void AppendScheme(G4int aPoint, const G4InterpolationScheme &aScheme)
G4double Interpolate(G4InterpolationScheme aScheme, G4double x, G4double x1, G4double x2, G4double y1, G4double y2) const
double G4double
Definition: G4Types.hh:76
G4InterpolationScheme GetScheme(G4int anIndex)

Here is the call graph for this function:

G4ParticleHPVector & G4ParticleHPVector::operator= ( const G4ParticleHPVector right)

Definition at line 125 of file G4ParticleHPVector.cc.

126  {
127  if(&right == this) return *this;
128 
129  G4int i;
130 
131  totalIntegral = right.totalIntegral;
132  if(right.theIntegral!=0) theIntegral = new G4double[right.nEntries];
133  for(i=0; i<right.nEntries; i++)
134  {
135  SetPoint(i, right.GetPoint(i)); // copy theData
136  if(right.theIntegral!=0) theIntegral[i] = right.theIntegral[i];
137  }
138  theManager = right.theManager;
139  label = right.label;
140 
141  Verbose = right.Verbose;
142  the15percentBorderCash = right.the15percentBorderCash;
143  the50percentBorderCash = right.the50percentBorderCash;
144  theHash = right.theHash;
145  return *this;
146  }
int G4int
Definition: G4Types.hh:78
void SetPoint(G4int i, const G4ParticleHPDataPoint &it)
double G4double
Definition: G4Types.hh:76
const G4ParticleHPDataPoint & GetPoint(G4int i) const

Here is the call graph for this function:

void G4ParticleHPVector::ReHash ( )
inline

Definition at line 150 of file G4ParticleHPVector.hh.

151  {
152  theHash.Clear();
153  Hash();
154  }

Here is the call graph for this function:

Here is the caller graph for this function:

G4double G4ParticleHPVector::Sample ( )

Definition at line 373 of file G4ParticleHPVector.cc.

374  {
376  G4int j;
377  for(j=0; j<GetVectorLength(); j++)
378  {
379  if(GetY(j)<0) SetY(j, 0);
380  }
381 
382  if(theBuffered.size() !=0 && G4UniformRand()<0.5)
383  {
384  result = theBuffered[0];
385  theBuffered.erase(theBuffered.begin());
386  if(result < GetX(GetVectorLength()-1) ) return result;
387  }
388  if(GetVectorLength()==1)
389  {
390  result = theData[0].GetX();
391  }
392  else
393  {
394  if(theIntegral==0) { IntegrateAndNormalise(); }
395  G4int icounter=0;
396  G4int icounter_max=1024;
397  do
398  {
399  icounter++;
400  if ( icounter > icounter_max ) {
401  G4cout << "Loop-counter exceeded the threshold value at " << __LINE__ << "th line of " << __FILE__ << "." << G4endl;
402  break;
403  }
404 //080808
405 /*
406  G4double rand;
407  G4double value, test, baseline;
408  baseline = theData[GetVectorLength()-1].GetX()-theData[0].GetX();
409  do
410  {
411  value = baseline*G4UniformRand();
412  value += theData[0].GetX();
413  test = GetY(value)/maxValue;
414  rand = G4UniformRand();
415  }
416  //while(test<rand);
417  while( test < rand && test > 0 );
418  result = value;
419 */
420  G4double rand;
421  G4double value, test;
422  G4int jcounter=0;
423  G4int jcounter_max=1024;
424  do
425  {
426  jcounter++;
427  if ( jcounter > jcounter_max ) {
428  G4cout << "Loop-counter exceeded the threshold value at " << __LINE__ << "th line of " << __FILE__ << "." << G4endl;
429  break;
430  }
431  rand = G4UniformRand();
432  G4int ibin = -1;
433  for ( G4int i = 0 ; i < GetVectorLength() ; i++ )
434  {
435  if ( rand < theIntegral[i] )
436  {
437  ibin = i;
438  break;
439  }
440  }
441  if ( ibin < 0 ) G4cout << "TKDB 080807 " << rand << G4endl;
442  // result
443  rand = G4UniformRand();
444  G4double x1, x2;
445  if ( ibin == 0 )
446  {
447  x1 = theData[ ibin ].GetX();
448  value = x1;
449  break;
450  }
451  else
452  {
453  x1 = theData[ ibin-1 ].GetX();
454  }
455 
456  x2 = theData[ ibin ].GetX();
457  value = rand * ( x2 - x1 ) + x1;
458  //***********************************************************************
459  /*
460  test = GetY ( value ) / std::max ( GetY( ibin-1 ) , GetY ( ibin ) );
461  */
462  //***********************************************************************
463  //EMendoza - Always linear interpolation:
464  G4double y1=theData[ ibin-1 ].GetY();
465  G4double y2=theData[ ibin ].GetY();
466  G4double mval=(y2-y1)/(x2-x1);
467  G4double bval=y1-mval*x1;
468  test =(mval*value+bval)/std::max ( GetY( ibin-1 ) , GetY ( ibin ) );
469  //***********************************************************************
470  }
471  while ( G4UniformRand() > test ); // Loop checking, 11.05.2015, T. Koi
472  result = value;
473 //080807
474  }
475  while(IsBlocked(result)); // Loop checking, 11.05.2015, T. Koi
476  }
477  return result;
478  }
G4double G4ParticleHPJENDLHEData::G4double result
G4int GetVectorLength() const
int G4int
Definition: G4Types.hh:78
void SetY(G4int i, G4double x)
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
const XML_Char int const XML_Char * value
Definition: expat.h:331
G4double GetX(G4int i) const
G4double GetY(G4double x)
T max(const T t1, const T t2)
brief Return the largest of the two arguments
#define G4endl
Definition: G4ios.hh:61
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

Here is the caller graph for this function:

G4double G4ParticleHPVector::SampleLin ( )
inline

Definition at line 324 of file G4ParticleHPVector.hh.

325  {
327  if(theIntegral==0) IntegrateAndNormalise();
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  }
G4double G4ParticleHPJENDLHEData::G4double result
G4int GetVectorLength() const
G4double Lin(G4double x, G4double x1, G4double x2, G4double y1, G4double y2)
int G4int
Definition: G4Types.hh:78
#define G4UniformRand()
Definition: Randomize.hh:97
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

void G4ParticleHPVector::SetData ( G4int  i,
G4double  x,
G4double  y 
)
inline

Definition at line 96 of file G4ParticleHPVector.hh.

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  }
void SetData(G4double e, G4double x)

Here is the call graph for this function:

Here is the caller graph for this function:

void G4ParticleHPVector::SetEnergy ( G4int  i,
G4double  e 
)
inline

Definition at line 108 of file G4ParticleHPVector.hh.

109  {
110  Check(i);
111  theData[i].SetX(e);
112  }

Here is the call graph for this function:

void G4ParticleHPVector::SetInterpolationManager ( const G4InterpolationManager aManager)
inline

Definition at line 483 of file G4ParticleHPVector.hh.

484  {
485  theManager = aManager;
486  }

Here is the caller graph for this function:

void G4ParticleHPVector::SetInterpolationManager ( G4InterpolationManager aMan)
inline

Definition at line 493 of file G4ParticleHPVector.hh.

494  {
495  theManager = aMan;
496  }
void G4ParticleHPVector::SetLabel ( G4double  aLabel)
inline

Definition at line 251 of file G4ParticleHPVector.hh.

252  {
253  label = aLabel;
254  }

Here is the caller graph for this function:

void G4ParticleHPVector::SetPoint ( G4int  i,
const G4ParticleHPDataPoint it 
)
inline

Definition at line 89 of file G4ParticleHPVector.hh.

90  {
91  G4double x = it.GetX();
92  G4double y = it.GetY();
93  SetData(i, x, y);
94  }
void SetData(G4int i, G4double x, G4double y)
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

Here is the caller graph for this function:

void G4ParticleHPVector::SetScheme ( G4int  aPoint,
const G4InterpolationScheme aScheme 
)
inline

Definition at line 498 of file G4ParticleHPVector.hh.

499  {
500  theManager.AppendScheme(aPoint, aScheme);
501  }
void AppendScheme(G4int aPoint, const G4InterpolationScheme &aScheme)

Here is the call graph for this function:

Here is the caller graph for this function:

void G4ParticleHPVector::SetVerbose ( G4int  ff)
inline

Definition at line 71 of file G4ParticleHPVector.hh.

72  {
73  Verbose = ff;
74  }
void G4ParticleHPVector::SetX ( G4int  i,
G4double  e 
)
inline

Definition at line 103 of file G4ParticleHPVector.hh.

104  {
105  Check(i);
106  theData[i].SetX(e);
107  }

Here is the call graph for this function:

Here is the caller graph for this function:

void G4ParticleHPVector::SetXsec ( G4int  i,
G4double  x 
)
inline

Definition at line 119 of file G4ParticleHPVector.hh.

120  {
121  Check(i);
122  if(x>maxValue) maxValue=x;
123  theData[i].SetY(x);
124  }

Here is the call graph for this function:

void G4ParticleHPVector::SetY ( G4int  i,
G4double  x 
)
inline

Definition at line 113 of file G4ParticleHPVector.hh.

114  {
115  Check(i);
116  if(x>maxValue) maxValue=x;
117  theData[i].SetY(x);
118  }

Here is the call graph for this function:

Here is the caller graph for this function:

void G4ParticleHPVector::ThinOut ( G4double  precision)

Definition at line 296 of file G4ParticleHPVector.cc.

297  {
298  // anything in there?
299  if(GetVectorLength()==0) return;
300  // make the new vector
301  G4ParticleHPDataPoint * aBuff = new G4ParticleHPDataPoint[nPoints];
302  G4double x, x1, x2, y, y1, y2;
303  G4int count = 0, current = 2, start = 1;
304 
305  // First element always goes and is never tested.
306  aBuff[0] = theData[0];
307 
308  // Find the rest
309  while(current < GetVectorLength()) // Loop checking, 11.05.2015, T. Koi
310  {
311  x1=aBuff[count].GetX();
312  y1=aBuff[count].GetY();
313  x2=theData[current].GetX();
314  y2=theData[current].GetY();
315 
316  if ( x1-x2 == 0 ) {
317  //Following block added for avoiding div 0 error on Release + G4FPE_DEBUG
318  for ( G4int j=start; j<current; j++ ) {
319  x = theData[j].GetX();
320  y = (y2+y1)/2.;
321  if ( std::abs( y-theData[j].GetY() ) > precision*y ) {
322  aBuff[++count] = theData[current-1]; // for this one, everything was fine
323  start = current; // the next candidate
324  break;
325  }
326  }
327  } else {
328  for(G4int j=start; j<current; j++)
329  {
330  x = theData[j].GetX();
331  if(x1-x2 == 0) y = (y2+y1)/2.;
332  else y = theInt.Lin(x, x1, x2, y1, y2);
333  if (std::abs(y-theData[j].GetY())>precision*y)
334  {
335  aBuff[++count] = theData[current-1]; // for this one, everything was fine
336  start = current; // the next candidate
337  break;
338  }
339  }
340  }
341  current++ ;
342  }
343  // The last one also always goes, and is never tested.
344  aBuff[++count] = theData[GetVectorLength()-1];
345  delete [] theData;
346  theData = aBuff;
347  nEntries = count+1;
348 
349  // Rebuild the Hash;
350  if(theHash.Prepared())
351  {
352  ReHash();
353  }
354  }
G4int GetVectorLength() const
G4bool Prepared() const
G4double Lin(G4double x, G4double x1, G4double x2, G4double y1, G4double y2)
int G4int
Definition: G4Types.hh:78
G4double GetY(G4double x)
double G4double
Definition: G4Types.hh:76

Here is the call graph for this function:

Here is the caller graph for this function:

void G4ParticleHPVector::Times ( G4double  factor)
inline

Definition at line 76 of file G4ParticleHPVector.hh.

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  }
int G4int
Definition: G4Types.hh:78
G4double GetY(G4double x)

Here is the call graph for this function:

Here is the caller graph for this function:

Friends And Related Function Documentation

G4ParticleHPVector& operator+ ( G4ParticleHPVector left,
G4ParticleHPVector right 
)
friend

Definition at line 40 of file G4ParticleHPVector.cc.

41  {
43  G4int j=0;
44  G4double x;
45  G4double y;
46  G4int running = 0;
47  for(G4int i=0; i<left.GetVectorLength(); i++)
48  {
49  while(j<right.GetVectorLength()) // Loop checking, 11.05.2015, T. Koi
50  {
51  if(right.GetX(j)<left.GetX(i)*1.001)
52  {
53  x = right.GetX(j);
54  y = right.GetY(j)+left.GetY(x);
55  result->SetData(running++, x, y);
56  j++;
57  }
58  //else if(std::abs((right.GetX(j)-left.GetX(i))/(left.GetX(i)+right.GetX(j)))>0.001)
59  else if( left.GetX(i)+right.GetX(j) == 0
60  || std::abs((right.GetX(j)-left.GetX(i))/(left.GetX(i)+right.GetX(j))) > 0.001 )
61  {
62  x = left.GetX(i);
63  y = left.GetY(i)+right.GetY(x);
64  result->SetData(running++, x, y);
65  break;
66  }
67  else
68  {
69  break;
70  }
71  }
72  if(j==right.GetVectorLength())
73  {
74  x = left.GetX(i);
75  y = left.GetY(i)+right.GetY(x);
76  result->SetData(running++, x, y);
77  }
78  }
79  result->ThinOut(0.02);
80  return *result;
81  }
G4double G4ParticleHPJENDLHEData::G4double result
G4int GetVectorLength() const
void SetData(G4int i, G4double x, G4double y)
int G4int
Definition: G4Types.hh:78
G4double GetX(G4int i) const
G4double GetY(G4double x)
double G4double
Definition: G4Types.hh:76
void ThinOut(G4double precision)

The documentation for this class was generated from the following files: