Changeset 1476 for trunk/WuerzburgSoft/Thomas/mphys/analp.C

Timestamp:

08/02/02 10:50:32 (22 years ago)

Author:

tbretz

Message:

*** empty log message ***

File:

• trunk/WuerzburgSoft/Thomas/mphys/analp.C (modified) (7 diffs)

trunk/WuerzburgSoft/Thomas/mphys/analp.C

@@ -20 +20 @@
 #include "MFit.h"
 #include "MPhoton.h"
+
+TH1D Smear(const TH1 &h, Double_t sigma)
+{
+    TH1D ret;
+
+    MH::SetBinning((TH1*)&ret, (TH1*)&h);
+
+    Int_t n = h.GetNbinsX();
+
+    for (int i=1; i<=n; i++)
+    {
+        Double_t xi = sqrt(h.GetBinLowEdge(i+1)*h.GetBinLowEdge(i));
+        for (int j=1; j<=n; j++)
+        {
+            Double_t xj = sqrt(h.GetBinLowEdge(j+1)*h.GetBinLowEdge(j));
+
+            Double_t dx = log10(xj/xi)/sigma;
+
+            Double_t gaus = exp(-dx*dx/2)/(sigma*sqrt(TMath::Pi()*2));
+
+            gaus *= h.GetBinContent(i);
+
+            ret.Fill(xj, gaus/sqrt(n));
+        }
+    }
+    return ret;
+}
 
 void GetRange(TChain *chain, const char *name, Int_t &nbins, Double_t &min, Double_t &max, Double_t conv=1, Bool_t zero=kTRUE)
@@ -110 +137 @@
     // -------------------------------------------------------------------
 
