#include <TH1F.h>
#include <TH1D.h>
#include <TH2F.h>
#include <TH2D.h>
#include <TFile.h>
#include <TF1.h>
#include <TProfile.h>
#include <TTree.h>
#include <unistd.h>

void tpcp0d_new(string filename){

  gStyle->SetLineScalePS(1);
  gStyle->SetOptFit(11111111); 
  gStyle->SetLineWidth(1);
  TString oFile("tpcp0dOut");
#define ANG 1.0
#define MAXHITS 2000
  // Define histograms
  
  TH1F *dy10 = new TH1F("dy10","dy, #Delta t < 10 ns",100.,-500.,500.);
  TH1F *dx10 = new TH1F("dx10","dx, #Delta t < 10 ns",100.,-500.,500.);
  TH1F *dy50 = new TH1F("dy50","dy, #Delta t < 50 ns",100.,-200.,200.);
  TH1F *dx50 = new TH1F("dx50","dx, #Delta t < 50 ns",100.,-200.,200.);
  TH1F *dy100 = new TH1F("dy100","dy, #Delta t < 100 ns",100.,-500.,500.);
  TH1F *dx100 = new TH1F("dx100","dx, #Delta t < 100 ns",100.,-500.,500.);
  TH1F *dy1000 = new TH1F("dy1000","dy, #Delta t < 1000 ns",100.,-500.,500.);
  TH1F *dx1000 = new TH1F("dx1000","dx, #Delta t < 1000 ns",100.,-500.,500.);
  
  formatHist(dy10);  formatHist(dy50);  formatHist(dy100);      formatHist(dy1000);
  formatHist(dx10);  formatHist(dx50);  formatHist(dx100);      formatHist(dx1000);

  // dj vs. k
  TH2F *dxvsx = new TH2F("dxvsx","dx vs x",10,-900.,900.,20,-20.,20.); 
  TH2F *dxvsy = new TH2F("dxvsy","dx vs y",10,-900.,900.,20,-20.,20.); 
  TH2F *dyvsx = new TH2F("dyvsx","dy vs x",10,-900.,900.,20,-20.,20.); 
  TH2F *dyvsy = new TH2F("dyvsy","dy vs y",10,-900.,900.,20,-20.,20.); 
  
  // dj vs. tk
  TH2F *dyvstx = new TH2F("dyvstx","  ",20,-ANG,ANG,20,-20.,20.);
  TH2F *dyvsty = new TH2F("dyvsty","  ",20,-ANG,ANG,20,-20.,20.);
  TH2F *dxvstx = new TH2F("dxvstx","  ",20,-ANG,ANG,20,-20.,20.);
  TH2F *dxvsty = new TH2F("dxvsty","  ",20,-ANG,ANG,20,-20.,20.);
  
  // tk
  TH1F *tyT = new TH1F("tyT","ty",100.,-ANG,ANG);
  TH1F *txT = new TH1F("txT","tx",100.,-ANG,ANG);
  TH1F *txN = new TH1F("txN","tx (x<0)",100.,-ANG,ANG);
  TH1F *txP = new TH1F("txP","tx (x>0)",100.,-ANG,ANG);
  
  // tk vs. j
  TH2F *txvsx = new TH2F("txvsx","  ",10,-900.,900.,20,-20.,20.); 
  TH2F *txvsy = new TH2F("txvsy","  ",10,-900.,900.,20,-20.,20.); 
  TH2F *tyvsx = new TH2F("tyvsx","  ",10,-900.,900.,20,-20.,20.); 
  TH2F *tyvsy = new TH2F("tyvsy","  ",10,-900.,900.,20,-20.,20.); 
  


  // Read in data files
  TChain *chain = new TChain("tpcp0dtree");
  chain->AddFile(filename.c_str(),0);
  chain->UseCurrentStyle();
  chain->GetEntry(0);
  gROOT->ForceStyle(1);
  double xtrack; 
  double ytrack; 
  double ztrack;
  double tx;
  double ty; 
  double rho;
  double time; 
  int nyhits; 
  double y[MAXHITS];
  double thity[MAXHITS];
  double yz[MAXHITS];
  double yx[MAXHITS];
  double dy[MAXHITS];
  double qy[MAXHITS];
  int nxhits; 
  double x[MAXHITS];
  double thitx[MAXHITS];
  double xz[MAXHITS];
  double xy[MAXHITS];
  double dx[MAXHITS];
  double qx[MAXHITS];

  // TTree *tree = (TTree*)chain->GetTree();

  chain->SetBranchAddress("xtrack",&xtrack);    // x coord of start of tpc track
  chain->SetBranchAddress("ytrack",&ytrack);    // y coord ...
  chain->SetBranchAddress("ztrack",&ztrack);    // z coord ...
  chain->SetBranchAddress("txtrack",&tx);       // dx/dz slope of start of track
  chain->SetBranchAddress("tytrack",&ty);       // dy/dz slope ...
  chain->SetBranchAddress("rho",&rho);          // curvature of track
  chain->SetBranchAddress("time",&time);        // initial time of track
  chain->SetBranchAddress("nyhit",&nyhits);     // number of p0d y hits
  chain->SetBranchAddress("thity",thity);               // y coord of P0D y hit
  chain->SetBranchAddress("y",y);               // y coord of P0D y hit
  chain->SetBranchAddress("yx",yx);             // x-coord of extrap TPC track
  chain->SetBranchAddress("yz",yz);             // z-coord of p0d y hit
  chain->SetBranchAddress("dy",dy);             // = yy - y (residual)
  chain->SetBranchAddress("nxhit",&nxhits);     // same for x
  chain->SetBranchAddress("x",x);               //
  chain->SetBranchAddress("thitx",thitx);               //
  chain->SetBranchAddress("xz",xz);             //
  chain->SetBranchAddress("xy",xy);             //
  chain->SetBranchAddress("dx",dx);             //
  
  Long64_t nentries = chain->GetEntries();
  
  std::cout << " ENTRIES : " << nentries << std::endl; 

  for (Long64_t iev = 0; iev < nentries; iev++){
    chain->GetEntry(iev);

    if( TMath::Abs(time) < 0.01 ) continue;

    for(int i = 0; i < nyhits; i++ ){ // Loop over y hit
      
      tyT->Fill(ty);
      if(sqrt(pow(time-thity[i],2)) < 10){
        dy10->Fill(dy[i]);
      }
      if( (time - thity[i] < 20) && (time - thity[i] > -15)){
      //if( (time - thity[i] < 60) && (time - thity[i] > -25)){
        dy50->Fill(dy[i]);
      }
      if(sqrt(pow(time-thity[i],2)) < 100){
        dy100->Fill(dy[i]);
      }
      if(sqrt(pow(time-thity[i],2)) < 1000){
        dy1000->Fill(dy[i]);
      }
      
      dyvsx->Fill(yx[i],dy[i]);
      dyvsy->Fill(y[i],dy[i]);

      dyvsty->Fill(ty,dy[i]);
      dyvstx->Fill(tx,dy[i]);
  
      tyvsx->Fill(yx[i],ty);
      tyvsy->Fill(y[i],ty);
    }
    
    for(int i = 0; i < nxhits; i++ ){  // loop over x hits


      txT->Fill(tx);
      
      if(sqrt(pow(time-thitx[i],2)) < 10){
        dx10->Fill(dx[i]);
      }
      if( time - thitx[i] < 20 && time - thitx[i] > -15){
      //if( time - thitx[i] < 40 && time - thitx[i] > -35){
        dx50->Fill(dx[i]);
      }
      if(sqrt(pow(time-thitx[i],2)) < 100){
        dx100->Fill(dx[i]);
      }
      if(sqrt(pow(time-thitx[i],2)) < 1000){
        dx1000->Fill(dx[i]);
      }
      dxvsx->Fill(x[i],dx[i]);
      dxvsy->Fill(xy[i],dx[i]);

      dxvsty->Fill(ty,dx[i]);
      dxvstx->Fill(tx,dx[i]);

      txvsx->Fill(yx[i],tx);
      txvsy->Fill(y[i],tx);
    }

  }


  // plot them ************************************************
  TCanvas *c0;
  TCanvas *c1;
  TCanvas *c2;

  c0 = new TCanvas("c0","dk vs j");
  c0->Divide(2,2);
  c0->cd(1);
  fitslices(dyvsx, 0)->Fit("pol1");
  c0->cd(2);
  fitslices(dyvsy, 0)->Fit("pol1");
  c0->cd(3);
  fitslices(dxvsx, 0)->Fit("pol1");   
  c0->cd(4);
  fitslices(dxvsy, 0)->Fit("pol1");  
  c0->Update();
 // c0->Print("tpcp0d00_"+oFile+".pdf");




  fFunk = getFitFunction();
  c1 = new TCanvas("c1","dy and dx");
  c1->Divide(2,3);
  //for(int k = 1; k<=6; k++) c1->cd(k)->SetGridy(1);
  c1->cd(1)->SetLogy(1);
  dy10->Fit(fFunk);
  c1->cd(2)->SetLogy(1); 
  dx10->Fit(fFunk);
  c1->cd(3)->SetLogy(1); 
  dy100->Fit(fFunk);
  c1->cd(4)->SetLogy(1); 
  dx100->Fit(fFunk);
  c1->cd(5)->SetLogy(1); 
  dy1000->Fit(fFunk);
  c1->cd(6)->SetLogy(1); 
  dx1000->Fit(fFunk);
  c1->Update();
  //c1->Print("tpcp0d04_"+oFile+".pdf");

  

  c2 = new TCanvas("c2","dx,dy vs tx,ty");
  c2->Divide(2,2);
  c2->cd(1);
  fitslices(dyvsty, 0)->Fit("pol1");
  c2->cd(2);
  fitslices(dyvstx, 0)->Fit("pol1");
  c2->cd(3);
  fitslices(dxvsty, 0)->Fit("pol1");
  c2->cd(4);
  fitslices(dxvstx, 0)->Fit("pol1");
  c2->Update(); 
  //c2->Print("tpcp0d09_"+oFile+".pdf");

  c3 = new TCanvas("c3", "dx and dy, 50 ns");
  cout << "drawing dx dy" << endl;
  chain->Draw("dx:yz>>h2x(20,-3500,-1000,20,-250,250)", "sqrt(pow(time-thitx,2)) < 50");
  TH2D* h2X = (TH2D*)gDirectory->Get("h2x")->Clone("h2X");
  cout << "drawing dx dy" << endl;
  chain->Draw("dy:yz>>h2y(20,-3500,-1000,20,-250,250)", "sqrt(pow(time-thity,2)) < 50");
  TH2D* h2Y = (TH2D*)gDirectory->Get("h2y")->Clone("h2Y");
  TH1D *hZY =  fitslices(h2Y);
  TH1D *hZX =  fitslices(h2X);
  c3->Clear();
  c3->Divide(2,2);
  c3->cd(1);
  dx50->Draw();
  p50x = getPaveText("p50x", dx50);
  p50x->Draw();
  //dx50->Fit(fFunk);
  c3->cd(2);
  dy50->Draw();
  p50y = getPaveText("p50y", dy50);
  p50y->Draw();
  c3->cd(3);
  h2X_1->SetTitle("dx vs. z");
  h2X_1->GetXaxis()->SetTitle("z [mm]");
  h2Y_1->GetXaxis()->SetTitle("z [mm]");
  h2X_1->GetYaxis()->SetTitle("dx [mm]");
  h2Y_1->GetYaxis()->SetTitle("dy [mm]");
  h2Y_1->SetTitle("dy vs. z");
  h2X_1->Draw();
  c3->cd(4);
  h2Y_1->Draw();
  c3->Update();

  

  cFinal = new TCanvas("c", "", 200,200);
/*
  TPaveText *pt = new TPaveText();
  pt->SetX1NDC(0.3);pt->SetX2NDC(0.7); pt->SetY1NDC(0.4);pt->SetY2NDC(0.6);
  pt->AddText("Click to Finish");
  pt->Draw();
  cFinal->Update();
  cFinal->WaitPrimitive();
*/

  c0->Print("djvsdk.pdf");
  c1->Print("dj.pdf");
  c2->Print("djvstk.pdf");
  




  delete c0;
  delete c1;
  delete c2;
  delete c3;
  delete cFinal;
  delete fFunk;
}


