Index: trunk/MagicSoft/Mars/mjobs/MJSpectrum.cc
===================================================================
--- trunk/MagicSoft/Mars/mjobs/MJSpectrum.cc	(revision 7091)
+++ trunk/MagicSoft/Mars/mjobs/MJSpectrum.cc	(revision 7094)
@@ -52,4 +52,5 @@
 #include "MBinning.h"
 #include "MDataSet.h"
+#include "MMcCorsikaRunHeader.h"
 
 // Spectrum
@@ -58,4 +59,5 @@
 #include "../mhflux/MHCollectionArea.h"
 #include "../mhflux/MHEnergyEst.h"
+#include "../mhflux/MMcSpectrumWeight.h"
 
 // Eventloop
@@ -155,4 +157,48 @@
 }
 
+// --------------------------------------------------------------------------
+//
+// Read the first MMcCorsikaRunHeader from the RunHeaders tree in
+// the dataset.
+// The simulated energy range and spectral slope is initialized from
+// there.
+// In the following eventloops the forced check in MMcSpectrumWeight
+// ensures, that the spectral slope and energy range doesn't change.
+//
+Bool_t  MJSpectrum::InitWeighting(const MDataSet &set, MMcSpectrumWeight &w) const
+{
+    fLog->Separator("Initialize energy weighting");
+
+    MParList l;
+    l.AddToList(&w);
+    if (!CheckEnv(l))
+    {
+        *fLog << err << "ERROR - Reading resources for MMcSpectrumWeight failed." << endl;
+        return kFALSE;
+    }
+
+    TChain chain("RunHeaders");
+    set.AddFilesOn(chain);
+
+    MMcCorsikaRunHeader *h=0;
+    chain.SetBranchAddress("MMcCorsikaRunHeader.", &h);
+    chain.GetEntry(1);
+
+    if (!h)
+    {
+        *fLog << err << "ERROR - Couldn't read MMcCorsikaRunHeader from DataSet." << endl;
+        return kFALSE;
+    }
+
+    if (!w.Set(*h))
+    {
+        *fLog << err << "ERROR - Initializing MMcSpectrumWeight failed." << endl;
+        return kFALSE;
+    }
+
+    w.Print();
+    return kTRUE;
+}
+
 Float_t MJSpectrum::ReadInput(MParList &plist, TH1D &h1, TH1D &h2)
 {
@@ -211,10 +257,15 @@
 }
 
-Bool_t MJSpectrum::ReadOrigMCDistribution(const MDataSet &set, TH1 &h) const
+Bool_t MJSpectrum::ReadOrigMCDistribution(const MDataSet &set, TH1 &h, MMcSpectrumWeight &weight) const
 {
     // Some debug output
     fLog->Separator("Compiling original MC distribution");
 
-    *fLog << inf << "Please stand by, this may take a while..." << flush;
+    weight.SetNameMcEvt("MMcEvtBasic");
+    const TString w(weight.GetFormulaWeights());
+    weight.SetNameMcEvt();
+
+    *fLog << inf << "Using weights: " << w << endl;
+    *fLog << "Please stand by, this may take a while..." << flush;
 
     if (fDisplay)
@@ -236,5 +287,5 @@
         h.SetYTitle("E [GeV]");
         h.SetZTitle("Counts");
-        chain.Draw("MMcEvtBasic.fEnergy:MMcEvtBasic.fTelescopeTheta*TMath::RadToDeg()>>ThetaEMC", "", "goff");
+        chain.Draw("MMcEvtBasic.fEnergy:MMcEvtBasic.fTelescopeTheta*TMath::RadToDeg()>>ThetaEMC", w, "goff");
     }
     else
@@ -243,5 +294,5 @@
         h.SetXTitle("\\Theta [\\circ]");
         h.SetYTitle("Counts");
-        chain.Draw("MMcEvtBasic.fTelescopeTheta*TMath::RadToDeg()>>ThetaMC", "", "goff");
+        chain.Draw("MMcEvtBasic.fTelescopeTheta*TMath::RadToDeg()>>ThetaMC", w, "goff");
     }
     h.SetDirectory(0);
@@ -413,5 +464,5 @@
 }
 
-Bool_t MJSpectrum::IntermediateLoop(MParList &plist, MH3 &mh1, TH1D &temp1, const MDataSet &set) const
+Bool_t MJSpectrum::IntermediateLoop(MParList &plist, MH3 &mh1, TH1D &temp1, const MDataSet &set, MMcSpectrumWeight &weight) const
 {
     MTaskList tlist1;
@@ -439,4 +490,5 @@
     }
     tlist1.AddToList(&readmc);
