/* ======================================================================== *\
!
! *
! * This file is part of MARS, the MAGIC Analysis and Reconstruction
! * Software. It is distributed to you in the hope that it can be a useful
! * and timesaving tool in analysing Data of imaging Cerenkov telescopes.
! * It is distributed WITHOUT ANY WARRANTY.
! *
! * Permission to use, copy, modify and distribute this software and its
! * documentation for any purpose is hereby granted without fee,
! * provided that the above copyright notice appear in all copies and
! * that both that copyright notice and this permission notice appear
! * in supporting documentation. It is provided "as is" without express
! * or implied warranty.
! *
!
!
!   Author(s): Thomas Bretz, 4/2005 <mailto:tbretz@astro.uni-wuerzburg.de>
!
!   Copyright: MAGIC Software Development, 2000-2005
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//
//  MJSpectrum
//
// Program to calculate spectrum
//
/////////////////////////////////////////////////////////////////////////////
#include "MJSpectrum.h"

// Root
#include <TF1.h>
#include <TH1.h>
#include <TH2.h>
#include <TLine.h>
#include <TFile.h>
#include <TChain.h>
#include <TLatex.h>
#include <TCanvas.h>
#include <TObjArray.h>

// Environment
#include "MLog.h"
#include "MLogManip.h"

#include "MStatusArray.h"
#include "MStatusDisplay.h"

// Container
#include "MH3.h"
#include "MBinning.h"
#include "MDataSet.h"
#include "MMcCorsikaRunHeader.h"

// Spectrum
#include "../mhflux/MAlphaFitter.h"
#include "../mhflux/MHAlpha.h"
#include "../mhflux/MHCollectionArea.h"
//#include "../mhflux/MHThreshold.h"
#include "../mhflux/MHEnergyEst.h"
#include "../mhflux/MMcSpectrumWeight.h"

// Eventloop
#include "MEvtLoop.h"
#include "MTaskList.h"
#include "MParList.h"

// Tasks/Filter
#include "MReadMarsFile.h"
#include "MReadMarsFile.h"
#include "MFEventSelector2.h"
#include "MFDataMember.h"
#include "MEnergyEstimate.h"
#include "MTaskEnv.h"
#include "MFillH.h"
#include "MHillasCalc.h"
//#include "MSrcPosCalc.h"
#include "MContinue.h"

ClassImp(MJSpectrum);

using namespace std;

MJSpectrum::MJSpectrum(const char *name, const char *title)
    : fCutQ(0), fCut0(0),fCut1(0), fCut2(0), fCut3(0), fCutS(0),
    fEstimateEnergy(0), fCalcHadronness(0), fRefill(kFALSE),
    fSimpleMode(kTRUE), fRawMc(kFALSE), fNoThetaWeights(kFALSE)
{
    fName  = name  ? name  : "MJSpectrum";
    fTitle = title ? title : "Standard program to calculate spectrum";
}

MJSpectrum::~MJSpectrum()
{
    if (fCutQ)
        delete fCutQ;
    if (fCut0)
        delete fCut0;
    if (fCut1)
        delete fCut1;
    if (fCut2)
        delete fCut2;
    if (fCut3)
        delete fCut3;
    if (fCutS)
        delete fCutS;
    if (fEstimateEnergy)
        delete fEstimateEnergy;
    if (fCalcHadronness)
        delete fCalcHadronness;
}

// --------------------------------------------------------------------------
//
// Setup a task estimating the energy. The given task is cloned.
//
void MJSpectrum::SetEnergyEstimator(const MTask *task)
{
    if (fEstimateEnergy)
        delete fEstimateEnergy;
    fEstimateEnergy = task ? (MTask*)task->Clone() : 0;
}

Bool_t MJSpectrum::ReadTask(MTask* &task, const char *name, Bool_t mustexist) const
{
    if (task)
    {
        delete task;
        task = 0;
    }

    task = (MTask*)gFile->Get(name);
    if (!task)
    {
        if (!mustexist)
            return kTRUE;
        *fLog << err << dbginf << "ERROR - " << name << " doen't exist in file!" << endl;
        return kFALSE;
    }
    if (!task->InheritsFrom(MTask::Class()))
    {
        *fLog << err << dbginf << "ERROR - " << name << " read doesn't inherit from MTask!" << endl;
        delete task;
        return kFALSE;
    }

    task->SetName(name);

    if (dynamic_cast<MContinue*>(task))
        dynamic_cast<MContinue*>(task)->SetAllowEmpty();

    return kTRUE;
}

void MJSpectrum::PrintSetup(const MAlphaFitter &fit) const
{
    fLog->Separator("Alpha Fitter");
    *fLog << all;
    fit.Print("result");

    fLog->Separator("Used Cuts");
    fCutQ->Print();
    fCut0->Print();
    fCut1->Print();
    fCutS->Print();
    fCut2->Print();
    fCut3->Print();
}

// --------------------------------------------------------------------------
//
// 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");

    // Read Resources
    if (!CheckEnv(w))
    {
        *fLog << err << "ERROR - Reading resources for MMcSpectrumWeight failed." << endl;
        return kFALSE;
    }

    // Read the first corsika RunHeader from the MC files
    TChain chain("RunHeaders");
    if (!set.AddFilesOn(chain))
        return kFALSE;

    MMcCorsikaRunHeader *h=0;
    chain.SetBranchAddress("MMcCorsikaRunHeader.", &h);
    chain.GetEntry(1);

    if (!h)
    {
        *fLog << err << "ERROR - Couldn't read MMcCorsikaRunHeader from DataSet." << endl;
        return kFALSE;
    }

    // Propagate the run header to MMcSpectrumWeight
    if (!w.Set(*h))
    {
        *fLog << err << "ERROR - Initializing MMcSpectrumWeight failed." << endl;
        return kFALSE;
    }

    // Print new setup of MMcSpectrumWeight
    w.Print();
    return kTRUE;
}

Float_t MJSpectrum::ReadInput(MParList &plist, TH1D &h1, TH1D &h2)
{
    *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 -1;
    }

    MStatusArray arr;
    if (arr.Read()<=0)
    {
        *fLog << "MStatusDisplay not found in file... abort." << endl;
        return -1;
    }

    TH1D *vstime = (TH1D*)arr.FindObjectInCanvas("Theta",  "TH1D", "OnTime");
    TH1D *size   = (TH1D*)arr.FindObjectInCanvas("Excess", "TH1D", "Hist");
    if (!vstime || !size)
        return -1;

    vstime->Copy(h1);
    size->Copy(h2);
    h1.SetDirectory(0);
    h2.SetDirectory(0);

    if (fDisplay)
        arr.DisplayIn(*fDisplay, "Hist");

    if (!ReadTask(fCutQ, "CutQ"))
        return -1;
    if (!ReadTask(fCut0, "Cut0"))
        return -1;
    if (!ReadTask(fCut1, "Cut1"))
        return -1;
    if (!ReadTask(fCut2, "Cut2"))
        return -1;
    if (!ReadTask(fCut3, "Cut3"))
        return -1;
    if (!ReadTask(fCalcHadronness, "CalcHadronness", kFALSE))
        return -1;

    TObjArray arrread;

    TIter Next(plist);
    TObject *o=0;
    while ((o=Next()))
        if (o->InheritsFrom(MBinning::Class()))
            arrread.Add(o);

    arrread.Add(plist.FindObject("MAlphaFitter"));

    if (!ReadContainer(arrread))
        return -1;

    return vstime->Integral();
}

Bool_t MJSpectrum::ReadOrigMCDistribution(const MDataSet &set, TH1 &h, MMcSpectrumWeight &weight) const
{
    // Some debug output
    fLog->Separator("Compiling original MC distribution");

    // The name of the input container is MMcEvtBasic
    weight.SetNameMcEvt("MMcEvtBasic");
    // Get the corresponding rule for the weights
    const TString w(weight.GetFormulaWeights());
    // Set back to MMcEvt
    weight.SetNameMcEvt();

    *fLog << inf << "Using weights: " << w << endl;
    *fLog << "Please stand by, this may take a while..." << flush;

    if (fDisplay)
        fDisplay->SetStatusLine1("Compiling MC distribution...");

    // Create chain
    TChain chain("OriginalMC");
    if (!set.AddFilesOn(chain))
        return kFALSE;

    // Prepare histogram
    h.Reset();
    h.Sumw2();

    // Fill histogram from chain
    h.SetDirectory(gROOT);
    if (h.InheritsFrom(TH2::Class()))
    {
        h.SetNameTitle("ThetaEMC", "Event-Distribution vs Theta and Energy for MC (produced)");
        h.SetXTitle("\\Theta [\\circ]");
        h.SetYTitle("E [GeV]");
        h.SetZTitle("Counts");
        chain.Draw("MMcEvtBasic.fEnergy:MMcEvtBasic.fTelescopeTheta*TMath::RadToDeg()>>ThetaEMC", w, "goff");
    }
    else
    {
        h.SetNameTitle("ThetaMC", "Event-Distribution vs Theta for MC (produced)");
        h.SetXTitle("\\Theta [\\circ]");
        h.SetYTitle("Counts");
        chain.Draw("MMcEvtBasic.fTelescopeTheta*TMath::RadToDeg()>>ThetaMC", w, "goff");
    }
    h.SetDirectory(0);

    *fLog << "done." << endl;
    if (fDisplay)
        fDisplay->SetStatusLine2("done.");

    if (h.GetEntries()>0)
        return kTRUE;

    *fLog << err << "ERROR - Histogram with original MC distribution empty..." << endl;

    return h.GetEntries()>0;
}

Bool_t MJSpectrum::GetThetaDistribution(TH1D &temp1, TH1D &temp2) const
{
    // Display some stuff
    if (fDisplay)
    {
        TCanvas &c = fDisplay->AddTab("ZdDist");
        c.Divide(2,2);

        // On-Time vs. Theta
        c.cd(1);
        gPad->SetBorderMode(0);
        temp1.DrawCopy();

        // Number of MC events (produced) vs Theta
        c.cd(2);
        gPad->SetBorderMode(0);
        temp2.SetName("NVsTheta");
        temp2.DrawCopy("hist");

        c.cd(4);
        gPad->SetBorderMode(0);

        c.cd(3);
        gPad->SetBorderMode(0);
    }

    // Calculate the Probability
    temp1.Divide(&temp2);
    temp1.Scale(fNoThetaWeights ? 1./temp1.GetMaximum() : 1./temp1.Integral());

    // Some cosmetics: Name, Axis, etc.
    temp1.SetName("ProbVsTheta");
    temp1.SetTitle("Probability vs. Zenith Angle to choose MC events");
    temp1.SetYTitle("Probability");
    if (fDisplay)
        temp1.DrawCopy("hist");

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// Display the final theta distribution.
//
Bool_t MJSpectrum::DisplayResult(const TH2D &h2) const
{
    if (!fDisplay || !fDisplay->CdCanvas("ZdDist"))
    {
        *fLog << err << "ERROR - Display or tab ZdDist vanished... abort." << endl;
        return kFALSE;
    }

    TH1D &proj = *h2.ProjectionX();
    proj.SetNameTitle("ThetaFinal", "Final Theta Distribution");
    proj.SetXTitle("\\Theta [\\circ]");
    proj.SetYTitle("Counts");
    proj.SetLineColor(kBlue);
    proj.SetDirectory(0);
    proj.SetBit(kCanDelete);

    TVirtualPad *pad = gPad;

    pad->cd(4);
    proj.DrawCopy();

    pad->cd(1);
    TH1D *theta = (TH1D*)gPad->FindObject("Theta");
    if (!theta)
    {
        *fLog << err << "ERROR - Theta-Histogram vanished... cannot proceed." << endl;
        return kFALSE;
    }

    if (proj.GetMaximum()==0)
    {
        *fLog << err;
        *fLog << "ERROR - The Zenith Angle distribution of your Monte Carlos doesn't overlap";
        *fLog << "        with the the Zenith Angle distribution of your observation.";
        theta->SetLineColor(kRed);
        return kFALSE;;
    }

    proj.Scale(theta->GetMaximum()/proj.GetMaximum());
    theta->SetMaximum(1.05*TMath::Max(theta->GetMaximum(), proj.GetMaximum()));

    proj.Draw("same");

    // Compare both histograms
    const Double_t prob = proj.Chi2Test(theta, "");
    if (prob==1)
        return kTRUE;

    if (prob>0.99)
    {
        *fLog << inf;
        *fLog << "The Zenith Angle distribution of your Monte Carlos fits well";
        *fLog << "with the the Zenith Angle distribution of your observation.";
        return kTRUE;
    }

    // <0.01: passen gar nicht
    // ==1: passen perfekt
    // >0.99: passer sehr gut

    if (prob<0.01)
    {
        *fLog << err;
        *fLog << "ERROR - The Zenith Angle distribution of your Monte Carlos does not fit";
        *fLog << "        with the the Zenith Angle distribution of your observation.";
        theta->SetLineColor(kRed);
        return kFALSE;
    }

    *fLog << warn;
    *fLog << "WARNING - The Zenith Angle distribution of your Monte Carlos doesn't fits well";
    *fLog << "          with the the Zenith Angle distribution of your observation.";
    return kTRUE;
}

// --------------------------------------------------------------------------
//
// 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*/
{
    // Try to find the class used to determin the signal!
    TString cls("MHAlpha");
    if (fDisplay)
    {
        TCanvas *c = fDisplay->GetCanvas("Hist");
        if (c)
        {
            TIter Next(c->GetListOfPrimitives());
            TObject *obj=0;
            while ((obj=Next()))
                if (obj->InheritsFrom(MHAlpha::Class()))
                    break;
            if (obj)
                cls = obj->ClassName();
        }
    }

    // Now fill the histogram
    *fLog << endl;
    fLog->Separator("Refill Excess");
    *fLog << endl;

    MTaskList tlist;
    plist.AddToList(&tlist);

    MReadTree read("Events");
    read.DisableAutoScheme();
    read.AddFile(fPathIn);

    MTaskEnv taskenv0("CalcHadronness");
    taskenv0.SetDefault(fCalcHadronness);

    MEnergyEstimate est;
    MTaskEnv taskenv1("EstimateEnergy");
    taskenv1.SetDefault(fEstimateEnergy ? fEstimateEnergy : &est);

    MContinue *cont = new MContinue("", "CutS");
    cont->SetAllowEmpty();

    if (fCutS)
        delete fCutS;
    fCutS = cont;

    // FIXME: Create HistE and HistEOff to be able to modify it from
    // the resource file.

    MFillH fill1(Form("HistEOff [%s]", cls.Data()), "", "FillHistEOff");
    MFillH fill2(Form("HistE    [%s]", cls.Data()), "", "FillHistE");

    MFDataMember f0("DataType.fVal", '<', 0.5, "FilterOffData");
    MFDataMember f1("DataType.fVal", '>', 0.5, "FilterOnData");

    fill1.SetFilter(&f0);
    fill2.SetFilter(&f1);

    tlist.AddToList(&read);
    //tlist.AddToList(&taskenv0); // not necessary, stored in file!
    tlist.AddToList(fCutS);
    tlist.AddToList(&taskenv1);
    tlist.AddToList(&f0);
    tlist.AddToList(&f1);
    tlist.AddToList(&fill1);
    tlist.AddToList(&fill2);

    // by setting it here it is distributed to all consecutive tasks
    tlist.SetAccelerator(MTask::kAccDontReset|MTask::kAccDontTime);

    MEvtLoop loop("RefillExcess"); // ***** fName *****
    loop.SetParList(&plist);
    loop.SetDisplay(fDisplay);
    loop.SetLogStream(fLog);

    if (!SetupEnv(loop))
        return kFALSE;

    if (!loop.Eventloop())
    {
        *fLog << err << GetDescriptor() << ": Refilling of data failed." << endl;
        return kFALSE;
    }

    if (!loop.GetDisplay())
    {
        *fLog << err << GetDescriptor() << ": Execution stopped by user." << endl;
        return kFALSE;
    }

    const MHAlpha *halpha = (MHAlpha *)plist.FindObject("HistE");
    if (!halpha)
    {
        *fLog << err << GetDescriptor() << ": HistE [MHAlpha] not found... abort." << endl;
        return kFALSE;
    }

    halpha->GetHEnergy().Copy(h2);
    h2.SetDirectory(0);

    return kTRUE;
}

