void draw_eff();
void draw_ptot();
void draw_hit_q();
void draw_clustersize();
void draw_reco_clustersize();
void draw_cluster_res();

void
mc_eval_plots()
{
  std::cout << "Try calling 'draw_all' for all plots" << std::endl;
}

void 
draw_all()
{
  draw_dndetach();
  draw_ptot();
  draw_hit_q();
  draw_clustersize();
  draw_reco_clustersize();
  draw_cluster_res();
  draw_dca_r();
  draw_dca_phi();
  draw_eff();
}

void
draw_dndetach()
{
  TTree* t = gROOT->FindObject("fkin_eval");
  if ( ! t ) { 
    std::cout << "Failed to find fkin_eval tree" << std::endl;
    return;
  }

  TTree* e = gROOT->FindObject("header_eval");
  if ( ! e ) { 
    std::cout << "Failed to find header_eval tree" << std::endl;
    return;
  }

  int nevts = e->GetEntries();

  TH1F* h_z = new TH1F("h_z","h_z",100,-100.0,100.0);
  TH1F* h_z1 = new TH1F("h_z1","h_z1",100,-100.0,100.0);

  TH1F* h_dndeta = new TH1F("h_dndeta","dn_{ch}/d#eta",80,-8.0,8.0);
  TH1F* h_dndeta1 = new TH1F("h_dndeta1","dn_{ch}/d#eta",80,-8.0,8.0);
  TH1F* h_dndeta2 = new TH1F("h_dndeta2","dn_{ch}/d#eta",80,-8.0,8.0);

  TCut ch_cut = "(idpart==2||idpart==3||idpart==5||idpart==6||idpart==8||idpart==9||idpart==11||idpart==12||idpart==14||idpart==15)";
  TCut v_cut = "(abs(z_vertex)<0.0001)";
  TCut acc_cut = "(nstationsHit>3)";
  TCut parent = "(itparent<0)";

  t->Project("h_dndeta","-log(tan(pthet/2*3.14159/180.0))",ch_cut);
  //t->Project("h_dndeta1","-log(tan(pthet/2*3.14159/180.0))","abs(z_vertex)<0.0001","same");
  t->Project("h_dndeta1","-log(tan(pthet/2*3.14159/180.0))",ch_cut && v_cut,"same");
  //t->Project("h_dndeta2","-log(tan(pthet/2*3.14159/180.0))","nstationsHit>3","same");
  t->Project("h_dndeta2","-log(tan(pthet/2*3.14159/180.0))",ch_cut && acc_cut,"same");

  t->Project("h_z","z_vertex",ch_cut);
  t->Project("h_z1","z_vertex",ch_cut&&parent);

  double deta = h_dndeta->GetXaxis()->GetBinWidth(1);

  h_dndeta->Scale(1.0/nevts/deta);
  h_dndeta1->Scale(1.0/nevts/deta);
  h_dndeta2->Scale(1.0/nevts/deta);

  TCanvas* c = gROOT->FindObject("dndetaCanvas");
  if ( ! c ) { 
    std::cout << "Creating dndetaCanvas" << std::endl;
    c = new TCanvas("dndetaCanvas","dndetaCanvas",1000,600);
  }
  c->Clear();
  c->Divide(2,1);
  
  c->cd(1);
  gPad->SetLogy(true);
  gPad->SetBottomMargin(0.12);
  gPad->SetLeftMargin(0.12);
  gPad->SetRightMargin(0.05);
  
  TLegend* legend0 = new TLegend(0.6,0.70,0.99,0.88,"","NDC");

  h_z->SetTitle("");
  h_z->SetStats(false);
  h_z->GetXaxis()->SetTitle("Z_{vtx} (cm)");
  h_z->GetXaxis()->SetTitleSize(0.05);
  h_z->GetYaxis()->SetTitle("Counts");
  h_z->GetYaxis()->SetTitleSize(0.05);
  //h_z->GetYaxis()->SetTitleOffset(1.1);
  //h_z->SetLineWidth(3);
  h_z->Draw();
  legend0->AddEntry(h_z,"All","F,BOX");

  //h_z1->SetFillColor(11);
  h_z1->SetLineColor(2);
  h_z1->SetFillColor(2);
  h_z1->SetFillStyle(3003);
  h_z1->SetLineWidth(3);
  h_z1->Draw("same");
  legend0->AddEntry(h_z1,"ITPARENT<0","F,BOX");

  legend0->Draw();

  TLegend* legend = new TLegend(0.75,0.70,0.99,0.95,"","NDC");

  c->cd(2);
  gPad->SetBottomMargin(0.12);
  gPad->SetLeftMargin(0.12);
  gPad->SetRightMargin(0.05);

  h_dndeta->SetTitle("");
  h_dndeta->SetStats(false);
  h_dndeta->GetXaxis()->SetTitle("#eta");
  h_dndeta->GetXaxis()->SetTitleSize(0.05);
  h_dndeta->GetYaxis()->SetTitle("dn_{ch}/d#eta");
  h_dndeta->GetYaxis()->SetTitleSize(0.05);
  h_dndeta->GetYaxis()->SetTitleOffset(1.2);
  h_dndeta->Draw();
  legend->AddEntry(h_dndeta,"All","F,BOX");

  h_dndeta1->SetLineColor(2);
  h_dndeta1->SetFillColor(2);
  h_dndeta1->SetFillStyle(3003);
  h_dndeta1->SetLineWidth(3);
  //h_dndeta1->SetFillColor(11);
  h_dndeta1->Draw("same");
  legend->AddEntry(h_dndeta1,"|Zvtx|<0.01","F,BOX");

  h_dndeta2->SetFillColor(13);
  h_dndeta2->SetLineWidth(3);
  h_dndeta2->Draw("same");
  legend->AddEntry(h_dndeta2,"4 Sta hit","F,BOX");

  legend->Draw();

  return;
}

void
draw_ptot()
{
  double z_vtx = 0.0; // nominal vtx position

  std::string tname = "match_eval";
  TTree* t = gROOT->FindObject(tname.c_str());
  if ( ! t ) { 
    std::cout << "Failed to find "  << tname << " tree" << std::endl;
    return;
  }

  char vz_sel[1024];
  sprintf(vz_sel,"arm==1&&abs(vtx_orig[2]-%f)<0.0001",z_vtx);

  int nbins = 50;
  double pmax = 15.0;
  TH1F* hptot = new TH1F("hptot","North Arm MC Tracks",nbins,0,pmax);
  TH1F* hptot_prim = new TH1F("hptot_prim","Momentum spectrum",nbins,0,pmax);
  TH1F* hptot_4sta = new TH1F("hptot_4sta","Momentum spectrum",nbins,0,pmax);
  TH1F* hptot_prim4sta = new TH1F("hptot_prim4sta","Momentum spectrum",nbins,0,pmax);

  t->Project("hptot","ptot_orig","arm==1");
  t->Project("hptot_prim","ptot_orig",vz_sel);

  sprintf(vz_sel,"arm==1&&abs(vtx_orig[2]-%f)<0.0001&&nstaHit>3",z_vtx);
  t->Project("hptot_prim4sta","ptot_orig",vz_sel);

  sprintf(vz_sel,"arm==1&&nstaHit>3",z_vtx);
  t->Project("hptot_4sta","ptot_orig",vz_sel);

  std::string cname = "ptotCanvas";
  TCanvas* c = gROOT->FindObject(cname.c_str());
  if ( ! c ) { 
    std::cout << "Creating " << cname << std::endl;
    c = new TCanvas(cname.c_str(),cname.c_str());
  }
  c->Clear();
  c->SetLogy();
  gPad->SetBottomMargin(0.12);
  gPad->SetLeftMargin(0.12);
  gPad->SetRightMargin(0.05);

  TLegend* legend0 = new TLegend(0.6,0.40,0.97,0.8,"","NDC");
  
  hptot->GetXaxis()->SetTitle("MC momentum (GeV/c)");
  hptot->GetXaxis()->SetTitleSize(0.05);
  hptot->GetXaxis()->SetLabelSize(0.06);
  hptot->GetYaxis()->SetTitle("counts");
  hptot->GetYaxis()->SetTitleSize(0.05);
  hptot->GetXaxis()->SetTitleOffset(1.1);
  hptot->GetYaxis()->SetTitleOffset(1.1);
  hptot->GetYaxis()->SetLabelSize(0.06);
  hptot->SetMarkerStyle(20);
  hptot->SetLineWidth(3);
  hptot->Draw("e");
  legend0->AddEntry(hptot,"All tracks","PL,BOX");
  
  hptot_prim->SetLineColor(2);
  hptot_prim->SetFillColor(2);
  hptot_prim->SetFillStyle(3003);
  hptot_prim->SetLineWidth(3);
  hptot_prim->Draw("same");
  //hptot->Draw("e,same");
  legend0->AddEntry(hptot_prim,"#DeltaZ<0.0001","F,BOX");
  
  hptot_4sta->SetLineColor(6);
  hptot_4sta->SetFillColor(6);
  hptot_4sta->SetFillStyle(3003);
  hptot_4sta->SetLineWidth(3);
  hptot_4sta->Draw("same");
  legend0->AddEntry(hptot_4sta,"4 Sta hit","F,BOX");

  hptot_prim4sta->SetLineColor(4);
  hptot_prim4sta->SetFillColor(4);
  hptot_prim4sta->SetFillStyle(3003);
  hptot_prim4sta->SetLineWidth(3);
  hptot_prim4sta->Draw("same");
  hptot->Draw("e,same");
  legend0->AddEntry(hptot_prim4sta,"#DeltaZ<0.0001 && 4 Sta","F,BOX");

  legend0->Draw();

  return;
}

