/* ======================================================================== *\
!
! *
! * This file is part of MARS, the MAGIC Analysis and Reconstruction
! * Software. It is distributed to you in the hope that it can be a useful
! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
! * It is distributed WITHOUT ANY WARRANTY.
! *
! * Permission to use, copy, modify and distribute this software and its
! * documentation for any purpose is hereby granted without fee,
! * provided that the above copyright notice appear in all copies and
! * that both that copyright notice and this permission notice appear
! * in supporting documentation. It is provided "as is" without express
! * or implied warranty.
! *
!
!
!   Author(s): Harald Kornmayer 1/2001 (harald@mppmu.mpg.de)
!   Author(s): Thomas Bretz  12/2000 (tbretz@uni-sw.gwdg.de)
!
!   Copyright: MAGIC Software Development, 2000-2001
!
!
\* ======================================================================== */

//////////////////////////////////////////////////////////////////////////////
//                                                                          //
//                                                                          //
//////////////////////////////////////////////////////////////////////////////
#include "MPhoton.h"

#include <iostream.h>

#include <TF1.h>
#include <TH1.h>
#include <TPad.h>
#include <TCanvas.h>

ClassImp(MPhoton);

Double_t MPhoton::Sigma_gg(Double_t *x, Double_t *k)
{
    const Double_t m2 = x[0]; // m2: (E0/sqrt(s))^2

    const Double_t r0 = 2.81794092e-15; // [m] = e^2/4/pi/m/eps0/c^2

    const Double_t beta2 = 1.-m2;
    const Double_t beta  = sqrt(beta2);

    const Double_t p1 = r0*r0*TMath::Pi()/2;

    // ----- Extreme Relativistic -------
    // return p1*2 * m*m*m* (log(2./m)-1);
    // ----------------------------------

    const Double_t p2 = 3.-beta2*beta2;
    const Double_t p3 = log((1.+beta)/(1.-beta));
    const Double_t p4 = beta*2*(1.+m2);

    const Double_t sigma = p1*m2*(p2*p3-p4); // [m^2]

    return sigma;
}

Double_t MPhoton::Int1(Double_t *x, Double_t *k)
{
    const Double_t costheta = x[0];

    const Double_t Eg = k[0];
    const Double_t Ep = k[1];

    const Double_t E0 = 511e-6; // [GeV]

    Double_t s = E0/Eg*E0/Ep/(1.-costheta)/2;
    if (s>1)   // Why is this necessary???
        return 0;

    const Double_t sigma = Sigma_gg(&s);  // [m^2]

    return sigma/2 * (1.-costheta); // [m^2]
}

Double_t MPhoton::Int2(Double_t *x, Double_t *k)
{
    const Double_t E0 = 511e-6; // [GeV]

    Double_t Ep = x[0];
    Double_t z  = k[1];

    const Double_t Eg = k[0];

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

    static TF1 f("int1", Int1, 0, 0, 2);

    const Double_t from   = -1.0;
    const Double_t to     = 1.-E0*E0/(2.*Eg*Ep);              // Originally Was: 1.
    const Double_t int1   = f.Integral(from, to, val, 1e-2);  // [m^2]
    const Double_t planck = MParticle::Planck(&Ep, &z);       // [GeV^2]

    const Double_t res = planck * int1;

    return res; // [GeV^2 m^2]
}

// --------------------------------------------------------------------------
//
// Returns 0 in case IL becomes (numerically) infinite.
//
Double_t MPhoton::InteractionLength(Double_t *x, Double_t *k)
{
    Double_t E0 = 511e-6;                  // [GeV]
    Double_t c  = 299792458;               // [m/s]
    Double_t e  = 1.602176462e-19;         // [C]
    Double_t h  = 1e-9/e*6.62606876e-34;   // [GeVs]
    Double_t hc = h*c;                     // [GeVm]
    Double_t pc = 1./3.258;                // [pc/ly]
    Double_t ly = 3600.*24.*365.*c;        // [m/ly]

    Double_t Eg = x[0];
    Double_t z  = k ? k[0] : 0;

    if (Eg<100)
        return 1e50;

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

    static TF1 f("int2", Int2, 0, 0, 2);

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

    if (Eg<5e4)
        //inf = 3e-11*(z+1)*pow(10, 4.7*0.5-log10(Eg)*0.5);
        inf = 3e-11*(z+1)*pow(10, 4.7-log10(Eg));
        //inf = 3e-11*(z+1)*pow(10, 7.0-log10(Eg)*1.5);
        //inf = 3e-11*(z+1)*pow(10, 8.2-log10(Eg)*1.75);
        //inf = 3e-11*(z+1)*pow(10, 9.4-log10(Eg)*2);

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

    if (int2==0)
    {
        //cout << "--->  Int2==0 <---" << endl;
        return 0;
    }

    /* Planck constants: konst */
    const Double_t konst = 4.*TMath::Pi() * 2. / (hc*hc*hc);
    int2 *= konst;           // [1 / m]

    Double_t res = 1./ int2; // [m]
    res *= pc/ly * 1e-3;     // [kpc]

    if (res > 1e50) return 1e50;
    if (res < 0)    return 1e35;

    return res;              //[kpc]
}

Double_t MPhoton::GetInteractionLength(Double_t energy, Double_t z)
{
    return InteractionLength(&energy, &z);
}

Double_t MPhoton::GetInteractionLength() const
{
    return InteractionLength((Double_t*)&fEnergy, (Double_t*)&fZ);
}

void MPhoton::DrawInteractionLength(Double_t z, Double_t from, Double_t to, Option_t *opt)
{
    if (!gPad)
        new TCanvas("ILPhoton", "Mean Interaction Length Photon");
    else
        gPad->GetVirtCanvas()->cd(4);

    TF1 f1("length", InteractionLength, from, to, 1);
    f1.SetParameter(0, z);

    gPad->SetLogx();
    gPad->SetLogy();
    gPad->SetGrid();
    f1.SetMinimum(1);
    f1.SetMaximum(1e9);
    f1.SetLineWidth(1);

    TH1 &h=*f1.DrawCopy(opt)->GetHistogram();

    h.SetTitle("Mean Interaction Length (Photon)");
    h.SetXTitle("E [GeV]");
    h.SetYTitle("x [kpc]");

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

void MPhoton::DrawInteractionLength() const
{
    DrawInteractionLength(fZ);
}

void MPhoton::Fill(TH1 &h, Double_t idx, Double_t w) const
{
    h.Fill(fEnergy, pow(fEnergy, idx)*w);
}