source: trunk/WuerzburgSoft/Thomas/mphys/MPhoton.cc@ 1350

/* ======================================================================== *\
!
! *
! * 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>

ClassImp(MPhoton);

Double_t MPhoton::Planck(Double_t *x, Double_t *k=NULL)
{
    //
    // Planck, per unit volume, per unit energy
    //
    // constants moved out of function, see below
    //
    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]
}

Double_t MPhoton::Sigma_gg(Double_t *x, Double_t *k=NULL)
{
    Double_t s = x[0]; // omega: CM mass

    Double_t E0 = 511e-6;         // [GeV]
    Double_t r0 = 2.81794092e-15; // [m] = e^2/4/pi/m/eps0/c^2

    Double_t m  = E0/s;

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

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

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

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

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

    return sigma;
}

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

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

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

    Double_t s = Eg*Ep/E0*(1.-costheta)*2;

    if (s<E0)   // Why is this necessary???
        return 0;

    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)
{
    Double_t E0 = 511e-6; // [GeV]

    Double_t Ep = x[0];
    Double_t Eg = k[0];

    Double_t z  = k[1];

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

    Double_t from = -1.0;
    Double_t to   = 1.-E0*E0/2./Eg/Ep; // Originally Was: 1.

    TF1 f("int1", Int1, from, to, 2);
    Double_t int1   = f.Integral(from, to, val);  // [m^2]
    Double_t planck = Planck(&Ep, &z);            // [GeV^2]

    Double_t res = planck * int1;

    res *= Eg/E0*1e-9; // FIXME!!!!!!!!!! WHICH FACTOR IS THIS????

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

// --------------------------------------------------------------------------
//
// Returns 0 in case IL becomes (numerically) infinite.
//
Double_t MPhoton::InteractionLength(Double_t *x, Double_t *k=NULL)
{
    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;

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

    Double_t lolim = E0*E0/Eg;
    Double_t inf   = 4e-12; //E0*E0/Eg * sqrt(Eg);

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

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

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

    /* Planck constants: konst */
    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);
}