#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 "save.C"

#include "../prt/PrtTrackInfo.h"
#include "../prt/PrtPhotonInfo.h"
#include "../prt/PrtAmbiguityInfo.h"

#include <iostream>
TString info, path;
Int_t saveflag = 1;
TCanvas* c = new TCanvas("c","c",0,0,400,600);
TPad *pad1 = new TPad("pad1", "The pad 80% of the height", 0,0.2, 1.0,1.0, 0);
TPad *pad2 = new TPad("pad2", "The pad 20% of the height", 0,0.0, 1.0,0.2, 0);
TF1 *fg = new TF1("fgaus","[0]*exp(-0.5*((x-[1])/[2])*(x-[1])/[2]) + [3] + [4]*x",-0.04,0.04);

void show(int ievent, TH1F *h1,TH1F *h2,TH1F *h3,  Double_t trackTheta, Double_t trackPhi){

  h2->SetTitle(Form("Track #%d with (#theta, #varphi) = (%2.1f^{o}, %2.1f^{o})", ievent, trackTheta, trackPhi));

  pad1->cd();

  fg->SetParLimits(0,0,10000);
  //h2->SetMaximum(h2->GetMaximum()*1.4);  
  fg->SetParameter(0,50);   // peak
  fg->SetParameter(1,0);   // peak
  fg->FixParameter(2,0.009); // width
  fg->SetParameter(3,5);   // p0
  fg->SetParameter(4,1);   // p1
  Int_t status = h2->Fit("fgaus","MQ","",-0.04,0.04);
  fg->ReleaseParameter(2); // width
  status = h2->Fit("fgaus","MQ","",-0.04,0.04);
  Double_t sigma = fg->GetParameter(2)*1000;
  Double_t mean = fg->GetParameter(1)*1000;

  h2->GetXaxis()->SetTitle("#theta_{track} - #theta_{LUT}, [rad]");
  h2->GetYaxis()->SetTitle("entries, [#]");
  h2->Draw();

  h3->SetLineColor(2);
  h3->GetXaxis()->SetTitle("#theta_{track} - #theta_{LUT}, [rad]");
  h3->GetYaxis()->SetTitle("entries, [#]");
  h3->Draw("same");

  pad2->cd();

  h1->SetFillColor(11);
  h1->SetFillStyle(3001);
  h1->GetXaxis()->SetTitle("#theta_{track} - #theta_{LUT}, [rad]");
  h1->GetYaxis()->SetTitle("entries, [#]");
  h1->Draw();
  h1->ShowPeaks(0.01,"",0.7);
  
  c->Modified();
  c->Update();
  //c->WaitPrimitive();
  Int_t thetaid = trackTheta;
  Int_t sigmaid = sigma+0.5;
  Int_t meanid = mean+20;
  save(c,path,Form("dist_%d_%d_%d",thetaid,sigmaid,meanid),info,saveflag);
  

  h1->Reset();
  h2->Reset();
  h3->Reset();
}


