source: trunk/MagicSoft/Mars/macros/calibration.C@ 4841

/* ======================================================================== *\
!
! *
! * 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, 11/2003 <mailto:markus@ifae.es>
!
!   Copyright: MAGIC Software Development, 2000-2004
!
!
\* ======================================================================== */
/////////////////////////////////////////////////////////////////////////////
//
//  calibration.C
//
//  Needs as arguments the run number of a calibration file ("*_C_*.root") and 
//  the run number of the corresponding pedestal file ("*_P_*.root"). 
//
//  The TString inpath has to be set correctly.
//
//  The macro searches for the pulser colour which corresponds to the calibration
//  run number. If the run number is smaller than 20000, pulser colour "CT1" 
//  is assumed, otherwise, it searches for the strings "green", "blue", "uv" or 
//  "ct1" in the filenames. If no colour or multiple colours are found, the 
//  execution is aborted.  
//
//  The container MBadPixelsCam is created and followed during the execution of the 
//  rest of the macro.
// 
//  A first loop over the pedestal file is performed using the class MJPedestal
//
//  The container MCalibrationQECam is created and followed during the execution of the 
//  rest of the macro.
//
//  A loop over the calibration files is performed using the class MJCalibration. 
//  The results are displayed using the MJCalibration::SetNormalDisplay() facility, 
//  but other displays can easily be uncommented. 
//  The call to MJCalibration skips the relative time calibration, which can be 
//  uncommented as well. 
// 
//  Last, a third loop is performed over the calibration file again in order to 
//  "calibrate" it and test the resulting outcome.
//
/////////////////////////////////////////////////////////////////////////////
#include "MJPedestal.h"
#include "MJCalibration.h"
#include "MJExtractSignal.h"
#include "MJExtractCalibTest.h"
#include "MExtractFixedWindowPeakSearch.h"
#include "MExtractSlidingWindow.h"
#include "MExtractFixedWindow.h"
#include "MExtractFixedWindowSpline.h"
#include "MExtractAmplitudeSpline.h"
#include "MExtractTimeHighestIntegral.h"
#include "MExtractTimeFastSpline.h"
#include "MRunIter.h"
#include "MStatusDisplay.h"
#include "MCalibrationQECam.h"
#include "MHCalibrationTestCam.h"
#include "MHCalibrationTestPix.h"
#include "MBadPixelsCam.h"
#include "MArgs.h"
#include "MArray.h"
#include "MLog.h"
#include "MParContainer.h"

#include "TStyle.h"
#include "TObject.h"
#include "TObjectTable.h"
#include "TSystem.h"

#include <fstream>

using namespace std;

