//-----------------------------------------------------------------------
//
// Convoluted Landau and Gaussian Fitting Function
//         (using ROOT's Landau and Gauss functions)
//
//  Based on a Fortran code by R.Fruehwirth (fruhwirth@hephy.oeaw.ac.at)
//  Adapted for C++/ROOT by H.Pernegger (Heinz.Pernegger@cern.ch) and
//   Markus Friedl (Markus.Friedl@cern.ch)
//
//  to execute this example, do:
//  root > .x langaus.C
// or
//  root > .x langaus.C++
//
//-----------------------------------------------------------------------

#include "../../Tools/Lib.C"
#include "../../Tools/lhcbStyle.C"
#include "../../Tools/Style.C"
#include "../../Tools/FindAndReplaceString.C"

Double_t langaufun(Double_t *x, Double_t *par) {

    //Fit parameters:
    //par[0]=Width (scale) parameter of Landau density
    //par[1]=Most Probable (MP, location) parameter of Landau density
    //par[2]=Total area (integral -inf to inf, normalization constant)
    //par[3]=Width (sigma) of convoluted Gaussian function
    //
    //In the Landau distribution (represented by the CERNLIB approximation),
    //the maximum is located at x=-0.22278298 with the location parameter=0.
    //This shift is corrected within this function, so that the actual
    //maximum is identical to the MP parameter.

    // Numeric constants
    Double_t invsq2pi = 0.3989422804014;   // (2 pi)^(-1/2)
    Double_t mpshift  = -0.22278298;       // Landau maximum location

    // Control constants
    Double_t np = 100.0;      // number of convolution steps
    Double_t sc =   5.0;      // convolution extends to +-sc Gaussian sigmas

    // Variables
    Double_t xx;
    Double_t mpc;
    Double_t fland;
    Double_t sum = 0.0;
    Double_t xlow,xupp;
    Double_t step;
    Double_t i;


    // MP shift correction
    mpc = par[1] - mpshift * par[0];

    // Range of convolution integral
    xlow = x[0] - sc * par[3];
    xupp = x[0] + sc * par[3];

    step = (xupp-xlow) / np;

    // Convolution integral of Landau and Gaussian by sum
    for(i=1.0; i<=np/2; i++) {
        xx = xlow + (i-.5) * step;
        fland = TMath::Landau(xx,mpc,par[0]) / par[0];
        sum += fland * TMath::Gaus(x[0],xx,par[3]);

        xx = xupp - (i-.5) * step;
        fland = TMath::Landau(xx,mpc,par[0]) / par[0];
        sum += fland * TMath::Gaus(x[0],xx,par[3]);
    }

    return (par[2] * step * sum * invsq2pi / par[3]);
}



TF1 *langaufit(TH1F *his, Double_t *fitrange, Double_t *startvalues, Double_t *parlimitslo, Double_t *parlimitshi, Double_t *fitparams, Double_t *fiterrors, Double_t *ChiSqr, Int_t *NDF)
{
    // Once again, here are the Landau * Gaussian parameters:
    //   par[0]=Width (scale) parameter of Landau density
    //   par[1]=Most Probable (MP, location) parameter of Landau density
    //   par[2]=Total area (integral -inf to inf, normalization constant)
    //   par[3]=Width (sigma) of convoluted Gaussian function
    //
    // Variables for langaufit call:
    //   his             histogram to fit
    //   fitrange[2]     lo and hi boundaries of fit range
    //   startvalues[4]  reasonable start values for the fit
    //   parlimitslo[4]  lower parameter limits
    //   parlimitshi[4]  upper parameter limits
    //   fitparams[4]    returns the final fit parameters
    //   fiterrors[4]    returns the final fit errors
    //   ChiSqr          returns the chi square
    //   NDF             returns ndf

    Int_t i;
    Char_t FunName[100];

    sprintf(FunName,"Fitfcn_%s",his->GetName());

    TF1 *ffitold = (TF1*)gROOT->GetListOfFunctions()->FindObject(FunName);
    if (ffitold) delete ffitold;

    TF1 *ffit = new TF1(FunName,langaufun,fitrange[0],fitrange[1],4);
    ffit->SetParameters(startvalues);
    ffit->SetParNames("Width","MP","Area","GSigma");

    for (i=0; i<4; i++) {
        ffit->SetParLimits(i, parlimitslo[i], parlimitshi[i]);
    }

    his->Fit(FunName,"RB0");   // fit within specified range, use ParLimits, do not plot

    ffit->GetParameters(fitparams);    // obtain fit parameters
    for (i=0; i<4; i++) {
        fiterrors[i] = ffit->GetParError(i);     // obtain fit parameter errors
    }
    ChiSqr[0] = ffit->GetChisquare();  // obtain chi^2
    NDF[0] = ffit->GetNDF();           // obtain ndf

    return (ffit);              // return fit function

}


