#include <G4Material.hh>
#include <G4SystemOfUnits.hh>

#include <RAT/Optics.hh>
#include <RAT/DB.hh>
#include <RAT/Log.hh>
using namespace RAT;

#include <CLHEP/Units/SystemOfUnits.h>
#include <CLHEP/Units/PhysicalConstants.h>
#include <sstream>
using namespace CLHEP;
using namespace std;

void
Optics::LoadOptics()
{
  DBLinkGroup optics = DB::Get()->GetLinkGroup( "OPTICS" );
  DBLinkGroup::iterator iMaterialOptics;
  for( iMaterialOptics = optics.begin(); iMaterialOptics != optics.end(); ++iMaterialOptics )
    {
      const string name = iMaterialOptics->first;
      DBLinkPtr dbTable = iMaterialOptics->second;
      Log::Assert( G4Material::GetMaterial( name ) != NULL, "Optics::LoadOptics: Optical properties defined for non-existant material." );
      LoadOpticalProperties( name, dbTable );
    }
}


void
Optics::LoadOpticalProperties( const std::string& name,
                               DBLinkPtr dbTable )
{
  vector<string> properties;
  try { properties = dbTable->GetSArray( "PROPERTY_LIST" ); }
  catch( DBNotFoundError& e ) { Log::Die( "Optics::LoadOpticalProperties: No PROPERTY_LIST for material " + name ); }
  G4Material* material = G4Material::GetMaterial( name );
  G4MaterialPropertiesTable* propertiesTable = material->GetMaterialPropertiesTable();
  // Now loop over and load in the properties
  bool absScaling = false, rsScaling = false, mieScaling = false;
  for( vector<string>::const_iterator iProperty = properties.begin(); iProperty != properties.end(); ++iProperty )
    {
      if( *iProperty == string( "LIGHT_YIELD" ) || iProperty->find( "REEMISSION_PROB" ) != string::npos ||
          *iProperty == string( "MIEHG_FORWARD" ) || *iProperty == string( "MIEHG_BACKWARD" ) ||
          *iProperty == string( "MIEHG_FORWARD_RATIO" ) || *iProperty == string( "RINDEX_MEAN" ) ||
          *iProperty == string( "RS_SCALE_FACTOR" ) )
        propertiesTable->AddConstProperty( iProperty->c_str(), dbTable->GetD( *iProperty ) );
      else if( *iProperty == string( "ISOTHERMAL_COMPRESSIBILITY" ) )
        propertiesTable->AddConstProperty( iProperty->c_str(), dbTable->GetD( *iProperty ) / pascal ); // Convert units
      else if( *iProperty == string( "NUM_COMP" ) )
        propertiesTable->AddConstProperty( iProperty->c_str(), dbTable->GetI( *iProperty ) );
      else if( *iProperty == string( "ABSLENGTH_SCALING" ) )
        absScaling = true;
      else if( *iProperty == string( "RSLENGTH_SCALING" ) )
        rsScaling = true;
      else if( *iProperty == string( "MIEHG_SCALING" ) )
        mieScaling = true;
      else if( *iProperty == string( "RINDEX" ) || *iProperty == string( "KINDEX" ) ||
               *iProperty == string( "SCINTILLATION" ) || iProperty->find( "SCINTILLATION_WLS" ) != string::npos ||
               iProperty->find( "ABSLENGTH" ) != string::npos || iProperty->find( "RSLENGTH") != string::npos ||
               iProperty->find( "MIEHG" ) != string::npos || *iProperty == string( "GROUPVEL" ))
        // Note this loads in optical component vectors, rescale occurs later if needed
        propertiesTable->AddProperty( iProperty->c_str(), LoadWavelengthPropertyVector( *iProperty, dbTable ) );
      else if( *iProperty == string( "SCINTWAVEFORM" ) || *iProperty == string( "SCINTWAVEFORMalpha" ) ||
               iProperty->find( "REEMITWAVEFORM" ) != string::npos || *iProperty == string( "SCINTMOD" ) ||
               *iProperty == string( "SCINTMODalpha" ) || *iProperty == string( "SCINTMODneutron" ) )
        propertiesTable->AddProperty( iProperty->c_str(), LoadPropertyVector( *iProperty, dbTable ) );
      else if (*iProperty == string("SCINT_RISE_TIME"))
        propertiesTable->AddConstProperty(iProperty->c_str(), dbTable->GetD(*iProperty));
      else
        Log::Die( "Optics::LoadOpticalProperties: Unknown property " + *iProperty );
    }
  try
    {
      const int numComp = dbTable->GetI( "NUM_COMP" );
      // Must duplicate the SCINTILLATION_WLS data (if it exists), this should change with #191
      G4MaterialPropertyVector* current = propertiesTable->GetProperty( "SCINTILLATION_WLS" );
      if( current != NULL )
        {
          for( int component = 0; component < numComp; component++ )
            {
              stringstream propertyName; propertyName << "SCINTILLATION_WLS" << component;
              if( propertiesTable->GetProperty( propertyName.str().c_str() ) != NULL )
                continue;
              // Copy the SCINTILLATION_WLS vector, this has to be done awkwardly...
              G4MaterialPropertyVector* copy = new G4MaterialPropertyVector();
              for( size_t bin = 0; bin < current->GetVectorLength(); bin++ )
                copy->InsertValues( current->Energy( bin ), (*current)[bin] );
              propertiesTable->AddProperty( propertyName.str().c_str(), copy );
            }
        }
    }
  catch( DBNotFoundError& e ) { /* Not a problem */ }
  // Now rescale ABSLENGTH, RSLENGTH and MIEHG if there are multiple components but Only rescale if a single property doesn't exist, i.e. the scaling is required
  if( absScaling )
    {
      if( propertiesTable->GetProperty( "ABSLENGTH" ) == NULL )
        RescaleProperty( dbTable, propertiesTable, "ABSLENGTH" );
      else
        Log::Die( "Optics::LoadOpticalProperties: ABSLENGTH is defined alongside ABSLENGTH_SCALING." );
    }
  if( rsScaling )
    {
      if( propertiesTable->GetProperty( "RSLENGTH" ) == NULL )
        RescaleProperty( dbTable, propertiesTable, "RSLENGTH" );
      else
        Log::Die( "Optics::LoadOpticalProperties: RSLENGTH is defined alongside RSLENGTH_SCALING." );
    }
  if( mieScaling )
    {
      if( propertiesTable->GetProperty( "MIEHG" ) == NULL )
        RescaleProperty( dbTable, propertiesTable, "MIEHG" );
      else
        Log::Die( "Optics::LoadOpticalProperties: MIEHG is defined alongside MIEHG_SCALING." );
    }
}

