source: trunk/MagicSoft/Mars/mhcalib/MHCalibrationPulseTimeCam.cc@ 6690

/* ======================================================================== *\
!
! *
! * 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   02/2004 <mailto:markus@ifae.es>
!
!   Copyright: MAGIC Software Development, 2000-2004
!
!
\* ======================================================================== */
/////////////////////////////////////////////////////////////////////////////
//                                                               
// MHCalibrationPulseTimeCam                                               
//
// Fills the extracted signals of MExtractedSignalCam into the MHCalibrationPix-classes 
// MHCalibrationPulseTimeHiGainPix and MHCalibrationPulseTimeLoGainPix for every:
//
// - Pixel, stored in the TOrdCollection's MHCalibrationCam::fHiGainArray and 
//   MHCalibrationCam::fLoGainArray
//
// - Average pixel per AREA index (e.g. inner and outer for the MAGIC camera), 
//   stored in the TOrdCollection's MHCalibrationCam::fAverageHiGainAreas and 
//   MHCalibrationCam::fAverageLoGainAreas
//
// - Average pixel per camera SECTOR (e.g. sectors 1-6 for the MAGIC camera), 
//   stored in the TOrdCollection's MHCalibrationCam::fAverageHiGainSectors and 
//   MHCalibrationCam::fAverageLoGainSectors
// 
// Every signal is taken from MExtractedSignalCam and filled into a histogram and 
// an array, in order to perform a Fourier analysis (see MHGausEvents). 
// The signals are moreover averaged on an event-by-event basis and written into 
// the corresponding average pixels.
//
// Additionally, the (FADC slice) position of the maximum is stored in an Absolute 
// Arrival Time histogram. This histogram serves for a rough cross-check if the 
// signal does not lie at or outside the edges of the extraction window. 
//
// The PulseTime histograms are fitted to a Gaussian, mean and sigma with its errors 
// and the fit probability are extracted. If none of these values are NaN's and 
// if the probability is bigger than MHGausEvents::fProbLimit (default: 0.5%), 
// the fit is declared valid.
// Otherwise, the fit is repeated within ranges of the previous mean 
// +- MHCalibrationPix::fPickupLimit (default: 5) sigma (see MHCalibrationPix::RepeatFit())
// In case this does not make the fit valid, the histogram means and RMS's are 
// taken directly (see MHCalibrationPix::BypassFit()) and the following flags are set:
// - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kHiGainNotFitted ) or  
// - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kLoGainNotFitted ) and 
// - MBadPixelsPix::SetUnsuitable(   MBadPixelsPix::kUnreliableRun   ) 
// 
// Outliers of more than MHCalibrationPix::fPickupLimit (default: 5) sigmas 
// from the mean are counted as Pickup events (stored in MHCalibrationPix::fPickup) 
//
// Unless more than fNumHiGainSaturationLimit (default: 1%) of the overall FADC 
// slices show saturation, the following flag is set:
// - MCalibrationPulseTimePix::SetHiGainSaturation();
// In that case, the calibration constants are derived from the low-gain results.
//
// If more than fNumLoGainSaturationLimit (default: 1%) of the overall 
// low-gain FADC slices saturate, the following flags are set:
// - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kLoGainSaturation ) and
// - MBadPixelsPix::SetUnsuitable(   MBadPixelsPix::kUnsuitableRun    )
// 
// The class also fills arrays with the signal vs. event number, creates a fourier 
// spectrum and investigates if the projected fourier components follow an exponential 
// distribution. In case that the probability of the exponential fit is less than 
// MHGausEvents::fProbLimit (default: 0.5%), the following flags are set:
// - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kHiGainOscillating ) or
// - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kLoGainOscillating ) and
// - MBadPixelsPix::SetUnsuitable(   MBadPixelsPix::kUnreliableRun     )
// 
// This same procedure is performed for the average pixels.
//
// The following results are written into MCalibrationPulseTimeCam:
//
// - MCalibrationPix::SetHiGainSaturation() 
// - MCalibrationPix::SetHiGainMean()
// - MCalibrationPix::SetHiGainMeanErr()
// - MCalibrationPix::SetHiGainSigma()
// - MCalibrationPix::SetHiGainSigmaErr()
// - MCalibrationPix::SetHiGainProb()
// - MCalibrationPix::SetHiGainNumPickup()
//
// For all averaged areas, the fitted sigma is multiplied with the square root of 
// the number involved pixels in order to be able to compare it to the average of 
// sigmas in the camera.
//
/////////////////////////////////////////////////////////////////////////////
#include "MHCalibrationPulseTimeCam.h"
#include "MHCalibrationCam.h"

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

