//////////////////////////////////////////////////////////////////////////////
//
//  TRolke
//
//  This class computes confidence intervals for the rate of a Poisson
//  process in the presence of uncertain background and/or efficiency.
//
//  The treatment and the resulting limits are fully frequentist. The
//  limit calculations make use of the profile likelihood method.
//
//      Author: Jan Conrad (CERN) 2004
//      Updated: Johan Lundberg (CERN) 2009
//
//      Copyright CERN 2004,2009           Jan.Conrad@cern.ch,
//                                     Johan.Lundberg@cern.ch
//
//  For a full list of methods and their syntax, and build instructions,
//  consult the header file TRolke.h.
//  --------------------------------
//
//  Examples/tutorials are found in the separate file Rolke.C
//  ---------------------------------------------------------
//
//////////////////////////////////////////////////////////////////////////////
//
//
//  TRolke implements the following Models
//  =======================================
//
//  The signal is always assumed to be Poisson, with the following
//  combinations of models of background and detection efficiency:
//
//  If unsure, first consider model 3, 4 or 5.
//
//       1: SetPoissonBkgBinomEff(x,y,z,tau,m)
//
//          Background: Poisson
//          Efficiency: Binomial
//
//          when the background is simultaneously measured
//          from sidebands (or MC), and
//          the signal efficiency was determined from Monte Carlo
//
//       2: SetPoissonBkgGaussEff(x,y,em,sde,tau)
//
//          Background: Poisson
//          Efficiency: Gaussian
//
//          when the background is simultaneously measured
//          from sidebands (or MC), and
//          the efficiency is modeled as Gaussian
//
//       3: SetGaussBkgGaussEff(x,bm,em,sde,sdb)
//
//          Background: Gaussian
//          Efficiency: Gaussian
//
//          when background and efficiency can both be
//          modeled as Gaussian.
//
//       4: SetPoissonBkgKnownEff(x,y,tau,e)
//
//          Background: Poisson
//          Efficiency: Known
//
//          when the background is simultaneously measured
//          from sidebands (or MC).
//
//       5: SetGaussBkgKnownEff(x,bm,sdb,e)
//
//          Background: Gaussian
//          Efficiency: Known
//
//          when background is Gaussian
//
//       6: SetKnownBkgBinomEff(x,z,b,m)
//
//          Background: Known
//          Efficiency: Binomial
//
//          when signal efficiency was determined from Monte Carlo
//
//       7: SetKnownBkgGaussEff(x,em,sde,b)
//
//          Background: Known
//          Efficiency: Gaussian
//
//          when background is known and efficiency Gaussian
//
//  Parameters and further explanation
//  ==================================
//
//  For all models:
//  ---------------
//
//    x = number of observed events in the experiment
//
//    Efficiency (e or em) is the detection probability for signal.
//    A low efficiency hence generally means weaker limits.
//    If the efficiency of an experiment (with analysis cuts) is
//    dealt with elsewhere, em or e can be set to one.
//
//  For Poisson background measurements (sideband or MC):
//  -----------------------------------------------------
//
//    y = number of observed events in background region
//    tau =
//        Either: the ratio between signal and background region
//        in case background is observed.
//        Or: the ratio between observed and simulated live-time
//        in case background is determined from MC.
//
//  For Gaussian efficiency or background:
//  --------------------------------------
//
//    bm  = estimate of the background
//    sdb = corresponding standard deviation
//
//    em  = estimate of the efficiency
//    sde = corresponding standard deviation
//
//        If the efficiency scale of dealt with elsewhere,
//        set em to 1 and sde to the relative uncertainty.
//
//  For Binomial signal efficiency:
//  -------------------------------
//
//     m = number of MC events generated
//     z = number of MC events observed
//
//  For the case of known background expectation or known efficiency:
//  -----------------------------------------------------------------
//
//     e = true efficiency (considered known)
//     b = background expectation value (considered known)
//
//
////////////////////////////////////////////////////////////////////
//
//
//  The confidence level (CL) is set either at construction
//  time or with either of SetCL or SetCLSigmas
//
//  The TRolke method is very similar to the one used in MINUIT (MINOS).
//
//  Two options are offered to deal with cases where the maximum likelihood
//  estimate (MLE) is not in the physical region. Version "bounded likelihood"
//  is the one used by MINOS if bounds for the physical region are chosen.
//  Unbounded likelihood (the default) allows the MLE to be in the
//  unphysical region. It has however better coverage.
//  For more details consult the reference (see below).
//
//  For a description of the method and its properties:
//
//  W.Rolke, A. Lopez, J. Conrad and Fred James
//  "Limits and Confidence Intervals in presence of nuisance parameters"
//   http://lanl.arxiv.org/abs/physics/0403059
//   Nucl.Instrum.Meth.A551:493-503,2005
//
//  Should I use TRolke, TFeldmanCousins, TLimit?
//  ----------------------------------------------
//
//  1. Does TRolke make TFeldmanCousins obsolete?
//
//  Certainly not. TFeldmanCousins is the fully frequentist construction and
//  should be used in case of no (or negligible) uncertainties. It is however
//  not capable of treating uncertainties in nuisance parameters. In other
//  words, it does not handle background expectations or signal efficiencies
//  which are known only with some limited accuracy.
//
//  TRolke is designed for this case and it is shown in the reference above
//  that it has good coverage properties for most cases, and can be used
//  where FeldmannCousins can't.
//
//  2. What are the advantages of TRolke over TLimit?
//
//  TRolke is fully frequentist. TLimit treats nuisance parameters Bayesian.
//
//  For a coverage study of a Bayesian method refer to
//  physics/0408039 (Tegenfeldt & J.C). However, this note studies
//  the coverage of Feldman&Cousins with Bayesian treatment of nuisance
//  parameters. To make a long story short: using the Bayesian method you
//  might introduce a small amount of over-coverage (though I haven't shown it
//  for TLimit). On the other hand, coverage of course is a not so interesting
//  when you consider yourself a Bayesian.
//
///////////////////////////////////////////////////////////////////////////

#include "TRolke.h"
#include "TMath.h"
#include "Riostream.h"

ClassImp(TRolke)

//__________________________________________________________________________
TRolke::TRolke(Double_t CL, Option_t * /*option*/)
:  fCL(CL),
   fUpperLimit(0.0),
   fLowerLimit(0.0),
   fBounding(false),  // true gives bounded likelihood
   fNumWarningsDeprecated1(0),
   fNumWarningsDeprecated2(0)
{
//constructor with optional Confidence Level argument.
//'option' is not used.

   SetModelParameters();
}

//___________________________________________________________________________
TRolke::~TRolke()
{
}

/* ______________________ new interface _____________________ */
void TRolke::SetPoissonBkgBinomEff(Int_t x, Int_t y, Int_t z, Double_t tau, Int_t m)
{
// Model 1: Background - Poisson, Efficiency - Binomial
//    x   : number of observed events in the experiment
//    y   : number of observed events in background region
//    z   : number of MC events observed
//    tau : ratio parameter (read TRolke.cxx for details)
//    m   : number of MC events generated
   SetModelParameters(
         x  ,       //   Int_t x,
         y  ,       //   Int_t y,
         z  ,       //   Int_t z,
         0  ,       //   Double_t bm,
         0  ,       //   Double_t em,
         0  ,       //   Double_t e,
         1  ,       //   Int_t mid,
         0  ,       //   Double_t sde,
         0  ,       //   Double_t sdb,
         tau,       //   Double_t tau,
         0  ,       //   Double_t b,
         m);        //   Int_t m
}

void TRolke::SetPoissonBkgGaussEff(Int_t x, Int_t y, Double_t em, Double_t tau, Double_t sde)
{
// Model 2: Background - Poisson, Efficiency - Gaussian
//    x   : number of observed events in the experiment
//    y   : number of observed events in background region
//    em  : estimate of the efficiency
//    tau : ratio parameter (read TRolke.cxx for details)
//    sde : efficiency estimate's standard deviation
   SetModelParameters(
         x  ,       //   Int_t x,
         y  ,       //   Int_t y,
         0  ,       //   Int_t z,
         0  ,       //   Double_t bm,
         em ,       //   Double_t em,
         0  ,       //   Double_t e,
         2  ,       //   Int_t mid,
         sde,       //   Double_t sde,
         0  ,       //   Double_t sdb,
         tau,       //   Double_t tau,
         0  ,       //   Double_t b,
         0);        //   Int_t m

}

void TRolke::SetGaussBkgGaussEff(Int_t x, Double_t bm, Double_t em, Double_t sde, Double_t sdb)
{
// Model 3: Background - Gaussian, Efficiency - Gaussian (x,bm,em,sde,sdb)
//    x   : number of observed events in the experiment
//    bm  : estimate of the background
//    em  : estimate of the efficiency
//    sde : efficiency estimate's standard deviation
//    sdb : background estimate's standard deviation
   SetModelParameters(
         x  ,       //   Int_t x,
         0  ,       //   Int_t y,
         0  ,       //   Int_t z,
         bm ,       //   Double_t bm,
         em ,       //   Double_t em,
         0  ,       //   Double_t e,
         3  ,       //   Int_t mid,
         sde,       //   Double_t sde,
         sdb,       //   Double_t sdb,
         0  ,       //   Double_t tau,
         0  ,       //   Double_t b,
         0);        //   Int_t m

}

void TRolke::SetPoissonBkgKnownEff(Int_t x, Int_t y, Double_t tau, Double_t e)
{
// Model 4: Background - Poisson, Efficiency - known     (x,y,tau,e)
//    x   : number of observed events in the experiment
//    y   : number of observed events in background region
//    tau : ratio parameter (read TRolke.cxx for details)
//    e   : true efficiency (considered known)
   SetModelParameters(
         x  ,       //   Int_t x,
         y  ,       //   Int_t y,
         0  ,       //   Int_t z,
         0  ,       //   Double_t bm,
         0  ,       //   Double_t em,
         e  ,       //   Double_t e,
         4  ,       //   Int_t mid,
         0  ,       //   Double_t sde,
         0  ,       //   Double_t sdb,
         tau,       //   Double_t tau,
         0  ,       //   Double_t b,
         0);        //   Int_t m

}

void TRolke::SetGaussBkgKnownEff(Int_t x, Double_t bm, Double_t sdb, Double_t e)
{
// Model 5: Background - Gaussian, Efficiency - known    (x,bm,sdb,e
//    x   : number of observed events in the experiment
//    bm  : estimate of the background
//    sdb : background estimate's standard deviation
//    e   : true efficiency (considered known)
   SetModelParameters(
         x  ,       //   Int_t x,
         0  ,       //   Int_t y,
         0  ,       //   Int_t z,
         bm ,       //   Double_t bm,
         0  ,       //   Double_t em,
         e  ,       //   Double_t e,
         5  ,       //   Int_t mid,
         0  ,       //   Double_t sde,
         sdb,       //   Double_t sdb,
         0  ,       //   Double_t tau,
         0  ,       //   Double_t b,
         0);        //   Int_t m

}

void TRolke::SetKnownBkgBinomEff(Int_t x, Int_t z, Int_t m, Double_t b)
{
// Model 6: Background - known, Efficiency - Binomial    (x,z,m,b)
//    x   : number of observed events in the experiment
//    z   : number of MC events observed
//    m   : number of MC events generated
//    b   : background expectation value (considered known)
   SetModelParameters(
         x  ,       //   Int_t x,
         0  ,       //   Int_t y
         z  ,       //   Int_t z,
         0  ,       //   Double_t bm,
         0  ,       //   Double_t em,
         0  ,       //   Double_t e,
         6  ,       //   Int_t mid,
         0  ,       //   Double_t sde,
         0  ,       //   Double_t sdb,
         0  ,       //   Double_t tau,
         b  ,       //   Double_t b,
         m);        //   Int_t m

}

