source: trunk/MagicSoft/Mars/mhist/MHCalibrationBlindPixel.cc@ 2544

/* ======================================================================== *\
!
! *
! * 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): Markus Gaug 11/2003 <mailto:markus@ifae.es>
!
!   Copyright: MAGIC Software Development, 2000-2002
!
!
\* ======================================================================== */

//////////////////////////////////////////////////////////////////////////////
//                                                                          //
//  MHCalibrationBlindPixel                                                 //
//                                                                          //
//  Performs all the necessary fits to extract the mean number of photons   //
//              out of the derived light flux                               //
//                                                                          //
//////////////////////////////////////////////////////////////////////////////
#include "MHCalibrationBlindPixel.h"
#include "MHCalibrationConfig.h"
#include "MCalibrationFits.h"

#include <TStyle.h>
#include <TMath.h>
#include <TPad.h>

#include <TMinuit.h>
#include <TFitter.h>

#include <TF1.h>
#include <TH2.h>
#include <TCanvas.h>
#include <TPaveText.h>
#include <TRandom.h>

#include "MBinning.h"
#include "MParList.h"

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

ClassImp(MHCalibrationBlindPixel);

using namespace std;
// --------------------------------------------------------------------------
//
// Default Constructor. 
//
MHCalibrationBlindPixel::MHCalibrationBlindPixel(const char *name, const char *title)
    : fSinglePheFit(NULL), fTimeGausFit(NULL)
{

    fName  = name  ? name  : "MHCalibrationBlindPixel";
    fTitle = title ? title : "Fill the accumulated charges and times all Blind Pixel events and perform fits";

    // Create a large number of bins, later we will rebin
    fBPQfirst = 0;
    fBPQlast  = gkStartBlindPixelBinNr;
    fBPQnbins = gkStartBlindPixelBinNr;

    fHBPQ = new TH1I("HBPQ","Distribution of Summed FADC Slices",fBPQnbins,fBPQfirst,fBPQlast);
    fHBPQ->SetXTitle("Sum FADC Slices");
    fHBPQ->SetYTitle("Nr. of events");
    fHBPQ->Sumw2();

    fErrBPQfirst = 0.;
    fErrBPQlast  = gkStartBlindPixelBinNr;
    fErrBPQnbins = gkStartBlindPixelBinNr;

    fHBPErrQ = new TH1F("HBPErrQ","Distribution of Variances of Summed FADC Slices",
                        fErrBPQnbins,fErrBPQfirst,fErrBPQlast);
    fHBPErrQ->SetXTitle("Variance Summed FADC Slices");
    fHBPErrQ->SetYTitle("Nr. of events");
    fHBPErrQ->Sumw2();

    Axis_t tfirst = -0.5;
    Axis_t tlast  = 15.5;
    Int_t nbins   = 16;

    fHBPT = new TH1I("HBPT","Distribution of Mean Arrival Times",nbins,tfirst,tlast);
    fHBPT->SetXTitle("Mean Arrival Times [FADC slice nr]");
    fHBPT->SetYTitle("Nr. of events");
    fHBPT->Sumw2();

    // We define a reasonable number and later enlarge it if necessary
    nbins = 20000;
    Axis_t nfirst = -0.5;
    Axis_t nlast  = (Axis_t)nbins - 0.5;

    fHBPQvsN = new TH1I("HBPQvsN","Sum of Charges vs. Event Number",nbins,nfirst,nlast);
    fHBPQvsN->SetXTitle("Event Nr.");
    fHBPQvsN->SetYTitle("Sum of FADC slices");

    fgSinglePheFitFunc = &gfKto8;
    fgSinglePheFitNPar = 5;
}

MHCalibrationBlindPixel::~MHCalibrationBlindPixel()
{

  delete fHBPQ;
  delete fHBPT;
  delete fHBPErrQ;
  
  if (fSinglePheFit)
    delete fSinglePheFit;
  if (fTimeGausFit)
    delete fTimeGausFit;
}



void MHCalibrationBlindPixel::ResetBin(Int_t i)
{
    fHBPQ->SetBinContent (i, 1.e-20);
    fHBPErrQ->SetBinContent  (i, 1.e-20);
    fHBPT->SetBinContent(i, 1.e-20);
}


