#include "Database.h"
#include "TMath.h"
#include <cstdlib>
#include <cstdio>
#include <iostream>
#include <fstream>
#include <string>
//
using namespace std;
static const double pi = TMath::Pi();

//_____________________________________________________________________________
int database::read_infile()
{ 
  char dummy[100];  

  FILE* fd=fopen(infil,"rb");
  if (fd == NULL) {
    cerr << "Can't find input file, " <<  infil << " , exiting" <<endl;
    exit (8);
  }
  fscanf(fd,"%s %s\n",dummy,dat_file);
  fscanf(fd,"%s %s\n",dummy,db_in);
  fscanf(fd,"%s %s\n",dummy,db_out);
  fscanf(fd,"%s %f\n",dummy,&E0);
  fscanf(fd,"%s %i\n",dummy,&spectro);
  fscanf(fd,"%s %f\n",dummy,&sp[spectro].th);
  fscanf(fd,"%s %f\n",dummy,&vert_D);
  fscanf(fd,"%s %f\n",dummy,&D);
  fscanf(fd,"%s %f\n",dummy,&target_angle);
  fscanf(fd,"%s %i\n",dummy,&nzpeak);
  for(int i=1; i<=nzpeak; i++) {
    fscanf(fd,"%s %f\n",dummy,&zpeak[i]);
  }
  fscanf(fd,"%s %i\n",dummy,&nysievpeak);
  for(int i=1; i<=nysievpeak; i++)   {
    fscanf(fd,"%s %f\n",dummy,&ysievpeak[i]);
  }
  fscanf(fd,"%s %i\n",dummy,&nxsievpeak);
  for(int i=1; i<=nxsievpeak; i++) {
    fscanf(fd,"%s %f\n",dummy,&xsievpeak[i]);
  }
  fscanf(fd,"%s %i\n",dummy,&npeak);
  for(int i=1; i<=npeak; i++) {
    fscanf(fd,"%i %i %i %i\n",&use_this[i], &pk_z[i],&pk_y[i],&pk_x[i]);
    peak_1[i]=zpeak[pk_z[i]];
    peak_2[i]=ysievpeak[pk_y[i]];
    peak_3[i]=xsievpeak[pk_x[i]];
  }
  printf("Selected peak locations are: \n\n");
  printf("peak_number      zreact     ysieve  xsieve\n\n"); 
  for(int i=1; i<=npeak; i++)  
    printf("%i           %f        %f         %f\n", 
	   i, peak_1[i],peak_2[i],peak_3[i]);
  //
  // Calculate zoffset: ie, the intersection of the spectroemter axis 
  // with the beam line
  //
  printf("\n");
  zoffset=-1*D/(pi*sp[spectro].th/180.0);
  //
  //
  // Write a .kumac files with the sieve hole grid and target peaks
  //
/*
  FILE* kumac=fopen("z_pattern.kumac","wb");
  for(int i=1; i<=nzpeak; i++){
    fprintf(kumac,"line  -100000.0 %f 100000.0 %f \n",zpeak[i], zpeak[i]);
  }
  fclose(kumac); 
  kumac=fopen("grid_pattern.kumac","wb");
  for(int i=1; i<=nysievpeak; i++){
    fprintf(kumac,"line  %f -0.5 %f 0.5 \n",ysievpeak[i], ysievpeak[i]);
  }
  for(int i=1; i<=nxsievpeak; i++){
    fprintf(kumac,"line -0.5 %f 0.5 %f \n",xsievpeak[i], xsievpeak[i]);
  }
  fclose(kumac); 
*/
  fclose(fd); 
  return 0;
}

