Index: trunk/WuerzburgSoft/Thomas/mphys/Changelog
===================================================================
--- trunk/WuerzburgSoft/Thomas/mphys/Changelog	(revision 1506)
+++ trunk/WuerzburgSoft/Thomas/mphys/Changelog	(revision 1507)
@@ -1,3 +1,20 @@
                                                                   -*-*- END -*-*-
+ 2002/08/19: Thomas Bretz
+
+   * MCascade.[h,cc]:
+     - removed argument of DoInvCompton
+     - delete GRandom only if set
+     - set gRandom to 0 after deletion
+
+   * MElectron.[h,cc]:
+     - added Sigma_ge
+     - removed argument from DoInvCompton
+     - Implemented Omega_sigmae and RandomThetaE
+     - Changed old wrong dertermination of interaction angle to 
+       correct simulation. 
+     - Reengenieered some formulas to be more correct in terms of numerics
+
+
+
  2002/08/02: Thomas Bretz
 
Index: trunk/WuerzburgSoft/Thomas/mphys/MCascade.cc
===================================================================
--- trunk/WuerzburgSoft/Thomas/mphys/MCascade.cc	(revision 1506)
+++ trunk/WuerzburgSoft/Thomas/mphys/MCascade.cc	(revision 1507)
@@ -50,5 +50,5 @@
 }
 
-Double_t RandomTheta(Double_t Eg, Double_t Ep)
+Double_t RandomThetaG(Double_t Eg, Double_t Ep)
 {
     Double_t E0 = 511e-6; // [GeV]
@@ -59,5 +59,5 @@
         return 0;
 
-    static TF1 func("RndTheta", Sbar_sigmas, 0, 0, 0);
+    static TF1 func("RndThetaG", Sbar_sigmas, 0, 0, 0);
 
     func.SetRange(0, log10(f));
@@ -139,6 +139,5 @@
         }
 
-        // WRONG!
-        MPhoton *p = e.DoInvCompton(0);
+        MPhoton *p = e.DoInvCompton();
 
         fBranchElectrons->GetTree()->Fill();
@@ -200,5 +199,5 @@
 
         Double_t Ep = pow(10, pe);
