//
// This file is basically a wrapper for the Fortran program
// contained in 2pi_event_gen_v7.f.
//
// Fortran generator code - Rory Miskimen
// C wrapper - David Lawrence
//

#include <stdio.h>
#include <unistd.h>
#include <stdlib.h>

#include "cppvars.h"

double genevents_(void);

void SelectTarget(double ztgt, double Atgt, double tgtlen, double radlen);
void SelectTarget2(double ztgt, double Atgt, double tgtlen, double radlen, double c_den, double a_den);
void SelectRhoPhotoProdCrossSection(double sigma_0, double At);

//----------------------------
// main
//----------------------------
int main(int narg, char* argv[])
{
	// Defaults
	int target_Z = 50;
	misc_.num_events = 100000;
	misc_.num_test_events = 1000000;
	misc_.rand_iseed = 3457;
	initialstate_.Ebeam = 5.5;
	initialstate_.polarization = 1.0;
	initialstate_.gamma_beam_rate = 1.0E7;
	initialstate_.pion_or_muon = 1;
	initialstate_.rho_background = 0;
	rhob_.rho_width_model = 1;

	// Parse the command line arguments
	for(int iarg = 1; iarg<narg; iarg++){
	
	}
	
	// Set target parameters
	switch(target_Z){
		case 1:  SelectTarget ( 1.,   1., .0338, 63.04); break;               // !based on 30 cm target
		case 6:  SelectTarget ( 6.,  12., .05  , 42.70); break;               // 
		case 14: SelectTarget2(14.,  28., .05  , 21.82, 3.14 , .537 ); break; // !units of c_den and a_den are fm
		case 20: SelectTarget2(20.,  40., .05  , 16.12, 3.51 , .563 ); break; // !target length is in units of RL
		case 26: SelectTarget2(26.,  56., .05  , 13.84, 3.971, .5935); break; // !RL is g/cm^2
		case 50: SelectTarget2(50., 116., .05  ,  8.82, 5.416, .552 ); break; //
		case 82: SelectTarget (82., 208., .05  ,  6.37); break;               // !90% cut is at 0.6 degree
	}
	
	// Set Rho photoproduction cross section parameters
	switch(target_Z){
//		case  1: SelectRhoPhotoProdCrossSection(109.e3,    6.0); break; // !Ballam on proton nb/GeV^2, Gev^-2
		case  1: SelectRhoPhotoProdCrossSection(130.e3,    7.9); break; // !ABBHHM on proton nb/GeV^2
		case  6: SelectRhoPhotoProdCrossSection(  0.  ,    0.0); break; // Carbon
		case 14: SelectRhoPhotoProdCrossSection( 37.9e6, 100. ); break; // !Alvensleben on silicon, extrapolated
		case 20: SelectRhoPhotoProdCrossSection( 67.0e6, 102. ); break; // !Alvensleben on 40Ca, extrapolated , should be 102
		case 26: SelectRhoPhotoProdCrossSection(  0.  ,    0.0); break; // Iron
		case 50: SelectRhoPhotoProdCrossSection(312.e6,  196. ); break; // !Alvensleben on 116Sn, extrapolated
		case 82: SelectRhoPhotoProdCrossSection(647.e6,  290. ); break; // !Alvensleben on lead nb/GeV^2
	}	

	// Call generator
	genevents_();
	
	// Inform user of result

	return 0;
}


//----------------------------
// SelectTarget
//----------------------------
void SelectTarget(double ztgt, double Atgt, double tgtlen, double radlen)
{
	SelectTarget2(ztgt, Atgt, tgtlen, radlen, 0.0, 0.0);
}

//----------------------------
// SelectTarget2
//----------------------------
void SelectTarget2(double ztgt, double Atgt, double tgtlen, double radlen, double c_den, double a_den)
{
	tgtvars_.z_target       = ztgt;
	tgtvars_.a_target       = Atgt;
	tgtvars_.target_len     = tgtlen;
	tgtvars_.radiation_len  = radlen;
	rhob_.c_density         = c_den;
	rhob_.a_density	        = a_den;
}

//----------------------------
// SelectRhoPhotoProdCrossSection
//----------------------------
void SelectRhoPhotoProdCrossSection(double sigma_0, double At)
{
	if(sigma_0==0.0 && At==0.0){
		printf("ERROR: Values of 0 passed into SelectRhoPhotoProdCrossSection().\n");
		printf("       This probably means an invalid target z was defined.\n");
		printf("       Currently, this must be one of: 1,14,20,50,82\n");
		exit(0);
	}

	rhob_.sigma_rho_0 = sigma_0;
	rhob_.A_t         = At;
}