//           Analysis for beta-alpha counter
//               Author: Eric Vazquez

gROOT->Reset();
#include "TCanvas.h"
#include "TMath.h"
#include "TGraph.h"

TGraph *gr;
Double_t Pi = 3.14159265358979323846;

Double_t myfcn4( Double_t *x, Double_t *par)  {
  Double_t fit;
  Double_t lambBi = 1/(60.55/60);
  Double_t lambPb = 1/10.64;
  Double_t lambRa = 1/(3.6319*24);
  Double_t lambTh = 1/(1.9116*365*24);
  Double_t ATh = par[0],ARa = par[1],APb = par[2],ABi = par[3];
  Double_t gamma = lambPb / (lambPb - lambRa);
  Double_t beta = lambBi / (lambBi - lambRa);
  Double_t alpha = lambBi / (lambBi - lambPb);
  Double_t eBi = TMath::Exp(-lambBi * x[0]);
  Double_t ePb = TMath::Exp(-lambPb * x[0]);
  Double_t eRa = TMath::Exp(-lambRa * x[0]);
  Double_t rate1 = ABi * eBi; 
  //cout << lambBi <<" "<<eBi<<endl;
  Double_t rate2 = APb * alpha * (ePb-eBi);
  Double_t rate3 = ARa * gamma * (beta*eRa - alpha*ePb + (alpha-beta)*eBi);
  Double_t rate4 = ATh * (1 - beta*gamma*eRa + (gamma-1)*alpha*ePb + (beta*gamma-alpha*gamma+alpha-1)*eBi);
  Double_t rate = rate1 + rate2 + rate3 + rate4;
  //  cout << ATh << ARa << APb << ABi << endl;
  //  cout <<"R "<<rate<<endl;
  return rate;
}

Double_t logpoisson( Double_t *x, Double_t *par)  {
  Double_t fit;
  Double_t logpoisson;
  //cout <<  myfcn4(x,par) <<endl;
  // cout <<"torate"<<endl;
  logpoisson = - myfcn4(x,par) + ( x[1] * TMath::Log( myfcn4(x,par) ) );
  //cout <<"fromrate"<<endl;
  fit = logpoisson;
  return fit;
}

