source: trunk/MagicSoft/Mars/mtemp/mifae/library/MCalibrateDC.cc@ 4061

/* ======================================================================== *\
!
! *
! * 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): Javier López , 5/2004 <mailto:jlopez@ifae.es>
!
!   Copyright: MAGIC Software Development, 2000-2004
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//
//  MCalibrateDC
//
/////////////////////////////////////////////////////////////////////////////
#include "MCalibrateDC.h"

#include <TString.h>
#include <TH1F.h>
#include <TF1.h>

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

#include "MGeomCam.h"
#include "MGeomPix.h"
#include "MCameraDC.h"
#include "MTime.h"

#include "MParList.h"

ClassImp(MCalibrateDC);
using namespace std;

MCalibrateDC::MCalibrateDC(TString filename, const char *name, const char *title)
{
  fName  = name  ? name  : "MCalibrateDC";
  fTitle = title ? title : "Taks to intercalibrate the DC of all pmts from a continuos light run";

  fStartingMissCalibration.Set(2004,3,8);
  fEndingMissCalibration.Set(2004,4,15);
  
  fDCMissCalibrationFactor = 30./40.96;
  fFileName = filename;
  
  Int_t   nbins = 120;
  Float_t min = 0;
  Float_t max = 30.;
  fCalHist = new TH1F("calhist","",nbins,min,max);

  fDCCalibration = 1.;
  fDCCalibrationRMS = 0.;
}
MCalibrateDC::~MCalibrateDC()
{
  delete  fCalHist;
}

Int_t MCalibrateDC::PreProcess(MParList *pList)
{

    fGeomCam = (MGeomCam*)pList->FindObject(AddSerialNumber("MGeomCam"));

    if (!fGeomCam)
    {
      *fLog << err << AddSerialNumber("MGeomCam") << " not found ... aborting" << endl;
      return kFALSE;
    }

    fCurr = (MCameraDC*)pList->FindObject(AddSerialNumber("MCameraDC"));

    if (!fCurr)
    {
      *fLog << err << AddSerialNumber("MCameraDC") << " not found ... aborting" << endl;
      return kFALSE;
    }

    fTimeCurr = (MTime*)pList->FindObject(AddSerialNumber("MTimeCurrents"));

    if (!fTimeCurr)
    {
      *fLog << err << AddSerialNumber("MTimeCurrents") << " not found ... aborting" << endl;
      return kFALSE;
    }

    // Run over the continuos light run to get the DC intercalibration factors
    fNumPixels = fGeomCam->GetNumPixels();
    fDCCalibrationFactor.Set(fNumPixels);
    fDCCalibrationFactor.Reset(1.);

    if ( fFileName != "" )
      {
        // FIXME::Here it should be done the loop over the continuos light run
        DCCalibrationCalc();
 
        for (UInt_t pix=1; pix<fNumPixels; pix++)
          fDCCalibrationFactor[pix] = fDCCalibration/fDisplay[pix];

      }

    return kTRUE;
}

Int_t MCalibrateDC::Process()
{

  if (*fTimeCurr >= fStartingMissCalibration && *fTimeCurr <= fEndingMissCalibration)
    {
      for (UInt_t pix=1; pix<fNumPixels; pix++)
        {
          MGeomPix& pixel = (*fGeomCam)[pix];
          if (pixel.GetSector() >=3 && pixel.GetSector() <=5)
            fDCCalibrationFactor[pix]*=fDCMissCalibrationFactor;
        }
    }

  for (UInt_t pix=1; pix<fNumPixels; pix++)
    (*fCurr)[pix]= (*fCurr)[pix]*fDCCalibrationFactor[pix];
          
  return kTRUE;
}

Bool_t MCalibrateDC::DCCalibrationCalc()
{

   for (UInt_t pix=1; pix<fNumPixels; pix++)
       fCalHist->Fill(fDisplay.GetBinContent(pix+1));

   Float_t nummaxprobdc = fCalHist->GetBinContent(fCalHist->GetMaximumBin());
   Float_t maxprobdc = fCalHist->GetBinCenter(fCalHist->GetMaximumBin());
   UInt_t bin = fCalHist->GetMaximumBin();
   do
   {
       bin++;
   }
   while(fCalHist->GetBinContent(bin)/nummaxprobdc > 0.5);
   Float_t halfmaxprobdc = fCalHist->GetBinCenter(bin);

   *fLog << dbg << " maxprobdc[high] " << maxprobdc << "[" << nummaxprobdc << "] ";
   *fLog << dbg << " halfmaxprobdc[high] " << halfmaxprobdc << "[" << fCalHist->GetBinContent(bin) << "]" << endl;

   Float_t rmsguess = TMath::Abs(maxprobdc-halfmaxprobdc);
   Float_t min = maxprobdc-3*rmsguess;
   min = (min<0.?0.:min);
   Float_t max = maxprobdc+3*rmsguess;

   *fLog << dbg << " maxprobdc " << maxprobdc << " rmsguess " << rmsguess << endl;

   TF1 func("func","gaus",min,max);
   func.SetParameters(nummaxprobdc, maxprobdc, rmsguess);
   
   fCalHist->Fit("func","QR0");

   UInt_t aproxnumdegrees = 6*(bin-fCalHist->GetMaximumBin());
   Float_t chiq = func.GetChisquare();
   fDCCalibration = func.GetParameter(1);
   fDCCalibrationRMS = func.GetParameter(2);

   *fLog << dbg << " fDCCalibration " << fDCCalibration << " fDCCalibrationRMS " << fDCCalibrationRMS << " chiq/ndof " << chiq << "/" << aproxnumdegrees << endl;

   Int_t numsigmas = 5;
   Axis_t minbin = fDCCalibration-numsigmas*fDCCalibrationRMS/fCalHist->GetBinWidth(1);
   minbin=minbin<1?1:minbin;
   Axis_t maxbin = fDCCalibration+numsigmas*fDCCalibrationRMS/fCalHist->GetBinWidth(1);
   *fLog << dbg << " Number of pixels with dc under " << numsigmas << " sigmas = " << fCalHist->Integral((int)minbin,(int)maxbin) << endl;

    //Check results from the fit are consistent
    if (TMath::Abs(fDCCalibration-maxprobdc) > rmsguess || fDCCalibrationRMS > rmsguess)
      {
        *fLog << warn << GetName() << " Calibration DC fit give non consistent results." << endl;
        *fLog << warn << " maxprobdc " << maxprobdc << " rmsguess " << rmsguess << endl;
        *fLog << warn << " fDCCalibration " << fDCCalibration << " fDCCalibrationRMS " << fDCCalibrationRMS << " chiq/ndof " << chiq << "/" << aproxnumdegrees << endl;
        fDCCalibration = maxprobdc;
        fDCCalibrationRMS = rmsguess/1.175; // FWHM=2.35*rms
      }

   return kTRUE;
}