Changeset 7147 for trunk/MagicSoft/Mars/mjobs

Timestamp:

06/13/05 09:41:57 (20 years ago)

Author:

tbretz

Message:

*** empty log message ***

File:

• trunk/MagicSoft/Mars/mjobs/MJSpectrum.cc (modified) (13 diffs)

trunk/MagicSoft/Mars/mjobs/MJSpectrum.cc

@@ -38 +38 @@
 #include <TFile.h>
 #include <TChain.h>
+#include <TLatex.h>
 #include <TCanvas.h>
 #include <TObjArray.h>
@@ -84 +85 @@
 MJSpectrum::MJSpectrum(const char *name, const char *title)
     : fCut0(0),fCut1(0), fCut2(0), fCut3(0), fEstimateEnergy(0),
-    fRefill(kFALSE), fSimpleMode(kTRUE), fRawMc(kFALSE)
+    fRefill(kFALSE), fSimpleMode(kTRUE), fRawMc(kFALSE),
+    fNoThetaWeights(kFALSE)
 {
     fName  = name  ? name  : "MJSpectrum";
@@ -336 +338 @@
     // Calculate the Probability
     temp1.Divide(&temp2);
-    temp1.Scale(1./temp1.GetMaximum());
+    temp1.Scale(fNoThetaWeights ? 1./temp1.GetMaximum() : 1./temp1.Integral());
 
     // Some cosmetics: Name, Axis, etc.
@@ -348 +350 @@
 }
 
+// --------------------------------------------------------------------------
+//
+// Display the final theta distribution.
+//
 void MJSpectrum::DisplayResult(const TH2D &h2) const
 {
@@ -376 +382 @@
 }
 
+// --------------------------------------------------------------------------
+//
+// Fills the excess histogram (vs E-est) from the events stored in the
+// ganymed result file and therefor estimates the energy.
+//
+// The resulting histogram excess-vs-energy ist copied into h2.
+//
 Bool_t MJSpectrum::Refill(MParList &plist, TH1D &h2) const
 {
@@ -531 +544 @@
 }
 