void
draw_hit_q()
{
  std::string tname = "mc_fvtx_hit";
  TTree* t = gROOT->FindObject(tname.c_str());
  if ( ! t ) { 
    std::cout << "Failed to find "  << tname << " tree" << std::endl;
    return;
  }

  TH1F* hitq = new TH1F("hitq","Hit Charge Distribution",100,0.0,50000.);

  t->Project("hitq","qhit","");

  std::string cname = "hitqCanvas";
  TCanvas* c = gROOT->FindObject(cname.c_str());
  if ( ! c ) { 
    std::cout << "Creating " << cname << std::endl;
    c = new TCanvas(cname.c_str(),cname.c_str(),533,400);
  }
  c->Clear();
  c->SetLogy();
  gPad->SetBottomMargin(0.12);
  gPad->SetLeftMargin(0.12);
  gPad->SetRightMargin(0.05);

  hitq->GetXaxis()->SetTitle("Hit Charge");
  hitq->GetXaxis()->SetTitleSize(0.05);
  hitq->GetXaxis()->SetLabelSize(0.03);
  hitq->GetYaxis()->SetTitle("counts");
  hitq->GetYaxis()->SetTitleSize(0.05);
  hitq->GetXaxis()->SetTitleOffset(1.1);
  hitq->GetYaxis()->SetTitleOffset(1.1);
  hitq->GetYaxis()->SetLabelSize(0.06);
  hitq->SetLineWidth(2);
  std::cout << "draw" << std::endl;
  hitq->Draw();

  gPad->Update();
//  ((TPaveStats*)hitq->GetFunction("stats"))->SetX1NDC(0.7);
//  ((TPaveStats*)hitq->GetFunction("stats"))->SetY1NDC(0.7);
  gPad->Modified();

  return;
}

void
draw_clustersize()
{
  std::string tname = "mc_radiograph";
  TTree* t = gROOT->FindObject(tname.c_str());
  if ( ! t ) { 
    std::cout << "Failed to find "  << tname << " tree" << std::endl;
    return;
  }

  TH1F* hclust = new TH1F("hclust","North Arm Cluster Sizes",61,-0.5,60.5);
  //TH1F* hclust_prim = new TH1F("hclust_prim","North Arm Cluster Sizes (Primaries)",26,-0.5,25.5);

  t->Project("hclust","size","arm==1");
  //t->Project("hptot_prim","ptot_orig","arm==1&&itparent<0");

  std::string cname = "clusCanvas";
  TCanvas* c = gROOT->FindObject(cname.c_str());
  if ( ! c ) { 
    std::cout << "Creating " << cname << std::endl;
    c = new TCanvas(cname.c_str(),cname.c_str(),533,400);
  }
  c->Clear();
  c->SetLogy();
  gPad->SetBottomMargin(0.12);
  gPad->SetLeftMargin(0.12);
  gPad->SetRightMargin(0.05);

  hclust->GetXaxis()->SetTitle("Cluster size");
  hclust->GetXaxis()->SetTitleSize(0.05);
  hclust->GetXaxis()->SetLabelSize(0.06);
  hclust->GetYaxis()->SetTitle("counts");
  hclust->GetYaxis()->SetTitleSize(0.05);
  hclust->GetXaxis()->SetTitleOffset(1.1);
  hclust->GetYaxis()->SetTitleOffset(1.1);
  hclust->GetYaxis()->SetLabelSize(0.06);
  //hclust->SetMarkerStyle(20);
  hclust->SetLineWidth(2);
  std::cout << "draw" << std::endl;
  hclust->Draw();

  gPad->Update();
//  ((TPaveStats*)hclust->GetFunction("stats"))->SetX1NDC(0.7);
//  ((TPaveStats*)hclust->GetFunction("stats"))->SetY1NDC(0.7);
  gPad->Modified();

  return;
}

void
draw_reco_clustersize()
{
  std::string tname = "mc_cluster";
  TTree* t = gROOT->FindObject(tname.c_str());
  if ( ! t ) { 
    std::cout << "Failed to find "  << tname << " tree" << std::endl;
    return;
  }

  TH1F* hclust2 = new TH1F("hclust2","North Arm Cluster Sizes",31,-0.5,30.5);

  t->Project("hclust2","size","");
  //t->Project("hptot_prim","ptot_orig","arm==1&&itparent<0");

  std::string cname = "RecoclusCanvas";
  TCanvas* c = gROOT->FindObject(cname.c_str());
  if ( ! c ) { 
    std::cout << "Creating " << cname << std::endl;
    c = new TCanvas(cname.c_str(),cname.c_str(),533,400);
  }
  c->Clear();
  c->SetLogy();
  gPad->SetBottomMargin(0.12);
  gPad->SetLeftMargin(0.12);
  gPad->SetRightMargin(0.05);

  hclust2->GetXaxis()->SetTitle("Reconstructed Cluster size");
  hclust2->GetXaxis()->SetTitleSize(0.05);
  hclust2->GetXaxis()->SetLabelSize(0.06);
  hclust2->GetYaxis()->SetTitle("counts");
  hclust2->GetYaxis()->SetTitleSize(0.05);
  hclust2->GetXaxis()->SetTitleOffset(1.1);
  hclust2->GetYaxis()->SetTitleOffset(1.1);
  hclust2->GetYaxis()->SetLabelSize(0.06);
  //hclust2->SetMarkerStyle(20);
  hclust2->SetLineWidth(2);
  std::cout << "draw" << std::endl;
  hclust2->Draw();

  gPad->Update();
//  ((TPaveStats*)hclust2->GetFunction("stats"))->SetX1NDC(0.7);
//  ((TPaveStats*)hclust2->GetFunction("stats"))->SetY1NDC(0.7);
  gPad->Modified();

  return;
}

void
draw_cluster_res()
{
  std::string tname = "mc_cluster";
  TTree* t = gROOT->FindObject(tname.c_str());
  if ( ! t ) { 
    std::cout << "Failed to find "  << tname << " tree" << std::endl;
    return;
  }

  TH1F* clusres = new TH1F("clusres","North Arm Cluster Sizes",50,0.0,0.005);

  t->Project("clusres","res","");
  //t->Project("hptot_prim","ptot_orig","arm==1&&itparent<0");

  std::string cname = "ClusResCanvas";
  TCanvas* c = gROOT->FindObject(cname.c_str());
  if ( ! c ) { 
    std::cout << "Creating " << cname << std::endl;
    c = new TCanvas(cname.c_str(),cname.c_str(),533,400);
  }
  c->Clear();
  c->SetLogy();
  gStyle->SetOptStat(1110); // RMS, mean, NumEntries, no hist name
  gPad->SetBottomMargin(0.12);
  gPad->SetLeftMargin(0.12);
  gPad->SetRightMargin(0.05);

  clusres->SetStats(true);
  clusres->GetXaxis()->SetTitle("Coordinate Resolution (cm)");
  clusres->GetXaxis()->SetTitleSize(0.05);
  clusres->GetXaxis()->SetLabelSize(0.04);
  clusres->GetYaxis()->SetTitle("counts");
  clusres->GetYaxis()->SetTitleSize(0.05);
  clusres->GetXaxis()->SetTitleOffset(1.1);
  clusres->GetYaxis()->SetTitleOffset(1.1);
  clusres->GetYaxis()->SetLabelSize(0.06);
  //clusres->SetMarkerStyle(20);
  clusres->SetLineWidth(2);
  std::cout << "draw" << std::endl;
  clusres->Draw();

  gPad->Update();
//  ((TPaveStats*)clusres->GetFunction("stats"))->SetX1NDC(0.7);
//  ((TPaveStats*)clusres->GetFunction("stats"))->SetY1NDC(0.7);
  gPad->Modified();

  return;
}

void
draw_dca_r()
{
  std::string tname2 = "mc_trk_eval";
  TTree* t2 = gROOT->FindObject(tname2.c_str());
  if ( ! t2 ) { 
    std::cout << "Failed to find "  << tname2 << " tree" << std::endl;
    return;
  }

  std::string cname2 = "DCArCanvas";
  TCanvas* c2 = gROOT->FindObject(cname2.c_str());
  if ( ! c2 ) { 
    std::cout << "Creating " << cname2 << std::endl;
    c2 = new TCanvas(cname2.c_str(),cname2.c_str(),533,400);
  }
  c2->Clear();
  gStyle->SetOptStat(1110); // RMS, mean, NumEntries, no hist name
  gPad->SetBottomMargin(0.12);
  gPad->SetLeftMargin(0.12);
  gPad->SetRightMargin(0.05);

  std::cout << "draw" << std::endl;
m1=new TH1F("m1","DCA resolution vs. Momentum",10,1,15);
e1=new TProfile("e1","DCA resolution vs. Momentum",10,1,15,"s");
TCut rres="abs((x0reco-x0mc)*(px0mc/sqrt(px0mc**2+py0mc**2))+(y0reco-y0mc)*(py0mc/sqrt(px0mc**2+py0mc**2)))<.1";
mc_trk_eval->Draw("(x0reco-x0mc)*(px0mc/sqrt(px0mc**2+py0mc**2))+(y0reco-y0mc)*(py0mc/sqrt(px0mc**2+py0mc**2)):sqrt(px0mc**2+py0mc**2+pz0mc**2)>>e1",rres&&"abs(x0reco)>0&&abs(y0reco)>0&&abs(z0mc)<10");
m1->Sumw2();
m1->SetBinContent(0,e1->GetBinError(0));
m1->SetBinContent(1,e1->GetBinError(1));
m1->SetBinContent(2,e1->GetBinError(2));
m1->SetBinContent(3,e1->GetBinError(3));
m1->SetBinContent(4,e1->GetBinError(4));
m1->SetBinContent(5,e1->GetBinError(5));
m1->SetBinContent(6,e1->GetBinError(6));
m1->SetBinContent(7,e1->GetBinError(7));
m1->SetBinContent(8,e1->GetBinError(8));
m1->SetBinContent(9,e1->GetBinError(9));
m1->SetBinContent(10,e1->GetBinError(10));
m1->SetMarkerStyle(2);
m1->SetMarkerColor(4);
m1->SetMinimum(0.0);
m1->GetXaxis()->SetTitle("Momentum (GeV/c)");
m1->GetYaxis()->SetTitle("DCA_r Sigma (cm)");
m1->GetYaxis()->SetTitleOffset(1.4);
m1->GetYaxis()->SetTitleSize(0.035);
m1->GetYaxis()->SetLabelSize(0.035);
m1->SetMaximum(0.050);
  m1->Draw("p");

  gPad->Update();
//  ((TPaveStats*)clusres->GetFunction("stats"))->SetX1NDC(0.7);
//  ((TPaveStats*)clusres->GetFunction("stats"))->SetY1NDC(0.7);
//  gPad->Modified();

  return;
}

