Geant4 9.6.0
Toolkit for the simulation of the passage of particles through matter
Loading...
Searching...
No Matches
G4NeutronHPVector Class Reference

#include <G4NeutronHPVector.hh>

Public Member Functions

 G4NeutronHPVector ()
 
 G4NeutronHPVector (G4int n)
 
 ~G4NeutronHPVector ()
 
G4NeutronHPVectoroperator= (const G4NeutronHPVector &right)
 
void SetVerbose (G4int ff)
 
void Times (G4double factor)
 
void SetPoint (G4int i, const G4NeutronHPDataPoint &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 G4NeutronHPDataPointGetPoint (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::ifstream &aDataFile)
 
void Init (std::ifstream &aDataFile, G4int total, G4double ux=1., G4double uy=1.)
 
void Init (std::ifstream &aDataFile, G4double ux=1., G4double uy=1.)
 
void ThinOut (G4double precision)
 
void SetLabel (G4double aLabel)
 
G4double GetLabel ()
 
void CleanUp ()
 
void Merge (G4NeutronHPVector *active, G4NeutronHPVector *passive)
 
void Merge (G4InterpolationScheme aScheme, G4double aValue, G4NeutronHPVector *active, G4NeutronHPVector *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

G4NeutronHPVectoroperator+ (G4NeutronHPVector &left, G4NeutronHPVector &right)
 

Detailed Description

Definition at line 50 of file G4NeutronHPVector.hh.

Constructor & Destructor Documentation

◆ G4NeutronHPVector() [1/2]

G4NeutronHPVector::G4NeutronHPVector ( )

Definition at line 80 of file G4NeutronHPVector.cc.

81 {
82 theData = new G4NeutronHPDataPoint[20];
83 nPoints=20;
84 nEntries=0;
85 Verbose=0;
86 theIntegral=0;
87 totalIntegral=-1;
88 isFreed = 0;
89 maxValue = -DBL_MAX;
90 the15percentBorderCash = -DBL_MAX;
91 the50percentBorderCash = -DBL_MAX;
92 label = -DBL_MAX;
93
94 }
#define DBL_MAX
Definition: templates.hh:83

◆ G4NeutronHPVector() [2/2]

G4NeutronHPVector::G4NeutronHPVector ( G4int  n)

Definition at line 96 of file G4NeutronHPVector.cc.

97 {
98 nPoints=std::max(n, 20);
99 theData = new G4NeutronHPDataPoint[nPoints];
100 nEntries=0;
101 Verbose=0;
102 theIntegral=0;
103 totalIntegral=-1;
104 isFreed = 0;
105 maxValue = -DBL_MAX;
106 the15percentBorderCash = -DBL_MAX;
107 the50percentBorderCash = -DBL_MAX;
108 }

◆ ~G4NeutronHPVector()

G4NeutronHPVector::~G4NeutronHPVector ( )

Definition at line 110 of file G4NeutronHPVector.cc.

111 {
112// if(Verbose==1)G4cout <<"G4NeutronHPVector::~G4NeutronHPVector"<<G4endl;
113 delete [] theData;
114// if(Verbose==1)G4cout <<"Vector: delete theData"<<G4endl;
115 delete [] theIntegral;
116// if(Verbose==1)G4cout <<"Vector: delete theIntegral"<<G4endl;
117 isFreed = 1;
118 }

Member Function Documentation

◆ CleanUp()

void G4NeutronHPVector::CleanUp ( )
inline

Definition at line 255 of file G4NeutronHPVector.hh.

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 }

Referenced by Merge().

◆ Debug()

G4double * G4NeutronHPVector::Debug ( )
inline

Definition at line 362 of file G4NeutronHPVector.hh.

363 {
364 return theIntegral;
365 }

◆ Dump()

void G4NeutronHPVector::Dump ( )

Definition at line 194 of file G4NeutronHPVector.cc.

195 {
196 G4cout << nEntries<<G4endl;
197 for(G4int i=0; i<nEntries; i++)
198 {
199 G4cout << theData[i].GetX()<<" ";
200 G4cout << theData[i].GetY()<<" ";
201// if (i!=1&&i==5*(i/5)) G4cout << G4endl;
202 G4cout << G4endl;
203 }
204 G4cout << G4endl;
205 }
int G4int
Definition: G4Types.hh:66
#define G4endl
Definition: G4ios.hh:52
G4DLLIMPORT std::ostream G4cout

◆ Get15percentBorder()

G4double G4NeutronHPVector::Get15percentBorder ( )

Definition at line 441 of file G4NeutronHPVector.cc.

442 {
443 if(the15percentBorderCash>-DBL_MAX/2.) return the15percentBorderCash;
444 G4double result;
445 if(GetVectorLength()==1)
446 {
447 result = theData[0].GetX();
448 the15percentBorderCash = result;
449 }
450 else
451 {
452 if(theIntegral==0) { IntegrateAndNormalise(); }
453 G4int i;
454 result = theData[GetVectorLength()-1].GetX();
455 for(i=0;i<GetVectorLength();i++)
456 {
457 if(theIntegral[i]/theIntegral[GetVectorLength()-1]>0.15)
458 {
459 result = theData[std::min(i+1, GetVectorLength()-1)].GetX();
460 the15percentBorderCash = result;
461 break;
462 }
463 }
464 the15percentBorderCash = result;
465 }
466 return result;
467 }
double G4double
Definition: G4Types.hh:64
G4int GetVectorLength() const

◆ Get50percentBorder()

G4double G4NeutronHPVector::Get50percentBorder ( )

Definition at line 469 of file G4NeutronHPVector.cc.

470 {
471 if(the50percentBorderCash>-DBL_MAX/2.) return the50percentBorderCash;
472 G4double result;
473 if(GetVectorLength()==1)
474 {
475 result = theData[0].GetX();
476 the50percentBorderCash = result;
477 }
478 else
479 {
480 if(theIntegral==0) { IntegrateAndNormalise(); }
481 G4int i;
482 G4double x = 0.5;
483 result = theData[GetVectorLength()-1].GetX();
484 for(i=0;i<GetVectorLength();i++)
485 {
486 if(theIntegral[i]/theIntegral[GetVectorLength()-1]>x)
487 {
488 G4int it;
489 it = i;
490 if(it == GetVectorLength()-1)
491 {
492 result = theData[GetVectorLength()-1].GetX();
493 }
494 else
495 {
496 G4double x1, x2, y1, y2;
497 x1 = theIntegral[i-1]/theIntegral[GetVectorLength()-1];
498 x2 = theIntegral[i]/theIntegral[GetVectorLength()-1];
499 y1 = theData[i-1].GetX();
500 y2 = theData[i].GetX();
501 result = theLin.Lin(x, x1, x2, y1, y2);
502 }
503 the50percentBorderCash = result;
504 break;
505 }
506 }
507 the50percentBorderCash = result;
508 }
509 return result;
510 }
G4double Lin(G4double x, G4double x1, G4double x2, G4double y1, G4double y2)

◆ GetBlocked()

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

Definition at line 532 of file G4NeutronHPVector.hh.

532{return theBlocked;}

◆ GetBuffered()

std::vector< G4double > G4NeutronHPVector::GetBuffered ( )
inline

Definition at line 533 of file G4NeutronHPVector.hh.

533{return theBuffered;}

◆ GetEnergy()

G4double G4NeutronHPVector::GetEnergy ( G4int  i) const
inline

◆ GetIntegral()

G4double G4NeutronHPVector::GetIntegral ( )
inline

Definition at line 471 of file G4NeutronHPVector.hh.

472 {
473 if(totalIntegral<-0.5) Integrate();
474 return totalIntegral;
475 }

Referenced by G4NeutronHPLabAngularEnergy::Sample().

◆ GetInterpolationManager()

const G4InterpolationManager & G4NeutronHPVector::GetInterpolationManager ( ) const
inline

Definition at line 482 of file G4NeutronHPVector.hh.

483 {
484 return theManager;
485 }

Referenced by G4NeutronHPLabAngularEnergy::Sample().

◆ GetLabel()

G4double G4NeutronHPVector::GetLabel ( )
inline

Definition at line 250 of file G4NeutronHPVector.hh.

251 {
252 return label;
253 }

Referenced by Merge(), G4NeutronHPArbitaryTab::Sample(), and G4NeutronHPLabAngularEnergy::Sample().

◆ GetMeanX()

G4double G4NeutronHPVector::GetMeanX ( )
inline

Definition at line 502 of file G4NeutronHPVector.hh.

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 }
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)

