source: trunk/MagicSoft/Mars/mjobs/MJExtractCalibTest.cc@ 4068

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

/////////////////////////////////////////////////////////////////////////////
//
//  MJExtractCalibTest
//
/////////////////////////////////////////////////////////////////////////////
#include "MJExtractCalibTest.h"

#include <TFile.h>
#include <TStyle.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 "MHCamera.h"

#include "MPedestalCam.h"
#include "MBadPixelsCam.h"
#include "MCerPhotEvt.h"
#include "MArrivalTime.h"
#include "MCalibrationChargeCam.h"
#include "MCalibrationRelTimeCam.h"
#include "MCalibrationQECam.h"
#include "MHCamEvent.h"

#include "MReadMarsFile.h"
#include "MGeomApply.h"
#include "MExtractSlidingWindow.h"
#include "MExtractor.h"
#include "MExtractTime.h"
#include "MExtractTimeFastSpline.h"
#include "MFCosmics.h"
#include "MContinue.h"
#include "MFillH.h"
#include "MCalibrate.h"
#include "MCalibrateRelTimes.h"
#include "MPedPhotCalc.h"

#include "MStatusDisplay.h"

ClassImp(MJExtractCalibTest);

using namespace std;
// --------------------------------------------------------------------------
//
// Default constructor. 
//
// Sets fRuns to 0, fExtractor to NULL, fTimeExtractor to NULL
//
MJExtractCalibTest::MJExtractCalibTest(const char *name, const char *title) 
    : fRuns(NULL), fExtractor(NULL), fTimeExtractor(NULL)
{
    fName  = name  ? name  : "MJExtractCalibTest";
    fTitle = title ? title : "Tool to extract, calibrate and test signals from a file";
}


void MJExtractCalibTest::DisplayResult(MParList &plist)
{
  if (!fDisplay)
    return;
  
  //
  // Update display
  //
  TString title = fDisplay->GetTitle();
  title += "--  Extraction-Calibration-Test ";
  title += fRuns->GetRunsAsString();
  title += "  --";
  fDisplay->SetTitle(title);

  //
  // Get container from list
  //
  MGeomCam &geomcam = *(MGeomCam*)plist.FindObject("MGeomCam");
  
  // Create histograms to display
  MHCamera disp1 (geomcam, "Test;Photons",           "Mean of calibrated Photons");
  MHCamera disp2 (geomcam, "Test;SigmaPhotons",      "Sigma of calibrated photons");
  MHCamera disp3 (geomcam, "Test;PhotonsPerArea",    "Equiv. Cherenkov Photons per Area");
  MHCamera disp4 (geomcam, "Test;SigmaPhotPerArea",  "Sigma equiv. Cher. Photons per Area");

  // Fitted charge means and sigmas
  disp1.SetCamContent(fTestCam,  0);
  disp1.SetCamError(  fTestCam,  1);
  disp2.SetCamContent(fTestCam,  2);
  disp2.SetCamError(  fTestCam,  3);
  disp3.SetCamContent(fTestCam,  7);
  disp3.SetCamError(  fTestCam,  8);
  disp4.SetCamContent(fTestCam,  9);
  disp4.SetCamError(  fTestCam,  10);

  disp1.SetYTitle("Photons");
  disp2.SetYTitle("\\sigma_{phot}");
  disp3.SetYTitle("Photons per Area [mm^{-2}]");
  disp4.SetYTitle("\\sigma_{phot} per Area [mm^{-2}]");
  
  gStyle->SetOptStat(1111);
  gStyle->SetOptFit();

  TCanvas &c = fDisplay->AddTab("TestCharges");
  c.Divide(4,4);

  CamDraw(c, 1, 4, disp1, 2, 1);
  CamDraw(c, 2, 4, disp2, 2, 1);        
  CamDraw(c, 3, 4, disp3, 1, 1);        
  CamDraw(c, 4, 4, disp4, 2, 1);        

  return;

}


void MJExtractCalibTest::DisplayResultT(MParList &plist)
{
  if (!fDisplay)
    return;
  
  //
  // Update display
  //
  TString title = fDisplay->GetTitle();
  title += "--  Extraction-Calibration-Test-Time";
  title += fRuns->GetRunsAsString();
  title += "  --";
  fDisplay->SetTitle(title);

  //
  // Get container from list
  //
  MGeomCam &geomcam = *(MGeomCam*)plist.FindObject("MGeomCam");
  
  // Create histograms to display
  MHCamera disp1 (geomcam, "Test;Arr.Times",           "Mean of calibrated Arr.Times");
  MHCamera disp2 (geomcam, "Test;SigmaArr.Times",      "Sigma of calibrated Arr.Times");

  // Fitted charge means and sigmas
  disp1.SetCamContent(fTestTimeCam,  0);
  disp1.SetCamError(  fTestTimeCam,  1);
  disp2.SetCamContent(fTestTimeCam,  2);
  disp2.SetCamError(  fTestTimeCam,  3);

  disp1.SetYTitle("Mean Arr.Times [FADC units]");
  disp2.SetYTitle("\\sigma_{t} [FADC units]");
  
  gStyle->SetOptStat(1111);
  gStyle->SetOptFit();

  TCanvas &c = fDisplay->AddTab("TestTimes");
  c.Divide(4,2);

  CamDraw(c, 1, 2, disp1, 5, 1);
  CamDraw(c, 2, 2, disp2, 5, 1);        

  return;

}


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

TString MJExtractCalibTest::GetOutputFile() const
{
  if (!fRuns)
    return "";
  
  return Form("%s/%s-Test.root", (const char*)fOutputPath, (const char*)fRuns->GetRunsAsFileName());
}


