#define analtof1_cxx
// The class definition in analtof1.h has been generated automatically
// by the ROOT utility TTree::MakeSelector(). This class is derived
// from the ROOT class TSelector. For more information on the TSelector
// framework see $ROOTSYS/README/README.SELECTOR or the ROOT User Manual.

// The following methods are defined in this file:
//    Begin():        called every time a loop on the tree starts,
//                    a convenient place to create your histograms.
//    SlaveBegin():   called after Begin(), when on PROOF called only on the
//                    slave servers.
//    Process():      called for each event, in this function you decide what
//                    to read and fill your histograms.
//    SlaveTerminate: called at the end of the loop on the tree, when on PROOF
//                    called only on the slave servers.
//    Terminate():    called at the end of the loop on the tree,
//                    a convenient place to draw/fit your histograms.
//
// To use this file, try the following session on your Tree T:
//
// Root > T->Process("analtof1.C")
// Root > T->Process("analtof1.C","some options")
// Root > T->Process("analtof1.C+")
//

#include "analtof1.h"
#include <TH2.h>
#include <TStyle.h>
#include <sstream>

int DEBUG = 0;

void analtof1::Begin(TTree * /*tree*/)
{
  // The Begin() function is called at the start of the query.
  // When running with PROOF Begin() is only called on the client.
  // The tree argument is deprecated (on PROOF 0 is passed).
  
  std::cout<<"Begin() Called!"<<endl;
  memset(Hit_Count,0,2*45*2*4);
  TRK_Count = 0;
  
  MCdata = 0;

  TString option = GetOption();
  stringstream ss(option.Data());
  ss>>RunNumber;
  cout<<"This is Run "<<RunNumber<<endl;
  ss>>MCdata;
  cout<<"Is This MC data? "<<MCdata<<endl;
   
  char fn[128];
  sprintf(fn,"EFF/hists_run%s.root",RunNumber.Data());
  if (MCdata){
    sprintf(fn,"EFF/hists_MC_run%s.root",RunNumber.Data());
  }
  OF = new TFile(fn,"RECREATE");
  
  
  zvertex = new TH1F("zvertex", "track z vertex position",140,30., 100.);
  zvertex->GetXaxis()->SetTitle("track vertex z [cm]");

  ntofpoints = new TH1F("ntofpoints","number of tof points",20, 0., 20.);

  char hnam[128];
  char htit[128];
  for (int k=0; k<45; k++) { // loop over paddles
    sprintf(hnam,"pvsdx%02d",k);
    sprintf(htit,"Momentum vs deltaX H-Plane Paddle%02d",k);
    pvsdx[k] = new TH2F(hnam,htit, 100, -50., 50.,40, 0., 10.);
    pvsdx[k]->GetXaxis()->SetTitle("Track_X-TOF_X [cm]");
    pvsdx[k]->GetYaxis()->SetTitle("Track momentum [GeV/c]");
    sprintf(hnam,"pvsdy%02d",k);
    sprintf(htit,"Momentum vs deltaY V-Plane Paddle%02d",k);
    pvsdy[k] = new TH2F(hnam,htit, 100, -50., 50.,40, 0., 10.); 
    pvsdy[k]->GetXaxis()->SetTitle("Track_Y-TOF_Y [cm]");
    pvsdy[k]->GetYaxis()->SetTitle("Track momentum [GeV/c]");
  }
  
  pvsTOFdx = new TH2F("pvsTOFdx","Momentum vs Nearest deltaX", 100, -50., 50.,40, 0., 10.);
  pvsTOFdx->GetXaxis()->SetTitle("Track_X-TOF_X");
  pvsTOFdx->GetYaxis()->SetTitle("Track momentum [GeV/c]");
  pvsTOFdy = new TH2F("pvsTOFdy","Momentum vs Nearest deltaY", 100, -50., 50.,40, 0., 10.);
  pvsTOFdy->GetXaxis()->SetTitle("Track_Y-TOF_Y");
  pvsTOFdy->GetYaxis()->SetTitle("Track momentum [GeV/c]");

  for (int k=0;k<20;k++){
    sprintf(hnam,"proj_vs_nearest_X%d",k);
    sprintf(htit,"Projected vs Nearest hit horizontal plane p = %4.1f GeV",k*0.5+0.5);    
    proj_vs_nearest_H[k] = new TH2F(hnam, htit, 45, 0., 45., 45, 0., 45.);
    proj_vs_nearest_H[k]->GetXaxis()->SetTitle("Horizontal Paddle with Hit [#]");
    proj_vs_nearest_H[k]->GetYaxis()->SetTitle("Projected Horizontal Paddle [#]");

    sprintf(hnam,"proj_vs_nearest_X4goodV%d",k);
    sprintf(htit,"Projected vs Nearest hit horizontal (given good Vhit) plane p = %4.1f GeV",k*0.5+0.5);    
    proj_vs_nearest_H4goodV[k] = new TH2F(hnam, htit, 45, 0., 45., 45, 0., 45.);
    proj_vs_nearest_H4goodV[k]->GetXaxis()->SetTitle("Horizontal Paddle with Hit [#]");
    proj_vs_nearest_H4goodV[k]->GetYaxis()->SetTitle("Projected Horizontal Paddle [#]");


    sprintf(hnam,"proj_vs_nearest_Y%d",k);
    sprintf(htit,"Projected vs Nearest hit vertical plane p = %4.1f GeV",k*0.5+0.5);    
    proj_vs_nearest_V[k] = new TH2F(hnam, htit, 45, 0., 45., 45, 0., 45.);
    proj_vs_nearest_V[k]->GetXaxis()->SetTitle("Horizontal Paddle with Hit [#]");
    proj_vs_nearest_V[k]->GetYaxis()->SetTitle("Projected Horizontal Paddle [#]");

    sprintf(hnam,"proj_vs_nearest_Y4goodH%d",k);
    sprintf(htit,"Projected vs Nearest hit vertical (given good Hhit) plane p = %4.1f GeV",k*0.5+0.5);    
    proj_vs_nearest_V4goodH[k] = new TH2F(hnam, htit, 45, 0., 45., 45, 0., 45.);
    proj_vs_nearest_V4goodH[k]->GetXaxis()->SetTitle("Horizontal Paddle with Hit [#]");
    proj_vs_nearest_V4goodH[k]->GetYaxis()->SetTitle("Projected Horizontal Paddle [#]");
  }
  
  xy = new TH2F("xy", "Hits in xy plane", 82, -126., 126., 82, -126., 126.);
  for (int j=0; j<84;j++){
    for (int k=0; k<84;k++){
      sprintf(hnam,"xy_p_vs_dx%02d%02d",j,k);
      sprintf(htit,"p vs dx bined-xy xbin%02d ybin%02d",j, k);    
      xy_p_vs_dx[j][k] = new TH2F(hnam, htit, 100, -10, 10., 50, 0., 10.);
 
      sprintf(hnam,"xy_p_vs_dy%02d%02d",j,k);
      sprintf(htit,"p vs dy bined-xy xbin%02d ybin%02d",j, k);    
      xy_p_vs_dy[j][k] = new TH2F(hnam, htit, 100, -10, 10., 50, 0., 10.);

      sprintf(hnam,"xyQ_H_vs_V%02d%02d",j,k);
      sprintf(htit,"BarStatus Horizontal vs Vertical xbin%02d ybin%02d",j, k);    
      xyQ_H_vs_V[j][k] = new TH2F(hnam, htit, 5, -0.5, 4.5, 5, -0.5, 4.5 ); 


    }
  }

}

void analtof1::SlaveBegin(TTree * /*tree*/)
{
   // The SlaveBegin() function is called after the Begin() function.
   // When running with PROOF SlaveBegin() is called on each slave server.
   // The tree argument is deprecated (on PROOF 0 is passed).

   std::cout<<"SlaveBegin() Called!"<<endl;

   TString option = GetOption();

}

Bool_t analtof1::Process(Long64_t entry)
{
  // The Process() function is called for each entry in the tree (or possibly
  // keyed object in the case of PROOF) to be processed. The entry argument
  // specifies which entry in the currently loaded tree is to be processed.
  // It can be passed to either analtof1::GetEntry() or TBranch::GetEntry()
  // to read either all or the required parts of the data. When processing
  // keyed objects with PROOF, the object is already loaded and is available
  // via the fObject pointer.
  //
  // This function should contain the "body" of the analysis. It can contain
  // simple or elaborate selection criteria, run algorithms on the data
  // of the event and typically fill histograms.
  //
  // The processing can be stopped by calling Abort().
  //
  // Use fStatus to set the return value of TTree::Process().
  //
  // The return value is currently not used.
  
  GetEntry(entry);
  //cout<<entry<<endl;

  zvertex->Fill(TrackVertexZ);
  ntofpoints->Fill((float)((unsigned int)NumTOFPoints));

  double p = TrackP3->Mag();
  if ((p>12.) || (p<0.5)){
    return kTRUE;
  }

  float R = 0;
  R = TMath::Sqrt(ProjectedTOFY*ProjectedTOFY + ProjectedTOFX*ProjectedTOFX);
  //cout<<R<<endl;
  //if (R>125.){
  //  return kTRUE;
  //}

  unsigned int FDCPlane4Hits = (TrackFDCPlanes>>18) & 0x3F;
  unsigned int FDCPlane3Hits = (TrackFDCPlanes>>12) & 0x3F;
  unsigned int FDCPlane2Hits = (TrackFDCPlanes>>6) & 0x3F;
  unsigned int FDCPlane1Hits = (TrackFDCPlanes) & 0x3F;
  if (!(FDCPlane4Hits & 0x38)){ // make sure the 4th FDC has at least one hit in the most downstream 3 planes
    return kTRUE;
  }

  if ((FDCPlane3Hits == 0) || (FDCPlane2Hits == 0) || (FDCPlane1Hits == 0) ){
    return kTRUE;
  }

  // horizontal plane hit 
  float paddleH = (float)((unsigned int)NearestTOFHit_Horizontal);
  float paddleHproj = (float)((unsigned int)ProjectedTOFBar_Horizontal);
  // 201-244 means single ended paddle hit South/bottom
  // 101-144 means single ended paddle hit North/top
  if (paddleH>200){
    paddleH -=200;
  } else if (paddleH>100){
    paddleH -=100;
  }
  
  float paddleV = (float)((unsigned int)NearestTOFHit_Vertical);
  float paddleVproj = (float)((unsigned int)ProjectedTOFBar_Vertical);
  if (paddleV>200){
    paddleV -=200;
  } else if (paddleV>100){
    paddleV -=100;
  }
  
  int pbin = (int)(p*2.);
  proj_vs_nearest_H[pbin]->Fill(paddleH, paddleHproj);
  proj_vs_nearest_V[pbin]->Fill(paddleV, paddleVproj);


  if (((unsigned int)NearestTOFHit_Vertical) == ((unsigned int)ProjectedTOFBar_Vertical)){
    proj_vs_nearest_H4goodV[pbin]->Fill(paddleH, paddleHproj);
  }
  if  (((unsigned int)NearestTOFHit_Horizontal) == ((unsigned int)ProjectedTOFBar_Horizontal)){
    proj_vs_nearest_V4goodH[pbin]->Fill(paddleV, paddleVproj);
  }


  if (DEBUG){
    if (paddleHproj<1){
      cout<<"H-projection failed:"<<endl;
      cout<<ProjectedTOFX<<" / "<<ProjectedTOFY<<endl;
    }
    if (paddleVproj<1){
      cout<<"V-projection failed:"<<endl;
      cout<<ProjectedTOFX<<" / "<<ProjectedTOFY<<endl;
    }
  }

  xy->Fill(ProjectedTOFX, ProjectedTOFY);
  int GlobalBinNumber = xy->FindBin(ProjectedTOFX, ProjectedTOFY);

  int binx = (ProjectedTOFX+123.)/3.;
  int biny = (ProjectedTOFY+123.)/3.;
  int binz = 0;

  xy->GetBinXYZ(GlobalBinNumber,binx,biny,binz);
  binx -=1;
  biny -=1;
  if (binx<0)
    return kTRUE;
  if (binx>83)
    return kTRUE;
  if (biny<0)
    return kTRUE;
  if (biny>83)
    return kTRUE;

  //cout<<"x/y = "<<ProjectedTOFX<<"/"<<ProjectedTOFY<<"           "<<binx<<" / "<<biny<<endl;

  xy_p_vs_dx[binx][biny]->Fill(NearestTOFPointDeltaX, p);
  xy_p_vs_dy[binx][biny]->Fill(NearestTOFPointDeltaY, p);

  UShort_t X =NearestTOFPointStatus;

  int VBStatus = 3;
  if ((X-8100*VBStatus)<0){
    VBStatus -= 1;
    if ((X-8100*VBStatus)<0){
      VBStatus -= 1;
    }
    if ((X-8100*VBStatus)<0){
      VBStatus -= 1;
    }
  }
  X -= 8100*VBStatus;

  int HBStatus = 3;
  if ((X-2025*HBStatus)<0){
    HBStatus -= 1;
    if ((X-2025*HBStatus)<0){
      HBStatus -= 1;
    }
    if ((X-2025*HBStatus)<0){
      HBStatus -= 1;
    }
  }
  X -= 2025*HBStatus;

  int VBar = 44;
  for (int k=44; k>0; k--){
    if (X-(45*k)<0){
      continue;
    } else {
      VBar = k;
      break;
    }
  }
  X -= 45*VBar;

  int HBar = X;

  //cout<<HBar<<"  "<<VBar<<"   "<<HBStatus<<"   "<<VBStatus<<endl;

  xyQ_H_vs_V[binx][biny]->Fill(VBStatus, HBStatus);

  pvsTOFdx->Fill(NearestTOFPointDeltaX,p); 
  pvsTOFdy->Fill(NearestTOFPointDeltaY,p); 

  // Vertical Plane
  pvsdy[(int)paddleV]->Fill(NearestTOFPointDeltaY,p);
  // Horizontal Plane
  pvsdx[(int)paddleH]->Fill(NearestTOFPointDeltaX,p);
  
  TRK_Count += 1.;
  
  if (TRK_Count<1) {
    cout<<ProjectedTOFX<<"  "<<ProjectedTOFY<<endl;
    cout<<(unsigned int)ProjectedTOFBar_Horizontal<<" / "<<(unsigned int)ProjectedTOFBar_Vertical<<"   "<<IsMatchedToTrack<<endl;
  }
    
  Hit_Count[1][(unsigned int)ProjectedTOFBar_Horizontal][0] +=1;
  Hit_Count[0][ProjectedTOFBar_Vertical][0] +=1;

  unsigned int NH = (unsigned int)NearestTOFHit_Horizontal;
  
  if (NH){
    if ((ProjectedTOFX<0) && NH>200){
      NH -= 200;
    }
    if ((ProjectedTOFX>0) && NH>100){
      NH -= 100;
    }
    
    //    if (((unsigned int)ProjectedTOFBar_Horizontal) == NH){
    int a = ((unsigned int)ProjectedTOFBar_Horizontal) - NH;
    if (TMath::Abs(a)<2){
      Hit_Count[1][(unsigned int)ProjectedTOFBar_Horizontal][1] +=1;
    }
  }

  unsigned int NV = (unsigned int)NearestTOFHit_Vertical;
  if (NV){
    if ((ProjectedTOFY<0) && NV>200){
      NV -= 200;
    }
    if ((ProjectedTOFY>0) && NV>100){
      NV -= 100;
    }
    
    
    //if (((unsigned int)ProjectedTOFBar_Vertical) == NV){
    int a = ((unsigned int)ProjectedTOFBar_Vertical) - NV;
    if (TMath::Abs(a)<2){
      Hit_Count[0][(unsigned int)ProjectedTOFBar_Vertical][1] +=1;
    }
  }
  
  return kTRUE;
}