#include "MParList.h"

#include "MHCalibrationPix.h"

#include "MCalibrationIntensityCam.h"
#include "MCalibrationChargeCam.h"
#include "MCalibrationChargePix.h"

#include "MGeomCam.h"
#include "MGeomPix.h"

#include "MBadPixelsIntensityCam.h"
#include "MBadPixelsCam.h"
#include "MBadPixelsPix.h"

#include "MRawEvtData.h"
#include "MRawRunHeader.h"
#include "MRawEvtPixelIter.h"

#include "MExtractedSignalCam.h"
#include "MExtractedSignalPix.h"

#include "MArrayI.h"
#include "MArrayD.h"

#include <TOrdCollection.h>
#include <TPad.h>
#include <TVirtualPad.h>
#include <TCanvas.h>
#include <TStyle.h>
#include <TF1.h>
#include <TLatex.h>
#include <TLegend.h>
#include <TGraph.h>
#include <TEnv.h>

ClassImp(MHCalibrationPulseTimeCam);

using namespace std;

const Int_t   MHCalibrationPulseTimeCam::fgHiGainNbins =  40;
const Axis_t  MHCalibrationPulseTimeCam::fgHiGainFirst =  -0.5;
const Axis_t  MHCalibrationPulseTimeCam::fgHiGainLast  =  19.5;
const Float_t MHCalibrationPulseTimeCam::fgProbLimit   = 0.0001;
const Byte_t  MHCalibrationPulseTimeCam::fgLowerSignalLimit = 100;
const TString MHCalibrationPulseTimeCam::gsHistName    = "PulseTime";
const TString MHCalibrationPulseTimeCam::gsHistTitle   = "Extracted Times";
const TString MHCalibrationPulseTimeCam::gsHistXTitle  = "Time [FADC slices]";
const TString MHCalibrationPulseTimeCam::gsHistYTitle  = "Nr. events";
const TString MHCalibrationPulseTimeCam::fgReferenceFile = "mjobs/calibrationref.rc";
// --------------------------------------------------------------------------
//
// Default Constructor. 
//
// Sets:
// - all pointers to NULL
//
// - fNbins to fgHiGainNbins
// - fFirst to fgHiGainFirst
// - fLast  to fgHiGainLast 
//
// - fHistName   to gsHistName  
// - fHistTitle  to gsHistTitle 
// - fHistXTitle to gsHistXTitle
// - fHistYTitle to gsHistYTitle
//
MHCalibrationPulseTimeCam::MHCalibrationPulseTimeCam(const char *name, const char *title)
    : fRawEvt(NULL)
{

  fName  = name  ? name  : "MHCalibrationPulseTimeCam";
  fTitle = title ? title : "Class to fill the extracted pulse times for cosmics ";
  
  SetNbins(fgHiGainNbins);
  SetFirst(fgHiGainFirst);
  SetLast (fgHiGainLast );

  SetProbLimit(fgProbLimit);

  SetHistName  (gsHistName  .Data());
  SetHistTitle (gsHistTitle .Data());
  SetHistXTitle(gsHistXTitle.Data());
  SetHistYTitle(gsHistYTitle.Data());

  SetReferenceFile();
  SetLoGain(kFALSE);
  SetOscillations(kFALSE);

  SetLowerSignalLimit();

  fInnerRefTime = 5.;
  fOuterRefTime = 5.;
}
// --------------------------------------------------------------------------
//
// Creates new MHCalibrationPulseTimeCam only with the averaged areas:
// the rest has to be retrieved directly, e.g. via: 
//     MHCalibrationPulseTimeCam *cam = MParList::FindObject("MHCalibrationPulseTimeCam");
//  -  cam->GetAverageSector(5).DrawClone();
//  -  (*cam)[100].DrawClone()
//
TObject *MHCalibrationPulseTimeCam::Clone(const char *) const
{

  MHCalibrationPulseTimeCam *cam = new MHCalibrationPulseTimeCam();

  //
  // Copy the data members
  //
  cam->fColor                  = fColor;
  cam->fRunNumbers             = fRunNumbers;
  cam->fPulserFrequency        = fPulserFrequency;
  cam->fFlags                  = fFlags;
  cam->fNbins                  = fNbins;
  cam->fFirst                  = fFirst;
  cam->fLast                   = fLast;
  cam->fReferenceFile          = fReferenceFile;

  if (!IsAverageing())
    return cam;

  const Int_t navhi   =  fAverageHiGainAreas->GetSize();

  for (int i=0; i<navhi; i++)
    cam->fAverageHiGainAreas->AddAt(GetAverageHiGainArea(i).Clone(),i);

  return cam;
}

