Changeset 1374 for trunk/WuerzburgSoft/Thomas/mphys

Timestamp:

06/21/02 08:28:43 (23 years ago)

Author:

tbretz

Message:

*** empty log message ***

Location:

trunk/WuerzburgSoft/Thomas/mphys

Files:

• Changelog (modified) (1 diff)
• analp.C (modified) (12 diffs)
• phys.C (modified) (18 diffs)

trunk/WuerzburgSoft/Thomas/mphys/Changelog

@@ -1 +1 @@
                                                                   -*-*- END -*-*-
+ 2002/06/20: Thomas Bretz
+
+   * analp.C:
+     - clean up the code
+     - added the point-spread histograms from different directions
+ 
+   * mphys.C:
+     - removed the angular histograms
+     - added the plot spectral index
+     - added the timer to get a usable ubdated window
+
+
+
  2002/06/20: Thomas Bretz
 

trunk/WuerzburgSoft/Thomas/mphys/analp.C

@@ -2 +2 @@
 
 #include <float.h>
+
 Double_t GetMin(TH1 *h)
 {
@@ -13 +14 @@
 }
 
-void GetRange(TChain *chain, const char *name, Double_t &min, Double_t &max)
+void GetRange(TChain *chain, const char *name, Int_t &nbins, Double_t &min, Double_t &max, Double_t conv=1)
 {
     TString str("MPhoton.MParticle.");
@@ -32 +33 @@
         if (!leaf)
             return 0;
+
         Double_t val = leaf->GetValue();
 
@@ -42 +44 @@
             max = val;
     }
+
+    min *= conv;
+    max *= conv;
+
+    chain.GetPlayer();
+    TTreePlayer &play = *chain.GetPlayer();
+
+    Int_t n;
+    play.FindGoodLimits(nbins, n, min, max, kFALSE);
+    nbins = n;
 }
 
