// Author: Andrew M. Rogers, NSCL 07/01/2006
// Original Design: Daniel Bazin, NSCL
//* Copyright (C) 2006-2008 Andrew M. Rogers

#include "TS800.h"
#include "TS800Map.h"

#include <iostream>
#include <math.h>

ClassImp(TS800Map);

//////////////////////////////////////////////////////////////////////////
//                                                                      //
// TS800Map                                                             //
//                                                                      //
//                                                                      //
//////////////////////////////////////////////////////////////////////////

//______________________________________________________________________________
TS800Map::TS800Map(const TS800Map & map) : TObject(map)
{
  // -- Copy Constructor.
  ((TS800Map&)map).Copy(*this);
}


//______________________________________________________________________________
void TS800Map::Copy(TObject &map) const
{
  // -- Copy this method.

  TObject::Copy((TObject&)map);
  for(Int_t i=0; i<TS800_TRACK_PARAMETERS; i++)       ((TS800Map&)map).maxcoefficient[i] = maxcoefficient[i];

  for(Int_t i=0; i<TS800_TRACK_PARAMETERS; i++){
    for(Int_t j=0; j<TS800_TRACK_COEFFICIENTS; j++) { ((TS800Map&)map).order[i][j]       = order[i][j];
                                                      ((TS800Map&)map).coefficient[i][j] = coefficient[i][j];}}

  for(Int_t i=0; i<TS800_TRACK_PARAMETERS; i++){
    for(Int_t j=0; j<TS800_TRACK_PARAMETERS; j++){
      for(Int_t k=0; k<TS800_TRACK_COEFFICIENTS; k++){((TS800Map&)map).exponent[i][j][k] = exponent[i][j][k];}}}

  ((TS800Map&)map).maxorder = maxorder;
  for(Int_t i=0; i<120; i++) ((TS800Map&)map).title[i] = title[i];
  
  // Clear();
}


//______________________________________________________________________________
void TS800Map::InitClass(TS800 *itop) 
{
  // --
  
  parent = itop;

  for (Int_t i=0; i<TS800_TRACK_PARAMETERS; i++) {
    maxcoefficient[i] = 0;
    for (Int_t j=0; j<TS800_TRACK_COEFFICIENTS; j++) {
      order[i][j] = 0;
      coefficient[i][j] = 0;
      for (Int_t k=0; k<TS800_TRACK_PARAMETERS; k++) {
	exponent[i][k][j] = 0;
      }
    }
  }
  maxorder = 0;
}

