Back to home page

EIC code displayed by LXR
 
 

    

File indexing completed on 2025-01-18 09:58:53 

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     G4double GetXsec(G4double e, G4int min)
0149     {
0150       G4int i;
0151       for (i = min; i < nEntries; i++) {
0152         if (theData[i].GetX() > e) break;
0153       }
0154       G4int low = i - 1;
0155       G4int high = i;
0156       if (i == 0) {
0157         low = 0;
0158         high = 1;
0159       }
0160       else if (i == nEntries) {
0161         low = nEntries - 2;
0162         high = nEntries - 1;
0163       }
0164       G4double y;
0165       if (e < theData[nEntries - 1].GetX()) {
0166         // Protect against doubled-up x values
0167         if ((theData[high].GetX() - theData[low].GetX()) / theData[high].GetX() < 0.000001) {
0168           y = theData[low].GetY();
0169         }
0170         else {
0171           y = theInt.Interpolate(theManager.GetScheme(high), e, theData[low].GetX(),
0172                                  theData[high].GetX(), theData[low].GetY(), theData[high].GetY());
0173         }
0174       }
0175       else {
0176         y = theData[nEntries - 1].GetY();
0177       }
0178       return y;
0179     }
0180 
0181     inline G4double GetY(G4double x) { return GetXsec(x); }
0182     inline G4int GetVectorLength() const { return nEntries; }
0183 
0184     inline G4double GetY(G4int i)
0185     {
0186       if (i < 0) i = 0;
0187       if (i >= GetVectorLength()) i = GetVectorLength() - 1;
0188       return theData[i].GetY();
0189     }
0190 
0191     inline G4double GetY(G4int i) const
0192     {
0193       if (i < 0) i = 0;
0194       if (i >= GetVectorLength()) i = GetVectorLength() - 1;
0195       return theData[i].GetY();
0196     }
0197     void Dump();
0198 
0199     inline void InitInterpolation(std::istream& aDataFile) { theManager.Init(aDataFile); }
0200 
0201     void Init(std::istream& aDataFile, G4int total, G4double ux = 1., G4double uy = 1.)
0202     {
0203       G4double x, y;
0204       for (G4int i = 0; i < total; i++) {
0205         aDataFile >> x >> y;
0206         x *= ux;
0207         y *= uy;
0208         SetData(i, x, y);
0209         if (0 == nEntries % 10) {
0210           theHash.SetData(nEntries - 1, x, y);
0211         }
0212       }
0213     }
0214 
0215     void Init(std::istream& aDataFile, G4double ux = 1., G4double uy = 1.)
0216     {
0217       G4int total;
0218       aDataFile >> total;
0219       delete[] theData;
0220       theData = new G4ParticleHPDataPoint[total];
0221       nPoints = total;
0222       nEntries = 0;
0223       theManager.Init(aDataFile);
0224       Init(aDataFile, total, ux, uy);
0225     }
0226 
0227     void ThinOut(G4double precision);
0228 
0229     inline void SetLabel(G4double aLabel) { label = aLabel; }
0230 
0231     inline G4double GetLabel() { return label; }
0232 
0233     inline void CleanUp()
0234     {
0235       nEntries = 0;
0236       theManager.CleanUp();
0237       maxValue = -DBL_MAX;
0238       theHash.Clear();
0239       // 080811 TK DB
0240       delete[] theIntegral;
0241       theIntegral = nullptr;
0242     }
0243 
0244     // merges the vectors active and passive into *this
0245     inline void Merge(G4ParticleHPVector* active, G4ParticleHPVector* passive)
0246     {
0247       CleanUp();
0248       G4int s_tmp = 0, n = 0, m_tmp = 0;
0249       G4ParticleHPVector* tmp;
0250       G4int a = s_tmp, p = n, t;
0251       while (a < active->GetVectorLength()
0252              && p < passive->GetVectorLength())  // Loop checking, 11.05.2015, T. Koi
0253       {
0254         if (active->GetEnergy(a) <= passive->GetEnergy(p)) {
0255           G4double xa = active->GetEnergy(a);
0256           G4double yy = active->GetXsec(a);
0257           SetData(m_tmp, xa, yy);
0258           theManager.AppendScheme(m_tmp, active->GetScheme(a));
0259           m_tmp++;
0260           a++;
0261           G4double xp = passive->GetEnergy(p);
0262 
0263           // 080409 TKDB
0264           // if( std::abs(std::abs(xp-xa)/xa)<0.001 ) p++;
0265           if (!(xa == 0) && std::abs(std::abs(xp - xa) / xa) < 0.001) p++;
0266         }
0267         else {
0268           tmp = active;
0269           t = a;
0270           active = passive;
0271           a = p;
0272           passive = tmp;
0273           p = t;
0274         }
0275       }
0276       while (a != active->GetVectorLength())  // Loop checking, 11.05.2015, T. Koi
0277       {
0278         SetData(m_tmp, active->GetEnergy(a), active->GetXsec(a));
0279         theManager.AppendScheme(m_tmp++, active->GetScheme(a));
0280         a++;
0281       }
0282       while (p != passive->GetVectorLength())  // Loop checking, 11.05.2015, T. Koi
0283       {
0284         if (std::abs(GetEnergy(m_tmp - 1) - passive->GetEnergy(p)) / passive->GetEnergy(p) > 0.001)
0285         // if(std::abs(GetEnergy(m)-passive->GetEnergy(p))/passive->GetEnergy(p)>0.001)
0286         {
0287           SetData(m_tmp, passive->GetEnergy(p), passive->GetXsec(p));
0288           theManager.AppendScheme(m_tmp++, active->GetScheme(p));
0289         }
0290         p++;
0291       }
0292     }
0293 
0294     void Merge(G4InterpolationScheme aScheme, G4double aValue, G4ParticleHPVector* active,
0295                G4ParticleHPVector* passive);
0296 
0297     G4double SampleLin()  // Samples X according to distribution Y, linear int
0298     {
0299       G4double result;
0300       if (theIntegral == nullptr) IntegrateAndNormalise();
0301       if (GetVectorLength() == 1) {
0302         result = theData[0].GetX();
0303       }
0304       else {
0305         G4int i;
0306         G4double rand = G4UniformRand();
0307 
0308         // this was replaced
0309         //      for(i=1;i<GetVectorLength();i++)
0310         //      {
0311         //  if(rand<theIntegral[i]/theIntegral[GetVectorLength()-1]) break;
0312         //      }
0313 
0314         // by this (begin)
0315         for (i = GetVectorLength() - 1; i >= 0; i--) {
0316           if (rand > theIntegral[i] / theIntegral[GetVectorLength() - 1]) break;
0317         }
0318         if (i != GetVectorLength() - 1) i++;
0319         // until this (end)
0320 
0321         G4double x1, x2, y1, y2;
0322         y1 = theData[i - 1].GetX();
0323         x1 = theIntegral[i - 1];
0324         y2 = theData[i].GetX();
0325         x2 = theIntegral[i];
0326         if (std::abs((y2 - y1) / y2)
0327             < 0.0000001)  // not really necessary, since the case is excluded by construction
0328         {
0329           y1 = theData[i - 2].GetX();
0330           x1 = theIntegral[i - 2];
0331         }
0332         result = theLin.Lin(rand, x1, x2, y1, y2);
0333       }
0334       return result;
0335     }
0336 
0337     G4double Sample();  // Samples X according to distribution Y
0338 
0339     G4double* Debug() { return theIntegral; }
0340 
0341     inline void IntegrateAndNormalise()
0342     {
0343       G4int i;
0344       if (theIntegral != nullptr) return;
0345       theIntegral = new G4double[nEntries];
0346       if (nEntries == 1) {
0347         theIntegral[0] = 1;
0348         return;
0349       }
0350       theIntegral[0] = 0;
0351       G4double sum = 0;
0352       G4double x1 = 0;
0353       G4double x0 = 0;
0354       for (i = 1; i < GetVectorLength(); i++) {
0355         x1 = theData[i].GetX();
0356         x0 = theData[i - 1].GetX();
0357         if (std::abs(x1 - x0) > std::abs(x1 * 0.0000001)) {
0358           //********************************************************************
0359           // EMendoza -> the interpolation scheme is not always lin-lin
0360           /*
0361                 sum+= 0.5*(theData[i].GetY()+theData[i-1].GetY())*
0362                           (x1-x0);
0363           */
0364           //********************************************************************
0365           G4InterpolationScheme aScheme = theManager.GetScheme(i);
0366           G4double y0 = theData[i - 1].GetY();
0367           G4double y1 = theData[i].GetY();
0368           G4double integ = theInt.GetBinIntegral(aScheme, x0, x1, y0, y1);
0369 #if defined WIN32 - VC
0370           if (!_finite(integ)) {
0371             integ = 0;
0372           }
0373 #elif defined __IBMCPP__
0374           if (isinf(integ) || isnan(integ)) {
0375             integ = 0;
0376           }
0377 #else
0378           if (std::isinf(integ) || std::isnan(integ)) {
0379             integ = 0;
0380           }
0381 #endif
0382           sum += integ;
0383           //********************************************************************
0384         }
0385         theIntegral[i] = sum;
0386       }
0387       G4double total = theIntegral[GetVectorLength() - 1];
0388       for (i = 1; i < GetVectorLength(); i++) {
0389         theIntegral[i] /= total;
0390       }
0391     }
0392 
0393     inline void Integrate()
0394     {
0395       G4int i;
0396       if (nEntries == 1) {
0397         totalIntegral = 0;
0398         return;
0399       }
0400       G4double sum = 0;
0401       for (i = 1; i < GetVectorLength(); i++) {
0402         if (std::abs((theData[i].GetX() - theData[i - 1].GetX()) / theData[i].GetX()) > 0.0000001) {
0403           G4double x1 = theData[i - 1].GetX();
0404           G4double x2 = theData[i].GetX();
0405           G4double y1 = theData[i - 1].GetY();
0406           G4double y2 = theData[i].GetY();
0407           G4InterpolationScheme aScheme = theManager.GetScheme(i);
0408           if (aScheme == LINLIN || aScheme == CLINLIN || aScheme == ULINLIN) {
0409             sum += 0.5 * (y2 + y1) * (x2 - x1);
0410           }
0411           else if (aScheme == LINLOG || aScheme == CLINLOG || aScheme == ULINLOG) {
0412             G4double a = y1;
0413             G4double b = (y2 - y1) / (G4Log(x2) - G4Log(x1));
0414             sum += (a - b) * (x2 - x1) + b * (x2 * G4Log(x2) - x1 * G4Log(x1));
0415           }
0416           else if (aScheme == LOGLIN || aScheme == CLOGLIN || aScheme == ULOGLIN) {
0417             G4double a = G4Log(y1);
0418             G4double b = (G4Log(y2) - G4Log(y1)) / (x2 - x1);
0419             sum += (G4Exp(a) / b) * (G4Exp(b * x2) - G4Exp(b * x1));
0420           }
0421           else if (aScheme == HISTO || aScheme == CHISTO || aScheme == UHISTO) {
0422             sum += y1 * (x2 - x1);
0423           }
0424           else if (aScheme == LOGLOG || aScheme == CLOGLOG || aScheme == ULOGLOG) {
0425             G4double a = G4Log(y1);
0426             G4double b = (G4Log(y2) - G4Log(y1)) / (G4Log(x2) - G4Log(x1));
0427             sum +=
0428               (G4Exp(a) / (b + 1))
0429               * (G4Pow::GetInstance()->powA(x2, b + 1) - G4Pow::GetInstance()->powA(x1, b + 1));
0430           }
0431           else {
0432             throw G4HadronicException(
0433               __FILE__, __LINE__, "Unknown interpolation scheme in G4ParticleHPVector::Integrate");
0434           }
0435         }
0436       }
0437       totalIntegral = sum;
0438     }
0439 
0440     inline G4double GetIntegral()  // linear interpolation; use with care
0441     {
0442       if (totalIntegral < -0.5) Integrate();
0443       return totalIntegral;
0444     }
0445 
0446     inline void SetInterpolationManager(const G4InterpolationManager& aManager)
0447     {
0448       theManager = aManager;
0449     }
0450 
0451     inline const G4InterpolationManager& GetInterpolationManager() const { return theManager; }
0452 
0453     inline void SetInterpolationManager(G4InterpolationManager& aMan) { theManager = aMan; }
0454 
0455     inline void SetScheme(G4int aPoint, const G4InterpolationScheme& aScheme)
0456     {
0457       theManager.AppendScheme(aPoint, aScheme);
0458     }
0459 
0460     inline G4InterpolationScheme GetScheme(G4int anIndex) { return theManager.GetScheme(anIndex); }
0461 
0462     G4double GetMeanX()
0463     {
0464       G4double result;
0465       G4double running = 0;
0466       G4double weighted = 0;
0467       for (G4int i = 1; i < nEntries; i++) {
0468         running +=
0469           theInt.GetBinIntegral(theManager.GetScheme(i - 1), theData[i - 1].GetX(),
0470                                 theData[i].GetX(), theData[i - 1].GetY(), theData[i].GetY());
0471         weighted += theInt.GetWeightedBinIntegral(theManager.GetScheme(i - 1),
0472                                                   theData[i - 1].GetX(), theData[i].GetX(),
0473                                                   theData[i - 1].GetY(), theData[i].GetY());
0474       }
0475       result = weighted / running;
0476       return result;
0477     }
0478 
0479     // Finds maximum cross section between two values of kinetic energy
0480     G4double GetMaxY(G4double emin, G4double emax);
0481 
0482     /*
0483       void Block(G4double aX)
0484       {
0485         theBlocked.push_back(aX);
0486       }
0487 
0488       void Buffer(G4double aX)
0489       {
0490         theBuffered.push_back(aX);
0491       }
0492     */
0493 
0494     std::vector<G4double> GetBlocked() { return theBlocked; }
0495     std::vector<G4double> GetBuffered() { return theBuffered; }
0496 
0497     //  void SetBlocked(const std::vector<G4double> &aBlocked) {theBlocked = aBlocked;}
0498     //  void SetBuffered(const std::vector<G4double> &aBuffer) {theBuffered = aBuffer;}
0499 
0500     G4double Get15percentBorder();
0501     G4double Get50percentBorder();
0502 
0503   private:
0504     void Check(G4int i);
0505 
0506     G4bool IsBlocked(G4double aX);
0507 
0508   private:
0509     G4ParticleHPInterpolator theLin;
0510 
0511   private:
0512     G4double totalIntegral;
0513 
0514     G4ParticleHPDataPoint* theData;  // the data
0515     G4InterpolationManager theManager;  // knows how to interpolate the data.
0516     G4double* theIntegral;
0517     G4int nEntries;
0518     G4int nPoints;
0519     G4double label;
0520 
0521     G4ParticleHPInterpolator theInt;
0522     G4int Verbose;
0523     // debug only
0524     G4int isFreed;
0525 
0526     G4ParticleHPHash theHash;
0527     G4double maxValue;
0528 
0529     std::vector<G4double> theBlocked;
0530     std::vector<G4double> theBuffered;
0531     G4double the15percentBorderCash;
0532     G4double the50percentBorderCash;
0533 };
0534 
0535 #endif
[ Source navigation ] [ Diff markup ] [ Identifier search ] [ general search ]

This page was automatically generated by the 2.3.7 LXR engine.

The LXR team
Valid CSS 2.1! Valid HTML 4.01!