@@ -47 +59 @@
 {
     TChain chain("Photons");
-    // chain.Add("cascade_100kpc_a_1e6.root");
-    // chain.Add("cascade_100kpc_b_1e6.root");
-    // chain.Add("cascade_500kpc_0*.root");
-    chain.Add("cascade_500kpc_21d_0*.root");
-    // chain.Add("cascade_100kpc_0*.root");
-    // chain.Add("cascade_100kpc_14*.root");
-    // chain.Add("cascade_100kpc_0*.root");
-    // chain.Add("cascade_500kpc_21d_*.root");
-    // chain.Add("cascade_100kpc_4*.root");
-    // chain.Add("cascade_500kpc_21d_B1e-18_*.root");
-
+    chain.Add("cascade_500kpc_*.root");
+
+    MBinning binsx;
+    binsx.SetEdgesLog(21, 1e4, 1e11);
+
+    Double_t alpha = -1.8;
+    Double_t plot  =  1.8;
+
+    Int_t nbins = 20;
+
+    // ------------------------------
 
     Int_t n = chain.GetEntries();
 
-    cout << "Found " << n << " entries." << endl;
+    cout << endl << "Found " << n << " entries." << endl;
 
     if (n==0)
         return;
-
-    MBinning binsx;
-    binsx.SetEdgesLog(21, 1e4, 1e11);
 
     MBinning binspolx;
@@ -74 +83 @@
     binspolx.SetEdges(16, -180, 180);
 
-    const Double_t ls = 299792458;         // [m/s]
-    const Double_t ly = 3600.*24.*365.*ls; // [m/ly]
-    const Double_t pc = 1./3.258;          // [pc/ly]
-
-    Double_t max;
-    Double_t min;
-    Int_t nbins=20;
-
-    TTreePlayer &play = *chain.GetPlayer();
-
-    GetRange(&chain, "fTheta", min, max);
-    Double_t conv = kRad2Deg*60;
-    min *= conv;
-    max *= conv;
-    cout << "fTheta: " << min  << " < " << max  << " [']" << endl;
-    play.FindGoodLimits(10, nbins, min, max, kFALSE);
-    cout << "    " << nbins << ": " << min << " < " << max << endl;
+    Double_t min, max;
+    GetRange(&chain, "fTheta", nbins, min, max, kRad2Deg*60);
     binspola.SetEdges(nbins, min, max);
-
-    GetRange(&chain, "fR", min, max);
-    max/=4;
-    cout << "fR:     " << min  << " < " << max  << " [kpc]" << endl;
-    play.FindGoodLimits(10, nbins, min, max, kFALSE);
-    cout << "    " << nbins << ": " << min << " < " << max << endl;
+    cout << "fTheta, " << nbins << ": " << min << " ... " << max << " [']" << endl;
+
+    GetRange(&chain, "fR", nbins, min, max);
     binspolr.SetEdges(nbins, min, max);
+    cout << "fR,     " << nbins << ": " << min << " ... " << max << " [kpc]" << endl;
 
     TH1D h;
-    h.SetName("Photons");
+    TH1D prim;
     MH::SetBinning(&h, &binsx);
-
-    TH1D prim;
     MH::SetBinning(&prim, &binsx);
 
@@ -115 +105 @@
     chain.SetBranchAddress("MPhoton.", &p);
 
-    Double_t weight =  0;
-    Double_t alpha  = -2;
-    Double_t plot   =  2;
-    Double_t E;
+    chain.GetEntry(0);
+    if (!p->IsPrimary())
+        return;
+
+    Double_t z = p->GetZ();
+    Double_t R = MParticle::RofZ(&z);
+
+    cout << "Z = " << z << endl;
+    cout << "R = " << R << "kpc" << endl;
+
+    Double_t weight = 0;
+
     for (int i=0; i<n; i++)
     {
@@ -125 +123 @@
         Double_t Ep = p->GetEnergy();
 
-        if (i==0)
-            weight = pow(Ep, alpha);
-
         if (p->IsPrimary())
         {
             weight = pow(Ep, alpha);
-            prim.Fill(Ep, pow(Ep,plot) * weight);
-            E = Ep;
+            prim.Fill(Ep, pow(Ep, plot+alpha));
             continue;
         }
 
-        //if (Ep==E)
-        //    continue;
-
         h.Fill(Ep, pow(Ep,plot) * weight);
-        Double_t v = p->GetPhi()*kRad2Deg;
-        r.Fill(fmod(v+180, 360)-180, p->GetR(), weight);
-        v = p->GetPsi()*kRad2Deg;
-        a.Fill(fmod(v+180, 360)-180, p->GetTheta()*kRad2Deg*60, weight);
+
+        Double_t phi = fmod(p->GetPhi()*kRad2Deg+180, 360)-180;
+        Double_t psi = fmod(p->GetPsi()*kRad2Deg+180, 360)-180;
+
+        r.Fill(phi, p->GetR(),                 weight);
+        a.Fill(psi, p->GetTheta()*kRad2Deg*60, weight);
     }
 
@@ -163 +156 @@
     c->Divide(2,2);
 
+    gStyle->SetPalette(1, 0);
+
     c->cd(1);
-    gPad->SetLogy();
-    gPad->SetLogz();
-    gPad->SetTheta(70);
-    gPad->SetPhi(90);
-    h.SetTitle(Form(" E^{%.1f}  Spectra ", alpha));
-    h.SetXTitle("E\\gamma [GeV]");
-    h.SetYTitle(Form("E^{%.1f} * Counts", plot));
     r.SetTitle(" Distance from Observer ");
     r.GetXaxis()->SetLabelOffset(-0.015);
@@ -177 +165 @@
     r.SetXTitle("\\Phi [\\circ]");
     r.SetYTitle("R [kpc]");
-    r.DrawCopy("surf1pol");
+    TPad &pad = *gPad;
+    pad.Divide(2,2);
+    pad.cd(1);
+    gPad->SetLogy();
+    gPad->SetLogz();
+    gPad->SetTheta(0);
+    gPad->SetPhi(90);
+    g = r.DrawCopy("surf1pol");
+    pad.cd(2);
+    gPad->SetLogy();
+    gPad->SetLogz();
+    gPad->SetTheta(70);
+    gPad->SetPhi(90);
+    g->Draw("surf1pol");
+    pad.cd(3);
+    gPad->SetLogy();
+    gPad->SetLogz();
+    gPad->SetTheta(90);
+    gPad->SetPhi(90);
+    g->Draw("surf1pol");
+    pad.cd(4);
+    gPad->SetLogy();
+    gPad->SetLogz();
+    gPad->SetTheta(20);
+    gPad->SetPhi(90);
+    g->Draw("surf1pol");
 
     c->cd(2);
-    gPad->SetLogy();
-    gPad->SetLogz();
-    gPad->SetTheta(70);
-    gPad->SetPhi(90);
     a.SetTitle(" Hit Angle for Observer ");
     a.GetXaxis()->SetLabelOffset(-0.015);
@@ -190 +199 @@
     a.SetXTitle("\\Psi [\\circ]");
     a.SetYTitle("\\Theta [']");
-    a.DrawCopy("surf1pol");
+    TPad &pad = *gPad;
+    pad.Divide(2,2);
+    pad.cd(1);
+    gPad->SetLogy();
+    gPad->SetLogz();
+    gPad->SetTheta(0);
+    gPad->SetPhi(90);
+    g = a.DrawCopy("surf1pol");
+    pad.cd(2);
+    gPad->SetLogy();
+    gPad->SetLogz();
+    gPad->SetTheta(70);
+    gPad->SetPhi(90);
+    g->Draw("surf1pol");
+    pad.cd(3);
+    gPad->SetLogy();
+    gPad->SetLogz();
+    gPad->SetTheta(90);
+    gPad->SetPhi(90);
+    g->Draw("surf1pol");
+    pad.cd(4);
+    gPad->SetLogy();
+    gPad->SetLogz();
+    gPad->SetTheta(20);
+    gPad->SetPhi(90);
+    g->Draw("surf1pol");
 
     c->cd(3);
@@ -227 +261 @@
     gPad->SetLogx();
     gPad->SetLogy();
-    prim.SetFillStyle(0);
     h.SetFillStyle(0);
+    h.SetName("Photons");
+    h.SetTitle(Form(" E^{%.1f}  Spectra ", alpha));
+    h.SetXTitle("E\\gamma [GeV]");
+    h.SetYTitle(Form("E^{%.1f} * Counts", plot));
     h.GetXaxis()->SetLabelOffset(-0.015);
     h.GetXaxis()->SetTitleOffset(1.1);
-//    h.GetXaxis()->SetRangeUser(1e4, 1e9);
+    h.SetMarkerStyle(kMultiply);
+    prim.SetFillStyle(0);
     prim.SetMarkerStyle(kPlus);
-    h.SetMarkerStyle(kMultiply);
     prim.SetMarkerColor(kRed);
     prim.SetLineColor(kRed);
-
-    h.DrawCopy("P");
-    prim.DrawCopy("Psame");
-    prim.DrawCopy("Csame");
-    h.DrawCopy("Csame");
+    TH1 &g1=*h.DrawCopy("P");
+    TH1 &g2=*prim.DrawCopy("Psame");
+    g2.Draw("Csame");
+    g1.Draw("Csame");
  
     c->cd(2);

trunk/WuerzburgSoft/Thomas/mphys/phys.C

@@ -10 +10 @@
 #include <TFile.h>
 #include <TTree.h>
+#include <TTimer.h>
+#include <TStyle.h>
 #include <TBranch.h>
-#include <TStyle.h>
 #include <TRandom.h>
 #include <TCanvas.h>
@@ -121 +122 @@
     Double_t startz = MParticle::ZofR(&R);
 
-    const char *filename = "data/cascade_delme.root";
+    const char *filename = "cascade_delme.root";
 
     const Double_t B = 0;
@@ -128 +129 @@
     cout << "Z = " << startz << endl;
 
-    static TRandom rand(0);
     MPairProduction pair;
 
-    Double_t runtime = 60; // [s]
+    Double_t runtime = 5*60; // [s]
 
     Double_t lo = 1e4;
     Double_t hi = 1e11;
     Double_t alpha = -2;
+    Double_t plot  =  2;
 
     Int_t nbins = 21; //4*log10(hi/lo);
@@ -142 +143 @@
     TH1D histsrc;
 
-    hist.SetMinimum(lo*lo * pow(lo, alpha)/10);
-    histsrc.SetMinimum(lo*lo * pow(lo, alpha)/10);
+    hist.SetMinimum(pow(lo, plot+alpha)/10);
+    histsrc.SetMinimum(pow(lo, plot+alpha)/10);
 
     TList listg;
@@ -152 +153 @@
 
     gStyle->SetOptStat(10);
-
-    TH2D position;
-    TH2D angle;
-    TH1D histpos;
-    TH1D hista;
-
-    MBinning binspolx;
-    MBinning binspoly1;
-    MBinning binspoly2;
-    binspolx.SetEdges(16, -180, 180);
-    binspoly1.SetEdges(20, 0, 5e-6);
-    binspoly2.SetEdges(20, 0, 1e-5);
-    MH::SetBinning(&angle,    &binspolx, &binspoly1);
-    MH::SetBinning(&position, &binspolx, &binspoly2);
-    MH::SetBinning(&hista,    &binspoly1);
-    MH::SetBinning(&histpos,  &binspoly2);
-
-    hista.SetTitle("Particle Angle");
-    angle.SetTitle("Particle Angle");
-    histpos.SetTitle("Particle Position");
-    position.SetTitle("Particle Position");
 
     //
@@ -190 +170 @@
     MH::SetBinning(&histsrc, &bins);
 
-    MH::MakeDefCanvas();
+    TCanvas *c=MH::MakeDefCanvas();
 
     gPad->SetGrid();
@@ -196 +176 @@
     gPad->SetLogy();
 
-    histsrc.SetName("Spectrum");
-    histsrc.SetXTitle("E [GeV]");
-    histsrc.SetYTitle("E^{2}\\dotCounts");
-    histsrc.GetXaxis()->SetLabelOffset(-0.015);
-    histsrc.GetXaxis()->SetTitleOffset(1.1);
-
-    histsrc.Draw("P");
-    hist.Draw("Psame");
-    histsrc.Draw("Csame");
-    hist.Draw("Csame");
+    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");
 
     // ------------------------------
 
-    MPhoton *p4file  = NULL;
-    MElectron *e4file  = NULL;
-
+    void *ptr = NULL;
     TFile file(filename, "RECREATE", "Intergalactic cascade", 9);
-    TTree *T1 = new TTree("Photons",   "Photons from Cascade");
-    TTree *T2 = new TTree("Electrons", "Electrons in the Cascade");
-    TBranch *B1 = T1->Branch("MPhoton.",   "MPhoton",   &p4file, 32000);
-    TBranch *B2 = T2->Branch("MElectron.", "MElectron", &e4file, 32000);
+    TTree   *T1 = new TTree ("Photons",    "Photons from Cascade");
+    TTree   *T2 = new TTree ("Electrons",  "Electrons in the Cascade");
+    TBranch *B1 = T1->Branch("MPhoton.",   "MPhoton",   &ptr);
+    TBranch *B2 = T2->Branch("MElectron.", "MElectron", &ptr);
 
     // ------------------------------
+
+    TTimer timer("gSystem->ProcessEvents();", 250, kFALSE);
+    timer.TurnOn();
 
     TStopwatch clock;
     clock.Start();
-
-    Int_t timer = 0;
 
     Int_t n=0;
@@ -234 +213 @@
         Double_t E = GetEnergy(nbins, lo, hi);
         Double_t weight = pow(E, alpha);
-        histsrc.Fill(E, E*E * weight);
+        histsrc.Fill(E, pow(E, plot) * weight);
 
         MPhoton *gamma=new MPhoton(E, startz);
@@ -248 +227 @@
         Double_t Esum=0;
         Double_t Emis=0;
-
         while (1)
         {
@@ -265 +243 @@
                     cout << "!" << flush;
 
-                    hist.Fill(Eg, Eg*Eg*weight);
-                    position.Fill(p->GetPhi()*kRad2Deg, p->GetR(), weight);
-                    angle.Fill(p->GetPsi()*kRad2Deg, p->GetTheta()*kRad2Deg, weight);
-                    histpos.Fill(p->GetR(), weight);
-                    hista.Fill(p->GetTheta()*kRad2Deg, weight);
+                    hist.Fill(Eg, pow(Eg, plot)*weight);
 
                     p->SetIsPrimary(kFALSE);
@@ -286 +260 @@
                 phot.SetParameter(0, Eg);
                 phot.SetParameter(1, p->GetZ());
-                /*
-                if (phot.Integral(-log10(Eg)-8, -10.5, (Double_t*)NULL, 1e-2)==0)
-                {
-                    cout << "z" << flush;
-                    continue;
-                }
-                */
+
                 Double_t Ep    = pow(10, phot.GetRandom());
                 if (Ep==pow(10,0))
@@ -327 +295 @@
             while ((e=(MElectron*)Next()))
             {
-                e->SetIsPrimary();
+                e->SetIsPrimary(kTRUE);
                 T2->Fill();
                 e->SetIsPrimary(kFALSE);
@@ -340 +308 @@
                     {
                         cout << "!" << flush;
-                        Esum += e->GetEnergy();
-                        Emis += e->GetEnergy();
                         break;
                     }
 
                     // WRONG!
-                    Double_t theta = 0; //rand.Uniform(TMath::Pi()/2)+TMath::Pi()*3/4;
-                    MPhoton *p = e->DoInvCompton(theta);
+                    MPhoton *p = e->DoInvCompton(0);
 
                     T2->Fill();
@@ -371 +336 @@
                     {
                         cout << "T" << flush;
-                        Esum += e->GetEnergy();
-                        Emis += e->GetEnergy();
                         break;
                     }
@@ -379 +342 @@
                     {
                         cout << "E" << flush;
-                        Esum += e->GetEnergy();
-                        Emis += e->GetEnergy();
                         break;
                     }
                 }
+                Esum += e->GetEnergy();
+                Emis += e->GetEnergy();
             }
             liste.Delete();
@@ -389 +352 @@
         cout << " ----> " << Form("%3.1f %3.1e / %3.1f %3.1e", Emis/E, Emis, Esum/E, Esum) << endl;
 
-        Int_t now = (int)(TStopwatch::GetRealTime()- starttime)/5;
-        if (now!=timer || n<10)
-        {
-            histsrc.SetTitle(Form("E^%.1f, z=%f, T=%d'%d\", N=%d", alpha, startz,
-                                  (int)runtime/60, (int)runtime%60, hist.GetEntries()));
-            gPad->Modified();
-            gPad->Update();
-            timer = now;
-        }
+        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);
 
     }