void TRolke::SetKnownBkgGaussEff(Int_t x, Double_t em, Double_t sde, Double_t b)
{
// Model 7: Background - known, Efficiency - Gaussian    (x,em,sde,b)
//    x   : number of observed events in the experiment
//    em  : estimate of the efficiency
//    sde : efficiency estimate's standard deviation
//    b   : background expectation value (considered known)
   SetModelParameters(
         x  ,       //   Int_t x,
         0  ,       //   Int_t y
         0  ,       //   Int_t z,
         0  ,       //   Double_t bm,
         em ,       //   Double_t em,
         0  ,       //   Double_t e,
         7  ,       //   Int_t mid,
         sde,       //   Double_t sde,
         0  ,       //   Double_t sdb,
         0  ,       //   Double_t tau,
         b  ,       //   Double_t b,
         0);        //   Int_t m

}

bool TRolke::GetLimits(Double_t& low, Double_t& high)
{
/* Calculate and get the upper and lower limits for the pre-specified model */
   if ((f_mid<1)||(f_mid>7)) {
      std::cerr << "TRolke - Error: Model id "<< f_mid<<std::endl;
      if (f_mid<1) {
         std::cerr << "TRolke - Please specify a model with e.g. 'SetGaussBkgGaussEff' (read the docs in Rolke.cxx )"<<std::endl;
      }
      return false;
   }

   ComputeInterval(f_x, f_y, f_z, f_bm, f_em, f_e, f_mid, f_sde, f_sdb, f_tau, f_b, f_m);
   low = fLowerLimit;
   high = fUpperLimit;
   if (low < high) {
      return true;
   }else{
      std::cerr << "TRolke - Warning: no limits found" <<std::endl;
      return false;
   }
}

Double_t TRolke::GetUpperLimit()
{
/* Calculate and get upper limit for the pre-specified model */
   Double_t low(0), high(0);
   GetLimits(low,high);
   return fUpperLimit;
}

Double_t TRolke::GetLowerLimit()
{
/* Calculate and get lower limit for the pre-specified model */
   Double_t low(0), high(0);
   GetLimits(low,high);
   return fLowerLimit;
}

Double_t TRolke::GetBackground()
{
/* Return a simple background value (estimate/truth) given the pre-specified model */
   Double_t background = 0;
   switch (f_mid) {
      case 1:
      case 2:
      case 4:
         background = f_y / f_tau;
         break;
      case 3:
      case 5:
         background = f_bm;
         break;
      case 6:
      case 7:
         background = f_b;
         break;
      default:
         std::cerr << "TRolke::GetBackground(): Model NR: " <<
         f_mid << " unknown"<<std::endl;
         return 0;
   }
   return background;
}

bool TRolke::GetSensitivity(Double_t& low, Double_t& high, Double_t pPrecision)
{
// get the upper and lower average limits based on the specified model.
// No uncertainties are considered for the Poisson weights in the averaging sum
   Double_t background = GetBackground();

   Double_t weight = 0;
   Double_t weightSum = 0;

   int loop_x = 0;

   while (1) {
      ComputeInterval(loop_x, f_y, f_z, f_bm, f_em, f_e, f_mid, f_sde, f_sdb, f_tau, f_b, f_m);
      weight = TMath::PoissonI(loop_x, background);
      low += fLowerLimit * weight;
      high += fUpperLimit * weight;
      weightSum += weight;
      if (loop_x > (background + 1)) { // don't stop too early
         if ((weightSum > (1 - pPrecision)) || (weight < 1e-12)) break;
      }
      loop_x++;
   }
   low /= weightSum;
   high /= weightSum;

   return (low < high); // could also add more detailed test
}


bool TRolke::GetLimitsQuantile(Double_t& low, Double_t& high, Int_t& out_x, Double_t integral)
{
/* get the upper and lower limits for the outcome corresponding to
   a given quantile.
   For integral=0.5 this gives the median limits
   in repeated experiments. The returned out_x is the corresponding
   (e.g. median) value of x.
   No uncertainties are considered for the Poisson weights when calculating
   the Poisson integral */
   Double_t background = GetBackground();
   Double_t weight = 0;
   Double_t weightSum = 0;
   Int_t loop_x = 0;

   while (1) {
      weight = TMath::PoissonI(loop_x, background);
      weightSum += weight;
      if (weightSum >= integral) {
         break;
      }
      loop_x++;
   }

   out_x = loop_x;

   ComputeInterval(loop_x, f_y, f_z, f_bm, f_em, f_e, f_mid, f_sde, f_sdb, f_tau, f_b, f_m);
   low = fLowerLimit;
   high = fUpperLimit;
   return (low < high); // could also add more detailed test

}

bool TRolke::GetLimitsML(Double_t& low, Double_t& high, Int_t& out_x)
{
/* get the upper and lower limits for the most likely outcome.
   The returned out_x is the corresponding value of x
   No uncertainties are considered for the Poisson weights when finding ML */
   Double_t background = GetBackground();

   Int_t loop_x = 0; // this can be optimized if needed.
   Int_t loop_max = 1000 + (Int_t)background; //     --||--

   Double_t max = TMath::PoissonI(loop_x, background);
   while (loop_x <= loop_max) {
      if (TMath::PoissonI(loop_x + 1, background) < max) {
         break;
      }
      loop_x++;
      max = TMath::PoissonI(loop_x, background);
   }
   if (loop_x >= loop_max) {
      std::cout << "internal error finding maximum of distribution" << std::endl;
      return false;
   }

   out_x = loop_x;

   ComputeInterval(loop_x, f_y, f_z, f_bm, f_em, f_e, f_mid, f_sde, f_sdb, f_tau, f_b, f_m);
   low = fLowerLimit;
   high = fUpperLimit;
   return (low < high); // could also add more detailed test

}

bool TRolke::GetCriticalNumber(Int_t& ncrit, Int_t maxtry)
{
// get the value of x corresponding to rejection of the null hypothesis.
// This means a lower limit >0 with the pre-specified Confidence Level.
// Optionally give maxtry; the maximum value of x to try. Of not, or if
// maxtry<0 an automatic mode is used.
   Double_t background = GetBackground();

   int j = 0;
   int rolke_ncrit = -1;
   int maxj =maxtry ;
   if(maxtry<1){
     maxj = 1000 + (Int_t)background; // max value to try
   }
   for (j = 0;j < maxj;j++) {
      Int_t rolke_x = j;
      ComputeInterval(rolke_x, f_y, f_z, f_bm, f_em, f_e, f_mid, f_sde, f_sdb, f_tau, f_b, f_m);
      double rolke_ll = fLowerLimit;
      if (rolke_ll > 0) {
         rolke_ncrit = j;
         break;
      }
   }

   if (rolke_ncrit == -1) {
     std::cerr << "TRolke GetCriticalNumber : Error: problem finding rolke inverse. Specify a larger maxtry value. maxtry was: " << maxj << ". highest x considered was j "<< j<< std::endl;
      ncrit = -1;
      return false;
   } else {
      ncrit = rolke_ncrit;
      return true;
   }
}

void TRolke::SetSwitch(bool bnd) {
/* Deprecated name for SetBounding. */
   if(fNumWarningsDeprecated1<2){
      std::cerr << "*******************************************" <<std::endl;
      std::cerr << "TRolke - Warning: 'SetSwitch' is depricated and may be removed from future releases:" <<std::endl;
      std::cerr << " - Use 'SetBounding' instead "<<std::endl;
      std::cerr << "*******************************************" <<std::endl;
      fNumWarningsDeprecated1++;
   }
   SetBounding(bnd);
}

void TRolke::Print(Option_t*) const {
/* Dump internals. Print members. */
   std::cout << "*******************************************" <<std::endl;
   std::cout << "* TRolke::Print() - dump of internals:                " <<std::endl;
   std::cout << "*"<<std::endl;
   std::cout << "* model id, mid = "<<f_mid <<std::endl;
   std::cout << "*"<<std::endl;
   std::cout << "*             x = "<<f_x   <<std::endl;
   std::cout << "*            bm = "<<f_bm  <<std::endl;
   std::cout << "*            em = "<<f_em  <<std::endl;
   std::cout << "*           sde = "<<f_sde <<std::endl;
   std::cout << "*           sdb = "<<f_sdb <<std::endl;
   std::cout << "*             y = "<<f_y   <<std::endl;
   std::cout << "*           tau = "<<f_tau <<std::endl;
   std::cout << "*             e = "<<f_e   <<std::endl;
   std::cout << "*             b = "<<f_b   <<std::endl;
   std::cout << "*             m = "<<f_m   <<std::endl;
   std::cout << "*             z = "<<f_z   <<std::endl;
   std::cout << "*"<<std::endl;
   std::cout << "*            CL = "<<fCL <<std::endl;
   std::cout << "*      Bounding = "<<fBounding <<std::endl;
   std::cout << "*"<<std::endl;
   std::cout << "* calculated on demand only:"<<std::endl;
   std::cout << "*   fUpperLimit = "<<fUpperLimit<<std::endl;
   std::cout << "*   fLowerLimit = "<<fLowerLimit<<std::endl;
   std::cout << "*******************************************" <<std::endl;
}


Double_t TRolke::CalculateInterval(Int_t x, Int_t y, Int_t z, Double_t bm, Double_t em, Double_t e, Int_t mid, Double_t sde, Double_t sdb, Double_t tau, Double_t b, Int_t m){
/* Deprecated and error prone model selection interface.
   It's use is trongly discouraged. 'mid' is the model ID (1 to 7).
   This method is provided for backwards compatibility/developer use only. */
//    x   : number of observed events in the experiment
//    y   : number of observed events in background region
//    z   : number of MC events observed
//    bm  : estimate of the background
//    em  : estimate of the efficiency
//    e   : true efficiency (considered known)
//    mid : internal model id (really, you should not use this method at all)
//    sde : efficiency estimate's standard deviation
//    sdb : background estimate's standard deviation
//    tau : ratio parameter (read TRolke.cxx for details)
//    b   : background expectation value (considered known)
//    m   : number of MC events generated
   if (fNumWarningsDeprecated2<2 ) {
      std::cerr << "*******************************************" <<std::endl;
      std::cerr << "TRolke - Warning: 'CalculateInterval' is depricated and may be removed from future releases:" <<std::endl;
      std::cerr << " - Use e.g. 'SetGaussBkgGaussEff' and 'GetLimits' instead (read the docs in Rolke.cxx )"<<std::endl;
      std::cerr << "*******************************************" <<std::endl;
      fNumWarningsDeprecated2++;
   }
   SetModelParameters(
         x,
         y,
         z,
         bm,
         em,
         e,
         mid,
         sde,
         sdb,
         tau,
         b,
         m);
   return ComputeInterval(f_x, f_y, f_z, f_bm, f_em, f_e, f_mid, f_sde, f_sdb, f_tau, f_b, f_m);
}


// --------------- Private methods ----------------
void TRolke::SetModelParameters(Int_t x, Int_t y, Int_t z, Double_t bm, Double_t em, Double_t e, Int_t mid, Double_t sde, Double_t sdb, Double_t tau, Double_t b, Int_t m)
{
//___________________________________________________________________________
//    x   : number of observed events in the experiment
//    y   : number of observed events in background region
//    z   : number of MC events observed
//    bm  : estimate of the background
//    em  : estimate of the efficiency
//    e   : true efficiency (considered known)
//    mid : internal model id
//    sde : efficiency estimate's standard deviation
//    sdb : background estimate's standard deviation
//    tau : ratio parameter (read TRolke.cxx for details)
//    b   : background expectation value (considered known)
//    m   : number of MC events generated
   f_x   = x   ;
   f_y   = y   ;
   f_z   = z   ;
   f_bm  = bm  ;
   f_em  = em  ;
   f_e   = e   ;
   f_mid = mid ;
   f_sde = sde ;
   f_sdb = sdb ;
   f_tau = tau ;
   f_b   = b   ;
   f_m   = m   ;
}

