tschum/LimitTool/Model.cc

#include "Model.h"
#include<vector>
#include<string>
#include<sstream>
#include "RooFit.h"
#include "RooAbsPdf.h"
#include "RooProdPdf.h"
#include "RooPoisson.h"
#include "RooGaussian.h"
#include "RooLognormal.h"
#include "RooRealVar.h"
#include "RooConstVar.h"
#include "RooFormulaVar.h"
#include "RooArgSet.h"
#include "RooPolynomial.h"
#include "RooArgList.h"
#include "RooDataSet.h"
#include "RooNLLVar.h"
#include "TMath.h"
#include "RooRandom.h"
#include "RooPlot.h"
#include "RooGamma.h"
#include "RooRealVar.h"
#include "RooAbsReal.h"
#include <iostream>
#include <fstream>
#include <istream>
#include <cstdlib>
#include <cstdarg>
#include <cstring>
#include <cstdio>
#include"TCanvas.h"
#include"TFile.h"
#include"TH1F.h"

#include "RooStats/BayesianCalculator.h"
#include "RooStats/ProfileLikelihoodCalculator.h"
#include "RooStats/LikelihoodIntervalPlot.h"
#include "Model.h"
#include "RooStats/HybridCalculatorOriginal.h"
#include "RooStats/HybridResult.h"
#include "RooStats/HybridPlot.h"
#include "RooStats/ModelConfig.h"
#include "RooStats/FeldmanCousins.h"
#include "RooStats/ToyMCSampler.h"
#include "RooStats/PointSetInterval.h"
#include "RooStats/ConfidenceBelt.h"
#include "RooWorkspace.h"
#include "RooStats/HypoTestInverter.h"
#include "RooStats/HypoTestInverterResult.h"
#include "RooStats/HypoTestPlot.h"

#include "RooStats/MCMCCalculator.h"
#include "RooStats/MCMCInterval.h"
#include "RooStats/MCMCIntervalPlot.h"


using namespace std;
using namespace RooFit;
using namespace RooStats;


string int_to_string(int i){
    std::stringstream out;
    out << i;
    return  out.str();
    
  };

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
Model::Model():n_channels(0),n_backgrounds(0),n_nuisances(0),range_nuisances(-1),range_POI(-1),range_n_obs(-1),use_nuisances(false),corrupted(true)
{
  POI_set=0;
  observable_set=0;
  nuisance_set=0;
  nuisance_prior_pdf=0;
  POI_prior_pdf=0;
  
  POI=0;
  signal=0;
  background=0;
  observables=0;
  nuisance_names=0;
  nuisance_upper_limit=0;
  nuisance_parameters=0;
  nuisance_widths=0;
  scaling_factors=0;
  nuisance_priors=0;
  sb_means=0;
  b_means=0;
  sb_poissons=0;
  b_poissons=0;
  sb_likelihood=0;
  b_likelihood=0;
  data=0;
}

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
Model::Model(ConfigFile* bkgConfig, ConfigFile* sigConfig, int selectChannel) {

  POI_set=0;
  observable_set=0;
  nuisance_set=0;
  nuisance_prior_pdf=0;
  POI_prior_pdf=0;
  
  POI=0;
  signal=0;
  background=0;
  observables=0;
  nuisance_names=0;
  nuisance_upper_limit=0;
  nuisance_parameters=0;
  nuisance_widths=0;
  scaling_factors=0;
  nuisance_priors=0;
  sb_means=0;
  b_means=0;
  sb_poissons=0;
  b_poissons=0;
  sb_likelihood=0;
  b_likelihood=0;
  data=0;
  nuisance_means=0;
  
  //store the input from the background and data file here
  vector<double> *signal_expectation =new vector<double>;
  vector< vector<double>* > *background_expectations =new vector< vector<double>* >;
  vector< vector< vector<double>* >* >* systematic_uncertainties=0;
  vector<int> *n_input=new vector<int>;

  //read [definitions]
  this->n_channels =   bkgConfig->Value("definitions","n_channels");
  this->n_backgrounds =   bkgConfig->Value("definitions","n_backgrounds");
  this->n_nuisances =   bkgConfig->Value("definitions","n_nuisance");
  this->range_nuisances=   bkgConfig->Value("definitions","nuisance_range");
  this->range_POI=   bkgConfig->Value("definitions","parameter_of_interest_range");
  this->range_n_obs=   bkgConfig->Value("definitions","observation_range");
  this->lumiScale= bkgConfig->Value("definitions","lumi_scale",1);

  bool sepStatUnc = bkgConfig->Value("definitions","seperate_stat_unc",false);
  bool sumSigUnc  = bkgConfig->Value("definitions","sum_signal_unc",false);


  if( n_channels      !=   sigConfig->Value("definitions","n_channels") ) {
    cout<<"n_channels don't match for background and signal!!!"<<endl;
    exit(-1);
  }
  if( n_nuisances     <   sigConfig->Value("definitions","n_nuisance") ) {
    cout<<"n_nuisances don't match for background and signal!!!"<<endl;
    exit(-1);
  }
  

  if(selectChannel > 0 ) this->n_channels = 1;
  
  //read [observation] 
  //fill the vector with the observations
  for(int i=1;i<=n_channels;i++){
    string value_id = (string)"n"+int_to_string(i);
    if(selectChannel > 0 ) value_id = (string)"n"+int_to_string(selectChannel);
    int value =   bkgConfig->Value("observation",value_id.c_str());
    //     cout<<value_id.c_str()<<" "<<value<<endl;
    n_input->push_back(int(value*lumiScale+0.5));
  }

  sumS=0;
  //read [signal yield] 
  //fill the vector with the signal expectations
  for(int i=1;i<=n_channels;i++){
    string value_id = (string)"s"+int_to_string(i);
    if(selectChannel > 0 ) value_id = (string)"s"+int_to_string(selectChannel);
    double value =   sigConfig->Value("signal_yield",value_id.c_str());
    //   cout<<value_id.c_str()<<" "<<value<<endl;
    signal_expectation->push_back(value*lumiScale);
    sumS += value*lumiScale;
  }

  sumB=0;
  //read [background_yield_n]
  //fill the vectors with the background expectations
  for(int i=1;i<=n_backgrounds;i++){
    string section_id = (string)"background_yield_" + int_to_string(i);
    vector<double>* background_yield = new vector<double>;      
    for(int j=1;j<=n_channels;j++){
      string value_id = (string)"b" + int_to_string(j);
      if(selectChannel > 0 ) value_id = (string)"b"+int_to_string(selectChannel);
      double value =   bkgConfig->Value(section_id.c_str(),value_id.c_str());
      //      cout<<value_id.c_str()<<" "<<value<<endl;
      background_yield->push_back(value*lumiScale);
      sumB += value*lumiScale;
    }
    background_expectations->push_back(background_yield);
  }
   

  sumBerr=0;
  //read [nuisance_n]
  //fill the object for the systematic uncertainties
  if(this->n_nuisances==0){
    systematic_uncertainties=0;
  }
  else{
    systematic_uncertainties=new vector< vector< vector<double>* >* >;
    this->nuisance_names=new vector<string>;
    this->nuisance_upper_limit=new vector<double>;
    //read [nuisance_n]
    vector<double>* sum_yield = new vector<double>;
    for(int w=0;w<n_channels;w++)
      sum_yield->push_back(0);

    for(int i=1;i<=n_nuisances;i++){  

      string section_id = (string)"nuisance_" + int_to_string(i);
      string name =  bkgConfig->Value(section_id.c_str(),"name");
      string name_compare =  sigConfig->Value(section_id.c_str(),"name","");
      
      if(name_compare== "") {
        cout<<"nuisance name "<<name<<" not found, set to 0 for signal!"<<endl;

      } else if( name != name_compare){
        cout<<"nuisance names don't match for background and signal!!!"<<endl;
        cout<<"background n"<<i<<" has name: "<< name<<endl;
        cout<<"signal n"<<i<<" has name: "<< name_compare<<endl;
        exit(-1);
      }

      if( !sepStatUnc || name.find("BKG_MC_") == string::npos ) {

        vector< vector<double>* >* nuisance_yields=new vector<vector<double>*>;
        vector<double>* nuisance_yield = new vector<double>;
        for(int j=1;j<=n_channels;j++){
          string value_id = (string)"sigma_s_" + int_to_string(j);
          if(selectChannel > 0 ) value_id = (string)"sigma_s_"+int_to_string(selectChannel);
          double value =   sigConfig->Value(section_id.c_str(),value_id.c_str(),0);
          if(name.find("_MC_")==string::npos)
            nuisance_yield->push_back(value);
          else 
            nuisance_yield->push_back(value/sqrt(lumiScale));
          //    cout<<value_id.c_str()<<" "<<value<<endl;
        }
        if(name.find("BKG_")==string::npos&&sumSigUnc) {
 
          for(int j=0;j<n_channels;j++){
            sum_yield->at(j) +=  pow(nuisance_yield->at(j),2);
          }
        } else  {
          this->nuisance_names->push_back(name);
          //upper range for uncertainties -> (X*sigma)
          double upper_limit=  bkgConfig->Value(section_id.c_str(),"upper_limit");
          this->nuisance_upper_limit->push_back(upper_limit);

          nuisance_yields->push_back(nuisance_yield);
          for(int k=1;k<=n_backgrounds;k++){
            vector<double>* nuisance_yield = new vector<double>;
            for(int j=1;j<=n_channels;j++){
              string value_id = (string)"sigma_b" + int_to_string(k) + "_" + int_to_string(j);
              if(selectChannel > 0 ) value_id = (string)"sigma_b" + int_to_string(k) + "_" + int_to_string(selectChannel);
              double value =   bkgConfig->Value(section_id.c_str(),value_id.c_str());
              if(name.find("_MC_")==string::npos) {
                nuisance_yield->push_back(value);
                sumBerr += pow(value,2);
              } else {
                nuisance_yield->push_back(value/sqrt(lumiScale));
                sumBerr += pow(value/sqrt(lumiScale),2);
              }
              //                 cout<<value_id.c_str()<<" "<<value<<endl;
            }
            nuisance_yields->push_back(nuisance_yield);
          }
          systematic_uncertainties->push_back(nuisance_yields);
        }
      } else {
        for(int t=1;t<=n_channels;t++){

          vector< vector<double>* >* nuisance_yields=new vector<vector<double>*>;

          this->nuisance_names->push_back(name+ "_" + int_to_string(t));
          //upper range for uncertainties -> (X*sigma)
          double upper_limit=  bkgConfig->Value(section_id.c_str(),"upper_limit");
          this->nuisance_upper_limit->push_back(upper_limit);

          vector<double>* nuisance_yield = new vector<double>;
          for(int j=1;j<=n_channels;j++){
            nuisance_yield->push_back(0);
          }
          nuisance_yields->push_back(nuisance_yield);
          for(int k=1;k<=n_backgrounds;k++){
            vector<double>* nuisance_yield = new vector<double>;
            for(int j=1;j<=n_channels;j++){
              string value_id = (string)"sigma_b" + int_to_string(k) + "_" + int_to_string(j);
              if(selectChannel > 0 ) value_id = (string)"sigma_b" + int_to_string(k) + "_" + int_to_string(selectChannel);
              double value =   bkgConfig->Value(section_id.c_str(),value_id.c_str());
              if(t==j)  {
                nuisance_yield->push_back(value/sqrt(lumiScale));
                sumBerr += pow(value/sqrt(lumiScale),2);
              } else nuisance_yield->push_back(0);
              //              cout<<value_id.c_str()<<" "<<value<<endl;
            }
            nuisance_yields->push_back(nuisance_yield);
          }
          systematic_uncertainties->push_back(nuisance_yields);
        }
      }
    }
  

    if(sumSigUnc) {

      this->nuisance_names->push_back("SUM_SIG_UNC");
      this->nuisance_upper_limit->push_back(5);
      vector< vector<double>* >* nuisance_yields=new vector<vector<double>*>;
      vector<double>* nuisance_yield = new vector<double>;
      for(int j=0;j<n_channels;j++){
        nuisance_yield->push_back( sqrt(sum_yield->at(j)) );
      }
      nuisance_yields->push_back(nuisance_yield);
      for(int k=1;k<=n_backgrounds;k++){
        vector<double>* nuisance_yield = new vector<double>;
        for(int j=1;j<=n_channels;j++){
          nuisance_yield->push_back(0);
        }
        nuisance_yields->push_back(nuisance_yield);
      }
      systematic_uncertainties->push_back(nuisance_yields);
    }

    //rescale num of nuisance parameters
    n_nuisances = systematic_uncertainties->size();

  }

  check_input(n_input,signal_expectation,background_expectations,systematic_uncertainties);

  if(systematic_uncertainties==0){
    use_nuisances=false;
  }
  else use_nuisances=true;

  create_POI();
  create_observables();
  create_signal_expectations(signal_expectation);
  create_background_expectations(background_expectations);
  create_nuisances(systematic_uncertainties);
  create_scaling_factors(systematic_uncertainties);
  create_POI_prior();

  create_nuisance_priors();
  create_sb_likelihood();
  create_b_likelihood();
  create_dataset(n_input);

  create_pseudo_datasets(signal_expectation,background_expectations,systematic_uncertainties);

  create_teststats();

}
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
Model::Model(std::vector<double>* signal_expectation,std::vector< std::vector<double>* >* background_expectation,std::vector< vector< vector<double>* >* >* systematic_uncertainties,std::vector<double>* nuisance_upper_limits,std::vector<int>* n_input,int nuisance_range,int POI_range,int observation_range):n_channels(signal_expectation->size()),n_backgrounds(background_expectation->size()),n_nuisances(systematic_uncertainties->size()),range_nuisances(nuisance_range),range_POI(POI_range),range_n_obs(observation_range){
  
  POI_set=0;
  observable_set=0;
  nuisance_set=0;
  nuisance_prior_pdf=0;
  POI_prior_pdf=0;
  
  POI=0;
  signal=0;
  background=0;
  observables=0;
  nuisance_names=0;
  nuisance_parameters=0;
  nuisance_widths=0;
  scaling_factors=0;
  nuisance_priors=0;
  sb_means=0;
  b_means=0;
  sb_poissons=0;
  b_poissons=0;
  sb_likelihood=0;
  b_likelihood=0;
  data=0;
  nuisance_means=0;


  nuisance_upper_limit=nuisance_upper_limits;
  
  check_input(n_input,signal_expectation,background_expectation,systematic_uncertainties);
  if(systematic_uncertainties==0){
    use_nuisances=false;
  }
  else use_nuisances=true;
  create_POI();
  create_observables();
  create_signal_expectations(signal_expectation);
  create_background_expectations(background_expectation);
  create_nuisances(systematic_uncertainties);
  create_scaling_factors(systematic_uncertainties);
  create_POI_prior();
  create_nuisance_priors();
  create_sb_likelihood();
  create_b_likelihood();
  if(n_input!=0){
    create_dataset(n_input);
  }


}

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
Model::~Model(){
  if(data!=0){
    delete data;    
  }
  //  cout<<"deleted"<<endl;
  if(b_likelihood!=0){
    delete b_likelihood;    
  }
  //  cout<<"deleted"<<endl;
  if(sb_likelihood!=0){
    delete sb_likelihood;    
  }
  //  cout<<"deleted"<<endl;
  if(b_poissons!=0){
    for(unsigned i=0;i<b_poissons->size();i++){
      delete b_poissons->at(i);   
    }
    delete b_poissons;
  }
  //  cout<<"deleted"<<endl;
  if(sb_poissons!=0){
    for(unsigned i=0;i<sb_poissons->size();i++){
      delete sb_poissons->at(i);   
    }
    delete sb_poissons;
  }
  //  cout<<"deleted"<<endl;
  if(b_means!=0){
    for(unsigned i=0;i<b_means->size();i++){
      delete b_means->at(i);    
    }
    delete b_means;
  }
  //  cout<<"deleted"<<endl;
  if(sb_means!=0){
    for(unsigned i=0;i<sb_means->size();i++){
      delete sb_means->at(i); 
    }
    delete sb_means;
  }
  //  cout<<"deleted"<<endl;
  if(nuisance_priors!=0){
    for(unsigned i=0;i<nuisance_priors->size();i++){
      delete nuisance_priors->at(i);  
    }
    delete nuisance_priors;
  }
  //  cout<<"deleted"<<endl;
  if(scaling_factors!=0){
    for(unsigned i=0;i<scaling_factors->size();i++){
      for(unsigned j=0;j<scaling_factors->at(i)->size();j++){
        for(unsigned k=0;k<scaling_factors->at(i)->at(j)->size();k++){
          delete scaling_factors->at(i)->at(j)->at(k);    
        }
        delete scaling_factors->at(i)->at(j);
      }
      delete scaling_factors->at(i);
    }
    delete scaling_factors;
  }
  //  cout<<"deleted"<<endl;
  if(nuisance_widths!=0){
    for(unsigned i=0;i<nuisance_widths->size();i++){
      delete nuisance_widths->at(i);  
    }
    delete nuisance_widths;
  }
  //  cout<<"deleted"<<endl;
  if(nuisance_parameters!=0){
    for(unsigned i=0;i<nuisance_parameters->size();i++){
      delete nuisance_parameters->at(i);  
    }
    delete nuisance_parameters;
  }
  //  cout<<"deleted"<<endl;
  if(observables!=0){
    for(unsigned i=0;i<observables->size();i++){
      delete observables->at(i);  
    }
    delete observables;
  }
  //  cout<<"deleted"<<endl;
  if(background!=0){
    for(unsigned i=0;i<background->size();i++){
      for(unsigned j=0;j<background->at(i)->size();j++){
        delete background->at(i)->at(j);
      }
      delete background->at(i);
    }
    delete background;
  }
  //  cout<<"deleted"<<endl;
  if(signal!=0){
    for(unsigned i=0;i<signal->size();i++){
      delete signal->at(i);    
    }
    delete signal;
  }
  //  cout<<"deleted"<<endl;
  if(POI!=0){
    delete POI;
  }
  //  cout<<"deleted"<<endl;
  if(POI_prior_pdf!=0){
    delete POI_prior_pdf;
  }
  //  cout<<"deleted"<<endl;
  if(nuisance_prior_pdf!=0){
    delete nuisance_prior_pdf;
  }
  //  cout<<"deleted"<<endl;
  if(nuisance_set!=0){
    delete nuisance_set;
  }
  //  cout<<"deleted"<<endl;
  if(observable_set!=0){
    delete observable_set;
  }
  //  cout<<"deleted"<<endl;
  if(POI_set!=0){
    delete POI_set;
  }
  cout<<"----------------------------"<<endl;
  cout<<endl;
}

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
bool Model::check_input(std::vector<int>* n_input,std::vector<double>* signal_expectation,std::vector< std::vector<double>* >* background_expectation,std::vector< std::vector< std::vector<double>* >* >* systematic_uncertainties){
  
  //  cout<<"Checking, if input is consistent\n"<<endl;
  unsigned n=signal_expectation->size();
  for(unsigned i=0;i<background_expectation->size();i++){
    if(background_expectation->at(i)->size()!=n){
      this->corrupted=true;
      cout<<"ERROR:             size of signal and background expectation "<<i<<" don't match"<<endl;
      return(false);
    }    
  }
  if(n_input!=0){
    if(n_input->size()!=n){
      this->corrupted=true;
      cout<<"ERROR:             size of signal expectation and observation don't match"<<endl;
      return(false);      
    }    
  }
  if(systematic_uncertainties!=0){
    for(unsigned i=0;i<systematic_uncertainties->size();i++){
      if(systematic_uncertainties->at(i)->size()!=background_expectation->size()+1){
        this->corrupted=true;
        cout<<"ERROR:   at "<<i<<"  :    systematic uncertainty matrix doesn't contain values for signal and all backgrounds"<<endl;
        return(false);
      }      
      for(unsigned j=0;j<systematic_uncertainties->at(i)->size();j++){
        if(systematic_uncertainties->at(i)->at(j)->size()!=n){
          this->corrupted=true;
          cout<<"ERROR:         size of signal expectation and systematic uncertainty matrix don't match"<<endl;
          return(false);
        }
      }
    }
    if(nuisance_upper_limit->size() != systematic_uncertainties->size() ) {
      this->corrupted=true;
      cout<<"ERROR:             size of nuisance upper limits and systematic uncertainties don't match"<<endl;
      return(false);
    }
  }


  cout<<"Checked, that input is consistent!\n"<<endl;

  this->corrupted=false;
  return(true);
}
// this overwrites existing dataset
void Model::set_dataset(RooDataSet* d){
  delete data;
  data=d;

}

