//-----------------------------------------------------------------------
// Fit ADC spectra with landau convoluted with Gaussian as a function of
// the distance from the pmt for both left and right side
// then fit these data points do determine attenuation length
// use adchists_run%d.root to get the 2d histograms x-position(TDC) vs. ADCIntegral
//
//
//-----------------------------------------------------------------------
#include "TTree.h"
#include "TFile.h"
#include "TH1.h"
#include "TH2.h"
#include "TF1.h"
#include "TROOT.h"
#include "TStyle.h"
#include "TGraphErrors.h"
#include "TMultiGraph.h"
#include "TMath.h"
#include "TCanvas.h"
#include "TText.h"

#include <iostream>
#include <fstream>

using namespace std;

int DEBUG = 2;
Double_t SNRPeak, SNRFWHM;
int RunNumber;
float findstart( TH1D *h);

Double_t langaufun(Double_t *x, Double_t *par) {

   //Fit parameters:
   //par[0]=Width (scale) parameter of Landau density
   //par[1]=Most Probable (MP, location) parameter of Landau density
   //par[2]=Total area (integral -inf to inf, normalization constant)
   //par[3]=Width (sigma) of convoluted Gaussian function
   //
   //In the Landau distribution (represented by the CERNLIB approximation), 
   //the maximum is located at x=-0.22278298 with the location parameter=0.
   //This shift is corrected within this function, so that the actual
   //maximum is identical to the MP parameter.

      // Numeric constants
      Double_t invsq2pi = 0.3989422804014;   // (2 pi)^(-1/2)
      Double_t mpshift  = -0.22278298;       // Landau maximum location

      // Control constants
      Double_t np = 100.0;      // number of convolution steps
      Double_t sc =   5.0;      // convolution extends to +-sc Gaussian sigmas

      // Variables
      Double_t xx;
      Double_t mpc;
      Double_t fland;
      Double_t sum = 0.0;
      Double_t xlow,xupp;
      Double_t step;
      Double_t i;


      // MP shift correction
      mpc = par[1] - mpshift * par[0]; 

      // Range of convolution integral
      xlow = x[0] - sc * par[3];
      xupp = x[0] + sc * par[3];

      step = (xupp-xlow) / np;

      // Convolution integral of Landau and Gaussian by sum
      for(i=1.0; i<=np/2; i++) {
         xx = xlow + (i-.5) * step;
         fland = TMath::Landau(xx,mpc,par[0]) / par[0];
         sum += fland * TMath::Gaus(x[0],xx,par[3]);

         xx = xupp - (i-.5) * step;
         fland = TMath::Landau(xx,mpc,par[0]) / par[0];
         sum += fland * TMath::Gaus(x[0],xx,par[3]);
      }

      return (par[2] * step * sum * invsq2pi / par[3]);
}



TF1 *langaufit(TH1D *his, Double_t *fitrange, Double_t *startvalues, Double_t *parlimitslo, Double_t *parlimitshi, Double_t *fitparams, Double_t *fiterrors, Double_t *ChiSqr, Int_t *NDF)
{
   // Once again, here are the Landau * Gaussian parameters:
   //   par[0]=Width (scale) parameter of Landau density
   //   par[1]=Most Probable (MP, location) parameter of Landau density
   //   par[2]=Total area (integral -inf to inf, normalization constant)
   //   par[3]=Width (sigma) of convoluted Gaussian function
   //
   // Variables for langaufit call:
   //   his             histogram to fit
   //   fitrange[2]     lo and hi boundaries of fit range
   //   startvalues[4]  reasonable start values for the fit
   //   parlimitslo[4]  lower parameter limits
   //   parlimitshi[4]  upper parameter limits
   //   fitparams[4]    returns the final fit parameters
   //   fiterrors[4]    returns the final fit errors
   //   ChiSqr          returns the chi square
   //   NDF             returns ndf

   Int_t i;
   Char_t FunName[100];

   sprintf(FunName,"Fitfcn_%s",his->GetName());

   TF1 *ffitold = (TF1*)gROOT->GetListOfFunctions()->FindObject(FunName);
   if (ffitold) delete ffitold;

   TF1 *ffit = new TF1(FunName,langaufun,fitrange[0],fitrange[1],4);
   ffit->SetParameters(startvalues);
   ffit->SetParNames("Width","MP","Area","GSigma");

   for (i=0; i<4; i++) {
      ffit->SetParLimits(i, parlimitslo[i], parlimitshi[i]);
   }

   ffit->FixParameter(3,0.001);
   his->Fit(FunName,"RB0");   // fit within specified range, use ParLimits, do not plot

   ffit->GetParameters(fitparams);    // obtain fit parameters
   for (i=0; i<4; i++) {
      fiterrors[i] = ffit->GetParError(i);     // obtain fit parameter errors
   }
   ChiSqr[0] = ffit->GetChisquare();  // obtain chi^2
   NDF[0] = ffit->GetNDF();           // obtain ndf

   return (ffit);              // return fit function

}