Int_t langaupro(Double_t *params, Double_t &maxx, Double_t &FWHM) {

    // Seaches for the location (x value) at the maximum of the
    // Landau-Gaussian convolute and its full width at half-maximum.
    //
    // The search is probably not very efficient, but it's a first try.

    Double_t p,x,fy,fxr,fxl;
    Double_t step;
    Double_t l,lold;
    Int_t i = 0;
    Int_t MAXCALLS = 10000;


    // Search for maximum

    p = params[1] - 0.1 * params[0];
    step = 0.05 * params[0];
    lold = -2.0;
    l    = -1.0;


    while ( (l != lold) && (i < MAXCALLS) ) {
        i++;

        lold = l;
        x = p + step;
        l = langaufun(&x,params);

        if (l < lold)
            step = -step/10;

        p += step;
    }

    if (i == MAXCALLS)
        return (-1);

    maxx = x;

    fy = l/2;


    // Search for right x location of fy

    p = maxx + params[0];
    step = params[0];
    lold = -2.0;
    l    = -1e300;
    i    = 0;


    while ( (l != lold) && (i < MAXCALLS) ) {
        i++;

        lold = l;
        x = p + step;
        l = TMath::Abs(langaufun(&x,params) - fy);

        if (l > lold)
            step = -step/10;

        p += step;
    }

    if (i == MAXCALLS)
        return (-2);

    fxr = x;


    // Search for left x location of fy

    p = maxx - 0.5 * params[0];
    step = -params[0];
    lold = -2.0;
    l    = -1e300;
    i    = 0;

    while ( (l != lold) && (i < MAXCALLS) ) {
        i++;

        lold = l;
        x = p + step;
        l = TMath::Abs(langaufun(&x,params) - fy);

        if (l > lold)
            step = -step/10;

        p += step;
    }

    if (i == MAXCALLS)
        return (-3);


    fxl = x;

    FWHM = fxr - fxl;
    return (0);
}

int N = 512;
// string filename = "~flionett/work/LHCb/UTTestBeam/Repository/KeplerDev_v3r0/output/Run_Bias_Scan-B6-A-299-8431_Tuple.root";
// string filename = "~flionett/work/LHCb/UTTestBeam/Repository/KeplerDev_v3r0/output/Plots/AnalysisOutput_A6_300_1_8431.root";
// int channel = -1;

void TbUTLandauFits(string filename, int channel);

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

    if ((argc ==2) && (string(argv[1]) == "--info"))
    {
        cout << "**************************************************" << endl;
        cout << "Put comments here." << endl;
        cout << "**************************************************" << endl;
        return 0;
    }
    else if (argc < 3)
    {
        cout << "**************************************************" << endl;
        cout << "Error! Arguments missing..." << endl;
        cout << "Please use the following format:" << endl;
        cout << "./TbUTLandauFits [1] [2]" << endl;
        cout << "with:" << endl;
        cout << "[1] = filename (full path);" << endl;
        cout << "[2] = channel." << endl;
        cout << "Type ./TbUTLandauFits --info for more information." << endl;
        cout << "**************************************************" << endl;
        return 0;
    }
    else if (argc == 3)
    {
        cout << "**************************************************" << endl;
        cout << "Filename: " << argv[1] << endl;
        cout << "Channel: " << argv[2] << endl;
        cout << "**************************************************" << endl;
        TbUTLandauFits(string(argv[1]),atoi(argv[2]));
        return 0;
    }
    else 
    {
        cout << "**************************************************" << endl;
        cout << "Error! Too many arguments provided..." << endl;
        cout << "**************************************************" << endl;
        return 0;
    }
}

