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

/////////////////////////////////////////////////////////////////////////////
//
// MMuonCalibPar
//
// Storage Container for muon
//
//  This class holds some information for a calibraion using muons. Muons 
// are identified by using the class of the MMuonSearchParCalc. You can fill 
// these information by using the MMuonCalibParCalc. See also these class
// manuals.
//
//  
//  Input Containers:
//   [MGeomCam]
//   [MCerPhotEvt]
//   [MMuonSearchPar]
//
/////////////////////////////////////////////////////////////////////////////
#include "MMuonCalibPar.h"

#include <fstream>

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

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

using namespace std;

ClassImp(MMuonCalibPar);

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

    fHistPhi.SetName("HistPhi");
    fHistPhi.SetTitle("HistPhi");
    fHistPhi.SetXTitle("phi [deg.]");
    fHistPhi.SetYTitle("sum of ADC");
    fHistPhi.SetDirectory(NULL);
    fHistPhi.SetFillStyle(4000);
    fHistPhi.UseCurrentStyle();

    fHistWid.SetName("HistWid");
    fHistWid.SetTitle("HistWid");
    fHistWid.SetXTitle("distance from the ring center [deg.]");
    fHistWid.SetYTitle("sum of ADC");
    fHistWid.SetDirectory(NULL);
    fHistWid.SetFillStyle(4000);
    fHistWid.UseCurrentStyle();

    fKeepHist = kFALSE;         // normally histgram is not kept
    fEnableImpactCalc = kFALSE; // By default the calculation of impact parameter is skipped.
    fDisablePreCuts = kFALSE;   // By default the pre cuts will be applied.
    fUseCleanForWid = kFALSE;   // By default all the pixels will be used for the histogram of arc width.

    fMargin = 60.;  // in mm
    fArcPhiThres  = 100.;
    fArcWidThres  = 100.;
    fArcPhiBinNum = 20;
    fArcPhiHistStartVal = -180.; // deg.
    fArcPhiHistEndVal   = 180.;  // deg.
    fArcWidBinNum = 28;
    fArcWidHistStartVal = 0.3;   // deg.
    fArcWidHistEndVal   = 1.7;   // deg.
}

// --------------------------------------------------------------------------
//
void MMuonCalibPar::Reset()
{
  fArcLen      = -1.;
  fArcPhi      =  0.;
  fArcWid      = -1.;
  fChiArcPhi   = -1.;
  fChiArcWid   = -1.;
  fMuonSize    =  0.;
  fEstImpact   = -1.;
  fUseUnmap    = kFALSE;
  fPeakPhi     =  0.;
}

