#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 "TGraph.h"

#include <vector>
#include <iostream>

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

const Double_t prismangle = 45;

TCanvas* c1 = new TCanvas("c1","c1",900,600);
// TH2F *hHits = new TH2F("hHist","",1000,-3,19,1000,-10,10 );
TH2F *hHits = new TH2F("hHist","",500,-40,340,500,-100,100 );

TH2F *hHitsa[140];
const Int_t nmcp = 12;
const Int_t nrow = 3;
const Int_t ncol = 6;

Double_t offset = -9;
Double_t shiftx=1;
Double_t shifty=4.5;

// const Int_t ncol = 5;

// Double_t offset = -15;
// Double_t shiftx=14;
// Double_t shifty=1;

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

Double_t ms=mcpSize[0];

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

Double_t prismShift = (prismSize[0])/2. -50/2.;
Double_t total = 0;

void drawMcp(Double_t x, Double_t y, Int_t color){
  gPad->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){
  gPad->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;

  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_t integrateMcpa(TH2F *fHits, Double_t x, Double_t y){
  x = x + prismShift - prismSize[0]/2.;
  y = y + 0;// - prismSize[1]/2.;

  Double_t integral = 0;

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

double getIntegral(const double *xa ){
  Double_t integral =0;
 
  Double_t tX, tY = -(mcpSize[0]+xa[4]);
  for(Int_t i=0; i<nrow; i++){
    tX = -mcpSize[0]/2.;
    for(Int_t j=0; j<ncol; j++){
      tX += mcpSize[0]+xa[j];
      integral += integrateMcp(tX,tY);
    }
    tY += mcpSize[1]+xa[4];
  }
  if(integral==0) integral = 1;
  return 1/integral;
}

void optimize(TString path = ".", TString name = ""){
  Int_t saveflag = -1;
  TString info = "MCP layout optimization.";
  path = createDir("rdata", info, saveflag); 

  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, [mm]");
  hHits->GetYaxis()->SetTitle("y, [mm]");
  hHits->GetYaxis()->SetTitleOffset(1.2);


  for(Int_t i=0; i<141; i++){
    hHitsa[i] = new TH2F( Form("%d degree", i),Form("%d", i),500,-20,320,500,-100,100 );
    hHitsa[i]->SetTitle(Form("path # %d",i));
    hHitsa[i]->GetXaxis()->SetTitle("x, [cm]");
    hHitsa[i]->GetYaxis()->SetTitle("y, [cm]");
    hHitsa[i]->GetXaxis()->SetTitleOffset(1.15);
    hHitsa[i]->GetYaxis()->SetTitleOffset(1.2); 
    hHitsa[i]->SetStats(0);
  }

  gSystem->Load("../build/libprtlib.so");
  TH1F *hTime = new TH1F(Form("id - %d", 0),Form("id%d", 0),500,0,40);
  
  Double_t time;
  PrtEvent* fEvent = 0;
  PrtHit hit;
  TChain *ch = new TChain("data");
  // ch->Add("../build/hits0.root");
  
  //ch->Add("/SAT/hera/had1/dervish/data/prt/opt1/hits*.root");
  //ch->Add("/SAT/hera/had1/dervish/data/prt/opt2/hits*.root");
  ch->Add("/SAT/hera/had1/dervish/data/prt/opt325/hits*.root"); //pi 1.7 //opt31
  // ch->Add("/SAT/hera/had1/dervish/data/prt/opt45/hits*.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);
    Int_t angle = fEvent->GetAngle() + 0.001;
    
    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.;
      Int_t nref = hit.GetNreflectionsInPrizm(); 
      //      if(time<30)
      if(nref<4){
	hHits->Fill(hit.GetGlobalPos().X(),hit.GetGlobalPos().Y());
	hHitsa[angle]->Fill(hit.GetGlobalPos().X(),hit.GetGlobalPos().Y());
      }
      
      hTime->Fill(time); 
    }
  }  
  total = hHits->GetEntries();
  c1->cd();
  hHits->Draw("colz");
  
  Double_t tX, tY = -(mcpSize[0]+7);
  Int_t sintegral = 0;
  for(Int_t i=0; i<nrow; i++){ 
    tX = -mcpSize[0]/2.;
    for(Int_t j=0; j<ncol; j++){
      tX += mcpSize[0]+7;
      drawMcp(tX,tY,2);
      sintegral += integrateMcp(tX,tY);
    }
    tY += mcpSize[1]+7;
  }
  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.0000001);
  min->SetPrintLevel(1);

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

  double step[] = {1,1,1,1,1};
  double a[] = {7,7,7,7,7};
  min->SetFunction(f);
 
  for(Int_t i=0; i<5; i++){
    min->SetVariable(i,Form("a%d",i),a[i], step[i]);
    min->SetVariableLimits(i,1,10000);
  }
  min->SetVariableLimits(0,0,10000);
  min->SetVariableLimits(4,4.5,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");


  tY = -(mcpSize[0]+xs[4]);
  Int_t fintegral = 0;
  for(Int_t i=0; i<nrow; i++){
    tX =  -mcpSize[0]/2.;
    for(Int_t j=0; j<ncol; j++){
      tX += mcpSize[0]+xs[j];
      drawMcp(tX,tY,3);
      fintegral += integrateMcp(tX,tY);
    }
    tY += mcpSize[1] + xs[4];
  }

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

  TPaveText *pave1 = new TPaveText(20,100,100,110);
  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,100,180,110);
  pave2->SetFillColor(0);
  pave2->SetBorderSize(0);
  TText *t2=pave2->AddText(Form("%d",fintegral));
  t2->SetTextColor(3);
  t2->SetTextFont(42);
  pave2->Draw();

  save(c1,path,"layout",info,saveflag);

  TCanvas* cTime = new TCanvas("cTime","cTime",900,600);
  hTime->Draw();
  save(cTime,path,"time",info,saveflag);

  for(Int_t i=0; i<5; i++) std::cout<<"a"<<i<<" "<< xs[i]<<std::endl;



  TCanvas* c4 = new TCanvas("c4","c4",900,600);
  // TPad *pad1 = new TPad("pad1", "The pad 60% of the height", 0,0.4, 1.0,1.0, 0);
  // TPad *pad2 = new TPad("pad2", "The pad 40% of the height", 0,0.0, 1.0,0.4, 0);
  // pad1->Draw();
  // pad2->Draw();

  for(Int_t ist=1; ist<2; ist++){ //0-50
  
    // pad1->cd();
    hHits->Draw("colz");
    Double_t fs[ncol+1];
    for(Int_t i=0; i<ncol; i++) fs[i] = ist*shiftx;
    fs[0]=offset;
    fs[ncol]=shifty;
  
    Int_t gpoints=0, s1integral = 0;
    TGraph *gDet1 = new TGraph();
    for(Int_t ia=22; ia<140; ia++){
      Int_t tintegral = 0;
      tY = -(mcpSize[0]+fs[ncol]);
      for(Int_t i=0; i<nrow; i++){
	tX =  -mcpSize[0]/2.;
	for(Int_t j=0; j<ncol; j++){
	  tX += mcpSize[0]+fs[j];
	  // if(ia==22) drawMcp(tX,tY,1);
	  tintegral += integrateMcpa(hHitsa[ia],tX,tY);
	}
	tY += mcpSize[1] + fs[ncol];
      }
      //std::cout<<"tintegral "<<tintegral <<std::endl;
      s1integral += tintegral;
      gDet1->SetPoint(gpoints++,ia,tintegral);
    }
 
    for(Int_t i=0; i<ncol; i++) fs[i] = ist*shiftx;
    fs[0]=offset;
    fs[ncol]=shifty;
  
    // TCanvas* ct = new TCanvas("ct","ct",900,600);
    gpoints=0; 
    Int_t s2integral = 0;
    TGraph *gDet2 = new TGraph();
    TGraph *gDet3 = new TGraph();
    for(Int_t ia=22; ia<140; ia++){
      // ct->cd();
      // hHitsa[ia]->Draw("colz");
      // ct->Update();
      // ct->Modified();
      // ct->WaitPrimitive();
      Int_t tintegral = 0;
      tY = -(mcpSize[0]+fs[ncol]);
      for(Int_t i=0; i<nrow; i++){
	tX =  -mcpSize[0]/2.;
	if(i==1)  tX -= offset;
	for(Int_t j=0; j<ncol; j++){
	  tX += mcpSize[0]+fs[j];
	  if(ia==22) drawMcp(tX,tY,2);
	  tintegral += integrateMcpa(hHitsa[ia],tX,tY);
	}
	tY += mcpSize[1] + fs[ncol];
      }
      //std::cout<<"tintegral "<<tintegral <<std::endl;
      s2integral += tintegral;
      gDet2->SetPoint(gpoints++,ia,tintegral);
      gDet3->SetPoint(gpoints-1,ia,hHitsa[ia]->GetEntries());
    }


    TPaveText *pave1 = new TPaveText(20,100,100,110);
    pave1->SetFillColor(0);
    pave1->SetBorderSize(0);
    TText *t1=pave1->AddText(Form("%d",s1integral));
    t1->SetTextFont(42);
    t1->SetTextColor(1);
    pave1->Draw();
    TPaveText *pave2 = new TPaveText(100,100,180,110);
    pave2->SetFillColor(0);
    pave2->SetBorderSize(0);
    TText *t2=pave2->AddText(Form("%d",s2integral));
    t2->SetTextColor(2);
    t2->SetTextFont(42);
    pave2->Draw();

    // TCanvas* cDet = new TCanvas("cDet","cDet",900,600);
    // pad2->cd();
    gDet1->SetLineColor(1);
    gDet1->GetXaxis()->SetTitle("x gap, [mm]");
    gDet1->GetYaxis()->SetTitle("entries, [#]");
    gDet1->GetXaxis()->SetTitleOffset(0.9);
    gDet1->GetYaxis()->SetTitleOffset(0.9); 
    gDet1->GetXaxis()->SetTitleSize(0.05);
    gDet1->GetYaxis()->SetTitleSize(0.05); 
    gDet1->GetXaxis()->SetLabelSize(0.05);
    gDet1->GetYaxis()->SetLabelSize(0.05); 

    // gDet1->Draw("APL");
    // gDet2->SetLineColor(2);
    // gDet2->Draw("same");
    // gDet3->SetLineColor(4);
    // gDet3->Draw("same");

    // c4->Update();
    // c4->Modified();
    // c4->WaitPrimitive();
    save(c4,path,Form("animtime_%d",ist),info,saveflag);
  }
}