// @(#)root/roostats:$Id$
// Author: Kyle Cranmer, Lorenzo Moneta, Gregory Schott, Wouter Verkerke
/*************************************************************************
 * Copyright (C) 1995-2008, Rene Brun and Fons Rademakers.               *
 * All rights reserved.                                                  *
 *                                                                       *
 * For the licensing terms see $ROOTSYS/LICENSE.                         *
 * For the list of contributors see $ROOTSYS/README/CREDITS.             *
 *************************************************************************/

/*****************************************************************************
 * Project: RooStats
 * Package: RooFit/RooStats  
 * @(#)root/roofit/roostats:$Id$
 * Original Author: Kyle Cranmer
 *   Kyle Cranmer, Lorenzo Moneta, Gregory Schott, Wouter Verkerke
 *
 *****************************************************************************/



//_________________________________________________________
//
// BEGIN_HTML
// ConfidenceBelt is a concrete implementation of the ConfInterval interface.  
// It implements simple general purpose interval of arbitrary dimensions and shape.
// It does not assume the interval is connected.
// It uses either a RooDataSet (eg. a list of parameter points in the interval) or
// a RooDataHist (eg. a Histogram-like object for small regions of the parameter space) to
// store the interval.  
// END_HTML
//
//

#ifndef RooStats_ConfidenceBelt
#include "RooStats/ConfidenceBelt.h"
#endif

#include "RooDataSet.h"
#include "RooDataHist.h"

#include "RooStats/RooStatsUtils.h"

ClassImp(RooStats::ConfidenceBelt) ;

using namespace RooStats;
using namespace std;

//____________________________________________________________________
ConfidenceBelt::ConfidenceBelt() : 
   TNamed(), fParameterPoints(0)
{
   // Default constructor
}

//____________________________________________________________________
ConfidenceBelt::ConfidenceBelt(const char* name) :
  TNamed(name,name), fParameterPoints(0)
{
   // Alternate constructor
}

//____________________________________________________________________
ConfidenceBelt::ConfidenceBelt(const char* name, const char* title) :
   TNamed(name,title), fParameterPoints(0)
{
   // Alternate constructor
}

//____________________________________________________________________
ConfidenceBelt::ConfidenceBelt(const char* name, RooAbsData& data) :
  TNamed(name,name), fParameterPoints((RooAbsData*)data.Clone("PointsToTestForBelt"))
{
   // Alternate constructor
}

//____________________________________________________________________
ConfidenceBelt::ConfidenceBelt(const char* name, const char* title, RooAbsData& data) :
   TNamed(name,title), fParameterPoints((RooAbsData*)data.Clone("PointsToTestForBelt"))
{
   // Alternate constructor
}



//____________________________________________________________________
ConfidenceBelt::~ConfidenceBelt()
{
   // Destructor

}


//____________________________________________________________________
Double_t ConfidenceBelt::GetAcceptanceRegionMin(RooArgSet& parameterPoint, Double_t cl, Double_t leftside) {

  if(cl>0 || leftside > 0) cout <<"using default cl, leftside for now" <<endl;
  AcceptanceRegion * region = GetAcceptanceRegion(parameterPoint, cl,leftside); 
  return (region) ? region->GetLowerLimit() : TMath::QuietNaN();
}

//____________________________________________________________________
Double_t ConfidenceBelt::GetAcceptanceRegionMax(RooArgSet& parameterPoint, Double_t cl, Double_t leftside) {
  if(cl>0 || leftside > 0) cout <<"using default cl, leftside for now" <<endl;
  AcceptanceRegion * region = GetAcceptanceRegion(parameterPoint, cl,leftside); 
  return (region) ? region->GetUpperLimit() : TMath::QuietNaN();
}

//____________________________________________________________________
vector<Double_t> ConfidenceBelt::ConfidenceLevels() const {
  vector<Double_t> levels;
  return levels;
}

