//____________________________________________________________________________
/*!

\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)) + fGenPar->MuonRangeTolerance;

  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->CtMin >= 0.) 	// we are generating only up-going events	
    fZmaxVol = fGenPar->Can2;
  else 
    fZmaxVol = min(fGenPar->Can2 + RmaxSW, fGenPar->SiteDepth);

  if(fGenPar->CtMax <= 0.){ 	// we are generating only down-going events
    fZminVol = fGenPar->Can1;
    fHSeaWater = fZmaxVol - fZminVol;
    fHRock = 0.;
  }  
  else if (fGenPar->Can1>=RmaxSW){ // the distance beetwen the can bottom and the seabed is larger than the range
    fZminVol = fGenPar->Can1 - RmaxSW;
    fHSeaWater = fZmaxVol - fZminVol;    
    fHRock = 0.;
  }
  else {
    fHSeaWater = fZmaxVol;
    fHRock = RmaxSR - fGenPar->Can1*fGenPar->RhoSW/fGenPar->RhoSR;  // we calculate the height of rock taking into account the layer of water below the can
    fZminVol = - fHRock;  
  }

  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->CtMin >= 0.) 	// we are generating only up-going events	
    fZmaxVol = fGenPar->Can2;
  else 
    fZmaxVol = min(fGenPar->Can2 + RmaxSW, fGenPar->SiteDepth);

  if(fGenPar->CtMax <= 0.){ 	// we are generating only down-going events
    fZminVol = fGenPar->Can1;
    fHSeaWater = fZmaxVol - fZminVol;
    fHRock = 0.;
  }  
  else if (fGenPar->Can1>=RmaxSW){ // the distance beetwen the can bottom and the seabed is larger than the range
    fZminVol = fGenPar->Can1 - RmaxSW;
    fHSeaWater = fZmaxVol - fZminVol;    
    fHRock = 0.;
  }
  else {
    fHSeaWater = fZmaxVol;
    fHRock = RmaxSR - fGenPar->Can1*fGenPar->RhoSW/fGenPar->RhoSR;  // we calculate the height of rock taking into account the layer of water below the can
    fZminVol = - fHRock;  
  }

  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->CtMin >= 0.) 	// we are generating only up-going events	
    fZmaxVol = fGenPar->Can2;
  else 
    fZmaxVol = min(fGenPar->Can2 + Rmax, fGenPar->SiteDepth);

  if(fGenPar->CtMax <= 0.){ 	// we are generating only down-going events
    fZminVol = fGenPar->Can1;
    fHSeaWater = fZmaxVol - fZminVol;
    fHRock = 0.;
  }  
  else if (fGenPar->Can1>=Rmax){ // the distance beetwen the can bottom and the seabed is larger than the range
    fZminVol = fGenPar->Can1 - Rmax;
    fHSeaWater = fZmaxVol - fZminVol;    
    fHRock = 0.;
  }
  else {
    fHSeaWater = fZmaxVol;
    fHRock = Rmax - fGenPar->Can1*fGenPar->RhoSW/fGenPar->RhoSR;  // we calculate the height of rock taking into account the layer of water below the can
    fZminVol = - fHRock;  
  }

  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"
// only if the can bottom is on the seabed we add the bedrock
  
  fHSeaWater = fGenPar->Can2-fGenPar->Can1;  
  fZmaxVol = fGenPar->Can2;  
  
  if(fGenPar->Can1>0){
    fZminVol = fGenPar->Can1;
    fHRock = 0.;
  }
  else {
    fZminVol = - fGenPar->HBedRock;
    fHRock = fGenPar->HBedRock; 
  }  	

  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   = max(fZmaxVol, fabs(fZminVol));
  double water_half_height = fHSeaWater*0.5;
  double rock_half_height  = fHRock*0.5;

  fRockCenter  = -rock_half_height;
  fSeaCenter   = fZmaxVol - water_half_height;

  DefineMedia();

  // define the volumes
  // top volume

  double RMax = fRVol*1.2;
  double HalfH = top_half_height*1.2;

  if (fGenPar->Can1==0 && fGenPar->Can2==0){ // to be checked

    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);
/*
    TGeoTranslation * TransCan = new TGeoTranslation(0.,0.,(fGenPar->Can2+fGenPar->Can1)/2.*1.E2-water_half_height*1.E2);
    TGeoVolume *CanVol = fGeoManager->MakeTube("CanVol",fSeaWater,0.,fGenPar->Can3*1.E2,(fGenPar->Can2-fGenPar->Can1)/2.*1.E2);
    CanVol->SetFillColor(kGreen);
    CanVol->SetLineColor(kGreen);
    WaterVol->AddNode(CanVol,3,TransCan);


*/
    }
/*


    //for visualization

    TGeoTranslation * TransCan = new TGeoTranslation(0.,0.,(fGenPar->Can2+fGenPar->Can1)/2.*1.E2);
    TGeoVolume *CanVol = fGeoManager->MakeTube("CanVol",fSeaWater,0.,fGenPar->Can3*1.E2,(fGenPar->Can2-fGenPar->Can1)/2.*1.E2);
    CanVol->SetFillColor(kGreen);
    CanVol->SetLineColor(kGreen);
    TopVolume->AddNodeOverlap(CanVol,3,TransCan);

*/


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

//_____________________________________________________________________________
void GSeaGeometry::CalcRadii(void)
{
  double fHSeaWpRock = fHSeaWater+fHRock;
  double sqr2fRVol = 4.* fRVol*fRVol;

  double HBedRock;	
  if(fGenPar->Can1>0) HBedRock=0.;	
  else HBedRock = fGenPar->HBedRock;

  fRL = sqrt(fHSeaWpRock*fHSeaWpRock + sqr2fRVol);

  fX0 = 0.;
  fY0 = 0.;

  if(fGenPar->RTOpt.compare("vol")==0){
    fZ0 = (fZminVol+fZmaxVol)*0.5;
    fRT = fRL*0.5;
    fRT += fRTPlus;
  }
  else if(fGenPar->RTOpt.compare("proj")==0){
    fZ0= fGenPar->Can2-(fGenPar->Can2-fGenPar->Can1+HBedRock)*0.5;
    fRT = -fRTPlus;
  }
  else if(fGenPar->RTOpt.compare("can")==0 || fGenPar->PropMode){
    fZ0= fGenPar->Can2-(fGenPar->Can2-fGenPar->Can1+HBedRock)*0.5;
    double Can2m1pBedrock= (fGenPar->Can2-fGenPar->Can1+HBedRock);
    fRT = sqrt(4.*(fGenPar->Can3)*(fGenPar->Can3)+Can2m1pBedrock*Can2m1pBedrock)*0.5;
    fRT += fRTPlus;
  }
  else {
    fZ0= fGenPar->Can2-(fGenPar->Can2-fGenPar->Can1+HBedRock)*0.5;
    fRT=atof(fGenPar->RTOpt.c_str());
  }

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

  return;

}