+// --------------------------------------------------------------------------
+//
+// Calculate the final spectrum from:
+//  - collection area
+//  - excess
+//  - correction coefficients
+//  - ontime
+// and display it
+//
 TArrayD MJSpectrum::DisplaySpectrum(MHCollectionArea &area, TH1D &excess, MHEnergyEst &hest, Double_t ontime) const
 {
@@ -573 +595 @@
     weights.DrawCopy();
 
-    //spectrum.Divide(&weights);
-    spectrum.Multiply(&weights);
-    spectrum.SetNameTitle("Flux", "N/TeVsm^{2} versus Energy (after unfolding)");
+    //spectrum.Multiply(&weights);
+    spectrum.SetNameTitle("Flux", "Spectrum (after unfolding)");
+    spectrum.SetBit(TH1::kNoStats);
 
     for (int i=0; i<excess.GetNbinsX(); i++)
     {
+        spectrum.SetBinContent(i+1, spectrum.GetBinContent(i+1)*weights.GetBinContent(i+1));
+        spectrum.SetBinError(i+1,   spectrum.GetBinError(i+1)  *weights.GetBinContent(i+1));
+
         spectrum.SetBinContent(i+1, spectrum.GetBinContent(i+1)/spectrum.GetBinWidth(i+1)*1000);
         spectrum.SetBinError(i+1,   spectrum.GetBinError(i+1)/  spectrum.GetBinWidth(i+1)*1000);
@@ -589 +614 @@
     gPad->SetGridx();
     gPad->SetGridy();
+    spectrum.SetMinimum(1e-12);
     spectrum.SetXTitle("E [GeV]");
     spectrum.SetYTitle("N/TeVsm^{2}");
     spectrum.DrawCopy();
 
-    TF1 f("f", "[1]*(x/1e3)^[0]", 50, 3e4);
+    TF1 f("f", "[1]*(x/1e3)^[0]", 10, 3e4);
     f.SetParameter(0, -2.87);
     f.SetParameter(1, 1.9e-6);
     f.SetLineColor(kGreen);
-    spectrum.Fit(&f, "NIM", "", 55, 2e4);
+    spectrum.Fit(&f, "NIM", "", 100, 4000);
     f.DrawCopy("same");
 
-    /*
-     TString str;
-     str += "(1.68#pm0.15)10^{-7}";
-     str += "(\\frac{E}{TeV})^{-2.59#pm0.06}";
-     str += "\\frac{ph}{TeVm^{2}s}";
-
-     TLatex tex;
-     tex.DrawLatex(2e2, 7e-5, str);
-     */
+    const Double_t p0 = f.GetParameter(0);
+    const Double_t p1 = f.GetParameter(1);
+
+    const Double_t e0 = f.GetParError(0);
+    const Double_t e1 = f.GetParError(1);
+
+    const Int_t    np  = TMath::Nint(TMath::Floor(TMath::Log10(p1)));
+    const Double_t exp = TMath::Power(10, np);
+
+    TString str;
+    str += Form("(%.2f#pm%.2f)10^{%d}", p1/exp, e1/exp, np);
+    str += Form("(\\frac{E}{TeV})^{-%.2f#pm%.2f}", p0, e0);
+    str += "\\frac{ph}{TeVm^{2}s}";
+
+    TLatex tex;
+    tex.SetTextSize(0.045);
+    tex.SetBit(TLatex::kTextNDC);
+    tex.DrawLatex(0.45, 0.935, str);
+
+    str = Form("\\chi^{2}/NDF=%.2f", f.GetChisquare()/f.GetNDF());
+    tex.DrawLatex(0.70, 0.83, str);
 
     TArrayD res(2);
     res[0] = f.GetParameter(0);
     res[1] = f.GetParameter(1);
+
+
+    // Plot other spectra  from Whipple
+    f.SetParameter(0, -2.45);
+    f.SetParameter(1, 3.3e-7);
+    f.SetRange(300, 8000);
+    f.SetLineColor(kBlack);
+    f.SetLineStyle(kDashed);
+    f.DrawCopy("same");
+
+    // Plot other spectra  from Cangaroo
+    f.SetParameter(0, -2.53);
+    f.SetParameter(1, 2.0e-7);
+    f.SetRange(7000, 50000);
+    f.SetLineColor(kBlack);
+    f.SetLineStyle(kDashed);
+    f.DrawCopy("same");
+
+    // Plot other spectra  from Robert
+    f.SetParameter(0, -2.59);
+    f.SetParameter(1, 2.58e-7);
+    f.SetRange(150, 1500);
+    f.SetLineColor(kBlack);
+    f.SetLineStyle(kDashed);
+    f.DrawCopy("same");
+
+    // Plot other spectra  from HEGRA
+    f.SetParameter(0, -2.61);
+    f.SetParameter(1, 2.7e-7);
+    f.SetRange(1000, 20000);
+    f.SetLineColor(kBlack);
+    f.SetLineStyle(kDashed);
+    f.DrawCopy("same");
+
     return res;
 }
 
+// --------------------------------------------------------------------------
+//
+// Scale some image parameter plots using the scale factor and plot them
+// together with the corresponding MC histograms.
+// Called from DisplaySize
+//
 Bool_t MJSpectrum::PlotSame(MStatusArray &arr, MParList &plist, const char *name, const char *tab, const char *plot, Double_t scale) const
 {
@@ -664 +742 @@
 }
 
+// --------------------------------------------------------------------------
+//
+// Take a lot of histograms and plot them together in one plot.
+// Calls PlotSame
+//
 Bool_t MJSpectrum::DisplaySize(MParList &plist, Double_t scale) const
 {
@@ -698 +781 @@
     gPad->SetGridy();
 
-    excess->SetTitle("Number of excess events vs Size (data, mc/blue)");
+    excess->SetTitle("Excess events vs Size (data/black, mc/blue)");
     excess = excess->DrawCopy("E2");
     // Don't do this on the original object!
@@ -727 +810 @@
             histsel->KolmogorovTest(excess, "DX");
             fLog->Separator("Chi^2 Test");
-            histsel->Chi2Test(excess, "P");
+            const Double_t p = histsel->Chi2Test(excess, "P");
+
+            TLatex tex;
+            tex.SetBit(TLatex::kTextNDC);
+            tex.DrawLatex(0.7, 0.93, Form("P(\\chi^{2})=%.0f", p*100));
         }
     }
@@ -814 +901 @@
     if (!GetThetaDistribution(temp1, temp2))
         return kFALSE;
+
+    if (!fNoThetaWeights)
+        weight.SetZdWeights(&temp1);
 
     TH1D excess;
@@ -930 +1020 @@
 
     tlist2.AddToList(&read);
-    if (!fRawMc)
+    if (!fRawMc && fNoThetaWeights)
         tlist2.AddToList(&contsel);
     tlist2.AddToList(&calc);