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

/* ======================================================================== *\
!
! *
! * 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 <math.h>     // for aqlphas
#include <iostream.h>

#include <TF1.h>
#include <TH1.h>
#include <TMatrixD.h>
#include <TVectorD.h>

#include <TRandom.h>

#include <TPad.h>
#include <TCanvas.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 MElectron::DiSum(Double_t *x, Double_t *k)
{
    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)
{
    //
    // 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) // 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 flim;
    if (epsilon<1e-14)
    {
        const Double_t d = E/(E0*E0);

        Double_t omega1 = 1e-13*d;
        Double_t omega2 = 1e-12*d;

        const Double_t f1 = Flim(&omega1);
        const Double_t f2 = Flim(&omega2);

        const Double_t m = log10(f2/f1);
        const Double_t t = pow(f2, 13)/pow(f1, 12);

        flim = pow(epsilon, m) * t;
    }
    else
    {
        Double_t omega  = epsilon*E/(E0*E0);
        flim = Flim(&omega);
    }

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

Double_t MElectron::InteractionLength(Double_t *E, Double_t *k)
{
    // 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     = 2e-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.5;

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

    fP.SetNpx(50);
    fQ.SetNpx(50);

    fP.SetRange(lolim, -10.6);
    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/E0)*(E/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)
{
    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 gammabeta = sqrt(gamma*gamma-1);

    const Double_t f = fEnergy/e;

    Double_t t=-10;
    Double_t arg;
    do
    {
        if (t>-10)
            cout << "~" << flush;

        //
        // Create an angle uniformly distributed in the range of possible
        // interaction angles due to the known energies.
        //
        // The distibution is a theta-function which is incorrect, but
        // should be correct in a first order...
        //
        t = gRandom->Uniform(TMath::Pi())+TMath::Pi()/2;
        Double_t er = epsilon*(gamma-gammabeta*cos(t)); // photon energy rest frame
        Double_t n  = sqrt(fEnergy*fEnergy-E0*E0)/e+1;
        arg = (f - E0/er - 1)/n;

        /*  ------------------ some hints ------------------------------
         -1 < arg < 1
         -n < f - r/er - 1 < n
         1-f-n < r/er < 1-f+n
         r/(1-f+n) < er < r/(1-f-n)
         r/(1-f+n) < gamma-gammabeta*cos(t) < r/(1-f-n)
         r/(1-f+n)-gamma < -gammabeta*cos(t) < r/(1-f-n)-gamma
         gamma-r/(1-f-n) < gammabeta*cos(t) < gamma-r/(1-f+n)
         (gamma-r/(1-f-n))/gammabeta < cos(t) < (gamma-r/(1-f+n))/gammabeta
         acos((gamma-r/(1-f+n))/gammabeta) < t < acos((gamma-r/(1-f-n))/gammabeta)


         (gamma+E0/(n*arg+1-f)/epsilon)/gammabeta = cos(t);
         1 = (epsilon/E0*cos(t)*gammabeta-gamma)*(n*arg+1-f);

         arg=1:   (gamma+E0/epsilon*(1-f+n))/gammabeta = cos(t);
         arg=-1 : (gamma+E0/epsilon*(1-f-n))/gammabeta = cos(t);

         (E0/(1-f-n)/epsilon+gamma)/gammabeta < cos(t) < (E0/(1-f+n)/epsilon+gamma)/gammabeta
         */

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

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

    const Double_t thetastar = thetas-theta1s;

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

    fEnergy -= e;

    const Double_t phi = gRandom->Uniform(TMath::Pi()*2);

    /*
    MPhoton &p = *new MPhoton(e, fZ);
    p = *this;
    */
    MPhoton &p = *new MPhoton(*this, e);
    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);
}