void
Optics::RescaleProperty( DBLinkPtr dbTable,
                         G4MaterialPropertiesTable* propertiesTable,
                         const std::string& property )
{
  vector<double> scalings;
  try
    {
      scalings = dbTable->GetDArray( property + "_SCALING" );
    }
  catch( DBNotFoundError& e )
    {
      warn << "No " + property + " for " + dbTable->GetIndex() << newline;
      return;
    }
  vector<double> reciprocalSum;
  vector<double> energies;
  for( size_t iScale = 0; iScale < scalings.size(); iScale++ )
    {
      stringstream propertyName; propertyName << property << iScale;
      G4MaterialPropertyVector* current = propertiesTable->GetProperty( propertyName.str().c_str() );
      if( current == NULL ) // No guarantee this component has this property
        continue;
      if( reciprocalSum.empty() )
        {
          reciprocalSum.resize( current->GetVectorLength(), 0.0 );
          energies.resize( current->GetVectorLength(), -1.0 );
        }
      else if( reciprocalSum.size() != current->GetVectorLength() )
        Log::Die( "Optics::RescaleProperty: " + property + " components are different length arrays." );

      current->ScaleVector( 1.0, 1.0 / scalings[iScale] ); // Scale the values (not energy)
      for( size_t bin = 0; bin < current->GetVectorLength(); bin++ )
        {
          reciprocalSum[bin] += 1.0 / (*current)[bin];
          energies[bin] = current->Energy( bin );
        }
    }
  // This is empty if no component data exists, in which case don't add the property to the properties table
  if( reciprocalSum.empty() )
    return;
  G4MaterialPropertyVector* total = new G4MaterialPropertyVector();
  for( size_t bin = 0; bin < energies.size(); bin++ )
    total->InsertValues( energies[bin], 1.0 / reciprocalSum[bin] );
  propertiesTable->AddProperty( property.c_str(), total );
}

