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

/* ======================================================================== *\
!
! *
! * 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 10/2004 <mailto:mase@mppmu.mpg.de>
!              Markus Meyer 10/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, the calculation
// is done here. 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"

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()
{
  fRadius  = -1.;
  fDeviation  = -1.;
  fCenterX = 0.;
  fCenterY = 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 = TMath::Tan(hillas.GetDelta());

  dx = a/TMath::Sqrt(tmp_r+a*a)/3.;
  dy = -tmp_r/TMath::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::CalcMinimumDeviation
( 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.
  CalcMinimumDeviation(geom, evt, xtmp1,ytmp1,hillas.GetMeanX(),
                       hillas.GetMeanY(), dev, opt_rad, new_sigma, 
                       newx, newy);

  fRadius = opt_rad;
  fDeviation = new_sigma;
  
  fCenterX = newx;
  fCenterY = newy;
  
  SetReadyToSave();
} 

void MMuonSearchPar::Print(Option_t *) const
{
    *fLog << all;
    *fLog << "Muon Parameters (" << GetName() << ")" << endl;
    *fLog << " - Est. Radius   [mm]  = " << fRadius  << endl;
    *fLog << " - Deviation     [mm]  = " << fDeviation  << endl;
    *fLog << " - Center Pos. X [mm]  = " << fCenterX << endl;
    *fLog << " - Center Pos. Y [mm]  = " << fCenterY << endl;
}

void MMuonSearchPar::Print(const MGeomCam &geom, Option_t *) const
{
    *fLog << all;
    *fLog << "Muon Parameters (" << GetName() << ")" << endl;
    *fLog << " - Est. Radius   [deg.]  = " << fRadius*geom.GetConvMm2Deg()   << endl;
    *fLog << " - Deviation     [deg.]  = " << fDeviation*geom.GetConvMm2Deg()   << endl;
    *fLog << " - Center Pos. X [deg.]  = " << fCenterX*geom.GetConvMm2Deg()  << endl;
    *fLog << " - Center Pos. Y [deg.]  = " << fCenterY*geom.GetConvMm2Deg()  << endl;
}