Referenced by G4NeutronHPLabAngularEnergy::Sample().

◆ GetPoint()

const G4NeutronHPDataPoint & G4NeutronHPVector::GetPoint ( G4int  i) const
inline

Definition at line 127 of file G4NeutronHPVector.hh.

127{ return theData[i]; }

Referenced by G4NeutronHPIsoData::FillChannelData(), and operator=().

◆ GetScheme()

G4InterpolationScheme G4NeutronHPVector::GetScheme ( G4int  anIndex)
inline

Definition at line 497 of file G4NeutronHPVector.hh.

498 {
499 return theManager.GetScheme(anIndex);
500 }

Referenced by Merge().

◆ GetVectorLength()

◆ GetX()

◆ GetXsec() [1/3]

G4double G4NeutronHPVector::GetXsec ( G4double  e)

Definition at line 145 of file G4NeutronHPVector.cc.

146 {
147 if(nEntries == 0) return 0;
148 if(!theHash.Prepared()) Hash();
149 G4int min = theHash.GetMinIndex(e);
150 G4int i;
151 for(i=min ; i<nEntries; i++)
152 {
153 //if(theData[i].GetX()>e) break;
154 if(theData[i].GetX() >= e) break;
155 }
156 G4int low = i-1;
157 G4int high = i;
158 if(i==0)
159 {
160 low = 0;
161 high = 1;
162 }
163 else if(i==nEntries)
164 {
165 low = nEntries-2;
166 high = nEntries-1;
167 }
168 G4double y;
169 if(e<theData[nEntries-1].GetX())
170 {
171 // Protect against doubled-up x values
172 //if( (theData[high].GetX()-theData[low].GetX())/theData[high].GetX() < 0.000001)
173 if ( theData[high].GetX() !=0
174 //080808 TKDB
175 //&&( theData[high].GetX()-theData[low].GetX())/theData[high].GetX() < 0.000001)
176 &&( std::abs( (theData[high].GetX()-theData[low].GetX())/theData[high].GetX() ) < 0.000001 ) )
177 {
178 y = theData[low].GetY();
179 }
180 else
181 {
182 y = theInt.Interpolate(theManager.GetScheme(high), e,
183 theData[low].GetX(), theData[high].GetX(),
184 theData[low].GetY(), theData[high].GetY());
185 }
186 }
187 else
188 {
189 y=theData[nEntries-1].GetY();
190 }
191 return y;
192 }
G4int GetMinIndex(G4double e) const
G4bool Prepared() const
G4double Interpolate(G4InterpolationScheme aScheme, G4double x, G4double x1, G4double x2, G4double y1, G4double y2) const
G4double GetX(G4int i) const

