#ifndef BACKGROUNDCOMPONENT_HH
#define BACKGROUNDCOMPONENT_HH

#include "Component.hh"
#include "dayfraction.hh"
#include "Mat.hh"    

/*! A background-component is meant for e.g. atmspheric neutrino  background.
    Such a contribution is typically taking from data, with maybe the normalisation
    a free parameter. */
  

struct BackgroundComponent : public Component
{  

	FullDetResponse& detres;

	map< channel, TH2D> hbkg; // histogram of all the background for this channel (sum over flavors)

	map< channel, TH1D> hbkg_erec; // 1d pdf of erec only.

	map< channel, TH3D >          _pdf_erec_dayfraction_decl; // pdf of Erec and time-of-day vs decl
	map< channel, vector<TH2D> >  _pdf_erec_dayfraction_decl_slices; 

	map< channel, TH1D > _dP_dOmega_decl;          // number of events found at decl (modulo norm)
	map< channel, TH2D > _dP_dOmega_dlogerec_decl; // integrated over time.

	vector<double> _dp_domega_dloge;               // search-event cache

	vector< channel> channels;


	
	BackgroundComponent( string name, FullDetResponse& detr ) : Component(name), detres(detr)
	{
		// add up all the background channels

		norm.value = 0;


		for ( auto chan : detr.channels )
		{
			for (auto flav : detr.flavors )
			{
				if( startswith(flavor_name(flav),"a") || startswith(flavor_name(flav),"n" )){
					if ( hbkg.find( chan ) == hbkg.end() ) {
						hbkg[chan] = detr.get( flav, chan ).h_BKGREC["ATMNU_conv"]; // copy
					}
					else {
						hbkg[chan].Add( & detr.get( flav, chan ).h_BKGREC["ATMNU_conv"] );
					}
					hbkg[chan].Add( & detr.get( flav, chan ).h_BKGREC["ATMNU_prompt"] );
				} else {
					if ( hbkg.find( chan ) == hbkg.end() ) {
						hbkg[chan] = detr.get( flav, chan ).h_BKGREC["ATMMUON"];
					}
					else {
						hbkg[chan].Add( & detr.get( flav, chan ).h_BKGREC["ATMMUON"] );
					}
				}
			}

			norm.value += hbkg[chan].Integral();
			normalize( hbkg[chan] );
			init_cos_zenith_decl(chan);
		}

		norm.fit_max = norm.value * 1.3;
		norm.fit_min = norm.value * 0.7;

		info();
	}


	void init_cos_zenith_decl(channel chan)
	{
		_pdf_erec_dayfraction_decl[chan] = TH3D("_pdf_erec_dayfraction_decl",
										  "_pdf_erec_dayfraction_decl",
										  hbkg[chan].GetNbinsX(), hbkg[chan].GetXaxis()->GetXmin(), hbkg[chan].GetXaxis()->GetXmax(),
										  100,0,1,
										  40,-1,1);

		//_rate_at_sindecl = TH1D( "_rate_at_sindecl","background rate_vs_decl", 40,-1,1 );
		
		print ("initializing background vs declination histogram");

		TH3D& H = _pdf_erec_dayfraction_decl[chan];

		for (int bx =1 ; bx <= H.GetXaxis()->GetNbins() ; bx ++ )
		{
			double logerec = H.GetXaxis()->GetBinCenter(bx);

			print (bx, "/", H.GetXaxis()->GetNbins() );


			for (int bz =1 ; bz <= H.GetZaxis()->GetNbins() ; bz ++ )
			{
				double sindecl = H.GetZaxis()->GetBinCenter( bz );


				for (int by =1 ; by <= H.GetYaxis()->GetNbins() ; by ++ )
				{
					double dayf = H.GetYaxis()->GetBinCenter( by );

					double costheta = dayfraction_to_coszenith( 0, asin(sindecl), detres.det(), dayf ); // ra = 0 on purpose
					
					double v = dP_dOmega_dlogE( chan, logerec, costheta ); 

					int bin = H.GetBin( bx, by, bz );
					H.SetBinContent( bin , v / H.GetYaxis()->GetNbins() );

				} // by
			}//bz
		}//bx

		_pdf_erec_dayfraction_decl_slices[chan] = split_z( _pdf_erec_dayfraction_decl[chan] );

		_dP_dOmega_dlogerec_decl[chan] = * (( TH2D*) H.Project3D("ZX"));
		//_dP_dOmega_dlogerec_decl[chan].Scale( H.GetXaxis()->GetBinWidth(1) );

		_dP_dOmega_decl[chan] = *(H.ProjectionZ());
		_dP_dOmega_decl[chan].Scale( H.GetXaxis()->GetBinWidth(1) );
	}


	TH2D& get_dayf_logerec_hist( channel chan, double sindecl )
	{
		int b = _pdf_erec_dayfraction_decl[chan].GetZaxis()->FindBin( sindecl);
		return _pdf_erec_dayfraction_decl_slices[chan][b];
	}


	double dP_dOmega_dlogErec_at_declination( channel chan, double logerec, double sindecl )
	{
		 return _dP_dOmega_dlogerec_decl[chan].Interpolate( logerec, sindecl );
	}

	double dN_dOmega_dlogErec_at_declination( channel chan, double logerec, double sindecl )
	{
		 return norm.value * dP_dOmega_dlogErec_at_declination( chan, logerec, sindecl );
	}