Int_t langaupro(Double_t *params, Double_t &maxx, Double_t &FWHM) {

   // Seaches for the location (x value) at the maximum of the 
   // Landau-Gaussian convolute and its full width at half-maximum.
   //
   // The search is probably not very efficient, but it's a first try.

   Double_t p,x,fy,fxr,fxl;
   Double_t step;
   Double_t l,lold;
   Int_t i = 0;
   Int_t MAXCALLS = 10000;


   // Search for maximum

   p = params[1] - 0.1 * params[0];
   step = 0.05 * params[0];
   lold = -2.0;
   l    = -1.0;


   while ( (l != lold) && (i < MAXCALLS) ) {
      i++;

      lold = l;
      x = p + step;
      l = langaufun(&x,params);
 
      if (l < lold)
         step = -step/10;
 
      p += step;
   }

   if (i == MAXCALLS)
      return (-1);

   maxx = x;

   fy = l/2;


   // Search for right x location of fy

   p = maxx + params[0];
   step = params[0];
   lold = -2.0;
   l    = -1e300;
   i    = 0;


   while ( (l != lold) && (i < MAXCALLS) ) {
      i++;

      lold = l;
      x = p + step;
      l = TMath::Abs(langaufun(&x,params) - fy);
 
      if (l > lold)
         step = -step/10;
 
      p += step;
   }

   if (i == MAXCALLS)
      return (-2);

   fxr = x;


   // Search for left x location of fy

   p = maxx - 0.5 * params[0];
   step = -params[0];
   lold = -2.0;
   l    = -1e300;
   i    = 0;

   while ( (l != lold) && (i < MAXCALLS) ) {
      i++;

      lold = l;
      x = p + step;
      l = TMath::Abs(langaufun(&x,params) - fy);
 
      if (l > lold)
         step = -step/10;
 
      p += step;
   }

   if (i == MAXCALLS)
      return (-3);


   fxl = x;

   FWHM = fxr - fxl;
   return (0);
}

TF1*  langaus(TH1D *hist) {

   // Fitting histogram
   printf("Fitting...\n");

   // Setting fit range and start values
   Double_t FitRange[2];
   Double_t StartValues[4], LimitLow[4], LimitHigh[4], fp[4], fpe[4];
   
   float x = findstart(hist);
   cout<<"start of signal around "<<x<<endl;

   hist->GetXaxis()->SetRangeUser(1500.,20000.);
   int maxbin = hist->GetMaximumBin();
   double max  = hist->GetBinCenter(maxbin);
   FitRange[0] = max-1000.;
   FitRange[1] = max+3000.;

   if (FitRange[1]>23000.){
     FitRange[1]=23000.;
   }

   cout<<"Fitrange:  "<<FitRange[0]<<"  "<<FitRange[1]<<endl;

   StartValues[0] = 100; //Width
   StartValues[1] = max; //MP
   StartValues[2] = hist->GetEntries(); //Area
   StartValues[3] = 10.0; //GSigma

   LimitHigh[0] = 1000.0;
   LimitLow[0]  = 100.;
   LimitHigh[1] = FitRange[1];
   LimitLow[1]  = FitRange[0];
   LimitHigh[2] = hist->GetEntries()*500.;
   LimitLow[2]  = hist->GetEntries()*0.01;
   LimitHigh[3] = 3000.0;
   LimitLow[3]  = 0.00001;

   Double_t chisqr;
   Int_t    ndf;
   TF1 *fitfunc = langaufit(hist,FitRange,StartValues,
			    LimitLow,LimitHigh,fp,fpe,&chisqr,&ndf);
   
   double w = fitfunc->GetParameter(0);
   double mp = fitfunc->GetParameter(1);
   FitRange[0] = mp-w*2.5;
   FitRange[1] = mp+w*12.;
   fitfunc = langaufit(hist,FitRange,StartValues,
		       LimitLow,LimitHigh,fp,fpe,&chisqr,&ndf);
   
   FitRange[0] = mp-w*2.5;
   FitRange[1] = mp+w*12.;
   fitfunc = langaufit(hist,FitRange,StartValues,
		       LimitLow,LimitHigh,fp,fpe,&chisqr,&ndf);
   

   
   langaupro(fp,SNRPeak,SNRFWHM);

   return fitfunc;
}

