#include <math.h>
#include <iostream.h>
#include <fstream.h>

#include <TStopwatch.h>
#include <TF1.h>
#include <TH1.h>
#include <TRandom.h>
#include <TGraph.h>
#include <TCanvas.h>
#include <TStyle.h>

#include "mbase/MParContainer.h"
#include "mphys/MPhoton.h"
#include "mphys/MElectron.h"
#include "mphys/MPairProduction.h"
#include "mhist/MBinning.h"
#include "mhist/MH.h"

// 2.96
// 2.87
// 2.73

// ================================================================
Double_t DiSum(Double_t *x, Double_t *k=NULL)
{
    Double_t t = x[0];

    Double_t disum = t;
    Double_t add = 0;
                   
    Double_t eps = fabs(t*1e-1);

    Int_t    n   = 2;
    Double_t pow = t*t;   // t^2

    do
    {
        add = pow/n/n;

        pow *= t;       // pow = t^n
        n++;

        disum += add;

    } while (fabs(add)>eps);

    return disum;
}

Double_t F(Double_t *x, Double_t *k=NULL)
{
    Double_t o = x[0];
    Double_t s = -2.*o;

//    if (o<1e-10)
//        return 2.125; //-3./8.;

    return -o/4. + (9./4. + 1./o + o/2.) * log(1.+2.*o) + 1./8.*(1.+2.*o) + MElectron::Li2(&s); //- 3./8.
}

void rkck(Double_t y[], Double_t dydx[], int n, Double_t x, Double_t h, Double_t yout[],
          Double_t yerr[], void (*derivs)(Double_t, Double_t[], Double_t[]))
{
    /*
     * ---------------------------------------------------------
     * Numerical recipes for C, Chapter 16.1, Runge-Kutta Method
     * ---------------------------------------------------------
     */
    const Double_t a2  = 0.2;
    const Double_t a3  = 0.3;
    const Double_t a4  = 0.6;
    const Double_t a5  = 1.0;
    const Double_t a6  = 0.875;
    const Double_t b21 = 0.2;
    const Double_t b31 = 3.0/40.0;
    const Double_t b32 = 9.0/40.0;
    const Double_t b41 = 0.3;
    const Double_t b42 = -0.9;
    const Double_t b43 = 1.2;
    const Double_t b51 = -11.0/54.0;
    const Double_t b52 = 2.5;
    const Double_t b53 = -70.0/27.0;
    const Double_t b54 = 35.0/27.0;
    const Double_t b61 = 1631.0/55296.0;
    const Double_t b62 = 175.0/512.0;
    const Double_t b63 = 575.0/13824.0;
    const Double_t b64 = 44275.0/110592.0;
    const Double_t b65 = 253.0/4096.0;
    const Double_t c1  = 37.0/378.0;
    const Double_t c3  = 250.0/621.0;
    const Double_t c4  = 125.0/594.0;
    const Double_t c6  = 512.0/1771.0;
    const Double_t dc5 = -277.00/14336.0;

    const Double_t dc1 = c1-2825.0/27648.0;
    const Double_t dc3 = c3-18575.0/48384.0;
    const Double_t dc4 = c4-13525.0/55296.0;
    const Double_t dc6 = c6-0.25;

    Double_t ak2[n];
    Double_t ak3[n];
    Double_t ak4[n];
    Double_t ak5[n];
    Double_t ak6[n];
    Double_t ytemp[n];

    for (int i=0; i<n; i++)
        ytemp[i] = y[i]+b21*h*dydx[i];

    (*derivs)(x+a2*h,ytemp,ak2);

    for (int i=0; i<n; i++)
        ytemp[i] = y[i]+h*(b31*dydx[i]+b32*ak2[i]);

    (*derivs)(x+a3*h,ytemp,ak3);

    for (int i=0; i<n; i++)
        ytemp[i] = y[i]+h*(b41*dydx[i]+b42*ak2[i]+b43*ak3[i]);

    (*derivs)(x+a4*h,ytemp,ak4);

    for (int i=0; i<n; i++)
        ytemp[i] = y[i]+h*(b51*dydx[i]+b52*ak2[i]+b53*ak3[i]+b54*ak4[i]);

    (*derivs)(x+a5*h,ytemp,ak5);

    for (int i=0; i<n; i++)
        ytemp[i] = y[i]+h*(b61*dydx[i]+b62*ak2[i]+b63*ak3[i]+b64*ak4[i]+b65*ak5[i]);

    (*derivs)(x+a6*h,ytemp,ak6);

    for (int i=0; i<n; i++)
        yout[i]=y[i]+h*(c1*dydx[i]+c3*ak3[i]+c4*ak4[i]+c6*ak6[i]);

    for (int i=0; i<n; i++)
        yerr[i]=h*(dc1*dydx[i]+dc3*ak3[i]+dc4*ak4[i]+dc5*ak5[i]+dc6*ak6[i]);
}