-        Double_t theta = RandomTheta(Eg, Ep);
+        Double_t theta = RandomThetaG(Eg, Ep);
         if (theta==0)
         {
@@ -316,6 +315,8 @@
 MCascade::MCascade()
 {
+    if (gRandom)
+        delete gRandom;
+
     TRandom r(0);
-    delete gRandom;
     gRandom = new TRandom3(r.GetSeed());
 
@@ -332,4 +333,5 @@
 {
     delete gRandom;
+    gRandom = 0;
 }
 
@@ -373,4 +375,6 @@
 
     // ------------------------------
+
+    cout << endl;
 
     cout << "R = " << fSrcR << "kpc" << endl;
Index: trunk/WuerzburgSoft/Thomas/mphys/MElectron.cc
===================================================================
--- trunk/WuerzburgSoft/Thomas/mphys/MElectron.cc	(revision 1506)
+++ trunk/WuerzburgSoft/Thomas/mphys/MElectron.cc	(revision 1507)
@@ -47,4 +47,28 @@
 ClassImp(MElectron);
 
+Double_t MElectron::Sigma_ge(Double_t *x, Double_t *k)
+{
+    const Double_t E0 = 511e-6;                 // [GeV]
+    const Double_t c  = 299792458;              // [m/s]
+    const Double_t e  = 1.602176462e-19;        // [C]
+    const Double_t h  = 1e-9/e*6.62606876e-34;  // [GeVs]
+    const Double_t a  = 1./137;                 // [1]
+
+    const Double_t o = x[0];
+
+    const Double_t o1  = o+1;
+    const Double_t o21 = o*2+1;
+
+    const Double_t s1 = TMath::Pi()*2;
+    const Double_t s2 = a*h*c/E0;             // [m]
+    const Double_t s3 = o1/(o*o*o);
+    const Double_t s4 = o*2*o1/o21;
+    const Double_t s5 = log(o21);
+    const Double_t s6 = s5/o/2;
+    const Double_t s7 = (o*3+1)/(o21*o21);
+
+    return s2*s2/s1*(s3*(s4-s5)+s6-s7); // [m^2]
+}
+
 Double_t MElectron::Li(Double_t *x, Double_t *k)
 {
@@ -310,5 +334,41 @@
 }
 
-MPhoton *MElectron::DoInvCompton(Double_t theta)
+
+Double_t Omega_sigmae(Double_t *x, Double_t *k)
+{
+    Double_t sbar = pow(10,x[0]);
+
+    Double_t omega = (sbar-1)/2;
+
+    Double_t sigma = MElectron::Sigma_ge(&omega);
+
+    return (sbar-1)*sigma*1e28;
+}
+
+Double_t RandomThetaE(Double_t Ee, Double_t Ep)
+{
+    Double_t E0 = 511e-6; // [GeV]
+
+    Double_t f = 2*Ee/E0*Ep/E0;
+
+    static TF1 func("RndThetaE", Omega_sigmae, 0, 0, 0);
+
+    Double_t beta = sqrt(1-E0/Ee*E0/Ee);
+
+    //func.SetRange(0, log10(1+f*(1+beta)));
+    func.SetRange(log10(1+f*(1-beta)), log10(1+f*(1+beta)));
+    func.SetNpx(50);
+
+    Double_t sbar = pow(10, func.GetRandom());
+
+    Double_t bcost = 1 - (sbar-1)/f;
+    return bcost;
+
+    Double_t theta = acos(bcost/beta);
+
+    return theta;
+}
+
+MPhoton *MElectron::DoInvCompton()
 {
     const Double_t E0 = 511e-6; //[GeV]
@@ -317,61 +377,38 @@
     const Double_t e = GetEnergyLoss(&epsilon);
 
-    // er: photon energy before interaction, rest frame
-    // e:  photon energy after interaction, lab
-
-    const Double_t gamma = fEnergy/E0;
-    // const Double_t beta  = sqrt(1.-1./(gamma*gamma));
+    // epsilon: photon energy befor interaction, lab
+    // e:       photon energy after interaction, lab
+
+    const Double_t gamma     = fEnergy/E0;
     const Double_t gammabeta = sqrt(gamma*gamma-1);
-
-    const Double_t f = fEnergy/e;
-
-    Double_t t=-10;
-    Double_t arg;
-    do
-    {
-        if (t>-10)
-            cout << "~" << flush;
-
-        //
-        // 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 beta      = gammabeta/gamma;
+
+    const Double_t bcost = RandomThetaE(fEnergy, epsilon);
+    const Double_t t     = acos(bcost/beta);
+
+    const Double_t f = epsilon/fEnergy;
+    const Double_t r = gamma*(1-bcost);
+
+    Double_t arg = (1 - 1/(gamma*r) - f)/(beta-f);
+
+    if (arg<-1 || arg>1)
+        cout << "<" << (int)(t*180/TMath::Pi()) << "°>" << flush;
+
+    //
+    // Due to numerical uncertanties arg can be something like:
+    // 1+1e-15 which we cannot allow
+    //
+    if (arg<-1)
+        arg = -1;
+    if (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)/(beta*r)/*(1-bcost)/gammabeta*/);
 
     const Double_t thetastar = thetas-theta1s;
 
-    const Double_t theta1 = atan(sin(thetastar)/(gamma*cos(thetastar)+gammabeta));
+    //    const Double_t theta1 = atan(sin(thetastar)/(gamma*cos(thetastar)+gammabeta));
+    const Double_t theta1 = atan(sin(thetastar)/(gamma*(cos(thetastar)+beta)));
 
     fEnergy -= e;
@@ -379,13 +416,14 @@
     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);
 
-    const Double_t beta2  = sqrt(1.-E0/fEnergy*E0/fEnergy);
-    const Double_t theta2 = asin((epsilon*sin(t)-e*sin(theta1))/fEnergy/beta2);
+    /*
+     const Double_t beta2  = sqrt(1.-E0/fEnergy*E0/fEnergy);
+     const Double_t theta2 = asin((epsilon*sin(t)-e*sin(theta1))/fEnergy/beta2);
+     */
+
+    const Double_t div    = sqrt(fEnergy*fEnergy-E0*E0);
+    const Double_t theta2 = asin((epsilon*sin(t)-e*sin(theta1))/div);
 
     SetNewDirection(theta2, phi);
Index: trunk/WuerzburgSoft/Thomas/mphys/MElectron.h
===================================================================
--- trunk/WuerzburgSoft/Thomas/mphys/MElectron.h	(revision 1506)
+++ trunk/WuerzburgSoft/Thomas/mphys/MElectron.h	(revision 1507)
@@ -24,4 +24,8 @@
     // ----------------------------------------------------------------
 
+    static Double_t Sigma_ge(Double_t *x, Double_t *k=NULL);
+
+    // ----------------------------------------------------------------
+
     static Double_t DiSum(Double_t *x, Double_t *k=NULL);
     static Double_t Li(Double_t *x, Double_t *k);
@@ -45,5 +49,5 @@
     // ----------------------------------------------------------------
 
-    MPhoton *DoInvCompton(Double_t theta);
+    MPhoton *DoInvCompton();
     Bool_t SetNewPositionB(Double_t B);
 
Index: trunk/WuerzburgSoft/Thomas/mphys/MPhoton.h
===================================================================
--- trunk/WuerzburgSoft/Thomas/mphys/MPhoton.h	(revision 1506)
+++ trunk/WuerzburgSoft/Thomas/mphys/MPhoton.h	(revision 1507)
@@ -19,6 +19,6 @@
     }
 
-        MPhoton(MParticle &p) : MParticle(p, MParticle::kEGamma) { }
-        MPhoton(MParticle &p, Double_t e) : MParticle(p, e, MParticle::kEGamma) { }
+    MPhoton(MParticle &p) : MParticle(p, MParticle::kEGamma) { }
+    MPhoton(MParticle &p, Double_t e) : MParticle(p, e, MParticle::kEGamma) { }
 
     void operator=(MParticle &p) { MParticle::operator=(p); }