// --------------------------------------------------------------------------
//
// Gets the pointers to:
// - MRawEvtData
//
Bool_t MHCalibrationPulseTimeCam::SetupHists(const MParList *pList)
{

  fRawEvt = (MRawEvtData*)pList->FindObject("MRawEvtData");
  if (!fRawEvt)
  {
      *fLog << err << dbginf << "MRawEvtData not found... aborting." << endl;
      return kFALSE;
  }

  return kTRUE;
}

// --------------------------------------------------------------------------
//
// Gets or creates the pointers to:
// - MExtractedSignalCam
// - MCalibrationPulseTimeCam or MCalibrationIntensityPulseTimeCam
// - MBadPixelsCam
//
// Initializes the number of used FADC slices from MExtractedSignalCam 
// into MCalibrationPulseTimeCam and test for changes in that variable
//
// Calls:
// - InitHiGainArrays()
// 
// Sets:
// - fSumhiarea   to nareas 
// - fSumloarea   to nareas 
// - fSumhisector to nsectors 
// - fSumlosector to nsectors 
//
Bool_t MHCalibrationPulseTimeCam::ReInitHists(MParList *pList)
{

  MExtractedSignalCam *signal = 
    (MExtractedSignalCam*)pList->FindObject(AddSerialNumber("MExtractedSignalCam"));
  if (!signal)
  {
      *fLog << err << "MExtractedSignalCam not found... abort." << endl;
      return kFALSE;
  }

  if (!InitCams(pList,"PulseTime"))
    return kFALSE;

  const Int_t npixels  = fGeom->GetNumPixels();
  const Int_t nsectors = fGeom->GetNumSectors();
  const Int_t nareas   = fGeom->GetNumAreas();

  InitHiGainArrays(npixels,nareas,nsectors);

  fSumhiarea  .Set(nareas); 
  fSumhisector.Set(nsectors); 

  return kTRUE;
}