Bool_t MJSpectrum::IntermediateLoop(MParList &plist, MH3 &mh1, TH1D &temp1, const MDataSet &set, MMcSpectrumWeight &weight) const
{
    MTaskList tlist1;
    plist.Replace(&tlist1);

    MReadMarsFile readmc("OriginalMC");
    //readmc.DisableAutoScheme();
    if (!set.AddFilesOn(readmc))
        return kFALSE;

    readmc.EnableBranch("MMcEvtBasic.fTelescopeTheta");
    readmc.EnableBranch("MMcEvtBasic.fEnergy");

    mh1.SetLogy();
    mh1.SetLogz();
    mh1.SetName("ThetaE");

    MFillH fill0(&mh1);
    //fill0.SetDrawOption("projx only");

    MBinning *bins2 = (MBinning*)plist.FindObject("BinningEnergyEst");
    MBinning *bins3 = (MBinning*)plist.FindObject("BinningTheta");
    if (bins2 && bins3)
    {
        bins2->SetName("BinningThetaEY");
        bins3->SetName("BinningThetaEX");
    }
    tlist1.AddToList(&readmc);
    tlist1.AddToList(&weight);

    temp1.SetXTitle("MMcEvtBasic.fTelescopeTheta*kRad2Deg");
    MH3 mh3mc(temp1);

    MFEventSelector2 sel1(mh3mc);
    sel1.SetHistIsProbability();

    fill0.SetFilter(&sel1);

    if (!fRawMc)
        tlist1.AddToList(&sel1);
    tlist1.AddToList(&fill0);

    // by setting it here it is distributed to all consecutive tasks
    tlist1.SetAccelerator(MTask::kAccDontReset|MTask::kAccDontTime);

    MEvtLoop loop1("IntermediateLoop"); // ***** fName *****
    loop1.SetParList(&plist);
    loop1.SetLogStream(fLog);
    loop1.SetDisplay(fDisplay);

    if (!SetupEnv(loop1))
        return kFALSE;

    if (!loop1.Eventloop(fMaxEvents))
    {
        *fLog << err << GetDescriptor() << ": Processing of MC-data failed." << endl;
        return kFALSE;
    }

    if (!loop1.GetDisplay())
    {
        *fLog << err << GetDescriptor() << ": Execution stopped by user." << endl;
        return kFALSE;
    }

    if (bins2 && bins3)
    {
        bins2->SetName("BinningEnergyEst");
        bins3->SetName("BinningTheta");
    }

    return kTRUE;
}

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");

    TString str = FormString(f);

    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;
}