///set nuisances to zero andconstant
void Model::set_nuisance_const(){
  
 RooArgSet* nuisances=nuisance_set;
 TIterator* it=nuisances->createIterator();
 //RooAbsArg* nui;
 // while((RooAbsReal*)it->Next()){
 RooRealVar* nui;
 while( (nui=(RooRealVar*) it->Next()) ){
   //(RooAbsReal*)it->setVal(0);
   RooRealVar* nui2=dynamic_cast<RooRealVar*>(nui);
   //cout<<"nui "<<nui2<<endl;
   nui2->setVal(0);
   nui2->setConstant(kTRUE);
   
   
 }

 
}
///set nuisances to zero andconstant
void Model::set_nuisance_var(){
  
 RooArgSet* nuisances=nuisance_set;
 TIterator* it=nuisances->createIterator();
 RooRealVar* nui;
 // while((RooAbsReal*)it->Next()){
 while( (nui=(RooRealVar*)it->Next()) ){
   //(RooAbsReal*)it->setVal(0);
   RooRealVar* nui2=dynamic_cast<RooRealVar*>(nui);
   //cout<<"nui "<<nui2<<endl;
   nui2->setVal(0);
   nui2->setConstant(kFALSE);
   
   
 }

}

void Model::set_poi_const(double p=0){
  RooRealVar* ss=get_POI();
  ss->setVal(p);
  ss->setConstant(kTRUE);


}
void Model::set_poi_var(double p=0){
  RooRealVar* ss=get_POI();  
  ss->setVal(p);
  ss->setConstant(kFALSE);

}


