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

void tpcfgd(void){

  gStyle->SetOptFit(11111111); 
TFile* f = new TFile("1kDSTrackerRecoOut.root", "recreate");

  TH2F *dyvsy1 = new TH2F("dyvsy1","  ",20,-1200.,1200.,20,-20.,20.); 
  TH2F *dyvsy2 = new TH2F("dyvsy2","  ",20,-1200.,1200.,20,-20.,20.); 
  TH2F *dyvsx1 = new TH2F("dyvsx1","  ",20,-1200.,1200.,20,-20.,20.); 
  TH2F *dyvsx2 = new TH2F("dyvsx2","  ",20,-1200.,1200.,20,-20.,20.); 
  
  TH2F *dxvsx = new TH2F("dxvsx","  ",10,-900.,900.,20,-20.,20.); 
  TH2F *dxvsy = new TH2F("dxvsy","  ",10,-900.,900.,20,-20.,20.); 
  TH2F *dyvsx = new TH2F("dyvsx","  ",10,-900.,900.,20,-20.,20.); 
  TH2F *dyvsy = new TH2F("dyvsy","  ",10,-900.,900.,20,-20.,20.); 

  TH2F *dxvsx1 = new TH2F("dxvsx1","  ",20,-1200.,1200.,20,-20.,20.); 
  TH2F *dxvsx2 = new TH2F("dxvsx2","  ",20,-1200.,1200.,20,-20.,20.);
  TH2F *dxvsy1 = new TH2F("dxvsy1","  ",20,-1200.,1200.,20,-20.,20.); 
  TH2F *dxvsy2 = new TH2F("dxvsy2","  ",20,-1200.,1200.,20,-20.,20.);

  TH2F *tdyvsy1 = new TH2F("tdyvsy1","  ",20,-1200.,1200.,20,-20.,20.); 
  TH2F *tdyvsy2 = new TH2F("tdyvsy2","  ",20,-1200.,1200.,20,-20.,20.); 
  TH2F *tdyvsy3 = new TH2F("tdyvsy3","  ",20,-1200.,1200.,20,-20.,20.); 
  TH2F *tdyvsx1 = new TH2F("tdyvsx1","  ",20,-1200.,1200.,20,-20.,20.); 
  TH2F *tdyvsx2 = new TH2F("tdyvsx2","  ",20,-1200.,1200.,20,-20.,20.); 
  TH2F *tdyvsx3 = new TH2F("tdyvsx3","  ",20,-1200.,1200.,20,-20.,20.); 
  
  TH2F *tdxvsx1 = new TH2F("tdxvsx1","  ",20,-1200.,1200.,20,-20.,20.); 
  TH2F *tdxvsx2 = new TH2F("tdxvsx2","  ",20,-1200.,1200.,20,-20.,20.);
  TH2F *tdxvsx3 = new TH2F("tdxvsx3","  ",20,-1200.,1200.,20,-20.,20.);
  TH2F *tdxvsy1 = new TH2F("tdxvsy1","  ",20,-1200.,1200.,20,-20.,20.); 
  TH2F *tdxvsy2 = new TH2F("tdxvsy2","  ",20,-1200.,1200.,20,-20.,20.);
  TH2F *tdxvsy3 = new TH2F("tdxvsy3","  ",20,-1200.,1200.,20,-20.,20.);

  TH1F *dyT = new TH1F("dyT","  ",100.,-100.,100.);
  TH1F *dxT = new TH1F("dxT","  ",100.,-100.,100.);

#define ANG 1.0
#define BIN 10

  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.);


  TH2F *dyvstyF1 = new TH2F("dyvstyF1","  ",BIN,-ANG,ANG,20,-20.,20.);
  TH2F *dyvstyF2 = new TH2F("dyvstyF2","  ",BIN,-ANG,ANG,20,-20.,20.);
  TH2F *dyvstyL1 = new TH2F("dyvstyL1","  ",BIN,-ANG,ANG,20,-20.,20.);
  TH2F *dyvstyL2 = new TH2F("dyvstyL2","  ",BIN,-ANG,ANG,20,-20.,20.);


  TH2F *dxvstxF1 = new TH2F("dxvstxF1","  ",BIN,-ANG,ANG,20,-20.,20.);
  TH2F *dxvstxF2 = new TH2F("dxvstxF2","  ",BIN,-ANG,ANG,20,-20.,20.);
  TH2F *dxvstxL1 = new TH2F("dxvstxL1","  ",BIN,-ANG,ANG,20,-20.,20.);
  TH2F *dxvstxL2 = new TH2F("dxvstxL2","  ",BIN,-ANG,ANG,20,-20.,20.);


  TH1F *dxN = new TH1F("dxN","  ",100.,-100.,100.);
  TH1F *dxP = new TH1F("dxP","  ",100.,-100.,100.);

  TH1F *dy1 = new TH1F("dy1","  ",100.,-100.,100.);
  TH1F *dy2 = new TH1F("dy2","  ",100.,-100.,100.);

  TH1F *dx1 = new TH1F("dx1","  ",100.,-100.,100.);
  TH1F *dx2 = new TH1F("dx2","  ",100.,-100.,100.);

  TH1F *tdy1 = new TH1F("tdy1","  ",100.,-100.,100.);
  TH1F *tdy2 = new TH1F("tdy2","  ",100.,-100.,100.);
  TH1F *tdy3 = new TH1F("tdy3","  ",100.,-100.,100.);

  TH1F *tdx1 = new TH1F("tdx1","  ",100.,-100.,100.);
  TH1F *tdx2 = new TH1F("tdx2","  ",100.,-100.,100.);
  TH1F *tdx3 = new TH1F("tdx3","  ",100.,-100.,100.);

  TCanvas *c0 = new TCanvas("c0"," Dx,y vs x,y  2 FGD   ");
  TCanvas *c1 = new TCanvas("c1"," Dx,y 2 FGD   ");
  TCanvas *c2 = new TCanvas("c2"," Dx,y vs x,y  3 TPC  ");
  TCanvas *c3 = new TCanvas("c3"," Dx,y 3 TPC  ");
  TCanvas *c4 = new TCanvas("c4"," Dx,y  ");
  TCanvas *c6 = new TCanvas("c6"," Dx,y vs tx,ty  ");
  TCanvas *c7 = new TCanvas("c7"," Dx vs tx for different layers  ");
  TCanvas *c8 = new TCanvas("c8"," Dy vs ty for different layers  ");

  TChain *chain = new TChain("tpcfgdtree");

  
      if( FileExists("../cmt/1kDSTrackerReco.root") ) {
  chain->AddFile("../cmt/1kDSTrackerReco.root",0);
}	
  chain->GetEntry(0);

  double xtrack; 
  double ytrack; 
  double ztrack;
  double tx;
  double ty; 
  double rho;
  double time; 
  int nyhits; 
  double y[100];
  double yz[100];
  double yx[100];
  double dy[100];
  double qy[100];
  int nxhits; 
  double x[100];
  double xz[100];
  double xy[100];
  double dx[100];
  double qx[100];

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

  chain->SetBranchAddress("xtrack",&xtrack);
  chain->SetBranchAddress("ytrack",&ytrack);
  chain->SetBranchAddress("ztrack",&ztrack);
  chain->SetBranchAddress("txtrack",&tx);
  chain->SetBranchAddress("tytrack",&ty);
  chain->SetBranchAddress("rho",&rho);
  chain->SetBranchAddress("time",&time);
  chain->SetBranchAddress("nyhit",&nyhits);
  chain->SetBranchAddress("y",y);
  chain->SetBranchAddress("yx",yx);
  chain->SetBranchAddress("yz",yz);
  chain->SetBranchAddress("dy",dy);
  chain->SetBranchAddress("nxhit",&nxhits);
  chain->SetBranchAddress("x",x);
  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++){
    
    if( iev == 5017 ) continue;
    if( iev == 5018 ) continue; 
    if( iev == 17958 ) continue; 
    if( iev == 17959 ) continue; 
    if( iev == 20630 ) continue; 
    if( iev == 20631 ) continue; 

    chain->GetEntry(iev);

    if( TMath::Abs(time) < 0.01 ) continue; 
     
    int tpc = TPC(ztrack);

    for(int i = 0; i < nyhits; i++ ){

      int fgd =  FGD(yz[i]);

      dyT->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]);
  
      if( fgd == -1 && tpc == 1) dyvstyF1->Fill(ty,dy[i]);
      if( fgd == 1 && tpc == 2 ) dyvstyL1->Fill(ty,dy[i]);
      if( fgd == -2 && tpc == 2) dyvstyF2->Fill(ty,dy[i]);
      if( fgd == 2 && tpc == 3) dyvstyL2->Fill(ty,dy[i]);

      if( abs(fgd) == 1 ) {
  dyvsy1->Fill(y[i],dy[i]);
  dyvsx1->Fill(yx[i],dy[i]);
  dy1->Fill(dy[i]);
      }
      else if ( abs(fgd) == 2 ){
  dyvsy2->Fill(y[i],dy[i]);
  dyvsx2->Fill(yx[i],dy[i]);
  dy2->Fill(dy[i]);
      }

      if( tpc == 1 ) {
  tdyvsy1->Fill(y[i],dy[i]);
  tdyvsx1->Fill(yx[i],dy[i]);
  tdy1->Fill(dy[i]);
      }
      else if ( tpc == 2 ){
  tdyvsy2->Fill(y[i],dy[i]);
  tdyvsx2->Fill(yx[i],dy[i]);
  tdy2->Fill(dy[i]);
      }
      else if ( tpc == 3 ){
  tdyvsy3->Fill(y[i],dy[i]);
  tdyvsx3->Fill(yx[i],dy[i]);
  tdy3->Fill(dy[i]);
      }

    }

    for(int i = 0; i < nxhits; i++ ){

      int fgd =  FGD(xz[i]);

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

      if( fgd == -1 && tpc == 1 ) dxvstxF1->Fill(tx,dx[i]);
      if( fgd == 1 && tpc == 2 ) dxvstxL1->Fill(tx,dx[i]);
      if( fgd == -2 && tpc == 2 ) dxvstxF2->Fill(tx,dx[i]);
      if( fgd == 2 && tpc == 3 ) dxvstxL2->Fill(tx,dx[i]);

      if( x[i] < 0 ) 
  dxN->Fill(dx[i]);
      else 
  dxP->Fill(dx[i]);

      dxvsx->Fill(x[i],dx[i]);
      dxvsy->Fill(xy[i],dx[i]);

      if( abs(fgd) == 1 ) {
  dxvsx1->Fill(x[i],dx[i]);
  dxvsy1->Fill(xy[i],dx[i]);
  dx1->Fill(dx[i]);
      }
      else if ( abs(fgd) == 2 ){
  dxvsx2->Fill(x[i],dx[i]);
  dxvsy2->Fill(xy[i],dx[i]);
  dx2->Fill(dx[i]);
      }

      if( tpc == 1 ) {
  tdxvsx1->Fill(x[i],dx[i]);
  tdxvsy1->Fill(xy[i],dx[i]);
  tdx1->Fill(dx[i]);
      }
      else if ( tpc == 2 ){
  tdxvsx2->Fill(x[i],dx[i]);
  tdxvsy2->Fill(xy[i],dx[i]);
  tdx2->Fill(dx[i]);
      }
      else if ( tpc == 3 ){
  tdxvsx3->Fill(x[i],dx[i]);
  tdxvsy3->Fill(xy[i],dx[i]);
  tdx3->Fill(dx[i]);
      }

    }
        
  }

  c0->Divide(2,4); 
  c0->cd(1);
  fitslices(dyvsx1)->Draw();
  
  c0->cd(2);
  fitslices(dyvsx2)->Draw(); 

  c0->cd(3);
  fitslices(dyvsy1)->Draw();
  
  c0->cd(4);
  fitslices(dyvsy2)->Draw(); 

  c0->cd(5);
  fitslices(dxvsx1)->Draw();
  
  c0->cd(6);
  fitslices(dxvsx2)->Draw(); 

  c0->cd(7);
  fitslices(dxvsy1)->Draw();
  
  c0->cd(8);
  fitslices(dxvsy2)->Draw(); 

  c0->Update();

   TF1 *f1 = new TF1("f1", "gaus", -10, 10);
  c1->Divide(2,2); 
  c1->cd(1);
  dy1->Fit("f1", "R");
  dy1->Draw();
  
  c1->cd(2);
  dy2->Fit("f1", "R");
  dy2->Draw();

  c1->cd(3);
  dx1->Fit("f1", "R");
  dx1->Draw();
  
  c1->cd(4);
  dx2->Fit("f1", "R");
  dx2->Draw(); 

  c1->Update();


  c2->Divide(3,4); 
  c2->cd(1);
  fitslices(tdyvsy1)->Draw();
  
  c2->cd(2);
  fitslices(tdyvsy2)->Draw(); 

  c2->cd(3);
  fitslices(tdyvsy3)->Draw(); 


  c2->cd(4);
  fitslices(tdyvsx1)->Draw();
  
  c2->cd(5);
  fitslices(tdyvsx2)->Draw(); 

  c2->cd(6); 
  fitslices(tdyvsx3)->Draw(); 


  c2->cd(7);
  fitslices(tdxvsx1)->Draw();
  
  c2->cd(8);
  fitslices(tdxvsx2)->Draw(); 
  
  c2->cd(9);
  fitslices(tdxvsx3)->Draw(); 

  c2->cd(10);
  fitslices(tdxvsy1)->Draw();
  
  c2->cd(11);
  fitslices(tdxvsy2)->Draw(); 

  c2->cd(12);
  fitslices(tdxvsy3)->Draw(); 

  c2->Update();


  c3->Divide(3,2); 

  c3->cd(1);
  tdy1->Fit("f1", "R");
  tdy1->Draw();
  
  c3->cd(2);
  tdy2->Fit("f1", "R");
  tdy2->Draw();
  
  c3->cd(3); 
  tdy3->Fit("f1", "R");
  tdy3->Draw();
  
  c3->cd(4);
  tdx1->Fit("f1", "R");
  tdx1->Draw();
  
  c3->cd(5);
  tdx2->Fit("f1", "R");
  tdx2->Draw(); 
  
  c3->cd(6); 
  tdx3->Fit("f1", "R");
  tdx3->Draw(); 


  c3->Update();

  c4->Divide(2,2); 
  c4->cd(1); 
  dyT->Fit("gaus");
  c4->cd(2); 
  dxT->Fit("gaus");
  c4->cd(3); 
  dxN->Fit("gaus");
  c4->cd(4); 
  dxP->Fit("gaus");

  c4->Update();


  c6->Divide(2,4);

  c6->cd(1);
  fitslices(dyvsty)->Fit("pol1");

  c6->cd(2);
  fitslices(dyvstx)->Fit("pol1");

  c6->cd(3);
  fitslices(dxvsty)->Fit("pol1");

  c6->cd(4);
  fitslices(dxvstx)->Fit("pol1");
  c6->cd(5);

  fitslices(dxvsx)->Draw();
  c6->cd(6);

  fitslices(dxvsy)->Draw();
  c6->cd(7);

  fitslices(dyvsx)->Draw();
  c6->cd(8);

  fitslices(dyvsy)->Draw();

  c6->Update(); 


  c7->Divide(2,2); 

  c7->cd(1); 
  fitslices(dxvstxF1)->Fit("pol1");

  c7->cd(2); 
  fitslices(dxvstxF2)->Fit("pol1");

  c7->cd(3); 
  fitslices(dxvstxL1)->Fit("pol1");

  c7->cd(4); 
  fitslices(dxvstxL2)->Fit("pol1");

  c7->Update(); 


  c8->Divide(2,2); 

  c8->cd(1); 
  fitslices(dyvstyF1)->Fit("pol1");

  c8->cd(2); 
  fitslices(dyvstyF2)->Fit("pol1");

  c8->cd(3); 
  fitslices(dyvstyL1)->Fit("pol1");

  c8->cd(4); 
  fitslices(dyvstyL2)->Fit("pol1");

  c8->Update(); 
  
  /*
  dyvsy1->Write();
dyvsy2->Write();
  dyvsx1->Write();
  dyvsx2->Write();
  
  dxvsx->Write();
  dxvsy->Write();
  dyvsx->Write();
  dyvsy ->Write();

  dxvsx1->Write();
  dxvsx2->Write();
  dxvsy1->Write();
  dxvsy2->Write();

  tdyvsy1->Write();
  tdyvsy2->Write();
  tdyvsy3->Write();
  tdyvsx1->Write();
  tdyvsx2->Write();
  tdyvsx3->Write();
  
  tdxvsx1->Write();
  tdxvsx2->Write();
  tdxvsx3->Write();
  tdxvsy1->Write();
  tdxvsy2->Write();
  tdxvsy3 ->Write();

  dyT->Write();
  dxT->Write();


  dyvstx ->Write();
  dyvsty ->Write();
  dxvstx ->Write();
  dxvsty ->Write();


  dyvstyF1->Write();
  dyvstyF2 ->Write();
  dyvstyL1 ->Write();
  dyvstyL2->Write();


  dxvstxF1->Write();
  dxvstxF2->Write();
  dxvstxL1->Write();
  dxvstxL2->Write();


  dxN->Write();
  dxP->Write();

  dy1 ->Write();
  dy2->Write();

  dx1->Write();
  dx2->Write();

  tdy1->Write();
  tdy2->Write();
  tdy3->Write();

  tdx1->Write();
  tdx2->Write();
  tdx3->Write();
  
*/
  f->Write();
}


int FGD(int z){

  if( z > 100. && z < 300. ) return -1;  // begin fgd 1 

  if( z > 300.  && z < 500 ) return 1;  // end fgd 1

  if( z > 1400. && z < 1600 ) return -2; // begin fgd 2 

  if( z > 1600. && z < 1900 ) return 2; // end fgd 2 

}

int TPC(int z){

  if( z < 100. ) return 1; 

  if( z < 1400. ) return 2; 

  return 3; 

}

TH1D  *fitslices(TH2F *hist){

  hist->FitSlicesY(); 

  char name[256];

  sprintf(name,"%s_1",hist->GetName()); 

  
  TH1D *hfit =  (TH1D*)gDirectory->Get(name);
  
  (hfit->GetYaxis())->SetRangeUser(-10.,10.); 

  return hfit;

}

bool FileExists(const char *name){
  struct stat stFileInfo;
  int intStat;

  // Attempt to get the file attributes
  intStat = stat(name,&stFileInfo); 
  if(intStat == 0) 
    return true;
  else{
    return false;
  }
}