+    tlist1.AddToList(&weight);
 
     temp1.SetXTitle("MMcEvtBasic.fTelescopeTheta*kRad2Deg");
@@ -477,8 +529,9 @@
         bins3->SetName("BinningTheta");
     }
+
     return kTRUE;
 }
 
-void MJSpectrum::DisplaySpectrum(MHCollectionArea &area, TH1D &excess, MHEnergyEst &hest, Double_t ontime) const
+TArrayD MJSpectrum::DisplaySpectrum(MHCollectionArea &area, TH1D &excess, MHEnergyEst &hest, Double_t ontime) const
 {
     TH1D collarea(area.GetHEnergy());
@@ -546,5 +599,5 @@
     f.SetParameter(1, 1.9e-6);
     f.SetLineColor(kGreen);
-    spectrum.Fit(&f, "NI", "", 55, 2e4);
+    spectrum.Fit(&f, "NIM", "", 55, 2e4);
     f.DrawCopy("same");
 
@@ -558,9 +611,13 @@
      tex.DrawLatex(2e2, 7e-5, str);
      */
-}
-
-Bool_t MJSpectrum::PlotSame(MStatusArray &arr, MParList &plist, const char *name, const char *tab, const char *plot) const
-{
-    cout << name << endl;
+
+    TArrayD res(2);
+    res[0] = f.GetParameter(0);
+    res[1] = f.GetParameter(1);
+    return res;
+}
+
+Bool_t MJSpectrum::PlotSame(MStatusArray &arr, MParList &plist, const char *name, const char *tab, const char *plot, Double_t scale) const
+{
     TString same(name);
     same += "Same";
@@ -587,13 +644,13 @@
 
     const MAlphaFitter *fit = (MAlphaFitter*)plist.FindObject("MAlphaFitter");
-    const Double_t scale = fit ? fit->GetScaleFactor() : 1;
+    const Double_t ascale = fit ? fit->GetScaleFactor() : 1;
 
     gPad->SetBorderMode(0);
     h2->SetLineColor(kBlack);
     h3->SetLineColor(kBlue);
-    h2->Add(h1, -scale);
-
-    h2->Scale(1./h2->Integral());
-    h3->Scale(1./h3->Integral());
+    h2->Add(h1, -ascale);
+
+    //h2->Scale(1./ontime);   //h2->Integral());
+    h3->Scale(scale);    //h3->Integral());
 
     h2->SetMaximum(1.05*TMath::Max(h2->GetMaximum(), h3->GetMaximum()));
@@ -609,5 +666,5 @@
 }
 
-Bool_t MJSpectrum::DisplaySize(MParList &plist) const
+Bool_t MJSpectrum::DisplaySize(MParList &plist, Double_t scale) const
 {
     *fLog << inf << "Reading from file: " << fPathIn << endl;
@@ -644,21 +701,33 @@
 
     excess->SetTitle("Number of excess events vs Size (data, mc/blue)");
-    excess->Scale(1./excess->Integral());
-    excess = excess->DrawCopy();
+    excess = excess->DrawCopy("E2");
     // Don't do this on the original object!
     excess->SetStats(kFALSE);
+    excess->SetMarkerStyle(kFullDotMedium);
+    excess->SetFillColor(kBlack);
+    excess->SetFillStyle(0);
+    excess->SetName("Excess  ");
+    excess->SetDirectory(0);
 
     TObject *o=0;
-    if ((o=plist.FindObject("ExcessSize")))
+    if ((o=plist.FindObject("ExcessMC")))
     {
         TH1 *histsel = (TH1F*)o->FindObject("");
         if (histsel)
         {
-            histsel->Scale(1./histsel->Integral());
+            if (scale<0)
+                scale = excess->Integral()/histsel->Integral();
+
+            histsel->Scale(scale);
             histsel->SetLineColor(kBlue);
             histsel->SetBit(kCanDelete);
-            histsel = histsel->DrawCopy("same");
+            histsel = histsel->DrawCopy("E1 same");
             // Don't do this on the original object!
             histsel->SetStats(kFALSE);
+
+            fLog->Separator("Kolmogorov Test");
+            histsel->KolmogorovTest(excess, "DX");
+            fLog->Separator("Chi^2 Test");
+            histsel->Chi2Test(excess, "P");
         }
     }
@@ -666,17 +735,17 @@
     // -------------- Comparison of Image Parameters --------------
     c.cd(2);
-    PlotSame(arr, plist, "Dist",   "HilSrc",  "MHHilSrcMCPost");
+    PlotSame(arr, plist, "Dist",   "HilSrc",  "MHHilSrcMCPost", scale);
 
     c.cd(3);
-    PlotSame(arr, plist, "Length", "PostCut", "MHHillasMCPost");
+    PlotSame(arr, plist, "Length", "PostCut", "MHHillasMCPost", scale);
 
     c.cd(4);
-    PlotSame(arr, plist, "M3l",    "HilExt",  "MHHilExtMCPost");
+    PlotSame(arr, plist, "M3l",    "HilExt",  "MHHilExtMCPost", scale);
 
     c.cd(5);
-    PlotSame(arr, plist, "Conc1",  "NewPar",  "MHNewParMCPost");
+    PlotSame(arr, plist, "Conc1",  "NewPar",  "MHNewParMCPost", scale);
 
     c.cd(6);
-    PlotSame(arr, plist, "Width",  "PostCut", "MHHillasMCPost");
+    PlotSame(arr, plist, "Width",  "PostCut", "MHHillasMCPost", scale);
 
     return kTRUE;
@@ -734,9 +803,13 @@
     }
 
