/***************************************************************************** 
  * Project: RooFit                                                           * 
  *                                                                           * 
  * Simple Poisson PDF
  * author: Kyle Cranmer <cranmer@cern.ch>
  *                                                                           * 
  *****************************************************************************/ 

//////////////////////////////////////////////////////////////////////////////
//
// BEGIN_HTML
// Poisson pdf
// END_HTML
//

#include <iostream> 

#include "RooPoisson.h" 
#include "RooAbsReal.h" 
#include "RooAbsCategory.h" 

#include "RooRandom.h"
#include "RooMath.h"
#include "TMath.h"
#include "Math/ProbFuncMathCore.h"

#include "TError.h"

using namespace std;

ClassImp(RooPoisson) 



//_____________________________________________________________________________
RooPoisson::RooPoisson(const char *name, const char *title, 
		       RooAbsReal& _x,
		       RooAbsReal& _mean,
		       Bool_t noRounding) :
  RooAbsPdf(name,title), 
  x("x","x",this,_x),
  mean("mean","mean",this,_mean),
  _noRounding(noRounding),
  _protectNegative(false)
{ 
  // Constructor  
} 



//_____________________________________________________________________________
 RooPoisson::RooPoisson(const RooPoisson& other, const char* name) :  
   RooAbsPdf(other,name), 
   x("x",this,other.x),
   mean("mean",this,other.mean),
   _noRounding(other._noRounding),
   _protectNegative(other._protectNegative)
{ 
   // Copy constructor
} 




//_____________________________________________________________________________
Double_t RooPoisson::evaluate() const 
{ 
  // Implementation in terms of the TMath Poisson function

  Double_t k = _noRounding ? x : floor(x);  
  if(_protectNegative && mean<0) 
    return 1e-3;
  return TMath::Poisson(k,mean) ;
} 



//_____________________________________________________________________________
Double_t RooPoisson::getLogVal(const RooArgSet* s) const 
{
  // calculate and return the negative log-likelihood of the Poisson                                                                                                                                    
  return RooAbsPdf::getLogVal(s) ;
//   Double_t prob = getVal(s) ;
//   return prob ;

  // Make inputs to naming conventions of RooAbsPdf::extendedTerm
  Double_t expected=mean ;
  Double_t observed=x ;

  // Explicitly handle case Nobs=Nexp=0
  if (fabs(expected)<1e-10 && fabs(observed)<1e-10) {
    return 0 ;
  }  

  // Explicitly handle case Nexp=0
  if (fabs(observed)<1e-10) {
    return -1*expected;
  }

  // Michaels code for log(poisson) in RooAbsPdf::extendedTer with an approximated log(observed!) term
  Double_t extra=0;
  if(observed<1000000) {
    extra = -observed*log(expected)+expected+TMath::LnGamma(observed+1.);    
  } else {
    //if many observed events, use Gauss approximation                                                                                                                                                 
    Double_t sigma_square=expected;
    Double_t diff=observed-expected;
    extra=-log(sigma_square)/2 + (diff*diff)/(2*sigma_square);
  }

//   if (fabs(extra)>100 || log(prob)>100) {
//     cout << "RooPoisson::getLogVal(" << GetName() << ") mu=" << expected << " x = " << x << " -log(P) = " << extra << " log(evaluate()) = " << log(prob) << endl ;
//   }
  
//   if (fabs(extra+log(prob))>1) {
//     cout << "RooPoisson::getLogVal(" << GetName() << ") WARNING mu=" << expected << " x = " << x << " -log(P) = " << extra << " log(evaluate()) = " << log(prob) << endl ;
//   }

  //return log(prob);
  return -extra-analyticalIntegral(1,0) ; //log(prob);

}


//_____________________________________________________________________________
Int_t RooPoisson::getAnalyticalIntegral(RooArgSet& allVars, RooArgSet& analVars, const char* /*rangeName*/) const 
{
  if (matchArgs(allVars,analVars,x)) return 1 ;
  if (matchArgs(allVars, analVars, mean)) return 2;
  return 0 ;
}



//_____________________________________________________________________________
Double_t RooPoisson::analyticalIntegral(Int_t code, const char* rangeName) const 
{
  R__ASSERT(code == 1 || code == 2) ;

  if(_protectNegative && mean<0) 
    return exp(-2*mean); // make it fall quickly

  if (code == 1) {
    // Implement integral over x as summation. Add special handling in case
    // range boundaries are not on integer values of x
    Double_t xmin = x.min(rangeName) ;
    Double_t xmax = x.max(rangeName) ;

    // Protect against negative lower boundaries
    if (xmin<0) xmin=0 ;

    Int_t ixmin = Int_t (xmin) ;
    Int_t ixmax = Int_t (xmax)+1 ;

    Double_t fracLoBin = 1-(xmin-ixmin) ;
    Double_t fracHiBin = 1-(ixmax-xmax) ;

    if (!x.hasMax()) {
      if (xmin<1e-6) {
        return 1 ;
      } else {
      
      // Return 1 minus integral from 0 to x.min() 

        if(ixmin == 0){ // first bin
	       return TMath::Poisson(0, mean)*(xmin-0);
        }      
        Double_t sum(0) ;
        sum += TMath::Poisson(0,mean)*fracLoBin ;
        sum+= ROOT::Math::poisson_cdf(ixmin-2, mean) - ROOT::Math::poisson_cdf(0,mean) ;
        sum += TMath::Poisson(ixmin-1,mean)*fracHiBin ;
        return 1-sum ;
      }
    }
  
    if(ixmin == ixmax-1){ // first bin
      return TMath::Poisson(ixmin, mean)*(xmax-xmin);
    }  

    Double_t sum(0) ;
    sum += TMath::Poisson(ixmin,mean)*fracLoBin ;
    if (RooNumber::isInfinite(xmax)){
      sum+= 1.-ROOT::Math::poisson_cdf(ixmin,mean) ;
    }  else {
      sum+= ROOT::Math::poisson_cdf(ixmax-2, mean) - ROOT::Math::poisson_cdf(ixmin,mean) ;
      sum += TMath::Poisson(ixmax-1,mean)*fracHiBin ;
    }
  
    return sum ;

  } else if(code == 2) {
   
    // the integral with respect to the mean is the integral of a gamma distribution 
    Double_t mean_min = mean.min(rangeName);
    Double_t mean_max = mean.max(rangeName);
     
    Double_t ix;
    if(_noRounding) ix = x + 1;
    else ix = Int_t(TMath::Floor(x)) + 1.0; // negative ix does not need protection (gamma returns 0.0)
    
    return ROOT::Math::gamma_cdf(mean_max, ix, 1.0) - ROOT::Math::gamma_cdf(mean_min, ix, 1.0);
  }

  return 0;

}







//_____________________________________________________________________________
Int_t RooPoisson::getGenerator(const RooArgSet& directVars, RooArgSet &generateVars, Bool_t /*staticInitOK*/) const
{
  // Advertise internal generator in x

  if (matchArgs(directVars,generateVars,x)) return 1 ;  
  return 0 ;
}



//_____________________________________________________________________________
void RooPoisson::generateEvent(Int_t code)
{
  // Implement internal generator using TRandom::Poisson 

  R__ASSERT(code==1) ;
  Double_t xgen ;
  while(1) {    
    xgen = RooRandom::randomGenerator()->Poisson(mean);
    if (xgen<=x.max() && xgen>=x.min()) {
      x = xgen ;
      break;
    }
  }
  return;
}


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