Index: /trunk/WuerzburgSoft/Thomas/mphys/Changelog
===================================================================
--- /trunk/WuerzburgSoft/Thomas/mphys/Changelog	(revision 1359)
+++ /trunk/WuerzburgSoft/Thomas/mphys/Changelog	(revision 1360)
@@ -1,3 +1,17 @@
                                                                   -*-*- END -*-*-
+ 2002/06/13: Thomas Bretz
+
+   * MParticle.h:
+     - Implemented Set/IsPrimary
+
+   * energyloss.C:
+     - changed output range for ZofR
+
+   * phys.C:
+     - changed energy randomizer
+     - added root file output
+
+
+
  2002/06/12: Thomas Bretz
 
Index: /trunk/WuerzburgSoft/Thomas/mphys/MParticle.h
===================================================================
--- /trunk/WuerzburgSoft/Thomas/mphys/MParticle.h	(revision 1359)
+++ /trunk/WuerzburgSoft/Thomas/mphys/MParticle.h	(revision 1360)
@@ -16,7 +16,8 @@
 public:
     typedef enum { kEGamma, kEElectron, kEPositron } ParticleType_t;
+    enum { kEIsPrimary = BIT(14) };
 
 private:
-    const ParticleType_t fPType;
+    const ParticleType_t fPType; //! Particle type
 
 protected:
@@ -45,4 +46,9 @@
     }
 
+    void SetIsPrimary(Bool_t is=kTRUE) { is ? SetBit(kEIsPrimary) : ResetBit(kEIsPrimary); }
+    Bool_t IsPrimary() const { return TestBit(kEIsPrimary); }
+
+    // ----------------------------------------------------------------
+
     virtual Double_t GetInteractionLength() const = 0;
 
Index: /trunk/WuerzburgSoft/Thomas/mphys/MPhoton.h
===================================================================
--- /trunk/WuerzburgSoft/Thomas/mphys/MPhoton.h	(revision 1359)
+++ /trunk/WuerzburgSoft/Thomas/mphys/MPhoton.h	(revision 1360)
@@ -9,5 +9,7 @@
 {
 public:
-    MPhoton(Double_t e=0, Double_t z=0) : MParticle(MParticle::kEGamma)
+
+    MPhoton(Double_t e=0, Double_t z=0)
+        : MParticle(MParticle::kEGamma)
     {
         fEnergy = e;
Index: /trunk/WuerzburgSoft/Thomas/mphys/energyloss.C
===================================================================
--- /trunk/WuerzburgSoft/Thomas/mphys/energyloss.C	(revision 1359)
+++ /trunk/WuerzburgSoft/Thomas/mphys/energyloss.C	(revision 1360)
@@ -449,5 +449,5 @@
     new TCanvas("RZ", "r and z");
 
-    TF1 f1("ZofR", MParticle::ZofR, 0, 4.5e5, 0);
+    TF1 f1("ZofR", MParticle::ZofR, 0, 7.1e6, 0);
     TF1 f2("RofZ", MParticle::RofZ, 0, 5,     0);
 
@@ -479,7 +479,5 @@
 void energyloss()
 {
-    EnergyLossRate();
-
-    MH::MakeDefCanvas();
+//    EnergyLossRate();
 
     DrawRZ();
Index: /trunk/WuerzburgSoft/Thomas/mphys/phys.C
===================================================================
--- /trunk/WuerzburgSoft/Thomas/mphys/phys.C	(revision 1359)
+++ /trunk/WuerzburgSoft/Thomas/mphys/phys.C	(revision 1360)
@@ -8,4 +8,7 @@
 #include <TH2.h>
 #include <TList.h>
+#include <TFile.h>
+#include <TTree.h>
+#include <TBranch.h>
 #include <TStyle.h>
 #include <TRandom.h>
@@ -64,4 +67,51 @@
 }
 
+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*bin);
+
+    //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()
 {
@@ -75,5 +125,5 @@
     MPairProduction pair;
 
-    Double_t runtime = 30*60; // [s]
+    Double_t runtime = 15*60; // [s]
 
     Double_t lo = 1e4; 
@@ -81,8 +131,5 @@
     Double_t alpha = -2;
 
-    Double_t nbins = 24; //4*log10(hi/lo);
-
-    TF1 src ("PrimarySpectrum", PrimSpect, log10(lo), log10(hi), 1); // 10GeV-1000PeV
-    src.SetParameter(0, 0);
+    Int_t nbins = 18; //4*log10(hi/lo);
 
     TH1D hist;
@@ -150,6 +197,14 @@
     hist.Draw("Csame");
 
-    ofstream fsrc("source.dat");
-    ofstream fcas("cascade.dat");
+    // ------------------------------
+
+    MPhoton *p4file  = NULL;
+    MElectron *e4file  = NULL;
+
+    TFile file("cascade.root", "RECREATE", "Intergalactic cascade", 9);
+    TTree *T  = new TTree("Photons",   "Photons from Cascade");
+    TTree *T2 = new TTree("Electrons", "Electrons in the Cascade");
+    TBranch *B =T ->Branch("MPhoton.",   "MPhoton",   &p4file, 32000);
+    TBranch *B2=T2->Branch("MElectron.", "MElectron", &e4file, 32000);
 
     // ------------------------------
@@ -163,18 +218,20 @@
     Double_t starttime = TStopwatch::GetRealTime(); // s
     while (TStopwatch::GetRealTime()<starttime+runtime)
+        for (int i=0; i<nbins; i++)
     {
         n++;
 
-        Double_t E = pow(10, src.GetRandom());
+        Double_t E = GetEnergy(nbins, lo, hi);
         Double_t weight = pow(E, alpha);
-
         histsrc.Fill(E, E*E * weight);
 
-        fsrc << E << endl;
+        MPhoton *gamma=new MPhoton(E, startz);
+        gamma->SetIsPrimary();
+        listg.Add(gamma);
+
+        B->SetAddress(&gamma);
+        T->Fill();
 
         cout << "--> " << n << ". " << Form("%d: %3.1e", (int)(starttime+runtime-TStopwatch::GetRealTime()), E) << flush;
-
-        MPhoton *gamma=new MPhoton(E, startz);
-        listg.Add(gamma);
 
         while (1)
@@ -185,4 +242,5 @@
             TIter NextP(&listg);
             MPhoton *p = NULL;
+            B->SetAddress(&p);
             while ((p=(MPhoton*)NextP()))
             {
@@ -199,5 +257,6 @@
                     hista.Fill(p->GetTheta()*kRad2Deg, weight);
 
-                    fcas << Eg << endl;
+                    T->Fill();
+
                     delete listg.Remove(p);
                     continue;
@@ -242,6 +301,11 @@
             TIter Next(&liste);
             MElectron *e = NULL;
+            B2->SetAddress(&e);
             while ((e=(MElectron*)Next()))
             {
+                e->SetIsPrimary();
+                T2->Fill();
+                e->SetIsPrimary(kFALSE);
+
                 if (e->GetEnergy()<lo)
                     continue;
@@ -259,4 +323,6 @@
                     Double_t theta = rand.Uniform(TMath::Pi()*2);
                     MPhoton *p = e->DoInvCompton(theta);
+
+                    T2->Fill();
 
                     if (fabs(p->GetTheta()*3437)<1 &&  // < 1min
@@ -300,4 +366,6 @@
     clock.Print();
 
+    file.Write();
+
     cout << "Created " << n << " gammas from source with E^" << alpha << endl;
     cout << "Processed " << Form("%.1f", n/(TStopwatch::GetRealTime()-starttime)) << " gammas/sec." << endl;