//____________________________________________________________________________
/*!

\class   genie::flux::GSaeRealPointFlux

\brief   A flux driver for the point-like fluxes

\author  Carla Distefano <distefano_c@lns.infn.it>
         LNS-INFN, Catania

\created February 13, 2015

\cpright  Copyright (c) 2015-2019, The KM3NeT Collaboration
          For the full text of the license see $GSEAGEN/LICENSE

*/
//____________________________________________________________________________


#include <fstream>
#include <cassert>
#include <cstdio>
#include <iostream>
#include <fstream>

#include <TH2D.h>
#include <TMath.h>
#include <TF2.h>

#include "SeaNuDrivers/GSeaRealPointFlux.h"

using namespace std;
using namespace genie;
using namespace genie::flux;


//___________________________________________________________________________
GSeaRealPointFlux::GSeaRealPointFlux(GenParam * GenPar) :
GSeaAtmoFlux(GenPar)
{
  LOG("Flux", pNOTICE)
       << "Instantiating the real point neutrino flux driver";

  this->LoadFluxFile();
}


//___________________________________________________________________________
GSeaRealPointFlux::~GSeaRealPointFlux()
{

}
//___________________________________________________________________________
void GSeaRealPointFlux::LoadFluxFile(void){


 for(unsigned i=0; i<fGenPar->FluxFiles.size(); i++){

    this->SetBinSizes(fGenPar->FluxFiles[i].FluxSimul,fGenPar->EvMin,fGenPar->EvMax);

    string hname=fGenPar->FluxFiles[i].FluxSimul;
    hname.append("_");
    hname.append(PDGLibrary::Instance()->Find(fGenPar->FluxFiles[i].NuType)->GetName());

    if(fGenPar->FluxFiles[i].FluxSimul.find("FUNC")!=string::npos){
      this->FillFuncFlux(i, hname);
    }

  }

  LOG("GSeaRealPointFlux", pNOTICE) << "Neutrino flux simulation data loaded!";

  return;
}

//___________________________________________________________________________
void GSeaRealPointFlux::SetBinSizes(string FluxSim, double Emin, double Emax)
{

  if(FluxSim=="FUNC1") {
    fNumCosThetaBins            = 40;
    fCosThetaMin                = -1.;
    fCosThetaMax                =  1.;
    fMinEv                      = Emin;
    fMaxEv                      = Emax;
    fNumLogEvBinsPerDecade      = 10;
    fDeltaLogEv = 1./fNumLogEvBinsPerDecade;
    fNumLogEvBins               = (int)(TMath::Log10(fMaxEv)-TMath::Log10(fMinEv))*fNumLogEvBinsPerDecade;
  }
  else {
    LOG("Flux", pFATAL) << "Uknonwn flux simulation: " << FluxSim <<" in point-mode generation";
    gAbortingInErr = true;
    exit(1);
  }

  return;
}