#include "L1AlgoInputSettings.h"
#include "L1Station.h"
#include <iostream>
using namespace std;
using std::ios;

#ifdef ROOT_IS_USED
#include "TFile.h"
#include "TProfile2D.h"
#endif

L1AlgoInputSettings::L1AlgoInputSettings( const L1AlgoInputSettings& a)
{
  SetSettings( a.GetSettings() );
}

const L1AlgoInputSettings& L1AlgoInputSettings::operator=( const L1AlgoInputSettings& a)
{
  SetSettings( a.GetSettings() );
  return a;
}

void L1AlgoInputSettings::SetSettings( const std::vector<fscal>& geo )
{
  for ( int iGeo = 0; iGeo < 10000; iGeo++ ) {
    fGeo[iGeo] = geo[iGeo];
  }
}

istream& L1AlgoInputSettings::eatwhite(istream& is) // skip spaces
{
  char c;
  while (is.get(c)) {
    if (isspace(c)==0) {
      is.putback(c);
      break;
    }
  }
  return is;
}

bool L1AlgoInputSettings::ReadSettingsFromFile(const string work_dir, const int iVerbose)
{
  string fname = work_dir + string("geo_algo.txt");
  ifstream GeoFile(fname.data());	
  if ( !GeoFile.is_open() ) return 0;
  
  if ( iVerbose > 0 ) cout << " -I- Read geometry from " << fname << endl;
  int i;
  for (i = 0; !GeoFile.eof(); i++){
    fscal tmp;
    GeoFile >> tmp >> eatwhite;
    fGeo[i] = tmp;
  };
  GeoFile.close();
  
#ifdef ROOT_IS_USED
  
  int ind = 12;
  int NStations = fGeo[ind++];
  int NMvdStations = fGeo[ind++];
  vector<L1Station> vStations(NStations);
  
  for( int iSt=0; iSt<NStations; iSt++ ){
    L1Station &st = vStations[iSt];
    st.z = fGeo[ind++];
    st.materialInfo.thick = fGeo[ind++];
    st.Rmin =  fGeo[ind++];
    st.Rmax =  fGeo[ind++];
    st.materialInfo.RL    = fGeo[ind++];
    st.materialInfo.RadThick = st.materialInfo.thick/st.materialInfo.RL;
    st.materialInfo.logRadThick = log(st.materialInfo.RadThick);
    
    ind+=4;
    int N= static_cast<int>( fGeo[ind++] );
    for( int iC=0; iC<N; iC++ ) st.fieldSlice.cx[iC] = fGeo[ind++];
    for( int iC=0; iC<N; iC++ ) st.fieldSlice.cy[iC] = fGeo[ind++];
    for( int iC=0; iC<N; iC++ ) st.fieldSlice.cz[iC] = fGeo[ind++];
  }

  fRadThick.resize(NStations);
    // MVD does not use map
  string dataDir = work_dir;
  string matBudgetFileName = dataDir;
  matBudgetFileName += "materialBudget.root";
    // Read STS Radiation Thickness table
  for( int iSta = 0; iSta < NMvdStations;iSta++ ) {
    fRadThick[iSta].SetBins(1, 100); // mvd should be in +-100 cm square
    fRadThick[iSta].table.resize(1);
    fRadThick[iSta].table[0].resize(1);
    fRadThick[iSta].table[0][0] = vStations[iSta].materialInfo.RadThick[0];
  }
  
  if ( access( matBudgetFileName.c_str(), F_OK ) != -1 ) {
    TFile* oldfile = gFile;
    TFile *rlFile = new TFile(matBudgetFileName.data());

    if( iVerbose > 0 )
      cout << " -I- STS Material budget file is " << matBudgetFileName << "." << endl;
    for( int j = 0, iSta = NMvdStations; iSta < NStations; iSta++, j++ ) {
      TString name = "Radiation Thickness [%]";
      name += ", Station";
      name += j+1; 
    
      TProfile2D* hStaRadLen = (TProfile2D*) rlFile->Get(name);
      if ( !hStaRadLen ) {
        cout << "L1: incorrect " << matBudgetFileName << " file. No " << name << endl; exit(1);
      }
      
      const int NBins = hStaRadLen->GetNbinsX(); // should be same in Y
      const float RMax = hStaRadLen->GetXaxis()->GetXmax(); // should be same as min
      fRadThick[iSta].SetBins(NBins,RMax); // TODO
      fRadThick[iSta].table.resize(NBins);
    
      for( int iB = 0; iB < NBins; iB++ ) {
        fRadThick[iSta].table[iB].resize(NBins);
        for( int iB2 = 0; iB2 < NBins; iB2++ ) {
          fRadThick[iSta].table[iB][iB2] = 0.01 * hStaRadLen->GetBinContent(iB,iB2); //0.0034;//0.003209;//
          if ( fRadThick[iSta].table[iB][iB2] < vStations[iSta].materialInfo.RadThick[0] )
          {
            fRadThick[iSta].table[iB][iB2] = vStations[iSta].materialInfo.RadThick[0];
          }
            // cout << iB << " " << iB2 << " " << fRadThick[iSta].table[iB][iB2] << endl; // dbg
        }
      }
    }
    rlFile->Close();
    rlFile->Delete();
    gFile = oldfile;
  }
  else{
    if( iVerbose > 0 )
    {
      cout  << " -W- No material budget file is found. Homogenious budget will be used" << endl;
    }
    for( int iSta = NMvdStations; iSta < NStations; iSta++ ) {
      fRadThick[iSta].SetBins(1, 100); // mvd should be in +-100 cm square
      fRadThick[iSta].table.resize(1);
      fRadThick[iSta].table[0].resize(1);
      fRadThick[iSta].table[0][0] = vStations[iSta].materialInfo.RadThick[0];
    }
  }
#endif
  if ( iVerbose > 0 ) cout << " -I- geometry has been read successfully." << endl;
  return 1;
}