source: trunk/MagicSoft/Mars/mjobs/MJPedestal.cc@ 3843

/* ======================================================================== *\
!
! *
! * 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, 1/2004 <mailto:tbretz@astro.uni-wuerzburg.de>
!              Markus Gaug ,04/2004 <mailto:markus@ifae.es>
!
!   Copyright: MAGIC Software Development, 2000-2004
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//
//  MJPedestal
//
/////////////////////////////////////////////////////////////////////////////
#include "MJPedestal.h"

#include <TF1.h>
#include <TFile.h>
#include <TStyle.h>
#include <TString.h>
#include <TCanvas.h>
#include <TSystem.h>

#include "MLog.h"
#include "MLogManip.h"

#include "MRunIter.h"
#include "MParList.h"
#include "MTaskList.h"
#include "MEvtLoop.h"

#include "MStatusDisplay.h"

#include "MGeomCam.h"
#include "MHCamera.h"
#include "MPedestalCam.h"
#include "MBadPixelsCam.h"

#include "MReadMarsFile.h"
#include "MRawFileRead.h"
#include "MGeomApply.h"
#include "MBadPixelsMerge.h"
#include "MPedCalcPedRun.h"

ClassImp(MJPedestal);

using namespace std;

MJPedestal::MJPedestal(const char *name, const char *title) 
    : fRuns(0), fDataCheck(kFALSE)

{
    fName  = name  ? name  : "MJPedestal";
    fTitle = title ? title : "Tool to create a pedestal file (MPedestalCam)";
}

TString MJPedestal::GetOutputFile() const
{
    if (!fRuns)
        return "";

    return Form("%s/%s-F0.root", (const char*)fOutputPath, (const char*)fRuns->GetRunsAsFileName());
}

Bool_t MJPedestal::ReadPedestalCam()
{
    const TString fname = GetOutputFile();

    if (gSystem->AccessPathName(fname, kFileExists))
    {
        *fLog << err << "Input file " << fname << " doesn't exist." << endl;
        return kFALSE;
    }

    *fLog << inf << "Reading from file: " << fname << endl;

    TFile file(fname, "READ");
    if (fPedestalCam.Read()<=0)
    {
        *fLog << err << "Unable to read MPedestalCam from " << fname << endl;
        return kFALSE;
    }

    if (file.FindKey("MBadPixelsCam"))
    {
        MBadPixelsCam bad;
        if (bad.Read()<=0)
        {
            *fLog << err << "Unable to read MBadPixelsCam from " << fname << endl;
            return kFALSE;
        }
        fBadPixels.Merge(bad);
    }

    if (fDisplay && !fDisplay->GetCanvas("Pedestals"))
        fDisplay->Read();

    return kTRUE;
}

void MJPedestal::DrawProjection(MHCamera *obj1, Int_t fit) const
{
    TH1D *obj2 = (TH1D*)obj1->Projection(obj1->GetName());
    obj2->SetDirectory(0);
    obj2->Draw();
    obj2->SetBit(kCanDelete);
    gPad->SetLogy();

    if (obj1->GetGeomCam().InheritsFrom("MGeomCamMagic"))
    {
        TArrayI s0(3);
        s0[0] = 6;
        s0[1] = 1;
        s0[2] = 2;

        TArrayI s1(3);
        s1[0] = 3;
        s1[1] = 4;
        s1[2] = 5;

        TArrayI inner(1);
        inner[0] = 0;

        TArrayI outer(1);
        outer[0] = 1;

        // Just to get the right (maximum) binning
        TH1D *half[4];
        half[0] = obj1->ProjectionS(s0, inner, "Sector 6-1-2 Inner");
        half[1] = obj1->ProjectionS(s1, inner, "Sector 3-4-5 Inner");
        half[2] = obj1->ProjectionS(s0, outer, "Sector 6-1-2 Outer");
        half[3] = obj1->ProjectionS(s1, outer, "Sector 3-4-5 Outer");

        for (int i=0; i<4; i++)
        {
            half[i]->SetLineColor(kRed+i);
            half[i]->SetDirectory(0);
            half[i]->SetBit(kCanDelete);
            half[i]->Draw("same");
        }
    }

    const Double_t min   = obj2->GetBinCenter(obj2->GetXaxis()->GetFirst());
    const Double_t max   = obj2->GetBinCenter(obj2->GetXaxis()->GetLast());
    const Double_t integ = obj2->Integral("width")/2.5066283;
    const Double_t mean  = obj2->GetMean();
    const Double_t rms   = obj2->GetRMS();
    const Double_t width = max-min;

    if (rms==0 || width==0)
        return;

    const TString dgausformula("([0]-[3])/[2]*exp(-0.5*(x-[1])*(x-[1])/[2]/[2])"
                               "+[3]/[5]*exp(-0.5*(x-[4])*(x-[4])/[5]/[5])");

    TF1 *f=0;
    switch (fit)
    {
        // FIXME: MAYBE add function to TH1?
    case 0:
        f = new TF1("sgaus", "gaus(0)", min, max);
        f->SetLineColor(kYellow);
        f->SetBit(kCanDelete);
        f->SetParNames("Area", "#mu", "#sigma");
        f->SetParameters(integ/rms, mean, rms);
        f->SetParLimits(0, 0,   integ);
        f->SetParLimits(1, min, max);
        f->SetParLimits(2, 0,   width/1.5);
        obj2->Fit(f, "QLRM");
        break;

    case 1:
        f = new TF1("dgaus", dgausformula, min, max);
        f->SetLineColor(kYellow);
        f->SetBit(kCanDelete);
        f->SetParNames("A_{tot}", "#mu_{1}", "#sigma_{1}", "A_{2}", "#mu_{2}", "#sigma_{2}");
        f->SetParameters(integ,         (min+mean)/2, width/4,
                         integ/width/2, (max+mean)/2, width/4);
        // The left-sided Gauss
        f->SetParLimits(0, integ-1.5,      integ+1.5);
        f->SetParLimits(1, min+(width/10), mean);
        f->SetParLimits(2, 0,              width/2);
        // The right-sided Gauss
        f->SetParLimits(3, 0,    integ);
        f->SetParLimits(4, mean, max-(width/10));
        f->SetParLimits(5, 0,    width/2);
        obj2->Fit(f, "QLRM");
        break;

    default:
        obj2->Fit("gaus", "Q");
        obj2->GetFunction("gaus")->SetLineColor(kYellow);
        break;
    }
}

void MJPedestal::CamDraw(TCanvas &c, Int_t x, Int_t y, MHCamera &cam1, Int_t fit)
{
    c.cd(x);
    gPad->SetBorderMode(0);
    MHCamera *obj1=(MHCamera*)cam1.DrawCopy("hist");

    obj1->AddNotify(&fPedestalCam);

    c.cd(x+y);
    gPad->SetBorderMode(0);
    obj1->Draw();

    c.cd(x+2*y);
    gPad->SetBorderMode(0);
    DrawProjection(obj1, fit);
}

void MJPedestal::DisplayResult(MParList &plist)
{

    if (!fDisplay)
        return;

    //
    // Update display
    //
    TString title = fDisplay->GetTitle();
    title += "--  Pedestal ";
    title += fRuns->GetRunsAsString();
    title += "  --";
    fDisplay->SetTitle(title);

    //
    // Get container from list
    //
    MGeomCam &geomcam = *(MGeomCam*)plist.FindObject("MGeomCam");

    //
    // Create container to display
    //
    MHCamera disp0(geomcam, "MPedestalCam;ped", "Pedestals");
    MHCamera disp1(geomcam, "MPedestalCam;var", "Sigma Pedestal");

    disp0.SetCamContent(fPedestalCam, 0);
    disp0.SetCamError  (fPedestalCam, 1);

    disp1.SetCamContent(fPedestalCam, 2);
    disp1.SetCamError  (fPedestalCam, 3);

    disp0.SetYTitle("P [fadc/slice]");
    disp1.SetYTitle("\\sigma_{P} [fadc/slice]");

    //
    // Display data
    //
    TCanvas &c3 = fDisplay->AddTab("Pedestals");
    c3.Divide(2,3);

    CamDraw(c3, 1, 2, disp0, 0);
    CamDraw(c3, 2, 2, disp1, 1);
}

Bool_t MJPedestal::WriteResult()
{
    if (fOutputPath.IsNull())
        return kTRUE;

    const TString oname(GetOutputFile());

    *fLog << inf << "Writing to file: " << oname << endl;

    TFile file(oname, "UPDATE");

    if (fDisplay && fDisplay->Write()<=0)
    {
        *fLog << err << "Unable to write MStatusDisplay to " << oname << endl;
        return kFALSE;
    }

    if (fPedestalCam.Write()<=0)
    {
        *fLog << err << "Unable to write MPedestalCam to " << oname << endl;
        return kFALSE;
    }

    if (fBadPixels.Write()<=0)
    {
        *fLog << err << "Unable to write MBadPixelsCam to " << oname << endl;
        return kFALSE;
    }

    return kTRUE;
}

void MJPedestal::SetOutputPath(const char *path)
{
    fOutputPath = path;
    if (fOutputPath.EndsWith("/"))
        fOutputPath = fOutputPath(0, fOutputPath.Length()-1);
}

Bool_t MJPedestal::Process()
{
    if (!ReadPedestalCam())
        return ProcessFile();

    return kTRUE;
}

Bool_t MJPedestal::ProcessFile()
{
    if (!fRuns)
    {
        *fLog << err << "No Runs choosen... abort." << endl;
        return kFALSE;
    }
    if (fRuns->GetNumRuns() != fRuns->GetNumEntries())
    {
        *fLog << err << "Number of files found doesn't match number of runs... abort." << endl;
        return kFALSE;
    }

    *fLog << inf;
    fLog->Separator(GetDescriptor());
    *fLog << "Calculate MPedestalCam from Runs " << fRuns->GetRunsAsString() << endl;
    *fLog << endl;

    MParList plist;
    MTaskList tlist;
    plist.AddToList(&tlist);

    MReadMarsFile read("Events");
    MRawFileRead rawread("");

    if (fDataCheck)
    {
//      rawread.AddFiles(*fRuns);
        tlist.AddToList(&rawread);
    }
    else
    {
        read.DisableAutoScheme();
        static_cast<MRead&>(read).AddFiles(*fRuns);
        tlist.AddToList(&read);
    }
    // Enable logging to file
    //*fLog.SetOutputFile(lname, kTRUE);

    // Setup Tasklist
    plist.AddToList(&fPedestalCam);
    plist.AddToList(&fBadPixels);

    MGeomApply      geomapl;
    // MBadPixelsMerge merge(&fBadPixels);
    //MExtractSignal sigcalc;
    MPedCalcPedRun  pedcalc;

    tlist.AddToList(&geomapl);
    //    tlist.AddToList(&merge);
    //tlist.AddToList(&sigcalc);
    tlist.AddToList(&pedcalc);

    //
    // Execute the eventloop
    //
    MEvtLoop evtloop(fName);
    evtloop.SetParList(&plist);
    evtloop.SetDisplay(fDisplay);
    evtloop.SetLogStream(fLog);

    // Execute first analysis
    if (!evtloop.Eventloop())
    {
        *fLog << err << GetDescriptor() << ": Failed." << endl;
        return kFALSE;
    }

    tlist.PrintStatistics();

    DisplayResult(plist);

    if (!WriteResult())
        return kFALSE;

    *fLog << inf << GetDescriptor() << ": Done." << endl;

    return kTRUE;
}