////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
void Model::create_observables(){
  // cout<<"Create observables\n"<<endl; 
 if(range_n_obs>0){
   observables=new vector<RooRealVar*>;
   observable_set=new RooArgSet();
   char* name=new char[1000];
   
   for(int i=0;i<this->n_channels;i++){
     sprintf (name, "observation_%d",i);
     RooRealVar* tmp_observables=new RooRealVar(name,name,0,range_n_obs);
     observables->push_back(tmp_observables); 
     observable_set->add(*observables->at(i));
   }      
   delete [] name;
 }
 else{
   cout<<"ERROR:        negative range for the observations\n"<<endl;
   this->corrupted=true;
 }
}

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
void Model::create_signal_expectations(std::vector<double>* signal_expectation){
  //  cout<<"Create signal expectations\n"<<endl;
  signal=new std::vector< RooConstVar*>;
  char* name=new char[1000];
  
  for(int i=0;i<this->n_channels;i++){
    sprintf (name, "signal_%d", i);
    RooConstVar* tmp_signal=new RooConstVar(name,name,signal_expectation->at(i));
    signal->push_back(tmp_signal);
  }
  delete [] name;
}

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
void Model::create_background_expectations(std::vector< std::vector<double>* >* background_expectation){
  //  cout<<"Create background expectations\n"<<endl;
  background =new vector< vector<RooConstVar*>* >;
  char* name=new char[1000];
  
  for(int i=0;i<this->n_backgrounds;i++){
    vector<RooConstVar*>* tmp_background_row=new vector<RooConstVar*>;
    for(int j=0;j<this->n_channels;j++){
      sprintf (name, "background_channel_%d_source_%d",j,i);
      RooConstVar* tmp_background=new RooConstVar(name,name,background_expectation->at(i)->at(j));
      tmp_background_row->push_back(tmp_background);
    }
    background->push_back(tmp_background_row);
  }
  delete [] name;
}

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
void Model::create_nuisances(std::vector< vector< vector<double>* >* >* systematic_uncertainties){
  //   cout<<"Create nuisance parameters\n"<<endl;
   if(range_nuisances>0){
     if(use_nuisances==true){   
       
       nuisance_set=new RooArgSet();
       nuisance_parameters=new vector<RooRealVar*>;
       nuisance_widths=new vector<RooConstVar*>;
       char* name=new char[1000];      
       char* title=new char[1000];
       
       for(unsigned i=0;i<systematic_uncertainties->size();i++){
         double greatest_error=0;
         for(unsigned j=0;j<systematic_uncertainties->at(i)->size();j++){
           for(unsigned k=0;k<systematic_uncertainties->at(i)->at(j)->size();k++){
             if(systematic_uncertainties->at(i)->at(j)->at(k) >greatest_error){
               greatest_error=systematic_uncertainties->at(i)->at(j)->at(k);      
             }
           }
         }
         if(nuisance_names==0){
           sprintf (name, "sigma_%d", i);
           sprintf (title, "sigma_%d", i);
         }
         else{
           sprintf (name, "sigma_%d", i);
           sprintf (title, "sigma_%s", this->nuisance_names->at(i).c_str());
         }
         RooConstVar* tmp_nuisance_width=new RooConstVar(name,title,greatest_error);
         nuisance_widths->push_back(tmp_nuisance_width);
         if(nuisance_names==0){
           sprintf (name, "delta_%d", i);
           sprintf (title, "delta_%d", i);
         }
         else{
           sprintf (name, "delta_%d", i);
           sprintf (title, "delta_%s", this->nuisance_names->at(i).c_str());
         }
         //RooRealVar* tmp_nuisance_parameters=new RooRealVar(name,title,0,-1*range_nuisances*nuisance_widths->at(i)->getVal(),1.0*range_nuisances*nuisance_widths->at(i)->getVal());
         //RooRealVar* tmp_nuisance_parameters=new RooRealVar(name,title,0,-1,1.0*range_nuisances*nuisance_widths->at(i)->getVal());
         
         double lower_limit=-1;
          std::string line = title;
         std::string::size_type pos=line.find("[LOGNORMAL]");
         if(pos != std::string::npos) lower_limit=-0.9999;
         RooRealVar* tmp_nuisance_parameters=new RooRealVar(name,title,0,lower_limit,1.0*nuisance_upper_limit->at(i)*nuisance_widths->at(i)->getVal());
         //      cout<<name<<" has upper limit "<<nuisance_upper_limit->at(i)<<endl;

         nuisance_parameters->push_back(tmp_nuisance_parameters);
         nuisance_set->add(*nuisance_parameters->at(i));              
       }
       delete [] name;
       delete [] title;
     }
     else cout<<"No systematic uncertainties specified, no nuisance parameters created\n"<<endl;
   }
   else{
     cout<<"ERROR:      negative range for the nuisance_parameters\n"<<endl;
     this->corrupted=true;
   }
}

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
void Model::create_scaling_factors(vector< vector< vector<double>* >* >* systematic_uncertainties){
  //  cout<<"Create scaling factors for nuisances\n"<<endl;
  if(use_nuisances==true &&nuisance_parameters!=0){     
    
    scaling_factors = new vector< vector< vector<RooConstVar*>* >* >;
    char* name=new char[1000]; 
  
    for(unsigned i=0;i<nuisance_parameters->size();i++){
      
      vector< vector<RooConstVar*>* >*  tmp_scaling_factors_for_one_nuisance = new vector< vector<RooConstVar*>* >;
      vector<RooConstVar*>* tmp_signal_scaling_for_one_nuisance = new vector<RooConstVar*>;

      for(int j=0;j<this->n_channels;j++){
        sprintf (name, "signal_scaling_factor_channel_%d_nuisance_%d", j,i);
        RooConstVar* tmp_signal_scaling=new RooConstVar(name,name,1.0*systematic_uncertainties->at(i)->at(0)->at(j)/nuisance_widths->at(i)->getVal());
        tmp_signal_scaling_for_one_nuisance->push_back(tmp_signal_scaling);
      }
      tmp_scaling_factors_for_one_nuisance->push_back(tmp_signal_scaling_for_one_nuisance);

      for(unsigned j=0;j<background->size();j++){
        
        vector<RooConstVar*>* tmp_background_scaling_for_one_nuisance_for_one_background_source = new vector<RooConstVar*>;
        
        for(int k=0;k<this->n_channels;k++){
          sprintf (name, "background_scaling_factor_channel_%d_source_%d_nuisance_%d",k,j,i);
          RooConstVar* tmp_background_scaling=new RooConstVar(name,name,1.0*systematic_uncertainties->at(i)->at(j+1)->at(k)/nuisance_widths->at(i)->getVal());
          tmp_background_scaling_for_one_nuisance_for_one_background_source->push_back(tmp_background_scaling);
        }
        tmp_scaling_factors_for_one_nuisance->push_back(tmp_background_scaling_for_one_nuisance_for_one_background_source);
      }
      scaling_factors->push_back(tmp_scaling_factors_for_one_nuisance);
    }
    delete [] name;
  }
  else{
    if(use_nuisances==true){
      cout<<"ERROR: systematics specified, but no nuisance parameters created yet; can't define scaling factors\n"<<endl;
      this->corrupted=true;
    }
    else cout<<"No systematics_specified, no scaling factors needed\n"<<endl;
  }
}


////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// Use this for Gaussian errors
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
void Model::create_nuisance_priors(){
  //void Model::create_nuisance_priors(){
  cout<<"Create nuisances priors\n"<<endl;  
  if(use_nuisances==true &&nuisance_parameters!=0){
    
    nuisance_priors=new vector<RooAbsPdf*>;
    nuisance_means=new vector<RooFormulaVar*>;
    RooArgList* nuisance_prior_list=new RooArgList();
    char* name=new char[1000]; 
    
    for(unsigned i=0;i<nuisance_parameters->size();i++){
      //here check if in name is gauss or logNormal
      //cout<<(*nuiscance_parameters)[i].GetTitle()<<endl;
      RooRealVar* dummy = nuisance_parameters->at(i);
      //cout<<dummy->GetTitle()<<endl;
      
      int distribution=0;
      std::string line = dummy->GetTitle();///now find in string gauss or logNormal
      std::string::size_type pos=line.find("[GAUSS]");
      if(pos != std::string::npos) distribution=0;
      pos=line.find("[LOGNORMAL]");
      if(pos != std::string::npos) distribution=1;
      
      if(distribution==0){
        ///using gaussian function as nuisance distribution
        cout<<" ---> "<<line<<" using Gaussian Distribuion"<<endl;
        sprintf (name, "nuisance_prior_%d", i);
        RooGaussian* tmp_nuisance_priors=new RooGaussian(name,name,*nuisance_parameters->at(i),RooFit::RooConst(0),*nuisance_widths->at(i));
        nuisance_priors->push_back(tmp_nuisance_priors);
        nuisance_prior_list->add(*nuisance_priors->at(i));
      }
      else if(distribution==1){ 
        ///using logNormal function as nuiscance distribution
        cout<<" ---> "<<line<<" using logNormal Distribuion"<<endl;
        char* name2=new char[1000]; 

        sprintf (name, "nuisance_mean_%d", i);
        sprintf(name2,"%s+1",nuisance_parameters->at(i)->GetName());
        RooFormulaVar* tmp_nuisance_means=new RooFormulaVar(name,name2,RooArgSet(*nuisance_parameters->at(i)));
        nuisance_means->push_back(tmp_nuisance_means);
        sprintf (name, "nuisance_prior_%d", i);
        //RooLognormal* tmp_nuisance_priors=new RooLognormal(name, name,*nuisance_means->at(i), RooFit::RooConst(1), RooFit::RooConst(TMath::Exp(nuisance_widths->at(i)->getVal())));
        //ok i have to adjust my mean parameter or each nusiance....
        ///by my handmade calculations this is mu=sigma^2
        double new_mean=nuisance_widths->at(i)->getVal()*nuisance_widths->at(i)->getVal();
        //and put it into the function as exponential...
        RooLognormal* tmp_nuisance_priors=new RooLognormal(name, name,*nuisance_means->at(i), RooFit::RooConst(TMath::Exp(new_mean)), RooFit::RooConst(TMath::Exp(nuisance_widths->at(i)->getVal())));
        //cout<<"----------------------------"<<name<<" "<<nuisance_widths->at(i)->getVal()<<" "<<TMath::Exp(nuisance_widths->at(i)->getVal())<<endl;
        nuisance_priors->push_back(tmp_nuisance_priors);
        nuisance_prior_list->add(*nuisance_priors->at(i));


      }
    }
    nuisance_prior_pdf=new RooProdPdf("nuisance_prior_pdf","nuisance_prior_pdf",*nuisance_prior_list);
    delete nuisance_prior_list;
    delete [] name;
  }
  else{
    if(use_nuisances==true){
      cout<<"ERROR: systematics specified, but no nuisance parameters created yet; can't create nuisance pdfs\n"<<endl;
      this->corrupted=true;
    }
    else cout<<"No systematics_specified, no nuisance_priors needed\n"<<endl;
  }

  cout<<endl;

}

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// Use this for Gaussian errors
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// void Model::create_nuisance_priors(){
//   //void Model::create_nuisance_priors(){
//   cout<<"Create nuisances priors\n"<<endl;  
//   if(use_nuisances==true &&nuisance_parameters!=0){
    
//     nuisance_priors=new vector<RooAbsPdf*>;
//     RooArgList* nuisance_prior_list=new RooArgList();
//     char* name=new char[1000]; 
    
//     for(int i=0;i<nuisance_parameters->size();i++){
//       sprintf (name, "nuisance_prior_%d", i);
//       RooGaussian* tmp_nuisance_priors=new RooGaussian(name,name,*nuisance_parameters->at(i),RooFit::RooConst(0),*nuisance_widths->at(i));
//       nuisance_priors->push_back(tmp_nuisance_priors);
//       nuisance_prior_list->add(*nuisance_priors->at(i));
//     }
//     nuisance_prior_pdf=new RooProdPdf("nuisance_prior_pdf","nuisance_prior_pdf",*nuisance_prior_list);
//     delete nuisance_prior_list;
//     delete [] name;
//   }
//   else{
//     if(use_nuisances==true){
//       cout<<"ERROR: systematics specified, but no nuisance parameters created yet; can't create nuisance pdfs\n"<<endl;
//       this->corrupted=true;
//     }
//     else cout<<"No systematics_specified, no nuisance_priors needed\n"<<endl;
//   }
// }

// ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// // Use this for LogNormal errors
// ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
void Model::create_nuisance_priors_2(){
     if(use_nuisances==true &&nuisance_parameters!=0){
     
     nuisance_priors=new vector<RooAbsPdf*>;
     nuisance_means=new vector<RooFormulaVar*>;
     RooArgList* nuisance_prior_list=new RooArgList();
     char* name=new char[1000]; 
     
     for(unsigned i=0;i<nuisance_parameters->size();i++){
       char* name2=new char[1000]; 
       sprintf (name, "nuisance_mean_%d", i);
       sprintf(name2,"%s+1",nuisance_parameters->at(i)->GetName());
       RooFormulaVar* tmp_nuisance_means=new RooFormulaVar(name,name2,RooArgSet(*nuisance_parameters->at(i)));
       nuisance_means->push_back(tmp_nuisance_means);
       sprintf (name, "nuisance_prior_%d", i);
       RooLognormal* tmp_nuisance_priors=new RooLognormal(name, name,*nuisance_means->at(i), RooFit::RooConst(1), RooFit::RooConst(TMath::Exp(nuisance_widths->at(i)->getVal())));

       nuisance_priors->push_back(tmp_nuisance_priors);
       nuisance_prior_list->add(*nuisance_priors->at(i));
       delete [] name2;
     }
     nuisance_prior_pdf=new RooProdPdf("nuisance_prior_pdf","nuisance_prior_pdf",*nuisance_prior_list);
     delete nuisance_prior_list;
     delete [] name;
   }
   else{
     if(use_nuisances==true){
       cout<<"ERROR: systematics specified, but no nuisance parameters created yet; can't create nuisance pdfs\n"<<endl;
       this->corrupted=true;
     }
     else cout<<"No systematics_specified, no nuisance_priors needed\n"<<endl;
   }
 }

// ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// // Use this for Gamma errors, does not work
// ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
// void Model::create_nuisance_priors(){
//   
//   if(use_nuisances==true &&nuisance_parameters!=0){
//     
//     nuisance_priors=new vector<RooAbsPdf*>;
//     nuisance_means=new vector<RooFormulaVar*>;
//     RooArgList* nuisance_prior_list=new RooArgList();
//     char* name=new char[1000]; 
//     RooConstVar gamma("gamma","gamma",20);
//     
//     for(int i=0;i<nuisance_parameters->size();i++){
//       char* name2=new char[1000]; 
//       sprintf (name, "beta_%d", i);
//       RooConstVar* beta=new RooConstVar(name,name,TMath::Sqrt(nuisance_widths->at(i)->getVal()*nuisance_widths->at(i)->getVal()/gamma.getVal()));
//       sprintf(name2,"%s+beta_%d*gamma",nuisance_parameters->at(i)->GetName(),i);
//       nuisance_parameters->at(i)->setRange(TMath::Max(-1.0,-beta->getVal()*gamma.getVal()),range_nuisances*nuisance_widths->at(i)->getVal());
//       sprintf (name, "nuisance_mean_%d", i);
//       RooFormulaVar* tmp_nuisance_means=new RooFormulaVar(name,name2,RooArgSet(*beta,gamma,*nuisance_parameters->at(i)));
//       nuisance_means->push_back(tmp_nuisance_means);
//       sprintf (name, "nuisance_prior_%d", i);
//       RooGamma* tmp_nuisance_priors=new RooGamma(name, name,*nuisance_means->at(i),RooFit::RooConst(gamma.getVal()),RooFit::RooConst(beta->getVal()),RooFit::RooConst(0));
//       tmp_nuisance_priors->Print();
//       nuisance_priors->push_back(tmp_nuisance_priors);
//       nuisance_prior_list->add(*nuisance_priors->at(i));
//       delete [] name2;
//       delete beta;
//     }
//     nuisance_prior_pdf=new RooProdPdf("nuisance_prior_pdf","nuisance_prior_pdf",*nuisance_prior_list);
//     nuisance_prior_pdf->Print();
//     delete nuisance_prior_list;
//     delete [] name;
//   }
//   else{
//     if(use_nuisances==true){
//       cout<<"ERROR: systematics specified, but no nuisance parameters created yet; can't create nuisance pdfs\n"<<endl;
//       this->corrupted=true;
//     }
//     else cout<<"No systematics_specified, no nuisance_priors needed\n"<<endl;
//   }
// }


////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
void Model::create_POI(){
  if(range_POI>0){
    //    cout<<"Create parameter of interest\n"<<endl;
    //POI=new RooRealVar("POI","POI",1.0*range_POI/2,0,range_POI);
    POI=new RooRealVar("POI","POI",1.0*range_POI/2,0,range_POI);
    POI_set=new RooArgSet(*POI);
  }
  else{
    cout<<"ERROR:       negative range specified, problem with input\n"<<endl;
    this->corrupted=true;
  }
}

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
void Model::create_POI_prior(){
  if(POI!=0){
    //    cout<<"Create flat prior for the parameter of interest\n"<<endl;
    RooAbsPdf* tmp_POI_prior= new RooPolynomial("tmp_POI_pior","tmp_POI_pior",*POI,RooArgList(RooFit::RooConst(0)));  
    //RooJeffreysPrior* tmp_POI_prior =new RooJeffreysPrior("jeffreys","jeffreys",*w.pdf("p"),*POI,*observable_set);
    //RooJeffreysPrior pi("jeffreys","jeffreys",*w.pdf("p"),*w.set("poi"),*w.set("obs"));
  
//     tmp->Print();
    POI_prior_pdf=new RooProdPdf("POI_prior","POI_prior",RooArgList(*tmp_POI_prior));
    POI_prior_pdf->getVal();
    delete tmp_POI_prior;
  }
  else{
    cout<<"ERROR:       no parameter of interest specified yet, can't create prior\n"<<endl;
    this->corrupted=true;
  }
}

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
void Model::create_sb_likelihood(){
  if(corrupted==false){
    cout<<"Create sb_likelihood\n"<<endl;
    sb_means = new vector<RooFormulaVar*>;
    sb_poissons = new vector<RooPoisson*>;
    vector<RooArgSet*>* sb_params =new vector<RooArgSet*>;
    char* name=new char[1000]; 
  
    for(int i=0;i<n_channels;i++){
      RooArgSet* tmp_sb_params=new RooArgSet(*POI,*signal->at(i));
      for(int j=0;j<n_backgrounds;j++){
        tmp_sb_params->add(*background->at(j)->at(i));
        if(use_nuisances==true){
          for(unsigned k=0;k<nuisance_parameters->size();k++){
            tmp_sb_params->add(*nuisance_parameters->at(k));
            tmp_sb_params->add(*scaling_factors->at(k)->at(0)->at(i));
            tmp_sb_params->add(*scaling_factors->at(k)->at(j+1)->at(i));
          }
        }
      }
      sb_params->push_back(tmp_sb_params);
    }
    //Poisson means and Poisson distributions
    stringstream* formula_string=new stringstream();  
    for(int i=0;i<n_channels;i++){
      formula_string->str("");
      *formula_string<<"1*(";
      *formula_string<<"POI";
      *formula_string<<"*signal_"<<i;
      if(use_nuisances==true){
        for(unsigned j=0;j<nuisance_parameters->size();j++){
          *formula_string<<"*(1+delta_"<<j<<"*signal_scaling_factor_channel_"<<i<<"_nuisance_"<<j<<")";
        }
      }
      for(int k=0;k<n_backgrounds;k++){
        *formula_string<<"+background_channel_"<<i<<"_source_"<<k;
        if(use_nuisances==true){        
          for(unsigned j=0;j<nuisance_parameters->size();j++){
            *formula_string<<"*(1+delta_"<<j<<"*background_scaling_factor_channel_"<<i<<"_source_"<<k<<"_nuisance_"<<j<<")";
          }
        }
      }
      *formula_string<<")";
      cout<<formula_string->str()<<endl;
      sprintf (name, "sb_mean_%d",i);
      RooFormulaVar* tmp_sb_means=new RooFormulaVar(name,formula_string->str().c_str(),*sb_params->at(i));
      sb_means->push_back(tmp_sb_means);
      sprintf (name, "sb_poisson_%d",i);
      RooPoisson* tmp_sb_poissons=new RooPoisson(name,name,*observables->at(i),*sb_means->at(i));
      sb_poissons->push_back(tmp_sb_poissons);
    }
  
    RooArgList *sb_pdfs=new RooArgList();
    for(int i=0;i<n_channels;i++)  {
      sb_pdfs->add(*sb_poissons->at(i));
    }
    sb_likelihood=new RooProdPdf("sb_likelihood","sb_likelihood",*sb_pdfs);
    
    delete sb_pdfs;
    for(unsigned i=0;i<sb_params->size();i++){
      delete sb_params->at(i);
    }
    delete sb_params;
    delete formula_string;
    delete [] name;
  }
}

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
void Model::create_b_likelihood(){
  
  if(corrupted==false){
    cout<<"Create b_likelihood\n"<<endl;
    b_means = new vector<RooFormulaVar*>;
    b_poissons = new vector<RooPoisson*>;
    vector<RooArgSet*>* b_params =new vector<RooArgSet*>;
    char* name=new char[1000];   
    for(int i=0;i<n_channels;i++){
      RooArgSet* tmp_b_params=new RooArgSet();
      for(int j=0;j<n_backgrounds;j++){
        tmp_b_params->add(*background->at(j)->at(i));
        if(use_nuisances==true){
          for(unsigned k=0;k<nuisance_parameters->size();k++){
            tmp_b_params->add(*nuisance_parameters->at(k));
            tmp_b_params->add(*scaling_factors->at(k)->at(j+1)->at(i));
          }
        }
      }
      b_params->push_back(tmp_b_params);
    }
    //Poisson means and Poisson distributions
    stringstream* formula_string=new stringstream();  
    for(int i=0;i<n_channels;i++){
      formula_string->str("");
      *formula_string<<"1*(";
      for(int k=0;k<n_backgrounds;k++){
        if(k==0){
          *formula_string<<"background_channel_"<<i<<"_source_"<<k;
        }
        else{
          *formula_string<<"+background_channel_"<<i<<"_source_"<<k;
        }
        if(use_nuisances==true){
          for(unsigned j=0;j<nuisance_parameters->size();j++){
            *formula_string<<"*(1+delta_"<<j<<"*background_scaling_factor_channel_"<<i<<"_source_"<<k<<"_nuisance_"<<j<<")";
          }
        }      
      }
      *formula_string<<")";
      sprintf (name, "b_mean_%d",i);
      cout<<formula_string->str()<<endl;
      //     b_params->at(i)->Print();
      RooFormulaVar* tmp_b_means=new RooFormulaVar(name,formula_string->str().c_str(),*b_params->at(i));
      b_means->push_back(tmp_b_means);
      sprintf (name, "b_poisson_%d",i);
      RooPoisson* tmp_b_poissons=new RooPoisson(name,name,*observables->at(i),*b_means->at(i));
      b_poissons->push_back(tmp_b_poissons);
    }
    RooArgList *b_pdfs=new RooArgList();
    for(int i=0;i<n_channels;i++)  {
      b_pdfs->add(*b_poissons->at(i));
    }
    b_likelihood=new RooProdPdf("b_likelihood","b_likelihood",*b_pdfs);
    delete b_pdfs;
    for(unsigned i=0;i<b_params->size();i++){
      delete b_params->at(i);
    }
    delete b_params;
    delete formula_string;
    delete [] name;
  }
}

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
bool Model::create_dataset(std::vector<int>* n_input){
  if(corrupted==false){  
    //    cout<<"Create data from input\n"<<endl;
    if(int(n_input->size())!=n_channels){
      cout<<"ERROR:     observation input doesn't match number of channels\n"<<endl;
      this->corrupted=true;
      return(false);
    }
    for(int i=0;i<n_channels;i++){
      observables->at(i)->setVal(n_input->at(i));  
    }
    if(data!=0){
      delete data;
    }
    data=new RooDataSet("data","data",*observable_set);
    data->add(*observable_set);
    return (true);
  }

  return(false);
}

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