//______________________________________________________________________________
Double_t* TS800Map::ReadMap(char* filename) 
{
  // -- Read the inverse map file generated by the S800 map server.
  // Returns:  Array with the Brho,A, and Z.
  //
  // The map file is formatted starting with a title line containing the Brho, A, Z, and
  // the experiment number.  There are four sections that follow each containing the 
  // coefficients used to calculate ata, yta, bta, dta in that order.
  //
  // Each line in a section contains,
  //      1.  Index:       Number used to index a parameter.
  //      2.  Coefficient: 
  //      3.  Order:     
  //      4.  Exponents:   6 exponents (4 used) for each input (xfp, afp, yfp, bfp).
  // 
  // In summary the format is,
  //
  // <Title> - Brho=##### - M=## - Q=## - Exp. #####
  //      <SECTION 1 (ata paramter)>     
  //      I  COEFFICIENT            ORDER EXPONENTS
  //      1  A1                       1    1 0  0 0  0 0
  //      2  A2                       1    0 1  0 0  0 0
  //      3  A3                       2    2 0  0 0  0 0
  //      .   .                       .          .
  //      .   .                       .          .
  //      N  AN                      nOrd  <xfp afp  |  yfp bfp  |  NONE NONE>
  //      <SECTION 2 (yta parameter)>
  //      I  COEFFICIENT            ORDER EXPONENTS
  //      1  B1                       1    1 0  0 0  0 0
  //      2  B2                       1    0 1  0 0  0 0
  //      3  B3                       2    2 0  0 0  0 0
  //      .   .                       .          .
  //      .   .                       .          .
  //      N  BN                      nOrd  <xfp afp  |  yfp bfp  |  NONE NONE>
  //      <SECTION 3 (bta parameter)>
  //      I  COEFFICIENT            ORDER EXPONENTS
  //      1  C1                       1    1 0  0 0  0 0
  //      2  C2                       1    0 1  0 0  0 0
  //      3  C3                       2    2 0  0 0  0 0
  //      .   .                       .          .
  //      .   .                       .          .
  //      N  CN                      nOrd  <xfp afp  |  yfp bfp  |  NONE NONE>
  //      <SECTION 4 (dta parameter)>
  //      I  COEFFICIENT            ORDER EXPONENTS
  //      1  D1                       1    1 0  0 0  0 0
  //      2  D2                       1    0 1  0 0  0 0
  //      3  D3                       2    2 0  0 0  0 0
  //      .   .                       .          .
  //      .   .                       .          .
  //      N  DN                      nOrd  <xfp afp  |  yfp bfp  |  NONE NONE>
  // 

  Double_t rVal[3] = {0,0,0};  // Return values.

  FILE*  file;
  Char_t line[80];
  Int_t  index, ord, par, exp[6];
  Double_t co;

  // Try to open the map file.
  file = fopen(filename, "r");
  if (file == NULL) {
    cout << "Sorry I couldn't find the file: " << filename << ". Aborting ..." << endl;
    return 0;
  }
  // Get the title and extract the Brho, A, and Z.
  fgets(title, 120, file);
  sscanf(title,"S800 inverse map - Brho=%lf - M=%lf - Q=%lf - Exp. %*i",&rVal[0],&rVal[1],&rVal[2]);

  // Read the map parameters.
  while (strstr(line, "COEFFICIENT") == NULL) fgets(line, 80, file);
  par = 0;
  while (!feof(file)) {
    fgets(line, 80, file);
    while (strstr(line, "------------------") == NULL) {
      sscanf(line, "%d %lf %d %d %d %d %d %d %d", &index, &co, &ord, &exp[0], &exp[1], &exp[2], &exp[3], &exp[4], &exp[5]);
      if (index > TS800_TRACK_COEFFICIENTS) {
	cout << "Too many coefficients in map.  Increase TS800_TRACK_COEFFICIENTS." << endl;
	break;
      }
      if (par > TS800_TRACK_PARAMETERS) {
	cout << "Too many parameters in map.  Increase TS800_TRACK_PARAMETERS." << endl;
	break;
      }
      maxcoefficient[par] = index;
      order[par][index-1] = ord;
      coefficient[par][index-1] = co;
      for (Int_t k=0; k<TS800_TRACK_PARAMETERS; k++) exponent[par][k][index-1] = exp[k];
      fgets(line, 80, file);
    }
    if (ord > maxorder) maxorder = ord;
    par++;
  }
  cout << " Done reading map from " << filename << "." << endl;
  cout << " Title: " << title;
  cout << " Order: " << maxorder << endl;
  fclose(file);

  // Set the fpTrack data members.
  //parent->fpTrack.brho   = rVal[0];
  parent->kBrho          = rVal[0];
  parent->fpTrack.mass   = (Int_t) rVal[1];
  parent->fpTrack.charge = (Int_t) rVal[2];

  return rVal;
}


//______________________________________________________________________________
Double_t TS800Map::Calculate(Int_t calcorder, Int_t parameter, Double_t *input) 
{
  // -- Calculate an S800 parameter using the invsere map.
  // 
  // The input is an array currently consisting of {xfp,afp,yfp,bfp}.
  //

  Double_t cumul=0;
  Double_t multiplicator;
  for (Int_t index=0; index<maxcoefficient[parameter]; index++) {
    if (calcorder < order[parameter][index]) break;
    multiplicator = 1;
    for (Int_t nex=0; nex<TS800_TRACK_PARAMETERS; nex++) {
      if (exponent[parameter][nex][index] != 0) {
	multiplicator *= pow(input[nex], exponent[parameter][nex][index]);
      }
    }
    cumul += multiplicator * coefficient[parameter][index];
  }
  return cumul;
}