source: trunk/Mars/hawc/plot_cleaned.C@ 19969

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

///////////////////////////////////////////////////////////////////////////
//
// plot_cleaned.C
// ==============
//
// This is a demonstration program showing how to do particular processing
// on a single event basis. Here we simply display uncleand und cleand
// images.
//
// Therefor MInteractiveTask is used, which gives you the possibility
// to develop new tasks without the need of compilation.
// The input is a merpp raw-data file.
//
///////////////////////////////////////////////////////////////////////////

// ==========================================================================
// ============ see plot_cleaned function at the end of the file ============
// ==========================================================================
Bool_t HandleInput()
{
    // This is a pure man's command line interface to wait for a key input
    // and allow exit but at the same time allow interaction with the GUI
    TTimer timer("gSystem->ProcessEvents();", 50, kFALSE);
    while (1)
    {
        // While reading the input process gui events asynchronously
        timer.TurnOn();
        const char *gl = Getline("Type 'q' to exit, <return> to go on: ");
        timer.TurnOff();

        TString input = gl;
        if (input=="q\n" || input==".q\n")
            return kFALSE;

        if (input=="\n")
            return kTRUE;
    };

    return kFALSE;
}

// Create four camera histograms
MHCamera display[4];

// All 'data members' that are required globally
TCanvas   *c;
MParList  *fParList;
MTaskList *fTaskList;

//
// Called like all PreProcess functions of tasks. Get the access to
// the containers necessary for you.
//
Int_t PreProcess(MParList *plist)
{
    // Get parameter and tasklist, see Process
    fParList = plist;
    fTaskList = (MTaskList*)plist->FindObject("MTaskList");

    // Get camera geoemtry
    MGeomCam *geomcam = (MGeomCam*)plist->FindObject("MGeomCam");

    // setup canvas and camera-histograms
    c = new TCanvas("Events", "Real Events", 600, 600);
    c->SetBorderMode(0);
    c->Divide(2,2);
    for (int i=0; i<4; i++)
    {
        // Tell all histograms which is the geometry to use
        display[i].SetGeometry(*geomcam);

        // Add the camera object to the pads
        c->cd(i+1);
        display[i].Draw();

        // Add the MHillas object (which pains the ellipse)
        gPad->cd(1);
        plist->FindObject("MHillas")->Draw();
    }
    return kTRUE;
}

//
// Called like all Process functions of tasks. Process a single
// event - here display it.
//
Int_t Process()
{
    // For simplicity we search in the Process function for the
    // objects. This is deprectaed! Store the pointers to the objects
    // as data member and get the pointers in PreProcess.
    MReadMarsFile *read  = (MReadMarsFile*)fTaskList->FindObject("MRead");
    MSignalCam    *cam   = (MSignalCam*)   fTaskList->FindObject("MSignalCam");
    MImgCleanStd  *clean = (MImgCleanStd*) fTaskList->FindObject("MImgCleanStd");
    MGeomCam      *geom  = (MGeomCam*)     fParList->FindObject("MGeomCam");

    // Ouput event number
    cout << "Event #" << read->GetNumEntry() << ":" << endl;

    // Fill the data into your camera-histograms
    display[0].SetCamContent(*(MSignalCam*)fParList->FindObject("MSignalCam"), 9);
    display[1].SetCamContent(*(MSignalCam*)fParList->FindObject("MSignalCam"), 3);
    display[2].SetCamContent(*(MSignalCam*)fParList->FindObject("MSignalCam"), 6);
    display[3].SetCamContent(*(MSignalCam*)fParList->FindObject("MSignalCam"), 8);

    // Index of island
    //display[3].SetCamContent(*(MSignalCam*)fParList->FindObject("MSignalCam"), 5);

    // Update the display
    for (int i=1; i<=4; i++)
    {
        c->GetPad(i)->Modified();
        c->GetPad(i)->Update();
    }

    // print the data on the console
    ((MHillas*)fParList->FindObject("MHillas"))->Print(*geom);
    ((MHillasExt*)fParList->FindObject("MHillasExt"))->Print(*geom);
    ((MNewImagePar*)fParList->FindObject("MNewImagePar"))->Print(*geom);

    cout << "\n===================================================================" << endl;
    bool rc = HandleInput();
    cout << "===================================================================" << endl;

    // wait for 'return'
    return rc;
}

//
// Called like all PostProcess functions of tasks. Delete
// instanciated objects.
//
Int_t PostProcess()
{
    delete c;
}

// ==========================================================================
//
// Run the macro with
//
//      root hawc/plot_cleaned.C\(\"00000001.003_Y_MonteCarlo003_Events.root\"\)
//
// The default file name is a Y-file (either as Monte Carlo-truth from ceres
// or from the callisto)
//
// From within root, the escape characters can be omitted, e.g.
//
//      root
//      [0] .x hawc/plot_cleaned.C("00000001.003_Y_MonteCarlo003_Events.root")
//
// ==========================================================================

void plot_cleaned(const char *fname)
{
    // Setup parameter- and tasklist
    MParList  plist;
    MTaskList tlist;
    plist.AddToList(&tlist);

    // Switch off blind pixels
    MBadPixelsCam badpixels;
    badpixels.InitSize(64);
    badpixels[63].SetUnsuitable(MBadPixelsPix::kUnsuitable);
    badpixels[62].SetUnsuitable(MBadPixelsPix::kUnsuitable);
    badpixels[61].SetUnsuitable(MBadPixelsPix::kUnsuitable);
    plist.AddToList(&badpixels);

    // setup reading task
    MReadMarsFile read("Events", fname);
    read.DisableAutoScheme();

    // setup a task making sure that all arrays are resized correctly
    MGeomApply geomapl;

    // Setup a print task calling TObject::Print
    MPrint print1("MMcEvt");
    MPrint print2("MRawEvtHeader");
    // Skip them if containers are not found
    print1.EnableSkip();
    print2.EnableSkip();

    // Switch off blind pixels
    MBadPixelsTreat treat;
    treat.SetProcessPedestalRun(kFALSE);
    treat.SetProcessPedestalEvt(kFALSE);

    // Setup image cleaning
    MImgCleanTime clean;
    clean.SetMinCount(0);
    clean.SetMinSize(25);
    clean.SetDeltaT(2.5);

    // Setup calculation of Image parameters
    MHillasCalc hcalc;

    // Setup intercative task calling the functions defined above
    MTaskInteractive mytask;
    mytask.SetPreProcess(PreProcess);
    mytask.SetProcess(Process);

    // Setup your tasklist
    tlist.AddToList(&read);
    tlist.AddToList(&geomapl);
    tlist.AddToList(&print1);
    tlist.AddToList(&print2);
    tlist.AddToList(&treat);
    tlist.AddToList(&clean);
    tlist.AddToList(&hcalc);
    tlist.AddToList(&mytask);

    // Run your analysis
    MEvtLoop evtloop;
    evtloop.SetParList(&plist);

    if (!evtloop.Eventloop())
        return;

    // Print statistics information about your loop
    tlist.PrintStatistics();
}