bool Model::create_pseudo_datasets(std::vector<double>* signal_expectation, std::vector< std::vector<double>* >* background_expectation, std::vector< std::vector< std::vector<double>* >* >* systematic_uncertainties){

  char* name=new char[1000]; 

  pseudo_b_obs_set=new RooArgSet();
  pseudo_bDN_obs_set=new RooArgSet();
  pseudo_bUP_obs_set=new RooArgSet();
  pseudo_sb_obs_set=new RooArgSet();

  for(int i=0;i<this->n_channels;i++){

    double b = 0;
    double bDN = 0;
    double bUP = 0;
    double sb = signal_expectation->at(i);

    for(int j=0;j<n_backgrounds;j++){

      b += background_expectation->at(j)->at(i);
      sb += background_expectation->at(j)->at(i);

      double uncsum2 = 0;
      for(unsigned k=0;k<nuisance_parameters->size();k++){

        
        uncsum2 += pow(systematic_uncertainties->at(k)->at(j+1)->at(i),2);
      }
      
      bDN += background_expectation->at(j)->at(i) * (1-sqrt(uncsum2));
      bUP += background_expectation->at(j)->at(i) * (1+sqrt(uncsum2));

    }
    
    sprintf (name, "observation_%d",i);

    RooRealVar* tmp_b=new RooRealVar(name,name,int(b+0.5),0,range_n_obs);
    tmp_b->Print();
    RooRealVar* tmp_bDN=new RooRealVar(name,name,int(bDN+0.5),0,range_n_obs);
    tmp_bDN->Print();
    RooRealVar* tmp_bUP=new RooRealVar(name,name,int(bUP+0.5),0,range_n_obs);
    tmp_bUP->Print();
    RooRealVar* tmp_sb=new RooRealVar(name,name,int(sb+0.5),0,range_n_obs);
    tmp_sb->Print();

 
    pseudo_b_obs_set->add(*tmp_b);
    pseudo_bDN_obs_set->add(*tmp_bDN);
    pseudo_bUP_obs_set->add(*tmp_bUP);
    pseudo_sb_obs_set->add(*tmp_sb);


  }      
  delete [] name;

  pseudo_data_b=new RooDataSet("pseudo_data_b","pseudo_data_b",*pseudo_b_obs_set);
  pseudo_data_b->add(*pseudo_b_obs_set);

  pseudo_data_b_minusSigma=new RooDataSet("pseudo_data_b_minusSigma","pseudo_data_b_minusSigma",*pseudo_bDN_obs_set);
  pseudo_data_b_minusSigma->add(*pseudo_bDN_obs_set);

  pseudo_data_b_plusSigma=new RooDataSet("pseudo_data_b_plusSigma","pseudo_data_b_plusSigma",*pseudo_bUP_obs_set);
  pseudo_data_b_plusSigma->add(*pseudo_bUP_obs_set);

  pseudo_data_sb=new RooDataSet("pseudo_data_sb","pseudo_data_sb",*pseudo_sb_obs_set);
  pseudo_data_sb->add(*pseudo_sb_obs_set);


  return true;
}
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
RooRealVar* Model::get_POI(){
 return(this->POI); 
}

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
RooArgSet* Model::get_POI_set(){
 return(this->POI_set); 
}

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
RooArgSet* Model::get_observable_set(){
 return(this->observable_set); 
}

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
RooArgSet* Model::get_nuisance_set(){
 return(this->nuisance_set); 
}
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
RooProdPdf* Model::get_POI_prior(){
  return(POI_prior_pdf);
}

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
RooProdPdf* Model::get_nuisance_prior_pdf(){
  return(nuisance_prior_pdf);
}

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
RooProdPdf* Model::get_sb_likelihood(){
  return(sb_likelihood);
}

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
RooProdPdf* Model::get_b_likelihood(){
  return(this->b_likelihood);
}

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
RooProdPdf* Model::get_complete_likelihood(){
  RooProdPdf* complete_likelihood=new RooProdPdf("complete_likelihood","complete_likelihood",RooArgList(*this->sb_likelihood,*this->nuisance_prior_pdf));
  return(complete_likelihood);
}

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
RooDataSet* Model::get_data(){
  return(this->data);
}
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
RooDataSet* Model::get_pseudo_data_b(){
  return(this->pseudo_data_b);
}
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
RooDataSet* Model::get_pseudo_data_bDN(){
  return(this->pseudo_data_b_minusSigma);
}
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
RooDataSet* Model::get_pseudo_data_bUP(){
  return(this->pseudo_data_b_plusSigma);
}
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
RooDataSet* Model::get_pseudo_data_sb(){
  return(this->pseudo_data_sb);
}
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
void Model::set_random_seed(int i){
 cout<<"Setting random seed to "<<i<<endl;
 RooRandom::randomGenerator()->SetSeed(i);
}
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
double Model::get_sum_signal(){
  return(this->sumS);
}
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
double Model::get_sum_background(){
  return(this->sumB);
}
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
double Model::get_total_bunc(){
  return(sqrt(this->sumBerr)*this->sumB);
}
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
double  Model::get_lumi_scale(){
  return(this->lumiScale);
}
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
void Model::Print(){
   ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
  //Some output of the model parameters
  ////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
  cout<<"---------------------------------------------------------------------------------------------------------------------"<<endl;
  cout<<"Model information"<<endl;
  cout<<"---------------------------------------------------------------------------------------------------------------------\n\n\n"<<endl;
  cout<<"Number of channels:    "<<this->n_channels<<"\nNumber of background sources per channel:       "<<this->n_backgrounds<<"\n\n\n"<<endl;
  cout<<"Channel                signal expectation                      background_expectation(s)\n"<<endl;
  for(int i=0;i<this->n_channels;i++){
   cout<<i<<"                   "<<this->signal->at(i)->getVal();
   for(int j=0;j<this->n_backgrounds;j++){
     cout<<"                            "<<this->background->at(j)->at(i)->getVal(); 
   }
   cout<<"\n";
  }
  
  if(use_nuisances==true){
    cout<<"\n\n---------------------------------------------------------------------------------------------------------------------"<<endl;
    cout<<"Systematic uncertainties"<<endl;
    cout<<"---------------------------------------------------------------------------------------------------------------------\n\n"<<endl;
    for(unsigned i=0;i<this->nuisance_parameters->size();i++){
      if(nuisance_names)    cout<<"Matrix for "<<nuisance_names->at(i).c_str()<<" uncertainty "<<"\n"<<endl; 
      else    cout<<"Matrix for uncertainty "<<i<<"\n"<<endl; 
      cout<<"                   sigma_signal            sigma_b\n"<<endl;
      for(int j=0;j<this->n_channels;j++){
        cout<<"Channel "<<j<<"          ";
        for(int k=0;k<this->n_backgrounds+1;k++){
          cout<<this->nuisance_widths->at(i)->getVal()*this->scaling_factors->at(i)->at(k)->at(j)->getVal()<<"                  ";
        }
        cout<<"\n";
      }
      cout<<"\n\n";
    }
    cout<<"\n\n---------------------------------------------------------------------------------------------------------------------"<<endl;
    cout<<"scaling factors"<<endl;
    cout<<"---------------------------------------------------------------------------------------------------------------------\n\n\n"<<endl;
    for(unsigned i=0;i<nuisance_parameters->size();i++){
      if(nuisance_names) cout<<"Scaling factors for "<<nuisance_names->at(i).c_str()<<" nuisance parameter "<<"\n"<<endl;
      else               cout<<"Scaling factors for nuisance parameter "<<i<<"\n"<<endl;
      cout<<"Width of nuisance distribution: "<<nuisance_widths->at(i)->getVal()<<"\n"<<endl;
      cout<<"           signal scaling          background scaling\n"<<endl;
      for(int j=0;j<n_channels;j++){
        cout<<"channel "<<j<<"  "<<scaling_factors->at(i)->at(0)->at(j)->getVal();
        for(unsigned k=0;k<background->size();k++){
          cout<<"               "<<scaling_factors->at(i)->at(k+1)->at(j)->getVal();
        }
        cout<<"\n"<<endl;
      }
      cout<<"\n\n\n"<<endl;
    }
  }
}


void Model::test(string file){
  ofstream myfile;
  myfile.open (file.c_str(), fstream::app);
  myfile<<endl;
  myfile<<"[observation]"<<endl;
  char* name=new char[1000]; 
  char* hmm=new char[1000];
        for(int i=0;i<this->n_channels;i++){
          sprintf (name, "observation_%d",i);
          //cout<<this->data->get()->getRealValue(name)<<endl;
          sprintf (hmm, "n%d",i+1);
          myfile<<hmm<<" = "<<this->data->get()->getRealValue(name)<<endl;
        }
 myfile.close();
 delete [] name;
}
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////

void Model::create_teststats() {


  RooRealVar* firstPOI = (RooRealVar*) this->get_POI_set()->first();
  firstPOI->setVal(1);
 
  RooNLLVar* sb_nll  = new RooNLLVar("sb_nll","sb_nll",*this->get_sb_likelihood(),*this->get_data());
  RooNLLVar* b_nll   = new RooNLLVar("b_nll","b_nll",*this->get_b_likelihood(),*this->get_data());
  test_Data = 2*(sb_nll->getVal()-b_nll->getVal());

  sb_nll  = new RooNLLVar ("sb_nll","sb_nll",*this->get_sb_likelihood(),*this->get_pseudo_data_b());
  b_nll   = new RooNLLVar ("b_nll","b_nll",*this->get_b_likelihood(),*this->get_pseudo_data_b());
  test_Bmean = 2*(sb_nll->getVal()-b_nll->getVal());

  sb_nll  = new RooNLLVar ("sb_nll","sb_nll",*this->get_sb_likelihood(),*this->get_pseudo_data_bDN());
  b_nll   = new RooNLLVar ("b_nll","b_nll",*this->get_b_likelihood(),*this->get_pseudo_data_bDN());
  test_Bdown = 2*(sb_nll->getVal()-b_nll->getVal());

  sb_nll  = new RooNLLVar ("sb_nll","sb_nll",*this->get_sb_likelihood(),*this->get_pseudo_data_bUP());
  b_nll   = new RooNLLVar ("b_nll","b_nll",*this->get_b_likelihood(),*this->get_pseudo_data_bUP());
  test_Bup = 2*(sb_nll->getVal()-b_nll->getVal());

  sb_nll  = new RooNLLVar ("sb_nll","sb_nll",*this->get_sb_likelihood(),*this->get_pseudo_data_sb());
  b_nll   = new RooNLLVar ("b_nll","b_nll",*this->get_b_likelihood(),*this->get_pseudo_data_sb());
  test_SBmean = 2*(sb_nll->getVal()-b_nll->getVal());

  cout<<"Test statistics for data = "<<test_Data<<endl;
  cout<<"Test statistics for b = "<<test_Bmean<<endl;
  cout<<"Test statistics for b - 1sigma = "<<test_Bdown<<endl;
  cout<<"Test statistics for b + 1sigma = "<<test_Bup<<endl;
  cout<<"Test statistics for sb = "<<test_SBmean<<endl;

}

