G4NeutronHPVector.hh

Go to the documentation of this file.
00001 //
00002 // ********************************************************************
00003 // * License and Disclaimer                                           *
00004 // *                                                                  *
00005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
00006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
00007 // * conditions of the Geant4 Software License,  included in the file *
00008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
00009 // * include a list of copyright holders.                             *
00010 // *                                                                  *
00011 // * Neither the authors of this software system, nor their employing *
00012 // * institutes,nor the agencies providing financial support for this *
00013 // * work  make  any representation or  warranty, express or implied, *
00014 // * regarding  this  software system or assume any liability for its *
00015 // * use.  Please see the license in the file  LICENSE  and URL above *
00016 // * for the full disclaimer and the limitation of liability.         *
00017 // *                                                                  *
00018 // * This  code  implementation is the result of  the  scientific and *
00019 // * technical work of the GEANT4 collaboration.                      *
00020 // * By using,  copying,  modifying or  distributing the software (or *
00021 // * any work based  on the software)  you  agree  to acknowledge its *
00022 // * use  in  resulting  scientific  publications,  and indicate your *
00023 // * acceptance of all terms of the Geant4 Software license.          *
00024 // ********************************************************************
00025 //
00026 // 070606 fix with Valgrind by T. Koi
00027 // 080409 Fix div0 error with G4FPE by T. Koi
00028 // 080811 Comment out unused method SetBlocked and SetBuffered
00029 //        Add required cleaning up in CleanUp by T. Koi
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 //    G4cout <<"G4NeutronHPVector::SetData called"<<nPoints<<" "<<nEntries<<G4endl;
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       // Protect against doubled-up x values
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 //080811 TK DB 
00262     delete[] theIntegral;
00263     theIntegral = NULL;
00264   }
00265 
00266   // merges the vectors active and passive into *this
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 //080409 TKDB 
00286         //if( std::abs(std::abs(xp-xa)/xa)<0.001 ) p++;
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       //if(std::abs(GetEnergy(m)-passive->GetEnergy(p))/passive->GetEnergy(p)>0.001)
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() // Samples X according to distribution Y, linear int
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       // this was replaced 
00332 //      for(i=1;i<GetVectorLength();i++)
00333 //      {
00334 //      if(rand<theIntegral[i]/theIntegral[GetVectorLength()-1]) break;
00335 //      }
00336 
00337 // by this (begin)
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 // until this (end)
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) // not really necessary, since the case is excluded by construction
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(); // Samples X according to distribution Y
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         //EMendoza -> the interpolation scheme is not always lin-lin
00389         /*
00390         sum+= 0.5*(theData[i].GetY()+theData[i-1].GetY())*
00391                   (x1-x0);
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() // linear interpolation; use with care
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   void Block(G4double aX)
00522   {
00523     theBlocked.push_back(aX);
00524   }
00525   
00526   void Buffer(G4double aX)
00527   {
00528     theBuffered.push_back(aX);
00529   }
00530 */
00531   
00532   std::vector<G4double> GetBlocked() {return theBlocked;}
00533   std::vector<G4double> GetBuffered() {return theBuffered;}
00534   
00535 //  void SetBlocked(const std::vector<G4double> &aBlocked) {theBlocked = aBlocked;}
00536 //  void SetBuffered(const std::vector<G4double> &aBuffer) {theBuffered = aBuffer;}
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; // the data
00556   G4InterpolationManager theManager; // knows how to interpolate the data.
00557   G4double * theIntegral;
00558   G4int nEntries;
00559   G4int nPoints;
00560   G4double label;
00561   
00562   G4NeutronHPInterpolator theInt;
00563   G4int Verbose;
00564   // debug only
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

Generated on Mon May 27 17:49:04 2013 for Geant4 by  doxygen 1.4.7