+    MMcSpectrumWeight weight;
+    if (!InitWeighting(set, weight))
+        return kFALSE;
+
     PrintSetup(fit);
     bins3.SetEdges(temp1, 'x');
 
     TH1D temp2(temp1);
-    if (!ReadOrigMCDistribution(set, temp2))
+    if (!ReadOrigMCDistribution(set, temp2, weight))
         return kFALSE;
 
@@ -754,5 +827,5 @@
         hist.UseCurrentStyle();
         MH::SetBinning(&hist, &bins3/*temp1.GetXaxis()*/, &bins2/*excess.GetXaxis()*/);
-        if (!ReadOrigMCDistribution(set, hist))
+        if (!ReadOrigMCDistribution(set, hist, weight))
             return kFALSE;
 
@@ -762,9 +835,14 @@
                 for (int x=0; x<hist.GetNbinsX(); x++)
                     hist.SetBinContent(x, y, hist.GetBinContent(x, y)*temp1.GetBinContent(x));
+            //hist.SetEntries(hist.Integral());
         }
     }
     else
-        if (!IntermediateLoop(plist, mh1, temp1, set))
+    {
+        weight.SetNameMcEvt("MMcEvtBasic");
+        if (!IntermediateLoop(plist, mh1, temp1, set, weight))
             return kFALSE;
+        weight.SetNameMcEvt();
+    }
 
     DisplayResult(fSimpleMode ? hist : (TH2D&)mh1.GetHist());
@@ -813,12 +891,13 @@
     MFillH fill3(&area, "", "FillCollectionArea");
     MFillH fill4(&hest, "", "FillEnergyEst");
+    fill3.SetWeight();
+    fill4.SetWeight();
 
     MH3 hsize("MHillas.fSize");
-    //MH3 henergy("MEnergyEst.fVal");
-    hsize.SetName("ExcessSize");
-    //henergy.SetName("EnergyEst");
-    MBinning bins(size, "BinningExcessSize");
+    hsize.SetName("ExcessMC");
+    hsize.Sumw2();
+
+    MBinning bins(size, "BinningExcessMC");
     plist.AddToList(&hsize);
-    //plist.AddToList(&henergy);
     plist.AddToList(&bins);
 
@@ -830,6 +909,5 @@
     MFillH fill6a("MHImgParMCPost [MHImagePar]",    "MImagePar",    "FillImgParPost");
     MFillH fill7a("MHNewParMCPost [MHNewImagePar]", "MNewImagePar", "FillNewParPost");
-    MFillH fill8a("ExcessSize     [MH3]",           "",             "FillExcessSize");
-    //MFillH fill9a("EnergyEst      [MH3]",           "",             "FillExcessEEst");
+    MFillH fill8a("ExcessMC       [MH3]",           "",             "FillExcessMC");
     fill1a.SetNameTab("PreCut");
     fill2a.SetNameTab("PostCut");
@@ -840,5 +918,12 @@
     fill7a.SetNameTab("NewPar");
     fill8a.SetBit(MFillH::kDoNotDisplay);
