/// \file
/// \ingroup tutorial_roofit
/// \notebook -js
/// Basic functionality: fitting, plotting, toy data generation on one-dimensional p.d.f
///
///  pdf = gauss(x,m,s)
///
/// \macro_image
/// \macro_output
/// \macro_code
/// \author 07/2008 - Wouter Verkerke

#include <fstream>
#include "RooRealVar.h"
#include "RooDataSet.h"
#include "RooGaussian.h"
#include "TCanvas.h"
#include "RooPlot.h"
#include "TAxis.h"
#include "TFile.h"
#include "TApplication.h"
#include "TROOT.h"
#include "RooAbsGenContext.h"
#include "RooRandom.h"
#include "RooAbsNumGenerator.h"
#include "RooErrorHandler.h"


using namespace RooFit;

using namespace std;

// Gaussian distribution constants//
const Double_t kC1 = 1.448242853;
const Double_t kC2 = 3.307147487;
const Double_t kC3 = 1.46754004;
const Double_t kD1 = 1.036467755;
const Double_t kD2 = 5.295844968;
const Double_t kD3 = 3.631288474;
const Double_t kHm = 0.483941449;
const Double_t kZm = 0.107981933;
const Double_t kHp = 4.132731354;
const Double_t kZp = 18.52161694;
const Double_t kPhln = 0.4515827053;
const Double_t kHm1 = 0.516058551;
const Double_t kHp1 = 3.132731354;
const Double_t kHzm = 0.375959516;
const Double_t kHzmp = 0.591923442;
/*zhm 0.967882898*/

const Double_t kAs = 0.8853395638;
const Double_t kBs = 0.2452635696;
const Double_t kCs = 0.2770276848;
const Double_t kB  = 0.5029324303;
const Double_t kX0 = 0.4571828819;
const Double_t kYm = 0.187308492 ;
const Double_t kS  = 0.7270572718 ;
const Double_t kT  = 0.03895759111;
RooDataSet *data;


RooRealVar x_gaus("x_gaus", "x_gaus", -10, 10);
RooRealVar mean("mean", "mean of gaussian", 1, -10, 10);
RooRealVar sigma("sigma", "width of gaussian", 1, 0.1, 10);
RooGaussian *gauss = new RooGaussian("gauss", "gaussian PDF", x_gaus, mean, sigma);


