#include "TMinuit.h"
#include "TMatrix.h"
#include "TString.h"
#include "TCanvas.h"
#include "TH1F.h"
#include "TF1.h"
#include "TLegend.h"
#include "TFile.h"
#include "TArrayF.h"
#include "TGraphErrors.h"
#include "MRawRunHeader.h"
#include "TPad.h"
#include "TPaveText.h"
#include "iostream.h"
#include "TObjArray.h"

#include "MHMatrix.h"

//--------------------------------------------------------------------------------------
// MAIN INPUT/OUTPUT

#include "IOMkn421.h"
//--------------------------------------------------------------------------------------

Double_t acut = 10; // alpha region of excess
Bool_t aplot  = 1;  // make only alpha plot (1) , do sig. optimization and alpha plot (0)

//--------------------------------------------------------------------------------------------------------------
//--------------------------------------------------------------------------------------------------------------

const TMatrix *matOn;
const TMatrix *matOff;

Int_t isize   = 0;
Int_t idist   = 1;
Int_t iwidth  = 2;
Int_t ilength = 3;
Int_t ialpha  = 4;

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------

Double_t signal(Double_t *x, Double_t *par){
    return par[0]/(TMath::Sqrt(2.*TMath::Pi())*par[1])*TMath::Exp(-(x[0]-par[2])*(x[0]-par[2])/2./par[1]/par[1]);
}

Double_t background(Double_t *x, Double_t *par){
    return par[0]+par[1]*x[0]+par[2]*pow(x[0],2.)+par[3]*x[0]*x[0]*x[0];
    //return par[0]+par[1]*x[0]+par[2]*x[0]*x[0]+par[3]*x[0]*x[0]*x[0];
}

Double_t fitFunction(Double_t *x, Double_t *par){
    return signal(x,par) + background(x,&par[3]);
}

