G4NeutronHPVector.cc

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 // neutron_hp -- source file
00027 // J.P. Wellisch, Nov-1996
00028 // A prototype of the low energy neutron transport model.
00029 //
00030 // 070523 bug fix for G4FPE_DEBUG on by A. Howard ( and T. Koi)
00031 // 080808 bug fix in Sample() and GetXsec() by T. Koi
00032 //
00033 #include "G4NeutronHPVector.hh"
00034 #include "G4SystemOfUnits.hh"
00035 
00036   // if the ranges do not match, constant extrapolation is used.
00037   G4NeutronHPVector & operator + (G4NeutronHPVector & left, G4NeutronHPVector & right)
00038   {
00039     G4NeutronHPVector * result = new G4NeutronHPVector;
00040     G4int j=0;
00041     G4double x;
00042     G4double y;
00043     G4int running = 0;
00044     for(G4int i=0; i<left.GetVectorLength(); i++)
00045     {
00046       while(j<right.GetVectorLength())
00047       {
00048         if(right.GetX(j)<left.GetX(i)*1.001)
00049         {
00050           x = right.GetX(j);
00051           y = right.GetY(j)+left.GetY(x);
00052           result->SetData(running++, x, y);
00053           j++;
00054         }
00055         //else if(std::abs((right.GetX(j)-left.GetX(i))/(left.GetX(i)+right.GetX(j)))>0.001)
00056         else if( left.GetX(i)+right.GetX(j) == 0 
00057               || std::abs((right.GetX(j)-left.GetX(i))/(left.GetX(i)+right.GetX(j))) > 0.001 )
00058         {
00059           x = left.GetX(i);
00060           y = left.GetY(i)+right.GetY(x);
00061           result->SetData(running++, x, y);
00062           break;
00063         }
00064         else
00065         {
00066           break;
00067         }
00068       }
00069       if(j==right.GetVectorLength())
00070       {
00071         x = left.GetX(i);
00072         y = left.GetY(i)+right.GetY(x);
00073         result->SetData(running++, x, y);     
00074       }
00075     }
00076     result->ThinOut(0.02);
00077     return *result;
00078   }
00079 
00080   G4NeutronHPVector::G4NeutronHPVector()
00081   {
00082     theData = new G4NeutronHPDataPoint[20]; 
00083     nPoints=20;
00084     nEntries=0;
00085     Verbose=0;
00086     theIntegral=0;
00087     totalIntegral=-1;
00088     isFreed = 0;
00089     maxValue = -DBL_MAX;
00090     the15percentBorderCash = -DBL_MAX;
00091     the50percentBorderCash = -DBL_MAX;
00092     label = -DBL_MAX;
00093 
00094   }
00095   
00096   G4NeutronHPVector::G4NeutronHPVector(G4int n)
00097   {
00098     nPoints=std::max(n, 20);
00099     theData = new G4NeutronHPDataPoint[nPoints]; 
00100     nEntries=0;
00101     Verbose=0;
00102     theIntegral=0;
00103     totalIntegral=-1;
00104     isFreed = 0;
00105     maxValue = -DBL_MAX;
00106     the15percentBorderCash = -DBL_MAX;
00107     the50percentBorderCash = -DBL_MAX;
00108   }
00109 
00110   G4NeutronHPVector::~G4NeutronHPVector()
00111   {
00112 //    if(Verbose==1)G4cout <<"G4NeutronHPVector::~G4NeutronHPVector"<<G4endl;
00113       delete [] theData;
00114 //    if(Verbose==1)G4cout <<"Vector: delete theData"<<G4endl;
00115       delete [] theIntegral;
00116 //    if(Verbose==1)G4cout <<"Vector: delete theIntegral"<<G4endl;
00117     isFreed = 1;
00118   }
00119   
00120   G4NeutronHPVector & G4NeutronHPVector::  
00121   operator = (const G4NeutronHPVector & right)
00122   {
00123     if(&right == this) return *this;
00124     
00125     G4int i;
00126     
00127     totalIntegral = right.totalIntegral;
00128     if(right.theIntegral!=0) theIntegral = new G4double[right.nEntries];
00129     for(i=0; i<right.nEntries; i++)
00130     {
00131       SetPoint(i, right.GetPoint(i)); // copy theData
00132       if(right.theIntegral!=0) theIntegral[i] = right.theIntegral[i];
00133     }
00134     theManager = right.theManager; 
00135     label = right.label;
00136   
00137     Verbose = right.Verbose;
00138     the15percentBorderCash = right.the15percentBorderCash;
00139     the50percentBorderCash = right.the50percentBorderCash;
00140     theHash = right.theHash;
00141    return *this;
00142   }
00143 
00144   
00145   G4double G4NeutronHPVector::GetXsec(G4double e) 
00146   {
00147     if(nEntries == 0) return 0;
00148     if(!theHash.Prepared()) Hash();
00149     G4int min = theHash.GetMinIndex(e);
00150     G4int i;
00151     for(i=min ; i<nEntries; i++)
00152     {
00153       //if(theData[i].GetX()>e) break;
00154       if(theData[i].GetX() >= e) break;
00155     }
00156     G4int low = i-1;
00157     G4int high = i;
00158     if(i==0)
00159     {
00160       low = 0;
00161       high = 1;
00162     }
00163     else if(i==nEntries)
00164     {
00165       low = nEntries-2;
00166       high = nEntries-1;
00167     }
00168     G4double y;
00169     if(e<theData[nEntries-1].GetX()) 
00170     {
00171       // Protect against doubled-up x values
00172       //if( (theData[high].GetX()-theData[low].GetX())/theData[high].GetX() < 0.000001)
00173       if ( theData[high].GetX() !=0 
00174        //080808 TKDB
00175        //&&( theData[high].GetX()-theData[low].GetX())/theData[high].GetX() < 0.000001)
00176        &&( std::abs( (theData[high].GetX()-theData[low].GetX())/theData[high].GetX() ) < 0.000001 ) )
00177       {
00178         y = theData[low].GetY();
00179       }
00180       else
00181       {
00182         y = theInt.Interpolate(theManager.GetScheme(high), e, 
00183                                theData[low].GetX(), theData[high].GetX(),
00184                                theData[low].GetY(), theData[high].GetY());
00185       }
00186     }
00187     else
00188     {
00189       y=theData[nEntries-1].GetY();
00190     }
00191     return y;
00192   }
00193 
00194   void G4NeutronHPVector::Dump()
00195   {
00196     G4cout << nEntries<<G4endl;
00197     for(G4int i=0; i<nEntries; i++)
00198     {
00199       G4cout << theData[i].GetX()<<" ";
00200       G4cout << theData[i].GetY()<<" ";
00201 //      if (i!=1&&i==5*(i/5)) G4cout << G4endl;
00202       G4cout << G4endl;
00203     }
00204     G4cout << G4endl;
00205   }
00206   
00207   void G4NeutronHPVector::Check(G4int i)
00208   {
00209     if(i>nEntries) throw G4HadronicException(__FILE__, __LINE__, "Skipped some index numbers in G4NeutronHPVector");
00210     if(i==nPoints)
00211     {
00212       nPoints = static_cast<G4int>(1.2*nPoints);
00213       G4NeutronHPDataPoint * buff = new G4NeutronHPDataPoint[nPoints];
00214       for (G4int j=0; j<nEntries; j++) buff[j] = theData[j];
00215       delete [] theData;
00216       theData = buff;
00217     }
00218     if(i==nEntries) nEntries=i+1;
00219   }
00220 
00221   void G4NeutronHPVector::
00222   Merge(G4InterpolationScheme aScheme, G4double aValue, 
00223         G4NeutronHPVector * active, G4NeutronHPVector * passive)
00224   { 
00225     // interpolate between labels according to aScheme, cut at aValue, 
00226     // continue in unknown areas by substraction of the last difference.
00227     
00228     CleanUp();
00229     G4int s_tmp = 0, n=0, m_tmp=0;
00230     G4NeutronHPVector * tmp;
00231     G4int a = s_tmp, p = n, t;
00232     while ( a<active->GetVectorLength() )
00233     {
00234       if(active->GetEnergy(a) <= passive->GetEnergy(p))
00235       {
00236         G4double xa  = active->GetEnergy(a);
00237         G4double yy = theInt.Interpolate(aScheme, aValue, active->GetLabel(), passive->GetLabel(),
00238                                                           active->GetXsec(a), passive->GetXsec(xa));
00239         SetData(m_tmp, xa, yy);
00240         theManager.AppendScheme(m_tmp, active->GetScheme(a));
00241         m_tmp++;
00242         a++;
00243         G4double xp = passive->GetEnergy(p);
00244         //if( std::abs(std::abs(xp-xa)/xa)<0.0000001&&a<active->GetVectorLength() ) 
00245         if ( xa != 0 
00246           && std::abs(std::abs(xp-xa)/xa) < 0.0000001 
00247           && a < active->GetVectorLength() )
00248         {
00249           p++;
00250           tmp = active; t=a;
00251           active = passive; a=p;
00252           passive = tmp; p=t;
00253         }
00254       } else {
00255         tmp = active; t=a;
00256         active = passive; a=p;
00257         passive = tmp; p=t;
00258       }
00259     }
00260     
00261     G4double deltaX = passive->GetXsec(GetEnergy(m_tmp-1)) - GetXsec(m_tmp-1);
00262     while (p!=passive->GetVectorLength()&&passive->GetEnergy(p)<=aValue)
00263     {
00264       G4double anX;
00265       anX = passive->GetXsec(p)-deltaX;
00266       if(anX>0)
00267       {
00268         //if(std::abs(GetEnergy(m-1)-passive->GetEnergy(p))/passive->GetEnergy(p)>0.0000001)
00269         if ( passive->GetEnergy(p) == 0 
00270           || std::abs(GetEnergy(m_tmp-1)-passive->GetEnergy(p))/passive->GetEnergy(p) > 0.0000001 )
00271         {
00272           SetData(m_tmp, passive->GetEnergy(p), anX);
00273           theManager.AppendScheme(m_tmp++, passive->GetScheme(p));
00274         }
00275       }
00276       p++;
00277     }
00278     // Rebuild the Hash;
00279     if(theHash.Prepared()) 
00280     {
00281       ReHash();
00282     }
00283   }
00284     
00285   void G4NeutronHPVector::ThinOut(G4double precision)
00286   {
00287     // anything in there?
00288     if(GetVectorLength()==0) return;
00289     // make the new vector
00290     G4NeutronHPDataPoint * aBuff = new G4NeutronHPDataPoint[nPoints];
00291     G4double x, x1, x2, y, y1, y2;
00292     G4int count = 0, current = 2, start = 1;
00293     
00294     // First element always goes and is never tested.
00295     aBuff[0] = theData[0];
00296     
00297     // Find the rest
00298     while(current < GetVectorLength())
00299     {
00300       x1=aBuff[count].GetX();
00301       y1=aBuff[count].GetY();
00302       x2=theData[current].GetX();
00303       y2=theData[current].GetY();
00304       for(G4int j=start; j<current; j++)
00305       {
00306         x = theData[j].GetX();
00307         if(x1-x2 == 0) y = (y2+y1)/2.;
00308         else y = theInt.Lin(x, x1, x2, y1, y2);
00309         if (std::abs(y-theData[j].GetY())>precision*y)
00310         {
00311           aBuff[++count] = theData[current-1]; // for this one, everything was fine
00312           start = current; // the next candidate
00313           break;
00314         }
00315       }
00316       current++ ;
00317     }
00318     // The last one also always goes, and is never tested.
00319     aBuff[++count] = theData[GetVectorLength()-1];
00320     delete [] theData;
00321     theData = aBuff;
00322     nEntries = count+1;
00323     
00324     // Rebuild the Hash;
00325     if(theHash.Prepared()) 
00326     {
00327       ReHash();
00328     }
00329   }
00330 
00331   G4bool G4NeutronHPVector::IsBlocked(G4double aX)
00332   {
00333     G4bool result = false;
00334     std::vector<G4double>::iterator i;
00335     for(i=theBlocked.begin(); i!=theBlocked.end(); i++)
00336     {
00337       G4double aBlock = *i;
00338       if(std::abs(aX-aBlock) < 0.1*MeV)
00339       {
00340         result = true;
00341         theBlocked.erase(i);
00342         break;
00343       }
00344     }
00345     return result;
00346   }
00347 
00348   G4double G4NeutronHPVector::Sample() // Samples X according to distribution Y
00349   {
00350     G4double result;
00351     G4int j;
00352     for(j=0; j<GetVectorLength(); j++)
00353     {
00354       if(GetY(j)<0) SetY(j, 0);
00355     }
00356     
00357     if(theBuffered.size() !=0 && G4UniformRand()<0.5) 
00358     {
00359       result = theBuffered[0];
00360       theBuffered.erase(theBuffered.begin());
00361       if(result < GetX(GetVectorLength()-1) ) return result;
00362     }
00363     if(GetVectorLength()==1)
00364     {
00365       result = theData[0].GetX();
00366     }
00367     else
00368     {
00369       if(theIntegral==0) { IntegrateAndNormalise(); }
00370       do
00371       {
00372 //080808
00373 /*
00374         G4double rand;
00375         G4double value, test, baseline;
00376         baseline = theData[GetVectorLength()-1].GetX()-theData[0].GetX();
00377         do
00378         {
00379           value = baseline*G4UniformRand();
00380           value += theData[0].GetX();
00381           test = GetY(value)/maxValue;
00382           rand = G4UniformRand();
00383         }
00384         //while(test<rand);
00385         while( test < rand && test > 0 );
00386         result = value;
00387 */
00388         G4double rand;
00389         G4double value, test;
00390         do 
00391         {
00392            rand = G4UniformRand();
00393            G4int ibin = -1;
00394            for ( G4int i = 0 ; i < GetVectorLength() ; i++ )
00395            {
00396               if ( rand < theIntegral[i] ) 
00397               {
00398                  ibin = i; 
00399                  break;
00400               }
00401            }
00402            if ( ibin < 0 ) G4cout << "TKDB 080807 " << rand << G4endl; 
00403            // result 
00404            rand = G4UniformRand();
00405            G4double x1, x2; 
00406            if ( ibin == 0 ) 
00407            {
00408               x1 = theData[ ibin ].GetX(); 
00409               value = x1; 
00410               break;
00411            }
00412            else 
00413            {
00414               x1 = theData[ ibin-1 ].GetX();
00415            }
00416            
00417            x2 = theData[ ibin ].GetX();
00418            value = rand * ( x2 - x1 ) + x1;
00419            //***********************************************************************
00420            /*
00421            test = GetY ( value ) / std::max ( GetY( ibin-1 ) , GetY ( ibin ) ); 
00422            */
00423            //***********************************************************************
00424            //EMendoza - Always linear interpolation:
00425            G4double y1=theData[ ibin-1 ].GetY();
00426            G4double y2=theData[ ibin ].GetY();
00427            G4double mval=(y2-y1)/(x2-x1);
00428            G4double bval=y1-mval*x1;
00429            test =(mval*value+bval)/std::max ( GetY( ibin-1 ) , GetY ( ibin ) ); 
00430            //***********************************************************************
00431         }
00432         while ( G4UniformRand() > test );
00433         result = value;
00434 //080807
00435       }
00436       while(IsBlocked(result));
00437     }
00438     return result;
00439   }
00440 
00441   G4double G4NeutronHPVector::Get15percentBorder()
00442   {    
00443     if(the15percentBorderCash>-DBL_MAX/2.) return the15percentBorderCash;
00444     G4double result;
00445     if(GetVectorLength()==1)
00446     {
00447       result = theData[0].GetX();
00448       the15percentBorderCash = result;
00449     }
00450     else
00451     {
00452       if(theIntegral==0) { IntegrateAndNormalise(); }
00453       G4int i;
00454       result = theData[GetVectorLength()-1].GetX();
00455       for(i=0;i<GetVectorLength();i++)
00456       {
00457         if(theIntegral[i]/theIntegral[GetVectorLength()-1]>0.15) 
00458         {
00459           result = theData[std::min(i+1, GetVectorLength()-1)].GetX();
00460           the15percentBorderCash = result;
00461           break;
00462         }
00463       }
00464       the15percentBorderCash = result;
00465     }
00466     return result;
00467   }
00468 
00469   G4double G4NeutronHPVector::Get50percentBorder()
00470   {    
00471     if(the50percentBorderCash>-DBL_MAX/2.) return the50percentBorderCash;
00472     G4double result;
00473     if(GetVectorLength()==1)
00474     {
00475       result = theData[0].GetX();
00476       the50percentBorderCash = result;
00477     }
00478     else
00479     {
00480       if(theIntegral==0) { IntegrateAndNormalise(); }
00481       G4int i;
00482       G4double x = 0.5;
00483       result = theData[GetVectorLength()-1].GetX();
00484       for(i=0;i<GetVectorLength();i++)
00485       {
00486         if(theIntegral[i]/theIntegral[GetVectorLength()-1]>x) 
00487         {
00488           G4int it;
00489           it = i;
00490           if(it == GetVectorLength()-1)
00491           {
00492             result = theData[GetVectorLength()-1].GetX();
00493           }
00494           else
00495           {
00496             G4double x1, x2, y1, y2;
00497             x1 = theIntegral[i-1]/theIntegral[GetVectorLength()-1];
00498             x2 = theIntegral[i]/theIntegral[GetVectorLength()-1];
00499             y1 = theData[i-1].GetX();
00500             y2 = theData[i].GetX();
00501             result = theLin.Lin(x, x1, x2, y1, y2);
00502           }
00503           the50percentBorderCash = result;
00504           break;
00505         }
00506       }
00507       the50percentBorderCash = result;
00508     }
00509     return result;
00510   }

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