// -------------------------------------------------------------------------
//
// Draw a legend with the fit results
//
void MHCalibrationBlindPixel::DrawLegend()
{

  fFitLegend = new TPaveText(0.05,0.05,0.95,0.95);

  if (fFitOK) 
      fFitLegend->SetFillColor(80);
  else
      fFitLegend->SetFillColor(2);    
  
  fFitLegend->SetLabel("Results of the single PhE Fit (to k=6):");
  fFitLegend->SetTextSize(0.05);

  char line1[32];
  sprintf(line1,"Mean: #lambda = %2.2f #pm %2.2f",GetLambda(),GetLambdaErr());
  fFitLegend->AddText(line1);

  char line2[32];
  sprintf(line2,"Pedestal: #mu_{0} = %2.2f #pm %2.2f",GetMu0(),GetMu0Err());
  fFitLegend->AddText(line2);

  char line3[32];
  sprintf(line3,"Width Pedestal: #sigma_{0} = %2.2f #pm %2.2f",GetSigma0(),GetSigma0Err());
  fFitLegend->AddText(line3);

  char line4[32];
  sprintf(line4,"1^{st} Phe-peak: #mu_{1} = %2.2f #pm %2.2f",GetMu1(),GetMu1Err());
  fFitLegend->AddText(line4);

  char line5[32];
  sprintf(line5,"Width 1^{st} Phe-peak: #sigma_{1} = %2.2f #pm %2.2f",GetSigma1(),GetSigma1Err());
  fFitLegend->AddText(line5);

  char line7[32];
  sprintf(line7,"#chi^{2} / N_{dof}: %4.2f / %3i",GetChiSquare(),GetNdf());
  fFitLegend->AddText(line7);

  char line8[32];
  sprintf(line8,"Probability: %4.2f ",GetProb());
  fFitLegend->AddText(line8);

  if (fFitOK)
    fFitLegend->AddText("Result of the Fit: OK");
  else
    fFitLegend->AddText("Result of the Fit: NOT OK");

  fFitLegend->SetBit(kCanDelete);
  fFitLegend->Draw();

}


// -------------------------------------------------------------------------
//
// Draw the histogram
//
void MHCalibrationBlindPixel::Draw(Option_t *opt) 
{

    gStyle->SetOptFit(0);
    gStyle->SetOptStat(1111111);

    TCanvas *c = MakeDefCanvas(this,550,700);

    c->Divide(2,2);

    gROOT->SetSelectedPad(NULL);

    c->cd(1);
    gPad->SetLogy(1);
    gPad->SetTicks();

    fHBPQ->DrawCopy(opt);
    
    if (fSinglePheFit)
      {
        if (fFitOK)
          fSinglePheFit->SetLineColor(kGreen);          
        else
          fSinglePheFit->SetLineColor(kRed);

        fSinglePheFit->DrawCopy("same");
        c->Modified();
        c->Update();
      }

    c->cd(2);
    DrawLegend();
    c->Modified();
    c->Update();

    c->cd(3);
    gPad->SetLogy(1);
    gPad->SetBorderMode(0);
    fHBPT->DrawCopy(opt);

    if (fHBPT->GetFunction("GausTime"))
      {
        TF1 *tfit = fHBPT->GetFunction("GausTime");
        if (tfit->GetProb() < 0.01)
          tfit->SetLineColor(kRed);
        else
          tfit->SetLineColor(kGreen);

        tfit->DrawCopy("same");
        c->Modified();
        c->Update();
      }
    
    c->cd(4);

    fHBPQvsN->DrawCopy(opt);

    c->Modified();
    c->Update();
}



Bool_t MHCalibrationBlindPixel::SimulateSinglePhe(Double_t lambda, Double_t mu0, Double_t mu1, Double_t sigma0, Double_t sigma1) 
{
  gRandom->SetSeed();

  if (fHBPQ->GetEntries() != 0)
    {
      *fLog << err << "Histogram " << fHBPQ->GetTitle() << " is already filled. " << endl;
      *fLog << err << "Create new class MHCalibrationBlindPixel for simulation! " << endl;
      return kFALSE;
    }
  
  TF1 *simulateSinglePhe = new TF1("simulateSinglePhe",fgSinglePheFitFunc,
                                   fBPQfirst,fBPQlast,fgSinglePheFitNPar);
  
  simulateSinglePhe->SetParameters(lambda,mu0,mu1,sigma0,sigma1);
  simulateSinglePhe->SetParNames("#lambda","#mu_0","#mu_1","#sigma_0","#sigma_1");
  simulateSinglePhe->SetNpx(fBPQnbins);  

  for (Int_t i=0;i<10000; i++) 
    {
      fHBPQ->Fill(simulateSinglePhe->GetRandom());
    }
  
  return kTRUE;
}