void analtof1::SlaveTerminate()
{
  // The SlaveTerminate() function is called after all entries or objects
  // have been processed. When running with PROOF SlaveTerminate() is called
  // on each slave server.
  std::cout<<"SlaveTerminate() Called!"<<endl;
  
}

void analtof1::Terminate()
{
  // The Terminate() function is the last function to be called during
  // a query. It always runs on the client, it can be used to present
  // the results graphically or save the results to file.
  
  //  std::cout<< "Terminate() Called!" << endl;
  // std:cout<<"Num of Trks: "<<TRK_Count<<endl;
  
  ofstream of;
  char ofnam[128];
  sprintf(ofnam,"EFF/paddle_efficiency_run%s.dat",RunNumber.Data());
  if (MCdata){
    sprintf(ofnam,"EFF/paddle_efficiency_MC_run%s.dat",RunNumber.Data());
  }
  of.open(ofnam);
  of<< TRK_Count<<endl;
  for (int k=1;k<45;k++){
    float r1 = (float)Hit_Count[0][k][1]/(float)Hit_Count[0][k][0];
    float r2 = (float)Hit_Count[1][k][1]/(float)Hit_Count[1][k][0];
    
    float dd = ((float)Hit_Count[0][k][0]-(float)Hit_Count[0][k][1]) + r1*r1/(float)Hit_Count[0][k][0];
    float dr1 = TMath::Sqrt(dd) / (float)Hit_Count[0][k][0];
    
    dd = ((float)Hit_Count[1][k][0] - (float)Hit_Count[1][k][1]) + r2*r2/(float)Hit_Count[1][k][0];
    float dr2 = TMath::Sqrt(dd) / (float)Hit_Count[1][k][0];    
    
    of<<k<<"  "<<r1<<"  "<<dr1<<"   "<<r2<<"  "<<dr2<<endl;
  }
  of.close();

  
  zvertex->Write();
  ntofpoints->Write();
  for (int k=0; k<44; k++) { // loop over paddles
    pvsdx[k]->Write();
    pvsdy[k]->Write();
  }
  

  pvsTOFdx->Write();
  pvsTOFdy->Write();

  for (int k=0;k<20;k++){
    proj_vs_nearest_H[k]->Write();
    proj_vs_nearest_H4goodV[k]->Write();
    proj_vs_nearest_V[k]->Write();
    proj_vs_nearest_V4goodH[k]->Write();
  }

  xy->Write();
  for (int j=0;j<84;j++){
    for (int k=0;k<84;k++){
      xy_p_vs_dx[j][k]->Write();
      xy_p_vs_dy[j][k]->Write();
      xyQ_H_vs_V[j][k]->Write();
    }
  }


  OF->Close();

}