//
//
#include "Event.h"
#include "Database.h"
#include "TFile.h" 
#include "TNtuple.h"
#include "TMath.h"
#include "opt_parameter.h"
#include <cstdlib>
#include <cstdio>
#include <iostream>
#include <fstream>

extern database db;

using namespace std;

void data::read_data(void)
{ 
  FILE* fd=fopen(db.dat_file,"rb");
  if (fd == NULL) {
    cerr << "Can't find data file " << db.dat_file << " , exiting" <<'\n';
    exit (8);
  }
  //
  // For single arm optimizations datafiles should have:
  // peak number the events belongs to, x_det, theta_det, y_det, phi_det,
  // x_tg, beamx and beamy
  // For emiss optimizaiton coordinates are read for both spectrometers.
  // In this case energy loss for the beam electron and the proton are 
  // also required.
  //
  tot_ev=0;
  for(int i=1; i<=NEVENT; i++) {
    int read_ok = 0;
    Event& evt = ev[i];
    if(db.fit_indx<=4) {
      sptr& spc = evt.sp[db.spectro];
      read_ok=fscanf(fd,"%i %lf %lf %lf %lf %lf %lf %lf", 
		     &evt.ipeak, &spc.xdet,
		     &spc.thdet, &spc.ydet,
		     &spc.phdet, &spc.xtg,
		     &spc.beamx, &spc.beamy);
    }
    else if(db.fit_indx==10) {
      sptr& sp1 = evt.sp[1];
      sptr& sp2 = evt.sp[2];
      read_ok=fscanf(fd,"%lf %lf %lf %lf %lf %lf %lf "
		     "%lf %lf %lf %lf %lf %lf %lf", 
		     &evt.pmiss,   &sp1.xdet,
		     &sp1.thdet,   &sp1.ydet,
		     &sp1.phdet,   &sp1.xtg,
		     &sp2.xdet,    &sp2.thdet,
		     &sp2.ydet,    &sp2.phdet,
		     &sp2.xtg,     &evt.pe0_loss,
		     &evt.pe_loss, &evt.px_loss);
      evt.ipeak = (read_ok > 0) ? 1 : 0;
    }
// index 5 is not used
//     else if(db.fit_indx==5) {
//       read_ok = fscanf(fd,"%i  %lf %lf %lf %lf %lf "
// 		       "%lf %lf %lf %lf %lf %lf %lf %lf %lf", 
// 		       &evt.ipeak, &evt.pmiss, &sp1.xdet,
// 		       &sp1.thdet, &sp1.ydet,
// 		       &sp1.phdet, &sp1.xtg,
// 		       &sp2.xdet,  &sp2.thdet,
// 		       &sp2.ydet,  &sp2.phdet,
// 		       &sp2.xtg,   &evt.pe0_loss,
// 		       &evt.pe_loss, &evt.px_loss);
//       evt.ipeak = (read_ok > 0) ? 1 : 0;
//     }
    if( read_ok > 0) {
      tot_ev++;
      evt.bad = 0;
    } else
      evt.bad = 1;
  }
  fclose(fd);     
}

