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

Timestamp:

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

Author:

tbretz

Message:

*** empty log message ***

File:

• trunk/WuerzburgSoft/Thomas/mphys/MElectron.cc (modified) (13 diffs)

trunk/WuerzburgSoft/Thomas/mphys/MElectron.cc

@@ -129 +129 @@
     const Double_t E = k[0];
 
-    Double_t omega  = epsilon*E/(E0*E0);
+    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(&omega)*n;
+    return flim*n;
 }
 
@@ -250 +270 @@
     const Double_t E0 = 511e-6; //[GeV]
 
-    const Double_t lolim = -log10(E)/7*4-13;
-
-    TF1 fP("p_e", p_e, lolim, -10.6, 2);
-    TF1 fQ("G",   G_q, 0, 1., 1);
-
+    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);
@@ -288 +312 @@
 MPhoton *MElectron::DoInvCompton(Double_t theta)
 {
-    static TRandom rand(0);
-
     const Double_t E0 = 511e-6; //[GeV]
 
@@ -304 +326 @@
     const Double_t f = fEnergy/e;
 
-    Double_t t=-1;
+    Double_t t=-10;
     Double_t arg;
     do
     {
-        if (t>0)
+        if (t>-10)
             cout << "~" << flush;
-        t = rand.Uniform(TMath::Pi())+TMath::Pi()/2;
-        Double_t er  = epsilon*(gamma-gammabeta*cos(t)); // photon energy rest frame
-        arg = (f - E0/er - 1)/(sqrt(fEnergy*fEnergy-E0*E0)/e+1);
+
+        //
+        // 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 thetas  = atan(sin(t)/(gamma*cos(t)-gammabeta));
 
     const Double_t thetastar = thetas-theta1s;
@@ -325 +377 @@
     fEnergy -= e;
 
-    const Double_t phi = rand.Uniform(TMath::Pi()*2);
-
+    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);
 
@@ -374 +429 @@
         return SetNewPosition();
 
-    static TRandom rand(0);
-    Double_t x = rand.Exp(GetInteractionLength());
+    Double_t x = gRandom->Exp(GetInteractionLength());
 
     // -----------------------------
@@ -381 +435 @@
     const Double_t E0 = 511e-6;            //[GeV]
 
-    Double_t B_theta = rand.Uniform(TMath::Pi());
-    Double_t B_phi   = rand.Uniform(TMath::Pi()*2);
-
-    TMatrixD M(3,3);
+    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);
@@ -398 +452 @@
     M(2, 2) =  0;
 
-    const Double_t beta = sqrt(1-E0/fEnergy*E0/fEnergy);
+    const Double_t beta = sqrt(1.-E0/fEnergy*E0/fEnergy);
 
     //
@@ -404 +458 @@
     // which the x-axis is identical with the B-Field
     //
-    TVectorD v(3);
+    static TVectorD v(3);
     v(0) = beta*sin(fTheta)*cos(fPsi);
     v(1) = beta*sin(fTheta)*sin(fPsi);
@@ -418 +472 @@
     // -----------------------------
 
-    TMatrixD N(3,3);
+    static TMatrixD N(3,3);
     N(0, 0) =  1;
     N(1, 0) =  0;
@@ -436 +490 @@
                                     // v(2) = 0
     // -----------------------------
-    TVectorD p(3);
+    static TVectorD p(3);
 
     Double_t rho = 0;    
@@ -469 +523 @@
     }
 
-    TVectorD s(3);
+    static TVectorD s(3);
     s(0) = fR*cos(fPhi);
     s(1) = fR*sin(fPhi);
@@ -494 +548 @@
     // -----------------------------
 
-    TVectorD w(3);
+    static TVectorD w(3);
     w(0) = beta_p/beta;
     w(1) = beta_o/beta*cos(rho);