#include "DetResponse.hh"
#include "Table.hh"

void DetResponse::dus() {}


Table DefaultDetResponse::compute_rates_table ( Flux& flux,                             
												double livetime, // seconds
												bool point, 
												double declination, 
												double cutoff_angle,
												double cutoff_logerec )
{
	Table T("log-E","Flux","Aeff","RawRate","EffAng","EffEnergy","Rate","Bg-Rate");

	T.title = string("  Event Rate for "+ channel_name(chan) + " " + flavor_name (flav) );
	if ( point ) T.title += " point source at decl=" + str( declination);
	T.title += "angle/E-cut = " + str(cutoff_angle) + " " + str( cutoff_logerec );

	auto res = compute_rates( flux, livetime, point, declination, cutoff_angle, cutoff_logerec, &T );
	double totrateeff = res[0];
	double totrate    = res[2];
	double bg         = res[1];
	
	T << "total" << "-" << "-" << totrate << "-" << "-" << totrateeff << bg ;
	return T;
}

vector<double> DefaultDetResponse::compute_rates(Flux& flux,                             
												double livetime, // seconds
												bool point, 
												double declination, 
												double cutoff_angle,
												double cutoff_logerec,
												Table* T  )
{
	double totrate(0), totrateeff(0);
	if (this->flav != muon) {
		EffectiveArea1D A;

		if ( point ) {
			A =  Aeff.GetAreaAtDecl( declination );
		}
		else {
			A = Aeff.GetAverage();
		}

		TH1D h = A.GetRate( flux, livetime );

		for (int bx = 1; bx <= h.GetNbinsX() ; bx++ ) 
		{
			const double loge = h.GetBinCenter(bx);
			const double flx  = flux.dNdE( flav, loge );
			const double aeff = A.H.GetBinContent(bx); // big H = Aeff

			const double rate = h.GetBinContent(bx);
			const double eff_ang  = psf.fraction_below_angle( loge, cutoff_angle );

			const double sindecl = sin(declination);
			const double eff_erec = erec_fraction_above_sindecl( loge, sindecl, cutoff_logerec);

			if (T) (*T) << loge << flx << aeff << rate << eff_ang << eff_erec << rate * eff_erec * eff_ang << "-";

			totrate += rate;
			totrateeff += rate * eff_erec * eff_ang;
		}
	}

	const double bg = background_rate( livetime, declination, cutoff_angle, cutoff_logerec );

	return {totrateeff, bg , totrate };
}

vector<double> DefaultDetResponse::compute_rates_angle_energy_range(Flux& flux,                             
												double livetime, // seconds
												bool point, 
												double declination, 
												double angle_min,
												double angle_max,
												double logerec_min,
												double logerec_max  )
{

	double totrate(0), totrateeff(0);

	if (this->flav != muon) {

		EffectiveArea1D A;

		if ( point ) {
			A =  Aeff.GetAreaAtDecl( declination );
		}
		else {
			A = Aeff.GetAverage();
		}

		TH1D h = A.GetRate( flux, livetime );


		for (int bx = 1; bx <= h.GetNbinsX() ; bx++ ) 
		{
			const double loge = h.GetBinCenter(bx);
			const double flx  = flux.dNdE( flav, loge );
			const double aeff = A.H.GetBinContent(bx); // big H = Aeff

			const double rate = h.GetBinContent(bx);
			const double eff_ang  = psf.fraction_between_angle( loge, angle_min, angle_max );

			const double sindecl = sin(declination);
			const double eff_erec = erec_fraction_range_above_sindecl( loge, sindecl, logerec_min, logerec_max);

			// cout << "eff_ang " << eff_ang << " eff_erec " << eff_erec << endl; 

			totrate += rate;
			totrateeff += rate * eff_erec * eff_ang;
		}

	}

	const double bg = background_rate_range( livetime, declination, angle_min, angle_max, logerec_min, logerec_max );


	return {totrateeff, bg , totrate };
}


void FullDetResponse::show_rates(Flux& flux, 
	                             double livetime,
	                             bool point, 
	                             double declination, 
	                             double cutoff_angle,
	                             double cutoff_logerec )
{
	for (auto flav : flavors )
	{
		for ( auto chan : channels )
		{
			Table T = get( flav, chan ).compute_rates_table( flux, livetime, point, declination, cutoff_angle, cutoff_logerec );
			print(T);
		}
	}

}

void FullDetResponse::compute_rates_summary(Flux& flux, bool point, double declination, double cutoff_angle )
{

	Table T( "flavour" );
	for (auto c : channels ) T.add_column( "Rate-" + channel_name(c) ) ;

	for ( auto flav : flavors) {
		if (flav != muon ) {

			T << flavor_name( flav );

			for (auto chan : channels ) {

				EffectiveArea1D A;

				if ( point ) {
					A =  Aeff( flav, chan ).GetAreaAtDecl( declination );
				}
				else {
					A = Aeff( flav, chan ).GetAverage();
				}

				double rate  = A.GetTotalRate( flux );
				T << rate;
			}
		}
	}
	print(T);
}