// --------------------------------------------------------------------------
//
// Retrieves from MExtractedSignalCam:
// - first used LoGain FADC slice
//
// Retrieves from MGeomCam:
// - number of pixels
// - number of pixel areas
// - number of sectors
//
// For all TOrdCollection's (including the averaged ones), the following steps are performed: 
//
// 1) Fill PulseTimes histograms (MHGausEvents::FillHistAndArray()) with:
// - MExtractedSignalPix::GetExtractedSignalHiGain();
// - MExtractedSignalPix::GetExtractedSignalLoGain();
//
Bool_t MHCalibrationPulseTimeCam::FillHists(const MParContainer *par, const Stat_t w)
{

  MExtractedSignalCam *signal = (MExtractedSignalCam*)par;
  if (!signal)
    {
      *fLog << err << "No argument in MExtractedSignalCam::Fill... abort." << endl;
      return kFALSE;
    }
  
  const UInt_t nareas   = fGeom->GetNumAreas();
  const UInt_t nsectors = fGeom->GetNumSectors();

  fSumhiarea  .Reset(); 
  fSumhisector.Reset(); 

  MRawEvtPixelIter pixel(fRawEvt);
  while (pixel.Next())
    {
      
      const Int_t i = pixel.GetPixelId();

      MHCalibrationPix &histhi = (*this)[i];

      if (histhi.IsExcluded())
        continue;

      const MExtractedSignalPix &pix = (*signal)[i];
      
      const Int_t sathi = (Int_t)pix.GetNumHiGainSaturated();

      if (sathi)
        continue;

      Byte_t *start = pixel.GetHiGainSamples();
      Byte_t *end   = start + pixel.GetNumHiGainSamples();
      Byte_t *p     = start;
      Byte_t max    = 0;
      Int_t  maxpos = 0;

      while (p < end)
        {
          if (*p > max)
            {
              max = *p;
              maxpos = p-start-1;
            }
          p++;
        }

      start   = pixel.GetLoGainSamples();
      end     = start + pixel.GetNumLoGainSamples();
      p       = start;

      while (p < end)
        {
          if (*p > max)
            {
              max = *p;
              maxpos = p-start+pixel.GetNumHiGainSamples() - 1;
            }
          p++;
        }

      if (max < fLowerSignalLimit)
        continue;

      const Float_t time = (Float_t)maxpos;
      histhi.FillHist(time);
        
      //      *fLog << inf << time << endl;

      const Int_t aidx   = (*fGeom)[i].GetAidx();
      const Int_t sector = (*fGeom)[i].GetSector();

      fSumhiarea[aidx]      += time;
      fSumhisector[sector]  += time;

    }
  
  for (UInt_t j=0; j<nareas; j++)
    {

      const Int_t npix = fAverageAreaNum[j];

      if (npix == 0)
        continue;

      MHCalibrationPix &hipix = GetAverageHiGainArea(j);
      hipix.FillHist(fSumhiarea[j]/npix);
      
    }

  for (UInt_t j=0; j<nsectors; j++)
    {

      const Int_t npix = fAverageSectorNum[j];

      if (npix == 0)
        continue;

      MHCalibrationPix &hipix = GetAverageHiGainSector(j);
      hipix.FillHist(fSumhisector [j]/npix);
    }

  return kTRUE;
}

// --------------------------------------------------------------------------
//
// For all TOrdCollection's (including the averaged ones), the following steps are performed: 
//
// 1) Returns if the pixel is excluded.
// 2) Tests saturation. In case yes, set the flag: MCalibrationPix::SetHiGainSaturation()
//    or the flag: MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kLoGainSaturated )
// 3) Store the absolute arrival times in the MCalibrationPulseTimePix's. If flag 
//    MCalibrationPix::IsHiGainSaturation() is set, the Low-Gain arrival times are stored, 
//    otherwise the Hi-Gain ones.
// 4) Calls to MHCalibrationCam::FitHiGainArrays() and MCalibrationCam::FitLoGainArrays() 
//    with the flags:
//    - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kHiGainNotFitted )
//    - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kLoGainNotFitted )
//    - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kHiGainOscillating )
//    - MBadPixelsPix::SetUncalibrated( MBadPixelsPix::kLoGainOscillating )
//
Bool_t MHCalibrationPulseTimeCam::FinalizeHists()
{

  *fLog << endl;

  MCalibrationCam *calcam  = fIntensCam ? fIntensCam->GetCam() : fCam;
  //
  // Perform the fitting for the High Gain (done in MHCalibrationCam)
  // 
  for (Int_t i=0; i<fHiGainArray->GetSize(); i++)
    {
      
      MHCalibrationPix &hist = (*this)[i];
      
      if (hist.IsExcluded())
        continue;
      
      MCalibrationPix &pix    = (*calcam)[i];
      CalcHists(hist,pix);
    }

  if (!IsAverageing())
    return kTRUE;
  
  for (Int_t j=0; j<fAverageHiGainAreas->GetSize(); j++)
    {
      
      MHCalibrationPix &hist = GetAverageHiGainArea(j);      
      MCalibrationPix  &pix  = calcam->GetAverageArea(j);
      CalcHists(hist,pix);
  }
  
  for (Int_t j=0; j<fAverageHiGainSectors->GetSize(); j++)
    {
      MHCalibrationPix &hist = GetAverageHiGainSector(j);      
      MCalibrationPix  &pix  = calcam->GetAverageSector(j);
      CalcHists(hist,pix);
    }
  
  return kTRUE;
}

