source: trunk/WuerzburgSoft/Thomas/mphys/MElectron.cc@ 1369

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

#include <iostream.h>

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

#include "MPhoton.h"

ClassImp(MElectron);

Double_t MElectron::Li(Double_t *x, Double_t *k)
{
    const Double_t t = x[0];
    return log(1.-t)/t;
}

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

    const Double_t eps = fabs(t*1e-2);

    Double_t disum = t;
    Double_t add = 0;

    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 MElectron::Li2(Double_t *x, Double_t *k=NULL)
{
    //
    // Dilog, Li2
    // ----------
    //
    //   Integral(0, 1) = konst;
    //   Double_t konst = 1./6*TMath::Pi()*TMath::Pi();
    //
    //   x[0]: z
    //
    const Double_t z = x[0];

    if (fabs(z)<1)
        return DiSum(x);

    // TF1 IntLi("Li", Li, 0, z, 0);
    static TF1 IntLi("Li", Li, 0, 0, 0);
    const Double_t integ = IntLi.Integral(0, z, (Double_t*)NULL, 1e-2);
    return -integ;
}

Double_t MElectron::Flim(Double_t *x, Double_t *k=NULL) // F(omegap)-F(omegam)  mit  b-->1   (Maple)
{
    const Double_t w = x[0];

    const Double_t w4   = w*4;
    const Double_t wsqr = w*w;

    const Double_t u1 = (w*wsqr*16 + wsqr*40 + w*17 + 2)*log(w4 + 1);
    const Double_t u2 = -w4*(wsqr*2 + w*9 + 2);
    const Double_t d  = w4*(w4 + 1);

    Double_t s = -w*2*(1+1); // -2*omega*(1+beta)
    const Double_t li2 = Li2(&s);

    const Double_t res = (u1+u2)/d + li2;

    return res; //<1e-10? 0 : res;
}

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

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

    const Double_t E = k[0];

    Double_t omega  = epsilon*E/(E0*E0);

    const Double_t n = MParticle::Planck(&epsilon, &z)/epsilon/epsilon; // [1]
    return Flim(&omega)*n;
}

Double_t MElectron::InteractionLength(Double_t *E, Double_t *k=NULL)
{
    // E    = electron energy,          ~  TeV(?)  1e12
    // e    = photon energy,            ~  meV(?)  1e-3
    // mc^2 = electron rest mass energy ~.5keV(?)  .5e3
    //
    // x^-1 = int( n(epsilon)/2beta * ((mc^2)^2/eE)^2 * int ( omega*sigma(omega), omega=o-..o+), epsilon=0..inf)
    //
    // o+/- = omage_0 (1 +- beta)
    //
    // omega_0 = eE/(mc^2)^2  ~1e12*1e-3/.25e6=4e3
    //
    // --> x^-1 = (alpha*hc)^2/4pibetaE^2 * int(n(epsilon)/epsilon^2 *( F(o+)-F(o-)), epsilon=0..inf)
    //
    // F(o) = -o/4 + (9/4 + 1/o + o/2) * ln(1+2o) + 1/8(1+2o) - 3/8 + Li2(-2o)
    //
    // Li2(x) = int(ln(1-t)/t, t=0..x)
    //
    // F(o+)~F(2o) = -o/2 + (9/4 + 1/2o + o) * ln(1+4o) + 1/8(1+4o) - 3/8 + Li2(-4o)
    // F(o-)~F(0)  = 14/8 = 1.75

    const Double_t E0     = 511e-6;                        // [GeV]
    const Double_t E02    = E0*E0;                         // [GeV^2]
    const Double_t c      = 299792458;                     // [m/s]
    const Double_t e      = 1.602176462e-19;               // [C]
    const Double_t h      = 1e-9/e*6.62606876e-34;         // [GeVs]
    const Double_t hc     = h*c;                           // [GeVm]
    const Double_t alpha  = 1./137.;                       // [1]

    const Double_t z      = k ? k[0] : 0;

    /* -------------- old ----------------
       Double_t from   = 1e-15;
       Double_t to     = 1e-11;
       eps = [default];
       -----------------------------------
    */
    static TF1 func("Compton", Compton, 0, 0, 2);         // [0, inf]

    const Double_t from   = 1e-17;
    const Double_t to     = 1e-11;

    Double_t val[3] = { E[0], z };                        // E[GeV]

    Double_t integ  = func.Integral(from, to, val, 1e-2); // [Gev] [0, inf]

    const Double_t aE     = alpha/E[0];                   // [1/GeV]

    const Double_t beta   = 1;

    const Double_t konst  = 2.*E02/hc/beta;               // [1 / GeV m]
    const Double_t ret    = konst * (aE*aE) * integ;      // [1 / m]

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

    return (1./ret)/ly*pc/1000;                           // [kpc]
}

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

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

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

