/***************************************************************************
 *
 * 2015/02/12 Kei Moriya
 *
 * Plotter for output of analyze_KsHunt.
 *
 ***************************************************************************/

#include <iostream>
#include <iomanip>
#include <fstream>
#include <cmath>
#include <vector>
#include <sstream>
#include <map>
#include "unistd.h" // to use optarg
#include <assert.h>

// ROOT header files
#include "TROOT.h"
#include "TFile.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TProfile.h"
#include "TTree.h"
#include "TMath.h"
#include "TLatex.h"
#include "TLine.h"
#include "TLegend.h"
#include "TPaveText.h"
#include "TCanvas.h"
#include "TPad.h"
#include "TLorentzVector.h"
#include "TProfile.h"
#include "TGraphErrors.h"

#include "rootStyles.h"
#include "/home/kmoriya/include/physicsConstants.h"

using namespace std;

void printUsage();

int main(int argc, char **argv){
  string infilename = "";
  extern char* optarg;
  // Check command line arguments
  int c;
  while((c = getopt(argc,argv,"hi:")) != -1){
    switch(c){
    case 'h':
      printUsage();
      exit(-1);
      break;
    case 'i':
      // Specify infile
      infilename = optarg;
      cout << "infile name : " << infilename << endl;
      break;
      // case 'c':
      //   cout << "create plots for each channel" << endl;
      //   doChannels = true;
      // break;
    default:
      break;
    }
  }
  //___________________________________________________________________________________________

  if(infilename==""){
    printUsage();
    return 0;
  }

  gStyle->SetOptStat(0); // 111111
  gStyle->SetHistFillStyle(0);
  gStyle->SetHistLineColor(kBlue-7);
  gStyle->SetTitleSize(0.060,"XY");
  gStyle->SetLabelSize(0.040,"XY");
  gStyle->SetTitleOffset(0.800,"X");
  gStyle->SetTitleOffset(0.900,"Y");
  gStyle->SetTitleFont(132,"XY");
  gStyle->SetLabelFont(132,"XY");
  gROOT->ForceStyle();

  char outdir[200];
  sprintf(outdir,"/home/kmoriya/halld_my/KsHunt/analysis/figures/%s",infilename.c_str());
  char command[200];
  sprintf(command,"mkdir -p %s",outdir);
  system(command);

  TFile *infile = new TFile(infilename.c_str());
  if(!infile){
    cout << "File " << infilename << " does not exist" << endl;
    abort();
  }

  char hname[200];
  char hnewname[200];

  // kfit info
  TH1F *hkfitCL = (TH1F*)infile->Get("hkfitCL")->Clone("_hkfitCL");

  // M(pi+,pi-)
  TH1F *hM2pi = (TH1F*)infile->Get("hM2pi")->Clone("_hM2pi");
  TH2F *hvrho_M2pi = (TH2F*)infile->Get("hvrho_M2pi")->Clone("_hvrho_M2pi");
  TH2F *hvz_M2pi = (TH2F*)infile->Get("hvz_M2pi")->Clone("_hvz_M2pi");
  TH1F *hM2pi_allVertexSelect = (TH1F*)infile->Get("hM2pi_allVertexSelect")->Clone("_hM2pi_allVertexSelect");

  // M(pi+,pi-,unused pi+-)
  TH1F *hM3pi = (TH1F*)infile->Get("hM3pi")->Clone("_hM3pi");
  TH1F *hM3pi_allVertexKsSelect = (TH1F*)infile->Get("hM3pi_allVertexKsSelect")->Clone("_hM3pi_allVertexKsSelect");

  // M(unused p,unused pi-)
  TH1F *hMppi = (TH1F*)infile->Get("hMppi")->Clone("_hMppi");
  TH1F *hMppi_allVertexKsSelect = (TH1F*)infile->Get("hMppi_allVertexKsSelect")->Clone("_hMppi_allVertexKsSelect");

  // MM from Ks vs M(pi+,pi-)
  TH2F *hMM_M2pi = (TH2F*)infile->Get("hMM_M2pi")->Clone("_hMM_M2pi");
  TH2F *hMM_M2pi_allVertexKsSelect = (TH2F*)infile->Get("hMM_M2pi_allVertexKsSelect")->Clone("_hMM_M2pi_allVertexKsSelect");

  // Momentum
  TH2F *hPpip_Ppim = (TH2F*)infile->Get("hPpip_Ppim")->Clone("_hPpip_Ppim");
  TH2F *hPpip_Ppim_allVertexKsSelect = (TH2F*)infile->Get("hPpip_Ppim_allVertexKsSelect")->Clone("_hPpip_Ppim_allVertexKsSelect");;
  TH1F *hPKs = (TH1F*)infile->Get("hPKs")->Clone("_hPKs");
  TH1F *hPKs_allVertexKsSelect = (TH1F*)infile->Get("hPKs_allVertexKsSelect")->Clone("_hPKs_allVertexKsSelect");

  // dE/dx
  TH2F *hdEdx_pip_CDC = (TH2F*)infile->Get("hdEdx_pip_CDC")->Clone("_hdEdx_pip_CDC");
  TH2F *hdEdx_pip_FDC = (TH2F*)infile->Get("hdEdx_pip_FDC")->Clone("_hdEdx_pip_FDC");
  TH2F *hdEdx_pip_ST = (TH2F*)infile->Get("hdEdx_pip_ST")->Clone("_hdEdx_pip_ST");

  TH2F *hdEdx_pim_CDC = (TH2F*)infile->Get("hdEdx_pim_CDC")->Clone("_hdEdx_pim_CDC");
  TH2F *hdEdx_pim_FDC = (TH2F*)infile->Get("hdEdx_pim_FDC")->Clone("_hdEdx_pim_FDC");
  TH2F *hdEdx_pim_ST = (TH2F*)infile->Get("hdEdx_pim_ST")->Clone("_hdEdx_pim_ST");

  TH2F *hdEdx_pip_CDC_selectKs = (TH2F*)infile->Get("hdEdx_pip_CDC_selectKs")->Clone("_hdEdx_pip_CDC_selectKs");
  TH2F *hdEdx_pip_FDC_selectKs = (TH2F*)infile->Get("hdEdx_pip_FDC_selectKs")->Clone("_hdEdx_pip_FDC_selectKs");
  TH2F *hdEdx_pip_ST_selectKs = (TH2F*)infile->Get("hdEdx_pip_ST_selectKs")->Clone("_hdEdx_pip_ST_selectKs");

  TH2F *hdEdx_pim_CDC_selectKs = (TH2F*)infile->Get("hdEdx_pim_CDC_selectKs")->Clone("_hdEdx_pim_CDC_selectKs");
  TH2F *hdEdx_pim_FDC_selectKs = (TH2F*)infile->Get("hdEdx_pim_FDC_selectKs")->Clone("_hdEdx_pim_FDC_selectKs");
  TH2F *hdEdx_pim_ST_selectKs = (TH2F*)infile->Get("hdEdx_pim_ST_selectKs")->Clone("_hdEdx_pim_ST_selectKs");

  // unused K+-
  TH2F *hdEdx_kp_CDC = (TH2F*)infile->Get("hdEdx_kp_CDC")->Clone("_hdEdx_kp_CDC");
  TH2F *hdEdx_kp_FDC = (TH2F*)infile->Get("hdEdx_kp_FDC")->Clone("_hdEdx_kp_FDC");
  TH2F *hdEdx_kp_ST = (TH2F*)infile->Get("hdEdx_kp_ST")->Clone("_hdEdx_kp_ST");

  TH2F *hdEdx_km_CDC = (TH2F*)infile->Get("hdEdx_km_CDC")->Clone("_hdEdx_km_CDC");
  TH2F *hdEdx_km_FDC = (TH2F*)infile->Get("hdEdx_km_FDC")->Clone("_hdEdx_km_FDC");
  TH2F *hdEdx_km_ST = (TH2F*)infile->Get("hdEdx_km_ST")->Clone("_hdEdx_km_ST");

  TH2F *hdEdx_kp_CDC_selectKs = (TH2F*)infile->Get("hdEdx_kp_CDC_selectKs")->Clone("_hdEdx_kp_CDC_selectKs");
  TH2F *hdEdx_kp_FDC_selectKs = (TH2F*)infile->Get("hdEdx_kp_FDC_selectKs")->Clone("_hdEdx_kp_FDC_selectKs");
  TH2F *hdEdx_kp_ST_selectKs = (TH2F*)infile->Get("hdEdx_kp_ST_selectKs")->Clone("_hdEdx_kp_ST_selectKs");

  TH2F *hdEdx_km_CDC_selectKs = (TH2F*)infile->Get("hdEdx_km_CDC_selectKs")->Clone("_hdEdx_km_CDC_selectKs");
  TH2F *hdEdx_km_FDC_selectKs = (TH2F*)infile->Get("hdEdx_km_FDC_selectKs")->Clone("_hdEdx_km_FDC_selectKs");
  TH2F *hdEdx_km_ST_selectKs = (TH2F*)infile->Get("hdEdx_km_ST_selectKs")->Clone("_hdEdx_km_ST_selectKs");

  // vertices
  TH2F *hvy_vx_pip_measured = (TH2F*)infile->Get("hvy_vx_pip_measured")->Clone("_hvy_vx_pip_measured");
  TH2F *hvy_vx_pip = (TH2F*)infile->Get("hvy_vx_pip")->Clone("_hvy_vx_pip");

  TH1F *hvz_pip_measured = (TH1F*)infile->Get("hvz_pip_measured")->Clone("_hvz_pip_measured");
  TH1F *hvz_pip = (TH1F*)infile->Get("hvz_pip")->Clone("_hvz_pip");

  TH2F *hvy_vx_pim_measured = (TH2F*)infile->Get("hvy_vx_pim_measured")->Clone("_hvy_vx_pim_measured");
  TH2F *hvy_vx_pim = (TH2F*)infile->Get("hvy_vx_pim")->Clone("_hvy_vx_pim");

  TH1F *hvz_pim_measured = (TH1F*)infile->Get("hvz_pim_measured")->Clone("_hvz_pim_measured");
  TH1F *hvz_pim = (TH1F*)infile->Get("hvz_pim")->Clone("_hvz_pim");

  TH2F *hvz_pip_pim_measured = (TH2F*)infile->Get("hvz_pip_pim_measured")->Clone("_hvz_pip_pim_measured");
  TH2F *hvz_pip_pim = (TH2F*)infile->Get("hvz_pip_pim")->Clone("_hvz_pip_pim");

  TCanvas *canvas = new TCanvas("canvas","canvas",1500,800);
  canvas->SetLeftMargin(0.12);
  canvas->SetRightMargin(0.10);
  canvas->SetTopMargin(0.03);
  canvas->SetBottomMargin(0.14);
  char filename[400];

  TLatex *latex = new TLatex();
  latex->SetTextColor(kBlue);
  latex->SetTextSize(0.055);
  latex->SetTextFont(132);
  latex->SetTextAlign(12);
  latex->SetNDC(1);
  char text[400];

  TLegend *legend = new TLegend(0.20,0.75,0.90,0.95);
  legend->SetNColumns(2);
  legend->SetFillStyle(0);
  legend->SetBorderSize(0);
  legend->SetTextSize(0.040);

  Int_t plot = 1;

  TLine *line_Ks = new TLine(m_KZero,0,m_KZero,0);
  line_Ks->SetLineColor(kRed);
  line_Ks->SetLineStyle(2);
  line_Ks->SetLineWidth(1);

  TLine *line_Kstar = new TLine(m_KstarPlus,0,m_KstarPlus,0);
  line_Kstar->SetLineColor(kRed);
  line_Kstar->SetLineStyle(2);
  line_Kstar->SetLineWidth(1);

  TLine *line_Lambda = new TLine(m_Lambda,0,m_Lambda,0);
  line_Lambda->SetLineColor(kRed);
  line_Lambda->SetLineStyle(2);
  line_Lambda->SetLineWidth(1);

  TLine *line_proton = new TLine(m_p,0,m_p,0);
  line_proton->SetLineColor(kRed);
  line_proton->SetLineStyle(2);
  line_proton->SetLineWidth(1);

  TLine *line_SigmaPlus = new TLine(m_SigmaPlus,0,m_SigmaPlus,0);
  line_SigmaPlus->SetLineColor(kRed);
  line_SigmaPlus->SetLineStyle(2);
  line_SigmaPlus->SetLineWidth(1);

  // kfit info
  canvas->SetLogy(1);
  hkfitCL->Draw();
  sprintf(filename,"%s/%2.2d___kfitCL.pdf",outdir,plot);
  canvas->SaveAs(filename);
  plot++;
  canvas->SetLogy(0);
  //_______________________________________________________________________________________________________

  // M(pi+,pi-)
  hM2pi->SetAxisRange(0.20,1.0);
  hM2pi->Draw();
  line_Ks->SetY2(hM2pi->GetMaximum() * 1.05);
  line_Ks->Draw("same");
  sprintf(filename,"%s/%2.2d___M2pi.pdf",outdir,plot);
  canvas->SaveAs(filename);
  plot++;
  //_______________________________________________________________________________________________________

  // Plot with cuts on vrho, vz, Ks mass
  hM2pi_allVertexSelect->SetAxisRange(0.20,1.0);
  hM2pi_allVertexSelect->Draw();

  line_Ks->SetY2(hM2pi_allVertexSelect->GetMaximum() * 1.05);
  line_Ks->Draw("same");

  // Fit
  TF1 *fM2pi_allVertexSelect = new TF1("fM2pi_allVertexSelect","[0]*exp(-pow(x-[1],2.)/2./pow([2],2.)) + [3] + [4]*x + [5]* pow(x,2.)",0.4,0.6);
  fM2pi_allVertexSelect->SetParameter(0,hM2pi_allVertexSelect->GetMaximum());
  fM2pi_allVertexSelect->SetParameter(1,m_KZero);
  fM2pi_allVertexSelect->SetParameter(2,0.050);
  hM2pi_allVertexSelect->Fit(fM2pi_allVertexSelect,"EQRN");
  fM2pi_allVertexSelect->SetLineColor(kRed);
  fM2pi_allVertexSelect->SetLineStyle(2);
  fM2pi_allVertexSelect->SetLineWidth(1);
  fM2pi_allVertexSelect->Draw("same");

  // Plot signal Gaussian
  TF1 *fSignalM2pi_allVertexSelect = new TF1("fSignalM2pi_allVertexSelect","[0]*exp(-pow(x-[1],2.)/2./pow([2],2.))",0.4,0.6);
  fSignalM2pi_allVertexSelect->SetParameter(0,fM2pi_allVertexSelect->GetParameter(0));
  fSignalM2pi_allVertexSelect->SetParameter(1,fM2pi_allVertexSelect->GetParameter(1));
  fSignalM2pi_allVertexSelect->SetParameter(2,fM2pi_allVertexSelect->GetParameter(2));
  fSignalM2pi_allVertexSelect->SetLineColor(kCyan+2);
  fSignalM2pi_allVertexSelect->SetLineStyle(2);
  fSignalM2pi_allVertexSelect->SetLineWidth(1);
  fSignalM2pi_allVertexSelect->Draw("same");

  // Plot bg
  TF1 *fBgM2pi_allVertexSelect = new TF1("fBgM2pi_allVertexSelect","[0] + [1] * x + [2] * pow(x,2.)",0.4,0.6);
  fBgM2pi_allVertexSelect->SetParameter(0,fM2pi_allVertexSelect->GetParameter(3));
  fBgM2pi_allVertexSelect->SetParameter(1,fM2pi_allVertexSelect->GetParameter(4));
  fBgM2pi_allVertexSelect->SetParameter(2,fM2pi_allVertexSelect->GetParameter(5));
  fBgM2pi_allVertexSelect->SetLineColor(kGreen+2);
  fBgM2pi_allVertexSelect->SetLineStyle(2);
  fBgM2pi_allVertexSelect->SetLineWidth(1);
  fBgM2pi_allVertexSelect->Draw("same");

  sprintf(text,"width: %4.1f #pm %4.1f MeV",fM2pi_allVertexSelect->GetParameter(2)* 1000.,fM2pi_allVertexSelect->GetParError(2)* 1000.);
  latex->DrawLatex(0.55,0.80,text);

  sprintf(text,"m_{fit} = %4.2f #pm %4.2f MeV/c^{2}",fM2pi_allVertexSelect->GetParameter(1)* 1000, fM2pi_allVertexSelect->GetParError(1)* 1000);
  latex->DrawLatex(0.55,0.75,text);

  sprintf(text,"m_{PDG} = %4.2f MeV/c^{2}",m_KZero * 1000.);
  latex->DrawLatex(0.55,0.70,text);

  sprintf(text,"counts: %5.1f #pm %5.1f",
	  fM2pi_allVertexSelect->GetParameter(0) * sqrt(2. * TMath::Pi()) * fM2pi_allVertexSelect->GetParameter(2) / hM2pi_allVertexSelect->GetBinWidth(1),
	  fM2pi_allVertexSelect->GetParError(0)  * sqrt(2. * TMath::Pi()) * fM2pi_allVertexSelect->GetParameter(2) / hM2pi_allVertexSelect->GetBinWidth(1));
  latex->DrawLatex(0.55,0.65,text);

  sprintf(filename,"%s/%2.2d___M2pi_allVertexSelect.pdf",outdir,plot);
  canvas->SaveAs(filename);
  plot++;
  //_______________________________________________________________________________________________________

  // vrho vs M(pi+,pi-)
  canvas->SetLogz(1);
  hvrho_M2pi->SetAxisRange(0.20,1.0);
  hvrho_M2pi->Draw("colz");
  line_Ks->SetY1(hvrho_M2pi->GetYaxis()->GetXmin());
  line_Ks->SetY2(hvrho_M2pi->GetYaxis()->GetXmax());
  line_Ks->Draw("same");
  sprintf(filename,"%s/%2.2d___vrho_M2pi.pdf",outdir,plot);
  canvas->SaveAs(filename);
  plot++;
  canvas->SetLogz(0);
  //_______________________________________________________________________________________________________

  // vrho vs M(pi+,pi-) (closeup)
  canvas->SetLogz(1);
  hvrho_M2pi->SetAxisRange(0,5,"Y");
  hvrho_M2pi->Draw("colz");
  line_Ks->SetY1(hvrho_M2pi->GetYaxis()->GetXmin());
  line_Ks->SetY2(hvrho_M2pi->GetYaxis()->GetXmax());
  line_Ks->Draw("same");
  sprintf(filename,"%s/%2.2d___vrho_M2pi_close.pdf",outdir,plot);
  canvas->SaveAs(filename);
  plot++;
  canvas->SetLogz(0);
  //_______________________________________________________________________________________________________

  // Project onto M(pi+,pi-) with cut on vrho
  Int_t firstbin, lastbin;
  Double_t min = 1.0; // cm
  Double_t y0 = hvrho_M2pi->GetYaxis()->GetBinLowEdge(1);
  Double_t biny = hvrho_M2pi->GetYaxis()->GetBinWidth(1);
  firstbin = (Int_t)((min - y0)/biny);
  lastbin  = hvrho_M2pi->GetNbinsY();
  cout << "y0 = " << y0 << " biny = " << biny << endl;
  cout << "firstbin = " << firstbin << endl;
  TH1F *hM2pi_vrhocut = (TH1F*)hvrho_M2pi->ProjectionX("hM2pi_vrhocut",firstbin, lastbin);
  hM2pi_vrhocut->Draw();
  line_Ks->SetY1(0);
  line_Ks->SetY2(hM2pi_vrhocut->GetMaximum() * 1.05);
  line_Ks->Draw("same");

  sprintf(text,"v_{#rho} > %5.3f cm",min);
  latex->DrawLatex(0.60,0.85,text);

  sprintf(filename,"%s/%2.2d___M2pi_vrhocut.pdf",outdir,plot);
  canvas->SaveAs(filename);
  plot++;
  //_______________________________________________________________________________________________________

  // Get Ks yield as function of vrho cut
  TGraphErrors *gKsYieldAgainstVrho = new TGraphErrors();
  TGraphErrors *gKsSBAgainstVrho = new TGraphErrors();
  TGraphErrors *gKsFOMAgainstVrho = new TGraphErrors();
  TH1F *hM2pi_vrhocut_vary;
  TF1 *fM2pi_vrhocut_vary;
  Double_t yield,error;
  Int_t first, last;
  Double_t area;
  Double_t sob, fom;
  for(Int_t ycutbin=0;ycutbin<hvrho_M2pi->GetNbinsY();ycutbin++){
    // cout << "fitting to min bin = " << ycutbin << endl;
    // Get minimum of bin
    Double_t ymin = hvrho_M2pi->GetYaxis()->GetBinLowEdge(ycutbin);
    hM2pi_vrhocut_vary = (TH1F*)hvrho_M2pi->ProjectionX("hM2pi_vrhocut_vary",ycutbin,lastbin);
    // Do fit
    fM2pi_vrhocut_vary = new TF1("fM2pi_vrhocut_vary","[0]*exp(-pow(x-[1],2.)/2./pow([2],2.)) + [3] + [4]*x + [5]* pow(x,2.)",0.4,0.6);
    fM2pi_vrhocut_vary->SetParameter(0,hM2pi_vrhocut_vary->GetMaximum() * 0.5);
    fM2pi_vrhocut_vary->SetParameter(1,m_KZero);
    fM2pi_vrhocut_vary->SetParameter(2,0.013);
    hM2pi_vrhocut_vary->Fit(fM2pi_vrhocut_vary,"EQRN");

    // Get yields, put in graph
    yield = fM2pi_vrhocut_vary->GetParameter(0) * sqrt(2. * TMath::Pi()) * fabs(fM2pi_vrhocut_vary->GetParameter(2)) / hM2pi_vrhocut_vary->GetBinWidth(1);
    error = fM2pi_vrhocut_vary->GetParError(0)  * sqrt(2. * TMath::Pi()) * fabs(fM2pi_vrhocut_vary->GetParameter(2)) / hM2pi_vrhocut_vary->GetBinWidth(1);
    gKsYieldAgainstVrho->SetPoint(gKsYieldAgainstVrho->GetN(),ymin,yield);
    gKsYieldAgainstVrho->SetPointError(gKsYieldAgainstVrho->GetN()-1,0,error);
    // cout << "fit to min bin = " << ycutbin << " ymin = " << ymin << " yield = " << yield << " +/- " << error << endl;
    
    // Set S/B
    // Defined as total yield of Gaussian divided by area under +/- 0.013 * 3 MeV of Ks
    sob   = yield / (fM2pi_vrhocut_vary->Integral(m_KZero - 0.013 * 3.,m_KZero + 0.013 * 3.) / hM2pi_vrhocut_vary->GetBinWidth(1));
    error = error / (fM2pi_vrhocut_vary->Integral(m_KZero - 0.013 * 3.,m_KZero + 0.013 * 3.) / hM2pi_vrhocut_vary->GetBinWidth(1));
    gKsSBAgainstVrho->SetPoint(gKsSBAgainstVrho->GetN(),ymin,sob);
    gKsSBAgainstVrho->SetPointError(gKsSBAgainstVrho->GetN()-1,0,error);
    // cout << "yield = " << yield
    //      << " total = " << fM2pi_vrhocut_vary->Integral(m_KZero - 0.013 * 3.,m_KZero + 0.013 * 3.) / hM2pi_vrhocut_vary->GetBinWidth(1)
    //      << " bin width = " << hM2pi_vrhocut_vary->GetBinWidth(1)
    //      << " S/B = " << sob << endl;

    // Set FOM
    fom   = yield / sqrt((fM2pi_vrhocut_vary->Integral(m_KZero - 0.013 * 3.,m_KZero + 0.013 * 3.) / hM2pi_vrhocut_vary->GetBinWidth(1)));
    gKsFOMAgainstVrho->SetPoint(gKsFOMAgainstVrho->GetN(),ymin,fom);

    hM2pi_vrhocut_vary = 0;
    fM2pi_vrhocut_vary = 0;
  }

  // Graph of yield of Ks against vrho cut
  gKsYieldAgainstVrho->SetLineColor(kBlue);
  gKsYieldAgainstVrho->SetLineWidth(1);
  gKsYieldAgainstVrho->SetLineStyle(1);
  gKsYieldAgainstVrho->SetMarkerColor(kBlue);
  gKsYieldAgainstVrho->SetMarkerSize(1.5);
  gKsYieldAgainstVrho->SetMarkerStyle(20);
  gKsYieldAgainstVrho->GetXaxis()->SetTitle("v_{#rho} min (cm)");
  gKsYieldAgainstVrho->GetYaxis()->SetTitle("K_{S} yield");

  gKsYieldAgainstVrho->GetXaxis()->SetLimits(0,10);
  gKsYieldAgainstVrho->SetMinimum(0);
  gKsYieldAgainstVrho->SetMaximum(2000.);
  gKsYieldAgainstVrho->Draw("ALP");
  sprintf(filename,"%s/%2.2d___gKsYieldAgainstVrho.pdf",outdir,plot);
  canvas->SaveAs(filename);
  plot++;
  canvas->SetLogz(0);
  //_______________________________________________________________________________________________________

  // Graph of S/B of Ks against vrho cut
  gKsSBAgainstVrho->SetLineColor(kBlue);
  gKsSBAgainstVrho->SetLineWidth(1);
  gKsSBAgainstVrho->SetLineStyle(1);
  gKsSBAgainstVrho->SetMarkerColor(kBlue);
  gKsSBAgainstVrho->SetMarkerSize(1.5);
  gKsSBAgainstVrho->SetMarkerStyle(20);
  gKsSBAgainstVrho->GetXaxis()->SetTitle("v_{#rho} min (cm)");
  gKsSBAgainstVrho->GetYaxis()->SetTitle("K_{S} S/B");

  gKsSBAgainstVrho->GetXaxis()->SetLimits(0,10);
  gKsSBAgainstVrho->SetMinimum(0);
  gKsSBAgainstVrho->SetMaximum(1.);
  gKsSBAgainstVrho->Draw("ALP");
  sprintf(filename,"%s/%2.2d___gKsSBAgainstVrho.pdf",outdir,plot);
  canvas->SaveAs(filename);
  plot++;
  canvas->SetLogz(0);
  //_______________________________________________________________________________________________________

  // Graph of FOM of Ks against vrho cut
  gKsFOMAgainstVrho->SetLineColor(kBlue);
  gKsFOMAgainstVrho->SetLineWidth(1);
  gKsFOMAgainstVrho->SetLineStyle(1);
  gKsFOMAgainstVrho->SetMarkerColor(kBlue);
  gKsFOMAgainstVrho->SetMarkerSize(1.5);
  gKsFOMAgainstVrho->SetMarkerStyle(20);
  gKsFOMAgainstVrho->GetXaxis()->SetTitle("v_{#rho} min (cm)");
  gKsFOMAgainstVrho->GetYaxis()->SetTitle("K_{S} FOM");

  gKsFOMAgainstVrho->GetXaxis()->SetLimits(0,10);
  gKsFOMAgainstVrho->SetMinimum(0);
  gKsFOMAgainstVrho->SetMaximum(100.);
  gKsFOMAgainstVrho->Draw("ALP");
  sprintf(filename,"%s/%2.2d___gKsFOMAgainstVrho.pdf",outdir,plot);
  canvas->SaveAs(filename);
  plot++;
  canvas->SetLogz(0);
  //_______________________________________________________________________________________________________

  // vz vs M(pi+,pi-)
  canvas->SetLogz(1);
  TH1F *hvz_M2pi_copy = (TH1F*)hvz_M2pi->Clone("hvz_M2pi_copy");
  hvz_M2pi_copy->SetAxisRange(0.20,1.0);
  hvz_M2pi_copy->Draw("colz");
  line_Ks->SetY1(hvz_M2pi_copy->GetYaxis()->GetXmin());
  line_Ks->SetY2(hvz_M2pi_copy->GetYaxis()->GetXmax());
  line_Ks->Draw("same");
  sprintf(filename,"%s/%2.2d___vz_M2pi.pdf",outdir,plot);
  canvas->SaveAs(filename);
  plot++;
  canvas->SetLogz(0);
  //_______________________________________________________________________________________________________

  // vz vs M(pi+,pi-) (closeup)
  canvas->SetLogz(1);
  hvz_M2pi_copy->SetAxisRange(40,80,"Y");
  hvz_M2pi_copy->Draw("colz");
  line_Ks->SetY1(40);
  line_Ks->SetY2(80);
  line_Ks->Draw("same");
  sprintf(filename,"%s/%2.2d___vz_M2pi_close.pdf",outdir,plot);
  canvas->SaveAs(filename);
  plot++;
  canvas->SetLogz(0);
  //_______________________________________________________________________________________________________

  // Project onto M(pi+,pi-) with cut on vz
  min = 70.0; // cm
  y0 = hvz_M2pi->GetYaxis()->GetBinLowEdge(1);
  biny = hvz_M2pi->GetYaxis()->GetBinWidth(1);
  firstbin = (Int_t)((min - y0)/biny);
  lastbin  = hvz_M2pi->GetNbinsY();
  cout << "y0 = " << y0 << " biny = " << biny << endl;
  cout << "firstbin = " << firstbin << endl;
  TH1F *hM2pi_vzcut = (TH1F*)hvz_M2pi->ProjectionX("hM2pi_vzcut",firstbin, lastbin);

  // This is for plots to send to Curtis
  canvas->SetLeftMargin(0.12);
  canvas->SetBottomMargin(0.16);
  canvas->SetRightMargin(0.03);
  canvas->SetTopMargin(0.03);

  legend->Clear();
  legend->SetTextSize(0.060);
  legend->SetX1(0.800);
  legend->SetX2(1.800);
  legend->SetY1(500);
  legend->SetY2(600);

  hM2pi_vzcut->GetXaxis()->SetTitleFont(132);
  hM2pi_vzcut->GetXaxis()->SetTitleSize(0.080);
  hM2pi_vzcut->GetXaxis()->SetLabelFont(132);
  hM2pi_vzcut->GetXaxis()->SetLabelSize(0.075);
  hM2pi_vzcut->GetXaxis()->SetTitle("#pi^{+} #pi^{-} Invariant Mass (GeV/c^{2})");
  hM2pi_vzcut->GetXaxis()->CenterTitle();

  hM2pi_vzcut->GetYaxis()->SetTitleFont(132);
  hM2pi_vzcut->GetYaxis()->SetTitleSize(0.080);
  hM2pi_vzcut->GetYaxis()->SetTitleOffset(0.700);
  hM2pi_vzcut->GetYaxis()->SetLabelFont(132);
  hM2pi_vzcut->GetYaxis()->SetLabelSize(0.075);
  hM2pi_vzcut->GetYaxis()->SetLabelOffset(0.018);
  hM2pi_vzcut->SetLineColor(kBlack);
  hM2pi_vzcut->SetLineWidth(1);
  hM2pi_vzcut->GetYaxis()->SetTitle("Events");
  hM2pi_vzcut->GetYaxis()->CenterTitle();

  hM2pi_vzcut->SetAxisRange(0.2,2);
  hM2pi_vzcut->SetLineWidth(2);
  hM2pi_vzcut->Draw();
  legend->AddEntry(hM2pi_vzcut,"Data","L");

  line_Ks->SetY1(0);
  line_Ks->SetY2(hM2pi_vzcut->GetMaximum() * 1.05);
  line_Ks->SetLineColor(kMagenta);
  line_Ks->Draw("same");

  sprintf(text,"v_{z} > %5.3f cm",min);
  latex->SetTextColor(kBlack);
  latex->SetTextSize(0.080);
  // latex->DrawLatex(0.60,0.85,text);

  TF1 *fM2pi_vzcut = new TF1("fM2pi_vzcut","[0]*exp(-pow(x-[1],2.)/2./pow([2],2.)) + [3] + [4]*x + [5]* pow(x,2.)",0.4,0.6);
  fM2pi_vzcut->SetParameter(0,hM2pi_vzcut->GetMaximum());
  fM2pi_vzcut->SetParameter(1,m_KZero);
  fM2pi_vzcut->SetParameter(2,0.050);
  hM2pi_vzcut->Fit(fM2pi_vzcut,"EQRN");
  fM2pi_vzcut->SetLineColor(kRed);
  fM2pi_vzcut->SetLineStyle(2);
  fM2pi_vzcut->SetLineWidth(3);
  fM2pi_vzcut->Draw("same");

  // Plot signal Gaussian
  TF1 *fSignalM2pi_vzcut = new TF1("fSignalM2pi_vzcut","[0]*exp(-pow(x-[1],2.)/2./pow([2],2.))",0.4,0.6);
  fSignalM2pi_vzcut->SetParameter(0,fM2pi_vzcut->GetParameter(0));
  fSignalM2pi_vzcut->SetParameter(1,fM2pi_vzcut->GetParameter(1));
  fSignalM2pi_vzcut->SetParameter(2,fM2pi_vzcut->GetParameter(2));
  fSignalM2pi_vzcut->SetLineColor(kBlue);
  fSignalM2pi_vzcut->SetLineStyle(1);
  fSignalM2pi_vzcut->SetLineWidth(3);
  fSignalM2pi_vzcut->Draw("same");
  legend->AddEntry(fSignalM2pi_vzcut,"K_{S} #rightarrow #pi^{+} #pi^{-}","L");

  // Plot bg
  TF1 *fBgM2pi_vzcut = new TF1("fBgM2pi_vzcut","[0] + [1] * x + [2] * pow(x,2.)",0.4,0.6);
  fBgM2pi_vzcut->SetParameter(0,fM2pi_vzcut->GetParameter(3));
  fBgM2pi_vzcut->SetParameter(1,fM2pi_vzcut->GetParameter(4));
  fBgM2pi_vzcut->SetParameter(2,fM2pi_vzcut->GetParameter(5));
  fBgM2pi_vzcut->SetLineColor(kGreen+2);
  fBgM2pi_vzcut->SetLineStyle(2);
  fBgM2pi_vzcut->SetLineWidth(3);
  // fBgM2pi_vzcut->Draw("same");
  // legend->Draw("same");

  sprintf(text,"GlueX Commissioning: Fall 2014");
  latex->DrawLatex(0.35,0.85,text);

  sprintf(text,"K_{S} #rightarrow #pi^{+} #pi^{-}");
  latex->DrawLatex(0.35,0.75,text);

  sprintf(text,"width: %4.1f #pm %4.1f MeV",fM2pi_vzcut->GetParameter(2)* 1000.,fM2pi_vzcut->GetParError(2)* 1000.);
  latex->DrawLatex(0.50,0.65,text);

  sprintf(text,"m_{fit}   = %5.1f #pm %5.1f MeV/c^{2}",fM2pi_vzcut->GetParameter(1)* 1000, fM2pi_vzcut->GetParError(1)* 1000);
  latex->DrawLatex(0.50,0.55,text);

  sprintf(text,"m_{PDG} = %4.2f MeV/c^{2}",m_KZero * 1000.);
  latex->DrawLatex(0.50,0.45,text);

  sprintf(text,"counts: %5.1f #pm %5.1f",
	  fM2pi_vzcut->GetParameter(0) * sqrt(2. * TMath::Pi()) * fM2pi_vzcut->GetParameter(2) / hM2pi_vzcut->GetBinWidth(1),
	  fM2pi_vzcut->GetParError(0)  * sqrt(2. * TMath::Pi()) * fM2pi_vzcut->GetParameter(2) / hM2pi_vzcut->GetBinWidth(1));
  latex->DrawLatex(0.50,0.35,text);

  sprintf(filename,"%s/%2.2d___M2pi_vzcut.pdf",outdir,plot);
  canvas->SaveAs(filename);
  plot++;
  //_______________________________________________________________________________________________________

  // Get Ks yield as function of vz cut
  TGraphErrors *gKsYieldAgainstVz = new TGraphErrors();
  TGraphErrors *gKsSBAgainstVz = new TGraphErrors();
  TGraphErrors *gKsFOMAgainstVz = new TGraphErrors();
  TH1F *hM2pi_vzcut_vary;
  TF1 *fM2pi_vzcut_vary;
  for(Int_t ycutbin=0;ycutbin<hvz_M2pi->GetNbinsY();ycutbin++){
    // cout << "fitting to min bin = " << ycutbin << endl;
    // Get minimum of bin
    Double_t ymin = hvz_M2pi->GetYaxis()->GetBinLowEdge(ycutbin);
    lastbin  = hvz_M2pi->GetNbinsY();
    hM2pi_vzcut_vary = (TH1F*)hvz_M2pi->ProjectionX("hM2pi_vzcut_vary",ycutbin,lastbin);
    // Do fit
    fM2pi_vzcut_vary = new TF1("fM2pi_vzcut_vary","[0]*exp(-pow(x-[1],2.)/2./pow([2],2.)) + [3] + [4]*x + [5]* pow(x,2.)",0.4,0.6);
    fM2pi_vzcut_vary->SetParameter(0,hM2pi_vzcut_vary->GetMaximum() * 0.5);
    fM2pi_vzcut_vary->SetParameter(1,m_KZero);
    fM2pi_vzcut_vary->SetParameter(2,0.013);
    hM2pi_vzcut_vary->Fit(fM2pi_vzcut_vary,"EQRN");

    // Get yields, put in graph
    yield = fM2pi_vzcut_vary->GetParameter(0) * sqrt(2. * TMath::Pi()) * fabs(fM2pi_vzcut_vary->GetParameter(2)) / hM2pi_vzcut_vary->GetBinWidth(1);
    error = fM2pi_vzcut_vary->GetParError(0)  * sqrt(2. * TMath::Pi()) * fabs(fM2pi_vzcut_vary->GetParameter(2)) / hM2pi_vzcut_vary->GetBinWidth(1);
    gKsYieldAgainstVz->SetPoint(gKsYieldAgainstVz->GetN(),ymin,yield);
    gKsYieldAgainstVz->SetPointError(gKsYieldAgainstVz->GetN()-1,0,error);
    // cout << "fit to min bin = " << ycutbin << " ymin = " << ymin << " yield = " << yield << " +/- " << error << endl;
    
    // Set S/B
    // Defined as total yield of Gaussian divided by area under +/- 0.013 * 3 MeV of Ks
    // cout << "hM2pi_vzcut_vary->GetBinLowEdge(1) = " << hM2pi_vzcut_vary->GetBinLowEdge(1)<< " hM2pi_vzcut_vary->GetBinWidth(1) = " << hM2pi_vzcut_vary->GetBinWidth(1) << endl;
    first = (Int_t)((m_KZero - 0.013 * 3. - hM2pi_vzcut_vary->GetBinLowEdge(1)) / hM2pi_vzcut_vary->GetBinWidth(1));
    last  = (Int_t)((m_KZero + 0.013 * 3. - hM2pi_vzcut_vary->GetBinLowEdge(1)) / hM2pi_vzcut_vary->GetBinWidth(1));
    // cout << "first = " << first << " last = " << last << endl;
    area  = hM2pi_vzcut_vary->Integral(first,last);
    sob   = yield / (fM2pi_vzcut_vary->Integral(m_KZero - 0.013 * 3.,m_KZero + 0.013 * 3.) / hM2pi_vzcut_vary->GetBinWidth(1));
    error = error / (fM2pi_vzcut_vary->Integral(m_KZero - 0.013 * 3.,m_KZero + 0.013 * 3.) / hM2pi_vzcut_vary->GetBinWidth(1));
    gKsSBAgainstVz->SetPoint(gKsSBAgainstVz->GetN(),ymin,sob);
    gKsSBAgainstVz->SetPointError(gKsSBAgainstVz->GetN()-1,0,error);

    // cout << "yield = " << yield
    //      << " total = " << fM2pi_vzcut_vary->Integral(m_KZero - 0.013 * 3.,m_KZero + 0.013 * 3.) / hM2pi_vzcut_vary->GetBinWidth(1)
    //      << " bin width = " << hM2pi_vzcut_vary->GetBinWidth(1) << endl
    //      << " area from histogram = " << area << endl
    //      << " total counts of histogram = " << hM2pi_vzcut_vary->Integral() << endl
    //      << " S/B = " << sob << endl;
    
    // Set FOM
    fom   = yield / sqrt((fM2pi_vzcut_vary->Integral(m_KZero - 0.013 * 3.,m_KZero + 0.013 * 3.) / hM2pi_vzcut_vary->GetBinWidth(1)));
    gKsFOMAgainstVz->SetPoint(gKsFOMAgainstVz->GetN(),ymin,fom);

    hM2pi_vzcut_vary = 0;
    fM2pi_vzcut_vary = 0;
  }

  // Graph of yield of Ks against vz cut
  gKsYieldAgainstVz->SetLineColor(kBlue);
  gKsYieldAgainstVz->SetLineWidth(1);
  gKsYieldAgainstVz->SetLineStyle(1);
  gKsYieldAgainstVz->SetMarkerColor(kBlue);
  gKsYieldAgainstVz->SetMarkerSize(1.5);
  gKsYieldAgainstVz->SetMarkerStyle(20);
  gKsYieldAgainstVz->GetXaxis()->SetTitle("v_{z} min (cm)");
  gKsYieldAgainstVz->GetYaxis()->SetTitle("K_{S} yield");

  gKsYieldAgainstVz->GetXaxis()->SetLimits(0,90);
  gKsYieldAgainstVz->SetMinimum(0);
  gKsYieldAgainstVz->SetMaximum(5000.);
  gKsYieldAgainstVz->Draw("ALP");
  sprintf(filename,"%s/%2.2d___gKsYieldAgainstVz.pdf",outdir,plot);
  canvas->SaveAs(filename);
  plot++;
  canvas->SetLogz(0);
  //_______________________________________________________________________________________________________

  // Graph of S/B of Ks against vz cut
  gKsSBAgainstVz->SetLineColor(kBlue);
  gKsSBAgainstVz->SetLineWidth(1);
  gKsSBAgainstVz->SetLineStyle(1);
  gKsSBAgainstVz->SetMarkerColor(kBlue);
  gKsSBAgainstVz->SetMarkerSize(1.5);
  gKsSBAgainstVz->SetMarkerStyle(20);
  gKsSBAgainstVz->GetXaxis()->SetTitle("v_{z} min (cm)");
  gKsSBAgainstVz->GetYaxis()->SetTitle("K_{S} S/B");

  gKsSBAgainstVz->GetXaxis()->SetLimits(0,90);
  gKsSBAgainstVz->SetMinimum(0);
  gKsSBAgainstVz->SetMaximum(1.);
  gKsSBAgainstVz->Draw("ALP");
  sprintf(filename,"%s/%2.2d___gKsSBAgainstVz.pdf",outdir,plot);
  canvas->SaveAs(filename);
  plot++;
  canvas->SetLogz(0);
  //_______________________________________________________________________________________________________

  // Graph of FOM of Ks against vz cut
  gKsFOMAgainstVz->SetLineColor(kBlue);
  gKsFOMAgainstVz->SetLineWidth(1);
  gKsFOMAgainstVz->SetLineStyle(1);
  gKsFOMAgainstVz->SetMarkerColor(kBlue);
  gKsFOMAgainstVz->SetMarkerSize(1.5);
  gKsFOMAgainstVz->SetMarkerStyle(20);
  gKsFOMAgainstVz->GetXaxis()->SetTitle("v_{z} min (cm)");
  gKsFOMAgainstVz->GetYaxis()->SetTitle("K_{S} FOM");

  gKsFOMAgainstVz->GetXaxis()->SetLimits(0,90);
  gKsFOMAgainstVz->SetMinimum(0);
  gKsFOMAgainstVz->SetMaximum(50.);
  gKsFOMAgainstVz->Draw("ALP");
  sprintf(filename,"%s/%2.2d___gKsFOMAgainstVz.pdf",outdir,plot);
  canvas->SaveAs(filename);
  plot++;
  canvas->SetLogz(0);
  //_______________________________________________________________________________________________________

  // |p| of pi+, pi-
  hPpip_Ppim->Draw("colz");
  sprintf(filename,"%s/%2.2d___Ppip_Ppim.pdf",outdir,plot);
  canvas->SaveAs(filename);
  plot++;
  //_______________________________________________________________________________________________________

  // |p| of pi+, pi- with cuts on vrho, vz, Ks
  hPpip_Ppim_allVertexKsSelect->Draw("colz");
  sprintf(filename,"%s/%2.2d___Ppip_Ppim_allVertexKsSelect.pdf",outdir,plot);
  canvas->SaveAs(filename);
  plot++;
  //_______________________________________________________________________________________________________

  // |p| of Ks
  legend->Clear();
  hPKs->Draw();
  legend->AddEntry(hPKs,"no cuts","L");
  hPKs_allVertexKsSelect->SetLineColor(kRed);
  hPKs_allVertexKsSelect->Scale(hPKs->GetMaximum() / hPKs_allVertexKsSelect->GetMaximum() / 2.);
  hPKs_allVertexKsSelect->Draw("same");
  legend->AddEntry(hPKs_allVertexKsSelect,"cut on v_{#rho},v_{z},K_{S} mass (scaled)","L");
  legend->Draw("same");
  sprintf(filename,"%s/%2.2d___PKs.pdf",outdir,plot);
  canvas->SaveAs(filename);
  plot++;
  //_______________________________________________________________________________________________________

  // M(pi+,pi-,pi)
  legend->Clear();
  // hM3pi->SetAxisRange(0.20,2.0);
  hM3pi->Draw();
  legend->AddEntry(hM3pi,"no cuts","L");
  line_Kstar->SetY2(hM3pi->GetMaximum() * 1.05);
  line_Kstar->Draw("same");

  hM3pi_allVertexKsSelect->SetLineColor(kRed);
  hM3pi_allVertexKsSelect->Scale(hM3pi->GetMaximum() / hM3pi_allVertexKsSelect->GetMaximum());
  hM3pi_allVertexKsSelect->Draw("same");
  legend->AddEntry(hM3pi_allVertexKsSelect,"cuts on v_{#rho},v_{z},K_{S} mass (scaled)","L");

  legend->Draw("same");
  sprintf(filename,"%s/%2.2d___M3pi.pdf",outdir,plot);
  canvas->SaveAs(filename);
  plot++;
  //_______________________________________________________________________________________________________

  // M(unused p,unused pi-)
  legend->Clear();
  // hMppi->SetAxisRange(1.0,2.0);
  legend->AddEntry(hMppi,"no cuts","L");
  hMppi->Draw();
  line_Lambda->SetY2(hMppi->GetMaximum() * 1.05);
  line_Lambda->Draw("same");
  // Cut on vrho, vz, Ks mass
  hMppi_allVertexKsSelect->SetLineColor(kRed);
  hMppi_allVertexKsSelect->Scale(hMppi->GetMaximum() / hMppi_allVertexKsSelect->GetMaximum() / 2.);
  hMppi_allVertexKsSelect->Draw("same");
  legend->AddEntry(hMppi_allVertexKsSelect,"cuts on v_{#rho},v_{z},K_{S} mass (scaled)","L");
  legend->Draw("same");
  sprintf(filename,"%s/%2.2d___Mppi.pdf",outdir,plot);
  canvas->SaveAs(filename);
  plot++;
  //_______________________________________________________________________________________________________

  // MM vs M(pi+,pi-)
  canvas->Clear();
  canvas->Divide(2,2,.002,.002);
  canvas->cd(1);
  hMM_M2pi->Draw("colz");

  canvas->cd(2);
  TH1F *hMM_M2pi___projY = (TH1F*)hMM_M2pi->ProjectionY("hMM_M2pi___projY");
  hMM_M2pi___projY->Draw();
  line_proton->SetY2(hMM_M2pi___projY->GetMaximum() * 1.05);
  line_proton->Draw("same");
  line_SigmaPlus->SetY2(hMM_M2pi___projY->GetMaximum() * 1.05);
  line_SigmaPlus->Draw("same");

  canvas->cd(3);
  TH1F *hMM_M2pi___projX = (TH1F*)hMM_M2pi->ProjectionX("hMM_M2pi___projX");
  hMM_M2pi___projX->Draw();
  line_Ks->SetY2(hMM_M2pi___projX->GetMaximum() * 1.05);
  line_Ks->Draw("same");

  sprintf(filename,"%s/%2.2d___MM_M2pi.pdf",outdir,plot);
  canvas->SaveAs(filename);
  plot++;
  canvas->Clear();
  //_______________________________________________________________________________________________________

  // MM vs M(pi+,pi-) with cut on vrho, vz, Ks mass
  canvas->Clear();
  canvas->Divide(2,2,.002,.002);
  canvas->cd(1);
  hMM_M2pi_allVertexKsSelect->Draw("colz");

  canvas->cd(2);
  TH1F *hMM_M2pi_allVertexKsSelect___projY = (TH1F*)hMM_M2pi_allVertexKsSelect->ProjectionY("hMM_M2pi_allVertexKsSelect___projY");
  hMM_M2pi_allVertexKsSelect___projY->Draw();
  line_proton->SetY2(hMM_M2pi_allVertexKsSelect___projY->GetMaximum() * 1.05);
  line_proton->Draw("same");
  line_SigmaPlus->SetY2(hMM_M2pi_allVertexKsSelect___projY->GetMaximum() * 1.05);
  line_SigmaPlus->Draw("same");

  canvas->cd(3);
  TH1F *hMM_M2pi_allVertexKsSelect___projX = (TH1F*)hMM_M2pi_allVertexKsSelect->ProjectionX("hMM_M2pi_allVertexKsSelect___projX");
  hMM_M2pi_allVertexKsSelect___projX->Draw();
  line_Ks->SetY2(hMM_M2pi_allVertexKsSelect___projX->GetMaximum() * 1.05);
  line_Ks->Draw("same");

  sprintf(filename,"%s/%2.2d___MM_M2pi_allVertexKsSelecty.pdf",outdir,plot);
  canvas->SaveAs(filename);
  plot++;
  canvas->Clear();
  //_______________________________________________________________________________________________________

  // dE/dx pi+ no cuts
  canvas->Divide(2,2,.002,.002);
  canvas->cd(1);
  hdEdx_pip_CDC->Draw("colz");
  latex->DrawLatex(0.50,0.60,"#pi^{+} CDC");
  canvas->cd(2);
  gPad->SetLogz(1);
  hdEdx_pip_FDC->Draw("colz");
  latex->DrawLatex(0.50,0.60,"#pi^{+} FDC");
  canvas->cd(3);
  hdEdx_pip_ST->Draw("colz");
  latex->DrawLatex(0.50,0.60,"#pi^{+} ST");
  sprintf(filename,"%s/%2.2d___dEdx_pip.pdf",outdir,plot);
  canvas->SaveAs(filename);
  plot++;
  canvas->Clear();
  //_______________________________________________________________________________________________________

  // dE/dx pi- no cuts
  canvas->Divide(2,2,.002,.002);
  canvas->cd(1);
  hdEdx_pim_CDC->Draw("colz");
  latex->DrawLatex(0.50,0.60,"#pi^{-} CDC");
  canvas->cd(2);
  gPad->SetLogz(1);
  hdEdx_pim_FDC->Draw("colz");
  latex->DrawLatex(0.50,0.60,"#pi^{-} FDC");
  canvas->cd(3);
  hdEdx_pim_ST->Draw("colz");
  latex->DrawLatex(0.50,0.60,"#pi^{-} ST");
  sprintf(filename,"%s/%2.2d___dEdx_pim.pdf",outdir,plot);
  canvas->SaveAs(filename);
  plot++;
  canvas->Clear();
  //_______________________________________________________________________________________________________

  // dE/dx pi+ Ks select
  canvas->Divide(2,2,.002,.002);
  canvas->cd(1);
  hdEdx_pip_CDC_selectKs->Draw("colz");
  latex->DrawLatex(0.50,0.60,"#pi^{+} CDC");
  canvas->cd(2);
  gPad->SetLogz(1);
  hdEdx_pip_FDC_selectKs->Draw("colz");
  latex->DrawLatex(0.50,0.60,"#pi^{+} FDC");
  canvas->cd(3);
  hdEdx_pip_ST_selectKs->Draw("colz");
  latex->DrawLatex(0.50,0.60,"#pi^{+} ST");
  sprintf(filename,"%s/%2.2d___dEdx_pip_selectKs.pdf",outdir,plot);
  canvas->SaveAs(filename);
  plot++;
  canvas->Clear();
  //_______________________________________________________________________________________________________

  // dE/dx pi- Ks select
  canvas->Divide(2,2,.002,.002);
  canvas->cd(1);
  hdEdx_pim_CDC_selectKs->Draw("colz");
  latex->DrawLatex(0.50,0.60,"#pi^{-} CDC");
  canvas->cd(2);
  gPad->SetLogz(1);
  hdEdx_pim_FDC_selectKs->Draw("colz");
  latex->DrawLatex(0.50,0.60,"#pi^{-} FDC");
  canvas->cd(3);
  hdEdx_pim_ST_selectKs->Draw("colz");
  latex->DrawLatex(0.50,0.60,"#pi^{-} ST");
  sprintf(filename,"%s/%2.2d___dEdx_pim_selectKs.pdf",outdir,plot);
  canvas->SaveAs(filename);
  plot++;
  canvas->Clear();
  //_______________________________________________________________________________________________________

  // dE/dx unused K+ no cuts
  canvas->Divide(2,2,.002,.002);
  canvas->cd(1);
  hdEdx_kp_CDC->Draw("colz");
  latex->DrawLatex(0.50,0.60,"unused K^{+} CDC");
  canvas->cd(2);
  gPad->SetLogz(1);
  hdEdx_kp_FDC->Draw("colz");
  latex->DrawLatex(0.50,0.60,"unused K^{+} FDC");
  canvas->cd(3);
  hdEdx_kp_ST->Draw("colz");
  latex->DrawLatex(0.50,0.60,"unused K^{+} ST");
  sprintf(filename,"%s/%2.2d___dEdx_kp.pdf",outdir,plot);
  canvas->SaveAs(filename);
  plot++;
  canvas->Clear();
  //_______________________________________________________________________________________________________

  // dE/dx unused K- no cuts
  canvas->Divide(2,2,.002,.002);
  canvas->cd(1);
  hdEdx_km_CDC->Draw("colz");
  latex->DrawLatex(0.50,0.60,"unused K^{-} CDC");
  canvas->cd(2);
  gPad->SetLogz(1);
  hdEdx_km_FDC->Draw("colz");
  latex->DrawLatex(0.50,0.60,"unused K^{-} FDC");
  canvas->cd(3);
  hdEdx_km_ST->Draw("colz");
  latex->DrawLatex(0.50,0.60,"unused K^{-} ST");
  sprintf(filename,"%s/%2.2d___dEdx_km.pdf",outdir,plot);
  canvas->SaveAs(filename);
  plot++;
  canvas->Clear();
  //_______________________________________________________________________________________________________

  // dE/dx unused K+ Lamba select
  canvas->Divide(2,2,.002,.002);
  canvas->cd(1);
  hdEdx_kp_CDC_selectKs->Draw("colz");
  latex->DrawLatex(0.50,0.60,"unused K^{+} CDC");
  canvas->cd(2);
  gPad->SetLogz(1);
  hdEdx_kp_FDC_selectKs->Draw("colz");
  latex->DrawLatex(0.50,0.60,"unused K^{+} FDC");
  canvas->cd(3);
  hdEdx_kp_ST_selectKs->Draw("colz");
  latex->DrawLatex(0.50,0.60,"unused K^{+} ST");
  sprintf(filename,"%s/%2.2d___dEdx_kp_selectKs.pdf",outdir,plot);
  canvas->SaveAs(filename);
  plot++;
  canvas->Clear();
  //_______________________________________________________________________________________________________

  // dE/dx unused K- Lamba select
  canvas->Divide(2,2,.002,.002);
  canvas->cd(1);
  hdEdx_km_CDC_selectKs->Draw("colz");
  latex->DrawLatex(0.50,0.60,"unused K^{-} CDC");
  canvas->cd(2);
  gPad->SetLogz(1);
  hdEdx_km_FDC_selectKs->Draw("colz");
  latex->DrawLatex(0.50,0.60,"unused K^{-} FDC");
  canvas->cd(3);
  hdEdx_km_ST_selectKs->Draw("colz");
  latex->DrawLatex(0.50,0.60,"unused K^{-} ST");
  sprintf(filename,"%s/%2.2d___dEdx_km_selectKs.pdf",outdir,plot);
  canvas->SaveAs(filename);
  plot++;
  canvas->Clear();
  //_______________________________________________________________________________________________________

  // vertex for pi+
  canvas->Divide(2,2,.002,.002);
  canvas->cd(1);
  gPad->SetLogz(1);
  hvy_vx_pip_measured->Draw("colz");
  latex->DrawLatex(0.70,0.85,"measured #pi^{+}");
  canvas->cd(2);
  gPad->SetLogz(1);
  hvy_vx_pip->Draw("colz");
  latex->DrawLatex(0.70,0.85,"fitted #pi^{+}");
  canvas->cd(3);
  legend->Clear();
  hvz_pip_measured->SetAxisRange(40,100);
  hvz_pip_measured->SetMaximum(hvz_pip->GetMaximum() * 1.05);
  hvz_pip_measured->Draw();
  legend->AddEntry(hvz_pip_measured,"measured #pi^{+}","L");
  hvz_pip->SetLineColor(kRed);
  hvz_pip->Draw("same");
  legend->AddEntry(hvz_pip,"fitted #pi^{+}","L");
  legend->Draw("same");
  canvas->cd(4);
  canvas->cd(4)->SetLogy(1);
  legend->Clear();
  hvz_pip_measured->SetMinimum(0.5);
  hvz_pip_measured->Draw();
  legend->AddEntry(hvz_pip_measured,"measured #pi^{+}","L");
  hvz_pip->SetLineColor(kRed);
  hvz_pip->Draw("same");
  legend->AddEntry(hvz_pip,"fitted #pi^{+}","L");
  legend->Draw("same");

  sprintf(filename,"%s/%2.2d___vertex_pip.pdf",outdir,plot);
  canvas->SaveAs(filename);
  plot++;
  canvas->Clear();
  //_______________________________________________________________________________________________________

  // vertex for pi-
  canvas->Divide(2,2,.002,.002);
  canvas->cd(1);
  gPad->SetLogz(1);
  hvy_vx_pim_measured->Draw("colz");
  latex->DrawLatex(0.70,0.85,"measured #pi^{-}");
  canvas->cd(2);
  gPad->SetLogz(1);
  hvy_vx_pim->Draw("colz");
  latex->DrawLatex(0.70,0.85,"fitted #pi^{-}");
  canvas->cd(3);
  legend->Clear();
  hvz_pim_measured->SetAxisRange(40,100);
  hvz_pim_measured->SetMaximum(hvz_pim->GetMaximum() * 1.05);
  hvz_pim_measured->Draw();
  legend->AddEntry(hvz_pim_measured,"measured #pi^{-}","L");
  hvz_pim->SetLineColor(kRed);
  hvz_pim->Draw("same");
  legend->AddEntry(hvz_pim,"fitted #pi^{-}","L");
  legend->Draw("same");
  canvas->cd(4);
  canvas->cd(4)->SetLogy(1);
  legend->Clear();
  hvz_pim_measured->SetMinimum(0.5);
  hvz_pim_measured->Draw();
  legend->AddEntry(hvz_pim_measured,"measured #pi^{-}","L");
  hvz_pim->SetLineColor(kRed);
  hvz_pim->Draw("same");
  legend->AddEntry(hvz_pim,"fitted #pi^{-}","L");
  legend->Draw("same");

  sprintf(filename,"%s/%2.2d___vertex_pim.pdf",outdir,plot);
  canvas->SaveAs(filename);
  plot++;
  canvas->Clear();
  //_______________________________________________________________________________________________________

  // vertex of pi+ vs pi-
  gPad->SetLogz(1);
  hvz_pip_pim_measured->Draw("colz");
  sprintf(filename,"%s/%2.2d___vz_pip_pim_measured.pdf",outdir,plot);
  canvas->SaveAs(filename);
  plot++;
  canvas->Clear();
  //_______________________________________________________________________________________________________

  // vertex of pi+ vs pi-
  gPad->SetLogz(1);
  hvz_pip_pim->Draw("colz");
  sprintf(filename,"%s/%2.2d___vz_pip_pim.pdf",outdir,plot);
  canvas->SaveAs(filename);
  plot++;
  canvas->Clear();
  //_______________________________________________________________________________________________________

  return 0;
}

void printUsage(){
      cout << "Usage: plotter -i [infile name] " << endl;
      cout << "Options:" << endl;
      cout << "\t-i [infile name]" << endl;
      // cout << "\t-c create plots for each channel" << endl;
}