#include "Math/Minimizer.h"
#include "Math/Factory.h"
#include "Math/Functor.h"
#include "TRandom2.h"
#include "TError.h"
#include "TPaveText.h"
#include "TText.h"
#include "TColor.h"

#include <iostream>


#include "TROOT.h"
#include "TSystem.h"
#include "TStyle.h"
#include "TCanvas.h"
#include "TPad.h"
#include "TH1.h"
#include "TH2.h"
#include "TF1.h"
#include "TFile.h"
#include "TTree.h"
#include "TClonesArray.h"
#include "TVector3.h"
#include "TMath.h"
#include "TChain.h"

#include <vector>
#include <iostream>

#include "../prt/PrtHit.h"
#include "../prt/PrtEvent.h" 

const Double_t prismangle = 45;

TCanvas* c1 = new TCanvas("c1","c1",500,600);
// TH2F *hHits = new TH2F("hHist","",1000,-30,190,1000,-100,100 );
TH2F *hHits = new TH2F("hHist","",1000,-20,320,1000,-100,100 );
const Int_t nmcp = 12;
const Int_t nrow = 3;
const Int_t ncol = 4;

Double_t mcpSize[] = {53,53,1};
Double_t mcpPos[nmcp][3];

Double_t ms=mcpSize[0];

Double_t prismSize[] = {30+300*tan(prismangle*TMath::Pi()/180.),170};

Double_t prismShift = (prismSize[0])/2. -30/2.;
// Double_t tX = prismShift + (prismSize[0]-mcpSize[0])/2.;
// Double_t tY = (-prismSize[1]+mcpSize[1])/2.;
Double_t total = 0;


void drawMcp(Double_t x, Double_t y, Int_t color){
  c1->cd();
  x = x + prismShift - prismSize[0]/2.;
  y = y + 0 - prismSize[1]/2.;
  TBox *box2 = new TBox(x-mcpSize[0]/2.,y-mcpSize[1]/2.,x+mcpSize[0]/2.,y+mcpSize[1]/2.);
  box2->SetFillStyle(0);
  box2->SetLineColor(color);
  box2->SetLineWidth(2);
  box2->Draw("same");
}

void drawPrism(Double_t x, Double_t y){
  c1->cd();
 
  TBox *box2 = new TBox(x-prismSize[0]/2.,y-prismSize[1]/2.,x+prismSize[0]/2.,y+prismSize[1]/2.);
  box2->SetFillStyle(0);
  box2->SetLineColor(4);
  box2->SetLineWidth(2);
  box2->Draw("same");
}

Double_t integrateMcp(Double_t x, Double_t y){
  x = x + prismShift - prismSize[0]/2.;
  y = y + 0 - prismSize[1]/2.;

  Double_t integral = 0;
  // TAxis *axis = hHits->GetXaxis();
  // int bmin = axis->FindBin(xmin);
  // int bmax = axis->FindBin(xmax);
  // double integral = h->Integral(bmin,bmax);
  // integral -= h->GetBinContent(bmin)*(xmin-axis->GetBinLowEdge(bmin))/axis->GetBinWidth(bmin);
  // integral -= h->GetBinContent(bmax)*(axis->GetBinUpEdge(bmax)-xmax)/axis->GetBinWidth(bmax);

  Int_t binx1 = hHits->GetXaxis()->FindBin(x-mcpSize[1]/2.);
  Int_t binx2 = hHits->GetXaxis()->FindBin(x+mcpSize[1]/2.);
 
  Int_t biny1 = hHits->GetYaxis()->FindBin(y-mcpSize[1]/2.);
  Int_t biny2 = hHits->GetYaxis()->FindBin(y+mcpSize[1]/2.);
  
  integral = hHits->Integral(binx1,binx2,biny1,biny2);
  
  return integral;
}



double getIntegral(const double *xa ){
  Double_t integral =0;
  for(Int_t i=0; i<nrow; i++){
    for(Int_t j=0; j<ncol; j++){
      integral += integrateMcp(xa[j],xa[ncol+i]);
    }
  }
  if(xa[5]-xa[4]<ms || xa[6]-xa[5]<ms )  integral = 1;
  //if(xa[2]-xa[1]<ms || xa[3]-xa[2]<ms )  integral = 1;
  
  if(integral==0) integral = 1;
  //std::cout<<"total "<< total <<std::endl;

  return 1/integral;
}

