source: trunk/MagicSoft/Mars/mmuon/MMuonSearchPar.cc@ 5173

/* ======================================================================== *\
!
! *
! * 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): Keiichi Mase 09/2004 <mailto:mase@mppmu.mpg.de>
!              Markus Meyer 09/2004 <mailto:meyer@astro.uni-wuerzburg.de>
!
!   Copyright: MAGIC Software Development, 2000-2004
!
!
\* ======================================================================== */

/////////////////////////////////////////////////////////////////////////////
//
// MMuonSearchPar
//
// Storage Container for muon
//
//  This class is the container for muon parameters. Actually, in this class, 
// muons are searched by fitting the image with a circle. (This function will
// be called by using the class of MMuonSearchParCalc.) This container
// especially holds the information of the results of the search (fit of 
// a image by a circle).
//
//  In order to use further information of muons such as the width of arcs,
// the arc length along it, the muons size. Use the infomation stored in 
// MMuonCalibPar. The information will be available by using the task of 
// MMuonCalibParCalc.
// 
// 
// --- How to search muons --- 
//  (This information is a little bit technical. You can skip to read if you 
//   don't need the technical information.)
// 
// 1. A temporal center position of a circle is determined by using 
//  the Hillas parameters. Assumed that the center position will be on the
//  line which is perpendicular to the longitudinal image axis and the 
//  distance from the gravity center of the image to the center position of
//  a ring is approximately 1 deg. (corresponding to the Cherenkov angle.).
//  Therefore, we will have two candidates of the center positions.
// 2. Find the ring radius which gives the minimum RMS between the camera
//  images and the estimated circle.
// 3. Select one temporal position which gives smaller RMS as a true temporal
//  center position.
// 4. Changing the center position of a circle on the camera plane from the 
// determined temporal center position, find the position which gives the 
// minimum RMS of the fit.
//
//
// --- Remark ---
//  This method to search for muons is not fully optimized yet. However,
// it is good idea to use the temporal estimated center position from
// the Hillas parameters in order to reduce time to search. In addition,
// This method is faster than the MINUIT.
//
//
//  Input Containers:
//   [MGeomCam]
//   [MHillas]
//   [MCerPhotEvt]
//
/////////////////////////////////////////////////////////////////////////////
#include "MMuonSearchPar.h"

#include <fstream>

#include "MLog.h"
#include "MLogManip.h"
#include "MHillas.h"
#include "MGeomCam.h"
#include "MGeomPix.h"
#include "MCerPhotEvt.h"
#include "MCerPhotPix.h"
#include "MBinning.h"

using namespace std;

ClassImp(MMuonSearchPar);

// --------------------------------------------------------------------------
//
// Default constructor.
//
MMuonSearchPar::MMuonSearchPar(const char *name, const char *title)
{
  fName  = name  ? name  : "MMuonSearchPar";
  fTitle = title ? title : "Muon search parameters";
}

// --------------------------------------------------------------------------
//
void MMuonSearchPar::Reset()
{
  fRad  = -1.;
  fDev  = -1.;
  fCenX = 0.;
  fCenY = 0.;
  fNoMuon = kFALSE;
}

// --------------------------------------------------------------------------
//
//  Get the tempolary center of a ring from the Hillas parameters.
//  Two candidates of the position is returened.
//
void MMuonSearchPar::CalcTempCenter(const MHillas &hillas, 
      Float_t *xtmp1, Float_t *ytmp1, Float_t *xtmp2, Float_t *ytmp2)
{
  Float_t a,dx,dy;
  Float_t tmp_r = 300.;  // assume that the temporal cherenkov angle is 1 deg. (300 mm).

  a = tan(hillas.GetDelta());

  dx = a/sqrt(tmp_r+a*a)/3.;
  dy = -tmp_r/sqrt(1+a*a)/3.;

  *xtmp1 = hillas.GetMeanX() + dx;
  *ytmp1 = hillas.GetMeanY() + dy;
  *xtmp2 = hillas.GetMeanX() - dx;
  *ytmp2 = hillas.GetMeanY() - dy;
}

// --------------------------------------------------------------------------
//
//  This function gives you the ring radius fitted best to the camera image
//  and its RMS for the input position.
//
Bool_t MMuonSearchPar::CalcRadius(const MGeomCam &geom, const MCerPhotEvt &evt,
      Float_t x, Float_t y, Float_t *r, Float_t *sigma)
{
  Float_t mean_r=0., dev_r=0., sums=0., tmp=0.;

  const Int_t entries = evt.GetNumPixels();

  for (Int_t i=0; i<entries; i++ ){
    const MCerPhotPix &pix = evt[i];

    if (!pix.IsPixelUsed())
      continue;

    const MGeomPix &gpix = geom[pix.GetPixId()];

    tmp=TMath::Sqrt((gpix.GetX()-x)*(gpix.GetX()-x)
                    +(gpix.GetY()-y)*(gpix.GetY()-y));

    mean_r += pix.GetNumPhotons()*tmp;
    dev_r  += pix.GetNumPhotons()*tmp*tmp;
    sums   += pix.GetNumPhotons();
  }

  if(sums<1.E-10)
      return kFALSE;

  *r = mean_r/sums;

  if(dev_r/sums-(*r)*(*r)<1.E-10)
      return kFALSE;

  *sigma = TMath::Sqrt(dev_r/sums-(*r)*(*r)); 

  return kTRUE;
} 

