// @(#)root/physics:$Id$
// Author: Adrian Bevan  2001

/*************************************************************************
 * Copyright (C) 1995-2004, Rene Brun and Fons Rademakers.               *
 * Copyright (C) 2001, Liverpool University.                             *
 * All rights reserved.                                                  *
 *                                                                       *
 * For the licensing terms see $ROOTSYS/LICENSE.                         *
 * For the list of contributors see $ROOTSYS/README/CREDITS.             *
 *************************************************************************/

////////////////////////////////////////////////////////////////////////////
// TFeldmanCousins
//
// class to calculate the CL upper limit using
// the Feldman-Cousins method as described in PRD V57 #7, p3873-3889
//
// The default confidence interval calvculated using this method is 90%
// This is set either by having a default the constructor, or using the
// appropriate fraction when instantiating an object of this class (e.g. 0.9)
//
// The simple extension to a gaussian resolution function bounded at zero
// has not been addressed as yet -> `time is of the essence' as they write
// on the wall of the maze in that classic game ...
//
//    VARIABLES THAT CAN BE ALTERED
//    -----------------------------
// => depending on your desired precision: The intial values of fMuMin,
// fMuMax, fMuStep and fNMax are those used in the PRD:
//   fMuMin = 0.0
//   fMuMax = 50.0
//   fMuStep= 0.005
// but there is total flexibility in changing this should you desire.
//
//
// see example of use in $ROOTSYS/tutorials/math/FeldmanCousins.C
//
// see note about: "Should I use TRolke, TFeldmanCousins, TLimit?"
//  in the TRolke class description.
//
// Author: Adrian Bevan, Liverpool University
//
// Copyright Liverpool University 2001       bevan@slac.stanford.edu
///////////////////////////////////////////////////////////////////////////

#include "Riostream.h"
#include "TMath.h"
#include "TFeldmanCousins.h"

ClassImp(TFeldmanCousins)

//______________________________________________________________________________
TFeldmanCousins::TFeldmanCousins(Double_t newFC, TString options)
{
   //constructor
   fCL          = newFC;
   fUpperLimit  = 0.0;
   fLowerLimit  = 0.0;
   fNobserved   = 0.0;
   fNbackground = 0.0;
   options.ToLower();
   if (options.Contains("q")) fQUICK = 1;
   else                       fQUICK = 0;

   fNMax   = 50;
   fMuStep = 0.005;
   SetMuMin();
   SetMuMax();
   SetMuStep();
}


//______________________________________________________________________________
TFeldmanCousins::~TFeldmanCousins()
{
}


//______________________________________________________________________________
Double_t TFeldmanCousins::CalculateLowerLimit(Double_t Nobserved, Double_t Nbackground)
{
////////////////////////////////////////////////////////////////////////////////////////////
// given Nobserved and Nbackground, try different values of mu that give lower limits that//
// are consistent with Nobserved.  The closed interval (plus any stragglers) corresponds  //
// to the F&C interval                                                                    //
////////////////////////////////////////////////////////////////////////////////////////////

   CalculateUpperLimit(Nobserved, Nbackground);
   return fLowerLimit;
}


//______________________________________________________________________________
Double_t TFeldmanCousins::CalculateUpperLimit(Double_t Nobserved, Double_t Nbackground)
{
////////////////////////////////////////////////////////////////////////////////////////////
// given Nobserved and Nbackground, try different values of mu that give upper limits that//
// are consistent with Nobserved.  The closed interval (plus any stragglers) corresponds  //
// to the F&C interval                                                                    //
////////////////////////////////////////////////////////////////////////////////////////////

   fNobserved   = Nobserved;
   fNbackground = Nbackground;

   Double_t mu = 0.0;

   // for each mu construct the ranked table of probabilities and test the
   // observed number of events with the upper limit
   Double_t min = -999.0;
   Double_t max = 0;
   Int_t iLower = 0;

   Int_t i;
   for(i = 0; i <= fNMuStep; i++) {
      mu = fMuMin + (Double_t)i*fMuStep;
      Int_t goodChoice = FindLimitsFromTable( mu );
      if( goodChoice ) {
         min = mu;
         iLower = i;
         break;
      }
   }

   //==================================================================
   // For quicker evaluation, assume that you get the same results when
   // you expect the uppper limit to be > Nobserved-Nbackground.
   // This is certainly true for all of the published tables in the PRD
   // and is a reasonable assumption in any case.
   //==================================================================

   Double_t quickJump = 0.0;
   if (fQUICK)          quickJump = Nobserved-Nbackground-fMuMin;
   if (quickJump < 0.0) quickJump = 0.0;

   for(i = iLower+1; i <= fNMuStep; i++) {
      mu = fMuMin + (Double_t)i*fMuStep + quickJump;
      Int_t goodChoice = FindLimitsFromTable( mu );
      if( !goodChoice ) {
         max = mu;
         break;
      }
   }

   fUpperLimit = max;
   fLowerLimit = min;

   return max;
}