// The <filename> variable is the filename with the complete path.
void TbUTLandauFits(string filename, int channel) {
    TString inFilename = filename.substr(0,filename.find(".root")); // No file extension.
    cout << inFilename << endl;
    TString outFilename = filename.substr(0,filename.find(".root"))+"_LandauFits"; // No file extension.
    cout << outFilename << endl;

    TString pathToInput = filename.substr(0,filename.find_last_of("/"));
    cout << pathToInput << endl;
    TString pathToOutput = filename.substr(0,filename.find_last_of("/"))+"/LandauFits";
    cout << pathToOutput << endl;
    TString pathToFigures = filename.substr(0,filename.find_last_of("/"))+"/LandauFits/Figures";

    // Create output folders if they do not exist.
    TString pathToOutputMake = "mkdir -p "+pathToOutput;
    cout << "Create folder " << pathToOutput << " where output will be saved" << endl;
    system(string(pathToOutputMake).c_str());
    TString path_to_make_figures = "mkdir -p "+pathToFigures;
    system(string(path_to_make_figures).c_str());
    cout << "Setting figures to: " << path_to_make_figures << endl;

    TFile *fIn = new TFile(inFilename+".root","READ");
    TFile *fOut = new TFile(outFilename+".root","RECREATE");

    TH1F *hS = 0;
    TH1F *hS1ch[N];

    if (channel == -1) {
        TList *list = new TList;
        for (int i=0; i<N; ++i) {
            // Convert the integer to string.
            stringstream s;
            s << i;
            string ch(s.str());

            hS1ch[i] = (TH1F *)fIn->Get("hlandau_"+TString(ch)); 
            InitHist(hS1ch[i],Form("Signal distribution on channel %d",i),"ADC counts","");
            list->Add(hS1ch[i]);
        }
        hS=(TH1F *)hS1ch[0]->Clone("hlandau_0");
        InitHist(hS,"Signal distribution","ADC counts","");
        hS->Reset();
        hS->Merge(list);
        /*
           TTree *t = (TTree *)fIn->Get("TbUT/Clusters");
           vector<double> clustersCharge(10);

           t->SetBranchAddress("clustersCharge",&(&clustersCharge[10]));
           if (fIn->GetListOfKeys()->Contains("TbUT/Clusters/clustersCharge")) {


           hS = (TH1F *)fIn->Get("TbUT/Clusters/clustersCharge");
           }
           else {
           cout << "ERROR! No input object found." << endl;
           return ;
           }
           */
    }
    else {
        // Convert the integer to string.
        stringstream s;
        s << channel;
        string ch(s.str());

        hS = (TH1F *)fIn->Get("hlandau_"+TString(ch));
        InitHist(hS,"Signal distribution","ADC counts","");
        if (hS->GetEntries() == 0) {
            cout << "ERROR! The selected histogram is empty." << endl;
            return;
        }
    }

    // Fitting S histo
    printf("Fitting...\n");

    // Setting fit range and start values
    Double_t fr[2];
    Double_t sv[4], pllo[4], plhi[4], fp[4], fpe[4];
    fr[0]=0.4*hS->GetMean();
    fr[1]=3.0*hS->GetMean();

    pllo[0]=0.; pllo[1]=100.0; pllo[2]=1.; pllo[3]=0.01;
    plhi[0]=50.; plhi[1]=350.0; plhi[2]=100000000.; plhi[3]=50.;
    sv[0]=30.; sv[1]=250.; sv[2]=500.; sv[3]=3.;

    Double_t chisqr;
    Int_t    ndf;
    TF1 *fitsnr = langaufit(hS,fr,sv,pllo,plhi,fp,fpe,&chisqr,&ndf);

    Double_t SPeak, SFWHM;
    langaupro(fp,SPeak,SFWHM);

    printf("Fitting done\nPlotting results...\n");

    // Global style settings
    gStyle->SetOptStat(1111);
    gStyle->SetOptFit(111);
    gStyle->SetLabelSize(0.03,"x");
    gStyle->SetLabelSize(0.03,"y");

    // hS->GetXaxis()->SetRange(0,70);
    // hS->Draw();
    // fitsnr->Draw("lsame");

    TCanvas *cS = new TCanvas("cS","",800,600);
    cS->cd();
    fitsnr->SetLineColor(kRed);
    DrawHistFunc(cS,hS,fitsnr,pathToFigures);

    fOut->cd();
    hS->Write();
    fIn->Close();
    fOut->Close();

    return;
}