Bool_t MJExtractCalibTest::ProcessD(MPedestalCam &pedcam, MCalibrationChargeCam &calcam, MCalibrationQECam &qecam)
{
    // const TString fname = GetOutputFile();
  
//  if (gSystem->AccessPathName(fname, kFileExists))
    return ProcessFileD(pedcam,calcam,qecam);
  
    // return kTRUE;
}

Bool_t MJExtractCalibTest::ProcessFileD(MPedestalCam &pedcam, MCalibrationChargeCam &calcam, MCalibrationQECam &qecam)
{
  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 MExtractedSignalCam from Runs " << fRuns->GetRunsAsString() << endl;
  *fLog << endl;
  
  MCerPhotEvt    cerphot;
  MPedPhotCam    pedphot;
  // Setup Lists
  MParList plist;
  plist.AddToList(&pedcam);
  plist.AddToList(&calcam);
  plist.AddToList(&qecam);
  plist.AddToList(&cerphot);
  plist.AddToList(&pedphot);
  plist.AddToList(&fTestCam);
  plist.AddToList(&fBadPixels);

  MTaskList tlist;
  plist.AddToList(&tlist);

  // Setup Task-lists
  MReadMarsFile read("Events");
  read.DisableAutoScheme();
  static_cast<MRead&>(read).AddFiles(*fRuns);
  
  MGeomApply            apply; // Only necessary to craete geometry
  MExtractSlidingWindow extract2;
  MCalibrate            photcalc;
  photcalc.SetCalibrationMode(MCalibrate::kFfactor);
  MPedPhotCalc    pedphotcalc;  

  MHCamEvent evt("ExtSignal");
  evt.SetType(0);
  MFillH fill(&evt, "MExtractedSignalCam");
  
  MFillH fillcam("MHCalibrationTestCam", "MCerPhotEvt");
  fillcam.SetNameTab("Test");

  MFCosmics cosmics;
  MContinue cont(&cosmics);
  
  tlist.AddToList(&read);
  tlist.AddToList(&apply);

  if (fExtractor)
    tlist.AddToList(fExtractor);
  else
  {
      *fLog << warn << GetDescriptor() 
            << ": No extractor has been chosen, take default MExtractSlidingWindow " << endl;
      tlist.AddToList(&extract2);
  }

  //  if (TestBit(kEnableGraphicalOutput))
  tlist.AddToList(&cont);
  tlist.AddToList(&fill);
  tlist.AddToList(&photcalc);
  tlist.AddToList(&pedphotcalc);
  tlist.AddToList(&fillcam);
  
  // Create and setup 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 (!WriteResultD())
    return kFALSE;

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

Bool_t MJExtractCalibTest::ProcessT(MPedestalCam &pedcam, MCalibrationRelTimeCam &relcam)
{

//  const TString fname = GetOutputFile();
  
//  if (gSystem->AccessPathName(fname, kFileExists))
    return ProcessFileT(pedcam,relcam);
  
//  return kTRUE;
}

Bool_t MJExtractCalibTest::ProcessFileT(MPedestalCam &pedcam, MCalibrationRelTimeCam &relcam)
{

  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 MExtractedSignalCam from Runs " << fRuns->GetRunsAsString() << endl;
  *fLog << endl;
  
  MArrivalTime     arrtime;

  // Setup Lists
  MParList plist;
  plist.AddToList(&pedcam);
  plist.AddToList(&relcam);
  plist.AddToList(&arrtime);
  plist.AddToList(&fTestTimeCam);
  plist.AddToList(&fBadPixels);
  
  MTaskList tlist;
  plist.AddToList(&tlist);

  // Setup Task-lists
  MReadMarsFile read("Events");
  read.DisableAutoScheme();
  static_cast<MRead&>(read).AddFiles(*fRuns);
  
  MGeomApply             apply; // Only necessary to craete geometry
  MExtractTimeFastSpline extract;
  MExtractSlidingWindow  extcharge; // Only for the cosmics filter
  MCalibrateRelTimes     timecalc;
  MFCosmics              cosmics;
  MContinue              cont(&cosmics);

  MHCamEvent evt("ExtTimes");
  evt.SetType(0);
  MFillH fill(&evt, "MArrivalTimeCam");
  
  MFillH fillcam("MHCalibrationTestTimeCam", "MArrivalTime");
  fillcam.SetNameTab("TestTime");

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

  if (fTimeExtractor)
    tlist.AddToList(fTimeExtractor);
  else
  {
      *fLog << warn << GetDescriptor() 
            << ": No extractor has been chosen, take default MExtractTimeFastSpline " << endl;
      tlist.AddToList(&extract);
  }

  tlist.AddToList(&extcharge);
  tlist.AddToList(&cont);
  tlist.AddToList(&fill);
  tlist.AddToList(&timecalc);
  tlist.AddToList(&fillcam);
  
  // Create and setup 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();
  
  DisplayResultT(plist);

  if (!WriteResultT())
    return kFALSE;

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


Bool_t MJExtractCalibTest::ReadPedPhotCam()
{
    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 (fPedPhotCam.Read()<=0)
    {
        *fLog << "Unable to read MPedPhotCam from " << fname << endl;
        return kFALSE;
    }

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

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

    return kTRUE;
}

Bool_t MJExtractCalibTest::WriteResultD()
{
    
    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 (fPedPhotCam.Write()<=0)
    {
        *fLog << err << "Unable to write MPedPhotCam to " << oname << endl;
        return kFALSE;
    }

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

    return kTRUE;

}

Bool_t MJExtractCalibTest::WriteResultT()
{
    
    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 (fTestTimeCam.Write()<=0)
    {
        *fLog << err << "Unable to write MHCalibrationTestTimeCam to " << oname << endl;
        return kFALSE;
    }

    return kTRUE;

}