//______________________________________________________________________________
Int_t TFeldmanCousins::FindLimitsFromTable( Double_t mu )
{
///////////////////////////////////////////////////////////////////
// calculate the probability table for a given mu for n = 0, NMAX//
// and return 1 if the number of observed events is consistent   //
// with the CL bad                                               //
///////////////////////////////////////////////////////////////////

   Double_t *p          = new Double_t[fNMax];   //the array of probabilities in the interval MUMIN-MUMAX
   Double_t *r          = new Double_t[fNMax];   //the ratio of likliehoods = P(Mu|Nobserved)/P(MuBest|Nobserved)
   Int_t    *rank       = new Int_t[fNMax];      //the ranked array corresponding to R (largest first)
   Double_t *muBest     = new Double_t[fNMax];
   Double_t *probMuBest = new Double_t[fNMax];

   //calculate P(i | mu) and P(i | mu)/P(i | mubest)
   Int_t i;
   for(i = 0; i < fNMax; i++) {
      muBest[i] = (Double_t)(i - fNbackground);
      if(muBest[i]<0.0) muBest[i] = 0.0;
      probMuBest[i] = Prob(i, muBest[i],  fNbackground);
      p[i]          = Prob(i, mu,  fNbackground);
      if(probMuBest[i] == 0.0) r[i] = 0.0;
      else                     r[i] = p[i]/probMuBest[i];
   }

   //rank the likelihood ratio
   TMath::Sort(fNMax, r, rank);

   //search through the probability table and get the i for the CL
   Double_t sum = 0.0;
   Int_t iMax = rank[0];
   Int_t iMin = rank[0];
   for(i = 0; i < fNMax; i++) {
      sum += p[rank[i]];
      if(iMax < rank[i]) iMax = rank[i];
      if(iMin > rank[i]) iMin = rank[i];
      if(sum >= fCL) break;
   }

   delete [] p;
   delete [] r;
   delete [] rank;
   delete [] muBest;
   delete [] probMuBest;

   if((fNobserved <= iMax) && (fNobserved >= iMin)) return 1;
   else return 0;
}


//______________________________________________________________________________
Double_t TFeldmanCousins::Prob(Int_t N, Double_t mu, Double_t B)
{
////////////////////////////////////////////////
// calculate the poissonian probability for   //
// a mean of mu+B events with a variance of N //
////////////////////////////////////////////////

   return TMath::Poisson( N, mu+B);
}

//______________________________________________________________________________
void TFeldmanCousins::SetMuMax(Double_t newMax)
{
   //set maximum value of signal to use in calculating the tables
   fMuMax   = newMax;
   fNMax    = (Int_t)newMax;
   SetMuStep(fMuStep);
}