////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////  CLS LIMIT  CALCULATION METHODS ///////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////////////////////////////////////////////


void significance_Hybrid(Model* model,int n_toys,int random_seed=0){
  cout<<"///////////////////////////////////////////////////////////////////////////////////////////"<<endl;
  cout<<"Calculating significance with the Hybrid method"<<endl;
  cout<<"///////////////////////////////////////////////////////////////////////////////////////////"<<endl;
  
  //set the random seed
  RooRandom::randomGenerator()->SetSeed(random_seed);
    
   //get the calculator
  HybridCalculatorOriginal myhc(*model->get_data(),*model->get_sb_likelihood(),*model->get_b_likelihood());
  
  //for numbercounting experiments
  myhc.PatchSetExtended(false);
  
  //set likelihood ratio as the test statistics
  myhc.SetTestStatistic(1);
  
  //define the systematics to be used
  if (model->get_nuisance_set()) {
    myhc.UseNuisance(true);
    myhc.SetNuisancePdf(*model->get_nuisance_prior_pdf());
    myhc.SetNuisanceParameters(*model->get_nuisance_set());
  } else {
    myhc.UseNuisance(false);                            
  }
  
  //define the number of toys to be done
  myhc.SetNumberOfToys(n_toys);

  //get the Hypotestresult
  HybridResult* hcResult = myhc.GetHypoTest();
  double significance = hcResult->Significance();
  //double CLS= hcResult->Cls();
  //double CLSerror= hcResult->ClsError(); 

  //  myhc.PrintMore("");
  cout <<"HybridCalculator significance:" << significance<< endl;
}

HybridResult* UL_significance_Hybrid(Model* model,int n_toys,int random_seed=0,double sig=1){
  cout<<"///////////////////////////////////////////////////////////////////////////////////////////"<<endl;
  cout<<"Calculating significance with the Hybrid method"<<endl;
  cout<<"///////////////////////////////////////////////////////////////////////////////////////////"<<endl;
  
  RooMsgService::instance().setGlobalKillBelow(RooFit::FATAL);

  //set the random seed
  RooRandom::randomGenerator()->SetSeed(random_seed);
    
   //get the calculator
  HybridCalculatorOriginal myhc(*model->get_data(),*model->get_sb_likelihood(),*model->get_b_likelihood());
  
  //for numbercounting experiments
  myhc.PatchSetExtended(false);
  
  //set likelihood ratio as the test statistics
  myhc.SetTestStatistic(1);


  //define the systematics to be used
  if (model->get_nuisance_set()) {
    myhc.UseNuisance(true);
    myhc.SetNuisancePdf(*model->get_nuisance_prior_pdf());
    myhc.SetNuisanceParameters(*model->get_nuisance_set());
  } else {
    myhc.UseNuisance(false);                            
  }

  //  myhc.UseNuisance(false);


  //RooArgSet* poi= model->get_POI_set();
  //  poi->first()->Set(sig);
  RooRealVar* firstPOI = (RooRealVar*) model->get_POI_set()->first();
  firstPOI->setVal(sig);


  double expSigmas = fabs(model->get_teststat_Bmean() - model->get_teststat_SBmean())/fabs(model->get_teststat_Bmean() - model->get_teststat_Bup());

  cout<<"Expected Sigmas Significance from Test statistics  = "<<expSigmas<<endl;


  if(! isfinite(expSigmas) || expSigmas > 5.5) {
    cout<<" => BIG signal: Set CLs to 0!"<<endl;
    return 0;
  }


  //Test if we have a very small signal

  myhc.SetNumberOfToys(200);
  HybridResult* hcResult = myhc.GetHypoTest();
  hcResult->SetDataTestStatistics( model->get_teststat_Bdown() );

  if(hcResult->CLs() > 0.5) {
    cout<<"Small signal: Stop running toys after 200!!"<<endl;
    return hcResult;
  }

  delete hcResult;


  //define the number of toys to be done
  myhc.SetNumberOfToys(n_toys);
    
  return myhc.GetHypoTest();
 
  // hcResult->PrintMore("");

  // //  HypoTestPlot *plot = new HypoTestPlot(*hcResult,100);
  // HybridPlot* plot=hcResult->GetPlot("hcPlot","p Values Plot",100);
  // TCanvas *c1=new TCanvas;
  // plot->Draw();
  // c1->SaveAs("hybrid_Result.eps");

  // double significance = hcResult->Significance();
  // double CLS= hcResult->CLs();
  // double CLB= hcResult->CLb();
  // double CLsplusb= hcResult->CLsplusb();
  // double CLSerror= hcResult->CLsError(); 
  

  // // compute the mean expected significance from toys
  // double mean_sb_toys_test_stat = plot->GetSBmean();
  // hcResult->SetDataTestStatistics(mean_sb_toys_test_stat);
  // double toys_significance = hcResult->Significance();

  // cout<<endl;
  // cout<<"POI: "<<sig<<endl;
  // cout <<"significance of data: " << significance<<endl;
  // cout <<"mean significance of toys: " << toys_significance<<endl;
  // cout<<"CLs: "<<CLS<<endl;
  // cout<<"CLb: "<<CLB<<endl;
  // cout<<"CLsplusb: "<<CLsplusb<<endl;
  // cout<<"CLserror: "<<CLSerror<<endl;

  
}


void upper_limit_Hybrid_auto_scan(Model* model,double confidence,int n_toys_per_point,double lower_scan_bound, double upper_scan_bound,int random_seed=0,bool useCLS=true,double accuracy=0.005){
  cout<<"///////////////////////////////////////////////////////////////////////////////////////////"<<endl;
  cout<<"Calculating upper limit with the Hybrid method"<<endl;
  cout<<"///////////////////////////////////////////////////////////////////////////////////////////"<<endl;
  
  //set the random seed
  RooRandom::randomGenerator()->SetSeed(random_seed);
    
   //get the calculator
  HybridCalculatorOriginal myhc(*model->get_data(),*model->get_sb_likelihood(),*model->get_b_likelihood());
  
  //for numbercounting experiments
  myhc.PatchSetExtended(false);
  
  //set likelihood ratio as the test statistics
  myhc.SetTestStatistic(1);
  
  //define the systematics to be used
  if (model->get_nuisance_set()) {
    myhc.UseNuisance(true);
    myhc.SetNuisancePdf(*model->get_nuisance_prior_pdf());
    myhc.SetNuisanceParameters(*model->get_nuisance_set());
  } else {
    myhc.UseNuisance(false);                            
  }
  
  //define the number of toys to be done
  myhc.SetNumberOfToys(n_toys_per_point);
  
  //get the interval calculator
  HypoTestInverter myInverter(myhc,*model->get_POI());
  
  //use the CLS method
  myInverter.UseCLs(useCLS);
  
  //set the confidence
  myInverter.SetTestSize((1-confidence)*2);
  
  //run the auto scan in range lower_scan_bound - upper_scan_bound with accuracy parameter accuracy
  myInverter.RunAutoScan(lower_scan_bound,upper_scan_bound,myInverter.Size()/2,accuracy);
 
  //get the result
  HypoTestInverterResult* results = myInverter.GetInterval();
  
  double upperLimit = results->UpperLimit();
  std::cout <<confidence<< "% upper limit is: " <<  upperLimit << std::endl;
}


void upper_limit_Hybrid_fixed_scan(Model* model,double confidence,int n_toys_per_point,double lower_scan_bound, double upper_scan_bound,int n_steps,int random_seed=0,bool useCLS=true){
  cout<<"///////////////////////////////////////////////////////////////////////////////////////////"<<endl;
  cout<<"Calculating upper limit with the Hybrid method"<<endl;
  cout<<"///////////////////////////////////////////////////////////////////////////////////////////"<<endl;
  
  //set the random seed
  RooRandom::randomGenerator()->SetSeed(random_seed);
    
   //get the calculator
  HybridCalculatorOriginal myhc(*model->get_data(),*model->get_sb_likelihood(),*model->get_b_likelihood());
  
  //for numbercounting experiments
  myhc.PatchSetExtended(false);
  
  //set likelihood ratio as the test statistics
  myhc.SetTestStatistic(1);
  
  //define the systematics to be used
  if (model->get_nuisance_set()) {
    myhc.UseNuisance(true);
    myhc.SetNuisancePdf(*model->get_nuisance_prior_pdf());
    myhc.SetNuisanceParameters(*model->get_nuisance_set());
  } else {
    myhc.UseNuisance(false);                            
  }
  
  //define the number of toys to be done
  myhc.SetNumberOfToys(n_toys_per_point);
  
  //get the interval calculator
  HypoTestInverter myInverter(myhc,*model->get_POI());
  
  //use the CLS method
  myInverter.UseCLs(useCLS);
  
  //set the confidence
  myInverter.SetTestSize((1-confidence)*2);
  
  //run the fixed scan in range lower_scan_bound - upper_scan_bound with n_steps steps
  myInverter.RunFixedScan(n_steps,lower_scan_bound,upper_scan_bound);
 
  //get the result
  HypoTestInverterResult* results = myInverter.GetInterval();
  
  double upperLimit = results->UpperLimit();
  std::cout <<confidence<< "% upper limit is: " <<  upperLimit << std::endl;
}