void TRolke::SetModelParameters()
{
/* Clear internal model */
   SetModelParameters(0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0);
   f_mid=0;
}


Double_t TRolke::ComputeInterval(Int_t x, Int_t y, Int_t z, Double_t bm, Double_t em, Double_t e, Int_t mid, Double_t sde, Double_t sdb, Double_t tau, Double_t b, Int_t m)
{
//___________________________________________________________________________
// ComputeInterval, the internals.
//    x   : number of observed events in the experiment
//    y   : number of observed events in background region
//    z   : number of MC events observed
//    bm  : estimate of the background
//    em  : estimate of the efficiency
//    e   : true efficiency (considered known)
//    mid : internal model id (really, you should not use this method at all)
//    sde : efficiency estimate's standard deviation
//    sdb : background estimate's standard deviation
//    tau : ratio parameter (read TRolke.cxx for details)
//    b   : background expectation value (considered known)
//    m   : number of MC events generated

   //calculate interval
   Int_t done = 0;
   Double_t limit[2];

   limit[1] = Interval(x, y, z, bm, em, e, mid, sde, sdb, tau, b, m);

   if (limit[1] > 0) {
      done = 1;
   }

   if (! fBounding) {

      Int_t trial_x = x;

      while (done == 0) {
         trial_x++;
         limit[1] = Interval(trial_x, y, z, bm, em, e, mid, sde, sdb, tau, b, m);
         if (limit[1] > 0) done = 1;
      }
   }

   return limit[1];
}

Double_t TRolke::Interval(Int_t x, Int_t y, Int_t z, Double_t bm, Double_t em, Double_t e, Int_t mid, Double_t sde, Double_t sdb, Double_t tau, Double_t b, Int_t m)
{
//_____________________________________________________________________
// Internal helper function 'Interval'
//
//    x   : number of observed events in the experiment
//    y   : number of observed events in background region
//    z   : number of MC events observed
//    bm  : estimate of the background
//    em  : estimate of the efficiency
//    e   : true efficiency (considered known)
//    mid : internal model id (really, you should not use this method at all)
//    sde : efficiency estimate's standard deviation
//    sdb : background estimate's standard deviation
//    tau : ratio parameter (read TRolke.cxx for details)
//    b   : background expectation value (considered known)
//    m   : number of MC events generated

   Double_t dchi2 = TMath::ChisquareQuantile(fCL, 1);
   Double_t tempxy[2], limits[2] = {0, 0};
   Double_t slope, fmid, low, flow, high, fhigh, test, ftest, mu0, maximum, target, l, f0;
   Double_t med = 0;
   Double_t maxiter = 1000, acc = 0.00001;
   Int_t i;
   Int_t bp = 0;

   if ((mid != 3) && (mid != 5)) bm = y;
   if ((mid == 3) || (mid == 5)) {
      if (bm == 0) bm = 0.00001;
   }

   if ((mid == 6) || (mid == 7)) {
      if (bm == 0) bm = 0.00001;
   }

   if ((mid <= 2) || (mid == 4)) bp = 1;


   if (bp == 1 && x == 0 && bm > 0) {
      for (i = 0; i < 2; i++) {
         x++;
         tempxy[i] = Interval(x, y, z, bm, em, e, mid, sde, sdb, tau, b, m);
      }

      slope = tempxy[1] - tempxy[0];
      limits[1] = tempxy[0] - slope;
      limits[0] = 0.0;
      if (limits[1] < 0) limits[1] = 0.0;
      goto done;
   }

   if (bp != 1 && x == 0) {

      for (i = 0; i < 2; i++) {
         x++;
         tempxy[i] = Interval(x, y, z, bm, em, e, mid, sde, sdb, tau, b, m);
      }
      slope = tempxy[1] - tempxy[0];
      limits[1] = tempxy[0] - slope;
      limits[0] = 0.0;
      if (limits[1] < 0) limits[1] = 0.0;
      goto done;
   }

   if (bp != 1  && bm == 0) {
      for (i = 0; i < 2; i++) {
         bm++;
         limits[1] = Interval(x, y, z, bm, em, e, mid, sde, sdb, tau, b, m);
         tempxy[i] = limits[1];
      }
      slope = tempxy[1] - tempxy[0];
      limits[1] = tempxy[0] - slope;
      if (limits[1] < 0) limits[1] = 0;
      goto done;
   }

   if (x == 0 && bm == 0) {
      x++;
      bm++;
      limits[1] = Interval(x, y, z, bm, em, e, mid, sde, sdb, tau, b, m);
      tempxy[0] = limits[1];
      x  = 1;
      bm = 2;
      limits[1] = Interval(x, y, z, bm, em, e, mid, sde, sdb, tau, b, m);
      tempxy[1] = limits[1];
      x  = 2;
      bm = 1;
      limits[1] = Interval(x, y, z, bm, em, e, mid, sde, sdb, tau, b, m);
      limits[1] = 3 * tempxy[0] - tempxy[1] - limits[1];
      if (limits[1] < 0) limits[1] = 0;
      goto done;
   }

   mu0 = Likelihood(0, x, y, z, bm, em, mid, sde, sdb, tau, b, m, 1);
   maximum = Likelihood(0, x, y, z, bm, em, mid, sde, sdb, tau, b, m, 2);
   test = 0;
   f0 = Likelihood(test, x, y, z, bm, em, mid, sde, sdb, tau, b, m, 3);
   if (fBounding) {
      if (mu0 < 0) maximum = f0;
   }

   target = maximum - dchi2;
   if (f0 > target) {
      limits[0] = 0;
   } else {
      if (mu0 < 0) {
         limits[0] = 0;
         limits[1] = 0;
      }

      low   = 0;
      flow  = f0;
      high  = mu0;
      fhigh = maximum;
      for (i = 0; i < maxiter; i++) {
         l = (target - fhigh) / (flow - fhigh);
         if (l < 0.2) l = 0.2;
         if (l > 0.8) l = 0.8;
         med = l * low + (1 - l) * high;
         if (med < 0.01) {
            limits[1] = 0.0;
            goto done;
         }
         fmid = Likelihood(med, x, y, z, bm, em, mid, sde, sdb, tau, b, m, 3);
         if (fmid > target) {
            high  = med;
            fhigh = fmid;
         } else {
            low  = med;
            flow = fmid;
         }
         if ((high - low) < acc*high) break;
      }
      limits[0] = med;
   }

   if (mu0 > 0) {
      low  = mu0;
      flow = maximum;
   } else {
      low  = 0;
      flow = f0;
   }

   test = low + 1 ;
   ftest = Likelihood(test, x, y, z, bm, em, mid, sde, sdb, tau, b, m, 3);
   if (ftest < target) {
      high  = test;
      fhigh = ftest;
   } else {
      slope = (ftest - flow) / (test - low);
      high  = test + (target - ftest) / slope;
      fhigh = Likelihood(high, x, y, z, bm, em, mid, sde, sdb, tau, b, m, 3);
   }

   for (i = 0; i < maxiter; i++) {
      l = (target - fhigh) / (flow - fhigh);
      if (l < 0.2) l = 0.2;
      if (l > 0.8) l = 0.8;
      med  = l * low + (1. - l) * high;
      fmid = Likelihood(med, x, y, z, bm, em, mid, sde, sdb, tau, b, m, 3);

      if (fmid < target) {
         high  = med;
         fhigh = fmid;
      } else {
         low  = med;
         flow = fmid;
      }

      if (high - low < acc*high) break;
   }

   limits[1] = med;

done:

   // apply known efficiency
   if ((mid == 4) || (mid == 5)) {
      limits[0] /= e;
      limits[1] /= e;
   }

   fUpperLimit = limits[1];
   fLowerLimit = TMath::Max(limits[0], 0.0);

   return limits[1];
}


Double_t TRolke::Likelihood(Double_t mu, Int_t x, Int_t y, Int_t z, Double_t bm, Double_t em, Int_t mid, Double_t sde, Double_t sdb, Double_t tau, Double_t b, Int_t m, Int_t what)
{
//___________________________________________________________________________
//Internal helper function
// Chooses between the different profile likelihood functions to use for the
// different models.
// Returns evaluation of the profile likelihood functions.

   switch (mid) {
      case 1:
         return EvalLikeMod1(mu, x, y, z, tau, m, what);
      case 2:
         return EvalLikeMod2(mu, x, y, em, sde, tau, what);
      case 3:
         return EvalLikeMod3(mu, x, bm, em, sde, sdb, what);
      case 4:
         return EvalLikeMod4(mu, x, y, tau, what);
      case 5:
         return EvalLikeMod5(mu, x, bm, sdb, what);
      case 6:
         return EvalLikeMod6(mu, x, z, b, m, what);
      case 7:
         return EvalLikeMod7(mu, x, em, sde, b, what);
      default:
         std::cerr << "TRolke::Likelihood(...): Model NR: " <<
         f_mid << " unknown"<<std::endl;
         return 0;
   }

   return 0;
}

Double_t TRolke::EvalLikeMod1(Double_t mu, Int_t x, Int_t y, Int_t z, Double_t tau, Int_t m, Int_t what)
{
//_________________________________________________________________________
//
// Calculates the Profile Likelihood for MODEL 1:
//  Poisson background/ Binomial Efficiency
// what = 1: Maximum likelihood estimate is returned
// what = 2: Profile Likelihood of Maximum Likelihood estimate is returned.
// what = 3: Profile Likelihood of Test hypothesis is returned
// otherwise parameters as described in the beginning of the class)
   Double_t f  = 0;
   Double_t zm = Double_t(z) / m;

   if (what == 1) {
      f = (x - y / tau) / zm;
   }

   if (what == 2) {
      mu = (x - y / tau) / zm;
      Double_t b  = y / tau;
      Double_t e = zm;
      f = LikeMod1(mu, b, e, x, y, z, tau, m);
   }

   if (what == 3) {
      if (mu == 0) {
         Double_t b = (x + y) / (1.0 + tau);
         Double_t e = zm;
         f = LikeMod1(mu, b, e, x, y, z, tau, m);
      } else {
         Double_t e = 0;
         Double_t b = 0;
         ProfLikeMod1(mu, b, e, x, y, z, tau, m);
         f = LikeMod1(mu, b, e, x, y, z, tau, m);
      }
   }

   return f;
}


Double_t TRolke::LikeMod1(Double_t mu, Double_t b, Double_t e, Int_t x, Int_t y, Int_t z, Double_t tau, Int_t m)
{
//________________________________________________________________________
// Profile Likelihood function for MODEL 1:
// Poisson background/ Binomial Efficiency

   double s = e*mu+b;
   double lls = - s;
   if (x > 0) lls = x*TMath::Log(s) - s - LogFactorial(x);
   double bg = tau*b;
   double llb =  -bg;
   if ( y > 0) llb =  y*TMath::Log( bg) - bg - LogFactorial(y);

   double lle = 0;  // binomial log-like
   if (z == 0)         lle = m * TMath::Log(1-e);
   else if ( z == m)   lle = m * TMath::Log(e);
   else                lle =   z * TMath::Log(e) + (m - z)*TMath::Log(1 - e) + LogFactorial(m) - LogFactorial(m-z) - LogFactorial(z);

   double f = 2*( lls + llb + lle);
   return f;
}


// this code is non-sense - // need to solve using Minuit
struct LikeFunction1 {
};