void
draw_dca_phi()
{
  std::string tname = "mc_trk_eval";
  TTree* t = gROOT->FindObject(tname.c_str());
  if ( ! t ) { 
    std::cout << "Failed to find "  << tname << " tree" << std::endl;
    return;
  }

  std::string cname = "DCAPhiCanvas";
  TCanvas* c = gROOT->FindObject(cname.c_str());
  if ( ! c ) { 
    std::cout << "Creating " << cname << std::endl;
    c = new TCanvas(cname.c_str(),cname.c_str(),533,400);
  }
  c->Clear();
  gStyle->SetOptStat(1110); // RMS, mean, NumEntries, no hist name
  gPad->SetBottomMargin(0.12);
  gPad->SetLeftMargin(0.12);
  gPad->SetRightMargin(0.05);

  std::cout << "draw" << std::endl;
m2=new TH1F("m2","DCA Phi resolution vs. Momentum",10,1,15);
e2=new TProfile("e2","DCA Phi resolution vs. Momentum",10,1,15,"s");
TCut phires="abs((y0reco-y0mc)*(px0mc/sqrt(px0mc**2+py0mc**2))-(x0reco-x0mc)*(py0mc/sqrt(px0mc**2+py0mc**2)))<.1";
mc_trk_eval->Draw("(y0reco-y0mc)*(px0mc/sqrt(px0mc**2+py0mc**2))-(x0reco-x0mc)*(py0mc/sqrt(px0mc**2+py0mc**2)):sqrt(px0mc**2+py0mc**2+pz0mc**2)>>e2",phires&&"abs(x0reco)>0&&abs(y0reco)>0&&abs(z0mc)<10");
m2->Sumw2();
m2->SetBinContent(0,e2->GetBinError(0));
m2->SetBinContent(1,e2->GetBinError(1));
m2->SetBinContent(2,e2->GetBinError(2));
m2->SetBinContent(3,e2->GetBinError(3));
m2->SetBinContent(4,e2->GetBinError(4));
m2->SetBinContent(5,e2->GetBinError(5));
m2->SetBinContent(6,e2->GetBinError(6));
m2->SetBinContent(7,e2->GetBinError(7));
m2->SetBinContent(8,e2->GetBinError(8));
m2->SetBinContent(9,e2->GetBinError(9));
m2->SetMarkerStyle(2);
m2->SetMarkerColor(2);
m2->SetMinimum(0.0);
m2->GetXaxis()->SetTitle("Momentum (GeV/c)");
m2->GetYaxis()->SetTitle("DCA_phi Sigma (cm)");
 m2->GetYaxis()->SetTitleOffset(1.0);   //1.4
 m2->GetYaxis()->SetTitleSize(0.050);    // 0.035
m2->GetYaxis()->SetLabelSize(0.035);
m2->SetMaximum(0.08);

  m2->Draw("p");

  gPad->Update();
//  ((TPaveStats*)clusres->GetFunction("stats"))->SetX1NDC(0.7);
//  ((TPaveStats*)clusres->GetFunction("stats"))->SetY1NDC(0.7);
//  gPad->Modified();

  return;
}

void
draw_resolution()
{
  int nbins = 100;
  double maxz =  0.5;

  double maxr =  1.0;

  std::string tname = "match_eval";
  TTree* t = gROOT->FindObject(tname.c_str());
  if ( ! t ) { 
    std::cout << "Failed to find "  << tname << " tree" << std::endl;
    return;
  }
 
  delete gROOT->FindObject("hz0");
  delete gROOT->FindObject("hz01d");
  delete gROOT->FindObject("hr0");

  TH1F* hz01d = new TH1F("hz01d","R(Z=0)",nbins,-maxz,maxz);

  TH2F* hz0 = new TH2F("hz0","R(Z=0) vs. p",100,0.0,10.0,nbins,-maxz,maxz);
  hz0->GetXaxis()->SetTitle("MC p (GeV/c)");
  hz0->GetYaxis()->SetTitle("R(Z=0) (cm)");
  hz0->GetXaxis()->SetTitleSize(0.05);
  hz0->GetYaxis()->SetTitleSize(0.05);
  hz0->GetXaxis()->SetTitleOffset(1.0);
  hz0->GetYaxis()->SetTitleOffset(1.1);

  TH2F* hr0 = new TH2F("hr0","Z(R=0)",100,0.0,10.0,nbins,-maxr,maxr);
  hr0->GetXaxis()->SetTitle("MC p (GeV/c)");
  hr0->GetYaxis()->SetTitle("Z(R=0) (cm)");
  hr0->GetXaxis()->SetTitleSize(0.05);
  hr0->GetYaxis()->SetTitleSize(0.05);
  hr0->GetYaxis()->SetTitleOffset(1.1);

  t->Project("hz0","reco_r_offsetFit:ptot_orig","arm==1");
  t->Project("hr0","-reco_r_offsetFit/reco_r_slopeFit:ptot_orig","arm==1");
 
  hz0->FitSlicesY(0,0,100,5); // Fit x bins 0 to 100 only, and only if containing min counts
  hr0->FitSlicesY(0,0,100,5); // Fit x bins 0 to 100 only, and only if containing min counts

  std::string cname = "resCanvas";
  TCanvas* c = gROOT->FindObject(cname.c_str());
  if ( ! c ) { 
    std::cout << "Creating " << cname << std::endl;
    c = new TCanvas(cname.c_str(),cname.c_str(),1000,600);
  }
  c->Clear();
  c->Divide(2,1);

  gStyle->SetOptFit(1111);  // prob, chi2/ndf, errors, par names
  gStyle->SetOptStat(1110); // RMS, mean, NumEntries, no hist name

  TLatex txt;
  char chtxt[1024];
  txt.SetNDC(true);

  c->cd(1);
  gPad->SetLogy(false);
  gPad->SetBottomMargin(0.11);
  gPad->SetLeftMargin(0.11);
  gPad->SetRightMargin(0.12);
  hz0->Draw("zcol");

  TH1D *hz0_1 = (TH1D*)gDirectory->Get("hz0_1"); // par 1 (mean)
  TH1D *hz0_2 = (TH1D*)gROOT->FindObject("hz0_2"); // par 2 (sigma)
  TH1D *hr0_1 = (TH1D*)gDirectory->Get("hr0_1"); // par 1 (mean)
  TH1D *hr0_2 = (TH1D*)gROOT->FindObject("hr0_2"); // par 2 (sigma)

  hz0_2->SetMarkerStyle(20);

  TPad* inset_1 = new TPad("sigmaZ0Pad","sigmaZ0Pad",0.38,0.6,0.76,0.89);
  inset_1->SetRightMargin(0.02);
  inset_1->SetTopMargin(0.04);
  inset_1->SetLeftMargin(0.15);
  inset_1->Draw();
  inset_1->cd();
  inset_1->SetGridy();
  hz0_2->SetTitle("");
  hz0_2->GetXaxis()->SetRangeUser(0.0,10.0);
  hz0_2->GetYaxis()->SetRangeUser(-0.01,0.2);
  hz0_2->GetXaxis()->SetLabelSize(0.075);
  hz0_2->GetXaxis()->SetTitle();
  hz0_2->GetYaxis()->SetLabelSize(0.075);
  hz0_2->SetStats(false);
  hz0_2->Draw();
  sprintf(chtxt,"#sigma");
  txt.SetTextSize(0.2);
  txt.DrawLatex(0.75,0.8,chtxt);


  c->cd(1);
  TPad* inset_2 = new TPad("meanZ0Pad","meanZ0Pad",0.38,0.15,0.76,0.45);
  inset_2->SetRightMargin(0.02);
  inset_2->SetTopMargin(0.04);
  inset_2->SetLeftMargin(0.15);
  inset_2->Draw();
  inset_2->cd();
  inset_2->SetGridy();
  hz0_1->SetTitle("");
  hz0_1->GetXaxis()->SetRangeUser(0.0,10.0);
  hz0_1->GetYaxis()->SetRangeUser(-0.03,0.03);
  hz0_1->GetXaxis()->SetLabelSize(0.075);
  hz0_1->GetXaxis()->SetTitle();
  hz0_1->GetYaxis()->SetLabelSize(0.075);
  hz0_1->SetStats(false);
  hz0_1->SetMarkerStyle(20);
  hz0_1->Draw();
  sprintf(chtxt,"Mean");
  txt.SetTextSize(0.15);
  txt.DrawLatex(0.65,0.8,chtxt);

  c->cd(2);
  gPad->SetLogy(false);
  gPad->SetBottomMargin(0.11);
  gPad->SetLeftMargin(0.11);
  gPad->SetRightMargin(0.12);
  hr0->Draw("zcol");

  TPad* inset_3 = new TPad("sigmaR0Pad","sigmaR0Pad",0.38,0.6,0.76,0.89);
  inset_3->SetRightMargin(0.02);
  inset_3->SetTopMargin(0.04);
  inset_3->SetLeftMargin(0.15);
  inset_3->Draw();
  inset_3->cd();
  inset_3->SetGridy();
  hr0_2->SetTitle("");
  hr0_2->GetXaxis()->SetRangeUser(0.0,10.0);
  hr0_2->GetYaxis()->SetRangeUser(-0.01,0.3);
  hr0_2->GetXaxis()->SetLabelSize(0.075);
  hr0_2->GetXaxis()->SetTitle();
  hr0_2->GetYaxis()->SetLabelSize(0.075);
  hr0_2->SetStats(false);
  hr0_2->SetMarkerStyle(20);
  hr0_2->Draw();
  sprintf(chtxt,"#sigma");
  txt.SetTextSize(0.2);
  txt.DrawLatex(0.75,0.8,chtxt);

  c->cd(2);
  TPad* inset_4 = new TPad("meanR0Pad","meanR0Pad",0.38,0.15,0.76,0.45);
  inset_4->SetRightMargin(0.02);
  inset_4->SetTopMargin(0.04);
  inset_4->SetLeftMargin(0.15);
  inset_4->Draw();
  inset_4->cd();
  inset_4->SetGridy();
  hr0_1->SetTitle("");
  hr0_1->GetXaxis()->SetRangeUser(0.0,10.0);
  hr0_1->GetYaxis()->SetRangeUser(-0.04,0.04);
  hr0_1->GetXaxis()->SetLabelSize(0.075);
  hr0_1->GetXaxis()->SetTitle();
  hr0_1->GetYaxis()->SetLabelSize(0.075);
  hr0_1->SetStats(false);
  hr0_1->SetMarkerStyle(20);
  hr0_1->Draw();
  sprintf(chtxt,"Mean");
  txt.SetTextSize(0.15);
  txt.DrawLatex(0.65,0.8,chtxt);

  return;
}