// --------------------------------------------------------------------------
//
// 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
{
    // Create copies of the histograms
    TH1D collarea(area.GetHEnergy());
    TH1D spectrum(excess);
    TH1D weights;

    // Get spill-over corrections from energy estimation
    hest.GetWeights(weights);

    // Print effective on-time
    cout << "Effective On time: " << ontime << "s" << endl;

    // Setup spectrum plot
    spectrum.SetNameTitle("Preliminary", "N/sm^{2} versus Energy (before unfolding)");
    spectrum.SetYTitle("N/sm^{2}");
    spectrum.SetDirectory(NULL);
    spectrum.SetBit(kCanDelete);

    // Divide by collection are and on-time
    spectrum.Scale(1./ontime);
    spectrum.Divide(&collarea);

    // Draw spectrum before applying spill-over corrections
    TCanvas &c1 = fDisplay->AddTab("Spectrum");
    c1.Divide(2,2);
    c1.cd(1);
    gPad->SetBorderMode(0);
    gPad->SetLogx();
    gPad->SetLogy();
    gPad->SetGridx();
    gPad->SetGridy();
    collarea.DrawCopy();

    c1.cd(2);
    gPad->SetBorderMode(0);
    gPad->SetLogx();
    gPad->SetLogy();
    gPad->SetGridx();
    gPad->SetGridy();
    TH1D *spec=(TH1D*)spectrum.DrawCopy();
    FitSpectrum(*spec);

    c1.cd(3);
    gPad->SetBorderMode(0);
    gPad->SetLogx();
    gPad->SetLogy();
    gPad->SetGridx();
    gPad->SetGridy();
    weights.DrawCopy();

    // Apply spill-over correction (done't take the errors into account)
    // They are supposed to be identical with the errors of the
    // collection area and cancel out.
    //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);
    }

    c1.cd(4);
    gPad->SetBorderMode(0);
    gPad->SetLogx();
    gPad->SetLogy();
    gPad->SetGridx();
    gPad->SetGridy();
    spectrum.SetMinimum(1e-12);
    spectrum.SetXTitle("E [GeV]");
    spectrum.SetYTitle("dN/dE [N/TeVsm^{2}]");
    spec = (TH1D*)spectrum.DrawCopy();

    // Calculate Spectrum * E^2
    for (int i=0; i<spectrum.GetNbinsX(); i++)
    {
        const Double_t e = TMath::Sqrt(spectrum.GetBinLowEdge(i+1)*spectrum.GetBinLowEdge(i+2))*1e-3;

        spectrum.SetBinContent(i+1, spectrum.GetBinContent(i+1)*e*e);
        spectrum.SetBinError(  i+1, spectrum.GetBinError(i+1)  *e*e);
    }

    // Display dN/dE*E^2 for conveinience
    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);
    spectrum.Fit(&f, "NIM", "", 100, 5000);
    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;
    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);

    return res;
  */