void cdrawrecoprj(TString inFile = "../build/reco.root")
{
  inFile = "/SAT/hera/had1/dervish/data/prt/proton27/reco_*.root";
  
  gSystem->Load("../build/libprtlib.so");

  info = "spr: Single photon resolution and mean. Default geometry 10. Ideal lens. Chromatic correction. LUT from averaged photon directions.";
  path = createDir("rdata",info, saveflag); 
  Bool_t bEffMap = true;
  Int_t naver=1;

  gStyle->SetOptFit(111);
  TString outFile = "res.root";
  pad1->Draw();
  pad2->Draw();

  TClonesArray* fTrackInfo = new TClonesArray("PrtTrackInfo");
  TChain *ch = new TChain("recodata");
  ch->Add(inFile);
  ch->SetBranchAddress("PrtTrackInfo", &fTrackInfo);
  
  const Int_t nphibins= 100;
  const Int_t nthetabins= 118;
  const Double_t phimin=0, phimax=16,thetamin=22, thetamax=140;
 
  Int_t phiid=0, thetaid=0, evReflections;
  Double_t angle, theta, phi, difftime, mcCherenkov, cherenkov, 
    hitTime, mcTimeInBar, barTime, evTime,time;
  PrtTrackInfo *track;
  PrtPhotonInfo photon;
  PrtAmbiguityInfo ambiguity;
  Bool_t reflected;

  Int_t aEntr[500000] = {0};


  Int_t ids[nphibins][nthetabins] = {{0}};
  
  Double_t dphi=(phimax-phimin)/(Double_t)nphibins;
  Double_t dtheta=(thetamax-thetamin)/(Double_t)nthetabins;

  
  TH1F *hist[nphibins][nthetabins], *hist1[nphibins][nthetabins], *histfull[nphibins][nthetabins];

  for(int i=0; i<nphibins;i++){
    for(int j=0; j<nthetabins;j++){
      hist[i][j] = new TH1F(Form("h_1p%d_t%d", i,j),Form("h_%d%d", i,j), 40,-0.07,0.07);
      if(!bEffMap){
	hist1[i][j] = new TH1F(Form("h_2p%d_t%d", i,j),Form("h_%d%d", i,j), 40,-0.07,0.07);
	histfull[i][j] = new TH1F(Form("h_f%d_t%d", i,j),Form("h_%d%d", i,j), 500,-0.8,0.8);
      }
    }
  }

  TH2F * hEffMap = new TH2F("Eff_map","Lut Info"  ,nthetabins,thetamin,thetamax,nphibins,phimin,phimax);
  TH2F * hMeanMap = new TH2F("Mean_map","Lut Info",nthetabins,thetamin,thetamax,nphibins,phimin,phimax);

  TH2F * hResol = new TH2F("Resol","SPR vs. Theta",nthetabins,thetamin,thetamax,100,0,30);
  TH2F * hMean = new TH2F("Mean","Mean vs. Theta" ,nthetabins,thetamin,thetamax,100,-20,20);

  std::cout<<" End of initialization "<<std::endl;
  std::cout<<"# of entries = "<<ch->GetEntries() <<std::endl;
  for (Int_t ievent=0; ievent<ch->GetEntries(); ievent++){
    ch->GetEntry(ievent);
    if(ievent%1000==0) cout<<"Event # "<<ievent<<endl;
    for (Int_t t=0; t<fTrackInfo->GetEntriesFast(); t++){ 
      track= (PrtTrackInfo*) fTrackInfo->At(t);

      //      theta = track->GetMcMomentum().Theta()*180/TMath::Pi();
      //      phi = track->GetMcMomentum().Phi()*180/TMath::Pi();
      theta = track->GetAngle()+0.5;
      std::cout<<"theta  "<<theta <<std::endl;
      
      phi = 0;

      //if(phi<0) phi +=360;
      phiid=(phi-phimin)/dphi;
      thetaid=(theta-thetamin)/dtheta;

      //mcCherenkov = track->GetMcCherenkov();
      mcTimeInBar = track->GetMcTimeInBar();
      
      for(Int_t p=0; p<track->GetPhotonSize(); p++){
	photon = track->GetPhoton(p);
	hitTime = photon.GetHitTime();
	reflected = photon.GetReflected();
	evReflections = photon.GetEvReflections();
	// theta = photon.GetMcPrimeMomentumInBar().Theta()*180/TMath::Pi();
	// phi = photon.GetMcPrimeMomentumInBar().Phi()*180/TMath::Pi();
	// phiid=(phi-phimin)/dphi;
	// thetaid=(theta-thetamin)/dtheta;
	// std::cout<<"phiid "<<phi <<"   thetaid  " << thetaid <<std::endl;
	
	mcCherenkov = photon.GetMcCherenkovInBar();
		  
	for(Int_t a=0; a<photon.GetAmbiguitySize(); a++){
	  ambiguity = photon.GetAmbiguity(a);
	  cherenkov = ambiguity.GetCherencov();
	  barTime = ambiguity.GetBarTime();
	  evTime = ambiguity.GetEvTime();
	  angle = cherenkov - mcCherenkov;
	  time = barTime + evTime;
	  difftime = time - (hitTime-mcTimeInBar);
 
	  // if(fabs(difftime)<2.0){
	  //   double b1 = -13.2+0.021*theta;
	  //   double b2 = 8.4-0.27*theta;
	  //   double corr = -b1*difftime;
	  //   if(theta>80) corr = -b2*difftime;
	  //   angle += corr/1000.;
	  // }
	  if(!bEffMap) histfull[phiid][thetaid]->Fill(angle);
	  hist[phiid][thetaid]->Fill(angle);
	  
	  // if(difftime>0) angle += 0.004;
	  // else  angle -= 0.004;
	
	   // if(fabs(difftime)<0.5){
	   //   if(difftime>0) angle += difftime*20/1000.;
	   //  else  angle -= difftime*10/1000.;
	   // }


	  // if(difftime<0) hist[phiid][thetaid]->Fill(angle);
	  // else  hist1[phiid][thetaid]->Fill(angle);
	}
      }
    }

    if(bEffMap){
      ids[phiid][thetaid]++;
      if(ids[phiid][thetaid]>=naver && hist[phiid][thetaid]->GetEntries()>0){
	fg->SetParLimits(0,0,5000);     
	fg->SetParameter(0,50);   // peak
	fg->SetParameter(1,0);   // peak
	fg->FixParameter(2,0.009); // width
	fg->SetParameter(3,5);   // p0
	fg->SetParameter(4,1);   // p1
	Int_t status = hist[phiid][thetaid]->Fit("fgaus","MQ","",-0.04,0.04);
	fg->ReleaseParameter(2); // width
	status = hist[phiid][thetaid]->Fit("fgaus","MQ","",-0.04,0.04);
	if(status==4000){
	  Double_t sigma = fg->GetParameter(2)*1000;
	  Double_t mean = fg->GetParameter(1)*1000;
	  if(sigma < 32 && sigma > 0 && fabs(mean)<20 )  {
	    if(sigma < 18){ 
	      hEffMap->Fill(theta,phi,sigma);
	      hMeanMap->Fill(theta,phi,mean);
	      aEntr[hEffMap->FindBin(theta,phi)]++;
	    }
	    hResol->Fill(theta,sigma);
	    hMean->Fill(theta,mean);
	  }
	}
	hist[phiid][thetaid]->Reset();
	ids[phiid][thetaid]=0;
      }
    }else{
      ids[phiid][thetaid]++;
      if(ids[phiid][thetaid]>=naver && hist[phiid][thetaid]->GetEntries()){
	show(ievent, histfull[phiid][thetaid], hist[phiid][thetaid], hist1[phiid][thetaid], theta, phi);
	ids[phiid][thetaid]=0;
      }
    }
  }

  for(int i=0; i<hEffMap->GetSize()-2;i++){
    if(hEffMap->GetBinContent(i)>0) hEffMap->SetBinContent(i, hEffMap->GetBinContent(i)/(Double_t)aEntr[i]);
    if(fabs(hMeanMap->GetBinContent(i))>0) hMeanMap->SetBinContent(i, hMeanMap->GetBinContent(i)/(Double_t)aEntr[i]);
  }
  if(bEffMap){
    c->Clear();
    c->Divide(1,2);
    // c->cd(1);
    // hEffMap->GetXaxis()->SetTitle("#theta, [degree]");
    // hEffMap->GetYaxis()->SetTitle("#varphi, [degree]");
    // hEffMap->SetMaximum(14);
    // hEffMap->SetMinimum(4);
    // hEffMap->Draw("colz");
    
    // c->cd(2);
    // hMeanMap->GetXaxis()->SetTitle("#theta, [degree]");
    // hMeanMap->GetYaxis()->SetTitle("#varphi, [degree]");
    // hMeanMap->Draw("colz");
   
    c->cd(1);
    hResol->GetXaxis()->SetTitle("#theta, [degree]");
    hResol->GetYaxis()->SetTitle("#sigma, [mrad]");
    hResol->Draw("colz");

    c->cd(2);
    hMean->GetXaxis()->SetTitle("#theta, [degree]");
    hMean->GetYaxis()->SetTitle("mean, [mrad]");
    hMean->Draw("colz");

    gROOT->SetBatch(kTRUE);
    TCanvas* cs = new TCanvas("cs","cs",0,0,400,600);    
    hEffMap->Draw("colz");
    save((TPad*)gPad,path,"spr",info,saveflag);
    hMeanMap->Draw("colz");
    save((TPad*)gPad,path,"mean",info,saveflag);
    hResol->Draw("colz"); 
    save((TPad*)gPad,path,"spr_prj",info,saveflag);
    hMean->Draw("colz");
    save((TPad*)gPad,path,"mean_prj",info,saveflag);
    gROOT->SetBatch(kFALSE);
  }
  
}