	/*! fill the vector _dp_domega_dloge with dP_dOmega_dlogE, which does not depend
	    on any of the Params */

	virtual void init( DataSet& dset )
	{
		_dp_domega_dloge.resize( dset.data.size() );

		for (unsigned i = 0 ; i < dset.data.size(); i++ ) {
			
			auto& sev = dset.data[i];
			if ( sev.index != i ) fatal ("dataset indices not correct");

			_dp_domega_dloge[i] = dP_dOmega_dlogE( (channel) sev.chan, sev.Eproxy, cos( sev.zenith ) );
		}
	}

	virtual void deinit() {
		_dp_domega_dloge.resize(0);
	}

	virtual void update() {
		// nothing to do: only norm can change.
	}


	virtual double dN_dOmega_dlogE(SearchEvent &sev )
	{
		TIMEFUNC	

		double p = 0;

		if (_dp_domega_dloge.size() == 0 ) {
			p = dP_dOmega_dlogE( (channel) sev.chan, sev.Eproxy, cos( sev.zenith ) );
		}
		else p = _dp_domega_dloge[ sev.index ];

		
		double r = norm.value * p;
		if (r<0) fatal( "illegal value", r, sev.str() );
		return r;
	}

	double dP_dOmega_dlogE( channel chan, double E_rec, double cos_zenith_rec )
	{
		return get_hist(chan).Interpolate( E_rec, cos_zenith_rec ) / (2*pi);
	}

	double check_integral( channel chan ) {

		TH2D& h = get_hist(chan);

		double n = 0;

		for ( int bx = 1 ; bx <= h.GetXaxis()->GetNbins() ; bx++ )
		{
			double dloge = h.GetXaxis()->GetBinWidth(bx);

			for ( int by = 1 ; by <= h.GetYaxis()->GetNbins() ; by++ )
			{
				double logerc = h.GetXaxis()->GetBinCenter(bx);
				double coszen = h.GetYaxis()->GetBinCenter(by);
				double sr = h.GetYaxis()->GetBinWidth(by) * 2*pi; // solid angle in this bin

				n += sr * dloge * dP_dOmega_dlogE( chan, logerc, coszen );

			}
		}
		return n;
	}

	double check_integral_declination( channel chan ) {

		TH2D& h = get_hist(chan);

		double n = 0;

		for ( int bx = 1 ; bx <= h.GetXaxis()->GetNbins() ; bx++ )
		{
			double dloge = h.GetXaxis()->GetBinWidth(bx);

			for ( int by = 1 ; by <= h.GetYaxis()->GetNbins() ; by++ )
			{
				double logerc = h.GetXaxis()->GetBinCenter(bx);
				double sindecl = h.GetYaxis()->GetBinCenter(by);
				double sr = h.GetYaxis()->GetBinWidth(by) * 2*pi; // solid angle in this bin
				n += sr * dloge * dP_dOmega_dlogErec_at_declination( chan, logerc, sindecl );
			}
		}
		return n;
	}


	TH2D& get_hist( channel chan ) {

		auto c = hbkg.find( chan );
		if ( c== hbkg.end() ) {
			fatal ("background called for a channel it does not have ", chan );
		}
		return (*c).second;
	}


	virtual SearchEvent generate_event() // nb : only track for the moment
	{
		TIMEFUNC

		double z, E;
		channel chan = track;

		hbkg[chan].GetRandom2( E, z );
		SearchEvent r( acos(z), 2*pi*gRandom->Rndm(), period.random(), E );
		r.compute_coordinates( detres.det() );
		r.chan = chan;
		//dbg ( r.chan , r.__str__() );
		return r;
	}

	/*! version that generates background events in a small region of the 
	 *  celestial sky. */


	vector<SearchEvent> generate_events_region( unsigned howmany, channel chan, EquatorialCoords center, double radius )
	{
		vector<SearchEvent> v(howmany);

		double sindecl = sin ( center.dec() );
		TH2D& h = get_dayf_logerec_hist( chan, sindecl );
		double dayf;
		Mat R( center.asVec() );

		for (unsigned i =0; i<howmany;i++) {
			SearchEvent& r = v[i];

			double x   = sqrt( gRandom->Rndm() ) *radius;
			double phi = 2* pi * gRandom->Rndm();
			Vec v; v.set_angles( x, phi );
		
			r.coordinates.fromVec( R* v );
			r.t = dayfraction_to_time( dayf );
			r.compute_source_direction( detres.det() );

			// lets just check the pdf is not zero
			double p = dN_dOmega_dlogE( r );
			if ( p <= 0 ) fatal(" background generation produced an impossible event ");
		}
		return v;
	}


	void info() {

		Table T ("channel","binsx","binsy","integral");
		T.title = "background-component info";

		for( auto p : hbkg ) {
			auto& h = p.second;
			auto x = p.first;
			T << int(x) << h.GetNbinsX() << h.GetNbinsY() << h.Integral();
		}

		print (T);


	}

	// -----boilerplate-----
	virtual BackgroundComponent *clone() const
	{
		return new BackgroundComponent(*this);
	}

	virtual ~BackgroundComponent() {}
	ClassDef( BackgroundComponent, 1)

};

#endif