// merge  the HTI results 
#include "TFile.h"
#include "TMath.h"
#include "TError.h"
#include "TCanvas.h"


#include "RooStats/HypoTestInverterResult.h"
#include "RooStats/HypoTestInverterPlot.h"
#include "RooStats/SamplingDistPlot.h"
#include "RooStats/RooStatsUtils.h"
#include "RooDataSet.h"
#include "TSystem.h"

#include <iostream>


//gSystem->Load("libRooStats");
using namespace RooStats;

// plot detailed output for one parameter + test statistics
bool PlotDetailedOutput(HypoTestResult * r, const char * poiName, double poiValue) { 

   RooDataSet * nullData = r->GetNullDetailedOutput();
   RooDataSet * altData = r->GetAltDetailedOutput();
   if (!nullData || !altData) { 
      Error("PlotDetailedOutput","Detailed output is not present in result file");
      return false;
   }
   
   TCanvas * c2 = new TCanvas(); 
   c2->SetTitle(TString::Format("Output for poi = %f",poiValue) );
   //RooArgList pars(*nullData->get(0));
   TString muName1 = TString("ModelConfig_TS0_fitUncond_")+TString(poiName);
   TString muName0 = TString("ModelConfig_with_poi_0_TS0_fitUncond_")+TString(poiName);
   TString tsName1 = TString("ModelConfig_TS0");
   TString tsName0 = TString("ModelConfig_with_poi_0_TS0");

   TH1 * h1 = new TH1D("h1","POI S+B",100,1,0);
   TH1 * h0 = new TH1D("h0","POI B",100,1,0);
   TH1 * ts1 = new TH1D("ts1","TS S+B",100,1,0);
   TH1 * ts0 = new TH1D("ts1","TS B",100,1,0);

   // fill the histograms
   for (int i = 0; i < nullData->numEntries() ; ++i) { 
      const RooArgSet * args = nullData->get(i);
      h1->Fill(args->getRealValue(muName1, -999 ) );
      ts1->Fill(args->getRealValue(tsName1, -999 ) );
   }
   for (int i = 0; i < altData->numEntries() ; ++i) { 
      const RooArgSet * args = altData->get(i);
      h0->Fill(args->getRealValue(muName0, -999 ) );
      ts0->Fill(args->getRealValue(tsName0, -999 ) );
   }

//    TH1 * h1 = nullData->createHistogram(muName1, (RooRealVar&)pars[1],Binning(100,0.5,1.5));
//    TH1 * h0 = altData->createHistogram(muName0, (RooRealVar&)pars[1],Binning(100,-.1,.1));
// //     TH1 * hy = nullData->createHistogram("muy", (RooRealVar&)pars[],Binning(100,-5,5));
//    TH1 * ts1 = nullData->createHistogram(tsName1, (RooRealVar&)pars[0],Binning(100,0,10) );
//    TH1 * ts0 = altData->createHistogram(tsName0, (RooRealVar&)pars[0],Binning(100,0,10) );
     c2->Divide(2,2);
     c2->cd(1);
     h1->Draw();
     //h1->Fit("gaus","L");
     c2->cd(2);
     h0->Draw();
     //h0->Fit("gaus","L");
     c2->cd(3);
     ts1->Draw();
     TF1 * f1 = new TF1("f1","[0]*ROOT::Math::chisquared_pdf(2*x,[1])");
     f1->SetParameters(1,2);
     //ts1->Fit(f1,"L");
     c2->cd(4);
     ts0->Draw();

     return true;
}

// plot histograms for all detailed outputs
bool PlotAllDetailedOutput(HypoTestResult * r, double poiValue, bool nullOutput) {
   RooDataSet * data = 0; 
   if (nullOutput) 
      data = r->GetNullDetailedOutput();
   else 
      data = r->GetAltDetailedOutput();

   if (!data) { 
      Error("PlotDetailedOutput","Detailed output is not present in result file");
      return false;
   }
   

   // loop on the data variables 

   const RooArgSet * vars = data->get(); 

   TCanvas * c2 = new TCanvas(); 
   c2->SetTitle(TString::Format("Output for poi = %f",poiValue) );
   int n = vars->getSize(); 
   int ny = TMath::CeilNint(TMath::Sqrt(n));
   int nx = TMath::CeilNint(double(n)/ny);
   c2->Divide( nx,ny);

   TFile * f = new TFile("detailedOutputHistograms.root","UPDATE");
   
   TString type = (nullOutput) ? "null" : "alt"; 
   
   RooArgList lvars(*vars);
   for (int i = 0; i < n; ++i) { 
      RooAbsArg & arg = lvars[i]; 
      TString hname = TString::Format("%s/poi_%f/h_%s",type.Data(),poiValue,arg.GetName() ); 
      RooAbsRealLValue * var = dynamic_cast<RooAbsRealLValue*>(&arg);
      //TH1 * h1 = data->createHistogram(hname, *var, RooFit::AutoBinning(100,0.05) );
      TH1 * h1 = new TH1D(hname,hname,100,1,0); 
      for (int j = 0; j < data->numEntries(); ++j) {
         const RooArgSet * ev = data->get(j); 
         h1->Fill(ev->getRealValue(arg.GetName() ) );
      }
      c2->cd(i+1);
      h1->Draw();
      h1->Write();
   }
   gPad->Update();

   // get also the fit information
   RooAbsData * fitData = r->GetFitInfo(); 
   if (!fitData) return false; 
   const RooArgSet * fitvars = fitData->get(0);
   // loop on all variable with null or all 

   RooArgList fitVarList(*fitvars);
   if (fitVarList.getSize() < 1) return false; 
   
   // make a label histogram with the fit data information
   TString hname1 = TString::Format("fitData_%s_poi_%f",type.Data(),poiValue ); 
   TH1D * hf = new TH1D(hname1,hname1,fitvars->getSize()/2,0,1);

   for (int i = 0; i < fitVarList.getSize(); ++i) { 
      RooAbsArg & arg = fitVarList[i];
      TString name = arg.GetName(); 
      if (nullOutput && name.Contains("fitNull") ) { 
         name.ReplaceAll("fitNull_","");
      }
      if (nullOutput && name.Contains("fitAlt") ) { 
         name.ReplaceAll("fitAlt_","");
      }
      double value = ( (RooRealVar &) arg).getVal();
      double error = ( (RooRealVar &) arg).getError();
      int bin = hf->Fill(name,value);
      hf->SetBinError(bin,error);
      hf->SetMinimum(-3); hf->SetMaximum(3);         
   }
   hf->Write();


   f->Close();
   return true;



}

// make a TTree with the detailed output
bool makeTree(HypoTestResult * r, double poiValue, bool nullOutput) {
   RooDataSet * data = 0; 
   if (nullOutput) 
      data = r->GetNullDetailedOutput();
   else 
      data = r->GetAltDetailedOutput();

   if (!data) { 
      Error("PlotDetailedOutput","Detailed output is not present in result file");
      return false;
   }

   TFile * f = new TFile("detailedOutputTree.root","UPDATE");
   
   TString type = (nullOutput) ? "null" : "alt"; 

   // save also all variables as a TTree
   TString tname1=TString::Format("tree_%s_poi_%f",type.Data(),poiValue);
   TString tname2=TString::Format("fitdata_%s_poi_%f",type.Data(),poiValue);

   TTree * tree1 = RooStats::GetAsTTree(tname1,tname1, *data); 

   tree1->Write();

   // get also the fit information
   RooDataSet * fitData = r->GetFitInfo(); 
   if (!fitData) return false; 

   TTree * tree2 = RooStats::GetAsTTree(tname2,tname2, *fitData); 

   tree2->Write();
   f->Close();

   return true; 
}
   
void readHypoTestInvDetailedOutput(const char * file1, const char *poiName = "s", bool plotAllPar = false, bool nullOutput = true) {


   TFile * f1 = new TFile(file1);
   TString resultName = "result_";
   resultName += poiName;
   HypoTestInverterResult * r1 = (HypoTestInverterResult*) f1->Get(resultName);

   if (!r1) { 
      Error("read_CLs","result with name %s is not found in file ",resultName.Data() );
      f1->ls();
      return;
   }


   if (plotAllPar) { 
      // remove old files 
      gSystem->Exec("rm detailedOutputHistograms.root");
      gSystem->Exec("rm detailedOutputTree.root");

   }

   // loop on the points
   for (int i = 0; i < r1->ArraySize(); ++i) {

      HypoTestResult * r = r1->GetResult(i); 

      bool ok = false; 
      if (!plotAllPar)
         ok = PlotDetailedOutput(r, poiName, r1->GetXValue(i)); 
      else 
         ok = PlotAllDetailedOutput(r, r1->GetXValue(i),nullOutput);          

      if (!ok) return;

      makeTree(r, r1->GetXValue(i),nullOutput); 


   }


   
}