#include <TROOT.h>
#include <TFile.h>
#include <fstream>
#include <TString.h>
#include <TH1.h>
#include <TObjString.h>
#include <TCut.h>
#include <TTree.h>
#include <TCanvas.h>
#include <iostream>
#include "TConfidenceLevel.h"
#include "TLimit.h"
#include "TLimitDataSource.h"

#define maxchan 20
#define maxmass 20
#define precision 50000
#define outtxtfile "hwwcombined.limits"
#define outrootfile "hwwcombined.root"

int Combine(TString MassFile="mass.list",TString ListFile="channels.list", Float_t xsfl=1., Float_t xsfh=1., Float_t xsfs=1.)
{
	// first open the ListFile to get the list of channels
	ifstream infile1((const char*)ListFile);
	if (!infile1) { cerr << "Fichier channels non existant" << endl; return -1; }
	TString filelist[maxchan];
	Int_t line = 0;
	while(infile1 >> filelist[line++])
		if(line>=maxchan) { cout << "Too many channels, using only " << maxchan << endl; break; }
	infile1.close();
	Int_t nch = line-1;
	// then open the MassFile to get the list of masses
	ifstream infile2((const char*)MassFile);
	if (!infile2) { cerr << "Fichier mass non existant" << endl; return -1; }
	Float_t mass[maxmass] = {0};
	line = 0;
	while(infile2 >> mass[line++])
		if(line>=maxmass) { cout << "Too many masses, using only " << maxmass << endl; break; }
	infile2.close();
	Int_t nma = line -1;
	// open the output text file and output ROOT file
	ofstream outfile(outtxtfile);
	TFile *rootfile = new TFile(outrootfile,"RECREATE");
	TTree *roottree = new TTree("hwwlimits","H->WW limits");
	Float_t cls;
	roottree->Branch("cls",&cls,"F");
	Float_t clb;
	roottree->Branch("clb",&clb,"F");
	Float_t ecls;
	roottree->Branch("ecls",&ecls,"F");
	Float_t eclb;
	roottree->Branch("eclb",&eclb,"F");
	Float_t amass;
	roottree->Branch("mass",&amass,"F");
	Float_t xsec;
	roottree->Branch("xsec_factor",&xsec,"F");
	TConfidenceLevel *myconfidence = new TConfidenceLevel;
	roottree->Branch("ConfidenceLevel","TConfidenceLevel",&myconfidence,32000,99);
	roottree->Branch("ConfidenceLevel_ns","TConfidenceLevel",&myconfidence,32000,0);
	delete myconfidence;
	// big loop over the masses
	for(Int_t midx=0;midx<nma;midx++)
	{
		// big loop over the cross-sections
		for(Float_t i=xsfl;i<=xsfh;i+=xsfs)
		{
			TLimitDataSource* mydatasource = new TLimitDataSource();
			mydatasource->SetOwner(true);
			// loop on the files and add the channels
			for(Int_t j = 0; j< nch; j++)
			{
				// open the file
				TFile theFile(filelist[j],"READ");
				theFile.cd();
				// get the data
				TH1F* tempo_data = (TH1F*)theFile.Get("data");
				tempo_data->SetDirectory(0);
				// get the background
				TH1F* tempo_bg = (TH1F*)theFile.Get("bg");
				tempo_bg->SetDirectory(0);
				// get the signal
				char sname[20];
				sprintf(sname,"signal%d",midx+1);
				TH1F* tempo_sig = (TH1F*)theFile.Get(sname);
				tempo_sig->SetDirectory(0);
				tempo_sig->Scale(i);
				tempo_sig->Sumw2();
				// add the channel
				mydatasource->AddChannel(tempo_sig,tempo_bg,tempo_data);
				// close the file
				theFile.Close();
			}
			rootfile->cd();
			myconfidence = TLimit::ComputeLimit(mydatasource,precision);
			ecls  = myconfidence->GetExpectedCLs_b();
	                eclb  = myconfidence->GetExpectedCLb_b();
			cls   = myconfidence->CLs();
	                clb   = myconfidence->CLb();
			amass = mass[midx];
			xsec  = i;
			outfile <<amass << " " << xsec << " " << cls << " " << clb << " " << ecls << " " << eclb << endl;
			roottree->Fill();
			delete mydatasource;
			delete myconfidence;
		}
	}
	// conclusion (close all what remains open)
	outfile.close();
	rootfile->cd();
	roottree->Write();
	rootfile->Close();
	return 0;
}