00001
00002
00003
00004
00005
00006
00007
00008
00009
00010
00011
00012
00013
00014
00015
00016
00017
00018
00019
00020
00021
00022
00023
00024
00025
00026
00027
00028
00029
00030
00031 #ifndef G4NeutronHPVector_h
00032 #define G4NeutronHPVector_h 1
00033
00034 #include "G4NeutronHPDataPoint.hh"
00035 #include "G4PhysicsVector.hh"
00036 #include "G4NeutronHPInterpolator.hh"
00037 #include "Randomize.hh"
00038 #include "G4ios.hh"
00039 #include <fstream>
00040 #include "G4InterpolationManager.hh"
00041 #include "G4NeutronHPInterpolator.hh"
00042 #include "G4NeutronHPHash.hh"
00043 #include <cmath>
00044 #include <vector>
00045
00046 #if defined WIN32-VC
00047 #include <float.h>
00048 #endif
00049
00050 class G4NeutronHPVector
00051 {
00052 friend G4NeutronHPVector & operator + (G4NeutronHPVector & left,
00053 G4NeutronHPVector & right);
00054
00055 public:
00056
00057 G4NeutronHPVector();
00058
00059 G4NeutronHPVector(G4int n);
00060
00061 ~G4NeutronHPVector();
00062
00063 G4NeutronHPVector & operator = (const G4NeutronHPVector & right);
00064
00065 inline void SetVerbose(G4int ff)
00066 {
00067 Verbose = ff;
00068 }
00069
00070 inline void Times(G4double factor)
00071 {
00072 G4int i;
00073 for(i=0; i<nEntries; i++)
00074 {
00075 theData[i].SetY(theData[i].GetY()*factor);
00076 }
00077 if(theIntegral!=0)
00078 {
00079 theIntegral[i] *= factor;
00080 }
00081 }
00082
00083 inline void SetPoint(G4int i, const G4NeutronHPDataPoint & it)
00084 {
00085 G4double x = it.GetX();
00086 G4double y = it.GetY();
00087 SetData(i, x, y);
00088 }
00089
00090 inline void SetData(G4int i, G4double x, G4double y)
00091 {
00092
00093 Check(i);
00094 if(y>maxValue) maxValue=y;
00095 theData[i].SetData(x, y);
00096 }
00097 inline void SetX(G4int i, G4double e)
00098 {
00099 Check(i);
00100 theData[i].SetX(e);
00101 }
00102 inline void SetEnergy(G4int i, G4double e)
00103 {
00104 Check(i);
00105 theData[i].SetX(e);
00106 }
00107 inline void SetY(G4int i, G4double x)
00108 {
00109 Check(i);
00110 if(x>maxValue) maxValue=x;
00111 theData[i].SetY(x);
00112 }
00113 inline void SetXsec(G4int i, G4double x)
00114 {
00115 Check(i);
00116 if(x>maxValue) maxValue=x;
00117 theData[i].SetY(x);
00118 }
00119 inline G4double GetEnergy(G4int i) const { return theData[i].GetX(); }
00120 inline G4double GetXsec(G4int i) { return theData[i].GetY(); }
00121 inline G4double GetX(G4int i) const
00122 {
00123 if (i<0) i=0;
00124 if(i>=GetVectorLength()) i=GetVectorLength()-1;
00125 return theData[i].GetX();
00126 }
00127 inline const G4NeutronHPDataPoint & GetPoint(G4int i) const { return theData[i]; }
00128
00129 void Hash()
00130 {
00131 G4int i;
00132 G4double x, y;
00133 for(i=0 ; i<nEntries; i++)
00134 {
00135 if(0 == (i+1)%10)
00136 {
00137 x = GetX(i);
00138 y = GetY(i);
00139 theHash.SetData(i, x, y);
00140 }
00141 }
00142 }
00143
00144 void ReHash()
00145 {
00146 theHash.Clear();
00147 Hash();
00148 }
00149
00150 G4double GetXsec(G4double e);
00151 G4double GetXsec(G4double e, G4int min)
00152 {
00153 G4int i;
00154 for(i=min ; i<nEntries; i++)
00155 {
00156 if(theData[i].GetX()>e) break;
00157 }
00158 G4int low = i-1;
00159 G4int high = i;
00160 if(i==0)
00161 {
00162 low = 0;
00163 high = 1;
00164 }
00165 else if(i==nEntries)
00166 {
00167 low = nEntries-2;
00168 high = nEntries-1;
00169 }
00170 G4double y;
00171 if(e<theData[nEntries-1].GetX())
00172 {
00173
00174 if( (theData[high].GetX()-theData[low].GetX())/theData[high].GetX() < 0.000001)
00175 {
00176 y = theData[low].GetY();
00177 }
00178 else
00179 {
00180 y = theInt.Interpolate(theManager.GetScheme(high), e,
00181 theData[low].GetX(), theData[high].GetX(),
00182 theData[low].GetY(), theData[high].GetY());
00183 }
00184 }
00185 else
00186 {
00187 y=theData[nEntries-1].GetY();
00188 }
00189 return y;
00190 }
00191
00192 inline G4double GetY(G4double x) {return GetXsec(x);}
00193 inline G4int GetVectorLength() const {return nEntries;}
00194
00195 inline G4double GetY(G4int i)
00196 {
00197 if (i<0) i=0;
00198 if(i>=GetVectorLength()) i=GetVectorLength()-1;
00199 return theData[i].GetY();
00200 }
00201
00202 inline G4double GetY(G4int i) const
00203 {
00204 if (i<0) i=0;
00205 if(i>=GetVectorLength()) i=GetVectorLength()-1;
00206 return theData[i].GetY();
00207 }
00208 void Dump();
00209
00210 inline void InitInterpolation(std::ifstream & aDataFile)
00211 {
00212 theManager.Init(aDataFile);
00213 }
00214
00215 void Init(std::ifstream & aDataFile, G4int total, G4double ux=1., G4double uy=1.)
00216 {
00217 G4double x,y;
00218 for (G4int i=0;i<total;i++)
00219 {
00220 aDataFile >> x >> y;
00221 x*=ux;
00222 y*=uy;
00223 SetData(i,x,y);
00224 if(0 == nEntries%10)
00225 {
00226 theHash.SetData(nEntries-1, x, y);
00227 }
00228 }
00229 }
00230
00231 void Init(std::ifstream & aDataFile,G4double ux=1., G4double uy=1.)
00232 {
00233 G4int total;
00234 aDataFile >> total;
00235 if(theData!=0) delete [] theData;
00236 theData = new G4NeutronHPDataPoint[total];
00237 nPoints=total;
00238 nEntries=0;
00239 theManager.Init(aDataFile);
00240 Init(aDataFile, total, ux, uy);
00241 }
00242
00243 void ThinOut(G4double precision);
00244
00245 inline void SetLabel(G4double aLabel)
00246 {
00247 label = aLabel;
00248 }
00249
00250 inline G4double GetLabel()
00251 {
00252 return label;
00253 }
00254
00255 inline void CleanUp()
00256 {
00257 nEntries=0;
00258 theManager.CleanUp();
00259 maxValue = -DBL_MAX;
00260 theHash.Clear();
00261
00262 delete[] theIntegral;
00263 theIntegral = NULL;
00264 }
00265
00266
00267 inline void Merge(G4NeutronHPVector * active, G4NeutronHPVector * passive)
00268 {
00269 CleanUp();
00270 G4int s_tmp = 0, n=0, m_tmp=0;
00271 G4NeutronHPVector * tmp;
00272 G4int a = s_tmp, p = n, t;
00273 while (a<active->GetVectorLength()&&p<passive->GetVectorLength())
00274 {
00275 if(active->GetEnergy(a) <= passive->GetEnergy(p))
00276 {
00277 G4double xa = active->GetEnergy(a);
00278 G4double yy = active->GetXsec(a);
00279 SetData(m_tmp, xa, yy);
00280 theManager.AppendScheme(m_tmp, active->GetScheme(a));
00281 m_tmp++;
00282 a++;
00283 G4double xp = passive->GetEnergy(p);
00284
00285
00286
00287 if ( !( xa == 0 ) && std::abs(std::abs(xp-xa)/xa)<0.001 ) p++;
00288 } else {
00289 tmp = active;
00290 t=a;
00291 active = passive;
00292 a=p;
00293 passive = tmp;
00294 p=t;
00295 }
00296 }
00297 while (a!=active->GetVectorLength())
00298 {
00299 SetData(m_tmp, active->GetEnergy(a), active->GetXsec(a));
00300 theManager.AppendScheme(m_tmp++, active->GetScheme(a));
00301 a++;
00302 }
00303 while (p!=passive->GetVectorLength())
00304 {
00305 if(std::abs(GetEnergy(m_tmp-1)-passive->GetEnergy(p))/passive->GetEnergy(p)>0.001)
00306
00307 {
00308 SetData(m_tmp, passive->GetEnergy(p), passive->GetXsec(p));
00309 theManager.AppendScheme(m_tmp++, active->GetScheme(p));
00310 }
00311 p++;
00312 }
00313 }
00314
00315 void Merge(G4InterpolationScheme aScheme, G4double aValue,
00316 G4NeutronHPVector * active, G4NeutronHPVector * passive);
00317
00318 G4double SampleLin()
00319 {
00320 G4double result;
00321 if(theIntegral==0) IntegrateAndNormalise();
00322 if(GetVectorLength()==1)
00323 {
00324 result = theData[0].GetX();
00325 }
00326 else
00327 {
00328 G4int i;
00329 G4double rand = G4UniformRand();
00330
00331
00332
00333
00334
00335
00336
00337
00338 for(i=GetVectorLength()-1; i>=0 ;i--)
00339 {
00340 if(rand>theIntegral[i]/theIntegral[GetVectorLength()-1]) break;
00341 }
00342 if(i!=GetVectorLength()-1) i++;
00343
00344
00345 G4double x1, x2, y1, y2;
00346 y1 = theData[i-1].GetX();
00347 x1 = theIntegral[i-1];
00348 y2 = theData[i].GetX();
00349 x2 = theIntegral[i];
00350 if(std::abs((y2-y1)/y2)<0.0000001)
00351 {
00352 y1 = theData[i-2].GetX();
00353 x1 = theIntegral[i-2];
00354 }
00355 result = theLin.Lin(rand, x1, x2, y1, y2);
00356 }
00357 return result;
00358 }
00359
00360 G4double Sample();
00361
00362 G4double * Debug()
00363 {
00364 return theIntegral;
00365 }
00366
00367 inline void IntegrateAndNormalise()
00368 {
00369 G4int i;
00370 if(theIntegral!=0) return;
00371 theIntegral = new G4double[nEntries];
00372 if(nEntries == 1)
00373 {
00374 theIntegral[0] = 1;
00375 return;
00376 }
00377 theIntegral[0] = 0;
00378 G4double sum = 0;
00379 G4double x1 = 0;
00380 G4double x0 = 0;
00381 for(i=1;i<GetVectorLength();i++)
00382 {
00383 x1 = theData[i].GetX();
00384 x0 = theData[i-1].GetX();
00385 if (std::abs(x1-x0) > std::abs(x1*0.0000001) )
00386 {
00387
00388
00389
00390
00391
00392
00393
00394 G4InterpolationScheme aScheme = theManager.GetScheme(i);
00395 G4double y0 = theData[i-1].GetY();
00396 G4double y1 = theData[i].GetY();
00397 G4double integ=theInt.GetBinIntegral(aScheme,x0,x1,y0,y1);
00398 #if defined WIN32-VC
00399 if(!_finite(integ)){integ=0;}
00400 #elif defined __IBMCPP__
00401 if(isinf(integ)||isnan(integ)){integ=0;}
00402 #else
00403 if(std::isinf(integ)||std::isnan(integ)){integ=0;}
00404 #endif
00405 sum+=integ;
00406
00407 }
00408 theIntegral[i] = sum;
00409 }
00410 G4double total = theIntegral[GetVectorLength()-1];
00411 for(i=1;i<GetVectorLength();i++)
00412 {
00413 theIntegral[i]/=total;
00414 }
00415 }
00416
00417 inline void Integrate()
00418 {
00419 G4int i;
00420 if(nEntries == 1)
00421 {
00422 totalIntegral = 0;
00423 return;
00424 }
00425 G4double sum = 0;
00426 for(i=1;i<GetVectorLength();i++)
00427 {
00428 if(std::abs((theData[i].GetX()-theData[i-1].GetX())/theData[i].GetX())>0.0000001)
00429 {
00430 G4double x1 = theData[i-1].GetX();
00431 G4double x2 = theData[i].GetX();
00432 G4double y1 = theData[i-1].GetY();
00433 G4double y2 = theData[i].GetY();
00434 G4InterpolationScheme aScheme = theManager.GetScheme(i);
00435 if(aScheme==LINLIN||aScheme==CLINLIN||aScheme==ULINLIN)
00436 {
00437 sum+= 0.5*(y2+y1)*(x2-x1);
00438 }
00439 else if(aScheme==LINLOG||aScheme==CLINLOG||aScheme==ULINLOG)
00440 {
00441 G4double a = y1;
00442 G4double b = (y2-y1)/(std::log(x2)-std::log(x1));
00443 sum+= (a-b)*(x2-x1) + b*(x2*std::log(x2)-x1*std::log(x1));
00444 }
00445 else if(aScheme==LOGLIN||aScheme==CLOGLIN||aScheme==ULOGLIN)
00446 {
00447 G4double a = std::log(y1);
00448 G4double b = (std::log(y2)-std::log(y1))/(x2-x1);
00449 sum += (std::exp(a)/b)*(std::exp(b*x2)-std::exp(b*x1));
00450 }
00451 else if(aScheme==HISTO||aScheme==CHISTO||aScheme==UHISTO)
00452 {
00453 sum+= y1*(x2-x1);
00454 }
00455 else if(aScheme==LOGLOG||aScheme==CLOGLOG||aScheme==ULOGLOG)
00456 {
00457 G4double a = std::log(y1);
00458 G4double b = (std::log(y2)-std::log(y1))/(std::log(x2)-std::log(x1));
00459 sum += (std::exp(a)/(b+1))*(std::pow(x2,b+1)-std::pow(x1,b+1));
00460 }
00461 else
00462 {
00463 throw G4HadronicException(__FILE__, __LINE__, "Unknown interpolation scheme in G4NeutronHPVector::Integrate");
00464 }
00465
00466 }
00467 }
00468 totalIntegral = sum;
00469 }
00470
00471 inline G4double GetIntegral()
00472 {
00473 if(totalIntegral<-0.5) Integrate();
00474 return totalIntegral;
00475 }
00476
00477 inline void SetInterpolationManager(const G4InterpolationManager & aManager)
00478 {
00479 theManager = aManager;
00480 }
00481
00482 inline const G4InterpolationManager & GetInterpolationManager() const
00483 {
00484 return theManager;
00485 }
00486
00487 inline void SetInterpolationManager(G4InterpolationManager & aMan)
00488 {
00489 theManager = aMan;
00490 }
00491
00492 inline void SetScheme(G4int aPoint, const G4InterpolationScheme & aScheme)
00493 {
00494 theManager.AppendScheme(aPoint, aScheme);
00495 }
00496
00497 inline G4InterpolationScheme GetScheme(G4int anIndex)
00498 {
00499 return theManager.GetScheme(anIndex);
00500 }
00501
00502 G4double GetMeanX()
00503 {
00504 G4double result;
00505 G4double running = 0;
00506 G4double weighted = 0;
00507 for(G4int i=1; i<nEntries; i++)
00508 {
00509 running += theInt.GetBinIntegral(theManager.GetScheme(i-1),
00510 theData[i-1].GetX(), theData[i].GetX(),
00511 theData[i-1].GetY(), theData[i].GetY());
00512 weighted += theInt.GetWeightedBinIntegral(theManager.GetScheme(i-1),
00513 theData[i-1].GetX(), theData[i].GetX(),
00514 theData[i-1].GetY(), theData[i].GetY());
00515 }
00516 result = weighted / running;
00517 return result;
00518 }
00519
00520
00521
00522
00523
00524
00525
00526
00527
00528
00529
00530
00531
00532 std::vector<G4double> GetBlocked() {return theBlocked;}
00533 std::vector<G4double> GetBuffered() {return theBuffered;}
00534
00535
00536
00537
00538 G4double Get15percentBorder();
00539 G4double Get50percentBorder();
00540
00541 private:
00542
00543 void Check(G4int i);
00544
00545 G4bool IsBlocked(G4double aX);
00546
00547 private:
00548
00549 G4NeutronHPInterpolator theLin;
00550
00551 private:
00552
00553 G4double totalIntegral;
00554
00555 G4NeutronHPDataPoint * theData;
00556 G4InterpolationManager theManager;
00557 G4double * theIntegral;
00558 G4int nEntries;
00559 G4int nPoints;
00560 G4double label;
00561
00562 G4NeutronHPInterpolator theInt;
00563 G4int Verbose;
00564
00565 G4int isFreed;
00566
00567 G4NeutronHPHash theHash;
00568 G4double maxValue;
00569
00570 std::vector<G4double> theBlocked;
00571 std::vector<G4double> theBuffered;
00572 G4double the15percentBorderCash;
00573 G4double the50percentBorderCash;
00574
00575 };
00576
00577 #endif