Impala Documentation

template<class HistT>

void Impala::Core::Histogram::FitWeibullMarginalOld< HistT >::DoFit	(	HistT *	orgHist,
		int	orgBins
	)			`[inline, private]`

Definition at line 63 of file FitWeibullMarginalOld.h.

Referenced by Impala::Core::Histogram::FitWeibullMarginalOld< HistT >::FitWeibullMarginalOld().

00064     {
00065         // find mu for folded histogram (2 bins together)
00066         double x;
00067         double mode = 0.0;
00068         int bin0 = 0;
00069         int i;
00070         mMu = 0.0;
00071         for (i=0, x = mMinval+0.5*mDx; i<orgBins; i++, x += mDx)
00072             mMu += orgHist[i]*x;
00073         bin0 = (int)((mMu-mMinval)/mDx);
00074 
00075         // fill new histogram with mirrored data
00076         mBins = bin0>orgBins/2 ? bin0 : orgBins-bin0;
00077         mHist = new HistT [mBins];
00078         for (i=0; i<mBins; i++)
00079             mHist[i] = 0;
00080         for (i=0; i<orgBins; i++)
00081         {
00082             int j = bin0 > i ? bin0-i-1 : i-bin0;
00083             mHist[j] += orgHist[i];
00084         }
00085 
00086         double sumy(0), sumysq(0);
00087         double gammal, gammah;
00088         double xgm, gfm, tol;
00089 
00090         mNorm = 0;
00091         for (i=0, x = mDx*0.5; i<mBins; i++, x += mDx)
00092         {
00093             if (mHist[i]>mPrecision)
00094             {
00095                 double logy = log(x);
00096                 sumy += mHist[i]*logy;
00097                 sumysq += mHist[i]*logy*logy;
00098                 mNorm += mHist[i];
00099             }
00100         }
00101 
00102         sumy /= mNorm;
00103         sumysq /= mNorm;
00104         mGamma = 1.28/sqrt((sumysq-sumy*sumy)*mBins/(mBins-1.0));
00105         xgm = exp(sumy);
00106         //    tol = 2.0*0.000001*gamma;
00107         tol = 2.0*0.001*mGamma;
00108         gfm = GFunct(mGamma,xgm);
00109 
00110         int its = 0;
00111 
00112         if (gfm >= 0.0)
00113         {
00114             while (gfm > 0.0)
00115             {
00116                 gammah = mGamma;
00117                 mGamma /= 2.0;
00118                 if (mGamma < tol)
00119                     break;
00120                 gfm = GFunct(mGamma,xgm);
00121                 if (its++ > 30)
00122                     break;
00123             }
00124             gammal = mGamma;
00125         }
00126         its = 0;
00127         while (gfm < 0.0)
00128         {
00129             gammal = mGamma;
00130             mGamma *= 2.0;
00131             gfm = GFunct(mGamma,xgm);
00132             if (its++ > 30)
00133                 break;
00134         }
00135         gammah = mGamma;
00136 
00137         its = 0;
00138         while ((gammah-gammal) > tol)
00139         {
00140             mGamma = (gammal+gammah)/2.0;
00141             gfm = GFunct(mGamma,xgm);
00142             if (gfm >= 0.0)
00143                 gammah = mGamma;
00144             else
00145                 gammal = mGamma;
00146             if (its++ > 30)
00147                 break;
00148         }
00149 
00150         mGamma = (gammal+gammah)/2.0;
00151         mBeta = BetaEst(mGamma,xgm);
00152     }

Here is the call graph for this function: