source: trunk/MagicSoft/Mars/manalysis/MSigmabar.cc@ 5200

/* ======================================================================== *\
!
! *
! * 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): Robert Wagner,   10/2002 <mailto:magicsoft@rwagner.de>
!   Author(s): Wolfgang Wittek, 01/2003 <mailto:wittek@mppmu.mpg.de>
!   Author(s): Thomas Bretz,    04/2003 <mailto:tbretz@astro.uni-wuerzburg.de>
!
!   Copyright: MAGIC Software Development, 2003
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//                                                                         //
// MSigmabar                                                               //
//                                                                         //
// This is the storage container to hold information about                 //
// "average" pedestal RMS                                                  //
//                                                                         //
// The "average" pedestal RMS is calculated as                             //
//     <pedRMS> = sqrt( sum_i( (pedRMS_i)^2/area_i ) / no.of pixels )      //
//                                                                         //
//     which is the sqrt of the average (pedRMS^2 per area)                //
//                                                                         //
/////////////////////////////////////////////////////////////////////////////
#include "MSigmabar.h"

#include <TMath.h>

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

#include "MParList.h"

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

#include "MCerPhotEvt.h"
#include "MCerPhotPix.h"

#include "MPedPhotCam.h"
#include "MPedPhotPix.h"

ClassImp(MSigmabar);

using namespace std;

// --------------------------------------------------------------------------
//
MSigmabar::MSigmabar(const char *name, const char *title)
{
    fName  = name  ? name  : "MSigmabar";
    fTitle = title ? title : "Storage container for Sigmabar";
}

// --------------------------------------------------------------------------
//
MSigmabar::~MSigmabar()
{
  // do nothing special.
}

void MSigmabar::Reset()
{
    fSigmabar      = -1;
    fInnerPixels   = -1;
    fOuterPixels   = -1;
    fSigmabarInner = -1;
    fSigmabarOuter = -1;

    memset(fSigmabarSector, 0, sizeof(fSigmabarSector));
    memset(fSigmabarInnerSector, 0, sizeof(fSigmabarInnerSector));
    memset(fSigmabarOuterSector, 0, sizeof(fSigmabarOuterSector));
}

// --------------------------------------------------------------------------
//
// Actual calculation of sigmabar. This is done for each of the six sectors
// separately due to their possibly different HV behavior. Also inner and
// outer pixels are treated separately.
//
// Preliminary! Works for CT1 test, for real MAGIC crosschecks still have
// to be done. Also implementation details will be updated, like 
// determination of sector to which a respective pixel belongs
//
Float_t MSigmabar::Calc(const MGeomCam &geom, const MPedPhotCam &ped, 
                        const MCerPhotEvt &evt)
{
    Int_t innerPixels[6];
    Int_t outerPixels[6];
    Float_t innerSum[6];
    Float_t outerSum[6];

    memset(innerPixels, 0, sizeof(innerPixels));
    memset(outerPixels, 0, sizeof(outerPixels));
    memset(innerSum, 0, sizeof(innerSum));
    memset(outerSum, 0, sizeof(outerSum));

    //
    // sum up sigma^2/area for each sector, 
    // separately for inner and outer region;
    //
    // consider all pixels with Cherenkov photon information
    // and require "Used"
    //

    const UInt_t npix = evt.GetNumPixels();

    //*fLog << "MSigmabar : npix = " << npix << endl;

    for (UInt_t i=0; i<npix; i++)
    {
        MCerPhotPix &cerpix = evt.operator[](i);
        if (!cerpix.IsPixelUsed())
            continue;

        /*
         if ( cerpix.GetNumPhotons() == 0 )
         {
         *fLog << "MSigmabar::Calc(); no.of photons is 0 for used pixel"
         << endl;
         continue;
         }
        */

        const Int_t idx = cerpix.GetPixId();
        if (idx == 0)
        {
          //*fLog << "MSigmabar : id = 0;  pixel '0' is used, ignore it" 
          //      << endl;
          continue;
        }

        // ratio is the area of pixel 0
        //          divided by the area of the current pixel
        const Double_t ratio = geom.GetPixRatio(idx);

        //*fLog << "pixel id, ratio = " << idx << ",  " << ratio << endl;

        const MGeomPix &gpix = geom[idx];

        // This is wrong : rounding has to be done on positive values
        //Int_t sector = (Int_t)(atan2(gpix.GetY(),gpix.GetX())*6 
        //               / (TMath::Pi()*2));
        //if (sector<0)
        //    sector+=6;

        Float_t sectorf = atan2(gpix.GetY(),gpix.GetX()) * 6.0 
                          / (TMath::Pi()*2);
        if (sectorf < 0.0)
          sectorf += 6.0;
        Int_t sector = (Int_t)sectorf;

        // count only those pixels which have a sigma != 0.0
        const Float_t sigma = ped[idx].GetRms();

        if ( sigma <= 0 )
            continue;

        //if (area < 1.5)
        if (ratio > 0.5)
        {
            innerPixels[sector]++;
            innerSum[sector]+= sigma*sigma * ratio;
        }
        else
        {
            outerPixels[sector]++;
            outerSum[sector]+= sigma*sigma * ratio;
        }
    }

    fInnerPixels   = 0;
    fOuterPixels   = 0;
    Double_t fSumInner = 0;
    Double_t fSumOuter = 0;
    for (UInt_t i=0; i<6; i++)
    {
        fSumInner      += innerSum[i];
        fInnerPixels   += innerPixels[i];
        fSumOuter      += outerSum[i];
        fOuterPixels   += outerPixels[i];
    }

    if (fInnerPixels > 0) fSigmabarInner = sqrt(fSumInner / fInnerPixels);
    if (fOuterPixels > 0) fSigmabarOuter = sqrt(fSumOuter / fOuterPixels);

    //
    // this is the sqrt of the average sigma^2/area 
    //
    fSigmabar = (fInnerPixels+fOuterPixels)<=0 ? 0:
                sqrt( (fSumInner+fSumOuter)/(fInnerPixels+fOuterPixels) );

    for (UInt_t i=0; i<6; i++)
    {
        const Double_t ip  = innerPixels[i];
        const Double_t op  = outerPixels[i];
        const Double_t iss = innerSum[i];
        const Double_t oss = outerSum[i];

        const Double_t sum = ip + op;
        fSigmabarInnerSector[i] = ip <=0 ? 0 :        sqrt(iss/ip);
        fSigmabarOuterSector[i] = op <=0 ? 0 :        sqrt(oss/op);
        fSigmabarSector[i]      = sum<=0 ? 0 : sqrt((iss+oss)/sum);
    }

    //TString opt = "";
    //Print(opt);

  return fSigmabarInner;
}

// --------------------------------------------------------------------------
//
void MSigmabar::Print(Option_t *) const
{
  *fLog << all << endl;
  *fLog << "Total number of inner pixels is " << fInnerPixels << endl;
  *fLog << "Total number of outer pixels is " << fOuterPixels << endl;
  *fLog << endl;

  *fLog << "Sigmabar     Overall : " << fSigmabar << "   ";
  *fLog << " Sectors: ";
  for (Int_t i=0;i<6;i++)
      *fLog << fSigmabarSector[i] << ", ";
  *fLog << endl;


  *fLog << "Sigmabar     Inner   : " << fSigmabarInner << "   ";
  *fLog << " Sectors: ";
  for (Int_t i=0;i<6;i++)
      *fLog << fSigmabarInnerSector[i] << ", ";
  *fLog << endl;


  *fLog << "Sigmabar     Outer   : " << fSigmabarOuter << "   ";
  *fLog << " Sectors: ";
  for (Int_t i=0;i<6;i++)
      *fLog << fSigmabarOuterSector[i] << ", ";
  *fLog << endl;

}