void upper_limit_Hybrid_one_point(Model* model,double value_to_run_for,int n_toys,int random_seed=0,bool useCLS=true){
  cout<<"///////////////////////////////////////////////////////////////////////////////////////////"<<endl;
  cout<<"Calculating upper limit with the Hybrid method"<<endl;
  cout<<"///////////////////////////////////////////////////////////////////////////////////////////"<<endl;
  
  //set the random seed
  RooRandom::randomGenerator()->SetSeed(random_seed);
    
   //get the calculator
  HybridCalculatorOriginal myhc(*model->get_data(),*model->get_sb_likelihood(),*model->get_b_likelihood());
  
  //for numbercounting experiments
  myhc.PatchSetExtended(false);
  
  //set likelihood ratio as the test statistics
  myhc.SetTestStatistic(1);
  
  //define the systematics to be used
  if (model->get_nuisance_set()) {
    myhc.UseNuisance(true);
    myhc.SetNuisancePdf(*model->get_nuisance_prior_pdf());
    myhc.SetNuisanceParameters(*model->get_nuisance_set());
  } else {
    myhc.UseNuisance(false);                            
  }
  
  //define the number of toys to be done
  myhc.SetNumberOfToys(n_toys);
  
  //get the interval calculator
  HypoTestInverter myInverter(myhc,*model->get_POI());
  
  //use the CLS method
  myInverter.UseCLs(useCLS);
  
  //run one point
  myInverter.RunOnePoint(value_to_run_for);
 }
 
 void central_interval_Hybrid_fixed_scan(Model* model,double confidence,int n_toys_per_point,double lower_scan_bound, double upper_scan_bound,int n_steps,int random_seed=0,bool useCLS=true){
  cout<<"///////////////////////////////////////////////////////////////////////////////////////////"<<endl;
  cout<<"Calculating central interval with the Hybrid method"<<endl;
  cout<<"///////////////////////////////////////////////////////////////////////////////////////////"<<endl;
  
  //set the random seed
  RooRandom::randomGenerator()->SetSeed(random_seed);
    
   //get the calculator
  HybridCalculatorOriginal myhc(*model->get_data(),*model->get_sb_likelihood(),*model->get_b_likelihood());
  
  //for numbercounting experiments
  myhc.PatchSetExtended(false);
  
  //set likelihood ratio as the test statistics
  myhc.SetTestStatistic(1);
  
  //define the systematics to be used
  if (model->get_nuisance_set()) {
    myhc.UseNuisance(true);
    myhc.SetNuisancePdf(*model->get_nuisance_prior_pdf());
    myhc.SetNuisanceParameters(*model->get_nuisance_set());
  } else {
    myhc.UseNuisance(false);                            
  }
  
  //define the number of toys to be done
  myhc.SetNumberOfToys(n_toys_per_point);
  
  //get the interval calculator
  HypoTestInverter myInverter(myhc,*model->get_POI());
  
  //use the CLS method
  myInverter.UseCLs(useCLS);
  
  //set the confidence
  myInverter.SetTestSize(1-confidence);
  
  //run the auto fixed in range lower_scan_bound - upper_scan_bound with n_steps steps
  myInverter.RunFixedScan(n_steps,lower_scan_bound,upper_scan_bound);

 
  //get the result
  HypoTestInverterResult* results = myInverter.GetInterval();
  
  double upperLimit = results->UpperLimit();
  double lowerLimit = results->LowerLimit();
  std::cout <<confidence<<"% interval is: " << lowerLimit<<" , "<< upperLimit << std::endl;
}
 
 void central_interval_Hybrid_auto_scan(Model* model,double confidence,int n_toys_per_point,double lower_scan_bound, double upper_scan_bound,int random_seed=0,bool useCLS=true,double accuracy=0.005){
  cout<<"///////////////////////////////////////////////////////////////////////////////////////////"<<endl;
  cout<<"Calculating central interval with the Hybrid method"<<endl;
  cout<<"///////////////////////////////////////////////////////////////////////////////////////////"<<endl;
  
   //set the random seed
   RooRandom::randomGenerator()->SetSeed(random_seed);
  
  //get the HybridCalculatorOriginal   
  HybridCalculatorOriginal myhc(*model->get_data(),*model->get_sb_likelihood(),*model->get_b_likelihood());
  
  //for numbercounting experiments
  myhc.PatchSetExtended(false);
  
  //set likelihood ratio as the test statistics
  myhc.SetTestStatistic(1);
  
//   define the systematics to be used
  if (model->get_nuisance_set()) {
    myhc.UseNuisance(true);
    myhc.SetNuisancePdf(*model->get_nuisance_prior_pdf());
    myhc.SetNuisanceParameters(*model->get_nuisance_set());
  } else {
    myhc.UseNuisance(false);                            
  }
  
  //define the number of toys to be done
  myhc.SetNumberOfToys(n_toys_per_point);
  
  //get the interval calculator
  HypoTestInverter myInverter(myhc,*model->get_POI());
  
  //use the CLS method
  myInverter.UseCLs(useCLS);
  
  //set the confidence
  myInverter.SetTestSize(1-confidence);
  //scan the values for X to find the interval
  //requires the search range and the desired accuracy
  myInverter.RunAutoScan(lower_scan_bound,upper_scan_bound,myInverter.Size(),accuracy);
  //get the result
  HypoTestInverterResult* results = myInverter.GetInterval();
  
  double upperLimit = results->UpperLimit();
  double lowerLimit = results->LowerLimit();
  std::cout <<confidence<<"% interval is: " << lowerLimit<<" , "<< upperLimit << std::endl;
  
}


HypoTestResult* significance_Profile_Likelihood(Model* model){
  cout<<"///////////////////////////////////////////////////////////////////////////////////////////"<<endl;
  cout<<"Calculating significance with the Profile Likelihood method"<<endl;
  cout<<"///////////////////////////////////////////////////////////////////////////////////////////"<<endl;
  
  
  //      //define paramsOfInterest
     RooArgSet *paramsOfInterest=new RooArgSet("paramsOfInterest");
     paramsOfInterest->addClone(*model->get_POI());     //clone because we need s for complete likelihood
     
     //get the calculator
     ProfileLikelihoodCalculator plc(*model->get_data(),*model->get_complete_likelihood(), *paramsOfInterest);

     //define null hypothesis
     paramsOfInterest->setRealValue(model->get_POI()->GetName(),0);
     plc.SetNullParameters(*paramsOfInterest);
     
     //get the hypotest
     //    HypoTestResult* lrhypo=plc.GetHypoTest();
     //     double significance=lrhypo->Significance();     
     //     cout<<"Profile Likelihood significance is: "<<significance<<endl;


     return plc.GetHypoTest();
}


