//____________________________________________________________________________
/*!
\class genie::flux::GSaeRealAtmoFlux
\brief A flux driver for the Atmospheric and Diffuse Neutrino Fluxes
\author Carla Distefano <distefano_c@lns.infn.it>
LNS-INFN, Catania
\created September 13, 2012
\cpright Copyright (c) 2012-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 <TH3D.h>
#include <TMath.h>
#include <TF1.h>
#include <TF3.h>
#include <TFile.h>
#ifdef _GENIEVERSIONGTEQ3__
#include "Framework/Conventions/Constants.h"
#include "Framework/Messenger/Messenger.h"
#include "Framework/ParticleData/PDGLibrary.h"
#else
#include "Conventions/Constants.h"
#include "Messenger/Messenger.h"
#include "PDG/PDGLibrary.h"
#endif
#include "SeaNuDrivers/GSeaRealAtmoFlux.h"
using namespace std;
using namespace genie;
using namespace genie::flux;
using namespace genie::constants;
//___________________________________________________________________________________________________
double GaisserFluxParam(double *x, double *param) {
double gamma=2.7;
double eps_pi=115;
double eps_K=850;
double Phi0=param[0];
double Apinu=param[1];
double Aknu=param[2];
double Bpinu=param[3];
double Bknu=param[4];
double energy=x[0];
double costheta=fabs(x[1]);
double PhiN=Phi0*pow(energy/1E3,-gamma);
double frac1=Apinu/(1.+ Bpinu*energy*costheta/eps_pi);
double frac2= Aknu/(1.+ Bknu*energy*costheta/eps_K);
double flux=PhiN*(frac1+frac2);
return flux;
}
//___________________________________________________________________________
GSeaRealAtmoFlux::GSeaRealAtmoFlux(GenParam * GenPar) :
GSeaAtmoFlux(GenPar)
{
LOG("GSeaRealAtmoFlux", pNOTICE)
<< "Instantiating the real atmospheric neutrino flux driver";
this->LoadFluxFile();
}
//___________________________________________________________________________
GSeaRealAtmoFlux::~GSeaRealAtmoFlux()
{
}
//___________________________________________________________________________
void GSeaRealAtmoFlux::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);
}
else if(fGenPar->FluxFiles[i].FluxSimul == "BARTOL" || fGenPar->FluxFiles[i].FluxSimul == "FLUKA"){
TH2D * histo = this->CreateFluxHisto2D(hname.c_str(),hname.c_str());
for(unsigned j=0; j<fGenPar->FluxFiles[i].FileNames.size(); j++)
this->FillFluxHisto2D(fGenPar->FluxFiles[i].FluxSimul,histo, fGenPar->FluxFiles[i].FileNames[j]);
FluxHisto2D flux;
flux.Flux2D=histo;
flux.Range[0]=fMinEv;
flux.Range[1]=fMaxEv;
flux.Range[2]=fCosThetaMin;
flux.Range[3]=fCosThetaMax;
map<int,vector<FluxHisto2D> >::iterator MapEntry = fFlux2DRange.find(fGenPar->FluxFiles[i].NuType);
MapEntry->second.push_back(flux);
if(fGenPar->FluxFiles[i].FluxSimul == "BARTOL"){ // no extrapolation for fluka at the moment
if(fGenPar->EvMax>fMaxEv){
FluxFunc2D fluxhe=this->ExtrapolateFlux(histo);
map<int,vector<FluxFunc2D> >::iterator MapEntryHe = fFlux2FRange.find(fGenPar->FluxFiles[i].NuType);
MapEntryHe->second.push_back(fluxhe);
fMaxEv=fGenPar->EvMax;
}
}
}
else if(fGenPar->FluxFiles[i].FluxSimul == "HONDA"){
TH2D * histo2D = this->CreateFluxHisto2D(hname.c_str(),hname.c_str());
TH3D * histo3D = this->CreateFluxHisto3D(hname.c_str(),hname.c_str());
// for(unsigned j=0; j<fGenPar->FluxFiles[i].FileNames.size(); j++)
// this->FillFluxHisto2D(fGenPar->FluxFiles[i].FluxSimul,histo, fGenPar->FluxFiles[i].FileNames[j]);
bool Only1PhiBin = this->FillHondaFlux(fGenPar->FluxFiles[i].NuType, histo2D, histo3D, fGenPar->FluxFiles[i].FileNames[0]); // only 1 file per nu type!!!
if(Only1PhiBin){
FluxHisto2D flux;
flux.Flux2D=histo2D;
flux.Range[0]=fMinEv;
flux.Range[1]=fMaxEv;
flux.Range[2]=fCosThetaMin;
flux.Range[3]=fCosThetaMax;
map<int,vector<FluxHisto2D> >::iterator MapEntry = fFlux2DRange.find(fGenPar->FluxFiles[i].NuType);
MapEntry->second.push_back(flux);
}
else {
FluxHisto3D flux;
flux.Flux3D=histo3D;
flux.Range[0]=fMinEv;
flux.Range[1]=fMaxEv;
flux.Range[2]=fCosThetaMin;
flux.Range[3]=fCosThetaMax;
flux.Range[4]=fPhiMin;
flux.Range[5]=fPhiMax;
map<int,vector<FluxHisto3D> >::iterator MapEntry = fFlux3DRange.find(fGenPar->FluxFiles[i].NuType);
MapEntry->second.push_back(flux);
}
if(fGenPar->EvMax>fMaxEv){
FluxFunc2D flux=this->ExtrapolateFlux(histo2D);
map<int,vector<FluxFunc2D> >::iterator MapEntry = fFlux2FRange.find(fGenPar->FluxFiles[i].NuType);
MapEntry->second.push_back(flux);
fMaxEv=fGenPar->EvMax;
}
if(Only1PhiBin)delete histo3D;
else delete histo2D;
}
else {
LOG("GSeaRealAtmoFlux", pFATAL) << "Uknonwn flux simulation: " << fGenPar->FluxFiles[i].FluxSimul;
gAbortingInErr = true;
exit(1);
}
}
// in case of atmospheric neutrinos we have to check is an extrapolation to highe energies is needed
LOG("GSeaAtmoFlux", pNOTICE) << "Neutrino flux simulation data loaded!";
//exit(1);
return;
}
//___________________________________________________________________________
void GSeaRealAtmoFlux::SetBinSizes(string FluxSim, double Emin, double Emax)
{
if(FluxSim == "BARTOL") {// BARTOL FLUX
fNumCosThetaBins = 20;
fCosThetaMin = -1.;
fCosThetaMax = 1.;
fNumPhiBins = 1;
fPhiMin = 0.;
fPhiMax = 2.*kPi;
fNumLogEvBins = 50;
fNumLogEvBinsPerDecade = 10;
fMinEv = 0.1; // GeV
fDeltaLogEv = 1./fNumLogEvBinsPerDecade;
fMaxEv = TMath::Power(10.,TMath::Log10(fMinEv)+fDeltaLogEv*fNumLogEvBins);
}
else if(FluxSim == "FLUKA") {// FLUKA3D FLUX
fNumCosThetaBins = 40;
fCosThetaMin = -1.;
fCosThetaMax = 1.;
fNumPhiBins = 1;
fPhiMin = 0.;
fPhiMax = 2.*kPi;
fNumLogEvBins = 61;
fNumLogEvBinsPerDecade = 20;
fMinEv = 0.1; // GeV
fDeltaLogEv = 1./fNumLogEvBinsPerDecade;
fMaxEv = TMath::Power(10.,TMath::Log10(fMinEv)+fDeltaLogEv*fNumLogEvBins);
}
else if(FluxSim == "HONDA"){// HONDA FLUX
fNumCosThetaBins = 20;
fCosThetaMin = -1.;
fCosThetaMax = 1.;
fNumPhiBins = 12;
fPhiMin = 0.;
fPhiMax = 2.*kPi;
fNumLogEvBins = 101;
fNumLogEvBinsPerDecade = 20;
fDeltaLogEv = 1./fNumLogEvBinsPerDecade;
fMinEv = TMath::Power(10.,TMath::Log10(0.1)-fDeltaLogEv/2.); // GeV
fMaxEv = TMath::Power(10.,TMath::Log10(fMinEv)+fDeltaLogEv*fNumLogEvBins);
}
else { // no histos created or histo read from files, we need just fMinEv and fMaxEv for GENIE
fNumCosThetaBins = 1;
fCosThetaMin = -1.;
fCosThetaMax = 1.;
fNumPhiBins = 1;
fPhiMin = 0.;
fPhiMax = 2.*kPi;
fMinEv = Emin;
fMaxEv = Emax;
fNumLogEvBinsPerDecade = 1;
fDeltaLogEv = 1./fNumLogEvBinsPerDecade;
fNumLogEvBins = (int)(TMath::Log10(fMaxEv)-TMath::Log10(fMinEv))*fNumLogEvBinsPerDecade;
}
return;
}
//___________________________________________________________________________
bool GSeaRealAtmoFlux::FillFluxHisto2D(string FluxSim, TH2D * histo, string fluxname)
{
if(FluxSim == "BARTOL")
this->FillBartolFlux(histo, fluxname);
else if (FluxSim == "FLUKA")
this->FillFlukaFlux(histo, fluxname);
else
return false;
return true;
}
//___________________________________________________________________________
bool GSeaRealAtmoFlux::FillBartolFlux(TH2D * histo, string filename)
{
LOG("GSeaRealAtmoFlux", pNOTICE) << "Loading: " << filename;
if(!histo) {
LOG("GSeaRealAtmoFlux", pERROR) << "Null flux histogram!";
return false;
}
ifstream flux_stream(filename.c_str(), ios::in);
if(!flux_stream) {
LOG("GSeaRealAtmoFlux", pERROR) << "Could not open file: " << filename;
return false;
}
double energy, costheta, flux;
double junkd; // throw away error estimates
// throw away comment line
flux_stream.ignore(99999, '\n');
while ( !flux_stream.eof() ) {
flux = 0.0;
flux_stream >> energy >> costheta >> flux >> junkd >> junkd;
if( flux>0.0 ){
// Compensate for logarithmic units - dlogE=dE/E
// [Note: should do this explicitly using bin widths]
flux /= energy;
LOG("GSeaRealAtmoFlux", pDEBUG)
<< "Flux[Ev = " << energy
<< ", cosTheta = " << costheta << "] = " << flux;
histo->Fill(energy,costheta,flux);
}
}
return true;
}
//___________________________________________________________________________
bool GSeaRealAtmoFlux::FillFlukaFlux(TH2D * histo, string filename)
{
LOG("GSeaRealAtmoFlux", pNOTICE) << "Loading: " << filename;
if(!histo) {
LOG("GSeaRealAtmoFlux", pERROR) << "Null flux histogram!";
return false;
}
ifstream flux_stream(filename.c_str(), ios::in);
if(!flux_stream) {
LOG("GSeaRealAtmoFlux", pERROR) << "Could not open file: " << filename;
return false;
}
double energy, costheta, flux;
char j1, j2;
while ( !flux_stream.eof() ) {
flux = 0.0;
flux_stream >> energy >> j1 >> costheta >> j2 >> flux;
if( flux>0.0 ){
LOG("GSeaRealAtmoFlux", pDEBUG)
<< "Flux[Ev = " << energy
<< ", cosTheta = " << costheta << "] = " << flux;
// note: reversing the Fluka sign convention for zenith angle
histo->Fill(energy,-costheta,flux);
}
}
return true;
}
//___________________________________________________________________________
bool GSeaRealAtmoFlux::FillHondaFlux(int NuType, TH2D * histo2D, TH3D * histo3D, string filename)
{
LOG("GSeaRealAtmoFlux", pNOTICE) << "Loading: " << filename;
ifstream file(filename.c_str());
string line;
// check if Azimuth angle dependent or Azimuth angle averaged flux Tables.
bool Only1PhiBin=true;
float cosZmin,cosZmax;
float PhiMin, PhiMax, PhiMin0=-1;
int NumPhiBins=0;
while(getline(file,line)){
size_t average = line.find("average");
if (average!=string::npos){
char c[10];
size_t found1 = line.rfind("=");
size_t found2 = line.rfind("]");
string str = line.substr(found1+1, found2-found1-1);
sscanf (str.c_str(), "%f %s %f ",&PhiMin, c,&PhiMax);
if(PhiMin>PhiMin0){
NumPhiBins++;
PhiMin0=PhiMin;
}
}
}
if(NumPhiBins==1)LOG("GSeaRealAtmoFlux", pNOTICE) << "HONDA flux: Azimuth angle averaged flux file";
else if(NumPhiBins==12){
Only1PhiBin=false;
LOG("GSeaRealAtmoFlux", pNOTICE) << "HONDA flux: Azimuth angle dependent flux file";
}
else {
LOG("GSeaRealAtmoFlux", pFATAL) << "HONDA flux: No standard file";
gAbortingInErr = true;
exit(1);
}
file.clear();
file.seekg(0, ios::beg);
float PhiMean=0.,CosZMean=0.;
float energy=0.,flux=0.,NuMuFlux,NuMubarFlux,NuEFlux,NuEbarFlux;
double PhiCenter[histo3D->GetZaxis()->GetNbins()];
histo3D->GetZaxis()->GetCenter(PhiCenter);
while(getline(file,line)){
size_t average = line.find("average");
if (average!=string::npos){
char c[10];
size_t found1 = line.find("=");
size_t found2 = line.find(",");
string str = line.substr(found1+1, found2-found1-1);
sscanf (str.c_str(), "%f %s %f ",&cosZmin, c,&cosZmax);
found1 = line.rfind("=");
found2 = line.rfind("]");
str = line.substr(found1+1, found2-found1-1);
sscanf (str.c_str(), "%f %s %f ",&PhiMin, c,&PhiMax);
CosZMean=(cosZmin+cosZmax)/2.;
PhiMean=(PhiMax+PhiMin+180.)/2.*kPi/180.; // PhiMax+PhiMin+180. --> phi angles are measured counter-clockwise from the south and not from north.
PhiMean=fmod(PhiMean+2*kPi,2*kPi); // normalization between 0 e 180 deg
}
else {
if(line.find("Enu")==string::npos){
sscanf (line.c_str(),"%f %f %f %f %f",&energy,&NuMuFlux,&NuMubarFlux,&NuEFlux,&NuEbarFlux);
if(NuType==14)flux=NuMuFlux;
else if(NuType==-14)flux=NuMubarFlux;
else if(NuType==12)flux=NuEFlux;
else if(NuType==-12)flux=NuEbarFlux;
if(Only1PhiBin){
for(int i=0; i<histo3D->GetZaxis()->GetNbins(); i++)
histo3D->Fill(energy,CosZMean,PhiCenter[i],flux);
histo2D->Fill(energy,CosZMean,flux);
}
else {
histo3D->Fill(energy,CosZMean,PhiMean,flux);
histo2D->Fill(energy,CosZMean,flux/(double)NumPhiBins);
}
}
}
}
file.close();
return Only1PhiBin;
}
//___________________________________________________________________________
FluxFunc2D GSeaRealAtmoFlux::ExtrapolateFlux(TH2D * histo2D){
const int npar = 5;
double GaisserParam[npar];
GaisserParam[0]=1e-4;
GaisserParam[1]=8.9E-3;
GaisserParam[2]=2.6E-3;
GaisserParam[3]=2.77;
GaisserParam[4]=1.18;
string name=histo2D->GetName();
name.append("_fit");
// ci facciamo restituire direttamente la funzione di fit
TF2 * f2 = new TF2(name.c_str(),GaisserFluxParam,1.E3,1.E8,-1.,1.,npar);
f2->SetParameters(GaisserParam);
f2->SetParLimits(0,0.5*1E-4,1.5*1E-4);
f2->SetParLimits(1,8.5E-3,10.5E-3);
f2->SetParLimits(2,2.E-3,3E-3);
f2->SetParLimits(3,2.5,3);
f2->SetParLimits(4,1.,2);
histo2D->Fit(name.c_str(),"LRQ","",1.E3,1.E4);
name.append(".root");
TFile OutFile( name.c_str(), "recreate");
histo2D->Write();
OutFile.Close();
FluxFunc2D flux;
flux.Flux2D=f2;
flux.Range[0]=histo2D->GetXaxis()->GetXmax();
flux.Range[1]=fGenPar->EvMax;
flux.Range[2]=fCosThetaMin;
flux.Range[3]=fCosThetaMax;
return flux;
}