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

Timestamp:

07/29/02 09:10:23 (22 years ago)

Author:

tbretz

Message:

*** empty log message ***

File:

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

trunk/WuerzburgSoft/Thomas/mphys/MParticle.cc

@@ -9 +9 @@
 #include <TRandom.h>
 #include <TMatrixD.h>
+#include <TRotation.h>
 
 ClassImp(MParticle);
@@ -39 +40 @@
     const Double_t r  = x[0] /pc*ly*1e3;  // [m]
 
-    const Double_t R = 1./(1-r*H0/c/2);   // [1]
+    const Double_t R = 1./(1.-r*H0/c/2);   // [1]
 
     return R*R - 1;
@@ -67 +68 @@
      const Double_t pc = 1./3.258;                    // [pc/ly]
 
-     const Double_t R = c/H0 * 2 * (1 - 1./sqrt(z1)); // [m]
+     const Double_t R = c/H0 * 2 * (1. - 1./sqrt(z1)); // [m]
 
      return R * pc/ly/1e3;                            // [kpc]
@@ -85 +86 @@
     if (!isloaded)
     {
-        Double_t c   = 299792458;             // [m/s]
-        Double_t e   = 1.602176462e-19;       // [C]
-        Double_t h   = 1e-9/e*6.62606876e-34; // [GeVs]
-        Double_t hc  = h*c;                   // [GeVm]
+        Double_t c     = 299792458;             // [m/s]
+        Double_t e     = 1.602176462e-19;       // [C]
+        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);
 
@@ -123 +124 @@
     }
 
+    static TAxis &axez = *hist2->GetXaxis();
+    static TAxis &axee = *hist2->GetYaxis();
+
     //
     // y = (y1-y0)/(x1-x0) * (x-x0) + y0
@@ -129 +133 @@
     Double_t E = x[0];
 
-    Int_t i = hist2->GetXaxis()->FindFixBin(z);
-    Int_t j = hist2->GetYaxis()->FindFixBin(E);
-
-    Double_t z1 = hist2->GetXaxis()->GetBinLowEdge(i+1);
-    Double_t z0 = hist2->GetXaxis()->GetBinLowEdge(i);
-
-    Double_t E1 = hist2->GetYaxis()->GetBinLowEdge(j+1);
-    Double_t E0 = hist2->GetYaxis()->GetBinLowEdge(j);
+    Int_t i = axez.FindFixBin(z);
+    Int_t j = axee.FindFixBin(E);
+
+    Double_t z1 = axez.GetBinLowEdge(i+1);
+    Double_t z0 = axez.GetBinLowEdge(i);
+
+    Double_t E1 = axee.GetBinLowEdge(j+1);
+    Double_t E0 = axee.GetBinLowEdge(j);
 
     Double_t n00 = hist2->GetBinContent(i,   j);
@@ -207 +211 @@
 void MParticle::InitRandom()
 {
-    static TRandom rnd(0);
-    fPhi = rnd.Uniform(TMath::Pi()*2);
-    fPsi = rnd.Uniform(TMath::Pi()*2);
+    fPhi = gRandom->Uniform(TMath::Pi()*2);
+    fPsi = gRandom->Uniform(TMath::Pi()*2);
 }
 
 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);
+    static TMatrixD B(3, 3);
+
+    const Double_t ct = cos(fTheta);
+    const Double_t st = sin(fTheta);
+    const Double_t cp = cos(fPsi);
+    const Double_t sp = sin(fPsi);
+
+    /*
+    TRotation B( ct*cp, ct*sp, -st,
+                -sp,    cp,     0,
+                 st*cp, st*sp,  ct);
+                 */
+
+    // first row
+    B(0, 0) =  ct*cp;
+    B(1, 0) =  ct*sp;
+    B(2, 0) = -st;
+
+    // second row
+    B(0, 1) = -sp;
+    B(1, 1) =  cp;
     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);
+    // third row
+    B(0, 2) = st*cp;
+    B(1, 2) = st*sp;
+    B(2, 2) = ct;
+
+    // ------------------------------
+
+    static TVectorD r(3);
+
+    const Double_t sint = sin(theta);
+
+    /*
+     TVector3 r(sint*cos(phi), sint*sin(phi), cos(theta));
+     */
+
+    r(0) = sint*cos(phi);
+    r(1) = sint*sin(phi);
     r(2) = cos(theta);
 
@@ -253 +276 @@
     Bool_t rc=kTRUE;
 
-    TVectorD x(3);
-
-    x(0) = sin(fTheta)*cos(fPsi);
-    x(1) = sin(fTheta)*sin(fPsi);
-    x(2) = cos(fTheta);
-
-    // ------------------------------
-
-    const Double_t R = RofZ(&fZ);
-
-    const Double_t dx = R - dr*x(2);
-
+    const Double_t st = sin(fTheta);
+
+    TVector3 d(st*cos(fPsi), st*sin(fPsi), cos(fTheta));
+
+    // ------------------------------
+
+    const Double_t R  = RofZ(&fZ);
+    const Double_t dx = R - dr*d.z();
     if (dx < 0)
     {
-        dr = R/x(2); // R>0 --> x(2)>0
+        dr = R/d.z(); // R>0 --> x(2)>0
         rc = kFALSE;
     }
-    else
-        fX += dr*(1.-x(2));
-
-    x  *= dr;
-
-    // ------------------------------
-
-    TVectorD r(3);
-
-    r(0) = fR*cos(fPhi);
-    r(1) = fR*sin(fPhi);
-    r(2) = R;
-
-    // ------------------------------
-
-    r -= x;
+
+    fX += dr*(1.-d(2));
+    d  *= dr;
+
+    // ------------------------------
+
+    TVector3 r(fR*cos(fPhi), fR*sin(fPhi), R);
+
+    r -= d;
+
+    // ------------------------------
 
     if (fR!=0)
-        fPhi = atan2(r(1), r(0));
-
-    fR = sqrt(r(0)*r(0)+r(1)*r(1));
+        fPhi = atan2(r.y(), r.x());
+    fR = sqrt(r.x()*r.x()+r.y()*r.y());
     fZ = ZofR(&r(2));
 
@@ -298 +311 @@
 Bool_t MParticle::SetNewPosition()
 {
-    static TRandom rand(0);
-    Double_t r = rand.Exp(GetInteractionLength());
+    Double_t r = gRandom->Exp(GetInteractionLength());
     return SetNewPosition(r);
 }