Double_t FMIN(Double_t a, Double_t b)
{
    return a<b ? a : b;
}

Double_t FMAX(Double_t a, Double_t b)
{
    return a>b ? a : b;
}

void SolvEq(Double_t y[], Double_t dydx[], int n, Double_t *x, Double_t htry, Double_t eps,
            Double_t yscal[], Double_t *hdid, Double_t *hnext,
            void (*derivs)(Double_t, Double_t[], Double_t[])) // rkqs
{
    /*
     * ---------------------------------------------------------
     * Numerical recipes for C, Chapter 16.1, Runge-Kutta Method
     * ---------------------------------------------------------
     */
    const Double_t SAFETY = 0.9;
    const Double_t PGROW  = -0.2;
    const Double_t PSHRNK = -0.25;
    const Double_t ERRCON = 1.89e-4;

    Double_t yerr[n];
    Double_t ytemp[n];

    Double_t h = htry;
    Double_t errmax;
    while (1)
    {
        rkck(y, dydx, n, *x, h, ytemp, yerr, derivs);

        errmax=0.0;

        for (int i=0; i<n; i++)
            errmax = FMAX(errmax, fabs(yerr[i]/yscal[i]) );

        errmax /= eps;

        if (errmax <= 1.0)
            break;

        Double_t htemp = SAFETY*h*pow(errmax,PSHRNK);

        h = (h >= 0.0 ? FMAX(htemp,0.1*h) : FMIN(htemp,0.1*h));

        Double_t xnew= (*x) + h;

        if (xnew != *x)
            continue;

        cout << "stepsize underflow in rkqs" << endl;
        break;
    }

    *hnext = errmax>ERRCON ? SAFETY*h*pow(errmax,PGROW) : 5.0*h;

    *x += (*hdid=h);

    for (int i=0; i<n; i++)
        y[i] = ytemp[i];
}

Double_t dEdt(Double_t E, Double_t t, Double_t z0=-1)
{
    /* ------------------------------------
     * Lower limit as differential Equation
     * ------------------------------------

     Double_t T     = 2.96;                  // [K]
     Double_t sigma = 6.653e-29;             // [m^2]
     Double_t E0    = 511e-6;                // [GeV]
     Double_t e     = 1.602176462e-19;       // [C]
     Double_t kB    = 1e-9/e*1.3806503e-23;  // [GeV/K]
     Double_t h     = 1e-9/e*6.62606876e-34; // [GeVs]
     Double_t hc    = h*c;                   // [GeVm]
     Double_t pi    = TMath::Pi();           // [1]
     Double_t khc   = pi*kB/hc;              // [1 / K m]
     Double_t a     = 8./15 * pi * khc*khc*khc*khc * hc;         // [Gev / K^4 / m^3]

     Double_t konst = 4./3. * sigma * a *T*T*T*T *c /E0 /E0;

     return -konst *E*E;
     */
    Double_t pc = 1./3.258;                // [pc/ly]
    Double_t c  = 299792458;               // [m/s]
    Double_t ly = 3600.*24.*365.*c;        // [m/ly]
    Double_t r  = t * c / ly * pc / 1000;  // [kpc]
    Double_t R  = MParticle::RofZ(&z0) - r;
    Double_t z  = z0>=0 ? MParticle::ZofR(&R) : 0;

    Double_t e     = 1.602176462e-19;       // [C]
    Double_t T     = 2.96*(z+1);            // [K]
    Double_t kB    = 1e-9/e*1.3806503e-23;  // [GeV/K]
    Double_t kBT   = kB*T;                  // [GeV]
    Double_t alpha = 1./137;                // [|]
    Double_t h     = 1e-9/e*6.62606876e-34; // [GeVs]
    Double_t E0    = 511e-6;                // [GeV]


    Double_t k = TMath::Pi() * alpha * kBT;

    Double_t ln = 4.*kBT/E0/E0;

    return -1./3/h* k*k * (log(ln*E)-1.9805);
}