inline Double_t MElectron::p_e(Double_t *x, Double_t *k)
{
    Double_t e = pow(10, x[0]);
    return Compton(&e, k);
    /*
     Double_t z  = k[1];

     const Double_t E  = k[0];

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

     Double_t omega = e*E/E02;

     const Double_t n = Planck(&e, &z);

     const Double_t F = Flim(&omega)/omega/omega;

     return n*F*1e26;
     */
}

Double_t MElectron::G_q(Double_t *x, Double_t *k)
{
    const Double_t q     = x[0];
    const Double_t Gamma = k[0];

    const Double_t Gq = Gamma*q;

    const Double_t s1 = 2.*q*log(q);
    const Double_t s2 = (1.+2.*q);
    const Double_t s3 = (Gq*Gq)/(1.+Gq)/2.;

    return s1+(s2+s3)*(1.-q);
}


Double_t MElectron::EnergyLoss(Double_t *x, Double_t *k, Double_t *ep)
{
    const Double_t E = x[0];
    const Double_t z = k ? k[0] : 0;

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

    const Double_t lolim = -log10(E)/7*4-13;

    TF1 fP("p_e", p_e, lolim, -10.8, 2);
    TF1 fQ("G",   G_q, 0, 1., 1);

    fP.SetParameter(0, E);
    fP.SetParameter(1, z);

    const Double_t e = pow(10, fP.GetRandom());

    if (ep)
        *ep = e;

    const Double_t omega = e*E/E0/E0;
    const Double_t Gamma = 4.*omega;

    fQ.SetParameter(0, Gamma);

    const Double_t q  = fQ.GetRandom();
    const Double_t Gq = Gamma*q;

    const Double_t e1 = Gq*E/(1.+Gq);

    return e1;
}

Double_t MElectron::GetEnergyLoss(Double_t E, Double_t z, Double_t *ep)
{
    return EnergyLoss(&E, &z);
}

Double_t MElectron::GetEnergyLoss(Double_t *ep) const
{
    return EnergyLoss((Double_t*)&fEnergy, (Double_t*)&fZ, ep);
}

MPhoton *MElectron::DoInvCompton(Double_t theta)
{
    static TRandom rand(0);

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

    Double_t epsilon;
    const Double_t e = GetEnergyLoss(&epsilon);

    // er: photon energy before interaction, rest frame
    // e:  photon energy after interaction, lab

    const Double_t gamma     = fEnergy/E0;
    const Double_t beta      = sqrt(1.-1./(gamma*gamma));

    const Double_t f = fEnergy/e;

    Double_t t=-1;
    Double_t arg;
    do
    {
        if (t>0)
            cout << "~" << flush;
        t = rand.Uniform(TMath::Pi()/2)+TMath::Pi()*3/4;
        Double_t er  = gamma*epsilon*(1.-beta*cos(t)); // photon energy rest frame
        arg = (f - E0/er - 1)/(f*beta+1);

    } while (arg<-1 || arg>1);

    const Double_t theta1s = acos(arg);
    const Double_t thetas = atan(sin(t)/(gamma*(cos(t)-beta)));

    const Double_t thetastar = thetas-theta1s;

    const Double_t theta1 = atan(sin(thetastar)/(gamma*(cos(thetastar)+beta)));

    fEnergy -= e;

    const Double_t phi = rand.Uniform(TMath::Pi()*2);

    MPhoton &p = *new MPhoton(e, fZ);
    p = *this;
    p.SetNewDirection(theta1, phi);

    const Double_t beta2  = sqrt(1.-E0/fEnergy*E0/fEnergy);
    const Double_t theta2 = asin((epsilon*sin(t)-e*sin(theta1))/fEnergy/beta2);

    SetNewDirection(theta2, phi);

    return &p;
}

void MElectron::DrawInteractionLength(Double_t z)
{
    if (!gPad)
        new TCanvas("IL_Electron", "Mean Interaction Length Electron");
    else
        gPad->GetVirtCanvas()->cd(3);

    TF1 f2("length", MElectron::InteractionLength, 1e3, 1e10, 0);
    f2.SetParameter(0, z);

    gPad->SetLogx();
    gPad->SetLogy();
    gPad->SetGrid();
    f2.SetLineWidth(1);

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

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

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

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