void MHCalibrationBlindPixel::ChangeFitFunc(BPFitFunc fitfunc, Int_t par)
{
  
  fgSinglePheFitFunc = fitfunc;
  fgSinglePheFitNPar = par;

}



Bool_t MHCalibrationBlindPixel::FitSinglePhe(Axis_t rmin, Axis_t rmax, Option_t *opt) 
{

  if (fSinglePheFit)
    return kFALSE;


  //
  // Get the fitting ranges
  //
  rmin = (rmin != 0.) ? rmin : fBPQfirst;
  rmax = (rmax != 0.) ? rmax : fBPQlast;

  //
  // First guesses for the fit (should be as close to reality as possible, 
  // otherwise the fit goes gaga because of high number of dimensions ...
  //
  const Stat_t   entries      = fHBPQ->GetSumOfWeights();
  const Double_t lambda_guess = 0.2;
  const Double_t mu_0_guess = fHBPQ->GetBinCenter(fHBPQ->GetMaximumBin());
  const Double_t si_0_guess = mu_0_guess/500.;
  const Double_t mu_1_guess = mu_0_guess + 2.;
  const Double_t si_1_guess = si_0_guess + si_0_guess;

  fSinglePheFit = new TF1("SinglePheFit",fgSinglePheFitFunc,rmin,rmax,fgSinglePheFitNPar+1);
  fSinglePheFit->SetParameters(lambda_guess,mu_0_guess,mu_1_guess,si_0_guess,si_1_guess,entries);
  fSinglePheFit->SetParNames("#lambda","#mu_0","#mu_1","#sigma_0","#sigma_1","area");
  fSinglePheFit->SetParLimits(0,0.,5.);
  fSinglePheFit->SetParLimits(1,rmin,rmax);
  fSinglePheFit->SetParLimits(2,rmin,rmax);
  fSinglePheFit->SetParLimits(3,1.0,rmax-rmin);
  fSinglePheFit->SetParLimits(4,1.7,rmax-rmin);
  fSinglePheFit->SetParLimits(5,0.,2.*entries);
  //
  // Normalize the histogram to facilitate faster fitting of the area
  // For speed reasons, FKto8 is normalized to Sqrt(2 pi).
  // Therefore also normalize the histogram to that value 
  //
  // ROOT gives us another nice example of user-unfriendly behavior:
  // Although the normalization of the function fSinglePhe and the 
  // Histogram fHBPQ agree (!!), the fit does not normalize correctly INTERNALLY
  // in the fitting procedure !!!
  // 
  // This has to do with the fact that the internal function histogramming 
  // uses 100 bins and does not adapt to the binning of the fitted histogram, unlike PAW does
  // (very important if you use Sumw2, see e.g. ROOTTALK: Mon May 26 1997 - 09:56:03 MEST)
  // 
  // So, WE have to adapt to that internal flaw of ROOT:
  //
  const Int_t  npx     = fSinglePheFit->GetNpx();
  const Int_t  bins    = fHBPQ->GetXaxis()->GetLast()-fHBPQ->GetXaxis()->GetFirst();
  //  fHBPQ->Scale(gkSq2Pi*(float)bins/npx/entries);

  // 
  // we need this, otherwise, ROOT does not calculate the area correctly
  // don't ask me why it does not behave correctly, it's one of the nasty
  // mysteries of ROOT which takes you a whole day to find out :-)
  //
  //  fSinglePheFit->SetNpx(fQnbins);  

  fHBPQ->Fit("SinglePheFit",opt);

  fLambda = fSinglePheFit->GetParameter(0);
  fMu0    = fSinglePheFit->GetParameter(1);
  fMu1    = fSinglePheFit->GetParameter(2);
  fSigma0 = fSinglePheFit->GetParameter(3);
  fSigma1 = fSinglePheFit->GetParameter(4);
  
  fLambdaErr = fSinglePheFit->GetParError(0);
  fMu0Err    = fSinglePheFit->GetParError(1);
  fMu1Err    = fSinglePheFit->GetParError(2);
  fSigma0Err = fSinglePheFit->GetParError(3);
  fSigma1Err = fSinglePheFit->GetParError(4);

  fProb      = fSinglePheFit->GetProb();
  fChisquare = fSinglePheFit->GetChisquare();
  fNdf       = fSinglePheFit->GetNDF();

  *fLog << "Results of the Blind Pixel Fit: " << endl;
  *fLog << "Chisquare: " << fChisquare << endl;
  *fLog << "DoF: " << fNdf << endl;
  *fLog << "Probability: " << fProb << endl;
  *fLog << "Integral: " << fSinglePheFit->Integral(rmin,rmax);

  //
  // The fit result is accepted under condition
  // The Probability is greater than gkProbLimit (default 0.01 == 99%)
  //
  if (fProb < gkProbLimit) 
    {
      *fLog << warn << "Prob: " << fProb << " is smaller than the allowed value: " << gkProbLimit << endl;
      fFitOK = kFALSE;
      return kFALSE;
    }
  
  
  fFitOK = kTRUE;
    
  return kTRUE;
}

 
void MHCalibrationBlindPixel::CutAllEdges()
{

  //
  // The number 100 is necessary because it is the internal binning
  // of ROOT functions. A call to SetNpx() does NOT help
  // If you find another solution which WORKS!!, please tell me!!
  //
  Int_t nbins = 100;

  *fLog << "New number of bins in HSinQ: " << CutEdges(fHBPQ,nbins) << endl;

  fBPQfirst = fHBPQ->GetBinLowEdge(fHBPQ->GetXaxis()->GetFirst());
  fBPQlast  = fHBPQ->GetBinLowEdge(fHBPQ->GetXaxis()->GetLast())+fHBPQ->GetBinWidth(0);
  fBPQnbins = nbins;

  *fLog << "New number of bins in HErrQ: " << CutEdges(fHBPErrQ,30) << endl;
  fErrBPQfirst = fHBPErrQ->GetBinLowEdge(fHBPErrQ->GetXaxis()->GetFirst());
  fErrBPQlast  = fHBPErrQ->GetBinLowEdge(fHBPErrQ->GetXaxis()->GetLast())+fHBPErrQ->GetBinWidth(0);
  fErrBPQnbins = nbins;

  CutEdges(fHBPQvsN,0);

}

