/* ======================================================================== *\
!
! *
! * 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 et al,  09/2002 <mailto:tbretz@astro.uni-wuerzburg.de>
!
!   Copyright: MAGIC Software Development, 2000-2002
!
!
\* ======================================================================== */


void estimate()
{
    //
    // This is a demonstration program which calculates the Hillas
    // parameter out of a Magic root file (raw data file).
    //

    //
    // Create a empty Parameter List and an empty Task List
    // The tasklist is identified in the eventloop by its name
    //
    MParList plist;


    MTaskList tlist;
    plist.AddToList(&tlist);

    //
    // Now setup the tasks and tasklist:
    // ---------------------------------
    //
    MReadMarsFile read("Events");
    read.DisableAutoScheme();
    /*
     read.AddFile("star.root");
     read.AddFile("star2.root");
     */
    read.AddFile("gammas.root");

    // Create a filter for the gamma events
    MFParticleId fgamma("MMcEvt", '=', kGAMMA);

    MTaskList tlist2;
    tlist2.SetFilter(&fgamma);

    MEnergyEstParam eest;
    eest.Add("MHillasSrc");

    //
    // Use this to change the binnign of the histograms to CT1-style
    //
    Bool_t usect1 = kTRUE;

    //
    // Here you set up the coefficients of the parametrization
    // (MEnergyEstParam)
    //
    TArrayD fA(5);
    fA[0] =  -2539; // [cm]
    fA[1] =    900; // [cm]
    fA[2] =   17.5; // [cm]
    fA[3] =      4;
    fA[4] =   58.3;

    TArrayD fB(7);
    fB[0] =    -8.64; // [GeV]
    fB[1] =   -0.069; // [GeV]
    fB[2] =  0.00066; // [GeV]
    fB[3] =    0.033; // [GeV]
    fB[4] = 0.000226; // [GeV]
    fB[5] =  4.14e-8; // [GeV]
    fB[6] =    -0.06;

    eest.SetCoeffA(fA);
    eest.SetCoeffB(fB);

    MH3 mh3e("MMcEvt.fEnergy",     "(MEnergyEst.fEnergy/MMcEvt.fEnergy-1)*(MEnergyEst.fEnergy/MMcEvt.fEnergy-1)");
    MH3 mh3i("MMcEvt.fImpact/100", "(MEnergyEst.fImpact/MMcEvt.fImpact-1)*(MEnergyEst.fImpact/MMcEvt.fImpact-1)");
    MH3 mh3eo("MMcEvt.fEnergy",     "MEnergyEst.fEnergy/MMcEvt.fEnergy-1");
    MH3 mh3io("MMcEvt.fImpact/100", "MEnergyEst.fImpact/MMcEvt.fImpact-1");

    MH3 mh3e2("MEnergyEst.fEnergy",     "(MEnergyEst.fEnergy/MMcEvt.fEnergy-1)*(MEnergyEst.fEnergy/MMcEvt.fEnergy-1)");
    MH3 mh3i2("MEnergyEst.fImpact/100", "(MEnergyEst.fImpact/MMcEvt.fImpact-1)*(MEnergyEst.fImpact/MMcEvt.fImpact-1)");
    MH3 mh3eo2("MEnergyEst.fEnergy",     "MEnergyEst.fEnergy/MMcEvt.fEnergy-1");
    MH3 mh3io2("MEnergyEst.fImpact/100", "MEnergyEst.fImpact/MMcEvt.fImpact-1");

    MH3 mhe("MMcEvt.fEnergy",     "MEnergyEst.fEnergy");
    MH3 mhi("MMcEvt.fImpact/100", "MEnergyEst.fImpact/100");

    mh3e.SetName("HistEnergy");
    mh3i.SetName("HistImpact");
    mh3eo.SetName("HistEnergyOffset");
    mh3io.SetName("HistImpactOffset");

    mh3e2.SetName("HistEnergy");
    mh3i2.SetName("HistImpact");
    mh3eo2.SetName("HistEnergyOffset");
    mh3io2.SetName("HistImpactOffset");

    mhe.SetName("HistEE");
    mhi.SetName("HistII");

    MFillH hfille(&mh3e);
    MFillH hfilli(&mh3i);
    MFillH hfilleo(&mh3eo);
    MFillH hfillio(&mh3io);

    MFillH hfille2(&mh3e2);
    MFillH hfilli2(&mh3i2);
    MFillH hfilleo2(&mh3eo2);
    MFillH hfillio2(&mh3io2);

    MFillH hfillee(&mhe);
    MFillH hfillii(&mhi);

    MBinning binsex("BinningHistEnergyX");
    MBinning binsey("BinningHistEnergyY");
    MBinning binsix("BinningHistImpactX");
    MBinning binsiy("BinningHistImpactY");
    MBinning binseox("BinningHistEnergyOffsetX");
    MBinning binseoy("BinningHistEnergyOffsetY");
    MBinning binsiox("BinningHistImpactOffsetX");
    MBinning binsioy("BinningHistImpactOffsetY");
    MBinning binseex("BinningHistEEX");
    MBinning binsiix("BinningHistIIX");
    MBinning binseey("BinningHistEEY");
    MBinning binsiiy("BinningHistIIY");

    binsex.SetEdgesLog(50, usect1 ? 300: 10, usect1 ? 50000 : 1e4);
    binsey.SetEdges(50, 0, 1.75);
    binseox.SetEdgesLog(50, usect1 ? 300 : 10, usect1 ? 50000 : 1e4);
    binseoy.SetEdges(50, -1.75, 1.75);

    binsix.SetEdges(50, 0, usect1 ? 450 : 300);
    binsiy.SetEdges(50, 0, 1.75);
    binsiox.SetEdges(50, 0, usect1 ? 450 : 300);
    binsioy.SetEdges(50, -1.75, 1.75);

    binseex.SetEdgesLog(50, usect1 ? 300 : 10, usect1 ? 50000 : 15e3);
    binseey.SetEdgesLog(50, usect1 ? 300 : 1,  usect1 ? 50000 : 2e3);
    binsiix.SetEdges(50, 0, usect1 ? 450 : 300);
    binsiiy.SetEdges(50, 0, usect1 ? 450 : 150);

    plist.AddToList(&binsex);
    plist.AddToList(&binsey);
    plist.AddToList(&binsix);
    plist.AddToList(&binsiy);
    plist.AddToList(&binseox);
    plist.AddToList(&binseoy);
    plist.AddToList(&binsiox);
    plist.AddToList(&binsioy);
    plist.AddToList(&binseex);
    plist.AddToList(&binseey);
    plist.AddToList(&binsiix);
    plist.AddToList(&binsiiy);

    //
    //  Setup tasklists
    //
    tlist.AddToList(&read);
    tlist.AddToList(&fgamma);
    tlist.AddToList(&tlist2);

    tlist2.AddToList(&eest);
    tlist2.AddToList(&hfille);
    tlist2.AddToList(&hfilli);
    tlist2.AddToList(&hfilleo);
    tlist2.AddToList(&hfillio);

    tlist2.AddToList(&hfille2);
    tlist2.AddToList(&hfilli2);
    tlist2.AddToList(&hfilleo2);
    tlist2.AddToList(&hfillio2);

    tlist2.AddToList(&hfillee);
    tlist2.AddToList(&hfillii);

    /*
     MPrint p("MEnergyEst");
     tlist2.AddToList(&p);
     */

    //
    // Create and setup the eventloop
    //
    MProgressBar bar;

    MEvtLoop evtloop;
    evtloop.SetProgressBar(&bar);
    evtloop.SetParList(&plist);

    //
    // Execute your analysis
    //
    if (!evtloop.Eventloop())
        return;

    tlist.PrintStatistics();

    const TString text = "\\sqrt{<y>}=%.0f%%";

    TCanvas *c=new TCanvas("Est1", "Estimates vs. E_{true}");
    c->Divide(2,2);
    c->cd(1);
    mh3i.DrawClone("PROFXnonew");
    TLatex *t = new TLatex(180, 1.6, Form(text, sqrt(mh3i.GetHist().GetMean(2))*100));
    t->Draw();
    c->cd(2);
    mh3e.DrawClone("PROFXnonew");
    t = new TLatex(2.7, 1.6, Form(text, sqrt(mh3e.GetHist().GetMean(2))*100));
    t->Draw();
    c->cd(3);
    mh3io.DrawClone("PROFXnonew");
    c->cd(4);
    mh3eo.DrawClone("PROFXnonew");

    c=new TCanvas("Est2", "Estimates vs. E_{est}");
    c->Divide(2,2);
    c->cd(1);
    mh3i2.DrawClone("PROFXnonew");
    c->cd(2);
    mh3e2.DrawClone("PROFXnonew");
    c->cd(3);
    mh3io2.DrawClone("PROFXnonew");
    c->cd(4);
    mh3eo2.DrawClone("PROFXnonew");

    mhe.DrawClone("PROFX");
    mhi.DrawClone("PROFX");
}