//____________________________________________________________________
void ConfidenceBelt::AddAcceptanceRegion(RooArgSet& parameterPoint, Int_t dsIndex,
					 Double_t lower, Double_t upper,
					 Double_t cl, Double_t leftside){
  
  if(cl>0 || leftside > 0) cout <<"using default cl, leftside for now" <<endl;

  RooDataSet*  tree = dynamic_cast<RooDataSet*>(  fParameterPoints );
  RooDataHist* hist = dynamic_cast<RooDataHist*>( fParameterPoints );  

  //  cout << "add: " << tree << " " << hist << endl;

  if( !this->CheckParameters(parameterPoint) )
    std::cout << "problem with parameters" << std::endl;

  Int_t luIndex = fSamplingSummaryLookup.GetLookupIndex(cl, leftside);
  //  cout << "lookup index = " << luIndex << endl;
  if(luIndex <0 ) {
    fSamplingSummaryLookup.Add(cl,leftside);
    luIndex = fSamplingSummaryLookup.GetLookupIndex(cl, leftside);
    cout << "lookup index = " << luIndex << endl;
  }
  AcceptanceRegion* thisRegion = new AcceptanceRegion(luIndex, lower, upper);
  
  if( hist ) {
    // need a way to get index for given point
    // Can do this by setting hist's internal parameters to desired values
    // need a better way
    //    RooStats::SetParameters(&parameterPoint, const_cast<RooArgSet*>(hist->get())); 
    //    int index = hist->calcTreeIndex(); // get index
    int index = hist->getIndex(parameterPoint); // get index
    //    cout << "hist index = " << index << endl;

    // allocate memory if necessary.  numEntries is overkill?
    if((Int_t)fSamplingSummaries.size() <= index) {
      fSamplingSummaries.reserve( hist->numEntries() ); 
      fSamplingSummaries.resize( hist->numEntries() ); 
    }

    // set the region for this point (check for duplicate?)
    fSamplingSummaries.at(index) = *thisRegion;
  }
  else if( tree ){
    //    int index = tree->getIndex(parameterPoint); 
    int index = dsIndex;
    //    cout << "tree index = " << index << endl;

    // allocate memory if necessary.  numEntries is overkill?
    if((Int_t)fSamplingSummaries.size() <= index){
      fSamplingSummaries.reserve( tree->numEntries()  ); 
      fSamplingSummaries.resize( tree->numEntries() ); 
    }

    // set the region for this point (check for duplicate?)
    fSamplingSummaries.at( index ) = *thisRegion;
  }
}

//____________________________________________________________________
void ConfidenceBelt::AddAcceptanceRegion(RooArgSet& parameterPoint, AcceptanceRegion region, 
					 Double_t cl, Double_t leftside){
  
  if(cl>0 || leftside > 0) cout <<"using default cl, leftside for now" <<endl;

  RooDataSet*  tree = dynamic_cast<RooDataSet*>(  fParameterPoints );
  RooDataHist* hist = dynamic_cast<RooDataHist*>( fParameterPoints );

  if( !this->CheckParameters(parameterPoint) )
    std::cout << "problem with parameters" << std::endl;
  
  
  if( hist ) {
    // need a way to get index for given point
    // Can do this by setting hist's internal parameters to desired values
    // need a better way
    //    RooStats::SetParameters(&parameterPoint, const_cast<RooArgSet*>(hist->get())); 
    //    int index = hist->calcTreeIndex(); // get index
    int index = hist->getIndex(parameterPoint); // get index

    // allocate memory if necessary.  numEntries is overkill?
    if((Int_t)fSamplingSummaries.size() < index) fSamplingSummaries.reserve( hist->numEntries() ); 

    // set the region for this point (check for duplicate?)
    fSamplingSummaries.at(index) = region;
  }
  else if( tree ){
    tree->add( parameterPoint ); // assume it's unique for now
    int index = tree->numEntries() - 1; //check that last point added has index nEntries -1
    // allocate memory if necessary.  numEntries is overkill?
    if((Int_t)fSamplingSummaries.size() < index) fSamplingSummaries.reserve( tree->numEntries()  ); 

    // set the region for this point (check for duplicate?)
    fSamplingSummaries.at( index ) = region;
  }
}

//____________________________________________________________________
AcceptanceRegion* ConfidenceBelt::GetAcceptanceRegion(RooArgSet &parameterPoint, Double_t cl, Double_t leftside) 
{  
   // Method to determine if a parameter point is in the interval

  if(cl>0 || leftside > 0) cout <<"using default cl, leftside for now" <<endl;

  RooDataSet*  tree = dynamic_cast<RooDataSet*>(  fParameterPoints );
  RooDataHist* hist = dynamic_cast<RooDataHist*>( fParameterPoints );
  
  if( !this->CheckParameters(parameterPoint) ){
    std::cout << "problem with parameters" << std::endl;
    return 0; 
  }
  
  if( hist ) {
    // need a way to get index for given point
    // Can do this by setting hist's internal parameters to desired values
    // need a better way
    //    RooStats::SetParameters(&parameterPoint, const_cast<RooArgSet*>(hist->get())); 
    //    Int_t index = hist->calcTreeIndex(); // get index
    int index = hist->getIndex(parameterPoint); // get index
    return &(fSamplingSummaries.at(index).GetAcceptanceRegion());
  }
  else if( tree ){
    // need a way to get index for given point
    //    RooStats::SetParameters(&parameterPoint, tree->get()); // set tree's parameters to desired values
    Int_t index = 0; //need something like tree->calcTreeIndex(); 
    const RooArgSet* thisPoint = 0;
    for(index=0; index<tree->numEntries(); ++index){
      thisPoint = tree->get(index); 
      bool samePoint = true;
      TIter it = parameterPoint.createIterator();
      RooRealVar *myarg; 
      while ( samePoint && (myarg = (RooRealVar *)it.Next())) { 
	if(myarg->getVal() != thisPoint->getRealValue(myarg->GetName()))
	  samePoint = false;
      }
      if(samePoint) 
	break;
    }

    return &(fSamplingSummaries.at(index).GetAcceptanceRegion());
  }
  else {
      std::cout << "dataset is not initialized properly" << std::endl;
  }

  return 0;
  
}