void anal(const char * infile1, int chn ,const char * chain, int hours, const char * outfile1)
{
  
  gStyle->SetPalette(1);
  TCanvas *c1 = new TCanvas("c1","Time",600,600);
  TCanvas *c2 = new TCanvas("c2","First and long",600,600);
  TCanvas *c3 = new TCanvas("c3","Tail",600,600);
  TCanvas *c4 = new TCanvas("c4","Counts per hour",600,600);
  c1->Divide(1,2);
  c2->Divide(2,2);
  c3->Divide(1,2);
  //c4->Divide(2,2);
  //c1->SetGrid();
  //c2->SetGrid();
  //c3->SetGrid();
  //c4->SetGrid();
 
  Float_t runtime,pulsetime;
  Float_t beta,alpha,tail;
  Int_t channel;
  TFile *f = new TFile(infile1);
  T->SetBranchAddress("channel",&channel);
  T->SetBranchAddress("pulsetime",&pulsetime);
  T->SetBranchAddress("runtime",&runtime);
  T->SetBranchAddress("beta",&beta);
  T->SetBranchAddress("alpha",&alpha);
  T->SetBranchAddress("tail",&tail);
  Int_t nentries = (Int_t)T->GetEntries();

  Double_t var1=hours;
  Double_t parms=4;
  const char * U;  
  const char * Th;
 
  Double_t mintdc,maxtdc;
  if (chain == "U") {
    mintdc=10000;
    maxtdc=700000;
  }
  if (chain == "Th") {
    mintdc=0;
    maxtdc=1500;
  }

  Double_t minbeta=50,maxbeta=4000;
  Double_t minalpha=0,maxalpha=3500;
  Double_t mintail=0.,maxtail=1.0;

  //  Double_t downalpha=50,upalpha=4000;
  Double_t downalpha=1300,upalpha=2500;
  Double_t downbeta=50,upbeta=4000;
  Double_t downtail=0.04,uptail=0.18;

  TH1F *rt = new TH1F("rt","time since start of run in hours",var1,0.,var1);
  TH1F *pt = new TH1F("pt","time between pulses for beta and alpha in ns",100,mintdc,maxtdc);
  TH1F *bet = new TH1F("bet","integral of beta candidate in ADC counts",100,minbeta,maxbeta);
  TH1F *alph = new TH1F("alph","long integral of alpha candidate in ADC counts",100,minalpha,maxalpha);
  TH1F *tailfr = new TH1F("tailfr","tail fraction",100,mintail,maxtail);

  TH2F *longtail = new TH2F("lovsta","tail fraction vs long integral alpha",400,minalpha,maxalpha,200,mintail,maxtail);
  TH2F *longtail2 = new TH2F("lovsta2","tail fraction vs long integral alpha",400,downalpha,upalpha,200,downtail,uptail);

  TH1F *timelota = new TH1F("timelota","counts per hour vs time",var1,0.,var1);
 
  longtail->SetMarkerStyle(2); 
  longtail->SetMarkerColor(kBlue); 
  longtail2->SetMarkerStyle(2); 
  longtail2->SetMarkerColor(kGreen); 

  //            Select events, cuts, ...
  //----------------------------------------------------------------------

  for (Int_t i=0;i<nentries;i++) {
    T->GetEntry(i);
    if (channel==chn)
      {
	if (pulsetime>mintdc && pulsetime<maxtdc){
	  if (beta>minbeta && beta<maxbeta){
	    if (alpha>minalpha && alpha<maxalpha){
	      rt->Fill(runtime);
	      pt->Fill(pulsetime);
	      bet->Fill(beta);
	      alph->Fill(alpha);
	      tailfr->Fill(tail);
	      longtail->Fill(alpha,tail);
	      if (tail > downtail && tail < uptail){
		if (alpha > downalpha && alpha < upalpha){
		  if (beta > downbeta && beta < upbeta){
		    longtail2->Fill(alpha,tail);
		    timelota->SetFillColor(kGreen);
		    timelota->SetLineColor(kGreen);
		    timelota->Fill(runtime);
		  }
		}
	      }
	    }
	  }
	}
      }
  }
    
  gr = new TGraph(var1);
  Double_t bins;
  for (Int_t i=0;i<var1;i++) {
    bins[i]= timelota->GetBinContent(i+1);
    gr->SetPoint(i,bins[i],i);
  }
  
 
  // evaluate deltas
  //----------------------------------------------------------

  TF2 *fitlogpoi = new TF2("fitlogpoi", logpoisson,0,0,0,0,4);
  TF1 *fitfunc = new TF1("fitfunc", myfcn4,0,0,4);

  Double_t *xbin = gr->GetX();

  Double_t fitres;
  Double_t fitlog=0;
  Double_t logsum;
  Double_t dlldp1=0,d2lldp21=0;
  Double_t dlldp2=0,d2lldp22=0;
  Double_t dlldp3=0,d2lldp23=0;
  Double_t dlldp4=0,d2lldp24=0;
  Double_t delta1;
  Double_t inifit1,inifit2,inifit3,inifit4;
  Double_t A1Bi,A1Pb,A1Ra,A1Th;
  Double_t totfit;
  Double_t bin;
  Double_t par1=1,par2=1,par3=1,par4=1;
  
  for (Int_t i=0;i<var1;i++){
    if (i==35) {xbin[i]=5;}
    if (i==239) {xbin[i]=-1;}
    if (i==290) {xbin[i]=-1;}
    if (xbin[i] > -1){
    //log sum
    //--------------------------------
    logsum=0;
    for (Int_t j=1;j<=xbin[i];j++) {
      logsum = logsum + TMath::Log(j);
    }
    
    bin=i+1;
    //fit poisson
    //-------------------------------
    fitlogpoi->SetParameters(par1,par2,par3,par4);
    fitres[i] = fitlogpoi->Eval(bin,xbin[i]);
    fitlog = fitlog + fitres[i] - logsum;
   
    //    cout << "BIN "<< bin << endl;
    //cout << logsum << endl;
    //cout << inifit1 << endl;
    //cout << fitlog << endl;


    // fit with all params
    //----------------------------------
    fitfunc->SetParameters(1,1,1,1);
    totfit[i] = fitfunc->Eval(bin);
    
    //    cout << "mufull "<< totfit[i] <<endl;

    //fit one by one params
    //----------------------------------
    fitfunc->SetParameters(par1,0,0,0);
    inifit1 = fitfunc->Eval(bin);
    dlldp1 = dlldp1 + ( inifit1 * ( 1 - ( xbin[i] / totfit[i] ) ) );
    d2lldp21 = d2lldp21 + ( inifit1 * inifit1 * xbin[i] / ( totfit[i] * totfit[i] ) );
    delta1 = - dlldp1 / d2lldp21;
    
    fitfunc->SetParameters(0,par2,0,0);
    inifit2 = fitfunc->Eval(bin);
    dlldp2 = dlldp2 + ( inifit2 * ( 1 - ( xbin[i] / totfit[i] ) ) );
    d2lldp22 = d2lldp22 + ( inifit2 * inifit2 * xbin[i] / ( totfit[i] * totfit[i] ) );
    delta2 = - dlldp2 / d2lldp22;
    
    fitfunc->SetParameters(0,0,par3,0);
    inifit3 = fitfunc->Eval(bin);
    dlldp3 = dlldp3 + ( inifit3 * ( 1 - ( xbin[i] / totfit[i] ) ) );
    d2lldp23 = d2lldp23 + ( inifit3 * inifit3 * xbin[i] / ( totfit[i] * totfit[i] ) );
    delta3 = - dlldp3 / d2lldp23;

    fitfunc->SetParameters(0,0,0,par4);
    inifit4 = fitfunc->Eval(bin);
    dlldp4 = dlldp4 + ( inifit4 * ( 1 - ( xbin[i] / totfit[i] ) ) );
    d2lldp24 = d2lldp24 + ( inifit4 * inifit4 * xbin[i] / ( totfit[i] * totfit[i] ) );
    delta4 = - dlldp4 / d2lldp24;
    
    //get func with new params

    par1=par1+delta1;
    par2=par2+delta2;
    par3=par3+delta3;
    par4=par4+delta4;

    //fitlogpoi->SetParameters(par1,par2,par3,par4);
    //fitres[i] = fitlogpoi->Eval(bin,xbin[i]);
    //    fitlog = fitlog + fitres[i] - logsum;

    //ut << "mu(4) "<< bin <<" "<< inifit4 <<endl;
    //ut << "mu(4)-1 "<< bin-1 <<" "<< inifit4 <<endl;
    //cout << "I "<< bin << endl;

    //cout << xbin[i]<<endl;
    //cout << fitfunc->Eval(bin)<<endl;  
    //cout << inifit2<<endl;  
    //cout << inifit3<<endl;  
    //cout << inifit4<<endl;  
    

      //dlldp[k] = dlldp[k] + ( inifit[k] * ( 1 - ( xbin[i] / totfit[i] ) ) );
      //d2lldp2[k] = d2lldp2[k] + ( inifit[k] * inifit[k] * xbin[i] / ( totfit[i] * totfit[i] ) );
  }
  
  //for (Int_t k=1;k<=4;k++){
    //    deltas[k] = - dlldp[k] / d2lldp2[k];
  //}  
  //cout << "--------"<<endl;
  //cout << delta1<<endl;  
  //cout << inifit1<<endl;
  //cout << inifit2<<endl;  
  //cout << inifit3<<endl;  
  // cout << inifit4<<endl;  
  //cout << deltas[2]<<endl;  
  //cout << inifit[2]<<endl;  
  //cout << deltas[3]<<endl;  
  //cout << inifit[3]<<endl;  
  //cout << deltas[4]<<endl;  
  //  cout << inifit[4]<<endl;  
  }

  cout << fitlog << endl;
  cout << delta1 << endl;
  cout << delta2 << endl;
  cout << delta3 << endl;
  cout << delta4 << endl;

 //             Fit function with Fitter
  //--------------------------------------------------------------
  
  //  TVirtualFitter *fitter = TVirtualFitter::GetFitter();  //the current fitter
  
  //TVirtualFitter::SetDefaultFitter("Minuit");
  //TVirtualFitter *fitter = TVirtualFitter::Fitter(0,4);
  //fitter->SetFCN(myfcn4);
  
  //fitter->FitLikelihood(myfcn4)
  //  fitter->SetFitMethod(loglikelihood);
  //  fitter->FitLikelihood("");

  //fitter->SetParameter(0, "ABi",   1, 0.01, 0,5);
  //fitter->SetParameter(1, "APb",   1, 0.01, 0,5);
  //fitter->SetParameter(2, "ARa",   1, 0.01, 0,5);
  //fitter->SetParameter(3, "ATh",   1, 0.01, 0,5);
  
   //Double_t arglist[1] = {0};
  //  fitter->ExecuteCommand("MIGRAD", arglist, 0);

  //double arglist[100];
  //arglist[0] = 0;
  // set print level
  //  fitter->ExecuteCommand("SHO FCN",arglist,2);
  //fitter->ExecuteCommand("SET PRINT",arglist,2);
  //arglist[0] = 5000; // number of function calls
  //arglist[1] = 0.5; // tolerance
  //fitter->ExecuteCommand("MIGRAD",arglist,2);
  
  //-----------------------------------------------------------------

  c1->cd(1);
  rt->Draw();
  c1->cd(2);
  pt->Draw();
  
  c2->cd(1);
  bet->Draw();
  c2->cd(2);
  alph->Draw();
  c2->cd(3);
  tailfr->Draw();
  
  c3->cd(1);
  longtail->Draw("P");
  longtail2->Draw("PSAME");

  c4->cd();
  timelota->Draw();

  c1->Print("calibration-time.eps");
  c2->Print("calibration-filo.eps");
  c3->Print("calibration-tail.eps");
  c4->Print(outfile1);
 


}