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

/* ======================================================================== *\
!
! *
! * 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 useful for
// the calibration of our telescope. This class store the information into the
// container of MMuonCalibPar. 
//
//  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;
//
//   mucalibcalc.SetMargin(60.);
//
//  You can retrieve the histogram (TH1F) using the function of GetHistPhi() 
// (also GetHistWid()). Therefore, you can draw the histogram such as
//
//   MParList  plist;
//   MMuonCalibPar muparcalib;
//   plist.AddToList(&muparcalib);.
//   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;
// // 
// //   mucalibcalc.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;
//  
//   mucalibcalc.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.
// 
//
//  Input Containers:
//   MGeomCam
//   MSignalCam
//   MMuonSearchPar
//
//  Output Containers:
//   MMuonCalibPar
//
//////////////////////////////////////////////////////////////////////////////

#include "MMuonCalibParCalc.h"

#include <TH1.h>
#include <TF1.h>
#include <TMinuit.h>

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

#include "MParList.h"

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

#include "MSignalCam.h"

#include "MMuonCalibPar.h"
#include "MMuonSetup.h"
#include "MMuonSearchPar.h"
#include "MHSingleMuon.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)
//    : fEnableImpactCalc(kFALSE)
{
    fName  = name  ? name  : gsDefName.Data();
    fTitle = title ? title : gsDefTitle.Data();
}

// -------------------------------------------------------------------------
//
Int_t MMuonCalibParCalc::PreProcess(MParList *pList)
{
    fGeomCam = (MGeomCam*)pList->FindObject("MGeomCam");
    if (!fGeomCam)
    {
        *fLog << err << "MGeomCam not found... abort." << endl;
        return kFALSE;
    }

    fHist = (MHSingleMuon*)pList->FindObject("MHSingleMuon");
    if (!fHist)
    {
        *fLog << err << "MHSingleMuon not found... abort." << endl;
        return kFALSE;
    }

    fMuonSetup = (MMuonSetup*)pList->FindObject("MMuonSetup");
    if (!fMuonSetup)
    {
        *fLog << err << "MMuonSetup not found... abort." << endl;
        return kFALSE;
    }

    fMuonCalibPar = (MMuonCalibPar*)pList->FindCreateObj("MMuonCalibPar");
    if (!fMuonCalibPar)
        return kFALSE;

    fMuonSearchPar = (MMuonSearchPar*)pList->FindCreateObj("MMuonSearchPar");
    if (!fMuonSearchPar)
        return kFALSE;

    return kTRUE;
}

// -------------------------------------------------------------------------
//
Int_t MMuonCalibParCalc::Process()
{
    // Calculation of Arc Phi etc...
    const Float_t thresphi   = fMuonSetup->GetThresholdArcPhi()  /fGeomCam->GetConvMm2Deg();
    const Float_t threswidth = fMuonSetup->GetThresholdArcWidth()/fGeomCam->GetConvMm2Deg();

    Double_t peakphi, arcphi;
    Double_t width, chi;

    if (!fHist->CalcPhi(thresphi, peakphi, arcphi))
        return kTRUE;

    if (!fHist->CalcWidth(threswidth, width, chi))
        return kTRUE;

    // Get Muon Size
    fMuonCalibPar->SetMuonSize(fHist->GetHistPhi().Integral());

    // Calculate Arc Length
    fMuonCalibPar->SetPeakPhi(peakphi);
    fMuonCalibPar->SetArcPhi(arcphi);

    const Float_t conv = TMath::RadToDeg()*fGeomCam->GetConvMm2Deg();
    fMuonCalibPar->SetArcLength(fMuonCalibPar->GetArcPhi()*fMuonSearchPar->GetRadius()*conv);

    // Calculation of Arc Width etc...
    fMuonCalibPar->SetChiArcWidth(chi);
    fMuonCalibPar->SetArcWidth(width);

    // Check if this is a 'Write-Out' candidate
    if (fMuonCalibPar->GetArcPhi()>160    && fMuonSearchPar->GetRadius()<400 &&
        fMuonSearchPar->GetDeviation()<50 && fMuonSearchPar->GetRadius()>170)
        fMuonCalibPar->SetReadyToSave();

    return kTRUE;
}