// --------------------------------------------------------------------------
// 
//  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.
// You can use the image after the cleaning by using the function of UseCleanForWid().
// See the manual of MMuonCalibParCalc.
// 
void MMuonCalibPar::FillHist
( const MGeomCam &geom, const MCerPhotEvt &evt, 
  const MMuonSearchPar &musearch)
{
  // preparation for a histgram
  MBinning binsphi;
  binsphi.SetEdges(fArcPhiBinNum, fArcPhiHistStartVal, fArcPhiHistEndVal);
  binsphi.Apply(fHistPhi);

  MBinning binswid;
  binswid.SetEdges(fArcWidBinNum, fArcWidHistStartVal, fArcWidHistEndVal);
  binswid.Apply(fHistWid);

 const Int_t entries = evt.GetNumPixels();

  // the position of the center of a muon ring
  const Float_t cenx = musearch.GetCenX();
  const Float_t ceny = musearch.GetCenY();
  
  for (Int_t i=0; i<entries; i++ )
    {
      MCerPhotPix &pix = (evt)[i];
      
      const MGeomPix &gpix = (geom)[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 < musearch.GetRad() + fMargin && dist > musearch.GetRad() - fMargin)
	{
	  // 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())
	    {
	      fUseUnmap = kTRUE;
	      continue;
	    }
	  fHistPhi.Fill(ang*kRad2Deg, pix.GetNumPhotons());
	  fMuonSize += pix.GetNumPhotons();
	}

      if(pix.GetPixId()>397)
	continue;  // use only the inner pixles
      
      if(fUseCleanForWid)
	{
	  if(!pix.IsPixelUsed())
	    continue;
	}

      fHistWid.Fill(dist*geom.GetConvMm2Deg(), pix.GetNumPhotons());
    }


  // 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<fArcPhiBinNum+1; i++)
      {
	Float_t rougherr  = TMath::Sqrt(TMath::Abs(fHistPhi.GetBinContent(i))*ADC2PhEl)/ADC2PhEl*Ffactor;
	{
	  fHistPhi.SetBinError(i, rougherr);
	}
      }
    for(Int_t i=0; i<fArcWidBinNum+1; i++)
      {
	Float_t rougherr = TMath::Sqrt(TMath::Abs(fHistWid.GetBinContent(i))*ADC2PhEl)/ADC2PhEl*Ffactor;
	{
	  fHistWid.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 MMuonCalibPar::CalcPhi
(const MGeomCam &geom, const MCerPhotEvt &evt,
 const MMuonSearchPar &musearch)
{
  Float_t convbin2val = (fArcPhiHistEndVal-fArcPhiHistStartVal)/
    (Float_t)fArcPhiBinNum;

    // adjust the peak to 0
    Float_t maxval = 0.;
    Int_t   maxbin = 0;
    maxval = fHistPhi.GetMaximum();
    maxbin = fHistPhi.GetMaximumBin();
    fPeakPhi = 180.-(Float_t)(maxbin-1)*convbin2val; 
    TArrayD tmp;
    tmp.Set(fArcPhiBinNum+1);
    for(Int_t i=1; i<fArcPhiBinNum+1; i++)
      {
	tmp[i] = fHistPhi.GetBinContent(i);
      }
    for(Int_t i=1; i<fArcPhiBinNum+1; i++)
      {
	Int_t id;
	id = i + (maxbin-(Int_t)((Float_t)fArcPhiBinNum/2.)-1);
	if(id>fArcPhiBinNum)
	  {
	    id-=(fArcPhiBinNum);
	  }
	if(id<=0)
	  {
	    id+=(fArcPhiBinNum);
	  }
	fHistPhi.SetBinContent(i,tmp[id]);
      }
    maxbin = (Int_t)((Float_t)fArcPhiBinNum/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<fArcPhiBinNum+1; i++)
    {
      Float_t content = fHistPhi.GetBinContent(i);
      Float_t content_pre = fHistPhi.GetBinContent(i-1);
      
      if(content > fArcPhiThres && content_pre < fArcPhiThres)
	{
	  startfitval = (Float_t)(i-1)*convbin2val+fArcPhiHistStartVal;
	}
      if(i==maxbin)
	IsInMaxim = kTRUE;
      
      if(content < fArcPhiThres && IsInMaxim == kTRUE)
	{
	  endfitval = (Float_t)(i-1)*convbin2val+fArcPhiHistStartVal;
	  break;
	}
      endfitval = fArcPhiHistEndVal;
    }
  
  effbinnum = (Int_t)((endfitval-startfitval)/convbin2val);
  
  fArcPhi = effbinnum*convbin2val;
  fArcLen = fArcPhi*3.1415926/180.*musearch.GetRad()*geom.GetConvMm2Deg();
  
  if(fEnableImpactCalc)
      CalcImpact(geom, musearch, 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 MMuonCalibPar::CalcImpact
( const MGeomCam &geom, const MMuonSearchPar &musearch, 
  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 = musearch.GetRad()*geom.GetConvMm2Deg()*3.1415926/180.;
  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.);
  
  if(fKeepHist)
    {
      fHistPhi.Fit("f1","RQEM");
    }
  else
    {
      fHistPhi.Fit("f1","RQEM0");
    }
  
  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 = musearch.GetRad()*geom.GetConvMm2Deg()*3.1415926/180.;
  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.);
  
  if(fKeepHist)
    {
      fHistPhi.Fit("f2","RQEM+");
    }
  else
    {
      fHistPhi.Fit("f2","RQEM0+");
    }
  
  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)
    {
      fEstImpact = -1.;
    }
  if(chi2 > chi1)
    {
      fEstImpact = estip1;
      fChiArcPhi = chi1;
    }
  else
    {
      fEstImpact = estip2;
      fChiArcPhi = 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 MMuonCalibPar::CalcWid
(const MGeomCam &geom, const MCerPhotEvt &evt,
 const MMuonSearchPar &musearch)
{
  Float_t convbin2val = (fArcWidHistEndVal - fArcWidHistStartVal)
    /(Float_t)fArcWidBinNum;

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

	if(content > fArcWidThres)
	  effbinnum++;

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

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

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

    if(fKeepHist)
      {
	fHistWid.Fit("f1","QR","",startfitval,endfitval);
      }
    else
      {
	fHistWid.Fit("f1","QR0","",startfitval,endfitval);
      }
    
    if(effbinnum>3)
      fChiArcWid = 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 MMuonCalibPar::Calc
(const MGeomCam &geom, const MCerPhotEvt &evt, 
 const MMuonSearchPar &musearch)
{
  // sanity check
  if(evt.GetNumPixels() < 3)
    return kCONTINUE;

  // if an event does not seem to be like muon, it will be skipped.
  if(musearch.IsNoMuon()) 
    return kCONTINUE; 

  // Pre Cuts 1
  if(!fDisablePreCuts)
    {
      if(musearch.GetRad() < 180. || musearch.GetRad() > 400.)
	return kCONTINUE;
      if(musearch.GetDev() > 50.)
	return kCONTINUE;
    }

  Reset();

  FillHist(geom,evt,musearch);

  CalcPhi(geom,evt,musearch);

  // Pre Cuts 2
  if(!fDisablePreCuts)
    {
      if(fMuonSize < 2000. || fArcPhi < 150.)
	return kCONTINUE;
    }

  fArcWid = CalcWid(geom,evt,musearch); 
  
  // Pre Cuts 3
  if(!fDisablePreCuts)
    {
      if(fArcWid < 0.16)
	return kCONTINUE;
    }

  SetReadyToSave();
  
  if(!fKeepHist){
    fHistPhi.Reset();
    fHistWid.Reset();
  }

  return kTRUE;
} 

void MMuonCalibPar::Print(Option_t *) const
{
    *fLog << all;
    *fLog << "Muon Parameters (" << GetName() << ")"    << endl;
    *fLog << " - Arc Len.    [deg.]  = " << fArcLen     << endl;
    *fLog << " - Arc Phi     [deg.]  = " << fArcPhi     << endl;
    *fLog << " - Arc Wid.    [deg.]  = " << fArcWid     << endl;
    *fLog << " - Chi Arc Phi [x2/ndf]= " << fChiArcPhi  << endl;
    *fLog << " - Chi Arc Wid [x2/ndf]= " << fChiArcWid  << endl;
    *fLog << " - Est. I. P.  [m]     = " << fEstImpact  << endl;
    *fLog << " - Size of muon        = " << fMuonSize   << endl;
    *fLog << " - Peak Phi    [deg.]  = " << fPeakPhi    << endl;
    *fLog << " - UseUnmap            = " << fUseUnmap   << endl;
}