//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
{
    Double_t w1=par[0];
    Double_t w2=par[1];
    Double_t l1=par[2];
    Double_t l2=par[3];
    Double_t d1=par[4];
    Double_t d2=par[5];

    Double_t non=0;
    Double_t noff=0;

    TH1F hon("hon","",18,0.,90.);//5 deg bins
    for(Int_t i=0;i<matOn->GetNrows();i++)
    {
        if( (*matOn)(i,isize)<SizeLo) continue;
        if(((*matOn)(i,idist)>d1)  ||((*matOn)(i,idist)<d2  )) continue;

        if(((*matOn)(i,iwidth)>w1) ||((*matOn)(i,iwidth)<w2))  continue;
        if(((*matOn)(i,ilength)>l1)||((*matOn)(i,ilength)<l2)) continue;

        Float_t alpha=(*matOn)(i,ialpha);
        hon.Fill(TMath::Abs(alpha));
    }

    for(Int_t i=1;i<=4;i++)
        non+=hon.GetBinContent(i);

    //-----------------------------------------------------------------
    /*TH1F hoff("hoff","",18,0.,90.);//5 deg bins
    for(Int_t i=0;i<matOff->GetNrows();i++)
    {
        if( (*matOff)(i,isize)<SizeLo) continue;
        if(((*matOff)(i,idist)>d1)  ||((*matOff)(i,idist)<d2  )) continue;

        if(((*matOff)(i,iwidth)>w1) ||((*matOff)(i,iwidth)<w2))  continue;
        if(((*matOff)(i,ilength)>l1)||((*matOff)(i,ilength)<l2)) continue;

        Float_t alpha=(*matOff)(i,ialpha);
        hoff.Fill(TMath::Abs(alpha));
    }

    for(Int_t i=1;i<=4;i++)
        noff+=hoff.GetBinContent(i);*/
    //-----------------------------------------------------------------

    // flat background region
    for(Int_t i=7;i<=15;i++)
        noff+=hon.GetBinContent(i);

    Float_t offscale=4./9.;
    noff*=offscale;

    if(non+noff<0.001)
        f=0.;
    else
        f = -(non-noff)/(TMath::Sqrt(non+noff));

    cout<<"Non="<<non<<" Noff="<<noff*offscale<<" Sig="<<-f<<endl;
}
//------------------------------------------------------------------------------
//------------------------------------------------------------------------------
void CutOptim()
{
    //-----------------------------------------------------------------------
    // read matrices

    MHMatrix hmatOn;
    TFile *fileOn=new TFile(fileMatOn.Data());
    hmatOn.Read("Matrix");
    delete fileOn;
    const TMatrix &mOn=hmatOn.GetM();
    matOn=(const TMatrix*)&mOn;

    MHMatrix hmatOff;
    TFile *fileOff=new TFile(fileMatOff.Data());
    hmatOff.Read("Matrix");
    delete fileOff;
    const TMatrix &mOff=hmatOff.GetM();
    matOff=(const TMatrix*)&mOff;

    //-----------------------------------------------------------------------

    Int_t nrowOn  = matOn->GetNrows();
    Int_t nrowOff = matOff->GetNrows();

    Int_t ncolOn  = matOn->GetNcols();
    Int_t ncolOff = matOff->GetNcols();

    cout<<"matrix On:  ncol="<<ncolOn<<" nrow="<<nrowOn<<endl;
    cout<<"matrix Off: ncol="<<ncolOff<<" nrow="<<nrowOff<<endl;

    TCanvas *c = new TCanvas("c1","Alpha distribution after g/h separation",200,10,800,500);
    c->SetFillColor(0);

    //-----------------------------------------------------------------------
    // optimize cuts

    Int_t iflag;
    Double_t val;
    Double_t vstart[6];

    Double_t *gin=NULL;

    Double_t arglist[10];
    Int_t ierflg = 0;

    TMinuit *gMinuit = new TMinuit(10);
    TGraph* gr=NULL;

    if(!aplot)
        c->Divide(2,1);
    else
        goto APlot;

    // start values
    vstart[0]=WidthUp;
    vstart[1]=WidthLo;
    vstart[2]=LengthUp;
    vstart[3]=LengthLo;
    vstart[4]=DistUp;
    vstart[5]=DistLo;

    gMinuit->SetFCN(fcn);

    arglist[0] = 1;
    gMinuit->mnexcm("SET ERR", arglist ,1,ierflg);

    // Set starting values and step sizes for parameters
    static Double_t step[6];
    const Double_t def_step=0.2;
    step[0]=def_step;
    step[1]=def_step;
    step[2]=def_step;
    step[3]=def_step;
    step[4]=def_step;
    step[5]=def_step;

    gMinuit->mnparm(0,"w1", vstart[0], step[0], 0.,0.,ierflg);
    gMinuit->mnparm(1,"w2", vstart[1], step[1], 0.,0.,ierflg);
    gMinuit->mnparm(2,"l1", vstart[2], step[2], 0.,0.,ierflg);
    gMinuit->mnparm(3,"l2", vstart[3], step[3], 0.,0.,ierflg);
    gMinuit->mnparm(4,"d1", vstart[4], step[4], 0.2,1.,ierflg);
    gMinuit->mnparm(5,"d2", vstart[5], step[5], 0.5,1.5,ierflg);

    //gMinuit->FixParameter(0);  // WidthUp
    gMinuit->FixParameter(1);   // WidthLo
    //gMinuit->FixParameter(2);  // LengthUp
    gMinuit->FixParameter(3);   // LengthLo
    gMinuit->FixParameter(4);   // DistUp
    gMinuit->FixParameter(5);   // DistLo

    gMinuit->Eval(6,gin,val,vstart,iflag);

    arglist[0] = 5000;
    arglist[1] = 1.;

    //gMinuit->mnexcm("MIGRAD", arglist ,2,ierflg);
    gMinuit->mnexcm("SIM", arglist ,2,ierflg);

    // Print results
    Double_t amin,edm,errdef;
    Int_t nvpar,nparx,icstat;
    gMinuit->mnstat(amin,edm,errdef,nvpar,nparx,icstat);
    gMinuit->mnprin(3,amin);

    Double_t fpar[6];
    Double_t fparerr[6];
    for(Int_t i=0;i<6;i++)
        gMinuit->GetParameter(i,fpar[i],fparerr[i]);

    WidthUp  = fpar[0];
    WidthLo  = fpar[1];
    LengthUp = fpar[2];
    LengthLo = fpar[3];
    DistUp   = fpar[4];
    DistLo   = fpar[5];


    // parameter's sigma plots
    c->cd(2);
    // 2 sigma contour
    /*
    gMinuit->SetErrorDef(4);
    gr = (TGraph*) gMinuit->Contour(100,0,2);
    gr->DrawClone("AL");*/

    // 1 sigma contour
    gMinuit->SetErrorDef(1);
    gr = (TGraph*) gMinuit->Contour(100,0,2);
    gr->DrawClone("AL");
    /*
    // contour of min. function
    c->cd(2);
    arglist[0] = 1; // WidthUp
    gMinuit->mnexcm("SCA", arglist ,1,ierflg);
    TGraph* gr = (TGraph*)gMinuit->GetPlot();
    gr->DrawClone("al");

    c->cd(4);
    arglist[0] = 4; // LengthLo
    gMinuit->mnexcm("SCA", arglist ,1,ierflg);
    gr = (TGraph*)gMinuit->GetPlot();
    gr->DrawClone("al");//*/

    c->cd(1);

APlot:
    //-----------------------------------------------------------------------
    // plot alpha histogram with current cuts and calc significance

    // alpha histogram
    TH1F *halpha=new TH1F("halpha","",18,0.,90.);
    TH1F *halphaOff=new TH1F("halphaOff","",18,0.,90.);

    halpha->SetStats(kFALSE);
    halphaOff->SetStats(kFALSE);

    for(Int_t i=0;i<nrowOn;i++)
    {
        if (mOn(i,isize)<SizeLo) continue;
        if((mOn(i,idist)<DistLo)||(mOn(i,idist)>DistUp)) continue;

        if((mOn(i,iwidth) < WidthLo) ||(mOn(i,iwidth) > WidthUp))  continue;
        if((mOn(i,ilength)< LengthLo)||(mOn(i,ilength)> LengthUp)) continue;

        halpha->Fill(TMath::Abs(mOn(i,ialpha)));
    }
    halpha->SetLineWidth(2);
    halpha->DrawClone();

    for(Int_t i=0;i<nrowOff;i++)
    {
        if (mOff(i,isize)<SizeLo) continue;
        if((mOff(i,idist)<DistLo)||(mOff(i,idist)>DistUp)) continue;

        if((mOff(i,iwidth) < WidthLo) ||(mOff(i,iwidth) > WidthUp))  continue;
        if((mOff(i,ilength)< LengthLo)||(mOff(i,ilength)> LengthUp)) continue;

        halphaOff->Fill(TMath::Abs(mOff(i,ialpha)));

    }

    // scale Off histogram to On data using bins 6-18
    Double_t nsample_on=0;
    Double_t nsample_off=0;
    for(Int_t i=6;i<18;i++)
        nsample_on+=halpha->GetBinContent(i);
    for(Int_t i=6;i<18;i++)
        nsample_off+=halphaOff->GetBinContent(i);

    Double_t scal=nsample_on/nsample_off;

    halphaOff->SetLineColor(2);
    halphaOff->Scale(scal);
    halphaOff->Draw("same");

    gPad->SetGridx();
    gPad->SetGridy();

    // define functions
    TF1 *backFcn = new TF1("backFcn",background,20.,89.5,4);
    backFcn->SetLineWidth(2);
    backFcn->SetLineColor(kRed);

    TF1 *signalFcn = new TF1("signalFcn",signal,0.,89.5,3);
    signalFcn->SetLineWidth(2);
    signalFcn->SetLineColor(kBlue);
    signalFcn->SetNpx(500);

    TF1 *fitFcn = new TF1("fitFcn",fitFunction,0.,89.5,7);
    fitFcn->SetNpx(500);
    fitFcn->SetLineWidth(2);
    fitFcn->SetLineColor(kMagenta);

    // start-parameters
    Double_t p3start=halpha->GetBinContent(halpha->GetNbinsX()-1);
    Double_t p4start=0.;
    Double_t p5start=0.;
    Double_t p6start=0.;

    Double_t p0start=(halpha->GetEntries()*halpha->GetBinWidth(1)-p3start)*0.05;
    Double_t p1start=5.;
    Double_t p2start=0.;

    cout<<"start values:"<<endl;

    cout<<" gaussian:"<<endl;
    cout<<"  p0 = "<<p0start<<endl;
    cout<<"  p1 = "<<p1start<<endl;
    cout<<"  p2 = "<<p2start<<endl;

    cout<<" pol3:"<<endl;
    cout<<"  p0 = "<<p3start<<endl;
    cout<<"  p1 = "<<p4start<<endl;
    cout<<"  p2 = "<<p5start<<endl;
    cout<<"  p3 = "<<p6start<<endl<<endl;


    // background fit
    Double_t parb[4];
    backFcn->SetParameters(p3start,p4start,p5start,p6start);
    backFcn->FixParameter(1,0);    // deriv. at zero equal zero
    //backFcn->FixParameter(2,0);  // no 2nd order term
    backFcn->FixParameter(3,0);    // no 3rd order term
    //halphaOff->Fit ("backFcn","RN");
    halpha->Fit ("backFcn","RN");
    backFcn->GetParameters(parb);

    // total fit
    Double_t par[7];
    fitFcn->SetParameters(p0start,p1start,p2start,
                          parb[0],parb[1],parb[2],parb[3]);

    fitFcn->FixParameter(2,0.);
    fitFcn->FixParameter(3,parb[0]);
    fitFcn->FixParameter(4,parb[1]);
    fitFcn->FixParameter(5,parb[2]);
    fitFcn->FixParameter(6,parb[3]);
    halpha->Fit ("fitFcn","RN");

    // draw the legend
    TLegend *legend=new TLegend(0.4,0.5,0.93,0.65);
    legend->SetTextFont(40);
    legend->SetTextSize(0.04);
    legend->SetFillColor(19);
    legend->AddEntry(halpha,"Data","lp");
    TString strpol("");
    strpol=Form("Polynom");
    legend->AddEntry(backFcn,strpol.Data(),"l");
    legend->AddEntry(fitFcn,"Global Fit","l");
    //legend->Draw();

    // drawing the functions
    fitFcn->GetParameters(par);
    fitFcn->SetRange(-90.,90.);
    fitFcn->SetLineWidth(2);
    fitFcn->Draw("same");

    signalFcn->SetParameters(par);
    signalFcn->SetLineWidth(2);
    //signalFcn->Draw("same");

    backFcn->SetParameters(&par[3]);
    backFcn->SetRange(-90,90);
    backFcn->SetLineWidth(2);
    backFcn->Draw("same");


    //-----------------------------------------------------------------------
    // calculate significance

    // analytical method
    Double_t nSig=signalFcn->Integral(0.,acut)/halpha->GetBinWidth(1);
    Double_t nOffScaled=backFcn->Integral(0.,acut)/halpha->GetBinWidth(1);

    Double_t nOn=nSig+nOffScaled;

    // significance after Li/Ma
    Double_t theta=1.; // should be = Ton/Toff ???
    nSig=nOn-nOffScaled;
    Double_t nOff=nOffScaled/theta;

    Double_t siglima = sqrt(2*(nOn*TMath::Log((1+theta)*nOn/(theta*(nOn+nOff)))
                               +nOff*TMath::Log((1+theta)*(nOff/(nOn+nOff)))));

    // "standard" significance
    Double_t sigstd  = nSig/sqrt(nOn+nOffScaled);


    TPaveText *pt = new TPaveText(0.25,0.75,0.56,0.98,"brNDC");
 
    pt->SetTextFont(60);
    pt->SetTextSize(0.03);
    pt->SetTextAlign(12);

    TString strcut1(Form("cuts:"));
    TString strcut3(Form(" |alpha| <  %4.2f ",acut));
    pt->AddText(strcut3.Data());
    TString str("");
    pt->AddText(str);

    TString strsig(Form("NExcess = %d",Int_t(nSig)));
    pt->AddText(strsig.Data());
    TString strbak(Form("NOff    = %d",Int_t(nOffScaled)));
    pt->AddText(strbak.Data());
    pt->AddText(str);

    TString strsg(Form("Significance (Li/Ma) = %4.2f",siglima));
    pt->AddText(strsg.Data());

    cout<<"results for events with alpha < "<<acut<<" degrees :"<<endl;
    cout<<" nbackground  = "<<nOffScaled<<endl;
    cout<<" theta        = "<<theta<<endl;
    cout<<" nsignal      = "<<nSig<<endl;
    cout<<" sig. (Li/Ma) = "<<siglima<<endl;
    cout<<" sig.-std     = "<<sigstd<<endl;
    pt->Draw();
    gPad->SetGridx();
    gPad->SetGridy();

    c->Update();

    cout<<endl<<"Optimized Cuts:"<<endl;
    cout<<" "<<DistLo<<" < dist <"<<DistUp<<endl;
    cout<<" "<<WidthLo<<" < width <"<<WidthUp<<endl;
    cout<<" "<<LengthLo<<" < length <"<<LengthUp<<endl;

    return;
}