void TRolke::ProfLikeMod1(Double_t mu, Double_t &b, Double_t &e, Int_t x, Int_t y, Int_t z, Double_t tau, Int_t m)
{
//________________________________________________________________________
// Helper for calculation of estimates of efficiency and background for model 1
//

   Double_t med = 0.0, fmid;
   Int_t maxiter = 1000;
   Double_t acc = 0.00001;
   Double_t emin = ((m + mu * tau) - TMath::Sqrt((m + mu * tau) * (m + mu * tau) - 4 * mu * tau * z)) / 2 / mu / tau;

   Double_t low  = TMath::Max(1e-10, emin + 1e-10);
   Double_t high = 1 - 1e-10;

   for (Int_t i = 0; i < maxiter; i++) {
      med = (low + high) / 2.;

      fmid = LikeGradMod1(med, mu, x, y, z, tau, m);

      if (high < 0.5) acc = 0.00001 * high;
      else           acc = 0.00001 * (1 - high);

      if ((high - low) < acc*high) break;

      if (fmid > 0) low  = med;
      else         high = med;
   }

   e = med;
   Double_t eta = Double_t(z) / e - Double_t(m - z) / (1 - e);

   b = Double_t(y) / (tau - eta / mu);
}

//___________________________________________________________________________
Double_t TRolke::LikeGradMod1(Double_t e, Double_t mu, Int_t x, Int_t y, Int_t z, Double_t tau, Int_t m)
{
//gradient model likelihood
   Double_t eta, etaprime, bprime, f;
   eta = static_cast<double>(z) / e - static_cast<double>(m - z) / (1.0 - e);
   etaprime = (-1) * (static_cast<double>(m - z) / ((1.0 - e) * (1.0 - e)) + static_cast<double>(z) / (e * e));
   Double_t b = y / (tau - eta / mu);
   bprime = (b * b * etaprime) / mu / y;
   f = (mu + bprime) * (x / (e * mu + b) - 1) + (y / b - tau) * bprime + eta;
   return f;
}

//___________________________________________________________________________
Double_t TRolke::EvalLikeMod2(Double_t mu, Int_t x, Int_t y, Double_t em, Double_t sde, Double_t tau, Int_t what)
{
// Calculates the Profile Likelihood for MODEL 2:
//  Poisson background/ Gauss Efficiency
// what = 1: Maximum likelihood estimate is returned
// what = 2: Profile Likelihood of Maximum Likelihood estimate is returned.
// what = 3: Profile Likelihood of Test hypothesis is returned
// otherwise parameters as described in the beginning of the class)
   Double_t v =  sde * sde;
   Double_t coef[4], roots[3];
   Double_t f = 0;

   if (what == 1) {
      f = (x - y / tau) / em;
   }

   if (what == 2) {
      mu = (x - y / tau) / em;
      Double_t b = y / tau;
      Double_t e = em;
      f = LikeMod2(mu, b, e, x, y, em, tau, v);
   }

   if (what == 3) {
      if (mu == 0) {
         Double_t b = (x + y) / (1 + tau);
         Double_t e = em ;
         f = LikeMod2(mu, b, e, x, y, em, tau, v);
      } else {
         coef[3] = mu;
         coef[2] = mu * mu * v - 2 * em * mu - mu * mu * v * tau;
         coef[1] = (- x) * mu * v - mu * mu * mu * v * v * tau - mu * mu * v * em + em * mu * mu * v * tau + em * em * mu - y * mu * v;
         coef[0] = x * mu * mu * v * v * tau + x * em * mu * v - y * mu * mu * v * v + y * em * mu * v;

         TMath::RootsCubic(coef, roots[0], roots[1], roots[2]);

         Double_t e = roots[1];
         Double_t b;
         if ( v > 0) b = y / (tau + (em - e) / mu / v);
         else b = y/tau;
         f = LikeMod2(mu, b, e, x, y, em, tau, v);
      }
   }

   return f;
}

//_________________________________________________________________________
Double_t TRolke::LikeMod2(Double_t mu, Double_t b, Double_t e, Int_t x, Int_t y, Double_t em, Double_t tau, Double_t v)
{
// Profile Likelihood function for MODEL 2:
// Poisson background/Gauss Efficiency

   double s = e*mu+b;
   double lls = - s;
   if (x > 0) lls = x*TMath::Log(s) - s - LogFactorial(x);
   double bg = tau*b;
   double llb =  -bg;
   if ( y > 0) llb =  y*TMath::Log( bg) - bg - LogFactorial(y);
   double lle = 0;
   if ( v > 0) lle = - 0.9189385 - TMath::Log(v) / 2 - (em - e)*(em - e) / v / 2;

   return 2*( lls + llb + lle);
}

//_____________________________________________________________________
Double_t TRolke::EvalLikeMod3(Double_t mu, Int_t x, Double_t bm, Double_t em, Double_t sde, Double_t sdb, Int_t what)
{
// Calculates the Profile Likelihood for MODEL 3:
// Gauss  background/ Gauss Efficiency
// what = 1: Maximum likelihood estimate is returned
// what = 2: Profile Likelihood of Maximum Likelihood estimate is returned.
// what = 3: Profile Likelihood of Test hypothesis is returned
// otherwise parameters as described in the beginning of the class)

   Double_t f = 0.;
   Double_t  v = sde * sde;
   Double_t  u = sdb * sdb;

   if (what == 1) {
      f = (x - bm) / em;
   }


   if (what == 2) {
      mu = (x - bm) / em;
      Double_t b  = bm;
      Double_t e  = em;
      f  = LikeMod3(mu, b, e, x, bm, em, u, v);
   }


   if (what == 3) {
      if (mu == 0.0) {
         Double_t b = ((bm - u) + TMath::Sqrt((bm - u) * (bm - u) + 4 * x * u)) / 2.;
         Double_t e = em;
         f = LikeMod3(mu, b, e, x, bm, em, u, v);
      } else {
         Double_t e = em;
         Double_t b = bm;
         if ( v > 0) {
            Double_t temp[3];
            temp[0] = mu * mu * v + u;
            temp[1] = mu * mu * mu * v * v + mu * v * u - mu * mu * v * em + mu * v * bm - 2 * u * em;
            temp[2] = mu * mu * v * v * bm - mu * v * u * em - mu * v * bm * em + u * em * em - mu * mu * v * v * x;
            e = (-temp[1] + TMath::Sqrt(temp[1] * temp[1] - 4 * temp[0] * temp[2])) / 2 / temp[0];
            b = bm - (u * (em - e)) / v / mu;
         }
         f = LikeMod3(mu, b, e, x, bm, em, u, v);
      }
   }

   return f;
}

//____________________________________________________________________
Double_t TRolke::LikeMod3(Double_t mu, Double_t b, Double_t e, Int_t x, Double_t bm, Double_t em, Double_t u, Double_t v)
{
// Profile Likelihood function for MODEL 3:
// Gauss background/Gauss Efficiency

   double s = e*mu+b;
   double lls = - s;
   if (x > 0) lls = x*TMath::Log(s) - s - LogFactorial(x);
   double llb =  0;
   if ( u > 0) llb = - 0.9189385 - TMath::Log(u) / 2 - (bm - b)*(bm - b) / u / 2;
   double lle = 0;
   if ( v > 0) lle = - 0.9189385 - TMath::Log(v) / 2 - (em - e)*(em - e) / v / 2;

   return 2*( lls + llb + lle);

}

//____________________________________________________________________
Double_t TRolke::EvalLikeMod4(Double_t mu, Int_t x, Int_t y, Double_t tau, Int_t what)
{
// Calculates the Profile Likelihood for MODEL 4:
// Poiss  background/Efficiency known
// what = 1: Maximum likelihood estimate is returned
// what = 2: Profile Likelihood of Maximum Likelihood estimate is returned.
// what = 3: Profile Likelihood of Test hypothesis is returned
// otherwise parameters as described in the beginning of the class)
   Double_t f = 0.0;

   if (what == 1) f = x - y / tau;
   if (what == 2) {
      mu = x - y / tau;
      Double_t b  = y / tau;
      f  = LikeMod4(mu, b, x, y, tau);
   }
   if (what == 3) {
      if (mu == 0.0) {
         Double_t b = Double_t(x + y) / (1 + tau);
         f = LikeMod4(mu, b, x, y, tau);
      } else {
         Double_t b = (x + y - (1 + tau) * mu + sqrt((x + y - (1 + tau) * mu) * (x + y - (1 + tau) * mu) + 4 * (1 + tau) * y * mu)) / 2 / (1 + tau);
         f = LikeMod4(mu, b, x, y, tau);
      }
   }
   return f;
}

//___________________________________________________________________
Double_t TRolke::LikeMod4(Double_t mu, Double_t b, Int_t x, Int_t y, Double_t tau)
{
// Profile Likelihood function for MODEL 4:
// Poiss background/Efficiency known

   double s = mu+b;
   double lls = - s;
   if (x > 0) lls = x*TMath::Log(s) - s - LogFactorial(x);
   double bg = tau*b;
   double llb =  -bg;
   if ( y > 0) llb =  y*TMath::Log( bg) - bg - LogFactorial(y);

   return 2*( lls + llb);
}

//___________________________________________________________________
Double_t TRolke::EvalLikeMod5(Double_t mu, Int_t x, Double_t bm, Double_t sdb, Int_t what)
{
// Calculates the Profile Likelihood for MODEL 5:
// Gauss  background/Efficiency known
// what = 1: Maximum likelihood estimate is returned
// what = 2: Profile Likelihood of Maximum Likelihood estimate is returned.
// what = 3: Profile Likelihood of Test hypothesis is returned
// otherwise parameters as described in the beginning of the class)

   Double_t u = sdb * sdb;
   Double_t f = 0;

   if (what == 1) {
      f = x - bm;
   }
   if (what == 2) {
      mu = x - bm;
      Double_t b  = bm;
      f  = LikeMod5(mu, b, x, bm, u);
   }

   if (what == 3) {
      Double_t b = ((bm - u - mu) + TMath::Sqrt((bm - u - mu) * (bm - u - mu) - 4 * (mu * u - mu * bm - u * x))) / 2;
      f = LikeMod5(mu, b, x, bm, u);
   }
   return f;
}

//_______________________________________________________________________
Double_t TRolke::LikeMod5(Double_t mu, Double_t b, Int_t x, Double_t bm, Double_t u)
{
// Profile Likelihood function for MODEL 5:
// Gauss background/Efficiency known

   double s = mu+b;
   double lls = - s;
   if (x > 0) lls = x*TMath::Log(s) - s - LogFactorial(x);
   double llb =  0;
   if ( u > 0) llb = - 0.9189385 - TMath::Log(u) / 2 - (bm - b)*(bm - b) / u / 2;

   return 2*( lls + llb);
}

//_______________________________________________________________________
Double_t TRolke::EvalLikeMod6(Double_t mu, Int_t x, Int_t z, Double_t b, Int_t m, Int_t what)
{
// Calculates the Profile Likelihood for MODEL 6:
// Background known/Efficiency binomial
// what = 1: Maximum likelihood estimate is returned
// what = 2: Profile Likelihood of Maximum Likelihood estimate is returned.
// what = 3: Profile Likelihood of Test hypothesis is returned
// otherwise parameters as described in the beginning of the class)

   Double_t coef[4], roots[3];
   Double_t f = 0.;
   Double_t zm = Double_t(z) / m;

   if (what == 1) {
      f = (x - b) / zm;
   }

   if (what == 2) {
      mu = (x - b) / zm;
      Double_t e  = zm;
      f  = LikeMod6(mu, b, e, x, z, m);
   }
   if (what == 3) {
      Double_t e;
      if (mu == 0) {
         e = zm;
      } else {
         coef[3] = mu * mu;
         coef[2] = mu * b - mu * x - mu * mu - mu * m;
         coef[1] = mu * x - mu * b + mu * z - m * b;
         coef[0] = b * z;
         TMath::RootsCubic(coef, roots[0], roots[1], roots[2]);
         e = roots[1];
      }
      f = LikeMod6(mu, b, e, x, z, m);
   }
   return f;
}

