#include "TF1.h"

Double_t fact(Int_t x)
{
  if(x<=1) return 1;
  else return x*fact(x-1);
} 

Double_t poissonpdf(Double_t mu, Int_t n)
{
  if( mu <= 0 )
    {
    if( n == 0 ) return 1;
    else return 0;
    }
  else return exp(-mu)*pow(mu,n)/fact(n);
}

Double_t gausspdf(Double_t m, Double_t s, Double_t x)
{
  //  return exp(-pow((x-m)/s,2)/2.)/sqrt(2*TMath::Pi())/s;
  return exp(-pow((x-m)/s,2)/2.)/s;
}

Double_t jointpoissonpdf(Double_t mu1, Double_t mu2, Int_t n)
{
  if( mu1 <= 0 || mu2 <= 0 || n <= 0 ) return poissonpdf(mu1, n);
  else 
    {
      Double_t tmpy = 0;
      for(Int_t i=1; i<=n; i++) 
	tmpy += poissonpdf(mu1,i)*poissonpdf(i*mu2,n-i);
      return tmpy;
    } 
}

Double_t poissonfunc(Double_t *x, Double_t *par)
{
  // par[0] = mu;
  // par[1] = intrinsic single photon signal width;
  // par[2] = scaling;
  // par[3] = signal photon amplitude;
  // par[4] = pedestal width;
  // par[5] = pedestal;
  // par[6] = cross-talk;
  Double_t nmax = 50;
  Double_t wn = par[4];
  Double_t tmpy = 0;
  for(Int_t i=0; i<nmax; i++)
    {
      wn = sqrt(pow(par[4],2)+i*pow(par[1],2));
      tmpy += par[2]*jointpoissonpdf(par[0],par[6],i)*gausspdf(par[5]+i*par[3],wn,x[0]);
    }
  return tmpy;
}

int fit(TH1F* ff, double par[7], double par_err[7], int xtalk = 0)
{
  // Define joint Poisson distribution.
  TF1* pf = new TF1("pf",poissonfunc,0,800,7);

  pf->SetNpx(1000);
  pf->SetLineColor(2);
  pf->SetLineWidth(1);
  pf->SetParName(0,"#mu");
  pf->SetParName(1,"#sigma_{signal}");
  pf->SetParName(2,"Amplitude");
  pf->SetParName(3,"#Delta_{p.e.}");
  pf->SetParName(4,"#sigma_{pedestal}");
  pf->SetParName(5,"Pedestal");
  pf->SetParName(6,"Cross-talk");

  // Find Pedestal
  double max = ff->GetMaximum();
  int nbin = ff->GetNbinsX();
  
  // Find the center of first peak: pedestal
  int ped = 50;
  for(int i = ped; i<=nbin; i++)
    {
      if(ff->GetBinContent(i)>max/10.)
	{
	  ped = i;
	  break;
	}
    }
  
  //  ped++;
  
  for(int i = ped; i<=nbin; i++)
    {
      if(ff->GetBinContent(i)<ff->GetBinContent(i-1))
	{
	  ped = i;
	  break;
	}
    }
  ped = ped-1;
  ped = par[5];

  // Calculate mu using Poisson distribution
  double chmax = 19.5;
  double chmin = 0.5;
  
  double del = 60;
  //  ff->SetAxisRange(ped-0.5*del,ped+chmax*del,"X");
  double mean = ff->GetMean();
  double rms = ff->GetRMS();
  double ntotal = ff->Integral();
  double mu = pow((mean-ped)/rms,2);

  pf->SetParameter(0,mu*1.2);  // mu
  pf->SetParameter(1,7);  // sigma
  pf->SetParameter(2,ntotal/2.5);  // amp
  pf->SetParameter(3,del);  // Delta
  pf->SetParameter(4,5);  // sigma_p
  pf->SetParameter(5,ped);  // ped
  pf->FixParameter(6,0);
  ff->Fit(pf,"QNL","",ped-chmin*del,ped+chmax*del);
  pf->GetParameters(par);

  del = par[3];
  //  ff->SetAxisRange(ped-chmin*del,ped+chmax*del,"X");
  mean = ff->GetMean();
  rms = ff->GetRMS();
  ntotal = ff->Integral();
  mu = pow((mean-ped)/rms,2);

  cout<<"estimated pedestal: "<<ped<<endl;
  cout<<"estimated mu: "<<mu<<endl;
  cout<<"estimated amplitude: "<<ntotal/2.5<<endl;

  pf->SetParameter(0,mu*1.2);  // mu
  pf->SetParameter(1,7);  // sigma
  pf->SetParameter(2,ntotal/2.5);  // amp
  pf->SetParameter(3,del);  // Delta
  pf->SetParameter(4,5);  // sigma_p
  pf->SetParameter(5,ped);  // ped
  pf->FixParameter(6,0);
  ff->Fit(pf,"L","",ped-chmin*del,ped+chmax*del);
  pf->GetParameters(par);
 
  if(xtalk)
    {
      pf->ReleaseParameter(6);
      pf->SetParameter(0,par[0]*.85);  // subtract X-talk
      pf->SetParLimits(6,0.05,0.35);
      pf->SetParameter(6,0.15);  // XT  
      ff->Fit(pf,"L","",par[5]-chmin*par[3],par[5]+chmax*par[3]);
      pf->GetParameters(par);
    }

  for(int i = 0; i<7; i++) par_err[i] = pf->GetParError(i);
 
 return 1;
}