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

Timestamp:

06/13/02 16:36:19 (22 years ago)

Author:

tbretz

Message:

*** empty log message ***

File:

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

trunk/WuerzburgSoft/Thomas/mphys/MElectron.cc

@@ -49 +49 @@
     // constants moved out of function
     //
-    Double_t E   = x[0];                    // [GeV]
-    Double_t z   = k ? k[0] : 0;
-
-    Double_t T   = 2.96*(z+1);              // [K]
-    Double_t e   = 1.602176462e-19;         // [C]
-    Double_t kB  = 1e-9/e*1.3806503e-23;    // [GeV/K]
-
-    Double_t EkT = E/kB/T;
+    const Double_t E   = x[0];                    // [GeV]
+    const Double_t z   = k ? k[0] : 0;
+
+    const Double_t T   = 2.96*(z+1);              // [K]
+    const Double_t e   = 1.602176462e-19;         // [C]
+    const Double_t kB  = 1e-9/e*1.3806503e-23;    // [GeV/K]
+
+    const Double_t EkT = E/kB/T;
 
     /*
@@ -71 +71 @@
 Double_t MElectron::Li(Double_t *x, Double_t *k)
 {
-    Double_t t = x[0];
+    const Double_t t = x[0];
 
     return log(1.-t)/t;
@@ -84 +84 @@
 
     TF1 IntLi("Li", Li, 0, z, 0);
-    Double_t integ = IntLi.Integral(0, z);
+    const Double_t integ = IntLi.Integral(0, z);
 
     /*
@@ -105 +105 @@
 Double_t MElectron::Flim(Double_t *x, Double_t *k=NULL) // F(omegap)-F(omegam)  mit  b-->1   (Maple)
 {
-    Double_t w = x[0];
-
-    Double_t w4   = w*4;
-    Double_t wsqr = w*w;
-
-    Double_t u1 = (w*wsqr*16 + wsqr*40 + w*17 + 2)*log(w4 + 1);
-    Double_t u2 = -w4*(wsqr*2 + w*9 + 2);
-    Double_t d  = w4*(w4 + 1);
+    const Double_t w = x[0];
+
+    const Double_t w4   = w*4;
+    const Double_t wsqr = w*w;
+
+    const Double_t u1 = (w*wsqr*16 + wsqr*40 + w*17 + 2)*log(w4 + 1);
+    const Double_t u2 = -w4*(wsqr*2 + w*9 + 2);
+    const Double_t d  = w4*(w4 + 1);
 
     Double_t s   = -w*2*(1+1); // -2*omega*(1+beta)
-    Double_t li2 = Li2(&s);
-
-    Double_t res = (u1+u2)/d + li2;
+    const Double_t li2 = Li2(&s);
+
+    const Double_t res = (u1+u2)/d + li2;
 
     return res; //<1e-10? 0 : res;
@@ -124 +124 @@
 Double_t MElectron::Compton(Double_t *x, Double_t *k)
 {
-    Double_t E0  = 511e-6; //[GeV]
-    Double_t E02 = E0*E0;
+    const Double_t E0  = 511e-6; //[GeV]
+    const Double_t E02 = E0*E0;
 
     Double_t epsilon = x[0];
@@ -134 +134 @@
     Double_t omega  = epsilon*E/E02;
 
-    Double_t n = Planck(&epsilon, &z)/epsilon/epsilon; // [1]
+    const Double_t n = Planck(&epsilon, &z)/epsilon/epsilon; // [1]
     return Flim(&omega)*n;
 }
@@ -160 +160 @@
     // F(o-)~F(0)  = 14/8 = 1.75
 
-    Double_t E0     = 511e-6;                        // [GeV]
-    Double_t E02    = E0*E0;                         // [GeV^2]
-    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 alpha  = 1./137.;                       // [1]
-
-    Double_t z      = k ? k[0] : 0;
+    const Double_t E0     = 511e-6;                        // [GeV]
+    const Double_t E02    = E0*E0;                         // [GeV^2]
+    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 hc     = h*c;                           // [GeVm]
+    const Double_t alpha  = 1./137.;                       // [1]
+
+    const Double_t z      = k ? k[0] : 0;
 
     // Double_t beta = sqrt(1-E0/E*E0/E);
-    Double_t beta   = 1; //0.999999999999999999999999999;
+    const Double_t beta   = 1; //0.999999999999999999999999999;
 
     Double_t val[3] = { E[0], beta, z };             // E[GeV]
 
-    Double_t from   = 1e-17;
-    Double_t to     = 1e-11;
+    const Double_t from   = 1e-17;
+    const Double_t to     = 1e-11;
 
     /* -------------- old ----------------
@@ -186 +186 @@
     TF1 func("Compton", Compton, from, to, 3);      // [0, inf]
 
-    Double_t integ  = func.Integral(from, to, val, 1e-15); // [Gev] [0, inf]
-
-    Double_t aE     = alpha/E[0];                   // [1/GeV]
-
-    Double_t konst  = 2.*E02/hc/beta;               // [1 / GeV m]
-    Double_t ret    = konst * (aE*aE) * integ;      // [1 / m]
-
-    Double_t ly     = 3600.*24.*365.*c;             // [m/ly]
-    Double_t pc     = 1./3.258;                     // [pc/ly]
+    const Double_t integ  = func.Integral(from, to, val, 1e-15); // [Gev] [0, inf]
+
+    const Double_t aE     = alpha/E[0];                   // [1/GeV]
+
+    const Double_t konst  = 2.*E02/hc/beta;               // [1 / GeV m]
+    const Double_t ret    = konst * (aE*aE) * integ;      // [1 / m]
+
+    const Double_t ly     = 3600.*24.*365.*c;             // [m/ly]
+    const Double_t pc     = 1./3.258;                     // [pc/ly]
 
     return (1./ret)/ly*pc/1000;                     // [kpc]
@@ -214 +214 @@
 {
     Double_t e  = pow(10, x[0]);
-    Double_t E  = k[0];
     Double_t z  = k[1];
 
-    Double_t E0  = 511e-6; //[GeV]
-    Double_t E02 = E0*E0;
+    const Double_t E  = k[0];
+
+    const Double_t E0  = 511e-6; //[GeV]
+    const Double_t E02 = E0*E0;
 
     Double_t omega = e*E/E02;
 
-    Double_t n = Planck(&e, &z);
-
-    Double_t F = Flim(&omega)/omega/omega;
+    const Double_t n = Planck(&e, &z);
+
+    const Double_t F = Flim(&omega)/omega/omega;
 
     return n*F*1e26;
@@ -231 +232 @@
 Double_t MElectron::G_q(Double_t *x, Double_t *k)
 {
-    Double_t q     = x[0];
-    Double_t Gamma = k[0];
-
-    Double_t Gq = Gamma*q;
-
-    Double_t s1 = 2.*q*log(q);
-    Double_t s2 = (1.+2.*q);
-    Double_t s3 = (Gq*Gq)/(1.+Gq)/2.;
+    const Double_t q     = x[0];
+    const Double_t Gamma = k[0];
+
+    const Double_t Gq = Gamma*q;
+
+    const Double_t s1 = 2.*q*log(q);
+    const Double_t s2 = (1.+2.*q);
+    const Double_t s3 = (Gq*Gq)/(1.+Gq)/2.;
 
     return s1+(s2+s3)*(1.-q);
@@ -246 +247 @@
 Double_t MElectron::EnergyLoss(Double_t *x, Double_t *k, Double_t *ep)
 {
-    Double_t E = x[0];
-    Double_t z = k ? k[0] : 0;
-
-    Double_t E0 = 511e-6; //[GeV]
-
-    Double_t lolim = -log10(E)/7*4-13;
+    const Double_t E = x[0];
+    const Double_t z = k ? k[0] : 0;
+
+    const Double_t E0 = 511e-6; //[GeV]
+
+    const Double_t lolim = -log10(E)/7*4-13;
 
     TF1 fP("p_e", p_e, lolim, -10.8, 2);
@@ -259 +260 @@
     fP.SetParameter(1, z);
 
-    Double_t e = pow(10, fP.GetRandom());
+    const Double_t e = pow(10, fP.GetRandom());
 
     if (ep)
         *ep = e;
 
-    Double_t omega = e*E/E0/E0;
-    Double_t Gamma = 4.*omega;
+    const Double_t omega = e*E/E0/E0;
+    const Double_t Gamma = 4.*omega;
 
     // --old-- fQ.SetRange(1e-6, 1./(1.+ 2.*Gamma));
     fQ.SetParameter(0, Gamma);
 
-    Double_t q  = fQ.GetRandom();
-    Double_t Gq = Gamma*q;
-
-    Double_t e1 = Gq*E/(1.+Gq);
+    const Double_t q  = fQ.GetRandom();
+    const Double_t Gq = Gamma*q;
+
+    const Double_t e1 = Gq*E/(1.+Gq);
 
     return e1;
@@ -292 +293 @@
     static TRandom rand(0);
 
-    Double_t E0 = 511e-6; //[GeV]
+    const Double_t E0 = 511e-6; //[GeV]
 
     Double_t epsilon;
-    Double_t e = GetEnergyLoss(&epsilon);
+    const Double_t e = GetEnergyLoss(&epsilon);
 
     // er: photon energy before interaction, rest frame
     // e:  photon energy after interaction, lab
 
-    Double_t gamma     = fEnergy/E0;
-    Double_t beta      = sqrt(1.-1./(gamma*gamma));
-    //Double_t gammabeta = sqrt(gamma*gamma-1);
-
-    Double_t f = fEnergy/e;
-
-    Double_t t;
+    const Double_t gamma     = fEnergy/E0;
+    const Double_t beta      = sqrt(1.-1./(gamma*gamma));
+
+    const Double_t f = fEnergy/e;
+
+    Double_t t=-1;
     Double_t arg;
     do
     {
-        t = rand.Uniform(TMath::Pi()*2);
+        if (t>0)
+            cout << "~" << flush;
+        t = rand.Uniform(TMath::Pi()/2)+TMath::Pi()*3/4;
         Double_t er  = gamma*epsilon*(1.-beta*cos(t)); // photon energy rest frame
         arg = (f - E0/er - 1)/(f*beta+1);
-        cout << "~" << flush;
 
     } while (arg<-1 || arg>1);
 
-    Double_t theta1s = acos(arg);
-    Double_t thetas = atan(sin(t)/(gamma*(cos(t)-beta)));
-
-    Double_t thetastar = thetas-theta1s;
-
-    Double_t theta1 = atan(sin(thetastar)/(gamma*(cos(thetastar)+beta)));
-
-    /*
-     cout << "(" << theta1/TMath::Pi()*180 << ",";
-     cout << theta1s/TMath::Pi()*180<< ",";
-     cout << arg << ")" << flush;
-     */
+    const Double_t theta1s = acos(arg);
+    const Double_t thetas = atan(sin(t)/(gamma*(cos(t)-beta)));
+
+    const Double_t thetastar = thetas-theta1s;
+
+    const Double_t theta1 = atan(sin(thetastar)/(gamma*(cos(thetastar)+beta)));
 
     fEnergy -= e;
@@ -335 +330 @@
     p = *this;
     p.SetNewDirection(theta1, rand.Uniform(TMath::Pi()*2));
-
-    // MISSING: Electron angle
-    //
-    // E1 = fEnergy  (after!)
-    //
-    // sin(t) = (epsilon sin(theta) - e sin(atan(tg)))/sqrt(E1*E1-E0*E0)
 
     return &p;