//_______________________________________________________________________
Double_t TRolke::LikeMod6(Double_t mu, Double_t b, Double_t e, Int_t x, Int_t z, Int_t m)
{
// Profile Likelihood function for MODEL 6:
// background known/ Efficiency binomial

   double s = e*mu+b;
   double lls = - s;
   if (x > 0) lls = x*TMath::Log(s) - s - LogFactorial(x);

   double lle = 0;
   if (z == 0)        lle = m * TMath::Log(1-e);
   else if ( z == m)  lle = m * TMath::Log(e);
   else               lle =   z * TMath::Log(e) + (m - z)*TMath::Log(1 - e) + LogFactorial(m) - LogFactorial(m-z) - LogFactorial(z);

   return 2* (lls + lle);
}


//___________________________________________________________________________
Double_t TRolke::EvalLikeMod7(Double_t mu, Int_t x, Double_t em, Double_t sde, Double_t b, Int_t what)
{
// Calculates the Profile Likelihood for MODEL 7:
// background known/Efficiency Gauss
// what = 1: Maximum likelihood estimate is returned
// what = 2: Profile Likelihood of Maximum Likelihood estimate is returned.
// what = 3: Profile Likelihood of Test hypothesis is returned
// otherwise parameters as described in the beginning of the class)

   Double_t v = sde * sde;
   Double_t f = 0.;

   if (what ==  1) {
      f = (x - b) / em;
   }

   if (what == 2) {
      mu = (x - b) / em;
      Double_t e  = em;
      f  = LikeMod7(mu, b, e, x, em, v);
   }

   if (what == 3) {
      Double_t e;
      if (mu == 0) {
         e = em;
      } else {
         e = (-(mu * em - b - mu * mu * v) - TMath::Sqrt((mu * em - b - mu * mu * v) * (mu * em - b - mu * mu * v) + 4 * mu * (x * mu * v - mu * b * v + b * em))) / (- mu) / 2;
      }
      f = LikeMod7(mu, b, e, x, em, v);
   }

   return f;
}

//___________________________________________________________________________
Double_t TRolke::LikeMod7(Double_t mu, Double_t b, Double_t e, Int_t x, Double_t em, Double_t v)
{
// Profile Likelihood function for MODEL 6:
// background known/ Efficiency gaussian

   double s = e*mu+b;
   double lls = - s;
   if (x > 0) lls = x*TMath::Log(s) - s - LogFactorial(x);

   double lle = 0;
   if ( v > 0) lle = - 0.9189385 - TMath::Log(v) / 2 - (em - e)*(em - e) / v / 2;

   return 2*( lls + lle);
}

//______________________________________________________________________
Double_t TRolke::EvalPolynomial(Double_t x, const Int_t  coef[], Int_t N)
{
// evaluate polynomial
   const Int_t   *p;
   p = coef;
   Double_t ans = *p++;
   Int_t i = N;

   do
      ans = ans * x  +  *p++;
   while (--i);

   return ans;
}

//______________________________________________________________________
Double_t TRolke::EvalMonomial(Double_t x, const Int_t coef[], Int_t N)
{
// evaluate mononomial
   Double_t ans;
   const Int_t   *p;

   p   = coef;
   ans = x + *p++;
   Int_t i = N - 1;

   do
      ans = ans * x  + *p++;
   while (--i);

   return ans;
}