/*
    // 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");
 */
}

// --------------------------------------------------------------------------
//
// 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
{
    TString same(name);
    same += "Same";

    TH1 *h1  = (TH1*)arr.FindObjectInCanvas(name, "TH1F", tab);
    TH1 *h2  = (TH1*)arr.FindObjectInCanvas(same, "TH1F", tab);
    if (!h1 || !h2)
        return kFALSE;

    TObject *obj = plist.FindObject(plot);
    if (!obj)
    {
        *fLog << warn << plot << " not in parameter list... skipping." << endl;
        return kFALSE;
    }

    TH1 *h3  = (TH1*)obj->FindObject(name);
    if (!h3)
    {
        *fLog << warn << name << " not found in " << plot << "... skipping." << endl;
        return kFALSE;
    }


    const MAlphaFitter *fit = (MAlphaFitter*)plist.FindObject("MAlphaFitter");
    const Double_t ascale = fit ? fit->GetScaleFactor() : 1;

    gPad->SetBorderMode(0);
    h2->SetLineColor(kBlack);
    h3->SetLineColor(kBlue);
    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()));

    h2 = h2->DrawCopy();
    h3 = h3->DrawCopy("same");

    // Don't do this on the original object!
    h2->SetStats(kFALSE);
    h3->SetStats(kFALSE);

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// Take a lot of histograms and plot them together in one plot.
// Calls PlotSame
//
Bool_t MJSpectrum::DisplaySize(MParList &plist, Double_t scale) 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;
    }

    TH1 *excess = (TH1D*)arr.FindObjectInCanvas("Excess", "TH1D", "Hist");
    if (!excess)
        return kFALSE;

    // ------------------- Plot excess versus size ------------------- 

    TCanvas &c = fDisplay->AddTab("Excess");
    c.Divide(3,2);
    c.cd(1);
    gPad->SetBorderMode(0);
    gPad->SetLogx();
    gPad->SetLogy();
    gPad->SetGridx();
    gPad->SetGridy();

    excess->SetTitle("Excess events vs Size (data/black, mc/blue)");
    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("ExcessMC")))
    {
        TH1 *histsel = (TH1F*)o->FindObject("");
        if (histsel)
        {
            if (scale<0)
                scale = excess->Integral()/histsel->Integral();

            histsel->Scale(scale);
            histsel->SetLineColor(kBlue);
            histsel->SetBit(kCanDelete);
            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");
            const Double_t p = histsel->Chi2Test(excess, "P");

            TLatex tex;
            tex.SetBit(TLatex::kTextNDC);
            tex.DrawLatex(0.75, 0.93, Form("P(\\chi^{2})=%.0f%%", p*100));
        }
    }

    // -------------- Comparison of Image Parameters --------------
    c.cd(2);
    PlotSame(arr, plist, "Dist",   "HilSrc",  "MHHilSrcMCPost", scale);

    c.cd(3);
    PlotSame(arr, plist, "Length", "PostCut", "MHHillasMCPost", scale);

    c.cd(4);
    PlotSame(arr, plist, "M3l",    "HilExt",  "MHHilExtMCPost", scale);

    c.cd(5);
    PlotSame(arr, plist, "Conc1",  "NewPar",  "MHNewParMCPost", scale);

    c.cd(6);
    PlotSame(arr, plist, "Width",  "PostCut", "MHHillasMCPost", scale);

    return kTRUE;
}

