Changeset 1352 for trunk/WuerzburgSoft/Thomas/mphys/MParticle.cc

Timestamp:

06/10/02 08:49:28 (22 years ago)

Author:

tbretz

Message:

*** empty log message ***

File:

• trunk/WuerzburgSoft/Thomas/mphys/MParticle.cc (modified) (2 diffs)

trunk/WuerzburgSoft/Thomas/mphys/MParticle.cc

@@ -7 +7 @@
 #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 ZofR(Double_t *x, Double_t *k=NULL)
+{
+    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 RofZ(Double_t *x, Double_t *k=NULL)
+{
+    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)
+    : fPType(t), fZ(0), fR(0), fPhi(0), fTheta(0), fPsi(0)
 {
     //
@@ -25 +70 @@
 }
 
+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);
+}