//____________________________________________________________________________
/*!

\class   genie::flux::GSaeGeometry

\brief   A geometry driver for an under-water neutrino telescope

\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 <cassert>
#include <cmath>
#include <vector>

#include <TSystem.h>

#ifdef _GENIEVERSIONGTEQ3__
#include "Framework/Conventions/Constants.h"
#include "Framework/Conventions/Units.h"
#include "Framework/Messenger/Messenger.h"
#include "Framework/Numerical/RandomGen.h"
#include "Framework/ParticleData/PDGCodeList.h"
#include "Framework/ParticleData/PDGCodes.h"
#include "Framework/ParticleData/PDGUtils.h"
#include "Framework/ParticleData/PDGLibrary.h"
#include "Framework/Utils/PrintUtils.h"
#else
#include "Conventions/Constants.h"
#include "Conventions/Units.h"
#include "Messenger/Messenger.h"
#include "Numerical/RandomGen.h"
#include "PDG/PDGCodeList.h"
#include "PDG/PDGCodes.h"
#include "PDG/PDGUtils.h"
#include "PDG/PDGLibrary.h"
#include "Utils/PrintUtils.h"
#endif

#include "GSeaGeometry.h"

using namespace std;
using namespace genie;
using namespace genie::geometry;
using namespace genie::constants;


//___________________________________________________________________________
GSeaGeometry::GSeaGeometry(GenParam * GenPar)
{

  fGenPar = GenPar;

  fRT = 0.;
  fRL = 0.;
  fRTPlus = 0.;
  fRVol = 0.;
  fHSeaWater = 0.;
  fHRock = 0.;
  fSeaCenter = 0.;
  fRockCenter = 0.;

  if(fGenPar->PropMode) fPEarth = new GPEarth(GenPar);
  else this->DefineMedia();

}

GSeaGeometry::~GSeaGeometry(){

  if(fVacuum){fVacuum=0; delete fVacuum;}
  if(fSeaWater){fSeaWater=0; delete fSeaWater;}
  if(fMixRock){fMixRock=0; delete fMixRock;}

  if(fMatVacuum){fMatVacuum=0; delete fMatVacuum;}
  if(fMixSeawater){fMixSeawater=0; delete fMixSeawater;}
  if(fMixRock){fMixRock=0; delete fMixRock;}

  if(fGeoManager){fGeoManager=0; delete fGeoManager;}

  if(fPEarth){fPEarth=0; delete fPEarth;}

}
//___________________________________________________________________________
bool GSeaGeometry::LoadMuonRange(void)
{
  //Load maximum muon ranges
  double energy,rangesr,rangesw;

  string filename = gSystem->Getenv("GSEAGEN");
  filename.append("/dat/muon_rmax_music.dat");

  if( gSystem->Getenv("GMURNGFL") ) {
    filename = gSystem->Getenv("GMURNGFL");
    SLOG("GSeaGeometry", pINFO) << "$GMURNGFL env.var = " << filename;
  }

  bool is_accessible = ! (gSystem->AccessPathName(filename.c_str()));

  if(is_accessible) {
    SLOG("GSeaGeometry", pINFO) << "Loading max muon ranges from file " << filename;
    ifstream range_stream(filename.c_str(), ios::in);
    int i=0;
    while ( range_stream >> energy >> rangesr >> rangesw ) {
      fRangeSW.SetPoint(i,energy,rangesw);
      fRangeSR.SetPoint(i,energy,rangesr);
      i++;
    }
    range_stream.close();
    return true;
  } else {
    LOG("GSeaGeometry", pFATAL)
       << "Max muon ranges file [" << filename << "] is not accessible!";
    gAbortingInErr = true;
    exit(1);
    return false;
  }


  return false;
}
//___________________________________________________________________________
#ifdef _HETAU_ENABLED__
bool GSeaGeometry::LoadTauRange(void)
{
  //Load maximum tau ranges
  double energy,rangesr,rangesw;

  string filename = gSystem->Getenv("GSEAGEN");
  filename.append("/dat/tau_rmax_music.dat");

  if( gSystem->Getenv("GTAURNGFL") ) {
    filename = gSystem->Getenv("GTAURNGFL");
    SLOG("GSeaGeometry", pINFO) << "$GTAURNGFL env.var = " << filename;
  }

  bool is_accessible = ! (gSystem->AccessPathName(filename.c_str()));

  if(is_accessible) {
    SLOG("GSeaGeometry", pINFO) << "Loading max tau ranges from file " << filename;
    ifstream range_stream(filename.c_str(), ios::in);
    int i=0;
    while ( range_stream >> energy >> rangesr >> rangesw ) {
      fRangeSW.SetPoint(i,energy,rangesw);
      fRangeSR.SetPoint(i,energy,rangesr);
      i++;
    }
    range_stream.close();
    return true;
  } else {
    LOG("GSeaGeometry", pFATAL)
       << "Max tau ranges file [" << filename << "] is not accessible!";
    gAbortingInErr = true;
    exit(1);
    return false;
  }


  return false;
}
#endif
//___________________________________________________________________________
double GSeaGeometry::GetRangeSW(double energy)
{
  double range = pow(10.,fRangeSW.Eval(log10(energy)))/(fGenPar->RhoSW);

  return range;
}
//___________________________________________________________________________
double GSeaGeometry::GetRangeSR(double energy)
{

  double range = pow(10.,fRangeSR.Eval(log10(energy)))/(fGenPar->RhoSR);

  return range;
}
//___________________________________________________________________________
void GSeaGeometry::DefineMedia(void)
{


  fGeoManager =  new TGeoManager("VolGenGeo", "generation volume geometry");

  // define some materials and mixture

  // vacuum
  fMatVacuum = new TGeoMaterial("MatVacuum");
  fMatVacuum->SetA(0.);
  fMatVacuum->SetZ(0.);
  fMatVacuum->SetDensity(0.);

  map<int,double>::iterator iter;

  // seawater
  fMixSeawater = new TGeoMixture("MatSeaWater",fGenPar->SeaWaterComp.size(),fGenPar->RhoSW);
  iter = fGenPar->SeaWaterComp.begin();
  for( ; iter != fGenPar->SeaWaterComp.end(); ++iter)
    fMixSeawater->AddElement(pdg::IonPdgCodeToA(iter->first), pdg::IonPdgCodeToZ(iter->first), iter->second);

  // rock
  fMixRock = new TGeoMixture("MatRock",fGenPar->RockComp.size(),fGenPar->RhoSR);
  iter = fGenPar->RockComp.begin();
  for( ; iter != fGenPar->RockComp.end(); ++iter)
    fMixRock->AddElement(pdg::IonPdgCodeToA(iter->first), pdg::IonPdgCodeToZ(iter->first), iter->second);


  // define the media
  fVacuum = new TGeoMedium("Vacuum", 0, fMatVacuum);
  fSeaWater = new TGeoMedium("SeaWater", 1, fMixSeawater);
  fRock = new TGeoMedium("Rock", 2, fMixRock);

  return;
}

//___________________________________________________________________________
void GSeaGeometry::BuildGeometryMuon(double EneMax, string filename)
{

  if(fRangeSW.GetN()==0){
      LoadMuonRange();
  }


  double RmaxSW = GetRangeSW(EneMax);
  double RmaxSR = GetRangeSR(EneMax);

  if(fGenPar->CtMax>0.)
    fHRock = RmaxSR;		// we are generating up-going events
  else
    fHRock = 0.;		// we are generating only down-going events

  if(fGenPar->CtMin >= 0.)		// we are generating only up-going events
    fHSeaWater = fGenPar->Can2-fGenPar->Can1;
  else
    fHSeaWater = TMath::Min((fGenPar->Can2-fGenPar->Can1)+RmaxSW,fGenPar->SiteDepth);

  fRVol = fGenPar->Can3+RmaxSW;
  fRTPlus = 100.;

  if(!fGenPar->PropMode)this->BuildGeometry(filename);

  return;
}
//___________________________________________________________________________
#ifdef _HETAU_ENABLED__
void GSeaGeometry::BuildGeometryTauTrack(double EneMax, string filename)
{

  if(fRangeSW.GetN()==0){
      LoadTauRange();
  }

  double RmaxSW = GetRangeSW(EneMax);
  double RmaxSR = GetRangeSR(EneMax);

  if(fGenPar->CtMax>0.)
    fHRock = RmaxSR;    // we are generating up-going events
  else
    fHRock = 0.;    // we are generating only down-going events

  if(fGenPar->CtMin >= 0.)    // we are generating only up-going events
    fHSeaWater = fGenPar->Can2-fGenPar->Can1;
  else
    fHSeaWater = TMath::Min((fGenPar->Can2-fGenPar->Can1)+RmaxSW,fGenPar->SiteDepth);

  fRVol = fGenPar->Can3+RmaxSW;
  fRTPlus = 100.;

  if(!fGenPar->PropMode)this->BuildGeometry(filename);

  return;

}
//___________________________________________________________________________
void GSeaGeometry::BuildGeometryTauShower(double EneMax, string filename)
{

  double Rmax = EneMax/kTauMass * 7.*2.91e-13 * kLightSpeed/(units::meter/units::second);

  if(fGenPar->CtMax>0.)
    fHRock = Rmax;    // we are generating up-going events
  else
    fHRock = 0.;    // we are generating only down-going events

  if(fGenPar->CtMin >= 0.)    // we are generating only up-going events
    fHSeaWater = fGenPar->Can2-fGenPar->Can1;
  else
    fHSeaWater = TMath::Min((fGenPar->Can2-fGenPar->Can1)+Rmax,fGenPar->SiteDepth);

  fRVol = fGenPar->Can3+Rmax;
  fRTPlus = 100.;

  if(!fGenPar->PropMode)this->BuildGeometry(filename);

  return;

}
#endif
//___________________________________________________________________________
void GSeaGeometry::BuildGeometryElectron(string filename)
{

  // the interaction volume is the detector "can"
  fHRock = fGenPar->HBedRock;	
  	
  fHSeaWater = fGenPar->Can2-fGenPar->Can1;
  fRVol = fGenPar->Can3;

  if(!fGenPar->PropMode)this->BuildGeometry(filename);

  return;
}
//___________________________________________________________________________
void GSeaGeometry::BuildGeometry(string filename)
{

  //  build Earth as ROOT geometry: lenght in cm, density in g/cm^3

  double top_half_height   = (fHSeaWater+fHRock)/2.;
  double water_half_height = fHSeaWater/2.;
  double rock_half_height  = fHRock/2.;

  fRockCenter  = -rock_half_height+fGenPar->Can1;
  fSeaCenter   = water_half_height+fGenPar->Can1;

  DefineMedia();

  // define the volumes
  // top volume

  double RMax = fRVol;
  double HalfH = top_half_height;

  if (fGenPar->Can1==0 && fGenPar->Can2==0){

    TGeoVolume* TopVolume = fGeoManager->MakeSphere("TopVolume",fVacuum,0.,RMax*1.E2);
    fGeoManager->SetTopVolume(TopVolume);
    TGeoVolume* WaterVol = fGeoManager->MakeSphere("WaterVol",fSeaWater,0.,fRVol*1.E2);
    WaterVol->SetFillColor(kBlue);
    WaterVol->SetLineColor(kBlue);
    WaterVol->SetLineWidth(2);
    WaterVol->SetLineStyle(2);
    TopVolume->AddNode(WaterVol,1);
    
    fGeoManager->CloseGeometry();
    cout << "exporting geometry" << endl;
    
    fGeoManager->Export(filename.c_str());

  }else{
  
    TGeoVolume *TopVolume = fGeoManager->MakeTube("TopVolume",fVacuum,0.,RMax*1.E2,2.*HalfH*1.E2);


    fGeoManager->SetTopVolume(TopVolume);

    // define the translations
    TGeoTranslation * TransRock = new TGeoTranslation(0.,0.,fRockCenter*1.E2);
    TGeoTranslation * TransWater = new TGeoTranslation(0.,0.,fSeaCenter*1.E2);


    // rock volume
    if(fHRock>0.){
      TGeoVolume *RockVol = fGeoManager->MakeTube("RockVol",fRock,0.,fRVol*1.E2,rock_half_height*1.E2);
      RockVol->SetFillColor(kOrange+3);
      RockVol->SetLineColor(kOrange+3);
      TopVolume->AddNode(RockVol,2,TransRock);
    }

    // seawater volume
    if(fHSeaWater>0.){
      TGeoVolume *WaterVol = fGeoManager->MakeTube("WaterVol",fSeaWater,0.,fRVol*1.E2,water_half_height*1.E2);
      WaterVol->SetFillColor(kBlue);
      WaterVol->SetLineColor(kBlue);
      TopVolume->AddNode(WaterVol,1,TransWater);
    }

    fGeoManager->CloseGeometry();
    fGeoManager->Export(filename.c_str());
  }
  return;
}

//_____________________________________________________________________________
void GSeaGeometry::CalcRadii(void)
{

  fRL = TMath::Power((fHSeaWater+fHRock),2)+TMath::Power(2.*fRVol,2);
  fRL = TMath::Sqrt(fRL);

  if(fGenPar->RTOpt.compare("vol")==0){
    fX0=0.;
    fY0=0.;
    fZ0= fHSeaWater+fGenPar->Can1-(fHSeaWater+fHRock)/2.;
    fRT = TMath::Sqrt(TMath::Power(2.*fRVol,2)+TMath::Power((fHSeaWater+fHRock),2))/2.;
    fRT += fRTPlus;
  }
  else if(fGenPar->RTOpt.compare("proj")==0){
    fX0=0.;
    fY0=0.;
    fZ0= fGenPar->Can2-(fGenPar->Can2-fGenPar->Can1+fGenPar->HBedRock)/2.;
    fRT = -fRTPlus;
  }
  else if(fGenPar->RTOpt.compare("can")==0 || fGenPar->PropMode){
    fX0=0.;
    fY0=0.;
    fZ0= fGenPar->Can2-(fGenPar->Can2-fGenPar->Can1+fGenPar->HBedRock)/2.;
    fRT = TMath::Sqrt(TMath::Power(2.*fGenPar->Can3,2)+TMath::Power((fGenPar->Can2-fGenPar->Can1+fGenPar->HBedRock),2))/2.;
    fRT += fRTPlus;
  }
  else {
    fX0=0.;
    fY0=0.;
    fZ0= fGenPar->Can2-(fGenPar->Can2-fGenPar->Can1+fGenPar->HBedRock)/2.;

    fRT=atof(fGenPar->RTOpt.c_str());
  }

  if(fGenPar->PropMode)fRL = 2*fPEarth->GetREarth()*1E3;

  return;

}