//
//_____________________________________________________________________________
//   
int database::read_dp()
{ 
  char dummy[100];
  
  FILE* fd=fopen(infil,"rb");
  if (fd == NULL){
    cerr << "Can't find input file, " << infil << " , exiting" <<endl;
    exit (8);
  }
  fscanf(fd,"%s %s\n",dummy,dat_file);
  fscanf(fd,"%s %s\n",dummy,db_in);
  fscanf(fd,"%s %s\n",dummy,db_out);
  fscanf(fd,"%s %f\n",dummy,&E0);
  fscanf(fd,"%s %i\n",dummy,&spectro);
  fscanf(fd,"%s %f\n",dummy,&sp[spectro].th);
  fscanf(fd,"%s %f\n",dummy,&vert_D);
  fscanf(fd,"%s %f\n",dummy,&D);
  //
  fscanf(fd,"%s %i\n",dummy,&npeak);
  for(int i=1; i<=npeak; i++) {
    fscanf(fd,"%i %s %lf %lf %f %f\n",&use_this[i],dummy,
	   &sp[spectro].B0[i],&m_A[i],&eloss1[i],&eloss2[i]);
  }
  fclose(fd);
  return 0;
}
//   
//_____________________________________________________________________________
//   
int database::dpk_peak()
{ 
  double p0[MAXPEAK];
  double f_rec,e_final;
  //
  // Calculate the final momentum of the scattered electrons
  //
  cout << " Calculated dp_kin values for each peak" << endl ; 
  cout << " peak #        dp_kin " << endl ;
  //
  for(int i=1; i<=npeak; i++) { 
    p0[i]=0.0;
    for(int ii=0; ii<=3; ii++) { 
      p0[i]= p0[i] + sp[spectro].Gamma[ii]*pow(sp[spectro].B0[i],ii);
    }
    f_rec= 1 + 2*E0*pow(sin(sp[spectro].th_0/2),2)/m_A[i];
    e_final=(E0-eloss1[i])/f_rec - eloss2[i];
    peak_0[i]=(e_final- p0[i])/p0[i];
    cout << i << "   " << peak_0[i] << endl ;
  }
  return 0;
}
//
//_____________________________________________________________________________
//   
int database::read_emiss()
{ 
  char dummy[100];
  
  FILE* fd=fopen(infil,"rb");
  if (fd == NULL) {
    cerr << "Can't find input file, " << infil << "  , exiting" <<'\n';
    exit (8);
  }
  fscanf(fd,"%s %s\n",dummy,dat_file);
  fscanf(fd,"%s %s\n",dummy,db_in);
  fscanf(fd,"%s %s\n",dummy,db_out);
  fscanf(fd,"%s %f\n",dummy,&E0);
  fscanf(fd,"%s %f\n",dummy,&sp[1].th);
  fscanf(fd,"%s %f\n",dummy,&sp[2].th);
  fscanf(fd,"%s %i\n",dummy,&npeak);
  for(int i=1; i<=npeak; i++){
    fscanf(fd,"%i %s %f %lf %lf %lf %f\n",
	   &use_this[i],dummy,&peak_0[i],&sp[1].B0[i],&sp[2].B0[i],
	   &m_A[i],&m_B[i]);
  }
  fclose(fd);
  return 0;
}

//_____________________________________________________________________________
int  database::init()
{
 //
 // Note that the spectrometer constants are hardwired in
 // The values given here are good to 1.10-3 level and
 // is accuracy is sufficient for the optimizations.
 //
 // Target to sieve collimator distances are from
 // http://hallaweb.jlab.org/news/minutes/sieves_2000.htm 
 //
 // Spectrometer 1 is the Right spectrometer 
 //
 //
 // sp[1].dsiev=1.1755;
 //
 // *****NL 05/07/03 modify the distance for septum for the right
 // spectrometer
 //
  sp[1].dsiev=0.8022; //0.81
  sp[2].dsiev=0.8025; // L for left arem
//sp[2].dsiev=1.1843;

  sp[1].lines = sp[2].lines = 0;

  // Read input file. dp and emiss input files have a special format.
  int status = 0;
  if( fit_indx == 2 )
    status = read_dp();
  else if( fit_indx == 10 )
    status = read_emiss();
  else
    status = read_infile();

  // Quit on error
  if( status != 0 )
    return status;

  //
  // Assign hit spectrometer: For y_tg, angle and dp (singles)
  // this is the one spectrometer hit. For emiss both spectrometers
  //
  if(fit_indx<=4){
    hit_arm[1]=spectro;
    hit_arm[2]=spectro;
  }
  else if(fit_indx==10){
    hit_arm[1]=1;
    hit_arm[2]=2;
  }

  // convert the spectrometer angles to radian
  for ( int i = hit_arm[1]; i <=hit_arm[2]; ++i)
    sp[i].th_0=pi*sp[i].th/180.0;

  return 0;
}

//_____________________________________________________________________________
const char* database::id() const
{
  return "$Id: Database.cpp,v 1.11 2004/02/25 22:15:22 ole Exp $";
}