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

/* ======================================================================== *\
!
! *
! * 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>
!
!   Copyright: MAGIC Software Development, 2002
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//                                                                         //
// MSigmabar                                                               //
//                                                                         //
// This is the storage container to hold information about the mean sigma  //
// (aka Sigmabar) of all pedestals                                         //
//                                                                         //
/////////////////////////////////////////////////////////////////////////////
#include "MSigmabar.h"

#include <TMath.h>

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

#include "MGeomCam.h"
#include "MPedestalCam.h"
#include "MGeomPix.h"
#include "MPedestalPix.h"

ClassImp(MSigmabar);

MSigmabar::MSigmabar(const char *name, const char *title)
{
    fName  = name  ? name  : "MSigmabar";
    fTitle = title ? title : "Storage container for Sigmabar";
    
    fSigmabar = 0.0;
    fSigmabarInner = 0.0;
    fSigmabarOuter = 0.0;
    fRatioA = 0.0;

    fCalcPixNum=kTRUE;
}

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


// --------------------------------------------------------------------------
//
// 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
//
Bool_t MSigmabar::Calc(const MGeomCam &geom, const MPedestalCam &ped)
{
  const UInt_t npix = ped.GetSize();
  Float_t innerSquaredSum[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
  Float_t outerSquaredSum[6] = {0.0, 0.0, 0.0, 0.0, 0.0, 0.0};
  Int_t innerPixels[6] = {0,0,0,0,0,0};
  Int_t outerPixels[6] = {0,0,0,0,0,0};
  Int_t currentSector;
  Float_t angle;
  
  for (UInt_t i=0; i<npix;i++)
    {
      const MGeomPix gpix = geom[i];
      angle=((6*atan2(gpix.GetX(),gpix.GetY())/(2*TMath::Pi()))-0.5);
      if (angle<0) angle+=6;
      currentSector=(Int_t)angle;          
      geom.GetPixRatio(i) == 1 ? innerPixels[currentSector]++ : outerPixels[currentSector]++;   
      
      const MPedestalPix pix = ped[i]; 
      geom.GetPixRatio(i) == 1 ? innerSquaredSum[currentSector]+=(pix.GetSigma()*pix.GetSigma()) : outerSquaredSum[currentSector]+=((pix.GetSigma()*pix.GetSigma()) / geom.GetPixRatio(i));

      // Get once and forever the ratio of areas outer/inner pixels
      if (fRatioA && (geom.GetPixRatio(i) != 1)) fRatioA=geom.GetPixRatio(i);
    }

  // Overall Sigma
  fSigmabarInner=0; fInnerPixels=0;
  fSigmabarOuter=0; fOuterPixels=0;
  for (UInt_t i=0; i<6; i++) {
    fSigmabarInner+=innerSquaredSum[i];
    fInnerPixels  +=innerPixels[i];
    fSigmabarOuter+=outerSquaredSum[i];
    fOuterPixels  +=outerPixels[i];
  }

  fSigmabarInner/=fInnerPixels;
  if (fSigmabarOuter != 0) fSigmabarOuter/=fOuterPixels; 
  fSigmabar=sqrt(fSigmabarInner + fSigmabarOuter/( fOuterPixels==0 ? 1 : fRatioA)); 
  fSigmabarInner=sqrt(fSigmabarInner);
  fSigmabarOuter=sqrt(fSigmabarOuter);
  
  for (UInt_t i=0; i<6; i++) {
    fSigmabarInnerSector[i]=innerSquaredSum[i]/innerPixels[i];
    fSigmabarOuterSector[i]=outerSquaredSum[i]/( outerSquaredSum[i]==0 ? 1: outerPixels[i] );

    fSigmabarSector[i]=sqrt(fSigmabarInnerSector[i] + fSigmabarOuterSector[i]*fRatioA);
    fSigmabarInnerSector[i]=sqrt(fSigmabarInnerSector[i]);
    fSigmabarOuterSector[i]=sqrt(fSigmabarOuterSector[i]);
  }
    
  // Did all our calculations work? fOuterPixels==0 could happen, however.
  return (fInnerPixels!=0);
}


void MSigmabar::Print(Option_t *) const
{
  *fLog << all;
  *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;
  *fLog << "Total number of inner pixels found is " << fInnerPixels << endl;
  *fLog << "Total number of outer pixels found is " << fOuterPixels << endl;
  *fLog << "Ratio of areas outer/inner pixels found is " << fRatioA << endl;
  *fLog << endl;
}