// --------------------------------------------------------------------------
//
void MJSpectrum::DisplayCutEfficiency(const MHCollectionArea &area0, const MHCollectionArea &area1) const
{
    if (!fDisplay)
        return;

    const TH1D &trig = area0.GetHEnergy();
    TH1D &cut = (TH1D&)*area1.GetHEnergy().Clone();

    fDisplay->AddTab("CutEff");

    gPad->SetBorderMode(0);
    gPad->SetFrameBorderMode(0);
    gPad->SetLogx();
    gPad->SetGridx();
    gPad->SetGridy();

    cut.Divide(&trig);
    cut.Scale(100);
    cut.SetNameTitle("CutEff", "Background Supression: Cut efficiency (after star)");
    cut.SetYTitle("\\eta [%]");
    cut.SetDirectory(0);
    cut.SetMinimum(0);
    cut.SetMaximum(100);
    cut.SetBit(kCanDelete);
    cut.Draw();

    TLine line;
    line.SetLineColor(kBlue);
    line.SetLineWidth(2);
    line.SetLineStyle(kDashed);
    line.DrawLine(cut.GetBinLowEdge(1), 50, cut.GetBinLowEdge(cut.GetNbinsX()+1), 50);
}

Bool_t MJSpectrum::Process(const MDataSet &set)
{
    if (!set.IsValid())
    {
        *fLog << err << "ERROR - DataSet invalid!" << endl;
        return kFALSE;
    }

    CheckEnv();

    // --------------------------------------------------------------------------------

    *fLog << inf;
    fLog->Separator(GetDescriptor());
    *fLog << "Compile Monte Carlo Sample (data set " << set.GetName() << ")" << endl;
    *fLog << endl;

    MBinning bins1("BinningAlpha");
    MBinning bins2("BinningEnergyEst");
    MBinning bins3("BinningTheta");
    MBinning bins4("BinningFalseSource");
    MBinning bins5("BinningWidth");
    MBinning bins6("BinningLength");
    MBinning bins7("BinningDist");
    MBinning bins8("BinningMaxDist");
    MBinning bins9("BinningM3Long");
    MBinning bins0("BinningConc1");

    MAlphaFitter fit;

    MParList plist;
    plist.AddToList(&bins1);
    plist.AddToList(&bins3);
    plist.AddToList(&bins4);
    plist.AddToList(&bins5);
    plist.AddToList(&bins6);
    plist.AddToList(&bins7);
    plist.AddToList(&bins8);
    plist.AddToList(&bins9);
    plist.AddToList(&bins0);
    plist.AddToList(&fit);

    TH1D temp1, size;
    const Float_t ontime = ReadInput(plist, temp1, size);
    if (ontime<0)
    {
        *fLog << err << GetDescriptor() << ": Could not determin effective on time..." << endl;
        return kFALSE;
    }

    plist.AddToList(&bins2);

    // Initialize weighting to a new spectru
    // m as defined in the resource file
    MMcSpectrumWeight weight;
    if (!InitWeighting(set, weight))
        return kFALSE;

    bins3.SetEdges(temp1, 'x');

    // Read the MC distribution as produced by corsika into
    // temp2 (1D) and apply the weights previously determined
    TH1D temp2(temp1);
    if (!ReadOrigMCDistribution(set, temp2, weight))
        return kFALSE;

    // Calculate the weights according to the zenith angle distribution
    // of the raw-data which have to be applied to the MC events
    if (!GetThetaDistribution(temp1, temp2))
        return kFALSE;

    // Tell the weighting task about the ZA-weights
    if (!fNoThetaWeights)
        weight.SetWeightsZd(&temp1);

    // Refill excess histogram to determin the excess events
    TH1D excess;
    if (!Refill(plist, excess))
        return kFALSE;

    // Print the setup and result of the MAlphaFitter
    // Print used cuts
    PrintSetup(fit);

    TH2D hist;
    MH3 mh1("MMcEvtBasic.fTelescopeTheta*kRad2Deg", "MMcEvtBasic.fEnergy");
    if (fSimpleMode)
    {
        // Now we read the MC distribution as produced by corsika again
        // with the spectral weights applied into a histogram vs.
        // zenith angle and energy
        hist.UseCurrentStyle();
        MH::SetBinning(&hist, &bins3/*temp1.GetXaxis()*/, &bins2/*excess.GetXaxis()*/);
        if (!ReadOrigMCDistribution(set, hist, weight))
            return kFALSE;

        // No we apply the the zenith-angle-weights to the corsika produced
        // MC distribution. Unfortunately this must be done manually
        // because we are multiplying column by column
        if (!fRawMc)
        {
            for (int y=0; y<hist.GetNbinsY(); y++)
                for (int x=0; x<hist.GetNbinsX(); x++)
                {
                    hist.SetBinContent(x, y, hist.GetBinContent(x, y)*temp1.GetBinContent(x));
                    hist.SetBinError(x, y,   hist.GetBinError(x, y)  *temp1.GetBinContent(x));
                }
        }
    }
    else
    {
        // This rereads the original MC spectrum and aplies both
        // weights, spectral weights and ZA-weights.
        weight.SetNameMcEvt("MMcEvtBasic");
        if (!IntermediateLoop(plist, mh1, temp1, set, weight))
            return kFALSE;
        weight.SetNameMcEvt();
    }

    if (!DisplayResult(fSimpleMode ? hist : (TH2D&)mh1.GetHist()))
        return kFALSE;

    // ------------------------- Final loop --------------------------

    *fLog << endl;
    fLog->Separator("Calculate efficiencies");
    *fLog << endl;

    MTaskList tlist2;
    plist.Replace(&tlist2);

    MReadMarsFile read("Events");
    read.DisableAutoScheme();
    if (!set.AddFilesOn(read))
        return kFALSE;

    // 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 area0("TriggerArea");
    MHCollectionArea area1;
    area0.SetHistAll(fSimpleMode ? hist : (TH2D&)mh1.GetHist());
    area1.SetHistAll(fSimpleMode ? hist : (TH2D&)mh1.GetHist());

    MHEnergyEst      hest;

    MFillH fill30(&area0, "", "FillTriggerArea");
    MFillH fill31(&area1, "", "FillCollectionArea");
    MFillH fill4(&hest,   "", "FillEnergyEst");
    MFillH fill5("MHThreshold", "", "FillThreshold");
    fill30.SetWeight();
    fill31.SetWeight();
    fill4.SetWeight();
    fill5.SetWeight();
    fill30.SetNameTab("TrigArea");
    fill31.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();

    MTaskEnv taskenv0("CalcHadronness");
    taskenv0.SetDefault(fCalcHadronness);

    MEnergyEstimate est;
    MTaskEnv taskenv1("EstimateEnergy");
    taskenv1.SetDefault(fEstimateEnergy ? fEstimateEnergy : &est);

    tlist2.AddToList(&read);
    // If no weighting should be applied but the events should
    // be thrown away according to the theta distribution
    // it is enabled here
    if (!fRawMc && fNoThetaWeights)
        tlist2.AddToList(&contsel);
    //tlist2.AddToList(&calc);
    //tlist2.AddToList(&hcalc1);
    //tlist2.AddToList(&hcalc2);
    tlist2.AddToList(&weight);
    tlist2.AddToList(&fill1a);
    tlist2.AddToList(&fill30);
    tlist2.AddToList(fCutQ);
    tlist2.AddToList(fCut0);
    tlist2.AddToList(&taskenv0);
    tlist2.AddToList(fCutS);
    tlist2.AddToList(fCut1);
    tlist2.AddToList(fCut2);
    tlist2.AddToList(fCut3);
    tlist2.AddToList(&taskenv1);
    tlist2.AddToList(&fill31);
    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);

    // by setting it here it is distributed to all consecutive tasks
    tlist2.SetAccelerator(MTask::kAccDontReset|MTask::kAccDontTime);

    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(area1, 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 *= area1.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/area1.GetEntries();

    // Print normalization constant
    cout << "MC normalization factor:  " << scale << endl;

    // Display cut efficiency
    DisplayCutEfficiency(area0, area1);

    // Overlay normalized plots
    DisplaySize(plist, scale);

    // check if output should be written
    if (fPathOut.IsNull())
        return kTRUE;

    // Write the output
    TObjArray cont;
    cont.Add((TObject*)GetEnv());

    if (fDisplay)
        cont.Add(fDisplay);

    return WriteContainer(cont, "", "RECREATE");
}
