#ifndef ISOTROPICPOWERLAWCOMPONENT
#define ISOTROPICPOWERLAWCOMPONENT

#include "Component.hh"


/*! IsotropicPowerLawComponent represents flux of neutrinos that is equal
    from all directions and follows a power law energy distribution. The
    likelihood computation uses the zenith angle, while the acceptance is
    implicitly assumed to be uniform for azimuth. */

struct IsotropicPowerLawComponent : public Component
{

	TH3D Hdens; // dP_dlogErec_dOmega vs logErec and cos_zenith_rec

	Parameter gamma = { "slope of the flux", 2., 0., 10., 0.1};

	virtual double dN_dOmega_dlogE(SearchEvent &sev)
	{
		double E_rec = sev.Eproxy;
		double cos_zenith_rec = cos( sev.zenith );
		return norm.value * dP_dOmega_dlogE( E_rec, cos_zenith_rec , gamma );
	}

	double dP_dOmega_dlogE( double E_rec, double cos_zenith_rec , double gamma )
	{
		return interpolate( 
			&Hdens ,                         
			E_rec, cos_zenith_rec, gamma ,   
			true ) / (2 * pi);                          // clamp to bincenter
	}

	IsotropicPowerLawComponent() {}

	IsotropicPowerLawComponent( string name, DetResponse& dres, int nbins=40, double fmin=1., double fmax=3. ) : Component(name)
	{
		gamma.fit_min = fmin;
		gamma.fit_max = fmax;

		TH1D hz("hz", "dummy for axis", nbins, gamma.fit_min, gamma.fit_max );

		vector<TH2D> stack;

		cout << "Inializing IsotropicPowerLawComponent" << flush;

		for (int b = 1; b <= hz.GetNbinsX() ; b++ )
		{
			double gam = hz.GetBinCenter( b );
		
			// get rate vs true Erec and theta 		
			IsotropicFlux *flux = new IsotropicFlux( "pow(E,"+str(-gam)+")" );
			TH2D h_true = dres.Aeff.GetRate(  flux ); // Nevnts per bin
			check( dres.Aeff.Haeff_costheta );

			normalize( h_true , true );
			h_true.SetName( "h_true"); 
			check( h_true );

			// smear with resolution function to get rate vs reco. quantities
			auto fun = [&]( int a,int b,int c,int d ) { return dres.dP_dlogErec_dOmega_zenith_bins(a,b,c,d);};
			
			TH2D h_rec  = smeared_bin( h_true  , fun );
			h_rec.SetName("h_rec");
			normalize( h_rec , true );
			check( h_rec );
			stack.push_back( h_rec ) ;
			cout << "." << flush;
		}

		Hdens = join( stack, "gamma", hz.GetXaxis()->GetXmin(), hz.GetXaxis()->GetXmax(), false  );
		Hdens.SetTitle("Hdens");
		check( Hdens );
	}

	virtual SearchEvent generate_event()
	{
		static double _gamma = 1e99;
		static TH2D hh;

		if ( _gamma != gamma.value )
		{
			_gamma = gamma.value;
			hh = slice_at_z( Hdens , _gamma );
		}

		static SearchEvent r;
		double z;
		hh.GetRandom2( r.Eproxy, z );
		r.zenith  = acos( z );
		r.azimuth = 2* pi * gRandom->Rndm();
		return r;
	}

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

	virtual void update() 
	{}


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

};


#endif