CMS 3D CMS Logo

/data/refman/pasoursint/CMSSW_5_3_3/src/FastSimulation/Utilities/src/GammaFunctionGenerator.cc

Go to the documentation of this file.
00001 //FAMOS headers
00002 #include "FastSimulation/Utilities/interface/GammaFunctionGenerator.h"
00003 #include "FastSimulation/Utilities/interface/RandomEngine.h"
00004 
00005 GammaFunctionGenerator::GammaFunctionGenerator(const RandomEngine* engine) :
00006   random(engine)
00007 {
00008 
00009   xmax = 30.;
00010 
00011   for(unsigned i=1;i<=12;++i)
00012     {
00013       // For all let's put the limit at 2.*(alpha-1)      (alpha-1 is the max of the dist)
00014       approxLimit.push_back(2*((double)i));
00015       myIncompleteGamma.a().setValue((double)i);
00016       integralToApproxLimit.push_back(myIncompleteGamma(approxLimit[i-1]));
00017       theGammas.push_back(
00018        GammaNumericalGenerator(random,(double)i,1.,0,approxLimit[i-1]+1.));
00019     }
00020   coreCoeff.push_back(0.);  // alpha=1 not used
00021   coreCoeff.push_back(1./8.24659e-01);
00022   coreCoeff.push_back(1./7.55976e-01);
00023   coreCoeff.push_back(1./7.12570e-01);
00024   coreCoeff.push_back(1./6.79062e-01);
00025   coreCoeff.push_back(1./6.65496e-01);
00026   coreCoeff.push_back(1./6.48736e-01);
00027   coreCoeff.push_back(1./6.25185e-01);
00028   coreCoeff.push_back(1./6.09188e-01);
00029   coreCoeff.push_back(1./6.06221e-01);
00030   coreCoeff.push_back(1./6.05057e-01);
00031 }
00032 
00033 GammaFunctionGenerator::~GammaFunctionGenerator() {}
00034 
00035 double GammaFunctionGenerator::shoot() const
00036 {
00037   if(alpha<0.) return -1.;
00038   if(badRange) return xmin/beta;
00039   if(alpha<12)
00040     {
00041 
00042       if (alpha == na)
00043         {
00044           return gammaInt ()/beta;
00045         }
00046       else if (na == 0)
00047         {
00048           return gammaFrac ()/beta;
00049         }
00050       else
00051         {
00052           double gi=gammaInt ();
00053           double gf=gammaFrac ();
00054           return (gi+gf)/beta;
00055         }
00056     }
00057   else
00058     {
00059       // an other generator has to be used in such a case
00060       return -1.;
00061     }
00062 }
00063 
00064 double GammaFunctionGenerator::gammaFrac () const
00065 {
00066   /* This is exercise 16 from Knuth; see page 135, and the solution is
00067      on page 551.  */
00068 
00069   double p, q, x, u, v;
00070   p = M_E / (frac + M_E);
00071   do
00072     {
00073       u = random->flatShoot();
00074       v = random->flatShoot();
00075 
00076       if (u < p)
00077         {
00078           x = exp ((1 / frac) * log (v));
00079           q = exp (-x);
00080         }
00081       else
00082         {
00083           x = 1 - log (v);
00084           q = exp ((frac - 1) * log (x));
00085         }
00086     }
00087   while (random->flatShoot() >= q);
00088 
00089   return x;
00090 }
00091 
00092 double GammaFunctionGenerator::gammaInt() const
00093 {
00094   // Exponential distribution : no approximation
00095   if(na==1)
00096     {
00097       return xmin-log(random->flatShoot());
00098     }
00099 
00100   unsigned gn=na-1;
00101 
00102   // are we sure to be in the tail 
00103   if(coreProba==0.)
00104     return xmin-coreCoeff[gn]*log(random->flatShoot());
00105 
00106   // core-tail interval
00107   if(random->flatShoot()<coreProba)
00108     {
00109       return theGammas[gn].gamma_lin();
00110     }
00111   //  std::cout << " Tail called " << std::endl;
00112   return approxLimit[gn]-coreCoeff[gn]*log(random->flatShoot());
00113 }
00114 
00115 void GammaFunctionGenerator::setParameters(double a,double b, double xm)
00116 {
00117   //  std::cout << "Setting parameters " << std::endl;
00118   alpha=a;
00119   beta=b;
00120   xmin=xm*beta;
00121   if(xm>xmax) 
00122     {
00123       badRange=true;
00124       return;
00125     }
00126   badRange=false;
00127   na=0;
00128   
00129   if(alpha>0.&&alpha<12) 
00130     na=(unsigned)floor(alpha);
00131 
00132   frac=alpha-na;  
00133   // Now calculate the probability to shoot between approxLimit and xmax
00134   // The Incomplete gamma is normalized to 1 
00135   if(na<=1) return;
00136   
00137   myIncompleteGamma.a().setValue((double)na);
00138   
00139   unsigned gn=na-1;
00140   //  std::cout << " na " << na << " xm " << xm << " beta " << beta << " xmin " << xmin << " approxLimit " << approxLimit[gn] << std::endl;
00141   if(xmin>approxLimit[gn]) 
00142     {
00143       coreProba=0.;
00144     }
00145   else
00146     {
00147       double tmp=(xmin!=0.) ?myIncompleteGamma(xmin) : 0.;
00148       coreProba=(integralToApproxLimit[gn]-tmp)/(1.-tmp);
00149       theGammas[gn].setSubInterval(xmin,approxLimit[gn]);
00150     }
00151   //  std::cout << " Proba " << coreProba << std::endl;
00152 }