#include <iostream>
#include <iomanip>
#include <fstream>
#include <iostream>
#include <sstream>
#include <numeric>
#include <vector>

#include "TCanvas.h"
#include "TSystem.h"
#include "TROOT.h"
#include "TStyle.h"
#include "TF1.h"
#include "TAxis.h"
#include "TGraph.h"

struct OverallFillFunction {
  OverallFillFunction(const std::vector<TF1 *> & flist) : fFuncList(flist) {}
  double operator() (const double * x, const double *p) {
    double result = 0;
    for (unsigned int i = 0; i < fFuncList.size(); ++i)
      result += fFuncList[i]->EvalPar(x,p);
    return result;
  }
  std::vector<TF1*> fFuncList;
};

void savePlot(TCanvas*, TString);
void SumOfNExp() {

  gROOT->SetStyle("Plain");
  gStyle->SetPalette(1,0);
  gStyle->SetOptStat(000000);
  gStyle->SetOptFit(000000);
  gStyle->SetOptTitle(0);
  gStyle->SetGridStyle(3);
  gROOT->UseCurrentStyle();

  // Fake Lumi data
  // const int nbx = 100;
  const int nbx = 150;
  std::vector<double> bxvec;
  for(int i=0; i<nbx; ++i) {
    bxvec.push_back(i);
  }

  // std::vector<double> lumivec = {0,0,0,0,0,0,0,0,0,0,99,98,97,96,95,94,93,92,91,90,0,0,0,0,0,99,98,97,96,95,94,93,92,91,90,0,0,0,0,0,
  //				 99,98,97,96,95,94,93,92,91,90,0,0,0,0,0,99,98,97,96,95,94,93,92,91,90,0,0,0,0,0,
  //                               99,98,97,96,95,94,93,92,91,90,0,0,0,0,0,99,98,97,96,95,94,93,92,91,90,0,0,0,0,0};

  std::vector<double> lumivec = {0,0,0,0,0,0,0,0,0,0,80,90,100,110,100,90,100,90,100,90,0,0,0,0,0,0,0,0,0,0,80,90,100,110,100,90,100,90,100,90,0,0,0,0,0,0,0,0,0,0,
				 80,90,100,110,100,90,100,90,100,90,0,0,0,0,0,0,0,0,0,0,80,90,100,110,100,90,100,90,100,90,0,0,0,0,0,0,0,0,0,0,80,90,100,110,100,90,100,90,100,90,
				 0,0,0,0,0,0,0,0,0,0,80,90,100,110,100,90,100,90,100,90,0,0,0,0,0,0,0,0,0,0,80,90,100,110,100,90,100,90,100,90,0,0,0,0,0,0,0,0,0,0};
  
  // Fake Lumi Graph
  TGraph * LumiGr = new TGraph(nbx, &bxvec[0], &lumivec[0]);
  LumiGr->SetMarkerColor(kBlue);  LumiGr->SetLineColor(kBlue);  LumiGr->SetMarkerSize(0.5); LumiGr->SetMarkerStyle(20);
  
  
  // ====================================================
  // Multiply each L_i[bx] with EXP(-T/TAU) [in bx units]
  // ====================================================
    
  std::vector<TF1 *> func_vec;
  for(int i=0; i<nbx; ++i) {
    
    std::string fncname = "f";
    fncname = fncname+std::to_string(i);
    
    TF1 * f = new TF1(fncname.c_str(),"(x<[0])*[1] + (x>=[0] && x<[0]+10)*[2]*exp(-1*(x-[0])*[3]) + (x>=[0]+10)*[1]",0,nbx-1);
    f->FixParameter(0,i);
    f->FixParameter(1,0);
    f->FixParameter(2,lumivec[i]);
    f->SetParameter(3,10); // 10 non funziona
    func_vec.push_back(f);
    
    std::cout<<"Creating function f"<<i<<" = (x<"<<f->GetParameter(0)<<") * "<<f->GetParameter(1)<<" + ";
    std::cout<<"(x>="<<f->GetParameter(0)<<" && x<"<<f->GetParameter(0)+10<<") * "<<f->GetParameter(2)<<" * exp( - ( x - "<<f->GetParameter(0)<<") *  "<<f->GetParameter(3)<<") + ";
    std::cout<<"(x>="<<f->GetParameter(0)+10<<") * "<<f->GetParameter(1)<<std::endl;
    std::cout<<"| f"<<i<<"("<<i-1<<") = "<<f->Eval(i-1)<<" | f"<<i<<"("<<i<<") = "<<f->Eval(i)<<" | f"<<i<<"("<<i+1<<") = "<<f->Eval(i+1)<<" | f"<<i<<"("<<i+2<<") = "<<f->Eval(i+2)<<" | "<<std::endl;
  }
  
  // sum all exponentials together 
  TF1 * fsum = new TF1("fsum",OverallFillFunction(func_vec),0,nbx-1,0);
  
  double cumulative_sum = 0;
  for(int i=0; i<nbx; ++i) {
    double funcval = func_vec[i]->Eval(0);
    cumulative_sum += funcval;
    std::cout<<" f"<<i<<"(0) = "<<funcval<<" cumulative sum = "<<cumulative_sum<<std::endl;
  }
  std::cout<<" f_sum(0) = "<<fsum->Eval(0)<<std::endl;
  
  // Make TGraph that samples fsum at each bx value
  std::vector<double> lumivec_exp;
  for(int i=0; i<nbx; ++i) {
    lumivec_exp.push_back(fsum->Eval(i));
  }
  TGraph * Gsum = new TGraph(nbx, &bxvec[0], &lumivec_exp[0]);
  Gsum->SetLineColor(kMagenta); Gsum->SetLineWidth(2); Gsum->SetMarkerColor(kMagenta); Gsum->SetMarkerSize(0.5); Gsum->SetMarkerStyle(20);

  // Modify Parameter in all fitting functions
  std::vector<double> lumivec_exp_mod;
  for(int j=0; j<nbx; ++j) {
    func_vec[j]->SetParameter(3,0.5);
  }
  for(int i=0; i<nbx; ++i) {
    lumivec_exp_mod.push_back(fsum->Eval(i));
  }
  TGraph * Hsum = new TGraph(nbx, &bxvec[0], &lumivec_exp_mod[0]);
  Hsum->SetLineColor(kRed); Hsum->SetLineWidth(2); Hsum->SetMarkerColor(kRed); Hsum->SetMarkerSize(0.5); Hsum->SetMarkerStyle(20);
  
  TCanvas *c1  = new TCanvas("c1", "c1", 0,0,1000,600);
  c1->SetFillStyle(4000);    c1->SetFillColor(4000); c1->SetFrameFillColor(4000); c1->SetFrameFillStyle(4000);
  c1->SetFrameBorderMode(0); c1->SetTicks(1,1);      c1->SetGrid(1,1);            // c1->SetLogy(1); // c1->SetLogx(1);
  c1->SetLeftMargin(0.025);   c1->SetRightMargin(0.025); c1->SetTopMargin(0.05);
  c1->Divide(1,3);
  
  c1->cd(1); // Show the Lumi
  LumiGr->Draw("APLE");
  LumiGr->GetXaxis()->SetRangeUser(0,nbx-1);
  
  c1->cd(2); // Show the sum of exp decays 
  Gsum->Draw("APLE");
  Gsum->GetXaxis()->SetRangeUser(0,nbx-1);

  c1->cd(3); // Show the modified sum of exp decays 
  Hsum->Draw("APLE");
  Hsum->GetXaxis()->SetRangeUser(0,nbx-1);

  
  // Gsum->Print();
  
  c1->Update();
  c1->cd(1);
  savePlot(c1,"SumOfNExp");

}

void savePlot(TCanvas * c, TString name) {
  c->SaveAs(name+".png");
  c->SaveAs(name+".eps");
  c->SaveAs(name+".pdf");
  std::cout<<"Plot "<<name<<" saved"<<std::endl;
}