source: trunk/WuerzburgSoft/Thomas/mphys/MParticle.cc@ 1365

///////////////////////////////////////////////////////////////////////
//
// MParticle
//
///////////////////////////////////////////////////////////////////////

#include "MParticle.h"

#include <TRandom.h>
#include <TMatrixD.h>

ClassImp(MParticle);

/**************************************************
 *
 * H0 = 50./3.0857e19; // [km / Mpc s] --> [s^-1]
 *
 **************************************************/

Double_t MParticle::ZofR(Double_t *x, Double_t *k)
{
    /*
     const Double_t c  = 299792458;        // [m/s]
     const Double_t H0 = 50./3.0857e19;    // [km / Mpc s] --> [s^-1]

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

     const Double_t R = r*H0/c;            // [1]

     return (R+1+sqrt(R*2+1))/2 - 1;
    */
    const Double_t c  = 299792458;        // [m/s]
    const Double_t H0 = 50./3.0857e19;    // [km / Mpc s] --> [s^-1]

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

    const Double_t R = 1./(1-r*H0/c/2);   // [1]

    return R*R - 1;
}

Double_t MParticle::RofZ(Double_t *x, Double_t *k)
{
    /*
     Double_t z1 = x[0] + 1;

     const Double_t c  = 299792458;                 // [m/s]
     const Double_t H0 = 50./3.0857e19;             // [km / Mpc s] --> [s^-1]

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

     const Double_t R = c/H0 * 2 * (z1 - sqrt(z1)); // [m]

     return  R * pc/ly/1e3;                         // [kpc]
     */
     Double_t z1 = x[0] + 1;

     const Double_t c  = 299792458;                   // [m/s]
     const Double_t H0 = 50./3.0857e19;               // [km / Mpc s] --> [s^-1]

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

     const Double_t R = c/H0 * 2 * (1 - 1./sqrt(z1)); // [m]

     return R * pc/ly/1e3;                            // [kpc]
}

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

MParticle::MParticle(ParticleType_t t, const char *name, const char *title)
    : fPType(t), fZ(0), fR(0), fPhi(0), fTheta(0), fPsi(0)
{
    //
    //  default constructor
    //
}

void MParticle::SetNewDirection(Double_t theta, Double_t phi)
{
    TMatrixD B(3, 3);

    B(0, 0) =  cos(fTheta)*cos(fPsi);
    B(1, 0) =  cos(fTheta)*sin(fPsi);
    B(2, 0) = -sin(fTheta);

    B(0, 1) = -sin(fPsi);
    B(1, 1) =  cos(fPsi);
    B(2, 1) =  0;

    B(0, 2) = sin(fTheta)*cos(fPsi);
    B(1, 2) = sin(fTheta)*sin(fPsi);
    B(2, 2) = cos(fTheta);

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

    TVectorD r(3);

    r(0) = sin(theta)*cos(phi);
    r(1) = sin(theta)*sin(phi);
    r(2) = cos(theta);

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

    r *= B;

    fTheta = sqrt(r(0)*r(0)+r(1)*r(1)); // Numerically bad: acos(r(2));
    fPsi   = atan2(r(1), r(0));

    /*
     if (fTheta*2 > TMath::Pi())
        fTheta = fabs(fTheta-TMath::Pi());
     */
}

Bool_t MParticle::SetNewPosition(Double_t dr)
{
    Bool_t rc=kTRUE;

    TVectorD x(3);

    x(0) = sin(fTheta)*cos(fPsi);
    x(1) = sin(fTheta)*sin(fPsi);
    x(2) = cos(fTheta);

    x *= dr;

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

    const Double_t R = RofZ(&fZ);

    if (x(2) > R*cos(fTheta))
    {
        x *= R/dr;
        rc = kFALSE;
    }

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

    TVectorD r(3);

    r(0) = fR*cos(fPhi);
    r(1) = fR*sin(fPhi);
    r(2) = R;

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

    r -= x;

    fR   = sqrt(r(0)*r(0)+r(1)*r(1));
    fPhi = atan2(r(1), r(0));
    fZ   = ZofR(&r(2));

    return rc;
}

Bool_t MParticle::SetNewPosition()
{
    static TRandom rand(0);
    Double_t r = rand.Exp(GetInteractionLength());
    return SetNewPosition(r);
}