/* ======================================================================== *\
!
! *
! * 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 "MGenIRPhoton.h"

#include <TMath.h>
#include <TF1.h>

#include "MParList.h"
#include "MPhoton.h"

ClassImp(MGenIRPhoton);

/*
Double_t Planck(Double_t *x, Double_t *k)
{
    const Double_t E   = x[0];             //[eV]

    const Double_t T   = 2.87;             //[K]

    const Double_t e   = 1.602176462e-19;  //[C]
    const Double_t kB  = 1.3806503e-23/e;  //[eV/K]

    const Double_t E3  = E*E*E;
    const Double_t EkT = E/kB/T;

    //const Double_t h     = 1e-9/e*6.62606876e-34; //[GeVs]
    //const Double_t h2    = h*h;
    //const Double_t c     = 299792458;             //[m/s]
    //const Double_t c3    = c*c*c;
    //const Double_t konst = 4/c3/h2;

    return E3/(exp(EkT)-1);
}
*/

Double_t MGenIRPhoton::Planck(Double_t *x, Double_t *k=NULL)
{
    //
    // Planck, per unit volume, per unit energy
    //
    // constants moved out of function
    //
    Double_t E   = x[0];                    // [GeV]
    Double_t z   = k ? k[0] : 0;

    Double_t T   = 2.96*(z+1);              // [K]
    Double_t e   = 1.602176462e-19;         // [C]
    Double_t kB  = 1e-9/e*1.3806503e-23;    // [GeV/K]

    Double_t EkT = E/kB/T;

    /*
     //Double_t c   = 299792458;             // [m/s]
     //Double_t h   = 1e-9/e*6.62606876e-34; // [GeVs]
     //Double_t hc  = h*c;                   // [GeVm]
     Double_t konst = 4.*TMath::Pi() * 2. / (hc*hc*hc);
     return konst * E*E / (exp(EkT)-1.); // [1 / GeV / m^3 ]
     */

    return E*E / (exp(EkT)-1.); // [GeV^2]
}


// --------------------------------------------------------------------------
MGenIRPhoton::MGenIRPhoton()
{
    fSrc = new TF1("Planck", Planck, 0, .5e-2, 0);
}

MGenIRPhoton::~MGenIRPhoton()
{
    delete fSrc;
}

// --------------------------------------------------------------------------
MParticle *MGenIRPhoton::GetRandom()
{
    const Double_t pi2 = TMath::Pi() * 2;

    MPhoton *p = new MPhoton;
    /*
     MParticle *p = new MParticle(MParticle::kEGamma);
     p->SetAngle(fRand.Uniform(pi2));
     p->SetEnergy(fSrc->GetRandom()*1e-9);
     */

    return p;
}