void
draw_kfresult()
{
  int nbins = 100;
  double maxz =  0.5;

  double maxr =  1.0;

  std::string tname = "mc_trk_eval";
  TTree* t = gROOT->FindObject(tname.c_str());
  if ( ! t ) { 
    std::cout << "Failed to find "  << tname << " tree" << std::endl;
    return;
  }
 
  delete gROOT->FindObject("hz0_kf");
  delete gROOT->FindObject("hz01d_kf");
  delete gROOT->FindObject("hr0_kf");

  TH1F* hz01d = new TH1F("hz01d_kf","R(Z=0)",nbins,-maxz,maxz);

  TH2F* hz0 = new TH2F("hz0_kf","R(Z=0) vs. MC p",100,0.0,10.0,nbins,-maxz,maxz);
  hz0->GetXaxis()->SetTitle("MC p (GeV/c)");
  hz0->GetYaxis()->SetTitle("Impact Parameter (Z=0) (cm)");
  hz0->GetXaxis()->SetTitleSize(0.05);
  hz0->GetYaxis()->SetTitleSize(0.05);
  hz0->GetXaxis()->SetTitleOffset(1.0);
  hz0->GetYaxis()->SetTitleOffset(1.1);

  TH2F* hz0_B = new TH2F("hz0_kf_B","R(0) vs. p_{MC} from B daughters",100,0.0,10.0,nbins,-maxz,maxz);
  hz0_B->GetXaxis()->SetTitle("MC p (GeV/c)");
  hz0_B->GetYaxis()->SetTitle("Impact Parameter (Z=0) (cm)");
  hz0_B->GetXaxis()->SetTitleSize(0.05);
  hz0_B->GetYaxis()->SetTitleSize(0.05);
  hz0_B->GetXaxis()->SetTitleOffset(1.0);
  hz0_B->GetYaxis()->SetTitleOffset(1.1);

  TH2F* hphi0 = new TH2F("hphi0_kf","PHI(Z=0) vs. MC p",100,0.0,10.0,nbins,-maxz,maxz);
  hphi0->GetXaxis()->SetTitle("MC p (GeV/c)");
  hphi0->GetYaxis()->SetTitle("(#phi) Impact Parameter (Z=0) (cm)");
  hphi0->GetXaxis()->SetTitleSize(0.05);
  hphi0->GetYaxis()->SetTitleSize(0.05);
  hphi0->GetXaxis()->SetTitleOffset(1.0);
  hphi0->GetYaxis()->SetTitleOffset(1.1);

  TH2F* hphi0mc = new TH2F("hphi0mc_kf","PHI(Z=0) vs. MC p",100,0.0,10.0,nbins,-maxz,maxz);
  hphi0mc->GetXaxis()->SetTitle("MC p (GeV/c)");
  hphi0mc->GetYaxis()->SetTitle("(#phi) Impact Parameter (Z=0) from p_{MC} (cm)");
  hphi0mc->GetXaxis()->SetTitleSize(0.05);
  hphi0mc->GetYaxis()->SetTitleSize(0.05);
  hphi0mc->GetXaxis()->SetTitleOffset(1.0);
  hphi0mc->GetYaxis()->SetTitleOffset(1.1);

  TH2F* hr0 = new TH2F("hr0_kf","Z(R=0) vs. MC p",100,0.0,10.0,nbins,-maxr,maxr);
  hr0->GetXaxis()->SetTitle("MC p (GeV/c)");
  hr0->GetYaxis()->SetTitle("Z(R=0) (cm)");
  hr0->GetXaxis()->SetTitleSize(0.05);
  hr0->GetYaxis()->SetTitleSize(0.05);
  hr0->GetYaxis()->SetTitleOffset(1.1);

  // Apply the cuts:
  // - number of mchits for the mc track is > 3
  // - number of coords found > 2
  // - number of stations hit > 3
  // - KF reconstruction succeeded
  //
  t->Project("hz0_kf",
	     "x0reco*(px0reco/sqrt(px0reco*px0reco+py0reco*py0reco))+y0reco*(py0reco/sqrt(px0reco*px0reco+py0reco*py0reco)):sqrt(px0mc*px0mc+py0mc*py0mc+pz0mc*pz0mc)",
	     "nhitst>3&&nhitsfoundt>2&&nstationst>3&&recosuc==1");
  //t->Project("hr0","-reco_r_offsetFit/reco_r_slopeFit:ptot_orig","arm==1");

  t->Project("hz0_kf_B",
	     "x0reco*(px0reco/sqrt(px0reco*px0reco+py0reco*py0reco))+y0reco*(py0reco/sqrt(px0reco*px0reco+py0reco*py0reco)):sqrt(px0mc*px0mc+py0mc*py0mc+pz0mc*pz0mc)",
	     "nhitst>3&&nhitsfoundt>2&&nstationst>3&&recosuc==1&&(abs(parent_pdgid)==511||abs(parent_pdgid==521))");
 
  t->Project("hphi0_kf",
	     "y0reco*(px0reco/sqrt(px0reco*px0reco+py0reco*py0reco))-x0reco*(py0reco/sqrt(px0reco*px0reco+py0reco*py0reco)):sqrt(px0mc*px0mc+py0mc*py0mc+pz0mc*pz0mc)",
	     "nhitst>3&&nhitsfoundt>2&&nstationst>3&&recosuc==1");

  t->Project("hphi0mc_kf",
	     "y0reco*(px0mc/sqrt(px0mc*px0mc+py0mc*py0mc))-x0reco*(py0mc/sqrt(px0mc*px0mc+py0mc*py0mc)):sqrt(px0mc*px0mc+py0mc*py0mc+pz0mc*pz0mc)",
	     "nhitst>3&&nhitsfoundt>2&&nstationst>3&&recosuc==1&&abs(z0reco)<0.01");

  hz0->FitSlicesY(0,0,100,5); // Fit x bins 0 to 100 only, and only if containing min counts
  hr0->FitSlicesY(0,0,100,5); // Fit x bins 0 to 100 only, and only if containing min counts

  hz0_B->FitSlicesY(0,0,100,5); // Fit x bins 0 to 100 only, and only if containing min counts

  std::string cname = "kfCanvasR";
  TCanvas* c = gROOT->FindObject(cname.c_str());
  if ( ! c ) { 
    std::cout << "Creating " << cname << std::endl;
    c = new TCanvas(cname.c_str(),cname.c_str(),900,540);
  }
  c->Clear();
  c->Divide(2,1);

  std::string cname = "kfCanvasPhi";
  TCanvas* c1 = gROOT->FindObject(cname.c_str());
  if ( ! c1 ) { 
    std::cout << "Creating " << cname << std::endl;
    c1 = new TCanvas(cname.c_str(),cname.c_str(),900,540);
  }
  c1->Clear();
  c1->Divide(2,1);

  gStyle->SetPalette(1);
  gStyle->SetOptFit(1111);  // prob, chi2/ndf, errors, par names
  gStyle->SetOptStat(1110); // RMS, mean, NumEntries, no hist name

  TLatex txt;
  char chtxt[1024];
  txt.SetNDC(true);

  c->cd(1);
  gPad->SetLogy(false);
  gPad->SetBottomMargin(0.11);
  gPad->SetLeftMargin(0.11);
  gPad->SetRightMargin(0.12);
  hz0->Draw("zcol");

  TH1D *hz0_1 = (TH1D*)gDirectory->Get("hz0_kf_1"); // par 1 (mean)
  TH1D *hz0_2 = (TH1D*)gROOT->FindObject("hz0_kf_2"); // par 2 (sigma)
  TH1D *hz0_B_1 = (TH1D*)gDirectory->Get("hz0_kf_B_1"); // par 1 (mean)
  TH1D *hz0_B_2 = (TH1D*)gROOT->FindObject("hz0_kf_B_2"); // par 2 (sigma)
  TH1D *hr0_1 = (TH1D*)gDirectory->Get("hr0_kf_1"); // par 1 (mean)
  TH1D *hr0_2 = (TH1D*)gROOT->FindObject("hr0_kf_2"); // par 2 (sigma)

  hz0_2->SetMarkerStyle(20);

  TPad* inset_1 = new TPad("sigmaZ0Pad_kf","sigmaZ0Pad_kf",0.38,0.6,0.76,0.89);
  inset_1->SetRightMargin(0.02);
  inset_1->SetTopMargin(0.04);
  inset_1->SetLeftMargin(0.15);
  inset_1->Draw();
  inset_1->cd();
  inset_1->SetGridy();
  hz0_2->SetTitle("");
  hz0_2->GetXaxis()->SetRangeUser(0.0,10.0);
  hz0_2->GetYaxis()->SetRangeUser(-0.01,0.2);
  hz0_2->GetXaxis()->SetLabelSize(0.075);
  hz0_2->GetXaxis()->SetTitle();
  hz0_2->GetYaxis()->SetLabelSize(0.075);
  hz0_2->SetStats(false);
  hz0_2->Draw();
  sprintf(chtxt,"#sigma");
  txt.SetTextSize(0.2);
  txt.DrawLatex(0.75,0.8,chtxt);


  c->cd(1);
  TPad* inset_2 = new TPad("meanZ0Pad_kf","meanZ0Pad_kf",0.38,0.15,0.76,0.45);
  inset_2->SetRightMargin(0.02);
  inset_2->SetTopMargin(0.04);
  inset_2->SetLeftMargin(0.15);
  inset_2->Draw();
  inset_2->cd();
  inset_2->SetGridy();
  hz0_1->SetTitle("");
  hz0_1->GetXaxis()->SetRangeUser(0.0,10.0);
  hz0_1->GetYaxis()->SetRangeUser(-0.03,0.03);
  hz0_1->GetXaxis()->SetLabelSize(0.075);
  hz0_1->GetXaxis()->SetTitle();
  hz0_1->GetYaxis()->SetLabelSize(0.075);
  hz0_1->SetStats(false);
  hz0_1->SetMarkerStyle(20);
  hz0_1->Draw();
  sprintf(chtxt,"Mean");
  txt.SetTextSize(0.15);
  txt.DrawLatex(0.65,0.8,chtxt);

  c->cd(2);
  gPad->SetLogy(false);
  gPad->SetBottomMargin(0.11);
  gPad->SetLeftMargin(0.11);
  gPad->SetRightMargin(0.12);
  //hr0->Draw("zcol");
  hz0_B->Draw("zcol");

  TPad* inset_3 = new TPad("sigmaR0Pad_kf","sigmaR0Pad_kf",0.38,0.6,0.76,0.89);
  inset_3->SetRightMargin(0.02);
  inset_3->SetTopMargin(0.04);
  inset_3->SetLeftMargin(0.15);
  inset_3->Draw();
  inset_3->cd();
  inset_3->SetGridy();
  hz0_B_2->SetTitle("");
  hz0_B_2->GetXaxis()->SetRangeUser(0.0,10.0);
  hz0_B_2->GetYaxis()->SetRangeUser(-0.01,0.3);
  hz0_B_2->GetXaxis()->SetLabelSize(0.075);
  hz0_B_2->GetXaxis()->SetTitle();
  hz0_B_2->GetYaxis()->SetLabelSize(0.075);
  hz0_B_2->SetStats(false);
  hz0_B_2->SetMarkerStyle(20);
  hz0_B_2->Draw();
  sprintf(chtxt,"#sigma");
  txt.SetTextSize(0.2);
  txt.DrawLatex(0.75,0.8,chtxt);

  c->cd(2);
  TPad* inset_4 = new TPad("meanR0Pad_kf","meanR0Pad_kf",0.38,0.15,0.76,0.45);
  inset_4->SetRightMargin(0.02);
  inset_4->SetTopMargin(0.04);
  inset_4->SetLeftMargin(0.15);
  inset_4->Draw();
  inset_4->cd();
  inset_4->SetGridy();
  hz0_B_1->SetTitle("");
  hz0_B_1->GetXaxis()->SetRangeUser(0.0,10.0);
  hz0_B_1->GetYaxis()->SetRangeUser(-0.04,0.04);
  hz0_B_1->GetXaxis()->SetLabelSize(0.075);
  hz0_B_1->GetXaxis()->SetTitle();
  hz0_B_1->GetYaxis()->SetLabelSize(0.075);
  hz0_B_1->SetStats(false);
  hz0_B_1->SetMarkerStyle(20);
  hz0_B_1->Draw();
  sprintf(chtxt,"Mean");
  txt.SetTextSize(0.15);
  txt.DrawLatex(0.65,0.8,chtxt);

  c1->cd(1);
  gPad->SetLogy(false);
  gPad->SetBottomMargin(0.11);
  gPad->SetLeftMargin(0.11);
  gPad->SetRightMargin(0.12);
  hphi0->Draw("zcol");

  c1->cd(2);
  gPad->SetLogy(false);
  gPad->SetBottomMargin(0.11);
  gPad->SetLeftMargin(0.11);
  gPad->SetRightMargin(0.12);
  hphi0mc->Draw("zcol");

  return;
}