Bool_t MElectron::SetNewPositionB(Double_t B)
{
    if (B==0)
        return SetNewPosition();

    Double_t x = gRandom->Exp(GetInteractionLength());

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

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

    Double_t B_theta = gRandom->Uniform(TMath::Pi());
    Double_t B_phi   = gRandom->Uniform(TMath::Pi()*2);

    static TMatrixD M(3,3);

    M(0, 0) =  sin(B_theta)*cos(B_phi);
    M(1, 0) =  cos(B_theta)*cos(B_phi);
    M(2, 0) = -sin(B_phi);

    M(0, 1) =  sin(B_theta)*sin(B_phi);
    M(1, 1) =  cos(B_theta)*sin(B_phi);
    M(2, 1) =  cos(B_phi);

    M(0, 2) =  cos(B_theta);
    M(1, 2) = -sin(B_theta);
    M(2, 2) =  0;

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

    //
    // rotate vector of velocity into a system in
    // which the x-axis is identical with the B-Field
    //
    static TVectorD v(3);
    v(0) = beta*sin(fTheta)*cos(fPsi);
    v(1) = beta*sin(fTheta)*sin(fPsi);
    v(2) = beta*cos(fTheta);
    v *= M;

    //
    // Now rotate the system so, that the velocity vector
    // lays in the y-z-plain
    //
    Double_t chi = atan(-v(2)/v(1));

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

    static TMatrixD N(3,3);
    N(0, 0) =  1;
    N(1, 0) =  0;
    N(2, 0) =  0;

    N(0, 1) =  0;
    N(1, 1) = -cos(chi);
    N(2, 1) = -sin(chi);

    N(0, 2) =  0;
    N(1, 2) =  sin(chi);
    N(2, 2) = -cos(chi);
    v *= N;

    Double_t beta_p = v(0);         // beta parallel
    Double_t beta_o = v(1);         // beta orthogonal
                                    // v(2) = 0
    // -----------------------------
    static TVectorD p(3);

    Double_t rho = 0;    
    if (B>0)
    {
        const Double_t c  = 299792458;               // [m/s]
        const Double_t ly = 3600.*24.*365.*c;        // [m/ly]
        const Double_t pc = 1./3.258;                // [pc/ly]

        Double_t r = (fEnergy*1e9)*beta_o/(c*B)*(pc/ly/1e3);  // [kpc]

        rho = beta_o/beta*x/r;                                // [2pi]

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

        if (fabs(rho*3437)>5*60)  // > 5*60min
        {
            cout << "r" << flush;
            return kFALSE;
        }

        p(0) = beta_p/beta*x;
        p(1) = GetPType()==kEElectron?r*sin(rho):-r*sin(rho);
        p(2) = r*(1.-cos(rho));
    }
    else
    {
        p(0) = beta_p/beta*x;
        p(1) = beta_o/beta*x;
        p(2) = 0;
        cout << "------------- HEY! --------------" << endl;
    }

    static TVectorD s(3);
    s(0) = fR*cos(fPhi);
    s(1) = fR*sin(fPhi);
    s(2) = RofZ(&fZ);

    TMatrixD N2(TMatrixD::kTransposed, N);
    TMatrixD M2(TMatrixD::kTransposed, M);
    p *= N2;
    p *= M2;

    if (p(2)<x) // happens sometimes in case B==0
    {
        cout << "----- HA: " << B << " " << x << " " << p(2) << " " << x-p(2) << endl;
        p(2)=x;
    }

    s -= p;

    fR   = sqrt(s(0)*s(0)+s(1)*s(1));
    fPhi = atan2(s(1), s(0));
    fZ   = ZofR(&s(2));
    fX  += x-p(2);

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

    static TVectorD w(3);
    w(0) = beta_p/beta;
    w(1) = beta_o/beta*cos(rho);
    w(2) = beta_o/beta*sin(rho);

    w *= N2;
    w *= M2;

    fPsi   = atan2(w(1), w(0));
    fTheta = asin(sqrt(w(0)*w(0)+w(1)*w(1)));

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

    return fZ<0 ? kFALSE : kTRUE;
}