g4simulation/g4main/EcoMug.h

0001 /////////////////////////////////////////////////////////////////////////////////////
0002 // EcoMug: Efficient COsmic MUon Generator                                         //
0003 // Copyright (C) 2021 Davide Pagano <davide.pagano@unibs.it>                       //
0004 // EcoMug is based on the following work:                                          //
0005 // D. Pagano, G. Bonomi, A. Donzella, A. Zenoni, G. Zumerle, N. Zurlo,             //
0006 // "EcoMug: an Efficient COsmic MUon Generator for cosmic-ray muons applications", //
0007 // doi:10.1016/j.nima.2021.165732                                                  //
0008 //                                                                                 //
0009 // This program is free software: you can redistribute it and/or modify            //
0010 // it under the terms of the GNU General Public License as published by            //
0011 // the Free Software Foundation, either version 3 of the License, or               //
0012 // (at your option) any later version.                                             //
0013 //                                                                                 //
0014 // This program is distributed in the hope that it will be useful,                 //
0015 // but WITHOUT ANY WARRANTY; without even the implied warranty of                  //
0016 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the                   //
0017 // GNU General Public License for more details.                                    //
0018 //                                                                                 //
0019 // You should have received a copy of the GNU General Public License               //
0020 // along with this program.  If not, see <https://www.gnu.org/licenses/>.          //
0021 /////////////////////////////////////////////////////////////////////////////////////
0022
0023 #ifndef EcoMug_H
0024 #define EcoMug_H
0025
0026 //#include <math.h>
0027 #include <array>
0028 #include <functional>
0029 #include <random>
0030
0031 //! Fast generation of random numbers
0032 //! This class is based on the xoroshiro128+ generator.
0033 //! https://prng.di.unimi.it/
0034 class EMRandom
0035 {
0036  public:
0037   EMRandom()
0038   {
0039     std::random_device rd;
0040     std::mt19937 gen(rd());
0041     std::uniform_int_distribution<uint64_t> dis(0, std::numeric_limits<uint64_t>::max());
0042     s[0] = 12345.;  // dis(gen);
0043     s[1] = 12345;   // dis(gen);
0044   };
0045
0046   void SetSeed(uint64_t seed)
0047   {
0048     s[0] = seed;
0049     s[1] = seed;
0050   };
0051
0052   double GenerateRandomDouble()
0053   {
0054     uint64_t x = next();
0055     return to_double(x);
0056   };
0057
0058   double GenerateRandomDouble(double x1, double x2)
0059   {
0060     return (x2 - x1) * GenerateRandomDouble() + x1;
0061   };
0062
0063   int64_t rotl(const uint64_t x, int k)
0064   {
0065     return (x << k) | (x >> (64 - k));
0066   };
0067
0068   uint64_t next()
0069   {
0070     const uint64_t s0 = s[0];
0071     uint64_t s1 = s[1];
0072     const uint64_t result = s0 + s1;
0073     s1 ^= s0;
0074     s[0] = rotl(s0, 55) ^ s1 ^ (s1 << 14);
0075     s[1] = rotl(s1, 36);
0076     return result;
0077   };
0078
0079   double to_double(uint64_t x)
0080   {
0081     union U
0082     {
0083       uint64_t i;
0084       double d;
0085     };
0086     U u = {UINT64_C(0x3FF) << 52 | x >> 12};
0087     return u.d - 1.0;
0088   };
0089
0090   uint64_t s[2];
0091 };
0092 ///////////////////////////////////////////////////////////////
0093 ///////////////////////////////////////////////////////////////
0094
0095 //! Class for maximization based on "Whale Optimization Algorithm"
0096 //! doi:10.1016/j.advengsoft.2016.01.008
0097 class EMMaximization
0098 {
0099  private:
0100   EMRandom mRandom;
0101   size_t mPopSize;
0102   size_t mNIter;
0103   int mGenMethod;  // 0 = sky, 1 = cylinder, 2 = hspere
0104   double m_a;
0105   double m_a2;
0106   std::vector<std::vector<double> > mRanges;
0107   std::vector<std::vector<double> > mPopulation;
0108   double mBestCost;
0109   std::vector<double> mBestSolution;
0110   std::function<double(double, double)> mFunc;
0111
0112  public:
0113   EMMaximization(const EMRandom& random, int genMethod)
0114     : mRandom(random)
0115     , mPopSize(200)
0116     , mNIter(500)
0117     , mGenMethod(genMethod)
0118     , m_a(0.)
0119     , m_a2(0.)
0120     , mBestCost(-1.)
0121   {
0122     // cppcheck-suppress [useInitializationList]
0123     mFunc = &DefaultJ;
0124   };
0125   ///////////////////////////////////////////////////////////////
0126
0127   static double DefaultJ(double p, double theta)
0128   {
0129     double n = std::max(0.1, 2.856 - 0.655 * log(p));
0130     return 1600 * pow(p, 0.279) * pow(cos(theta), n);
0131   }
0132   ///////////////////////////////////////////////////////////////
0133
0134   void SetParameters(double minP, double maxP, double minTheta, double maxTheta)
0135   {
0136     mRanges.push_back({minP, maxP});
0137     mRanges.push_back({minTheta, maxTheta});
0138   }
0139   ///////////////////////////////////////////////////////////////
0140
0141   void SetParameters(double minP, double maxP, double minTheta, double maxTheta, double minPhi, double maxPhi)
0142   {
0143     mRanges.push_back({minP, maxP});
0144     mRanges.push_back({minTheta, maxTheta});
0145     mRanges.push_back({minPhi, maxPhi});
0146     mRanges.push_back({0, M_PI / 2.});
0147   }
0148   ///////////////////////////////////////////////////////////////
0149
0150   void SetFunction(std::function<double(double, double)> func)
0151   {
0152     mFunc = func;
0153   }
0154   ///////////////////////////////////////////////////////////////
0155
0156   double SkyFunc(double p, double theta)
0157   {
0158     return mFunc(p, theta) * cos(theta) * sin(theta);
0159   }
0160   ///////////////////////////////////////////////////////////////
0161
0162   double CylFunc(double p, double theta)
0163   {
0164     return mFunc(p, theta) * pow(sin(theta), 2);
0165   }
0166   ///////////////////////////////////////////////////////////////
0167
0168   double HSFunc(double p, double theta, double phi, double theta0)
0169   {
0170     return mFunc(p, theta) * (sin(theta0) * sin(theta) * cos(phi) + cos(theta0) * cos(theta)) * sin(theta);
0171   }
0172   ///////////////////////////////////////////////////////////////
0173
0174   double Evaluate(std::vector<double>& v)
0175   {
0176     if (mGenMethod == 0)
0177     {
0178       return SkyFunc(v[0], v[1]);
0179     }
0180     else if (mGenMethod == 1)
0181     {
0182       return CylFunc(v[0], v[1]);
0183     }
0184     else
0185     {
0186       return HSFunc(v[0], v[1], v[2], v[3]);
0187     }
0188     return -1;
0189   }
0190   ///////////////////////////////////////////////////////////////
0191
0192   void Evaluate()
0193   {
0194     double value;
0195     for (size_t i = 0; i < mPopSize; ++i)
0196     {
0197       value = Evaluate(mPopulation[i]);
0198       if (value > mBestCost)
0199       {
0200         mBestCost = value;
0201         mBestSolution = mPopulation[i];
0202       }
0203     }
0204   }
0205   ///////////////////////////////////////////////////////////////
0206
0207   void Init()
0208   {
0209     size_t dim = mRanges.size();
0210     mPopulation.resize(mPopSize);
0211     for (size_t i = 0; i < mPopSize; ++i)
0212     {
0213       mPopulation[i].resize(dim);
0214       for (size_t j = 0; j < dim; ++j)
0215       {
0216         mPopulation[i][j] = mRandom.GenerateRandomDouble(mRanges[j][0], mRanges[j][1]);
0217       }
0218     }
0219   }
0220   ///////////////////////////////////////////////////////////////
0221
0222   void UpdateParameters(size_t t)
0223   {
0224     m_a = 2. - t * (2. / mNIter);
0225     m_a2 = -1. + t * ((-1.) / mNIter);
0226   }
0227   ///////////////////////////////////////////////////////////////
0228
0229   void Move()
0230   {
0231     double r1, r2, A, C, b, l, rw, p, D_tmp, D_best, distance;
0232     std::vector<double> tmp;
0233     for (size_t i = 0; i < mPopulation.size(); ++i)
0234     {
0235       r1 = mRandom.GenerateRandomDouble();
0236       r2 = mRandom.GenerateRandomDouble();
0237       A = 2 * m_a * r1 - m_a;
0238       C = 2 * r2;
0239       b = 1.;
0240       l = (m_a2 - 1) * mRandom.GenerateRandomDouble() + 1;
0241       p = mRandom.GenerateRandomDouble();
0242
0243       for (size_t j = 0; j < mPopulation[0].size(); ++j)
0244       {
0245         if (p < 0.5)
0246         {
0247           if (fabs(A) >= 1)
0248           {
0249             rw = floor(mRandom.GenerateRandomDouble() * mPopulation.size());
0250             tmp = mPopulation[rw];
0251             D_tmp = fabs(C * tmp[j] - mPopulation[i][j]);
0252             mPopulation[i][j] = tmp[j] - A * D_tmp;
0253           }
0254           else
0255           {
0256             D_best = fabs(C * mBestSolution[j] - mPopulation[i][j]);
0257             mPopulation[i][j] = mBestSolution[j] - A * D_best;
0258           }
0259         }
0260         else
0261         {
0262           distance = fabs(mBestSolution[j] - mPopulation[i][j]);
0263           mPopulation[i][j] = distance * exp(b * l) * cos(l * 2 * M_PI) + mBestSolution[j];
0264         }
0265         if (mPopulation[i][j] < mRanges[j][0]) mPopulation[i][j] = mRanges[j][0];
0266         if (mPopulation[i][j] > mRanges[j][1]) mPopulation[i][j] = mRanges[j][1];
0267       }
0268     }
0269   }
0270   ///////////////////////////////////////////////////////////////
0271
0272   double Maximize()
0273   {
0274     Init();
0275     Evaluate();
0276     for (size_t iter = 1; iter < mNIter; ++iter)
0277     {
0278       UpdateParameters(iter);
0279       Move();
0280       Evaluate();
0281     }
0282     return mBestCost;
0283   }
0284 };
0285 ///////////////////////////////////////////////////////////////
0286 ///////////////////////////////////////////////////////////////
0287
0288 //! Class for the generation of cosmic muons
0289 class EcoMug
0290 {
0291   friend class EMRandom;
0292
0293  public:
0294   /// Possible generation methods
0295   enum EMGeometry
0296   {
0297     Sky,       ///< generation from a plane (flat sky)
0298     Cylinder,  ///< generation from a cylinder
0299     HSphere    ///< generation from a half-sphere
0300   };
0301
0302  private:
0303   EMGeometry mGenMethod;
0304   std::array<double, 3> mGenerationPosition;
0305   double mGenerationTheta;
0306   double mGenerationPhi;
0307   double mGenerationMomentum;
0308   double mMinimumMomentum;
0309   double mMaximumMomentum;
0310   double mMinimumTheta;
0311   double mMaximumTheta;
0312   double mMinimumPhi;
0313   double mMaximumPhi;
0314   int mCharge;
0315   double mCylinderMinPositionPhi;
0316   double mCylinderMaxPositionPhi;
0317   double mHSphereMinPositionPhi;
0318   double mHSphereMaxPositionPhi;
0319   double mHSphereMinPositionTheta;
0320   double mHSphereMaxPositionTheta;
0321   double mHSphereCosMinPositionTheta;
0322   double mHSphereCosMaxPositionTheta;
0323   double mJPrime;
0324   double mN;
0325   double mRandAccRej;
0326   double mPhi0;
0327   double mTheta0;
0328   bool mAccepted;
0329   std::array<double, 2> mSkySize;
0330   std::array<double, 3> mSkyCenterPosition;
0331   double mCylinderHeight;
0332   double mCylinderRadius;
0333   std::array<double, 3> mCylinderCenterPosition;
0334   double mHSphereRadius;
0335   //  double mMaxFuncSkyCylinder;
0336   std::array<double, 3> mHSphereCenterPosition;
0337   EMRandom mRandom;
0338   std::array<double, 3> mMaxJ;
0339   std::array<double, 3> mMaxCustomJ;
0340   std::function<double(double, double)> mJ;
0341
0342  public:
0343   // Default constructor
0344   EcoMug()
0345     : mGenMethod(Sky)
0346     , mGenerationPosition({{0., 0., 0.}})
0347     , mGenerationTheta(0.)
0348     , mGenerationPhi(0.)
0349     , mGenerationMomentum(0.)
0350     , mMinimumMomentum(0.01)
0351     , mMaximumMomentum(1000.)
0352     , mMinimumTheta(0.)
0353     , mMaximumTheta(M_PI / 2.)
0354     , mMinimumPhi(0.)
0355     , mMaximumPhi(2. * M_PI)
0356     , mCharge(1)
0357     , mCylinderMinPositionPhi(0.)
0358     , mCylinderMaxPositionPhi(2. * M_PI)
0359     , mHSphereMinPositionPhi(0.)
0360     , mHSphereMaxPositionPhi(2. * M_PI)
0361     , mHSphereMinPositionTheta(0.)
0362     , mHSphereMaxPositionTheta(M_PI / 2.)
0363     , mHSphereCosMinPositionTheta(1.)
0364     , mHSphereCosMaxPositionTheta(0.)
0365     , mJPrime(0.)
0366     , mN(0.)
0367     , mRandAccRej(0.)
0368     , mPhi0(0.)
0369     , mTheta0(0.)
0370     , mAccepted(false)
0371     , mSkySize({{0., 0.}})
0372     , mSkyCenterPosition({{0., 0., 0.}})
0373     , mCylinderHeight(0.)
0374     , mCylinderRadius(0.)
0375     , mCylinderCenterPosition({{0., 0., 0.}})
0376     , mHSphereRadius(0.)
0377     ,
0378     /* mMaxFuncSkyCylinder(5.3176), */ mHSphereCenterPosition({{0., 0., 0.}})
0379   {
0380     // cppcheck-suppress [useInitializationList]
0381     mMaxJ = {-1., -1., -1.};
0382     // cppcheck-suppress [useInitializationList]
0383     mMaxCustomJ = {-1., -1., -1.};
0384   };
0385
0386   ///////////////////////////////////////////////////////////////
0387   // Methods to access the parameters of the generated muon
0388   ///////////////////////////////////////////////////////////////
0389   /// Get the generation position
0390   const std::array<double, 3>& GetGenerationPosition() const
0391   {
0392     return mGenerationPosition;
0393   };
0394   /// Get the generation momentum
0395   double GetGenerationMomentum() const
0396   {
0397     return mGenerationMomentum;
0398   };
0399   /// Get the generation momentum
0400   void GetGenerationMomentum(std::array<double, 3>& momentum) const
0401   {
0402     momentum = {
0403         mGenerationMomentum * sin(mGenerationTheta) * cos(mGenerationPhi),
0404         mGenerationMomentum * sin(mGenerationTheta) * sin(mGenerationPhi),
0405         mGenerationMomentum * cos(mGenerationTheta)};
0406   };
0407   /// Get the generation theta
0408   double GetGenerationTheta() const
0409   {
0410     return mGenerationTheta;
0411   };
0412   /// Get the generation phi
0413   double GetGenerationPhi() const
0414   {
0415     return mGenerationPhi;
0416   };
0417   /// Get charge
0418   int GetCharge() const
0419   {
0420     return mCharge;
0421   };
0422   ///////////////////////////////////////////////////////////////
0423
0424   ///////////////////////////////////////////////////////////////
0425   // Methods for the geometry of the generation
0426   ///////////////////////////////////////////////////////////////
0427   /// Set generation from sky
0428   void SetUseSky()
0429   {
0430     mGenMethod = Sky;
0431   };
0432   /// Set cylindrical generation
0433   void SetUseCylinder()
0434   {
0435     mGenMethod = Cylinder;
0436   };
0437   /// Set half-sphere generation
0438   void SetUseHSphere()
0439   {
0440     mGenMethod = HSphere;
0441   };
0442   /// Set the generation method (Sky, Cylinder or HSphere)
0443   void SetGenerationMethod(EMGeometry genM)
0444   {
0445     mGenMethod = genM;
0446   };
0447
0448   /// Get the generation method (Sky, Cylinder or HSphere)
0449   EMGeometry GetGenerationMethod() const
0450   {
0451     return mGenMethod;
0452   };
0453   ///////////////////////////////////////////////////////////////
0454
0455   ///////////////////////////////////////////////////////////////
0456   // Common methods to all geometries
0457   ///////////////////////////////////////////////////////////////
0458   /// Set the differential flux J. Accepted functions are like
0459   /// double J(double momentum, double theta)
0460   /// momentum has to be in GeV/c and theta in radians
0461   void SetDifferentialFlux(std::function<double(double, double)> J)
0462   {
0463     mJ = J;
0464   };
0465   /// Set the seed for the internal PRNG (if 0 a random seed is used)
0466   void SetSeed(uint64_t seed)
0467   {
0468     if (seed > 0) mRandom.SetSeed(seed);
0469   };
0470   /// Set minimum generation Momentum
0471   void SetMinimumMomentum(double momentum)
0472   {
0473     mMinimumMomentum = momentum;
0474   };
0475   /// Set maximum generation Momentum
0476   void SetMaximumMomentum(double momentum)
0477   {
0478     mMaximumMomentum = momentum;
0479   };
0480   /// Set minimum generation Theta
0481   void SetMinimumTheta(double theta)
0482   {
0483     mMinimumTheta = theta;
0484   };
0485   /// Set maximum generation Theta
0486   void SetMaximumTheta(double theta)
0487   {
0488     mMaximumTheta = theta;
0489   };
0490   /// Set minimum generation Phi
0491   void SetMinimumPhi(double phi)
0492   {
0493     mMinimumPhi = phi;
0494   };
0495   /// Set maximum generation Phi
0496   void SetMaximumPhi(double phi)
0497   {
0498     mMaximumPhi = phi;
0499   };
0500
0501   /// Get minimum generation Momentum
0502   double GetMinimumMomentum() const
0503   {
0504     return mMinimumMomentum;
0505   };
0506   /// Get maximum generation Momentum
0507   double GetMaximumMomentum() const
0508   {
0509     return mMaximumMomentum;
0510   };
0511   /// Get minimum generation Theta
0512   double GetMinimumTheta() const
0513   {
0514     return mMinimumTheta;
0515   };
0516   /// Get maximum generation Theta
0517   double GetMaximumTheta() const
0518   {
0519     return mMaximumTheta;
0520   };
0521   /// Get minimum generation Phi
0522   double GetMinimumPhi() const
0523   {
0524     return mMinimumPhi;
0525   };
0526   /// Get maximum generation Phi
0527   double GetMaximumPhi() const
0528   {
0529     return mMaximumPhi;
0530   };
0531   ///////////////////////////////////////////////////////////////
0532
0533   ///////////////////////////////////////////////////////////////
0534   // Methods for the plane-based generation
0535   ///////////////////////////////////////////////////////////////
0536   /// Set sky size
0537   void SetSkySize(const std::array<double, 2>& size)
0538   {
0539     mSkySize = size;
0540   };
0541
0542   /// Set sky center position
0543   void SetSkyCenterPosition(const std::array<double, 3>& position)
0544   {
0545     mSkyCenterPosition = position;
0546   };
0547   ///////////////////////////////////////////////////////////////
0548
0549   ///////////////////////////////////////////////////////////////
0550   // Methods for the cylinder-based generation
0551   ///////////////////////////////////////////////////////////////
0552   /// Set cylinder radius
0553   void SetCylinderRadius(double radius)
0554   {
0555     mCylinderRadius = radius;
0556   };
0557   /// Set cylinder height
0558   void SetCylinderHeight(double height)
0559   {
0560     mCylinderHeight = height;
0561   };
0562   /// Set cylinder center position
0563   void SetCylinderCenterPosition(const std::array<double, 3>& position)
0564   {
0565     mCylinderCenterPosition = position;
0566   };
0567   void SetCylinderMinPositionPhi(double phi)
0568   {
0569     mCylinderMinPositionPhi = phi;
0570   };
0571   void SetCylinderMaxPositionPhi(double phi)
0572   {
0573     mCylinderMaxPositionPhi = phi;
0574   };
0575   /// Get cylinder radius
0576   double GetCylinderRadius() const
0577   {
0578     return mCylinderRadius;
0579   };
0580   /// Get cylinder height
0581   double GetCylinderHeight() const
0582   {
0583     return mCylinderHeight;
0584   };
0585   /// Get cylinder center position
0586   const std::array<double, 3>& GetCylinderCenterPosition() const
0587   {
0588     return mCylinderCenterPosition;
0589   };
0590   ///////////////////////////////////////////////////////////////
0591
0592   ///////////////////////////////////////////////////////////////
0593   // Methods for the half sphere-based generation
0594   ///////////////////////////////////////////////////////////////
0595   /// Set half-sphere radius
0596   void SetHSphereRadius(double radius)
0597   {
0598     mHSphereRadius = radius;
0599   };
0600   /// Set half-sphere center position
0601   void SetHSphereCenterPosition(const std::array<double, 3>& position)
0602   {
0603     mHSphereCenterPosition = position;
0604   };
0605   void SetHSphereMinPositionPhi(double phi)
0606   {
0607     mHSphereMinPositionPhi = phi;
0608   };
0609   void SetHSphereMaxPositionPhi(double phi)
0610   {
0611     mHSphereMaxPositionPhi = phi;
0612   };
0613   void SetHSphereMinPositionTheta(double theta)
0614   {
0615     mHSphereMinPositionTheta = theta;
0616     mHSphereCosMinPositionTheta = cos(mHSphereMinPositionTheta);
0617   };
0618   void SetHSphereMaxPositionTheta(double theta)
0619   {
0620     mHSphereMaxPositionTheta = theta;
0621     mHSphereCosMaxPositionTheta = cos(mHSphereMaxPositionTheta);
0622   };
0623   /// Get half-sphere radius
0624   double GetHSphereRadius() const
0625   {
0626     return mHSphereRadius;
0627   };
0628   /// Get half-sphere center position
0629   const std::array<double, 3>& GetHSphereCenterPosition() const
0630   {
0631     return mHSphereCenterPosition;
0632   };
0633   ///////////////////////////////////////////////////////////////
0634
0635  private:
0636   double F1Cumulative(double x)
0637   {
0638     return 1. - 8.534790171171021 / pow(x + 2.68, 87. / 40.);
0639   };
0640
0641   double F1Inverse(double x)
0642   {
0643     return (2.68 - 2.68 * pow(1. - x, 40. / 87.)) / pow(1. - x, 40. / 87.);
0644   };
0645
0646   double maxSkyJFunc()
0647   {
0648     return 1600 * pow(mMaximumMomentum, 0.279) * pow(cos(0.76158), 1.1) * sin(0.76158);
0649   };
0650
0651   double maxCylJFunc()
0652   {
0653     return 1600 * pow(mMaximumMomentum, 0.279) * pow(cos(1.35081), 0.1) * pow(sin(1.35081), 2);
0654   };
0655
0656   double maxHSJFunc()
0657   {
0658     return 1600 * pow(mMaximumMomentum, 0.279) * pow(cos(1.26452), 0.1) * (sin(1.26452) * sin(1.26452) + cos(1.26452) * cos(1.26452)) * sin(1.26452);
0659   };
0660
0661   double GenerateMomentumF1()
0662   {
0663     double z = mRandom.GenerateRandomDouble(F1Cumulative(mMinimumMomentum), F1Cumulative(mMaximumMomentum));
0664     return F1Inverse(z);
0665   };
0666
0667   void GeneratePositionSky()
0668   {
0669     mGenerationPosition[0] = mRandom.GenerateRandomDouble(mSkyCenterPosition[0] - mSkySize[0] / 2., mSkyCenterPosition[0] + mSkySize[0] / 2.);
0670     mGenerationPosition[1] = mRandom.GenerateRandomDouble(mSkyCenterPosition[1] - mSkySize[1] / 2., mSkyCenterPosition[1] + mSkySize[1] / 2.);
0671     mGenerationPosition[2] = mSkyCenterPosition[2];
0672   };
0673
0674   void GeneratePositionCylinder()
0675   {
0676     mPhi0 = mRandom.GenerateRandomDouble(mCylinderMinPositionPhi, mCylinderMaxPositionPhi);
0677     mGenerationPosition[0] = mCylinderCenterPosition[0] + mCylinderRadius * cos(mPhi0);
0678     mGenerationPosition[1] = mCylinderCenterPosition[1] + mCylinderRadius * sin(mPhi0);
0679     mGenerationPosition[2] = mRandom.GenerateRandomDouble(mCylinderCenterPosition[2] - mCylinderHeight / 2., mCylinderCenterPosition[2] + mCylinderHeight / 2.);
0680   };
0681
0682   void ComputeMaximumCustomJ()
0683   {
0684     EMMaximization maximizer(mRandom, mGenMethod);
0685     maximizer.SetFunction(mJ);
0686     if (mGenMethod == 0 || mGenMethod == 1)
0687     {
0688       maximizer.SetParameters(mMinimumMomentum, mMaximumMomentum, mMinimumTheta, mMaximumTheta);
0689     }
0690     else
0691     {
0692       maximizer.SetParameters(mMinimumMomentum, mMaximumMomentum, mMinimumTheta, mMaximumTheta, mMinimumPhi, mMaximumPhi);
0693     }
0694     mMaxCustomJ[mGenMethod] = maximizer.Maximize();
0695   };
0696
0697   void ComputeMaximum()
0698   {
0699     EMMaximization maximizer(mRandom, mGenMethod);
0700     if (mGenMethod == 0 || mGenMethod == 1)
0701     {
0702       maximizer.SetParameters(mMinimumMomentum, mMaximumMomentum, mMinimumTheta, mMaximumTheta);
0703     }
0704     else
0705     {
0706       maximizer.SetParameters(mMinimumMomentum, mMaximumMomentum, mMinimumTheta, mMaximumTheta, mMinimumPhi, mMaximumPhi);
0707     }
0708     mMaxJ[mGenMethod] = maximizer.Maximize();
0709   };
0710
0711  public:
0712   ///////////////////////////////////////////////////////////////
0713   /// Generate a cosmic muon from the pre-defined J
0714   ///////////////////////////////////////////////////////////////
0715   void Generate()
0716   {
0717     mAccepted = false;
0718
0719     if (mMaxJ[mGenMethod] < 0) ComputeMaximum();
0720
0721     // Sky or cylinder generation
0722     if (mGenMethod == Sky || mGenMethod == Cylinder)
0723     {
0724       // Generation of the momentum and theta angle
0725       while (!mAccepted)
0726       {
0727         mRandAccRej = mRandom.GenerateRandomDouble();
0728         mGenerationTheta = mRandom.GenerateRandomDouble(mMinimumTheta, mMaximumTheta);
0729         mGenerationMomentum = GenerateMomentumF1();
0730         mN = 2.856 - 0.655 * log(mGenerationMomentum);
0731         if (mN < 0.1) mN = 0.1;
0732
0733         if (mGenMethod == Sky)
0734         {
0735           mJPrime = 1600 * pow(mGenerationMomentum, 0.279) * pow(cos(mGenerationTheta), mN + 1) * sin(mGenerationTheta);
0736           if (mMaxJ[mGenMethod] * mRandAccRej < mJPrime) mAccepted = true;
0737         }
0738
0739         if (mGenMethod == Cylinder)
0740         {
0741           mJPrime = 1600 * pow(mGenerationMomentum, 0.279) * pow(cos(mGenerationTheta), mN) * pow(sin(mGenerationTheta), 2);
0742           if (mMaxJ[mGenMethod] * mRandAccRej < mJPrime) mAccepted = true;
0743         }
0744       }
0745       mGenerationTheta = M_PI - mGenerationTheta;
0746
0747       // Generation of the position and phi angle
0748       if (mGenMethod == Sky)
0749       {
0750         GeneratePositionSky();
0751         mGenerationPhi = mRandom.GenerateRandomDouble(mMinimumPhi, mMaximumPhi);
0752       }
0753       if (mGenMethod == Cylinder)
0754       {
0755         mAccepted = false;
0756         GeneratePositionCylinder();
0757         while (!mAccepted)
0758         {
0759           mRandAccRej = mRandom.GenerateRandomDouble();
0760           mGenerationPhi = mRandom.GenerateRandomDouble(mMinimumPhi, mMaximumPhi);
0761           if (mRandAccRej < fabs(cos(mGenerationPhi))) mAccepted = true;
0762         }
0763         mGenerationPhi = mGenerationPhi + mPhi0;
0764         if (mGenerationPhi >= 2. * M_PI) mGenerationPhi -= 2. * M_PI;
0765
0766         // Check if the muon is inward
0767         if (sin(mGenerationTheta) * cos(mGenerationPhi) * mGenerationPosition[0] + sin(mGenerationTheta) * sin(mGenerationPhi) * mGenerationPosition[1] > 0) Generate();
0768       }
0769     }
0770
0771     // Half-sphere generation
0772     if (mGenMethod == HSphere)
0773     {
0774       // Generation point on the half-sphere
0775       mPhi0 = mRandom.GenerateRandomDouble(mHSphereMinPositionPhi, mHSphereMaxPositionPhi);
0776       while (!mAccepted)
0777       {
0778         mRandAccRej = mRandom.GenerateRandomDouble();
0779         mTheta0 = acos(mRandom.GenerateRandomDouble(mHSphereCosMaxPositionTheta, mHSphereCosMinPositionTheta));
0780         mGenerationTheta = mRandom.GenerateRandomDouble(mMinimumTheta, mMaximumTheta);
0781         mGenerationPhi = mRandom.GenerateRandomDouble(mMinimumPhi, mMaximumPhi);
0782         mGenerationMomentum = GenerateMomentumF1();
0783         mN = 2.856 - 0.655 * log(mGenerationMomentum);
0784         if (mN < 0.1) mN = 0.1;
0785
0786         mJPrime = 1600 * pow(mGenerationMomentum, 0.279) * pow(cos(mGenerationTheta), mN) * (sin(mGenerationTheta) * sin(mTheta0) * cos(mGenerationPhi) + cos(mGenerationTheta) * cos(mTheta0)) * sin(mGenerationTheta);
0787         if (mMaxJ[mGenMethod] * mRandAccRej < mJPrime) mAccepted = true;
0788       }
0789
0790       mGenerationPosition[0] = mHSphereRadius * sin(mTheta0) * cos(mPhi0) + mHSphereCenterPosition[0];
0791       mGenerationPosition[1] = mHSphereRadius * sin(mTheta0) * sin(mPhi0) + mHSphereCenterPosition[1];
0792       mGenerationPosition[2] = mHSphereRadius * cos(mTheta0) + mHSphereCenterPosition[2];
0793
0794       mGenerationTheta = M_PI - mGenerationTheta;
0795       mGenerationPhi = mGenerationPhi + mPhi0;
0796       if (mGenerationPhi >= 2 * M_PI) mGenerationPhi -= 2 * M_PI;
0797
0798       mGenerationPhi += M_PI;
0799       if (mGenerationPhi >= 2 * M_PI) mGenerationPhi -= 2 * M_PI;
0800     }
0801
0802     // Generate the charge
0803     if (mRandom.GenerateRandomDouble(-100, 128) >= 0)
0804       mCharge = 1;
0805     else
0806       mCharge = -1;
0807   };
0808   ///////////////////////////////////////////////////////////////
0809
0810   ///////////////////////////////////////////////////////////////
0811   /// Generate a cosmic muon for the user-defined J
0812   ///////////////////////////////////////////////////////////////
0813   void GenerateFromCustomJ()
0814   {
0815     mAccepted = false;
0816
0817     if (mMaxCustomJ[mGenMethod] < 0) ComputeMaximumCustomJ();
0818
0819     // Sky or cylinder generation
0820     if (mGenMethod == Sky || mGenMethod == Cylinder)
0821     {
0822       // Generation of the momentum and theta angle
0823       while (!mAccepted)
0824       {
0825         mRandAccRej = mRandom.GenerateRandomDouble();
0826         mGenerationTheta = mRandom.GenerateRandomDouble(mMinimumTheta, mMaximumTheta);
0827         mGenerationMomentum = mRandom.GenerateRandomDouble(mMinimumMomentum, mMaximumMomentum);
0828
0829         if (mGenMethod == Sky)
0830         {
0831           mJPrime = mJ(mGenerationMomentum, mGenerationTheta) * cos(mGenerationTheta) * sin(mGenerationTheta);
0832           if (mMaxCustomJ[mGenMethod] * mRandAccRej < mJPrime) mAccepted = true;
0833         }
0834
0835         if (mGenMethod == Cylinder)
0836         {
0837           mJPrime = mJ(mGenerationMomentum, mGenerationTheta) * pow(sin(mGenerationTheta), 2) * cos(mGenerationPhi);
0838           if (mMaxCustomJ[mGenMethod] * mRandAccRej < mJPrime) mAccepted = true;
0839         }
0840       }
0841       mGenerationTheta = M_PI - mGenerationTheta;
0842
0843       // Generation of the position and phi angle
0844       if (mGenMethod == Sky)
0845       {
0846         GeneratePositionSky();
0847         mGenerationPhi = mRandom.GenerateRandomDouble(mMinimumPhi, mMaximumPhi);
0848       }
0849       if (mGenMethod == Cylinder)
0850       {
0851         mAccepted = false;
0852         GeneratePositionCylinder();
0853         while (!mAccepted)
0854         {
0855           mRandAccRej = mRandom.GenerateRandomDouble();
0856           mGenerationPhi = mRandom.GenerateRandomDouble(mMinimumPhi, mMaximumPhi);
0857           if (mRandAccRej < fabs(cos(mGenerationPhi))) mAccepted = true;
0858         }
0859         mGenerationPhi = mGenerationPhi + mPhi0;
0860         if (mGenerationPhi >= 2. * M_PI) mGenerationPhi -= 2. * M_PI;
0861
0862         // Check if the muon is inward
0863         if (sin(mGenerationTheta) * cos(mGenerationPhi) * mGenerationPosition[0] + sin(mGenerationTheta) * sin(mGenerationPhi) * mGenerationPosition[1] > 0) Generate();
0864       }
0865     }
0866
0867     // Half-sphere generation
0868     if (mGenMethod == HSphere)
0869     {
0870       // Generation point on the half-sphere
0871       mPhi0 = mRandom.GenerateRandomDouble(mHSphereMinPositionPhi, mHSphereMaxPositionPhi);
0872       while (!mAccepted)
0873       {
0874         mRandAccRej = mRandom.GenerateRandomDouble();
0875         mTheta0 = acos(mRandom.GenerateRandomDouble(mHSphereCosMaxPositionTheta, mHSphereCosMinPositionTheta));
0876         mGenerationTheta = mRandom.GenerateRandomDouble(mMinimumTheta, mMaximumTheta);
0877         mGenerationPhi = mRandom.GenerateRandomDouble(mMinimumPhi, mMaximumPhi);
0878         mGenerationMomentum = mRandom.GenerateRandomDouble(mMinimumMomentum, mMaximumMomentum);
0879
0880         mJPrime = mJ(mGenerationMomentum, mGenerationTheta) * (sin(mTheta0) * sin(mGenerationTheta) * cos(mGenerationPhi) + cos(mTheta0) * cos(mGenerationTheta)) * sin(mGenerationTheta);
0881         if (mMaxCustomJ[mGenMethod] * mRandAccRej < mJPrime) mAccepted = true;
0882       }
0883
0884       mGenerationPosition[0] = mHSphereRadius * sin(mTheta0) * cos(mPhi0) + mHSphereCenterPosition[0];
0885       mGenerationPosition[1] = mHSphereRadius * sin(mTheta0) * sin(mPhi0) + mHSphereCenterPosition[1];
0886       mGenerationPosition[2] = mHSphereRadius * cos(mTheta0) + mHSphereCenterPosition[2];
0887
0888       mGenerationTheta = M_PI - mGenerationTheta;
0889       mGenerationPhi = mGenerationPhi + mPhi0;
0890       if (mGenerationPhi >= 2 * M_PI) mGenerationPhi -= 2 * M_PI;
0891
0892       mGenerationPhi += M_PI;
0893       if (mGenerationPhi >= 2 * M_PI) mGenerationPhi -= 2 * M_PI;
0894     }
0895
0896     // Generate the charge
0897     if (mRandom.GenerateRandomDouble(-100, 128) >= 0)
0898       mCharge = 1;
0899     else
0900       mCharge = -1;
0901   };
0902   ///////////////////////////////////////////////////////////////
0903 };
0904 ///////////////////////////////////////////////////////////////
0905 ///////////////////////////////////////////////////////////////
0906
0907 #endif