void
draw_eff()
{
  int nbins = 25;
  double pmin = 0.0;
  double pmax = 15.0;
  std::string tname = "match_eval";
  std::string cname = "matchEffCanvas";

  TTree* t = gROOT->FindObject(tname.c_str());
  if ( ! t ) { 
    std::cout << "draw_eff: Failed to find " << tname << " tree" << std::endl;
    return;
  }

  TCanvas* c_tmp = gROOT->FindObject(cname.c_str());
  if ( ! c_tmp ) { 
    std::cout << "Creating " << cname.c_str() << std::endl;
    //c_tmp = new TCanvas(cname.c_str(),cname.c_str(),470,360);
    c_tmp = new TCanvas(cname.c_str(),cname.c_str(),533,400);
  } 
  c_tmp->cd();
  gPad->SetBottomMargin(0.12);
  gPad->SetRightMargin(0.1);

  TH1F* hptotal = new TH1F("hptotal","ptot spectrum",nbins,0.0,pmax);
  TH1F* eff = new TH1F("eff","eff of finding >=1 match",nbins,0.0,pmax);
  TH1F* eff1 = new TH1F("eff1","rate of finding =1 match",nbins,0.0,pmax);
  TH1F* eff2 = new TH1F("eff2","rate of finding >1 match",nbins,0.0,pmax);
  //TH1F* matchEq0 = new TH1F("matchEq0","fraction of tracks with nmatch==0",nbins,0.0,pmax);
  TH1F* matchGT0 = new TH1F("matchGT0","fraction of tracks with nmatch>0",nbins,0.0,pmax);
  TH1F* matchEq1 = new TH1F("matchEq1","fraction of tracks with nmatch==1",nbins,0.0,pmax);
  TH1F* matchGT1 = new TH1F("matchGT1","fraction of tracks with nmatch>1",nbins,0.0,pmax);

  t->Project("hptotal","ptot_orig","arm==1&&nstaHit==4");
  t->Project("matchGT0","ptot_orig","arm==1&&nstaHit==4&&nreco>0");
  t->Project("matchEq1","ptot_orig","arm==1&&nstaHit==4&&nreco==1");
  t->Project("matchGT1","ptot_orig","arm==1&&nstaHit==4&&nreco>1");

//   t->Project("hptotal","ptot_orig","nstaHit==4");
//   t->Project("matchGT0","ptot_orig","nstaHit==4&&nreco>0");
//   t->Project("matchEq1","ptot_orig","nstaHit==4&&nreco==1");
//   t->Project("matchGT1","ptot_orig","nstaHit==4&&nreco>1");

  eff->Divide(matchGT0,hptotal,1.0,1.0);
  eff1->Divide(matchEq1,hptotal,1.0,1.0);
  eff2->Divide(matchGT1,hptotal,1.0,1.0);

  TLegend* leg = new TLegend(0.19,0.30,0.55,0.65,"","NDC");
  leg->SetName("Legend");

  TH1F* hframe = new TH1F("hframe","hframe",eff->GetXaxis()->GetNbins(),0,pmax);
  hframe->SetMaximum(1.05);
  //hframe->SetMinimum(0.25);
  hframe->SetTitle("North Arm, 4 Stations Required");
  hframe->SetStats(false);
  hframe->GetXaxis()->SetLabelSize(0.05);
  hframe->GetXaxis()->SetTitleSize(0.06);
  hframe->GetXaxis()->SetTitleOffset(0.82);
  hframe->GetXaxis()->SetTitle("MC p (GeV/c)");
  hframe->GetYaxis()->SetLabelSize(0.05);
  hframe->GetYaxis()->SetTitle("Rate of matching MC to Reco (fraction)");
  hframe->GetYaxis()->SetTitleSize(0.05);
  hframe->Draw();

  eff->SetMarkerStyle(20);
  eff->Draw("pl,same");
  leg->AddEntry(eff,"Matches>0","pl");

  eff1->SetMarkerStyle(20);
  eff1->SetMarkerColor(2);
  eff1->SetLineColor(2);
  eff1->Draw("pl,same");
  leg->AddEntry(eff1,"Matches=1","pl");

  gPad->Update();

  float rightmax = 1.1*eff2->GetMaximum();
  if ( rightmax > 0.0 ) {
    float scale = gPad->GetUymax()/rightmax;
    eff2->Scale(scale);
  }

  eff2->SetMarkerStyle(20);
  eff2->SetMarkerColor(4);
  eff2->SetLineColor(4);
  eff2->Draw("pl,same");
  leg->AddEntry(eff2,"Matches>1","pl");

  if ( rightmax == 0.0 ) rightmax = 1.0;
  TGaxis *axis = new TGaxis(gPad->GetUxmax(),gPad->GetUymin(),
			    gPad->GetUxmax(), gPad->GetUymax(),0,rightmax,505,"+L");
  axis->SetLabelSize(0.05);
  axis->SetLineColor(4);
  axis->SetLabelColor(4);
  axis->Draw();

  leg->Draw();

//   TLatex txt;
//   char chtxt[1024];
//   txt.SetNDC(true);
//   sprintf(chtxt,"4 Sta Required");
//   txt.DrawLatex(0.25,0.65,chtxt);

//   TPad* inset = new TPad("matchGt1","matchGt1",0.5,0.15,0.93,0.6);
//   inset->SetLeftMargin(0.15);
//   inset->SetRightMargin(0.01);
//   inset->SetTopMargin(0.02);
//   inset->Draw();
//   inset->cd();
//   TH1* h_inset = hframe->Clone();
//   double max = eff2->GetMaximum();
//   h_inset->SetMaximum(max+0.1*max);
//   h_inset->SetMinimum(0.0);
//   h_inset->SetTitle("");
//   h_inset->GetXaxis()->SetTitle("");
//   h_inset->GetXaxis()->SetLabelSize(0.07);
//   h_inset->GetYaxis()->SetLabelSize(0.07);
//   h_inset->GetYaxis()->SetTitle("");
//   h_inset->SetStats(false);
//   h_inset->Draw();
  
//   eff2->Draw("pl,same");

  return;
}

