Index: trunk/MagicSoft/Mars/mjobs/MJSpectrum.cc
===================================================================
--- trunk/MagicSoft/Mars/mjobs/MJSpectrum.cc	(revision 7390)
+++ trunk/MagicSoft/Mars/mjobs/MJSpectrum.cc	(revision 7403)
@@ -76,5 +76,5 @@
 #include "MTaskEnv.h"
 #include "MFillH.h"
-#include "MHillasCalc.h"
+//#include "MHillasCalc.h"
 //#include "MSrcPosCalc.h"
 #include "MContinue.h"
@@ -87,5 +87,5 @@
     : fCut0(0),fCut1(0), fCut2(0), fCut3(0), fEstimateEnergy(0),
     fRefill(kFALSE), fSimpleMode(kTRUE), fRawMc(kFALSE),
-    fNoThetaWeights(kFALSE)
+    fNoThetaWeights(kFALSE), fWeighting(kTRUE)
 {
     fName  = name  ? name  : "MJSpectrum";
@@ -551,4 +551,88 @@
 }
 
+TString MJSpectrum::FormString(const TF1 &f, Byte_t type)
+{
+    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);
+
+    const Float_t fe0 = TMath::Log10(e0);
+    const Float_t fe1 = TMath::Log10(e1);
+
+    // calc number of (gueltige ziffern) digits to be displayed
+    const char *f0  = fe0-TMath::Floor(fe0)>0.3 ? "3.1" : "1.0";
+    const char *f1  = fe1-TMath::Floor(fe1)>0.3 ? "3.1" : "1.0";
+
+    TString str;
+    switch (type)
+    {
+    case 0:
+        {
+            const TString fmt0 = Form("(%%%sf#pm%%%sf)\\bullet10^{%%d}",   f1, f1);
+            const TString fmt1 = Form("(\\frac{E}{TeV})^{%%%sf#pm%%%sf}", f0, f0);
+
+            str  = Form(fmt0.Data(), p1/exp, e1/exp, np);
+            str += Form(fmt1.Data(), p0, e0);
+            str += "\\frac{ph}{TeVm^{2}s}";
+        }
+        break;
+    case 1:
+        str = Form("\\chi^{2}/NDF=%.2f/%d", f.GetChisquare(),f.GetNDF());
+        break;
+    case 2:
+        str = Form("P=%.0f%%", 100*TMath::Prob(f.GetChisquare(),  f.GetNDF()));
+        break;
+    }
+    return str;
+}
+
+TArrayD MJSpectrum::FitSpectrum(TH1D &spectrum) const
+{
+    TF1 f("f", "[1]*(x/1e3)^[0]", 10, 3e4);
+    f.SetParameter(0, -2.5);
+    f.SetParameter(1, spectrum.GetBinContent(spectrum.GetXaxis()->FindFixBin(1000)));
+    f.SetLineColor(kBlue);
+    spectrum.Fit(&f, "NI", "", 80, 1150); // M skipped
+    f.DrawCopy("same");
+
+    /*
+    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 = FormString(f);
+    /*
+    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.SetTextAlign(31);
+    tex.DrawLatex(0.89, 0.935, str);
+
+    str = FormString(f, 1);
+    tex.DrawLatex(0.89, 0.83, str);
+
+    str = FormString(f, 2);
+    tex.DrawLatex(0.89, 0.735, str);
+
+    TArrayD res(2);
+    res[0] = f.GetParameter(0);
+    res[1] = f.GetParameter(1);
+    return res;
+}
+
 // --------------------------------------------------------------------------
 //
@@ -568,4 +652,7 @@
 
     cout << "Effective On time: " << ontime << "s" << endl;
+
+//    for (int i=0; i<collarea.GetNbinsX(); i++)
+//        collarea.SetBinError(i+1, 0);//collarea.GetBinError(i+1)/2);
 
     spectrum.SetDirectory(NULL);
@@ -602,5 +689,5 @@
     weights.DrawCopy();
 
-    //spectrum.Multiply(&weights);
+    //spectrum.Multiply(&weights); // Currentyl we skip the error distribution!
     spectrum.SetNameTitle("Flux", "Spectrum (after unfolding)");
     spectrum.SetBit(TH1::kNoStats);
@@ -608,9 +695,16 @@
     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));
+        if (i+1>2)
+        {
+            //Float_t E = spectrum.GetXaxis()->GetBinWidth(i+1);
+            //Float_t c = log10(E)*3.012e-1+1.905e-1;
+            //spectrum.SetBinContent(i+1, spectrum.GetBinContent(i+1)*c);
+            //spectrum.SetBinContent(i+1, spectrum.GetBinContent(i+1)*c);
+            spectrum.SetBinError(i+1,   spectrum.GetBinError(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);
+        spectrum.SetBinError(i+1,   spectrum.GetBinError(i+1)  /spectrum.GetBinWidth(i+1)*1000);
     }
 
@@ -623,11 +717,33 @@
     spectrum.SetMinimum(1e-12);
     spectrum.SetXTitle("E [GeV]");
-    spectrum.SetYTitle("N/TeVsm^{2}");
+    spectrum.SetYTitle("dN/dE [N/TeVsm^{2}]");
+    TH1D *spec = (TH1D*)spectrum.DrawCopy();
+
+    // Calculate Spectrum * E^2
+    for (int i=0; i<spectrum.GetNbinsX(); i++)
+    {
+        spectrum.SetBinContent(i+1, spectrum.GetBinContent(i+1)*spectrum.GetBinWidth(i+1)*spectrum.GetBinWidth(i+1));
+        spectrum.SetBinError(  i+1, spectrum.GetBinError(i+1)  *spectrum.GetBinWidth(i+1)*spectrum.GetBinWidth(i+1));
+    }
+    fDisplay->AddTab("Check");
+    gPad->SetLogx();
+    gPad->SetLogy();
+    gPad->SetGrid();
+
+    spectrum.SetName("FluxStd");
+    spectrum.SetMarkerStyle(kFullDotMedium);
+    spectrum.SetTitle("Differential flux times E^{2}");
+    spectrum.SetYTitle("E^{2}dN/dE [N TeV/sm^{2}]");
+    spectrum.SetDirectory(0);
     spectrum.DrawCopy();
 
+    // Fit Spectrum
+    c1.cd(4);
+    return FitSpectrum(*spec);
+    /*
     TF1 f("f", "[1]*(x/1e3)^[0]", 10, 3e4);
     f.SetParameter(0, -2.87);
     f.SetParameter(1, 1.9e-6);
-    f.SetLineColor(kGreen);
+    f.SetLineColor(kBlue);
     spectrum.Fit(&f, "NIM", "", 100, 5000);
     f.DrawCopy("same");
@@ -657,5 +773,5 @@
     TArrayD res(2);
     res[0] = f.GetParameter(0);
-    res[1] = f.GetParameter(1);
+    res[1] = f.GetParameter(1);*/
 
 /*
@@ -691,6 +807,6 @@
     f.SetLineStyle(kDashed);
     f.DrawCopy("same");
+    return res;
  */
-    return res;
 }
 