◆ GetXsec() [2/3]

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

Definition at line 151 of file G4NeutronHPVector.hh.

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 }

◆ GetXsec() [3/3]

◆ GetY() [1/3]

G4double G4NeutronHPVector::GetY ( G4double  x)
inline

◆ GetY() [2/3]

G4double G4NeutronHPVector::GetY ( G4int  i)
inline

Definition at line 195 of file G4NeutronHPVector.hh.

196 {
197 if (i<0) i=0;
198 if(i>=GetVectorLength()) i=GetVectorLength()-1;
199 return theData[i].GetY();
200 }

◆ GetY() [3/3]

G4double G4NeutronHPVector::GetY ( G4int  i) const
inline

Definition at line 202 of file G4NeutronHPVector.hh.

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

◆ Hash()

void G4NeutronHPVector::Hash ( )
inline

Definition at line 129 of file G4NeutronHPVector.hh.

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 }
void SetData(G4int index, G4double x, G4double y)
G4double GetY(G4double x)

Referenced by GetXsec(), and ReHash().

◆ Init() [1/2]

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

Definition at line 231 of file G4NeutronHPVector.hh.

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 }
void Init(G4int aScheme, G4int aRange)
void Init(std::ifstream &aDataFile, G4int total, G4double ux=1., G4double uy=1.)

◆ Init() [2/2]

◆ InitInterpolation()

void G4NeutronHPVector::InitInterpolation ( std::ifstream &  aDataFile)
inline

Definition at line 210 of file G4NeutronHPVector.hh.

211 {
212 theManager.Init(aDataFile);
213 }

Referenced by G4NeutronHPPartial::InitData(), and G4NeutronHPPartial::InitInterpolation().

◆ Integrate()

void G4NeutronHPVector::Integrate ( )
inline

Definition at line 417 of file G4NeutronHPVector.hh.

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 }
G4InterpolationScheme