#ifndef M_PI
# define M_PI		3.14159265358979323846	/* pi */
#endif

// Normalizes angle to (pi,-pi]
double
angle_normalize(const double x)
{
  double ret;
  if ( x > -3 * M_PI && x <= 3 * M_PI )
    {
      ret = x;
      if ( ret > M_PI )        ret -= 2 * M_PI;
      else if ( ret <= -M_PI ) ret += 2 * M_PI;
    }
  else 
    {
      const double x_pi_mod = fmod(x + M_PI, 2 * M_PI);      
      ret = x_pi_mod > 0 ? x_pi_mod - M_PI : x_pi_mod + M_PI;
    }

  return ret;
}

void
draw_matching()
{
  std::string mctrk_tree_name = "mc_trk_eval";
  std::string trk_tree_name = "trk_eval";
  std::string matchingCanvas = "matchCanvas";

  TTree* mctrks = gROOT->FindObject(mctrk_tree_name.c_str());
  if ( ! mctrks ) { 
    std::cout << "Failed to find tree " << mctrk_tree_name << std::endl;
    return;
  }

  TTree* trks = gROOT->FindObject(trk_tree_name.c_str());
  if ( ! trks ) { 
    std::cout << "Failed to find tree " << trk_tree_name << std::endl;
    return;
  }

  double x1mc, y1mc, z1mc;
  double x4mc, y4mc, z4mc;
  int nhits;
  int hitIds[100];
  mctrks->SetBranchAddress("x1mc",&x1mc);
  mctrks->SetBranchAddress("y1mc",&x1mc);
  mctrks->SetBranchAddress("z1mc",&x1mc);
  mctrks->SetBranchAddress("x4mc",&x1mc);
  mctrks->SetBranchAddress("y4mc",&x1mc);
  mctrks->SetBranchAddress("z4mc",&x1mc);
  mctrks->SetBranchAddress("nhitst",&nhits);
  mctrks->SetBranchAddress("hitstIds",hitIds);

  double dslope = 0.001;
  double dr = 0.1;

  TH2F* hMatch = new TH2F("hMatch","Matching Dist",100,-dslope,dslope,100,-dr,dr);

  TCanvas* cMatch = gROOT->FindObject(matchingCanvas.c_str());
  if ( ! cMatch ) { 
    std::cout << "Creating canvas " << matchingCanvas << std::endl;
    cMatch = new TCanvas(matchingCanvas.c_str(),matchingCanvas.c_str(),470,360);
  } 
  cMatch->Clear();
  cMatch->cd();
  gPad->SetBottomMargin(0.12);
  gPad->SetRightMargin(0.05);

}

void
draw_compareKF()
{
  draw_resolution();
  draw_kfresult();

  TH1D *hz0_kf_2 = (TH1D*)gROOT->FindObject("hz0_kf_2"); // par 2 (sigma)
  TH1D *hz0_2 = (TH1D*)gROOT->FindObject("hz0_2"); // par 2 (sigma)

  std::string cname = "compareKFCanvas";
  TCanvas* c = gROOT->FindObject(cname.c_str());
  if ( ! c ) { 
    std::cout << "Creating " << cname << std::endl;
    c = new TCanvas(cname.c_str(),cname.c_str(),533,400);
  }
  c->Clear();
  gPad->SetRightMargin(0.04);
  gPad->SetTopMargin(0.04);
  gPad->SetLeftMargin(0.15);
  gPad->SetGridy();
  gPad->SetGridx();
  gPad->SetBottomMargin(0.15);
  //gPad->SetLeftMargin(0.12);
  //gPad->SetRightMargin(0.05);

  c->cd();

  hz0_2->GetXaxis()->SetTitle("MC p (GeV/c)");
  hz0_2->GetXaxis()->SetTitleSize(0.07);
  hz0_2->GetYaxis()->SetTitle("#sigma (cm)");
  hz0_2->GetYaxis()->SetTitleSize(0.075);
  hz0_2->Draw();

  TH1* tmp = hz0_kf_2->Clone("hz0_kf_2_copy");
  tmp->SetMarkerColor(2);
  tmp->SetLineColor(2);
  tmp->Draw("same");

  TLegend* legend0 = new TLegend(0.6,0.70,0.99,0.88,"","NDC");
  legend0->AddEntry(hz0_2,"R-Z linear fit","PL");
  legend0->AddEntry(tmp,"KF fit & extrap","PL");
  legend0->Draw();
}

void
draw_pres()
{
  int nbins = 100;
  double maxz =  5.0;
  double maxr =  1.0;

  std::string tname = "mc_trk_eval";
  TTree* t = gROOT->FindObject(tname.c_str());
  if ( ! t ) { 
    std::cout << "Failed to find "  << tname << " tree" << std::endl;
    return;
  }

  TH2F* hpres = new TH2F("hpres","Momentum Resolution",100,0.0,10.0,nbins,-maxz,maxz);
  hpres->GetXaxis()->SetTitle("p_{MC} (GeV/c)");
  hpres->GetYaxis()->SetTitle("(p - p_{MC}) / p_{MC}");
//   hpres->GetXaxis()->SetTitleSize(0.05);
//   hpres->GetYaxis()->SetTitleSize(0.05);
//   hpres->GetXaxis()->SetTitleOffset(1.0);
//   hpres->GetYaxis()->SetTitleOffset(1.1);

  t->Project("hpres",
	     "(sqrt(px0reco*px0reco+py0reco*py0reco+pz0reco*pz0reco)-sqrt(px0mc*px0mc+py0mc*py0mc+pz0mc*pz0mc))/sqrt(px0mc*px0mc+py0mc*py0mc+pz0mc*pz0mc):sqrt(px0mc*px0mc+py0mc*py0mc+pz0mc*pz0mc)",
	     "nhitst>3&&nhitsfoundt>2&&nstationst>3&&recosuc==1");

  std::string cname = "presCanvas";
  TCanvas* c = gROOT->FindObject(cname.c_str());
  if ( ! c ) { 
    std::cout << "Creating " << cname << std::endl;
    c = new TCanvas(cname.c_str(),cname.c_str(),470,360);
  }
  c->Clear();
  c->SetLogz();
  gPad->SetBottomMargin(0.12);
  gPad->SetRightMargin(0.1);

  gStyle->SetPalette(1);

  hpres->GetXaxis()->SetLabelSize(0.05);
  hpres->GetXaxis()->SetTitleSize(0.06);
  hpres->GetXaxis()->SetTitleOffset(0.82);
  //hpres->GetXaxis()->SetTitle("MC p (GeV/c)");
  hpres->GetYaxis()->SetLabelSize(0.05);
  hpres->GetYaxis()->SetLabelOffset(0.01);
  //hpres->GetYaxis()->SetTitle("Rate of matching MC to Reco (fraction)");
  hpres->GetYaxis()->SetTitleSize(0.05);
  hpres->GetYaxis()->SetTitleOffset(0.82);

  hpres->Draw("zcol");
  
  gPad->Update();
  TPaveStats* s = (TPaveStats*) gPad->GetPrimitive("stats");
  s->SetX1NDC(0.7);
  s->SetX2NDC(0.91);
  gPad->Modified();
}