void MHCalibrationPulseTimeCam::CalcHists(MHCalibrationPix &hist, MCalibrationPix &pix)
{


  if (hist.IsEmpty() || hist.IsOnlyOverflow() || hist.IsOnlyUnderflow())
    {
      *fLog << warn << GetDescriptor()
            << ": Only overflow or underflow in hi-gain pixel: " << pix.GetPixId() << endl;
      return;
    }  

  hist.BypassFit();

  pix.SetHiGainMean       ( hist.GetMean()      );
  pix.SetHiGainMeanVar    ( hist.GetMeanErr() * hist.GetMeanErr()   );
  pix.SetHiGainRms        ( hist.GetHistRms()   );
  pix.SetHiGainSigma      ( hist.GetSigma()     );
  pix.SetHiGainSigmaVar   ( hist.GetSigmaErr()* hist.GetSigmaErr()  );

  if (IsDebug())
    {
      *fLog << dbginf << GetDescriptor() << ": ID " << GetName() 
            << " Mean: "         << hist.GetMean    ()
            << " MeanErr: "      << hist.GetMeanErr ()
            << " MeanSigma: "    << hist.GetSigma   ()
            << " MeanSigmaErr: " << hist.GetSigmaErr()
            << " Prob: "         << hist.GetProb    ()
            << endl;
    }
  return;
}


// --------------------------------------------------------------------------
//
// Calls MHCalibrationPix::DrawClone() for pixel idx
//
void MHCalibrationPulseTimeCam::DrawPixelContent(Int_t idx) const
{
  (*this)[idx].DrawClone();
}


// -----------------------------------------------------------------------------
// 
// Default draw:
//
// Displays the averaged areas, both High Gain and Low Gain 
//
// Calls the Draw of the fAverageHiGainAreas and fAverageLoGainAreas objects with options
//
void MHCalibrationPulseTimeCam::Draw(const Option_t *opt)
{

  const Int_t nareas = fAverageHiGainAreas->GetSize();
  if (nareas == 0)
    return;

  TString option(opt);
  option.ToLower();

  // 
  // From here on , the datacheck - Draw
  //
  TVirtualPad *pad = gPad ? gPad : MH::MakeDefCanvas(this);  
  pad->SetBorderMode(0);
  pad->Divide(1,nareas);

  //
  // Loop over inner and outer pixels
  //  
  for (Int_t i=0; i<nareas;i++) 
     {
       
       pad->cd(i+1);
       
       MHCalibrationPix &hipix = GetAverageHiGainArea(i);
       DrawDataCheckPixel(hipix,i ? fOuterRefTime : fInnerRefTime);
    }      
}

