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

/* ======================================================================== *\
!
! *
! * 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::Planck(Double_t *x, Double_t *k=NULL)
{
    //
    // Planck, per unit volume, per unit energy
    //
    // constants moved out of function, see below
    //
    const Double_t E   = x[0];                    // [GeV]
    const Double_t z   = k ? k[0] : 0;

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

    const 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)
{
    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=NULL)
{
    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 };

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

    TF1 f("int1", Int1, from, to, 2);

    const Double_t int1   = f.Integral(from, to, val);  // [m^2]
    const Double_t planck = 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=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 > 1e-8 ? E0*E0/Eg : 1e-8/Eg;
    Double_t inf   = 3e-11; 

    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);
}

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

    TF1 f1("length", InteractionLength, 1e4, 1e11, 1);
    f1.SetParameter(0, z);

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

    TH1 &h=*f1.DrawCopy()->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);
}