//
// TAGMboundaries.h - class to compute the boundaries in focal plane position
//                    and energy of the tagger microscope fiber columns.
//
// author: richard.t.jones at uconn.edu
// version: june 27, 2014
//
// usage:
//           root[] .L TAGMboundaries.C+O
//           root[] TAGMboundaries tagm;
//           root[] tagm.Write(); 

#include <TAGMboundaries.h>
#include <TROOT.h>
#include <TTree.h>
#include <math.h>

#include <fstream>
#include <sstream>
#include <iostream>

ClassImp(TAGMboundaries);

const char *TAGMboundaries::fFPtable_struct =
           "FPx_cm/D:Ephoton_GeV/D:beta_deg/D";

TAGMboundaries::TAGMboundaries(const char *rootfile)
 : fRootfile(0)
{
   fXorigin_cm = 336.22;
   fYorigin_cm = 0;
   fEendpoint_GeV = 12.0;

   fTAGM_sums_ncolumns = 102;
   fTAGM_singles_ncolumns = 4;
   fTAGM_singles_column = new int[4];
   fTAGM_singles_column[0] = 4;
   fTAGM_singles_column[1] = 24;
   fTAGM_singles_column[2] = 48;
   fTAGM_singles_column[3] = 96;

   fTAGM_fiber_width_cm = 0.1880;
   fTAGM_fiber_pitch_cm = 0.2150;
   fTAGM_fiber_beta_deg = new double[fTAGM_sums_ncolumns];
   for (int col=0; col < fTAGM_sums_ncolumns; ++col) {
      fTAGM_fiber_beta_deg[col] = 12.50-(col/6)*0.09;
   }
   if (rootfile != 0)
      fRootfile = new TFile(rootfile,"UPDATE");
   else
      fRootfile = new TFile("TAGMboundaries.root","UPDATE");
}

TAGMboundaries::~TAGMboundaries()
{
   delete fTAGM_fiber_beta_deg;
   if (fRootfile != 0)
      fRootfile->Close();
}

int TAGMboundaries::loadFPtable_from_textfile(const char *textfile)
{
   ifstream tfile(textfile);
   if (!tfile.is_open()) {
      std::cerr << "Error in TAGMboundaries::loadFPtable_from_textfile:"
                << std::endl << "  input text file " << textfile
                << " not found." << std::endl;
      return 0;
   }

   fFPtable.clear();
   int entry=0;
   while (tfile.good()) {
      char tline[99];
      tfile.getline(tline,99);
      std::stringstream tstr(tline);
      double Eelectron, FPx, beta;
      tstr >> Eelectron >> FPx >> beta;
      if (Eelectron == 0 || FPx == 0 || beta == 0)
         continue;
      fFPtable_entry.Ephoton_GeV = fEendpoint_GeV - Eelectron;
      fFPtable_entry.FPx_cm = FPx;
      fFPtable_entry.beta_deg = beta;
      fFPtable.push_back(fFPtable_entry);
      ++entry;
   }
   return fFPtable.size();
}

int TAGMboundaries::loadFPtable_from_rootfile(const char *rootfile)
{
   TFile *rfile;
   if (rootfile != 0)
      rfile = new TFile(rootfile,"UPDATE");

   TTree *fptree = (TTree*)gROOT->FindObject("FPtree");
   if (fptree == 0) {
      std::cerr << "Error in TAGMboundaries::loadFPtable_from_rootfile: "
                << std::endl << "  no object FPtree found in input file."
                << std::endl;
      return 0;
   }

   fFPtable.clear();
   TBranch *b_FPtable;
   fptree->SetBranchAddress("FPtable", &fFPtable_entry.FPx_cm, &b_FPtable);
   int len = fptree->GetEntries();
   for (int entry=0; entry < len; ++entry) {
      b_FPtable->GetEntry(entry);
      fFPtable.push_back(fFPtable_entry);
   }
   return fFPtable.size();
}

int TAGMboundaries::saveFPtable_to_rootfile(const char *rootfile)
{
   TFile *rfile;
   if (rootfile != 0)
      rfile = new TFile(rootfile,"UPDATE");
   else if (fRootfile != 0)
      rfile = fRootfile;
   else {
      std::cerr << "Error in TAGMboundaries::saveFPtable_to_rootfile: "
                << std::endl << "  no output root file specified."
                << std::endl;
      return 0;
   }

   TTree *fptree = new TTree("FPtree","tagger focal plane map table");
   fptree->Branch("FPtable",&fFPtable_entry,fFPtable_struct);
   unsigned int entry;
   for (entry=0; entry < fFPtable.size(); ++entry) {
      fFPtable_entry = fFPtable[entry];
      fptree->Fill();
   }
   fptree->Write();
   return entry;
}