//______________________________________________________________________________
void TFeldmanCousins::SetMuStep(Double_t newMuStep)
{
   //set the step in signal to use when generating tables
   if(newMuStep == 0.0) {
      std::cout << "TFeldmanCousins::SetMuStep ERROR New step size is zero - unable to change value"<< std::endl;
      return;
   } else {
      fMuStep = newMuStep;
      fNMuStep = (Int_t)((fMuMax - fMuMin)/fMuStep);
   }
}
 TFeldmanCousins.cxx:1
 TFeldmanCousins.cxx:2
 TFeldmanCousins.cxx:3
 TFeldmanCousins.cxx:4
 TFeldmanCousins.cxx:5
 TFeldmanCousins.cxx:6
 TFeldmanCousins.cxx:7
 TFeldmanCousins.cxx:8
 TFeldmanCousins.cxx:9
 TFeldmanCousins.cxx:10
 TFeldmanCousins.cxx:11
 TFeldmanCousins.cxx:12
 TFeldmanCousins.cxx:13
 TFeldmanCousins.cxx:14
 TFeldmanCousins.cxx:15
 TFeldmanCousins.cxx:16
 TFeldmanCousins.cxx:17
 TFeldmanCousins.cxx:18
 TFeldmanCousins.cxx:19
 TFeldmanCousins.cxx:20
 TFeldmanCousins.cxx:21
 TFeldmanCousins.cxx:22
 TFeldmanCousins.cxx:23
 TFeldmanCousins.cxx:24
 TFeldmanCousins.cxx:25
 TFeldmanCousins.cxx:26
 TFeldmanCousins.cxx:27
 TFeldmanCousins.cxx:28
 TFeldmanCousins.cxx:29
 TFeldmanCousins.cxx:30
 TFeldmanCousins.cxx:31
 TFeldmanCousins.cxx:32
 TFeldmanCousins.cxx:33
 TFeldmanCousins.cxx:34
 TFeldmanCousins.cxx:35
 TFeldmanCousins.cxx:36
 TFeldmanCousins.cxx:37
 TFeldmanCousins.cxx:38
 TFeldmanCousins.cxx:39
 TFeldmanCousins.cxx:40
 TFeldmanCousins.cxx:41
 TFeldmanCousins.cxx:42
 TFeldmanCousins.cxx:43
 TFeldmanCousins.cxx:44
 TFeldmanCousins.cxx:45
 TFeldmanCousins.cxx:46
 TFeldmanCousins.cxx:47
 TFeldmanCousins.cxx:48
 TFeldmanCousins.cxx:49
 TFeldmanCousins.cxx:50
 TFeldmanCousins.cxx:51
 TFeldmanCousins.cxx:52
 TFeldmanCousins.cxx:53
 TFeldmanCousins.cxx:54
 TFeldmanCousins.cxx:55
 TFeldmanCousins.cxx:56
 TFeldmanCousins.cxx:57
 TFeldmanCousins.cxx:58
 TFeldmanCousins.cxx:59
 TFeldmanCousins.cxx:60
 TFeldmanCousins.cxx:61
 TFeldmanCousins.cxx:62
 TFeldmanCousins.cxx:63
 TFeldmanCousins.cxx:64
 TFeldmanCousins.cxx:65
 TFeldmanCousins.cxx:66
 TFeldmanCousins.cxx:67
 TFeldmanCousins.cxx:68
 TFeldmanCousins.cxx:69
 TFeldmanCousins.cxx:70
 TFeldmanCousins.cxx:71
 TFeldmanCousins.cxx:72
 TFeldmanCousins.cxx:73
 TFeldmanCousins.cxx:74
 TFeldmanCousins.cxx:75
 TFeldmanCousins.cxx:76
 TFeldmanCousins.cxx:77
 TFeldmanCousins.cxx:78
 TFeldmanCousins.cxx:79
 TFeldmanCousins.cxx:80
 TFeldmanCousins.cxx:81
 TFeldmanCousins.cxx:82
 TFeldmanCousins.cxx:83
 TFeldmanCousins.cxx:84
 TFeldmanCousins.cxx:85
 TFeldmanCousins.cxx:86
 TFeldmanCousins.cxx:87
 TFeldmanCousins.cxx:88
 TFeldmanCousins.cxx:89
 TFeldmanCousins.cxx:90
 TFeldmanCousins.cxx:91
 TFeldmanCousins.cxx:92
 TFeldmanCousins.cxx:93
 TFeldmanCousins.cxx:94
 TFeldmanCousins.cxx:95
 TFeldmanCousins.cxx:96
 TFeldmanCousins.cxx:97
 TFeldmanCousins.cxx:98
 TFeldmanCousins.cxx:99
 TFeldmanCousins.cxx:100
 TFeldmanCousins.cxx:101
 TFeldmanCousins.cxx:102
 TFeldmanCousins.cxx:103
 TFeldmanCousins.cxx:104
 TFeldmanCousins.cxx:105
 TFeldmanCousins.cxx:106
 TFeldmanCousins.cxx:107
 TFeldmanCousins.cxx:108
 TFeldmanCousins.cxx:109
 TFeldmanCousins.cxx:110
 TFeldmanCousins.cxx:111
 TFeldmanCousins.cxx:112
 TFeldmanCousins.cxx:113
 TFeldmanCousins.cxx:114
 TFeldmanCousins.cxx:115
 TFeldmanCousins.cxx:116
 TFeldmanCousins.cxx:117
 TFeldmanCousins.cxx:118
 TFeldmanCousins.cxx:119
 TFeldmanCousins.cxx:120
 TFeldmanCousins.cxx:121
 TFeldmanCousins.cxx:122
 TFeldmanCousins.cxx:123
 TFeldmanCousins.cxx:124
 TFeldmanCousins.cxx:125
 TFeldmanCousins.cxx:126
 TFeldmanCousins.cxx:127
 TFeldmanCousins.cxx:128
 TFeldmanCousins.cxx:129
 TFeldmanCousins.cxx:130
 TFeldmanCousins.cxx:131
 TFeldmanCousins.cxx:132
 TFeldmanCousins.cxx:133
 TFeldmanCousins.cxx:134
 TFeldmanCousins.cxx:135
 TFeldmanCousins.cxx:136
 TFeldmanCousins.cxx:137
 TFeldmanCousins.cxx:138
 TFeldmanCousins.cxx:139
 TFeldmanCousins.cxx:140
 TFeldmanCousins.cxx:141
 TFeldmanCousins.cxx:142
 TFeldmanCousins.cxx:143
 TFeldmanCousins.cxx:144
 TFeldmanCousins.cxx:145
 TFeldmanCousins.cxx:146
 TFeldmanCousins.cxx:147
 TFeldmanCousins.cxx:148
 TFeldmanCousins.cxx:149
 TFeldmanCousins.cxx:150
 TFeldmanCousins.cxx:151
 TFeldmanCousins.cxx:152
 TFeldmanCousins.cxx:153
 TFeldmanCousins.cxx:154
 TFeldmanCousins.cxx:155
 TFeldmanCousins.cxx:156
 TFeldmanCousins.cxx:157
 TFeldmanCousins.cxx:158
 TFeldmanCousins.cxx:159
 TFeldmanCousins.cxx:160
 TFeldmanCousins.cxx:161
 TFeldmanCousins.cxx:162
 TFeldmanCousins.cxx:163
 TFeldmanCousins.cxx:164
 TFeldmanCousins.cxx:165
 TFeldmanCousins.cxx:166
 TFeldmanCousins.cxx:167
 TFeldmanCousins.cxx:168
 TFeldmanCousins.cxx:169
 TFeldmanCousins.cxx:170
 TFeldmanCousins.cxx:171
 TFeldmanCousins.cxx:172
 TFeldmanCousins.cxx:173
 TFeldmanCousins.cxx:174
 TFeldmanCousins.cxx:175
 TFeldmanCousins.cxx:176
 TFeldmanCousins.cxx:177
 TFeldmanCousins.cxx:178
 TFeldmanCousins.cxx:179
 TFeldmanCousins.cxx:180
 TFeldmanCousins.cxx:181
 TFeldmanCousins.cxx:182
 TFeldmanCousins.cxx:183
 TFeldmanCousins.cxx:184
 TFeldmanCousins.cxx:185
 TFeldmanCousins.cxx:186
 TFeldmanCousins.cxx:187
 TFeldmanCousins.cxx:188
 TFeldmanCousins.cxx:189
 TFeldmanCousins.cxx:190
 TFeldmanCousins.cxx:191
 TFeldmanCousins.cxx:192
 TFeldmanCousins.cxx:193
 TFeldmanCousins.cxx:194
 TFeldmanCousins.cxx:195
 TFeldmanCousins.cxx:196
 TFeldmanCousins.cxx:197
 TFeldmanCousins.cxx:198
 TFeldmanCousins.cxx:199
 TFeldmanCousins.cxx:200
 TFeldmanCousins.cxx:201
 TFeldmanCousins.cxx:202
 TFeldmanCousins.cxx:203
 TFeldmanCousins.cxx:204
 TFeldmanCousins.cxx:205
 TFeldmanCousins.cxx:206
 TFeldmanCousins.cxx:207
 TFeldmanCousins.cxx:208
 TFeldmanCousins.cxx:209
 TFeldmanCousins.cxx:210
 TFeldmanCousins.cxx:211
 TFeldmanCousins.cxx:212
 TFeldmanCousins.cxx:213
 TFeldmanCousins.cxx:214
 TFeldmanCousins.cxx:215
 TFeldmanCousins.cxx:216
 TFeldmanCousins.cxx:217
 TFeldmanCousins.cxx:218
 TFeldmanCousins.cxx:219
 TFeldmanCousins.cxx:220
 TFeldmanCousins.cxx:221
 TFeldmanCousins.cxx:222
 TFeldmanCousins.cxx:223
 TFeldmanCousins.cxx:224
 TFeldmanCousins.cxx:225
 TFeldmanCousins.cxx:226
 TFeldmanCousins.cxx:227
 TFeldmanCousins.cxx:228
 TFeldmanCousins.cxx:229
 TFeldmanCousins.cxx:230
 TFeldmanCousins.cxx:231
 TFeldmanCousins.cxx:232
 TFeldmanCousins.cxx:233
 TFeldmanCousins.cxx:234