void
draw_eta_z_accept()
{
  std::string tname = "mc_trk_eval";
  TTree* t = gROOT->FindObject(tname.c_str());
  if ( ! t ) { 
    std::cout << "Failed to find "  << tname << " tree" << std::endl;
    return;
  }

  std::string cname = "etazCanvas";
  TCanvas* c = gROOT->FindObject(cname.c_str());
  if ( ! c ) { 
    std::cout << "Creating " << cname << std::endl;
    c = new TCanvas(cname.c_str(),cname.c_str(),470,360);
  }
  c->Clear();
  gPad->SetBottomMargin(0.12);
  gPad->SetLeftMargin(0.12);
  gPad->SetRightMargin(0.1);
  gPad->SetGridx(true);
  gPad->SetGridy(true);

  std::string cname = "thetazCanvas";
  TCanvas* c1 = gROOT->FindObject(cname.c_str());
  if ( ! c1 ) { 
    std::cout << "Creating " << cname << std::endl;
    c1 = new TCanvas(cname.c_str(),cname.c_str(),470,360);
  }
  c1->Clear();
  gPad->SetBottomMargin(0.12);
  gPad->SetLeftMargin(0.12);
  gPad->SetRightMargin(0.1);
  gPad->SetGridx(true);
  gPad->SetGridy(true);

  TH2F* h = new TH2F("hetaz","#eta-z acceptance, North Arm",
		     60,-20.0,20.0,50,0.0,3.5);
  h->GetXaxis()->SetTitle("MC Track Z0 (cm)");
  h->GetYaxis()->SetTitle("#eta");
  h->GetXaxis()->SetLabelSize(0.05);
  h->GetXaxis()->SetTitleSize(0.06);
  h->GetXaxis()->SetTitleOffset(0.82);
  h->GetYaxis()->SetLabelSize(0.05);
  h->GetYaxis()->SetLabelOffset(0.01);
  h->GetYaxis()->SetTitleSize(0.07);
  h->GetYaxis()->SetTitleOffset(0.82);
  h->SetStats(0);
  h->SetFillColor(1);

  TH2F* hh = new TH2F("hthetaz","#theta-z acceptance, North Arm",
		      h->GetXaxis()->GetNbins(),h->GetXaxis()->GetXmin(),h->GetXaxis()->GetXmax(),
		      30,0.0,60);
  hh->GetXaxis()->SetTitle("MC Track Z0 (cm)");
  hh->GetYaxis()->SetTitle("#theta (degrees)");
  hh->GetXaxis()->SetLabelSize(0.05);
  hh->GetXaxis()->SetTitleSize(0.06);
  hh->GetXaxis()->SetTitleOffset(0.82);
  hh->GetYaxis()->SetLabelSize(0.05);
  hh->GetYaxis()->SetLabelOffset(0.01);
  hh->GetYaxis()->SetTitleSize(0.07);
  hh->GetYaxis()->SetTitleOffset(0.82);
  hh->SetStats(0);
  hh->SetFillColor(1);

  TH1* h1 = h->Clone("hetaz_vtx1"); 
  h1->SetFillColor(2);

  TH1* h2 = h->Clone("hetaz_vtx2"); 
  h2->SetFillColor(3);

  TH1* hh1 = hh->Clone("hthetaz_vtx1"); 
  hh1->SetFillColor(2);

  TH1* hh2 = hh->Clone("hthetaz_vtx2"); 
  hh2->SetFillColor(3);

  // Plot inclusive tracks
  TCut a = "abs(atan2(py0mc,px0mc))<TMath::Pi()/4";
  int n = 0;
  n = t->Project("hetaz","-log(tan(atan2(sqrt(px0mc**2+py0mc**2),pz0mc)/2.0)):z0mc",
		 "arm==1&&primId==0&&pz0mc>0&&nstationst==4");

  // Plot 4-station hits with 1 barrel hits
  n = t->Project("hetaz_vtx1","-log(tan(atan2(sqrt(px0mc**2+py0mc**2),pz0mc)/2.0)):z0mc",
		 "arm==1&&primId==0&&pz0mc>0&&nstationst==4&&nhitsb==1");


  // Plot 4-station hits with 2 barrel hits
  n = t->Project("hetaz_vtx2","-log(tan(atan2(sqrt(px0mc**2+py0mc**2),pz0mc)/2.0)):z0mc",
		 "arm==1&&primId==0&&pz0mc>0&&nstationst==4&&nhitsb==2");

  // do the same for the theta plot
  //
  n = t->Project("hthetaz","atan2(sqrt(px0mc**2+py0mc**2),pz0mc)*180.0/TMath::Pi():z0mc",
		 "arm==1&&primId==0&&pz0mc>0&&nstationst==4");

  n = t->Project("hthetaz_vtx1","atan2(sqrt(px0mc**2+py0mc**2),pz0mc)*180.0/TMath::Pi():z0mc",
		 "arm==1&&primId==0&&pz0mc>0&&nstationst==4&&nhitsb==1");

  n = t->Project("hthetaz_vtx2","atan2(sqrt(px0mc**2+py0mc**2),pz0mc)*180.0/TMath::Pi():z0mc",
		 "arm==1&&primId==0&&pz0mc>0&&nstationst==4&&nhitsb==2");

  // Draw the same legend on each canvas
  TLegend* l = new TLegend(0.6,0.6,0.89,0.89,"4 FVTX Sta Hit","NDC");
  l->AddEntry(h,"No VTX req","F");
  l->AddEntry(h1,"VTX hits = 1","F");
  l->AddEntry(h2,"VTX hits = 2","F");

  c->cd();
  h->Draw("box");
  //h1->Scale(4); // To make the box plot easier to see
  h1->Draw("box same");
  //h2->Scale(4); // To make the box plot easier to see
  h2->Draw("box same");
  l->Draw();

  c1->cd();
  hh->Draw("box");
  hh1->Draw("box same");
  hh2->Draw("box same");
  l->Draw();

  return;

}

void
draw_chi2slndf()
{
  std::string tname = "mc_trk_eval";
  TTree* t = gROOT->FindObject(tname.c_str());
  if ( ! t ) { 
    std::cout << "Failed to find "  << tname << " tree" << std::endl;
    return;
  }

  std::string cname = "chi2slndfCanvas";
  TCanvas* c = gROOT->FindObject(cname.c_str());
  if ( ! c ) { 
    std::cout << "Creating " << cname << std::endl;
    c = new TCanvas(cname.c_str(),cname.c_str(),800,600);
  }
  c->Clear();
  c->Divide(4,2);

  TCut fourStaTrue = "nstationst>=4";
  TCut chi2slCut = "chi2sl/ndfsl<4";

  for(int i=0; i<8; i++)
    {
      c->cd(i+1);
      gPad->SetBottomMargin(0.12);
      gPad->SetLeftMargin(0.15);
      gPad->SetRightMargin(0.1);
      gPad->SetGridx(true);
      gPad->SetGridy(true);

      
      char cut[100];
      sprintf(cut,"ndfsl==%d",2*i+4);
      char name[100];
      //char title[100];
      sprintf(name,"hchi2ndf_%d",2*i+4);

      TH1F* h = new TH1F(name,"",90,0.0,3.0);
      h->SetStats(0);
      h->GetXaxis()->SetTitle("#chi^{2}/n.d.f.");
      h->GetXaxis()->SetLabelSize(0.05);
      h->GetXaxis()->SetTitleSize(0.1);
      h->GetXaxis()->SetTitleOffset(0.5);

      t->Project(name,"chi2sl/ndfsl", fourStaTrue && cut);
      
      h->Draw();
      
      TLatex txt;
      txt.SetNDC(true);
      txt.SetTextAlign(22);
      txt.SetTextSize(0.1);
      sprintf(cut,"n.d.f. = %d",2*i+4);
      txt.DrawLatex(0.5,0.85,cut);
    }

  std::string cname = "chi2KFndfCanvas";
  TCanvas* c2 = gROOT->FindObject(cname.c_str());
  if ( ! c2 ) { 
    std::cout << "Creating " << cname << std::endl;
    c2 = new TCanvas(cname.c_str(),cname.c_str(),800,600);
  }
  c2->Clear();
  c2->Divide(4,2);

  for(int i=0; i<8; i++)
    {
      c2->cd(i+1);
      gPad->SetBottomMargin(0.12);
      gPad->SetLeftMargin(0.15);
      gPad->SetRightMargin(0.1);
      gPad->SetGridx(true);
      gPad->SetGridy(true);

      
      char cut[100];
      sprintf(cut,"chisqreco/chisqrecopdf==%d",2*(i+4)-5);
      char name[100];
      //char title[100];
      sprintf(name,"hchi2ndfKF_%d",2*(i+4)-5);

      TH1F* h = new TH1F(name,"",90,0.0,3.0);
      h->SetStats(0);
      h->GetXaxis()->SetTitle("#chi^{2}/n.d.f.");
      h->GetXaxis()->SetLabelSize(0.05);
      h->GetXaxis()->SetTitleSize(0.1);
      h->GetXaxis()->SetTitleOffset(0.5);

      t->Project(name,"chisqrecopdf", fourStaTrue && cut);
      
      h->Draw();
      
      TLatex txt;
      txt.SetNDC(true);
      txt.SetTextAlign(22);
      txt.SetTextSize(0.1);
      sprintf(cut,"n.d.f. = %d",2*(i+4)-5);
      txt.DrawLatex(0.5,0.85,cut);
    }

}