Referenced by GetIntegral().

◆ IntegrateAndNormalise()

void G4NeutronHPVector::IntegrateAndNormalise ( )
inline

Definition at line 367 of file G4NeutronHPVector.hh.

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 }

Referenced by Get15percentBorder(), Get50percentBorder(), Sample(), and SampleLin().

◆ Merge() [1/2]

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

Definition at line 221 of file G4NeutronHPVector.cc.

224 {
225 // interpolate between labels according to aScheme, cut at aValue,
226 // continue in unknown areas by substraction of the last difference.
227
228 CleanUp();
229 G4int s_tmp = 0, n=0, m_tmp=0;
230 G4NeutronHPVector * tmp;
231 G4int a = s_tmp, p = n, t;
232 while ( a<active->GetVectorLength() )
233 {
234 if(active->GetEnergy(a) <= passive->GetEnergy(p))
235 {
236 G4double xa = active->GetEnergy(a);
237 G4double yy = theInt.Interpolate(aScheme, aValue, active->GetLabel(), passive->GetLabel(),
238 active->GetXsec(a), passive->GetXsec(xa));
239 SetData(m_tmp, xa, yy);
240 theManager.AppendScheme(m_tmp, active->GetScheme(a));
241 m_tmp++;
242 a++;
243 G4double xp = passive->GetEnergy(p);
244 //if( std::abs(std::abs(xp-xa)/xa)<0.0000001&&a<active->GetVectorLength() )
245 if ( xa != 0
246 && std::abs(std::abs(xp-xa)/xa) < 0.0000001
247 && a < active->GetVectorLength() )
248 {
249 p++;
250 tmp = active; t=a;
251 active = passive; a=p;
252 passive = tmp; p=t;
253 }
254 } else {
255 tmp = active; t=a;
256 active = passive; a=p;
257 passive = tmp; p=t;
258 }
259 }
260
261 G4double deltaX = passive->GetXsec(GetEnergy(m_tmp-1)) - GetXsec(m_tmp-1);
262 while (p!=passive->GetVectorLength()&&passive->GetEnergy(p)<=aValue)
263 {
264 G4double anX;
265 anX = passive->GetXsec(p)-deltaX;
266 if(anX>0)
267 {
268 //if(std::abs(GetEnergy(m-1)-passive->GetEnergy(p))/passive->GetEnergy(p)>0.0000001)
269 if ( passive->GetEnergy(p) == 0
270 || std::abs(GetEnergy(m_tmp-1)-passive->GetEnergy(p))/passive->GetEnergy(p) > 0.0000001 )
271 {
272 SetData(m_tmp, passive->GetEnergy(p), anX);
273 theManager.AppendScheme(m_tmp++, passive->GetScheme(p));
274 }
275 }
276 p++;
277 }
278 // Rebuild the Hash;
279 if(theHash.Prepared())
280 {
281 ReHash();
282 }
283 }
void AppendScheme(G4int aPoint, const G4InterpolationScheme &aScheme)
G4InterpolationScheme GetScheme(G4int anIndex)
G4double GetEnergy(G4int i) const

◆ Merge() [2/2]

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

Definition at line 267 of file G4NeutronHPVector.hh.

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 }

Referenced by G4NeutronHPDiscreteTwoBody::Sample(), and G4NeutronHPLabAngularEnergy::Sample().

◆ operator=()

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

Definition at line 120 of file G4NeutronHPVector.cc.

122 {
123 if(&right == this) return *this;
124
125 G4int i;
126
127 totalIntegral = right.totalIntegral;
128 if(right.theIntegral!=0) theIntegral = new G4double[right.nEntries];
129 for(i=0; i<right.nEntries; i++)
130 {
131 SetPoint(i, right.GetPoint(i)); // copy theData
132 if(right.theIntegral!=0) theIntegral[i] = right.theIntegral[i];
133 }
134 theManager = right.theManager;
135 label = right.label;
136
137 Verbose = right.Verbose;
138 the15percentBorderCash = right.the15percentBorderCash;
139 the50percentBorderCash = right.the50percentBorderCash;
140 theHash = right.theHash;
141 return *this;
142 }
void SetPoint(G4int i, const G4NeutronHPDataPoint &it)
const G4NeutronHPDataPoint & GetPoint(G4int i) const

◆ ReHash()

void G4NeutronHPVector::ReHash ( )
inline