// -----------------------------------------------------------------------------
// 
// Draw the average pixel for the datacheck:
//
// Displays the averaged areas, both High Gain and Low Gain 
//
// Calls the Draw of the fAverageHiGainAreas and fAverageLoGainAreas objects with options
//
void MHCalibrationPulseTimeCam::DrawDataCheckPixel(MHCalibrationPix &pix, const Float_t refline)
{
  
  gPad->SetTicks();
  if (!pix.IsEmpty() && !pix.IsOnlyOverflow() && !pix.IsOnlyUnderflow())
    gPad->SetLogy();

  TH1F *hist = pix.GetHGausHist();

  TH1F *null = new TH1F("Null",hist->GetTitle(),100,
                        pix.GetFirst() > 0. ? pix.GetFirst() : 0.,
                        pix.GetLast() > pix.GetFirst()
                        ? ( pix.GetLast() > 450. 
                            ? ( fColor == MCalibrationCam::kBLUE ? 800. : 450. )
                            : pix.GetLast() )
                        : pix.GetFirst()*2.);

  null->SetMaximum(1.1*hist->GetMaximum());
  null->SetDirectory(NULL);
  null->SetBit(kCanDelete);
  null->SetStats(kFALSE);
  //
  // set the labels bigger
  //
  TAxis *xaxe = null->GetXaxis();
  TAxis *yaxe = null->GetYaxis();
  xaxe->CenterTitle();
  yaxe->CenterTitle();    
  xaxe->SetTitleSize(0.07);
  yaxe->SetTitleSize(0.07);    
  xaxe->SetTitleOffset(0.7);
  yaxe->SetTitleOffset(0.55);    
  xaxe->SetLabelSize(0.06);
  yaxe->SetLabelSize(0.06);    
  xaxe->SetTitle(hist->GetXaxis()->GetTitle());
  yaxe->SetTitle(hist->GetYaxis()->GetTitle());  

  null->Draw();
  hist->Draw("same");

  gStyle->SetOptFit();

  TF1 *fit = pix.GetFGausFit();

  if (fit)
  {
    switch ( fColor )
      {
      case MCalibrationCam::kGREEN:
        fit->SetLineColor(kGreen);
        break;
      case MCalibrationCam::kBLUE:
        fit->SetLineColor(kBlue);
        break;
      case MCalibrationCam::kUV:  
        fit->SetLineColor(106);
        break;
      case MCalibrationCam::kCT1: 
        fit->SetLineColor(006);
        break;
      default:
        fit->SetLineColor(kRed);
      }
    fit->Draw("same");
  }

  DisplayRefLines(null,refline);

  return;
  
}

void  MHCalibrationPulseTimeCam::DisplayRefLines(const TH1F *hist, const Float_t refline) const
{

  TGraph *uv10    = new TGraph(2);
  uv10->SetPoint(0,refline,0.);
  uv10->SetPoint(1,refline,hist->GetMaximum());
  uv10->SetBit(kCanDelete);
  uv10->SetLineColor(106);
  uv10->SetLineStyle(2);
  uv10->SetLineWidth(3);
  uv10->Draw("L");

  TLegend *leg = new TLegend(0.8,0.55,0.99,0.99);
  leg->SetBit(kCanDelete);
  leg->AddEntry(uv10,"10 Leds UV","l");

  leg->Draw();
}

Int_t MHCalibrationPulseTimeCam::ReadEnv(const TEnv &env, TString prefix, Bool_t print)
{

  Bool_t rc = kFALSE;

  if (MHCalibrationCam::ReadEnv(env,prefix,print))
    rc = kTRUE;
   
  if (IsEnvDefined(env, prefix, "LowerSignalLimit", print))
    {
      SetLowerSignalLimit(GetEnvValue(env,prefix,"LowerSignalLimit",fLowerSignalLimit));
      rc = kTRUE;
    }

  if (IsEnvDefined(env, prefix, "ReferenceFile", print))
    {
      SetReferenceFile(GetEnvValue(env,prefix,"ReferenceFile",fReferenceFile.Data()));
      rc = kTRUE;
    }

  TEnv refenv(fReferenceFile);

  fInnerRefTime = refenv.GetValue("InnerRefTime",fInnerRefTime);
  fOuterRefTime = refenv.GetValue("OuterRefTime",fOuterRefTime);

  return rc;
}