source: trunk/MagicSoft/Mars/mmuon/MMuonCalibParCalc.cc@ 5210

/* ======================================================================== *\
!
! *
! * 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
!
!
\* ======================================================================== */

//////////////////////////////////////////////////////////////////////////////
//
// MMuonCalibParCalc
//
// Task to calculate the muon parameters
//
//  This class allows you to get more muon information especially usefule for
// the calibration of our telescope. This class store the information into the
// container of MMuonCalibPar. The actual calculation is done in the class of 
// MMuonCalibPar. Please see the manual.
//
//  In order to make this class work, we need the information of the arc
// center and the radius. Therefore, we need to use the task of 
// MMuonSearchParCalc.
//
//  You can use this class such as the followings;
//
//   MTaskList tlist;
//   MMuonSearchParCalc musearchcalc;
//   MMuonCalibParCalc mucalibcalc;
//   tlist.AddToList(&musearchcalc);
//   tlist.AddToList(&mucalibcalc);.
//
//  You may change the allowed region to estimate muon parameters such as 
// Muon SIZE and ARC LENGTH. The default value is 60 mm (0.2 deg.). If the
// estimated radius of the arc is 1.0 degree, we take the photons in the 
// radius range from 0.8 to 1.2 degrees. You can change this value such as
// the followings;
//
//   MParList  plist;
//   MMuonCalibPar muparcalib;
//   plist.AddToList(&muparcalib);.
//   muparcalib.SetMargin(60.);
//
//  You can retrieve the histogram (TH1F) using the function of GetHistPhi() 
// (also GetHistWid()). Therefore, you can draw the histogram such as
//
//   muparcalib.GetHistPhi().Draw();.
//
//  In order to use another information of muons such as the center position 
// of the estimated circle, the radius of the circle. Use the infomation 
// stored in MMuonSearchPar. 
//
// 
//  For the faster computation, by default, the calculation of impact
// parameter is suppressed. If you want to calculate the impact parameter
// from the muon image, you can use the function of EnableImpactCalc(),
// namely;
//
//   muparcalib.EnableImpactCalc();.
//
//  In addition, for the faster comutation, pre cuts to select the candidates
// of muons for the calibration is done. You can set the values using the
// function of SetPreCuts. This function takes 5 variables. They correspond
// to the cur for the Arc Radius (low and high), the deviation of the fit
// (high), the Muon Size (low) and Arc Phi (low). You can set them such as 
//
//   mucalibcalc.SetPreCuts(180., 400., 50., 2000., 180.);
// 
//  If you want to disable the pre cuts, you can disable it by using the 
// function of DisablePreCuts(), namely;
//  
//   muparcalib.DisablePreCuts();.
//
// 
// ### TODO ###
//  Up to now, in the histogram the error of the signal is estimated from
// the signal using a rough conversion factor and a F-factor and this values
// are global for all pixels. This is not the case for the real data. This
// value should be taken from some containers. In addition, the error of 
// the pedestal is not taken into accout. The error treatment should be
// done correctly.
// 
//  Because of this rough treatment of the error, it may be better to use
// the image after the cleaning for the evaluation of the Arc Width. There
// is a possibility to use the image after the cleaning for the Arc Width.
// You can use the function of UseCleanForWid(). However, this treatment
// should be temporal and if the problem of the error is fixed, this function
// will be removed.
//
//
//  Input Containers:
//   [MGeomCam]
//   [MCerPhotEvt]
//   [MMuonSearchPar]
//
//  Output Containers:
//   [MMuonCalibPar]
//
//////////////////////////////////////////////////////////////////////////////

#include "MMuonCalibParCalc.h"

#include <fstream>

#include "MParList.h"

#include "MGeomCam.h"
#include "MSrcPosCam.h"
#include "MCerPhotEvt.h"
#include "MMuonSearchPar.h"
#include "MMuonCalibPar.h"
#include "MLog.h"
#include "MLogManip.h"

using namespace std;

ClassImp(MMuonCalibParCalc);

static const TString gsDefName  = "MMuonCalibParCalc";
static const TString gsDefTitle = "Calculate new image parameters";

// -------------------------------------------------------------------------
//
// Default constructor. 
//
MMuonCalibParCalc::MMuonCalibParCalc(const char *name, const char *title)
  : fCerPhotName("MCerPhotEvt")
{
    fName  = name  ? name  : gsDefName.Data();
    fTitle = title ? title : gsDefTitle.Data();

    fPreCuts[0]     = 180.;
    fPreCuts[1]     = 400.;
    fPreCuts[2]     = 50.;
    fPreCuts[3]     = 2000.;
    fPreCuts[4]     = 150.;
}

// -------------------------------------------------------------------------
//
Int_t MMuonCalibParCalc::PreProcess(MParList *pList)
{
  fCerPhotEvt = (MCerPhotEvt*)pList->FindObject(AddSerialNumber(fCerPhotName), "MCerPhotEvt");
  if (!fCerPhotEvt)
    {
      *fLog << dbginf << "MCerPhotEvt not found... aborting." << endl;
      return kFALSE;
    }
  
  fGeomCam = (MGeomCam*)pList->FindObject("MGeomCam");
  if (!fGeomCam)
    {
      *fLog << dbginf << "MGeomCam (Camera Geometry) missing in Parameter List... aborting." << endl;
      return kFALSE;
    }
  
  fMuonCalibPar = (MMuonCalibPar*)pList->FindCreateObj("MMuonCalibPar", "MMuonCalibPar");
  if (!fMuonCalibPar)
    {
      *fLog << dbginf << "MMuonCalibPar missing in Parameter List... aborting." << endl;
      return kFALSE;
    }
  
  fMuonSearchPar = (MMuonSearchPar*)pList->FindCreateObj("MMuonSearchPar", "MMuonSearchPar");
  if (!fMuonSearchPar)
    {
      *fLog << dbginf << "MMuonSearchPar missing in Parameter List... aborting." << endl;
      return kFALSE;
    }
  
  return kTRUE;
}

// -------------------------------------------------------------------------
//
Int_t MMuonCalibParCalc::Process()
{

  if(!fMuonCalibPar->Calc(*fGeomCam, *fCerPhotEvt, *fMuonSearchPar, fPreCuts))
    return kCONTINUE;

  return kTRUE;
}

void MMuonCalibParCalc::SetMargin(Float_t margin)
{
  fMuonCalibPar->SetMargin(margin);
}

void MMuonCalibParCalc::EnableImpactCalc()
{
  fMuonCalibPar->EnableImpactCalc();
}

void MMuonCalibParCalc::DisablePreCuts()
{
  fMuonCalibPar->DisablePreCuts();
}

void MMuonCalibParCalc::SetPreCuts
(Float_t radcutlow, Float_t radcuthigh, Float_t devcuthigh,
 Float_t musizecutlow, Float_t arcphicutlow)
{
  fPreCuts[0] = radcutlow;
  fPreCuts[1] = radcuthigh;
  fPreCuts[2] = devcuthigh;
  fPreCuts[3] = musizecutlow;
  fPreCuts[4] = arcphicutlow;
}