void dEdt(Double_t t, Double_t E[], Double_t dedt[])
{
    dedt[0] = dEdt(E[0], t);
    //cout << t << "\t" << E[0] << "\t" << dedt[0] << endl;
}

void DrawDevelopmentHiLim(Double_t E0, Double_t z, Option_t *opt="")
{
    Double_t t = 0;
    Double_t E[1] = { E0 };
    Double_t yscal[1] = { 1 };

    Double_t dedt[1];

    Double_t eps;// = 1e5;
    Double_t step = 5e6;

    Double_t hdid;
    Double_t hnext;

    cout << "Start: " << dedt[0] << endl;

    Int_t n = 15000;
    Double_t tres[n];
    Double_t Eres[n];
    int i;
    for (i=0; i<n; i++)
    {
        tres[i] = t;

        eps = E[0]/1e9; //9;

        dedt[0] = dEdt(E[0], t);

        SolvEq(E, dedt, 1, &t, step, eps, yscal, &hdid, &hnext, dEdt);

        step = hnext;

        Eres[i] = E[0];

        if (i==0) cout << "Did: " << hdid << endl;

        if (t>1.5e14 || E[0]<5e6)
            break;
//        cout << tres[i] << "\t" << Eres[i] << "\t(" << step << ")" << endl;
    }

    cout << i << endl;

    TGraph grp(i<n?i:n, tres, Eres);
    grp.Clone()->Draw(opt);

}

Double_t EnergyLossRateLoLim(Double_t *x, Double_t *k=NULL)
{
    Double_t t  = x[0];
    Double_t E  = k[0];
    Double_t t0 = k[1];

    Double_t c   = 299792458;               // [m/s]
    Double_t ly  = 3600.*24.*365.*c;        // [m/ly]
    Double_t pc  = 1./3.258;                // [pc/ly]

    Double_t r   = t * c / ly * pc / 1000;  // [kpc]
    Double_t R   = MParticle::RofZ(&k[2]) - r;
    Double_t z   = k[2]>=0 ? MParticle::ZofR(&R) : 0;

    Double_t T     = 2.96*(z+1);            // [K]
    //    Double_t alpha = 1./137;                // [1]
    Double_t sigma = 6.653e-29;             // [m^2]
    Double_t E0    = 511e-6;                // [GeV]
    Double_t e     = 1.602176462e-19;       // [C]
    Double_t kB    = 1e-9/e*1.3806503e-23;  // [GeV/K]
    Double_t h     = 1e-9/e*6.62606876e-34; // [GeVs]
    Double_t hc    = h*c;                   // [GeVm]
    Double_t pi    = TMath::Pi();           // [1]
    Double_t khc   = pi*kB/hc;              // [1 / K m]
    Double_t a     = 8./15 * pi * khc*khc*khc*khc * hc;         // [Gev / K^4 / m^3]
    Double_t konst = 4./3 * sigma * a * T*T*T*T * c / (E0* E0); // [1 / GeV s]

    Double_t ret = 1./(konst*(t-t0) + 1./E);

    return ret;
}

void DrawDevelopmentLoLim(Double_t t0, Double_t E0, Double_t z=-1, Option_t *opt="")
{
    // 8.7

    Double_t val[] = { E0, t0, z };
    Double_t t = 1.5e14;
    while (EnergyLossRateLoLim(&t, val)<1e4)
        t -= 0.01e14;

    TF1 *f1=new TF1("LoLim", EnergyLossRateLoLim, t, 1.5e14, 2);
    f1->SetParameter(0, E0);
    f1->SetParameter(1, t0);

    f1->Draw(opt);
    f1->SetBit(kCanDelete);
}

