#include <iostream>
#include <cmath>
#include "TString.h"
#include "TFile.h"
#include "TCut.h"

int Slices(TString inputFile = TString())
{
//     gROOT->SetBatch(kTRUE);
        
    if (inputFile.IsNull()) { cout << "No input file is given" << endl; return 1;}

    TFile *f = gROOT->GetListOfFiles()->FindObject(inputFile);
    if (!f) { f = new TFile(inputFile); }

    TCanvas *c1 = gROOT->GetListOfCanvases()->FindObject("c1");
    if (!c1) { c1 = new TCanvas("c1","c1",800,600); }

    ofstream fout("Slices.dat");

	TPostScript *ps = new TPostScript("Slices.ps",112);

    const int ibins = 50;
    const int NLayers = 4;
    const int NSlices = 8;
    float minZetLim[NLayers] = {-120.0, -120.0, -120.0, -120.0};
    float maxZetLim[NLayers] = {+120.0, +120.0, +120.0, +120.0};
    
    float sliceWidth[NLayers]; for (int i=0; i<NLayers; i++) { sliceWidth[i]=(maxZetLim[i]-minZetLim[i])/NSlices; }

    float minSlice[NLayers][NSlices];
    float maxSlice[NLayers][NSlices];
    float midSlice[NLayers][NSlices];

    char VrblSliceZETA[NLayers][NSlices][300], CutsSliceZETA[NLayers][NSlices][300], TitleSliceZETA[NLayers][NSlices][300];
    char VrblSliceRPHI[NLayers][NSlices][300], CutsSliceRPHI[NLayers][NSlices][300], TitleSliceRPHI[NLayers][NSlices][300];

    float minhistZETA[NLayers] = {-20.0, -10.0, -5.0, -15.0};
    float maxhistZETA[NLayers] = {+20.0, +10.0, +5.0, +15.0};

    float minhistRPHI[NLayers] = {-20.0, -10.0, -5.0, -15.0};
    float maxhistRPHI[NLayers] = {+20.0, +10.0, +5.0, +15.0};

    float ZETAmean[NLayers][NSlices], ZETAsigma[NLayers][NSlices], sigZETAmean[NLayers][NSlices], sigZETAsigma[NLayers][NSlices];
    float RPHImean[NLayers][NSlices], RPHIsigma[NLayers][NSlices], sigRPHImean[NLayers][NSlices], sigRPHIsigma[NLayers][NSlices];

    //  CC Data        
    float RCCFLR[5]  = { 72.9600, 85.3442, 87.1001, 91.4303, 99.2269 };     // Radii of the floors PS, EM1-EM4 for CC.                  (.tz files)
    
    for (int k=0; k<NLayers; k++)
    {
        for (int i=0; i<NSlices; i++)
        {
            minSlice[k][i] = minZetLim[k]+i*sliceWidth[k];
            maxSlice[k][i] = minZetLim[k]+(i+1)*sliceWidth[k];
            midSlice[k][i] = 0.5*(minSlice[k][i]+maxSlice[k][i]);

			// This part is for the residuals on ZETA.			
            ps->NewPage();
            sprintf(VrblSliceZETA[k][i],"(floor_z[0][%i]-hits_z[0][%i])>>histZETA",k+1,k+1);            
            sprintf(CutsSliceZETA[k][i],"floor_z[0][%i]>%f && floor_z[0][%i]<%f && i_id1==0 &&  i_id2==1 && hasTrackMatch[0]==1 && hits_z[0][%i]!=9999.0",k+1,minSlice[k][i],k+1,maxSlice[k][i],k+1);            
            sprintf(TitleSliceZETA[k][i],"floor_z[0][%i]>%f && floor_z[0][%i]<%f",k+1,minSlice[k][i],k+1,maxSlice[k][i]);

	TH1F *histZETA = new TH1F("histZETA",TitleSliceZETA[k][i],ibins,minhistZETA[k],maxhistZETA[k]);

	myTree->Draw(VrblSliceZETA[k][i],CutsSliceZETA[k][i]);	
            
            histZETA->Fit("gaus","ME"); c1->Update();
           
            ZETAmean[k][i]  = histZETA->GetFunction("gaus")->GetParameter(1);
            ZETAsigma[k][i] = histZETA->GetFunction("gaus")->GetParameter(2);
            
            sigZETAmean[k][i]  = histZETA->GetFunction("gaus")->GetParError(1);
            sigZETAsigma[k][i] = histZETA->GetFunction("gaus")->GetParError(2);
            
			// This part is for the residuals on RPHI.			
            ps->NewPage();
            sprintf(VrblSliceRPHI[k][i],"%f*(atan(floor_x[0][%i],floor_y[0][%i])-atan(hits_x[0][%i]-hits_y[0][%i]))>>histRPHI",RCCFLR[k+1],k+1,k+1,k+1,k+1);
            sprintf(CutsSliceRPHI[k][i],"floor_z[0][%i]>%f && floor_z[0][%i]<%f && i_id1==0 &&  i_id2==1 && hasTrackMatch[0]==1 && hits_z[0][%i]!=9999.0",k+1,minSlice[k][i],k+1,maxSlice[k][i],k+1);            
            sprintf(TitleSliceRPHI[k][i],"floor_z[0][%i]>%f && floor_z[0][%i]<%f",k+1,minSlice[k][i],k+1,maxSlice[k][i]);
            
	TH1F *histRPHI = new TH1F("histRPHI",TitleSliceRPHI[k][i],ibins,minhistRPHI[k],maxhistRPHI[k]);
            
	myTree->Draw(VrblSliceRPHI[k][i],CutsSliceRPHI[k][i]);	
            
            histRPHI->Fit("gaus", "ME"); c1->Update();

            RPHImean[k][i]  = histRPHI->GetFunction("gaus")->GetParameter(1);
            RPHIsigma[k][i] = histRPHI->GetFunction("gaus")->GetParameter(2);
            
            sigRPHImean[k][i]  = histRPHI->GetFunction("gaus")->GetParError(1);
            sigRPHIsigma[k][i] = histRPHI->GetFunction("gaus")->GetParError(2);
            
            fout << midSlice[k][i] << " " << ZETAmean[k][i] << " " << sigZETAmean[k][i] << endl;
            fout << midSlice[k][i] << " " << ZETAsigma[k][i] << " " << sigZETAsigma[k][i] << endl;

            fout << midSlice[k][i] << " " << RPHImean[k][i] << " " << sigRPHImean[k][i] << endl;
            fout << midSlice[k][i] << " " << RPHIsigma[k][i] << " " << sigRPHIsigma[k][i] << endl;
            
        }
    }

    ps->Close();

}