Changeset 1449 for trunk/WuerzburgSoft

Timestamp:

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

Author:

tbretz

Message:

*** empty log message ***

Location:

trunk/WuerzburgSoft/Thomas/mphys

Files:

• Changelog (modified) (1 diff)
• MCascade.cc (added)
• MCascade.h (added)
• MElectron.cc (modified) (13 diffs)
• MElectron.h (modified) (1 diff)
• MFit.cc (added)
• MFit.h (added)
• MPairProduction.cc (modified) (1 diff)
• MParticle.cc (modified) (9 diffs)
• MParticle.h (modified) (2 diffs)
• MPhoton.cc (modified) (1 diff)
• MPhoton.h (modified) (2 diffs)
• Makefile (modified) (2 diffs)
• PhysLinkDef.h (modified) (1 diff)
• analp.C (modified) (20 diffs)
• phys.C (modified) (1 diff)

trunk/WuerzburgSoft/Thomas/mphys/Changelog

@@ -1 +1 @@
                                                                   -*-*- END -*-*-
+ 2002/06/29: Thomas Bretz
+
+   * MElectron.[h,cc]:
+     - if the photon energy is to small return a straight line in
+       Compton
+     - changed lolimit of p_e integration -0.5
+     - changed number of calculated points to 50 for fP, fQ
+     - changed usage of TRandom to gRandom
+     - changed some TVector and TMatrix objects to static
+     - added new copy constructors
+
+   * MPairProduction.cc:
+     - use the new copy constructor of MElectron
+
+   * MParticle.[h,cc]:
+     - added new copy constructors
+     - Added a TAxis object in Planck
+     - replaced usage of TRandom by gRandom
+     - replaced TVector, TMatrix by TVector3, TRotation
+     - fixed a bug in SetNewPosition (fX was not incremented if the
+       particle reached the observers plain)
+     - incremented version-number
+
+   * MPhoton.[h,cc]:
+     - added Fill member function
+     - added new copy constructors
+
+   * MFit.[h,cc], MCascade.[h,cc]:
+     - added
+   
+   * phys.C:
+     - moved code to MCascade
+
+   * Makefile, PhysLinkDef.h:
+     - added MFit, MCascade
+
+   * analp.C:
+     - simplified
+     - added fit for cutoff
+
+
+
  2002/06/26: Thomas Bretz
 

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);

trunk/WuerzburgSoft/Thomas/mphys/MElectron.h

@@ -16 +16 @@
         fZ      = z;
     }
+
+    MElectron(MParticle &p, Bool_t type=kFALSE) : MParticle(p, type?MParticle::kEPositron:MParticle::kEElectron) { }
+    MElectron(MParticle &p, Double_t e, Bool_t type=kFALSE) : MParticle(p, e, type?MParticle::kEPositron:MParticle::kEElectron) {}
 
     void operator=(MParticle &p) { MParticle::operator=(p); }

trunk/WuerzburgSoft/Thomas/mphys/MPairProduction.cc

@@ -119 +119 @@
     // cout << " {" << f << "," << E << "," << atan(tan1) << "," << atan(-tan2) << "} " << flush;
 
-    MElectron &p0 = *new MElectron(E*(1.-f), gamma->GetZ(), kTRUE);
-    MElectron &p1 = *new MElectron(E*(1.+f), gamma->GetZ(), kFALSE);
-    p0 = *gamma; // Set Position and direction
-    p1 = *gamma; // Set Position and direction
+    MElectron &p0 = *new MElectron(*gamma, E*(1.-f), kTRUE);
+    MElectron &p1 = *new MElectron(*gamma, E*(1.+f), kFALSE);
 
-    static TRandom rand(0);
-    Double_t rnd = rand.Uniform(TMath::Pi() * 2);
+    Double_t rnd = gRandom->Uniform(TMath::Pi() * 2);
     p0.SetNewDirection(atan(+tan1), rnd);
     p1.SetNewDirection(atan(-tan2), rnd+TMath::Pi());

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);
 }

trunk/WuerzburgSoft/Thomas/mphys/MParticle.h

@@ -37 +37 @@
 public:
     MParticle(ParticleType_t t=kENone, const char *name=NULL, const char *title=NULL);
+    MParticle(MParticle &p, ParticleType_t t=kENone) : fPType(t) { *this = p; }
+    MParticle(MParticle &p, Double_t e, ParticleType_t t=kENone) : fPType(t) { *this = p; fEnergy = e; }
 
     void InitRandom();
@@ -89 +91 @@
     Double_t GetPsi() const    { return fPsi; }
 
-    ClassDef(MParticle, 1) // Container which holds hostograms for the Hillas parameters
+    ClassDef(MParticle, 2) // Container which holds hostograms for the Hillas parameters
 };
 

trunk/WuerzburgSoft/Thomas/mphys/MPhoton.cc

@@ -200 +200 @@
     DrawInteractionLength(fZ);
 }
+
+void MPhoton::Fill(TH1 &h, Double_t idx, Double_t w) const
+{
+    h.Fill(fEnergy, pow(fEnergy, idx)*w);
+}

trunk/WuerzburgSoft/Thomas/mphys/MPhoton.h

@@ -5 +5 @@
 #include "MParticle.h"
 #endif
+
+class TH1;
 
 class MPhoton : public MParticle
@@ -17 +19 @@
     }
 
+        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); }
+
+    void Fill(TH1 &h, Double_t idx, Double_t w) const;
 
     // ----------------------------------------------------------------

trunk/WuerzburgSoft/Thomas/mphys/Makefile

@@ -22 +22 @@
 #  connect the include files defined in the config.mk file
 #
-INCLUDES = -I. -I../mbase
+INCLUDES = -I. -I../mbase -I../mhist
 
 #------------------------------------------------------------------------------
@@ -29 +29 @@
 
 SRCFILES = MParticle.cc \
+           MPairProduction.cc \
            MPhoton.cc \
-           MElectron.cc \
-           MPairProduction.cc
+           MFit.cc \
+           MCascade.cc \
+           MElectron.cc
 
 SRCS    = $(SRCFILES)

trunk/WuerzburgSoft/Thomas/mphys/PhysLinkDef.h

