using namespace RooStats;
void poissonRuns_HipoInverter(){
  const char* infile =  "poissonRuns.out_separate.root"; 
  const char* workspaceName = "w";
  const char* modelConfigName = "Smodel";
  const char* dataName = "data";

  /////////////////////////////////////////////////////////////
  // First part is just to access the workspace file 
  ////////////////////////////////////////////////////////////

  // Check if example input file exists
  TFile *file = TFile::Open(infile);

  // get the workspace out of the file
  RooWorkspace* w = (RooWorkspace*) file->Get(workspaceName);

  // get the modelConfig out of the file and the data
  ModelConfig* sbModel = (ModelConfig*) w->obj(modelConfigName);
  RooAbsData* data = w->data(dataName) ;

  //construct the B hypo (s=0)
  RooStats::ModelConfig* bModel = (RooStats::ModelConfig*) sbModel->Clone("Bmodel") ;
  RooRealVar* poi = (RooRealVar*) bModel->GetParametersOfInterest()->first();
  poi->setVal(0);
  bModel->SetSnapshot(*poi);
   

  RooStats::AsymptoticCalculator  asympCalc(*data, *bModel, *sbModel);
  asympCalc.SetOneSided(true);

  RooStats::HypoTestInverter inverter(asympCalc);
  inverter.SetConfidenceLevel(0.90);
  inverter.UseCLs(true);

  inverter.SetVerbose(false);
  inverter.SetFixedScan(500,1.0,200.0); // set number of points , xmin and xmax

  RooStats::HypoTestInverterResult* result =  inverter.GetInterval();

  



  const char* infile_tot =  "poissonRuns.out_tot.root"; 
  const char* workspaceName_tot = "w2";
  const char* modelConfigName_tot = "Smodel_tot";
  const char* dataName_tot = "data_tot";

  /////////////////////////////////////////////////////////////
  // First part is just to access the workspace file 
  ////////////////////////////////////////////////////////////

  // Check if example input file exists
  TFile *file_tot = TFile::Open(infile_tot);

  // get the workspace out of the file
  RooWorkspace* w2 = (RooWorkspace*) file_tot->Get(workspaceName_tot);

  // get the modelConfig out of the file and the data
  ModelConfig* sbModel_tot = (ModelConfig*) w2->obj(modelConfigName_tot);
  RooAbsData* data_tot = w2->data(dataName_tot) ;

  //construct the B hypo (s=0)
  RooStats::ModelConfig* bModel_tot = (RooStats::ModelConfig*) sbModel_tot->Clone("Bmodel_tot") ;
  RooRealVar* poi_tot = (RooRealVar*) bModel_tot->GetParametersOfInterest()->first();
  poi_tot->setVal(0);
  bModel_tot->SetSnapshot(*poi_tot);
   

  RooStats::AsymptoticCalculator  asympCalc_tot(*data_tot, *bModel_tot, *sbModel_tot);
  asympCalc_tot.SetOneSided(true);

  RooStats::HypoTestInverter inverter_tot(asympCalc_tot);
  inverter_tot.SetConfidenceLevel(0.90);
  inverter_tot.UseCLs(true);

  inverter_tot.SetVerbose(false);
  inverter_tot.SetFixedScan(500,1.0,200.0); // set number of points , xmin and xmax

  RooStats::HypoTestInverterResult* result_tot =  inverter_tot.GetInterval();





  
  /*
  RooStats::FrequentistCalculator  freqCalc(*data, *bModel, *sbModel);
  RooStats::ProfileLikelihoodTestStat* plr = new RooStats::ProfileLikelihoodTestStat(*sbModel->GetPdf());
  plr->SetOneSided(true);

  RooStats::ToyMCSampler* toymcs = (RooStats::ToyMCSampler*) freqCalc.GetTestStatSampler();
  toymcs->SetTestStatistic(plr);
  if (!sbModel->GetPdf()->canBeExtended()) {
    toymcs->SetNEventsPerToy(1);
  }
  freqCalc.SetToys(40000,40000);
  RooStats::HypoTestInverter inverter2(freqCalc);
  inverter2.SetConfidenceLevel(0.90);
  inverter2.UseCLs(false);
  inverter2.SetVerbose(false);
  inverter2.SetFixedScan(20,0.1,14.0); // set number of points , xmin and xmax
  RooStats::HypoTestInverterResult* result2 =  inverter2.GetInterval();*/
  

  cout<<"RUN_BY_RUN"<<endl;
  cout << 100*inverter.ConfidenceLevel() << "%  upper limit : " << result->UpperLimit() << endl;
  std::cout << "Expected upper limits, using the B (alternate) model : " << std::endl;
  std::cout << " expected limit (median) " << result->GetExpectedUpperLimit(0) << std::endl;
  std::cout << " expected limit (-1 sig) " << result->GetExpectedUpperLimit(-1) << std::endl;
  std::cout << " expected limit (+1 sig) " << result->GetExpectedUpperLimit(1) << std::endl;
  std::cout << " expected limit (-2 sig) " << result->GetExpectedUpperLimit(-2) << std::endl;
  std::cout << " expected limit (+2 sig) " << result->GetExpectedUpperLimit(2) << std::endl;

  cout<<"TOTAL"<<endl;
  cout << 100*inverter_tot.ConfidenceLevel() << "%  upper limit : " << result_tot->UpperLimit() << endl;
  std::cout << "Expected upper limits, using the B (alternate) model : " << std::endl;
  std::cout << " expected limit (median) " << result_tot->GetExpectedUpperLimit(0) << std::endl;
  std::cout << " expected limit (-1 sig) " << result_tot->GetExpectedUpperLimit(-1) << std::endl;
  std::cout << " expected limit (+1 sig) " << result_tot->GetExpectedUpperLimit(1) << std::endl;
  std::cout << " expected limit (-2 sig) " << result_tot->GetExpectedUpperLimit(-2) << std::endl;
  std::cout << " expected limit (+2 sig) " << result_tot->GetExpectedUpperLimit(2) << std::endl;

  
  /* cout << 100*inverter2.ConfidenceLevel() << "%  upper limit : " << result2->UpperLimit() << endl;
  std::cout << "Expected upper limits, using the B (alternate) model : " << std::endl;
  std::cout << " expected limit (median) " << result2->GetExpectedUpperLimit(0) << std::endl;
  std::cout << " expected limit (-1 sig) " << result2->GetExpectedUpperLimit(-1) << std::endl;
  std::cout << " expected limit (+1 sig) " << result2->GetExpectedUpperLimit(1) << std::endl;
  std::cout << " expected limit (-2 sig) " << result2->GetExpectedUpperLimit(-2) << std::endl;
  std::cout << " expected limit (+2 sig) " << result2->GetExpectedUpperLimit(2) << std::endl;*/
  

  TCanvas* c1 = new TCanvas();
  RooStats::HypoTestInverterPlot* plot = new RooStats::HypoTestInverterPlot("HTI_Result_Plot_runBYrun","HypoTest Scan Result run-by-run",result);
  plot->Draw("CLb 2CL");  // plot also CLb and CLs+b
  c1->Draw();
 
  /*TCanvas* c2 = new TCanvas();
  RooStats::HypoTestInverterPlot* plot2 = new RooStats::HypoTestInverterPlot("HTI_Result_Plot","HypoTest Scan Result",result2);
  plot2->Draw("CLb 2CL");  // plot also CLb and CLs+b
  c2->Draw();*/

  TCanvas* c1_tot = new TCanvas();
  RooStats::HypoTestInverterPlot* plot_tot = new RooStats::HypoTestInverterPlot("HTI_Result_Plot_tot","HypoTest Scan Result tot",result_tot);
  plot_tot->Draw("CLb 2CL");  // plot also CLb and CLs+b
  c1_tot->Draw();
}