-    //fill9a.SetBit(MFillH::kDoNotDisplay);
+    fill1a.SetWeight();
+    fill2a.SetWeight();
+    fill3a.SetWeight();
+    fill4a.SetWeight();
+    fill5a.SetWeight();
+    fill6a.SetWeight();
+    fill7a.SetWeight();
+    fill8a.SetWeight();
 
     MEnergyEstimate est;
@@ -852,4 +937,5 @@
     tlist2.AddToList(&hcalc1);
     tlist2.AddToList(&hcalc2);
+    tlist2.AddToList(&weight);
     tlist2.AddToList(&fill1a);
     tlist2.AddToList(fCut0);
@@ -891,7 +977,29 @@
     // -------------------------- Spectrum ----------------------------
 
-    DisplaySpectrum(area, excess, hest, ontime);
-    DisplaySize(plist);
-
+    // Calculate and display spectrum (N/TeVsm^2 at 1TeV)
+    TArrayD res(DisplaySpectrum(area, excess, hest, ontime));
+
+    // Spectrum fitted (convert res[1] from TeV to GeV)
+    TF1 flx("flx", Form("%e*pow(x/1000, %f)", res[1]/1000, res[0]));
+
+    // Number of events this spectrum would produce per s and m^2
+    Double_t n = flx.Integral(weight.GetEnergyMin(), weight.GetEnergyMax());
+
+    // scale with effective collection area to get the event rate (N/s)
+    // scale with the effective observation time to absolute observed events
+    n *= area.GetCollectionAreaAbs()*ontime; // N
+
+    // Now calculate the scale factor from the number of events
+    // produced and the number of events which should have been
+    // observed with our telescope in the time ontime
+    const Double_t scale = n/area.GetEntries();
+
+    // Print normalization constant
+    cout << "MC normalization factor:  " << scale << endl;
+
+    // Overlay normalized plots
+    DisplaySize(plist, scale);
+
+    // check if output should be written
     if (!fPathOut.IsNull())
         fDisplay->SaveAsRoot(fPathOut);
Index: trunk/MagicSoft/Mars/mjobs/MJSpectrum.h
===================================================================
--- trunk/MagicSoft/Mars/mjobs/MJSpectrum.h	(revision 7091)
+++ trunk/MagicSoft/Mars/mjobs/MJSpectrum.h	(revision 7094)
@@ -18,4 +18,5 @@
 class MStatusArray;
 class MHCollectionArea;
+class MMcSpectrumWeight;
 
 class MJSpectrum : public MJob
@@ -35,15 +36,16 @@
     Bool_t  ReadTask(MTask* &task, const char *name) const;
     Float_t ReadInput(MParList &plist, TH1D &h1, TH1D &size);
-    Bool_t  ReadOrigMCDistribution(const MDataSet &set, TH1 &h) const;
+    Bool_t  ReadOrigMCDistribution(const MDataSet &set, TH1 &h, MMcSpectrumWeight &w) const;
     Bool_t  GetThetaDistribution(TH1D &temp1, TH1D &temp2) const;
     Bool_t  Refill(MParList &plist, TH1D &h) const;
+    Bool_t  InitWeighting(const MDataSet &set, MMcSpectrumWeight &w) const;
 
     // Display Output
     void    PrintSetup(const MAlphaFitter &fit) const;
     void    DisplayResult(const TH2D &mh1) const;
-    Bool_t  IntermediateLoop(MParList &plist, MH3 &h1, TH1D &temp1, const MDataSet &set) const;
-    void    DisplaySpectrum(MHCollectionArea &area, TH1D &excess, MHEnergyEst &hest, Double_t ontime) const;
-    Bool_t  DisplaySize(MParList &plist) const;
-    Bool_t  PlotSame(MStatusArray &arr, MParList &plist, const char *name, const char *tab, const char *plot) const;
+    Bool_t  IntermediateLoop(MParList &plist, MH3 &h1, TH1D &temp1, const MDataSet &set, MMcSpectrumWeight &w) const;
+    TArrayD DisplaySpectrum(MHCollectionArea &area, TH1D &excess, MHEnergyEst &hest, Double_t ontime) const;
+    Bool_t  DisplaySize(MParList &plist, Double_t scale) const;
+    Bool_t  PlotSame(MStatusArray &arr, MParList &plist, const char *name, const char *tab, const char *plot, Double_t scale) const;
 
 public: