#include <time.h>
#include <stdio.h>
#include "rnu.h"

Rnu::Rnu()
{
  //
  //  Default creator
  //
  printf("Creating rnu\n");
  for(Int_t i=0;i<NBIN;i++)
    {
      the_spectrum[i]=0.;
      the_energy[i]=((Double_t)i/10.)+0.05;
    }
}

void Rnu::CrossNuE(Double_t Enu, Double_t ga, Double_t gv, Double_t Munue)
{
  //
  // Load nu-bar - e cross section 
  //
  the_type=1;
  for(Int_t i=0;i<NBIN;i++)
    {
      the_spectrum[i]=dNuEdT(Enu,the_energy[i],ga,gv,Munue);
      printf("i %d the_energy %f the_spectrum %e\n",i,the_energy[i],the_spectrum[i]);
    }

}

void Rnu::NeutrinoSpectrum(Double_t* f)
{
  //
  // Loads neutrino spectrum.  f is a nine element array of weights
  // corresponding to columns in Table 1 of Vogel paper.
  //
  for(Int_t i=0;i<9;i++)
    {
      printf(" %f ",f[i]);
    }
  printf("\n");
  Double_t x[NNUBIN];
  Double_t y[NNUBIN];
  for(Int_t i=0;i<NNUBIN;i++)
    {
      x[i]=GetNeutrinoSpectrumEnergy(i);
      y[i]=0.;
      for(Int_t j=0;j<NUCOMP;j++)
	{
	  y[i]=y[i]+f[j]*GetNeutrinoSpectrumComponent(i,j);
	}
      printf("i %d e %f nu %f\n",i,x[i],y[i]);
    }
  TGraph* g=new TGraph(NNUBIN,x,y);
  TSpline3* s= new TSpline3("x",g);
  for(Int_t i=0;i<NBIN;i++)
    {
      Double_t energy=the_energy[i];
      if(energy<x[NNUBIN-1])the_spectrum[i]=s->Eval(energy);
      if(the_spectrum[i]<0.)the_spectrum[i]=0.;
      printf("energy %f spect %f %f\n",the_energy[i],the_spectrum[i],s->Eval(energy));
    }
}

Double_t Rnu::dNuEdT(Double_t Enu, Double_t T, Double_t ga, Double_t gv, Double_t Munue)
{
  //
  //  Routine to compute nu-bar - e spectrum
  //
  Double_t t1=Gf*Gf*ME*pow(HBARC,2.)/(2.*PI);
  Double_t t21=pow((gv+ga),2.);
  Double_t t22=pow((gv-ga),2.)*pow((1.-T/Enu),2.);
  Double_t t23=(pow(ga,2.)-pow(gv,2.))*ME*T/pow(Enu,2.);
  Double_t t3=0.;
  printf("t1 %e\n",t1);
  printf("t21 %e t22 %e t23 %e\n",t21,t22,t23);
  printf("t3 %e\n",t3);
  if(T!=0.)t3=(PI*pow(ALPHA,2.)*pow(Munue,2.)/pow(ME,2.))*(1-T/Enu)/T;
  Double_t Tmax=ME*(Enu+ME)/sqrt((2.*Enu+ME)*ME);
  printf("Tmax %f\n",Tmax);
  if(T<Tmax)
    {
      return t1*(t21+t22+t23)+t3;
    }else
      {
	return 0.;
      }
}

Double_t Rnu::GetNeutrinoSpectrumEnergy(Int_t i)
{
  //
  // Returns ith energy value for spectrum energy
  //
  Double_t e[NNUBIN]={1.0,1.5,2.0,2.5,3.0,3.5,4.0,4.5,5.0,5.5,6.0,6.5,7.0,7.5,8.0,8.5};
  Double_t out=0.;
  if(i>-1 && i<NNUBIN)
    {
      out=e[i];
    }
  return out;
}

Double_t Rnu::GetNeutrinoSpectrumComponent(Int_t i, Int_t comp)
{
  //
  //  Returns ith value of jth component.
  //
  Double_t U235therm[NNUBIN]={2.44, 1.73, 1.28,0.860,0.580,0.410,0.282,0.181,0.118,0.0720,0.0408,0.0195,0.0104,0.00506,0.00204,0.};
  Double_t U23505[NNUBIN]={2.44,1.72,1.28,0.850,0.590,0.410,0.284,0.184,0.0120,0.0740,0.0423,0.0207,0.0112,0.00556,0.00240,0.};
  Double_t U23514[NNUBIN]={2.17,1.43,1.04,0.680,0.460,0.320,0.213,0.134,0.0860,0.0520,0.0304,0.0148,0.00780,0.00353,0.00121,0.};
  Double_t U23805[NNUBIN]={2.77,1.99,1.51,1.06,0.760,0.550,0.386,0.258,0.171,0.106,0.0605,0.0315,0.0174,0.00928,0.00469,0.};
  Double_t U23814[NNUBIN]={2.59,1.80,1.35,0.940,0.670,0.480,0.337,0.223,0.147,0.0920,0.0545,0.0290,0.0165,0.00896,0.00457,0.};
  Double_t PU239therm[NNUBIN]={2.34,1.50,1.08,0.710,0.477,0.314,0.203,0.121,0.0746,0.0437,0.0237,0.0101,0.00493,0.00220,0.000778,0.};
  Double_t PU23905[NNUBIN]={2.35,1.50,1.09,0.714,0.483,0.319,0.207,0.124,0.0768,0.0452,0.0249,0.0111,0.00554,0.00255,0.000926,0.};
  Double_t PU24005[NNUBIN]={2.47,1.61,1.19,0.798,0.542,0.361,0.236,0.142,0.0882,0.0519,0.0280,0.0123,0.00619,0.00285,0.00111,0.};
  Double_t PU241therm[NNUBIN]={2.64,1.76,1.33,0.913,0.632,0.431,0.288,0.179,0.114,0.0682,0.0371,0.0173,0.00887,0.00449,0.00216,0.};
  Double_t out=-1.;
  if(i>-1 && i<NNUBIN)
    {
      if(comp==0) out=U235therm[i];
      if(comp==1) out=U23505[i];
      if(comp==2) out=U23514[i];
      if(comp==3) out=U23805[i];
      if(comp==4) out=U23814[i];
      if(comp==5) out=PU239therm[i];
      if(comp==6) out=PU23905[i];
      if(comp==7) out=PU24005[i];
      if(comp==8) out=PU241therm[i];
    }else
      {
	printf("GetNeutrinoSpectrumComp - bin out of range, %d\n",i);
      };
  if(out<0.)printf("GetNeutrinoSepectrumComp - component not found, %d\n",comp);
  return out;
}

void Rnu::MakeHist(char* name)
{
  //
  // Makes a TH1D histogram from array.
  //
  printf("GetHist\n");
  printf("%f %f\n",the_energy[0],the_energy[NBIN-1]);
  TH1D* h_x=new TH1D("h_x","",NBIN,the_energy[0],the_energy[NBIN-1]);
  for(Int_t i=0;i<NBIN;i++)
    {
      h_x->SetBinContent(i,the_spectrum[i]);
      //      printf("i %d e %f s %e\n",i,h_x->GetBinCenter(i),h_x->GetBinContent(i));
    }
  TH1D* h=(TH1D*) h_x->Clone();
  h->SetName(name);
  gDirectory->ls();
  h->SetXTitle("T(MeV)");
  h->GetXaxis()->SetTitleOffset(0.);
  if(the_type==1) h->SetYTitle("#sigma (cm^{2})");
}