include/Geant4/G4ParticleHPVector.hh

0001 //
0002 // ********************************************************************
0003 // * License and Disclaimer                                           *
0004 // *                                                                  *
0005 // * The  Geant4 software  is  copyright of the Copyright Holders  of *
0006 // * the Geant4 Collaboration.  It is provided  under  the terms  and *
0007 // * conditions of the Geant4 Software License,  included in the file *
0008 // * LICENSE and available at  http://cern.ch/geant4/license .  These *
0009 // * include a list of copyright holders.                             *
0010 // *                                                                  *
0011 // * Neither the authors of this software system, nor their employing *
0012 // * institutes,nor the agencies providing financial support for this *
0013 // * work  make  any representation or  warranty, express or implied, *
0014 // * regarding  this  software system or assume any liability for its *
0015 // * use.  Please see the license in the file  LICENSE  and URL above *
0016 // * for the full disclaimer and the limitation of liability.         *
0017 // *                                                                  *
0018 // * This  code  implementation is the result of  the  scientific and *
0019 // * technical work of the GEANT4 collaboration.                      *
0020 // * By using,  copying,  modifying or  distributing the software (or *
0021 // * any work based  on the software)  you  agree  to acknowledge its *
0022 // * use  in  resulting  scientific  publications,  and indicate your *
0023 // * acceptance of all terms of the Geant4 Software license.          *
0024 // ********************************************************************
0025 //
0026 //
0027 // 070606 fix with Valgrind by T. Koi
0028 // 080409 Fix div0 error with G4FPE by T. Koi
0029 // 080811 Comment out unused method SetBlocked and SetBuffered
0030 //        Add required cleaning up in CleanUp by T. Koi
0031 //
0032 // P. Arce, June-2014 Conversion neutron_hp to particle_hp
0033 //
0034 #ifndef G4ParticleHPVector_h
0035 #define G4ParticleHPVector_h 1
0036
0037 #include "G4Exp.hh"
0038 #include "G4InterpolationManager.hh"
0039 #include "G4Log.hh"
0040 #include "G4ParticleHPDataPoint.hh"
0041 #include "G4ParticleHPHash.hh"
0042 #include "G4ParticleHPInterpolator.hh"
0043 #include "G4PhysicsVector.hh"
0044 #include "G4Pow.hh"
0045 #include "G4ios.hh"
0046 #include "Randomize.hh"
0047
0048 #include <cmath>
0049 #include <fstream>
0050 #include <vector>
0051
0052 #if defined WIN32 - VC
0053 #  include <float.h>
0054 #endif
0055
0056 class G4ParticleHPVector
0057 {
0058     friend G4ParticleHPVector& operator+(G4ParticleHPVector& left, G4ParticleHPVector& right);
0059
0060   public:
0061     G4ParticleHPVector();
0062
0063     G4ParticleHPVector(G4int n);
0064
0065     ~G4ParticleHPVector();
0066
0067     G4ParticleHPVector& operator=(const G4ParticleHPVector& right);
0068
0069     inline void SetVerbose(G4int ff) { Verbose = ff; }
0070
0071     inline void Times(G4double factor)
0072     {
0073       G4int i;
0074       for (i = 0; i < nEntries; i++) {
0075         theData[i].SetY(theData[i].GetY() * factor);
0076       }
0077       if (theIntegral != nullptr) {
0078         theIntegral[i] *= factor;
0079       }
0080     }
0081
0082     inline void SetPoint(G4int i, const G4ParticleHPDataPoint& it)
0083     {
0084       G4double x = it.GetX();
0085       G4double y = it.GetY();
0086       SetData(i, x, y);
0087     }
0088
0089     inline void SetData(G4int i, G4double x, G4double y)
0090     {
0091       //    G4cout <<"G4ParticleHPVector::SetData called"<<nPoints<<" "<<nEntries<<G4endl;
0092       Check(i);
0093       if (y > maxValue) maxValue = y;
0094       theData[i].SetData(x, y);
0095     }
0096     inline void SetX(G4int i, G4double e)
0097     {
0098       Check(i);
0099       theData[i].SetX(e);
0100     }
0101     inline void SetEnergy(G4int i, G4double e)
0102     {
0103       Check(i);
0104       theData[i].SetX(e);
0105     }
0106     inline void SetY(G4int i, G4double x)
0107     {
0108       Check(i);
0109       if (x > maxValue) maxValue = x;
0110       theData[i].SetY(x);
0111     }
0112     inline void SetXsec(G4int i, G4double x)
0113     {
0114       Check(i);
0115       if (x > maxValue) maxValue = x;
0116       theData[i].SetY(x);
0117     }
0118     inline G4double GetEnergy(G4int i) const { return theData[i].GetX(); }
0119     inline G4double GetXsec(G4int i) { return theData[i].GetY(); }
0120     inline G4double GetX(G4int i) const
0121     {
0122       if (i < 0) i = 0;
0123       if (i >= GetVectorLength()) i = GetVectorLength() - 1;
0124       return theData[i].GetX();
0125     }
0126     inline const G4ParticleHPDataPoint& GetPoint(G4int i) const { return theData[i]; }
0127
0128     void Hash()
0129     {
0130       G4int i;
0131       G4double x, y;
0132       for (i = 0; i < nEntries; i++) {
0133         if (0 == (i + 1) % 10) {
0134           x = GetX(i);
0135           y = GetY(i);
0136           theHash.SetData(i, x, y);
0137         }
0138       }
0139     }
0140
0141     void ReHash()
0142     {
0143       theHash.Clear();
0144       Hash();
0145     }
0146
0147     G4double GetXsec(G4double e);
0148
0149     G4int GetEnergyIndex(G4double &e);  // method added by M. Zmeskal 02/2024
0150
0151     G4double GetXsec(G4double e, G4int min)
0152     {
0153       G4int i;
0154       for (i = min; i < nEntries; i++) {
0155         if (theData[i].GetX() > e) break;
0156       }
0157       G4int low = i - 1;
0158       G4int high = i;
0159       if (i == 0) {
0160         low = 0;
0161         high = 1;
0162       }
0163       else if (i == nEntries) {
0164         low = nEntries - 2;
0165         high = nEntries - 1;
0166       }
0167       G4double y;
0168       if (e < theData[nEntries - 1].GetX()) {
0169         // Protect against doubled-up x values
0170         if ((theData[high].GetX() - theData[low].GetX()) / theData[high].GetX() < 0.000001) {
0171           y = theData[low].GetY();
0172         }
0173         else {
0174           y = theInt.Interpolate(theManager.GetScheme(high), e, theData[low].GetX(),
0175                                  theData[high].GetX(), theData[low].GetY(), theData[high].GetY());
0176         }
0177       }
0178       else {
0179         y = theData[nEntries - 1].GetY();
0180       }
0181       return y;
0182     }
0183
0184     inline G4double GetY(G4double x) { return GetXsec(x); }
0185     inline G4int GetVectorLength() const { return nEntries; }
0186
0187     inline G4double GetY(G4int i)
0188     {
0189       if (i < 0) i = 0;
0190       if (i >= GetVectorLength()) i = GetVectorLength() - 1;
0191       return theData[i].GetY();
0192     }
0193
0194     inline G4double GetY(G4int i) const
0195     {
0196       if (i < 0) i = 0;
0197       if (i >= GetVectorLength()) i = GetVectorLength() - 1;
0198       return theData[i].GetY();
0199     }
0200     void Dump();
0201
0202     inline void InitInterpolation(std::istream& aDataFile) { theManager.Init(aDataFile); }
0203
0204     void Init(std::istream& aDataFile, G4int total, G4double ux = 1., G4double uy = 1.)
0205     {
0206       G4double x, y;
0207       for (G4int i = 0; i < total; i++) {
0208         aDataFile >> x >> y;
0209         x *= ux;
0210         y *= uy;
0211         SetData(i, x, y);
0212         if (0 == nEntries % 10) {
0213           theHash.SetData(nEntries - 1, x, y);
0214         }
0215       }
0216     }
0217
0218     void Init(std::istream& aDataFile, G4double ux = 1., G4double uy = 1.)
0219     {
0220       G4int total;
0221       aDataFile >> total;
0222       delete[] theData;
0223       theData = new G4ParticleHPDataPoint[total];
0224       nPoints = total;
0225       nEntries = 0;
0226       theManager.Init(aDataFile);
0227       Init(aDataFile, total, ux, uy);
0228     }
0229
0230     void ThinOut(G4double precision);
0231
0232     inline void SetLabel(G4double aLabel) { label = aLabel; }
0233
0234     inline G4double GetLabel() { return label; }
0235
0236     inline void CleanUp()
0237     {
0238       nEntries = 0;
0239       theManager.CleanUp();
0240       maxValue = -DBL_MAX;
0241       theHash.Clear();
0242       // 080811 TK DB
0243       delete[] theIntegral;
0244       theIntegral = nullptr;
0245     }
0246
0247     // merges the vectors active and passive into *this
0248     inline void Merge(G4ParticleHPVector* active, G4ParticleHPVector* passive)
0249     {
0250       CleanUp();
0251       G4int s_tmp = 0, n = 0, m_tmp = 0;
0252       G4ParticleHPVector* tmp;
0253       G4int a = s_tmp, p = n, t;
0254       while (a < active->GetVectorLength()
0255              && p < passive->GetVectorLength())  // Loop checking, 11.05.2015, T. Koi
0256       {
0257         if (active->GetEnergy(a) <= passive->GetEnergy(p)) {
0258           G4double xa = active->GetEnergy(a);
0259           G4double yy = active->GetXsec(a);
0260           SetData(m_tmp, xa, yy);
0261           theManager.AppendScheme(m_tmp, active->GetScheme(a));
0262           m_tmp++;
0263           a++;
0264           G4double xp = passive->GetEnergy(p);
0265
0266           // 080409 TKDB
0267           // if( std::abs(std::abs(xp-xa)/xa)<0.001 ) p++;
0268           if (!(xa == 0) && std::abs(std::abs(xp - xa) / xa) < 0.001) p++;
0269         }
0270         else {
0271           tmp = active;
0272           t = a;
0273           active = passive;
0274           a = p;
0275           passive = tmp;
0276           p = t;
0277         }
0278       }
0279       while (a != active->GetVectorLength())  // Loop checking, 11.05.2015, T. Koi
0280       {
0281         SetData(m_tmp, active->GetEnergy(a), active->GetXsec(a));
0282         theManager.AppendScheme(m_tmp++, active->GetScheme(a));
0283         a++;
0284       }
0285       while (p != passive->GetVectorLength())  // Loop checking, 11.05.2015, T. Koi
0286       {
0287         if (std::abs(GetEnergy(m_tmp - 1) - passive->GetEnergy(p)) / passive->GetEnergy(p) > 0.001)
0288         // if(std::abs(GetEnergy(m)-passive->GetEnergy(p))/passive->GetEnergy(p)>0.001)
0289         {
0290           SetData(m_tmp, passive->GetEnergy(p), passive->GetXsec(p));
0291           theManager.AppendScheme(m_tmp++, active->GetScheme(p));
0292         }
0293         p++;
0294       }
0295     }
0296
0297     void Merge(G4InterpolationScheme aScheme, G4double aValue, G4ParticleHPVector* active,
0298                G4ParticleHPVector* passive);
0299
0300     G4double SampleLin()  // Samples X according to distribution Y, linear int
0301     {
0302       G4double result;
0303       if (theIntegral == nullptr) IntegrateAndNormalise();
0304       if (GetVectorLength() == 1) {
0305         result = theData[0].GetX();
0306       }
0307       else {
0308         G4int i;
0309         G4double rand = G4UniformRand();
0310
0311         // this was replaced
0312         //      for(i=1;i<GetVectorLength();i++)
0313         //      {
0314         //  if(rand<theIntegral[i]/theIntegral[GetVectorLength()-1]) break;
0315         //      }
0316
0317         // by this (begin)
0318         for (i = GetVectorLength() - 1; i >= 0; i--) {
0319           if (rand > theIntegral[i] / theIntegral[GetVectorLength() - 1]) break;
0320         }
0321         if (i != GetVectorLength() - 1) i++;
0322         // until this (end)
0323
0324         G4double x1, x2, y1, y2;
0325         y1 = theData[i - 1].GetX();
0326         x1 = theIntegral[i - 1];
0327         y2 = theData[i].GetX();
0328         x2 = theIntegral[i];
0329         if (std::abs((y2 - y1) / y2)
0330             < 0.0000001)  // not really necessary, since the case is excluded by construction
0331         {
0332           y1 = theData[i - 2].GetX();
0333           x1 = theIntegral[i - 2];
0334         }
0335         result = theLin.Lin(rand, x1, x2, y1, y2);
0336       }
0337       return result;
0338     }
0339
0340     G4double Sample();  // Samples X according to distribution Y
0341
0342     G4double* Debug() { return theIntegral; }
0343
0344     inline void IntegrateAndNormalise()
0345     {
0346       G4int i;
0347       if (theIntegral != nullptr) return;
0348       theIntegral = new G4double[nEntries];
0349       if (nEntries == 1) {
0350         theIntegral[0] = 1;
0351         return;
0352       }
0353       theIntegral[0] = 0;
0354       G4double sum = 0;
0355       G4double x1 = 0;
0356       G4double x0 = 0;
0357       for (i = 1; i < GetVectorLength(); i++) {
0358         x1 = theData[i].GetX();
0359         x0 = theData[i - 1].GetX();
0360         if (std::abs(x1 - x0) > std::abs(x1 * 0.0000001)) {
0361           //********************************************************************
0362           // EMendoza -> the interpolation scheme is not always lin-lin
0363           /*
0364                 sum+= 0.5*(theData[i].GetY()+theData[i-1].GetY())*
0365                           (x1-x0);
0366           */
0367           //********************************************************************
0368           G4InterpolationScheme aScheme = theManager.GetScheme(i);
0369           G4double y0 = theData[i - 1].GetY();
0370           G4double y1 = theData[i].GetY();
0371           G4double integ = theInt.GetBinIntegral(aScheme, x0, x1, y0, y1);
0372 #if defined WIN32 - VC
0373           if (!_finite(integ)) {
0374             integ = 0;
0375           }
0376 #elif defined __IBMCPP__
0377           if (isinf(integ) || isnan(integ)) {
0378             integ = 0;
0379           }
0380 #else
0381           if (std::isinf(integ) || std::isnan(integ)) {
0382             integ = 0;
0383           }
0384 #endif
0385           sum += integ;
0386           //********************************************************************
0387         }
0388         theIntegral[i] = sum;
0389       }
0390       G4double total = theIntegral[GetVectorLength() - 1];
0391       for (i = 1; i < GetVectorLength(); i++) {
0392         theIntegral[i] /= total;
0393       }
0394     }
0395
0396     inline void Integrate()
0397     {
0398       G4int i;
0399       if (nEntries == 1) {
0400         totalIntegral = 0;
0401         return;
0402       }
0403       G4double sum = 0;
0404       for (i = 1; i < GetVectorLength(); i++) {
0405         if (std::abs((theData[i].GetX() - theData[i - 1].GetX()) / theData[i].GetX()) > 0.0000001) {
0406           G4double x1 = theData[i - 1].GetX();
0407           G4double x2 = theData[i].GetX();
0408           G4double y1 = theData[i - 1].GetY();
0409           G4double y2 = theData[i].GetY();
0410           G4InterpolationScheme aScheme = theManager.GetScheme(i);
0411           if (aScheme == LINLIN || aScheme == CLINLIN || aScheme == ULINLIN) {
0412             sum += 0.5 * (y2 + y1) * (x2 - x1);
0413           }
0414           else if (aScheme == LINLOG || aScheme == CLINLOG || aScheme == ULINLOG) {
0415             G4double a = y1;
0416             G4double b = (y2 - y1) / (G4Log(x2) - G4Log(x1));
0417             sum += (a - b) * (x2 - x1) + b * (x2 * G4Log(x2) - x1 * G4Log(x1));
0418           }
0419           else if (aScheme == LOGLIN || aScheme == CLOGLIN || aScheme == ULOGLIN) {
0420             G4double a = G4Log(y1);
0421             G4double b = (G4Log(y2) - G4Log(y1)) / (x2 - x1);
0422             sum += (G4Exp(a) / b) * (G4Exp(b * x2) - G4Exp(b * x1));
0423           }
0424           else if (aScheme == HISTO || aScheme == CHISTO || aScheme == UHISTO) {
0425             sum += y1 * (x2 - x1);
0426           }
0427           else if (aScheme == LOGLOG || aScheme == CLOGLOG || aScheme == ULOGLOG) {
0428             G4double a = G4Log(y1);
0429             G4double b = (G4Log(y2) - G4Log(y1)) / (G4Log(x2) - G4Log(x1));
0430             sum +=
0431               (G4Exp(a) / (b + 1))
0432               * (G4Pow::GetInstance()->powA(x2, b + 1) - G4Pow::GetInstance()->powA(x1, b + 1));
0433           }
0434           else {
0435             throw G4HadronicException(
0436               __FILE__, __LINE__, "Unknown interpolation scheme in G4ParticleHPVector::Integrate");
0437           }
0438         }
0439       }
0440       totalIntegral = sum;
0441     }
0442
0443     inline G4double GetIntegral()  // linear interpolation; use with care
0444     {
0445       if (totalIntegral < -0.5) Integrate();
0446       return totalIntegral;
0447     }
0448
0449     inline void SetInterpolationManager(const G4InterpolationManager& aManager)
0450     {
0451       theManager = aManager;
0452     }
0453
0454     inline const G4InterpolationManager& GetInterpolationManager() const { return theManager; }
0455
0456     inline void SetInterpolationManager(G4InterpolationManager& aMan) { theManager = aMan; }
0457
0458     inline void SetScheme(G4int aPoint, const G4InterpolationScheme& aScheme)
0459     {
0460       theManager.AppendScheme(aPoint, aScheme);
0461     }
0462
0463     inline G4InterpolationScheme GetScheme(G4int anIndex) { return theManager.GetScheme(anIndex); }
0464
0465     G4double GetMeanX()
0466     {
0467       G4double result;
0468       G4double running = 0;
0469       G4double weighted = 0;
0470       for (G4int i = 1; i < nEntries; i++) {
0471         running +=
0472           theInt.GetBinIntegral(theManager.GetScheme(i - 1), theData[i - 1].GetX(),
0473                                 theData[i].GetX(), theData[i - 1].GetY(), theData[i].GetY());
0474         weighted += theInt.GetWeightedBinIntegral(theManager.GetScheme(i - 1),
0475                                                   theData[i - 1].GetX(), theData[i].GetX(),
0476                                                   theData[i - 1].GetY(), theData[i].GetY());
0477       }
0478       result = weighted / running;
0479       return result;
0480     }
0481
0482     // Finds maximum cross section between two values of kinetic energy
0483     G4double GetMaxY(G4double emin, G4double emax);
0484
0485     /*
0486       void Block(G4double aX)
0487       {
0488         theBlocked.push_back(aX);
0489       }
0490
0491       void Buffer(G4double aX)
0492       {
0493         theBuffered.push_back(aX);
0494       }
0495     */
0496
0497     std::vector<G4double> GetBlocked() { return theBlocked; }
0498     std::vector<G4double> GetBuffered() { return theBuffered; }
0499
0500     //  void SetBlocked(const std::vector<G4double> &aBlocked) {theBlocked = aBlocked;}
0501     //  void SetBuffered(const std::vector<G4double> &aBuffer) {theBuffered = aBuffer;}
0502
0503     G4double Get15percentBorder();
0504     G4double Get50percentBorder();
0505
0506   private:
0507     void Check(G4int i);
0508
0509     G4bool IsBlocked(G4double aX);
0510
0511   private:
0512     G4ParticleHPInterpolator theLin;
0513
0514   private:
0515     G4double totalIntegral;
0516
0517     G4ParticleHPDataPoint* theData;  // the data
0518     G4InterpolationManager theManager;  // knows how to interpolate the data.
0519     G4double* theIntegral;
0520     G4int nEntries;
0521     G4int nPoints;
0522     G4double label;
0523
0524     G4ParticleHPInterpolator theInt;
0525     G4int Verbose;
0526     // debug only
0527     G4int isFreed;
0528
0529     G4ParticleHPHash theHash;
0530     G4double maxValue;
0531
0532     std::vector<G4double> theBlocked;
0533     std::vector<G4double> theBuffered;
0534     G4double the15percentBorderCash;
0535     G4double the50percentBorderCash;
0536 };
0537
0538 #endif