/* ======================================================================== *\
!
! *
! * 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 estct1()
{
    //
    // 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:
    // ---------------------------------
    //
    MReadTree read("Events", "~/ct1/MC_ON2.root");

    //MReadMarsFile read("Events");
    read.DisableAutoScheme();
    /*
     read.AddFile("star.root");
     read.AddFile("star2.root");
     */
    //read.AddFile("~/ct1/MC_ON2.root");

    MEnergyEstParam eest("Hillas");
    eest.Add("HillasSrc");

    //
    // 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] = -3907.74; //4916.4;     //-2414.75;
    fA[1] = 1162.3; //149.549;    // 1134.28;
    fA[2] = 199.351;//-558.209;   // 132.932;
    fA[3] = 0.403192;//0.270725;   //0.292845;
    fA[4] = 121.921;//107.001;    // 107.001;

    TArrayD fB(7);
    fB[0] =  677.821;//-8234.12;  //-4282.25;
    fB[1] =  11.3302;//23.2153;   //  18.892;
    fB[2] =  -0.0211177;//0.416372;  //0.193373;
    fB[3] =  -23.0217;//332.42;    //203.803;
    fB[4] =  -0.785242;//-0.701764; //-0.534876;
    fB[5] = 5.6413e-5;//-0.0131774; //-0.00789539;
    fB[6] =  -0.146595;//-0.162687; //  0.111913;

    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, usect1 ? 0.8 : 1.75);
    binseox.SetEdgesLog(50, usect1 ? 300 : 10, usect1 ? 50000 : 1e4);
    binseoy.SetEdges(50, usect1 ? -0.75 : -1.75, usect1 ? 0.75 : 1.75);

    binsix.SetEdges(50, 0, usect1 ? 275 : 300);
    binsiy.SetEdges(50, 0, usect1 ? 0.2 : 1.75);
    binsiox.SetEdges(50, 0, usect1 ? 275 : 300);
    binsioy.SetEdges(50, usect1 ? -0.75 : -1.75, usect1 ? 0.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 ? 275 : 300);
    binsiiy.SetEdges(50, 0, usect1 ? 275 : 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(&eest);
    tlist.AddToList(&hfille);
    tlist.AddToList(&hfilli);
    tlist.AddToList(&hfilleo);
    tlist.AddToList(&hfillio);

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

    tlist.AddToList(&hfillee);
    tlist.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%%";

    char txt[1000];

    TCanvas *c=new TCanvas("Est1", "Estimates vs. E_{true}");
    c->Divide(2,2);
    c->cd(1);
    mh3i.DrawClone("PROFXnonew");
    sprintf(txt, text.Data(), sqrt(mh3i.GetHist().GetMean(2))*100);
    TLatex *t = new TLatex(180, 0.15, txt);
    t->Draw();
    c->cd(2);
    mh3e.DrawClone("PROFXnonew");
    sprintf(txt, text.Data(), sqrt(mh3e.GetHist().GetMean(2))*100);
    t = new TLatex(3.5, 0.6, txt);
    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");
}
