Geant4  10.03.p02
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
G4ParticleHPVector.cc
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 // neutron_hp -- source file
27 // J.P. Wellisch, Nov-1996
28 // A prototype of the low energy neutron transport model.
29 //
30 // 070523 bug fix for G4FPE_DEBUG on by A. Howard ( and T. Koi)
31 // 080808 bug fix in Sample() and GetXsec() by T. Koi
32 //
33 // P. Arce, June-2014 Conversion neutron_hp to particle_hp
34 //
35 #include "G4ParticleHPVector.hh"
36 #include "G4SystemOfUnits.hh"
37 #include "G4Threading.hh"
38 
39  // if the ranges do not match, constant extrapolation is used.
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  }
82 
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  }
97 
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  }
112 
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  }
123 
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  }
147 
148 
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  }
204 
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  }
217 
218  void G4ParticleHPVector::Check(G4int i)
219  {
220  if(i>nEntries) throw G4HadronicException(__FILE__, __LINE__, "Skipped some index numbers in G4ParticleHPVector");
221  if(i==nPoints)
222  {
223  nPoints = static_cast<G4int>(1.2*nPoints);
224  G4ParticleHPDataPoint * buff = new G4ParticleHPDataPoint[nPoints];
225  for (G4int j=0; j<nEntries; j++) buff[j] = theData[j];
226  delete [] theData;
227  theData = buff;
228  }
229  if(i==nEntries) nEntries=i+1;
230  }
231 
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  }
295 
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  y = (y2+y1)/2.;
320  if ( std::abs( y-theData[j].GetY() ) > precision*y ) {
321  aBuff[++count] = theData[current-1]; // for this one, everything was fine
322  start = current; // the next candidate
323  break;
324  }
325  }
326  } else {
327  for(G4int j=start; j<current; j++)
328  {
329  x = theData[j].GetX();
330  if(x1-x2 == 0) y = (y2+y1)/2.;
331  else y = theInt.Lin(x, x1, x2, y1, y2);
332  if (std::abs(y-theData[j].GetY())>precision*y)
333  {
334  aBuff[++count] = theData[current-1]; // for this one, everything was fine
335  start = current; // the next candidate
336  break;
337  }
338  }
339  }
340  current++ ;
341  }
342  // The last one also always goes, and is never tested.
343  aBuff[++count] = theData[GetVectorLength()-1];
344  delete [] theData;
345  theData = aBuff;
346  nEntries = count+1;
347 
348  // Rebuild the Hash;
349  if(theHash.Prepared())
350  {
351  ReHash();
352  }
353  }
354 
355  G4bool G4ParticleHPVector::IsBlocked(G4double aX)
356  {
357  G4bool result = false;
358  std::vector<G4double>::iterator i;
359  for(i=theBlocked.begin(); i!=theBlocked.end(); i++)
360  {
361  G4double aBlock = *i;
362  if(std::abs(aX-aBlock) < 0.1*MeV)
363  {
364  result = true;
365  theBlocked.erase(i);
366  break;
367  }
368  }
369  return result;
370  }
371 
372  G4double G4ParticleHPVector::Sample() // Samples X according to distribution Y
373  {
375  G4int j;
376  for(j=0; j<GetVectorLength(); j++)
377  {
378  if(GetY(j)<0) SetY(j, 0);
379  }
380 
381  if(theBuffered.size() !=0 && G4UniformRand()<0.5)
382  {
383  result = theBuffered[0];
384  theBuffered.erase(theBuffered.begin());
385  if(result < GetX(GetVectorLength()-1) ) return result;
386  }
387  if(GetVectorLength()==1)
388  {
389  result = theData[0].GetX();
390  }
391  else
392  {
393  if(theIntegral==0) { IntegrateAndNormalise(); }
394  G4int icounter=0;
395  G4int icounter_max=1024;
396  do
397  {
398  icounter++;
399  if ( icounter > icounter_max ) {
400  G4cout << "Loop-counter exceeded the threshold value at " << __LINE__ << "th line of " << __FILE__ << "." << G4endl;
401  break;
402  }
403 //080808
404 /*
405  G4double rand;
406  G4double value, test, baseline;
407  baseline = theData[GetVectorLength()-1].GetX()-theData[0].GetX();
408  do
409  {
410  value = baseline*G4UniformRand();
411  value += theData[0].GetX();
412  test = GetY(value)/maxValue;
413  rand = G4UniformRand();
414  }
415  //while(test<rand);
416  while( test < rand && test > 0 );
417  result = value;
418 */
419  G4double rand;
421  G4int jcounter=0;
422  G4int jcounter_max=1024;
423  do
424  {
425  jcounter++;
426  if ( jcounter > jcounter_max ) {
427  G4cout << "Loop-counter exceeded the threshold value at " << __LINE__ << "th line of " << __FILE__ << "." << G4endl;
428  break;
429  }
430  rand = G4UniformRand();
431  G4int ibin = -1;
432  for ( G4int i = 0 ; i < GetVectorLength() ; i++ )
433  {
434  if ( rand < theIntegral[i] )
435  {
436  ibin = i;
437  break;
438  }
439  }
440  if ( ibin < 0 ) G4cout << "TKDB 080807 " << rand << G4endl;
441  // result
442  rand = G4UniformRand();
443  G4double x1, x2;
444  if ( ibin == 0 )
445  {
446  x1 = theData[ ibin ].GetX();
447  value = x1;
448  break;
449  }
450  else
451  {
452  x1 = theData[ ibin-1 ].GetX();
453  }
454 
455  x2 = theData[ ibin ].GetX();
456  value = rand * ( x2 - x1 ) + x1;
457  //***********************************************************************
458  /*
459  test = GetY ( value ) / std::max ( GetY( ibin-1 ) , GetY ( ibin ) );
460  */
461  //***********************************************************************
462  //EMendoza - Always linear interpolation:
463  G4double y1=theData[ ibin-1 ].GetY();
464  G4double y2=theData[ ibin ].GetY();
465  G4double mval=(y2-y1)/(x2-x1);
466  G4double bval=y1-mval*x1;
467  test =(mval*value+bval)/std::max ( GetY( ibin-1 ) , GetY ( ibin ) );
468  //***********************************************************************
469  }
470  while ( G4UniformRand() > test ); // Loop checking, 11.05.2015, T. Koi
471  result = value;
472 //080807
473  }
474  while(IsBlocked(result)); // Loop checking, 11.05.2015, T. Koi
475  }
476  return result;
477  }
478 
480  {
481  if(the15percentBorderCash>-DBL_MAX/2.) return the15percentBorderCash;
483  if(GetVectorLength()==1)
484  {
485  result = theData[0].GetX();
486  the15percentBorderCash = result;
487  }
488  else
489  {
490  if(theIntegral==0) { IntegrateAndNormalise(); }
491  G4int i;
492  result = theData[GetVectorLength()-1].GetX();
493  for(i=0;i<GetVectorLength();i++)
494  {
495  if(theIntegral[i]/theIntegral[GetVectorLength()-1]>0.15)
496  {
497  result = theData[std::min(i+1, GetVectorLength()-1)].GetX();
498  the15percentBorderCash = result;
499  break;
500  }
501  }
502  the15percentBorderCash = result;
503  }
504  return result;
505  }
506 
508  {
509  if(the50percentBorderCash>-DBL_MAX/2.) return the50percentBorderCash;
511  if(GetVectorLength()==1)
512  {
513  result = theData[0].GetX();
514  the50percentBorderCash = result;
515  }
516  else
517  {
518  if(theIntegral==0) { IntegrateAndNormalise(); }
519  G4int i;
520  G4double x = 0.5;
521  result = theData[GetVectorLength()-1].GetX();
522  for(i=0;i<GetVectorLength();i++)
523  {
524  if(theIntegral[i]/theIntegral[GetVectorLength()-1]>x)
525  {
526  G4int it;
527  it = i;
528  if(it == GetVectorLength()-1)
529  {
530  result = theData[GetVectorLength()-1].GetX();
531  }
532  else
533  {
534  G4double x1, x2, y1, y2;
535  x1 = theIntegral[i-1]/theIntegral[GetVectorLength()-1];
536  x2 = theIntegral[i]/theIntegral[GetVectorLength()-1];
537  y1 = theData[i-1].GetX();
538  y2 = theData[i].GetX();
539  result = theLin.Lin(x, x1, x2, y1, y2);
540  }
541  the50percentBorderCash = result;
542  break;
543  }
544  }
545  the50percentBorderCash = result;
546  }
547  return result;
548  }
G4double G4ParticleHPJENDLHEData::G4double result
G4double GetEnergy(G4int i) const
G4int GetVectorLength() const
G4bool Prepared() const
G4double Lin(G4double x, G4double x1, G4double x2, G4double y1, G4double y2)
std::vector< ExP01TrackerHit * > a
Definition: ExP01Classes.hh:33
void SetData(G4int i, G4double x, G4double y)
const char * p
Definition: xmltok.h:285
BasicVector3D< float > operator+(const BasicVector3D< float > &v)
tuple x
Definition: test.py:50
G4ParticleHPVector & operator=(const G4ParticleHPVector &right)
int G4int
Definition: G4Types.hh:78
G4double GetXsec(G4int i)
void SetY(G4int i, G4double x)
void AppendScheme(G4int aPoint, const G4InterpolationScheme &aScheme)
#define G4UniformRand()
Definition: Randomize.hh:97
G4GLOB_DLL std::ostream G4cout
const XML_Char int const XML_Char * value
Definition: expat.h:331
bool G4bool
Definition: G4Types.hh:79
G4InterpolationScheme GetScheme(G4int index) const
void Merge(G4ParticleHPVector *active, G4ParticleHPVector *passive)
G4double GetX(G4int i) const
void SetPoint(G4int i, const G4ParticleHPDataPoint &it)
G4double GetY(G4double x)
const G4int n
G4double Interpolate(G4InterpolationScheme aScheme, G4double x, G4double x1, G4double x2, G4double y1, G4double y2) const
G4bool IsWorkerThread()
Definition: G4Threading.cc:145
G4InterpolationScheme
T max(const T t1, const T t2)
brief Return the largest of the two arguments
T min(const T t1, const T t2)
brief Return the smallest of the two arguments
#define G4endl
Definition: G4ios.hh:61
static constexpr double MeV
Definition: G4SIunits.hh:214
double G4double
Definition: G4Types.hh:76
G4int GetMinIndex(G4double e) const
G4InterpolationScheme GetScheme(G4int anIndex)
const G4ParticleHPDataPoint & GetPoint(G4int i) const
program test
Definition: Main_HIJING.f:1
#define DBL_MAX
Definition: templates.hh:83
void ThinOut(G4double precision)