// @(#)root/roostats:$Id$
// Author: Kyle Cranmer, Lorenzo Moneta, Gregory Schott, Wouter Verkerke, Nils Ruthmann
/*************************************************************************
 * 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.             *
 *************************************************************************/




//////////////////////////////////////////////////////////////////////////////
// 
// BEGIN_HTML
// This modules allow to compute in the ToyMcStudy framework the ProfileLikelihood upper limit for each toy-MC sample generated

// END_HTML
//
//

#include "Riostream.h"

#include "RooDataSet.h"
//#include "RooRealVar.h"
#include "TString.h"
//#include "RooFit.h"
#include "RooFitResult.h"
#include "RooStats/UpperLimitMCSModule.h"
#include "RooMsgService.h"
#include "RooStats/ConfInterval.h"
#include "RooStats/PointSetInterval.h"
#include "RooStats/LikelihoodInterval.h"
#include "RooStats/LikelihoodIntervalPlot.h"
#include "RooStats/ProfileLikelihoodCalculator.h"
#include "TCanvas.h"
#include "RooMinuit.h"
#include "RooNLLVar.h"
#include "RooCmdArg.h"
#include "RooRealVar.h"



using namespace std;

ClassImp(RooStats::UpperLimitMCSModule);


using namespace RooStats ;



//_____________________________________________________________________________
UpperLimitMCSModule::UpperLimitMCSModule(const RooArgSet* poi, Double_t CL) : 
  RooAbsMCStudyModule(Form("UpperLimitMCSModule_%s",poi->first()->GetName()),Form("UpperLimitMCSModule_%s",poi->first()->GetName())),
  _parName(poi->first()->GetName()), 
  _plc(0),_ul(0),_poi(0), _data(0),_cl(CL), _model(0)
{
  std::cout<<"RooUpperLimitConstructor ParName:"<<_parName<<std::endl;
  std::cout<<"RooUpperLimitConstructor CL:"<<_cl<<std::endl;
  // Constructor of module with parameter to be interpreted as nSignal and the value of the
  // null hypothesis for nSignal (usually zero)
}



//_____________________________________________________________________________
UpperLimitMCSModule::UpperLimitMCSModule(const UpperLimitMCSModule& other) : 
  RooAbsMCStudyModule(other), 
  _parName(other._poi->first()->GetName()),
  _plc(0),_ul(0),_poi(other._poi), _data(0), _cl(other._cl), _model(other._model)
{
  // Copy constructor
}



//_____________________________________________________________________________
UpperLimitMCSModule:: ~UpperLimitMCSModule() 
{
  // Destructor

 
  if (_plc) {
    delete _plc ;
  }
  if (_data) {
    delete _data ;
  }
  if(_ul){
    delete _ul;
  }
  if(_poi){
     delete _poi;
  }
  if (_model){
    delete _model;
  }
}



//_____________________________________________________________________________
Bool_t UpperLimitMCSModule::initializeInstance()
{
  // Initialize module after attachment to RooMCStudy object

  // Check that parameter is also present in fit parameter list of RooMCStudy object
  if (!fitParams()->find(_parName.c_str())) {
    coutE(InputArguments) << "UpperLimitMCSModule::initializeInstance:: ERROR: No parameter named " << _parName << " in RooMCStudy!" << endl ;
    return kFALSE ;
  }
  
  //Construct the ProfileLikelihoodCalculator
  _poi=new RooArgSet(*(fitParams()->find(_parName.c_str())));
  std::cout<<"RooUpperLimit Initialize Instance: POI Set:"<<std::endl;
  _poi->Print("v");
  std::cout<<"RooUpperLimit Initialize Instance: End:"<<std::endl;


  
  TString ulName = Form("ul_%s",_parName.c_str()) ;
  TString ulTitle = Form("UL for parameter %s",_parName.c_str()) ;
  _ul = new RooRealVar(ulName.Data(),ulTitle.Data(),0) ;


  // Create new dataset to be merged with RooMCStudy::fitParDataSet
  _data = new RooDataSet("ULSigData","Additional data for UL study",RooArgSet(*_ul)) ;

  return kTRUE ;
}



//_____________________________________________________________________________
Bool_t UpperLimitMCSModule::initializeRun(Int_t /*numSamples*/) 
{
  // Initialize module at beginning of RooCMStudy run

  _data->reset() ;
  return kTRUE ;
}



//_____________________________________________________________________________
RooDataSet* UpperLimitMCSModule::finalizeRun() 
{
  // Return auxiliary dataset with results of delta(-log(L))
  // calculations of this module so that it is merged with
  // RooMCStudy::fitParDataSet() by RooMCStudy

  return _data ;
}



//_____________________________________________________________________________

// Bool_t UpperLimitMCSModule::processAfterFit(Int_t /*sampleNum*/)  
// {
//   // Save likelihood from nominal fit, fix chosen parameter to its
//   // null hypothesis value and rerun fit Save difference in likelihood
//   // and associated Gaussian significance in auxilary dataset

//   RooRealVar* par = static_cast<RooRealVar*>(fitParams()->find(_parName.c_str())) ;
//   par->setVal(_nullValue) ;
//   par->setConstant(kTRUE) ;
//   RooFitResult* frnull = refit() ;
//   par->setConstant(kFALSE) ;
  
//   _nll0h->setVal(frnull->minNll()) ;

//   Double_t deltaLL = (frnull->minNll() - nllVar()->getVal()) ;
//   Double_t signif = deltaLL>0 ? sqrt(2*deltaLL) : -sqrt(-2*deltaLL) ;
//   _sig0h->setVal(signif) ;
//   _dll0h->setVal(deltaLL) ;


//   _data->add(RooArgSet(*_nll0h,*_dll0h,*_sig0h)) ;

//   delete frnull ;
//   return kTRUE ;

// }


//_____________________________________________________________________________
Bool_t UpperLimitMCSModule::processBetweenGenAndFit(Int_t /*sampleNum*/) {

  std::cout<<"after generation Test"<<std::endl;

  if (!fitInitParams() || !genSample() || !fitParams() || !fitModel() ) return kFALSE; 

  static_cast<RooRealVar*>(_poi->first())->setVal(static_cast<RooRealVar*>(fitInitParams()->find(_parName.c_str()))->getVal());
  
  //_poi->first()->Print();
  static_cast<RooRealVar*>(_poi->first())->setBins(1000);
  //fitModel()->Print("v");

  std::cout<<"generated Entries:"<<genSample()->numEntries()<<std::endl;

  RooStats::ProfileLikelihoodCalculator plc( *(genSample()), *(fitModel()), *_poi);
   
  //PLC calculates intervals. for one sided ul multiply testsize by two
  plc.SetTestSize(2*(1-_cl));
  RooStats::ConfInterval* pllint=plc.GetInterval();

  if (!pllint) return kFALSE;

  std::cout<<"poi value: "<<((RooRealVar*)( _poi->first()))->getVal()<<std::endl;
  std::cout<<(static_cast<RooRealVar*>((fitParams()->find(_parName.c_str()))))->getVal()<<std::endl;
  std::cout<<((RooStats::LikelihoodInterval*)pllint)->UpperLimit((RooRealVar&)*(_poi->first()))<<std::endl;


  //Go to the fit Value for zour POI to make sure upperlimit works correct.
  //fitModel()->fitTo(*genSample());
  


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