// --------------------------------------------------------------------------
//
//  This function finds the center position of the circle which gives minimum 
// RMS of the fit, changing the center position of the circle.
//
void MMuonSearchPar::CalcMinimumDev
( const MGeomCam &geom, const MCerPhotEvt &evt, Float_t x, Float_t y, 
 Float_t xcog, Float_t ycog, Float_t sigma, Float_t *opt_rad, 
 Float_t *new_sigma, Float_t *newx, Float_t *newy )
{
  Float_t delta = 3.;  // 3 mm (1/10 of an inner pixel size) Step to move.
  Float_t rad_tmp,sig_tmp;
  Float_t r2;

  while(1) 
  {
      r2=(xcog-x)*(xcog-x)+(ycog-y)*(ycog-y);
      // Exit if the new estimated radius is above 2 deg. (600 mm).
      if(r2 > 360000.)
      { 
          *new_sigma=sigma;
          *opt_rad=rad_tmp;
          *newx=x; *newy=y;
          fNoMuon = kTRUE;
          break;
      }
      if(CalcRadius(geom,evt,x,y+delta,&rad_tmp,&sig_tmp))
      {
          if(sig_tmp<sigma)
          {
              sigma=sig_tmp; *opt_rad=rad_tmp;
              y=y+delta;
              continue;
          }
      }
      if(CalcRadius(geom,evt,x-delta,y,&rad_tmp,&sig_tmp))
      {
          if(sig_tmp<sigma)
          {
              sigma=sig_tmp; *opt_rad=rad_tmp;
              x=x-delta;
              continue;
          }
      }
      if(CalcRadius(geom,evt,x+delta,y,&rad_tmp,&sig_tmp))
      {
          if(sig_tmp<sigma)
          {
              sigma=sig_tmp; *opt_rad=rad_tmp;
              x=x+delta;
              continue;
          }
      }
      if(CalcRadius(geom,evt,x,y-delta,&rad_tmp,&sig_tmp))
      {
          if(sig_tmp<sigma)
          {
              sigma=sig_tmp; *opt_rad=rad_tmp;
              y=y-delta;
              continue;
          }
      }
      *new_sigma=sigma;
      *newx=x; *newy=y;
      break;
    }
}

// --------------------------------------------------------------------------
//
//  Calculation of muon parameters
//
void MMuonSearchPar::Calc
(const MGeomCam &geom, const MCerPhotEvt &evt, const MHillas &hillas)
{
  Reset();
  
  Float_t xtmp1,xtmp2,ytmp1,ytmp2;
  Float_t rad,dev,rad2,dev2;
  Float_t opt_rad,new_sigma,newx,newy;
    
  // gets temporaly center
  CalcTempCenter(hillas,&xtmp1,&ytmp1,&xtmp2,&ytmp2);
  
  // determine which position will be the true position. Here mainly
  // the curvature of a muon arc is relied on.
  CalcRadius(geom, evt, xtmp1,ytmp1,&rad,&dev);
  CalcRadius(geom, evt, xtmp2,ytmp2,&rad2,&dev2);
  if(dev2<dev){
    xtmp1=xtmp2; ytmp1=ytmp2; dev=dev2; rad=rad2;
  }
  
  // find the best fit.
  CalcMinimumDev(geom, evt, xtmp1,ytmp1,hillas.GetMeanX(),hillas.GetMeanY(), 
                 dev, &opt_rad, &new_sigma, &newx, &newy);

  fRad = opt_rad;
  fDev = new_sigma;
  
  fCenX = newx;
  fCenY = newy;
  
  SetReadyToSave();
} 

void MMuonSearchPar::Print(Option_t *) const
{
    *fLog << all;
    *fLog << "Muon Parameters (" << GetName() << ")" << endl;
    *fLog << " - Est. Rad.   [mm]  = " << fRad  << endl;
    *fLog << " - Dev.        [mm]  = " << fDev  << endl;
    *fLog << " - Cen. Pos. X [mm]  = " << fCenX << endl;
    *fLog << " - Cen. Pos. Y [mm]  = " << fCenY << endl;
}

void MMuonSearchPar::Print(const MGeomCam &geom, Option_t *) const
{
    *fLog << all;
    *fLog << "Muon Parameters (" << GetName() << ")" << endl;
    *fLog << " - Est. Rad.   [deg.]  = " << fRad*geom.GetConvMm2Deg()   << endl;
    *fLog << " - Dev.        [deg.]  = " << fDev*geom.GetConvMm2Deg()   << endl;
    *fLog << " - Cen. Pos. X [deg.]  = " << fCenX*geom.GetConvMm2Deg()  << endl;
    *fLog << " - Cen. Pos. Y [deg.]  = " << fCenY*geom.GetConvMm2Deg()  << endl;
}