#if 0
class track:public TObject{

	public:
		track(){}

		TVector3 pthrown;
		TVector3 pfit;
		double likelihood;
		double chisq;
		double pt_pull;
		double theta_pull;
		double phi_pull;
		bool isreconstructable;
		int Nstereo;
		int Ncdc;
		int Nfdc;

	private:
		ClassDef(track,1);

};
#endif
//--------------------
// efficiency
//--------------------
void efficiency(void)
{
	gROOT->Reset();
	gStyle->SetPalette(1,NULL);

	// Chain together root files
	TFile *ifile = new TFile("hd_root.root");
	ifile->cd("TRACKING");
	TTree *trkeff = (TTree*)gROOT->FindObject("trkeff");
	
	TCanvas *c1 = new TCanvas("c1");
	c1->SetTickx();
	c1->SetTicky();
	
	// Create output file to store the histograms in
	TFile *f = new TFile("efficiency.root","RECREATE");
	
	// It turns out that the resolutions used to determine the pulls
	// could be off which in turn, leads to an inaccurate chisq
	// and an inaccurate efficiency. To avoid this, we first fit the
	// pulls to gaussians in order to determine thier sigmas (which
	// really should be 1). These are then used to scale the event
	// by event pulls such that they are equivalent to sigma=1 and
	// the chisq re-calculated for use in determining the efficiency.
	double sigma_pt_pull = 1.0;
	double sigma_theta_pull = 1.0;
	double sigma_phi_pull = 1.0;
	
	TCut cut1 = "isreconstructable==1";
	
	// sigma pt_pull
	TH1D *pt_pull = new TH1D("pt_pull","",200, -10.0, 10.0);
	trkeff->Project("pt_pull","pt_pull", cut1);
	pt_pull->Fit("gaus","0Q");
	sigma_pt_pull = pt_pull->GetFunction("gaus")->GetParameter(2);
	
	// sigma sigma_theta_pull
	TH1D *theta_pull = new TH1D("theta_pull","",200, -10.0, 10.0);
	trkeff->Project("theta_pull","theta_pull", cut1);
	theta_pull->Fit("gaus","0Q");
	sigma_theta_pull = theta_pull->GetFunction("gaus")->GetParameter(2);
	
	// sigma pt_pull
	TH1D *phi_pull = new TH1D("phi_pull","",200, -10.0, 10.0);
	trkeff->Project("phi_pull","phi_pull", cut1);
	phi_pull->Fit("gaus","0Q");
	sigma_phi_pull = phi_pull->GetFunction("gaus")->GetParameter(2);
	
	cout<<"sigma_pt_pull="<<sigma_pt_pull<<" sigma_theta_pull="<<sigma_theta_pull<<" sigma_phi_pull="<<sigma_phi_pull<<endl;

	char chisq_str[256];
	char cut_str[256];
	sprintf(chisq_str, "((pow(pt_pull/%f,2.0) + pow(theta_pull/%f,2.0) + pow(phi_pull/%f,2.0))/3.0)", sigma_pt_pull, sigma_theta_pull, sigma_phi_pull);	
	sprintf(cut_str, "%s<1000.0", chisq_str);	
	TCut cut2 = cut_str;
	
	TH1D *chisq_corrected = new TH1D("chisq_corrected","#chi^{2} corrected", 10001, 0.0, 1000.0);
	trkeff->Project("chisq_corrected", chisq_str, cut1);

	TH2D *found_vs_pt_vs_theta = new TH2D("found_vs_pt_vs_theta","", 160, 0.0, 160.0, 50, 0.0, 2.0);
	TH2D *thrown_vs_pt_vs_theta = found_vs_pt_vs_theta->Clone("thrown_vs_pt_vs_theta");
	TH2D *eff_vs_pt_vs_theta = found_vs_pt_vs_theta->Clone("eff_vs_pt_vs_theta");
	
	trkeff->Project("found_vs_pt_vs_theta", "E.pthrown.Perp():E.pthrown.Theta()*57.3", cut1 && cut2);
	trkeff->Project("thrown_vs_pt_vs_theta", "E.pthrown.Perp():E.pthrown.Theta()*57.3", cut1);
	eff_vs_pt_vs_theta->Divide(found_vs_pt_vs_theta, thrown_vs_pt_vs_theta);
	
	eff_vs_pt_vs_theta->SetStats(0);
	eff_vs_pt_vs_theta->SetTitle("CDC Tracking Efficiency vs. p_{T} and #theta");
	eff_vs_pt_vs_theta->SetXTitle("#theta angle (degrees)");
	eff_vs_pt_vs_theta->SetYTitle("Transverse momentum (GeV/c)");
	eff_vs_pt_vs_theta->GetZaxis()->SetRangeUser(0.0, 1.0);
	eff_vs_pt_vs_theta->Draw("colz");
	
	AddStandardLabels(eff_vs_pt_vs_theta);
	
	c1->SaveAs("efficiency.gif");
	c1->SaveAs("efficiency.pdf");
	
	c1->Clear();
	c1->Divide(2,2);
	c1->cd(1);
	pt_pull->Fit("gaus");
	c1->cd(2);
	theta_pull->Fit("gaus");
	c1->cd(3);
	phi_pull->Fit("gaus");
	c1->cd(4);
	c1->GetPad(4)->SetLogy();
	chisq_corrected->GetXaxis()->SetRangeUser(0.0, 700.0);
	chisq_corrected->Draw();

	c1->SaveAs("pulls.gif");
	c1->SaveAs("pulls.pdf");
	
	int maxbin = eff_vs_pt_vs_theta->GetMaximumBin();
	cout<<"Maximum efficiency:"<<eff_vs_pt_vs_theta->GetBinContent(maxbin)<<endl;
	
	f->Write();
	delete f;
	
}


