Index: /trunk/WuerzburgSoft/Thomas/mphys/Changelog
===================================================================
--- /trunk/WuerzburgSoft/Thomas/mphys/Changelog	(revision 1362)
+++ /trunk/WuerzburgSoft/Thomas/mphys/Changelog	(revision 1363)
@@ -4,4 +4,17 @@
    * MParticle.h:
      - Implemented Set/IsPrimary
+     - removed pure abstract method (root can't read this)
+     - added default constructor
+
+   * MElectron.cc:
+     - changed many variables to const
+     - changed the output
+
+   * MPairProduction.cc:
+     - changed many variables to const
+
+   * MPhoton.[h,cc]:
+     - changed the upper and lower limit of the integration in Int2
+     - added upper and lower limit of draw to DrawInteractionLength
 
    * energyloss.C:
@@ -11,4 +24,11 @@
      - changed energy randomizer
      - added root file output
+     - removed log scale from point spread plots
+
+   * analp.C:
+     - added macro to analize the photons from the root file
+
+   * anale.C:
+     - added macro to analize the electrons from the root file
 
 
Index: /trunk/WuerzburgSoft/Thomas/mphys/MElectron.cc
===================================================================
--- /trunk/WuerzburgSoft/Thomas/mphys/MElectron.cc	(revision 1362)
+++ /trunk/WuerzburgSoft/Thomas/mphys/MElectron.cc	(revision 1363)
@@ -49,12 +49,12 @@
     // 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,5 +71,5 @@
 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,5 +84,5 @@
 
     TF1 IntLi("Li", Li, 0, z, 0);
-    Double_t integ = IntLi.Integral(0, z);
+    const Double_t integ = IntLi.Integral(0, z);
 
     /*
@@ -105,17 +105,17 @@
 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,6 +124,6 @@
 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,5 +134,5 @@
     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,21 +160,21 @@
     // 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,13 +186,13 @@
     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,15 +214,16 @@
 {
     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,12 +232,12 @@
 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,10 +247,10 @@
 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,19 +260,19 @@
     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,41 +293,35 @@
     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,10 +330,4 @@
     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;
Index: /trunk/WuerzburgSoft/Thomas/mphys/MPairProduction.cc
===================================================================
--- /trunk/WuerzburgSoft/Thomas/mphys/MPairProduction.cc	(revision 1362)
+++ /trunk/WuerzburgSoft/Thomas/mphys/MPairProduction.cc	(revision 1363)
@@ -79,6 +79,4 @@
     //  theta: interaction angle [rad]
     //
-
-
     const Double_t E0     = 511e-6;              // [GeV]
     const Double_t Eg     = gamma->GetEnergy();  // [GeV]
@@ -98,5 +96,5 @@
     const Double_t GammaH = (Eg+Ep)/sqrt(s);
 
-    Double_t psi = stheta/(GammaH*(ctheta-sqrt(sqrbetah)));
+    const Double_t psi = stheta/(GammaH*(ctheta-sqrt(sqrbetah)));
 
     fAngle->SetParameter(0, sqrt(sqrbetae));
@@ -104,16 +102,16 @@
     const Double_t alpha = psi-acos(fAngle->GetRandom());
 
-    Double_t salpha = sin(alpha);
-    Double_t calpha = cos(alpha);
+    const Double_t salpha = sin(alpha);
+    const Double_t calpha = cos(alpha);
 
     const Double_t tphi = stheta/(Eg/Ep+ctheta); // tan(phi)
 
-    Double_t bb = sqrt(sqrbetah/sqrbetae);
+    const Double_t bb = sqrt(sqrbetah/sqrbetae);
 
-    Double_t s1 = calpha/GammaH;
-    Double_t s2 = tphi*s1 - salpha - bb;
+    const Double_t s1 = calpha/GammaH;
+    const Double_t s2 = tphi*s1 - salpha - bb;
 
-    Double_t tan1 = ((salpha+bb)*tphi+s1)/s2;
-    Double_t tan2 = ((salpha-bb)*tphi+s1)/s2;
+    const Double_t tan1 = ((salpha+bb)*tphi+s1)/s2;
+    const Double_t tan2 = ((salpha-bb)*tphi+s1)/s2;
 
     const Double_t E = (Eg+Ep)/2;;
Index: /trunk/WuerzburgSoft/Thomas/mphys/MPhoton.cc
===================================================================
--- /trunk/WuerzburgSoft/Thomas/mphys/MPhoton.cc	(revision 1362)
+++ /trunk/WuerzburgSoft/Thomas/mphys/MPhoton.cc	(revision 1363)
@@ -152,6 +152,6 @@
     Double_t val[2] = { Eg, z };
 
-    Double_t lolim = E0*E0 > 1e-8 ? E0*E0/Eg : 1e-8/Eg;
-    Double_t inf   = 3e-11; 
+    Double_t lolim = E0*E0/Eg;
+    Double_t inf   = Eg<1e6 ? 3e-11*(z+1) : 3e-12*(z+1);
 
     TF1 f("int2", Int2, lolim, inf, 2);
@@ -188,5 +188,5 @@
 }
 
-void MPhoton::DrawInteractionLength(Double_t z)
+void MPhoton::DrawInteractionLength(Double_t z, Double_t from, Double_t to)
 {
     if (!gPad)
@@ -195,5 +195,5 @@
         gPad->GetVirtCanvas()->cd(4);
 
-    TF1 f1("length", InteractionLength, 1e4, 1e11, 1);
+    TF1 f1("length", InteractionLength, from, to, 1);
     f1.SetParameter(0, z);
 
Index: /trunk/WuerzburgSoft/Thomas/mphys/MPhoton.h
===================================================================
--- /trunk/WuerzburgSoft/Thomas/mphys/MPhoton.h	(revision 1362)
+++ /trunk/WuerzburgSoft/Thomas/mphys/MPhoton.h	(revision 1363)
@@ -32,5 +32,5 @@
     // ----------------------------------------------------------------
 
-    static void DrawInteractionLength(Double_t z);
+    static void DrawInteractionLength(Double_t z, Double_t from=1e4, Double_t to=1e11);
     void DrawInteractionLength() const;
 
Index: /trunk/WuerzburgSoft/Thomas/mphys/anale.C
===================================================================
--- /trunk/WuerzburgSoft/Thomas/mphys/anale.C	(revision 1362)
+++ /trunk/WuerzburgSoft/Thomas/mphys/anale.C	(revision 1363)
@@ -4,6 +4,7 @@
 {
     TChain chain("Electrons");
-    chain.Add("cascade01.root");
-    chain.Add("cascade02.root");
+    chain.Add("cascade_100kpc_01.root");
+    chain.Add("cascade_100kpc_02.root");
+    chain.Add("cascade_100kpc_03.root");
 
     MElectron *e = new MElectron;
Index: /trunk/WuerzburgSoft/Thomas/mphys/analp.C
===================================================================
--- /trunk/WuerzburgSoft/Thomas/mphys/analp.C	(revision 1362)
+++ /trunk/WuerzburgSoft/Thomas/mphys/analp.C	(revision 1363)
@@ -4,6 +4,7 @@
 {
     TChain chain("Photons");
-    chain.Add("cascade01.root");
-    chain.Add("cascade02.root");
+    chain.Add("cascade_100kpc_01.root");
+    chain.Add("cascade_100kpc_02.root");
+    chain.Add("cascade_100kpc_03.root");
 
     MPhoton *p = new MPhoton;
Index: /trunk/WuerzburgSoft/Thomas/mphys/energyloss.C
===================================================================
--- /trunk/WuerzburgSoft/Thomas/mphys/energyloss.C	(revision 1362)
+++ /trunk/WuerzburgSoft/Thomas/mphys/energyloss.C	(revision 1363)
@@ -477,6 +477,49 @@
 // -------------------------------------------------------------------
 
+Double_t func(Double_t *x, Double_t *k)
+{
+    return MPhoton::Int2(x, k)*1e68;
+}
+
 void energyloss()
 {
+/*    Double_t E0    = 511e-6;                // [GeV]
+
+    Double_t Eg = 1e4; //3.6e4;
+    Double_t z  = 5;
+
+    Double_t val[2] = { Eg, z };
+
+    Double_t lolim = E0*E0/Eg;
+    Double_t inf   = Eg<1e6 ? 3e-11*z : 3e-12*z;
+
+    cout << Eg << " " << z << " " << lolim << " " << inf << endl;
+
+    TF1 f("int2", func, lolim, inf, 2);
+
+    Double_t int2 = f.Integral(lolim, inf, val); //[GeV^3 m^2]
+
+    f.SetParameter(0, Eg);
+    f.SetParameter(1, z);
+
+    cout << int2 << endl;
+
+    new TCanvas("ILPhoton", "Mean Interaction Length Photon");
+
+    gPad->SetLogx();
+//    gPad->SetLogy();
+    gPad->SetGrid();
+
+    f.SetLineWidth(1);
+    f.DrawCopy();
+
+    gPad->Modified();
+    gPad->Update();
+
+    //    return;
+  */  MPhoton p;
+    p.DrawInteractionLength();
+
+    return;
 //    EnergyLossRate();
 
Index: /trunk/WuerzburgSoft/Thomas/mphys/phys.C
===================================================================
--- /trunk/WuerzburgSoft/Thomas/mphys/phys.C	(revision 1362)
+++ /trunk/WuerzburgSoft/Thomas/mphys/phys.C	(revision 1363)
@@ -116,8 +116,8 @@
 void DoIt()
 {
-    //    Double_t R = 1798; // [kpc]
-    Double_t startz = 0.003; //MParticle::ZofR(&R);
-
-    //    cout << "R = " << R << "kpc" << endl;
+    Double_t R = 500; // [kpc]
+    Double_t startz = MParticle::ZofR(&R);
+
+    cout << "R = " << R << "kpc" << endl;
     cout << "Z = " << startz << endl;
 
@@ -131,5 +131,5 @@
     Double_t alpha = -2;
 
-    Int_t nbins = 18; //4*log10(hi/lo);
+    Int_t nbins = 21; //4*log10(hi/lo);
 
     TH1D hist;
@@ -153,11 +153,13 @@
 
     MBinning binspolx;
-    MBinning binspoly;
+    MBinning binspoly1;
+    MBinning binspoly2;
     binspolx.SetEdges(16, -180, 180);
-    binspoly.SetEdgesLog(20, 1e-10, 1e-3);
-    MH::SetBinning(&position, &binspolx, &binspoly);
-    MH::SetBinning(&angle,    &binspolx, &binspoly);
-    MH::SetBinning(&histpos,  &binspoly);
-    MH::SetBinning(&hista,    &binspoly);
+    binspoly1.SetEdges(20, 0, 5e-9);
+    binspoly2.SetEdges(20, 0, 2e-7);
+    MH::SetBinning(&angle,    &binspolx, &binspoly1);
+    MH::SetBinning(&position, &binspolx, &binspoly2);
+    MH::SetBinning(&hista,    &binspoly1);
+    MH::SetBinning(&histpos,  &binspoly2);
 
     hista.SetTitle("Particle Angle");
@@ -257,4 +259,5 @@
                     hista.Fill(p->GetTheta()*kRad2Deg, weight);
 
+                    p->SetIsPrimary(kFALSE);
                     T->Fill();
 
@@ -393,5 +396,5 @@
     position.SetYTitle("Pos: R [kpc]");
     position.DrawCopy("surf1pol");
-    gPad->SetLogy();
+    //gPad->SetLogy();
 
     MH::MakeDefCanvas();
@@ -399,5 +402,5 @@
     angle.SetYTitle("Angle: \\Theta [\\circ]");
     angle.DrawCopy("surf1pol");
-    gPad->SetLogy();
+    //gPad->SetLogy();
 
     MH::MakeDefCanvas();