int main(int argc, char* argv[]){

  if (argc<=1){ 
    std::cerr<<"Error: Expected configs!"<<std::endl;
    exit(1);
  } else if (argc==2){ 
    std::cerr<<"Error: Expected background AND min 1 signal config!"<<std::endl;
    exit(1);
  } 


  bool writePlots = true;

  string bkgFileName=argv[1];
  std::cout << "opening background and data config file: "<<bkgFileName<<"\n"<<std::endl;
  cout<<"----------------------------\n"<<endl;
  //open the config file
  ConfigFile* bkgConfig = new ConfigFile(bkgFileName);

  int n_toys = bkgConfig->Value("definitions","n_toys");
  bool use_seperate_channels = bkgConfig->Value("definitions","use_seperate_channels",0);
  int n_sep_channels = 1;
  if(use_seperate_channels) n_sep_channels = bkgConfig->Value("definitions","n_channels");

  for (int file=2; file<argc; ++file){ 

    string fileName=argv[file];
    std::cout << "opening scan signal "<<file-1<<" from "<<argc-2<<": "<<fileName<<"\n"<<std::endl;


    for (int iChannnel=0; iChannnel<n_sep_channels; ++iChannnel){ 

      if(use_seperate_channels) std::cout << "Use seperated channels, now channel "<<iChannnel+1<<"/"<<n_sep_channels<<".\n"<<std::endl;

      ConfigFile* sigConfig = new ConfigFile(fileName);

      Model* model = 0;
      if(use_seperate_channels)    model = new Model(bkgConfig, sigConfig, iChannnel+1);
      else                         model = new Model(bkgConfig, sigConfig);

      //    model->Print();


      //       HypoTestResult* plcResult = significance_Profile_Likelihood(model);

      //       if(plcResult) {
      //        cout<<"\n ++++++++++"<<endl;
      //        plcResult->Print();
      //        cout<<"CL_s_Error: "<<plcResult->CLsError()<<endl;
      //        cout<<"plc significance: "<<plcResult->Significance()<<endl;
      //        cout<<"\n ++++++++++"<<endl;
      
      //       }

      cout<<"\n ++++++++++"<<endl;
      if(  model->get_lumi_scale() != 1) cout<<"Scale lumi with "<<model->get_lumi_scale()<<endl;
      cout<<"Signal= "<<model->get_sum_signal()<<endl;
      double soverb =model->get_sum_signal()/sqrt(model->get_sum_background());
      double soverberr = model->get_sum_signal()/sqrt( model->get_sum_background() + model->get_total_bunc() );
      cout<<"S / sqrt(B) = "<<soverb<<endl;
      cout<<"S / sqrt(B+unc) = "<<soverberr<<endl;


      int random_integer = rand();
      //      random_integer = 321;  //set to fixed value only for comparisions
      cout<<"random_integer: "<<random_integer<<endl;


      HybridResult* hcResult = UL_significance_Hybrid(model,n_toys,random_integer);

      double obsCLs = 0;
      double obsCLsError = 0;
      double obsSignificance = 1000;
      double expSignificance = 1000;
      double estimateSignificance = 1000;
      double expCLs = 0;
      double expCLsError = 0;
      double expCLsRMSDN = 0;
      double expCLsRMSUP = 0;

      HybridPlot* plot = 0;
      TCanvas *c1 = 0;


      if(hcResult) {


        vector<double> b  = hcResult->GetTestStat_b();
        vector<double> sb  = hcResult->GetTestStat_sb();

        double max = 0;
        double min = 0;

        for(unsigned int v=0;v<b.size()&&v<sb.size();v++) {
          if( max < b.at(v)  && isfinite(b.at(v)) )     max = b.at(v);
          if( min > sb.at(v) && isfinite(sb.at(v)) )    min = sb.at(v);
        }

        for(unsigned int v=0;v<b.size()&&v<sb.size();v++) {
          if(! isfinite(b.at(v)))       b.at(v) = max;
          if(! isfinite(sb.at(v)))     sb.at(v) = min;
        }

        HybridResult* corrResult = new HybridResult("HybridResult",sb,b);
        corrResult->SetDataTestStatistics( hcResult->GetTestStat_data() );


        plot=corrResult->GetPlot("hcPlot","p Values Plot",200);
        c1=new TCanvas;
        plot->Draw();


        cout<<"\n ++++++++++"<<endl;
        cout<<"Get center of B and SB model from the plot: "<<endl;
        double sbcenter = plot->GetSBCenter();
        double bcenter = plot->GetBCenter();
        double bUP = bcenter + plot->GetBrms();
        double bDN = bcenter - plot->GetBrms();


        if( sbcenter==0 ) sbcenter = model->get_teststat_SBmean();
        if( bcenter==0 ) {
          bcenter = model->get_teststat_Bmean();
          bUP = model->get_teststat_Bup();
          bDN =  model->get_teststat_Bdown();
        }

        estimateSignificance = fabs(bcenter-sbcenter)/fabs(bcenter-bUP);

        cout<<"SB center = "<<sbcenter<<" +/- "<<plot->GetSBrms()<<endl;
        cout<<"B center = "<<bcenter<<" +/- "<<plot->GetBrms()<<endl;
        cout<<"estimation of Significance: "<<estimateSignificance<<" sigma" <<endl;

        corrResult->PrintMore("");
        cout<<"   +/-  "<<corrResult->CLsError()<<endl;
        cout<<"data significance: "<<corrResult->Significance()<< " sigma" <<endl;

        // compute the mean expected significance from toys
        obsCLs =  corrResult->CLs();
        obsCLsError = corrResult->CLsError();
        obsSignificance = corrResult->Significance();


        corrResult->SetDataTestStatistics(sbcenter);

        if( corrResult->CLb() < 1 )
          expSignificance = corrResult->Significance();
        else
          expSignificance = 1000;

        //      cout<<"\n ++++++++++"<<endl;
        //      cout<<"Results from mean of SB toys:"<<endl;
        //      corrResult->PrintMore("");
        cout<<"Exp significance: "<<corrResult->Significance()<< " sigma" <<endl;

        corrResult->SetDataTestStatistics(bcenter);

        expCLs =  corrResult->CLs();
        expCLsError = corrResult->CLsError();

        //      cout<<"\n ++++++++++"<<endl;
        //      cout<<"Results from mean of B toys:"<<endl;
        //      corrResult->PrintMore("");

        cout<<"CL_s_exp: "<<expCLs<<endl;
        cout<<"    +/-   "<<expCLsError<<endl;


        corrResult->SetDataTestStatistics(bDN);

        expCLsRMSDN =  corrResult->CLs();
        cout<<"CL_s_expRMSdn: "<<expCLsRMSDN<<endl;
        cout<<"      +/-      "<< corrResult->CLsError()<<endl;

        corrResult->SetDataTestStatistics(bUP);

        expCLsRMSUP =  corrResult->CLs();
        cout<<"CL_s_expRMSup: "<<expCLsRMSUP<<endl;
        cout<<"      +/-      "<< corrResult->CLsError()<<endl;


        cout<<"\n ++++++++++"<<endl;

        if(writePlots) {
          size_t found = fileName.find_last_of("/");
          if(found == string::npos) found = 0;
          else found++;
          string shortFileName = fileName.substr(found,fileName.length());
    
          if(use_seperate_channels) c1->SaveAs(string(shortFileName.substr(0,shortFileName.find(".dat"))+"_channel"+int_to_string(iChannnel+1)+"_Result.eps").c_str());
          else                      c1->SaveAs(string(shortFileName.substr(0,shortFileName.find(".dat"))+"_Result.eps").c_str());
        }
        
        delete corrResult;
      }

      ///Write results to config

      ostringstream otext;


      if(!use_seperate_channels || iChannnel==0) {
        otext <<"[results_Hybrid]\n" 
              <<"BkgrConfig = "<<bkgFileName<<"\n"
              <<"n_toys = "<<n_toys<<"\n"
              <<"SoverB = "<<soverb<<"\n"
              <<"SoverBunc = "<<soverberr<<"\n"
              <<"SignalSum = "<<model->get_sum_signal()<<"\n"
              <<"BkgSum = "<<model->get_sum_background()<<"\n"
              <<"BkgUnc = "<<model->get_total_bunc()<<"\n";
      }


      if(!use_seperate_channels) {
        otext <<"obsCLs = "<<obsCLs<<"\n"
              <<"obsCLsError = "<<obsCLsError<<"\n"
              <<"obsSignificance = "<<obsSignificance<<"\n"
              <<"expCLs = "<<expCLs<<"\n"
              <<"expCLsError = "<<expCLsError<<"\n"
              <<"expSignificance = "<<expSignificance<<"\n"
              <<"expCLs_rmsDN = "<<expCLsRMSDN<<"\n"
              <<"expCLs_rmsUP = "<<expCLsRMSUP<<"\n"
              <<"estimateSignificance = "<<estimateSignificance<<"\n";
      } else {
        otext <<"obsCLs_ch"<<iChannnel+1<<" = "<<obsCLs<<"\n"
              <<"obsCLsError_ch"<<iChannnel+1<<" = "<<obsCLsError<<"\n"
              <<"obsSignificance_ch"<<iChannnel+1<<" = "<<obsSignificance<<"\n"
              <<"expCLs_ch"<<iChannnel+1<<" = "<<expCLs<<"\n"
              <<"expCLsError_ch"<<iChannnel+1<<" = "<<expCLsError<<"\n"
              <<"expSignificance_ch"<<iChannnel+1<<" = "<<expSignificance<<"\n"
              <<"expCLs_rmsDN_ch"<<iChannnel+1<<" = "<<expCLsRMSDN<<"\n"
              <<"expCLs_rmsUP_ch"<<iChannnel+1<<" = "<<expCLsRMSUP<<"\n"
              <<"estimateSignificance_ch"<<iChannnel+1<<" = "<<estimateSignificance<<"\n";
      }


      time_t rawtime;
      struct tm * timeinfo;
      time ( &rawtime );
      timeinfo = localtime ( &rawtime );
      ofstream ofile;
      ofile.open (fileName.c_str(),ios::app);
      if (ofile.good())
        std::cout << "adding results to '"<<fileName.c_str()<<"'"<<std::endl;
      else if ( (ofile.rdstate() & ifstream::failbit ) != 0 )
        std::cerr << "ERROR opening '"<<fileName.c_str()<<"'! Does the directory exist?"<<std::endl;


      if(!use_seperate_channels || iChannnel==0) ofile << "\n"<< asctime (timeinfo) <<"\n"<< otext.str()<<"\n"; 
      else                                       ofile << otext.str()<<"\n";


      ofile.close();


      if(plot) delete plot;
      if(c1)   delete c1;
      delete hcResult;
      delete sigConfig;
      delete model;
    }
  }
  

  //########################

  //   vector<double> *signal_expectation = new vector<double>;
  //   vector< vector<double>* > *background_expectation = new vector< vector<double>* >;
  //   vector< vector< vector<double>* >* >* systematic_uncertainties= new vector< vector< vector<double>* >* >;
  //   vector<double> *nuisance_upper_limit = new vector<double>;
  //   vector<int> *n_input= new vector<int>;

  //   int nuisance_range=10;
  //   int POI_range=20;
  //   int observation_range=500;


  //   //  signal_expectation->push_back(1.85562);
  //   signal_expectation->push_back(9.97288);
  //   signal_expectation->push_back(1.33057);
  //   signal_expectation->push_back(12.8879);


  //   std::vector<double>* bkg1 = new vector<double>;
  //   //  bkg1->push_back(68.0);
  //   bkg1->push_back(38.9);
  //   bkg1->push_back(8.0);
  //   bkg1->push_back(5.9);
  //   background_expectation->push_back(bkg1);


  //   std::vector<double>* unc1_sig = new vector<double>;
  //   std::vector<double>* unc2_sig = new vector<double>;
  //   std::vector<double>* unc1_bkg1 = new vector<double>;
  //   std::vector<double>* unc2_bkg1 = new vector<double>;

  //   std::vector< std::vector<double>* >* unc_type1 = new vector< vector<double>* >;
  //   std::vector< std::vector<double>* >* unc_type2 = new vector< vector<double>* >;

  //   //  unc1_sig->push_back(0.318846/1.85562);
  //   unc1_sig->push_back(1.55237/9.97288);
  //   unc1_sig->push_back(0.34926/1.33057);
  //   unc1_sig->push_back(2.2063/12.8879);


  //   unc_type1->push_back(unc1_sig);

  //   //  unc2_sig->push_back(0.0);
  //   unc2_sig->push_back(0.0);
  //   unc2_sig->push_back(0.0);
  //   unc2_sig->push_back(0.0);


  //   unc_type2->push_back(unc2_sig);


  //   // unc1_bkg1->push_back(13.6/68.0);
  //   // unc1_bkg1->push_back(8.6/38.9);
  //   // unc1_bkg1->push_back(1.2/8.0);
  //   // unc1_bkg1->push_back(1.4/5.9);

  //   //  unc1_bkg1->push_back(0.);
  //   unc1_bkg1->push_back(0.);
  //   unc1_bkg1->push_back(0.);
  //   unc1_bkg1->push_back(0.);

  //   unc_type1->push_back(unc1_bkg1);


  //   //  unc2_bkg1->push_back(13.6/68.0);
  //   unc2_bkg1->push_back(8.6/38.9);
  //   unc2_bkg1->push_back(1.2/8.0);
  //   unc2_bkg1->push_back(1.4/5.9);


  //   unc_type2->push_back(unc2_bkg1);

  //   systematic_uncertainties->push_back(unc_type1);

  //   systematic_uncertainties->push_back(unc_type2);


  //   nuisance_upper_limit->push_back(5);
  //   nuisance_upper_limit->push_back(5);


  //   //  n_input->push_back(65);
  //   n_input->push_back(31);
  //   n_input->push_back(6);
  //   n_input->push_back(9);


  //   Model* test=new Model(signal_expectation,background_expectation,systematic_uncertainties,nuisance_upper_limit,n_input,nuisance_range,POI_range,observation_range);

  //   //  Model* test=new Model(signal_expectation,background_expectation,systematic_uncertainties,nuisance_upper_limit,0,nuisance_range,POI_range,observation_range);

  //   test->Print();


  //   //  Model* test=new Model("RA7/test_60_230_COV.txt");


  //   test->get_b_likelihood()->getVariables()->Print();


  // //   TFile file("Model_results.root","recreate");

  // //   TH1F* hist_b_likelihood_0 =  (TH1F*) test->get_b_likelihood()->createHistogram("observation_0",100);
  // //   TH1F* hist_b_likelihood_1 =  (TH1F*) test->get_b_likelihood()->createHistogram("observation_1",100);
  // //   TH1F* hist_b_likelihood_2 =  (TH1F*) test->get_b_likelihood()->createHistogram("observation_2",100);
  // //   TH1F* hist_b_likelihood_3 =  (TH1F*) test->get_b_likelihood()->createHistogram("observation_3",100);

  // //   TH1F* hist_sb_likelihood_0 =  (TH1F*) test->get_sb_likelihood()->createHistogram("observation_0",100);
  // //   TH1F* hist_sb_likelihood_1 =  (TH1F*) test->get_sb_likelihood()->createHistogram("observation_1",100);
  // //   TH1F* hist_sb_likelihood_2 =  (TH1F*) test->get_sb_likelihood()->createHistogram("observation_2",100);
  // //   TH1F* hist_sb_likelihood_3 =  (TH1F*) test->get_sb_likelihood()->createHistogram("observation_3",100);

  // //   hist_b_likelihood_0->Write();
  // //   hist_b_likelihood_1->Write();
  // //   hist_b_likelihood_2->Write();
  // //   hist_b_likelihood_3->Write();

  // //   hist_sb_likelihood_0->Write();
  // //   hist_sb_likelihood_1->Write();
  // //   hist_sb_likelihood_2->Write();
  // //   hist_sb_likelihood_3->Write();

  // //   file.Close();

  // //  significance_Hybrid(test,5000,321);


  //   UL_significance_Hybrid(test,10000,321);
  //   //  upper_limit_Hybrid_auto_scan(test,0.95,50,0.01,5);


  // //   systematic_uncertainties

  // //   for (int file=1; file<argc; ++file){ 

  // //     string fileName=argv[file];
  // //     std::cout << "opening: "<<fileName<<std::endl;

  // //     ConfigFile config(argv[file]);


  // //     Model* test=new Model(signal_expectation,background_expectation,systematic_uncertainties,n_input,nuisance_range,POI_range,observation_range);
 

  // //     significance_Hybrid(test,1000,321);


  // //     delete test;

  // //     string out="out_"+fileName;
  // //     ofstream outfile(out.c_str());
  // //     outfile<<"##outputfile created by Model.cc from: "<<fileName<<endl;

  // //     outfile.close();

  // //     std::cout << "##outputfile created: "<<out<<std::endl;

  // //   }


  //   //  Model* test=new Model(file.c_str());
  //   //Model* test=new Model("test_60_230.txt");
  //   //Model* test=new Model("test_140_200.txt");
  //   //Model* test=new Model("MultiLepton_Model_Jt200_CU_DAVIS_ML01.txt");
  
  
  //   //   srand((unsigned)time(0));
  //   //   int random_integer = rand();
  //    //cout<<random_integer<<endl;
  //    //   cout<<"random_number: "<<"XXXseedXXX"<<endl;

  //   //   RooRandom::randomGenerator()->SetSeed(random_integer);


  return(0);
  
}
Revision:	1.1.1.1 (vendor branch)
Committed:	Fri Jun 17 13:09:26 2011 UTC (13 years, 10 months ago) by tschum
Content type:	text/plain
Branch:	tschum, MAIN
CVS Tags:	start, HEAD
Changes since 1.1:	+0 -0 lines
Log Message:	Evolved standalone RA7 limit calculation tool