G4ParticleHPVector Class Reference

#include <G4ParticleHPVector.hh>

Public Member Functions
	G4ParticleHPVector ()
	G4ParticleHPVector (G4int n)
	~G4ParticleHPVector ()
G4ParticleHPVector &	operator= (const G4ParticleHPVector &right)
void	SetVerbose (G4int ff)
void	Times (G4double factor)
void	SetPoint (G4int i, const G4ParticleHPDataPoint &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 G4ParticleHPDataPoint &	GetPoint (G4int i) const
void	Hash ()
void	ReHash ()
G4double	GetXsec (G4double e)
G4int	GetEnergyIndex (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::istream &aDataFile)
void	Init (std::istream &aDataFile, G4int total, G4double ux=1., G4double uy=1.)
void	Init (std::istream &aDataFile, G4double ux=1., G4double uy=1.)
void	ThinOut (G4double precision)
void	SetLabel (G4double aLabel)
G4double	GetLabel ()
void	CleanUp ()
void	Merge (G4ParticleHPVector *active, G4ParticleHPVector *passive)
void	Merge (G4InterpolationScheme aScheme, G4double aValue, G4ParticleHPVector *active, G4ParticleHPVector *passive)
G4double	SampleLin ()
G4double	Sample ()
G4double *	Debug ()
void	IntegrateAndNormalise ()
void	Integrate ()
G4double	GetIntegral ()
void	SetInterpolationManager (const G4InterpolationManager &aManager)
const G4InterpolationManager &	GetInterpolationManager () const
void	SetInterpolationManager (G4InterpolationManager &aMan)
void	SetScheme (G4int aPoint, const G4InterpolationScheme &aScheme)
G4InterpolationScheme	GetScheme (G4int anIndex)
G4double	GetMeanX ()
G4double	GetMaxY (G4double emin, G4double emax)
std::vector< G4double >	GetBlocked ()
std::vector< G4double >	GetBuffered ()
G4double	Get15percentBorder ()
G4double	Get50percentBorder ()

Friends
G4ParticleHPVector &	operator+ (G4ParticleHPVector &left, G4ParticleHPVector &right)

Detailed Description

Definition at line 56 of file G4ParticleHPVector.hh.

Constructor & Destructor Documentation

G4ParticleHPVector() [1/2]

G4ParticleHPVector::G4ParticleHPVector

(

)

Definition at line 80 of file G4ParticleHPVector.cc.

{
  theData = new G4ParticleHPDataPoint[20];
  nPoints = 20;
  nEntries = 0;
  Verbose = 0;
  theIntegral = nullptr;
  totalIntegral = -1;
  isFreed = 0;
  maxValue = -DBL_MAX;
  the15percentBorderCash = -DBL_MAX;
  the50percentBorderCash = -DBL_MAX;
  label = -DBL_MAX;
}

Referenced by Merge(), Merge(), operator+, and operator=().

G4ParticleHPVector() [2/2]

G4ParticleHPVector::G4ParticleHPVector

(

G4int

n

)

Definition at line 95 of file G4ParticleHPVector.cc.

{
  nPoints = std::max(n, 20);
  theData = new G4ParticleHPDataPoint[nPoints];
  nEntries = 0;
  Verbose = 0;
  theIntegral = nullptr;
  totalIntegral = -1;
  isFreed = 0;
  maxValue = -DBL_MAX;
  the15percentBorderCash = -DBL_MAX;
  the50percentBorderCash = -DBL_MAX;
  label = -DBL_MAX;
}

~G4ParticleHPVector()

G4ParticleHPVector::~G4ParticleHPVector

(

)

Definition at line 110 of file G4ParticleHPVector.cc.

{
  //    if(Verbose==1)G4cout <<"G4ParticleHPVector::~G4ParticleHPVector"<<G4endl;
  delete[] theData;
  //    if(Verbose==1)G4cout <<"Vector: delete theData"<<G4endl;
  delete[] theIntegral;
  //    if(Verbose==1)G4cout <<"Vector: delete theIntegral"<<G4endl;
  theHash.Clear();
  isFreed = 1;
}

Member Function Documentation

CleanUp()

void G4ParticleHPVector::CleanUp ( )

inline

Definition at line 236 of file G4ParticleHPVector.hh.

    {
      nEntries = 0;
      theManager.CleanUp();
      maxValue = -DBL_MAX;
      theHash.Clear();
      // 080811 TK DB
      delete[] theIntegral;
      theIntegral = nullptr;
    }

Referenced by Merge(), and Merge().

Debug()

G4double * G4ParticleHPVector::Debug ( )

inline

Definition at line 342 of file G4ParticleHPVector.hh.

{ return theIntegral; }

Dump()

void G4ParticleHPVector::Dump

(

)

Definition at line 206 of file G4ParticleHPVector.cc.

{
  G4cout << nEntries << G4endl;
  for (G4int i = 0; i < nEntries; i++) {
    G4cout << theData[i].GetX() << " ";
    G4cout << theData[i].GetY() << " ";
    //      if (i!=1&&i==5*(i/5)) G4cout << G4endl;
    G4cout << G4endl;
  }
  G4cout << G4endl;
}

Get15percentBorder()

G4double G4ParticleHPVector::Get15percentBorder

(

)

Definition at line 472 of file G4ParticleHPVector.cc.

{
  if (the15percentBorderCash > -DBL_MAX / 2.) return the15percentBorderCash;
  G4double result;
  if (GetVectorLength() == 1) {
    result = theData[0].GetX();
    the15percentBorderCash = result;
  }
  else {
    if (theIntegral == nullptr) {
      IntegrateAndNormalise();
    }
    G4int i;
    result = theData[GetVectorLength() - 1].GetX();
    for (i = 0; i < GetVectorLength(); i++) {
      if (theIntegral[i] / theIntegral[GetVectorLength() - 1] > 0.15) {
        result = theData[std::min(i + 1, GetVectorLength() - 1)].GetX();
        the15percentBorderCash = result;
        break;
      }
    }
    the15percentBorderCash = result;
  }
  return result;
}

Get50percentBorder()

G4double G4ParticleHPVector::Get50percentBorder

(

)

Definition at line 498 of file G4ParticleHPVector.cc.

{
  if (the50percentBorderCash > -DBL_MAX / 2.) return the50percentBorderCash;
  G4double result;
  if (GetVectorLength() == 1) {
    result = theData[0].GetX();
    the50percentBorderCash = result;
  }
  else {
    if (theIntegral == nullptr) {
      IntegrateAndNormalise();
    }
    G4int i;
    G4double x = 0.5;
    result = theData[GetVectorLength() - 1].GetX();
    for (i = 0; i < GetVectorLength(); i++) {
      if (theIntegral[i] / theIntegral[GetVectorLength() - 1] > x) {
        G4int it;
        it = i;
        if (it == GetVectorLength() - 1) {
          result = theData[GetVectorLength() - 1].GetX();
        }
        else {
          G4double x1, x2, y1, y2;
          x1 = theIntegral[i - 1] / theIntegral[GetVectorLength() - 1];
          x2 = theIntegral[i] / theIntegral[GetVectorLength() - 1];
          y1 = theData[i - 1].GetX();
          y2 = theData[i].GetX();
          result = theLin.Lin(x, x1, x2, y1, y2);
        }
        the50percentBorderCash = result;
        break;
      }
    }
    the50percentBorderCash = result;
  }
  return result;
}

GetBlocked()

std::vector< G4double > G4ParticleHPVector::GetBlocked ( )

inline

Definition at line 497 of file G4ParticleHPVector.hh.

{ return theBlocked; }

GetBuffered()

std::vector< G4double > G4ParticleHPVector::GetBuffered ( )

inline

Definition at line 498 of file G4ParticleHPVector.hh.

{ return theBuffered; }

GetEnergy()

G4double G4ParticleHPVector::GetEnergy ( G4int i ) const

inline

Definition at line 118 of file G4ParticleHPVector.hh.

{ return theData[i].GetX(); }

Referenced by G4ParticleHPChannel::Harmonise(), G4ParticleHPElementData::Harmonise(), Merge(), and Merge().

GetEnergyIndex()

G4int G4ParticleHPVector::GetEnergyIndex

(

G4double &

e

)

Definition at line 194 of file G4ParticleHPVector.cc.

{
  // returns energy index of the bin in which ekin is lower then emax and higher than emin for CALENDF
  // returns energy index of emax for NJOY
  if (nEntries == 0) return 0;
  if ( ( ! theHash.Prepared() ) && ( ! G4Threading::IsWorkerThread() ) ) Hash();
  G4int min = theHash.GetMinIndex(e);
  G4int i = 0;
  for (i=min ; i < nEntries; i++) if (theData[i].GetX() >= e) break;
  return i;
}

GetIntegral()

G4double G4ParticleHPVector::GetIntegral ( )

inline

Definition at line 443 of file G4ParticleHPVector.hh.

    {
      if (totalIntegral < -0.5) Integrate();
      return totalIntegral;
    }

GetInterpolationManager()

const G4InterpolationManager & G4ParticleHPVector::GetInterpolationManager ( ) const

inline

Definition at line 454 of file G4ParticleHPVector.hh.

{ return theManager; }

Referenced by G4ParticleHPLabAngularEnergy::Sample().

GetLabel()

G4double G4ParticleHPVector::GetLabel ( )

inline

Definition at line 234 of file G4ParticleHPVector.hh.

{ return label; }

Referenced by Merge().

GetMaxY()

G4double G4ParticleHPVector::GetMaxY	(	G4double	emin,
		G4double	emax )

Definition at line 537 of file G4ParticleHPVector.cc.

{
  G4double xsmax = 0.;
  if (emin > emax || nEntries == 0) return xsmax;
  if (emin >= theData[nEntries - 1].GetX()) {
    xsmax = theData[nEntries - 1].GetY();
    return xsmax;
  }
  if (emax <= theData[0].GetX()) {
    xsmax = theData[0].GetY();
    return xsmax;
  }
  if (!theHash.Prepared() && !G4Threading::IsWorkerThread()) Hash();
  // Find the lowest index, low, where x(energy) is higher than emin
  G4int i = theHash.GetMinIndex(emin);
  for (; i < nEntries; ++i) {
    if (theData[i].GetX() >= emin) break;
  }
  G4int low = i;
  // Find the lowest index, high, where x(energy) is higher than emax
  i = theHash.GetMinIndex(emax);
  for (; i < nEntries; ++i) {
    if (theData[i].GetX() >= emax) break;
  }
  G4int high = i;
  xsmax = GetXsec(emin);  // Set xsmax as low border
  // Find the highest cross-section
  for (i = low; i < high; ++i) {
    if (xsmax < theData[i].GetY()) xsmax = theData[i].GetY();
  }
  // Check if it is smaller than the high border (e.g. as for a monotonic increasing cross-section)
  G4double highborder = GetXsec(emax);
  if (xsmax < highborder) xsmax = highborder;
 
  if (xsmax == 0.) {
    throw G4HadronicException(__FILE__, __LINE__,
                              "G4ParticleHPVector::GetMaxY : called "
                              "G4Nucleus::GetBiasedThermalNucleus for DBRC, xsmax==0.");
  }
 
  return xsmax;
}

GetMeanX()

G4double G4ParticleHPVector::GetMeanX ( )

inline

Definition at line 465 of file G4ParticleHPVector.hh.

    {
      G4double result;
      G4double running = 0;
      G4double weighted = 0;
      for (G4int i = 1; i < nEntries; i++) {
        running +=
          theInt.GetBinIntegral(theManager.GetScheme(i - 1), theData[i - 1].GetX(),
                                theData[i].GetX(), theData[i - 1].GetY(), theData[i].GetY());
        weighted += theInt.GetWeightedBinIntegral(theManager.GetScheme(i - 1),
                                                  theData[i - 1].GetX(), theData[i].GetX(),
                                                  theData[i - 1].GetY(), theData[i].GetY());
      }
      result = weighted / running;
      return result;
    }

Referenced by G4ParticleHPLabAngularEnergy::Sample().

GetPoint()

const G4ParticleHPDataPoint & G4ParticleHPVector::GetPoint ( G4int i ) const

inline

Definition at line 126 of file G4ParticleHPVector.hh.

{ return theData[i]; }

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

GetScheme()

G4InterpolationScheme G4ParticleHPVector::GetScheme ( G4int anIndex )

inline

Definition at line 463 of file G4ParticleHPVector.hh.

{ return theManager.GetScheme(anIndex); }

Referenced by Merge(), and Merge().

GetVectorLength()

G4int G4ParticleHPVector::GetVectorLength ( ) const

inline

Definition at line 185 of file G4ParticleHPVector.hh.

{ return nEntries; }

Referenced by G4ParticleHPData::DoPhysicsVector(), G4ParticleHPIsoData::FillChannelData(), Get15percentBorder(), Get50percentBorder(), G4ParticleHPPhotonDist::GetPhotons(), GetX(), GetY(), GetY(), G4ParticleHPChannel::Harmonise(), G4ParticleHPElementData::Harmonise(), Integrate(), IntegrateAndNormalise(), Merge(), Merge(), operator+, G4ParticleHPContAngularPar::Sample(), G4ParticleHPDiscreteTwoBody::Sample(), G4ParticleHPLabAngularEnergy::Sample(), Sample(), SampleLin(), and ThinOut().

GetX()

G4double G4ParticleHPVector::GetX ( G4int i ) const

inline

Definition at line 120 of file G4ParticleHPVector.hh.

    {
      if (i < 0) i = 0;
      if (i >= GetVectorLength()) i = GetVectorLength() - 1;
      return theData[i].GetX();
    }

Referenced by G4ParticleHPData::DoPhysicsVector(), GetEnergyIndex(), GetMaxY(), G4ParticleHPPhotonDist::GetPhotons(), GetXsec(), GetXsec(), Hash(), Integrate(), operator+, G4ParticleHPContAngularPar::Sample(), G4ParticleHPDiscreteTwoBody::Sample(), G4ParticleHPLabAngularEnergy::Sample(), G4ParticleHPLegendreStore::Sample(), and Sample().

GetXsec() [1/3]

G4double G4ParticleHPVector::GetXsec

(

G4double

e

)

Definition at line 143 of file G4ParticleHPVector.cc.

{
  if (nEntries == 0) return 0;
  // if(!theHash.Prepared()) Hash();
  if (!theHash.Prepared()) {
    if (G4Threading::IsWorkerThread()) {
      ;
    }
    else {
      Hash();
    }
  }
  G4int min = theHash.GetMinIndex(e);
  G4int i;
  for (i = min; i < nEntries; i++) {
    // if(theData[i].GetX()>e) break;
    if (theData[i].GetX() >= e) break;
  }
  G4int low = i - 1;
  G4int high = i;
  if (i == 0) {
    low = 0;
    high = 1;
  }
  else if (i == nEntries) {
    low = nEntries - 2;
    high = nEntries - 1;
  }
  G4double y;
  if (e < theData[nEntries - 1].GetX()) {
    // Protect against doubled-up x values
    // if( (theData[high].GetX()-theData[low].GetX())/theData[high].GetX() < 0.000001)
    if (theData[high].GetX() != 0
        // 080808 TKDB
        //&&( theData[high].GetX()-theData[low].GetX())/theData[high].GetX() < 0.000001)
        && (std::abs((theData[high].GetX() - theData[low].GetX()) / theData[high].GetX())
            < 0.000001))
    {
      y = theData[low].GetY();
    }
    else {
      y = theInt.Interpolate(theManager.GetScheme(high), e, theData[low].GetX(),
                             theData[high].GetX(), theData[low].GetY(), theData[high].GetY());
    }
  }
  else {
    y = theData[nEntries - 1].GetY();
  }
  return y;
}

GetXsec() [2/3]

G4double G4ParticleHPVector::GetXsec ( G4double e, G4int min )

inline

Definition at line 151 of file G4ParticleHPVector.hh.

    {
      G4int i;
      for (i = min; i < nEntries; i++) {
        if (theData[i].GetX() > e) break;
      }
      G4int low = i - 1;
      G4int high = i;
      if (i == 0) {
        low = 0;
        high = 1;
      }
      else if (i == nEntries) {
        low = nEntries - 2;
        high = nEntries - 1;
      }
      G4double y;
      if (e < theData[nEntries - 1].GetX()) {
        // Protect against doubled-up x values
        if ((theData[high].GetX() - theData[low].GetX()) / theData[high].GetX() < 0.000001) {
          y = theData[low].GetY();
        }
        else {
          y = theInt.Interpolate(theManager.GetScheme(high), e, theData[low].GetX(),
                                 theData[high].GetX(), theData[low].GetY(), theData[high].GetY());
        }
      }
      else {
        y = theData[nEntries - 1].GetY();
      }
      return y;
    }

GetXsec() [3/3]

G4double G4ParticleHPVector::GetXsec ( G4int i )

inline

Definition at line 119 of file G4ParticleHPVector.hh.

{ return theData[i].GetY(); }

Referenced by GetMaxY(), GetY(), G4ParticleHPChannel::Harmonise(), G4ParticleHPElementData::Harmonise(), Merge(), and Merge().

GetY() [1/3]

G4double G4ParticleHPVector::GetY ( G4double x )

inline

Definition at line 184 of file G4ParticleHPVector.hh.

{ return GetXsec(x); }

Referenced by G4ParticleHPData::DoPhysicsVector(), GetMaxY(), G4ParticleHPPhotonDist::GetPhotons(), Hash(), operator+, G4ParticleHPContAngularPar::Sample(), G4ParticleHPDiscreteTwoBody::Sample(), G4ParticleHPLabAngularEnergy::Sample(), G4ParticleHPLegendreStore::Sample(), Sample(), ThinOut(), and Times().

GetY() [2/3]

G4double G4ParticleHPVector::GetY ( G4int i )

inline

Definition at line 187 of file G4ParticleHPVector.hh.

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

GetY() [3/3]

G4double G4ParticleHPVector::GetY ( G4int i ) const

inline

Definition at line 194 of file G4ParticleHPVector.hh.

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

Hash()

void G4ParticleHPVector::Hash ( )

inline

Definition at line 128 of file G4ParticleHPVector.hh.

    {
      G4int i;
      G4double x, y;
      for (i = 0; i < nEntries; i++) {
        if (0 == (i + 1) % 10) {
          x = GetX(i);
          y = GetY(i);
          theHash.SetData(i, x, y);
        }
      }
    }

Referenced by G4ParticleHPIsoData::FillChannelData(), GetEnergyIndex(), GetMaxY(), GetXsec(), and ReHash().

Init() [1/2]

void G4ParticleHPVector::Init ( std::istream & aDataFile, G4double ux = 1., G4double uy = 1. )

inline

Definition at line 218 of file G4ParticleHPVector.hh.

    {
      G4int total;
      aDataFile >> total;
      delete[] theData;
      theData = new G4ParticleHPDataPoint[total];
      nPoints = total;
      nEntries = 0;
      theManager.Init(aDataFile);
      Init(aDataFile, total, ux, uy);
    }

Init() [2/2]

void G4ParticleHPVector::Init ( std::istream & aDataFile, G4int total, G4double ux = 1., G4double uy = 1. )

inline

Definition at line 204 of file G4ParticleHPVector.hh.

    {
      G4double x, y;
      for (G4int i = 0; i < total; i++) {
        aDataFile >> x >> y;
        x *= ux;
        y *= uy;
        SetData(i, x, y);
        if (0 == nEntries % 10) {
          theHash.SetData(nEntries - 1, x, y);
        }
      }
    }

Referenced by G4ParticleHPArbitaryTab::Init(), and Init().

InitInterpolation()

void G4ParticleHPVector::InitInterpolation ( std::istream & aDataFile )

inline

Definition at line 202 of file G4ParticleHPVector.hh.

{ theManager.Init(aDataFile); }

Integrate()

void G4ParticleHPVector::Integrate ( )

inline

Definition at line 396 of file G4ParticleHPVector.hh.

    {
      G4int i;
      if (nEntries == 1) {
        totalIntegral = 0;
        return;
      }
      G4double sum = 0;
      for (i = 1; i < GetVectorLength(); i++) {
        if (std::abs((theData[i].GetX() - theData[i - 1].GetX()) / theData[i].GetX()) > 0.0000001) {
          G4double x1 = theData[i - 1].GetX();
          G4double x2 = theData[i].GetX();
          G4double y1 = theData[i - 1].GetY();
          G4double y2 = theData[i].GetY();
          G4InterpolationScheme aScheme = theManager.GetScheme(i);
          if (aScheme == LINLIN || aScheme == CLINLIN || aScheme == ULINLIN) {
            sum += 0.5 * (y2 + y1) * (x2 - x1);
          }
          else if (aScheme == LINLOG || aScheme == CLINLOG || aScheme == ULINLOG) {
            G4double a = y1;
            G4double b = (y2 - y1) / (G4Log(x2) - G4Log(x1));
            sum += (a - b) * (x2 - x1) + b * (x2 * G4Log(x2) - x1 * G4Log(x1));
          }
          else if (aScheme == LOGLIN || aScheme == CLOGLIN || aScheme == ULOGLIN) {
            G4double a = G4Log(y1);
            G4double b = (G4Log(y2) - G4Log(y1)) / (x2 - x1);
            sum += (G4Exp(a) / b) * (G4Exp(b * x2) - G4Exp(b * x1));
          }
          else if (aScheme == HISTO || aScheme == CHISTO || aScheme == UHISTO) {
            sum += y1 * (x2 - x1);
          }
          else if (aScheme == LOGLOG || aScheme == CLOGLOG || aScheme == ULOGLOG) {
            G4double a = G4Log(y1);
            G4double b = (G4Log(y2) - G4Log(y1)) / (G4Log(x2) - G4Log(x1));
            sum +=
              (G4Exp(a) / (b + 1))
              * (G4Pow::GetInstance()->powA(x2, b + 1) - G4Pow::GetInstance()->powA(x1, b + 1));
          }
          else {
            throw G4HadronicException(
              __FILE__, __LINE__, "Unknown interpolation scheme in G4ParticleHPVector::Integrate");
          }
        }
      }
      totalIntegral = sum;
    }

Referenced by GetIntegral().

IntegrateAndNormalise()

void G4ParticleHPVector::IntegrateAndNormalise ( )

inline

Definition at line 344 of file G4ParticleHPVector.hh.

    {
      G4int i;
      if (theIntegral != nullptr) return;
      theIntegral = new G4double[nEntries];
      if (nEntries == 1) {
        theIntegral[0] = 1;
        return;
      }
      theIntegral[0] = 0;
      G4double sum = 0;
      G4double x1 = 0;
      G4double x0 = 0;
      for (i = 1; i < GetVectorLength(); i++) {
        x1 = theData[i].GetX();
        x0 = theData[i - 1].GetX();
        if (std::abs(x1 - x0) > std::abs(x1 * 0.0000001)) {
          //********************************************************************
          // EMendoza -> the interpolation scheme is not always lin-lin
          /*
                sum+= 0.5*(theData[i].GetY()+theData[i-1].GetY())*
                          (x1-x0);
          */
          //********************************************************************
          G4InterpolationScheme aScheme = theManager.GetScheme(i);
          G4double y0 = theData[i - 1].GetY();
          G4double y1 = theData[i].GetY();
          G4double integ = theInt.GetBinIntegral(aScheme, x0, x1, y0, y1);
#if defined WIN32 - VC
          if (!_finite(integ)) {
            integ = 0;
          }
#elif defined __IBMCPP__
          if (isinf(integ) || isnan(integ)) {
            integ = 0;
          }
#else
          if (std::isinf(integ) || std::isnan(integ)) {
            integ = 0;
          }
#endif
          sum += integ;
          //********************************************************************
        }
        theIntegral[i] = sum;
      }
      G4double total = theIntegral[GetVectorLength() - 1];
      for (i = 1; i < GetVectorLength(); i++) {
        theIntegral[i] /= total;
      }
    }

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

Merge() [1/2]

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

Definition at line 234 of file G4ParticleHPVector.cc.

{
  // interpolate between labels according to aScheme, cut at aValue,
  // continue in unknown areas by substraction of the last difference.
 
  CleanUp();
  G4int s_tmp = 0, n = 0, m_tmp = 0;
  G4ParticleHPVector* tmp;
  G4int a = s_tmp, p = n, t;
  while (a < active->GetVectorLength())  // Loop checking, 11.05.2015, T. Koi
  {
    if (active->GetEnergy(a) <= passive->GetEnergy(p)) {
      G4double xa = active->GetEnergy(a);
      G4double yy = theInt.Interpolate(aScheme, aValue, active->GetLabel(), passive->GetLabel(),
                                       active->GetXsec(a), passive->GetXsec(xa));
      SetData(m_tmp, xa, yy);
      theManager.AppendScheme(m_tmp, active->GetScheme(a));
      m_tmp++;
      a++;
      G4double xp = passive->GetEnergy(p);
      // if( std::abs(std::abs(xp-xa)/xa)<0.0000001&&a<active->GetVectorLength() )
      if (xa != 0 && std::abs(std::abs(xp - xa) / xa) < 0.0000001 && a < active->GetVectorLength())
      {
        p++;
        tmp = active;
        t = a;
        active = passive;
        a = p;
        passive = tmp;
        p = t;
      }
    }
    else {
      tmp = active;
      t = a;
      active = passive;
      a = p;
      passive = tmp;
      p = t;
    }
  }
 
  G4double deltaX = passive->GetXsec(GetEnergy(m_tmp - 1)) - GetXsec(m_tmp - 1);
  while (p != passive->GetVectorLength()
         && passive->GetEnergy(p) <= aValue)  // Loop checking, 11.05.2015, T. Koi
  {
    G4double anX;
    anX = passive->GetXsec(p) - deltaX;
    if (anX > 0) {
      // if(std::abs(GetEnergy(m-1)-passive->GetEnergy(p))/passive->GetEnergy(p)>0.0000001)
      if (passive->GetEnergy(p) == 0
          || std::abs(GetEnergy(m_tmp - 1) - passive->GetEnergy(p)) / passive->GetEnergy(p)
               > 0.0000001)
      {
        SetData(m_tmp, passive->GetEnergy(p), anX);
        theManager.AppendScheme(m_tmp++, passive->GetScheme(p));
      }
    }
    p++;
  }
  // Rebuild the Hash;
  if (theHash.Prepared()) {
    ReHash();
  }
}

Merge() [2/2]

void G4ParticleHPVector::Merge ( G4ParticleHPVector * active, G4ParticleHPVector * passive )

inline

Definition at line 248 of file G4ParticleHPVector.hh.

    {
      CleanUp();
      G4int s_tmp = 0, n = 0, m_tmp = 0;
      G4ParticleHPVector* tmp;
      G4int a = s_tmp, p = n, t;
      while (a < active->GetVectorLength()
             && p < passive->GetVectorLength())  // Loop checking, 11.05.2015, T. Koi
      {
        if (active->GetEnergy(a) <= passive->GetEnergy(p)) {
          G4double xa = active->GetEnergy(a);
          G4double yy = active->GetXsec(a);
          SetData(m_tmp, xa, yy);
          theManager.AppendScheme(m_tmp, active->GetScheme(a));
          m_tmp++;
          a++;
          G4double xp = passive->GetEnergy(p);
 
          // 080409 TKDB
          // if( std::abs(std::abs(xp-xa)/xa)<0.001 ) p++;
          if (!(xa == 0) && std::abs(std::abs(xp - xa) / xa) < 0.001) p++;
        }
        else {
          tmp = active;
          t = a;
          active = passive;
          a = p;
          passive = tmp;
          p = t;
        }
      }
      while (a != active->GetVectorLength())  // Loop checking, 11.05.2015, T. Koi
      {
        SetData(m_tmp, active->GetEnergy(a), active->GetXsec(a));
        theManager.AppendScheme(m_tmp++, active->GetScheme(a));
        a++;
      }
      while (p != passive->GetVectorLength())  // Loop checking, 11.05.2015, T. Koi
      {
        if (std::abs(GetEnergy(m_tmp - 1) - passive->GetEnergy(p)) / passive->GetEnergy(p) > 0.001)
        // if(std::abs(GetEnergy(m)-passive->GetEnergy(p))/passive->GetEnergy(p)>0.001)
        {
          SetData(m_tmp, passive->GetEnergy(p), passive->GetXsec(p));
          theManager.AppendScheme(m_tmp++, active->GetScheme(p));
        }
        p++;
      }
    }

Referenced by G4ParticleHPDiscreteTwoBody::Sample(), and G4ParticleHPLabAngularEnergy::Sample().

operator=()

G4ParticleHPVector & G4ParticleHPVector::operator=

(

const G4ParticleHPVector &

right

)

Definition at line 121 of file G4ParticleHPVector.cc.

{
  if (&right == this) return *this;
 
  G4int i;
 
  totalIntegral = right.totalIntegral;
  if (right.theIntegral != nullptr) theIntegral = new G4double[right.nEntries];
  for (i = 0; i < right.nEntries; i++) {
    SetPoint(i, right.GetPoint(i));  // copy theData
    if (right.theIntegral != nullptr) theIntegral[i] = right.theIntegral[i];
  }
  theManager = right.theManager;
  label = right.label;
 
  Verbose = right.Verbose;
  the15percentBorderCash = right.the15percentBorderCash;
  the50percentBorderCash = right.the50percentBorderCash;
  theHash = right.theHash;
  return *this;
}

ReHash()

void G4ParticleHPVector::ReHash ( )

inline

Definition at line 141 of file G4ParticleHPVector.hh.

    {
      theHash.Clear();
      Hash();
    }

Referenced by Merge(), and ThinOut().

Sample()

G4double G4ParticleHPVector::Sample

(

)

Definition at line 375 of file G4ParticleHPVector.cc.

{
  G4double result = 0.;
  G4int j;
  for (j = 0; j < GetVectorLength(); j++) {
    if (GetY(j) < 0) SetY(j, 0);
  }
 
  if (!theBuffered.empty() && G4UniformRand() < 0.5) {
    result = theBuffered[0];
    theBuffered.erase(theBuffered.begin());
    if (result < GetX(GetVectorLength() - 1)) return result;
  }
  if (GetVectorLength() == 1) {
    result = theData[0].GetX();
  }
  else {
    if (theIntegral == nullptr) {
      IntegrateAndNormalise();
    }
    G4int icounter = 0;
    G4int icounter_max = 1024;
    do {
      icounter++;
      if (icounter > icounter_max) {
        G4cout << "Loop-counter exceeded the threshold value at " << __LINE__ << "th line of "
               << __FILE__ << "." << G4endl;
        break;
      }
      // 080808
      /*
              G4double rand;
              G4double value, test, baseline;
              baseline = theData[GetVectorLength()-1].GetX()-theData[0].GetX();
              do
              {
                value = baseline*G4UniformRand();
                value += theData[0].GetX();
                test = GetY(value)/maxValue;
                rand = G4UniformRand();
              }
              //while(test<rand);
              while( test < rand && test > 0 );
              result = value;
      */
      G4double rand;
      G4double value = 0., test;
      G4int jcounter = 0;
      G4int jcounter_max = 1024;
      do {
        jcounter++;
        if (jcounter > jcounter_max) {
          G4cout << "Loop-counter exceeded the threshold value at " << __LINE__ << "th line of "
                 << __FILE__ << "." << G4endl;
          break;
        }
        rand = G4UniformRand();
        G4int ibin = -1;
        for (G4int i = 0; i < GetVectorLength(); i++) {
          if (rand < theIntegral[i]) {
            ibin = i;
            break;
          }
        }
        if (ibin < 0) G4cout << "TKDB 080807 " << rand << G4endl;
        // result
        rand = G4UniformRand();
        G4double x1, x2;
        if (ibin == 0) {
          x1 = theData[ibin].GetX();
          value = x1;
          break;
        }
        x1 = theData[ibin - 1].GetX();
 
        x2 = theData[ibin].GetX();
        value = rand * (x2 - x1) + x1;
        //***********************************************************************
        /*
              test = GetY ( value ) / std::max ( GetY( ibin-1 ) , GetY ( ibin ) );
        */
        //***********************************************************************
        // EMendoza - Always linear interpolation:
        G4double y1 = theData[ibin - 1].GetY();
        G4double y2 = theData[ibin].GetY();
        G4double mval = (y2 - y1) / (x2 - x1);
        G4double bval = y1 - mval * x1;
        test = (mval * value + bval) / std::max(GetY(ibin - 1), GetY(ibin));
        //***********************************************************************
      } while (G4UniformRand() > test);  // Loop checking, 11.05.2015, T. Koi
      result = value;
      // 080807
    } while (IsBlocked(result));  // Loop checking, 11.05.2015, T. Koi
  }
  return result;
}

Referenced by G4ParticleHPPhotonDist::GetPhotons(), G4ParticleHPContAngularPar::Sample(), G4ParticleHPDiscreteTwoBody::Sample(), G4ParticleHPLabAngularEnergy::Sample(), and G4ParticleHPPartial::Sample().

SampleLin()

G4double G4ParticleHPVector::SampleLin ( )

inline

Definition at line 300 of file G4ParticleHPVector.hh.

    {
      G4double result;
      if (theIntegral == nullptr) IntegrateAndNormalise();
      if (GetVectorLength() == 1) {
        result = theData[0].GetX();
      }
      else {
        G4int i;
        G4double rand = G4UniformRand();
 
        // this was replaced
        //      for(i=1;i<GetVectorLength();i++)
        //      {
        //  if(rand<theIntegral[i]/theIntegral[GetVectorLength()-1]) break;
        //      }
 
        // by this (begin)
        for (i = GetVectorLength() - 1; i >= 0; i--) {
          if (rand > theIntegral[i] / theIntegral[GetVectorLength() - 1]) break;
        }
        if (i != GetVectorLength() - 1) i++;
        // until this (end)
 
        G4double x1, x2, y1, y2;
        y1 = theData[i - 1].GetX();
        x1 = theIntegral[i - 1];
        y2 = theData[i].GetX();
        x2 = theIntegral[i];
        if (std::abs((y2 - y1) / y2)
            < 0.0000001)  // not really necessary, since the case is excluded by construction
        {
          y1 = theData[i - 2].GetX();
          x1 = theIntegral[i - 2];
        }
        result = theLin.Lin(rand, x1, x2, y1, y2);
      }
      return result;
    }

SetData()

void G4ParticleHPVector::SetData ( G4int i, G4double x, G4double y )

inline

Definition at line 89 of file G4ParticleHPVector.hh.

    {
      //    G4cout <<"G4ParticleHPVector::SetData called"<<nPoints<<" "<<nEntries<<G4endl;
      Check(i);
      if (y > maxValue) maxValue = y;
      theData[i].SetData(x, y);
    }

Referenced by Init(), Merge(), Merge(), G4ParticleHPDiscreteTwoBody::Sample(), G4ParticleHPLabAngularEnergy::Sample(), G4ParticleHPLegendreStore::Sample(), and SetPoint().

SetEnergy()

void G4ParticleHPVector::SetEnergy ( G4int i, G4double e )

inline

Definition at line 101 of file G4ParticleHPVector.hh.

    {
      Check(i);
      theData[i].SetX(e);
    }

SetInterpolationManager() [1/2]

void G4ParticleHPVector::SetInterpolationManager ( const G4InterpolationManager & aManager )

inline

Definition at line 449 of file G4ParticleHPVector.hh.

    {
      theManager = aManager;
    }

Referenced by G4ParticleHPContAngularPar::Sample(), G4ParticleHPDiscreteTwoBody::Sample(), G4ParticleHPLabAngularEnergy::Sample(), and G4ParticleHPPartial::Sample().

SetInterpolationManager() [2/2]

void G4ParticleHPVector::SetInterpolationManager ( G4InterpolationManager & aMan )

inline

Definition at line 456 of file G4ParticleHPVector.hh.

{ theManager = aMan; }

SetLabel()

void G4ParticleHPVector::SetLabel ( G4double aLabel )

inline

Definition at line 232 of file G4ParticleHPVector.hh.

{ label = aLabel; }

SetPoint()

void G4ParticleHPVector::SetPoint ( G4int i, const G4ParticleHPDataPoint & it )

inline

Definition at line 82 of file G4ParticleHPVector.hh.

    {
      G4double x = it.GetX();
      G4double y = it.GetY();
      SetData(i, x, y);
    }

Referenced by operator=().

SetScheme()

void G4ParticleHPVector::SetScheme ( G4int aPoint, const G4InterpolationScheme & aScheme )

inline

Definition at line 458 of file G4ParticleHPVector.hh.

    {
      theManager.AppendScheme(aPoint, aScheme);
    }

Referenced by G4ParticleHPPartial::Sample().

SetVerbose()

void G4ParticleHPVector::SetVerbose ( G4int ff )

inline

Definition at line 69 of file G4ParticleHPVector.hh.

{ Verbose = ff; }

SetX()

void G4ParticleHPVector::SetX ( G4int i, G4double e )

inline

Definition at line 96 of file G4ParticleHPVector.hh.

    {
      Check(i);
      theData[i].SetX(e);
    }

Referenced by G4ParticleHPContAngularPar::Sample(), G4ParticleHPDiscreteTwoBody::Sample(), G4ParticleHPLabAngularEnergy::Sample(), and G4ParticleHPPartial::Sample().

SetXsec()

void G4ParticleHPVector::SetXsec ( G4int i, G4double x )

inline

Definition at line 112 of file G4ParticleHPVector.hh.

    {
      Check(i);
      if (x > maxValue) maxValue = x;
      theData[i].SetY(x);
    }

SetY()

void G4ParticleHPVector::SetY ( G4int i, G4double x )

inline

Definition at line 106 of file G4ParticleHPVector.hh.

    {
      Check(i);
      if (x > maxValue) maxValue = x;
      theData[i].SetY(x);
    }

Referenced by G4ParticleHPContAngularPar::Sample(), G4ParticleHPDiscreteTwoBody::Sample(), G4ParticleHPLabAngularEnergy::Sample(), G4ParticleHPPartial::Sample(), and Sample().

ThinOut()

void G4ParticleHPVector::ThinOut

(

G4double

precision

)

Definition at line 301 of file G4ParticleHPVector.cc.

{
  // anything in there?
  if (GetVectorLength() == 0) return;
  // make the new vector
  auto aBuff = new G4ParticleHPDataPoint[nPoints];
  G4double x, x1, x2, y, y1, y2;
  G4int count = 0, current = 2, start = 1;
 
  // First element always goes and is never tested.
  aBuff[0] = theData[0];
 
  // Find the rest
  while (current < GetVectorLength())  // Loop checking, 11.05.2015, T. Koi
  {
    x1 = aBuff[count].GetX();
    y1 = aBuff[count].GetY();
    x2 = theData[current].GetX();
    y2 = theData[current].GetY();
 
    if (x1 - x2 == 0) {
      // Following block added for avoiding div 0 error on Release + G4FPE_DEBUG
      for (G4int j = start; j < current; j++) {
        y = (y2 + y1) / 2.;
        if (std::abs(y - theData[j].GetY()) > precision * y) {
          aBuff[++count] = theData[current - 1];  // for this one, everything was fine
          start = current;  // the next candidate
          break;
        }
      }
    }
    else {
      for (G4int j = start; j < current; j++) {
        x = theData[j].GetX();
        if (x1 - x2 == 0)
          y = (y2 + y1) / 2.;
        else
          y = theInt.Lin(x, x1, x2, y1, y2);
        if (std::abs(y - theData[j].GetY()) > precision * y) {
          aBuff[++count] = theData[current - 1];  // for this one, everything was fine
          start = current;  // the next candidate
          break;
        }
      }
    }
    current++;
  }
  // The last one also always goes, and is never tested.
  aBuff[++count] = theData[GetVectorLength() - 1];
  delete[] theData;
  theData = aBuff;
  nEntries = count + 1;
 
  // Rebuild the Hash;
  if (theHash.Prepared()) {
    ReHash();
  }
}

Times()

void G4ParticleHPVector::Times ( G4double factor )

inline

Definition at line 71 of file G4ParticleHPVector.hh.

    {
      G4int i;
      for (i = 0; i < nEntries; i++) {
        theData[i].SetY(theData[i].GetY() * factor);
      }
      if (theIntegral != nullptr) {
        theIntegral[i] *= factor;
      }
    }

Friends And Related Symbol Documentation

operator+

G4ParticleHPVector & operator+ ( G4ParticleHPVector & left, G4ParticleHPVector & right )

friend

Definition at line 41 of file G4ParticleHPVector.cc.

{
  auto result = new G4ParticleHPVector;
  G4int j = 0;
  G4double x;
  G4double y;
  G4int running = 0;
  for (G4int i = 0; i < left.GetVectorLength(); i++) {
    while (j < right.GetVectorLength())  // Loop checking, 11.05.2015, T. Koi
    {
      if (right.GetX(j) < left.GetX(i) * 1.001) {
        x = right.GetX(j);
        y = right.GetY(j) + left.GetY(x);
        result->SetData(running++, x, y);
        j++;
      }
      // else if(std::abs((right.GetX(j)-left.GetX(i))/(left.GetX(i)+right.GetX(j)))>0.001)
      else if (left.GetX(i) + right.GetX(j) == 0
               || std::abs((right.GetX(j) - left.GetX(i)) / (left.GetX(i) + right.GetX(j))) > 0.001)
      {
        x = left.GetX(i);
        y = left.GetY(i) + right.GetY(x);
        result->SetData(running++, x, y);
        break;
      }
      else {
        break;
      }
    }
    if (j == right.GetVectorLength()) {
      x = left.GetX(i);
      y = left.GetY(i) + right.GetY(x);
      result->SetData(running++, x, y);
    }
  }
  result->ThinOut(0.02);
  return *result;
}

---

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

• G4ParticleHPVector.hh
• G4ParticleHPVector.cc