static TString outpath = "./";
static TString inpath  = "/home/rootdata/Calib/2004_07_06/";
static TString badfile = "";
//
// the default pedestal run for the calibration
//
static const Int_t   pedrun  = 31741;
//
// the default start calibration run 
//
static const Int_t   calrun1 = 31742;
//
// the default last calibration run (if 0, only one run is taken, otherwise consecutive runs 
// between calrun1 and calrun2)
//
static  const Int_t calrun2 = 0;
//
// A switch to output debugging information about Objects use
//
static Bool_t debug = kFALSE;
//
// A switch to use the Blind Pixel
//
static Bool_t blindpix = kTRUE;
// A switch to use the PIN Diode 
//
static Bool_t pindiode = kFALSE;
//
// Tell if you want to calibrate times:
//
static Bool_t useTimes = kFALSE;
//
// Tell if you want to use the display:
//
static Bool_t useDisplay = kTRUE;
//
// Tell if you want to test the result afterwards
//
static Bool_t useTest = kTRUE;
//
Int_t calibration(const Int_t prun=pedrun, 
                  const Int_t crun1=calrun1, const Int_t crun2=calrun2)
{

  if (debug)
    TObject::SetObjectStat(kTRUE);

  //
  // Choose the signal Extractor:
  //
  //  MExtractFixedWindowPeakSearch extractor;
  //  MExtractSlidingWindow  extractor;   
  //  MExtractFixedWindowSpline extractor;
  MExtractFixedWindow extractor;
  //  MExtractAmplitudeSpline extractor;
  //  MExtractTimeAndChargeSpline extractor;
  //
  // Set Ranges or Windows
  //
  extractor.SetRange(3,14,4,13);
  //  extractor.SetWindows(8,8);

  //
  // Choose the arrival time Extractor:
  //
  //  MExtractTimeHighestIntegral timeext;
  MExtractTimeFastSpline timeext;
  //
  // Set Ranges or Windows
  //
  timeext.SetRange(0,7,3,8);

  MRunIter pruns;
  MRunIter cruns;

  pruns.AddRun(prun,inpath);

  if (crun2==0)
    cruns.AddRun(crun1,inpath);
  else
    cruns.AddRuns(crun1,crun2,inpath);

  gStyle->SetOptStat(1111);
  gStyle->SetOptFit();
  gStyle->SetTitleSize(0.1,"u");
  gStyle->SetLineWidth(1);

  MStatusDisplay *display = NULL;

  if (useDisplay)
    {
      display = new MStatusDisplay;
      display->SetUpdateTime(3000);
      display->Resize(850,700);
    }
  
  /************************************/
  /* FIRST LOOP: PEDESTAL COMPUTATION */
  /************************************/

  MCalibrationQECam qecam;
  MBadPixelsCam     badcam;
  //
  // If you want to exclude pixels from the beginning, read 
  // an ascii-file with the corr. pixel numbers (see MBadPixelsCam)
  //
  if (!badfile.IsNull())
    {
      ifstream f(badfile.Data());
      badcam.AsciiRead(f);  
      f.close();
    }

  MJPedestal pedloop;
  pedloop.SetExtractor(&extractor);
  pedloop.SetInput(&pruns);
  pedloop.SetPathOut(outpath.Data());
  if (useDisplay)
    {
      pedloop.SetDisplay(display);
      pedloop.SetDataCheckDisplay();
    }
  pedloop.SetBadPixels(badcam);

  if (!pedloop.Process())
    return 1;

  /****************************************/
  /* SECOND LOOP: CALIBRATION COMPUTATION */
  /****************************************/

  MJCalibration calloop;

  if (debug)
    calloop.SetDebug();
  //
  // If you want to run the data-check on RAW DATA!!!, choose:
  //  calloop.SetDataCheck();
  // 
  // If you want to see the data-check plots only, choose:
  calloop.SetDataCheckDisplay();
  // 
  // For everything, you have ever dreamed of, choose:
  // calloop.SetFullDisplay();

  //
  // If you want to calibrate the times as well, choose:
  //
  calloop.SetRelTimeCalibration(useTimes);
  calloop.SetExtractor(&extractor);
  calloop.SetTimeExtractor(&timeext);
  calloop.SetInput(&cruns);
  calloop.SetPathOut(outpath.Data());
  if (useDisplay)
    calloop.SetDisplay(display);
  calloop.SetUseBlindPixel(blindpix);
  calloop.SetUsePINDiode(pindiode);
  calloop.SetQECam(qecam);
  calloop.SetBadPixels(pedloop.GetBadPixels());

  if (!calloop.Process(pedloop.GetPedestalCam()))
    return 2;

  //
  // The next lines are the use the Print() function and have 
  // all the results as ascii-tables:
  //
  if (debug)
    {
      MCalibrationChargeCam &chargecam   = calloop.GetCalibrationCam();
      MCalibrationQECam     &nqecam      = calloop.GetQECam();
      MBadPixelsCam         &badbad      = calloop.GetBadPixels();
      chargecam.Print();
      nqecam.Print();
      badbad.Print();
    }

  gLog << endl;
  gLog << "Mean number of photons from pulser Inner pixels (F-Factor Method): " 
       << calloop.GetCalibrationCam().GetNumPhotonsFFactorMethod() 
       << " +- " << calloop.GetCalibrationCam().GetNumPhotonsFFactorMethodErr() << endl;
  gLog << endl;

  /********************************************************************/
  /* THIRD LOOP: APPLY CALIBRATION TO THE CALIBRATION FILES FOR TESTS */
  /********************************************************************/

  if (useTest)
    {
      
      MJExtractCalibTest testloop;
      
      // If you want to see the data-check plots only, choose:
      testloop.SetDataCheckDisplay();

      testloop.SetExtractor(&extractor);
      testloop.SetTimeExtractor(&timeext);
      testloop.SetInput(&cruns);
      testloop.SetPathOut(outpath);
      if (useDisplay)
        testloop.SetDisplay(display);
      testloop.SetBadPixels(calloop.GetBadPixels());
      
      if (!testloop.ProcessD(pedloop.GetPedestalCam(),calloop.GetCalibrationCam(),calloop.GetQECam()))
        return 3;
      
      if (useTimes)
        if (!testloop.ProcessT(pedloop.GetPedestalCam(),calloop.GetRelTimeCam()))
      return 4;
      
    }
  
  /********************************************************************/
  /* FOURTH LOOP: APPLY CALIBRATION TO THE PEDESTAL FILES             */
  /********************************************************************/

  MJExtractSignal pedphotloop;

  pedphotloop.SetExtractor(&extractor);
  pedphotloop.SetTimeExtractor(&timeext);
  pedphotloop.SetInput(&pruns);
  pedphotloop.SetOutputPath(outpath);
  if (useDisplay)
    pedphotloop.SetDisplay(display);
  pedphotloop.SetBadPixels(calloop.GetBadPixels());
  
  if (!pedphotloop.ProcessP(pedloop.GetPedestalCam(),calloop.GetCalibrationCam(),calloop.GetQECam()))
    return 5;


  if (debug)
    TObject::SetObjectStat(kFALSE);

  //
  // Debugging at the end:
  // 
  if (debug)
    gObjectTable->Print();

  //
  // List of useful containers:
  // 
/*
  MPedestalCam          &pedcam      = pedloop.GetPedestalCam();
  MCalibrationChargeCam &chargecam   = calloop.GetCalibrationCam();
  MCalibrationQECam     &qecam       = calloop.GetCalibrationCam();
  MBadPixelsCam         &badcam      = calloop.GetBadPixels();
  MCalibrationTestCam   &testcam     = testloop.GetTestCam();
  MHCalibrationTestTimeCam &testtime = testloop.GetTestTimeCam();
*/

  return 0;

}