@@ -8 +8 @@
 #pragma link C++ class MElectron+;
 #pragma link C++ class MPhoton+;
+#pragma link C++ class MCascade+;
+#pragma link C++ class MFit+;
 #pragma link C++ class MPairProduction+;
 

trunk/WuerzburgSoft/Thomas/mphys/analp.C

@@ -2 +2 @@
 
 #include <float.h>
-
-Double_t GetMin(TH1 *h)
-{
-    Double_t min = FLT_MAX;
-    for (int i=1; i<=h->GetXaxis()->GetNbins(); i++)
-    {
-        if (h->GetBinContent(i)<min && h->GetBinContent(i)!=0)
-            min = h->GetBinContent(i);
-    }
-    return min;
-}
-
-void GetRange(TChain *chain, const char *name, Int_t &nbins, Double_t &min, Double_t &max, Double_t conv=1)
+#include <iomanip.h>
+
+#include <TH1.h>
+#include <TH2.h>
+#include <TLine.h>
+#include <TLeaf.h>
+#include <TChain.h>
+#include <TGaxis.h>
+#include <TStyle.h>
+#include <TCanvas.h>
+#include <TPolyLine.h>
+#include <TPolyMarker.h>
+
+#include "MH.h"
+#include "MBinning.h"
+
+#include "MFit.h"
+#include "MPhoton.h"
+
+void GetRange(TChain *chain, const char *name, Int_t &nbins, Double_t &min, Double_t &max, Double_t conv=1, Bool_t zero=kTRUE)
 {
     TString str("MPhoton.MParticle.");
+    str += name;
 
     MPhoton *p = new MPhoton;
 
     chain->SetBranchAddress("MPhoton.", &p);
+    chain->SetBranchStatus("*", 0);
+    chain->SetBranchStatus(str, 1);
 
     min =  FLT_MAX;
     max = -FLT_MAX;
 
-    Int_t n = chain->GetEntries();
+    Int_t  numt = chain->GetTreeNumber();
+    TLeaf *leaf = chain->GetLeaf(str);
+
+    if (!leaf && numt>=0)
+        return;
+
+    Stat_t n = chain->GetEntries();
     for (int i=0; i<n; i++)
     {
         chain->GetEntry(i);
 
-        TLeaf *leaf = chain->GetLeaf(str+name);
-        if (!leaf)
-            return 0;
+        if (numt != chain->GetTreeNumber())
+        {
+            if (!(leaf = chain->GetLeaf(str)))
+                return;
+
+            numt = chain->GetTreeNumber();
+        }
 
         Double_t val = leaf->GetValue();
@@ -39 +59 @@
             continue;
 
-        if (val < min)
+        if (val < min && (zero || val!=0))
             min = val;
         if (val > max)
@@ -48 +68 @@
     max *= conv;
 
-    chain.GetPlayer();
-    TTreePlayer &play = *chain.GetPlayer();
-
-    Int_t n;
-    play.FindGoodLimits(nbins, n, min, max, kFALSE);
-    nbins = n;
+    Int_t newn=0;
+    MH::FindGoodLimits(nbins, newn, min, max, kFALSE);
+    nbins = newn;
 }
 
 void draw1h1426(Double_t E, Double_t plot=1)
 {
-    plot=2;
-    Double_t level = log10(E);
-    level -= log10(5.73e-0*pow( 0.39e3, plot));
-
-    cout << "Level: " << level << endl;
-
-    TPolyMarker *m = new TPolyMarker(7);
-    /*
-     m->SetPoint(0,  0.28e3, 1.86e-0*pow( 0.28e3, plot)*level);
-     m->SetPoint(1,  0.39e3, 5.73e-0*pow( 0.39e3, plot)*level);
-     m->SetPoint(2,  0.80e3, 1.22e-0*pow( 0.80e3, plot)*level);
-     m->SetPoint(3,  1.55e3, 5.10e-2*pow( 1.55e3, plot)*level);
-     m->SetPoint(4,  2.82e3, 1.50e-2*pow( 2.82e3, plot)*level);
-     m->SetPoint(5,  5.33e3, 7.70e-3*pow( 5.33e3, plot)*level);
-     m->SetPoint(6, 10.00e3, 1.20e-4*pow(10.00e3, plot)*level);
-     */
-    m->SetPoint(0, log10( 0.28e3), log10(1.86e-0*pow( 0.28e3, plot))+level);
-    m->SetPoint(1, log10( 0.39e3), log10(5.73e-0*pow( 0.39e3, plot))+level);
-    m->SetPoint(2, log10( 0.80e3), log10(1.22e-0*pow( 0.80e3, plot))+level);
-    m->SetPoint(3, log10( 1.55e3), log10(5.10e-2*pow( 1.55e3, plot))+level);
-    m->SetPoint(4, log10( 2.82e3), log10(1.50e-2*pow( 2.82e3, plot))+level);
-    m->SetPoint(5, log10( 5.33e3), log10(7.70e-3*pow( 5.33e3, plot))+level);
-    m->SetPoint(6, log10(10.00e3), log10(1.20e-4*pow(10.00e3, plot))+level);
+    plot += 1;
+
+    const Double_t level = E / (5.73e-0*pow( 0.39e3, plot));
+
+    TPolyLine   *l = new TPolyLine(6);
+    TPolyMarker *m = new TPolyMarker(6);
+
+    //  m->SetPoint(0, log10( 0.28e3), log10(1.86e-0*pow( 0.28e3, plot)*level));
+    m->SetPoint(0, log10( 0.39e3)-0.3, log10(5.73e-0*pow( 0.39e3, plot)*level));
+    m->SetPoint(1, log10( 0.80e3)-0.3, log10(1.22e-0*pow( 0.80e3, plot)*level));
+    m->SetPoint(2, log10( 1.55e3)-0.3, log10(5.10e-2*pow( 1.55e3, plot)*level));
+    m->SetPoint(3, log10( 2.82e3)-0.3, log10(1.50e-2*pow( 2.82e3, plot)*level));
+    m->SetPoint(4, log10( 5.33e3)-0.3, log10(7.70e-3*pow( 5.33e3, plot)*level));
+    m->SetPoint(5, log10(10.00e3)-0.3, log10(1.20e-4*pow(10.00e3, plot)*level));
     m->SetMarkerStyle(kOpenStar);
-//    m->SetMarkerSize(1);
+    m->SetMarkerColor(kMagenta);
+    m->SetBit(kCanDelete);
     m->Draw();
+
+    //  l->SetPoint(0, log10( 0.28e3), log10(1.86e-0*pow( 0.28e3, plot)*level));
+    l->SetPoint(0, log10( 0.39e3)-0.3, log10(5.73e-0*pow( 0.39e3, plot)*level));
+    l->SetPoint(1, log10( 0.80e3)-0.3, log10(1.22e-0*pow( 0.80e3, plot)*level));
+    l->SetPoint(2, log10( 1.55e3)-0.3, log10(5.10e-2*pow( 1.55e3, plot)*level));
+    l->SetPoint(3, log10( 2.82e3)-0.3, log10(1.50e-2*pow( 2.82e3, plot)*level));
+    l->SetPoint(4, log10( 5.33e3)-0.3, log10(7.70e-3*pow( 5.33e3, plot)*level));
+    l->SetPoint(5, log10(10.00e3)-0.3, log10(1.20e-4*pow(10.00e3, plot)*level));
+    l->SetLineColor(kMagenta);
+    l->SetBit(kCanDelete);
+    l->Draw();
 }
 
 void analp()
 {
+    // -------------------------------------------------------------------
+
     TChain chain("Photons");
+    chain.Add("cascade_0.1_18_1e2_1e5_B0_256_01.root");
+    chain.Add("cascade_0.1_18_1e2_1e5_B0_256_02.root");
+    chain.Add("cascade_0.1_18_1e2_1e5_B0_256_03.root");
+    chain.Add("cascade_0.1_18_1e2_1e5_B0_256_04.root");
+    chain.Add("cascade_0.1_18_1e2_1e5_B0_512_01.root");
+    chain.Add("cascade_0.1_18_1e2_1e5_B0_512_03.root");
+    chain.Add("cascade_0.1_18_1e2_1e5_B0_512_02.root");
     /*
-     chain.Add("delme3.root");
-     chain.Add("delme3B.root");
-     chain.Add("delme3C.root");
-     chain.Add("delme3D.root");
-     chain.Add("delme3E.root");
+     chain.Add("cascade_0.03_18_1e2_1e5_B0_256_3.root");
+     chain.Add("cascade_0.03_18_1e2_1e5_B0_256_4.root");
+     chain.Add("cascade_0.03_18_1e2_1e5_B0_256_5.root");
+     chain.Add("cascade_0.03_18_1e2_1e5_B0_256_6.root");
+     chain.Add("cascade_0.03_18_1e2_1e5_B0_256_7.root");
+     chain.Add("cascade_0.03_18_1e2_1e5_B0_256_8.root");
+     chain.Add("cascade_0.03_18_1e2_1e5_B0_256_9.root");
+
+     chain.Add("cascade_0.03_18_1e2_1e5_B0_512_01.root");
+     chain.Add("cascade_0.03_18_1e2_1e5_B0_512_02.root");
      */
-    chain.Add("delme3H.root");
-
-    Int_t nbins =  24;
-    Double_t lo = 1e2;
-    Double_t hi = 1e6;
-
-    MBinning binsx;
-    binsx.SetEdgesLog(nbins, lo, hi);
-
-    // ------------------------------
-
-    Double_t cutoff = 1e12;
-
-    Double_t alpha = -1;
-    Double_t plot  =  1;
-
-    Int_t nbins = 16;
-
-    // ------------------------------
-
-    Int_t n = chain.GetEntries();
+
+    Int_t nbinse = 18;      // number of bins in your histogram
+    Double_t lo  = 1e2/1e3; // lower limit of the energy histogram
+    Double_t hi  = 1e5;     // upper limit of the energy histogram
+
+    Double_t minE = 1e2;    // minimum value produced in the simulation
+
+    Double_t cutoff = 1e12; // arbitrary energy cutoff (throw away the events)
+
+    Double_t alpha = -1;    // alpha of the energy spectrum (-1 is E^-2)
+    Double_t plot  =  1;    // plot of the enrgy spectrum (1 is E^2)
+
+    Int_t nbins = 16;       // number of bins you want to have in psi/phi
+
+    // -------------------------------------------------------------------
+
+    MBinning binspolx;
+    binspolx.SetEdges(16, -180, 180);
+
+    // -------------------------------------------------------------------
+
+    Stat_t n = chain.GetEntries();
 
     cout << endl << "Found " << n << " entries." << endl;
@@ -123 +158 @@
         return;
 
-    MBinning binspolx;
-    MBinning binspola;
-    MBinning binspolr;
-    binspolx.SetEdges(16, -180, 180);
+    // -----------------------------------
+
+    MBinning bins;
+
+    cout << "Determin min/max..." << endl;
 
     Double_t min, max;
     GetRange(&chain, "fTheta", nbins, min, max, kRad2Deg*60);
-    binspola.SetEdges(nbins, min, max*1.1);
     cout << "fTheta, " << nbins << ": " << min << " ... " << max << " [']" << endl;
 
+    bins.SetEdges(nbins, min, max*1.1);
+
+    TH2D a, a2, a3;
+    MH::SetBinning(&a,  &binspolx, &bins);
+    MH::SetBinning(&a2, &binspolx, &bins);
+    MH::SetBinning(&a3, &binspolx, &bins);
+
+    // -----------------------------------
+
     GetRange(&chain, "fR", nbins, min, max);
-    binspolr.SetEdges(nbins, min, max*1.1);
     cout << "fR,     " << nbins << ": " << min << " ... " << max << " [kpc]" << endl;
 
-    TH1D h, h2, h3;
-    TH1D prim;
-    MH::SetBinning(&h, &binsx);
-    MH::SetBinning(&h2, &binsx);
-    MH::SetBinning(&h3, &binsx);
+    bins.SetEdges(nbins, min, max*1.1);
+
+    TH2D r, r2, r3;
+    MH::SetBinning(&r,  &binspolx, &bins);
+    MH::SetBinning(&r2, &binspolx, &bins);
+    MH::SetBinning(&r3, &binspolx, &bins);
+
+    // -----------------------------------
+
+    MBinning binsx;
+    binsx.SetEdgesLog(nbinse, lo, hi);
+
+    TH1D prim, h, h2, h3;
+    MH::SetBinning(&h,    &binsx);
+    MH::SetBinning(&h2,   &binsx);
+    MH::SetBinning(&h3,   &binsx);
     MH::SetBinning(&prim, &binsx);
-
-    TH2D r, r2, r3;
-    TH2D a, a2, a3;
-    MH::SetBinning(&a, &binspolx, &binspola);
-    MH::SetBinning(&r, &binspolx, &binspolr);
-    MH::SetBinning(&a2, &binspolx, &binspola);
-    MH::SetBinning(&r2, &binspolx, &binspolr);
-    MH::SetBinning(&a3, &binspolx, &binspola);
-    MH::SetBinning(&r3, &binspolx, &binspolr);
+//    prim.Sumw2();
+    h.Sumw2();
+    h2.Sumw2();
+    h3.Sumw2();
+
+    // -----------------------------------
 
     MPhoton *p = new MPhoton;
     chain.SetBranchAddress("MPhoton.", &p);
+    chain.SetBranchStatus("*", 1);
 
     chain.GetEntry(0);
@@ -165 +217 @@
     cout << "Z = " << z << endl;
     cout << "R = " << R << "kpc" << endl;
+
+    // -----------------------------------
+
+    const Double_t skpc = 3600.*24.*365.*3.258*1e3;  // [s/kpc]
+    GetRange(&chain, "fX", nbins, min, max, 1, kFALSE);
+    min *= skpc;
+    max *= skpc;
+    if (min<1e-2)
+        min=1e-2;
+    bins.SetEdgesLog(nbins, min, max);
+    cout << "dX,     " << nbins << ": " << min << " ... " << max << " [s]" << endl;
+
+    TH1D t;
+    MH::SetBinning(&t, &bins);
+    t.Sumw2();
+
+    // ----------------------------------------------------------------------
+
+    chain.SetBranchAddress("MPhoton.", &p);
+    chain.SetBranchStatus("*", 1);
+
+    cout << "Filling: " << nbinse << " bins @ " << minE << " - " << hi << "... " << flush;
 
     Double_t weight = 0;
@@ -193 +267 @@
             continue;
 
-        h.Fill(Ep, pow(Ep,plot) * weight);
-
         Double_t phi = fmod(p->GetPhi()*kRad2Deg+180, 360)-180;
         Double_t psi = fmod(p->GetPsi()*kRad2Deg+180, 360)-180;
-
-        r.Fill(phi, isprim?0:p->GetR(),                 weight);
-        a.Fill(psi, isprim?0:p->GetTheta()*kRad2Deg*60, weight);
 
         if (isprim)
         {
             h3.Fill(Ep, pow(Ep,plot) * weight);
-            r3.Fill(phi, 0,                 weight);
+            r3.Fill(phi, 0, weight);
             a3.Fill(psi, 0, weight);
             isprim = kFALSE;
@@ -211 +280 @@
 
         h2.Fill(Ep, pow(Ep,plot) * weight);
-        r2.Fill(phi, p->GetR(),                 weight);
+        r2.Fill(phi, p->GetR(), weight);
         a2.Fill(psi, p->GetTheta()*kRad2Deg*60, weight);
+
+        Double_t tm1 = sqrt(R*R+p->GetR()*p->GetR());
+        Double_t tm2 = (R+p->GetX())/tm1 - 1;
+        t.Fill(p->GetX()<=0 ? 0 : R*tm2*skpc, weight);
     }
 
+    // ----------------------------------------------------------------------
+
     delete p;
 
-    TH1 *g=r.ProjectionX();
-    cout << "Mean fPhi: " << g->GetMean() << " +- " << g->GetRMS() << endl;
-    delete g;
-    g=a.ProjectionX();
-    cout << "Mean fPsi: " << g->GetMean() << " +- " << g->GetRMS() <<  endl;
-    delete g;
+    cout << "Done." << endl;
+
+    Double_t N = prim.GetBinContent(nbinse);
+    Double_t Nerr = prim.GetBinError(nbinse);
+    for (int i=1; i<=nbinse; i++)
+    {
+        Double_t E = prim.GetXaxis()->GetBinCenter(i);
+        if (E>minE)
+            continue;
+
+        Double_t E1 = prim.GetXaxis()->GetBinLowEdge(i+1);
+        Double_t E0 = prim.GetXaxis()->GetBinLowEdge(i);
+
+        E = sqrt(E1*E0);
+
+        prim.SetBinContent(i, N);
+        h3.SetBinContent(i, N);
+        prim.SetBinError(i, Nerr);
+        h3.SetBinError(i, Nerr);
+    }
+
+    h.Add(&h2, &h3);
+    r.Add(&r2, &r3);
+    a.Add(&a2, &a3);
+
+    if (alpha==-plot)
+    {
+        cout << "Observed Energy: " << h.Integral() << "GeV" << endl;
+        cout << "Emitted Energy:  " << prim.Integral() << "GeV" << endl;
+        cout << "Ratio: " << h.Integral()/prim.Integral() << endl;
+    }
+
+    cout << "Mean fPhi: " << r.GetMean() << " +- " << r.GetRMS() << endl;
+    cout << "Mean fPsi: " << a.GetMean() << " +- " << a.GetRMS() <<  endl;
 
     gStyle->SetOptStat(10);
+    gStyle->SetPalette(1, 0);
 
     TLine line;
 
+    // ----------------------------------------------------------------------
+
 //    delete gROOT->FindObject("Analysis Arrival");
-    c = MH::MakeDefCanvas("Analysis Arrival", "");
+    TCanvas *c = MH::MakeDefCanvas("Analysis Arrival", "");
     c->Divide(2,2);
-
-    gStyle->SetPalette(1, 0);
 
     c->cd(1);
@@ -241 +345 @@
     r.SetXTitle("\\Phi [\\circ]");
     r.SetYTitle("R [kpc]");
-    TPad &pad = *gPad;
+    TVirtualPad &pad = *gPad;
     pad.Divide(2,2);
     pad.cd(1);
@@ -248 +352 @@
     gPad->SetTheta(0);
     gPad->SetPhi(90);
-    g = r.DrawCopy("surf1pol");
+    TH1 *g = r.DrawCopy("surf1pol");
     pad.cd(2);
     gPad->SetLogy();
@@ -275 +379 @@
     a.SetXTitle("\\Psi [\\circ]");
     a.SetYTitle("\\Theta [']");
-    TPad &pad = *gPad;
-    pad.Divide(2,2);
-    pad.cd(1);
+    TVirtualPad &pad2 = *gPad;
+    pad2.Divide(2,2);
+    pad2.cd(1);
     gPad->SetLogy();
     gPad->SetLogz();
@@ -283 +387 @@
     gPad->SetPhi(90);
     g = a.DrawCopy("surf1pol");
-    pad.cd(2);
+    pad2.cd(2);
     gPad->SetLogy();
     gPad->SetLogz();
@@ -289 +393 @@
     gPad->SetPhi(90);
     g->Draw("surf1pol");
-    pad.cd(3);
+    pad2.cd(3);
     gPad->SetLogy();
     gPad->SetLogz();
@@ -295 +399 @@
     gPad->SetPhi(90);
     g->Draw("surf1pol");
-    pad.cd(4);
+    pad2.cd(4);
     gPad->SetLogy();
     gPad->SetLogz();
@@ -307 +411 @@
     g->SetDirectory(NULL);
     TH1 *g2=r2.ProjectionY();
-    g->SetMinimum(GetMin(g2)/2);
+    g->SetMinimum(MH::GetMinimumGT(*g2)/2);
     g->SetBit(kCanDelete);
     g->SetTitle(" Hit Observers Plain ");
@@ -329 +433 @@
     gPad->SetLogy();
     g=a.ProjectionY();
-    g->SetMinimum(GetMin(g)/2);
+    g->SetMinimum(MH::GetMinimumGT(*g)/2);
     g->SetDirectory(NULL);
     g->SetBit(kCanDelete);
@@ -350 +454 @@
     g->Draw("same");
 
-  //  delete gROOT->FindObject("Analysis Photons");
-    TCanvas *c = MH::MakeDefCanvas("Analysis Photons", "", 580, 870);
+    // ----------------------------------------------------------------------
+
+  //  delete gROOT->FindObject("Analysis Photons");
+    c = MH::MakeDefCanvas("Analysis Photons", "", 580, 870);
     c->Divide(1,2);
 
@@ -364 +470 @@
     h.GetXaxis()->SetLabelOffset(-0.015);
     h.GetXaxis()->SetTitleOffset(1.1);
-    h.SetMarkerStyle(kMultiply);
+    h.SetMarkerStyle(kPlus);
     prim.SetFillStyle(0);
     prim.SetMarkerStyle(kPlus);
     prim.SetMarkerColor(kRed);
     prim.SetLineColor(kRed);
-    TH1 &g1=*h.DrawCopy("P");
-    g2=*prim.DrawCopy("Psame");
+    TH1 &g1=*h.DrawCopy("P E4");
+    g2=prim.DrawCopy("Psame E0");
     g2->Draw("Csame");
     g1.Draw("Csame");
@@ -378 +484 @@
     h2.SetMarkerColor(kBlue);
     h2.SetLineColor(kBlue);
-    TH1 &g3=*h2.DrawCopy("Psame");
+    TH1 &g3=*h2.DrawCopy("Psame E2");
     g3.Draw("Csame");
     h3.SetFillStyle(0);
@@ -385 +491 @@
     h3.SetMarkerColor(kGreen);
     h3.SetLineColor(kGreen);
-    TH1 &g4=*h3.DrawCopy("Psame");
+    TH1 &g4=*h3.DrawCopy("Psame E1");
     g4.Draw("Csame");
 
@@ -394 +500 @@
     TH1D div;
     MH::SetBinning(&div, &binsx);
+    div.Sumw2();
     div.Divide(&h, &prim);
-    div.SetTitle("Spectrum / Primary Spectrum");
+    div.SetTitle(" Spectrum / Primary Spectrum ");
     div.GetXaxis()->SetLabelOffset(-0.015);
     div.GetXaxis()->SetTitleOffset(1.1);
     div.SetXTitle("E\\gamma [GeV]");
     div.SetYTitle("Ratio");
+    TH1 *gHist = div.DrawCopy("E4");
+
+    line.SetLineWidth(1);
+    line.SetLineColor(kGreen);
+    line.DrawLine(log10(lo), exp(-1), log10(hi), exp(-1));
     line.SetLineColor(kBlue);
-    line.SetLineWidth(1);
-    div.DrawCopy();
     line.DrawLine(log10(lo), 1, log10(hi), 1);
+
+    MFit fit("exp(-1)-[1]*log10(x/[0])");
+    for (int i=0; i<gHist->GetEntries(); i++)
+    {
+        if (gHist->GetBinContent(i+1)<exp(-1))
+        {
+            fit.SetRange(gHist->GetBinCenter(i-2), gHist->GetBinCenter(i+2));
+            cout << "Fitting from " << gHist->GetBinLowEdge(i-1) << " to ";
+            cout << gHist->GetBinLowEdge(i+1) << endl;
+            break;
+        }
+    }
+    fit.SetParameter(0, "t", 750,  10, 1e5);
+    fit.SetParameter(1, "m", 0.5, 0.1,  10);
+    fit.FitLog(gHist);
+    fit.Print();
+    fit.DrawCopy("same");
+
+    cout << "Cutoff:  " << setprecision(2) << fit[0]/1e3 << "TeV +- ";
+    cout << fit.GetParError(0) << endl;
+
+     // ----------------------------------------------------------------------
+
+    c = MH::MakeDefCanvas("Time Analysis", "", 580, 435);
+    gPad->SetLogx();
+    //    gPad->SetLogy();
+    t.SetName("Times");
+    t.SetTitle(" Arrival time distribution ");
+    t.SetXTitle("\\Delta  t [s]");
+    t.GetXaxis()->SetLabelOffset(-0.015);
+    t.GetXaxis()->SetTitleOffset(1.1);
+    t.DrawCopy("E4");
+    line.DrawLine(log10(1),           0, log10(1),           t.GetMaximum()*1.05);
+    line.DrawLine(log10(3600),        0, log10(3600),        t.GetMaximum()*1.05);
+    line.DrawLine(log10(3600*24),     0, log10(3600*24),     t.GetMaximum()*1.05);
+    line.DrawLine(log10(3600*24*365), 0, log10(3600*24*365), t.GetMaximum()*1.05);
 }

trunk/WuerzburgSoft/Thomas/mphys/phys.C

@@ -1 @@
-#include <iostream.h>
+#include "MCascade.h"
 
-#include <fstream.h>
+void phys()
+{
+    MCascade cascade;
 
-#include <TStopwatch.h>
-#include <TF1.h>
-#include <TH1.h>
-#include <TH2.h>
-#include <TList.h>
-#include <TFile.h>
-#include <TTree.h>
-#include <TTimer.h>
-#include <TStyle.h>
-#include <TBranch.h>
-#include <TRandom.h>
-#include <TCanvas.h>
+    cascade.SetSourceZ(0.1);             // Readshift of source
+    cascade.SetB(0/*1e-6*/);             // [G] mean magnetic field
+    cascade.SetBradius(50);              // [Mpc]
+    cascade.SetRuntime(8*60);            // [min] maximum time to run the simulation
+    cascade.SetEnergyBins(18, 1e2, 1e5); // [GeV]
+    cascade.SetMaxInvCompton(512);       // [#] maximum number of inv. Compton (-1 means infinite)
+    cascade.SetRatioInvCompton(0);       // [%] allowed Emis of electron (0 means disabled)
+    cascade.SetSpectralIndex(-1);        // -1 means a E^2 plot
+    cascade.SetDisplayIndex(1);          //  1 means a E^-2 spectrum
 
-#include "mbase/MParContainer.h"
-#include "mphys/MPhoton.h"
-#include "mphys/MElectron.h"
-#include "mphys/MPairProduction.h"
-#include "mhist/MBinning.h"
-#include "mhist/MH.h"
-
-// 2.96
-// 2.87
-// 2.73
-Double_t PrimSpect(Double_t *x, Double_t *k)
-{
-    return pow(pow(10, x[0]), k[0]);
+    //
+    // Run the simulation: filename, eventdisply (on/off)
+    //
+    //cascade.Run("cascade_0.03_18_1e2_1e5_B1e-6_50Mpc_256_1.root", kTRUE);
+    cascade.Run("cascade_0.1_18_1e2_1e5_B0_512_02.root", kFALSE);
 }
 
-Double_t PhotonSpect(Double_t *x, Double_t *k=NULL)
-{
-    Double_t Ep = pow(10, x[0]);
-
-    Double_t res = MPhoton::Int2(&Ep, k);
-    return res*1e55; //65/k[0];
-
-    //return MPhoton::Planck(&Ep, &k[1]);
-}
-
-Double_t Sbar_sigmas(Double_t *x, Double_t *k)
-{
-    Double_t sbar = pow(10, x[0]);
-
-    Double_t s = 1./(sbar*4);
-
-    Double_t sigma = MPhoton::Sigma_gg(&s);
-
-    return sigma*sbar*1e28;
-}
-
-Double_t RandomTheta(Double_t Eg, Double_t Ep)
-{
-    Double_t E0 = 511e-6; // [GeV]
-
-    Double_t f = Eg/E0*Ep/E0;
-
-    if (f<1)
-        return 0;
-
-    TF1 func("RndTheta", Sbar_sigmas, 0, log10(f), 0);
-
-    func.SetNpx(50);
-
-    Double_t sbar  = pow(10, func.GetRandom());
-    Double_t theta = acos(1.-sbar*2/f);
-
-    return theta;
-}
-
-Double_t GetEnergy(Int_t n, Double_t lo, Double_t hi)
-{
-    static int bin=0;
-    static TRandom rnd(0);
-
-    Double_t w = log10(hi/lo)/n;
-
-    Double_t E = lo*pow(10, rnd.Uniform(w) + w*(n-bin-1));
-
-    //cout << endl << w << " " << n << " " << w*bin << " " << E << endl;
-
-    if (++bin==n)
-        bin=0;
-
-    return E;
-
-    /*
-    // TF1 src ("PrimarySpectrum", PrimSpect, log10(lo), log10(hi), 1); // 10GeV-1000PeV
-    // src.SetParameter(0, 0);
-    // Double_t E = pow(10, src.GetRandom());
-
-    //
-    // 0: 11111111111111|22222222222222   2   2^0 + 1    1
-    // 1: 11111111|22222222222|33333333   3   2^1 + 1    2
-    // 2: 11111|22222|33333|44444|55555   5   2^2 + 1    7
-    // 3:  |     | |   | |  |  |   |                    15
-    //
-
-    static int run=0;
-    static int num=1;
-
-    Double_t w = log10(hi/lo)/((1<<run) + 1);
-
-    Double_t E = lo*pow(10, num*w);
-
-    if (num == (1<<run))
-    {
-        run++;
-        num=1;
-    }
-    else
-        num++;
-
-        return E;
-        */
-}
-
-void DoIt()
-{
-    // a -->  5
-    // b --> 10
-    // c -->  2
-
-    // delme, delmeB starting integration at lolim
-    // delme2, delme2B starting integration at -log10(Eg)-8
-    // delme3,  lolim, 24, 1e2-1e6,   2x inv.Compton
-    // delme3B, lolim, 24, 1e2-1e6,   4x inv.Compton
-    // delme3C, lolim, 24, 1e2-1e6,   8x inv.Compton
-    // delme3D, lolim, 24, 1e2-1e6,  16x inv.Compton
-    // delme3E, lolim, 24, 1e2-1e6,  32x inv.Compton
-
-    // delme3H, lolim, 24, 1e2-1e6, 256x inv.Compton 8*60
-
-    Double_t startz = 0.03; //MParticle::ZofR(&R);
-    Double_t R      = MParticle::RofZ(&startz); // [kpc]
-
-    const char *filename = "cascade_0.03_18_1e2_1e5_B0_256_1.root";
-    //const char *filename = "test.root";
-
-    const Double_t B = 0;//1e-6*1e-4; // [T=1e4G] mean magnetic field
-
-    Double_t runtime = 45*60;     // [s] maximum time to run the simulation
-
-    Int_t nbins = 18;      // number of bins produced in energy spectrum
-
-    Double_t lo = 1e2;     // lower limit of displayed spectrum
-    Double_t hi = 1e5;     // upper limit of spectrum (cutoff)
-
-    Int_t counter = 256;   // maximum number of inv. Compton (-1 means infinite)
-
-    Double_t alpha = -1;   // -1 means a E^2 plot
-    Double_t plot  =  1;   //  1 means a E^-2 spectrum
-
-    Double_t bubbler = R-1e3*10; // [kpc]
-    Double_t bubblez = MParticle::ZofR(&bubbler);
-
-    // --------------------------------------------------------------------
-
-    cout << "R = " << R << "kpc" << endl;
-    cout << "Z = " << startz << endl;
-
-    cout << "Setting up: Histograms... " << flush;
-
-    MPairProduction pair;
-
-    TH1D hist;
-    TH1D histsrc;
-
-    hist.SetMinimum(pow(lo, plot+alpha)/100);
-    histsrc.SetMinimum(pow(lo, plot+alpha)/100);
-
-    TList listg;
-    TList liste;
-
-    listg.SetOwner();
-    liste.SetOwner();
-
-    gStyle->SetOptStat(10);
-
-    //
-    // Don't change the order!!!
-    //
-    hist.SetFillStyle(0);
-    hist.SetMarkerStyle(kPlus);
-    histsrc.SetFillStyle(0);
-    histsrc.SetMarkerStyle(kMultiply);
-    histsrc.SetMarkerColor(kRed);
-    histsrc.SetLineColor(kRed);
-
-    MBinning bins;
-    bins.SetEdgesLog(nbins, lo, hi);
-
-    MH::SetBinning(&hist,    &bins);
-    MH::SetBinning(&histsrc, &bins);
-
-    TCanvas *c=MH::MakeDefCanvas();
-
-    gPad->SetGrid();
-    gPad->SetLogx();
-    gPad->SetLogy();
-
-    hist.SetName("Spectrum");
-    hist.SetXTitle("E [GeV]");
-    hist.SetYTitle(Form("E^{%.1f} Counts", plot));
-    hist.GetXaxis()->SetLabelOffset(-0.015);
-    hist.GetXaxis()->SetTitleOffset(1.1);
-
-    hist.Draw("P");
-    histsrc.Draw("Psame");
-    histsrc.Draw("same");
-    hist.Draw("same");
-
-    // ------------------------------
-
-    cout << "Output File... " << flush;
-
-    TFile file(filename, "RECREATE", "Intergalactic cascade", 9);
-    cout << "Trees... " << flush;
-    TTree   *T1 = new TTree ("Photons",    "Photons from Cascade");
-    TTree   *T2 = new TTree ("Electrons",  "Electrons in the Cascade");
-    cout << "Branches... " << flush;
-    MPhoton dummyp;
-    void *ptr = &dummyp;
-    TBranch *B1 = T1->Branch("MPhoton.",   "MPhoton",   &ptr);
-    MPhoton dummye;
-    ptr = &dummye;
-    TBranch *B2 = T2->Branch("MElectron.", "MElectron", &ptr);
-
-    // ------------------------------
-
-    cout << "Timers... " << flush;
-
-    TTimer timer("gSystem->ProcessEvents();", 250, kFALSE);
-    timer.TurnOn();
-
-    TStopwatch clock;
-    clock.Start();
-
-    cout << "Done. " << endl;
-
-    Int_t n=0;
-    Double_t starttime = TStopwatch::GetRealTime(); // s
-    while (TStopwatch::GetRealTime()<starttime+runtime)
-        for (int i=0; i<nbins; i++)
-    {
-        n++;
-
-        Double_t E = GetEnergy(nbins, lo, hi);
-        Double_t weight = pow(E, alpha);
-        histsrc.Fill(E, pow(E, plot) * weight);
-
-        MPhoton *gamma=new MPhoton(E, startz);
-        gamma->SetIsPrimary();
-        gamma->SetSrcR(R);
-        gamma->InitRandom();
-        listg.Add(gamma);
-
-        B1->SetAddress(&gamma);
-        T1->Fill();
-
-        cout << "--> " << n << ". " << Form("%d: %3.1e", (int)(starttime+runtime-TStopwatch::GetRealTime()), E) << flush;
-
-        Double_t Esum=0;
-        Double_t Emis=0;
-        while (1)
-        {
-            if (listg.GetSize()!=0)
-                cout << " |P" << flush;
-
-            TIter NextP(&listg);
-            MPhoton *p = NULL;
-            B1->SetAddress(&p);
-            while ((p=(MPhoton*)NextP()))
-            {
-                cout << ":" << flush;
-                Double_t Eg = p->GetEnergy();
-
-                Double_t E0 = 511e-6;
-                Double_t z  = p->GetZ();
-                Double_t lolim = E0*E0/Eg;
-                Double_t inf   = (Eg<1e6 ? 3e-11*(z+1) : 3e-12*(z+1));
-                if (Eg<5e4)
-                    inf = 3e-11*(z+1)*pow(10, 9.4-log10(Eg)*2);
-
-                //TF1 phot("PhotonSpectrum", PhotonSpect, -log10(Eg)-8, -10.5, 2);
-                //Double_t E0 = 511e-6; // [GeV]
-                TF1 phot("PhotonSpectrum", PhotonSpect, log10(lolim), log10(inf), 2);
-                phot.SetParameter(0, Eg);
-                while (1)
-                {
-                    if (!p->SetNewPosition())
-                    {
-                        cout << "!" << flush;
-
-                        hist.Fill(Eg, pow(Eg, plot)*weight);
-
-                        p->SetIsPrimary(kFALSE);
-                        T1->Fill();
-
-                        Esum += Eg;
-
-                        break;
-                    }
-
-                    //
-                    // Sample phtoton from background and interaction angle
-                    //
-                    phot.SetParameter(1, p->GetZ());
-                    Double_t pe = phot.GetRandom();
-                    if (pe==0)
-                    {
-                        cout << "z" << flush;
-                        continue;
-                    }
-
-                    Double_t Ep = pow(10, pe);
-                    Double_t theta = RandomTheta(Eg, Ep);
-                    if (theta==0)
-                    {
-                        cout << "t" << flush;
-                        continue;
-                    }
-
-                    if (!pair.Process(p, Ep, theta, &liste))
-                    {
-                        cout << "0" << flush;
-                        continue;
-                    }
-
-                    cout << "." << flush;
-                    break;
-                }
-                delete listg.Remove(p);
-            }
-
-            if (liste.GetSize()==0 && listg.GetSize()==0)
-                break;
-
-            if (liste.GetSize()!=0)
-                cout << " |E" << flush;
-
-            TIter Next(&liste);
-            MElectron *e = NULL;
-            B2->SetAddress(&e);
-            while ((e=(MElectron*)Next()))
-            {
-                e->SetIsPrimary(kTRUE);
-                T2->Fill();
-                e->SetIsPrimary(kFALSE);
-
-                Double_t Ee = e->GetEnergy();
-
-                if (Ee<lo)
-                    continue;
-
-                cout << ":" << flush;
-                int test = counter<0 ? -1 : 0;
-                while (test<0 ? true : (test++<counter))
-                {
-                    if (!e->SetNewPositionB(e->GetZ()>bubblez ? B : 0))
-                    {
-                        cout << "!" << flush;
-                        break;
-                    }
-
-                    // WRONG!
-                    MPhoton *p = e->DoInvCompton(0);
-
-                    T2->Fill();
-
-                    cout << "." << flush;
-                    listg.Add(p);
-                    /*
-                    if (fabs(p->GetTheta()*3437)<60 &&  // < 60min
-                        p->GetEnergy()>lo)
-                    {
-                        cout << "." << flush;
-                        listg.Add(p);
-                    }
-                    else
-                    {
-                        if (p->GetEnergy()<E*pow(10, alpha)/5 || p->GetEnergy()<=lo)
-                            cout << "e" << flush;
-                        else
-                            cout << "t" << flush; // ignored
-                        delete p;
-                        Esum += p->GetEnergy();
-                        Emis += p->GetEnergy();
-                    }
-                    */
-                    if (fabs(e->GetTheta()*3437)>60)  // < 60min
-                    {
-                        cout << "T" << flush;
-                        break;
-                    }
-
-                    if (e->GetEnergy()<Ee*1e-3) // <2e3
-                    {
-                        cout << "E" << flush;
-                        break;
-                    }
-                }
-                Esum += e->GetEnergy();
-                Emis += e->GetEnergy();
-            }
-            liste.Delete();
-        }
-        cout << " ----> " << Form("%3.1f %3.1e / %3.1f %3.1e", Emis/E, Emis, Esum/E, Esum) << endl;
-
-        hist.SetTitle(Form("E^{%.1f}  z=%f  T=%d'%d\"  N=%d", alpha, startz,
-                              (int)runtime/60, (int)runtime%60, n));
-
-        timer.Stop();
-        c->Update();
-        timer.Start(250);
-
-    }
-    cout << endl;
-
-    clock.Stop();
-    clock.Print();
-
-    timer.Stop();
-
-    file.Write();
-
-    cout << "Wrote: " << filename << endl;
-    cout << "Created " << n << " gammas from source with E^" << alpha << endl;
-    cout << "Processing time: " << Form("%.1f", (TStopwatch::GetRealTime()-starttime)/n) << " sec/gamma." << endl;
-
-    // ------------------------------
-
-    gPad->Clear();
-
-    hist.SetTitle(Form("E^{%.1f}  z=%f  T=%d'%d\"  N=%d", alpha, startz, (int)runtime/60, (int)runtime%60, n));
-
-    TH1 &h1 = *hist.DrawCopy("P");
-    TH1 &h2 = *histsrc.DrawCopy("Psame");
-    h2.Draw("Csame");
-    h1.Draw("Csame");
-}
-
-// -------------------------------------------------------------------
-void phys()
-{
-    /*
-    Double_t Eg = 10;
-    Double_t E0 = 511e-6;
-    Double_t z  = 0.03;
-    Double_t lolim = E0*E0/Eg;
-    Double_t inf   = (Eg<1e6 ? 3e-11*(z+1) : 3e-12*(z+1));
-    if (Eg<5e4)
-        inf = 3e-11*(z+1)*pow(10, 9.4-log10(Eg)*2);
-    TF1 phot("PhotonSpectrum", PhotonSpect, log10(lolim), log10(inf), 2);
-    phot.SetParameter(0, Eg);
-    phot.SetParameter(1, z);
-
-    Double_t val[2] = {Eg,z};
-    cout << phot.Integral(log10(lolim), log10(inf), val, 1e-2) << endl;
-
-    cout << lolim << " " << inf << endl;
-
-    phot.DrawCopy();
-    return;
-*/
-
-    DoIt();
-
-    /*
-     Double_t E = 1e10;
-     TF1 phot("PhotonSpectrum", PhotonSpect, -log10(E)-8, -10.5, 2);
-     phot.SetParameter(0, E);
-     phot.SetParameter(1, 5);
-     phot.DrawCopy();
-     return;
-     */
-
-    /*
-     Double_t Eg = 1e5;
-
-     Double_t E0    = 511e-6;                  // [GeV]
-     Double_t lolim = E0*E0/Eg;
-     Double_t inf   = 4e-12; //E0*E0/Eg * sqrt(Eg);
-
-     // 1e5:   5e-13, 4e-12
-     // 1e6:   5e-14, 2e-12
-     // 1e8:   5e-16, 1e-12
-     // 1e10:  5e-18, 1e-12
-
-     // 1e5:   2.6e-12, 4e-12
-     // 1e6:   2.6e-13, 2e-12
-     // 1e8:   2.6e-15, 1e-12
-     // 1e10:  1.6e-17, 1e-12
-
-     TF1 f("int2", MPhoton::Int2, lolim, inf, 2);
-     f.SetParameter(0, Eg);
-     f.SetParameter(1, 0);
-
-     MH::MakeDefCanvas();
-     gPad->SetLogx();
-     f.DrawCopy();
-    */
-}
-