int TPC(int z){
  if( z < 100. ) return 1; 
  if( z < 1400. ) return 2; 
  return 3; 
}

TH1D  *fitslices(TH2 *hist, int i = 0){
  hist->FitSlicesY(); 
  char name[256];
  sprintf(name,"%s_1",hist->GetName()); 
  TH1D *hfit =  (TH1D*)gDirectory->Get(name);
  hfit->SetTitle(name);
  hfit->SetLineColor(i+1);
  //(hfit->GetYaxis())->SetRangeUser(-10.,10.); 
  return hfit;
}

TF1 *getFitFunction(){
  // fit function
  TF1 *fFunk = new TF1("fFunk", "[0]*(1.0/sqrt(2*3.14159*[2]*[2]))*exp(-0.5*pow((x-[1])/[2],2))", -5000,5000);
  fFunk->SetParLimits(2, 0, 50);  // set the width of the gaussian
  fFunk->SetParLimits(5, 100, 1e5);  // set the width of the gaussian
  fFunk->SetParameter(0, 10);
  fFunk->SetParameter(0, 100);
  fFunk->SetParameter(2, 20);
  fFunk->SetParName(0, "N");
  fFunk->SetParName(1, "#mu");
  fFunk->SetParName(2, "#sigma");
  fFunk->SetNpx(1000);
  return fFunk;

}

void formatHist(TH1 *h1){
  h1->GetXaxis()->SetTitle("[mm]");
  h1->SetFillColor(kGray+2);
}

TPaveText* getPaveText(string name, TH1 *h1){
  pt = new TPaveText();
  pt->SetName(name.c_str());
  pt->AddText(Form("Mean: %.2f",h1->GetMean()));
  pt->AddText(Form("RMS: %.2f",h1->GetRMS()));
  pt->SetX1NDC(0.5);
  pt->SetX2NDC(0.7);
  pt->SetY1NDC(0.5);
  pt->SetY2NDC(0.7);
  return pt;
}