static void Usage()
{
    gLog << endl;
    gLog << "Usage:" << endl;
    gLog << endl;
    gLog << "   calibration [ped.run nr.] [first cal.run nr.] [last cal.run nr.]" << endl ;
    gLog << endl;
    gLog << "   ped.run.nr:        Run number of the pedestal file." << endl;
    gLog << "   first cal.run nr.: Run number of the first calibration file." << endl;
    gLog << "   last  cal.run nr.: Run number of the last  calibration file." << endl;
    gLog << endl;
    gLog << "All calibration runs between (first cal.run nr.) and (last  cal.run nr.) will be used" << endl;
    gLog << "If last.cal.run.nr is 0 (default), only one calibration run is taken"                  << endl;
    gLog << endl;
    gLog << "Additional Options: " << endl;
    gLog << "     --inpath=#          Find the data in inpath"                      << endl;
    gLog << "     --outpath=#         Write the output containers to outpath"       << endl;
    gLog << "     --badfile=#         Use the file # to exclude pixels from the beginning" << endl;
    gLog << "     --debug             Use the TObjectTable for debugging    "       << endl;
    gLog << "                             and write out the pixels as ascii tables" << endl;
    gLog << "     --useTimes          Calibrate the relative arrival times"         << endl;
    gLog << "     --useTest           Use the class MJExtractCalibTest to test the calibration on itself" << endl;
    gLog << "     --skipBlindPix      Skip the blind pixel calibration"             << endl;
    gLog << "     --skipPINDiode      Skip the PIN Diode   calibration"             << endl;
}


int main(int argc, char **argv)
{
  //
  // Evaluate arguments
  //
  MArgs arg(argc, argv);
  
  if (arg.HasOnly("-?") || arg.HasOnly("-h") || arg.HasOnly("--help"))
    {
      Usage();
      return -1;
    }
  
  debug    = arg.HasOnlyAndRemove("--debug")    || arg.HasOnlyAndRemove("-d");
  useTimes = arg.HasOnlyAndRemove("--useTimes") || arg.HasOnlyAndRemove("-t");
  useTest  = arg.HasOnlyAndRemove("--useTest")  || arg.HasOnlyAndRemove("-e");
  blindpix = !(arg.HasOnlyAndRemove("--skipBlindPix"));
  pindiode = !(arg.HasOnlyAndRemove("--skipPINDiode"));

  if (arg.HasOption("--inpath="))
    inpath = arg.GetStringAndRemove("--inpath=");

  if (arg.HasOption("--outpath="))
    outpath = arg.GetStringAndRemove("--outpath=");

  if (arg.HasOption("--badfile="))
    badfile = arg.GetStringAndRemove("--badfile=");

  if (gSystem->AccessPathName(badfile,kFileExists))
  {
    gLog << "WARNING: the bad pixels file '" << badfile.Data() << "' doesn't exist." << endl;
    badfile = "";
  }

  // check for the right usage of the program
  //
  if (arg.GetNumArguments()>4)
    {
      Usage();
      return -1;
    }

  //
  // Initialize Non-GUI (batch) mode
  //
  gROOT->SetBatch();
  
  //
  // Switch off the display
  //
  useDisplay = kFALSE;


  //
  // check for the arguments
  //
  Int_t pedr  = 0;
  Int_t calr1 = 0;
  Int_t calr2 = 0;

  const Int_t nargs = arg.GetNumArguments();

  if (nargs>=3)
    {
      pedr = arg.GetArgumentInt(0);
      calr1 = arg.GetArgumentInt(1);
      calr2 = arg.GetArgumentInt(2);
      return calibration(pedr,calr1,calr2);
    }

  if (nargs>=2)
    {
      pedr = arg.GetArgumentInt(0);
      calr1 = arg.GetArgumentInt(1);
      return calibration(pedr,calr1);
    }

  if (nargs>=1)
    {
      pedr = arg.GetArgumentInt(0);
      gLog << "PEDR: " << pedr << endl;
      return calibration(pedr);
    }

  return calibration();
}