//
//_____________________________________________________________________________
//
double sptr::zreact(void)
{ 
  double z_react;
  const database::spectrometer& spc = db.sp[db.spectro];
  z_react=-1*(ytg-db.zoffset*sin(spc.th_0))
    /(cos(spc.th_0)*phtg+sin(spc.th_0));
  z_react=z_react+beamx*(1.0-phtg*tan(spc.th_0))/(phtg+tan(spc.th_0));
  return(z_react);
}
//
//_____________________________________________________________________________
//
double sptr::dpkin(int peak)
{ 
  /*
  const double m_Ta = 181000.000;
  const double m_O = 14903.000;
  const double m_C = 11177.928;
  const double m_H = 938.3;
  const double m_Be = 8394.793;
  */
  double dp_kin;
  double P0;
  double mom, mom_kin;
  const database::spectrometer& spc = db.sp[db.spectro];
  P0=0.0;
  for(int ii=0; ii<=3; ii++) {
    P0= P0 + spc.Gamma[ii] * pow(spc.B0[peak],ii);
  }
  //
  mom=P0*(1+dp);
  //
  // Perform kinematic corrections
  //
  // Nilanga's original calculation
  //  double kin_fact=1+(mom*sin(spc.th_0)*phtg)/db.m_A[peak];
  //  mom_kin=mom*kin_fact;
  //
  // Yi's new calculation
  /*
  double m_A = pow(db.E0-db.eloss1[peak],2)*(1-cos(spc.th_0))/(db.E0-db.eloss1[peak]-mom);
  if(abs(m_A-m_Ta)<12000) m_A = m_Ta;
  else if(abs(m_A-m_C)<1000) m_A = m_C;
  else if(abs(m_A-m_Be)<800) m_A = m_Be;
  else if(abs(m_A-m_H)<100) m_A = m_H;
  */
  double m_A = db.m_A[peak]; 
 //  double d_tg = phtg + thtg*thtg/(2*spc.th_0);
  //  mom_kin = mom + m_A*sin(spc.th_0)*d_tg/pow(1-cos(spc.th_0)+m_A/(db.E0-db.eloss1[peak]),2);
  //  precise calculation
  double theta = atan(sqrt(pow(tan(spc.th_0+atan(phtg)),2)+pow(thtg,2)));
  double e_in = db.E0 - db.eloss1[peak];
  double p_corr = m_A*e_in/(m_A+e_in*(1-cos(spc.th_0))) - m_A*e_in/(m_A+e_in*(1-cos(theta)));
  mom_kin = mom + p_corr;

  dp_kin=(mom_kin-P0)/P0;
  return(dp_kin);
}
//
//_____________________________________________________________________________
//
double sptr::ytg_yi(int peak)
{
  // for ytg optimization without sieve
  const database::spectrometer& spc = db.sp[db.spectro];
  double y = beamx*(1/cos(spc.th_0)-sin(spc.th_0)*phtg)-db.zpeak[peak]*(cos(spc.th_0)*phtg+sin(spc.th_0))-db.D;
  return(y);
}
//
//_____________________________________________________________________________
//
double Event::emiss(int peak)
{ 
  double P0_1=0.0, P0_2=0.0;
  for(int ii=0; ii<=3; ii++) {
    P0_1= P0_1 + db.sp[1].Gamma[ii]  * pow(db.sp[1].B0[peak],ii);
    P0_2= P0_2 + db.sp[2].Gamma[ii] * pow(db.sp[2].B0[peak],ii); 
  }
  //
  double mom_1=P0_1*(1+sp[1].dp);
  double mom_2=P0_2*(1+sp[2].dp);
  //
  // Energy loss corrections
  //
  double E0_c=db.E0 - pe0_loss;
  mom_1 += pe_loss;
  mom_2 += px_loss;
  //
  double omega = E0_c-mom_1;
  double Tp=sqrt(pow(mom_2,2)+pow(938.27,2))-938.27;
  double TB=pow(pmiss,2)/2/db.m_B[peak];
  double e_miss=omega-Tp-TB;
  return(e_miss);
}
//
//_____________________________________________________________________________
//
double sptr::xsiev(void)
{ 
  double x_siev;
  x_siev=(thtg)*db.sp[db.spectro].dsiev+xtg;
  return(x_siev);
}
//
//_____________________________________________________________________________
//
double sptr::ysiev(void)
{
  double y_siev;
  y_siev=phtg*db.sp[db.spectro].dsiev+ytg;
  return(y_siev);
}
//
//
//_____________________________________________________________________________
//
// double sptr::delth()
// {
//   return(xtg*THT_XTG_COR);
// }
//
//_____________________________________________________________________________
//
// double sptr::deldp()
// {
//   return(xtg/DEL_XTG_COR);
// }
//
//_____________________________________________________________________________
//
void data::plot_data(void)
{
  TFile* hfile = 
    new TFile("hsimple.root","RECREATE","Demo");
  TNtuple* ntuple = 
    new TNtuple("ntuple","Demo ntuple",
		"xfp:thfp:yfp:phfp:thtg:ytg:phtg:dpkin:xsiev:"
		"ysiev:zreact");
  for(int i=1; i<=NEVENT; i++){
    sptr& sp2 = ev[i].sp[2];
    if(!ev[i].bad){
      ntuple->Fill(sp2.xfp,sp2.thfp,sp2.yfp,
		   sp2.phfp,
		   sp2.thtg,sp2.ytg,sp2.phtg,
		   sp2.dpkin(ev[i].ipeak),
		   sp2.xsiev(),sp2.ysiev(),
		   sp2.zreact());
    }
  }
  hfile->Write();
}

//_____________________________________________________________________________
const char* Event::id() const
{
  return "$Id: Event.cpp,v 1.6 2004/02/04 22:51:02 ole Exp $";
}