void optimizebb(TString path = ".", TString name = ""){
  const Int_t NRGBs = 5, NCont = 255;
  Double_t stops[NRGBs] = { 0.00, 0.34, 0.61, 0.84, 1.00 };
  Double_t red[NRGBs]   = { 0.00, 0.00, 0.87, 1.00, 0.51 };
  Double_t green[NRGBs] = { 0.00, 0.81, 1.00, 0.20, 0.00 };
  Double_t blue[NRGBs]  = { 0.51, 1.00, 0.12, 0.00, 0.00 };
  TColor::CreateGradientColorTable(NRGBs, stops, red, green, blue, NCont);
  gStyle->SetNumberContours(NCont);
  hHits->SetStats(0);
  hHits->GetXaxis()->SetTitle("x, [cm]");
  hHits->GetYaxis()->SetTitle("y, [cm]");
  hHits->GetYaxis()->SetTitleOffset(1.2);


  gSystem->Load("../build/libprtlib.so");
  TH1F *hTime = new TH1F(Form("id - %d", 0),Form("id%d", 0),500,0,25);
  
  Double_t time;
  PrtEvent* fEvent = 0;
  PrtHit hit;
  TChain *ch = new TChain("data");
  // ch->Add("../build/hits.root");
  ch->Add("/SAT/hera/had1/dervish/data/prt/iii2/hits*1_1.root");
  ch->SetBranchAddress("PrtEvent", &fEvent);
  std::cout<<"ch->GetEntries()  "<< ch->GetEntries()<<std::endl;
 
  for (Int_t ievent=0; ievent<ch->GetEntries(); ievent++){
    ch->GetEntry(ievent);
    if(ievent%100==0) cout<<"Event # "<<ievent<<endl;
    for(Int_t h=0; h<fEvent->GetHitSize(); h++){
      hit = fEvent->GetHit(h);
      time = hit.GetLeadTime(0)/1000.;
      if(time<50) hHits->Fill(hit.GetGlobalPos().X(),hit.GetGlobalPos().Y());
      hTime->Fill(time); 
    }
  }  
  total = hHits->GetEntries();
  c1->cd();
  hHits->Draw("colz");
  
  Double_t tX = mcpSize[0]/2.,tY = mcpSize[0]/2.;
  Int_t sintegral = 0;
  for(Int_t i=0; i<nrow; i++){
    tX = mcpSize[0]/2.;
    for(Int_t j=0; j<ncol; j++){
      drawMcp(tX,tY,2);
      sintegral += integrateMcp(tX,tY);
      tX += mcpSize[0];
    }
    tY += mcpSize[1];
  }
  std::cout<<"Sintegral  "<< sintegral <<std::endl;

  // algoName Migrad, Simplex,Combined,Scan  (default is Migrad)
  ROOT::Math::Minimizer* min = ROOT::Math::Factory::CreateMinimizer("Minuit", "Migrad");
  min->SetMaxFunctionCalls(100000); // for Minuit/Minuit2 
  min->SetMaxIterations(10000);  // for GSL 
  min->SetTolerance(0.000001);
  min->SetPrintLevel(1);

  ROOT::Math::Functor f(&getIntegral,7); 

  double step[] = {1,1,1,1,1,1,1};
  double a[] = {ms/2.,ms+ms/2.,2*ms+ms/2.,3*ms+ms/2.,ms/2.,ms+ms/2.,2*ms+ms/2.};
  min->SetFunction(f);
 
  for(Int_t i=0; i<7; i++){
    min->SetVariable(i,Form("a%d",i),a[i], step[i]);
    if(i<4) min->SetVariableLimits(i,i*ms+ms/2.,10000);

  }

  min->Minimize(); 
  min->SetStrategy(2);
  min->Minimize(); 

  const double *xs = min->X();

  std::cout << "Minimum: f(" << xs[0] << "," << xs[1] << "): " << min->MinValue()  << std::endl;
  if ( min->MinValue()  < 8 ) std::cout << "Converged" << std::endl;
  else  Error("NumericalMinimization","fail to converge");

  Int_t fintegral = 0;
  for(Int_t i=0; i<nrow; i++){
    for(Int_t j=0; j<ncol; j++){
      drawMcp(xs[j],xs[ncol+i],3);
      fintegral += integrateMcp(xs[j], xs[ncol+i]);
    }
  }
  
  for(Int_t i=0; i<7; i++) std::cout<<"a"<<i<<" "<< xs[i]<<std::endl;
    

  std::cout<<"Fintegral  "<< fintegral <<std::endl;
  
  drawPrism(prismShift,0);

  TPaveText *pave1 = new TPaveText(20,90,100,98);
  pave1->SetFillColor(0);
  pave1->SetBorderSize(0);
  TText *t1=pave1->AddText(Form("%d",sintegral));
  t1->SetTextFont(42);
  t1->SetTextColor(2);
  pave1->Draw();
  TPaveText *pave2 = new TPaveText(100,90,180,98);
  pave2->SetFillColor(0);
  pave2->SetBorderSize(0);
  TText *t2=pave2->AddText(Form("%d",fintegral));
  t2->SetTextColor(3);
  t2->SetTextFont(42);
  pave2->Draw();
}