//______________________________________________________________________
Double_t TRolke::LogFactorial(Int_t n)
{
// LogFactorial function (use the logGamma function via the relation Gamma(n+1) = n!

   return TMath::LnGamma(n+1);
}


 TRolke.cxx:1
 TRolke.cxx:2
 TRolke.cxx:3
 TRolke.cxx:4
 TRolke.cxx:5
 TRolke.cxx:6
 TRolke.cxx:7
 TRolke.cxx:8
 TRolke.cxx:9
 TRolke.cxx:10
 TRolke.cxx:11
 TRolke.cxx:12
 TRolke.cxx:13
 TRolke.cxx:14
 TRolke.cxx:15
 TRolke.cxx:16
 TRolke.cxx:17
 TRolke.cxx:18
 TRolke.cxx:19
 TRolke.cxx:20
 TRolke.cxx:21
 TRolke.cxx:22
 TRolke.cxx:23
 TRolke.cxx:24
 TRolke.cxx:25
 TRolke.cxx:26
 TRolke.cxx:27
 TRolke.cxx:28
 TRolke.cxx:29
 TRolke.cxx:30
 TRolke.cxx:31
 TRolke.cxx:32
 TRolke.cxx:33
 TRolke.cxx:34
 TRolke.cxx:35
 TRolke.cxx:36
 TRolke.cxx:37
 TRolke.cxx:38
 TRolke.cxx:39
 TRolke.cxx:40
 TRolke.cxx:41
 TRolke.cxx:42
 TRolke.cxx:43
 TRolke.cxx:44
 TRolke.cxx:45
 TRolke.cxx:46
 TRolke.cxx:47
 TRolke.cxx:48
 TRolke.cxx:49
 TRolke.cxx:50
 TRolke.cxx:51
 TRolke.cxx:52
 TRolke.cxx:53
 TRolke.cxx:54
 TRolke.cxx:55
 TRolke.cxx:56
 TRolke.cxx:57
 TRolke.cxx:58
 TRolke.cxx:59
 TRolke.cxx:60
 TRolke.cxx:61
 TRolke.cxx:62
 TRolke.cxx:63
 TRolke.cxx:64
 TRolke.cxx:65
 TRolke.cxx:66
 TRolke.cxx:67
 TRolke.cxx:68
 TRolke.cxx:69
 TRolke.cxx:70
 TRolke.cxx:71
 TRolke.cxx:72
 TRolke.cxx:73
 TRolke.cxx:74
 TRolke.cxx:75
 TRolke.cxx:76
 TRolke.cxx:77
 TRolke.cxx:78
 TRolke.cxx:79
 TRolke.cxx:80
 TRolke.cxx:81
 TRolke.cxx:82
 TRolke.cxx:83
 TRolke.cxx:84
 TRolke.cxx:85
 TRolke.cxx:86
 TRolke.cxx:87
 TRolke.cxx:88
 TRolke.cxx:89
 TRolke.cxx:90
 TRolke.cxx:91
 TRolke.cxx:92
 TRolke.cxx:93
 TRolke.cxx:94
 TRolke.cxx:95
 TRolke.cxx:96
 TRolke.cxx:97
 TRolke.cxx:98
 TRolke.cxx:99
 TRolke.cxx:100
 TRolke.cxx:101
 TRolke.cxx:102
 TRolke.cxx:103
 TRolke.cxx:104
 TRolke.cxx:105
 TRolke.cxx:106
 TRolke.cxx:107
 TRolke.cxx:108
 TRolke.cxx:109
 TRolke.cxx:110
 TRolke.cxx:111
 TRolke.cxx:112
 TRolke.cxx:113
 TRolke.cxx:114
 TRolke.cxx:115
 TRolke.cxx:116
 TRolke.cxx:117
 TRolke.cxx:118
 TRolke.cxx:119
 TRolke.cxx:120
 TRolke.cxx:121
 TRolke.cxx:122
 TRolke.cxx:123
 TRolke.cxx:124
 TRolke.cxx:125
 TRolke.cxx:126
 TRolke.cxx:127
 TRolke.cxx:128
 TRolke.cxx:129
 TRolke.cxx:130
 TRolke.cxx:131
 TRolke.cxx:132
 TRolke.cxx:133
 TRolke.cxx:134
 TRolke.cxx:135
 TRolke.cxx:136
 TRolke.cxx:137
 TRolke.cxx:138
 TRolke.cxx:139
 TRolke.cxx:140
 TRolke.cxx:141
 TRolke.cxx:142
 TRolke.cxx:143
 TRolke.cxx:144
 TRolke.cxx:145
 TRolke.cxx:146
 TRolke.cxx:147
 TRolke.cxx:148
 TRolke.cxx:149
 TRolke.cxx:150
 TRolke.cxx:151
 TRolke.cxx:152
 TRolke.cxx:153
 TRolke.cxx:154
 TRolke.cxx:155
 TRolke.cxx:156
 TRolke.cxx:157
 TRolke.cxx:158
 TRolke.cxx:159
 TRolke.cxx:160
 TRolke.cxx:161
 TRolke.cxx:162
 TRolke.cxx:163
 TRolke.cxx:164
 TRolke.cxx:165
 TRolke.cxx:166
 TRolke.cxx:167
 TRolke.cxx:168
 TRolke.cxx:169
 TRolke.cxx:170
 TRolke.cxx:171
 TRolke.cxx:172
 TRolke.cxx:173
 TRolke.cxx:174
 TRolke.cxx:175
 TRolke.cxx:176
 TRolke.cxx:177
 TRolke.cxx:178
 TRolke.cxx:179
 TRolke.cxx:180
 TRolke.cxx:181
 TRolke.cxx:182
 TRolke.cxx:183
 TRolke.cxx:184
 TRolke.cxx:185
 TRolke.cxx:186
 TRolke.cxx:187
 TRolke.cxx:188
 TRolke.cxx:189
 TRolke.cxx:190
 TRolke.cxx:191
 TRolke.cxx:192
 TRolke.cxx:193
 TRolke.cxx:194
 TRolke.cxx:195
 TRolke.cxx:196
 TRolke.cxx:197
 TRolke.cxx:198
 TRolke.cxx:199
 TRolke.cxx:200
 TRolke.cxx:201
 TRolke.cxx:202
 TRolke.cxx:203
 TRolke.cxx:204
 TRolke.cxx:205
 TRolke.cxx:206
 TRolke.cxx:207
 TRolke.cxx:208
 TRolke.cxx:209
 TRolke.cxx:210
 TRolke.cxx:211
 TRolke.cxx:212
 TRolke.cxx:213
 TRolke.cxx:214
 TRolke.cxx:215
 TRolke.cxx:216
 TRolke.cxx:217
 TRolke.cxx:218
 TRolke.cxx:219
 TRolke.cxx:220
 TRolke.cxx:221
 TRolke.cxx:222
 TRolke.cxx:223
 TRolke.cxx:224
 TRolke.cxx:225
 TRolke.cxx:226
 TRolke.cxx:227
 TRolke.cxx:228
 TRolke.cxx:229
 TRolke.cxx:230
 TRolke.cxx:231
 TRolke.cxx:232
 TRolke.cxx:233
 TRolke.cxx:234
 TRolke.cxx:235
 TRolke.cxx:236
 TRolke.cxx:237
 TRolke.cxx:238
 TRolke.cxx:239
 TRolke.cxx:240
 TRolke.cxx:241
 TRolke.cxx:242
 TRolke.cxx:243
 TRolke.cxx:244
 TRolke.cxx:245
 TRolke.cxx:246
 TRolke.cxx:247
 TRolke.cxx:248
 TRolke.cxx:249
 TRolke.cxx:250
 TRolke.cxx:251
 TRolke.cxx:252
 TRolke.cxx:253
 TRolke.cxx:254
 TRolke.cxx:255
 TRolke.cxx:256
 TRolke.cxx:257
 TRolke.cxx:258
 TRolke.cxx:259
 TRolke.cxx:260
 TRolke.cxx:261
 TRolke.cxx:262
 TRolke.cxx:263
 TRolke.cxx:264
 TRolke.cxx:265
 TRolke.cxx:266
 TRolke.cxx:267
 TRolke.cxx:268
 TRolke.cxx:269
 TRolke.cxx:270
 TRolke.cxx:271
 TRolke.cxx:272
 TRolke.cxx:273
 TRolke.cxx:274
 TRolke.cxx:275
 TRolke.cxx:276
 TRolke.cxx:277
 TRolke.cxx:278
 TRolke.cxx:279
 TRolke.cxx:280
 TRolke.cxx:281
 TRolke.cxx:282
 TRolke.cxx:283
 TRolke.cxx:284
 TRolke.cxx:285
 TRolke.cxx:286
 TRolke.cxx:287
 TRolke.cxx:288
 TRolke.cxx:289
 TRolke.cxx:290
 TRolke.cxx:291
 TRolke.cxx:292
 TRolke.cxx:293
 TRolke.cxx:294
 TRolke.cxx:295
 TRolke.cxx:296
 TRolke.cxx:297
 TRolke.cxx:298
 TRolke.cxx:299
 TRolke.cxx:300
 TRolke.cxx:301
 TRolke.cxx:302
 TRolke.cxx:303
 TRolke.cxx:304
 TRolke.cxx:305
 TRolke.cxx:306
 TRolke.cxx:307
 TRolke.cxx:308
 TRolke.cxx:309
 TRolke.cxx:310
 TRolke.cxx:311
 TRolke.cxx:312
 TRolke.cxx:313
 TRolke.cxx:314
 TRolke.cxx:315
 TRolke.cxx:316
 TRolke.cxx:317
 TRolke.cxx:318
 TRolke.cxx:319
 TRolke.cxx:320
 TRolke.cxx:321
 TRolke.cxx:322
 TRolke.cxx:323
 TRolke.cxx:324
 TRolke.cxx:325
 TRolke.cxx:326
 TRolke.cxx:327
 TRolke.cxx:328
 TRolke.cxx:329
 TRolke.cxx:330
 TRolke.cxx:331
 TRolke.cxx:332
 TRolke.cxx:333
 TRolke.cxx:334
 TRolke.cxx:335
 TRolke.cxx:336
 TRolke.cxx:337
 TRolke.cxx:338
 TRolke.cxx:339
 TRolke.cxx:340
 TRolke.cxx:341
 TRolke.cxx:342
 TRolke.cxx:343
 TRolke.cxx:344
 TRolke.cxx:345
 TRolke.cxx:346
 TRolke.cxx:347
 TRolke.cxx:348
 TRolke.cxx:349
 TRolke.cxx:350
 TRolke.cxx:351
 TRolke.cxx:352
 TRolke.cxx:353
 TRolke.cxx:354
 TRolke.cxx:355
 TRolke.cxx:356
 TRolke.cxx:357
 TRolke.cxx:358
 TRolke.cxx:359
 TRolke.cxx:360
 TRolke.cxx:361
 TRolke.cxx:362
 TRolke.cxx:363
 TRolke.cxx:364
 TRolke.cxx:365
 TRolke.cxx:366
 TRolke.cxx:367
 TRolke.cxx:368
 TRolke.cxx:369
 TRolke.cxx:370
 TRolke.cxx:371
 TRolke.cxx:372
 TRolke.cxx:373
 TRolke.cxx:374
 TRolke.cxx:375
 TRolke.cxx:376
 TRolke.cxx:377
 TRolke.cxx:378
 TRolke.cxx:379
 TRolke.cxx:380
 TRolke.cxx:381
 TRolke.cxx:382
 TRolke.cxx:383
 TRolke.cxx:384
 TRolke.cxx:385
 TRolke.cxx:386
 TRolke.cxx:387
 TRolke.cxx:388
 TRolke.cxx:389
 TRolke.cxx:390
 TRolke.cxx:391
 TRolke.cxx:392
 TRolke.cxx:393
 TRolke.cxx:394
 TRolke.cxx:395
 TRolke.cxx:396
 TRolke.cxx:397
 TRolke.cxx:398
 TRolke.cxx:399
 TRolke.cxx:400
 TRolke.cxx:401
 TRolke.cxx:402
 TRolke.cxx:403
 TRolke.cxx:404
 TRolke.cxx:405
 TRolke.cxx:406
 TRolke.cxx:407
 TRolke.cxx:408
 TRolke.cxx:409
 TRolke.cxx:410
 TRolke.cxx:411
 TRolke.cxx:412
 TRolke.cxx:413
 TRolke.cxx:414
 TRolke.cxx:415
 TRolke.cxx:416
 TRolke.cxx:417
 TRolke.cxx:418
 TRolke.cxx:419
 TRolke.cxx:420
 TRolke.cxx:421
 TRolke.cxx:422
 TRolke.cxx:423
 TRolke.cxx:424
 TRolke.cxx:425
 TRolke.cxx:426
 TRolke.cxx:427
 TRolke.cxx:428
 TRolke.cxx:429
 TRolke.cxx:430
 TRolke.cxx:431
 TRolke.cxx:432
 TRolke.cxx:433
 TRolke.cxx:434
 TRolke.cxx:435
 TRolke.cxx:436
 TRolke.cxx:437
 TRolke.cxx:438
 TRolke.cxx:439
 TRolke.cxx:440
 TRolke.cxx:441
 TRolke.cxx:442
 TRolke.cxx:443
 TRolke.cxx:444
 TRolke.cxx:445
 TRolke.cxx:446
 TRolke.cxx:447
 TRolke.cxx:448
 TRolke.cxx:449
 TRolke.cxx:450
 TRolke.cxx:451
 TRolke.cxx:452
 TRolke.cxx:453
 TRolke.cxx:454
 TRolke.cxx:455
 TRolke.cxx:456
 TRolke.cxx:457
 TRolke.cxx:458
 TRolke.cxx:459
 TRolke.cxx:460
 TRolke.cxx:461
 TRolke.cxx:462
 TRolke.cxx:463
 TRolke.cxx:464
 TRolke.cxx:465
 TRolke.cxx:466
 TRolke.cxx:467
 TRolke.cxx:468
 TRolke.cxx:469
 TRolke.cxx:470
 TRolke.cxx:471
 TRolke.cxx:472
 TRolke.cxx:473
 TRolke.cxx:474
 TRolke.cxx:475
 TRolke.cxx:476
 TRolke.cxx:477
 TRolke.cxx:478
 TRolke.cxx:479
 TRolke.cxx:480
 TRolke.cxx:481
 TRolke.cxx:482
 TRolke.cxx:483
 TRolke.cxx:484
 TRolke.cxx:485
 TRolke.cxx:486
 TRolke.cxx:487
 TRolke.cxx:488
 TRolke.cxx:489
 TRolke.cxx:490
 TRolke.cxx:491
 TRolke.cxx:492
 TRolke.cxx:493
 TRolke.cxx:494
 TRolke.cxx:495
 TRolke.cxx:496
 TRolke.cxx:497
 TRolke.cxx:498
 TRolke.cxx:499
 TRolke.cxx:500
 TRolke.cxx:501
 TRolke.cxx:502
 TRolke.cxx:503
 TRolke.cxx:504
 TRolke.cxx:505
 TRolke.cxx:506
 TRolke.cxx:507
 TRolke.cxx:508
 TRolke.cxx:509
 TRolke.cxx:510
 TRolke.cxx:511
 TRolke.cxx:512
 TRolke.cxx:513
 TRolke.cxx:514
 TRolke.cxx:515
 TRolke.cxx:516
 TRolke.cxx:517
 TRolke.cxx:518
 TRolke.cxx:519
 TRolke.cxx:520
 TRolke.cxx:521
 TRolke.cxx:522
 TRolke.cxx:523
 TRolke.cxx:524
 TRolke.cxx:525
 TRolke.cxx:526
 TRolke.cxx:527
 TRolke.cxx:528
 TRolke.cxx:529
 TRolke.cxx:530
 TRolke.cxx:531
 TRolke.cxx:532
 TRolke.cxx:533
 TRolke.cxx:534
 TRolke.cxx:535
 TRolke.cxx:536
 TRolke.cxx:537
 TRolke.cxx:538
 TRolke.cxx:539
 TRolke.cxx:540
 TRolke.cxx:541
 TRolke.cxx:542
 TRolke.cxx:543
 TRolke.cxx:544
 TRolke.cxx:545
 TRolke.cxx:546
 TRolke.cxx:547
 TRolke.cxx:548
 TRolke.cxx:549
 TRolke.cxx:550
 TRolke.cxx:551
 TRolke.cxx:552
 TRolke.cxx:553
 TRolke.cxx:554
 TRolke.cxx:555
 TRolke.cxx:556
 TRolke.cxx:557
 TRolke.cxx:558
 TRolke.cxx:559
 TRolke.cxx:560
 TRolke.cxx:561
 TRolke.cxx:562
 TRolke.cxx:563
 TRolke.cxx:564
 TRolke.cxx:565
 TRolke.cxx:566
 TRolke.cxx:567
 TRolke.cxx:568
 TRolke.cxx:569
 TRolke.cxx:570
 TRolke.cxx:571
 TRolke.cxx:572
 TRolke.cxx:573
 TRolke.cxx:574
 TRolke.cxx:575
 TRolke.cxx:576
 TRolke.cxx:577
 TRolke.cxx:578
 TRolke.cxx:579
 TRolke.cxx:580
 TRolke.cxx:581
 TRolke.cxx:582
 TRolke.cxx:583
 TRolke.cxx:584
 TRolke.cxx:585
 TRolke.cxx:586
 TRolke.cxx:587
 TRolke.cxx:588
 TRolke.cxx:589
 TRolke.cxx:590
 TRolke.cxx:591
 TRolke.cxx:592
 TRolke.cxx:593
 TRolke.cxx:594
 TRolke.cxx:595
 TRolke.cxx:596
 TRolke.cxx:597
 TRolke.cxx:598
 TRolke.cxx:599
 TRolke.cxx:600
 TRolke.cxx:601
 TRolke.cxx:602
 TRolke.cxx:603
 TRolke.cxx:604
 TRolke.cxx:605
 TRolke.cxx:606
 TRolke.cxx:607
 TRolke.cxx:608
 TRolke.cxx:609
 TRolke.cxx:610
 TRolke.cxx:611
 TRolke.cxx:612
 TRolke.cxx:613
 TRolke.cxx:614
 TRolke.cxx:615
 TRolke.cxx:616
 TRolke.cxx:617
 TRolke.cxx:618
 TRolke.cxx:619
 TRolke.cxx:620
 TRolke.cxx:621
 TRolke.cxx:622
 TRolke.cxx:623
 TRolke.cxx:624
 TRolke.cxx:625
 TRolke.cxx:626
 TRolke.cxx:627
 TRolke.cxx:628
 TRolke.cxx:629
 TRolke.cxx:630
 TRolke.cxx:631
 TRolke.cxx:632
 TRolke.cxx:633
 TRolke.cxx:634
 TRolke.cxx:635
 TRolke.cxx:636
 TRolke.cxx:637
 TRolke.cxx:638
 TRolke.cxx:639
 TRolke.cxx:640
 TRolke.cxx:641
 TRolke.cxx:642
 TRolke.cxx:643
 TRolke.cxx:644
 TRolke.cxx:645
 TRolke.cxx:646
 TRolke.cxx:647
 TRolke.cxx:648
 TRolke.cxx:649
 TRolke.cxx:650
 TRolke.cxx:651
 TRolke.cxx:652
 TRolke.cxx:653
 TRolke.cxx:654
 TRolke.cxx:655
 TRolke.cxx:656
 TRolke.cxx:657
 TRolke.cxx:658
 TRolke.cxx:659
 TRolke.cxx:660
 TRolke.cxx:661
 TRolke.cxx:662
 TRolke.cxx:663
 TRolke.cxx:664
 TRolke.cxx:665
 TRolke.cxx:666
 TRolke.cxx:667
 TRolke.cxx:668
 TRolke.cxx:669
 TRolke.cxx:670
 TRolke.cxx:671
 TRolke.cxx:672
 TRolke.cxx:673
 TRolke.cxx:674
 TRolke.cxx:675
 TRolke.cxx:676
 TRolke.cxx:677
 TRolke.cxx:678
 TRolke.cxx:679
 TRolke.cxx:680
 TRolke.cxx:681
 TRolke.cxx:682
 TRolke.cxx:683
 TRolke.cxx:684
 TRolke.cxx:685
 TRolke.cxx:686
 TRolke.cxx:687
 TRolke.cxx:688
 TRolke.cxx:689
 TRolke.cxx:690
 TRolke.cxx:691
 TRolke.cxx:692
 TRolke.cxx:693
 TRolke.cxx:694
 TRolke.cxx:695
 TRolke.cxx:696
 TRolke.cxx:697
 TRolke.cxx:698
 TRolke.cxx:699
 TRolke.cxx:700
 TRolke.cxx:701
 TRolke.cxx:702
 TRolke.cxx:703
 TRolke.cxx:704
 TRolke.cxx:705
 TRolke.cxx:706
 TRolke.cxx:707
 TRolke.cxx:708
 TRolke.cxx:709
 TRolke.cxx:710
 TRolke.cxx:711
 TRolke.cxx:712
 TRolke.cxx:713
 TRolke.cxx:714
 TRolke.cxx:715
 TRolke.cxx:716
 TRolke.cxx:717
 TRolke.cxx:718
 TRolke.cxx:719
 TRolke.cxx:720
 TRolke.cxx:721
 TRolke.cxx:722
 TRolke.cxx:723
 TRolke.cxx:724
 TRolke.cxx:725
 TRolke.cxx:726
 TRolke.cxx:727
 TRolke.cxx:728
 TRolke.cxx:729
 TRolke.cxx:730
 TRolke.cxx:731
 TRolke.cxx:732
 TRolke.cxx:733
 TRolke.cxx:734
 TRolke.cxx:735
 TRolke.cxx:736
 TRolke.cxx:737
 TRolke.cxx:738
 TRolke.cxx:739
 TRolke.cxx:740
 TRolke.cxx:741
 TRolke.cxx:742
 TRolke.cxx:743
 TRolke.cxx:744
 TRolke.cxx:745
 TRolke.cxx:746
 TRolke.cxx:747
 TRolke.cxx:748
 TRolke.cxx:749
 TRolke.cxx:750
 TRolke.cxx:751
 TRolke.cxx:752
 TRolke.cxx:753
 TRolke.cxx:754
 TRolke.cxx:755
 TRolke.cxx:756
 TRolke.cxx:757
 TRolke.cxx:758
 TRolke.cxx:759
 TRolke.cxx:760
 TRolke.cxx:761
 TRolke.cxx:762
 TRolke.cxx:763
 TRolke.cxx:764
 TRolke.cxx:765
 TRolke.cxx:766
 TRolke.cxx:767
 TRolke.cxx:768
 TRolke.cxx:769
 TRolke.cxx:770
 TRolke.cxx:771
 TRolke.cxx:772
 TRolke.cxx:773
 TRolke.cxx:774
 TRolke.cxx:775
 TRolke.cxx:776
 TRolke.cxx:777
 TRolke.cxx:778
 TRolke.cxx:779
 TRolke.cxx:780
 TRolke.cxx:781
 TRolke.cxx:782
 TRolke.cxx:783
 TRolke.cxx:784
 TRolke.cxx:785
 TRolke.cxx:786
 TRolke.cxx:787
 TRolke.cxx:788
 TRolke.cxx:789
 TRolke.cxx:790
 TRolke.cxx:791
 TRolke.cxx:792
 TRolke.cxx:793
 TRolke.cxx:794
 TRolke.cxx:795
 TRolke.cxx:796
 TRolke.cxx:797
 TRolke.cxx:798
 TRolke.cxx:799
 TRolke.cxx:800
 TRolke.cxx:801
 TRolke.cxx:802
 TRolke.cxx:803
 TRolke.cxx:804
 TRolke.cxx:805
 TRolke.cxx:806
 TRolke.cxx:807
 TRolke.cxx:808
 TRolke.cxx:809
 TRolke.cxx:810
 TRolke.cxx:811
 TRolke.cxx:812
 TRolke.cxx:813
 TRolke.cxx:814
 TRolke.cxx:815
 TRolke.cxx:816
 TRolke.cxx:817
 TRolke.cxx:818
 TRolke.cxx:819
 TRolke.cxx:820
 TRolke.cxx:821
 TRolke.cxx:822
 TRolke.cxx:823
 TRolke.cxx:824
 TRolke.cxx:825
 TRolke.cxx:826
 TRolke.cxx:827
 TRolke.cxx:828
 TRolke.cxx:829
 TRolke.cxx:830
 TRolke.cxx:831
 TRolke.cxx:832
 TRolke.cxx:833
 TRolke.cxx:834
 TRolke.cxx:835
 TRolke.cxx:836
 TRolke.cxx:837
 TRolke.cxx:838
 TRolke.cxx:839
 TRolke.cxx:840
 TRolke.cxx:841
 TRolke.cxx:842
 TRolke.cxx:843
 TRolke.cxx:844
 TRolke.cxx:845
 TRolke.cxx:846
 TRolke.cxx:847
 TRolke.cxx:848
 TRolke.cxx:849
 TRolke.cxx:850
 TRolke.cxx:851
 TRolke.cxx:852
 TRolke.cxx:853
 TRolke.cxx:854
 TRolke.cxx:855
 TRolke.cxx:856
 TRolke.cxx:857
 TRolke.cxx:858
 TRolke.cxx:859
 TRolke.cxx:860
 TRolke.cxx:861
 TRolke.cxx:862
 TRolke.cxx:863
 TRolke.cxx:864
 TRolke.cxx:865
 TRolke.cxx:866
 TRolke.cxx:867
 TRolke.cxx:868
 TRolke.cxx:869
 TRolke.cxx:870
 TRolke.cxx:871
 TRolke.cxx:872
 TRolke.cxx:873
 TRolke.cxx:874
 TRolke.cxx:875
 TRolke.cxx:876
 TRolke.cxx:877
 TRolke.cxx:878
 TRolke.cxx:879
 TRolke.cxx:880
 TRolke.cxx:881
 TRolke.cxx:882
 TRolke.cxx:883
 TRolke.cxx:884
 TRolke.cxx:885
 TRolke.cxx:886
 TRolke.cxx:887
 TRolke.cxx:888
 TRolke.cxx:889
 TRolke.cxx:890
 TRolke.cxx:891
 TRolke.cxx:892
 TRolke.cxx:893
 TRolke.cxx:894
 TRolke.cxx:895
 TRolke.cxx:896
 TRolke.cxx:897
 TRolke.cxx:898
 TRolke.cxx:899
 TRolke.cxx:900
 TRolke.cxx:901
 TRolke.cxx:902
 TRolke.cxx:903
 TRolke.cxx:904
 TRolke.cxx:905
 TRolke.cxx:906
 TRolke.cxx:907
 TRolke.cxx:908
 TRolke.cxx:909
 TRolke.cxx:910
 TRolke.cxx:911
 TRolke.cxx:912
 TRolke.cxx:913
 TRolke.cxx:914
 TRolke.cxx:915
 TRolke.cxx:916
 TRolke.cxx:917
 TRolke.cxx:918
 TRolke.cxx:919
 TRolke.cxx:920
 TRolke.cxx:921
 TRolke.cxx:922
 TRolke.cxx:923
 TRolke.cxx:924
 TRolke.cxx:925
 TRolke.cxx:926
 TRolke.cxx:927
 TRolke.cxx:928
 TRolke.cxx:929
 TRolke.cxx:930
 TRolke.cxx:931
 TRolke.cxx:932
 TRolke.cxx:933
 TRolke.cxx:934
 TRolke.cxx:935
 TRolke.cxx:936
 TRolke.cxx:937
 TRolke.cxx:938
 TRolke.cxx:939
 TRolke.cxx:940
 TRolke.cxx:941
 TRolke.cxx:942
 TRolke.cxx:943
 TRolke.cxx:944
 TRolke.cxx:945
 TRolke.cxx:946
 TRolke.cxx:947
 TRolke.cxx:948
 TRolke.cxx:949
 TRolke.cxx:950
 TRolke.cxx:951
 TRolke.cxx:952
 TRolke.cxx:953
 TRolke.cxx:954
 TRolke.cxx:955
 TRolke.cxx:956
 TRolke.cxx:957
 TRolke.cxx:958
 TRolke.cxx:959
 TRolke.cxx:960
 TRolke.cxx:961
 TRolke.cxx:962
 TRolke.cxx:963
 TRolke.cxx:964
 TRolke.cxx:965
 TRolke.cxx:966
 TRolke.cxx:967
 TRolke.cxx:968
 TRolke.cxx:969
 TRolke.cxx:970
 TRolke.cxx:971
 TRolke.cxx:972
 TRolke.cxx:973
 TRolke.cxx:974
 TRolke.cxx:975
 TRolke.cxx:976
 TRolke.cxx:977
 TRolke.cxx:978
 TRolke.cxx:979
 TRolke.cxx:980
 TRolke.cxx:981
 TRolke.cxx:982
 TRolke.cxx:983
 TRolke.cxx:984
 TRolke.cxx:985
 TRolke.cxx:986
 TRolke.cxx:987
 TRolke.cxx:988
 TRolke.cxx:989
 TRolke.cxx:990
 TRolke.cxx:991
 TRolke.cxx:992
 TRolke.cxx:993
 TRolke.cxx:994
 TRolke.cxx:995
 TRolke.cxx:996
 TRolke.cxx:997
 TRolke.cxx:998
 TRolke.cxx:999
 TRolke.cxx:1000
 TRolke.cxx:1001
 TRolke.cxx:1002
 TRolke.cxx:1003
 TRolke.cxx:1004
 TRolke.cxx:1005
 TRolke.cxx:1006
 TRolke.cxx:1007
 TRolke.cxx:1008
 TRolke.cxx:1009
 TRolke.cxx:1010
 TRolke.cxx:1011
 TRolke.cxx:1012
 TRolke.cxx:1013
 TRolke.cxx:1014
 TRolke.cxx:1015
 TRolke.cxx:1016
 TRolke.cxx:1017
 TRolke.cxx:1018
 TRolke.cxx:1019
 TRolke.cxx:1020
 TRolke.cxx:1021
 TRolke.cxx:1022
 TRolke.cxx:1023
 TRolke.cxx:1024
 TRolke.cxx:1025
 TRolke.cxx:1026
 TRolke.cxx:1027
 TRolke.cxx:1028
 TRolke.cxx:1029
 TRolke.cxx:1030
 TRolke.cxx:1031
 TRolke.cxx:1032
 TRolke.cxx:1033
 TRolke.cxx:1034
 TRolke.cxx:1035
 TRolke.cxx:1036
 TRolke.cxx:1037
 TRolke.cxx:1038
 TRolke.cxx:1039
 TRolke.cxx:1040
 TRolke.cxx:1041
 TRolke.cxx:1042
 TRolke.cxx:1043
 TRolke.cxx:1044
 TRolke.cxx:1045
 TRolke.cxx:1046
 TRolke.cxx:1047
 TRolke.cxx:1048
 TRolke.cxx:1049
 TRolke.cxx:1050
 TRolke.cxx:1051
 TRolke.cxx:1052
 TRolke.cxx:1053
 TRolke.cxx:1054
 TRolke.cxx:1055
 TRolke.cxx:1056
 TRolke.cxx:1057
 TRolke.cxx:1058
 TRolke.cxx:1059
 TRolke.cxx:1060
 TRolke.cxx:1061
 TRolke.cxx:1062
 TRolke.cxx:1063
 TRolke.cxx:1064
 TRolke.cxx:1065
 TRolke.cxx:1066
 TRolke.cxx:1067
 TRolke.cxx:1068
 TRolke.cxx:1069
 TRolke.cxx:1070
 TRolke.cxx:1071
 TRolke.cxx:1072
 TRolke.cxx:1073
 TRolke.cxx:1074
 TRolke.cxx:1075
 TRolke.cxx:1076
 TRolke.cxx:1077
 TRolke.cxx:1078
 TRolke.cxx:1079
 TRolke.cxx:1080
 TRolke.cxx:1081
 TRolke.cxx:1082
 TRolke.cxx:1083
 TRolke.cxx:1084
 TRolke.cxx:1085
 TRolke.cxx:1086
 TRolke.cxx:1087
 TRolke.cxx:1088
 TRolke.cxx:1089
 TRolke.cxx:1090
 TRolke.cxx:1091
 TRolke.cxx:1092
 TRolke.cxx:1093
 TRolke.cxx:1094
 TRolke.cxx:1095
 TRolke.cxx:1096
 TRolke.cxx:1097
 TRolke.cxx:1098
 TRolke.cxx:1099
 TRolke.cxx:1100
 TRolke.cxx:1101
 TRolke.cxx:1102
 TRolke.cxx:1103
 TRolke.cxx:1104
 TRolke.cxx:1105
 TRolke.cxx:1106
 TRolke.cxx:1107
 TRolke.cxx:1108
 TRolke.cxx:1109
 TRolke.cxx:1110
 TRolke.cxx:1111
 TRolke.cxx:1112
 TRolke.cxx:1113
 TRolke.cxx:1114
 TRolke.cxx:1115
 TRolke.cxx:1116
 TRolke.cxx:1117
 TRolke.cxx:1118
 TRolke.cxx:1119
 TRolke.cxx:1120
 TRolke.cxx:1121
 TRolke.cxx:1122
 TRolke.cxx:1123
 TRolke.cxx:1124
 TRolke.cxx:1125
 TRolke.cxx:1126
 TRolke.cxx:1127
 TRolke.cxx:1128
 TRolke.cxx:1129
 TRolke.cxx:1130
 TRolke.cxx:1131
 TRolke.cxx:1132
 TRolke.cxx:1133
 TRolke.cxx:1134
 TRolke.cxx:1135
 TRolke.cxx:1136
 TRolke.cxx:1137
 TRolke.cxx:1138
 TRolke.cxx:1139
 TRolke.cxx:1140
 TRolke.cxx:1141
 TRolke.cxx:1142
 TRolke.cxx:1143
 TRolke.cxx:1144
 TRolke.cxx:1145
 TRolke.cxx:1146
 TRolke.cxx:1147
 TRolke.cxx:1148
 TRolke.cxx:1149
 TRolke.cxx:1150
 TRolke.cxx:1151
 TRolke.cxx:1152
 TRolke.cxx:1153
 TRolke.cxx:1154
 TRolke.cxx:1155
 TRolke.cxx:1156
 TRolke.cxx:1157
 TRolke.cxx:1158
 TRolke.cxx:1159
 TRolke.cxx:1160
 TRolke.cxx:1161
 TRolke.cxx:1162
 TRolke.cxx:1163
 TRolke.cxx:1164
 TRolke.cxx:1165
 TRolke.cxx:1166
 TRolke.cxx:1167
 TRolke.cxx:1168
 TRolke.cxx:1169
 TRolke.cxx:1170
 TRolke.cxx:1171
 TRolke.cxx:1172
 TRolke.cxx:1173
 TRolke.cxx:1174
 TRolke.cxx:1175
 TRolke.cxx:1176
 TRolke.cxx:1177
 TRolke.cxx:1178
 TRolke.cxx:1179
 TRolke.cxx:1180
 TRolke.cxx:1181
 TRolke.cxx:1182
 TRolke.cxx:1183
 TRolke.cxx:1184
 TRolke.cxx:1185
 TRolke.cxx:1186
 TRolke.cxx:1187
 TRolke.cxx:1188
 TRolke.cxx:1189
 TRolke.cxx:1190
 TRolke.cxx:1191
 TRolke.cxx:1192
 TRolke.cxx:1193
 TRolke.cxx:1194
 TRolke.cxx:1195
 TRolke.cxx:1196
 TRolke.cxx:1197
 TRolke.cxx:1198
 TRolke.cxx:1199
 TRolke.cxx:1200
 TRolke.cxx:1201
 TRolke.cxx:1202
 TRolke.cxx:1203
 TRolke.cxx:1204
 TRolke.cxx:1205
 TRolke.cxx:1206
 TRolke.cxx:1207
 TRolke.cxx:1208
 TRolke.cxx:1209
 TRolke.cxx:1210
 TRolke.cxx:1211
 TRolke.cxx:1212
 TRolke.cxx:1213
 TRolke.cxx:1214
 TRolke.cxx:1215
 TRolke.cxx:1216
 TRolke.cxx:1217
 TRolke.cxx:1218
 TRolke.cxx:1219
 TRolke.cxx:1220
 TRolke.cxx:1221
 TRolke.cxx:1222
 TRolke.cxx:1223
 TRolke.cxx:1224
 TRolke.cxx:1225
 TRolke.cxx:1226
 TRolke.cxx:1227
 TRolke.cxx:1228
 TRolke.cxx:1229
 TRolke.cxx:1230
 TRolke.cxx:1231
 TRolke.cxx:1232
 TRolke.cxx:1233
 TRolke.cxx:1234
 TRolke.cxx:1235
 TRolke.cxx:1236
 TRolke.cxx:1237
 TRolke.cxx:1238
 TRolke.cxx:1239
 TRolke.cxx:1240
 TRolke.cxx:1241
 TRolke.cxx:1242
 TRolke.cxx:1243
 TRolke.cxx:1244
 TRolke.cxx:1245
 TRolke.cxx:1246
 TRolke.cxx:1247
 TRolke.cxx:1248
 TRolke.cxx:1249
 TRolke.cxx:1250
 TRolke.cxx:1251
 TRolke.cxx:1252
 TRolke.cxx:1253
 TRolke.cxx:1254
 TRolke.cxx:1255
 TRolke.cxx:1256
 TRolke.cxx:1257
 TRolke.cxx:1258
 TRolke.cxx:1259
 TRolke.cxx:1260
 TRolke.cxx:1261
 TRolke.cxx:1262
 TRolke.cxx:1263
 TRolke.cxx:1264
 TRolke.cxx:1265
 TRolke.cxx:1266
 TRolke.cxx:1267
 TRolke.cxx:1268
 TRolke.cxx:1269
 TRolke.cxx:1270
 TRolke.cxx:1271
 TRolke.cxx:1272
 TRolke.cxx:1273
 TRolke.cxx:1274
 TRolke.cxx:1275
 TRolke.cxx:1276
 TRolke.cxx:1277
 TRolke.cxx:1278
 TRolke.cxx:1279
 TRolke.cxx:1280
 TRolke.cxx:1281
 TRolke.cxx:1282
 TRolke.cxx:1283
 TRolke.cxx:1284
 TRolke.cxx:1285
 TRolke.cxx:1286
 TRolke.cxx:1287
 TRolke.cxx:1288
 TRolke.cxx:1289
 TRolke.cxx:1290
 TRolke.cxx:1291
 TRolke.cxx:1292
 TRolke.cxx:1293
 TRolke.cxx:1294
 TRolke.cxx:1295
 TRolke.cxx:1296
 TRolke.cxx:1297
 TRolke.cxx:1298
 TRolke.cxx:1299
 TRolke.cxx:1300
 TRolke.cxx:1301
 TRolke.cxx:1302
 TRolke.cxx:1303
 TRolke.cxx:1304
 TRolke.cxx:1305
 TRolke.cxx:1306
 TRolke.cxx:1307
 TRolke.cxx:1308
 TRolke.cxx:1309
 TRolke.cxx:1310
 TRolke.cxx:1311
 TRolke.cxx:1312
 TRolke.cxx:1313
 TRolke.cxx:1314
 TRolke.cxx:1315
 TRolke.cxx:1316
 TRolke.cxx:1317
 TRolke.cxx:1318
 TRolke.cxx:1319
 TRolke.cxx:1320
 TRolke.cxx:1321
 TRolke.cxx:1322
 TRolke.cxx:1323
 TRolke.cxx:1324
 TRolke.cxx:1325
 TRolke.cxx:1326
 TRolke.cxx:1327
 TRolke.cxx:1328
 TRolke.cxx:1329
 TRolke.cxx:1330
 TRolke.cxx:1331
 TRolke.cxx:1332
 TRolke.cxx:1333
 TRolke.cxx:1334
 TRolke.cxx:1335
 TRolke.cxx:1336
 TRolke.cxx:1337
 TRolke.cxx:1338
 TRolke.cxx:1339
 TRolke.cxx:1340
 TRolke.cxx:1341
 TRolke.cxx:1342
 TRolke.cxx:1343
 TRolke.cxx:1344
 TRolke.cxx:1345
 TRolke.cxx:1346
 TRolke.cxx:1347
 TRolke.cxx:1348
 TRolke.cxx:1349
 TRolke.cxx:1350
 TRolke.cxx:1351
 TRolke.cxx:1352
 TRolke.cxx:1353
 TRolke.cxx:1354
 TRolke.cxx:1355
 TRolke.cxx:1356
 TRolke.cxx:1357
 TRolke.cxx:1358
 TRolke.cxx:1359
 TRolke.cxx:1360
 TRolke.cxx:1361
 TRolke.cxx:1362
 TRolke.cxx:1363
 TRolke.cxx:1364
 TRolke.cxx:1365
 TRolke.cxx:1366
 TRolke.cxx:1367
 TRolke.cxx:1368
 TRolke.cxx:1369
 TRolke.cxx:1370
 TRolke.cxx:1371
 TRolke.cxx:1372
 TRolke.cxx:1373
 TRolke.cxx:1374
 TRolke.cxx:1375
 TRolke.cxx:1376
 TRolke.cxx:1377
 TRolke.cxx:1378
 TRolke.cxx:1379
 TRolke.cxx:1380
 TRolke.cxx:1381
 TRolke.cxx:1382
 TRolke.cxx:1383
 TRolke.cxx:1384
 TRolke.cxx:1385
 TRolke.cxx:1386
 TRolke.cxx:1387
 TRolke.cxx:1388
 TRolke.cxx:1389
 TRolke.cxx:1390
 TRolke.cxx:1391
 TRolke.cxx:1392
 TRolke.cxx:1393
 TRolke.cxx:1394
 TRolke.cxx:1395
 TRolke.cxx:1396
 TRolke.cxx:1397
 TRolke.cxx:1398
 TRolke.cxx:1399
 TRolke.cxx:1400
 TRolke.cxx:1401
 TRolke.cxx:1402
 TRolke.cxx:1403
 TRolke.cxx:1404
 TRolke.cxx:1405
 TRolke.cxx:1406
 TRolke.cxx:1407
 TRolke.cxx:1408
 TRolke.cxx:1409
 TRolke.cxx:1410
 TRolke.cxx:1411
 TRolke.cxx:1412
 TRolke.cxx:1413
 TRolke.cxx:1414
 TRolke.cxx:1415
 TRolke.cxx:1416
 TRolke.cxx:1417
 TRolke.cxx:1418
 TRolke.cxx:1419
 TRolke.cxx:1420
 TRolke.cxx:1421
 TRolke.cxx:1422
 TRolke.cxx:1423
 TRolke.cxx:1424
 TRolke.cxx:1425
 TRolke.cxx:1426
 TRolke.cxx:1427
 TRolke.cxx:1428
 TRolke.cxx:1429
 TRolke.cxx:1430
 TRolke.cxx:1431
 TRolke.cxx:1432
 TRolke.cxx:1433
 TRolke.cxx:1434
 TRolke.cxx:1435
 TRolke.cxx:1436
 TRolke.cxx:1437
 TRolke.cxx:1438
 TRolke.cxx:1439
 TRolke.cxx:1440
 TRolke.cxx:1441
 TRolke.cxx:1442
 TRolke.cxx:1443
 TRolke.cxx:1444