void
draw_slopecmp()
{
  std::string tname = "mc_trk_eval";
  TTree* t = gROOT->FindObject(tname.c_str());
  if ( ! t ) { 
    std::cout << "Failed to find "  << tname << " tree" << std::endl;
    return;
  }

  std::string cname = "slopecmpCanvas";
  TCanvas* c = gROOT->FindObject(cname.c_str());
  if ( ! c ) { 
    std::cout << "Creating " << cname << std::endl;
    c = new TCanvas(cname.c_str(),cname.c_str(),800,600);
  }
  c->Clear();
  c->Divide(2,2);

  h0x = new TH2F("h0x","Both Arms",100,-0.5,0.5,100,-0.5,0.5);
  h0y = new TH2F("h0y","Both Arms",100,-0.5,0.5,100,-0.5,0.5);
  h1x = new TH2F("h1x","South Arm",100,-0.5,0.5,100,-0.5,0.5);
  h1y = new TH2F("h1y","South Arm",100,-0.5,0.5,100,-0.5,0.5);

  hDiffx = new TH1F("hDiffx","MC - Reco X slopes",200,-0.02,0.02);
  hDiffy = new TH1F("hDiffy","MC - Reco Y slopes",200,-0.02,0.02);

  gStyle->SetOptStat("e");

  h0x->SetXTitle("p_{x,MC}/p_{z,MC}");
  h0x->SetYTitle("M_{x}");
  h0x->GetXaxis()->SetTitleSize(0.1);
  h0x->GetYaxis()->SetTitleSize(0.1);
  h0x->GetXaxis()->SetTitleOffset(0.7);
  h0x->GetYaxis()->SetTitleOffset(0.7);
  h0x->GetXaxis()->SetNdivisions(505);
  h0x->GetYaxis()->SetNdivisions(505);
  h0x->GetXaxis()->SetLabelSize(0.05);
  h0x->GetXaxis()->SetLabelOffset(0.02);
  h0x->GetYaxis()->SetLabelSize(0.05);

  h0y->SetXTitle("p_{y,MC}/p_{z,MC}");
  h0y->SetYTitle("M_{y}");
  h0y->GetXaxis()->SetTitleSize(0.1);
  h0y->GetYaxis()->SetTitleSize(0.1);
  h0y->GetXaxis()->SetTitleOffset(0.7);
  h0y->GetYaxis()->SetTitleOffset(0.7);
  h0y->GetXaxis()->SetNdivisions(505);
  h0y->GetYaxis()->SetNdivisions(505);
  h0y->GetXaxis()->SetLabelSize(0.05);
  h0y->GetXaxis()->SetLabelOffset(0.02);
  h0y->GetYaxis()->SetLabelSize(0.05);

  h1x->SetXTitle("p_{x,MC}/p_{z,MC}");
  h1x->SetYTitle("M_{x}");
  h1x->GetXaxis()->SetTitleSize(0.1);
  h1x->GetYaxis()->SetTitleSize(0.1);
  h1x->GetXaxis()->SetTitleOffset(0.7);
  h1x->GetYaxis()->SetTitleOffset(0.7);

  h1y->SetXTitle("p_{y,MC}/p_{z,MC}");
  h1y->SetYTitle("M_{y}");
  h1y->GetXaxis()->SetTitleSize(0.1);
  h1y->GetYaxis()->SetTitleSize(0.1);
  h1y->GetXaxis()->SetTitleOffset(0.7);
  h1y->GetYaxis()->SetTitleOffset(0.7);

  hDiffx->SetLineWidth(2);
  hDiffx->GetXaxis()->SetTitle("p_{x,MC}/p_{z,MC} - M_{x}"); 
  hDiffx->GetXaxis()->SetTitleSize(0.1);
  hDiffx->GetXaxis()->SetTitleOffset(0.85);
  hDiffx->GetXaxis()->SetNdivisions(505);
  hDiffx->GetYaxis()->SetNdivisions(505);
  hDiffx->GetXaxis()->SetLabelSize(0.05);
  hDiffx->GetXaxis()->SetLabelOffset(0.02);
  hDiffx->GetYaxis()->SetLabelSize(0.05);

  hDiffy->SetLineWidth(2);
  hDiffy->GetXaxis()->SetTitle("p_{y,MC}/p_{z,MC} - M_{y}"); 
  hDiffy->GetXaxis()->SetTitleSize(0.1);
  hDiffy->GetXaxis()->SetTitleOffset(0.85);
  hDiffy->GetXaxis()->SetNdivisions(505);
  hDiffy->GetYaxis()->SetNdivisions(505);
  hDiffy->GetXaxis()->SetLabelSize(0.05);
  hDiffy->GetXaxis()->SetLabelOffset(0.02);
  hDiffy->GetYaxis()->SetLabelSize(0.05);

  TH1F* hDiffx1 = hDiffx->Clone("hDiffx1");
  hDiffx1->SetLineWidth(2);
  hDiffx1->SetLineColor(kGreen+1);

  TH1F* hDiffx2 = hDiffx->Clone("hDiffx2");
  hDiffx2->SetLineWidth(2);
  hDiffx2->SetLineColor(2);

  TH1F* hDiffx3 = hDiffx->Clone("hDiffx3");
  hDiffx3->SetLineWidth(2);
  hDiffx3->SetLineColor(36);

  TH1F* hDiffy1 = hDiffy->Clone("hDiffy1");
  hDiffy1->SetLineWidth(2);
  hDiffy1->SetLineColor(kGreen+1);

  TH1F* hDiffy2 = hDiffy->Clone("hDiffy2");
  hDiffy2->SetLineWidth(2);
  hDiffy2->SetLineColor(2);

  TH1F* hDiffy3 = hDiffy->Clone("hDiffy3");
  hDiffy3->SetLineWidth(2);
  hDiffy3->SetLineColor(36);

  TCut goodChi2 = "chi2sl>0&&(chi2sl/ndfsl)<4";
  TCut withVtxLayers = "nhitsb>0";
  TCut noVtxLayers = "nhitsb==0";
  TCut largeNDF = "ndfsl>6";

  // Note: a cut on the ndf can get rid of grass in the reco vs. mc dists
  //
  gStyle->SetStatW();
  gStyle->SetStatH();

  c->cd(1);
  gPad->SetLeftMargin(0.15);
  gPad->SetBottomMargin(0.2);
  t->Project("h0x","mxsl:px0mc/pz0mc",goodChi2);
  h0x->Draw("zcol");

  c->cd(2);
  gPad->SetLeftMargin(0.15);
  gPad->SetBottomMargin(0.2);
  t->Project("h0y","mysl:py0mc/pz0mc",goodChi2);
  h0y->Draw("zcol");

  //gStyle->SetOptStat();

  gStyle->SetStatW(0.4);
  gStyle->SetStatH(0.2);
  gStyle->SetOptStat("emr");

  gStyle->SetOptStat(0);

  c->cd(3);
  gPad->SetLeftMargin(0.15);
  gPad->SetBottomMargin(0.2);
  gPad->SetLogy();
  t->Project("hDiffx","px0mc/pz0mc-mxsl",goodChi2);
  hDiffx->Draw();
  //t->Draw("px0mc/pz0mc-mxsl",goodChi2&&"ndfsl>6","same");
  //htemp->SetFillColor(42);
  //htemp->SetLineColor(kGreen);
  //htemp->SetLineWidth(2);
  t->Draw("px0mc/pz0mc-mxsl>>hDiffx1",goodChi2&&largeNDF,"same");
  t->Draw("px0mc/pz0mc-mxsl>>hDiffx2",goodChi2&&noVtxLayers,"same");
  t->Draw("px0mc/pz0mc-mxsl>>hDiffx3",goodChi2&&withVtxLayers,"same");

  TLegend* leg1 = new TLegend(0.6,0.7,0.95,0.99,"#chi^{2} / ndf < 4");
  leg1->SetFillColor(0);
  leg1->AddEntry(hDiffx,"no cuts");
  leg1->AddEntry(hDiffx3,"with VTX hits");
  leg1->AddEntry(hDiffx2,"no VTX hits");
  leg1->AddEntry(hDiffx1,"ndf > 6");
  leg1->Draw();

  c->cd(4);
  gPad->SetLeftMargin(0.15);
  gPad->SetBottomMargin(0.2);
  gPad->SetLogy();
  t->Draw("py0mc/pz0mc-mysl>>hDiffy",goodChi2);
  //hDiffy->Draw();
  t->Draw("py0mc/pz0mc-mysl>>hDiffy1",goodChi2&&largeNDF,"same");
  t->Draw("py0mc/pz0mc-mysl>>hDiffy2",goodChi2&&noVtxLayers,"same");
  t->Draw("py0mc/pz0mc-mysl>>hDiffy3",goodChi2&&withVtxLayers,"same");
  //htemp->SetFillColor(42);

  TLegend* leg2 = (TLegend*)leg1->Clone();
  leg2->Clear();
  leg2->SetHeader("#chi^{2} / ndf < 4");
  leg2->AddEntry(hDiffy,"no cuts");
  leg2->AddEntry(hDiffy3,"with VTX hits");
  leg2->AddEntry(hDiffy2,"no VTX hits");
  leg2->AddEntry(hDiffy1,"ndf > 6");
  leg2->Draw();

  //l.DrawLatex(0.6,0.5,"ndf > 6");
}

void
draw_intcmp()
{
  std::string tname = "mc_trk_eval";
  TTree* t = gROOT->FindObject(tname.c_str());
  if ( ! t ) { 
    std::cout << "Failed to find "  << tname << " tree" << std::endl;
    return;
  }

  std::string cname = "intcmpCanvas";
  TCanvas* c = gROOT->FindObject(cname.c_str());
  if ( ! c ) { 
    std::cout << "Creating " << cname << std::endl;
    c = new TCanvas(cname.c_str(),cname.c_str(),800,300);
  }
  c->Clear();
  c->Divide(2,1);

  hDiffx = new TH1F("hDiffx","MC - Reco X ints",200,-0.1,0.1);
  hDiffy = new TH1F("hDiffy","MC - Reco Y ints",200,-0.1,0.1);

  gStyle->SetOptStat("e");

  hDiffx->GetXaxis()->SetTitle("b_{x} (cm)"); 
  hDiffx->GetXaxis()->SetTitleSize(0.1);
  hDiffx->GetXaxis()->SetTitleOffset(0.9);
  hDiffx->GetXaxis()->SetNdivisions(505);
  hDiffx->GetXaxis()->SetLabelSize(0.07);
  hDiffx->GetYaxis()->SetNdivisions(505);
  hDiffx->GetYaxis()->SetLabelSize(0.07);

  hDiffy->GetXaxis()->SetTitle("b_{y} (cm)"); 
  hDiffy->GetXaxis()->SetLabelSize(0.07);
  hDiffy->GetXaxis()->SetTitleSize(0.1);
  hDiffy->GetXaxis()->SetTitleOffset(0.9);
  hDiffy->GetXaxis()->SetNdivisions(505);
  hDiffy->GetYaxis()->SetLabelSize(0.07);
  hDiffy->GetYaxis()->SetNdivisions(505);

  // Note: a cut on the ndf can get rid of grass in the reco vs. mc dists
  //

  gStyle->SetOptStat();

  c->cd(1);
  gPad->SetLeftMargin(0.15);
  gPad->SetBottomMargin(0.2);
  t->Project("hDiffx","x0sl","chi2sl/ndfsl<4");
  hDiffx->Draw();
  t->Draw("x0sl","chi2sl/ndfsl<4&&ndfsl>6","same");
  htemp->SetFillColor(42);

  TLatex l;
  l.SetTextColor(42);
  l.SetNDC(true);
  l.SetTextSize(0.1);
  l.DrawLatex(0.6,0.5,"ndf > 6");

  c->cd(2);
  gPad->SetLeftMargin(0.15);
  gPad->SetBottomMargin(0.2);
  t->Project("hDiffy","y0sl","chi2sl/ndfsl<4");
  hDiffy->Draw();
  t->Draw("y0sl","chi2sl/ndfsl<4&&ndfsl>6","same");
  htemp->SetFillColor(42);
  l.DrawLatex(0.6,0.5,"ndf > 6");

}