G4MaterialPropertyVector*
Optics::LoadPropertyVector( const std::string& propertyName,
                            DBLinkPtr dbTable )
{
  vector<double> val1, val2;
  try
    {
      val1 = dbTable->GetDArray( propertyName + "_value1" );
      val2 = dbTable->GetDArray( propertyName + "_value2" );
    }
  catch( DBNotFoundError& e ) { Log::Die( "Optics::LoadOpticalProperties: Values missing for " + propertyName ); }
  Log::Assert( val1.size() == val2.size() && !val1.empty(), "Optics::LoadOpticalProperties: Value arrays are miss-sized for " + propertyName );
  G4MaterialPropertyVector* propertyVector = new G4MaterialPropertyVector();
  int start = 0, dir = +1;
  if( val1.front() > val1.back() ) // Must reverse order
    {
      start = val1.size() - 1;
      dir = -1;
    }
  for(int index = start; index >=0 && index < int(val1.size()); index += dir)
    propertyVector->InsertValues( val1[index], val2[index] );
  return propertyVector;
}

G4MaterialPropertyVector*
Optics::LoadWavelengthPropertyVector( const std::string& propertyName,
                                      DBLinkPtr dbTable )
{
  enum EWavelength { eWavelength, eDyWavelength, eEnergy };
  int wavelengthOption = eEnergy;
  try
    {
      const string option = dbTable->GetS( propertyName + "_option" );
      if( option == string( "wavelength" ) )
        wavelengthOption = eWavelength;
      else if( option == string( "dy_dwavelength" ) )
        wavelengthOption = eDyWavelength;
      else if( option == string( "energy" ) )
        wavelengthOption = eEnergy;
      else
        Log::Die( "Optics::LoadPropertyVector: Unknown wavelength option for " + propertyName );
    }
  catch( DBNotFoundError& e ) { Log::Die( "Optics::LoadPropertyVector: Unknown wavelength option for " + propertyName ); }
  vector<double> val1, val2;
  try
    {
      val1 = dbTable->GetDArray( propertyName + "_value1" );
      val2 = dbTable->GetDArray( propertyName + "_value2" );
    }
  catch( DBNotFoundError& e ) { Log::Die( "Optics::LoadPropertyVector: Values missing for " + propertyName ); }
  Log::Assert( val1.size() == val2.size() && !val1.empty(), "Optics::LoadPropertyVector: Value arrays are miss-sized for " + propertyName );
  // Now must ensure created G4MaterailPropertyVector is ordered ascending in energy. By default go forward over values.
  int start = 0, dir = +1;
  if( wavelengthOption != eEnergy && val1.front() < val1.back() ) // Must invert the order as energy is 1/wavelength
    {
      start = val1.size() - 1;
      dir = -1;
    }
  else if( wavelengthOption == eEnergy && val1.front() > val1.back() ) // Must invert the energy order
    {
      start = val1.size() - 1;
      dir = -1;
    }
  G4MaterialPropertyVector* propertyVector = new G4MaterialPropertyVector();
  for( int index = start; index >= 0 && index < static_cast<int>( val1.size() ); index += dir )
    {
      double energy = val1[index];
      double value = val2[index];
      switch( wavelengthOption )
        {
        case eDyWavelength:
          energy = twopi * hbarc / ( val1[index] * nm );
          value *= val1[index] / energy; // Times value by the
          break;
        case eWavelength:
          energy = twopi * hbarc / ( val1[index] * nm );
          break;
        }
      propertyVector->InsertValues( energy, value );
    }
  return propertyVector;
}