//
// (3) Energy loss rate of electrons and 'high energy particle'
//
Double_t DrawDevelopment(Double_t E, Double_t z, Option_t *opt="", TH1 *hist=NULL)
{
    Double_t ly = 3600.*24.*365.; // [s/ly]
    Double_t pc = 1./3.258;       // [pc/ly]

    TGraph *graph = new TGraph;
    graph->SetPoint(0, 0, E);
    graph->SetMaximum(E*3); // *MENU*

    cout << "------ " << endl;

    static TRandom rand;

    Double_t x = 0;
    for (int i=1; i<10; i++)
    {
        Double_t l = rand.Exp(MElectron::InteractionLength(&E, &z));

        if (z>=0)
        {
            Double_t r = MParticle::RofZ(&z) - l;

            cout << " " << i << ". R=" << MParticle::RofZ(&z) << " l=" << l << " z=" << MParticle::ZofR(&r) << endl;

            z = r>=0 ? MParticle::ZofR(&r) : 0;

            if (z==0)
                cout << "z<0" << endl;
        }

        x += l;

        Double_t t = x/pc*ly*1000;
        graph->SetPoint(i*2-1, t, E);

        Double_t e1 = MElectron::GetEnergyLoss(E, z<0?0:z);

        E -= e1;

        if (hist)
            hist->Fill(e1);

        cout << " Ep=" << e1 << flush;

        graph->SetPoint(i*2, t, E);
    }

    graph->SetMinimum(E/3); // *MENU*
    graph->Draw(opt);
    graph->SetBit(kCanDelete);

    //if (E<31500)
    cout << "t=" << x*ly/pc*1000 << "\tE=" << E << "\tz=" << z << endl;

    return E;
}

void EnergyLossRate()
{
    if (gPad)
        gPad->Clear();

    Double_t E = 1.5e9; //  [GeV]
    Double_t z = 0.03;

    MBinning bins;
    bins.SetEdgesLog(18, 0.1, 1e9);

    TH1D hist;
    hist.SetName("Phot");
    hist.SetTitle("Photons from inverse Compton");
    MH::SetBinning(&hist, &bins);

    cout << "Working..." << flush;

    for (int i=0; i<50; i++)
    {
        cout << i << "." << flush;
        DrawDevelopment(E, z, i?"L":"AL", &hist);
    }

    //DrawDevelopmentLoLim(2e14, 1.64e2, "Lsame"); // seen
    DrawDevelopmentLoLim(1.78e14, 280, z, "Lsame"); // seen
    DrawDevelopmentHiLim(E, z, "L");
    gPad->SetLogy();

    new TCanvas("Photons", "Photons created in inverse Compton");
    hist.DrawCopy();
    gPad->SetLogx();

    cout << "...done." << endl;
}

void DrawRZ()
{
    new TCanvas("RZ", "r and z");

    TF1 f1("ZofR", MParticle::ZofR, 0, 7.1e6, 0);
    TF1 f2("RofZ", MParticle::RofZ, 0, 5,     0);

    gPad->Divide(2,2);

    gPad->GetVirtCanvas()->cd(1);
    TH1 *h = f1.DrawCopy()->GetHistogram();

    h->SetTitle("z(r)");
    h->SetXTitle("r [kpc]");
    h->SetYTitle("z");

    gPad->Modified();
    gPad->Update();

    gPad->GetVirtCanvas()->cd(2);
    h = f2.DrawCopy()->GetHistogram();

    h->SetTitle("r(z)");
    h->SetXTitle("z");
    h->SetYTitle("r [kpc]");

    gPad->Modified();
    gPad->Update();
}

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

Double_t func(Double_t *x, Double_t *k)
{
    return MPhoton::Int2(x, k)*1e68;
}

void energyloss()
{
/*    Double_t E0    = 511e-6;                // [GeV]

    Double_t Eg = 1e4; //3.6e4;
    Double_t z  = 5;

    Double_t val[2] = { Eg, z };

    Double_t lolim = E0*E0/Eg;
    Double_t inf   = Eg<1e6 ? 3e-11*z : 3e-12*z;

    cout << Eg << " " << z << " " << lolim << " " << inf << endl;

    TF1 f("int2", func, lolim, inf, 2);

    Double_t int2 = f.Integral(lolim, inf, val); //[GeV^3 m^2]

    f.SetParameter(0, Eg);
    f.SetParameter(1, z);

    cout << int2 << endl;

    new TCanvas("ILPhoton", "Mean Interaction Length Photon");

    gPad->SetLogx();
//    gPad->SetLogy();
    gPad->SetGrid();

    f.SetLineWidth(1);
    f.DrawCopy();

    gPad->Modified();
    gPad->Update();

    //    return;
  */  MPhoton p;
    p.DrawInteractionLength();

    return;
//    EnergyLossRate();

    DrawRZ();

    return;

}