double TAGMboundaries::PhotonEnergy_GeV(double x_cm, double y_cm)
{
   if (fFPtable.size() == 0) {
      if (loadFPtable_from_rootfile() == 0) {
         if (loadFPtable_from_textfile("12GeV-10MeVrays-6-2014.txt") == 0) {
            std::cerr << "Error in TAGMboundaries::fitFP_Ephoton_vs_FPx: "
                      << std::endl << "  unable to load FP tables"
                      << " from any known source." << std::endl;
            return 0;
         }
         saveFPtable_to_rootfile();
      }
   }

   // interpolate the table, ordered in increasing FPx_cm
 
   unsigned int entry;
   double x[2] = {0,0};
   for (entry=0; entry < fFPtable.size(); ++entry) {
      double beta = fFPtable[entry].beta_deg * M_PI/180;
      x[0] = x[1];
      x[1] = fFPtable[entry].FPx_cm + y_cm/tan(beta);
      if (x[1] > x_cm)
         break;
   }
   if (entry == 0 || entry == fFPtable.size()) {
      return 0;
   }
   double f = (x[1] - x_cm)/(x[1] - x[0]);
   return f*fFPtable[entry-1].Ephoton_GeV +
          (1-f)*fFPtable[entry].Ephoton_GeV;
}

int TAGMboundaries::Write(const char *textfile)
{
   ofstream tfile;
   if (textfile == 0)
      tfile.open("TAGMboundaries.txt");
   else
      tfile.open(textfile);
   if (!tfile.is_open()) {
      std::cerr << "Error in TAGMboundaries::Write: "
                << std::endl << "  unable to open output text file."
                << std::endl;
      return 0;
   }

   tfile << 
   "column xlow(cm)  xhigh(cm)   y(cm)  Elow(GeV)  Ehigh(GeV)   Flow   Fhigh"
         << std::endl <<
   "--------------------------------------------------------------------------"
         << std::endl;
  
   double x = fXorigin_cm;
   double y = fYorigin_cm;
   int chan = 1;
   for (int col=0; col < fTAGM_sums_ncolumns; ++col) {
      double tanbeta = tan(fTAGM_fiber_beta_deg[col] * M_PI/180);
      double xlow = x + (fTAGM_fiber_pitch_cm + fTAGM_fiber_width_cm)/
                        (2*tanbeta);
      double xhigh = xlow - fTAGM_fiber_width_cm/tanbeta;
      double Elow = PhotonEnergy_GeV(xlow,y);
      double Ehigh = PhotonEnergy_GeV(xhigh,y);
      double Flow = Elow/fEendpoint_GeV;
      double Fhigh = Ehigh/fEendpoint_GeV;
      x += fTAGM_fiber_pitch_cm/tanbeta;
      char line[80];
      sprintf(line,"%4d   %8.4f   %8.4f  %7.4f %8.4f   %8.4f   %7.4f %7.4f",
                    chan, xlow, xhigh, y, Elow, Ehigh, Flow, Fhigh);
      tfile << line << std::endl;
      ++chan;
   }   
#if APPEND_SINGLE_CHANNELS_TABLE
   x = fXorigin_cm;
   y = fYorigin_cm;
   int scol = 0;
   for (int col=0; col < fTAGM_sums_ncolumns; ++col) {
      double tanbeta = tan(fTAGM_fiber_beta_deg[col] * M_PI/180);
      double xlow = x + (fTAGM_fiber_pitch_cm + fTAGM_fiber_width_cm)/
                        (2*tanbeta);
      double xhigh = xlow - fTAGM_fiber_width_cm/tanbeta;
      double Elow = PhotonEnergy_GeV(xlow,y);
      double Ehigh = PhotonEnergy_GeV(xhigh,y);
      double Flow = Elow/fEendpoint_GeV;
      double Fhigh = Ehigh/fEendpoint_GeV;
      x += fTAGM_fiber_pitch_cm/tanbeta;
      if (col < fTAGM_singles_column[scol])
         continue;
      for (int j=0; j < 5; ++j) {
         char line[80];
         sprintf(line,"%4d   %8.4f   %8.4f  %7.4f %8.4f   %8.4f   %7.4f %7.4f",
                       chan, xlow, xhigh, y, Elow, Ehigh, Flow, Fhigh);
         tfile << line << std::endl;
         ++chan;
      }
      if (++scol == fTAGM_singles_ncolumns)
         break;
   }
#endif
   tfile.close();
   return chan;
}