Bool_t MHCalibrationBlindPixel::FitT(Axis_t rmin, Axis_t rmax, Option_t *opt) 
{
  
  if (fTimeGausFit)
    return kFALSE;

  rmin = (rmin != 0.) ? rmin : 0.;
  rmax = (rmax != 0.) ? rmax : 16.;

  const Stat_t   entries     = fHBPT->GetEntries();
  const Double_t mu_guess    = fHBPT->GetBinCenter(fHBPT->GetMaximumBin());
  const Double_t sigma_guess = (rmax - rmin)/2.;
  const Double_t area_guess  = entries/gkSq2Pi;

  fTimeGausFit = new TF1("GausTime","gaus",rmin,rmax);  
  fTimeGausFit->SetParameters(area_guess,mu_guess,sigma_guess);
  fTimeGausFit->SetParNames("Area","#mu","#sigma");
  fTimeGausFit->SetParLimits(0,0.,entries);
  fTimeGausFit->SetParLimits(1,rmin,rmax);
  fTimeGausFit->SetParLimits(2,0.,rmax-rmin);

  fHBPT->Fit("GausTime",opt);

  fMeanT     = fTimeGausFit->GetParameter(2);
  fSigmaT    = fTimeGausFit->GetParameter(3);
  fMeanTErr  = fTimeGausFit->GetParError(2);
  fSigmaTErr = fTimeGausFit->GetParError(3);

  Float_t prob = fTimeGausFit->GetProb();

  *fLog << "Results of the Times Fit: " << endl;
  *fLog << "Chisquare: "   << fTimeGausFit->GetChisquare() << endl;
  *fLog << "Ndf: "         << fTimeGausFit->GetNDF() << endl;
  *fLog << "Probability: " << prob << endl;

  if (prob < gkProbLimit) 
    {
      *fLog << warn << "Fit of the Arrival times failed ! " << endl;
      return kFALSE;
    }
  
  return kTRUE;

}