/* ======================================================================== *\
!
! *
! * 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]
//   [MCerPhotEvt]
//   [MMuonSearchPar]
//
//  Output Containers:
//   [MMuonCalibPar]
//
//////////////////////////////////////////////////////////////////////////////

#include "MMuonCalibParCalc.h"

#include <fstream>

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

#include "MParList.h"

#include "MGeomCam.h"
#include "MGeomPix.h"
#include "MSrcPosCam.h"
#include "MCerPhotEvt.h"
#include "MMuonSearchPar.h"
#include "MMuonCalibPar.h"
#include "MLog.h"
#include "MLogManip.h"
#include "MBinning.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)
  : fNameCerPhot("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.;

    fMargin = 60.;
    fArcPhiThres = 100.;
    fArcWidthThres = 100.;

    fEnableImpactCalc = kFALSE; // By default the calculation of impact parameter is skipped.
    fDisablePreCuts =   kFALSE; // By default the pre cuts will be applied.
}

// -------------------------------------------------------------------------
//
Int_t MMuonCalibParCalc::PreProcess(MParList *pList)
{
  fCerPhotEvt = (MCerPhotEvt*)pList->FindObject(AddSerialNumber(fNameCerPhot), "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;
}

// --------------------------------------------------------------------------
// 
//  This function fill the histograms in order to get muon parameters.
// For the evaluation of the Arc Width, we use only the signals in the inner 
// part. 
// 
void MMuonCalibParCalc::FillHist()
{
  Float_t MuonSize = 0.;

  Int_t binnumphi = fMuonCalibPar->fArcPhiBinNum;
  Int_t binnumwid = fMuonCalibPar->fArcWidthBinNum;

  // preparation for a histgram
  MBinning binsphi;
  binsphi.SetEdges(binnumphi, 
		   fMuonCalibPar->fArcPhiHistStartVal, 
		   fMuonCalibPar->fArcPhiHistEndVal);
  binsphi.Apply(*(fMuonCalibPar->fHistPhi));

  MBinning binswid;
  binswid.SetEdges(binnumwid, 
		   fMuonCalibPar->fArcWidthHistStartVal, 
		   fMuonCalibPar->fArcWidthHistEndVal);
  binswid.Apply(*(fMuonCalibPar->fHistWidth));

  const Int_t entries = (*fCerPhotEvt).GetNumPixels();

  // the position of the center of a muon ring
  const Float_t cenx = (*fMuonSearchPar).GetCenterX();
  const Float_t ceny = (*fMuonSearchPar).GetCenterY();
  
  for (Int_t i=0; i<entries; i++ )
    {
      MCerPhotPix &pix = (*fCerPhotEvt)[i];
      
      const MGeomPix &gpix = (*fGeomCam)[pix.GetPixId()];
      
      const Float_t dx = gpix.GetX() - cenx;
      const Float_t dy = gpix.GetY() - ceny;
      
     const Float_t dist = TMath::Sqrt(dx*dx+dy*dy);
      
      Float_t ang = TMath::ACos(dx/dist);
      if(dy>0)
	ang *= -1.0;
      
      // if the signal is not near the estimated circle, it is ignored.
      if(dist < (*fMuonSearchPar).GetRadius() + fMuonCalibPar->GetMargin() 
	 && dist > (*fMuonSearchPar).GetRadius() - fMuonCalibPar->GetMargin())
	{
	  // check whether ummapped pixel is used or not.
	  // if it is so, ingnore the pixel information since the pixels totally deteriorate the muon information.
	  if(pix.IsPixelUnmapped())
	    {
	      fMuonCalibPar->SetUseUnmap();
	      continue;
	    }
	  fMuonCalibPar->fHistPhi->Fill(ang*kRad2Deg, pix.GetNumPhotons());
	  MuonSize += pix.GetNumPhotons();
	}

      // use only the inner pixles. This is geometry dependent. This has to
      // be fixed!
      if(pix.GetPixId()>397)
	continue;  
      
      fMuonCalibPar->fHistWidth
	->Fill(dist*(*fGeomCam).GetConvMm2Deg(), pix.GetNumPhotons());
    }

  fMuonCalibPar->SetMuonSize(MuonSize);

  // error estimation (temporaly)
  //  The error is estimated from the signal. In order to do so, we have to
  // once convert the signal from ADC to photo-electron. Then we can get
  // the fluctuation such as F-factor*sqrt(phe). 
  //  Up to now, the error of pedestal is not taken into accout. This is not 
  // of course correct. We will include this soon.
    Double_t ADC2PhEl = 0.14;
    Double_t Ffactor = 1.4;
    for(Int_t i=0; i<binnumphi+1; i++)
      {
	Float_t rougherr  = TMath::Sqrt(TMath::Abs(fMuonCalibPar->fHistPhi->GetBinContent(i))*ADC2PhEl)/ADC2PhEl*Ffactor;
	{
	  fMuonCalibPar->fHistPhi->SetBinError(i, rougherr);
	}
      }
    for(Int_t i=0; i<binnumwid+1; i++)
      {
	Float_t rougherr = TMath::Sqrt(TMath::Abs(fMuonCalibPar->fHistWidth->GetBinContent(i))*ADC2PhEl)/ADC2PhEl*Ffactor;
	{
	  fMuonCalibPar->fHistWidth->SetBinError(i, rougherr);
	}
      }
}

// --------------------------------------------------------------------------
//
//  Photon distribution along the estimated circle is fitted with theoritical
// function in order to get some more information such as Arc Phi and Arc 
// Length.
//
void MMuonCalibParCalc::CalcPhi()
{
  Float_t thres = fMuonCalibPar->GetArcPhiThres();
  Float_t startval = fMuonCalibPar->fArcPhiHistStartVal;
  Float_t endval = fMuonCalibPar->fArcPhiHistEndVal;
  Int_t   binnum = fMuonCalibPar->fArcPhiBinNum;

  Float_t convbin2val = (endval-startval)/(Float_t)binnum;

    // adjust the peak to 0
    Float_t maxval = 0.;
    Int_t   maxbin = 0;
    maxval = fMuonCalibPar->fHistPhi->GetMaximum();
    maxbin = fMuonCalibPar->fHistPhi->GetMaximumBin();
    fMuonCalibPar->SetPeakPhi(180.-(Float_t)(maxbin-1.)*convbin2val); 
    TArrayD tmp;
    tmp.Set(binnum+1);
    for(Int_t i=1; i<binnum+1; i++)
      {
	tmp[i] = fMuonCalibPar->fHistPhi->GetBinContent(i);
      }
    for(Int_t i=1; i<binnum+1; i++)
      {
	Int_t id;
	id = i + (maxbin-(Int_t)((Float_t)binnum/2.)-1);
	if(id>binnum)
	  {
	    id-=(binnum);
	  }
	if(id<=0)
	  {
	    id+=(binnum);
	  }
	fMuonCalibPar->fHistPhi->SetBinContent(i,tmp[id]);
      }
    maxbin = (Int_t)((Float_t)binnum/2.)+1;

  // Determination of fitting region
  // The threshold is fixed with 100 [photons or ADC] in a bin. Therefore,
  // if you change the bin number, YOU HAVE TO CHANGE THIS VALUE!!!
  Float_t startfitval = 0.;
  Float_t endfitval   = 0.;
  Bool_t  IsInMaxim = kFALSE;
  Int_t effbinnum = 0;
  for(Int_t i=1; i<binnum+1; i++)
    {
      Float_t content = fMuonCalibPar->fHistPhi->GetBinContent(i);
      Float_t content_pre = fMuonCalibPar->fHistPhi->GetBinContent(i-1);
      
      if(content > thres && content_pre < thres)
	{
	  startfitval = (Float_t)(i-1)*convbin2val+startval;
	}
      if(i==maxbin)
	IsInMaxim = kTRUE;
      
      if(content < thres && IsInMaxim == kTRUE)
	{
	  endfitval = (Float_t)(i-1)*convbin2val+startval;
	  break;
	}
      endfitval = endval;
    }
  
  effbinnum = (Int_t)((endfitval-startfitval)/convbin2val);

  fMuonCalibPar->SetArcPhi(endfitval-startfitval);

  fMuonCalibPar->SetArcLength( fMuonCalibPar->GetArcPhi()*TMath::DegToRad()*(*fMuonSearchPar).GetRadius()*(*fGeomCam).GetConvMm2Deg());
  
  if(fEnableImpactCalc)
      CalcImpact(effbinnum, startfitval, endfitval);
}

// --------------------------------------------------------------------------
//
//  An impact parameter is calculated by fitting the histogram of photon
// distribution along the circle with a theoritical model. 
// (See G. Vacanti et. al., Astroparticle Physics 2, 1994, 1-11. 
// The function (6) is used here.) 
//
//  By default this calculation is suppressed because this calculation is 
// very time consuming. If you want to calculate an impact parameter,
// you can call the function of EnableImpactCalc().
//
void MMuonCalibParCalc::CalcImpact
(Int_t effbinnum, Float_t startfitval, Float_t endfitval)
{
  // Fit the distribution with Vacanti function. The function is different
  // for the impact parameter of inside or outside of our reflector. 
  // Then two different functions are applied to the photon distribution,
  // and the one which give us smaller chisquare value is taken as a 
  // proper one.
  Double_t val1,err1,val2,err2;
  // impact parameter inside mirror radius (8.5m)
  TString func1;
  Float_t tmpval = (*fMuonSearchPar).GetRadius()*(*fGeomCam).GetConvMm2Deg()*TMath::DegToRad();
  tmpval = sin(2.*tmpval)*8.5;
  func1 += "[0]*";
  func1 += tmpval;
  func1 += "*(sqrt(1.-([1]/8.5)**2*sin((x-[2])*3.1415926/180.)**2)+([1]/8.5)*cos((x-[2])*3.1415926/180.))";
  
  TF1 f1("f1",func1,startfitval,endfitval);
  f1.SetParameters(2000,3,0);
  f1.SetParLimits(1,0,8.5);
  f1.SetParLimits(2,-180.,180.);
  
  fMuonCalibPar->fHistPhi->Fit("f1","RQEM");
  
  Float_t chi1 = -1;
  Float_t chi2 = -1;
  if(effbinnum>3)
    chi1 = f1.GetChisquare()/((Float_t)(effbinnum-3));
  
  gMinuit->GetParameter(1,val1,err1);  // get the estimated IP
  Float_t estip1 = val1;
  
  // impact parameter beyond mirror area (8.5m)
  TString func2;
  Float_t tmpval2 = (*fMuonSearchPar).GetRadius()*(*fGeomCam).GetConvMm2Deg()*TMath::DegToRad();
  tmpval2 = sin(2.*tmpval2)*8.5*2.;
  func2 += "[0]*";
  func2 += tmpval2;
  func2 += "*sqrt(1.-(([1]/8.5)*sin((x-[2])*3.1415926/180.))**2)";
  
  TF1 f2("f2",func2,startfitval,endfitval);
  f2.SetParameters(2000,20,0);
  f2.SetParLimits(1,8.5,300.);
  f2.SetParLimits(2,-180.,180.);
  
  fMuonCalibPar->fHistPhi->Fit("f2","RQEM+");
  
  if(effbinnum>3)
    chi2 = f2.GetChisquare()/((Float_t)(effbinnum-3));
  
  gMinuit->GetParameter(1,val2,err2);  // get the estimated IP
  Float_t estip2 = val2;
  
  if(effbinnum<=3)
    {
      fMuonCalibPar->SetEstImpact(-1.);
    }
  if(chi2 > chi1)
    {
      fMuonCalibPar->SetEstImpact(estip1);
      fMuonCalibPar->SetChiArcPhi(chi1);
    }
  else
    {
      fMuonCalibPar->SetEstImpact(estip2);
      fMuonCalibPar->SetChiArcPhi(chi2);
    }
}

// --------------------------------------------------------------------------
//
//  Photon distribution of distance from the center of estimated ring is 
// fitted in order to get some more information such as ARC WIDTH which 
// can represent to the PSF of our reflector.
//
Float_t MMuonCalibParCalc::CalcWidth()
{
  Float_t startval = fMuonCalibPar->fArcWidthHistStartVal;
  Float_t endval = fMuonCalibPar->fArcWidthHistEndVal;
  Int_t   binnum = fMuonCalibPar->fArcWidthBinNum;
  Float_t thres = fMuonCalibPar->GetArcWidthThres();

  Float_t convbin2val = (endval - startval)
    /(Float_t)binnum;

    // determination of fitting region
    Int_t maxbin = fMuonCalibPar->fHistWidth->GetMaximumBin();
    Float_t startfitval = 0.;
    Float_t endfitval   = 0.;
    Bool_t   IsInMaxim = kFALSE;
    Int_t    effbinnum = 0;
    for(Int_t i=1; i<binnum+1; i++)
      {
	Float_t content = fMuonCalibPar->fHistWidth->GetBinContent(i);
	Float_t content_pre = fMuonCalibPar->fHistWidth->GetBinContent(i-1);

	if(content > thres)
	  effbinnum++;

	if(content > thres && content_pre < thres)
	  {
	    startfitval = (Float_t)(i-4)*convbin2val + startval;
	    if(startfitval<0.) startfitval = 0.;
	  }
	if(i==maxbin)
	  IsInMaxim = kTRUE;

	if(content < thres && IsInMaxim == kTRUE)
	  {
	    endfitval = (Float_t)(i+2)*convbin2val + startval;
	    if(endfitval>180.) endfitval = 180.;
	    break;
	  }
	endfitval = endval;
      }
    effbinnum = (Int_t)((endfitval-startfitval)/convbin2val);

    TF1 f1("f1","gaus",startfitval,endfitval);

    fMuonCalibPar->fHistWidth->Fit("f1","QR","",startfitval,endfitval);
    
    if(effbinnum>3)
      fMuonCalibPar->SetChiArcWidth(f1.GetChisquare()/((Float_t)(effbinnum-3)));

    Double_t val,err;
    gMinuit->GetParameter(2,val,err); // get the sigma value

    return val;
}

// --------------------------------------------------------------------------
//
//  Calculation of muon parameters
//
Int_t MMuonCalibParCalc::Calc(const Float_t *cuts)
{
  // sanity check
  if((*fCerPhotEvt).GetNumPixels() < 3)
    return kCONTINUE;

  // If an event does not seem to be like muon, the calculation will be skipped.
  if((*fMuonSearchPar).IsNoMuon()) 
    return kCONTINUE; 

  // Pre Cuts 1
  if(!fDisablePreCuts)
    {
      if((*fMuonSearchPar).GetRadius() < cuts[0] || (*fMuonSearchPar).GetRadius() > cuts[1])
	{
	  (*fMuonSearchPar).SetNoMuon();
	  return kCONTINUE;
	}
      if((*fMuonSearchPar).GetDeviation() > cuts[2])
	{
	  (*fMuonSearchPar).SetNoMuon();
	  return kCONTINUE;
	}
    }

  // initialization
  (*fMuonCalibPar).Reset();

  // Fill signals to histograms
  FillHist();

  // Calculation of Arc Phi etc...
  CalcPhi();

  // Pre Cuts 2
  if(!fDisablePreCuts)
    {
      if(fMuonCalibPar->GetMuonSize() < cuts[3] 
	 || fMuonCalibPar->GetArcPhi() < cuts[4])
	{
	  (*fMuonSearchPar).SetNoMuon();
	  return kCONTINUE;
	}
    }

  // Calculation of Arc Width etc...
  fMuonCalibPar->SetArcWidth(CalcWidth()); 
  
  return kTRUE;
} 


// -------------------------------------------------------------------------
//
Int_t MMuonCalibParCalc::Process()
{
  fMuonCalibPar->SetMargin(fMargin);
  fMuonCalibPar->SetArcPhiThres(fArcPhiThres);
  fMuonCalibPar->SetArcWidthThres(fArcWidthThres);

  if(!Calc(fPreCuts))
    return kCONTINUE;

  return kTRUE;
}

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;
}