Definition at line 144 of file G4NeutronHPVector.hh.

145 {
146 theHash.Clear();
147 Hash();
148 }

Referenced by Merge(), and ThinOut().

◆ Sample()

G4double G4NeutronHPVector::Sample ( )

Definition at line 348 of file G4NeutronHPVector.cc.

349 {
350 G4double result;
351 G4int j;
352 for(j=0; j<GetVectorLength(); j++)
353 {
354 if(GetY(j)<0) SetY(j, 0);
355 }
356
357 if(theBuffered.size() !=0 && G4UniformRand()<0.5)
358 {
359 result = theBuffered[0];
360 theBuffered.erase(theBuffered.begin());
361 if(result < GetX(GetVectorLength()-1) ) return result;
362 }
363 if(GetVectorLength()==1)
364 {
365 result = theData[0].GetX();
366 }
367 else
368 {
369 if(theIntegral==0) { IntegrateAndNormalise(); }
370 do
371 {
372//080808
373/*
374 G4double rand;
375 G4double value, test, baseline;
376 baseline = theData[GetVectorLength()-1].GetX()-theData[0].GetX();
377 do
378 {
379 value = baseline*G4UniformRand();
380 value += theData[0].GetX();
381 test = GetY(value)/maxValue;
382 rand = G4UniformRand();
383 }
384 //while(test<rand);
385 while( test < rand && test > 0 );
386 result = value;
387*/
388 G4double rand;
389 G4double value, test;
390 do
391 {
392 rand = G4UniformRand();
393 G4int ibin = -1;
394 for ( G4int i = 0 ; i < GetVectorLength() ; i++ )
395 {
396 if ( rand < theIntegral[i] )
397 {
398 ibin = i;
399 break;
400 }
401 }
402 if ( ibin < 0 ) G4cout << "TKDB 080807 " << rand << G4endl;
403 // result
404 rand = G4UniformRand();
405 G4double x1, x2;
406 if ( ibin == 0 )
407 {
408 x1 = theData[ ibin ].GetX();
409 value = x1;
410 break;
411 }
412 else
413 {
414 x1 = theData[ ibin-1 ].GetX();
415 }
416
417 x2 = theData[ ibin ].GetX();
418 value = rand * ( x2 - x1 ) + x1;
419 //***********************************************************************
420 /*
421 test = GetY ( value ) / std::max ( GetY( ibin-1 ) , GetY ( ibin ) );
422 */
423 //***********************************************************************
424 //EMendoza - Always linear interpolation:
425 G4double y1=theData[ ibin-1 ].GetY();
426 G4double y2=theData[ ibin ].GetY();
427 G4double mval=(y2-y1)/(x2-x1);
428 G4double bval=y1-mval*x1;
429 test =(mval*value+bval)/std::max ( GetY( ibin-1 ) , GetY ( ibin ) );
430 //***********************************************************************
431 }
432 while ( G4UniformRand() > test );
433 result = value;
434//080807
435 }
436 while(IsBlocked(result));
437 }
438 return result;
439 }
#define G4UniformRand()
Definition: Randomize.hh:53
void SetY(G4int i, G4double x)

Referenced by G4NeutronHPPhotonDist::GetPhotons(), G4NeutronHPArbitaryTab::Sample(), G4NeutronHPEvapSpectrum::Sample(), G4NeutronHPDiscreteTwoBody::Sample(), G4NeutronHPLabAngularEnergy::Sample(), G4NeutronHPContAngularPar::Sample(), and G4NeutronHPPartial::Sample().

◆ SampleLin()

G4double G4NeutronHPVector::SampleLin ( )
inline

Definition at line 318 of file G4NeutronHPVector.hh.

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 }

◆ SetData()

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

Definition at line 90 of file G4NeutronHPVector.hh.

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

Referenced by G4NeutronHPPartial::GetY(), G4NeutronHPChannel::Harmonise(), G4NeutronHPElementData::Harmonise(), G4NeutronIsoIsoCrossSections::Init(), Init(), Merge(), G4NeutronHPDiscreteTwoBody::Sample(), G4NeutronHPLabAngularEnergy::Sample(), G4NeutronHPLegendreStore::Sample(), and SetPoint().

◆ SetEnergy()

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

Definition at line 102 of file G4NeutronHPVector.hh.

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

◆ SetInterpolationManager() [1/2]

