//Double_t kRad2Deg = 180./TMath::Pi();

#include <float.h>
#include <iomanip.h>

#include <TH1.h>
#include <TH2.h>
#include <TLine.h>
#include <TLeaf.h>
#include <TChain.h>
#include <TGaxis.h>
#include <TStyle.h>
#include <TCanvas.h>
#include <TPolyLine.h>
#include <TPolyMarker.h>

#include "MH.h"
#include "MBinning.h"

#include "MFit.h"
#include "MPhoton.h"

TH1D Smear(const TH1 &h, Double_t sigma)
{
    TH1D ret;

    MH::SetBinning((TH1*)&ret, (TH1*)&h);

    Int_t n = h.GetNbinsX();

    for (int i=1; i<=n; i++)
    {
        Double_t xi = sqrt(h.GetBinLowEdge(i+1)*h.GetBinLowEdge(i));
        for (int j=1; j<=n; j++)
        {
            Double_t xj = sqrt(h.GetBinLowEdge(j+1)*h.GetBinLowEdge(j));

            Double_t dx = log10(xj/xi)/sigma;

            Double_t gaus = exp(-dx*dx/2)/(sigma*sqrt(TMath::Pi()*2));

            gaus *= h.GetBinContent(i);

            ret.Fill(xj, gaus/sqrt(n));
        }
    }
    return ret;
}

void GetRange(TChain *chain, const char *name, Int_t &nbins, Double_t &min, Double_t &max, Double_t conv=1, Bool_t zero=kTRUE)
{
    TString str("MPhoton.MParticle.");
    str += name;

    MPhoton *p = new MPhoton;

    chain->SetBranchAddress("MPhoton.", &p);
    chain->SetBranchStatus("*", 0);
    chain->SetBranchStatus(str, 1);

    min =  FLT_MAX;
    max = -FLT_MAX;

    Int_t  numt = chain->GetTreeNumber();
    TLeaf *leaf = chain->GetLeaf(str);

    if (!leaf && numt>=0)
        return;

    Stat_t n = chain->GetEntries();
    for (int i=0; i<n; i++)
    {
        chain->GetEntry(i);

        if (numt != chain->GetTreeNumber())
        {
            if (!(leaf = chain->GetLeaf(str)))
                return;

            numt = chain->GetTreeNumber();
        }

        Double_t val = leaf->GetValue();

        if (p->IsPrimary())
            continue;

        if (val < min && (zero || val!=0))
            min = val;
        if (val > max)
            max = val;
    }

    min *= conv;
    max *= conv;

    Int_t newn=0;
    MH::FindGoodLimits(nbins, newn, min, max, kFALSE);
    nbins = newn;
}

void draw1h1426(Double_t E, Double_t plot=1)
{
    plot += 1;

    const Double_t level = E / (5.73e-0*pow( 0.39e3, plot));

    TPolyLine   *l = new TPolyLine(6);
    TPolyMarker *m = new TPolyMarker(6);

    //  m->SetPoint(0, log10( 0.28e3), log10(1.86e-0*pow( 0.28e3, plot)*level));
    m->SetPoint(0, log10( 0.39e3)-0.3, log10(5.73e-0*pow( 0.39e3, plot)*level));
    m->SetPoint(1, log10( 0.80e3)-0.3, log10(1.22e-0*pow( 0.80e3, plot)*level));
    m->SetPoint(2, log10( 1.55e3)-0.3, log10(5.10e-2*pow( 1.55e3, plot)*level));
    m->SetPoint(3, log10( 2.82e3)-0.3, log10(1.50e-2*pow( 2.82e3, plot)*level));
    m->SetPoint(4, log10( 5.33e3)-0.3, log10(7.70e-3*pow( 5.33e3, plot)*level));
    m->SetPoint(5, log10(10.00e3)-0.3, log10(1.20e-4*pow(10.00e3, plot)*level));
    m->SetMarkerStyle(kOpenStar);
    m->SetMarkerColor(kMagenta);
    m->SetBit(kCanDelete);
    m->Draw();

    //  l->SetPoint(0, log10( 0.28e3), log10(1.86e-0*pow( 0.28e3, plot)*level));
    l->SetPoint(0, log10( 0.39e3)-0.3, log10(5.73e-0*pow( 0.39e3, plot)*level));
    l->SetPoint(1, log10( 0.80e3)-0.3, log10(1.22e-0*pow( 0.80e3, plot)*level));
    l->SetPoint(2, log10( 1.55e3)-0.3, log10(5.10e-2*pow( 1.55e3, plot)*level));
    l->SetPoint(3, log10( 2.82e3)-0.3, log10(1.50e-2*pow( 2.82e3, plot)*level));
    l->SetPoint(4, log10( 5.33e3)-0.3, log10(7.70e-3*pow( 5.33e3, plot)*level));
    l->SetPoint(5, log10(10.00e3)-0.3, log10(1.20e-4*pow(10.00e3, plot)*level));
    l->SetLineColor(kMagenta);
    l->SetBit(kCanDelete);
    l->Draw();
}

