//----------------------------------------
//
// macro to plot Beamline validation plots
//
// TODO: save information into TTree
// author: lysak@fnal.gov
// updates: 22 Jan 2009: Willis
//----------------------------------------

int create_beamline_plots(int MinRun=252836,int MaxRun=252910,char *Pass="17A")
{
  //-------------------------------
  //LOAD necesary modules, macros
  //-------------------------------

  TString outdir( gSystem->Getenv( "OUTPUTDIR"  ) );
  TString reldir( gSystem->Getenv( "RELEASEDIR" ) );


  cout << "-----------------------------------------------------------" << endl;
  cout << "Release directory: " << reldir << endl;
  cout << "Output  directory: " << outdir << endl;
  cout << "Database beamlines for period of runs: " 
       << MinRun << " - " << MaxRun << endl;
  cout << "-----------------------------------------------------------" << endl;

  fill_beamlines(MinRun,MaxRun,Pass);

  gSystem->ChangeDirectory( outdir );


  //-----------------------------------------------------------------
  // fill the arrays for ploting 
  // select only values which are fitted (fflag0=20) and 
  // not interpolated (fflag0=10) for a given detector
  //-----------------------------------------------------------------
  double *cotBeamX         = new double[npt];
  double *cotBeamY         = new double[npt];
  double *cotBeamXErr      = new double[npt];
  double *cotBeamYErr      = new double[npt];
  double *cotSlopeX        = new double[npt];
  double *cotSlopeY        = new double[npt];
  double *cotSlopeXErr     = new double[npt];
  double *cotSlopeYErr     = new double[npt];
  double *svxBeamX         = new double[npt];
  double *svxBeamY         = new double[npt];
  double *svxBeamXErr      = new double[npt];
  double *svxBeamYErr      = new double[npt];
  double *svxSlopeX        = new double[npt];
  double *svxSlopeY        = new double[npt];
  double *svxSlopeXErr     = new double[npt];
  double *svxSlopeYErr     = new double[npt];
  double *svxRun           = new double[npt];
  double *cotRun           = new double[npt];
  //number of runs with the fitted beamlines for SVX/COT
  int n_svx = 0;
  int n_cot = 0;
  
  cout << "-------------------------------------------" << endl;
  cout << "getting values for plotting                " << endl;
  cout << "-------------------------------------------" << endl;
  int svx_flag,cot_flag; 
  for (Int_t i=0;i<npt;i++)
    {
      svx_flag = svx_flag0[i] >> 16; 
      if ( (svx_flag==20) && (svx_status_bit[i]==1) && (cot_status_bit[i]==1) ) 
	{
	  svxBeamX[n_svx]    =svx_x[i];
	  svxBeamXErr[n_svx] =sqrt(svx_x_err[i]);
	  svxBeamY[n_svx]    =svx_y[i];
	  svxBeamYErr[n_svx] =sqrt(svx_y_err[i]);
	  svxSlopeX[n_svx]   =svx_slope_x[i];
	  svxSlopeXErr[n_svx]=sqrt(svx_slope_x_err[i]);
	  svxSlopeY[n_svx]   =svx_slope_y[i];
	  svxSlopeYErr[n_svx]=sqrt(svx_slope_y_err[i]);
	  svxRun[n_svx]      =run[i];
	  if ( (svx_y[i] == 0.16) || (svx_y[i] == 0.17) )
	    std::cout<<"Dummy value for run "<<run<<std::endl;
	  n_svx++;
	}
      
      cot_flag = cot_flag0[i] >> 16;
      if ( (cot_flag==20) && (cot_status_bit[i]==1) )
	{
	  cotBeamX[n_cot]    =cot_x[i];
	  cotBeamXErr[n_cot] =sqrt(cot_x_err[i]);
	  cotBeamY[n_cot]    =cot_y[i];
	  cotBeamYErr[n_cot] =sqrt(cot_y_err[i]);
	  cotSlopeX[n_cot]   =cot_slope_x[i];
	  cotSlopeXErr[n_cot]=sqrt(cot_slope_x_err[i]);
	  cotSlopeY[n_cot]   =cot_slope_y[i];
	  cotSlopeYErr[n_cot]=sqrt(cot_slope_y_err[i]);
	  cotRun[n_cot]      =run[i];
	  if ( (cot_y[i] == 0.16000) || (cot_y[i] == 0.17) )
	    std::cout<<"Dummy value for run "<<run<<std::endl;
	  n_cot++;
	}
    }

  //---------------------------------------------------------
  // calculate differences and pulls
  //---------------------------------------------------------
  double *diffX             = new double[npt];
  double *diffY             = new double[npt];
  double *diffXErr          = new double[npt];
  double *diffYErr          = new double[npt];
  double *diffSlopeX        = new double[npt];
  double *diffSlopeY        = new double[npt];
  double *diffSlopeXErr     = new double[npt];
  double *diffSlopeYErr     = new double[npt];
  double *pullX             = new double[npt];
  double *pullY             = new double[npt];
  double *pullXErr          = new double[npt];
  double *pullYErr          = new double[npt];
  double *pullSlopeX        = new double[npt];
  double *pullSlopeY        = new double[npt];
  double *diffRun           = new double[npt];
  
  TH1F *h_pullX      = new TH1F("h_pullX", "pulls x", 20, -5, 10);
  TH1F *h_pullY      = new TH1F("h_pullY", "pulls y", 20, -10, 10);
  TH1F *h_pullSlopeX = new TH1F("h_pullSlopeX", "pulls slopex", 20, -5, 5);
  TH1F *h_pullSlopeY = new TH1F("h_pullSlopeY", "pulls slopey", 20, -5, 5);

  int n_cot_svx=0;
  for (Int_t i = 0; i<n_svx; i++)
    for (int j=0; j<n_cot; j++)
      {
	if (svxRun[i]== cotRun[j]){
	  
	  diffX[n_cot_svx] = svxBeamX[i]-cotBeamX[j];
	  diffY[n_cot_svx] = svxBeamY[i]-cotBeamY[j];
	  diffSlopeX[n_cot_svx] = svxSlopeX[i]-cotSlopeX[j];
	  diffSlopeY[n_cot_svx] = svxSlopeY[i]-cotSlopeY[j];

	  if (diffY[n_cot_svx]>0.002) 
	    std::cout << svxRun[i] <<"  DeltaY="<< diffY[n_cot_svx] << endl;
	  
	  diffXErr[n_cot_svx] = sqrt(svxBeamXErr[i]*svxBeamXErr[i]+cotBeamXErr[j]*cotBeamXErr[j]);
	  diffYErr[n_cot_svx] = sqrt(svxBeamYErr[i]*svxBeamYErr[i]+cotBeamYErr[j]*cotBeamYErr[j]);
	  diffSlopeXErr[n_cot_svx] = sqrt(svxSlopeXErr[i]*svxSlopeXErr[i]+cotSlopeXErr[j]*cotSlopeXErr[j]);
	  diffSlopeYErr[n_cot_svx] = sqrt(svxSlopeYErr[i]*svxSlopeYErr[i]+cotSlopeYErr[j]*cotSlopeYErr[j]);
	  
	  diffRun[n_cot_svx] = svxRun[i];
	  pullX[n_cot_svx] = diffX[n_cot_svx]/diffXErr[n_cot_svx];
	  pullY[n_cot_svx] = diffY[n_cot_svx]/diffYErr[n_cot_svx];
	  pullSlopeX[n_cot_svx] = diffSlopeX[n_cot_svx]/diffSlopeXErr[n_cot_svx];
	  pullSlopeY[n_cot_svx] = diffSlopeY[n_cot_svx]/diffSlopeYErr[n_cot_svx];

	  if(pullY[n_cot_svx]>6) 
	    std::cout<<svxRun[i]<<"  pull="<<pullY[n_cot_svx]<<"  diffY="<<diffY[n_cot_svx]<<"   errdiffY="<<diffYErr[n_cot_svx]<<std::endl; 

	  h_pullX->Fill(pullX[n_cot_svx]);
	  h_pullY->Fill(pullY[n_cot_svx]);
	  h_pullSlopeX->Fill(pullSlopeX[n_cot_svx]);
	  h_pullSlopeY->Fill(pullSlopeY[n_cot_svx]);

	  n_cot_svx++;
	}
      }

  
  cout << "TOTAL RUNS: " << npt << endl;
  cout << "COT FITTED RUNS: " << n_cot << endl;
  cout << "SVX FITTED RUNS: " << n_svx << endl;
  cout << "COT AND SVX FITTED RUNS: " << n_cot_svx << endl;


  //------------------------------------
  // get the plots for fit code
  //------------------------------------
  plot_beamlines("svx_fit_code");
  TGraph *g_svx_fit_code = (TGraph*) c1->FindObject("Graph");
  c1->Close();
  
  plot_beamlines("cot_fit_code");
  TGraph *g_cot_fit_code = (TGraph*) c1->FindObject("Graph");
  c1->Close();
  
  /*
  plot_beamlines("svx_x");
  plot_beamlines("svx_y");
  plot_beamlines("svx_slope_x");
  plot_beamlines("svx_slope_y");
  plot_beamlines("cot_x");
  plot_beamlines("cot_y");
  plot_beamlines("cot_slope_x");
  plot_beamlines("cot_slope_y");
  plot_beamlines("dx");
  plot_beamlines("dy");
  plot_beamlines("dx_slope");
  plot_beamlines("dy_slope");
  */

  //------------------------------------------------------------
  // make and save the plots in .eps files (batch mode)
  //-----------------------------------------------------------
  int id=0;
  TCanvas *cc[100];

  // FIT CODES
  cc[id] = new TCanvas(Form("cc_%d",id),"COT fit code");
  cc[id]->cd();
  g_cot_fit_code->SetTitle("COT Beamline fit code");
  //g_cot_fit_code->SetMarkerStyle(20);
  //g_cot_fit_code->SetMarkerSize(1);
  //g_cot_fit_code->GetXaxis()->SetLabelSize(0.02);
  g_cot_fit_code->Draw("AP");
  cc[id]->SaveAs(Form("fitcode_cot_%d_%d.eps",MinRun,MaxRun));

  id++;
  cc[id] = new TCanvas(Form("cc_%d",id),"SVX fit code");
  cc[id]->cd();
  g_svx_fit_code->SetTitle("SVX Beamline fit code");
  g_svx_fit_code->Draw("AP");
  cc[id]->SaveAs(Form("fitcode_svx_%d_%d.eps",MinRun,MaxRun));

  // SVX
  id++;
  cc[id] = new TCanvas(Form("cc_%d",id),"SVX",800,800);
  cc[id]->Divide(2,2);
  cc[id]->cd(1);
  TGraphErrors* svx_xplot = new TGraphErrors(n_svx, svxRun, svxBeamX, 0, svxBeamXErr);
  svx_xplot->SetMarkerColor(2);
  svx_xplot->SetMarkerStyle(7);
  svx_xplot->SetTitle("beamX: SVX");
  svx_xplot->Draw("AP");
  cc[id]->cd(2);
  TGraphErrors* svx_yplot = new TGraphErrors(n_svx, svxRun, svxBeamY, 0, svxBeamYErr);
  svx_yplot->SetMarkerColor(2);
  svx_yplot->SetMarkerStyle(7);
  svx_yplot->SetTitle("beamY: SVX");
  svx_yplot->Draw("AP");
  cc[id]->cd(3);
  TGraphErrors* svx_slopexplot = new TGraphErrors(n_svx, svxRun, svxSlopeX, 0, svxSlopeXErr);
  svx_slopexplot->SetMarkerColor(2);
  svx_slopexplot->SetMarkerStyle(7);
  svx_slopexplot->SetTitle("beam slopeX: SVX");
  svx_slopexplot->Draw("AP");
  cc[id]->cd(4);
  TGraphErrors* svx_slopeyplot = new TGraphErrors(n_svx, svxRun, svxSlopeY, 0, svxSlopeYErr);
  svx_slopeyplot->SetMarkerColor(2);
  svx_slopeyplot->SetMarkerStyle(7);
  svx_slopeyplot->SetTitle("beam slopeY: SVX");
  svx_slopeyplot->Draw("AP");
  cc[id]->SaveAs(Form("svx_%d_%d.eps",MinRun,MaxRun));

  // COT
  id++;
  cc[id] = new TCanvas(Form("cc_%d",id),"COT",800,800);
  cc[id]->Divide(2,2);
  cc[id]->cd(1);
  TGraphErrors* cot_xplot = new TGraphErrors(n_cot, cotRun, cotBeamX, 0, cotBeamXErr);
  cot_xplot->SetMarkerColor(2);
  cot_xplot->SetMarkerStyle(7);
  cot_xplot->SetTitle("beamX: COT");
  cot_xplot->Draw("AP");
  cc[id]->cd(2);
  TGraphErrors* cot_yplot = new TGraphErrors(n_cot, cotRun, cotBeamY, 0, cotBeamYErr);
  cot_yplot->SetMarkerColor(2);
  cot_yplot->SetMarkerStyle(7);
  cot_yplot->SetTitle("beamY: COT");
  cot_yplot->Draw("AP");
  cc[id]->cd(3);
  TGraphErrors* cot_slopexplot = new TGraphErrors(n_cot, cotRun, cotSlopeX, 0, cotSlopeXErr);
  cot_slopexplot->SetMarkerColor(2);
  cot_slopexplot->SetMarkerStyle(7);
  cot_slopexplot->SetTitle("beam slopeX: COT");
  cot_slopexplot->Draw("AP");
  cc[id]->cd(4);
  TGraphErrors* cot_slopeyplot = new TGraphErrors(n_cot, cotRun, cotSlopeY, 0, cotSlopeYErr);
  cot_slopeyplot->SetMarkerColor(2);
  cot_slopeyplot->SetMarkerStyle(7);
  cot_slopeyplot->SetTitle("beam slopeY: COT");
  cot_slopeyplot->Draw("AP");
  cc[id]->SaveAs(Form("cot_%d_%d.eps",MinRun,MaxRun));

  // difference between SVX and COT
  id++;
  cc[id] = new TCanvas(Form("cc_%d",id),"Diff",800,800);
  cc[id]->Divide(2,2);
  cc[id]->cd(1);
  diff_xplot = new TGraphErrors(n_cot_svx, diffRun, diffX, 0, diffXErr);
  diff_xplot->SetMarkerColor(103);
  diff_xplot->SetMarkerStyle(7);
  diff_xplot->SetTitle("delta beamX: SVX-COT");
  diff_xplot->Draw("AP");
  cc[id]->cd(2);
  diff_yplot = new TGraphErrors(n_cot_svx, diffRun, diffY, 0, diffYErr);
  diff_yplot->SetMarkerColor(103);
  diff_yplot->SetMarkerStyle(7);
  diff_yplot->SetTitle("delta beamY: SVX-COT");
  diff_yplot->Draw("AP");
  cc[id]->cd(3);
  diff_slopexplot = new TGraphErrors(n_cot_svx, diffRun, diffSlopeX, 0, diffSlopeXErr);
  diff_slopexplot->SetMarkerColor(103);
  diff_slopexplot->SetMarkerStyle(7);
  diff_slopexplot->SetTitle("delta beam slopeX: SVX-COT");
  diff_slopexplot->Draw("AP");
  cc[id]->cd(4);
  diff_slopeyplot = new TGraphErrors(n_cot_svx, diffRun, diffSlopeY, 0, diffSlopeYErr);
  diff_slopeyplot->SetMarkerColor(103);
  diff_slopeyplot->SetMarkerStyle(7);
  diff_slopeyplot->SetTitle("delta beam slopeY: SVX-COT");
  diff_slopeyplot->Draw("AP");
  cc[id]->SaveAs(Form("diff_%d_%d.eps",MinRun,MaxRun));

  // PULLS
  id++;
  cc[id] = new TCanvas(Form("cc_%d",id),"Pull", 800, 800);
  cc[id]->Divide(2,2);
  cc[id]->cd(1);
  pull_xplot = new TGraphErrors(n_cot_svx, diffRun, pullX, 0, 0);
  pull_xplot->SetMarkerStyle(7);
  pull_xplot->SetTitle("pull beamX: SVX-COT");
  pull_xplot->Draw("AP");
  cc[id]->cd(2);
  pull_yplot = new TGraphErrors(n_cot_svx, diffRun, pullY, 0, 0);
  pull_yplot->SetMarkerStyle(7);
  pull_yplot->SetTitle("pull beamY: SVX-COT");
  pull_yplot->Draw("AP");
  cc[id]->cd(3);
  pull_slopexplot = new TGraphErrors(n_cot_svx, diffRun, pullSlopeX, 0, 0);
  pull_slopexplot->SetMarkerStyle(7);
  pull_slopexplot->SetTitle("pull slopeX: SVX-COT");
  pull_slopexplot->Draw("AP");
  cc[id]->cd(4);
  pull_slopeyplot = new TGraphErrors(n_cot_svx, diffRun, pullSlopeY, 0, 0);
  pull_slopeyplot->SetMarkerStyle(7);
  pull_slopeyplot->SetTitle("pull beam: SVX-COT");
  pull_slopeyplot->Draw("AP");
  cc[id]->SaveAs(Form("pulls_%d_%d.eps",MinRun,MaxRun));

  // PULL DISTRIBUTIONS
  id++;
  cc[id] = new TCanvas(Form("cc_%d",id), "pulls", 800, 800); 
  cc[id]->Divide(2,2);
  cc[id]->cd(1);
  h_pullX->Draw();
  h_pullX->SetFillColor(38);
  h_pullX->Fit("gaus");
  cc[id]->cd(2);
  h_pullY->Draw();
  h_pullY->SetFillColor(38);
  h_pullY->Fit("gaus");
  cc[id]->cd(3);
  h_pullSlopeX->Draw();
  h_pullSlopeX->SetFillColor(38);
  h_pullSlopeX->Fit("gaus");
  cc[id]->cd(4);
  h_pullSlopeY->Draw();
  h_pullSlopeY->SetFillColor(38);
  h_pullSlopeY->Fit("gaus");
  cc[id]->SaveAs(Form("pulls_distr_%d_%d.eps",MinRun,MaxRun));

  return 0;
}