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

#include <TStyle.h>
#include <TCanvas.h>
#include <TSystem.h>
#include <TF1.h>
#include <TProfile.h>
#include <TProfile2D.h>
#include <TMath.h>
#include <TGraph.h>
#include <TLine.h>
#include <TFitResultPtr.h>
#include <TFitResult.h>
#include <TFile.h>
#include <TLine.h>

#include <cstdio>
#include <stdio.h>
#include <stdint.h>

#include "Getline.h"
#include "MH.h"
#include "MArrayI.h"
#include "MDirIter.h"
#include "MFillH.h"
#include "MEvtLoop.h"
#include "MCamEvent.h"
#include "MGeomApply.h"
#include "MTaskList.h"
#include "MParList.h"
#include "MContinue.h"
#include "MBinning.h"
#include "MDrsCalibApply.h"
#include "MDrsCalibration.h"
#include "MRawFitsRead.h"
#include "MStatusDisplay.h"
#include "MTaskInteractive.h"
#include "MPedestalSubtractedEvt.h"
#include "MGeomCamFAMOUS.h"
#include "MRawRunHeader.h"
#include "MPedestalCam.h"
#include "MPedestalPix.h"
#include "MParameters.h"

// ==========================================================================
// ============ see plot_callisto function at the end of the file ===========
// ==========================================================================
int 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 reading the input process gui events asynchronously
    timer.TurnOn();
    const char *gl = Getline("Type 'q' to exit, or event number and <return> to go on: ");
    timer.TurnOff();

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

    return kTRUE;
}

// Create a histogram for the trace
//TH1D hist("Trace", "Waveform", 1024, -0.25, 511.75); //range = [-0.25,511.75]ns
TH1D hist("Trace", "Waveform (Pixel 5, Event 1)", 1024, -0.5, 1024-0.5);

// All 'data members' that are required globally
TCanvas *c;
MPedestalSubtractedEvt *fEvt = 0;

//
// Called like all PreProcess functions of tasks. Get the access to
// the containers necessary for you.
//
Int_t PreProcess(MParList *plist)
{
    fEvt = (MPedestalSubtractedEvt*)plist->FindObject("MPedestalSubtractedEvt");

    //if (!fEvt)
    //{
    //    gfLog << err << "MPedestalSubtractedEvt not found... aborting." << endl;
    //    return kFALSE;
    //}

    // Create a canvas and plot the histogram into it
    c = new TCanvas();
    hist.SetStats(kFALSE);
    hist.Draw();
    hist.SetXTitle("Samples");
    hist.SetYTitle("Amplitude [a.u.]");

    //TF1 constant("zero", "[0]", -0.25, 511.75);
    TF1 constant("zero", "[0]", -0.5, 1024);
    constant.SetParameter(0, 0);
    constant.SetLineColor(kBlack);
    constant.SetLineStyle(kDashed);
    constant.DrawCopy("same");

    // One p.e.
    constant.SetParameter(0, 20);
    constant.SetLineStyle(kDotted);
    constant.DrawCopy("same");

    constant.SetParameter(0, -20);
    constant.DrawCopy("same");


    return kTRUE;
}

Int_t Process()
{
    int pixel_index = 6; //plot trace only of this pixel
    int skipstart = 5;
    int skipend = 10;

    // Number of samples per pixel (usually 1024)
    int numsamples = fEvt->GetNumSamples();
    Float_t *cal = fEvt->GetSamples(pixel_index);

    // reset contents of histogram and fill new contents
    hist.Reset();
    for (int i=0; i<numsamples; i++)
        hist.Fill(i, cal[i]); // 0.5ns = 1 sample

    // Reset min/max range
    hist.SetMinimum(-50);
    hist.SetMaximum(50);

    //// Set a reasonable range
    //if (hist.GetMinimum()>-30)
    //    hist.SetMinimum(-20);
    //if (hist.GetMaximum()<60)
    //    hist.SetMaximum(50);



    TLine *line = new TLine(skipstart, 50, skipstart, -50);
    line->SetLineStyle(kSolid);
    line->SetLineColor(kRed);
    line->Draw();

    TLine *line2 = new TLine(1024-skipend, 50, 1024-skipend, -50);
    line2->SetLineStyle(kSolid);
    line2->SetLineColor(kRed);
    line2->Draw();

    // Signal root to update the canvas/pad
    c->Modified();
    c->Update();
    //c->SaveAs("/home/giangdo/Documents/Master/MA/pulse/data_single_pe/traces/pixel3/event_0.png");


    bool rc = HandleInput();

    // wait for 'return'
    return rc;
}


/*****************************************************************

 This is an example how to access and display calibrated raw data

 datafile:
    A data file written by the fadctrl or ceres, e.g.
    20170727_006.fits.fz or
    00000015.003_D_MonteCarlo019_Events.fits.fz

 drsfile:
    Usually the third of the three .drs.fits files from a
    DRS calibration sequence, e.g. 20170727_004.drs.fits

 To run the macro from the command line (assuming you are in a directory
 Mars/build where you have built your Mars environment) ou can do

    root ../hawc/plot_trace.C\(\"20191002_018.fits.fz\",\"20191002_013.drs.fits\"\)

 or from within root

    [0] .x ../hawc/plot_trace.C("20191002_018.fits.fz", "20191002_013.drs.fits")

******************************************************************/
int plot_trace(const char *datafile, const char *drsfile)
{
    // ======================================================

    // true:  Display correctly mapped pixels in the camera displays
    //        but the value-vs-index plot is in software/spiral indices
    // false: Display pixels in hardware/linear indices,
    //        but the order is the camera display is distorted.
    //        false is assumed automatically for files with the ISMC flag set.
    bool usemap = true;

    // mapping file (found in Mars/hawc)
    const char *mmap = usemap ? "/home/giangdo/Documents/Master/MA/Software/Mars/hawc/HAWCsEyemap181214.txt" : NULL;

    // ======================================================

    // Check if the requested mapping file is available
    if (mmap && gSystem->AccessPathName(mmap, kFileExists))
    {
        gLog << err << "ERROR - Cannot access mapping file '" << mmap << "'" << endl;
        return 11;
    }

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

    // Check and read the DRS calibration constants
    MDrsCalibration drscalib300;
    if (drsfile && !drscalib300.ReadFits(drsfile))
        return 31;

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

    // Setup camera geometry
    MGeomCamFAMOUS geom;

    // ======================================================

    // Setup task for reading the file
    MRawFitsRead read;
    read.LoadMap(mmap);
    read.AddFile(datafile);

    // Setup task which ensures sizes of containers consistent with camera geometry
    MGeomApply apply;

    // Setup task for DRS calibration
    MDrsCalibApply drsapply;
    drsapply.SetMaxNumPrevEvents(0); // Switched off step-correction -> crashes due to the requirement of N*9 ch

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

    // ======================================================

    MTaskList tlist;

    MParList plist;
    plist.AddToList(&tlist);
    plist.AddToList(&drscalib300);
    plist.AddToList(&geom);

    // ------------------ Setup eventloop and run analysis ---------------

    tlist.AddToList(&read);
    tlist.AddToList(&apply);
    tlist.AddToList(&drsapply);
    tlist.AddToList(&mytask);

    // ------------------ Setup and run eventloop ---------------

    MEvtLoop loop(gSystem->BaseName(datafile));
    loop.SetParList(&plist);

    if (!loop.Eventloop(100))
        return 4;

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

    return 0;
}