#include <iostream>
using namespace std;

#include <TFile.h>
#include <TH2.h>

// The following was needed to avoid undefined link errors on Mac OS X
extern "C" {void _gfortran_transfer_character_write(){cerr<<"AAHHH! We actually got called!"<<endl;}}

// Declare FORTRAN routine
extern "C" {void rezest_fdc_cdc_(float* P, float* LAMBDA, float* M, float* DP_OVER_P, float* DPHI_TOT, float* DTHETA_TOT);}

//------------------------
// rezest_fdc_cdc (wrapper routine
//------------------------
void rezest_fdc_cdc(float p, float lambda, float m, float &dp_over_p, float &dphi_tot, float &dtheta_tot)
{
	// Call FORTRAN routine
	rezest_fdc_cdc_(&p, &lambda, &m, &dp_over_p, &dphi_tot, &dtheta_tot);
}

//------------------------
// main
//------------------------
int main(int narg, char *argv[])
{
	TFile *f = new TFile("rezest.root", "RECREATE");
	
	TH1D *pres_vs_theta = new TH1D("pres_vs_theta", "#deltap/p vs #theta", 1000, 0.0, 90.0);
	TH2D *pres_vs_theta_vs_p = new TH2D("pres_vs_theta_vs_p", "#deltap/p vs #theta", 1000, 0.0, 90.0, 10, 0.0, 10.0);
	for(int pbin=1; pbin<=pres_vs_theta_vs_p->GetNbinsY(); pbin++){

		double p = pres_vs_theta_vs_p->GetYaxis()->GetBinCenter(pbin);
		
		if(pbin==2){
			char title[256];
			sprintf(title, "#deltap/p vs #theta for p=%3.1fGeV/c", p);
			pres_vs_theta->SetTitle(title);
		}
		
		for(int ibin=1; ibin<=pres_vs_theta->GetNbinsX(); ibin++){
			double theta = pres_vs_theta->GetXaxis()->GetBinCenter(ibin);
			float lambda = 3.1415927/2.0 - theta/57.3;
			
			float dp_over_p, dphi_tot, dtheta_tot;
			rezest_fdc_cdc(p, lambda, 0.139, dp_over_p, dphi_tot, dtheta_tot);
			
			if(pbin==2)pres_vs_theta->SetBinContent(ibin, dp_over_p);
			pres_vs_theta_vs_p->SetBinContent(ibin, pbin, dp_over_p);
		}
	}
	
	pres_vs_theta->SetStats(0);
	pres_vs_theta->GetYaxis()->SetRangeUser(0.0, 0.3);
	pres_vs_theta_vs_p->SetStats(0);
	pres_vs_theta_vs_p->GetZaxis()->SetRangeUser(0.0, 0.3);

	f->Write();
	f->Close();
	delete f;

	return 0;
}