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

#include "MParticle.h"

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

ClassImp(MParticle);

/***********************
 *
 * calc r from z:
 *
 *        R = 2*c/H0 *(1+z - sqrt(1+z))
 *
 *        z = -0.5 * (3 + R' +/- sqrt(R'+1))   R' = R*H0/c
 *
 *        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;
}

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]
}

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
    //  creates an a list of histograms for all pixels and both gain channels
    //

    //
    //   set the name and title of this object
    //
    
    fName  = name  ? name  : "MParticle";
    fTitle = title ? title : "Container for a particle";
}

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

#include <iostream.h>
#include <iomanip.h>

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