source: trunk/MagicSoft/Mars/macros/estimate.C@ 1574

/* ======================================================================== *\
!
! *
! * 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("HillasSource");

    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;
    /* fb5=size*width, ir 29.8% (RMS=23%)
     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;*/
    /* size*dist, ir 30%
     fB[0] =   -32; // [GeV]
     fB[1] =  -0.0295; // [GeV]
     fB[2] =   0.00112; // [GeV]
     fB[3] =   0.0644; // [GeV]
     fB[4] = 1e-5; // [GeV]
     fB[5] =  6.3e-6; // [GeV]
     fB[6] =   -0.06;
     */
    /* size*width, ir 30%
     fB[0] =   -17; // [GeV]
     fB[1] =  -0.06; // [GeV]
     fB[2] =   0.0013; // [GeV]
     fB[3] =   0.039; // [GeV]
     fB[4] = 0.00036; // [GeV]
     fB[5] =  5.6e-6; // [GeV]
     fB[6] =   -0.06;
     */
    /* fb5=width, ir 31%
     fB[0] =   0.14; // [GeV]
     fB[1] =  -0.078; // [GeV]
     fB[2] =   0.0003; // [GeV]
     fB[3] =  -0.171; // [GeV]
     fB[4] = 0.00038; // [GeV]
     fB[5] =  3e-5; // [GeV]
     fB[6] =   -0.06;
     */
    /* size*length, ir  32%
     fB[0] =   -1.14; // [GeV]
     fB[1] =  -0.075; // [GeV]
     fB[2] =   0.0008; // [GeV]
     fB[3] =  -0.148; // [GeV]
     fB[4] = 0.00037; // [GeV]
     fB[5] =  8e-7; // [GeV]
     fB[6] =   -0.06;
     */
    /* ir/dist 35%
     fB[0] =     7.7; // [GeV]
     fB[1] =  -0.076; // [GeV]
     fB[2] =   0.011; // [GeV]
     fB[3] =  -0.175; // [GeV]
     fB[4] = 0.00038; // [GeV]
     fB[5] = -0.0001; // [GeV]
     fB[6] =   -0.06;
     */


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

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

    mh3e.SetName("HistEnergy");
    mh3i.SetName("HistImpact");

    MFillH hfille(&mh3e);
    MFillH hfilli(&mh3i);

    MBinning binsex("BinningHistEnergyX");
    MBinning binsey("BinningHistEnergyY");
    MBinning binsix("BinningHistImpactX");
    MBinning binsiy("BinningHistImpactY");

    binsex.SetEdgesLog(50, 10, 1e4);
    binsey.SetEdges(50, 0, 7);

    binsix.SetEdges(50, 0, 300);
    binsiy.SetEdges(50, 0, 3);

    plist.AddToList(&binsex);
    plist.AddToList(&binsey);
    plist.AddToList(&binsix);
    plist.AddToList(&binsiy);

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

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

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

    //
    // Create and setup the eventloop
    //
    MEvtLoop evtloop;
    evtloop.SetParList(&plist);

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

    tlist.PrintStatistics();

    mh3i.DrawClone("PROFX");
    mh3e.DrawClone("PROFX");
}