//____________________________________________________________________
RooArgSet* ConfidenceBelt::GetParameters() const
{  
   // returns list of parameters
   return new RooArgSet(*(fParameterPoints->get()));
}

//____________________________________________________________________
Bool_t ConfidenceBelt::CheckParameters(RooArgSet &parameterPoint) const
{  

   if (parameterPoint.getSize() != fParameterPoints->get()->getSize() ) {
      std::cout << "size is wrong, parameters don't match" << std::endl;
      return false;
   }
   if ( ! parameterPoint.equals( *(fParameterPoints->get() ) ) ) {
      std::cout << "size is ok, but parameters don't match" << std::endl;
      return false;
   }
   return true;
}
 ConfidenceBelt.cxx:1
 ConfidenceBelt.cxx:2
 ConfidenceBelt.cxx:3
 ConfidenceBelt.cxx:4
 ConfidenceBelt.cxx:5
 ConfidenceBelt.cxx:6
 ConfidenceBelt.cxx:7
 ConfidenceBelt.cxx:8
 ConfidenceBelt.cxx:9
 ConfidenceBelt.cxx:10
 ConfidenceBelt.cxx:11
 ConfidenceBelt.cxx:12
 ConfidenceBelt.cxx:13
 ConfidenceBelt.cxx:14
 ConfidenceBelt.cxx:15
 ConfidenceBelt.cxx:16
 ConfidenceBelt.cxx:17
 ConfidenceBelt.cxx:18
 ConfidenceBelt.cxx:19
 ConfidenceBelt.cxx:20
 ConfidenceBelt.cxx:21
 ConfidenceBelt.cxx:22
 ConfidenceBelt.cxx:23
 ConfidenceBelt.cxx:24
 ConfidenceBelt.cxx:25
 ConfidenceBelt.cxx:26
 ConfidenceBelt.cxx:27
 ConfidenceBelt.cxx:28
 ConfidenceBelt.cxx:29
 ConfidenceBelt.cxx:30
 ConfidenceBelt.cxx:31
 ConfidenceBelt.cxx:32
 ConfidenceBelt.cxx:33
 ConfidenceBelt.cxx:34
 ConfidenceBelt.cxx:35
 ConfidenceBelt.cxx:36
 ConfidenceBelt.cxx:37
 ConfidenceBelt.cxx:38
 ConfidenceBelt.cxx:39
 ConfidenceBelt.cxx:40
 ConfidenceBelt.cxx:41
 ConfidenceBelt.cxx:42
 ConfidenceBelt.cxx:43
 ConfidenceBelt.cxx:44
 ConfidenceBelt.cxx:45
 ConfidenceBelt.cxx:46
 ConfidenceBelt.cxx:47
 ConfidenceBelt.cxx:48
 ConfidenceBelt.cxx:49
 ConfidenceBelt.cxx:50
 ConfidenceBelt.cxx:51
 ConfidenceBelt.cxx:52
 ConfidenceBelt.cxx:53
 ConfidenceBelt.cxx:54
 ConfidenceBelt.cxx:55
 ConfidenceBelt.cxx:56
 ConfidenceBelt.cxx:57
 ConfidenceBelt.cxx:58
 ConfidenceBelt.cxx:59
 ConfidenceBelt.cxx:60
 ConfidenceBelt.cxx:61
 ConfidenceBelt.cxx:62
 ConfidenceBelt.cxx:63
 ConfidenceBelt.cxx:64
 ConfidenceBelt.cxx:65
 ConfidenceBelt.cxx:66
 ConfidenceBelt.cxx:67
 ConfidenceBelt.cxx:68
 ConfidenceBelt.cxx:69
 ConfidenceBelt.cxx:70
 ConfidenceBelt.cxx:71
 ConfidenceBelt.cxx:72
 ConfidenceBelt.cxx:73
 ConfidenceBelt.cxx:74
 ConfidenceBelt.cxx:75
 ConfidenceBelt.cxx:76
 ConfidenceBelt.cxx:77
 ConfidenceBelt.cxx:78
 ConfidenceBelt.cxx:79
 ConfidenceBelt.cxx:80
 ConfidenceBelt.cxx:81
 ConfidenceBelt.cxx:82
 ConfidenceBelt.cxx:83
 ConfidenceBelt.cxx:84
 ConfidenceBelt.cxx:85
 ConfidenceBelt.cxx:86
 ConfidenceBelt.cxx:87
 ConfidenceBelt.cxx:88
 ConfidenceBelt.cxx:89
 ConfidenceBelt.cxx:90
 ConfidenceBelt.cxx:91
 ConfidenceBelt.cxx:92
 ConfidenceBelt.cxx:93
 ConfidenceBelt.cxx:94
 ConfidenceBelt.cxx:95
 ConfidenceBelt.cxx:96
 ConfidenceBelt.cxx:97
 ConfidenceBelt.cxx:98
 ConfidenceBelt.cxx:99
 ConfidenceBelt.cxx:100
 ConfidenceBelt.cxx:101
 ConfidenceBelt.cxx:102
 ConfidenceBelt.cxx:103
 ConfidenceBelt.cxx:104
 ConfidenceBelt.cxx:105
 ConfidenceBelt.cxx:106
 ConfidenceBelt.cxx:107
 ConfidenceBelt.cxx:108
 ConfidenceBelt.cxx:109
 ConfidenceBelt.cxx:110
 ConfidenceBelt.cxx:111
 ConfidenceBelt.cxx:112
 ConfidenceBelt.cxx:113
 ConfidenceBelt.cxx:114
 ConfidenceBelt.cxx:115
 ConfidenceBelt.cxx:116
 ConfidenceBelt.cxx:117
 ConfidenceBelt.cxx:118
 ConfidenceBelt.cxx:119
 ConfidenceBelt.cxx:120
 ConfidenceBelt.cxx:121
 ConfidenceBelt.cxx:122
 ConfidenceBelt.cxx:123
 ConfidenceBelt.cxx:124
 ConfidenceBelt.cxx:125
 ConfidenceBelt.cxx:126
 ConfidenceBelt.cxx:127
 ConfidenceBelt.cxx:128
 ConfidenceBelt.cxx:129
 ConfidenceBelt.cxx:130
 ConfidenceBelt.cxx:131
 ConfidenceBelt.cxx:132
 ConfidenceBelt.cxx:133
 ConfidenceBelt.cxx:134
 ConfidenceBelt.cxx:135
 ConfidenceBelt.cxx:136
 ConfidenceBelt.cxx:137
 ConfidenceBelt.cxx:138
 ConfidenceBelt.cxx:139
 ConfidenceBelt.cxx:140
 ConfidenceBelt.cxx:141
 ConfidenceBelt.cxx:142
 ConfidenceBelt.cxx:143
 ConfidenceBelt.cxx:144
 ConfidenceBelt.cxx:145
 ConfidenceBelt.cxx:146
 ConfidenceBelt.cxx:147
 ConfidenceBelt.cxx:148
 ConfidenceBelt.cxx:149
 ConfidenceBelt.cxx:150
 ConfidenceBelt.cxx:151
 ConfidenceBelt.cxx:152
 ConfidenceBelt.cxx:153
 ConfidenceBelt.cxx:154
 ConfidenceBelt.cxx:155
 ConfidenceBelt.cxx:156
 ConfidenceBelt.cxx:157
 ConfidenceBelt.cxx:158
 ConfidenceBelt.cxx:159
 ConfidenceBelt.cxx:160
 ConfidenceBelt.cxx:161
 ConfidenceBelt.cxx:162
 ConfidenceBelt.cxx:163
 ConfidenceBelt.cxx:164
 ConfidenceBelt.cxx:165
 ConfidenceBelt.cxx:166
 ConfidenceBelt.cxx:167
 ConfidenceBelt.cxx:168
 ConfidenceBelt.cxx:169
 ConfidenceBelt.cxx:170
 ConfidenceBelt.cxx:171
 ConfidenceBelt.cxx:172
 ConfidenceBelt.cxx:173
 ConfidenceBelt.cxx:174
 ConfidenceBelt.cxx:175
 ConfidenceBelt.cxx:176
 ConfidenceBelt.cxx:177
 ConfidenceBelt.cxx:178
 ConfidenceBelt.cxx:179
 ConfidenceBelt.cxx:180
 ConfidenceBelt.cxx:181
 ConfidenceBelt.cxx:182
 ConfidenceBelt.cxx:183
 ConfidenceBelt.cxx:184
 ConfidenceBelt.cxx:185
 ConfidenceBelt.cxx:186
 ConfidenceBelt.cxx:187
 ConfidenceBelt.cxx:188
 ConfidenceBelt.cxx:189
 ConfidenceBelt.cxx:190
 ConfidenceBelt.cxx:191
 ConfidenceBelt.cxx:192
 ConfidenceBelt.cxx:193
 ConfidenceBelt.cxx:194
 ConfidenceBelt.cxx:195
 ConfidenceBelt.cxx:196
 ConfidenceBelt.cxx:197
 ConfidenceBelt.cxx:198
 ConfidenceBelt.cxx:199
 ConfidenceBelt.cxx:200
 ConfidenceBelt.cxx:201
 ConfidenceBelt.cxx:202
 ConfidenceBelt.cxx:203
 ConfidenceBelt.cxx:204
 ConfidenceBelt.cxx:205
 ConfidenceBelt.cxx:206
 ConfidenceBelt.cxx:207
 ConfidenceBelt.cxx:208
 ConfidenceBelt.cxx:209
 ConfidenceBelt.cxx:210
 ConfidenceBelt.cxx:211
 ConfidenceBelt.cxx:212
 ConfidenceBelt.cxx:213
 ConfidenceBelt.cxx:214
 ConfidenceBelt.cxx:215
 ConfidenceBelt.cxx:216
 ConfidenceBelt.cxx:217
 ConfidenceBelt.cxx:218
 ConfidenceBelt.cxx:219
 ConfidenceBelt.cxx:220
 ConfidenceBelt.cxx:221
 ConfidenceBelt.cxx:222
 ConfidenceBelt.cxx:223
 ConfidenceBelt.cxx:224
 ConfidenceBelt.cxx:225
 ConfidenceBelt.cxx:226
 ConfidenceBelt.cxx:227
 ConfidenceBelt.cxx:228
 ConfidenceBelt.cxx:229
 ConfidenceBelt.cxx:230
 ConfidenceBelt.cxx:231
 ConfidenceBelt.cxx:232
 ConfidenceBelt.cxx:233
 ConfidenceBelt.cxx:234
 ConfidenceBelt.cxx:235
 ConfidenceBelt.cxx:236
 ConfidenceBelt.cxx:237
 ConfidenceBelt.cxx:238
 ConfidenceBelt.cxx:239
 ConfidenceBelt.cxx:240
 ConfidenceBelt.cxx:241
 ConfidenceBelt.cxx:242
 ConfidenceBelt.cxx:243
 ConfidenceBelt.cxx:244
 ConfidenceBelt.cxx:245
 ConfidenceBelt.cxx:246
 ConfidenceBelt.cxx:247
 ConfidenceBelt.cxx:248
 ConfidenceBelt.cxx:249
 ConfidenceBelt.cxx:250
 ConfidenceBelt.cxx:251
 ConfidenceBelt.cxx:252
 ConfidenceBelt.cxx:253
 ConfidenceBelt.cxx:254
 ConfidenceBelt.cxx:255
 ConfidenceBelt.cxx:256
 ConfidenceBelt.cxx:257
 ConfidenceBelt.cxx:258
 ConfidenceBelt.cxx:259
 ConfidenceBelt.cxx:260
 ConfidenceBelt.cxx:261
 ConfidenceBelt.cxx:262
 ConfidenceBelt.cxx:263
 ConfidenceBelt.cxx:264
 ConfidenceBelt.cxx:265
 ConfidenceBelt.cxx:266
 ConfidenceBelt.cxx:267
 ConfidenceBelt.cxx:268
 ConfidenceBelt.cxx:269
 ConfidenceBelt.cxx:270
 ConfidenceBelt.cxx:271
 ConfidenceBelt.cxx:272
 ConfidenceBelt.cxx:273
 ConfidenceBelt.cxx:274
 ConfidenceBelt.cxx:275
 ConfidenceBelt.cxx:276
 ConfidenceBelt.cxx:277
 ConfidenceBelt.cxx:278
 ConfidenceBelt.cxx:279
 ConfidenceBelt.cxx:280
 ConfidenceBelt.cxx:281
 ConfidenceBelt.cxx:282
 ConfidenceBelt.cxx:283