void getcenterMPV(int R){

  RunNumber = R;
  char fnam[128];

  sprintf(fnam,"localdir/run%d/tofdata_run%d.root", RunNumber, RunNumber);
  //sprintf(fnam,"calibration%d/adchists_run%d.root", RunNumber, RunNumber);
  TFile *f = new TFile(fnam,"READ");
  //f->ls();

  TCanvas *c2 = new TCanvas("c2","Attenuation length",800.,600.);

  double CenterMPV[2][44][2];

  for (int Plane=0;Plane<2;Plane++){
    
    for (int Paddle=0; Paddle<44;Paddle++) {

      int idx[2];
      idx[0] = Paddle + Plane*88;
      idx[1] = 44 + Paddle + Plane*88;
      
      for (int n=0;n<2;n++){
	CenterMPV[Plane][Paddle][n] = 0.;

	char hnam[128];
	//sprintf(hnam,"ADCHists%d",idx[n]);
	sprintf(hnam,"adc_hist%d", idx[n]);
	
	TH2D *h2d = (TH2D*)f->Get(hnam);
	TH1D *h;
	if ((Paddle==21 || Paddle == 22)){
	  h = h2d->ProjectionX("h",20,20);
	} else {
	  // h = h2d->ProjectionX("h",24,25);
	  h = h2d->ProjectionX("h",20,20);
	}

	if (h == NULL){
	  cout<<"ERROR not histogram:"<<hnam<<endl;
	  return;
	}

	if (h->GetEntries() < 50){
	  continue;
	}

	char htit[128];
	sprintf(htit,"Plane %d Paddle %d End %d",Plane, Paddle, n);
	h->SetTitle(htit);

	
	TH1D *toh;
	toh = h;
	
	TF1* TheFitFunction = langaus(toh);	
	double mpv = TheFitFunction->GetParameter(1);
	CenterMPV[Plane][Paddle][n] = mpv;
	
	if (DEBUG>1) {
	  h->Draw();
	  TheFitFunction->Draw("lsame");
	  gPad->Update();
	  //getchar();
	}
      }
    }
  }


  ofstream OUTF;
  char nam[128];
  sprintf(nam,"mpv_center_run%d.dat",R);
  OUTF.open(nam);
  
  for (int Plane=0;Plane<2;Plane++){
    for (int Paddle=0; Paddle<44;Paddle++) {
      for (int End=0; End<2;End++) {

	OUTF<<Plane<<"   "<<Paddle<<"   "<<End<<"   "<<	CenterMPV[Plane][Paddle][End]<<endl;

      }
    }
  }

  OUTF.close();

}


float findstart( TH1D *h){

  int NBins = h->GetNbinsX();

  float slopes[5];
  float values[6];
  int found = 0;
  int foundat = 0;
  for (int k=NBins; k>35; k--) {
    for (int n=0; n<6; n++){
      values[n] = h->GetBinContent(k-n);
    }
    int nf = 0;
    for (int n=0; n<5; n++){
      slopes[n] = values[n]-values[n+1];
      if (slopes[n]>0){
	nf++;
      }
    }
    if (nf>4){
      foundat = k;
      break;
    }
  }

  float x = h->GetBinCenter(foundat);

  return x;
}