//----------------
// AddStandardLabels
//----------------
void AddStandardLabels(TH2D *axes=NULL)
{
	// This will draw a label or two on the
	// current plot using the NDC coordinates.
	// It is put here to make sure all plots have
	// a consistent labeling.
	
	// Get Date
	TDatime *t = new TDatime();
	const char * month[]={"Jan.", "Feb.", "March", "April", "May", "June", "July", "Aug.", "Sep.", "Oct.", "Nov.", "Dec."};
	char date[256];
	sprintf(date, "%s %d, %d DL", month[t->GetMonth()-1], t->GetDay(), t->GetYear());
	
	// Date, Author
	TLatex *lab = new TLatex(0.7, 0.7, date);
	ConvertFromNDC(lab, axes);
	lab->SetTextSize(0.03);
	lab->SetTextAlign(33);
	lab->Draw();

	// SVN Revision
	lab = new TLatex(0.7, 0.645, "svn revision: current");
	ConvertFromNDC(lab, axes);
	lab->SetTextSize(0.02);
	lab->SetTextAlign(31);
	lab->Draw();

	// Candidate type
	lab = new TLatex(0.7, 0.635, "Candidates from finder");
	ConvertFromNDC(lab, axes);
	lab->SetTextSize(0.02);
	lab->SetTextAlign(32);
	lab->Draw();
	
	// Event type
	lab = new TLatex(0.7, 0.59, "Single #pi^{+} with MULS, LOSS, DRAY off");
	ConvertFromNDC(lab, axes);
	lab->SetTextSize(0.03);
	lab->SetTextAlign(31);
	lab->Draw();
}

//----------------
// ConvertFromNDC
//----------------
void ConvertFromNDC(TLatex *obj, TH2D *h=NULL)
{
	// Bugs in ROOT make it hard to plot labels consistently.
	// For 1D plots, the histogram axes define the coordinate
	// system. For 2D plots, we seem to be forced to use the
	// NDC. There does not seem to be an obvious way to tell
	// which we're using so we pass the information in in the
	// form of the "axes" histogram. If it is not NULL, then
	// we use it to define the limits. Otherwise, we do nothing.

	if(h==NULL)return;

	TAxis *xaxis = h->GetXaxis();
	int Nbinsx = xaxis->GetNbins();
	double xmin = xaxis->GetBinLowEdge(1);
	double xmax = xmin + xaxis->GetBinLowEdge(Nbinsx);

	TAxis *yaxis = h->GetYaxis();
	int Nbinsy = yaxis->GetNbins();
	double ymin = yaxis->GetBinLowEdge(1);
	double ymax = ymin + yaxis->GetBinLowEdge(Nbinsy);
	
	
	double x = obj->GetX();
	double y = obj->GetY();
	cout<<" in:  x="<<x<<"  y="<<y<<endl;
	
	x = xmin + (xmax-xmin)*(0.5+x/1.15);
	y = ymin + (ymax-ymin)*(0.5+y/1.15);
	
	obj->SetX(x);
	obj->SetY(y);

	cout<<"out:  x="<<x<<"  y="<<y<<endl;
}