@@ -405 +365 @@
     clock.Print();
 
+    timer.Stop();
+
     file.Write();
 
@@ -414 +376 @@
     gPad->Clear();
 
-    hist.SetTitle(Form("E^%.1f, z=%f, T=%d'%d\", N=%d", alpha, startz, (int)runtime/60, (int)runtime%60, n));
-    gPad->Modified();
-    gPad->Update();
-
-    hist.DrawCopy("P");
-    histsrc.DrawCopy("Psame");
-    hist.DrawCopy("Csame");
-    histsrc.DrawCopy("Csame");
-
-    gPad->SetGrid();
-    gPad->SetLogx();
-    gPad->SetLogy();
-
-    MH::MakeDefCanvas();
-    position.SetXTitle("Pos: \\Phi [\\circ]");
-    position.SetYTitle("Pos: R [kpc]");
-    position.DrawCopy("surf1pol");
-    //gPad->SetLogy();
-
-    MH::MakeDefCanvas();
-    angle.SetXTitle("Angle: \\Psi [\\circ]");
-    angle.SetYTitle("Angle: \\Theta [\\circ]");
-    angle.DrawCopy("surf1pol");
-    //gPad->SetLogy();
-
-    MH::MakeDefCanvas();
-    histpos.SetXTitle("R [kpc]");
-    histpos.SetYTitle("Counts");
-    histpos.DrawCopy();
-    gPad->SetLogx();
-
-    MH::MakeDefCanvas();
-    hista.SetXTitle("\\Theta [\\circ]");
-    hista.SetYTitle("Counts");
-    hista.DrawCopy();
-    gPad->SetLogx();
+    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");
 }