void analp()
{
    // -------------------------------------------------------------------

    TString dir = "~/data/";

    TChain chain("Photons");
    /*
     chain.Add(dir+"cascade_0.1_18_1e2_1e5_B0_256_01.root");
     chain.Add(dir+"cascade_0.1_18_1e2_1e5_B0_256_02.root");
     chain.Add(dir+"cascade_0.1_18_1e2_1e5_B0_256_03.root");
     chain.Add(dir+"cascade_0.1_18_1e2_1e5_B0_256_04.root");

     chain.Add(dir+"cascade_0.1_18_1e2_1e5_B0_512_01.root");
     chain.Add(dir+"cascade_0.1_18_1e2_1e5_B0_512_02.root");
     chain.Add(dir+"cascade_0.1_18_1e2_1e5_B0_512_03.root");
     chain.Add(dir+"cascade_0.1_18_1e2_1e5_B0_512_05.root");
     */

    //chain.Add(dir+"cascade_0.5_18_1e2_1e5_B0_512_01.root");
    //chain.Add(dir+"cascade_0.5_18_1e2_1e5_B1e-7_10Mpc_512_01.root");
    //chain.Add(dir+"cascade_0.5_18_1e2_1e5_B1e-7_50Mpc_512_01.root");
    chain.Add(dir+"cascade_0.5_18_1e2_1e5_B1e-7_100Mpc_512_01.root");

    //chain.Add(dir+"cascade_0.1_18_1e2_1e5_B1e-6_10Mpc_512_01.root");
    //chain.Add(dir+"cascade_0.1_18_1e2_1e5_B1e-6_50Mpc_512_01.root");
    //chain.Add(dir+"cascade_0.1_18_1e2_1e5_B1e-7_100Mpc_512_01.root");

    /*
     chain.Add(dir+"cascade_0.03_18_1e2_1e5_B0_512_01.root");
     chain.Add(dir+"cascade_0.03_18_1e2_1e5_B0_512_02.root");
     chain.Add(dir+"cascade_0.03_18_1e2_1e5_B0_512_03.root");
     chain.Add(dir+"cascade_0.03_18_1e2_1e5_B0_512_04.root");
     */

    //chain.Add(dir+"cascade_0.03_18_1e2_1e5_B1e-6_10Mpc_512_01.root");
    //chain.Add(dir+"cascade_0.03_18_1e2_1e5_B1e-6_50Mpc_512_01.root");
    //chain.Add(dir+"cascade_0.03_18_1e2_1e5_B1e-7_100Mpc_512_01.root");

    Int_t nbinse = 18;      // number of bins in your histogram
    Double_t lo  = 1e2/1e3; // lower limit of the energy histogram
    Double_t hi  = 1e5;     // upper limit of the energy histogram

    Double_t minE = 1e2;    // minimum value produced in the simulation

    Double_t cutoff = 1e12; // arbitrary energy cutoff (throw away the events)

    Double_t alpha = -1;    // alpha of the energy spectrum (-1 is E^-2)
    Double_t plot  =  1;    // plot of the enrgy spectrum (1 is E^2)

    Int_t nbins = 16;       // number of bins you want to have in psi/phi

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

    MBinning binspolx;
    binspolx.SetEdges(16, -180, 180);

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

    Stat_t n = chain.GetEntries();

    cout << endl << "Found " << n << " entries." << endl;

    if (n==0)
        return;

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

    MBinning bins;

    cout << "Determin min/max..." << endl;

    Double_t min, max;
    GetRange(&chain, "fTheta", nbins, min, max, kRad2Deg*60);
    cout << "fTheta, " << nbins << ": " << min << " ... " << max << " [']" << endl;

    bins.SetEdges(nbins, min, max*1.1);

    TH2D a, a2, a3;
    MH::SetBinning(&a,  &binspolx, &bins);
    MH::SetBinning(&a2, &binspolx, &bins);
    MH::SetBinning(&a3, &binspolx, &bins);

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

    GetRange(&chain, "fR", nbins, min, max);
    cout << "fR,     " << nbins << ": " << min << " ... " << max << " [kpc]" << endl;

    bins.SetEdges(nbins, min, max*1.1);

    TH2D r, r2, r3;
    MH::SetBinning(&r,  &binspolx, &bins);
    MH::SetBinning(&r2, &binspolx, &bins);
    MH::SetBinning(&r3, &binspolx, &bins);

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

    MBinning binsx;
    binsx.SetEdgesLog(nbinse, lo, hi);

    TH1D prim, h, h2, h3;
    MH::SetBinning(&h,    &binsx);
    MH::SetBinning(&h2,   &binsx);
    MH::SetBinning(&h3,   &binsx);
    MH::SetBinning(&prim, &binsx);
//    prim.Sumw2();
    h.Sumw2();
    h2.Sumw2();
    h3.Sumw2();

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

    MPhoton *p = new MPhoton;
    chain.SetBranchAddress("MPhoton.", &p);
    chain.SetBranchStatus("*", 1);

    chain.GetEntry(0);
    if (!p->IsPrimary())
        return;

    Double_t z = p->GetZ();
    Double_t R = MParticle::RofZ(&z);

    cout << "Z = " << z << endl;
    cout << "R = " << R << "kpc" << endl;

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

    const Double_t skpc = 3600.*24.*365.*3.258*1e3;  // [s/kpc]
    GetRange(&chain, "fX", nbins, min, max, 1, kFALSE);
    min *= skpc;
    max *= skpc;
    if (min<1e-2)
        min=1e-2;
    bins.SetEdgesLog(nbins, min, max);
    cout << "dX,     " << nbins << ": " << min << " ... " << max << " [s]" << endl;

    TH1D t;
    MH::SetBinning(&t, &bins);
    t.Sumw2();

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

    chain.SetBranchAddress("MPhoton.", &p);
    chain.SetBranchStatus("*", 1);

    cout << "Filling: " << nbinse << " bins @ " << minE << " - " << hi << "... " << flush;

    Double_t weight = 0;

    Bool_t isprim = kFALSE;
    Bool_t ignore = kFALSE;
    for (int i=0; i<n; i++)
    {
        chain.GetEntry(i);

        Double_t Ep = p->GetEnergy();

        if (p->IsPrimary())
        {
            if (Ep>cutoff)
            {
                ignore = kTRUE;
                continue;
            }
            ignore = kFALSE;
            weight = pow(Ep, alpha);
            prim.Fill(Ep, pow(Ep, plot+alpha));
            isprim = kTRUE;
            continue;
        }

        if (ignore)
            continue;

        Double_t phi = fmod(p->GetPhi()*kRad2Deg+180, 360)-180;
        Double_t psi = fmod(p->GetPsi()*kRad2Deg+180, 360)-180;

        if (isprim)
        {
            h3.Fill(Ep, pow(Ep,plot) * weight);
            r3.Fill(phi, 0.0, weight);
            a3.Fill(psi, 0.0, weight);
            isprim = kFALSE;
            continue;
        }

        h2.Fill(Ep, pow(Ep,plot) * weight);
        r2.Fill(phi, p->GetR(), weight);
        a2.Fill(psi, p->GetTheta()*kRad2Deg*60, weight);

        Double_t tm1 = sqrt(R*R+p->GetR()*p->GetR());
        Double_t tm2 = (R+p->GetX())/tm1 - 1;
        t.Fill(p->GetX()<=0 ? 0 : R*tm2*skpc, weight);
    }

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

    delete p;

    cout << "Done." << endl;

    Double_t N = prim.GetBinContent(nbinse);
    Double_t Nerr = prim.GetBinError(nbinse);
    for (int i=1; i<=nbinse; i++)
    {
        Double_t E = prim.GetXaxis()->GetBinCenter(i);
        if (E>minE)
            continue;

        Double_t E1 = prim.GetXaxis()->GetBinLowEdge(i+1);
        Double_t E0 = prim.GetXaxis()->GetBinLowEdge(i);

        E = sqrt(E1*E0);

        prim.SetBinContent(i, N);
        h3.SetBinContent(i, N);
        prim.SetBinError(i, Nerr);
        h3.SetBinError(i, Nerr);
    }

    h.Add(&h2, &h3);
    r.Add(&r2, &r3);
    a.Add(&a2, &a3);

    if (alpha==-plot)
    {
        cout << "Observed Energy: " << h.Integral() << "GeV" << endl;
        cout << "Emitted Energy:  " << prim.Integral() << "GeV" << endl;
        cout << "Ratio: " << h.Integral()/prim.Integral() << endl;
    }

    cout << "Mean fPhi: " << r.GetMean() << " +- " << r.GetRMS() << endl;
    cout << "Mean fPsi: " << a.GetMean() << " +- " << a.GetRMS() <<  endl;

    gStyle->SetOptStat(10);
    gStyle->SetPalette(1, 0);

    TLine line;

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

//    delete gROOT->FindObject("Analysis Arrival");
    TCanvas *c = MH::MakeDefCanvas("Analysis Arrival", "");
    c->Divide(2,2);

    c->cd(1);
    r.SetTitle(" Distance from Observer ");
    r.GetXaxis()->SetLabelOffset(-0.015);
    r.GetXaxis()->SetTitleOffset(1.1);
    r.GetXaxis()->SetRangeUser(1e4, 1e9);
    r.SetXTitle("\\Phi [\\circ]");
    r.SetYTitle("R [kpc]");
    TVirtualPad &pad = *gPad;
    pad.Divide(2,2);
    pad.cd(1);
    gPad->SetLogy();
    gPad->SetLogz();
    gPad->SetTheta(0);
    gPad->SetPhi(90);
    TH1 *g = r.DrawCopy("surf1pol");
    pad.cd(2);
    gPad->SetLogy();
    gPad->SetLogz();
    gPad->SetTheta(70);
    gPad->SetPhi(90);
    g->Draw("surf1pol");
    pad.cd(3);
    gPad->SetLogy();
    gPad->SetLogz();
    gPad->SetTheta(90);
    gPad->SetPhi(90);
    g->Draw("surf1pol");
    pad.cd(4);
    gPad->SetLogy();
    gPad->SetLogz();
    gPad->SetTheta(20);
    gPad->SetPhi(90);
    g->Draw("surf1pol");

    c->cd(2);
    a.SetTitle(" Hit Angle for Observer ");
    a.GetXaxis()->SetLabelOffset(-0.015);
    a.GetXaxis()->SetTitleOffset(1.1);
    a.GetXaxis()->SetRangeUser(1e4, 1e9);
    a.SetXTitle("\\Psi [\\circ]");
    a.SetYTitle("\\Theta [']");
    TVirtualPad &pad2 = *gPad;
    pad2.Divide(2,2);
    pad2.cd(1);
    gPad->SetLogy();
    gPad->SetLogz();
    gPad->SetTheta(0);
    gPad->SetPhi(90);
    g = a.DrawCopy("surf1pol");
    pad2.cd(2);
    gPad->SetLogy();
    gPad->SetLogz();
    gPad->SetTheta(70);
    gPad->SetPhi(90);
    g->Draw("surf1pol");
    pad2.cd(3);
    gPad->SetLogy();
    gPad->SetLogz();
    gPad->SetTheta(90);
    gPad->SetPhi(90);
    g->Draw("surf1pol");
    pad2.cd(4);
    gPad->SetLogy();
    gPad->SetLogz();
    gPad->SetTheta(20);
    gPad->SetPhi(90);
    g->Draw("surf1pol");

    c->cd(3);
    gPad->SetLogy();
    g=r.ProjectionY();
    g->SetDirectory(NULL);
    TH1 *g2=r2.ProjectionY();
    g->SetMinimum(MH::GetMinimumGT(*g2)/2);
    g->SetBit(kCanDelete);
    g->SetTitle(" Hit Observers Plain ");
    g->GetXaxis()->SetLabelOffset(0);
    g->GetXaxis()->SetTitleOffset(1.1);
    g->GetYaxis()->SetTitleOffset(1.3);
    g->SetXTitle("R [kpc]");
    g->SetYTitle("Counts");
    g->Draw();
    g2->SetDirectory(NULL);
    g2->SetBit(kCanDelete);
    g2->SetLineColor(kBlue);
    g2->Draw("same");
    g=r3.ProjectionY();
    g->SetDirectory(NULL);
    g->SetBit(kCanDelete);
    g->SetLineColor(kGreen);
    g->Draw("same");

    c->cd(4);
    gPad->SetLogy();
    g=a.ProjectionY();
    g->SetMinimum(MH::GetMinimumGT(*g)/2);
    g->SetDirectory(NULL);
    g->SetBit(kCanDelete);
    g->SetTitle("Hit Angle");
    g->GetXaxis()->SetLabelOffset(0);
    g->GetXaxis()->SetTitleOffset(1.1);
    g->GetYaxis()->SetTitleOffset(1.3);
    g->SetXTitle("\\Phi [']");
    g->SetYTitle("Counts");
    g->Draw();
    g=a2.ProjectionY();
    g->SetDirectory(NULL);
    g->SetBit(kCanDelete);
    g->SetLineColor(kBlue);
    g->Draw("same");
    g=a3.ProjectionY();
    g->SetDirectory(NULL);
    g->SetBit(kCanDelete);
    g->SetLineColor(kGreen);
    g->Draw("same");

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

    //  delete gROOT->FindObject("Analysis Photons");
    c = MH::MakeDefCanvas("Analysis Photons", "", 580, 870);
    c->Divide(1,2);

    c->cd(1);
    gPad->SetLogx();
    gPad->SetLogy();
    h3.SetTitle(Form(" E^{%.1f}  Spectra @ z=%.1e (R=%.0fkpc) ", alpha, z, R));
    h3.SetXTitle("E\\gamma [GeV]");
    h3.SetYTitle(Form("E^{%.1f} * Counts", plot));
    h3.GetXaxis()->SetLabelOffset(-0.015);
    h3.GetXaxis()->SetTitleOffset(1.1);
    h3.SetFillStyle(0);
    h3.SetName("PrimPhotons");
    h3.SetMarkerStyle(kFullCircle);
    h3.SetMarkerColor(kRed);
    h3.SetMarkerSize(0.8);
    h3.SetLineColor(kRed);
    h3.DrawCopy("C"); //E1

    TH1D hs = Smear(h, 0.2);
    hs.SetLineColor(kGreen);
    hs.SetFillStyle(0);
    hs.DrawCopy("Csame");

    h.SetFillStyle(0);
    h.SetName("Photons");
    h.SetMarkerStyle(kFullCircle);
    h.SetMarkerSize(0.8);
    h.DrawCopy("Csame E4");
    prim.SetFillStyle(0);
    prim.SetMarkerStyle(kPlus);
    prim.SetMarkerColor(kRed);
    prim.SetLineColor(kRed);
    prim.SetMarkerSize(0.8);
    prim.DrawCopy("Csame");
    h2.SetFillStyle(0);
    h2.SetName("SecPhotons");
    h2.SetMarkerStyle(kFullCircle);
    h2.SetMarkerColor(kBlue);
    h2.SetMarkerSize(0.8);
    h2.SetLineColor(kBlue);
    h2.DrawCopy("Csame E4"); //E2

    draw1h1426(prim.GetBinContent(2),plot);

    c->cd(2);
    gPad->SetLogx();
    TH1D div;
    MH::SetBinning(&div, &binsx);
    div.Sumw2();
    div.Divide(&h, &prim);
    div.SetTitle(" Spectrum / Primary Spectrum ");
    div.GetXaxis()->SetLabelOffset(-0.015);
    div.GetXaxis()->SetTitleOffset(1.1);
    div.SetXTitle("E\\gamma [GeV]");
    div.SetYTitle("Ratio");
    div.SetMarkerStyle(kPlus);
    TH1 *gHist = div.DrawCopy("E4");

    line.SetLineWidth(1);
    line.SetLineColor(kGreen);
    line.DrawLine(log10(lo), exp(-1), log10(hi), exp(-1));
    line.SetLineColor(kBlue);
    line.DrawLine(log10(lo), 1, log10(hi), 1);

    MFit fit("exp(-1)-[1]*log10(x/[0])");
    for (int i=0; i<gHist->GetEntries(); i++)
    {
        if (gHist->GetBinContent(i+1)<exp(-1))
        {
            fit.SetRange(gHist->GetBinCenter(i-2), gHist->GetBinCenter(i+2));
            cout << "Fitting from " << gHist->GetBinLowEdge(i-1) << " to ";
            cout << gHist->GetBinLowEdge(i+1) << endl;
            break;
        }
    }
    fit.SetParameter(0, "t", 750,   10, 1e5);
    fit.SetParameter(1, "m", 0.5,  0.1,  10);
    fit.FitLog(gHist);
    fit.Print();
    fit.DrawCopy("same");

    cout << "Cutoff:  " << setprecision(2) << fit[0]/1e3 << " +- ";
    cout << fit.GetParError(0)/1e3 << " TeV" << endl;

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

    c = MH::MakeDefCanvas("Time Analysis", "", 580, 435);
    gPad->SetLogx();
    //    gPad->SetLogy();
    t.SetName("Times");
    t.SetTitle(" Arrival time distribution ");
    t.SetXTitle("\\Delta  t [s]");
    t.GetXaxis()->SetLabelOffset(-0.015);
    t.GetXaxis()->SetTitleOffset(1.1);
    t.DrawCopy("E4");
    line.DrawLine(log10(1),           0, log10(1),           t.GetMaximum()*1.05);
    line.DrawLine(log10(3600),        0, log10(3600),        t.GetMaximum()*1.05);
    line.DrawLine(log10(3600*24),     0, log10(3600*24),     t.GetMaximum()*1.05);
    line.DrawLine(log10(3600*24*365), 0, log10(3600*24*365), t.GetMaximum()*1.05);
}