#include "PointSourceComponent.hh"



TH3D& PointSourceComponent::_pdf_erec_dayfraction_loga ( flavor flav, channel chan )
{
	static TTimeStamp T0 (2020, 1, 1, 0, 0, 0 );

	{	
	TH3D& r = HFC[flav][chan];
	if ( r.GetEntries() > 0 ) return r; 
    }
    
	print(name, "building lookup histogram.", flav, chan );

	string nam = "HDR" + flavor_name(flav) + channel_name(chan);
	string tit = "pdf of erec, dayfraction, loga for " + flavor_name(flav) + channel_name(chan);


	// we use h_ANGERR to derive the binnin of the nes histogram
	DefaultDetResponse& detres = detector_response.get( flav, chan );

	TAxis* XX = detres.h_ANGERR.GetXaxis();
	TAxis* ZZ = detres.h_ANGERR.GetZaxis();

	HFC[flav][chan] = TH3D( nam.c_str(), tit.c_str(),
							XX->GetNbins(), XX->GetXmin(), XX->GetXmax(),
							48            , 0            , 1            , // dayfraction
							ZZ->GetNbins(), ZZ->GetXmin(), ZZ->GetXmax() );

	TH3D& H = HFC[flav][chan];
	H.SetXTitle("Erec (GeV)");
	H.SetYTitle("Dayfraction");	
	H.SetZTitle("alpha (rad)");	

	


	// fill it.

	TTimeStamp T;

	


	for ( int bx = 1; bx <= H.GetNbinsX() ; bx++ )          // Erec slice
	{
		print ( " erec slice : ", bx, "/", H.GetNbinsX() );
		for ( int by = 1; by <= H.GetNbinsY() ; by++ )        // time of day from 0 to 1
		{
			double erec = H.GetXaxis()->GetBinCenter(bx);
			double dayf = H.GetYaxis()->GetBinCenter(by);

			// compute the true zenith angle for this time of day
			T.SetSec( T0.GetSec() + dayf * siderial_day() ) ;             // todo : use period
			EquatorialCoords eq( 0 , dec.value );                    // nb: ra set to zero on purpose
			Vec dtrue = -eq.track_direction( detres.det , T );
			double theta_true = dtrue.theta();

			if ( bx == 1 ) {
				print ("DAYF", dayf, theta_true );
			}


			for ( int bz = 1; bz <= H.GetNbinsZ() ; bz++ )         // log angular error
			{
				double loga = H.GetZaxis()->GetBinCenter(bz);

				double J = detres.psf.d_loga_d_omega ( loga );

				double alpha = pow(10, loga) * pi / 180 ;

				double v;
				if ( alpha > pi ) v = 0;
				else v = dN_dOmega_dlogErec_alpha( track, alpha , erec, theta_true , false ) / J ;

				if ( v < 0 )
				{
					print ("warning : negative value", v, "for", flav, chan , "in N_dOmega_dlogErec_alpha at", bz, "apha=" , alpha, "erec=",erec, "costhetatrue=",cos(theta_true) , "J=", J );
					v = 0;
				}

				int b = H.GetBin(bx, by, bz);
				H.SetBinContent( b, v );
			}
		}
	}


	return H;
}




Table PointSourceComponent::rate_info( flavor flav, channel chan ) 
{

	Table T("Ebin","logE","Flux (/GeV m2 s)","Aeff@decl (m)","Rate (/yr)");

	auto v = detector_response.Aeff(flav, chan ).GetAreaAtDecl( dec.value ).GetRateInfo( *flux );

	double c =0;
	for(unsigned i = 0; i<v[0].size() ; i++ ) 
	{
	  T << i << v[0][i] << v[1][i] << v[2][i] << v[3][i];
	  c += v[3][i];
	}

	T << "total" << "-" << "-" << "-" << c ;
	T << "xcheck" << "-" << "-" << "-" << Ntot(flav,chan);
	return T;
}


Table PointSourceComponent::rate_info() 
{
	Table R;

	for ( auto f : flavors ) {
	  for ( auto c : channels ) {
	      string n = flavor_name( f ) +"/"+ channel_name( c );

	      auto T = rate_info( f, c );

	      Table_row top;
	      top.resize( T.ncols() );
	      top[0] = n;
	      T.insert(T.begin(), top );

	      if ( R.ncols() == 0 ) 
	        {
	          R = T;
	          continue;
	        }

	      for ( unsigned i =0; i< R.size(); i++ ) {
	        R[i].insert( R[i].end(), T[i].begin(), T[i].end() );
	      }

	  }
	}

	return R;
}