{

  gROOT->Reset();

  cout << "plot title: ";
  TString title;
  cin >> title;
  cout << "title = /" << title << "/\n";
  pt = new TPaveText(0.05,0.85,0.95,0.95, "NDC"); // NDC sets coords
                                              // relative to pad dimensions
  pt->SetFillColor(0); // text is black on white
  pt->SetTextSize(0.04); 
  pt->SetTextAlign(12);
  text = pt->AddText(title);

 ifstream in, in_hd;
// we assume a file basic.dat in the current directory
// this file has 3 columns of float data
   in.open("thetascan.txt");
   in_hd.open("hdget_thetascan.txt");

   const Int_t n = 1000;
   Float_t lambda[n], theta[n]; // 1-2
   Float_t dp_p[n], dph[n], dth[n]; // 3-5
   Float_t dp_p_ms_c[n], dp_p_ms_f[n], dp_p_res_c[n], dp_p_res_f[n],
     dp_p_f[n], dp_p_c[n]; // 6-11
   Float_t dph_ms_c[n], dph_ms_f[n], dph_res_c[n], dph_res_f[n],
     dph_k_c[n], dph_k_f[n], dph_c[n], dph_f[n]; // 12-19
   Float_t dth_ms_c[n], dth_ms_f[n], dth_res_c[n], dth_res_f[n],
     dth_f[n], dth_c[n]; // 20-25
   Float_t p_hd[n], theta_hd[n], dp_p_hd[n], dth_hd[n], dph_hd[n];

   Int_t nl = 0;

   while (1) {
     in >> lambda[nl] >> theta[nl] >> dp_p[nl] >> dph[nl] >> dth[nl]
	>> dp_p_ms_c[nl] >> dp_p_ms_f[nl] >> dp_p_res_c[nl] >> dp_p_res_f[nl]
	>> dp_p_c[nl] >> dp_p_f[nl]
	>> dph_ms_c[nl] >> dph_ms_f[nl] >> dph_res_c[nl] >> dph_res_f[nl]
	>> dph_k_c[nl] >> dph_k_f[nl] >> dph_c[nl] >> dph_f[nl]
	>> dth_ms_c[nl] >> dth_ms_f[nl] >> dth_res_c[nl] >> dth_res_f[nl]
	>> dth_c[nl] >> dth_f[nl];
     if (nl%20 == 0) {
       printf("%d %f %f %f %f\n", nl, theta[nl], dp_p[nl], dph[nl], dth[nl]);
     }
     if (!in.good()) break;
     nl++;
   }
   printf(" found %d points\n",nl);

   in.close();

   Int_t nl_hd = 0;
   
   while (1) {
     in_hd >> p_hd[nl_hd] >> theta_hd[nl_hd] >> dp_p_hd[nl_hd]
	   >> dth_hd[nl_hd] >> dph_hd[nl_hd];
     dth_hd[nl_hd] *= 1.e-3;
     dph_hd[nl_hd] *= 1.e-3;
     if (nl_hd%20 == 0) {
       printf("%d %f %f %f %f %f %f %f\n", nl_hd, p_hd[nl_hd], dp_p_hd[nl_hd], dph_hd[nl_hd], dth_hd[nl_hd]);
     }
     if (!in_hd.good()) break;
     nl_hd++;
   }
   printf(" found %d points, hd\n",nl_hd);
   
   in_hd.close();
   
   c1 = new TCanvas("c1","relative transverse momentum resolution",900,700);

   pad1 = new TPad("pad1","The pad with the function",0.0,0.5,0.50,1.0,21);
   pad1->SetLeftMargin(0.15);
   pad1->SetBottomMargin(0.15);
   pad1->SetRightMargin(0.05);
   pad2 = new TPad("pad2","The pad with the histogram",0.5,0.5,1.0,1.0,21);
   pad2->SetLeftMargin(0.15);
   pad2->SetBottomMargin(0.15);
   pad2->SetRightMargin(0.05);
   pad3 = new TPad("pad3","The pad with the histogram",0.0,0.0,0.50,0.50,21);
   pad3->SetLeftMargin(0.15);
   pad3->SetBottomMargin(0.15);
   pad3->SetRightMargin(0.05);
   pad4 = new TPad("pad3","The pad with the histogram",0.5,0.0,1.0,0.50,21);

   pad1->Draw();
   pad2->Draw();
   pad3->Draw();
   pad4->Draw();

   pad1->cd();

   gr1 = new TGraph(nl, theta, dp_p_c);
   gr1->SetLineColor(2);
   gr1->SetLineWidth(3);

   gr2 = new TGraph(nl, theta, dp_p_ms_c);
   gr2->SetLineColor(4);
   gr2->SetLineWidth(3);

   gr3 = new TGraph(nl, theta, dp_p_res_c);
   gr3->SetLineColor(5);
   gr3->SetLineWidth(3);
   gr3->SetTitle("CDC");
   gr3->GetXaxis()->SetTitle("theta (degrees)");
   gr3->GetYaxis()->SetTitle("dpt/pt");
   gr3->GetXaxis()->SetRangeUser(0.0, 60.0);
   gr3->GetYaxis()->SetRangeUser(0, 0.08);
   gr3->GetXaxis()->SetTitleSize(0.06);
   gr3->GetYaxis()->SetTitleSize(0.06);
   gr3->GetXaxis()->SetLabelSize(0.06);
   gr3->GetYaxis()->SetLabelSize(0.06);

   gr3->Draw("AL");
   gr2->Draw("L");
   gr1->Draw("L");

   pad2->cd();

   gr4 = new TGraph(nl, theta, dp_p_f);
   gr4->SetLineColor(3);
   gr4->SetLineWidth(3);

   gr5 = new TGraph(nl, theta, dp_p_ms_f);
   gr5->SetLineColor(4);
   gr5->SetLineWidth(3);

   gr6 = new TGraph(nl, theta, dp_p_res_f);
   gr6->SetLineColor(5);
   gr6->SetLineWidth(3);
   gr6->SetTitle("FDC");
   gr6->GetXaxis()->SetTitle("theta (degrees)");
   gr6->GetYaxis()->SetTitle("dpt/pt");
   gr6->GetXaxis()->SetRangeUser(0.0, 60.0);
   gr6->GetYaxis()->SetRangeUser(0, 0.08);
   gr6->GetXaxis()->SetTitleSize(0.06);
   gr6->GetYaxis()->SetTitleSize(0.06);
   gr6->GetXaxis()->SetLabelSize(0.06);
   gr6->GetYaxis()->SetLabelSize(0.06);

   gr6->Draw("AL");
   gr5->Draw("L");
   gr4->Draw("L");

   pad3->cd();

   gr7 = new TGraph(nl, theta, dp_p);
   gr7->SetLineColor(1);
   gr7->SetLineWidth(3);

   gr8 = new TGraph(nl, theta, dp_p_c);
   gr8->SetLineColor(2);
   gr8->SetLineWidth(3);

   gr9 = new TGraph(nl, theta, dp_p_f);
   gr9->SetLineColor(3);
   gr9->SetLineWidth(3);
   gr9->SetTitle("Combined CDC and FDC");
   gr9->GetXaxis()->SetTitle("theta (degrees)");
   gr9->GetYaxis()->SetTitle("dpt/pt");
   gr9->GetXaxis()->SetRangeUser(0.0, 60.0);
   gr9->GetYaxis()->SetRangeUser(0.0, 0.08);
   gr9->GetXaxis()->SetTitleSize(0.06);
   gr9->GetYaxis()->SetTitleSize(0.06);
   gr9->GetXaxis()->SetLabelSize(0.06);
   gr9->GetYaxis()->SetLabelSize(0.06);

   gr10 = new TGraph(nl_hd, theta_hd, dp_p_hd);
   gr10->SetLineColor(7);
   gr10->SetLineWidth(3);

   gr9->Draw("AL");
   gr8->Draw("L");
   gr7->Draw("L");
   gr10->Draw("L");

   pad4->cd();

   pt->Draw();

   leg = new TLegend(0.05,0.05,0.95,0.8);
   leg->SetHeader("relative transverse momentum resolution");
   leg->AddEntry(gr3,"due to position resolution","l");
   leg->AddEntry(gr2,"due to multiple scattering","l");
   leg->AddEntry(gr8,"CDC only, pos+ms","l");
   leg->AddEntry(gr9,"FDC only, pos+ms","l");
   leg->AddEntry(gr7,"combined","l");
   leg->AddEntry(gr10,"HDGEANT","l");
   leg->Draw();

   c1->Print("thetascan_pt.png", "png");

   c2 = new TCanvas("c2","azimuthal angular resolution",900,700);

   pad1 = new TPad("pad1","The pad with the function",0.0,0.5,0.50,1.0,21);
   pad1->SetLeftMargin(0.15);
   pad1->SetBottomMargin(0.15);
   pad1->SetRightMargin(0.05);
   pad2 = new TPad("pad2","The pad with the histogram",0.5,0.5,1.0,1.0,21);
   pad2->SetLeftMargin(0.15);
   pad2->SetBottomMargin(0.15);
   pad2->SetRightMargin(0.05);
   pad3 = new TPad("pad3","The pad with the histogram",0.0,0.0,0.50,0.50,21);
   pad3->SetLeftMargin(0.15);
   pad3->SetBottomMargin(0.15);
   pad3->SetRightMargin(0.05);
   pad4 = new TPad("pad3","The pad with the histogram",0.5,0.0,1.0,0.50,21);
   pad1->Draw();
   pad2->Draw();
   pad3->Draw();
   pad4->Draw();

   pad1->cd();
   pad1->SetLogy(1);

   gr1 = new TGraph(nl, theta, dph_c);
   gr1->SetLineColor(2);
   gr1->SetLineWidth(2);

   gr2 = new TGraph(nl, theta, dph_ms_c);
   gr2->SetLineColor(4);
   gr2->SetLineWidth(3);

   gr3 = new TGraph(nl, theta, dph_res_c);
   gr3->SetLineColor(5);
   gr3->SetLineWidth(3);

   gr4 = new TGraph(nl, theta, dph_k_c);
   gr4->SetLineColor(6);
   gr4->SetLineWidth(3);
   gr4->SetTitle("CDC");
   gr4->GetXaxis()->SetTitle("theta (degrees)");
   gr4->GetYaxis()->SetTitle("dphi (rad)");
   gr4->GetYaxis()->SetRangeUser(1.e-4, 1.e-0);
   gr4->GetXaxis()->SetTitleSize(0.06);
   gr4->GetYaxis()->SetTitleSize(0.06);
   gr4->GetXaxis()->SetLabelSize(0.06);
   gr4->GetYaxis()->SetLabelSize(0.06);

   gr4->Draw("AL");
   gr3->Draw("L");
   gr2->Draw("L");
   gr1->Draw("L");

   pad2->cd();
   pad2->SetLogy(1);

   gr5 = new TGraph(nl, theta, dph_f);
   gr5->SetLineColor(3);
   gr5->SetLineWidth(2);

   gr6 = new TGraph(nl, theta, dph_ms_f);
   gr6->SetLineColor(4);
   gr6->SetLineWidth(3);

   gr7 = new TGraph(nl, theta, dph_res_f);
   gr7->SetLineColor(5);
   gr7->SetLineWidth(3);

   gr8 = new TGraph(nl, theta, dph_k_f);
   gr8->SetLineColor(6);
   gr8->SetLineWidth(3);
   gr8->SetTitle("FDC");
   gr8->GetXaxis()->SetTitle("theta (degrees)");
   gr8->GetYaxis()->SetTitle("dphi (rad)");
   gr8->GetYaxis()->SetRangeUser(1.e-4, 1.e-0);
   gr8->GetXaxis()->SetTitleSize(0.06);
   gr8->GetYaxis()->SetTitleSize(0.06);
   gr8->GetXaxis()->SetLabelSize(0.06);
   gr8->GetYaxis()->SetLabelSize(0.06);

   gr8->Draw("AL");
   gr7->Draw("L");
   gr6->Draw("L");
   gr5->Draw("L");

   pad3->cd();
   pad3->SetLogy(1);

   gr9 = new TGraph(nl, theta, dph);
   gr9->SetLineColor(1);
   gr9->SetLineWidth(3);

   gr10 = new TGraph(nl, theta, dph_c);
   gr10->SetLineColor(2);
   gr10->SetLineWidth(3);

   gr11 = new TGraph(nl, theta, dph_f);
   gr11->SetLineColor(3);
   gr11->SetLineWidth(3);
   gr11->SetTitle("Combined CDC and FDC");
   gr11->GetXaxis()->SetTitle("theta (degrees)");
   gr11->GetYaxis()->SetTitle("dphi (rad)");
   gr11->GetYaxis()->SetRangeUser(1.e-3, 1.e-0);
   gr11->GetXaxis()->SetTitleSize(0.06);
   gr11->GetYaxis()->SetTitleSize(0.06);
   gr11->GetXaxis()->SetLabelSize(0.06);
   gr11->GetYaxis()->SetLabelSize(0.06);

   gr12 = new TGraph(nl_hd, theta_hd, dph_hd);
   gr12->SetLineColor(7);
   gr12->SetLineWidth(3);

   gr11->Draw("AL");
   gr10->Draw("L");
   gr9->Draw("L");
   gr12->Draw("L");

   pad4->cd();

   pt->Draw();

   leg = new TLegend(0.05,0.05,0.95,0.8);
   leg->SetHeader("azimuthal angular resolution");
   leg->AddEntry(gr3,"due to position resolution","l");
   leg->AddEntry(gr2,"due to multiple scattering","l");
   leg->AddEntry(gr4,"due to curvature resolution","l");
   leg->AddEntry(gr10,"CDC only, pos+ms","l");
   leg->AddEntry(gr11,"FDC only, pos+ms","l");
   leg->AddEntry(gr9,"combined","l");
   leg->AddEntry(gr12,"HDGEANT","l");
   leg->Draw();

   c2->Print("thetascan_azimuth.png", "png");

   c3 = new TCanvas("c3","polar angle resolution",900,700);

   pad1 = new TPad("pad1","CDC",0.0,0.5,0.50,1.0,21);
   pad1->SetLeftMargin(0.15);
   pad1->SetBottomMargin(0.15);
   pad1->SetRightMargin(0.05);
   pad2 = new TPad("pad2","FDC",0.5,0.5,1.0,1.0,21);
   pad2->SetLeftMargin(0.15);
   pad2->SetBottomMargin(0.15);
   pad2->SetRightMargin(0.05);
   pad3 = new TPad("pad3","Combined",0.0,0.0,0.50,0.50,21);
   pad3->SetLeftMargin(0.15);
   pad3->SetBottomMargin(0.15);
   pad3->SetRightMargin(0.05);
   pad4 = new TPad("pad3","The pad with the histogram",0.5,0.0,1.0,0.50,21);
   pad1->Draw();
   pad2->Draw();
   pad3->Draw();
   pad4->Draw();

   pad1->cd();
   pad1->SetLogy(1);

   gr1 = new TGraph(nl, theta, dth_c);
   gr1->SetLineColor(2);
   gr1->SetLineWidth(3);

   gr2 = new TGraph(nl, theta, dth_ms_c);
   gr2->SetLineColor(4);
   gr2->SetLineWidth(3);

   gr3 = new TGraph(nl, theta, dth_res_c);
   gr3->SetLineColor(5);
   gr3->SetLineWidth(3);
   gr3->SetTitle("CDC");
   gr3->GetXaxis()->SetTitle("theta (degrees)");
   gr3->GetYaxis()->SetTitle("dtheta (rad)");
   gr3->GetYaxis()->SetRangeUser(1.e-4, 1.e-2);
   gr3->GetXaxis()->SetTitleSize(0.06);
   gr3->GetYaxis()->SetTitleSize(0.06);
   gr3->GetXaxis()->SetLabelSize(0.06);
   gr3->GetYaxis()->SetLabelSize(0.06);

   gr3->Draw("AL");
   gr2->Draw("L");
   gr1->Draw("L");

   pad2->cd();
   pad2->SetLogy(1);

   gr4 = new TGraph(nl, theta, dth_f);
   gr4->SetLineColor(3);
   gr4->SetLineWidth(3);

   gr5 = new TGraph(nl, theta, dth_ms_f);
   gr5->SetLineColor(4);
   gr5->SetLineWidth(3);

   gr6 = new TGraph(nl, theta, dth_res_f);
   gr6->SetLineColor(5);
   gr6->SetLineWidth(3);
   gr6->SetTitle("FDC");
   gr6->GetXaxis()->SetTitle("theta (degrees)");
   gr6->GetYaxis()->SetTitle("dtheta (rad)");
   gr6->GetYaxis()->SetRangeUser(1.e-4, 1.e-2);
   gr6->GetXaxis()->SetTitleSize(0.06);
   gr6->GetYaxis()->SetTitleSize(0.06);
   gr6->GetXaxis()->SetLabelSize(0.06);
   gr6->GetYaxis()->SetLabelSize(0.06);

   gr6->Draw("AL");
   gr5->Draw("L");
   gr4->Draw("L");

   pad3->cd();
   pad3->SetLogy(1);

   gr7 = new TGraph(nl, theta, dth);
   gr7->SetLineColor(1);
   gr7->SetLineWidth(3);

   gr8 = new TGraph(nl, theta, dth_c);
   gr8->SetLineColor(2);
   gr8->SetLineWidth(3);

   gr9 = new TGraph(nl, theta, dth_f);
   gr9->SetLineColor(3);
   gr9->SetLineWidth(3);
   gr9->SetTitle("Combined CDC and FDC");
   gr9->GetXaxis()->SetTitle("theta (degrees)");
   gr9->GetYaxis()->SetTitle("dtheta (rad)");
   gr9->GetYaxis()->SetRangeUser(1.e-4, 1.e-2);
   gr9->GetXaxis()->SetTitleSize(0.06);
   gr9->GetYaxis()->SetTitleSize(0.06);
   gr9->GetXaxis()->SetLabelSize(0.06);
   gr9->GetYaxis()->SetLabelSize(0.06);

   gr10 = new TGraph(nl_hd, theta_hd, dth_hd);
   gr10->SetLineColor(7);
   gr10->SetLineWidth(3);

   gr9->Draw("AL");
   gr8->Draw("L");
   gr7->Draw("L");
   gr10->Draw("L");

   pad4->cd();

   pt->Draw();

   leg = new TLegend(0.05,0.05,0.95,0.8);
   leg->SetHeader("polar angular resolution");
   leg->AddEntry(gr3,"due to position resolution","l");
   leg->AddEntry(gr2,"due to multiple scattering","l");
   leg->AddEntry(gr8,"CDC only, pos+ms","l");
   leg->AddEntry(gr9,"FDC only, pos+ms","l");
   leg->AddEntry(gr7,"combined","l");
   leg->AddEntry(gr10,"HDGEANT","l");
   leg->Draw();

   c3->Print("thetascan_polar.png", "png");

}