source: trunk/MagicSoft/Mars/manalysis/MMcCalibrationUpdate.cc@ 3534

/* ======================================================================== *\
!
! *
! * 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): Abelardo Moralejo, 12/2003 <mailto:moralejo@pd.infn.it>
!
!   Copyright: MAGIC Software Development, 2000-2003
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//
//  MMcCalibrationUpdate
//
//  This task looks for the ìnformation about FADC pedestals in
//  MMcFadcHeader and translates it to the pedestal mean and rms (in adc counts).
//  If not already existing in the parameter list, an MCalibrationCam object
//  is created, with the conversion factor between photons and ADC counts is 
//  set to 1 to allow the analysis to proceed.
//
//  Then it creates and fills also the MPedPhotCam object containing the pedestal
//  mean and rms in units of photons.
//
//  Input Containers:
//   MMcFadcHeader
//   MRawRunHeader
//  [MCalibrationCam] (if it existed previously)
//
//  Output Containers:
//   MPedPhotCam
//  [MCalibrationCam] (if it did not exist previously)
//
/////////////////////////////////////////////////////////////////////////////
#include "MMcCalibrationUpdate.h"

#include "MParList.h"

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

#include "MCalibrationChargePix.h"
#include "MCalibrationChargeCam.h"
#include "MExtractedSignalCam.h"
#include "MExtractedSignalPix.h"
#include "MGeomCam.h"
#include "MPedPhotCam.h"
#include "MPedPhotPix.h"

#include "MRawRunHeader.h"
#include "MMcRunHeader.hxx"
#include "MMcFadcHeader.hxx"

ClassImp(MMcCalibrationUpdate);

using namespace std;

MMcCalibrationUpdate::MMcCalibrationUpdate(const char *name, const char *title)
{
    fName  = name  ? name  : "MMcCalibrationUpdate";
    fTitle = title ? title : "Write MC pedestals and conversion factors into MCalibration Container";

    fADC2PhInner = 1.;
    fADC2PhOuter = 1.;

    fAmplitude = -1.;
    fAmplitudeOuter = -1.;
    fConversionHiLo = -1.;

    fFillCalibrationCam = kTRUE;
    fOuterPixelsGainScaling = kTRUE;
}

// --------------------------------------------------------------------------
//
// Check for the run type. Return kTRUE if it is a MC run or if there
// is no MC run header (old camera files) kFALSE in case of a different
// run type
//
Bool_t MMcCalibrationUpdate::CheckRunType(MParList *pList) const
{
    const MRawRunHeader *run = (MRawRunHeader*)pList->FindObject("MRawRunHeader");
    if (!run)
    {
        *fLog << warn << dbginf << "Warning - cannot check file type, MRawRunHeader not found." << endl;
        return kTRUE;
    }

    return run->IsMonteCarloRun();
}

// --------------------------------------------------------------------------
//
// Make sure, that there is an MCalibrationCam Object in the Parameter List.
//
Int_t MMcCalibrationUpdate::PreProcess(MParList *pList)
{
    fCalCam = (MCalibrationChargeCam*) pList->FindObject(AddSerialNumber("MCalibrationChargeCam"));
    if (!fCalCam)
    {
        fCalCam = (MCalibrationChargeCam*) pList->FindCreateObj(AddSerialNumber("MCalibrationChargeCam"));
        if (!fCalCam)
            return kFALSE;
    }
    else
    {
        fFillCalibrationCam = kFALSE;
        *fLog << inf << AddSerialNumber("MCalibrationChargeCam") << " already exists... " << endl;
    }

    fPedPhotCam = (MPedPhotCam*) pList->FindCreateObj(AddSerialNumber("MPedPhotCam"));
    if (!fPedPhotCam)
        return kFALSE;

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

    return kTRUE;
}

// --------------------------------------------------------------------------
//
// Check for the runtype.
// Search for MGeomCam and MMcFadcHeader.
// Fill the MCalibrationCam object.
//
Bool_t MMcCalibrationUpdate::ReInit(MParList *pList)
{
    //
    // If it is no MC file skip this function...
    //
    if (!CheckRunType(pList))
    {
        *fLog << inf << "This is no MC file... skipping." << endl;
        return kTRUE;
    }
        
    //
    // Now check the existence of all necessary containers.
    //
    fGeom = (MGeomCam*) pList->FindObject(AddSerialNumber("MGeomCam"));
    if (!fGeom)
    {
        *fLog << err << AddSerialNumber("MGeomCam") << " not found... aborting." << endl;
        return kFALSE;
    }

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

    MMcRunHeader* mcrunh = (MMcRunHeader*) pList->FindObject("MMcRunHeader");

    //
    // Initialize Fadc simulation parameters:
    //
    if (fAmplitude < 0)
    {
        fAmplitude = fHeaderFadc->GetAmplitud();
        if (mcrunh->GetCamVersion() > 60)
          {
            fAmplitudeOuter = fHeaderFadc->GetAmplitudOuter();
            fConversionHiLo = fHeaderFadc->GetLow2HighGain();
          }
        else // old MC files, camera  < v0.7
          {
            fAmplitudeOuter = fAmplitude;
            fConversionHiLo = 1;  // dummy value
          }

    }
    else   // Check that following files have all the same FADC parameters
    {
      if ( fabs(fHeaderFadc->GetAmplitud()-fAmplitude) > 1.e-6 )
        {
          *fLog << err << "Parameters of MMcFadcHeader are not the same for all files... aborting." << endl;
          return kFALSE;
        }

      if (mcrunh->GetCamVersion() > 60) // old MC files, camera < v0.7
        {
          if( fabs(fHeaderFadc->GetAmplitudOuter()-fAmplitudeOuter) > 1.e-6  ||
              fabs(fConversionHiLo-fHeaderFadc->GetLow2HighGain())  > 1.e-6 )
            {
              *fLog << err << "Parameters of MMcFadcHeader are not the same for all files... aborting." << endl;
              return kFALSE;
            }
        }
    }



    //
    // If MCalibrationCam already existed in the parameter list before
    // MMcCalibrationUpdate::PreProcess was executed (from a 
    // previous calibration loop) we must not fill it, hence nothing
    // else has to be done in ReInit:
    //
    if (!fFillCalibrationCam)
        return kTRUE;

    //
    // Set the ADC to photons conversion factor for outer pixels.
    // One can choose not to apply the known (in MC) gain factor between
    // inner and outer pixels, (fOuterPixelsGainScaling = kFALSE),
    // which may be useful for display purposes.
    //

    if (fOuterPixelsGainScaling)
      fADC2PhOuter = fADC2PhInner * (fAmplitude / fAmplitudeOuter);
    else
      fADC2PhOuter = fADC2PhInner;


    const int num = fCalCam->GetSize();

    fCalCam->SetBlindPixelMethodValid(kTRUE);

    for (int i=0; i<num; i++)
    {
        MCalibrationChargePix &calpix = (*fCalCam)[i];

        calpix.SetBlindPixelMethodValid();

        calpix.SetConversionHiLo(fConversionHiLo);
        calpix.SetConversionHiLoErr(0.);         // FIXME ?

        //
        // Write conversion factor ADC to photons (different for inner
        // and outer pixels).
        //

        Float_t adc2phot = (fGeom->GetPixRatio(i) < fGeom->GetPixRatio(0))?
          fADC2PhOuter : fADC2PhInner;

        calpix.SetConversionBlindPixelMethod(adc2phot, 0., 0.);
    }


    return kTRUE;
}


// --------------------------------------------------------------------------
//
// Fill the MCerPhotPed object
//
Int_t MMcCalibrationUpdate::Process()
{
    const int num = fCalCam->GetSize();

    for (int i=0; i<num; i++)
    {
        MExtractedSignalPix &sigpix = (*fSignalCam)[i];

        //
        // ped mean and rms per pixel, in ADC counts, according to signal 
        // calculation (hi or low gain and number of integrated slices):
        //
        const Float_t pedestmean = sigpix.IsLoGainUsed()? 
          fSignalCam->GetNumUsedLoGainFADCSlices()*fHeaderFadc->GetPedestal(i) : 
          fSignalCam->GetNumUsedHiGainFADCSlices()*fHeaderFadc->GetPedestal(i);

        //
        // In some cases, depending on the camera simulation parameters, one can
        // have very little or no noise in the FADC. In the case the rms of  
        // pedestal is zero, the pixel will be cleaned out later in the image
        // cleaning. To avoid this problem,we set a default value of 0.01 ADC
        // counts for the RMS per slice:
        //
        const Double_t used = (Double_t)(sigpix.IsLoGainUsed() ?
                                         fSignalCam->GetNumUsedLoGainFADCSlices() :
                                         fSignalCam->GetNumUsedHiGainFADCSlices());

        const Float_t rms0 = sigpix.IsLoGainUsed() ?
            fHeaderFadc->GetPedestalRmsLow(i) :
            fHeaderFadc->GetPedestalRmsHigh(i);

        const Float_t pedestrms = TMath::Sqrt(used) * (rms0>0 ? rms0 : 0.01);

        //
        // Write mean pedestal and pedestal rms per pixel
        // in number of photons:
        //
        MPedPhotPix &pedpix = (*fPedPhotCam)[i];

        MCalibrationChargePix &calpix = (*fCalCam)[i];
        Float_t adc2phot = calpix.GetMeanConversionBlindPixelMethod();
        Float_t hi2lo    = calpix.GetConversionHiLo();

        if (sigpix.IsLoGainUsed())
            pedpix.Set(adc2phot*hi2lo*pedestmean, adc2phot*hi2lo*pedestrms);
        else
            pedpix.Set(adc2phot*pedestmean, adc2phot*pedestrms);

    }

    return kTRUE;
}