int main(int argc, char **argv) {

    if(argc==1){
        printf("Usage  : %s  <inputfile-data> <inputfile-source>  <imagefile> \n",argv[0]);
        printf("Assumes there is a soft link called DATA to these files \n");
        printf("To exit, select File->Quit from the menu in the ROOT window. \n");
        exit(0);
    }

    // Extract input file name 1 from the argument list
    char input_file_1[100];   // input file
    strcpy(input_file_1,argv[1]);
    // Extract input file name 1 from the argument list
    char input_file_2[100];   // input file
    strcpy(input_file_2,argv[2]);

    // Extract ROOT output file name from the argument list
    char image_file[64];   // root file
    strcpy(image_file,argv[3]);

    gROOT->Reset();
    TApplication *TheApp = new TApplication("DAQanalisys",&argc,argv);

    // Make a soft link in the project directory to the data called DATA
    TString PathToData = "DATA/";
    TString input_1 = PathToData+input_file_1;
    TString input_2 = PathToData+input_file_2;
    cout << "\n"<<input_1<<" and "<<input_2<<" and "<<image_file<<" confirmed "<<endl;

    // Construct plot frame in 'x'
    RooPlot *xframe = x_gaus.frame(Title("Gaussian p.d.f. with while(1)"));

    // P l o t   m o d e l   a n d   c h a n g e   p a r a m e t e r   v a l u e s
    // ---------------------------------------------------------------------------

    TFile outputFile("filename.root", "RECREATE");
//    RooAbsData::setDefaultStorageType(RooAbsData::Tree);

    // Plot gauss in frame (i.e. in x_gaus)
    gauss->plotOn(xframe);

    // Change the value of sigma to 3
    sigma.setVal(3);

    // Plot gauss in frame (i.e. in x_gaus) and draw frame on canvas
    gauss->plotOn(xframe, LineColor(kRed));

    // G e n e r a t e   e v e n t s
    // -----------------------------
    bool verbose = 0;
    bool autoBinned = 0;
    const char* binnedTag = "";

    data = new RooDataSet("Gauss data","Generated from gaus data",x_gaus) ;
    data->convertToTreeStore();
    Double_t xgen,residual_y, residual_x,residual_z,residual_rn;

    // specify number of events to generate
    long long nEvents = 1000000000;

    // arrays containing random numbers
//    double y_arr[nEvents+148];
//    double x_arr[nEvents];
//    double yw_arr[nEvents];
//    double rn_arr[nEvents];

    // variables used to generate gaussian random numbers
    Double_t result;
    Double_t rn,x,y,z,x1,y1;
    Double_t repeat_x;
    Double_t repeat_y;
    int countin_y;
    int countin_x;
    int countin_z;
    int end;


    //variables to read random numbers into text files
    ofstream y_txt,x_txt,rn_txt,yw_txt,c3d2e3_txt,c5_txt,c5d1_txt,c5d2_txt,c5d3_txt,c5d4_txt
    ,c5d4e1_txt,c3d2e2_txt,c3d2e1_txt,c3d2_txt,c3d1_txt,c3_txt,c2_txt,c1_txt,c0_txt;
//    c1_txt.open("/Users/cashcrusaders/Documents/UJ_LimitSetting/bin/check.txt");
//    x_txt.open("/home/xola1/UJ_LimitSetting/src/x.txt");
//    rn_txt.open("/home/xola1/UJ_LimitSetting/src/rn.txt");
//    yw_txt.open("/home/xola1/UJ_LimitSetting/src/yw.txt");
//    y_txt.open("/home/xola1/UJ_LimitSetting/src/y.txt");

    // variables to count number of times code goes into if statements
    Int_t c3d2e3,c5,c5d1,c5d2,c5d3,c5d4,c5d4e1,c3d2e2,c3d2e1,c3d2,c3d1,c4,c3,c2,c1,c0,c5d4e2,c5d5,c3d2e4;

    // Generate a dataset of 1000 events in x_gaus from gauss
    for(long long evt =0;evt<nEvents;evt++)
    {   while(1)
        {
            do
            {
                y = gRandom->Rndm();
//                y_txt<<y<<endl;
//                y= y_arr[evt];
                if (y>kHm1)
                {
                    //************c1*****************//
                    c1++;
                    rn = 0;
//                    rn_txt<<rn<<endl;
                    result = kHp*y-kHp1;
                    break;
                }

                else if (y<kZm)
                {
                    //************c2*****************//
                    c2++;
                    rn = 0;
//                    rn_txt<<rn<<endl;
                    rn = kZp*y-1;
                    result = (rn>0) ? (1+rn) : (-1+rn);
                    break;
                }
                else if (y<kHm)
                {
                    //************c3*****************//
                    c3++;
                    rn = gRandom->Rndm();
//                    rn_arr[evt] = rn;
//                    rn_txt<<rn<<endl;
                    rn = rn-1+rn;
                    z = (rn>0) ? 2-rn : -2-rn;
                    if ((kC1-y)*(kC3+TMath::Abs(z))<kC2)
                    {
                        //************c3d1*****************//
                        c3d1++;
                        result = z; break; }
                    else
                    {
                        //************c3d2*****************//
                        c3d2++;
                        x = rn*rn;
                        if ((y+kD1)*(kD3+x)<kD2)
                        {
                            //************c3d2e1*****************//
                            c3d2e1++;
                            result = rn; break; }
                        else if (kHzmp-y<exp(-(z*z+kPhln)/2))
                        {
                            //************c3d2e2*****************//
                            c3d2e2++;
                            result = z; break; }
                        else if (y+kHzm<exp(-(x+kPhln)/2))
                        {
                            //************c3d2e3*****************//
                            c3d2e3++;
                            result = rn;
                            break;
                        }
                        else
                        {
                            c3d2e4++;
                        }
                    }
                }
                else
                {
                    c4++;
                    rn = 0;
//                    rn_txt<<rn<<endl;
                }
                while (1)
                {
                    //************c5*****************//
                    c5++;
                    x = gRandom->Rndm();
                    x1 = x;
                    y = kYm * gRandom->Rndm();
                    y1 = y;
                    z = kX0 - kS*x - y;
                    if (z>0)
                    {
                        //************c5d1*****************//
                        c5d1++;
                        rn = 2 + y / x;
                    }
                    else
                    {
                        //************c5d2*****************//
                        c5d2++;
                        x = 1-x;
                        y = kYm-y;
                        rn = -(2+y/x);
                    }
                    if ((y-kAs+x)*(kCs+x)+kBs<0)
                    {
                        //************c5d3*****************//
                        c5d3++;
                        cout << "extracted number "<<x1<<"and the event number "<<evt<<endl;
                        result = rn; break; }
                    else if (y<x+kT)
                    {
                        //************c5d4*****************//
                        c5d4++;
                        if (rn * rn < 4 * (kB - log(x)))
                        {
                            //************c5d4e1*****************//
                            c5d4e1++;
                            result = rn;
//                            cout << "extracted number "<<x1<<"and the event number "<<evt<<endl;
                            break;
                        }
                        else
                        {
                            c5d4e2++;
                        }
                    }
                    else
                    {
                        c5d5++;
                    }
//                    x_arr[evt] = x1;
//                    yw_arr[evt] = y1;
//                    x_txt<<x1<<endl;
//                    yw_txt<<y1<<endl;
//                    repeat_y = yw_arr[evt];
//                    repeat_x = x_arr[evt];
                    countin_x++;
                    countin_y++;
                }
//                x_arr[evt] = x1;
//                yw_arr[evt] = y1;
//                x_txt<<x1<<endl;
//                yw_txt<<y1<<endl;
//                repeat_y = yw_arr[evt];
//                repeat_x = x_arr[evt];
                countin_x++;
                countin_z++;
            } while(0);
            //meant to end loop as soon as it ends

            xgen = gauss->mean + gauss->sigma * result;
            if (xgen < gauss->x.max() && xgen > gauss->x.min()) {
                gauss->x = xgen;
//                end++;
                break;
            }else{
                c0++;
//                break;
            }
        }

        // WVE add check that event is in normRange

        data->addFast(x_gaus);
    }

    // Make a second plot frame in x_gaus and draw both the
    // data1 and the p.d.f in the frame
    cout << "the is the number of times c0 is called "<<end<< endl;
    cout << "the is the number of times c1 is called "<<c1<< endl;
    cout << "the is the number of times c2 is called "<<c2<< endl;
    cout << "the is the number of times c3 is called "<<c3<< endl;
    cout << "the is the number of times c3d1 is called "<<c3d1<< endl;
    cout << "the is the number of times c3d2 is called "<<c3d2<< endl;
    cout << "the is the number of times c3d2e1 is called "<<c3d2e1<< endl;
    cout << "the is the number of times c3d2e2 is called "<<c3d2e2<< endl;
    cout << "the is the number of times c3d2e3 is called "<<c3d2e3<< endl;
    cout << "the is the number of times c3d2e4 is called "<<c3d2e4<< endl;
    cout << "the is the number of times c4 is called "<<c4<< endl;
    cout << "the is the number of times c5 is called "<<c5<< endl;
    cout << "the is the number of times c5d1 is called "<<c5d1<< endl;
    cout << "the is the number of times c5d2 is called "<<c5d2<< endl;
    cout << "the is the number of times c5d3 is called "<<c5d3<< endl;
    cout << "the is the number of times c5d4 is called "<<c5d4<< endl;
    cout << "the is the number of times c5d5 is called "<<c5d5<< endl;
    cout << "the is the number of times c5d4e1 is called "<<c5d4e1<< endl;
    cout << "the is the number of times c5d4e2 is called "<<c5d4e2<< endl;
    cout << "the is the number of times repetitions in while loop is called "<<countin_y<< endl;
    cout << "the is the number of times repetitions in while loop is called "<<countin_z<< endl;


    RooPlot *xframe2 = x_gaus.frame(Title("Gaussian p.d.f. without while(1)"));
    data->plotOn(xframe2);

    // F i t   m o d e l   t o   d a t a
    // -----------------------------
    // Fit pdf to data1
//    gauss->fitTo(*data);

    // Print values of mean and sigma (that now reflect fitted values and errors)
    mean.Print();
    sigma.Print();

    // Draw all frames on a canvas
    TCanvas *c = new TCanvas("rf101_basics", "rf101_basics", 800, 400);
//   c->Divide(2);
//   c->cd(1);
//   gPad->SetLeftMargin(0.15);
//   xframe->GetYaxis()->SetTitleOffset(1.6);
//   xframe->Draw();
//   c->cd(2);
//    y_txt.close();
//    yw_txt.close();
//    x_txt.close();
//    rn_txt.close();
    data->statOn(xframe2);

    gPad->SetLeftMargin(0.15);
    xframe2->GetYaxis()->SetTitleOffset(1.6);
    xframe2->Draw();
    TFile f("gausplot","recreate");
    c->Update();
    f.Write();
    TheApp->Run();
}