void G4NeutronHPVector::SetInterpolationManager ( const G4InterpolationManager aManager)
inline

◆ SetInterpolationManager() [2/2]

void G4NeutronHPVector::SetInterpolationManager ( G4InterpolationManager aMan)
inline

Definition at line 487 of file G4NeutronHPVector.hh.

488 {
489 theManager = aMan;
490 }

◆ SetLabel()

void G4NeutronHPVector::SetLabel ( G4double  aLabel)
inline

Definition at line 245 of file G4NeutronHPVector.hh.

246 {
247 label = aLabel;
248 }

Referenced by G4NeutronHPLabAngularEnergy::Init(), and G4NeutronHPArbitaryTab::Init().

◆ SetPoint()

void G4NeutronHPVector::SetPoint ( G4int  i,
const G4NeutronHPDataPoint it 
)
inline

Definition at line 83 of file G4NeutronHPVector.hh.

84 {
85 G4double x = it.GetX();
86 G4double y = it.GetY();
87 SetData(i, x, y);
88 }

Referenced by G4NeutronHPIsoData::FillChannelData(), and operator=().

◆ SetScheme()

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

Definition at line 492 of file G4NeutronHPVector.hh.

493 {
494 theManager.AppendScheme(aPoint, aScheme);
495 }

Referenced by G4NeutronHPPartial::GetY(), and G4NeutronHPPartial::Sample().

◆ SetVerbose()

void G4NeutronHPVector::SetVerbose ( G4int  ff)
inline

Definition at line 65 of file G4NeutronHPVector.hh.

66 {
67 Verbose = ff;
68 }

◆ SetX()

void G4NeutronHPVector::SetX ( G4int  i,
G4double  e 
)
inline

◆ SetXsec()

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

Definition at line 113 of file G4NeutronHPVector.hh.

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

◆ SetY()

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

Definition at line 107 of file G4NeutronHPVector.hh.

108 {
109 Check(i);
110 if(x>maxValue) maxValue=x;
111 theData[i].SetY(x);
112 }

Referenced by Sample(), G4NeutronHPDiscreteTwoBody::Sample(), G4NeutronHPLabAngularEnergy::Sample(), G4NeutronHPContAngularPar::Sample(), G4NeutronHPPartial::Sample(), and G4NeutronHPPartial::SetY().

◆ ThinOut()

void G4NeutronHPVector::ThinOut ( G4double  precision)

Definition at line 285 of file G4NeutronHPVector.cc.

286 {
287 // anything in there?
288 if(GetVectorLength()==0) return;
289 // make the new vector
290 G4NeutronHPDataPoint * aBuff = new G4NeutronHPDataPoint[nPoints];
291 G4double x, x1, x2, y, y1, y2;
292 G4int count = 0, current = 2, start = 1;
293
294 // First element always goes and is never tested.
295 aBuff[0] = theData[0];
296
297 // Find the rest
298 while(current < GetVectorLength())
299 {
300 x1=aBuff[count].GetX();
301 y1=aBuff[count].GetY();
302 x2=theData[current].GetX();
303 y2=theData[current].GetY();
304 for(G4int j=start; j<current; j++)
305 {
306 x = theData[j].GetX();
307 if(x1-x2 == 0) y = (y2+y1)/2.;
308 else y = theInt.Lin(x, x1, x2, y1, y2);
309 if (std::abs(y-theData[j].GetY())>precision*y)
310 {
311 aBuff[++count] = theData[current-1]; // for this one, everything was fine
312 start = current; // the next candidate
313 break;
314 }
315 }
316 current++ ;
317 }
318 // The last one also always goes, and is never tested.
319 aBuff[++count] = theData[GetVectorLength()-1];
320 delete [] theData;
321 theData = aBuff;
322 nEntries = count+1;
323
324 // Rebuild the Hash;
325 if(theHash.Prepared())
326 {
327 ReHash();
328 }
329 }

Referenced by G4NeutronHPElementData::Init(), and G4NeutronHPIsoData::ThinOut().

◆ Times()

void G4NeutronHPVector::Times ( G4double  factor)
inline

Definition at line 70 of file G4NeutronHPVector.hh.

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 }

Referenced by G4NeutronIsoIsoCrossSections::Init(), and G4NeutronHPChannel::UpdateData().

Friends And Related Function Documentation

◆ operator+

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

Definition at line 37 of file G4NeutronHPVector.cc.

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

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