+    TString dir = "~/data/";
+
     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("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(dir+"cascade_0.1_18_1e2_1e5_B0_256_01.root");
+     chain.Add(dir+"cascade_0.1_18_1e2_1e5_B0_256_02.root");
+     chain.Add(dir+"cascade_0.1_18_1e2_1e5_B0_256_03.root");
+     chain.Add(dir+"cascade_0.1_18_1e2_1e5_B0_256_04.root");
+
+     chain.Add(dir+"cascade_0.1_18_1e2_1e5_B0_512_01.root");
+     chain.Add(dir+"cascade_0.1_18_1e2_1e5_B0_512_02.root");
+     chain.Add(dir+"cascade_0.1_18_1e2_1e5_B0_512_03.root");
+     chain.Add(dir+"cascade_0.1_18_1e2_1e5_B0_512_05.root");
      */
+
+    //chain.Add(dir+"cascade_0.5_18_1e2_1e5_B0_512_01.root");
+    //chain.Add(dir+"cascade_0.5_18_1e2_1e5_B1e-7_10Mpc_512_01.root");
+    //chain.Add(dir+"cascade_0.5_18_1e2_1e5_B1e-7_50Mpc_512_01.root");
+    chain.Add(dir+"cascade_0.5_18_1e2_1e5_B1e-7_100Mpc_512_01.root");
+
+    //chain.Add(dir+"cascade_0.1_18_1e2_1e5_B1e-6_10Mpc_512_01.root");
+    //chain.Add(dir+"cascade_0.1_18_1e2_1e5_B1e-6_50Mpc_512_01.root");
+    //chain.Add(dir+"cascade_0.1_18_1e2_1e5_B1e-7_100Mpc_512_01.root");
+
+    /*
+     chain.Add(dir+"cascade_0.03_18_1e2_1e5_B0_512_01.root");
+     chain.Add(dir+"cascade_0.03_18_1e2_1e5_B0_512_02.root");
+     chain.Add(dir+"cascade_0.03_18_1e2_1e5_B0_512_03.root");
+     chain.Add(dir+"cascade_0.03_18_1e2_1e5_B0_512_04.root");
+     */
+
+    //chain.Add(dir+"cascade_0.03_18_1e2_1e5_B1e-6_10Mpc_512_01.root");
+    //chain.Add(dir+"cascade_0.03_18_1e2_1e5_B1e-6_50Mpc_512_01.root");
+    //chain.Add(dir+"cascade_0.03_18_1e2_1e5_B1e-7_100Mpc_512_01.root");
 
     Int_t nbinse = 18;      // number of bins in your histogram
@@ -273 +314 @@
         {
             h3.Fill(Ep, pow(Ep,plot) * weight);
-            r3.Fill(phi, 0, weight);
-            a3.Fill(psi, 0, weight);
+            r3.Fill(phi, 0.0, weight);
+            a3.Fill(psi, 0.0, weight);
             isprim = kFALSE;
             continue;
@@ -456 +497 @@
     // ----------------------------------------------------------------------
 
-  //  delete gROOT->FindObject("Analysis Photons");
+    //  delete gROOT->FindObject("Analysis Photons");
     c = MH::MakeDefCanvas("Analysis Photons", "", 580, 870);
     c->Divide(1,2);
@@ -463 +504 @@
     gPad->SetLogx();
     gPad->SetLogy();
+    h3.SetTitle(Form(" E^{%.1f}  Spectra @ z=%.1e (R=%.0fkpc) ", alpha, z, R));
+    h3.SetXTitle("E\\gamma [GeV]");
+    h3.SetYTitle(Form("E^{%.1f} * Counts", plot));
+    h3.GetXaxis()->SetLabelOffset(-0.015);
+    h3.GetXaxis()->SetTitleOffset(1.1);
+    h3.SetFillStyle(0);
+    h3.SetName("PrimPhotons");
+    h3.SetMarkerStyle(kFullCircle);
+    h3.SetMarkerColor(kRed);
+    h3.SetMarkerSize(0.8);
+    h3.SetLineColor(kRed);
+    h3.DrawCopy("C"); //E1
+
+    TH1D hs = Smear(h, 0.2);
+    hs.SetLineColor(kGreen);
+    hs.SetFillStyle(0);
+    hs.DrawCopy("Csame");
+
     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.SetMarkerStyle(kPlus);
+    h.SetMarkerStyle(kFullCircle);
+    h.SetMarkerSize(0.8);
+    h.DrawCopy("Csame E4");
     prim.SetFillStyle(0);
     prim.SetMarkerStyle(kPlus);
     prim.SetMarkerColor(kRed);
     prim.SetLineColor(kRed);
-    TH1 &g1=*h.DrawCopy("P E4");
-    g2=prim.DrawCopy("Psame E0");
-    g2->Draw("Csame");
-    g1.Draw("Csame");
+    prim.SetMarkerSize(0.8);
+    prim.DrawCopy("Csame");
     h2.SetFillStyle(0);
     h2.SetName("SecPhotons");
-    h2.SetMarkerStyle(kMultiply);
+    h2.SetMarkerStyle(kFullCircle);
     h2.SetMarkerColor(kBlue);
+    h2.SetMarkerSize(0.8);
     h2.SetLineColor(kBlue);
-    TH1 &g3=*h2.DrawCopy("Psame E2");
-    g3.Draw("Csame");
-    h3.SetFillStyle(0);
-    h3.SetName("PrimPhotons");
-    h3.SetMarkerStyle(kMultiply);
-    h3.SetMarkerColor(kGreen);
-    h3.SetLineColor(kGreen);
-    TH1 &g4=*h3.DrawCopy("Psame E1");
-    g4.Draw("Csame");
+    h2.DrawCopy("Csame E4"); //E2
 
     draw1h1426(prim.GetBinContent(2),plot);
@@ -507 +554 @@
     div.SetXTitle("E\\gamma [GeV]");
     div.SetYTitle("Ratio");
+    div.SetMarkerStyle(kPlus);
     TH1 *gHist = div.DrawCopy("E4");
 
@@ -526 +574 @@
         }
     }
-    fit.SetParameter(0, "t", 750,  10, 1e5);
-    fit.SetParameter(1, "m", 0.5, 0.1,  10);
+    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;
+    cout << "Cutoff:  " << setprecision(2) << fit[0]/1e3 << " +- ";
+    cout << fit.GetParError(0)/1e3 << " TeV" << endl;
 
      // ----------------------------------------------------------------------