@@ -754,5 +870,6 @@
 // Calls PlotSame
 //
-Bool_t MJSpectrum::DisplaySize(MParList &plist, Double_t scale) const
+//#include <TSpline.h>
+Bool_t MJSpectrum::DisplaySize(MParList &plist, Double_t scale, TH1D &result) const
 {
     *fLog << inf << "Reading from file: " << fPathIn << endl;
@@ -773,5 +890,5 @@
     }
 
-    TH1 *excess = (TH1D*)arr.FindObjectInCanvas("Excess", "TH1D", "Hist");
+    TH1D *excess = (TH1D*)arr.FindObjectInCanvas("Excess", "TH1D", "Hist");
     if (!excess)
         return kFALSE;
@@ -789,5 +906,5 @@
 
     excess->SetTitle("Excess events vs Size (data/black, mc/blue)");
-    excess = excess->DrawCopy("E2");
+    excess = (TH1D*)excess->DrawCopy("E2");
     // Don't do this on the original object!
     excess->SetStats(kFALSE);
@@ -801,7 +918,72 @@
     if ((o=plist.FindObject("ExcessMC")))
     {
-        TH1 *histsel = (TH1F*)o->FindObject("");
+        TH1 *histsel = (TH1D*)o->FindObject("");
         if (histsel)
         {
+            // Calculate the weights such that the number of total
+            // MC events after cuts will stay constant
+
+            /*
+            Int_t first = -1;
+            for (int i=0; i<excess->GetNbinsX(); i++)
+                if (excess->GetBinContent(i+1)>0)
+                {
+                    first = i+1;
+                    break;
+                }
+
+            Int_t last = -1;
+            for (int i=excess->GetNbinsX(); i>0; i--)
+                if (excess->GetBinContent(i)>0)
+                {
+                    last = i;
+                    break;
+                    }*/
+            /*
+            Int_t first = 1;
+            Int_t last  = excess->GetNbinsX();
+
+            const Double_t xmin = TMath::Log10(excess->GetBinLowEdge(first));
+            const Double_t xmax = TMath::Log10(excess->GetBinLowEdge(last+1));
+            const Int_t    xnum = last-first+1;
+
+            //            Double_t xmin = TMath::Log10(excess->GetXaxis()->GetXmin());
+            //            Double_t xmax = TMath::Log10(excess->GetXaxis()->GetXmax());
+            //            Int_t    xnum = excess->GetNbinsX();
+
+            TArrayD arr1(xnum);
+            TArrayD arr2(xnum);
+            for (int i=0; i<xnum; i++)
+            {
+                Double_t v1 = excess->GetBinContent(i+first);
+                Double_t v2 = histsel->GetBinContent(i+first);
+                arr1[i] = v1<1 ? -1 : TMath::Log10(v1);
+                arr2[i] = v2<1 ? -1 : TMath::Log10(v2);
+            }
+            cout << endl << endl;
+
+            TSpline5 sp1("Excess", xmin, xmax, arr1.GetArray(), xnum+1);
+            TSpline5 sp2("MC",     xmin, xmax, arr2.GetArray(), xnum+1);
+
+            MBinning bins(xnum*5, TMath::Power(10, xmin), TMath::Power(10, xmax), "", "log");
+            bins.Apply(result);
+
+            for (int i=0; i<result.GetNbinsX(); i++)
+            {
+                Double_t lo = result.GetXaxis()->GetBinLowEdge(i+1);
+                Double_t up = result.GetXaxis()->GetBinUpEdge(i+1);
+
+                Double_t x  = TMath::Log10(lo*up)/2;
+                Double_t v1 = sp1.Eval(x);
+                Double_t v2 = sp2.Eval(x);
+                Double_t rc = v1<-0.9 || v2<-0.9 ? 0 : TMath::Power(10, v1-v2);
+                result.SetBinContent(i+1, rc);
+            }
+            */
+            excess->Copy(result);
+            result.Divide(histsel);
+            result.SetDirectory(0);
+
+            // Apply additional absolute scale factor
             if (scale<0)
                 scale = excess->Integral()/histsel->Integral();
@@ -814,4 +996,5 @@
             histsel->SetStats(kFALSE);
 
+            // Do some Chi^2 and Kolmogorov tests
             fLog->Separator("Kolmogorov Test");
             histsel->KolmogorovTest(excess, "DX");
@@ -819,4 +1002,5 @@
             const Double_t p = histsel->Chi2Test(excess, "P");
 
+            // Plot result
             TLatex tex;
             tex.SetBit(TLatex::kTextNDC);
@@ -840,4 +1024,179 @@
     c.cd(6);
     PlotSame(arr, plist, "Width",  "PostCut", "MHHillasMCPost", scale);
+
+    return kTRUE;
+}
+
+Bool_t MJSpectrum::DisplaySpectrumW(TH1D &energy0, TH1 &energy, TH1 &hsize, TH1D &weights, TH2D &hall, Double_t scale, TArrayD &res) const
+{
+    *fLog << inf << "Reading from file: " << fPathIn << endl;
+
+    TFile file(fPathIn, "READ");
+    if (!file.IsOpen())
+    {
+        *fLog << err << dbginf << "ERROR - Could not open file " << fPathIn << endl;
+        return kFALSE;
+    }
+
+    file.cd();
+    MStatusArray arr;
+    if (arr.Read()<=0)
+    {
+        *fLog << "MStatusDisplay not found in file... abort." << endl;
+        return kFALSE;
+    }
+
+    TH1D *exc = (TH1D*)arr.FindObjectInCanvas("Excess", "TH1D", "Hist");
+    if (!exc)
+        return kFALSE;
+
+    // energy: energy distribution after weighting
+    // energy0: energy distribution before weighting
+    // weights: weights to weight events depending on their size
+
+    TCanvas &c = fDisplay->AddTab("SpectrumW");
+    c.Divide(3,2);
+    c.cd(1);
+    gPad->SetBorderMode(0);
+    gPad->SetLogx();
+    gPad->SetLogy();
+    gPad->SetGrid();
+
+    TH1D *e = (TH1D*)exc->DrawCopy("E2");
+    // Don't do this on the original object!
+    e->SetName("ExcessData");
+    e->SetTitle("Excess events vs Size (data/black, mc/blue)");
+    e->SetDirectory(0);
+    e->SetStats(kFALSE);
+    e->SetBit(TH1::kNoStats);
+    e->ResetBit(TH1::kNoTitle);
+    e->SetMarkerStyle(kFullDotMedium);
+    e->SetFillColor(kBlack);
+    e->SetFillStyle(0);
+    e->SetName("Excess  ");
+    e->SetDirectory(0);
+
+    e = (TH1D*)hsize.DrawCopy("E1 same");
+    // Don't do this on the original object!
+    e->SetName("ExcessMC");
+    e->SetDirectory(0);
+    e->SetBit(kCanDelete);
+    e->SetBit(TH1::kNoStats);
+    e->ResetBit(TH1::kNoTitle);
+    e->SetLineColor(kBlue);
+    e->SetStats(kFALSE);
+
+    c.cd(2);
+    gPad->SetBorderMode(0);
+    gPad->SetLogx();
+    gPad->SetGrid();
+
+    weights.SetName("WeightsSize");
+    weights.SetDirectory(0);
+    weights.SetBit(TH1::kNoStats);
+    weights.ResetBit(TH1::kNoTitle);
+    weights.SetTitle("Size weights: size(data)/size(mc)");
+    weights.DrawCopy();
+
+    c.cd(4);
+    gPad->SetBorderMode(0);
+    gPad->SetLogx();
+    gPad->SetLogy();
+    gPad->SetGrid();
+    energy.SetBit(TH1::kNoStats);
+    energy.ResetBit(TH1::kNoTitle);
+    energy.SetTitle("Normalized energy distribution before (black) and after (blue) weighting");
+    energy.SetLineColor(kBlue);
+    TH1 *he = energy.DrawCopy();
+    TH1 *hf = energy0.DrawCopy("same");
+    he->Scale(1./he->Integral());
+    hf->Scale(1./hf->Integral());
+    he->SetName("EnergyOrig");
+    hf->SetName("EnergyWeighted");
+    he->SetDirectory(0);
+    hf->SetDirectory(0);
+
+    c.cd(5);
+    gPad->SetBorderMode(0);
+    gPad->SetLogx();
+    gPad->SetGrid();
+    energy0.SetBit(TH1::kNoStats);
+    energy0.ResetBit(TH1::kNoTitle);
+    energy0.SetName("WeightsEnergy");
+    energy0.SetTitle("Energy weights: E(before weighting)/E(after weighting)");
+    energy0.SetDirectory(0);
+    energy0.Divide(&energy);
+    energy0.DrawCopy();
+
+    c.cd(3);
+    gPad->SetBorderMode(0);
+    gPad->SetLogx();
+    gPad->SetLogy();
+    gPad->SetGrid();
+
+    TH1D *h1 = (TH1D*)hall.ProjectionY("ProjAllY", -1, 9999, "E");
+    h1->SetTitle("Normalized spectrum before and after Energy weights applied");
+    h1->SetDirectory(NULL);
+    h1->SetXTitle("E [GeV]");
+    h1->SetBit(kCanDelete);
+    h1->SetBit(TH1::kNoStats);
+    he = h1->DrawCopy();
+    h1->SetLineColor(kBlue);
+    // Convert the mc spectrum before weights to the mc spectrum after weights
+    // Therefor weights after/before are applied
+    h1->Divide(&energy0);
+    hf = h1->DrawCopy("same");
+    he->Scale(1./he->GetBinContent(he->GetNbinsX()/2));
+    hf->Scale(1./hf->GetBinContent(hf->GetNbinsX()/2));
+    he->SetName("SpectrumOrig");
+    hf->SetName("SpectrumWeighted");
+    he->SetDirectory(0);
+    hf->SetDirectory(0);
+
+    c.cd(6);
+    gPad->SetBorderMode(0);
+    gPad->SetLogx();
+    gPad->SetLogy();
+    gPad->SetGrid();
+
+    h1->SetName("Spectrum");
+    h1->SetTitle("Spectrum");
+    h1->SetXTitle("E [GeV]");
+    h1->SetYTitle("dN/dE [N/TeVsm^{2}]");
+    h1->SetDirectory(0);
+    h1->Scale(1./scale);
+
+    for (int i=0; i<h1->GetNbinsX(); i++)
+    {
+        h1->SetBinContent(i+1, h1->GetBinContent(i+1)/h1->GetBinWidth(i+1)*1000);
+        h1->SetBinError(  i+1, h1->GetBinError(i+1)  /h1->GetBinWidth(i+1)*1000);
+    }
+
+    TH1D *h2 = (TH1D*)h1->DrawCopy();
+    res = FitSpectrum(*h2);
+
+    // Calculate and Display Spectrum * E^2
+    TCanvas *c1 = fDisplay->GetCanvas("Check");
+    if (c1)
+    {
+        c1->cd();
+        gPad->SetLogx();
+        gPad->SetLogy();
+        gPad->SetGrid();
+
+        for (int i=0; i<h1->GetNbinsX(); i++)
+        {
+            h1->SetBinContent(i+1, h1->GetBinContent(i+1)*h1->GetBinWidth(i+1)*h1->GetBinWidth(i+1));
+            h1->SetBinError(  i+1, h1->GetBinError(i+1)  *h1->GetBinWidth(i+1)*h1->GetBinWidth(i+1));
+        }
+        h1->SetName("FluxW");
+        h1->SetTitle("Differential flux times E^{2}");
+        h1->SetDirectory(0);
+        h1->SetYTitle("E^{2}dN/dE [N TeV/sm^{2}]");
+        h1->SetLineColor(kBlue);
+        h1->SetMarkerColor(kBlue);
+        h1->SetMarkerStyle(kFullDotMedium);
+        h1->Draw("same");
+    }
 
     return kTRUE;
@@ -946,168 +1305,314 @@
     // ------------------------- Final loop --------------------------
 
-    *fLog << endl;
-    fLog->Separator("Calculate efficiencies");
-    *fLog << endl;
-
-    MTaskList tlist2;
-    plist.Replace(&tlist2);
-
-    MReadMarsFile read("Events");
-    read.DisableAutoScheme();
-    set.AddFilesOn(read);
-
-    // Selector to get correct (final) theta-distribution
-    temp1.SetXTitle("MPointingPos.fZd");
-
-    MH3 mh3(temp1);
-
-    MFEventSelector2 sel2(mh3);
-    sel2.SetHistIsProbability();
-
-    MContinue contsel(&sel2);
-    contsel.SetInverted();
-
-    // Get correct source position
-    //MSrcPosCalc calc;
-
-    // Calculate corresponding Hillas parameters
-    /*
-     MHillasCalc hcalc1;
-     MHillasCalc hcalc2("MHillasCalcAnti");
-     hcalc1.SetFlags(MHillasCalc::kCalcHillasSrc);
-     hcalc2.SetFlags(MHillasCalc::kCalcHillasSrc);
-     hcalc2.SetNameHillasSrc("MHillasSrcAnti");
-     hcalc2.SetNameSrcPosCam("MSrcPosAnti");
-     */
-
-    // Fill collection area and energy estimator (unfolding)
-    // Make sure to use the same binning for MHCollectionArea and MHEnergyEst
-    MHCollectionArea area;
-    area.SetHistAll(fSimpleMode ? hist : (TH2D&)mh1.GetHist());
-    MHEnergyEst      hest;
-
-    MFillH fill3(&area, "", "FillCollectionArea");
-    MFillH fill4(&hest, "", "FillEnergyEst");
-    MFillH fill5("MHThreshold", "", "FillThreshold");
-    fill3.SetWeight();
-    fill4.SetWeight();
-    fill5.SetWeight();
-    fill3.SetNameTab("ColArea");
-    fill4.SetNameTab("E-Est");
-    fill5.SetNameTab("Threshold");
-
-    MH3 hsize("MHillas.fSize");
-    hsize.SetName("ExcessMC");
-    hsize.Sumw2();
-
-    MBinning bins(size, "BinningExcessMC");
-    plist.AddToList(&hsize);
-    plist.AddToList(&bins);
-
-    MFillH fill1a("MHHillasMCPre  [MHHillas]",      "MHillas",      "FillHillasPre");
-    MFillH fill2a("MHHillasMCPost [MHHillas]",      "MHillas",      "FillHillasPost");
-    MFillH fill3a("MHVsSizeMCPost [MHVsSize]",      "MHillasSrc",   "FillVsSizePost");
-    MFillH fill4a("MHHilExtMCPost [MHHillasExt]",   "MHillasSrc",   "FillHilExtPost");
-    MFillH fill5a("MHHilSrcMCPost [MHHillasSrc]",   "MHillasSrc",   "FillHilSrcPost");
-    MFillH fill6a("MHImgParMCPost [MHImagePar]",    "MImagePar",    "FillImgParPost");
-    MFillH fill7a("MHNewParMCPost [MHNewImagePar]", "MNewImagePar", "FillNewParPost");
-    MFillH fill8a("ExcessMC       [MH3]",           "",             "FillExcessMC");
-    fill1a.SetNameTab("PreCut");
-    fill2a.SetNameTab("PostCut");
-    fill3a.SetNameTab("VsSize");
-    fill4a.SetNameTab("HilExt");
-    fill5a.SetNameTab("HilSrc");
-    fill6a.SetNameTab("ImgPar");
-    fill7a.SetNameTab("NewPar");
-    fill8a.SetBit(MFillH::kDoNotDisplay);
-    fill1a.SetWeight();
-    fill2a.SetWeight();
-    fill3a.SetWeight();
-    fill4a.SetWeight();
-    fill5a.SetWeight();
-    fill6a.SetWeight();
-    fill7a.SetWeight();
-    fill8a.SetWeight();
-
-    MEnergyEstimate est;
-    MTaskEnv taskenv1("EstimateEnergy");
-    taskenv1.SetDefault(fEstimateEnergy ? fEstimateEnergy : &est);
-
-    tlist2.AddToList(&read);
-    if (!fRawMc && fNoThetaWeights)
-        tlist2.AddToList(&contsel);
-    //tlist2.AddToList(&calc);
-    //tlist2.AddToList(&hcalc1);
-    //tlist2.AddToList(&hcalc2);
-    tlist2.AddToList(&weight);
-    tlist2.AddToList(&fill1a);
-    tlist2.AddToList(fCut0);
-    tlist2.AddToList(fCut1);
-    tlist2.AddToList(fCut2);
-    tlist2.AddToList(fCut3);
-    tlist2.AddToList(&taskenv1);
-    tlist2.AddToList(&fill3);
-    tlist2.AddToList(&fill4);
-    tlist2.AddToList(&fill5);
-    tlist2.AddToList(&fill2a);
-    tlist2.AddToList(&fill3a);
-    tlist2.AddToList(&fill4a);
-    tlist2.AddToList(&fill5a);
-    tlist2.AddToList(&fill6a);
-    tlist2.AddToList(&fill7a);
-    tlist2.AddToList(&fill8a);
-    //tlist2.AddToList(&fill9a);
-
-    MEvtLoop loop2("FillMonteCarlo"); // ***** fName *****
-    loop2.SetParList(&plist);
-    loop2.SetDisplay(fDisplay);
-    loop2.SetLogStream(fLog);
-
-    if (!SetupEnv(loop2))
-        return kFALSE;
-
-    if (!loop2.Eventloop(fMaxEvents))
-    {
-        *fLog << err << GetDescriptor() << ": Processing of MC-data failed." << endl;
-        return kFALSE;
-    }
-
-    if (!loop2.GetDisplay())
-    {
-        *fLog << err << GetDescriptor() << ": Execution stopped by user." << endl;
-        return kFALSE;
-    }
-
-    gLog.Separator("Energy Estimator");
-    if (plist.FindObject("EstimateEnergy"))
-        plist.FindObject("EstimateEnergy")->Print();
-
-    gLog.Separator("Spectrum");
-
-    // -------------------------- Spectrum ----------------------------
-
-    // 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);
+    TH1D weights;
+    TH1D energy0;
+    weights.SetDirectory(0);
+
+    const Int_t n = fWeighting ? 2 : 1;
+    for (int i=0; i<n; i++)
+    {
+        *fLog << endl;
+        fLog->Separator("Calculate efficiencies");
+        *fLog << endl;
+
+        if (i==1)
+            weight.SetWeightsSize(&weights);
+
+
+        MTaskList tlist2;
+        plist.Replace(&tlist2);
+
+        MReadMarsFile read("Events");
+        read.DisableAutoScheme();
+        set.AddFilesOn(read);
+
+        // Selector to get correct (final) theta-distribution
+        temp1.SetXTitle("MPointingPos.fZd");
+
+        MH3 mh3(temp1);
+
+        MFEventSelector2 sel2(mh3);
+        sel2.SetHistIsProbability();
+
+        MContinue contsel(&sel2);
+        contsel.SetInverted();
+
+//***        // Get correct source position
+//***        MSrcPosCalc calc;
+//***
+//***        // Calculate corresponding Hillas parameters
+//***        MHillasCalc hcalc1;
+//***        MHillasCalc hcalc2("MHillasCalcAnti");
+//***        hcalc1.SetFlags(MHillasCalc::kCalcHillasSrc);
+//***        hcalc2.SetFlags(MHillasCalc::kCalcHillasSrc);
+//***        hcalc2.SetNameHillasSrc("MHillasSrcAnti");
+//***        hcalc2.SetNameSrcPosCam("MSrcPosAnti");
+
+        // Fill collection area and energy estimator (unfolding)
+        // Make sure to use the same binning for MHCollectionArea and MHEnergyEst
+        MHCollectionArea area;
+        area.SetHistAll(fSimpleMode ? hist : (TH2D&)mh1.GetHist());
+        MHEnergyEst      hest;
+
+        MFillH fill3(&area, "", "FillCollectionArea");
+        MFillH fill4(&hest, "", "FillEnergyEst");
+        MFillH fill5("MHThreshold", "", "FillThreshold");
+        fill3.SetWeight();
+        fill4.SetWeight();
+        fill5.SetWeight();
+        fill3.SetNameTab("ColArea");
+        fill4.SetNameTab("E-Est");
+        fill5.SetNameTab("Threshold");
+
+        MH3 hsize("MHillas.fSize");
+        hsize.SetName("ExcessMC");
+        hsize.Sumw2();
+
+        MH3 energy("MMcEvt.fEnergy");
+        energy.SetName("EnergyEst"); // FIXME: Should be EnergyMC, but binning has name EnergyEst
+        energy.SetLogy();
+        energy.Sumw2();
+
+        MBinning bins(size, "BinningExcessMC");
+        plist.AddToList(&hsize);
+        plist.AddToList(&energy);
+        plist.AddToList(&bins);
+
+        MFillH fill1a("MHHillasMCPre  [MHHillas]",      "MHillas",      "FillHillasPre");
+        MFillH fill2a("MHHillasMCPost [MHHillas]",      "MHillas",      "FillHillasPost");
+        MFillH fill3a("MHVsSizeMCPost [MHVsSize]",      "MHillasSrc",   "FillVsSizePost");
+        MFillH fill4a("MHHilExtMCPost [MHHillasExt]",   "MHillasSrc",   "FillHilExtPost");
+        MFillH fill5a("MHHilSrcMCPost [MHHillasSrc]",   "MHillasSrc",   "FillHilSrcPost");
+        MFillH fill6a("MHImgParMCPost [MHImagePar]",    "MImagePar",    "FillImgParPost");
+        MFillH fill7a("MHNewParMCPost [MHNewImagePar]", "MNewImagePar", "FillNewParPost");
+        MFillH fill8a("ExcessMC       [MH3]",           "",             "FillExcessMC");
+        MFillH fill9a("EnergyEst      [MH3]",           "",             "FillEnergy");
+        fill1a.SetNameTab("PreCut");
+        fill2a.SetNameTab("PostCut");
+        fill3a.SetNameTab("VsSize");
+        fill4a.SetNameTab("HilExt");
+        fill5a.SetNameTab("HilSrc");
+        fill6a.SetNameTab("ImgPar");
+        fill7a.SetNameTab("NewPar");
+        fill9a.SetNameTab("Energy");
+        fill8a.SetBit(MFillH::kDoNotDisplay);
+        fill9a.SetBit(MFillH::kDoNotDisplay);
+        fill1a.SetWeight();
+        fill2a.SetWeight();
+        fill3a.SetWeight();
+        fill4a.SetWeight();
+        fill5a.SetWeight();
+        fill6a.SetWeight();
+        fill7a.SetWeight();
+        fill8a.SetWeight();
+        fill9a.SetWeight();
+
+        MEnergyEstimate est;
+        MTaskEnv taskenv1("EstimateEnergy");
+        taskenv1.SetDefault(fEstimateEnergy ? fEstimateEnergy : &est);
+
+        if (i==1)
+        {
+            est.SetName("SetEnergyMC");
+            est.SetRule("MMcEvt.fEnergy");
+        }
+
+        tlist2.AddToList(&read);
+        if (!fRawMc && fNoThetaWeights)
+            tlist2.AddToList(&contsel);
+//***        tlist2.AddToList(&calc);
+//***        tlist2.AddToList(&hcalc1);
+//***        tlist2.AddToList(&hcalc2);
+        tlist2.AddToList(&weight);
+        if (i==n-1)
+            tlist2.AddToList(&fill1a);
+        tlist2.AddToList(fCut0);
+        tlist2.AddToList(fCut1);
+        tlist2.AddToList(fCut2);
+        tlist2.AddToList(fCut3);
+        tlist2.AddToList(i==1 ? (MTask*)&est : (MTask*)&taskenv1);
+        tlist2.AddToList(&fill3);
+        if (i==0)
+            tlist2.AddToList(&fill4);
+        tlist2.AddToList(&fill5);
+        if (i==n-1)
+        {
+            tlist2.AddToList(&fill2a);
+            tlist2.AddToList(&fill3a);
+            tlist2.AddToList(&fill4a);
+            tlist2.AddToList(&fill5a);
+            tlist2.AddToList(&fill6a);
+            tlist2.AddToList(&fill7a);
+        }
+        tlist2.AddToList(&fill8a);
+        tlist2.AddToList(&fill9a);
+
+        MEvtLoop loop2("FillMonteCarlo"); // ***** fName *****
+        loop2.SetParList(&plist);
+        loop2.SetDisplay(fDisplay);
+        loop2.SetLogStream(fLog);
+
+//***<<<<<<< MJSpectrum.cc
+        if (!SetupEnv(loop2))
+            return kFALSE;
+//***=======
+//***    MFEventSelector2 sel2(mh3);
+//***    sel2.SetHistIsProbability();
+//***
+//***    MContinue contsel(&sel2);
+//***    contsel.SetInverted();
+//***
+//***    // Get correct source position
+//***    //MSrcPosCalc calc;
+//***
+//***    // Calculate corresponding Hillas parameters
+//***    /*
+//***     MHillasCalc hcalc1;
+//***     MHillasCalc hcalc2("MHillasCalcAnti");
+//***     hcalc1.SetFlags(MHillasCalc::kCalcHillasSrc);
+//***     hcalc2.SetFlags(MHillasCalc::kCalcHillasSrc);
+//***     hcalc2.SetNameHillasSrc("MHillasSrcAnti");
+//***     hcalc2.SetNameSrcPosCam("MSrcPosAnti");
+//***     */
+//***
+//***    // Fill collection area and energy estimator (unfolding)
+//***    // Make sure to use the same binning for MHCollectionArea and MHEnergyEst
+//***    MHCollectionArea area;
+//***    area.SetHistAll(fSimpleMode ? hist : (TH2D&)mh1.GetHist());
+//***    MHEnergyEst      hest;
+//***
+//***    MFillH fill3(&area, "", "FillCollectionArea");
+//***    MFillH fill4(&hest, "", "FillEnergyEst");
+//***    MFillH fill5("MHThreshold", "", "FillThreshold");
+//***    fill3.SetWeight();
+//***    fill4.SetWeight();
+//***    fill5.SetWeight();
+//***    fill3.SetNameTab("ColArea");
+//***    fill4.SetNameTab("E-Est");
+//***    fill5.SetNameTab("Threshold");
+//***
+//***    MH3 hsize("MHillas.fSize");
+//***    hsize.SetName("ExcessMC");
+//***    hsize.Sumw2();
+//***
+//***    MBinning bins(size, "BinningExcessMC");
+//***    plist.AddToList(&hsize);
+//***    plist.AddToList(&bins);
+//***
+//***    MFillH fill1a("MHHillasMCPre  [MHHillas]",      "MHillas",      "FillHillasPre");
+//***    MFillH fill2a("MHHillasMCPost [MHHillas]",      "MHillas",      "FillHillasPost");
+//***    MFillH fill3a("MHVsSizeMCPost [MHVsSize]",      "MHillasSrc",   "FillVsSizePost");
+//***    MFillH fill4a("MHHilExtMCPost [MHHillasExt]",   "MHillasSrc",   "FillHilExtPost");
+//***    MFillH fill5a("MHHilSrcMCPost [MHHillasSrc]",   "MHillasSrc",   "FillHilSrcPost");
+//***    MFillH fill6a("MHImgParMCPost [MHImagePar]",    "MImagePar",    "FillImgParPost");
+//***    MFillH fill7a("MHNewParMCPost [MHNewImagePar]", "MNewImagePar", "FillNewParPost");
+//***    MFillH fill8a("ExcessMC       [MH3]",           "",             "FillExcessMC");
+//***    fill1a.SetNameTab("PreCut");
+//***    fill2a.SetNameTab("PostCut");
+//***    fill3a.SetNameTab("VsSize");
+//***    fill4a.SetNameTab("HilExt");
+//***    fill5a.SetNameTab("HilSrc");
+//***    fill6a.SetNameTab("ImgPar");
+//***    fill7a.SetNameTab("NewPar");
+//***    fill8a.SetBit(MFillH::kDoNotDisplay);
+//***    fill1a.SetWeight();
+//***    fill2a.SetWeight();
+//***    fill3a.SetWeight();
+//***    fill4a.SetWeight();
+//***    fill5a.SetWeight();
+//***    fill6a.SetWeight();
+//***    fill7a.SetWeight();
+//***    fill8a.SetWeight();
+//***>>>>>>> 1.17
+
+        if (!loop2.Eventloop(fMaxEvents))
+        {
+            *fLog << err << GetDescriptor() << ": Processing of MC-data failed." << endl;
+            return kFALSE;
+        }
+
+//***<<<<<<< MJSpectrum.cc
+        if (!loop2.GetDisplay())
+        {
+            *fLog << err << GetDescriptor() << ": Execution stopped by user." << endl;
+            return kFALSE;
+        }
+//***=======
+//***    tlist2.AddToList(&read);
+//***    if (!fRawMc && fNoThetaWeights)
+//***        tlist2.AddToList(&contsel);
+//***    //tlist2.AddToList(&calc);
+//***    //tlist2.AddToList(&hcalc1);
+//***    //tlist2.AddToList(&hcalc2);
+//***    tlist2.AddToList(&weight);
+//***    tlist2.AddToList(&fill1a);
+//***    tlist2.AddToList(fCut0);
+//***    tlist2.AddToList(fCut1);
+//***    tlist2.AddToList(fCut2);
+//***    tlist2.AddToList(fCut3);
+//***    tlist2.AddToList(&taskenv1);
+//***    tlist2.AddToList(&fill3);
+//***    tlist2.AddToList(&fill4);
+//***    tlist2.AddToList(&fill5);
+//***    tlist2.AddToList(&fill2a);
+//***    tlist2.AddToList(&fill3a);
+//***    tlist2.AddToList(&fill4a);
+//***    tlist2.AddToList(&fill5a);
+//***    tlist2.AddToList(&fill6a);
+//***    tlist2.AddToList(&fill7a);
+//***    tlist2.AddToList(&fill8a);
+//***    //tlist2.AddToList(&fill9a);
+//***
+//***    MEvtLoop loop2("FillMonteCarlo"); // ***** fName *****
+//***    loop2.SetParList(&plist);
+//***    loop2.SetDisplay(fDisplay);
+//***    loop2.SetLogStream(fLog);
+//***>>>>>>> 1.17
+
+        gLog.Separator("Energy Estimator");
+        if (plist.FindObject("EstimateEnergy"))
+            plist.FindObject("EstimateEnergy")->Print();
+
+        gLog.Separator("Spectrum");
+
+        // -------------------------- Spectrum ----------------------------
+
+
+        TArrayD res;
+        if (i==0)
+        {
+            // Calculate and display spectrum (N/TeVsm^2 at 1TeV)
+            res = DisplaySpectrum(area, excess, hest, ontime);
+
+            ((TH1D&)energy.GetHist()).Copy(energy0);
+        }
+        else
+        {
+            const Double_t scale = ontime*area.GetCollectionAreaAbs();
+            if (!DisplaySpectrumW(energy0, energy.GetHist(), hsize.GetHist(), weights, fSimpleMode ? hist : (TH2D&)mh1.GetHist(), scale, res))
+                return kFALSE;
+        }
+
+        // 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 = i==0 ? n/area.GetEntries() : 1;
+
+        // Print normalization constant
+        cout << "MC normalization factor:  " << scale << endl;
+
+        // Overlay normalized plots, and calculate weights for second loop
+        DisplaySize(plist, scale, weights);
+    }
 
     // check if output should be written
Index: trunk/MagicSoft/Mars/mjobs/MJSpectrum.h
===================================================================
--- trunk/MagicSoft/Mars/mjobs/MJSpectrum.h	(revision 7390)
+++ trunk/MagicSoft/Mars/mjobs/MJSpectrum.h	(revision 7403)
@@ -6,4 +6,5 @@
 #endif
 
+class TF1;
 class TH1;
 class TH1D;
@@ -33,4 +34,5 @@
     Bool_t fRawMc;
     Bool_t fNoThetaWeights;
+    Bool_t fWeighting;
 
     // Read Input
@@ -46,6 +48,8 @@
     void    DisplayResult(const TH2D &mh1) const;
     Bool_t  IntermediateLoop(MParList &plist, MH3 &h1, TH1D &temp1, const MDataSet &set, MMcSpectrumWeight &w) const;
+    TArrayD FitSpectrum(TH1D &spectrum) const;
     TArrayD DisplaySpectrum(MHCollectionArea &area, TH1D &excess, MHEnergyEst &hest, Double_t ontime) const;
-    Bool_t  DisplaySize(MParList &plist, Double_t scale) const;
+    Bool_t  DisplaySpectrumW(TH1D &energy0, TH1 &energy, TH1 &hsize, TH1D &weights, TH2D &hall, Double_t scale, TArrayD &res) const;
+    Bool_t  DisplaySize(MParList &plist, Double_t scale, TH1D &result) const;
     Bool_t  PlotSame(MStatusArray &arr, MParList &plist, const char *name, const char *tab, const char *plot, Double_t scale) const;
 
@@ -59,6 +63,9 @@
     void EnableSimpleMode(Bool_t b=kTRUE) { fSimpleMode=b; }
     void EnableRawMc(Bool_t b=kTRUE)      { fRawMc=b; }
+    void EnableWeighting(Bool_t b=kTRUE)  { fWeighting=b; }
 
     void SetEnergyEstimator(const MTask *task);
+
+    static TString FormString(const TF1 &f, Byte_t type=0);
 
     ClassDef(MJSpectrum, 0) // Proh'gram to calculate spectrum
