#define speed_calib_cxx
#include "speed_calib.h"

#include "TCanvas.h"
#include "TGraph.h"
#include "TH2.h"
#include "TMath.h"
#include "TNtuple.h"
#include "TStyle.h"
#include <fstream.h>
#include <iomanip.h>       
#include <iostream>
#include <math.h>
#include <string.h>
#include <TMinuit.h>
 
const double c      =    29.9792;
const double L      =    279.500;
const double kmass  =   0.493677;
const double pimass = 0.13956995;
const double pmass  = 0.93827231;
const double r2pi   = 2.50662827479464932329;
const int    maxevt =      17000;
const int    NumBar =        216;
const int    npid   =          3;

//const int max_t0_trks = 50;

// For attenuation length
int          evts, jchan;
double       length[maxevt], lenght[maxevt], zps[maxevt];
int          nvt[NumBar];
double       path[NumBar][maxevt];
double att[maxevt],zpp[maxevt], attp[NumBar][maxevt], zpos[NumBar][maxevt],attpp[NumBar][maxevt];
double q_east[NumBar][maxevt], q_west[NumBar][maxevt];
double var_att[maxevt][NumBar],var_zpos[maxevt][NumBar];

double t0ge, t0gw, t0gei, t0gwi, t0g;
double tof, etof, tofe, tofw, tof_k, tof_p, tof_pi, tof_m;

char filesl[50], titlesl[50];

//--------------------------------------------------------------------------
// Function to be minimized by Minuit for the attenuation length
//--------------------------------------------------------------------------

void fcn_att(int    &npar,
             double *gin,
             double &f,
             double *par,
             int     iflag)
{
  const double r2pi = 2.50662827479464932329;
  double lambda  = par[0];
  double of = par[1];
  double rn = par[2];
  double ff = par[3];
  double rw = par[4];
  double lambdab  = par[5];
  double ofb = par[6];
  double logl = 0.;
  for ( int indx=0; indx<evts; indx++ ) {
        double dy  = zps[indx]-(lambda*lenght[indx]/2.+of);
	double dyb = zps[indx]-(lambdab*lenght[indx]/2.+ofb);
        double arg = ff*(1/(r2pi*rn))*exp(-dy*dy/(2*rn*rn))
                     + ((1-ff)/(r2pi*rw))*exp(-dyb*dyb/(2*rw*rw));
           if(arg>0)logl += log(arg);      
  }
  f = -2*logl;
}

void fcnatt(int &npar,double *gin,double &f,double *par,int iflag){

  double q = par[0];
  double a = par[2];
  double AttLen = par[1];
  double sigmas ;
  double logl  = 0.;
  double arg;
  double polz, z, t;
  for ( int i=0; i<evts; i++ ) {

    t     = att[i]; // ACD/(path-length*atl)
//     if (t>4000.) continue;
    if (jchan<216) {
      z = L/2-zpp[i];
    } else {
      z = L/2+zpp[i];
    }

    polz  = exp(-(z/AttLen) + a*z*z);
    sigmas = polz*par[3];
//     arg = (1.-par[4])*(TMath::Landau(t,q*polz,sigmas))/sigmas + par[4]/5250.;
    arg = (TMath::Landau(t,q*polz,sigmas))/sigmas;

//     if(arg>0.)logl += log(arg);

    double limit = TMath::Landau(5000.,0.,sigmas)/sigmas;
    if(arg>limit)logl += log(arg);
    else logl += log(limit);
    
  }  

 f = -2*logl;
}

// Global variables for the t0 calculator
// ======================================
int    Nmatches_sl, _skip_sl; double _skip_p;
int    Npulses[200];
double _m_sl[3] = {0.13956995, 0.493677, 0.93827231};
double _f_sl[3] = {0.8, 0.1, 0.1};
double _s[200][4], _p[200][4], _t[200][4], _et[200][4];
int tof_ntt, tof_ttpBar[200][4], tof_nttp[200], tof_nttpTrk[200], tof_inBar[200], zv_nzv;
double zv_zv[20], tof_hbInZ[200][4], tof_hb_s[200][4], tof_pp[200], tof_pt[200];
double tof_pulseQ[200][4], tp_tRaw[2][216], tof_ctz0bc[200], tof_hbPathLen[200][4], tof_ttpMip[200][4];

// Methods used to calculate the t0 of the event:
// getTofSL computes the time of flight
// ==============================================
float getTofSL(float s_cm, float p_gev, float m_gev) {

  float c_cm_per_ns = 29.9792458;
  float y;
  if (p_gev){
    y = m_gev / p_gev;
    return (s_cm / c_cm_per_ns)*sqrt(1+y*y);
  }
  return 0.;

} // end getTofSL(float s_cm, float p_gev, float m_gev)

// getMeanSL
// =========
void getMeanSL(double &mean, double &sigma, int size, double value[4], double error[4]) {

  int i;
  double w,sum_x,sum_w;
    
  sum_x = sum_w = 0.;
  for (i=0;i<size;i++){
    if (error[i])
      w = 1./(error[i]*error[i]);
    else
      w = 1.;
    sum_w += w;
    sum_x += w*value[i];
  }
  if (sum_w > 0. ){
    mean  = sum_x / sum_w;
    sigma = 1. / sqrt(sum_w);
  }

} // void getMeanSL(double &mean, double &sigma, int size, double value[4], double error[4])

// Function for t0 fit
// ===================
void pid_likelihood_fcn_sl(int &npar, double *gin, double &f, double *par, int iflag) {

  static const double norm = 2.5066283;
  double diff, sig;
  double outer, inner, prod;

// Calculate the likelihood:
//
//  L = Prod(tracks)[Sum(pids)[f(pid) Prod(pmts)[G(t0-t0(s,p,t,pid),sigma(z))]]]
//
// or (pijk = prob from jth match, kth pmt, that its particle i):
//
//  L = (f1*p111*p112*...*p11M + f2*p211*p212*...*p21M + f3*p311*p312*...*p31M)
//    * (f1*p121*p122*...*p12M + f2*p221*p222*...*p22M + f3*p321*p322*...*p32M)
//    ...
//

  outer = 0;
  for (int i = 0; i < Nmatches_sl; i++)
    {
      //don't use the match corresponding to skip, for an unbiased t0:
      if (i == _skip_sl) continue;
      inner = 0;
      for (int pid = 0; pid<3; pid++){
        prod = _f_sl[pid];
        for (int j = 0; j < Npulses[i]; j++){
          diff  = _t[i][j] - getTofSL(_s[i][j],_p[i][j],_m_sl[pid]) - par[0];
          sig   = _et[i][j];
          if (sig!=0)
            prod *= exp(-diff*diff/(2*sig*sig))/(norm*sig);
        }
        inner += prod;
      }
      if (inner>0.)
        outer += log(inner);
    }

  f = -2*outer;

} // void pid_likelihood_fcn_sl(int &npar, double *gin, double &f, double *par, int iflag)

// Fill parameters for track correction
double parc[NumBar*2][20], parr[NumBar*2][2], parq[NumBar*2][10];
void FillParameters(char   *directory, char   *file_name, int flag) {
  
  ifstream twe_result_par, tww_result_par, speed_result_par, offe_result_par, offw_result_par;
  ifstream zdepe_result, zdepw_result, rese_result, resw_result;
  sprintf(filesl,"%s/files/%s_timew-calib_e.txt",directory,file_name); twe_result_par.open(filesl,ios::in);
  sprintf(filesl,"%s/files/%s_timew-calib_w.txt",directory,file_name); tww_result_par.open(filesl,ios::in);
  sprintf(filesl,"%s/files/%s_spdtw-calib.txt",directory,file_name); speed_result_par.open(filesl,ios::in);
  sprintf(filesl,"%s/files/%s_T0-calib_e.txt",directory,file_name); offe_result_par.open(filesl,ios::in);
  sprintf(filesl,"%s/files/%s_T0-calib_w.txt",directory,file_name); offw_result_par.open(filesl,ios::in);
  sprintf(filesl,"%s/files/%s_Zres-calib_e.txt",directory,file_name); zdepe_result.open(filesl,ios::in);
  sprintf(filesl,"%s/files/%s_Zres-calib_w.txt",directory,file_name); zdepw_result.open(filesl,ios::in);
//   sprintf(filesl,"/home/sluca/TofCalibrations/533validation/IvanFitOut.txt",directory,file_name); zdepw_result.open(filesl,ios::in);
  sprintf(filesl,"%s/files/%s_PmtRes-calib_e.txt",directory,file_name); rese_result.open(filesl,ios::in);
  sprintf(filesl,"%s/files/%s_PmtRes-calib_w.txt",directory,file_name); resw_result.open(filesl,ios::in);
  ifstream atte_result, attw_result;
  sprintf(filesl,"%s/files/%s_Att-calib_e.txt",directory,file_name); atte_result.open(filesl,ios::in);
  sprintf(filesl,"%s/files/%s_Att-calib_w.txt",directory,file_name); attw_result.open(filesl,ios::in);
  for(int j=0;j<NumBar;j++) {

    double   resultzdepe[NumBar][16], resultzdepw[NumBar][16], resultrese[NumBar][2], resultresw[NumBar][2];
    double   resulte[NumBar][4], resultw[NumBar][4], result[NumBar][2], resultoffe[NumBar], resultoffw[NumBar];

    double   data[9],datap[9],datat[9],datoffe[5],datoffw[5],datzdepe[33],datzdepw[33], datrese[5], datresw[5];

    for(int i=0;i<4;i++) {
      resulte[j][i]=0.;
      resultw[j][i]=0.;
    }

    if (flag>=1) {
      for(int i=0;i<3;i++) twe_result_par >> data[i]; 
      for(int i=0;i<1;i++) resulte[j][i]= data[1+i*2];
      for(int i=0;i<3;i++) tww_result_par >> datap[i];
      for(int i=0;i<1;i++) resultw[j][i]= datap[1+i*2];
    }

    result[j][0]=10000000000000000.;
    result[j][1]=10000000000000000.;

    if (flag>=2) {
  
      for(int i=0;i<9;i++) speed_result_par >> datat[i];
      result[j][0]=datat[1];
      result[j][1]=datat[3];  

    }

    resultoffe[j] = 0.;
    resultoffw[j] = 0.;

    if (flag>=3) {

      for(int i=0;i<5;i++) offe_result_par >> datoffe[i];
      resultoffe[j] = datoffe[3];
      for(int i=0;i<5;i++) offw_result_par >> datoffw[i]; 
      resultoffw[j] = datoffw[3];

    }

    for(int i=0;i<16;i++) resultzdepe[j][i] = 0.0;
    for(int i=0;i<16;i++) resultzdepw[j][i] = 0.0;

    if (flag>=4) {

      for(int i=0;i<9;i++) zdepe_result >> datzdepe[i];
      for(int i=0;i<8;i++) zdepe_result >> datzdepe[i+9];
      for(int i=0;i<8;i++) zdepe_result >> datzdepe[i+17];
      for(int i=0;i<8;i++) zdepe_result >> datzdepe[i+25];
      for(int i=0;i<16;i++) resultzdepe[j][i] = datzdepe[1+i*2];

      for(int i=0;i<9;i++) zdepw_result >> datzdepw[i];
      for(int i=0;i<8;i++) zdepw_result >> datzdepw[i+9];
      for(int i=0;i<8;i++) zdepw_result >> datzdepw[i+17];
      for(int i=0;i<8;i++) zdepw_result >> datzdepw[i+25];
      for(int i=0;i<16;i++) resultzdepw[j][i] = datzdepw[1+i*2];
      
    }

    parr[j][0] = 0.100+(L/2)*0.0004;
    parr[j][1] = -0.0004; 
    parr[j+216][0] = 0.100+(L/2)*0.0004;
    parr[j+216][1] = 0.0004;

    if (flag>=5) {

      for(int i=0;i<5;i++) rese_result >> datrese[i];
      resultrese[j][0] = datrese[1]; resultrese[j][1] = datrese[3]; 
      for(int i=0;i<5;i++) resw_result >> datresw[i];
      resultresw[j][0] = datresw[1]; resultresw[j][1] = datresw[3]; 
      
      parr[j][0] = resultrese[j][0]+(L/2)*resultrese[j][1]; 
      parr[j][1] = -resultrese[j][1]; 
      parr[j+216][0] = resultresw[j][0]+(L/2)*resultresw[j][1]; 
      parr[j+216][1] = resultresw[j][1]; 
      
    }
    
    parc[j][0] = resultoffe[j] - resultzdepe[j][0] - resultzdepe[j][1]*L/2 - resultzdepe[j][2]*L*L/4 -
      resultzdepe[j][3]*L*L*L/8 - L/(2.*result[j][0]);
    parc[j][1] = + 1./result[j][0] - resultzdepe[j][1] - resultzdepe[j][2]*L - (3./4.)*resultzdepe[j][3]*L*L; 
    parc[j][2] = - resulte[j][0];
    parc[j][3] = - resultzdepe[j][2] - (3./2.)*L*resultzdepe[j][3];
    parc[j][4] = 0.;//- resulte[j][1]; 
    parc[j][5] = 0.;//- resulte[j][2];
    parc[j][6] = - resultzdepe[j][3];
    parc[j][7] = 0.;//- resulte[j][3];
    for (int y = 0; y<12; y++) parc[j][7+y+1] = - resultzdepe[j][4+y];

    parc[j+216][0] = resultoffw[j] - resultzdepw[j][0] - resultzdepw[j][1]*L/2 - resultzdepw[j][2]*L*L/4 -
      resultzdepw[j][3]*L*L*L/8 - L/(2.*result[j][0]);
    parc[j+216][1] = - 1./result[j][0] - resultzdepw[j][1] - resultzdepw[j][2]*L - (3./4.)*resultzdepw[j][3]*L*L; 
    parc[j+216][2] = - resultw[j][0];
    parc[j+216][3] = - resultzdepw[j][2] - (3./2.)*L*resultzdepw[j][3];
    parc[j+216][4] = 0.;//- resultw[j][1]; 
    parc[j+216][5] = 0.;//- resultw[j][2];
    parc[j+216][6] = - resultzdepw[j][3];
    parc[j+216][7] = 0.;//- resultw[j][3];
    for (int y = 0; y<12; y++) parc[j+216][7+y+1] = - resultzdepw[j][4+y];

    for(int i=0;i<10;i++) atte_result >> parq[j][i];
    for(int i=0;i<10;i++) attw_result >> parq[j+216][i];

  }
    
  twe_result_par.close();
  tww_result_par.close();
  speed_result_par.close();
  offe_result_par.close();
  offw_result_par.close();
  zdepe_result.close(); 
  zdepw_result.close();
  rese_result.close();
  resw_result.close();
  atte_result.close();
  attw_result.close();

}

// Fill parameters for track correction
double t0parc[NumBar*2][20], t0parr[NumBar*2][2], t0parq[NumBar*2][10];
void FillParametersT0(char   *directory, char   *file_name, int flag) {
  
  ifstream twe_result_par, tww_result_par, speed_result_par, offe_result_par, offw_result_par;
  ifstream zdepe_result, zdepw_result, rese_result, resw_result;
  sprintf(filesl,"%s/files/%s_timew-calib_e.txt",directory,file_name); twe_result_par.open(filesl,ios::in);
  sprintf(filesl,"%s/files/%s_timew-calib_w.txt",directory,file_name); tww_result_par.open(filesl,ios::in);
  sprintf(filesl,"%s/files/%s_spdtw-calib.txt",directory,file_name); speed_result_par.open(filesl,ios::in);
  sprintf(filesl,"%s/files/%s_T0-calib_e.txt",directory,file_name); offe_result_par.open(filesl,ios::in);
  sprintf(filesl,"%s/files/%s_T0-calib_w.txt",directory,file_name); offw_result_par.open(filesl,ios::in);
  sprintf(filesl,"%s/files/%s_Zres-calib_e.txt",directory,file_name); zdepe_result.open(filesl,ios::in);
  sprintf(filesl,"%s/files/%s_Zres-calib_w.txt",directory,file_name); zdepw_result.open(filesl,ios::in);
  sprintf(filesl,"%s/files/%s_PmtRes-calib_e.txt",directory,file_name); rese_result.open(filesl,ios::in);
  sprintf(filesl,"%s/files/%s_PmtRes-calib_w.txt",directory,file_name); resw_result.open(filesl,ios::in);
  ifstream atte_result, attw_result;
  sprintf(filesl,"%s/files/%s_Att-calib_e.txt",directory,file_name); atte_result.open(filesl,ios::in);
  sprintf(filesl,"%s/files/%s_Att-calib_w.txt",directory,file_name); attw_result.open(filesl,ios::in);
  for(int j=0;j<NumBar;j++) {

    double   resultzdepe[NumBar][16], resultzdepw[NumBar][16], resultrese[NumBar][2], resultresw[NumBar][2];
    double   resulte[NumBar][4], resultw[NumBar][4], result[NumBar][2], resultoffe[NumBar], resultoffw[NumBar];

    double   data[9],datap[9],datat[9],datoffe[5],datoffw[5],datzdepe[33],datzdepw[33], datrese[5], datresw[5];

    for(int i=0;i<4;i++) {
      resulte[j][i]=0.;
      resultw[j][i]=0.;
    }

    if (flag>=1) {
      for(int i=0;i<3;i++) twe_result_par >> data[i]; 
      for(int i=0;i<1;i++) resulte[j][i]=data[1+i*2];
      for(int i=0;i<3;i++) tww_result_par >> datap[i];
      for(int i=0;i<1;i++) resultw[j][i]=datap[1+i*2];
    }

    result[j][0]=10000000000000000.;
    result[j][1]=10000000000000000.;

    if (flag>=2) {
  
      for(int i=0;i<9;i++) speed_result_par >> datat[i];
      result[j][0]=datat[1];
      result[j][1]=datat[3];  

    }

    resultoffe[j] = 0.;
    resultoffw[j] = 0.;

    if (flag>=3) {

      for(int i=0;i<5;i++) offe_result_par >> datoffe[i];
      resultoffe[j] = datoffe[3];
      for(int i=0;i<5;i++) offw_result_par >> datoffw[i]; 
      resultoffw[j] = datoffw[3];

    }

    for(int i=0;i<16;i++) resultzdepe[j][i] = 0.0;
    for(int i=0;i<16;i++) resultzdepw[j][i] = 0.0;

    if (flag>=4) {

      for(int i=0;i<9;i++) zdepe_result >> datzdepe[i];
      for(int i=0;i<8;i++) zdepe_result >> datzdepe[i+9];
      for(int i=0;i<8;i++) zdepe_result >> datzdepe[i+17];
      for(int i=0;i<8;i++) zdepe_result >> datzdepe[i+25];
      for(int i=0;i<16;i++) resultzdepe[j][i] = datzdepe[1+i*2];

      for(int i=0;i<9;i++) zdepw_result >> datzdepw[i];
      for(int i=0;i<8;i++) zdepw_result >> datzdepw[i+9];
      for(int i=0;i<8;i++) zdepw_result >> datzdepw[i+17];
      for(int i=0;i<8;i++) zdepw_result >> datzdepw[i+25];
      for(int i=0;i<16;i++) resultzdepw[j][i] = datzdepw[1+i*2];
      
    }

    t0parr[j][0] = 0.100+(L/2)*0.0004;
    t0parr[j][1] = -0.0004; 
    t0parr[j+216][0] = 0.100+(L/2)*0.0004;
    t0parr[j+216][1] = 0.0004;

    if (flag>=5) {

      for(int i=0;i<5;i++) rese_result >> datrese[i];
      resultrese[j][0] = datrese[1]; resultrese[j][1] = datrese[3]; 
      for(int i=0;i<5;i++) resw_result >> datresw[i];
      resultresw[j][0] = datresw[1]; resultresw[j][1] = datresw[3]; 
      
      t0parr[j][0] = resultrese[j][0]+(L/2)*resultrese[j][1]; 
      t0parr[j][1] = -resultrese[j][1]; 
      t0parr[j+216][0] = resultresw[j][0]+(L/2)*resultresw[j][1]; 
      t0parr[j+216][1] = resultresw[j][1]; 
      
    }

    t0parc[j][0] = resultoffe[j] - resultzdepe[j][0] - resultzdepe[j][1]*L/2 - resultzdepe[j][2]*L*L/4 -
      resultzdepe[j][3]*L*L*L/8 - L/(2.*result[j][0]);
    t0parc[j][1] = + 1./result[j][0] - resultzdepe[j][1] - resultzdepe[j][2]*L - (3./4.)*resultzdepe[j][3]*L*L; 
    t0parc[j][2] = - resulte[j][0];
    t0parc[j][3] = - resultzdepe[j][2] - (3./2.)*L*resultzdepe[j][3];
    t0parc[j][4] = 0.;//- resulte[j][1]; 
    t0parc[j][5] = 0.;//- resulte[j][2];
    t0parc[j][6] = - resultzdepe[j][3];
    t0parc[j][7] = 0.;//- resulte[j][3];
    for (int y = 0; y<12; y++) t0parc[j][7+y+1] = - resultzdepe[j][4+y];

    t0parc[j+216][0] = resultoffw[j] - resultzdepw[j][0] - resultzdepw[j][1]*L/2 - resultzdepw[j][2]*L*L/4 -
      resultzdepw[j][3]*L*L*L/8 - L/(2.*result[j][0]);
    t0parc[j+216][1] = - 1./result[j][0] - resultzdepw[j][1] - resultzdepw[j][2]*L - (3./4.)*resultzdepw[j][3]*L*L; 
    t0parc[j+216][2] = - resultw[j][0];
    t0parc[j+216][3] = - resultzdepw[j][2] - (3./2.)*L*resultzdepw[j][3];
    t0parc[j+216][4] = 0.;//- resultw[j][1]; 
    t0parc[j+216][5] = 0.;//- resultw[j][2];
    t0parc[j+216][6] = - resultzdepw[j][3];
    t0parc[j+216][7] = 0.;//- resultw[j][3];
    for (int y = 0; y<12; y++) t0parc[j+216][7+y+1] = - resultzdepw[j][4+y];

    for(int i=0;i<10;i++) atte_result >> t0parq[j][i];
    for(int i=0;i<10;i++) attw_result >> t0parq[j+216][i];

  }

  twe_result_par.close();
  tww_result_par.close();
  speed_result_par.close();
  offe_result_par.close();
  offw_result_par.close();
  zdepe_result.close(); 
  zdepw_result.close();
  rese_result.close();
  resw_result.close();
  atte_result.close();
  attw_result.close();

}

double TimeCorr(double traw, double z, double q, int channel) {

  double zs = L/2. + z;
//   double zs = z;

  double fq = 1./sqrt(q);
  
  double tcorr =  -1.*traw;
  tcorr += parc[channel][0] + parc[channel][1]*z + parc[channel][2]*fq  + 
	  parc[channel][3]*z*z + parc[channel][4]*fq*fq + parc[channel][5]*z*fq + 
	  parc[channel][6]*z*z*z + parc[channel][7]*z*fq*fq +
          (parc[channel][8] +parc[channel][9]*zs +parc[channel][10]*zs*zs+parc[channel][11]*zs*zs*zs)*fq      +
          (parc[channel][12]+parc[channel][13]*zs+parc[channel][14]*zs*zs+parc[channel][15]*zs*zs*zs)*fq*fq   +
          (parc[channel][16]+parc[channel][17]*zs+parc[channel][18]*zs*zs+parc[channel][19]*zs*zs*zs)*fq*fq*fq;
//            parc[channel][8] +parc[channel][13]/sqrt(q) + 
//           (parc[channel][9] +parc[channel][14]/sqrt(q))*exp(-(parc[channel][10]+parc[channel][15]/sqrt(q))*zs) +
//           (parc[channel][11]+parc[channel][16]/sqrt(q))*exp(-(parc[channel][12]+parc[channel][17]/sqrt(q))*zs);
//  if (channel>215) cout << channel-216 << " " << tcorr << " " << parc[channel][8] << endl;
  return tcorr;

}

double TimeCorrT0(double traw, double z, double q, int channel) {

  double zs = L/2. + z;
//   double zs = z;

  double fq = 1./sqrt(q);

  double tcorr =  -1.*traw;
  tcorr += t0parc[channel][0] + t0parc[channel][1]*z + t0parc[channel][2]*fq  + 
	  t0parc[channel][3]*z*z + t0parc[channel][4]*fq*fq + t0parc[channel][5]*z*fq + 
	  t0parc[channel][6]*z*z*z + t0parc[channel][7]*z*fq*fq +
          (t0parc[channel][8] +t0parc[channel][9]*zs +t0parc[channel][10]*zs*zs+t0parc[channel][11]*zs*zs*zs)*fq      +
          (t0parc[channel][12]+t0parc[channel][13]*zs+t0parc[channel][14]*zs*zs+t0parc[channel][15]*zs*zs*zs)*fq*fq   +
          (t0parc[channel][16]+t0parc[channel][17]*zs+t0parc[channel][18]*zs*zs+t0parc[channel][19]*zs*zs*zs)*fq*fq*fq;
//            t0parc[channel][8] +t0parc[channel][13]/sqrt(q) + 
//           (t0parc[channel][9] +t0parc[channel][14]/sqrt(q))*exp(-(t0parc[channel][10]+t0parc[channel][15]/sqrt(q))*zs) +
//           (t0parc[channel][11]+t0parc[channel][16]/sqrt(q))*exp(-(t0parc[channel][12]+t0parc[channel][17]/sqrt(q))*zs);

//   int zbin = zs/30.;
//   if (zbin<0 && zbin>9) cout << "???" << endl;
//   else {
//     double cres = binpar[zbin][0] + binpar[zbin][1]*fq + binpar[zbin][2]*fq*fq + binpar[zbin][3]*fq*fq*fq;
//     tcorr -= cres;
//   }

  return tcorr;

}

// Select track as in standard tof-track match
double Chi2(double te, double tw, double qe, double qw, double zi, int br, double pl) {
//   double fqe = sqrt(1/qe);  
  double ctme = TimeCorr(te, zi, qe, br);
  double rte = parr[br][0] + parr[br][1]*zi;
  
//   double fqw = sqrt(1/qw); 
  double ctmw = TimeCorr(tw, zi, qw, br+216);
  double rtw = parr[br+216][0] + parr[br+216][1]*zi;
//   cout << ctme << " " << ctmw << " " << rte << " " << rtw << endl;
  double tdr = sqrt(rte*rte + rtw*rtw);
  double c2t = pow( (ctme - ctmw)/tdr , 2);

  double AF[2];
  AF[0] = exp( -(L/2.- zi)/parq[br][3]     + parq[br][5]*pow(L/2.- zi,2)      );
  AF[1] = exp( -(L/2.+ zi)/parq[br+216][3] + parq[br+216][5]*pow(L/2.+ zi,2)  );
  double mpe = parq[br][1]*AF[0];
  double mpw = parq[br+216][1]*AF[1];

  double qd = (qe/mpe - qw/mpw)/parq[br][9];
  qd *= 4.0/pl;
  double c2q = pow( qd , 2);

  return c2t + c2q;

}

double ParTable[432][40];
void ReadParTable(TString ParTableName) {

  ifstream ParTableFile;
  ParTableFile.open(ParTableName);

  for (int c = 0; c<432; c++)
    for (int p = 0; p<41; p++) {
      int channel;
      if (p==0) ParTableFile >> channel;
      else ParTableFile >> ParTable[c][p-1];
    }

}

double CorrectedTime(double t, double z, double q, int channel) {
    
  double fq;
  if   (q <= 0) fq = 0;
  else fq = sqrt(1/q);

  double par[40]; 
  for (int p = 0; p<40; p++) par[p] = ParTable[channel][p];
  
  double t_cor = -t;
  if(par[0])
      { 
	t_cor += par[1] + par[2]*z + par[3]*fq
	  + par[4]*z*z + par[5]*fq*fq + par[6]*z*fq + par[7]*z*z*z +
	  ( par[28] + par[29]*(L/2 + z)+ par[30]*(L/2 + z)*(L/2 + z)+ par[31]*(L/2 + z)*(L/2 + z)*(L/2 + z) )*fq + 
	  ( par[32] + par[33]*(L/2 + z)+ par[34]*(L/2 + z)*(L/2 + z)+ par[35]*(L/2 + z)*(L/2 + z)*(L/2 + z) )*fq*fq + 
	  ( par[36] + par[37]*(L/2 + z)+ par[38]*(L/2 + z)*(L/2 + z)+ par[39]*(L/2 + z)*(L/2 + z)*(L/2 + z) )*fq*fq*fq;
      } 

  return t_cor;
}

//--------------------------------------------------------------------------
// Event counter
//--------------------------------------------------------------------------

void speed_calib::EventCounter()
{
  if (fChain == 0) return;
  Int_t nentries = Int_t(fChain->GetEntries());
  Int_t nbytes = 0, nb = 0;

  cout << "\n LoopSpd: Number of input tracks = " << nentries << " \n" << endl;

  int currentEvent = -999;
  int eventCounter =    0;
  int kentry	   =    0;

  for (Int_t jentry=0; jentry<nentries; jentry++) {
    Int_t ientry = LoadTree(jentry);
    nb = fChain->GetEntry(jentry);   
    nbytes += nb;
    
    if ((jentry%trkPrescale) < 10)
    {
     // Prescaled track count: kentry
     ++kentry;

    //if (jentry%10000 == 0)cout <<  lightspeed_d0wrt << ", "<<  lightspeed_inBar << ", " 
    //                          << lightspeed_run << ", " <<  lightspeed_event << endl;
    if (kentry%200000 == 0) printf(" EventCounter:track %8d (run %d) \n",kentry,lightspeed_run);

    if (lightspeed_event != currentEvent) {
      currentEvent = lightspeed_event;
      ++eventCounter;
    }

    } // trkPrescale
    
  } // for..jentry

  cout << " Number of input events = " << eventCounter << endl;

}

//--------------------------------------------------------------------------
// Speed of light 
//--------------------------------------------------------------------------
       
void speed_calib::LoopSpd(int     NumEvt,
                          int     flag,
                          double  cut,
                          char   *file_name,
                          char   *directory)
{
  if (fChain == 0) return;
  Int_t nentries = Int_t(fChain->GetEntries());
  Int_t nbytes = 0, nb = 0;

  cout << "\n LoopSpd: Number of input tracks = " << nentries << " \n" << endl;

  if (flag < 10) {sprintf(filesl,"%s/root/%s_speed-calib.root",directory,file_name);} 
  else           {sprintf(filesl,"%s/root/%s_spdtw-calib.root",directory,file_name);}

  TFile *hfile = new TFile(filesl,"recreate");

  gStyle->SetOptFit(111);

  TH1F *speedm    = new TH1F("speedm",   "speedm",   100, 10,20  ); // z vs deltaT
  TH1F *errspeedm = new TH1F("errspeedm","errspeedm",100,  0, 0.2);
  TH1F *residm    = new TH1F("residm",   "residm",   100, -1, 1  );
  TH1F *chim      = new TH1F("chim",     "chim",     100,  0,10  );

  TH2F *speed [NumBar]; // z vs deltaT
  TH2F *speed_bef_cuts [NumBar]; // z vs deltaT
  TH1F *resid [NumBar]; // residuals
  TH2F *residz[NumBar]; // residuals (deltaT - 2z/v - off) vs z 
  TH2F *residt[NumBar]; // residuals (deltaT - 2z/v - off) vs deltaT    

  TH1F *tracks_before = new TH1F("tracks_before", "tracks_before", 216,0,215);
  TH1F *tracks_after  = new TH1F("tracks_after",  "tracks_after",  216,0,215);

  for (int i=0; i<NumBar; i++) {
    sprintf(titlesl,"speed_%d",i);
    speed[i]  = new TH2F(titlesl,titlesl,150,-150,150,150,-20,20);
    sprintf(titlesl,"speed_bef_cuts_%d",i);
    speed_bef_cuts[i]  = new TH2F(titlesl,titlesl,150,-150,150,150,-20,20);
    sprintf(titlesl,"resid_%d",i);
    resid[i]  = new TH1F(titlesl,titlesl,1000,-20,20);
    sprintf(titlesl,"resid_vs_z_%d",i);
    residz[i] = new TH2F(titlesl,titlesl,300,-150,150,300,-15,15); 
    sprintf(titlesl,"resid_vs_deltaT_%d",i);
    residt[i] = new TH2F(titlesl,titlesl,300,-15,15,300,-15,15); 
  }                     

  int cont[NumBar]; for (int i=0; i<NumBar; i++) cont[i] = 0;

  double data[9], result[NumBar][2];
  //double twe[NumBar], tww[NumBar], datae[9], dataw[9];
 
  for (int i=0; i<9; i++) {
    data [i] = 0.0;
    //datae[i] = 0.0;
    //dataw[i] = 0.0;
  }

  for (int i=0; i<NumBar; i++) {
    //twe[i] = 0.0;
    //tww[i] = 0.0;
    for (int j=0; j<2; j++) result[i][j] = 0.0;
  }

  ifstream speed_result, twe_result, tww_result;
  ofstream speed_calib;

  //------------------------------------------------------------------------
  // Open the speed of light and time walk calibration files
  // flag <= 0       --> No input file
  // 0 < flag < 10   --> Speed of light input file
  // 10 < flag < 100 --> Speed of light + TW input files
  // flag > 100      --> Speed of light after TW + TW input files
  //------------------------------------------------------------------------

  if (flag > 0) { 
    if (flag < 100) {
      sprintf(filesl,"%s/files/%s_speed-calib.txt",directory,file_name);
      speed_result.open(filesl,ios::in);
    } else {
      sprintf(filesl,"%s/files/%s_spdtw-calib.txt",directory,file_name);
      speed_result.open(filesl,ios::in);
    }

    if (flag > 10) {
      FillParameters(directory, file_name, 1);
//       sprintf(filesl,"%s/files/%s_timew-calib_e.txt",directory,file_name);
//       twe_result.open(filesl,ios::in);   
//       sprintf(filesl,"%s/files/%s_timew-calib_w.txt",directory,file_name);
//       tww_result.open(filesl,ios::in);
    }

    for(int j=0; j<NumBar; j++) {
      for(int i=0; i<9; i++) speed_result >> data[i];
      
      result[j][0] = data[1]; // Is saved the speed : 2/fit value
      result[j][1] = data[3];

//       if (flag > 10) {
//         for(int i=0; i<9; i++) {twe_result >> datae[i]; tww_result >> dataw[i];}
//         twe[j] = datae[1];
//         tww[j] = dataw[1];
//       }  
    }

    speed_result.close();

    if (flag > 10) {twe_result.close(); tww_result.close();} 
  }

  Int_t kentry = 0;

  // Loop over tracks
  if (NumEvt > nentries) NumEvt = nentries;
  for (Int_t jentry=0; jentry<NumEvt; jentry++) {
    Int_t ientry = LoadTree(jentry);
    nb = fChain->GetEntry(jentry);   
    nbytes += nb;  

    if ((jentry%trkPrescale) < 10)
    {
    // Prescaled track count: kentry
    ++kentry;
    
    int bar_number = int(lightspeed_bar);

    //if (kentry%200000 == 0) printf(" LoopSpd:track %8d (run %d) \n",kentry,lightspeed_run);

    tracks_before->Fill(bar_number);
     
    // Track selection
    if (   //(fabs(lightspeed_d0wrt) <  1.5) 
         (lightspeed_track_type  !=   2)         
        && (lightspeed_pt          >  0.4)         
        && (lightspeed_SiPhiHits   >    2)         
        && (lightspeed_cotAxHits   >   10)      
        && (lightspeed_cotStHits   >   10)      
        && (lightspeed_sl8Hits     >    4)      
        && (lightspeed_inBar       ==   1)
       ) {  

      tracks_after->Fill(bar_number);

      double diff = -999;

      // Bar number i     => east
      // Bar number i+216 => west

      if (flag > 10 && lightspeed_qe > 0 && lightspeed_qw > 0) {
	double tce = TimeCorr(lightspeed_te, lightspeed_z, lightspeed_qe, lightspeed_bar);
	double tcw = TimeCorr(lightspeed_tw, lightspeed_z, lightspeed_qw, lightspeed_bar+216);
	diff = -(tce - tcw);
//         diff = lightspeed_te + (twe[bar_number]/sqrt(lightspeed_qe)) -
//                lightspeed_tw - (tww[bar_number]/sqrt(lightspeed_qw));
      } else { 
        diff = lightspeed_te - lightspeed_tw; // When no TW info available
      }

      float tprim = 0.0;

      if (result[bar_number][0] > 0) 
        tprim = 2.0*lightspeed_z/result[bar_number][0]+result[bar_number][1];

      speed_bef_cuts[bar_number]->Fill(lightspeed_z,diff);
       // deltaT - (2z/v + offset) < 3.5 ns
      if (fabs(diff-tprim) < cut) speed[bar_number]->Fill(lightspeed_z,diff);

      resid [bar_number]->Fill(             diff-tprim);
      residz[bar_number]->Fill(lightspeed_z,diff-tprim);
      residt[bar_number]->Fill(diff,        diff-tprim);

      ++cont[bar_number];

    } // track selection

    } // trkPrescale
  } // loop over tracks
  
  printf(" LoopSpd: Tracks used %8d (run %d) \n",kentry,lightspeed_run);
  
  // Perform linear fits
  cout << "\n Performing linear fits for speed of light... \n";

  if (flag < 10) { 
    sprintf(filesl,"%s/files/%s_speed-calib.txt",directory,file_name);
    speed_calib.open(filesl);
  } else {
    sprintf(filesl,"%s/files/%s_spdtw-calib.txt",directory,file_name);
    speed_calib.open(filesl);
  }

  // Loop over bars
  for (int i=0; i<NumBar; i++) {

    double csq    = 999;
    double err    = 999;
    double errsig = 999;
    double sigma  = 999;
    double cov[4], par[4]; for (int j=0; j<4; j++) {cov[j] = par[j] = 999;}

    if (cont[i] > 25) {
      ajuste.gaus1(resid[i],-0.75,0.75,err,sigma,errsig,csq);
      par[0] = ajuste.fitlin(speed[i],-100,100,par[1],par[2],par[3],csq,cov[0],cov[1],cov[2]);
    } else {
      cout << " Warning! In bar " << i << " there are only " << cont[i] << " tracks" << endl;
    }

    residm->Fill(sigma); 
    speed_calib << setw(13) << i 
                << setw(13) << 2.0/par[0]
                << setw(13) << 2.0*par[1]/(par[0]*par[0])
                << setw(13) << par[2]
                << setw(13) << par[3]   
                << setw(13) << csq
                << setw(13) << cov[0]
                << setw(13) << cov[1]
                << setw(13) << cov[2] << endl; 

    if(par[0] > 0.0) {
      speedm->Fill(2.0/par[0]);
      errspeedm->Fill(par[1]);
      chim->Fill(csq);
    }
  } // loop over bars

  speed_calib.close();
  hfile->Write();
  hfile->Close();     
                                 
} // speed_calib::LoopSpd

//--------------------------------------------------------------------------
// Channel to channel offsets 
//--------------------------------------------------------------------------
       
void speed_calib::LoopT0Abs(int     NumEvt,
                            double  cut,
                            char   *file_name,
                            char   *directory)
{
  if (fChain == 0) return;
  Int_t nentries = Int_t(fChain->GetEntries());
  Int_t nbytes = 0, nb = 0;

  cout << "\n LoopT0Abs: Number of input tracks = " << nentries << " \n" << endl;

  sprintf(filesl,"%s/root/%s_t0-calib.root",directory,file_name); 

  TFile *hfile = new TFile(filesl,"recreate");

  gStyle->SetOptFit(111);

  TH1F *cabe   = new TH1F("cabe",  "T0_offe", 200,40,80);
  TH1F *cabw   = new TH1F("cabw",  "T0_offe", 200,40,80);
  TH1F *cabe2  = new TH1F("cabe2", "T0_offe2",200,40,80);
  TH1F *cabw2  = new TH1F("cabw2", "T0_offe2",200,40,80);
  TH1F *speede = new TH1F("speede","speede",   40,13,17);
  TH1F *speedw = new TH1F("speedw","speedw",   40,13,17);
  TH2F *parb [NumBar];
  TH2F *parb1[NumBar];
  TH1F *med  [NumBar];

  // New plots for tE+tW check 022403
  TH1F *sum_check_tE[NumBar];
  TH1F *sum_check_tW[NumBar];
  TH1F *sum_check_qE[NumBar];
  TH1F *sum_check_qW[NumBar];
  TH1F *sum_check_tof_pi = new TH1F("sum_check_tof_pi","sum_check_tof_pi",200,0, 20);
  TH1F *sum_check_p      = new TH1F("sum_check_p",     "sum_check_p",     200,0, 20);
  TH1F *sum_check_length = new TH1F("sum_check_length","sum_check_length",200,0,300);

  for (int i=0; i<NumBar; i++) {
    sprintf(titlesl,"diff_bar_%d",i);
    parb [i] = new TH2F(titlesl,titlesl,150,-150,150,200,-20, 20);
    sprintf(titlesl,"sum_bar_%d",i);
    parb1[i] = new TH2F(titlesl,titlesl,150,-150,150,250, 70,120); 
    sprintf(titlesl,"sum_med_%d",i);
    med  [i] = new TH1F(titlesl,titlesl,100,  50,150);

    sprintf(titlesl,"sum_check_tE_%d",i);  
    sum_check_tE[i] = new TH1F(titlesl,titlesl,200,0,60);
    sprintf(titlesl,"sum_check_tW_%d",i);  
    sum_check_tW[i] = new TH1F(titlesl,titlesl,200,0,60);
    sprintf(titlesl,"sum_check_qE_%d",i);  
    sum_check_qE[i] = new TH1F(titlesl,titlesl,200,0,10000);
    sprintf(titlesl,"sum_check_qW_%d",i);  
    sum_check_qW[i] = new TH1F(titlesl,titlesl,200,0,10000);
  }

  FillParameters(directory, file_name, 1);

  double diff[NumBar], sum[NumBar];

  double data[9];//, datae[9], dataw[9];
  double result[NumBar][2];//, resulte[NumBar], resultw[NumBar];

  for (int i=0; i<9; i++) {
    data [i] = 0.0;
//     datae[i] = 0.0;
//     dataw[i] = 0.0;
  }

  for (int i=0; i<NumBar; i++) {
//     resulte[i] = 0.0;
//     resultw[i] = 0.0;
    for (int j=0; j<2; j++) result[i][j] = 0.0;
  }

//   ifstream speed_result, twe_result, tww_result, ofe_result, ofw_result; 
  ifstream speed_result, ofe_result, ofw_result; 
  ofstream ofe_calib, ofw_calib; 

  sprintf(filesl,"%s/files/%s_spdtw-calib.txt",directory,file_name);
  speed_result.open(filesl,ios::in);
//   sprintf(filesl,"%s/files/%s_timew-calib_e.txt",directory,file_name);
//   twe_result.open(filesl,ios::in);
//   sprintf(filesl,"%s/files/%s_timew-calib_w.txt",directory,file_name);
//   tww_result.open(filesl,ios::in);

  for(int j=0; j<NumBar; j++) {
    for(int i=0; i<9; i++) {
      speed_result >> data [i];
//       twe_result   >> datae[i];
//       tww_result   >> dataw[i];
    }

    result [j][0] = data [1];
    result [j][1] = data [3];
//     resulte[j]    = datae[1];
//     resultw[j]    = dataw[1];
  }

  speed_result.close();
//   twe_result.close();
//   tww_result.close();
 
  for(int i=0; i<NumBar; i++) {
    sum [i] =    0;
    diff[i] = -999;
  }

  Int_t kentry = 0;
  
  // Loop over tracks
  if (NumEvt > nentries) NumEvt = nentries;
  for (Int_t jentry=0; jentry<NumEvt; jentry++) {
    Int_t ientry = LoadTree(jentry);
    nb = fChain->GetEntry(jentry);
    nbytes += nb;

    if ((jentry%trkPrescale)<10)
    {
    // Prescaled track count: kentry
    ++kentry;
 
    int bar_number = int(lightspeed_bar);

    //if (kentry%200000 == 0) printf(" LoopT0Abs:track %8d (run %d) \n",kentry,lightspeed_run);

    // Track selection
    if (   //(fabs(lightspeed_d0wrt) <  1.5) 
         (lightspeed_track_type  !=   2)         
        && (lightspeed_pt          >  0.4)         
        && (lightspeed_SiPhiHits   >    2)      
        && (lightspeed_cotAxHits   >   10)      
        && (lightspeed_cotStHits   >   10)      
        && (lightspeed_sl8Hits     >    4)
        && (lightspeed_inBar       ==   1)    
       ) {        

      if (lightspeed_qe > 0  && lightspeed_qw > 0) {
	double tce = TimeCorr(lightspeed_te, lightspeed_z, lightspeed_qe, lightspeed_bar);
	double tcw = TimeCorr(lightspeed_tw, lightspeed_z, lightspeed_qw, lightspeed_bar+216);
	diff[bar_number] = -(tce - tcw);
//         diff[bar_number] = lightspeed_te + resulte[bar_number]/sqrt(lightspeed_qe) -
//                            lightspeed_tw - resultw[bar_number]/sqrt(lightspeed_qw);
      } //else { 
  //      diff[bar_number] = lightspeed_te - lightspeed_tw; // When no TW info available
  //    }

//       tof_pi = (lightspeed_p/pimass)/sqrt(1+((lightspeed_p/pimass)*(lightspeed_p/pimass))); 
//       tof_pi = lightspeed_tofLength/(tof_pi*c);
      tof_pi = getTofSL(lightspeed_tofLength, lightspeed_p, pimass);

      sum_check_p     ->Fill(lightspeed_p);
      sum_check_length->Fill(lightspeed_tofLength);
      sum_check_tof_pi->Fill(tof_pi);
      sum_check_tE[bar_number]->Fill(lightspeed_te);
      sum_check_tW[bar_number]->Fill(lightspeed_tw);
      sum_check_qE[bar_number]->Fill(lightspeed_qe);
      sum_check_qW[bar_number]->Fill(lightspeed_qw);

      if (lightspeed_qe > 0 && lightspeed_qw > 0) {
	double tce = TimeCorr(lightspeed_te, lightspeed_z, lightspeed_qe, lightspeed_bar);
	double tcw = TimeCorr(lightspeed_tw, lightspeed_z, lightspeed_qw, lightspeed_bar+216);
	sum[bar_number] = -(tce + tcw) + 2.0*tof_pi;
//         sum[bar_number] = lightspeed_te + resulte[bar_number]/sqrt(lightspeed_qe) +
//                           lightspeed_tw + resultw[bar_number]/sqrt(lightspeed_qw) +
//                           2.0*tof_pi;
                        //2.0*lightspeed_tof_pi;
      } //else { 
   //     sum[bar_number] = lightspeed_te + 
   //                       lightspeed_tw +
   //                       2.0*tof_pi; 
   //                     //2.0*lightspeed_tof_pi; // When no TW info available
   //   }

      if (fabs(diff[bar_number]-2.0*lightspeed_z/result[bar_number][0]-result[bar_number][1]) < cut) 
        parb[bar_number]->Fill(lightspeed_z,diff[bar_number]);

      parb1[bar_number]->Fill(lightspeed_z,sum[bar_number]);
      med  [bar_number]->Fill(sum[bar_number]);

      sum [bar_number] =    0;
      diff[bar_number] = -999;

    } // track selection

    } // trkPrescale
    
  } // loop over tracks
  
  printf(" LoopT0Abs: Tracks used %8d (run %d) \n",kentry,lightspeed_run);

  cout << "\n Performing linear fits for t0 offsets... \n";

  sprintf(filesl,"%s/files/%s_T0-calib_e.txt",directory,file_name); ofe_calib.open(filesl);
  sprintf(filesl,"%s/files/%s_T0-calib_w.txt",directory,file_name); ofw_calib.open(filesl); 

  // Loop over bars
  for (int i=0; i<NumBar; i++) {

    double csq    = 999;
    double err    = 999;
    double errsig = 999;
    double sigma  = 999;
    double summe  = 999;

    double alpha    = 999;
    double offalpha = 999;
    double beta     = 999;
    double offbeta  = 999;

    double speedtmp = 999;
    double cabtmp   = 999;
    double cabtmp2  = 999;

    beta    = ajuste.fit(parb1[i],-100,100,alpha   ); // linear fit to the time sum vs z
    offbeta = ajuste.fit(parb [i],-100,100,offalpha); // linear fit to the time dif vz z
    summe   = ajuste.gaus1(med[i],50,150,err,sigma,errsig,csq); // average te+tw
 
    if (beta != 0 && offbeta != 0) {
      speedtmp = 2.0/(beta+offbeta);
      speede->Fill(speedtmp);
      if (speedtmp > 0) {
        cabtmp  = (alpha+offalpha)/2.0 + L/(2.0*speedtmp);
        //cabtmp2 = (summe+(L/(2.0*result[i][0]))+offalpha)/2.0;
        cabtmp2 = (summe+result[i][1])/2.0 + (L/2.0)*(1./result[i][0]);
        cabe ->Fill(cabtmp );     
        cabe2->Fill(cabtmp2);     
      }

      ofe_calib << setw(13) << i
                << setw(13) << cabtmp
                << setw(13) << speedtmp
                << setw(13) << cabtmp2
                << setw(13) << 2.0/offbeta
                << endl;

      speedtmp = 2.0/(offbeta-beta);
      speedw->Fill(speedtmp);
      if(speedtmp > 0) {
        cabtmp  = (alpha-offalpha)/2.0 + L/(2.0*speedtmp);
        //cabtmp2 = (summe+(L/(2.0*result[i][0]))-offalpha)/2.0;
        cabtmp2 = (summe-result[i][1])/2.0 + (L/2.0)*(1./result[i][0]);
        cabw ->Fill(cabtmp );     
        cabw2->Fill(cabtmp2);
      }

      ofw_calib << setw(13) << i+216
                << setw(13) << cabtmp
                << setw(13) << speedtmp
                << setw(13) << cabtmp2
                << setw(13) << 2.0/offbeta
                << endl;

    } else {

      ofe_calib << setw(13) << i
                << setw(13) << 0.0
                << setw(13) << 0.0
                << setw(13) << 0.0
                << setw(13) << 0.0
                << endl;

      ofw_calib << setw(13) << i+216
                << setw(13) << 0.0
                << setw(13) << 0.0
                << setw(13) << 0.0
                << setw(13) << 0.0
                << endl;

    } // end if beta != 0
  } // end loop over bars

  hfile->Write();                                        
  hfile->Close();
  ofe_calib.close(); 
  ofw_calib.close();

} // speed_calib::LoopT0Abs

//--------------------------------------------------------------------------
// Attenuation length
//--------------------------------------------------------------------------

void speed_calib::LoopAtt(int     NumEvt,
                          double  cut,
                          char   *file_name,
                          char   *directory)
{
  Int_t nentries = Int_t(fChain->GetEntries());
  Int_t nbytes = 0, nb = 0; 

  sprintf(filesl,"%s/root/%s_att-calib.root",directory,file_name); 
  TFile *hfile = new TFile(filesl,"recreate","");

  // Booking the histograms
  gStyle->SetOptFit(111);
  TH1F *attlenm        = new TH1F("attlenm",       "attlen cut",      120,150.,450.);
  TH1F *attlenm2       = new TH1F("attlenm2",   "attlenm ",   120,150.,450.);

  TH1F *resid [NumBar];
  TH2F *residz[NumBar];
  TH2F *attlen[NumBar];
  TH1F *resid2 [NumBar];
  TH2F *residz2[NumBar];
  TH2F *attlen2[NumBar];
 
  int    cont[NumBar]; 
  double data[15], result[NumBar][2],result2[NumBar][2];
  double Adc_e[NumBar], Adc_w[NumBar], Z[NumBar];
  double att;
//   float var_att2[maxevt][NumBar],var_zpos2[maxevt][NumBar];
  int suc[NumBar];//,suc2[NumBar]; 
  ofstream att_calib, att_calib2;
  ifstream att_result,att_result2;

  sprintf(filesl,"%s/files/%s_Att-calib.txt",directory,file_name); 
  att_result.open(filesl,ios::in);
  sprintf(filesl,"%s/files/%s_Att-calib_s2.txt",directory,file_name); 
  att_result2.open(filesl,ios::in);

  for(int j=0; j<NumBar; j++)
  {
    for(int i=0; i<15; i++) att_result >> data[i];
    result[j][0] = data[1];
    result[j][1] = data[3];
    
    for(int i=0; i<15; i++) att_result2 >> data[i];
    result2[j][0] = data[1];
    result2[j][1] = data[3];
    cout<<result[j][0]<<", "<<result[j][1]<<", "<<result2[j][0]<<", "<<result2[j][1]<<endl;
  }
  att_result.close();

  // Initialize histograms 
  for (int i=0; i<NumBar; i++)
  {
    sprintf(titlesl,"resid_bar_%d",i);
    resid[i] = new TH1F(titlesl,titlesl,100,-100.,100.);
    sprintf(titlesl,"attlen_%d",i);
    attlen[i] = new TH2F(titlesl,titlesl,600,-150.,150.,150,-3.,3.);
    sprintf(titlesl,"attlen_z_%d",i);
    residz[i] = new TH2F(titlesl,titlesl,600,-150.,150.,100,-100.,100.);
    sprintf(titlesl,"resid2_bar_%d",i);
    resid2[i] = new TH1F(titlesl,titlesl,100,-100.,100.);
    sprintf(titlesl,"attlen2_%d",i);
    attlen2[i] = new TH2F(titlesl,titlesl,600,-150.,150.,150,-3.,3.);
    sprintf(titlesl,"attlen2_z_%d",i);
    residz2[i] = new TH2F(titlesl,titlesl,600,-150.,150.,100,-100.,100.);
  }
 
  TMinuit *gMinuit = new TMinuit(7);
  gMinuit->SetFCN(fcn_att);
  double arglis[10];
  int ierr;
  arglis[0] = 1;
  gMinuit->mnexcm("SET ERR",arglis,1,ierr);
  arglis[0] = -1;
  gMinuit->mnexcm("SET PRINT",arglis,1,ierr);
//   gMinuit->mnexcm("SET SETNOWARNINGS",arglis,1,ierr);

  // Loop over tracks
  for(int i=0; i<NumBar; i++) {
    suc[i] = 0;
    //suc2[i] = 0;
  }

  int current_event = -1;
  int kentry = 0;

  if (NumEvt > nentries) NumEvt = nentries;
  for (Int_t jentry=0; jentry<NumEvt; jentry++) {
    Int_t ientry = LoadTree(jentry);
    nb = fChain->GetEntry(jentry);
    nbytes += nb;

    if ((jentry%trkPrescale)<10)
    {
    // Prescaled track count: kentry
    ++kentry;
     
    //if (kentry%200000 == 0) printf("LoopAtt:track %8d (run %d) \n",kentry,lightspeed_run);

    int bar_number = int(lightspeed_bar);

    if (current_event==-1) current_event = lightspeed_event;
    if (lightspeed_event != current_event) {
      current_event = lightspeed_event;

      // Loop over bars
      for (int i=0; i<NumBar; i++) {

        if(cont[i] != 1) continue;// Require 1 and only 1 hit in the bar

        (Adc_w[i] > 0) ? att = log(Adc_e[i]/Adc_w[i]) : att = -999;

	  var_att[suc[i]][i]  = att;
	  var_zpos[suc[i]][i] = Z[i]; 
	  suc[i]++;
	  
	  attlen2[i]->Fill(Z[i],att); 
	  if (fabs(Z[i])<=90.) attlen[i]->Fill(Z[i],att); 
	  
        resid2 [i]->Fill(Z[i]-(att*result2[i][0]+result2[i][1]));
        residz2[i]->Fill(Z[i],Z[i]-(att*result2[i][0]+result2[i][1])); 

	if (fabs(Z[i])>90.) continue;

        resid [i]->Fill(Z[i]-(att*result[i][0]+result[i][1]));
        residz[i]->Fill(Z[i],Z[i]-(att*result[i][0]+result[i][1]));

      } // loop over bars                  

      for (int i=0; i<NumBar; i++) {
        cont [i] =    0;
        Adc_e[i] = -999;
        Adc_w[i] = -999;   
        Z    [i] = -999;
      }

    } // if event_counter > num_events

    // Track selection
    if (   //(fabs(lightspeed_d0wrt) <  1.5) 
         (lightspeed_track_type  !=   2)         
        && (lightspeed_pt          >  0.4)         
        && (lightspeed_SiPhiHits   >    2)      
        && (lightspeed_cotAxHits   >   10)      
        && (lightspeed_cotStHits   >   10)      
        && (lightspeed_sl8Hits     >    4)
        && (lightspeed_inBar       ==   1)    
       ) {        

      ++cont[bar_number];

      Adc_e[bar_number] = lightspeed_qe;
      Adc_w[bar_number] = lightspeed_qw;
      Z    [bar_number] = lightspeed_z;
         
    } // track selection

    } // trkPrescale

  } // loop over tracks 

   printf("LoopAtt: Tracks used %8d (run %d) \n",kentry,lightspeed_run);

  // Perform the linear fits

  cout << endl << "Performing linear fits for att length..." << endl;

  sprintf(filesl,"%s/files/%s_Att-calib.txt",directory,file_name); att_calib.open(filesl);
  sprintf(filesl,"%s/files/%s_Att-calib_s2.txt",directory,file_name); att_calib2.open(filesl);

  double par[10], epar[10];
  //double err, csq, errsig, sigma;

  // Loop over bars
  for(int i=0;i<NumBar;i++) {

    //Loop for all the events
   evts=suc[i];
   for( int k =0; k<evts;k++){
     lenght[k] = var_att[k][i]; 
     zps[k] = var_zpos[k][i];}
    
    gMinuit->mncler();
    gMinuit->SetFCN(fcn_att);

    gMinuit->mnparm(0,"at",result2[i][0],0.1,0.,1000.,ierr);
    gMinuit->mnparm(1,"of",result2[i][1],0.1,0,0,ierr);
    gMinuit->mnparm(2,"rn",5.,0.001,0.,1000.,ierr);
    gMinuit->mnparm(3,"ff",.2,0.01,0.,1.,ierr);
    gMinuit->mnparm(4,"rw",100.,0.001,0.,1000.,ierr);
    gMinuit->mnparm(5,"at2",result2[i][0]/2.,0.1,0.,result2[i][0],ierr);
    gMinuit->mnparm(6,"of2",result2[i][1],0.1,0,0,ierr);
   

    arglis[0] = 1000;
    arglis[1] = 1;

    gMinuit->mnexcm("MIGRAD",arglis,1,ierr);

    if ( ierr != 0 ) { cout<<" bar " << i << " fails for: "<< ierr << endl;
    for(int k=0;k<6;k++){par[k] = 0.0;
    epar[k] = 0.0;}
    }   
     gMinuit ->GetParameter(0,par[0],epar[0]);
     gMinuit ->GetParameter(1,par[1],epar[1]);
     gMinuit ->GetParameter(2,par[2],epar[2]);
     gMinuit ->GetParameter(3,par[3],epar[3]);
     gMinuit ->GetParameter(4,par[4],epar[4]);
     gMinuit ->GetParameter(5,par[5],epar[5]);
     gMinuit ->GetParameter(6,par[6],epar[6]);


    att_calib2 << setw(13) << i        
              << setw(13) << par[0]   
              << setw(13) << epar[0]   
              << setw(13) << par[1]   
              << setw(13) << epar[1]   
              << setw(13) << par[2]     
              << setw(13) << epar[2]   
              << setw(13) << par[3]   
              << setw(13) << epar[3]   
              << setw(13) << par[4]   
              << setw(13) << epar[4]
              << setw(13) << par[5]   
              << setw(13) << epar[5]
              << setw(13) << par[6]   
              << setw(13) << epar[6]        << endl;    
    if(par[0] > 0 && par[0] < 1000000) {
      attlenm2   ->Fill(par[0]);
    }


    //Loop for events that are not at the ends of the bars
//   cout << "Performing linear fits for att length...(cutting edges)" << endl;
   evts=0;
   for( int k =0; k<suc[i];k++)
     if (fabs(var_zpos[k][i])<=90.) {
       lenght[k] = var_att[k][i]; 
       zps[k] = var_zpos[k][i];
       evts++;
     }
    
    gMinuit->mncler();
    gMinuit->SetFCN(fcn_att);

    gMinuit->mnparm(0,"at",par[0],0.1,0.,1000.,ierr);
    gMinuit->mnparm(1,"of",par[1],0.1,0,0,ierr);
    gMinuit->mnparm(2,"rn",5.,0.001,0.,100.,ierr);
    gMinuit->mnparm(3,"ff",.8,0.01,0.,1.,ierr);
    gMinuit->mnparm(4,"rw",100.,0.001,0.,1000.,ierr);
    gMinuit->mnparm(5,"at2",par[5]/2.,0.1,0.,par[5],ierr);
    gMinuit->mnparm(6,"of2",par[6],0.1,0,0,ierr);
   

    arglis[0] = 1000;
    arglis[1] = 1;

    gMinuit->mnexcm("MIGRAD",arglis,1,ierr);

    if ( ierr != 0 ) { cout<<" bar " << i << " fails for : "<< ierr
			   << " (z cut applied)" << endl;
       for(int k=0;k<6;k++){par[k] = 0.0;
          epar[k] = 0.0;}
     
    }   

     gMinuit ->GetParameter(0,par[0],epar[0]);
     gMinuit ->GetParameter(1,par[1],epar[1]);
     gMinuit ->GetParameter(2,par[2],epar[2]);
     gMinuit ->GetParameter(3,par[3],epar[3]);
     gMinuit ->GetParameter(4,par[4],epar[4]);
     gMinuit ->GetParameter(5,par[5],epar[5]);
     gMinuit ->GetParameter(6,par[6],epar[6]);

    att_calib << setw(13) << i        
              << setw(13) << par[0]   
              << setw(13) << epar[0]   
              << setw(13) << par[1]   
              << setw(13) << epar[1]   
              << setw(13) << par[2]     
              << setw(13) << epar[2]   
              << setw(13) << par[3]   
              << setw(13) << epar[3]   
              << setw(13) << par[4]   
              << setw(13) << epar[4]      
	      << setw(13) << par[5]   
              << setw(13) << epar[5]
              << setw(13) << par[6]   
              << setw(13) << epar[6]       << endl;    
    if(par[0] > 0 && par[0] < 1000000) {
      attlenm->Fill(par[0]);
    }


  } // loop over bars
   
  att_calib2.close();
  att_calib.close();   

  hfile->Write(); 
  hfile->Close();                          

} // speed_calib::LoopAtt


//--------------------------------------------------------------------------
// Attenuation lengths
//--------------------------------------------------------------------------

void speed_calib::LoopAttl(int NumEvt,
			   double cut,
			   char  *file_name,
			   char   *directory,
			   int   Flag) {
  
  cout << "Start Attenuation Length Analysis With Flag = " << Flag << endl;

  // open output file
  // ================
  if (Flag>1) 
    sprintf(filesl,"%s/root/%s_att-calib_res.root",directory,file_name); 
  else 
    sprintf(filesl,"%s/root/%s_att-calibt.root",directory,file_name); 
  TFile *hfile = new TFile(filesl,"recreate","");

  // Booking the histograms
  // ======================
  gStyle->SetOptFit(111);
  TH1F *attlenm        = new TH1F("attlenm",       "attlenm Manual",      120,150,450);
  TH1F *qmip    = new TH1F("qmip",   "qmip", 200,0.,5000.);
  TH1F *aterm     = new TH1F("aterm",    "aterm",  200,  -0.1,  0.1);
  TH1F *norma = new TH1F("norma","norma",500,  0,  4000.);
  TH1F *residm         = new TH1F("residm",        "residm",              300, -5,  5);
  TH1F *attlene[NumBar],*attlenecp[NumBar],*attlenec[NumBar], *attlenec_A[NumBar],*attlenec_N[NumBar];
  TH1F *attlenw[NumBar],*attlenwcp[NumBar],*attlenwc[NumBar],*attlenwc_A[NumBar],*attlenwc_N[NumBar];
  TH1F * attleneb1[216],* attleneb2[216],* attleneb3[216],* attleneb4[216];
  TH1F * attlenwb1[216],* attlenwb2[216],* attlenwb3[216],* attlenwb4[216];
  TF1 * lan = new TF1 ("lan","landau");
  TGraph* attbare = new TGraph(216) ;
  TGraph* attbarw = new TGraph(216) ;
  TH2F* adcreshe[216] ;
  TH2F* adcreshw[216] ;
  TH2F* adcreshec[216] ;
  TH2F* adcreshwc[216] ;
  TH2F* residuale[216];
  TH2F* residualw[216];
  TH1F* intresiduale[216];
  TH1F* intresidualw[216];
  TH2F* sigmaQ[NumBar];
  TH1F * sigmaQs[NumBar], *sigmaQsr[NumBar], *sigmaQz[NumBar][4];
  TH1F * charge,*charge_pn,*charge_p, *chargem, * pathl;

  // declare useful variables
  // ========================
  double Adc_e[NumBar], Adc_w[NumBar], Z[NumBar], PathB[NumBar];
  int    cont[NumBar]; 
  double data[15], result[NumBar][6];
  double date[216][11],datw[216][11];
  double adc1, adc2;

 for(int i=0; i<NumBar; i++) {
   nvt[i] = 0 ;
   Adc_e[i] = 0.;
   Adc_w[i] = 0.;
   Z[i] = 0.;
   PathB[i] = 0.;
   //for (int k = 0; k<11; k++) { date[i][k]; datw[i][k];}
 }
 
  // declare output-input text files
  // ===============================
  ofstream att_calib, lin_calib, Q0_calib;
  ifstream att_result;

  // Initialize histograms 
  // =====================
  TH1F * adchis = new TH1F("adchis","adchis",100,0.,5000.);

  if(Flag>1){
    sprintf(titlesl,"charge");
    charge = new TH1F(titlesl,titlesl,100,0.,2000.);
    sprintf(titlesl,"charge_pn");
    charge_pn = new TH1F(titlesl,titlesl,100,0.,6000.);
    sprintf(titlesl,"charge_p");
    charge_p = new TH1F(titlesl,titlesl,100,0.,6000.);
    sprintf(titlesl,"charge_m");
    chargem = new TH1F(titlesl,titlesl,100,0.,2000.);
    sprintf(titlesl,"path_le");
    pathl = new TH1F(titlesl,titlesl,100,0.,10.);
  }
  for (int i=0; i<NumBar; i++)
  {
    if(Flag>1){
      sprintf(titlesl,"resid_e%d",i);
      intresiduale[i] = new TH1F(titlesl,titlesl,100,-2000.,3000.);
      sprintf(titlesl,"resid_w%d",i);
      intresidualw[i] = new TH1F(titlesl,titlesl,100,-2000.,3000.);
      sprintf(titlesl,"resid_z_e%d",i);
      residuale[i] = new TH2F(titlesl,titlesl,100,-150.,150.,100,-2000.,3000.);
      sprintf(titlesl,"resid_z_w%d",i);
      residualw[i] = new TH2F(titlesl,titlesl,100,-150.,150.,100,-2000.,3000.);
      sprintf(titlesl,"attlenc_A_e%d",i);

      attlenec_A[i] = new TH1F(titlesl,titlesl,200,0.,10000.);      
      sprintf(titlesl,"attlenc_A_w%d",i);
      attlenwc_A[i] = new TH1F(titlesl,titlesl,200,0.,10000.);
      sprintf(titlesl,"attlenc_N_e%d",i);
      attlenec_N[i] = new TH1F(titlesl,titlesl,200,0.,5.);
      sprintf(titlesl,"attlenc_N_w%d",i);
      attlenwc_N[i] = new TH1F(titlesl,titlesl,200,0.,5.);
      sprintf(titlesl,"sigmaQ_z_%d",i);
      sigmaQ[i] = new TH2F(titlesl,titlesl,100,-150.,150.,150,-2.,2.);
      sprintf(titlesl,"sigmaQ_%d",i);
      sigmaQs[i] = new TH1F(titlesl,titlesl,200,-2.,2.);
      sprintf(titlesl,"sigmaQ_er_%d",i);
      sigmaQsr[i] = new TH1F(titlesl,titlesl,200,-0.005,0.005);
      for (int k = 0; k<4; k++) {
	sprintf(titlesl,"sigmaQz%d_%d",k,i);
	sigmaQz[i][k] = new TH1F(titlesl,titlesl,200,-2.,2.);
      }
    }
    sprintf(titlesl,"attlen_e%d",i);
    attlene[i] = new TH1F(titlesl,titlesl,50,0.,5000.);
    sprintf(titlesl,"attlen_w%d",i);
    attlenw[i] = new TH1F(titlesl,titlesl,50,0.,5000.);
    sprintf(titlesl,"attlencp_e%d",i);
    attlenecp[i] = new TH1F(titlesl,titlesl,100,0.,5000.);
    sprintf(titlesl,"attlencp_w%d",i);
    attlenwcp[i] = new TH1F(titlesl,titlesl,100,0.,5000.);
    sprintf(titlesl,"attlenc_e%d",i);
    attlenec[i] = new TH1F(titlesl,titlesl,100,0.,5000.);
    sprintf(titlesl,"attlenc_e_b1_%d",i);
    attleneb1[i] = new TH1F(titlesl,titlesl,100,0.,5000.);
    sprintf(titlesl,"attlenc_e_b2_%d",i);
    attleneb2[i] = new TH1F(titlesl,titlesl,100,0.,5000.);
    sprintf(titlesl,"attlenc_e_b3_%d",i);
    attleneb3[i] = new TH1F(titlesl,titlesl,100,0.,5000.);
    sprintf(titlesl,"attlenc_e_b4_%d",i);
    attleneb4[i] = new TH1F(titlesl,titlesl,100,0.,5000.);
    sprintf(titlesl,"attlenc_w%d",i);
    attlenwc[i] = new TH1F(titlesl,titlesl,100,0.,5000.);
    sprintf(titlesl,"attlenc_w_b1_%d",i);
    attlenwb1[i] = new TH1F(titlesl,titlesl,100,0.,5000.);
    sprintf(titlesl,"attlenc_w_b2_%d",i);
    attlenwb2[i] = new TH1F(titlesl,titlesl,100,0.,5000.);
    sprintf(titlesl,"attlenc_w_b3_%d",i);
    attlenwb3[i] = new TH1F(titlesl,titlesl,100,0.,5000.);
    sprintf(titlesl,"attlenc_w_b4_%d",i);
    attlenwb4[i] = new TH1F(titlesl,titlesl,100,0.,5000.);
    sprintf(titlesl,"adcres_z_e%d",i);
    adcreshe[i] = new TH2F(titlesl,titlesl,100,-150.,150.,150,0.,5000.);
    sprintf(titlesl,"adcres_z_w%d",i);
    adcreshw[i] = new TH2F(titlesl,titlesl,100,-150.,150.,150,0.,5000.);
    sprintf(titlesl,"adcresc_z_e%d",i);
    adcreshec[i] = new TH2F(titlesl,titlesl,100,-150.,150.,150,0.,5000.);
    sprintf(titlesl,"adcresc_z_w%d",i);
    adcreshwc[i] = new TH2F(titlesl,titlesl,100,-150.,150.,150,0.,5000.);
  }

  // initializing Minuit
  // ===================
  double arglis[4];
  //double ierr; // anachronism
  int ierr;
    //ierr = 0;
  TMinuit *gMinuit = new TMinuit(7);
  gMinuit->SetFCN(fcnatt);
  arglis[0] = 1;
  gMinuit->mnexcm("SET ERR",arglis,1,ierr);
  arglis[0] = -1;
  gMinuit->mnexcm("SET PRINT",arglis,1,ierr);
  gMinuit->mnexcm("SET SETNOWARNINGS",arglis,1,ierr);
 
  // read results from LoopAtt()
  // ===========================
  if (Flag==1) 
    sprintf(filesl,"%s/files/%s_Att-calib.txt",directory,file_name);
  else if (Flag==0) 
    sprintf(filesl,"%s/files/%s_Att-calib_s2.txt",directory,file_name); 
  att_result.open(filesl,ios::in);
  
  for(int j=0; j<NumBar; j++) {
    for(int i=0; i<15; i++) att_result >> data[i];
    result[j][0] = data[1];
    result[j][1] = data[3];
    //cout << result[j][0] <<", " << result[j][1] <<endl;
  }
  
  att_result.close();
  
  // Read results from LoopAttl()
  // ============================
  if(Flag>1){

    if (Flag==2)
      sprintf(filesl,"%s/files/%s_Att-calib_t_s2.txt",directory,file_name); 
    else if (Flag==3)
      sprintf(filesl,"%s/files/%s_Att-calib_t.txt",directory,file_name);
    ifstream attresult;
    attresult.open(filesl,ios::in);
   
    double date0, date1, date2, date3, date4, date5, date6, date7, date8;// date9, date10;
    for(int j=0; j<432; j++)
      {
	//if (!(j<5 || (j>215 && j<221))) continue;
	if(j<216){
	  attresult >> date0 >> date1 >> date2 
		     >> date3 >> date4 >> date5  
		    >> date6 >> date7 >> date8;
// 		     >> date9 >> date10; 
	  date[j][0] = date0;       
	  date[j][1] = date1;
	  date[j][2] = date2;
	  date[j][3] = date3;
	  date[j][4] = date4;
	  date[j][5] = date5;
	  date[j][6] = date6;
	  date[j][7] = date7;
	  date[j][8] = date8;
// 	  date[j][9] = date9;
// 	  date[j][10] = date10;
	  cout<<" corrc east = " << j << ", " << date[j][1] << ", " 
	      << date[j][3]<< ", " << date[j][5] << ", " << date[j][7] << endl;
	}
	if(j>215){
	  attresult >> date0 >> date1 >> date2 
		     >> date3 >> date4 >> date5  
		    >> date6 >> date7 >> date8;
// 		     >> date9 >> date10;  
	  datw[j-216][0] = date0;       
	  datw[j-216][1] = date1;
	  datw[j-216][2] = date2;
	  datw[j-216][3] = date3;
	  datw[j-216][4] = date4;
	  datw[j-216][5] = date5;
	  datw[j-216][6] = date6;
	  datw[j-216][7] = date7;
	  datw[j-216][8] = date8;
// 	  datw[j-216][9] = date9;
// 	  datw[j-216][10] = date10;
	  cout<<" corrc west = " << j-216 << ", " << datw[j-216][1] << ", " 
	      << datw[j-216][3]<< ", " << datw[j-216][5] << ", " 
	      << datw[j-216][7] << endl;
	}
      }

    attresult.close();
	
  }

  // Loop over tracks
  // ================
  int  current_event = -1;
  Int_t nentries = Int_t(fChain->GetEntries());
  Int_t nbytes = 0, nb = 0; 

  Int_t kentry = 0;
  
  if (NumEvt > nentries) NumEvt = nentries;
  for (Int_t jentry=0; jentry<NumEvt; jentry++) {

    // Keep the track info
    // ===================
    Int_t ientry = LoadTree(jentry);
    nb = fChain->GetEntry(jentry);
    nbytes += nb;

    if ((jentry%trkPrescale)<10)
    {
    // Prescaled track count: kentry
    ++kentry;

    // Print 
    // =====
    //if (kentry%200000 == 0) printf("LoopAtt:track %8d (run %d) \n",kentry,lightspeed_run);
    
    // If this is a new event, fill histos and variables for the previous one
    // ======================================================================
    if (current_event==-1) current_event = lightspeed_event;
    if (lightspeed_event != current_event) {
      current_event = lightspeed_event;

      // Loop over bars and fill global variables and histos
      // ===================================================
      double l,ll;
      for (int i=0; i<NumBar; i++) {
	
        if(cont[i] != 1) continue;// Require 1 and only 1 hit in the bar

        adc1 = Adc_e[i];
        adc2 = Adc_w[i];
	
        if(adc1 > 0 && adc2 > 0){ 

	  path[i][nvt[i]] = PathB[i];

	  adchis->Fill(adc1);
	  zpos[i][nvt[i]] = Z[i];
	  l = L/2 - Z[i];
	  ll = L/2 + Z[i];
	  attlene[i]->Fill(adc1); 
	  attlenw[i]->Fill(adc2); 
	  adcreshw[i]->Fill(Z[i],adc2/(path[i][nvt[i]]));
	  adcreshe[i]->Fill(Z[i],adc1/(path[i][nvt[i]]));
	  attp[i][nvt[i]]  = adc1/(path[i][nvt[i]]);
	  attpp[i][nvt[i]] = adc2/(path[i][nvt[i]]); 
	  
	  attlenecp[i]->Fill(adc1/path[i][nvt[i]]); 
	  attlenwcp[i]->Fill(adc2/path[i][nvt[i]]); 
	  
	  if(fabs(zpos[i][nvt[i]])<90.0)attlenec[i]->Fill(adc1/(path[i][nvt[i]]*exp(-(l/result[i][0]))));
	  
	  
	  if(-130<zpos[i][nvt[i]]&& zpos[i][nvt[i]]<-65.0)attleneb1[i]->Fill(adc1/(path[i][nvt[i]]*exp(-(l/result[i][0]))));
	  if(-65<zpos[i][nvt[i]]&& zpos[i][nvt[i]] <0.0)attleneb2[i]->Fill(adc1/(path[i][nvt[i]]*exp(-(l/result[i][0]))));
	  if(0.<zpos[i][nvt[i]] && zpos[i][nvt[i]]<65.0)attleneb3[i]->Fill(adc1/(path[i][nvt[i]]*exp(-(l/result[i][0]))));
	  if(65.<zpos[i][nvt[i]] && zpos[i][nvt[i]]<130.0)attleneb4[i]->Fill(adc1/(path[i][nvt[i]]*exp(-(l/result[i][0]))));
	  
	  
	  adcreshec[i]->Fill(Z[i],adc1/(path[i][nvt[i]]*exp(-(l/result[i][0]))));
	  if(fabs(zpos[i][nvt[i]])<90.0)attlenwc[i]->Fill(adc2/(path[i][nvt[i]]*exp(-(ll/result[i][0]))));
	  if(-130<zpos[i][nvt[i]] && zpos[i][nvt[i]]<-65.0)attlenwb1[i]->Fill(adc2/(path[i][nvt[i]]*exp(-(ll/result[i][0]))));
	  if(-65<zpos[i][nvt[i]] && zpos[i][nvt[i]]<0.0)attlenwb2[i]->Fill(adc2/(path[i][nvt[i]]*exp(-(ll/result[i][0]))));
	  if(0.<zpos[i][nvt[i]] && zpos[i][nvt[i]]<65.0)attlenwb3[i]->Fill(adc2/(path[i][nvt[i]]*exp(-(ll/result[i][0]))));
	  if(65.<zpos[i][nvt[i]] && zpos[i][nvt[i]]<130.0)attlenwb4[i]->Fill(adc2/(path[i][nvt[i]]*exp(-(ll/result[i][0]))));
	  
	  adcreshwc[i]->Fill(Z[i],adc2/((path[i][nvt[i]])*exp(-(ll/result[i][0]))));
	  //to get the sigma (Qe-Qw)
	  if(Flag>1){

	    intresiduale[i]->Fill(adc1/path[i][nvt[i]] - 
		    date[i][1]*exp(-(l/date[i][3])+date[i][5]*l*l));
	    intresidualw[i]->Fill(adc2/path[i][nvt[i]] - 
		    datw[i][1]*exp(-(ll/datw[i][3])+datw[i][5]*ll*ll));
	    residuale[i]->Fill(Z[i], adc1/path[i][nvt[i]] - 
		    date[i][1]*exp(-(l/date[i][3])+date[i][5]*l*l));
	    residualw[i]->Fill(Z[i], adc2/path[i][nvt[i]] - 
		    datw[i][1]*exp(-(ll/datw[i][3])+datw[i][5]*ll*ll));

// 	    if (attp[i][nvt[i]] < 4000.) 
	    attlenec_A[i]->Fill(adc1/(path[i][nvt[i]]*exp(-(l/date[i][3]))*exp(date[i][5]*l*l)));
// 	    if (attpp[i][nvt[i]] < 4000.) 
	    attlenwc_A[i]->Fill(adc2/(path[i][nvt[i]]*exp(-(ll/datw[i][3]))*exp(datw[i][5]*ll*ll)));
	    
	    attlenec_N[i]->Fill(adc1/(path[i][nvt[i]]*exp(-(l/date[i][3]))*exp(date[i][5]*l*l)*date[i][1]));
	    attlenwc_N[i]->Fill(adc2/(path[i][nvt[i]]*exp(-(ll/datw[i][3]))*exp(datw[i][5]*ll*ll)*datw[i][1]));
	    
                                                          	    charge->Fill((exp(-(l/date[i][3]))*exp(date[i][5]*l*l)*date[i][1]));
	    charge_pn->Fill(path[i][nvt[i]]*(exp(-(l/date[i][3]))*exp(date[i][5]*l*l)*date[i][1]));
	    charge_p->Fill(path[i][nvt[i]]*4.*(exp(-(l/date[i][3]))*exp(date[i][5]*l*l)*date[i][1]));
	    charge->Fill((exp(-(ll/datw[i][3]))*exp(datw[i][5]*ll*ll)*datw[i][1]));
	    charge_pn->Fill(path[i][nvt[i]]*(exp(-(ll/datw[i][3]))*exp(datw[i][5]*ll*ll)*datw[i][1]));
	    charge_p->Fill(path[i][nvt[i]]*4.*(exp(-(ll/datw[i][3]))*exp(datw[i][5]*ll*ll)*datw[i][1]));
	    chargem->Fill(adc1);
	    chargem->Fill(adc2);
	    pathl->Fill(path[i][nvt[i]]*4.);
	    
	    double tmp = (adc1/(path[i][nvt[i]]*exp(-(l/date[i][3]))*exp(date[i][5]*l*l)*date[i][1])) - (adc2/(path[i][nvt[i]]*exp(-(ll/datw[i][3]))*exp(datw[i][5]*ll*ll)*datw[i][1]));
	    sigmaQ[i]->Fill(zpos[i][nvt[i]],tmp);
	    sigmaQs[i]->Fill(tmp);
	    sigmaQsr[i]->Fill(tmp/date[i][1]);
	    if (Z[i]>65.) sigmaQz[i][0]->Fill(tmp);
	    else if (Z[i]>0.) sigmaQz[i][1]->Fill(tmp);
	    else if (Z[i]>-65.) sigmaQz[i][2]->Fill(tmp);
	    else sigmaQz[i][3]->Fill(tmp);
	    
	  }
	  nvt[i]++;
	}
	
      } // loop over bars                  

      // Now initialize the local variables for the incoming event
      // =========================================================
      for (int i=0; i<NumBar; i++) {
        cont [i] =    0;
        Adc_e[i] = -999;
        Adc_w[i] = -999;   
        Z    [i] = -999;
	PathB[i] = -999.;
      }
      
    } 
    
    // Select the track and save its information
    // =========================================
    int bar_number = int(lightspeed_bar);
    if (  // (fabs(lightspeed_d0wrt) <  1.5) 
         (lightspeed_track_type  !=   2)         
        && (lightspeed_pt          >  0.4)         
        && (lightspeed_SiPhiHits   >    2)      
        && (lightspeed_cotAxHits   >   10)      
        && (lightspeed_cotStHits   >   10)      
        && (lightspeed_sl8Hits     >    4)
        && (lightspeed_inBar       ==   1)    
       ) {              
      
      ++cont[bar_number];
      
      Adc_e[bar_number] = lightspeed_qe;
      Adc_w[bar_number] = lightspeed_qw;
      Z    [bar_number] = lightspeed_z;	  
      PathB[bar_number] = lightspeed_pathinBar/4.;
      
    } // track selection

    } // trkPrescale

  } // loop over tracks

  printf("LoopAtt: Tracks used %8d (run %d) \n",kentry,lightspeed_run);

  // Declare variables for the final fits
  // ==================================== 
  double par[7], epar[7];
  //double cov[3], datae[7];
  double par1;
  //double epar1, epar3, epar2, par3, par2;
  TF1 *fun;


  if (Flag>1) {

    ofstream parresulte, parresultw;    
    sprintf(filesl,"%s/files/%s_Att-calib_e.txt",directory,file_name);
    parresulte.open(filesl);    
    sprintf(filesl,"%s/files/%s_Att-calib_w.txt",directory,file_name);
    parresultw.open(filesl);

    TF1 *lan = new TF1("lan","landau+[3]",0.,10000.);

    double resolution[216][5];

    for (int i = 0; i<432; i++) {
      //if (!(i<5 || (i>215 && i<221))) continue;
    
      if (i<216)
	lan->SetParameters(800.,date[i][1],date[i][7],1.);
      else if (i>215)
	lan->SetParameters(800.,datw[i][1],datw[i][7],1.);
      
      if (i<216)
	attlenec_A[i]->Fit("lan","","",0.,7000.);
      else if (i>215)
	attlenwc_A[i-216]->Fit("lan","","",0.,7000.);

      double newsigma = lan->GetParameter(2); 
      double enewsigma = lan->GetParameter(3); 

/////////////////////////////////////////////////////////////////
// Correct bad fits
/////////////////////////////////////////////////////////////////

      if(newsigma<40 || newsigma > 249){
         TF1 *lan1 = new TF1("lan1","landau",0.,10000.);
         if (i<216)      attlenec_A[i]    ->Fit("lan1","qz","",0.,6000.);
         else if (i>215) attlenwc_A[i-216]->Fit("lan1","qz","",0.,6000.);

         TF1 *lan2 = new TF1("lan2","landau+[3]",0.,10000.);
         double par0_temp = lan1->GetParameter(0);
         double par1_temp = lan1->GetParameter(1);
         double par2_temp = lan1->GetParameter(2);
         lan2->SetParameters(par0_temp,par1_temp,par2_temp,5.);
    
         if (i<216)      attlenec_A[i]    ->Fit("lan2","qz","",0.,8000.);
         else if (i>215) attlenwc_A[i-216]->Fit("lan2","qz","",0.,8000.);

         newsigma  = lan2->GetParameter(2); 
         enewsigma = lan2->GetParError(2); 
      
      } // if(newsigma<40 || newsigma > 249)

/////////////////////////////////////////////////////////////////

      if (i<216) {
	parresulte << i; 
	parresulte << setw(13) << date[i][1]; 
	parresulte << setw(13) << date[i][2]; 
	parresulte << setw(13) << date[i][3]; 
	parresulte << setw(13) << date[i][4]; 
	parresulte << setw(13) << date[i][5]; 
	parresulte << setw(13) << date[i][6]; 
	parresulte << setw(13) << newsigma << setw(13) << enewsigma; 
// 	parresulte << setw(13) << date[i][9]; 
// 	parresulte << setw(13) << date[i][10]; 
      } else if (i>215) {
	parresultw << i-216;
	parresultw << setw(13) << datw[i-216][1]; 
	parresultw << setw(13) << datw[i-216][2]; 
	parresultw << setw(13) << datw[i-216][3]; 
	parresultw << setw(13) << datw[i-216][4]; 
	parresultw << setw(13) << datw[i-216][5]; 
	parresultw << setw(13) << datw[i-216][6]; 
	parresultw << setw(13) << newsigma << setw(13) << enewsigma; 
// 	parresultw << setw(13) << datw[i-216][9];  
// 	parresultw << setw(13) << datw[i-216][10]; 
      } 

      if (i<216) {

	TF1 *gau = new TF1("gau","gaus",-2.,2.);
	gau->SetParameters(100., 0., 0.07);
	
	for (int k = 0; k<4; k++) {
	  
	  float thrEntry = 0.15*sigmaQz[i][k]->GetMaximum();
	  int bin1 = 1;
	  while (sigmaQz[i][k]->GetBinContent(bin1)<thrEntry)
	    bin1++;
	  int bin2 = 200;
	  while (sigmaQz[i][k]->GetBinContent(bin2)<thrEntry)
	    bin2--;

	  sigmaQz[i][k]->Fit("gau","","",-2.+0.02*(bin1-1),-2.+0.02*bin2);
	  
	  resolution[i][k] = gau->GetParameter(2);
	  
	}

	float thrEntry = 0.15*sigmaQs[i]->GetMaximum();
	int bin1 = 1;
	while (sigmaQs[i]->GetBinContent(bin1)<thrEntry)
	  bin1++;
	int bin2 = 200;
	while (sigmaQs[i]->GetBinContent(bin2)<thrEntry)
	  bin2--;

	sigmaQs[i]->Fit("gau","","",-2.+0.02*(bin1-1),-2.+0.02*bin2);
	  
	resolution[i][4] = gau->GetParameter(2);

	parresulte << setw(13) << resolution[i][4];
	parresulte << endl;
	
      } else {
	
	parresultw << setw(13) << resolution[i-216][4];
	parresultw << endl;
	
      }
      
    }
      
    parresulte.close();
    parresultw.close();

    float meanrmsres = 0.;
    for (int l = 0; l<216; l++) {

      float meanres = 0., rmsres = 0.;
      for (int k = 0; k<4; k++) 
	meanres += resolution[l][k];
      meanres = meanres/4.;
      for (int k = 0; k<4; k++) 
	rmsres += (resolution[l][k]-meanres)*(resolution[l][k]-meanres);
      rmsres = sqrt(rmsres);
      meanrmsres += rmsres;

    }
    cout << meanrmsres/216 << endl;

  }
  
  if (Flag<2) {

    // Open the output text files
    // ==========================
    if (Flag==1)
      sprintf(filesl,"%s/files/%s_Att-calib_t.txt",directory,file_name); 
    else if (Flag==0)
      sprintf(filesl,"%s/files/%s_Att-calib_t_s2.txt",directory,file_name);
    att_calib.open(filesl);
    
    // Loop over bars
    // ==============
    for(int i=0;i<432;i++) {
      //if (!(i<5 || (i>215 && i<221))) continue;
      
      if ((i/100)*100==i)
	cout << "Now fitting channel " << i << endl;
      
      if(i<216){
	evts  = nvt[i];
	////attlenec[i]->Fit("landau","Q");   
	////fun = attlenec[i]->GetFunction("landau");
	//par3  = fun->GetParameter(2);
	////par1  = fun->GetParameter(1); 
	//par2  = fun->GetParameter(2); 
	//epar3  = fun->GetParError(2);
	//epar1  = fun->GetParError(1); 
	//epar2  = fun->GetParError(2); 
      }
      else if(i>215){
	evts  = nvt[i-216];
	////attlenwc[i-216]->Fit("landau","Q");   
	////fun = attlenwc[i-216]->GetFunction("landau");
	//par3  = fun->GetParameter(2);
	////par1  = fun->GetParameter(1); 
	//par2  = fun->GetParameter(2); 
	//epar3  = fun->GetParError(2);
	//epar1  = fun->GetParError(1); 
	//epar2  = fun->GetParError(2); 
      }
      
      
      jchan = i;
      for (int k= 0; k<evts;k++){
	
	if(i<216) {
	  att[k] = attp[i][k];
	  zpp[k] = zpos[i][k];
	} else if(i>215){
	  att[k] = attpp[i-216][k];
	  zpp[k] = zpos[i-216][k];
	}
	
      }
      
      
      gMinuit->mncler();
      gMinuit->SetFCN(fcnatt);
      //gMinuit->mnparm(0,"Q",par1,1.,100.,4000.,ierr);
      gMinuit->mnparm(0,"Q",500.,1.,100.,4000.,ierr);
      //gMinuit->mnparm(1,"AttLen",result[i][0],1.,100.,400.,ierr);
      gMinuit->mnparm(1,"AttLen",180.,1.,100.,400.,ierr);
      gMinuit->mnparm(2,"a",0.000016,0.000001,0.,0.001,ierr);
      gMinuit->mnparm(3,"sigma",250.0,0.01,50.,1000.,ierr);
      gMinuit->mnparm(4,"back",0.99,0.01,0.,1.,ierr);
      arglis[0] = 1000;
      arglis[1] = 1;
      
      gMinuit->mnexcm("MIGRAD",arglis,1,ierr);
      
      
      if ( ierr != 0 ) { cout<<" esto fallo por: "<< ierr<<" bar="<< i<<endl;
      for(int k=0;k<6;k++){par[k] = 0.0;
      epar[k] = 0.0;}}
      
      gMinuit->GetParameter(0,par[0],epar[0]);
      gMinuit->GetParameter(1,par[1],epar[1]);
      gMinuit->GetParameter(2,par[2],epar[2]);
      gMinuit->GetParameter(3,par[3],epar[3]);
//       gMinuit->GetParameter(4,par[4],epar[4]);
//       gMinuit->GetParameter(5,par[5],epar[5]);
      
      
      att_calib << i;
      att_calib<< setw(13) << par[0] 
	       << setw(13) << epar[0]   
	       << setw(13) << par[1]   
	       << setw(13) << epar[1]   
	       << setw(13) << par[2]
	       << setw(13) << epar[2]   
	       << setw(13) << par[3]   
	       << setw(13) << epar[3]    
// 	       << setw(13) << par[4]   
// 	       << setw(13) << epar[4]   
	       << endl;  
      
    } // loop over bars
    
    // close output files  
    att_calib.close();

  }
  
  hfile->Write(); 
  hfile->Close();                                       

} // final_calib::LoopAttl


// Variables for the residual corrections on z and PMT resolutions ---------

bool verbose_timeres = false;

int cont[NumBar], cont1[NumBar], no_cut_cont[NumBar];
int ibar, ibar_trkhi, ibar_trk[50];
int ADC, ichan, indx, flag_t, jj, max_tdc, nhit, track_used;

double datasp [9];
double datatwe[9], datatze[9], datat0e[5];
double datatww[9], datatzw[9], datat0w[5];

double spdlight[NumBar][2];
double spdlight_e[NumBar], t0_off_e[NumBar], tz_e[NumBar][4];
double spdlight_w[NumBar], t0_off_w[NumBar], tz_w[NumBar][4];
double adc_e[NumBar], t_e[NumBar], twalk_e[NumBar];
double adc_w[NumBar], t_w[NumBar], twalk_w[NumBar];

double curv[NumBar], trk_p[NumBar], trk_pt[NumBar], z[NumBar], z0bc[NumBar], pathbar[NumBar], inbar[NumBar];

double v, ve, vw, gammat, gammae, gammaw, beta, diff, p;
double masa, masa1, masae, masaw, p_trkhi, path_trkhi, pt_hi, z_trkhi;
double check_highest_pt;

double ADCe[50], ADCw[50], TDCe[50], TDCw[50], pt[50], z0diff[50];
double curv_trk[50], p_trk[50], path_trk[50], z_trk[50], z0bc_trk[50], pathbar_trk[50], inbar_trk[50];

double timeofflight[maxevt], timeofflight_p[maxevt][3];
double difft[maxevt], zp[maxevt], pt_trk[maxevt], et0_evt[maxevt], qzfit[maxevt];

//==========================================================================
//
// FCN(NPAR,GRAD,FVAL,XVAL,IFLAG) 
//
// npar  ==  NPAR  Number of currently variable parameters.
// gin   ==  GRAD  The (optional) vector of first derivatives.
// f     ==  FVAL  The calculated function value.
// par   ==  XVAL  Vector of (constant and variable) parameters.
// iflag ==  IFLAG Indicates what is to be calculated.
//
//==========================================================================

void fcn(int    &npar,
         double *gin,
         double &f,
         double *par,
         int     iflag)
{
  const double mpid[npid] = {pimass,kmass,pmass};
  const double fpid[npid] = {   0.8,  0.1,  0.1};

  t0g = par[0];
  double sig  =    0.1; // resolution at the ends of the bar (0.1 ns)
  double lin  = 0.0004; // slope of the resolution degradation along the bar
  double logl =    0.0;

  for (int i=0; i<nhit; i++) {
    double arg = 0.0;
    if (i == track_used || pt[i] < 0.0) continue;
    int indx = ibar_trk[i];
    if (cont[indx] != 1) continue; // only one track per bar

    for (int j=0; j<npid; j++) {  // npid = 0 (pion), 1 (kaon), 2 (proton)

      tof = (path_trk[i]/c)*sqrt(1+mpid[j]*mpid[j]/(p_trk[i]*p_trk[i]));
  
      // Difference East = de
      // (t0_off_e = east cable offset)
      
      double de = t0_off_e[indx] - TDCe[i] - tof - t0g;
      if (ADCe[i] > 0) de = de - twalk_e[indx]/sqrt(ADCe[i]);
      de = de - ((L/2)-z_trk[i])/spdlight[indx][0];

      // If z residuals correction is calculated (flag_t=0) we include it in "de"

      if (flag_t == 0) {
        de = de - ( tz_e[indx][0]
                   +tz_e[indx][1]*(z_trk[i]+L/2)
                   +tz_e[indx][2]*(z_trk[i]+L/2)*(z_trk[i]+L/2)
                   +tz_e[indx][3]*(z_trk[i]+L/2)*(z_trk[i]+L/2)*(z_trk[i]+L/2));
      }

      double ee = sig + lin*(L/2-z_trk[i]); // resolution equation along the east side of the bar 
      double pe = exp(-de*de/(2*ee*ee))/(r2pi*ee);

      // Difference West = de
      // (t0_off_w = west cable offset)
             
      double dw = t0_off_w[indx] - TDCw[i] - tof - t0g;
      if (ADCw[i] > 0) dw = dw - twalk_w[indx]/sqrt(ADCw[i]);
      dw = dw - ((L/2)+z_trk[i])/spdlight[indx][0];

      if (flag_t == 0) {
        dw = dw - ( tz_w[indx][0]
                   +tz_w[indx][1]*(z_trk[i]+L/2)
                   +tz_w[indx][2]*(z_trk[i]+L/2)*(z_trk[i]+L/2)
                   +tz_w[indx][3]*(z_trk[i]+L/2)*(z_trk[i]+L/2)*(z_trk[i]+L/2));
      }
     	
      double ew = sig + lin*(L/2+z_trk[i]); // resolution equation along the west side of the bar
      double pw = exp(-dw*dw/(2*ew*ew))/(r2pi*ew);

      double pp =  sqrt(1 + kmass *kmass /(p_trk[i]*p_trk[i])) 
                  -sqrt(1 + pimass*pimass/(p_trk[i]*p_trk[i]));
      pp   = exp(-pp*pp);

      arg += fpid[j]*pe*pw*pp; // weight for the p of the track j

    } // for..npid 

    if (arg > 0) logl += log(arg);

  } // for..nhit

  f = -2.0*logl;

} // void fcn

// void fcnZres(int    &npar,
//              double *gin,
//              double &f,
//              double *par,
//              int     iflag)
// {
//   double a    = par[0];
//   double b    = par[1];
//   double c    = par[2];
//   double d    = par[3];
//   double sig  =    0.1;
//   double lin  = 0.0004;
//   double logl =    0.0;
//   double dy, t, z, polz, tmp, l;

//   for (int i=0; i<evts; i++) {
//     t = difft[i];
//     if (fabs(t) > 0.7) continue;
//     z = zp[i] + L/2;
//     if (ichan < 216) l = L/2 - zp[i];
//     else             l = L/2 + zp[i];
//     tmp   = sig + lin*l;
//     polz  = a + b*z + c*z*z + d*z*z*z; 
//     dy    = (t-polz)/tmp;
//     logl += dy*dy;
//   }

//   f = logl;

// } // void fcnZres

int npar = 16; float cfitvalue;
void fcnZres(int    &npar,
             double *gin,
             double &f,
             double *par,
             int     iflag)
{
//   double a    = par[0];
//   double b    = par[1];
//   double c    = par[2];
//   double d    = par[3];
  double sig  =    0.1;
  double lin  = 0.0004;
  double logl =    0.0;
  double dy, t, z, polz, tmp, l, q;
  for (int i=0; i<evts; i++) {
    t = difft[i];
    q = qzfit[i];
    //if (fabs(t) > 0.7) continue;
    if (fabs(t-cfitvalue) > 0.7) continue;
    //if (q < 300.) continue;
    z = zp[i] + L/2;
    if (ichan < 216) l = L/2 - zp[i];
    else             l = L/2 + zp[i];
    tmp   = sig + lin*l;

//     polz =  par[0] + 
//       (par[1])*exp((par[2])*z) +
//       (par[3])*exp((par[4])*z);

//     polz =  par[0]+par[5]/sqrt(q) + 
//            (par[1]+par[6]/sqrt(q))*exp((par[2]+par[7]/sqrt(q))*z) +
//            (par[3]+par[8]/sqrt(q))*exp((par[4]+par[9]/sqrt(q))*z);

//     polz  = (par[0]  + par[1]*z  + par[2]*z*z  + par[3]*z*z*z )         +
//             (par[4]  + par[5]*z  + par[6]*z*z  + par[7]*z*z*z )/sqrt(q) +
//             (par[8]  + par[9]*z  + par[10]*z*z + par[11]*z*z*z)/q       +
//             (par[12] + par[13]*z + par[14]*z*z + par[15]*z*z*z)/(q*sqrt(q)); 

    polz  = (par[0]  + par[2]*z  + par[3]*z*z) +
            (par[1]  + par[4]*z)/sqrt(q) +
            (par[5])/q;

    dy    = (t-polz)/tmp;
    logl += dy*dy;
  }
  //cout << " " << logl << endl;
  f = logl;

} // void fcnZres

//==========================================================================
//
// FCN(NPAR,GRAD,FVAL,XVAL,IFLAG) 
//
// npar  ==  NPAR  Number of currently variable parameters.
// gin   ==  GRAD  The (optional) vector of first derivatives.
// f     ==  FVAL  The calculated function value.
// par   ==  XVAL  Vector of (constant and variable) parameters.
// iflag ==  IFLAG Indicates what is to be calculated.
//
//==========================================================================

void timeres(int    &npar,
             double *gin,
             double &f,
             double *par,
             int     iflag)
{
  double fpid[npid];
  fpid[0] = par[2];             // par[2] = f_pi
  fpid[1] = (1-par[2])*par[3];  // par[3] = kop = f_k / (f_k + f_p)
  fpid[2] = (1-par[2])*(1-par[3]);

  // fpid[0] = f_pi
  // fpid[1] = f_k
  // fpid[2] = f_p

  double sigma = par[0];
  double lin   = par[1];
  double frac  = par[4];
  double logl  = 0.0;

  for (int i=0; i<evts; i++) {
    double arg = 0;

    if (timeofflight[i] < 0) continue;

    double eet = sigma+lin*zp[i];    
    double ee  = sqrt(eet*eet + et0_evt[i]*et0_evt[i]);
    double eeb = par[5];
    double de  = timeofflight_p[i][0] - timeofflight[i];
    double pe1 = fpid[0]*exp(-de*de/(2.0*ee*ee))/(r2pi*ee);   
    double pb1 = exp(-(de-par[6])*(de-par[6])/(2.0*eeb*eeb))/(r2pi*eeb);

    de = timeofflight_p[i][1] - timeofflight[i];
    double pe2 = fpid[1]*exp(-de*de/(2.0*ee*ee))/(r2pi*ee);   

    de = timeofflight_p[i][2] - timeofflight[i];
    double pe3 = fpid[2]*exp(-de*de/(2.0*ee*ee))/(r2pi*ee);

    arg = frac*(pe1+pe2+pe3)+(1-frac)*pb1;
        if (verbose_timeres) cout << " timeres:arg=" << arg << endl;
    if (arg > 0.0) logl += log(arg);

  } // for..evts

  f = -2.0*logl;

  if (verbose_timeres) cout << " timeres:after  loop, fpid[0]=" << fpid[0]
                            << " fpid[1]=" << fpid[1]
                            << " fpid[2]=" << fpid[2]
                            << " sigma="   << sigma
                            << " lin="     << lin
                            << " frac="    << frac 
                            << " f="       << f << endl;

} // void timeres

double timeres2(double par[6],
                int    i)
{
  double fpid[npid];
  fpid[0] = par[2];
  fpid[1] = (1-par[2])*par[3];
  fpid[2] = (1-par[2])*(1-par[3]);

  double sigma  = par[0];
  double lin    = par[1];
  double frac   = par[4];
  double ee     = sigma+lin*zp[i];   
  double eeb    = ee*par[5];
  double arg    = 0.0;
  double delta1 = timeofflight_p[i][1] - timeofflight_p[i][0];
  double delta2 = timeofflight_p[i][2] - timeofflight_p[i][0];
  double de     = timeofflight[i] - timeofflight_p[i][0];
  double pe1    = fpid[0]*exp(-de*de/(2.0*ee*ee))/(r2pi*ee);
  double pb     = exp(-de*de/(2.0*eeb*eeb))/(r2pi*eeb);

  de = timeofflight_p[i][1] - timeofflight[i] - delta1;
  double pe2 = fpid[1]*exp(-de*de/(2.0*ee*ee))/(r2pi*ee);

  de = timeofflight_p[i][2] - timeofflight[i] - delta2;
  double pe3 = fpid[2]*exp(-de*de/(2.0*ee*ee))/(r2pi*ee);

  arg = frac*(pe1+pe2+pe3)+(1-frac)*pb;

  return arg;

} // double timeres2

//--------------------------------------------------------------------------
// Residual corrections on z and PMT resolutions
//--------------------------------------------------------------------------

void speed_calib::LoopRes(int     NumEvt,
                          int     flag,  
                          int     maxtrk,
                          int     ntrkt0,
                          double  pt_cut,
                          double  cut,
                          char   *file_name,
                          char   *directory)
{
  if (fChain == 0) return;
  Int_t nentries = Int_t(fChain->GetEntries());
  Int_t nbytes = 0, nb = 0;

  if      (flag == 0) cout << "\n Loop pmt resolutions: Number of tracks = " << nentries << " \n" << endl;
  else if (flag == 2) cout << "\n Loop z residuals results: Number of tracks = "     << nentries << " \n" << endl;
  else if (flag == 1) cout << "\n Loop z residuals: Number of tracks = "     << nentries << " \n" << endl;

  // modified _St1493
  if      (flag == 0) sprintf(filesl,"%s/root/%s_pmtres-calib.root",directory,file_name);
  else if (flag == 1) sprintf(filesl,"%s/root/%s_zres-calib.root",  directory,file_name);
  else if (flag == 2) sprintf(filesl,"%s/root/%s_zres-calib_results.root",  directory,file_name);

//   ofstream t0bias; t0bias.open("t0bias_unc.txt");

  TFile *hfile = new TFile(filesl,"recreate");

  gStyle->SetOptFit(111);

  TH1F *mass1       = new TH1F("mass1",    "mass1",    120, 0,2);  
  TH1F *mass        = new TH1F("mass",     "mass",     120, 0,2); 
  TH1F *masse       = new TH1F("masse",    "masse",    120, 0,2); // mass east
  TH1F *massw       = new TH1F("massw",    "massw",    120, 0,2); // mass west
  TH1F *masslo      = new TH1F("masslo",   "masslo",   120, 0,2); 
  TH1F *mass_diff   = new TH1F("mass_diff","mass_diff",100,-1,1);   

  TH1F *d0                     = new TH1F("d0",                    "d0",                    100,-2,  2);
  TH1F *bar_occ                = new TH1F("bar_occ",               "bar_occ",               216, 0,216);
  TH1F *bar_occ_after          = new TH1F("bar_occ_after",         "bar_occ_after",         216, 0,216);
  TH1F *only1hit               = new TH1F("only1hit",              "only1hit",              216, 0,216);
  TH1F *trk_evt_his            = new TH1F("trk_evt_his",           "trk_evt_his",            35, 0, 34);
  TH1F *trk_evt_after_his      = new TH1F("trk_evt_after_his",     "trk_evt_after_his",      35, 0, 34);
  TH1F *trk_bar_his            = new TH1F("trk_bar_his",           "trk_bar_his",             4, 0,  4);
  TH1F *trk_bar_after_his      = new TH1F("trk_bar_after_his",     "trk_bar_after_his",       4, 0,  4);
  TH1F *not_trk_pt_1           = new TH1F("not_trk_pt_1",          "not_trk_pt_1",          100, 0, 10);
  TH1F *sel_trk_pt_1           = new TH1F("sel_trk_pt_1",          "sel_trk_pt_1",          100, 0, 10);
  TH1F *sel_trk_pt_2           = new TH1F("sel_trk_pt_2",          "sel_trk_pt_2",          100, 0, 10);
  TH1F *sel_trk_pt_hi          = new TH1F("sel_trk_pt_hi",         "sel_trk_pt_hi",         100, 0, 10);
  TH1F *check_highest_pt_track = new TH1F("check_highest_pt_track","check_highest_pt_track",100, 0, 10);
  TH1F *nhit_before_index      = new TH1F("nhit_before_index",     "nhit_before_index",      20, 0, 20);
  TH1F *nhit_after_index       = new TH1F("nhit_after_index",      "nhit_after_index",       20, 0, 20);
  TH1F *off_diff               = new TH1F("off_diff",              "off_diff",              100,-9,  9);

  // Track plots (1)
  TH1F *check_length = new TH1F("check_length","check_length",100,   100,  400);
  TH1F *check_z      = new TH1F("check_z",     "check_z",     100,  -150,  150);
  TH1F *check_curv   = new TH1F("check_curv",  "check_curv",  100,-0.008,0.008);
  TH1F *check_phi0   = new TH1F("check_phi0",  "check_phi0",  100,     0,  6.3);
  TH1F *check_length_qual = new TH1F("check_length_qual","check_length_qual",100,   100,  400);
  TH1F *check_z_qual      = new TH1F("check_z_qual",     "check_z_qual",     100,  -150,  150);
  TH1F *check_curv_qual   = new TH1F("check_curv_qual",  "check_curv_qual",  100,-0.008,0.008);
  TH1F *check_phi0_qual   = new TH1F("check_phi0_qual",  "check_phi0_qual",  100,     0,  6.3);

  // Track selection plots (2)
  TH1F *check_d0wrt_before   = new TH1F("check_d0wrt_before",  "check_d0wrt_before",300,-1.5,1.5);
  TH1F *check_d0wrt_after    = new TH1F("check_d0wrt_after",   "check_d0wrt_after", 300,-1.5,1.5);
  TH1F *check_pt_before      = new TH1F("check_pt_before",     "check_pt_before",   100, 0,5);
  TH1F *check_pt_after       = new TH1F("check_pt_after",      "check_pt_after",    100, 0,5);
  TH1F *check_chi2_before    = new TH1F("check_chi2_before",   "check_chi2_before", 100,0,500);
  TH1F *check_chi2_after     = new TH1F("check_chi2_after",    "check_chi2_after",  100,0,500);
  TH1F *check_cotHits_before = new TH1F("check_cotHits_before","check_cotHits_before",97,0,97);
  TH1F *check_cotHits_after  = new TH1F("check_cotHits_after", "check_cotHits_after", 97,0,97);
  TH1F *check_sl8Hits_before = new TH1F("check_sl8Hits_before","check_sl8Hits_before",13,0,13);
  TH1F *check_sl8Hits_after  = new TH1F("check_sl8Hits_after", "check_sl8Hits_after", 13,0,13);

  TH2F *check_phi_d0_before  = new TH2F("check_phi_d0_before", "check_phi_d0_before",100,0,6.3,100,-2,2);
  TH2F *check_phi_d0_after   = new TH2F("check_phi_d0_after",  "check_phi_d0_after", 100,0,6.3,100,-2,2);

  TH2F *test_charge_e   = new TH2F("test_charge_e", "test_charge_e", 5000,0,5000,150,-20,20);
  TH2F *test_charge_w   = new TH2F("test_charge_w", "test_charge_w", 5000,0,5000,150,-20,20);

  // Track selection plots (3)
  TH1F *check_occ_after_d0wrt_cut   = new TH1F("check_occ_after_d0wrt_cut",  "check_occ_after_d0wrt_cut",  216, 0,216);
  TH1F *check_occ_after_pt_cut      = new TH1F("check_occ_after_pt_cut",     "check_occ_after_pt_cut",     216, 0,216);
  TH1F *check_occ_after_chi2_cut    = new TH1F("check_occ_after_chi2_cut",   "check_occ_after_chi2_cut",   216, 0,216);
  TH1F *check_occ_after_cotHits_cut = new TH1F("check_occ_after_cotHits_cut","check_occ_after_cotHits_cut",216, 0,216);
  TH1F *check_occ_after_sl8Hits_cut = new TH1F("check_occ_after_sl8Hits_cut","check_occ_after_sl8Hits_cut",216, 0,216);
  TH1F *check_occ_after_inBar_cut   = new TH1F("check_occ_after_inBar_cut",  "check_occ_after_inBar_cut",  216, 0,216);

  TH1F *check_occ_after_diff_cut    = new TH1F("check_occ_after_diff_cut",   "check_occ_after_diff_cut",   216, 0,216);
  TH1F *check_occ_after_hi_pt_cut   = new TH1F("check_occ_after_hi_pt_cut",  "check_occ_after_hi_pt_cut",  216, 0,216);
  TH1F *check_occ_after_ntrkt0_cut  = new TH1F("check_occ_after_ntrkt0_cut", "check_occ_after_ntrkt0_cut", 216, 0,216);
  TH1F *check_occ_in_t0_histo       = new TH1F("check_occ_in_t0_histo",      "check_occ_in_t0_histo",      216, 0,216);

  TH1F *diff_for_zres_histo = new TH1F("diff_for_zres_histo","diff_for_zres_histo",100,-100,100);

  TH1F *tracks_per_event_conv = new TH1F("tracks_per_event_conv","tracks_per_event_conv",20,0,20);
  TH1F *tracks_per_event_good_t0 = new TH1F("tracks_per_event_good_t0","tracks_per_event_good_t0",20,0,20);

  TH1F *tf    = new TH1F("tf",   "tf",   400,  0,20); // time of flight
  TH1F *tfe   = new TH1F("tfe",  "tfe",  800,-10,30); // time of flight east
  TH1F *tfw   = new TH1F("tfw",  "tfw",  800,-10,30); // time of flight west 
  TH1F *tf_pi = new TH1F("tf_pi","tf_pi",400,  0,10); // tof, fastest pion
  TH1F *tf_k  = new TH1F("tf_k", "tf_k", 400,  0,10); // tof, fastest kaon 
  TH1F *tf_p  = new TH1F("tf_p", "tf_p", 400,  0,10); // tof, fastest proton

  TH1F *t0    = new TH1F("t0",   "t0",  600,-30.0,30    );     
  TH1F *et0   = new TH1F("et0",  "et0", 100,  0.0, 0.3  ); 
  TH1F *sigme = new TH1F("sigme","sigm", 50,  0.0, 0.5  );     
  TH1F *sigma = new TH1F("sigma","siga", 50,  0.0, 0.5  );     
  TH1F *depe  = new TH1F("depe", "dep",  50,  0.0, 0.001);     
  TH1F *sigmw = new TH1F("sigmw","sigm", 50,  0.0, 0.5  );     
  TH1F *depw  = new TH1F("depw", "dep",  50,  0.0, 0.001); 
  
  TH1F *mean_t0_th1f = new TH1F("mean_t0_th1f","mean_t0_th1f",200,-1,1);  

  TH2F *ADCew      = new TH2F("ADCew",     "ADCew",     600, 0,6000,600,0,6000  );
  TH2F *tf_vs_tfpi = new TH2F("tf_vs_tfpi","tf_vs_tfpi",100, 0,  10,100,0,  10  );
  TH2F *mass_vs_p  = new TH2F("mass_vs_p", "mass_vs_p", 200,-2,   2, 75,0,   1.5);
  TH2F *mass2_vs_p = new TH2F("mass2_vs_p","mass2_vs_p",200,-2,   2, 75,0,   1.5);

  TH1F *de_his   [NumBar]; // tof_E - tof_pi
  TH1F *dw_his   [NumBar]; // tof_W - tof_pi
  TH1F *d_his    [NumBar]; // tof   - tof_pi
  TH1F *diffe_lik[NumBar];
  TH1F *diffw_lik[NumBar];
  TH1F *timediff [NumBar];
  //TH1F *timediff1[NumBar];
  //TH1F *timediff2[NumBar];
  //TH1F *timediff3[NumBar];
  //TH1F *timediff4[NumBar];
  //TH1F *residp   [NumBar]; // residuals for each bar

  //TH1F *check_timediff4_after_diff_cut[NumBar];

  TH2F* reside[NumBar];
  TH2F* residw[NumBar];

  TH2F* residqe[NumBar];
  TH2F* residqw[NumBar];

  TH2F* residiqe_z[NumBar][10];
  TH2F* residiqw_z[NumBar][10];

  TH2F* residgepi;
  TH2F* residgeka;
  TH2F* residgepr;

  TH2F* residgwpi;
  TH2F* residgwka;
  TH2F* residgwpr;

  TH2F* residgpi;
  TH2F* residgka;
  TH2F* residgpr;

  TH2F* residqgepi;
  TH2F* residqgeka;
  TH2F* residqgepr;

  TH2F* residqgwpi;
  TH2F* residqgwka;
  TH2F* residqgwpr;

  TH2F* residqgpi;
  TH2F* residqgka;
  TH2F* residqgpr;

  TH2F* residptgepi;
  TH2F* residptgeka;
  TH2F* residptgepr;

  TH2F* residptgwpi;
  TH2F* residptgwka;
  TH2F* residptgwpr;

  TH2F* residptgpi;
  TH2F* residptgka;
  TH2F* residptgpr;

  TH2F* residge;
  TH2F* residgw;

  TH2F* residgqe;
  TH2F* residgqw;

  TH2F* residgiqe;
  TH2F* residgiqw;

  TH2F* residge_q[10];
  TH2F* residgw_q[10];

  TH2F* residge_iq[10];
  TH2F* residgw_iq[10];

  TH2F* residgqe_z[10];
  TH2F* residgqw_z[10];

  TH2F* residgiqe_z[10];
  TH2F* residgiqw_z[10];

  int nevts[NumBar];

  double *zpos[NumBar], *pt_evt[NumBar], *q_east[NumBar], *q_west[NumBar];
  double *tof_e[NumBar], *tof_w[NumBar], *et0g[NumBar];
  double *tofk[NumBar], *tofp[NumBar], *tofpi[NumBar];

  double  mpar[20], empar[20]; for (int j=0; j<20; j++) {mpar[j] = -777; empar[j] = -777;}

  for(int i=0; i<NumBar; i++) { 

    sprintf(titlesl,"de_%d",i);
    de_his   [i] = new TH1F (titlesl,titlesl,200,-5,5);
    sprintf(titlesl,"dw_%d",i);
    dw_his   [i] = new TH1F (titlesl,titlesl,200,-5,5);
    sprintf(titlesl,"d_%d",i);
    d_his    [i] = new TH1F (titlesl,titlesl,200,-5,5);
    sprintf(titlesl,"like_%d",i);
    diffe_lik[i] = new TH1F (titlesl,titlesl,200,-5,5);
    sprintf(titlesl,"likw_%d",i);
    diffw_lik[i] = new TH1F (titlesl,titlesl,200,-5,5);
    sprintf(titlesl,"diff_%d",i);
    timediff [i] = new TH1F (titlesl,titlesl,250,-5,5);
    //sprintf(titlesl,"diff1_%d",i);
    //timediff1[i] = new TH1F (titlesl,titlesl,250,-5,5);
    //sprintf(titlesl,"diff2_%d",i);
    //timediff2[i] = new TH1F (titlesl,titlesl,250,-5,5);
    //sprintf(titlesl,"diff3_%d",i);
    //timediff3[i] = new TH1F (titlesl,titlesl,250,-5,5);
    //sprintf(titlesl,"diff4_%d",i);
    //timediff4[i] = new TH1F (titlesl,titlesl,250,-5,5);
    //sprintf(titlesl,"inv_resid_%d",i);
    //residp   [i] = new TH1F (titlesl,titlesl,250,-5,5);

    //sprintf(titlesl,"check_diff4_after_diff_cut_%d",i);
    //check_timediff4_after_diff_cut[i] = new TH1F (titlesl,titlesl,250,-5,5);

    sprintf(titlesl,"resid_east_%d",i);
    reside[i] = new TH2F (titlesl,titlesl ,10,0,300,80,-2,2);
    sprintf(titlesl,"resid_west_%d",i);
    residw[i] = new TH2F (titlesl,titlesl ,10,0,300,80,-2,2);

    sprintf(titlesl,"residq_east_%d",i);
    residqe[i] = new TH2F (titlesl,titlesl ,20,0,4000,80,-2,2);
    sprintf(titlesl,"residq_west_%d",i);
    residqw[i] = new TH2F (titlesl,titlesl ,20,0,4000,80,-2,2);

    for (int y = 0; y<10; y++) {
    
      sprintf(titlesl,"residiq_east_bar%d_z%d",i,y);
      residiqe_z[i][y] = new TH2F (titlesl,titlesl ,40,0.015,0.1,80,-2,2);
      sprintf(titlesl,"residiq_west_bar%d_z%d",i,y);
      residiqw_z[i][y] = new TH2F (titlesl,titlesl ,40,0.015,0.1,80,-2,2);

    }

  }
  
  sprintf(titlesl,"resid_east_pi");
  residgepi = new TH2F (titlesl,titlesl ,30,-150,150,80,-2,2);
  sprintf(titlesl,"resid_east_ka");
  residgeka = new TH2F (titlesl,titlesl ,30,-150,150,80,-2,2);
  sprintf(titlesl,"resid_east_pr");
  residgepr = new TH2F (titlesl,titlesl ,30,-150,150,80,-2,2);

  sprintf(titlesl,"residq_east_pi");
  residqgepi = new TH2F (titlesl,titlesl ,40,0,4000,80,-2,2);
  sprintf(titlesl,"residq_east_ka");
  residqgeka = new TH2F (titlesl,titlesl ,40,0,4000,80,-2,2);
  sprintf(titlesl,"residq_east_pr");
  residqgepr = new TH2F (titlesl,titlesl ,40,0,4000,80,-2,2);

  sprintf(titlesl,"residpt_east_pi");
  residptgepi = new TH2F (titlesl,titlesl ,40,0,1.,80,-2,2);
  sprintf(titlesl,"residpt_east_ka");
  residptgeka = new TH2F (titlesl,titlesl ,40,0,1.,80,-2,2);
  sprintf(titlesl,"residpt_east_pr");
  residptgepr = new TH2F (titlesl,titlesl ,40,0,1.,80,-2,2);

  sprintf(titlesl,"resid_west_pi");
  residgwpi = new TH2F (titlesl,titlesl ,30,-150,150,80,-2,2);
  sprintf(titlesl,"resid_west_ka");
  residgwka = new TH2F (titlesl,titlesl ,30,-150,150,80,-2,2);
  sprintf(titlesl,"resid_west_pr");
  residgwpr = new TH2F (titlesl,titlesl ,30,-150,150,80,-2,2);

  sprintf(titlesl,"residq_west_pi");
  residqgwpi = new TH2F (titlesl,titlesl ,40,0,4000,80,-2,2);
  sprintf(titlesl,"residq_west_ka");
  residqgwka = new TH2F (titlesl,titlesl ,40,0,4000,80,-2,2);
  sprintf(titlesl,"residq_west_pr");
  residqgwpr = new TH2F (titlesl,titlesl ,40,0,4000,80,-2,2);

  sprintf(titlesl,"residpt_west_pi");
  residptgwpi = new TH2F (titlesl,titlesl ,40,0,1.,80,-2,2);
  sprintf(titlesl,"residpt_west_ka");
  residptgwka = new TH2F (titlesl,titlesl ,40,0,1.,80,-2,2);
  sprintf(titlesl,"residpt_west_pr");
  residptgwpr = new TH2F (titlesl,titlesl ,40,0,1.,80,-2,2);

  sprintf(titlesl,"resid_mean_pi");
  residgpi = new TH2F (titlesl,titlesl ,30,-150,150,80,-2,2);
  sprintf(titlesl,"resid_mean_ka");
  residgka = new TH2F (titlesl,titlesl ,30,-150,150,80,-2,2);
  sprintf(titlesl,"resid_mean_pr");
  residgpr = new TH2F (titlesl,titlesl ,30,-150,150,80,-2,2);

  sprintf(titlesl,"residq_mean_pi");
  residqgpi = new TH2F (titlesl,titlesl ,40,0,4000,80,-2,2);
  sprintf(titlesl,"residq_mean_ka");
  residqgka = new TH2F (titlesl,titlesl ,40,0,4000,80,-2,2);
  sprintf(titlesl,"residq_mean_pr");
  residqgpr = new TH2F (titlesl,titlesl ,40,0,4000,80,-2,2);

  sprintf(titlesl,"residpt_mean_pi");
  residptgpi = new TH2F (titlesl,titlesl ,40,0,1.,80,-2,2);
  sprintf(titlesl,"residpt_mean_ka");
  residptgka = new TH2F (titlesl,titlesl ,40,0,1.,80,-2,2);
  sprintf(titlesl,"residpt_mean_pr");
  residptgpr = new TH2F (titlesl,titlesl ,40,0,1.,80,-2,2);

  sprintf(titlesl,"resid_east");
  residge = new TH2F (titlesl,titlesl ,30,0,300,80,-2,2);
  sprintf(titlesl,"resid_west");
  residgw = new TH2F (titlesl,titlesl ,30,0,300,80,-2,2);
  
  sprintf(titlesl,"residq_east");
  residgqe = new TH2F (titlesl,titlesl ,40,0,4000,80,-2,2);
  sprintf(titlesl,"residq_west");
  residgqw = new TH2F (titlesl,titlesl ,40,0,4000,80,-2,2);
  
  sprintf(titlesl,"residiq_east");
  residgiqe = new TH2F (titlesl,titlesl ,40,0.015,0.1,80,-2,2);
  sprintf(titlesl,"residiq_west");
  residgiqw = new TH2F (titlesl,titlesl ,40,0.015,0.1,80,-2,2);

  for (int y = 0; y<10; y++) {

    sprintf(titlesl,"resid_east_q%d",y);
    residge_q[y] = new TH2F (titlesl,titlesl ,30,0,300,80,-2,2);
    sprintf(titlesl,"resid_west_q%d",y);
    residgw_q[y] = new TH2F (titlesl,titlesl ,30,0,300,80,-2,2);

    sprintf(titlesl,"resid_east_iq%d",y);
    residge_iq[y] = new TH2F (titlesl,titlesl ,30,0,300,80,-2,2);
    sprintf(titlesl,"resid_west_iq%d",y);
    residgw_iq[y] = new TH2F (titlesl,titlesl ,30,0,300,80,-2,2);
    
    sprintf(titlesl,"residq_east_z%d",y);
    residgqe_z[y] = new TH2F (titlesl,titlesl ,40,0,4000,80,-2,2);
    sprintf(titlesl,"residq_west_z%d",y);
    residgqw_z[y] = new TH2F (titlesl,titlesl ,40,0,4000,80,-2,2);
    
    sprintf(titlesl,"residiq_east_z%d",y);
    residgiqe_z[y] = new TH2F (titlesl,titlesl ,40,0.015,0.1,80,-2,2);
    sprintf(titlesl,"residiq_west_z%d",y);
    residgiqw_z[y] = new TH2F (titlesl,titlesl ,40,0.015,0.1,80,-2,2);

  }

  if (flag==0) {//pmt
    FillParameters(directory, file_name, 4);
    FillParametersT0(directory, file_name, 4);
  } else if (flag==1) {//zres
    FillParameters(directory, file_name, 3);
    FillParametersT0(directory, file_name, 3);
  } else if (flag==2) {//zres results
    FillParameters(directory, file_name, 4);
    FillParametersT0(directory, file_name, 3);
//     ReadParTable("/home/sluca/TofCalibrations/533validation/caf/r6_12.txt");
  }

//   ifstream att_result; 
//   ifstream speed_result;
//   ifstream twe_result, t0e_result; 
//   ifstream tww_result, t0w_result;

// if (true) {

//   sprintf(filesl,"%s/files/%s_timew-calib_e.txt",directory,file_name);
//   twe_result.open(filesl,ios::in);   
//   sprintf(filesl,"%s/files/%s_timew-calib_w.txt",directory,file_name);
//   tww_result.open(filesl,ios::in);
//   sprintf(filesl,"%s/files/%s_T0-calib_e.txt",directory,file_name);
//   t0e_result.open(filesl,ios::in);   
//   sprintf(filesl,"%s/files/%s_T0-calib_w.txt",directory,file_name);
//   t0w_result.open(filesl,ios::in);
//   sprintf(filesl,"%s/files/%s_spdtw-calib.txt",directory,file_name);
//   speed_result.open(filesl,ios::in);

// } else {

//   sprintf(filesl,"/cdf/data03/s1/tof/rocio/Calibration/files/St1493_timew-calib_e.txt");
//   twe_result.open(filesl,ios::in);   
//   sprintf(filesl,"/cdf/data03/s1/tof/rocio/Calibration/files/St1493_timew-calib_w.txt");
//   tww_result.open(filesl,ios::in);
//   sprintf(filesl,"/cdf/data03/s1/tof/rocio/Calibration/files/St1493_T0-calib_e.txt");
//   t0e_result.open(filesl,ios::in);   
//   sprintf(filesl,"/cdf/data03/s1/tof/rocio/Calibration/files/St1493_T0-calib_w.txt");
//   t0w_result.open(filesl,ios::in);
//   sprintf(filesl,"/cdf/data03/s1/tof/rocio/Calibration/files/St1493_spdtw-calib_p.txt");
//   speed_result.open(filesl,ios::in);

// } // modified

//   for (int i=0; i<9; i++) {
//     datatwe[i] = 0.0;
//     datatww[i] = 0.0;
//     datatze[i] = 0.0;
//     datatzw[i] = 0.0;
//     datasp [i] = 0.0;
//     if (i < 5) datat0e[i] = 0.0;
//     if (i < 5) datat0w[i] = 0.0;
//   }

//   for (int j=0; j<NumBar; j++) {
//     spdlight  [j][0] = 0.0;
//     spdlight  [j][1] = 0.0;
//     twalk_e   [j]    = 0.0;
//     twalk_w   [j]    = 0.0;
//     spdlight_e[j]    = 0.0;
//     spdlight_w[j]    = 0.0;
//     t0_off_e  [j]    = 0.0;
//     t0_off_w  [j]    = 0.0;
//     for (int k=0; k<4; k++) {tz_e[j][k] = 0.0; tz_w[j][k] = 0.0;}
//   }

  for (int j=0; j<NumBar; j++) {
//     for (int i=0; i<9; i++) twe_result   >> datatwe[i];
//     for (int i=0; i<9; i++) tww_result   >> datatww[i];
//     for (int i=0; i<5; i++) t0e_result   >> datat0e[i];
//     for (int i=0; i<5; i++) t0w_result   >> datat0w[i];
//     for (int i=0; i<9; i++) speed_result >> datasp [i];

//     twalk_e   [j]    = datatwe[1];
//     twalk_w   [j]    = datatww[1];    
//     //t0_off_e  [j]    = datat0e[1]; // offset if s east and west are different
//     //t0_off_w  [j]    = datat0w[1]; // offset if s east and west are different
//     t0_off_e  [j]    = datat0e[3];   // offset if s east and west are the same
//     t0_off_w  [j]    = datat0w[3];   // offset if s east and west are the same

//     // modified
//     if (true) {spdlight[j][0] = datasp [1];} else {spdlight[j][0] = 2.0/datasp[1];}

//     spdlight  [j][1] = datasp [3];
//     //spdlight_e[j]    = datat0e[2]; // s if s east and west are different
//     //spdlight_w[j]    = datat0w[2]; // s if s east and west are different
//     spdlight_e[j]    = datat0e[4];   // s if s east and west are the same
//     spdlight_w[j]    = datat0w[4];   // s if s east and west are the same

    nevts [j] = 0;
    zpos  [j] = new double[maxevt];
    q_east[j] = new double[maxevt];
    q_west[j] = new double[maxevt];
    tof_e [j] = new double[maxevt];
    tof_w [j] = new double[maxevt];
    tofk  [j] = new double[maxevt];
    tofp  [j] = new double[maxevt];
    tofpi [j] = new double[maxevt];
    pt_evt[j] = new double[maxevt];
    et0g  [j] = new double[maxevt];

    for (int k=0; k<maxevt; k++) {
      zpos  [j][k] = 0.0;
      q_east[j][k] = 0.0;
      q_west[j][k] = 0.0;
      tof_e [j][k] = 0.0;
      tof_w [j][k] = 0.0;
      tofk  [j][k] = 0.0;
      tofp  [j][k] = 0.0;
      tofpi [j][k] = 0.0;
      pt_evt[j][k] = 0.0;
      et0g  [j][k] = 0.0;
    }
  } // for..NumBar

//   twe_result.close();
//   tww_result.close();
//   t0e_result.close();
//   t0w_result.close();
//   speed_result.close();

  // New 031203
  //int currentEventForT0Count = -999;
  //int countingEvents         =    0;
//ofstream event_t0_t0Error_File;
//sprintf(filesl,"%s/files/%s_event-t0.txt",directory,file_name); 
//event_t0_t0Error_File.open(filesl,ios::in);

//   if(flag == 0) {

//     ifstream tze_result, tzw_result;

// if (true) {

//     sprintf(filesl,"%s/files/%s_Zres-calib_e.txt",directory,file_name);
//     tze_result.open(filesl,ios::in);
//     sprintf(filesl,"%s/files/%s_Zres-calib_w.txt",directory,file_name);
//     tzw_result.open(filesl,ios::in);

// } else {

//     sprintf(filesl,"/cdf/data03/s1/tof/rocio/Calibration/files/St1493_Zres-calib_e.txt");
//     tze_result.open(filesl,ios::in);
//     sprintf(filesl,"/cdf/data03/s1/tof/rocio/Calibration/files/St1493_Zres-calib_w.txt");
//     tzw_result.open(filesl,ios::in);

// } // modified

//     for (int j=0; j<NumBar; j++) {   
//       for (int i=0; i<9; i++) tze_result >> datatze[i];
//       for (int i=0; i<9; i++) tzw_result >> datatzw[i];

//       tz_e[j][0] = datatze[1];
//       tz_e[j][1] = datatze[3];
//       tz_e[j][2] = datatze[5];
//       tz_e[j][3] = datatze[7];
//       tz_w[j][0] = datatzw[1];
//       tz_w[j][1] = datatzw[3];
//       tz_w[j][2] = datatzw[5];
//       tz_w[j][3] = datatzw[7];
 
//     } // for..NumBar

//     tze_result.close();
//     tzw_result.close();

//   } // if..flag == 0

  int    ierr     =    0;
  int    ierre    =    0;
  //double t0i      = -999;
  //double et0g_tmp = -999;

  double arglis[10]; for (int i=0; i<10; i++) arglis[i] = -999;

  TMinuit *gMinuit = new TMinuit(npar);
  gMinuit->SetFCN(fcn);
  arglis[0] = 1;
  gMinuit->mnexcm("SET ERR",          arglis,1,ierre);
  arglis[0] = -1;
  gMinuit->mnexcm("SET PRINT",        arglis,1,ierre);
  gMinuit->mnexcm("SET SETNOWARNINGS",arglis,1,ierre);

  TMinuit *gMinuit_r = new TMinuit(8);
  gMinuit_r->SetFCN(timeres);
  arglis[0] = 1;
  gMinuit_r->mnexcm("SET ERR",          arglis,1,ierr);
  arglis[0] = -1;
  gMinuit_r->mnexcm("SET PRINT",        arglis,1,ierr);
  gMinuit_r->mnexcm("SET SETNOWARNINGS",arglis,1,ierr);

// Begin: initialize some variables every event ----------------------------
      
  check_highest_pt = -999;
  pt_hi            =  2.0; // Set the high pt cut
  ibar_trkhi       = -999;
  t0g	           = -999;
  nhit             =    0;

  for (int i=0; i<50; i++) {
    pt      [i] = -999;
    ibar_trk[i] = -999;
  }

  for (int i=0; i<NumBar; i++) {
    no_cut_cont[i] =	0;
    cont       [i] =	0;
    cont1      [i] =	0;
    adc_e      [i] = -999;
    adc_w      [i] = -999;
    curv       [i] = -999;
    length     [i] = -999;
    t_e        [i] = -999;
    t_w        [i] = -999;
    trk_p      [i] = -999;
    trk_pt     [i] = -999;
    z	       [i] = -999;
  }

  jj = 0; // Low pt track counter

// End: initialize some variables every event ------------------------------

  int current_event  = -999;
  //int previous_event = -999;
  int event_counter  =    0;
  int num_events     =    1;

  int ntrk_event                 = 0;
  int ntrk_event_after_quality   = 0;
  //int events_with_nhit_ge_ntrkt0 = 0;
  //int events_with_nhit_lt_ntrkt0 = 0;

  //int total_track = 0;     // Initial number of tracks used to calculate t0

  // Fill parameters for track corrections
  
  Int_t kentry = 0;
  
  // Loop over tracks
  if (NumEvt > nentries) NumEvt = nentries;
  for (Int_t jentry=0; jentry<NumEvt; jentry++) {
    Int_t ientry = LoadTree(jentry);
    nb = fChain->GetEntry(jentry);
    nbytes += nb;

    if ((jentry%trkPrescale)<10)
    {
    // Prescaled track count: kentry
    ++kentry;

    int bar_number = int(lightspeed_bar);

    //if (jentry>5200000/10) continue;

    /*
    if (kentry%200000 == 0) {
      if (flag == 0) printf("Loop pmt resolutions: track %8d (run %d) \n",kentry,lightspeed_run);
      if (flag == 1) printf("Loop z residuals: track %8d (run %d) \n",    kentry,lightspeed_run);
      if (flag == 2) printf("Loop z residuals results:track %8d (run %d) \n",    kentry,lightspeed_run);
    }
    */

    d0->Fill(lightspeed_d0wrt);
    bar_occ->Fill(bar_number);
    ++no_cut_cont[bar_number];

    if (lightspeed_event != current_event) {
      //previous_event = current_event;
      current_event = lightspeed_event;
      ++event_counter;

      if (event_counter > 1) {
        trk_evt_his      ->Fill(ntrk_event              );
        trk_evt_after_his->Fill(ntrk_event_after_quality);

        for (int i=0; i<NumBar; i++) {
          trk_bar_his      ->Fill(no_cut_cont[i]);
          trk_bar_after_his->Fill(       cont[i]);
        }
      }

      ntrk_event               = 0;
      ntrk_event_after_quality = 0;
    }

    ntrk_event++;

    if (event_counter > num_events) {

      num_events = event_counter;

      // Loop over bars
      for (int i=0; i<NumBar; i++) {
        if (cont[i] != 1 || cont1[i] != 0) continue;

        only1hit->Fill(i);

//         if (adc_e[i] > 0) tofe = t_e[i] + twalk_e[i]/sqrt(adc_e[i]); else tofe = t_e[i];
//         if (adc_w[i] > 0) tofw = t_w[i] + twalk_w[i]/sqrt(adc_w[i]); else tofw = t_w[i];
//         if ((tofe > 0) && (tofw > 0)) diff = tofe - tofw; else diff = -999;
//         residp[i]->Fill(diff-(z[i]*2.0/spdlight[i][0]+spdlight[i][1]));

//         tofe = -999;
//         tofw = -999;

//         tofe = t_e[i] + (L/2 - z[i])/spdlight[i][0];
//         tofw = t_w[i] + (L/2 + z[i])/spdlight[i][0];
//         timediff1[i]->Fill(tofe-tofw);

//         if (adc_e[i] > 0 && adc_w[i] > 0) {
//            tofe += twalk_e[i]/sqrt(adc_e[i]);
//            tofw += twalk_w[i]/sqrt(adc_w[i]);
//         }
//         timediff2[i]->Fill(tofe-tofw);
 
//         double diff_for_zres = t_e[i] - t_w[i];
//         if (adc_e[i] > 0 && adc_w[i] > 0) {
//            diff_for_zres += twalk_e[i]/sqrt(adc_e[i]);
//            diff_for_zres -= twalk_w[i]/sqrt(adc_w[i]);
//         }
//         diff_for_zres = (diff_for_zres - spdlight[i][1])*spdlight[i][0]/2.0 - z[i];

//         diff_for_zres_histo->Fill(diff_for_zres);

//         tofe -= t0_off_e[i];
//         tofw -= t0_off_w[i];
//         timediff3[i]->Fill(tofe-tofw);       
//         off_diff    ->Fill(t0_off_e[i]-t0_off_w[i]);

//         if (flag == 0)  tofe +=  tz_e[i][0]
//                        +tz_e[i][1]*(z[i]+L/2)
//                        +tz_e[i][2]*(z[i]+L/2)*(z[i]+L/2)
//                        +tz_e[i][3]*(z[i]+L/2)*(z[i]+L/2)*(z[i]+L/2);

//         if (flag == 0)  tofw +=  tz_w[i][0]
//                        +tz_w[i][1]*(z[i]+L/2)
//                        +tz_w[i][2]*(z[i]+L/2)*(z[i]+L/2)
//                        +tz_w[i][3]*(z[i]+L/2)*(z[i]+L/2)*(z[i]+L/2);

//         diff = tofe - tofw;
//         timediff4[i]->Fill(diff);
	
//         if ((flag == 0) && (fabs(diff) > cut)) continue;          // pmt resolutions (cut = 0.5)
//         if ((flag >  0) && (fabs(diff_for_zres) > cut)) continue; // z residuals     (cut =  10)

//         check_timediff4_after_diff_cut[i]->Fill(diff);
//         check_occ_after_diff_cut->Fill(i);

	if ((trk_pt[i] > pt_cut) && (jj < maxtrk)) { 

          curv_trk[jj] = curv  [i];
	  pt      [jj] = trk_pt[i];
          p_trk   [jj] = trk_p [i];
          z_trk   [jj] = z     [i];      
          path_trk[jj] = length[i]; 
          ibar_trk[jj] = i;  
          ADCe    [jj] = adc_e [i];
          TDCe    [jj] = t_e   [i];
          ADCw    [jj] = adc_w [i];
          TDCw    [jj] = t_w   [i];
          z0bc_trk[jj] = z0bc  [i];
          pathbar_trk[jj] = pathbar[i];

          if (pt[jj] > check_highest_pt) check_highest_pt = pt[jj];

          sel_trk_pt_2->Fill(trk_pt[i]);
          jj++;
	}

        // Store any high pt track for initial t0 calculation
       	if (trk_pt[i] > pt_hi) {
          pt_hi      = trk_pt[i];
          p_trkhi    = trk_p [i];
          path_trkhi = length[i];
          ibar_trkhi = i;
	}
      } // loop over bars

      check_highest_pt_track->Fill(check_highest_pt);
      if (pt_hi > 2.0) sel_trk_pt_hi->Fill(pt_hi);

      nhit = jj;
      nhit_before_index->Fill(nhit);

      // Initial t0 calculation
      //indx = ibar_trkhi; this was for old t0 calculation
      indx = nhit;
      
      if (indx < 0) {
        nhit_after_index->Fill(0);
        tracks_per_event_conv->Fill(0);
        tracks_per_event_good_t0->Fill(0);
      }	
      else {

        nhit_after_index->Fill(nhit);

	// New T0 Calculation
	double CorrectedT0[20], ErrorCorrT0[20]; 
	for (int nv = 0; nv<zv_nzv; nv++) {
	  CorrectedT0[nv] = -999.;
	  ErrorCorrT0[nv] = -999.;
	}
	if (flag>=1) {
	  _skip_p = -999;
	  for (int nv = 0; nv<zv_nzv; nv++) 
	    CorrectedT0[nv] = speed_calib::CalculateT0(nv, &ErrorCorrT0[nv]); 
	}
	
//         if (nhit >= ntrkt0) ++events_with_nhit_ge_ntrkt0;
//         else ++events_with_nhit_lt_ntrkt0;

// //         tof_pi = (p_trkhi/pimass)/sqrt(1+(p_trkhi/pimass)*(p_trkhi/pimass)); 
// //         tof_pi = path_trkhi/(tof_pi*c);
// 	tof_pi = getTofSL(path_trkhi, p_trkhi, pimass);

//         t0gei = t0_off_e[indx] - t_e[indx];
//         t0gei = t0gei - (L/2-z[indx])/spdlight[indx][0];
//         if (adc_e[indx] > 0) t0gei = t0gei - twalk_e[indx]/sqrt(adc_e[indx]);

//         if (flag == 0) t0gei = t0gei - ( tz_e[indx][0]
//                                         +tz_e[indx][1]*(z[indx]+L/2)
//                                         +tz_e[indx][2]*(z[indx]+L/2)*(z[indx]+L/2)
//                                         +tz_e[indx][3]*(z[indx]+L/2)*(z[indx]+L/2)*(z[indx]+L/2));   

//         t0gei = t0gei - tof_pi;

//         t0gwi = t0_off_w[indx] - t_w[indx];
//         t0gwi = t0gwi - (L/2+z[indx])/spdlight[indx][0];
//         if (adc_w[indx] > 0) t0gwi = t0gwi - twalk_w[indx]/sqrt(adc_w[indx]);

//         if (flag == 0) t0gwi = t0gwi - ( tz_w[indx][0]
//                                         +tz_w[indx][1]*(z[indx]+L/2)
//                                         +tz_w[indx][2]*(z[indx]+L/2)*(z[indx]+L/2)
//                                         +tz_w[indx][3]*(z[indx]+L/2)*(z[indx]+L/2)*(z[indx]+L/2));   

//         t0gwi = t0gwi - tof_pi;

//         indx = -999;
//         track_used = -999;

// // ====== New for t0 cleaning =========	012804

// 	int total_track_t0 = 0;

// 	int   track_id[max_t0_trks];
//         float track_t0[max_t0_trks];
// 	float track_e0[max_t0_trks];
	
// 	for (int jjj=0; jjj<max_t0_trks; jjj++) {
// 	  track_id[jjj] = -999;
// 	  track_t0[jjj] = -999;
// 	  track_e0[jjj] = -999;
//         }
	
// // ====== New for t0 cleaning =========

//         // Calculate TOF for all low pt tracks
//         for (int i=0; i<nhit; i++ ) {
//           indx = ibar_trk[i];    

//           check_occ_after_hi_pt_cut->Fill(indx);

//           if (pt[i]    == 0)    break; // Skip event if no more low pt tracks
//           if (nhit < ntrkt0) continue; // Number of tracks for t0

//           if (cont[indx] != 1 || cont1[indx] != 0) continue;

//           check_occ_after_ntrkt0_cut->Fill(indx);

//           // Minuit fit starts here
//           flag_t = flag;

//           track_used = i;
//           t0i = (t0gwi + t0gei)/2;

//           gMinuit->mncler();
//           gMinuit->mnparm(0,"t0g",t0i,0.05,0,0,ierr); 

// //==========================================================================
// //
// // mnparm(NUM,CHNAM,STVAL,STEP,BND1,BND2,IERFLG*)
// //
// // NUM    Parameter number as referenced by user in FCN.
// // CHNAM  Character string of up to 10 characters containing
// //        the name which the user assigned to the given parameter
// // STVAL  Starting value.
// // STEP   Starting step size or approximate parameter error.
// // BND1   Lower bound (limit) on parameter value, if any.
// // BND2   Upper bound (limit) on parameter value, if any.
// // IERFLG Error return code: 0 if no error, >0 if request failed.
// //
// //==========================================================================

//           arglis[0] = 1000; // maxcalls
//           arglis[1] =  1.0; // tolerance   
//           gMinuit->mnexcm("MIGRAD",arglis,2,ierr);
//           gMinuit->GetParameter(0,t0g,et0g_tmp);

//           if (ierr != 0) {
//             cout << " t0 fit failed because ierr = " << ierr 
//                  << " (event " << lightspeed_event << "), et0g = " << et0g_tmp << endl;
//             continue;
//           } 
     
//           // New 012604
//           //if (abs(t0g) > 5.0) {
// 	  //  cout << "Rejected track because t0=" << t0g << ", > 5.0 ns" << endl; 
// 	  //  continue;
// 	  //}

// 	  track_t0[total_track_t0] = t0g;
// 	  track_e0[total_track_t0] = et0g_tmp;
//           track_id[total_track_t0] = i;
//           ++total_track_t0;

//         } // for..nhit. Loop closed here to reject tracks with bad t0
//         tracks_per_event_conv->Fill(total_track_t0);    // number of events with n tracks which converge

// // ====== New for t0 cleaning =========	
       
// 	int total_track_good_t0 = 0; // Total number of tracks with good t0

// 	int   track_good_id[max_t0_trks];
//         float track_good_t0[max_t0_trks];
// 	float track_good_e0[max_t0_trks];

// 	for (int jjj=0; jjj<max_t0_trks; jjj++) {
// 	  track_good_id[jjj] = -999;
// 	  track_good_t0[jjj] = -999;
// 	  track_good_e0[jjj] = -999;
//         }	 
	 
//         for (int jjj=0; jjj<total_track_t0; jjj++) {
// 	  float inv2_e0 = 0.0;
//           float mean_t0 = 0.0;
//           for (int kkk=0; kkk<total_track_t0; kkk++) {
//             if (kkk == jjj) continue;   // We don't consider the track used to calculate the mean
// 	    if (track_e0[kkk] > 0.0) {
// 	      mean_t0 += (track_t0[kkk] / (track_e0[kkk]*track_e0[kkk]));
// 	      inv2_e0 += (1.0 / (track_e0[kkk]*track_e0[kkk]));
// 	    }
// 	  }

//           if (total_track_t0 == 1) mean_t0 = track_t0[jjj];	  
// 	  else mean_t0 = mean_t0 / inv2_e0;
 	  
// 	  mean_t0_th1f->Fill(mean_t0 - track_t0[jjj]);
// 	  if (fabs(mean_t0 - track_t0[jjj]) < 10) {
// 	    track_good_id[total_track_good_t0] = track_id[jjj];
// 	    track_good_t0[total_track_good_t0] = track_t0[jjj];
// 	    track_good_e0[total_track_good_t0] = track_e0[jjj];
// 	    ++total_track_good_t0;
//           }
//         } // for..jjj
//         tracks_per_event_good_t0->Fill(total_track_good_t0);    // number of events with n tracks which pass the t0 cut
	
// // ====== New for t0 cleaning =========	

        for (int kkk=0; kkk<nhit; kkk++ ) {
	
// 	  int i = -999;
// 	  for (int jjj=0; jjj<total_track_good_t0; jjj++) {
// 	    if (kkk == track_good_id[jjj]) i = kkk;
// 	  }
	
// 	  if (i < 0) continue;
	
	  double CorrectedTrackT0 = -999.;
	  double ErrorCorrTrackT0 = -999.;
	  for (int nv = 0; nv<zv_nzv; nv++) 
	    if (fabs(z0bc_trk[kkk]-zv_zv[nv])<2.0) {
	      if (flag>=1) {
		CorrectedTrackT0 = CorrectedT0[nv];
		ErrorCorrTrackT0 = ErrorCorrT0[nv];
	      } else if (flag==0) {
		_skip_p = p_trk[kkk];
		CorrectedTrackT0 = speed_calib::CalculateT0(nv, &ErrorCorrTrackT0); 
	      } 
	    }
	  if (CorrectedTrackT0==-999. || 
	    !(CorrectedTrackT0<5. && CorrectedTrackT0>-5) ||
	    ErrorCorrTrackT0>0.1) continue;
	  int i = kkk;

// 	  float t0g_clean = track_good_t0[i];
// 	  float et0g_tmp_clean = track_good_e0[i];
	  float t0g_clean = CorrectedTrackT0;
	  float et0g_tmp_clean = ErrorCorrTrackT0;
	
          indx = ibar_trk[i];

          check_occ_in_t0_histo->Fill(indx);

//           if (previous_event != currentEventForT0Count) {
//             currentEventForT0Count = previous_event;
//             ++countingEvents;
//           }

          t0 ->Fill(t0g_clean);
          et0->Fill(et0g_tmp_clean);

          // Predicted TOF for different particles
//           tof_pi = (p_trk[i]/pimass)/sqrt(1+(p_trk[i]/pimass)*(p_trk[i]/pimass)); 
//           tof_pi = path_trk[i]/(tof_pi*c);
//           tof_k  = (p_trk[i]/kmass)/sqrt(1+(p_trk[i]/kmass)*(p_trk[i]/kmass)); 
//           tof_k  = path_trk[i]/(tof_k*c);
//           tof_p  = (p_trk[i]/pmass)/sqrt(1+(p_trk[i]/pmass)*(p_trk[i]/pmass)); 
//           tof_p  = path_trk[i]/(tof_p*c);
	  tof_pi = getTofSL(path_trk[i], p_trk[i], pimass);
	  tof_k  = getTofSL(path_trk[i], p_trk[i], kmass );
	  tof_p  = getTofSL(path_trk[i], p_trk[i], pmass );

          // TOF calculated using east pmt's
	  tofe = TimeCorr(TDCe[i], z_trk[i], ADCe[i], indx);
// 	  tofe = CorrectedTime(TDCe[i], z_trk[i], ADCe[i], indx);
//           tofe  = t0_off_e[indx] - TDCe[i];
//           tofe -= (L/2-z_trk[i])/spdlight[indx][0];
//           if (ADCe[i] > 0) tofe -= twalk_e[indx]/sqrt(ADCe[i]);

//           if (flag == 0) tofe -= ( tz_e[indx][0]
//                                   +tz_e[indx][1]*(z_trk[i]+L/2)
//                                   +tz_e[indx][2]*(z_trk[i]+L/2)*(z_trk[i]+L/2)
//                                   +tz_e[indx][3]*(z_trk[i]+L/2)*(z_trk[i]+L/2)*(z_trk[i]+L/2));

          tofe -= t0g_clean; 

          // TOF calculated using west pmt's
	  tofw = TimeCorr(TDCw[i], z_trk[i], ADCw[i], indx+216);
// 	  tofw = CorrectedTime(TDCw[i], z_trk[i], ADCw[i], indx+216);
//           tofw  = t0_off_w[indx] - TDCw[i];
//           tofw -= (L/2+z_trk[i])/spdlight[indx][0];
//           if (ADCw[i] > 0) tofw -= twalk_w[indx]/sqrt(ADCw[i]);

//           if (flag == 0) tofw -= ( tz_w[indx][0]
//                                   +tz_w[indx][1]*(z_trk[i]+L/2)
//                                   +tz_w[indx][2]*(z_trk[i]+L/2)*(z_trk[i]+L/2)
//                                   +tz_w[indx][3]*(z_trk[i]+L/2)*(z_trk[i]+L/2)*(z_trk[i]+L/2));

          tofw -= t0g_clean; 

          //tof = (tofe + tofw)/2.0;
	  double errore = parr[indx][0] + parr[indx][1]*z_trk[i];
	  double errorw = parr[indx+216][0] + parr[indx+216][1]*z_trk[i];
	  double terror = 1./(errore*errore) + 1./(errorw*errorw);
	  tof = tofe/(errore*errore) + tofw/(errorw*errorw);
	  tof = tof/terror;

          tf   ->Fill(tof   );
          tfe  ->Fill(tofe  );   
          tfw  ->Fill(tofw  );
          tf_pi->Fill(tof_pi);   
          tf_k ->Fill(tof_k );
          tf_p ->Fill(tof_p );   

          timediff[indx]->Fill(tofe - tofw);

          masae = 0;
          masaw = 0;
          masa  = 0;
          masa1 = 0;

          if (tofe > -10000 && tofw > -10000 && fabs(tofe - tofw) < 1000000) {
            v      = path_trk[i]/tof;  // speed using both pmts
            ve     = path_trk[i]/tofe; // speed using east pmts
            vw     = path_trk[i]/tofw; // speed using west pmts
            gammae = 0.0;              // gamma east pmts
            gammaw = 0.0;              // gamma west pmts  
            gammat = 0.0;              // gamma both pmts

            if (v  < c) gammat = 1.0/sqrt(1-((v *v )/(c*c)));
            if (ve < c) gammae = 1.0/sqrt(1-((ve*ve)/(c*c)));
            if (vw < c) gammaw = 1.0/sqrt(1-((vw*vw)/(c*c)));

           // Mass calculation
            if (gammae > 0) masae = p_trk[i]/(gammae*(ve/c)); // mass east 
            if (gammaw > 0) masaw = p_trk[i]/(gammaw*(vw/c)); // mass west
            if (gammat > 0) masa1 = p_trk[i]/(gammat*(v /c)); // mass both 

            masa = (masaw+masae)/2.0;
            if (masae > 0)              masse    ->Fill(masae        );
            if (masaw > 0)              massw    ->Fill(masaw        );
            if (masae > 0 && masaw > 0) mass_diff->Fill(masae - masaw);
            if (masa1 > 0)              mass     ->Fill(masa1        );
            if (masa  > 0)              mass1    ->Fill(masa         );
            if (masa1 > 0) {
              if (pt[i] < 1.2) masslo->Fill(masa1);
              if (curv_trk[i] > 0) {
                mass_vs_p ->Fill(p_trk[i],masa1      );
                mass2_vs_p->Fill(p_trk[i],masa1*masa1);
              } else {
                mass_vs_p ->Fill(-p_trk[i],masa1);
                mass2_vs_p->Fill(-p_trk[i],masa1*masa1);
              }
            } // if..masa1 > 0
	    
            // Resolution and stuff
            if (pt[i] > pt_cut && indx<216) {

// 	      if (fabs(tofe - tof_pi)<0.7) 
// 		t0bias << jentry << " " << tofe+t0g_clean << " " << t0g_clean << " " << tof_pi << endl; 

	      if (fabs(tofe - tof_pi)<0.7 && pt[i]<0.8) residgepi->Fill((z_trk[i]),(tofe-tof_pi));
	      if (fabs(tofe - tof_k )<0.7 && pt[i]<0.8) residgeka->Fill((z_trk[i]),(tofe-tof_k ));
	      if (fabs(tofe - tof_p )<0.7 && pt[i]<0.8) residgepr->Fill((z_trk[i]),(tofe-tof_p ));

	      if (fabs(tofe - tof_pi)<0.7 && pt[i]<0.8) residqgepi->Fill(ADCe[i],(tofe-tof_pi));
	      if (fabs(tofe - tof_k )<0.7 && pt[i]<0.8) residqgeka->Fill(ADCe[i],(tofe-tof_k ));
	      if (fabs(tofe - tof_p )<0.7 && pt[i]<0.8) residqgepr->Fill(ADCe[i],(tofe-tof_p ));

	      if (fabs(tofe - tof_pi)<0.7) residptgepi->Fill(pt[i],(tofe-tof_pi));
	      if (fabs(tofe - tof_k )<0.7) residptgeka->Fill(pt[i],(tofe-tof_k ));
	      if (fabs(tofe - tof_p )<0.7) residptgepr->Fill(pt[i],(tofe-tof_p ));

	      if (fabs(tofw - tof_pi)<0.7 && pt[i]<0.8) residgwpi->Fill((z_trk[i]),(tofw-tof_pi));
	      if (fabs(tofw - tof_k )<0.7 && pt[i]<0.8) residgwka->Fill((z_trk[i]),(tofw-tof_k ));
	      if (fabs(tofw - tof_p )<0.7 && pt[i]<0.8) residgwpr->Fill((z_trk[i]),(tofw-tof_p ));

	      if (fabs(tofw - tof_pi)<0.7 && pt[i]<0.8) residqgwpi->Fill(ADCw[i],(tofw-tof_pi));
	      if (fabs(tofw - tof_k )<0.7 && pt[i]<0.8) residqgwka->Fill(ADCw[i],(tofw-tof_k ));
	      if (fabs(tofw - tof_p )<0.7 && pt[i]<0.8) residqgwpr->Fill(ADCw[i],(tofw-tof_p ));
	      
	      if (fabs(tofw - tof_pi)<0.7) residptgwpi->Fill(pt[i],(tofw-tof_pi));
	      if (fabs(tofw - tof_k )<0.7) residptgwka->Fill(pt[i],(tofw-tof_k ));
	      if (fabs(tofw - tof_p )<0.7) residptgwpr->Fill(pt[i],(tofw-tof_p ));

	      if (fabs(tof - tof_pi)<0.7 && pt[i]<0.8) residgpi->Fill((z_trk[i]),(tof-tof_pi));
	      if (fabs(tof - tof_k )<0.7 && pt[i]<0.8) residgka->Fill((z_trk[i]),(tof-tof_k ));
	      if (fabs(tof - tof_p )<0.7 && pt[i]<0.8) residgpr->Fill((z_trk[i]),(tof-tof_p ));

	      if (fabs(tof - tof_pi)<0.7 && pt[i]<0.8) residqgpi->Fill((ADCe[i]+ADCw[i])/2.,(tof-tof_pi));
	      if (fabs(tof - tof_k )<0.7 && pt[i]<0.8) residqgka->Fill((ADCe[i]+ADCw[i])/2.,(tof-tof_k ));
	      if (fabs(tof - tof_p )<0.7 && pt[i]<0.8) residqgpr->Fill((ADCe[i]+ADCw[i])/2.,(tof-tof_p ));

	      if (fabs(tof - tof_pi)<0.7) residptgpi->Fill(pt[i],(tof-tof_pi));
	      if (fabs(tof - tof_k )<0.7) residptgka->Fill(pt[i],(tof-tof_k ));
	      if (fabs(tof - tof_p )<0.7) residptgpr->Fill(pt[i],(tof-tof_p ));

	    }
	    
            //if (pt[i] > pt_cut && indx<216 && fabs(tofe - tof_pi)<0.7 && 
	    //    fabs(tofw - tof_pi)<0.7) {
            if (pt[i] > pt_cut && 
		(flag!=0 || (fabs(tofe - tof_pi)<0.7 && 
			     fabs(tofw - tof_pi)<0.7)) 
		&& pt[i]<0.8) {

	      if (pt[i]<0.8) {
                 de_his[indx]->Fill(tofe - tof_pi);
                 dw_his[indx]->Fill(tofw - tof_pi);
                 d_his [indx]->Fill(tof - tof_pi );
              }
	      
              reside[indx]->Fill((z_trk[i]+L/2),(tofe-tof_pi));
              residw[indx]->Fill((z_trk[i]+L/2),(tofw-tof_pi));

              residqe[indx]->Fill(ADCe[i],(tofe-tof_pi));
              residqw[indx]->Fill(ADCw[i],(tofw-tof_pi));

	      residge->Fill((z_trk[i]+L/2),(tofe-tof_pi));
	      residgw->Fill((z_trk[i]+L/2),(tofw-tof_pi));

              residgqe->Fill(ADCe[i],(tofe-tof_pi));
              residgqw->Fill(ADCw[i],(tofw-tof_pi));

              residgiqe->Fill(1./sqrt(ADCe[i]),(tofe-tof_pi));
              residgiqw->Fill(1./sqrt(ADCw[i]),(tofw-tof_pi));

	      int qbin = (ADCe[i])/400;
	      if (qbin>=0 && qbin<=9) 
		residge_q[qbin]->Fill(z_trk[i]+L/2,(tofe-tof_pi));
	      qbin = (ADCw[i])/400;
	      if (qbin>=0 && qbin<=9) 
		residgw_q[qbin]->Fill(z_trk[i]+L/2,(tofw-tof_pi));
	      
	      int iqbin = 1./sqrt(ADCe[i])/0.01;
	      if (iqbin>=1 && iqbin<=9) 
		residge_iq[iqbin]->Fill(z_trk[i]+L/2,(tofe-tof_pi));
	      iqbin = 1./sqrt(ADCw[i])/0.01;
	      if (iqbin>=1 && iqbin<=9) 
		residgw_iq[iqbin]->Fill(z_trk[i]+L/2,(tofw-tof_pi));
	       
	      int zbin = (z_trk[i]+L/2)/30;
	      if (zbin>=0 && zbin<=9) {
		residgqe_z[zbin]->Fill(ADCe[i],(tofe-tof_pi));
		residgqw_z[zbin]->Fill(ADCw[i],(tofw-tof_pi));
		residgiqe_z[zbin]->Fill(1./sqrt(ADCe[i]),(tofe-tof_pi));
		residgiqw_z[zbin]->Fill(1./sqrt(ADCw[i]),(tofw-tof_pi));
		residiqe_z[indx][zbin]->Fill(1./sqrt(ADCe[i]),(tofe-tof_pi));
		residiqw_z[indx][zbin]->Fill(1./sqrt(ADCw[i]),(tofw-tof_pi));
	      }

            test_charge_e->Fill(tofe - tof_pi, lightspeed_qe);
            test_charge_w->Fill(tofw - tof_pi, lightspeed_qw);

 
              zpos  [indx][nevts[indx]] = z_trk[i];
	      q_east[indx][nevts[indx]] = ADCe[i];
	      q_west[indx][nevts[indx]] = ADCw[i];
              tofpi [indx][nevts[indx]] = tof_pi;      
              tofk  [indx][nevts[indx]] = tof_k;
              tofp  [indx][nevts[indx]] = tof_p;
              tof_e [indx][nevts[indx]] = tofe;
              tof_w [indx][nevts[indx]] = tofw;
              pt_evt[indx][nevts[indx]] = pt[i];
              et0g  [indx][nevts[indx]] = et0g_tmp_clean;
              nevts[indx]++;  

            } // if..pt > pt_cut....  
          } // if..tofe > 0
        } // close loop over low pt tracks (for..nhit)
      } // else if..indx

// Begin: initialize some variables every event ----------------------------

      check_highest_pt = -999;
      pt_hi            =  2.0; // Set the high pt cut
      ibar_trkhi       = -999;
      t0g              = -999;
      nhit             =    0;

      for (int i=0; i<50; i++) {
        pt      [i] = -999;
        ibar_trk[i] = -999;
      }

      for (int i=0; i<NumBar; i++) {
        no_cut_cont[i] =    0;
        cont       [i] =    0;
        cont1      [i] =    0;
        adc_e      [i] = -999;
        adc_w      [i] = -999;
        curv       [i] = -999;
        length     [i] = -999;
        t_e        [i] = -999;
        t_w        [i] = -999;
        trk_p      [i] = -999;
        trk_pt     [i] = -999;
        z          [i] = -999;
        z0bc       [i] = -999;
        pathbar    [i] = -999;
      }

      tof_ntt =    0; // Number of Tracks for T0 Calculation
      for (int tt = 0; tt<200; tt++) {
	tof_nttp[tt]    = 0;
	tof_nttpTrk[tt] = 0;
      }

      jj = 0; // Low pt track counter

// End: initialize some variables every event ------------------------------

    } // if event_counter > num_events

    // Track plots (1)
    check_length->Fill(lightspeed_tofLength);
    check_z	->Fill(lightspeed_z     );
    check_curv  ->Fill(lightspeed_curv  );
    check_phi0  ->Fill(lightspeed_phi0  );

    // Track selection plots (2,3)
    check_d0wrt_before  ->Fill(lightspeed_d0wrt  );
    check_pt_before     ->Fill(lightspeed_pt     );
    check_chi2_before   ->Fill(lightspeed_chi2   );
    check_cotHits_before->Fill(lightspeed_cotHits);
    check_sl8Hits_before->Fill(lightspeed_sl8Hits);
    check_phi_d0_before ->Fill(lightspeed_phi0,lightspeed_d0wrt);

    if (fabs(lightspeed_d0wrt) < 0.076) {
      check_d0wrt_after->Fill(lightspeed_d0wrt);
      check_occ_after_d0wrt_cut->Fill(lightspeed_bar);
      check_phi_d0_after->Fill(lightspeed_phi0,lightspeed_d0wrt);
    }
    if (lightspeed_pt > 0.4) {
      check_pt_after->Fill(lightspeed_pt);
      check_occ_after_pt_cut->Fill(lightspeed_bar);
    }
    if (lightspeed_chi2 < 1000) {
      check_chi2_after->Fill(lightspeed_chi2);
      check_occ_after_chi2_cut->Fill(lightspeed_bar);
    }       
    if (lightspeed_cotHits > 50) {
      check_cotHits_after->Fill(lightspeed_cotHits);
      check_occ_after_cotHits_cut->Fill(lightspeed_bar);
    }
    if (lightspeed_sl8Hits >  4) {
      check_sl8Hits_after->Fill(lightspeed_sl8Hits);
      check_occ_after_sl8Hits_cut->Fill(lightspeed_bar);
    }
    if (lightspeed_inBar == 1) {
      check_occ_after_inBar_cut->Fill(lightspeed_bar);
    }

    //float d0_cut;
   // if (lightspeed_track_type == 2) d0_cut = 0.076;
  //  else d0_cut = 0.025;

    // Fill Variables for T0 Calculation 
    zv_nzv = lightspeed_num_vtx;
    for (int iz = 0; iz<zv_nzv; iz++) zv_zv[iz] = lightspeed_zvtx[iz];
    if (Chi2(lightspeed_te, lightspeed_tw, lightspeed_qe, lightspeed_qw, 
	     lightspeed_z, int(lightspeed_bar), lightspeed_pathinBar) <14. && lightspeed_bar<216) {
      int indx = tof_ntt;
      for (int tt = 0; tt<tof_ntt; tt++)
	if (lightspeed_p==tof_pp[tt]) indx = tt;
      tof_pp[indx]                          = lightspeed_p;
      tof_ctz0bc[indx]                      = lightspeed_z0corr;
      tof_inBar[indx]                       = lightspeed_inBar;
      if (tof_nttp[indx]<4) {
	tof_ttpBar[indx][tof_nttp[indx]]      = int(lightspeed_bar);
	tof_ttpBar[indx][tof_nttp[indx]+1]    = int(lightspeed_bar);
	tof_hbInZ[indx][tof_nttp[indx]/2]     = lightspeed_z;
	tof_hb_s[indx][tof_nttp[indx]/2]      = lightspeed_tofLength;
	tof_hbPathLen[indx][tof_nttp[indx]/2] = lightspeed_pathinBar;
	tof_pulseQ[indx][tof_nttp[indx]]      = lightspeed_qe;
	tof_pulseQ[indx][tof_nttp[indx]+1]    = lightspeed_qw;
	tof_nttp[indx]    +=2;
      }
      tof_nttpTrk[indx] +=2;
      if (indx==tof_ntt) tof_ntt++;
      tp_tRaw[0][int(lightspeed_bar)] = lightspeed_tw;
      tp_tRaw[1][int(lightspeed_bar)] = lightspeed_te;
    }

    // Track selection
    if (lightspeed_num_vtx < 4  && 
//  (fabs(lightspeed_d0wrt) < 0.076) 
//         && (lightspeed_track_type  !=   2)         
//         && (lightspeed_pt          >  0.4)         
//         && (lightspeed_SiPhiHits   >    2)      
//         && (lightspeed_cotAxHits   >   10)      
//         && (lightspeed_cotStHits   >   10)      
//         && (lightspeed_sl8Hits     >    4)
//         && ((lightspeed_inBar == 1) || (flag > 0))
//        ) {     
	Chi2(lightspeed_te, lightspeed_tw, lightspeed_qe, lightspeed_qw, 
	     lightspeed_z, int(lightspeed_bar), lightspeed_pathinBar) <14. &&
//        ((lightspeed_pt<0.8 && lightspeed_pt>0.4 && lightspeed_qe>0. && lightspeed_qw>0.) || 
//	 (flag==2 || flag==0)) && 
   	lightspeed_bar<216) {
      
      ++ntrk_event_after_quality;
      ++cont[bar_number];
      bar_occ_after->Fill(bar_number);

      check_length_qual->Fill(lightspeed_tofLength);
      check_z_qual     ->Fill(lightspeed_z        );
      check_curv_qual  ->Fill(lightspeed_curv     );
      check_phi0_qual  ->Fill(lightspeed_phi0     );

      adc_e [bar_number] = lightspeed_qe;
      adc_w [bar_number] = lightspeed_qw;
      curv  [bar_number] = lightspeed_curv;
      length[bar_number] = lightspeed_tofLength;
      t_e   [bar_number] = lightspeed_te;
      t_w   [bar_number] = lightspeed_tw;
      trk_p [bar_number] = lightspeed_p;
      trk_pt[bar_number] = lightspeed_pt;
      z     [bar_number] = lightspeed_z;
      z0bc  [bar_number] = lightspeed_z0corr;
      pathbar[bar_number]= lightspeed_pathinBar;

      sel_trk_pt_1->Fill(trk_pt[bar_number]);

    } else {

      ++cont1[bar_number];
      not_trk_pt_1->Fill(lightspeed_pt);

    } // track selection

    } // trkPrescale

  } // loop over tracks

  if (flag == 0) printf("Loop pmt resolutions: Tracks %8d (run %d) \n",
	                kentry,lightspeed_run);
  if (flag == 1) printf("Loop z residuals: Tracks %8d (run %d) \n",
	                kentry,lightspeed_run);
  if (flag == 2) printf("Loop z residuals results: Tracks %8d (run %d) \n",
	                kentry,lightspeed_run);

  cout << "\n Number of events = " << event_counter << endl;

  // Residual corrections on z
  if (flag == 1) {
 
    ofstream tze_calib, tzw_calib, tz_calib_log;

    sprintf(filesl,"%s/files/%s_Zres-calib_e.txt",directory,file_name); 
    tze_calib.open(filesl);
    sprintf(filesl,"%s/files/%s_Zres-calib_w.txt",directory,file_name); 
    tzw_calib.open(filesl);
    sprintf(filesl,"%s/files/%s_Zres-calib.log",directory,file_name); 
    tz_calib_log.open(filesl);

//     double InParE[10] = {2.94043e-01, -1.10401572260608871e-01,5.43229e-03, -4.85164859488651179e-01,-1.25796e-01,
// 			 2.94043e-01*20, -1.10401572260608871e-01*20,5.43229e-03*20, -4.85164859488651179e-01*20,
// 			 -1.25796e-01*20};

//     double InParW[10] = {2.94043e-01, -1.10401572260608871e-01,5.43229e-03, -4.85164859488651179e-01,-1.25796e-01,
// 			 2.94043e-01*20, -1.10401572260608871e-01*20,5.43229e-03*20, -4.85164859488651179e-01*20,
// 			 -1.25796e-01*20};
    double InParE[6] = { 1.16336e-01, -4.90497e+00, -5.69027e-04, -7.04072e-06 ,  2.16620e-02, 7.43892e+01};
    double InParW[6] = { 1.16336e-01, -4.90497e+00, -5.69027e-04, -7.04072e-06 ,  2.16620e-02, 7.43892e+01};
    TString ParName[10] = {"PAR0", "PAR1", "PAR2", "PAR3", "PAR4", "PAR5", "PAR6", "PAR7", "PAR8", "PAR9"};
     
    for (int j=0 ;j<NumBar; j++) {
      evts = nevts[j];

      // evts = events that had one (and only one) track hitting the bar #j
      
//       char histoname[50];
//       char pipponame[50];
//       char pipponame1[50];
//       sprintf(histoname,"resid_east_%d",j);
//       TH2F* ehistotemp =  (TH2F*) gDirectory->Get(histoname);
//       sprintf(pipponame,"epippo%d",j);
//       ehistotemp->SetName(pipponame);
//       ehistotemp->FitSlicesY();
//       sprintf(pipponame1,"epippo%d_1",j);
//       TH1F* ehistotemp_1 =  (TH1F*) gDirectory->Get(pipponame1);
//       TF1 *f1 = new TF1("f1","[0]+[1]*exp([2]*x)+[3]*exp([4]*x)",0.,300.);
//       f1->SetParameters( 2.94043e-01, -1.10401572260608871e-01,5.43229e-03, -4.85164859488651179e-01,-1.25796e-01);
//       ehistotemp_1->Fit("f1","","",0.,300.);
//       double InParE[10];
//       for (int h = 0; h<5; h++)
// 	InParE[h] = f1->GetParameter(h);
//       sprintf(histoname,"resid_west_%d",j);
//       TH2F* whistotemp =  (TH2F*) gDirectory->Get(histoname);
//       sprintf(pipponame,"wpippo%d",j);
//       whistotemp->SetName(pipponame);
//       whistotemp->FitSlicesY();
//       sprintf(pipponame1,"wpippo%d_1",j);
//       TH1F* whistotemp_1 =  (TH1F*) gDirectory->Get(pipponame1);
//       f1->SetParameters( 2.94043e-01, -1.10401572260608871e-01,5.43229e-03, -4.85164859488651179e-01,-1.25796e-01);
//       whistotemp_1->Fit("f1","","",0.,300.);
//       double InParW[10];
//       for (int h = 0; h<5; h++)
// 	InParW[h] = f1->GetParameter(h);

      // East side
      ichan = j;  // ichan = photomultiplier
      for (int k=0; k<evts; k++) {
        qzfit[k] = q_east[j][k];
        zp   [k] = zpos[j][k];
        difft[k] = tof_e[j][k] - tofpi[j][k];
      }

      int mbinfite = de_his[j]->GetMaximumBin();
      cfitvalue = de_his[j]->GetBinLowEdge(mbinfite)+de_his[j]->GetBinWidth(mbinfite)/2.;

      for (int p = 0; p<npar; p++) mpar[p] = 0.;
      for (int p = 0; p<npar; p++) empar[p] = 0.;

      gMinuit->mncler();
      gMinuit->SetFCN(fcnZres);
//       for (int p = 0; p<npar-11; p++)
// 	gMinuit->mnparm(p,ParName[p],InParE[p],fabs(InParE[p])/1.,-fabs(InParE[p])*20,fabs(InParE[p])*20,ierr);
//       for (int p = 0; p<npar-11; p++)
// 	gMinuit->mnparm(p+5,ParName[p+5],
// 			InParE[p]*30,fabs(InParE[p]*30)/1.,-fabs(InParE[p]*30)*20,fabs(InParE[p]*30)*20,ierr);
      for (int p = 0; p<6; p++)
	gMinuit->mnparm(p,ParName[p],InParE[p],fabs(InParE[p])/1.,-fabs(InParE[p])*20,fabs(InParE[p])*20,ierr);

      arglis[0] = 1000;
      arglis[1] = 1;
      gMinuit->mnexcm("MIGRAD",arglis,1,ierr);
 
      if (ierr != 0 ) {
	tz_calib_log << ichan << " east side residual corrections on z: this failed because ierr = " << ierr << endl;
// 	for (int p = 0; p<npar-11; p++) mpar[p+4] = InParE[p];
      } // else
// 	for (int p = 0; p<npar-6; p++) gMinuit->GetParameter(p,mpar[p+4],empar[p+4]);
      gMinuit->GetParameter(0,mpar[0],empar[0]);
      gMinuit->GetParameter(1,mpar[4],empar[4]);
      gMinuit->GetParameter(2,mpar[1],empar[1]);
      gMinuit->GetParameter(3,mpar[2],empar[2]);
      gMinuit->GetParameter(4,mpar[5],empar[5]);
      gMinuit->GetParameter(5,mpar[8],empar[8]);


//       TF1 *f1e = new TF1("f1e","pol3",0.,1.);
//       f1e->SetParameters(-6.68537e-01, 6.07106e+01, -1.60566e+03, 1.29766e+04);
      
//       TH1F *PARE[4];
//       PARE[0] = new TH1F("PARE1","",10,0.,300.);
//       PARE[1] = new TH1F("PARE2","",10,0.,300.);
//       PARE[2] = new TH1F("PARE3","",10,0.,300.);
//       PARE[3] = new TH1F("PARE4","",10,0.,300.);

//       for (int zb = 0; zb<10; zb++) {

// 	TH1F *hn_1 = new TH1F("hn_1","hn_1",40,0.015,0.1);
// 	for (int iq = 1; iq<=40; iq++) {
// 	  TH1F *hn = new TH1F("hn","hn",80,-2,2);
// 	  for (int it = 0; it<=81; it++) 
// 	    hn->SetBinContent(iq, residiqe_z[j][zb]->GetBinContent(it+iq*82));
// 	  TF1 *f2 = new TF1("f2","gaus",0.,80.);
// 	  hn->Fit("f2");
// 	  hn_1->SetBinContent(iq, f2->GetParameter(1));
// 	  hn_1->SetBinError  (iq, f2->GetParError (1)); 
// 	}
// 	hn_1->Fit("f1e","","",0.,1./sqrt(300.));
// 	for (int y = 0; y<4; y++) {
// 	  double par  = f1e->GetParameter(y);
// 	  double epar = f1e->GetParError(y);
// 	  PARE[y]->SetBinContent(zb, par);
// 	  PARE[y]->SetBinError(zb, epar);
// 	}

//       }

//       for (int y = 0; y<4; y++) {
// 	PARE[y]->Fit("f1e","","",0.,300.);
// 	for (int p = 0; p<4; p++) {
// 	  mpar[p+y*4] = f1e->GetParameter(p);
// 	  empar[p+y*4] = f1e->GetParError(p);
// 	}
//       }

      tze_calib << j;
      for (int p = 0; p<4; p++) 
	tze_calib << setw(13) << mpar[0+p*4] << setw(13) << empar[0+p*4] 
		  << setw(13) << mpar[1+p*4] << setw(13) << empar[1+p*4] 
		  << setw(13) << mpar[2+p*4] << setw(13) << empar[2+p*4] 
		  << setw(13) << mpar[3+p*4] << setw(13) << empar[3+p*4] << endl;

      // West side
      ichan = j + 216;
      for (int k=0; k<evts; k++) {
        qzfit[k] = q_west[j][k];
        zp   [k] = zpos[j][k];
        difft[k] = tof_w[j][k]-tofpi[j][k];
      }

      int mbinfitw = dw_his[j]->GetMaximumBin();
      cfitvalue = dw_his[j]->GetBinLowEdge(mbinfitw)+dw_his[j]->GetBinWidth(mbinfitw)/2.;

      for (int p = 0; p<npar; p++) mpar[p] = 0.;
      for (int p = 0; p<npar; p++) empar[p] = 0.;

      gMinuit->mncler();
      gMinuit->SetFCN(fcnZres);
//       for (int p = 0; p<npar-11; p++)
// 	gMinuit->mnparm(p,ParName[p],InParW[p],fabs(InParW[p])/1.,-fabs(InParW[p])*20,fabs(InParW[p])*20,ierr);
//       for (int p = 0; p<npar-11; p++)
// 	gMinuit->mnparm(p+5,ParName[p+5],
// 			InParW[p]*30,fabs(InParW[p]*30)/1.,-fabs(InParW[p]*30)*20,fabs(InParW[p]*30)*20,ierr);
      for (int p = 0; p<6; p++)
	gMinuit->mnparm(p,ParName[p],InParW[p],fabs(InParW[p])/1.,-fabs(InParW[p])*20,fabs(InParW[p])*20,ierr);

      arglis[0] = 1000;
      arglis[1] = 1;
      gMinuit->mnexcm("MIGRAD",arglis,1,ierr);
 
      if (ierr != 0) {
	tz_calib_log << ichan << " west side residual corrections on z: this failed because ierr = " << ierr << endl;
// 	for (int p = 0; p<npar-11; p++) mpar[p+4] = InParW[p];
      } // else
// 	for (int p = 0; p<npar-6; p++) gMinuit->GetParameter(p,mpar[p+4],empar[p+4]);
      gMinuit->GetParameter(0,mpar[0],empar[0]);
      gMinuit->GetParameter(1,mpar[4],empar[4]);
      gMinuit->GetParameter(2,mpar[1],empar[1]);
      gMinuit->GetParameter(3,mpar[2],empar[2]);
      gMinuit->GetParameter(4,mpar[5],empar[5]);
      gMinuit->GetParameter(5,mpar[8],empar[8]);

//       TF1 *f1w = new TF1("f1w","pol3",0.,1.);
//       f1w->SetParameters(-6.68537e-01, 6.07106e+01, -1.60566e+03, 1.29766e+04);
      
//       TH1F *PARW[4];
//       PARW[0] = new TH1F("PARW1","",10,0.,300.);
//       PARW[1] = new TH1F("PARW2","",10,0.,300.);
//       PARW[2] = new TH1F("PARW3","",10,0.,300.);
//       PARW[3] = new TH1F("PARW4","",10,0.,300.);

//       for (int zb = 0; zb<10; zb++) {

// 	TH1F *hn_1 = new TH1F("hn_1","hn_1",40,0.015,0.1);
// 	for (int iq = 1; iq<=40; iq++) {
// 	  TH1F *hn = new TH1F("hn","hn",80,-2,2);
// 	  for (int it = 0; it<=81; it++) 
// 	    hn->SetBinContent(iq, residiqw_z[j][zb]->GetBinContent(it+iq*82));
// 	  TF1 *f2 = new TF1("f2","gaus",0.,80.);
// 	  hn->Fit("f2");
// 	  hn_1->SetBinContent(iq, f2->GetParameter(1));
// 	  hn_1->SetBinError  (iq, f2->GetParError (1)); 
// 	}
// 	hn_1->Fit("f1w","","",0.,1./sqrt(300.));
// 	for (int y = 0; y<4; y++) {
// 	  double par  = f1w->GetParameter(y);
// 	  double epar = f1w->GetParError(y);
// 	  PARW[y]->SetBinContent(zb, par);
// 	  PARW[y]->SetBinError(zb, epar);
// 	}

//       }

//       for (int y = 0; y<4; y++) {
// 	PARW[y]->Fit("f1w","","",0.,300.);
// 	for (int p = 0; p<4; p++) {
// 	  mpar[p+y*4] = f1w->GetParameter(p);
// 	  empar[p+y*4] = f1w->GetParError(p);
// 	}
//       }

      
      tzw_calib << j;
      for (int p = 0; p<4; p++) 
	tzw_calib << setw(13) << mpar[0+p*4] << setw(13) << empar[0+p*4] 
		  << setw(13) << mpar[1+p*4] << setw(13) << empar[1+p*4] 
		  << setw(13) << mpar[2+p*4] << setw(13) << empar[2+p*4] 
		  << setw(13) << mpar[3+p*4] << setw(13) << empar[3+p*4] << endl;

    }

    tze_calib.close();
    tzw_calib.close();
  }

  // PMT resolutions
  if (flag == 0) {

    ofstream rese_calib, resw_calib, rese1_calib;

if (true) {

    sprintf(filesl,"%s/files/%s_PmtRes-calib.txt",directory,file_name); 
    rese1_calib.open(filesl);
    sprintf(filesl,"%s/files/%s_PmtRes-calib_e.txt",directory,file_name); 
    rese_calib.open(filesl);
    sprintf(filesl,"%s/files/%s_PmtRes-calib_w.txt",directory,file_name); 
    resw_calib.open(filesl);

} else {

    sprintf(filesl,"%s/files/%s_PmtRes-calib_St1493.txt",directory,file_name); 
    rese1_calib.open(filesl,ios::in);
    sprintf(filesl,"%s/files/%s_PmtRes-calib_e_St1493.txt",directory,file_name); 
    rese_calib.open(filesl,ios::in);
    sprintf(filesl,"%s/files/%s_PmtRes-calib_w_St1493.txt",directory,file_name); 
    resw_calib.open(filesl,ios::in);

} // modified

    for (int j=0 ; j<NumBar; j++) {
      evts = nevts[j];
 
      // evts = events that had one (and only one) track hitting the bar #j

      ichan = j;  // ichan = photomultiplier
      
      for (int k=0; k<evts; k++) {
        zp            [k]    = L/2 - zpos[j][k];
        timeofflight  [k]    = (tof_e[j][k] + tof_w[j][k])/2;
        timeofflight_p[k][0] = tofpi [j][k];
        timeofflight_p[k][1] = tofk  [j][k];
        timeofflight_p[k][2] = tofp  [j][k];
        pt_trk        [k]    = pt_evt[j][k];
        et0_evt       [k]    = et0g  [j][k]; 
      }

      ierr = 0;
      gMinuit_r->mncler();
      gMinuit_r->SetFCN(timeres);
    //gMinuit_r->mnparm(NUM,CHNAM,STVAL,STEP,BND1,BND2,IERRFLG);
      gMinuit_r->mnparm(0,"sigma",0.12,0.05,0.0001,0.5,ierr);
      gMinuit_r->mnparm(1,"lin",0.0,0.01,0,0,ierr);
      gMinuit_r->mnparm(2,"f_pi",0.8,0.05,0,1,ierr);
      gMinuit_r->mnparm(3,"kop",0.5,0.05,0,1,ierr);  // f_k / (f_k + f_p)
      gMinuit_r->mnparm(4,"frac_b",0.8,0.05,0,1,ierr); // signal fraction
      gMinuit_r->mnparm(5,"sigma_b",2,0.1,0,100,ierr); // background sigma
      gMinuit_r->mnparm(6,"b_Ct",0,0.05,0,0,ierr);
      gMinuit_r->FixParameter(1);

      arglis[0] = 1000;
      arglis[1] =  1.0;
      gMinuit_r->mnexcm("MIGRAD",arglis,2,ierr);
     
      gMinuit_r->GetParameter(0,mpar[0],empar[0]);
      gMinuit_r->GetParameter(1,mpar[1],empar[1]);
      gMinuit_r->GetParameter(2,mpar[2],empar[2]);
      gMinuit_r->GetParameter(3,mpar[3],empar[3]);
      gMinuit_r->GetParameter(4,mpar[4],empar[4]);
      gMinuit_r->GetParameter(5,mpar[5],empar[5]);
      gMinuit_r->GetParameter(6,mpar[6],empar[6]);

      // if a bar doesn't work
      if (fabs(mpar[0]-0.12) < 1e-12 && fabs(mpar[1]-0.0004) < 1e-14) {
         mpar[0]  = (-0.999);
	 mpar[1]  = (-0.999);
	 empar[0] = (-0.999);
	 empar[1] = (-0.999);	
      }		     

      // if a bar doesn't converge
      if (ierr != 0) {
        //cout << " Mean PMT resolutions: failed because ierr = " 
             //<< ierr << " for pmt " << ichan << endl;
	mpar[0] -= 1.0;
	mpar[1] -= 1.0;
	empar[0] -= 1.0;
	empar[1] -= 1.0;			
      }

      sigma->Fill(mpar[0]);     

      // rese1_calib corresponds to St%s_PmtRes-calib.txt

      rese1_calib << setw(13) << j       
                  << setw(13) << mpar [0]
                  << setw(13) << empar[0] 
                  << setw(13) << mpar [1]
                  << setw(13) << empar[1] <<endl; 
      
      for (int k=0; k<evts; k++) {
        zp            [k]    = L/2 - zpos[j][k];
        timeofflight  [k]    = tof_e [j][k];
        timeofflight_p[k][0] = tofpi [j][k];
        timeofflight_p[k][1] = tofk  [j][k];
        timeofflight_p[k][2] = tofp  [j][k];
        pt_trk        [k]    = pt_evt[j][k];
        et0_evt       [k]    = et0g  [j][k]; 
      }

      ierr = 0;
      gMinuit_r->mncler();
      gMinuit_r->SetFCN(timeres);
      gMinuit_r->mnparm(0,"sigma",0.12,0.05,0.0001,0.5,ierr);
      gMinuit_r->mnparm(1,"lin",0.0004,0.01,0,1,ierr);
      gMinuit_r->mnparm(2,"f_pi",0.8,0.05,0,1,ierr);
      gMinuit_r->mnparm(3,"kop",0.5,0.05,0,1,ierr);
      gMinuit_r->mnparm(4,"frac_b",0.8,0.05,0,1,ierr);
      gMinuit_r->mnparm(5,"sigma_b",2,0.1,0,100,ierr);
      gMinuit_r->mnparm(6,"b_Ct",0,0.05,0,0,ierr);

      arglis[0] = 1000;
      arglis[1] =  1.0;
      gMinuit_r->mnexcm("MIGRAD",arglis,2,ierr);
     
      gMinuit_r->GetParameter(0,mpar[0],empar[0]);
      gMinuit_r->GetParameter(1,mpar[1],empar[1]);
      gMinuit_r->GetParameter(2,mpar[2],empar[2]);
      gMinuit_r->GetParameter(3,mpar[3],empar[3]);
      gMinuit_r->GetParameter(4,mpar[4],empar[4]);
      gMinuit_r->GetParameter(5,mpar[5],empar[5]);
      gMinuit_r->GetParameter(6,mpar[6],empar[6]);

      // if a bar doesn't work
      if (fabs(mpar[0]-0.12) < 1e-12 && fabs(mpar[1]-0.0004) < 1e-14) {
         cout << " Warning!!! PMT "<< j << " (east) doesn't work "  
	      << evts << endl;
         mpar[0]  = (-0.999);
	 mpar[1]  = (-0.999);
	 empar[0] = (-0.999);
	 empar[1] = (-0.999);	
      }		     

      // if a bar doesn't converge
      if (ierr != 0) {
        cout << " East PMT resolutions: failed because ierr = " 
             << ierr << " for pmt " << ichan << endl;
	mpar[0] -= 1.0;
	mpar[1] -= 1.0;
	empar[0] -= 1.0;
	empar[1] -= 1.0;			

      }

      rese_calib << setw(13) << j
                 << setw(13) << mpar [0]
                 << setw(13) << empar[0] 
                 << setw(13) << mpar [1]
                 << setw(13) << empar[1] <<endl; 

      sigme->Fill(mpar[0]);
      depe ->Fill(mpar[1]);
 
      for (int i=0; i<evts; i++) 
        diffe_lik[j]->Fill((timeofflight[i] - timeofflight_p[i][0]),timeres2(mpar,i));

      // West side       
      ichan = j + 216;
      for (int k=0; k<evts; k++) {
        zp            [k]    = L/2+zpos[j][k];
        timeofflight  [k]    = tof_w[j][k];
        timeofflight_p[k][0] = tofpi[j][k];
        timeofflight_p[k][1] = tofk[j][k];
        timeofflight_p[k][2] = tofp[j][k];
        pt_trk        [k]    = pt_evt[j][k];    
        et0_evt       [k]    = et0g[j][k];
      }

      ierr = 0;
      gMinuit_r->mncler();
      gMinuit_r->SetFCN(timeres);
      gMinuit_r->mnparm(0,"sigma",0.12,0.05,0.0001,0.5,ierr);
      gMinuit_r->mnparm(1,"lin",0.0004,0.01,0,1,ierr);
      gMinuit_r->mnparm(2,"f_pi",0.8,0.05,0,1,ierr);
      gMinuit_r->mnparm(3,"kop",0.5,0.05,0,1,ierr);
      gMinuit_r->mnparm(4,"frac_b",0.8,0.05,0,1,ierr);
      gMinuit_r->mnparm(5,"sigma_b",2,0.1,0,100,ierr);
      gMinuit_r->mnparm(6,"b_Ct",0,0.05,0,0,ierr);

      arglis[0] = 1000;
      arglis[1] =  1.0;
      gMinuit_r->mnexcm("MIGRAD",arglis,2,ierr);
      
      gMinuit_r->GetParameter(0,mpar[0],empar[0]);
      gMinuit_r->GetParameter(1,mpar[1],empar[1]);
      gMinuit_r->GetParameter(2,mpar[2],empar[2]);
      gMinuit_r->GetParameter(3,mpar[3],empar[3]);
      gMinuit_r->GetParameter(4,mpar[4],empar[4]);
      gMinuit_r->GetParameter(5,mpar[5],empar[5]);
      gMinuit_r->GetParameter(6,mpar[6],empar[6]);

      // if a bar doesn't work
      if (fabs(mpar[0]-0.12) < 1e-12 && fabs(mpar[1]-0.0004) < 1e-14) {
         cout << " Warning!!! PMT "<< j << " (west) doesn't work "
	      << evts << endl;
         mpar[0]  = (-0.999);
	 mpar[1]  = (-0.999);
	 empar[0] = (-0.999);
	 empar[1] = (-0.999);	
      }		     

      // if a bar doesn't converge
      if (ierr != 0) {
        cout << " West PMT resolutions: failed because ierr = " 
             << ierr << " for pmt " << ichan << endl;
	mpar[0] -= 1.0;
	mpar[1] -= 1.0;
	empar[0] -= 1.0;
	empar[1] -= 1.0;			
      }

      resw_calib << setw(13) << j       
                 << setw(13) << mpar [0]
                 << setw(13) << empar[0] 
                 << setw(13) << mpar [1]
                 << setw(13) << empar[1] <<endl; 

      sigmw->Fill(mpar[0]);
      depw ->Fill(mpar[1]);

      for (int i=0; i<evts; i++)
        diffw_lik[j]->Fill(timeofflight[i] - timeofflight_p[i][0],timeres2(mpar,i));
  
    } // for..NumBar

    rese_calib.close();
    resw_calib.close();
    rese1_calib.close();

  } // if..flag == 0

  // New 031203
//event_t0_t0Error_File.close();
//   cout << " Number of events with t0 results = " << countingEvents << endl;

  hfile->Write();
  hfile->Close();

} // speed_calib::LoopRes

//--------------------------------------------------------------------------
// Residual corrections on z and PMT resolutions
//--------------------------------------------------------------------------

void speed_calib::LoopVal(int     NumEvt,
                          char   *file_name,
                          char   *directory)
{
  if (fChain == 0) return;
  Int_t nentries = Int_t(fChain->GetEntries());
  Int_t nbytes = 0, nb = 0;

  cout << "\n Loop validation: Number of tracks = " << nentries << " \n" << endl;

  // modified _St1493
  sprintf(filesl,"%s/root/%s_valid-calib_newplots_zres.root",directory,file_name);

  TFile *hfile = new TFile(filesl,"recreate");

  gStyle->SetOptFit(111);

  // Book histos
  TH1F *t0    = new TH1F("t0",   "t0",  600,-30.0,30    );
  TH1D* ZRES = new TH1D("ZRES","ZRES",4000,-10.,10.); 

  TH1D* PTDIS = new TH1D("PTDIS","PTDIS",500,0.,100.);
  TH1D* PDIS = new TH1D("PDIS","PDIS",50,0.,10.);

  TH1D* QDIS1 = new TH1D("QDIS1","QDIS1",400,0.,4000.);
  TH1D* QNDIS1 = new TH1D("QNDIS1","QNDIS1",400,0,4000);
  TH1D* DTP1 = new TH1D("DTP1","DTP1",4000,-10.,10.);
  TH1D* DTK1 = new TH1D("DTK1","DTK1",4000,-10.,10.);
  TH2D* TQP1 = new TH2D("TQP1","TQP1",20,-0.4,0.4,100,0,4000);
  TH2D* TQK1 = new TH2D("TQK1","TQK1",20,-0.4,0.4,100,0,4000);
  TH2D* TQR1 = new TH2D("TQR1","TQR1",20,-0.4,0.4,100,0,4000);
  TH2D* TQP1N = new TH2D("TQP1N","TQP1N",20,-0.4,0.4,100,0,4);
  TH2D* TQK1N = new TH2D("TQK1N","TQK1N",20,-0.4,0.4,100,0,4);
  TH2D* TQR1N = new TH2D("TQR1N","TQR1N",20,-0.4,0.4,100,0,4);

  TH1D* QDIS2 = new TH1D("QDIS2","QDIS2",400,0.,4000.);
  TH1D* QNDIS2 = new TH1D("QNDIS2","QNDIS2",400,0,4000);
  TH1D* DTP2 = new TH1D("DTP2","DTP2",4000,-10.,10.);
  TH1D* DTK2 = new TH1D("DTK2","DTK2",4000,-10.,10.);
  TH2D* TQP2 = new TH2D("TQP2","TQP2",20,-0.4,0.4,100,0,4000);
  TH2D* TQK2 = new TH2D("TQK2","TQK2",20,-0.4,0.4,100,0,4000);
  TH2D* TQR2 = new TH2D("TQR2","TQR2",20,-0.4,0.4,100,0,4000);
  TH2D* TQP2N = new TH2D("TQP2N","TQP2N",20,-0.4,0.4,100,0,4);
  TH2D* TQK2N = new TH2D("TQK2N","TQK2N",20,-0.4,0.4,100,0,4);
  TH2D* TQR2N = new TH2D("TQR2N","TQR2N",20,-0.4,0.4,100,0,4);

  TH1D* QDIS3 = new TH1D("QDIS3","QDIS3",400,0.,4000.);
  TH1D* QNDIS3 = new TH1D("QNDIS3","QNDIS3",400,0,4000);
  TH1D* DTP3 = new TH1D("DTP3","DTP3",4000,-10.,10.);
  TH1D* DTK3 = new TH1D("DTK3","DTK3",4000,-10.,10.);
  TH2D* TQP3 = new TH2D("TQP3","TQP3",20,-0.4,0.4,100,0,4000);
  TH2D* TQK3 = new TH2D("TQK3","TQK3",20,-0.4,0.4,100,0,4000);
  TH2D* TQR3 = new TH2D("TQR3","TQR3",20,-0.4,0.4,100,0,4000);
  TH2D* TQP3N = new TH2D("TQP3N","TQP3N",20,-0.4,0.4,100,0,4);
  TH2D* TQK3N = new TH2D("TQK3N","TQK3N",20,-0.4,0.4,100,0,4);
  TH2D* TQR3N = new TH2D("TQR3N","TQR3N",20,-0.4,0.4,100,0,4);

  TH1D* QDIS4 = new TH1D("QDIS4","QDIS4",400,0.,4000.);
  TH1D* QNDIS4 = new TH1D("QNDIS4","QNDIS4",400,0,4000);
  TH1D* DTP4 = new TH1D("DTP4","DTP4",4000,-10.,10.);
  TH1D* DTK4 = new TH1D("DTK4","DTK4",4000,-10.,10.);
  TH2D* TQP4 = new TH2D("TQP4","TQP4",20,-0.4,0.4,100,0,4000);
  TH2D* TQK4 = new TH2D("TQK4","TQK4",20,-0.4,0.4,100,0,4000);
  TH2D* TQR4 = new TH2D("TQR4","TQR4",20,-0.4,0.4,100,0,4000);
  TH2D* TQP4N = new TH2D("TQP4N","TQP4N",20,-0.4,0.4,100,0,4);
  TH2D* TQK4N = new TH2D("TQK4N","TQK4N",20,-0.4,0.4,100,0,4);
  TH2D* TQR4N = new TH2D("TQR4N","TQR4N",20,-0.4,0.4,100,0,4);

  TH1D* QDIS5 = new TH1D("QDIS5","QDIS5",400,0.,4000.);
  TH1D* QNDIS5 = new TH1D("QNDIS5","QNDIS5",400,0,4000);
  TH1D* DTP5 = new TH1D("DTP5","DTP5",4000,-10.,10.);
  TH1D* DTK5 = new TH1D("DTK5","DTK5",4000,-10.,10.);
  TH2D* TQP5 = new TH2D("TQP5","TQP5",20,-0.4,0.4,100,0,4000);
  TH2D* TQK5 = new TH2D("TQK5","TQK5",20,-0.4,0.4,100,0,4000);
  TH2D* TQR5 = new TH2D("TQR5","TQR5",20,-0.4,0.4,100,0,4000);
  TH2D* TQP5N = new TH2D("TQP5N","TQP5N",20,-0.4,0.4,100,0,4);
  TH2D* TQK5N = new TH2D("TQK5N","TQK5N",20,-0.4,0.4,100,0,4);
  TH2D* TQR5N = new TH2D("TQR5N","TQR5N",20,-0.4,0.4,100,0,4);

  TH1D* ZZZ = new TH1D("ZZZ","ZZZ",1000,-150,150);
  TH2D* QZ  = new TH2D("QZ","QZ",100,0.,4000.,100,-150,150);
  TH2D* PZ  = new TH2D("PZ","PZ",50,0.,10.,100,-150,150);
  TH2D* ZQ0  = new TH2D("ZQ0","ZQ0",100,-150,150,100,0.,4000.);
  TH2D* PQ0  = new TH2D("PQ0","PQ0",50,0,10,100,0.,4000.);
  TH2D* ZQ1  = new TH2D("ZQ1","ZQ1",100,-150,150,100,0.,4000.);
  TH2D* PQ1  = new TH2D("PQ1","PQ1",50,0,10,100,0.,4000.);
  TH2D* PTQ0  = new TH2D("PTQ0","PTQ0",50,0.,10.,100,0.,4000.);
  TH2D* PTQ1  = new TH2D("PTQ1","PTQ1",50,0.,10.,100,0.,4000.);
  TH2D* ZQN = new TH2D("ZQN","ZQN",100,-150,150,100,0.,4.);
  TH2D* QNZ = new TH2D("QNZ","QNZ",100,0.,4.,100,-150,150);
  TH2D* QMZ = new TH2D("QMZ","QMZ",100,0.,4.,100,-150,150);
  TH2D* DTZ = new TH2D("DTZ","DTZ",800,-0.4,0.4,100,-150,150);
  TH2D* DTZE = new TH2D("DTZE","DTZE",800,-0.4,0.4,100,-150,150);
  TH2D* DTZW = new TH2D("DTZW","DTZW",800,-0.4,0.4,100,-150,150);
  TH2D* EWZ  = new TH2D("EWZ", "EWZ" ,100,-150,150,800,-0.4,0.4);
  TH2D* EWQ  = new TH2D("EWQ", "EWQ" ,100,0.,4000.,800,-0.4,0.4);
  TH2D* EWP  = new TH2D("EWP", "EWP" ,100,0.,10,800,-0.4,0.4);
  TH2D* EVW  = new TH2D("EVW", "EVW" ,160,-0.4,0.4,160,-0.4,0.4);
  TH2D* QEVQW= new TH2D("QEVQW", "QEVQW",100,0.,4000.,100,0.,4000.);

  TH2D* TQP5M = new TH2D("TQP5M","TQP5M",20,-0.4,0.4,100,0,4);
  TH2D* TQP4M = new TH2D("TQP4M","TQP4M",20,-0.4,0.4,100,0,4);
  TH2D* TQP3M = new TH2D("TQP3M","TQP3M",20,-0.4,0.4,100,0,4);
  TH2D* TQP2M = new TH2D("TQP2M","TQP2M",20,-0.4,0.4,100,0,4);
  TH2D* TQP1M = new TH2D("TQP1M","TQP1M",20,-0.4,0.4,100,0,4);

  TH1D* DTM1 = new TH1D("DTM1","DTM1",4000,-10.,10.);
  TH1D* DTM2 = new TH1D("DTM2","DTM2",4000,-10.,10.);
  TH1D* DTM3 = new TH1D("DTM3","DTM3",4000,-10.,10.);
  TH1D* DTM4 = new TH1D("DTM4","DTM4",4000,-10.,10.);
  TH1D* DTM5 = new TH1D("DTM5","DTM5",4000,-10.,10.);

  TH2D* TQM5M = new TH2D("TQM5M","TQM5M",20,-0.4,0.4,100,0,4);
  TH2D* TQM4M = new TH2D("TQM4M","TQM4M",20,-0.4,0.4,100,0,4);
  TH2D* TQM3M = new TH2D("TQM3M","TQM3M",20,-0.4,0.4,100,0,4);
  TH2D* TQM2M = new TH2D("TQM2M","TQM2M",20,-0.4,0.4,100,0,4);
  TH2D* TQM1M = new TH2D("TQM1M","TQM1M",20,-0.4,0.4,100,0,4);

  TH2D* TQM5N = new TH2D("TQM5N","TQM5N",20,-0.4,0.4,100,0,4);
  TH2D* TQM4N = new TH2D("TQM4N","TQM4N",20,-0.4,0.4,100,0,4);
  TH2D* TQM3N = new TH2D("TQM3N","TQM3N",20,-0.4,0.4,100,0,4);
  TH2D* TQM2N = new TH2D("TQM2N","TQM2N",20,-0.4,0.4,100,0,4);
  TH2D* TQM1N = new TH2D("TQM1N","TQM1N",20,-0.4,0.4,100,0,4);

  TH2D* TQP1W = new TH2D("TQP1W","TQP1W",20,-0.4,0.4,100,0,4000);
  TH2D* TQP2W = new TH2D("TQP2W","TQP2W",20,-0.4,0.4,100,0,4000);
  TH2D* TQP3W = new TH2D("TQP3W","TQP3W",20,-0.4,0.4,100,0,4000);
  TH2D* TQP4W = new TH2D("TQP4W","TQP4W",20,-0.4,0.4,100,0,4000);
  TH2D* TQP5W = new TH2D("TQP5W","TQP5W",20,-0.4,0.4,100,0,4000);

  TH2D* TQPE[5][3];
  TH2D* TPQE[5][3];
  TH2D* TQZE[5][3];
  TH2D* TZQE[5][3];
  TH2D* TQME[5][3];
  TH2D* TPZE[5][3];
  TH2D* TZPE[5][3];
  TH2D* TIQZE[5][3];
  TH2D* TIRQZE[5][3];
  
  for (int i=0; i<5; i++) {
    for (int j=0; j<3; j++) {
      sprintf(titlesl,"TQPE%dpt%d",i,j);
      TQPE[i][j]  = new TH2D(titlesl,titlesl,160,-0.4,0.4,50,0,4000);
      sprintf(titlesl,"TPQE%dpt%d",i,j);
      TPQE[i][j]  = new TH2D(titlesl,titlesl,160,-0.4,0.4,25,0,10);
      sprintf(titlesl,"TQZE%dpt%d",i,j);
      TQZE[i][j]  = new TH2D(titlesl,titlesl,160,-0.4,0.4,50,0,4000);
      sprintf(titlesl,"TZQE%dpt%d",i,j);
      TZQE[i][j]  = new TH2D(titlesl,titlesl,160,-0.4,0.4,30,-150.,150.);
      sprintf(titlesl,"TQME%dpt%d",i,j);
      TQME[i][j]  = new TH2D(titlesl,titlesl,160,-0.4,0.4,50,0,4000);
      sprintf(titlesl,"TPZE%dpt%d",i,j);
      TPZE[i][j]  = new TH2D(titlesl,titlesl,160,-0.4,0.4,25,0,10);
      sprintf(titlesl,"TZPE%dpt%d",i,j);
      TZPE[i][j]  = new TH2D(titlesl,titlesl,160,-0.4,0.4,30,-150.,150.);
      sprintf(titlesl,"TIQZE%dpt%d",i,j);
      TIQZE[i][j]  = new TH2D(titlesl,titlesl,160,-0.4,0.4,50,0.00025,0.01);
      sprintf(titlesl,"TIRQZE%dpt%d",i,j);
      TIRQZE[i][j]  = new TH2D(titlesl,titlesl,160,-0.4,0.4,50,0.01,0.1);
    }
  }

 

  // Read parameters
  FillParameters(directory, file_name, 5);
//   ifstream att_result; 
//   ifstream speed_result;
//   ifstream twe_result, t0e_result; 
//   ifstream tww_result, t0w_result;

// if (true) {

//   sprintf(filesl,"%s/files/%s_timew-calib_e.txt",directory,file_name);
//   twe_result.open(filesl,ios::in);   
//   sprintf(filesl,"%s/files/%s_timew-calib_w.txt",directory,file_name);
//   tww_result.open(filesl,ios::in);
//   sprintf(filesl,"%s/files/%s_T0-calib_e.txt",directory,file_name);
//   t0e_result.open(filesl,ios::in);   
//   sprintf(filesl,"%s/files/%s_T0-calib_w.txt",directory,file_name);
//   t0w_result.open(filesl,ios::in);
//   sprintf(filesl,"%s/files/%s_spdtw-calib.txt",directory,file_name);
//   speed_result.open(filesl,ios::in);

// } else {

//   sprintf(filesl,"/cdf/data03/s1/tof/rocio/Calibration/files/St1493_timew-calib_e.txt");
//   twe_result.open(filesl,ios::in);   
//   sprintf(filesl,"/cdf/data03/s1/tof/rocio/Calibration/files/St1493_timew-calib_w.txt");
//   tww_result.open(filesl,ios::in);
//   sprintf(filesl,"/cdf/data03/s1/tof/rocio/Calibration/files/St1493_T0-calib_e.txt");
//   t0e_result.open(filesl,ios::in);   
//   sprintf(filesl,"/cdf/data03/s1/tof/rocio/Calibration/files/St1493_T0-calib_w.txt");
//   t0w_result.open(filesl,ios::in);
//   sprintf(filesl,"/cdf/data03/s1/tof/rocio/Calibration/files/St1493_spdtw-calib_p.txt");
//   speed_result.open(filesl,ios::in);

// } // modified

//   for (int i=0; i<9; i++) {
//     datatwe[i] = 0.0;
//     datatww[i] = 0.0;
//     datatze[i] = 0.0;
//     datatzw[i] = 0.0;
//     datasp [i] = 0.0;
//     if (i < 5) datat0e[i] = 0.0;
//     if (i < 5) datat0w[i] = 0.0;
//   }

//   for (int j=0; j<NumBar; j++) {
//     spdlight  [j][0] = 0.0;
//     spdlight  [j][1] = 0.0;
//     twalk_e   [j]    = 0.0;
//     twalk_w   [j]    = 0.0;
//     spdlight_e[j]    = 0.0;
//     spdlight_w[j]    = 0.0;
//     t0_off_e  [j]    = 0.0;
//     t0_off_w  [j]    = 0.0;
//     for (int k=0; k<4; k++) {tz_e[j][k] = 0.0; tz_w[j][k] = 0.0;}
//   }

//   for (int j=0; j<NumBar; j++) {
//     for (int i=0; i<9; i++) twe_result   >> datatwe[i];
//     for (int i=0; i<9; i++) tww_result   >> datatww[i];
//     for (int i=0; i<5; i++) t0e_result   >> datat0e[i];
//     for (int i=0; i<5; i++) t0w_result   >> datat0w[i];
//     for (int i=0; i<9; i++) speed_result >> datasp [i];

//     twalk_e   [j]    = datatwe[1];
//     twalk_w   [j]    = datatww[1];    
//     //t0_off_e  [j]    = datat0e[1]; // offset if s east and west are different
//     //t0_off_w  [j]    = datat0w[1]; // offset if s east and west are different
//     t0_off_e  [j]    = datat0e[3];   // offset if s east and west are the same
//     t0_off_w  [j]    = datat0w[3];   // offset if s east and west are the same

//     // modified
//     if (true) {spdlight[j][0] = datasp [1];} else {spdlight[j][0] = 2.0/datasp[1];}

//     spdlight  [j][1] = datasp [3];
//     //spdlight_e[j]    = datat0e[2]; // s if s east and west are different
//     //spdlight_w[j]    = datat0w[2]; // s if s east and west are different
//     spdlight_e[j]    = datat0e[4];   // s if s east and west are the same
//     spdlight_w[j]    = datat0w[4];   // s if s east and west are the same


// //   int nevts[NumBar];

// //   double *zpos[NumBar], *pt_evt[NumBar];
// //   double *tof_e[NumBar], *tof_w[NumBar], *et0g[NumBar];
// //   double *tofk[NumBar], *tofp[NumBar], *tofpi[NumBar];

// //     nevts [j] = 0;
// //     zpos  [j] = new double[maxevt];
// //     tof_e [j] = new double[maxevt];
// //     tof_w [j] = new double[maxevt];
// //     tofk  [j] = new double[maxevt];
// //     tofp  [j] = new double[maxevt];
// //     tofpi [j] = new double[maxevt];
// //     pt_evt[j] = new double[maxevt];
// //     et0g  [j] = new double[maxevt];

// //     for (int k=0; k<maxevt; k++) {
// //       zpos  [j][k] = 0.0;
// //       tof_e [j][k] = 0.0;
// //       tof_w [j][k] = 0.0;
// //       tofk  [j][k] = 0.0;
// //       tofp  [j][k] = 0.0;
// //       tofpi [j][k] = 0.0;
// //       pt_evt[j][k] = 0.0;
// //       et0g  [j][k] = 0.0;
// //     }

//   } // for..NumBar

//   twe_result.close();
//   tww_result.close();
//   t0e_result.close();
//   t0w_result.close();
//   speed_result.close();

  // New 031203
  //int currentEventForT0Count = -999;
  //int countingEvents         =    0;
//ofstream event_t0_t0Error_File;
//sprintf(filesl,"%s/files/%s_event-t0.txt",directory,file_name); 
//event_t0_t0Error_File.open(filesl,ios::in);

//   ifstream tze_result, tzw_result;

// if (true) {

//     sprintf(filesl,"%s/files/%s_Zres-calib_e.txt",directory,file_name);
//     tze_result.open(filesl,ios::in);
//     sprintf(filesl,"%s/files/%s_Zres-calib_w.txt",directory,file_name);
//     tzw_result.open(filesl,ios::in);

// } else {

//     sprintf(filesl,"/cdf/data03/s1/tof/rocio/Calibration/files/St1493_Zres-calib_e.txt");
//     tze_result.open(filesl,ios::in);
//     sprintf(filesl,"/cdf/data03/s1/tof/rocio/Calibration/files/St1493_Zres-calib_w.txt");
//     tzw_result.open(filesl,ios::in);

// } // modified

//     for (int j=0; j<NumBar; j++) {   
//       for (int i=0; i<9; i++) tze_result >> datatze[i];
//       for (int i=0; i<9; i++) tzw_result >> datatzw[i];

//       tz_e[j][0] = datatze[1];
//       tz_e[j][1] = datatze[3];
//       tz_e[j][2] = datatze[5];
//       tz_e[j][3] = datatze[7];
//       tz_w[j][0] = datatzw[1];
//       tz_w[j][1] = datatzw[3];
//       tz_w[j][2] = datatzw[5];
//       tz_w[j][3] = datatzw[7];
 
//     } // for..NumBar

//     tze_result.close();
//     tzw_result.close();

  //int    ierr     =    0;
  //int    ierre    =    0;
  //double t0i      = -999;
  //double et0g_tmp = -999;

  jj = 0; // Low pt track counter

// End: initialize some variables every event ------------------------------

  int current_event  = -999;
  //int previous_event = -999;
  int event_counter  =    0;
  int num_events     =    1;

  int ntrk_event                 = 0;
  int ntrk_event_after_quality   = 0;
  //int events_with_nhit_ge_ntrkt0 = 0;
  //int events_with_nhit_lt_ntrkt0 = 0;

  //int total_track = 0;     // Initial number of tracks used to calculate t0

  // Fill parameters for track corrections
  FillParameters(directory, file_name, 2);

  Int_t kentry = 0;
  
  // Loop over tracks
  if (NumEvt > nentries) NumEvt = nentries;
  for (Int_t jentry=0; jentry<NumEvt; jentry++) {
    Int_t ientry = LoadTree(jentry);
    nb = fChain->GetEntry(jentry);
    nbytes += nb;

    if ((jentry%trkPrescale)<10)
    {
    // Prescaled track count: kentry
    ++kentry;

    int bar_number = int(lightspeed_bar);

    /*
    if (kentry%200000 == 0) {
      printf("Loop validation:track %8d (run %d) \n",kentry,lightspeed_run);
    }
    */
    
    ++no_cut_cont[bar_number];

    if (lightspeed_event != current_event) {
      //previous_event = current_event;
      current_event = lightspeed_event;
      ++event_counter;

      ntrk_event               = 0;
      ntrk_event_after_quality = 0;
    }

    ntrk_event++;

    if (event_counter > num_events) {

      num_events = event_counter;

      // Loop over bars
      for (int i=0; i<NumBar; i++) {
//         if (cont[i] != 1 || cont1[i] != 0) continue;


//         if (adc_e[i] > 0) tofe = t_e[i] + twalk_e[i]/sqrt(adc_e[i]); else tofe = t_e[i];
//         if (adc_w[i] > 0) tofw = t_w[i] + twalk_w[i]/sqrt(adc_w[i]); else tofw = t_w[i];
//         if ((tofe > 0) && (tofw > 0)) diff = tofe - tofw; else diff = -999;

//         tofe = -999;
//         tofw = -999;

//         tofe = t_e[i] + (L/2 - z[i])/spdlight[i][0];
//         tofw = t_w[i] + (L/2 + z[i])/spdlight[i][0];

//         if (adc_e[i] > 0 && adc_w[i] > 0) {
//            tofe += twalk_e[i]/sqrt(adc_e[i]);
//            tofw += twalk_w[i]/sqrt(adc_w[i]);
//         }
 
//         double diff_for_zres = t_e[i] - t_w[i];
//         if (adc_e[i] > 0 && adc_w[i] > 0) {
//            diff_for_zres += twalk_e[i]/sqrt(adc_e[i]);
//            diff_for_zres -= twalk_w[i]/sqrt(adc_w[i]);
//         }
//         diff_for_zres = (diff_for_zres - spdlight[i][1])*spdlight[i][0]/2.0 - z[i];


//         tofe -= t0_off_e[i];
//         tofw -= t0_off_w[i];

//         tofe +=  tz_e[i][0]
// 	  +tz_e[i][1]*(z[i]+L/2)
// 	  +tz_e[i][2]*(z[i]+L/2)*(z[i]+L/2)
// 	  +tz_e[i][3]*(z[i]+L/2)*(z[i]+L/2)*(z[i]+L/2);
	
//         tofw +=  tz_w[i][0]
// 	  +tz_w[i][1]*(z[i]+L/2)
// 	  +tz_w[i][2]*(z[i]+L/2)*(z[i]+L/2)
// 	  +tz_w[i][3]*(z[i]+L/2)*(z[i]+L/2)*(z[i]+L/2);
	
//         diff = tofe - tofw;
	
// 	float cut = 0.5;
//         if (fabs(diff) > cut) continue;          // pmt resolutions (cut = 0.5)

	float pt_cut = 0.; int maxtrk = 50;
	if ((trk_pt[i] > pt_cut) && (jj < maxtrk)) { 

          curv_trk[jj] = curv  [i];
	  pt      [jj] = trk_pt[i];
          p_trk   [jj] = trk_p [i];
          z_trk   [jj] = z     [i];      
          path_trk[jj] = length[i]; 
          ibar_trk[jj] = i;  
          ADCe    [jj] = adc_e [i];
          TDCe    [jj] = t_e   [i];
          ADCw    [jj] = adc_w [i];
          TDCw    [jj] = t_w   [i];
          z0bc_trk[jj] = z0bc  [i];
          pathbar_trk[jj] = pathbar  [i];
          inbar_trk[jj]= inbar  [i];

          if (pt[jj] > check_highest_pt) check_highest_pt = pt[jj];

          jj++;
	}

        // Store any high pt track for initial t0 calculation
       	if (trk_pt[i] > pt_hi) {
          pt_hi      = trk_pt[i];
          p_trkhi    = trk_p [i];
          path_trkhi = length[i];
          ibar_trkhi = i;
	}
      } // loop over bars

      nhit = jj;

      // Initial t0 calculation
      //indx = ibar_trkhi; this was for old t0 calculation
      indx = nhit;
      
      if (indx < 0) {
      }	
      else {

	// New T0 Calculation
	double CorrectedT0[20], ErrorCorrT0[20]; 
	for (int nv = 0; nv<zv_nzv; nv++) {
	  CorrectedT0[nv] = -999.;
	  ErrorCorrT0[nv] = -999.;
	}
	_skip_p = -999;
	for (int nv = 0; nv<zv_nzv; nv++) {
	  CorrectedT0[nv] = speed_calib::CalculateT0(nv, &ErrorCorrT0[nv]); 
	  t0->Fill(CorrectedT0[nv]);
	}

        for (int kkk=0; kkk<nhit; kkk++ ) {
	
	  double CorrectedTrackT0 = -999.;
	  //double ErrorCorrTrackT0 = -999.;
	  for (int nv = 0; nv<zv_nzv; nv++) 
	    if (fabs(z0bc_trk[kkk]-zv_zv[nv])<2.0) {
	      CorrectedTrackT0 = CorrectedT0[nv];
	      //ErrorCorrTrackT0 = ErrorCorrT0[nv];
// 	      _skip_p = p_trk[kkk];
// 	      CorrectedTrackT0 = speed_calib::CalculateT0(nv, &ErrorCorrTrackT0); 
	    }
	  if (CorrectedTrackT0==-999.) continue;
	  int i = kkk;

	  float t0g_clean = CorrectedTrackT0;
	  //float et0g_tmp_clean = ErrorCorrTrackT0;
	
          indx = ibar_trk[i];

          // Predicted TOF for different particles
	  tof_m = getTofSL(path_trk[i], p_trk[i], 0.105);
	  tof_pi = getTofSL(path_trk[i], p_trk[i], pimass);
	  tof_k  = getTofSL(path_trk[i], p_trk[i], kmass );
	  tof_p  = getTofSL(path_trk[i], p_trk[i], pmass );

          // TOF calculated using east pmt's
	  tofe = TimeCorr(TDCe[i], z_trk[i], ADCe[i], indx);
//           tofe  = t0_off_e[indx] - TDCe[i];
//           tofe -= (L/2-z_trk[i])/spdlight[indx][0];
//           if (ADCe[i] > 0) tofe -= twalk_e[indx]/sqrt(ADCe[i]);

// 	  ZRES->Fill(tz_e[indx][0]
// 		     +tz_e[indx][1]*(z_trk[i]+L/2)
// 		     +tz_e[indx][2]*(z_trk[i]+L/2)*(z_trk[i]+L/2)
// 		     +tz_e[indx][3]*(z_trk[i]+L/2)*(z_trk[i]+L/2)*(z_trk[i]+L/2));

//           tofe -= ( tz_e[indx][0]
// 		    +tz_e[indx][1]*(z_trk[i]+L/2)
// 		    +tz_e[indx][2]*(z_trk[i]+L/2)*(z_trk[i]+L/2)
// 		    +tz_e[indx][3]*(z_trk[i]+L/2)*(z_trk[i]+L/2)*(z_trk[i]+L/2));
	  
          tofe -= t0g_clean; 

          // TOF calculated using west pmt's
	  tofw = TimeCorr(TDCw[i], z_trk[i], ADCw[i], indx+216);
//           tofw  = t0_off_w[indx] - TDCw[i];
//           tofw -= (L/2+z_trk[i])/spdlight[indx][0];
//           if (ADCw[i] > 0) tofw -= twalk_w[indx]/sqrt(ADCw[i]);

// 	  ZRES->Fill(tz_w[indx][0]
// 		     +tz_w[indx][1]*(z_trk[i]+L/2)
// 		     +tz_w[indx][2]*(z_trk[i]+L/2)*(z_trk[i]+L/2)
// 		     +tz_w[indx][3]*(z_trk[i]+L/2)*(z_trk[i]+L/2)*(z_trk[i]+L/2));

//           tofw -= ( tz_w[indx][0]
// 		    +tz_w[indx][1]*(z_trk[i]+L/2)
// 		    +tz_w[indx][2]*(z_trk[i]+L/2)*(z_trk[i]+L/2)
// 		    +tz_w[indx][3]*(z_trk[i]+L/2)*(z_trk[i]+L/2)*(z_trk[i]+L/2));
	  
          tofw -= t0g_clean; 

	  double etimee = parr[indx][0] + parr[indx][1]*z_trk[i];
	  double etimew = parr[indx+216][0] + parr[indx+216][1]*z_trk[i];
	  double w = 1./(etimee*etimee) + 1./(etimew*etimew);
	  double tof = tofe*1./(etimee*etimee) + tofw*1./(etimew*etimew);
	  tof /= w;
//           tof = (tofe + tofw)/2.0;

	  double Charge = ADCe[i];
	  double Chargw = ADCw[i];
	  double Chargen = 1.;// (tof_pulseQ[i][0]/tof_ttpMip[i][0] +
// 			    tof_pulseQ[i][1]/tof_ttpMip[i][1])/2.;
	  double Chargem = 1.;// (tof_pulseQ[i][0]/tof_ttpMip[i][0] +
// 			    tof_pulseQ[i][1]/tof_ttpMip[i][1])/2.*(4./tof_hbPathLen[i][0]);

	  EWZ->Fill(z_trk[i], tofe-tofw);
	  EWQ->Fill(Charge, tofe-tofw);
	  QNZ->Fill(Chargen, z_trk[i]);	
	  QMZ->Fill(Chargem, z_trk[i]);	
	  QZ->Fill(Charge, z_trk[i]);	
	  PZ->Fill(pathbar_trk[i], z_trk[i]);

	  if (inbar_trk[i]==0) {
	    if (fabs(z_trk[i])<L/2.) ZQ0->Fill(z_trk[i], Charge);
	    PQ0->Fill(pathbar_trk[i], Charge);
	    PTQ0->Fill(pt[i], Charge);
	  } else if (inbar_trk[i]==1) {
	    if (fabs(z_trk[i])<L/2.) ZQ1->Fill(z_trk[i], Charge);
	    PQ1->Fill(pathbar_trk[i], Charge);
	    PTQ1->Fill(pt[i], Charge);
	  }
	  ZQN->Fill(z_trk[i], Chargen);	

	  if (z_trk[i]<-90.) {
	    EVW->Fill(tofe, tofw);
	    QEVQW->Fill(Charge, Chargw);	
	  }

	  ZZZ->Fill(z_trk[i]);

	  PTDIS->Fill(pt[i]);
	  PDIS->Fill(pathbar_trk[i]);
	  if (pt[i]>0.4) {
	    if (fabs(tofe-tof_pi)<0.2) TQP1->Fill(tofe-tof_pi, Charge);
	    if (fabs(tofw-tof_pi)<0.2) TQP1W->Fill(tofw-tof_pi, Chargw);
	    if (fabs(tofe-tof_k)<0.2) TQK1->Fill(tofe-tof_k, Charge);
	    if (fabs(tof-tof_pi)<0.2) TQP1N->Fill(tof-tof_pi, Chargen);
	    if (fabs(tof-tof_k)<0.2) TQK1N->Fill(tof-tof_k, Chargen);
	    if (fabs(tof-tof_pi)<0.2) TQP1M->Fill(tof-tof_pi, Chargem);
// 	    if (fabs(tof-tof_pi)<=0.000001) 
// 	      cout << tof << " "<< tof_pi << " " << getTofSL(path_trk[i], p_trk[i], 0.13956995) << endl; 
	    DTP1->Fill(tof-tof_pi);  
	    DTK1->Fill(tof-tof_k); 
	    if (inbar_trk[i]==0) QDIS1->Fill(Charge); 
	    if (inbar_trk[i]==1) QDIS2->Fill(Charge); 
// 	    QDIS1->Fill((tof_pulseQ[i][0]+tof_pulseQ[i][1])/2.);
	    if (fabs(tof-tof_m)<0.2) TQM1N->Fill(tof-tof_m, Chargen);
	    if (fabs(tof-tof_m)<0.2) TQM1M->Fill(tof-tof_m, Chargem);
	    DTM1->Fill(tof-tof_m);  
	  }
	  if (pt[i]<0.6) {
	    if (fabs(tofw-tof_pi)<0.2) TQP2W->Fill(tofw-tof_pi, Chargw);
	    if (fabs(tof-tof_pi)<0.2) DTZ->Fill(tof-tof_pi, z_trk[i]);
	    if (fabs(tofe-tof_pi)<0.2) DTZE->Fill(tofe-tof_pi, z_trk[i]);	
	    if (fabs(tofw-tof_pi)<0.2) DTZW->Fill(tofw-tof_pi, z_trk[i]);
	    if (fabs(tofe-tof_pi)<0.2) TQP2->Fill(tofe-tof_pi, Charge);
	    if (fabs(tofe-tof_k)<0.2) TQK2->Fill(tofe-tof_k, Charge);
	    if (fabs(tof-tof_pi)<0.2) TQP2N->Fill(tof-tof_pi, Chargen);
	    if (fabs(tof-tof_pi)<0.2) TQP2M->Fill(tof-tof_pi, Chargem);
	    if (fabs(tof-tof_k)<0.2) TQK2N->Fill(tof-tof_k, Chargen);
	    DTP2->Fill(tof-tof_pi);  
	    DTK2->Fill(tof-tof_k); 
	    QNDIS2->Fill(Charge); 
// 	    QDIS2->Fill((tof_pulseQ[i][0]+tof_pulseQ[i][1])/2.);
	    if (fabs(tof-tof_m)<0.2) TQM2N->Fill(tof-tof_m, Chargen);
	    if (fabs(tof-tof_m)<0.2) TQM2M->Fill(tof-tof_m, Chargem);
	    DTM2->Fill(tof-tof_m);  
	  }
	  if (pt[i]<0.8) {
	    if (fabs(tofw-tof_pi)<0.2) TQP3W->Fill(tofw-tof_pi, Chargw);
	    if (fabs(tofe-tof_pi)<0.2) TQP3->Fill(tofe-tof_pi, Charge);
	    if (fabs(tofe-tof_k)<0.2) TQK3->Fill(tofe-tof_k, Charge);
	    if (fabs(tof-tof_pi)<0.2) TQP3N->Fill(tof-tof_pi, Chargen);
	    if (fabs(tof-tof_k)<0.2) TQK3N->Fill(tof-tof_k, Chargen);
	    if (fabs(tof-tof_pi)<0.2) TQP3M->Fill(tof-tof_pi, Chargem);
	    DTP3->Fill(tof-tof_pi);  
	    DTK3->Fill(tof-tof_k); 
	    QNDIS3->Fill(Charge); 
	    QDIS3->Fill((tof_pulseQ[i][0]+tof_pulseQ[i][1])/2.);
	    if (fabs(tof-tof_m)<0.2) TQM3N->Fill(tof-tof_m, Chargen);
	    if (fabs(tof-tof_m)<0.2) TQM3M->Fill(tof-tof_m, Chargem);
	    DTM3->Fill(tof-tof_m);  
	  }
	  if (pt[i]<1.5 && pt[i]>0.8) {
	    if (fabs(tofw-tof_pi)<0.2) TQP4W->Fill(tofw-tof_pi, Chargw);
	    if (fabs(tofe-tof_pi)<0.2) TQP4->Fill(tofe-tof_pi, Charge);
	    if (fabs(tofe-tof_k)<0.2) TQK4->Fill(tofe-tof_k, Charge);
	    if (fabs(tof-tof_pi)<0.2) TQP4N->Fill(tof-tof_pi, Chargen);
	    if (fabs(tof-tof_k)<0.2) TQK4N->Fill(tof-tof_k, Chargen);
	    if (fabs(tof-tof_pi)<0.2) TQP4M->Fill(tof-tof_pi, Chargem);
	    DTP4->Fill(tof-tof_pi);  
	    DTK4->Fill(tof-tof_k); 
	    QNDIS4->Fill(Charge); 
	    QDIS4->Fill((tof_pulseQ[i][0]+tof_pulseQ[i][1])/2.);
	    if (fabs(tof-tof_m)<0.2) TQM4N->Fill(tof-tof_m, Chargen);
	    if (fabs(tof-tof_m)<0.2) TQM4M->Fill(tof-tof_m, Chargem);
	    DTM4->Fill(tof-tof_m);  
	  }
	  if (pt[i]>1.5) {
	    if (fabs(tofw-tof_pi)<0.2) TQP5W->Fill(tofw-tof_pi, Chargw);
	    if (fabs(tofe-tof_pi)<0.2) TQP5->Fill(tofe-tof_pi, Charge);
	    if (fabs(tofe-tof_k)<0.2) TQK5->Fill(tofe-tof_k, Charge);
	    if (fabs(tof-tof_pi)<0.2) TQP5N->Fill(tof-tof_pi, Chargen);
	    if (fabs(tof-tof_k)<0.2) TQK5N->Fill(tof-tof_k, Chargen);
	    if (fabs(tof-tof_pi)<0.2) TQP5M->Fill(tof-tof_pi, Chargem);
	    DTP5->Fill(tof-tof_pi);  
	    DTK5->Fill(tof-tof_k); 
	    QNDIS5->Fill(Charge); 
	    QDIS5->Fill((tof_pulseQ[i][0]+tof_pulseQ[i][1])/2.);
	    if (fabs(tof-tof_m)<0.2) TQM5N->Fill(tof-tof_m, Chargen);
	    if (fabs(tof-tof_m)<0.2) TQM5M->Fill(tof-tof_m, Chargem);
	    DTM5->Fill(tof-tof_m);  
	  }

	  int ptbin[3] = {0,0,0};
	  if (pt[i]<0.6) ptbin[0] = 1;
	  if (pt[i]<0.8) ptbin[1] = 1;
	  if (pt[i]>0.8) ptbin[2] = 1;
	  int qbin = -1;
	  if (Charge<800) qbin = 0;
	  else if (Charge<1000) qbin = 1;
	  else if (Charge<1250) qbin = 2;
	  else if (Charge<1500) qbin = 3;
	  else if (Charge>1500) qbin = 4;
	  int pbin = -1;
	  if (pathbar_trk[i]<2) pbin = 0;
	  else if (pathbar_trk[i]<3) pbin = 1;
	  else if (pathbar_trk[i]<4) pbin = 2;
	  else if (pathbar_trk[i]<5) pbin = 3;
	  else if (pathbar_trk[i]>5) pbin = 4;
	  int zbin = int((z_trk[i]+150.)/60.);
	  for (int j = 0; j<3; j++)
	    if (ptbin[j]==1 && qbin>=0 && qbin<=4 && pbin>=0 && pbin<=4 && zbin>=0 && zbin<=4 && fabs(tofe-tof_pi)<0.2) {
	      TQPE[pbin][j]->Fill(tofe-tof_pi, Charge);
	      TPQE[qbin][j]->Fill(tofe-tof_pi, pathbar_trk[i]);
	      TQZE[zbin][j]->Fill(tofe-tof_pi, Charge);
	      TZQE[qbin][j]->Fill(tofe-tof_pi, z_trk[i]);
	      TZPE[pbin][j]->Fill(tofe-tof_pi, z_trk[i]);
	      TPZE[zbin][j]->Fill(tofe-tof_pi, pathbar_trk[i]);
	      if (Charge!=0.) TIQZE[zbin][j]->Fill(tofe-tof_pi, 1./Charge);
	      if (Charge!=0.) TIQZE[zbin][j]->Fill(tofe-tof_pi, 1./sqrt(Charge));
	    }

	  int mbin = -1;
	  if (pt[i]<0.6) mbin = 0;
	  else if (pt[i]<0.8) mbin = 1;
	  else if (pt[i]<1.0) mbin = 2;
	  else if (pt[i]<1.4) mbin = 3;
	  else if (pt[i]>1.4) mbin = 4;
	  if (mbin>=0 && fabs(tofe-tof_pi)<0.2)  TQME[mbin][0]->Fill(tofe-tof_pi, Charge);
	  if (mbin>=0 && fabs(tofe-tof_pi)<0.2)  TQME[mbin][1]->Fill(tofe-tof_pi, Charge);
	  if (mbin>=0 && fabs(tofe-tof_pi)<0.2)  TQME[mbin][2]->Fill(tofe-tof_pi, Charge);


        } // close loop over low pt tracks (for..nhit)
      } // else if..indx

// Begin: initialize some variables every event ----------------------------

      check_highest_pt = -999;
      pt_hi            =  2.0; // Set the high pt cut
      ibar_trkhi       = -999;
      t0g              = -999;
      nhit             =    0;

      for (int i=0; i<50; i++) {
        pt      [i] = -999;
        ibar_trk[i] = -999;
      }

      for (int i=0; i<NumBar; i++) {
        no_cut_cont[i] =    0;
        cont       [i] =    0;
        cont1      [i] =    0;
        adc_e      [i] = -999;
        adc_w      [i] = -999;
        curv       [i] = -999;
        length     [i] = -999;
        t_e        [i] = -999;
        t_w        [i] = -999;
        trk_p      [i] = -999;
        trk_pt     [i] = -999;
        z          [i] = -999;
        z0bc       [i] = -999;
        pathbar    [i] = -999;
        inbar      [i] = -999;
      }

      tof_ntt =    0; // Number of Tracks for T0 Calculation
      for (int tt = 0; tt<200; tt++) {
	tof_nttp[tt]    = 0;
	tof_nttpTrk[tt] = 0;
      }

      jj = 0; // Low pt track counter

// End: initialize some variables every event ------------------------------

    } // if event_counter > num_events

    //float d0_cut;
   // if (lightspeed_track_type == 2) d0_cut = 0.076;
  //  else d0_cut = 0.025;

    // Fill Variables for T0 Calculation 
    zv_nzv = lightspeed_num_vtx;
    for (int iz = 0; iz<zv_nzv; iz++) zv_zv[iz] = lightspeed_zvtx[iz];
    if (Chi2(lightspeed_te, lightspeed_tw, lightspeed_qe, lightspeed_qw, 
	     lightspeed_z, int(lightspeed_bar), lightspeed_pathinBar) <14.) {
      int indx = tof_ntt;
      for (int tt = 0; tt<tof_ntt; tt++)
	if (lightspeed_p==tof_pp[tt]) indx = tt;
      tof_pp[indx]                          = lightspeed_p;
      tof_ctz0bc[indx]                      = lightspeed_z0corr;
      tof_inBar[indx]                       = lightspeed_inBar;
      if (tof_nttp[indx]<4) {
	tof_ttpBar[indx][tof_nttp[indx]]      = int(lightspeed_bar);
	tof_ttpBar[indx][tof_nttp[indx]+1]    = int(lightspeed_bar);
	tof_hbInZ[indx][tof_nttp[indx]/2]     = lightspeed_z;
	tof_hb_s[indx][tof_nttp[indx]/2]      = lightspeed_tofLength;
	tof_hbPathLen[indx][tof_nttp[indx]/2] = lightspeed_pathinBar;
	tof_pulseQ[indx][tof_nttp[indx]]      = lightspeed_qe;
	tof_pulseQ[indx][tof_nttp[indx]+1]    = lightspeed_qw;
	tof_nttp[indx]    +=2;
      }
      tof_nttpTrk[indx] +=2;
      if (indx==tof_ntt) tof_ntt++;
      tp_tRaw[0][int(lightspeed_bar)] = lightspeed_tw;
      tp_tRaw[1][int(lightspeed_bar)] = lightspeed_te;
    }

    // Track selection
    if (  //  1==1
//         && (fabs(lightspeed_d0wrt) < 0.076) 
// 	   && (lightspeed_track_type  !=   2)         
//         && (lightspeed_pt          >  0.4)         
//         && (lightspeed_SiPhiHits   >    2)      
//         && (lightspeed_cotAxHits   >   10)      
//         && (lightspeed_cotStHits   >   10)      
//         && (lightspeed_sl8Hits     >    4)
//         && (lightspeed_inBar == 1) 
//        ) {   
	Chi2(lightspeed_te, lightspeed_tw, lightspeed_qe, lightspeed_qw, 
	     lightspeed_z, int(lightspeed_bar), lightspeed_pathinBar) <14.) {     
  
      ++ntrk_event_after_quality;
      ++cont[bar_number];

      adc_e [bar_number] = lightspeed_qe;
      adc_w [bar_number] = lightspeed_qw;
      curv  [bar_number] = lightspeed_curv;
      length[bar_number] = lightspeed_tofLength;
      t_e   [bar_number] = lightspeed_te;
      t_w   [bar_number] = lightspeed_tw;
      trk_p [bar_number] = lightspeed_p;
      trk_pt[bar_number] = lightspeed_pt;
      z     [bar_number] = lightspeed_z;
      z0bc  [bar_number] = lightspeed_z0corr;
      pathbar[bar_number]= lightspeed_pathinBar;
      inbar[bar_number]  = lightspeed_inBar;
      

    } else {

      ++cont1[bar_number];

    } // track selection

    } // trkPrescale

  } // loop over tracks

  printf("Loop validation: Tracks used %8d (run %d) \n",kentry,lightspeed_run);

  cout << "\n Number of events = " << event_counter << endl;

  TQP1->FitSlicesX();
  TH1D *TQP1_1 = (TH1D*)gDirectory->Get("TQP1_1");
  TQP1_1->SetLineColor(kGreen);

  TQK1->FitSlicesX();
  TQR1->FitSlicesX();
  TH1D* TQK1_1 = (TH1D*)gDirectory->Get("TQK1_1");
  TH1D* TQR1_1 = (TH1D*)gDirectory->Get("TQR1_1");  
  TQK1_1->SetLineColor(kRed);  
  TQR1_1->SetLineColor(kBlue);

  TQP1N->FitSlicesX();
  TH1D *TQP1N_1 = (TH1D*)gDirectory->Get("TQP1N_1");
  TQK1N->FitSlicesX();
  TH1D *TQK1N_1 = (TH1D*)gDirectory->Get("TQK1N_1");
  TQR1N->FitSlicesX();
  TH1D *TQR1N_1 = (TH1D*)gDirectory->Get("TQR1N_1");

  TQP2->FitSlicesX();
  TH1D *TQP2_1 = (TH1D*)gDirectory->Get("TQP2_1");
  TQP2_1->SetLineColor(kGreen);

  TQK2->FitSlicesX();
  TQR2->FitSlicesX();
  TH1D* TQK2_1 = (TH1D*)gDirectory->Get("TQK2_1");
  TH1D* TQR2_1 = (TH1D*)gDirectory->Get("TQR2_1");  
  TQK2_1->SetLineColor(kRed);  
  TQR2_1->SetLineColor(kBlue);

  TQP2N->FitSlicesX();
  TH1D *TQP2N_1 = (TH1D*)gDirectory->Get("TQP2N_1");
  TQK2N->FitSlicesX();
  TH1D *TQK2N_1 = (TH1D*)gDirectory->Get("TQK2N_1");
  TQR2N->FitSlicesX();
  TH1D *TQR2N_1 = (TH1D*)gDirectory->Get("TQR2N_1");

  TQP3->FitSlicesX();
  TH1D *TQP3_1 = (TH1D*)gDirectory->Get("TQP3_1");
  TQP3_1->SetLineColor(kGreen);

  TQK3->FitSlicesX();
  TQR3->FitSlicesX();
  TH1D* TQK3_1 = (TH1D*)gDirectory->Get("TQK3_1");
  TH1D* TQR3_1 = (TH1D*)gDirectory->Get("TQR3_1");  
  TQK3_1->SetLineColor(kRed);  
  TQR3_1->SetLineColor(kBlue);

  TQP3N->FitSlicesX();
  TH1D *TQP3N_1 = (TH1D*)gDirectory->Get("TQP3N_1");
  TQK3N->FitSlicesX();
  TH1D *TQK3N_1 = (TH1D*)gDirectory->Get("TQK3N_1");
  TQR3N->FitSlicesX();
  TH1D *TQR3N_1 = (TH1D*)gDirectory->Get("TQR3N_1");

  TQP4->FitSlicesX();
  TH1D *TQP4_1 = (TH1D*)gDirectory->Get("TQP4_1");
  TQP4_1->SetLineColor(kGreen);

  TQK4->FitSlicesX();
  TQR4->FitSlicesX();
  TH1D* TQK4_1 = (TH1D*)gDirectory->Get("TQK4_1");
  TH1D* TQR4_1 = (TH1D*)gDirectory->Get("TQR4_1");  
  TQK4_1->SetLineColor(kRed);  
  TQR4_1->SetLineColor(kBlue);

  TQP4N->FitSlicesX();
  TH1D *TQP4N_1 = (TH1D*)gDirectory->Get("TQP4N_1");
  TQK4N->FitSlicesX();
  TH1D *TQK4N_1 = (TH1D*)gDirectory->Get("TQK4N_1");
  TQR4N->FitSlicesX();
  TH1D *TQR4N_1 = (TH1D*)gDirectory->Get("TQR4N_1");

  TQP5->FitSlicesX();
  TH1D *TQP5_1 = (TH1D*)gDirectory->Get("TQP5_1");
  TQP5_1->SetLineColor(kGreen);

  TQK5->FitSlicesX();
  TQR5->FitSlicesX();
  TH1D* TQK5_1 = (TH1D*)gDirectory->Get("TQK5_1");
  TH1D* TQR5_1 = (TH1D*)gDirectory->Get("TQR5_1");  
  TQK5_1->SetLineColor(kRed);  
  TQR5_1->SetLineColor(kBlue);

  TQP5N->FitSlicesX();
  TH1D *TQP5N_1 = (TH1D*)gDirectory->Get("TQP5N_1");
  TQK5N->FitSlicesX();
  TH1D *TQK5N_1 = (TH1D*)gDirectory->Get("TQK5N_1");
  TQR5N->FitSlicesX();
  TH1D *TQR5N_1 = (TH1D*)gDirectory->Get("TQR5N_1");

  DTZ->FitSlicesX();
  TH1D* DTZ_1 = (TH1D*)gDirectory->Get("DTZ_1");
  DTZE->FitSlicesX();
  TH1D* DTZE_1 = (TH1D*)gDirectory->Get("DTZE_1");
  DTZW->FitSlicesX();
  TH1D* DTZW_1 = (TH1D*)gDirectory->Get("DTZW_1");
  EWZ->FitSlicesY();
  TH1D* EWZ_1 = (TH1D*)gDirectory->Get("EWZ_1");
  EWQ->FitSlicesY();
  TH1D* EWQ_1 = (TH1D*)gDirectory->Get("EWQ_1");

  QNZ->FitSlicesX();
  TH1D* QNZ_1 = (TH1D*)gDirectory->Get("QNZ_1");
  ZQN->FitSlicesX();
  TH1D* ZQN_1 = (TH1D*)gDirectory->Get("ZQN_1");
  QMZ->FitSlicesX();
  TH1D* QMZ_1 = (TH1D*)gDirectory->Get("QMZ_1");
  QZ->FitSlicesX();
  TH1D* QZ_1 = (TH1D*)gDirectory->Get("QZ_1");
  PZ->FitSlicesX();
  TH1D* PZ_1 = (TH1D*)gDirectory->Get("PZ_1");

  ZQ0->FitSlicesX();
  TH1D* ZQ0_1 = (TH1D*)gDirectory->Get("ZQ0_1");
  PQ0->FitSlicesX();
  TH1D* PQ0_1 = (TH1D*)gDirectory->Get("PQ0_1");
  PTQ0->FitSlicesX();
  TH1D* PTQ0_1 = (TH1D*)gDirectory->Get("PTQ0_1");

  ZQ1->FitSlicesX();
  TH1D* ZQ1_1 = (TH1D*)gDirectory->Get("ZQ1_1");
  PQ1->FitSlicesX();
  TH1D* PQ1_1 = (TH1D*)gDirectory->Get("PQ1_1");
  PTQ1->FitSlicesX();
  TH1D* PTQ1_1 = (TH1D*)gDirectory->Get("PTQ1_1");

  TQP1M->FitSlicesX();
  TH1D *TQP1M_1 = (TH1D*)gDirectory->Get("TQP1M_1");
  TQP2M->FitSlicesX();
  TH1D *TQP2M_1 = (TH1D*)gDirectory->Get("TQP2M_1");
  TQP3M->FitSlicesX();
  TH1D *TQP3M_1 = (TH1D*)gDirectory->Get("TQP3M_1");
  TQP4M->FitSlicesX();
  TH1D *TQP4M_1 = (TH1D*)gDirectory->Get("TQP4M_1");
  TQP5M->FitSlicesX();
  TH1D *TQP5M_1 = (TH1D*)gDirectory->Get("TQP5M_1");

  TQM1M->FitSlicesX();
  TH1D *TQM1M_1 = (TH1D*)gDirectory->Get("TQM1M_1");
  TQM2M->FitSlicesX();
  TH1D *TQM2M_1 = (TH1D*)gDirectory->Get("TQM2M_1");
  TQM3M->FitSlicesX();
  TH1D *TQM3M_1 = (TH1D*)gDirectory->Get("TQM3M_1");
  TQM4M->FitSlicesX();
  TH1D *TQM4M_1 = (TH1D*)gDirectory->Get("TQM4M_1");
  TQM5M->FitSlicesX();
  TH1D *TQM5M_1 = (TH1D*)gDirectory->Get("TQM5M_1");

  TQM1N->FitSlicesX();
  TH1D *TQM1N_1 = (TH1D*)gDirectory->Get("TQM1N_1");
  TQM2N->FitSlicesX();
  TH1D *TQM2N_1 = (TH1D*)gDirectory->Get("TQM2N_1");
  TQM3N->FitSlicesX();
  TH1D *TQM3N_1 = (TH1D*)gDirectory->Get("TQM3N_1");
  TQM4N->FitSlicesX();
  TH1D *TQM4N_1 = (TH1D*)gDirectory->Get("TQM4N_1");
  TQM5N->FitSlicesX();
  TH1D *TQM5N_1 = (TH1D*)gDirectory->Get("TQM5N_1");
 
  TQP1W->FitSlicesX();
  TH1D *TQP1W_1 = (TH1D*)gDirectory->Get("TQP1W_1");
  TQP2W->FitSlicesX();
  TH1D *TQP2W_1 = (TH1D*)gDirectory->Get("TQP2W_1");
  TQP3W->FitSlicesX();
  TH1D *TQP3W_1 = (TH1D*)gDirectory->Get("TQP3W_1");
  TQP4W->FitSlicesX();
  TH1D *TQP4W_1 = (TH1D*)gDirectory->Get("TQP4W_1");
  TQP5W->FitSlicesX();
  TH1D *TQP5W_1 = (TH1D*)gDirectory->Get("TQP5W_1");

  TQPE[0][0]->FitSlicesX();
  TH1D *TQPE0pt0_1 = (TH1D*)gDirectory->Get("TQPE0pt0_1");
  TQPE[0][1]->FitSlicesX();
  TH1D *TQPE0pt1_1 = (TH1D*)gDirectory->Get("TQPE0pt1_1");
  TQPE[0][2]->FitSlicesX();
  TH1D *TQPE0pt2_1 = (TH1D*)gDirectory->Get("TQPE0pt2_1");
  TQPE[1][0]->FitSlicesX();
  TH1D *TQPE1pt0_1 = (TH1D*)gDirectory->Get("TQPE1pt0_1");
  TQPE[1][1]->FitSlicesX();
  TH1D *TQPE1pt1_1 = (TH1D*)gDirectory->Get("TQPE1pt1_1");
  TQPE[1][2]->FitSlicesX();
  TH1D *TQPE1pt2_1 = (TH1D*)gDirectory->Get("TQPE1pt2_1");
  TQPE[2][0]->FitSlicesX();
  TH1D *TQPE2pt0_1 = (TH1D*)gDirectory->Get("TQPE2pt0_1");
  TQPE[2][1]->FitSlicesX();
  TH1D *TQPE2pt1_1 = (TH1D*)gDirectory->Get("TQPE2pt1_1");
  TQPE[2][2]->FitSlicesX();
  TH1D *TQPE2pt2_1 = (TH1D*)gDirectory->Get("TQPE2pt2_1");
  TQPE[3][0]->FitSlicesX();
  TH1D *TQPE3pt0_1 = (TH1D*)gDirectory->Get("TQPE3pt0_1");
  TQPE[3][1]->FitSlicesX();
  TH1D *TQPE3pt1_1 = (TH1D*)gDirectory->Get("TQPE3pt1_1");
  TQPE[3][2]->FitSlicesX();
  TH1D *TQPE3pt2_1 = (TH1D*)gDirectory->Get("TQPE3pt2_1");
  TQPE[4][0]->FitSlicesX();
  TH1D *TQPE4pt0_1 = (TH1D*)gDirectory->Get("TQPE4pt0_1");
  TQPE[4][1]->FitSlicesX();
  TH1D *TQPE4pt1_1 = (TH1D*)gDirectory->Get("TQPE4pt1_1");
  TQPE[4][2]->FitSlicesX();
  TH1D *TQPE4pt2_1 = (TH1D*)gDirectory->Get("TQPE4pt2_1");

  TPQE[0][0]->FitSlicesX();
  TH1D *TPQE0pt0_1 = (TH1D*)gDirectory->Get("TPQE0pt0_1");
  TPQE[0][1]->FitSlicesX();
  TH1D *TPQE0pt1_1 = (TH1D*)gDirectory->Get("TPQE0pt1_1");
  TPQE[0][2]->FitSlicesX();
  TH1D *TPQE0pt2_1 = (TH1D*)gDirectory->Get("TPQE0pt2_1");
  TPQE[1][0]->FitSlicesX();
  TH1D *TPQE1pt0_1 = (TH1D*)gDirectory->Get("TPQE1pt0_1");
  TPQE[1][1]->FitSlicesX();
  TH1D *TPQE1pt1_1 = (TH1D*)gDirectory->Get("TPQE1pt1_1");
  TPQE[1][2]->FitSlicesX();
  TH1D *TPQE1pt2_1 = (TH1D*)gDirectory->Get("TPQE1pt2_1");
  TPQE[2][0]->FitSlicesX();
  TH1D *TPQE2pt0_1 = (TH1D*)gDirectory->Get("TPQE2pt0_1");
  TPQE[2][1]->FitSlicesX();
  TH1D *TPQE2pt1_1 = (TH1D*)gDirectory->Get("TPQE2pt1_1");
  TPQE[2][2]->FitSlicesX();
  TH1D *TPQE2pt2_1 = (TH1D*)gDirectory->Get("TPQE2pt2_1");
  TPQE[3][0]->FitSlicesX();
  TH1D *TPQE3pt0_1 = (TH1D*)gDirectory->Get("TPQE3pt0_1");
  TPQE[3][1]->FitSlicesX();
  TH1D *TPQE3pt1_1 = (TH1D*)gDirectory->Get("TPQE3pt1_1");
  TPQE[3][2]->FitSlicesX();
  TH1D *TPQE3pt2_1 = (TH1D*)gDirectory->Get("TPQE3pt2_1");
  TPQE[4][0]->FitSlicesX();
  TH1D *TPQE4pt0_1 = (TH1D*)gDirectory->Get("TPQE4pt0_1");
  TPQE[4][1]->FitSlicesX();
  TH1D *TPQE4pt1_1 = (TH1D*)gDirectory->Get("TPQE4pt1_1");
  TPQE[4][2]->FitSlicesX();
  TH1D *TPQE4pt2_1 = (TH1D*)gDirectory->Get("TPQE4pt2_1");

  TQZE[0][0]->FitSlicesX();
  TH1D *TQZE0pt0_1 = (TH1D*)gDirectory->Get("TQZE0pt0_1");
  TQZE[0][1]->FitSlicesX();
  TH1D *TQZE0pt1_1 = (TH1D*)gDirectory->Get("TQZE0pt1_1");
  TQZE[0][2]->FitSlicesX();
  TH1D *TQZE0pt2_1 = (TH1D*)gDirectory->Get("TQZE0pt2_1");
  TQZE[1][0]->FitSlicesX();
  TH1D *TQZE1pt0_1 = (TH1D*)gDirectory->Get("TQZE1pt0_1");
  TQZE[1][1]->FitSlicesX();
  TH1D *TQZE1pt1_1 = (TH1D*)gDirectory->Get("TQZE1pt1_1");
  TQZE[1][2]->FitSlicesX();
  TH1D *TQZE1pt2_1 = (TH1D*)gDirectory->Get("TQZE1pt2_1");
  TQZE[2][0]->FitSlicesX();
  TH1D *TQZE2pt0_1 = (TH1D*)gDirectory->Get("TQZE2pt0_1");
  TQZE[2][1]->FitSlicesX();
  TH1D *TQZE2pt1_1 = (TH1D*)gDirectory->Get("TQZE2pt1_1");
  TQZE[2][2]->FitSlicesX();
  TH1D *TQZE2pt2_1 = (TH1D*)gDirectory->Get("TQZE2pt2_1");
  TQZE[3][0]->FitSlicesX();
  TH1D *TQZE3pt0_1 = (TH1D*)gDirectory->Get("TQZE3pt0_1");
  TQZE[3][1]->FitSlicesX();
  TH1D *TQZE3pt1_1 = (TH1D*)gDirectory->Get("TQZE3pt1_1");
  TQZE[3][2]->FitSlicesX();
  TH1D *TQZE3pt2_1 = (TH1D*)gDirectory->Get("TQZE3pt2_1");
  TQZE[4][0]->FitSlicesX();
  TH1D *TQZE4pt0_1 = (TH1D*)gDirectory->Get("TQZE4pt0_1");
  TQZE[4][1]->FitSlicesX();
  TH1D *TQZE4pt1_1 = (TH1D*)gDirectory->Get("TQZE4pt1_1");
  TQZE[4][2]->FitSlicesX();
  TH1D *TQZE4pt2_1 = (TH1D*)gDirectory->Get("TQZE4pt2_1");

  TZQE[0][0]->FitSlicesX();
  TH1D *TZQE0pt0_1 = (TH1D*)gDirectory->Get("TZQE0pt0_1");
  TZQE[0][1]->FitSlicesX();
  TH1D *TZQE0pt1_1 = (TH1D*)gDirectory->Get("TZQE0pt1_1");
  TZQE[0][2]->FitSlicesX();
  TH1D *TZQE0pt2_1 = (TH1D*)gDirectory->Get("TZQE0pt2_1");
  TZQE[1][0]->FitSlicesX();
  TH1D *TZQE1pt0_1 = (TH1D*)gDirectory->Get("TZQE1pt0_1");
  TZQE[1][1]->FitSlicesX();
  TH1D *TZQE1pt1_1 = (TH1D*)gDirectory->Get("TZQE1pt1_1");
  TZQE[1][2]->FitSlicesX();
  TH1D *TZQE1pt2_1 = (TH1D*)gDirectory->Get("TZQE1pt2_1");
  TZQE[2][0]->FitSlicesX();
  TH1D *TZQE2pt0_1 = (TH1D*)gDirectory->Get("TZQE2pt0_1");
  TZQE[2][1]->FitSlicesX();
  TH1D *TZQE2pt1_1 = (TH1D*)gDirectory->Get("TZQE2pt1_1");
  TZQE[2][2]->FitSlicesX();
  TH1D *TZQE2pt2_1 = (TH1D*)gDirectory->Get("TZQE2pt2_1");
  TZQE[3][0]->FitSlicesX();
  TH1D *TZQE3pt0_1 = (TH1D*)gDirectory->Get("TZQE3pt0_1");
  TZQE[3][1]->FitSlicesX();
  TH1D *TZQE3pt1_1 = (TH1D*)gDirectory->Get("TZQE3pt1_1");
  TZQE[3][2]->FitSlicesX();
  TH1D *TZQE3pt2_1 = (TH1D*)gDirectory->Get("TZQE3pt2_1");
  TZQE[4][0]->FitSlicesX();
  TH1D *TZQE4pt0_1 = (TH1D*)gDirectory->Get("TZQE4pt0_1");
  TZQE[4][1]->FitSlicesX();
  TH1D *TZQE4pt1_1 = (TH1D*)gDirectory->Get("TZQE4pt1_1");
  TZQE[4][2]->FitSlicesX();
  TH1D *TZQE4pt2_1 = (TH1D*)gDirectory->Get("TZQE4pt2_1");

  TQME[0][0]->FitSlicesX();
  TH1D *TQME0pt0_1 = (TH1D*)gDirectory->Get("TQME0pt0_1");
  TQME[0][1]->FitSlicesX();
  TH1D *TQME0pt1_1 = (TH1D*)gDirectory->Get("TQME0pt1_1");
  TQME[0][2]->FitSlicesX();
  TH1D *TQME0pt2_1 = (TH1D*)gDirectory->Get("TQME0pt2_1");
  TQME[1][0]->FitSlicesX();
  TH1D *TQME1pt0_1 = (TH1D*)gDirectory->Get("TQME1pt0_1");
  TQME[1][1]->FitSlicesX();
  TH1D *TQME1pt1_1 = (TH1D*)gDirectory->Get("TQME1pt1_1");
  TQME[1][2]->FitSlicesX();
  TH1D *TQME1pt2_1 = (TH1D*)gDirectory->Get("TQME1pt2_1");
  TQME[2][0]->FitSlicesX();
  TH1D *TQME2pt0_1 = (TH1D*)gDirectory->Get("TQME2pt0_1");
  TQME[2][1]->FitSlicesX();
  TH1D *TQME2pt1_1 = (TH1D*)gDirectory->Get("TQME2pt1_1");
  TQME[2][2]->FitSlicesX();
  TH1D *TQME2pt2_1 = (TH1D*)gDirectory->Get("TQME2pt2_1");
  TQME[3][0]->FitSlicesX();
  TH1D *TQME3pt0_1 = (TH1D*)gDirectory->Get("TQME3pt0_1");
  TQME[3][1]->FitSlicesX();
  TH1D *TQME3pt1_1 = (TH1D*)gDirectory->Get("TQME3pt1_1");
  TQME[3][2]->FitSlicesX();
  TH1D *TQME3pt2_1 = (TH1D*)gDirectory->Get("TQME3pt2_1");
  TQME[4][0]->FitSlicesX();
  TH1D *TQME4pt0_1 = (TH1D*)gDirectory->Get("TQME4pt0_1");
  TQME[4][1]->FitSlicesX();
  TH1D *TQME4pt1_1 = (TH1D*)gDirectory->Get("TQME4pt1_1");
  TQME[4][2]->FitSlicesX();
  TH1D *TQME4pt2_1 = (TH1D*)gDirectory->Get("TQME4pt2_1");

  TZPE[0][0]->FitSlicesX();
  TH1D *TZPE0pt0_1 = (TH1D*)gDirectory->Get("TZPE0pt0_1");
  TZPE[0][1]->FitSlicesX();
  TH1D *TZPE0pt1_1 = (TH1D*)gDirectory->Get("TZPE0pt1_1");
  TZPE[0][2]->FitSlicesX();
  TH1D *TZPE0pt2_1 = (TH1D*)gDirectory->Get("TZPE0pt2_1");
  TZPE[1][0]->FitSlicesX();
  TH1D *TZPE1pt0_1 = (TH1D*)gDirectory->Get("TZPE1pt0_1");
  TZPE[1][1]->FitSlicesX();
  TH1D *TZPE1pt1_1 = (TH1D*)gDirectory->Get("TZPE1pt1_1");
  TZPE[1][2]->FitSlicesX();
  TH1D *TZPE1pt2_1 = (TH1D*)gDirectory->Get("TZPE1pt2_1");
  TZPE[2][0]->FitSlicesX();
  TH1D *TZPE2pt0_1 = (TH1D*)gDirectory->Get("TZPE2pt0_1");
  TZPE[2][1]->FitSlicesX();
  TH1D *TZPE2pt1_1 = (TH1D*)gDirectory->Get("TZPE2pt1_1");
  TZPE[2][2]->FitSlicesX();
  TH1D *TZPE2pt2_1 = (TH1D*)gDirectory->Get("TZPE2pt2_1");
  TZPE[3][0]->FitSlicesX();
  TH1D *TZPE3pt0_1 = (TH1D*)gDirectory->Get("TZPE3pt0_1");
  TZPE[3][1]->FitSlicesX();
  TH1D *TZPE3pt1_1 = (TH1D*)gDirectory->Get("TZPE3pt1_1");
  TZPE[3][2]->FitSlicesX();
  TH1D *TZPE3pt2_1 = (TH1D*)gDirectory->Get("TZPE3pt2_1");
  TZPE[4][0]->FitSlicesX();
  TH1D *TZPE4pt0_1 = (TH1D*)gDirectory->Get("TZPE4pt0_1");
  TZPE[4][1]->FitSlicesX();
  TH1D *TZPE4pt1_1 = (TH1D*)gDirectory->Get("TZPE4pt1_1");
  TZPE[4][2]->FitSlicesX();
  TH1D *TZPE4pt2_1 = (TH1D*)gDirectory->Get("TZPE4pt2_1");

  TPZE[0][0]->FitSlicesX();
  TH1D *TPZE0pt0_1 = (TH1D*)gDirectory->Get("TPZE0pt0_1");
  TPZE[0][1]->FitSlicesX();
  TH1D *TPZE0pt1_1 = (TH1D*)gDirectory->Get("TPZE0pt1_1");
  TPZE[0][2]->FitSlicesX();
  TH1D *TPZE0pt2_1 = (TH1D*)gDirectory->Get("TPZE0pt2_1");
  TPZE[1][0]->FitSlicesX();
  TH1D *TPZE1pt0_1 = (TH1D*)gDirectory->Get("TPZE1pt0_1");
  TPZE[1][1]->FitSlicesX();
  TH1D *TPZE1pt1_1 = (TH1D*)gDirectory->Get("TPZE1pt1_1");
  TPZE[1][2]->FitSlicesX();
  TH1D *TPZE1pt2_1 = (TH1D*)gDirectory->Get("TPZE1pt2_1");
  TPZE[2][0]->FitSlicesX();
  TH1D *TPZE2pt0_1 = (TH1D*)gDirectory->Get("TPZE2pt0_1");
  TPZE[2][1]->FitSlicesX();
  TH1D *TPZE2pt1_1 = (TH1D*)gDirectory->Get("TPZE2pt1_1");
  TPZE[2][2]->FitSlicesX();
  TH1D *TPZE2pt2_1 = (TH1D*)gDirectory->Get("TPZE2pt2_1");
  TPZE[3][0]->FitSlicesX();
  TH1D *TPZE3pt0_1 = (TH1D*)gDirectory->Get("TPZE3pt0_1");
  TPZE[3][1]->FitSlicesX();
  TH1D *TPZE3pt1_1 = (TH1D*)gDirectory->Get("TPZE3pt1_1");
  TPZE[3][2]->FitSlicesX();
  TH1D *TPZE3pt2_1 = (TH1D*)gDirectory->Get("TPZE3pt2_1");
  TPZE[4][0]->FitSlicesX();
  TH1D *TPZE4pt0_1 = (TH1D*)gDirectory->Get("TPZE4pt0_1");
  TPZE[4][1]->FitSlicesX();
  TH1D *TPZE4pt1_1 = (TH1D*)gDirectory->Get("TPZE4pt1_1");
  TPZE[4][2]->FitSlicesX();
  TH1D *TPZE4pt2_1 = (TH1D*)gDirectory->Get("TPZE4pt2_1");

  TIQZE[0][0]->FitSlicesX();
  TH1D *TIQZE0pt0_1 = (TH1D*)gDirectory->Get("TIQZE0pt0_1");
  TIQZE[0][1]->FitSlicesX();
  TH1D *TIQZE0pt1_1 = (TH1D*)gDirectory->Get("TIQZE0pt1_1");
  TIQZE[0][2]->FitSlicesX();
  TH1D *TIQZE0pt2_1 = (TH1D*)gDirectory->Get("TIQZE0pt2_1");
  TIQZE[1][0]->FitSlicesX();
  TH1D *TIQZE1pt0_1 = (TH1D*)gDirectory->Get("TIQZE1pt0_1");
  TIQZE[1][1]->FitSlicesX();
  TH1D *TIQZE1pt1_1 = (TH1D*)gDirectory->Get("TIQZE1pt1_1");
  TIQZE[1][2]->FitSlicesX();
  TH1D *TIQZE1pt2_1 = (TH1D*)gDirectory->Get("TIQZE1pt2_1");
  TIQZE[2][0]->FitSlicesX();
  TH1D *TIQZE2pt0_1 = (TH1D*)gDirectory->Get("TIQZE2pt0_1");
  TIQZE[2][1]->FitSlicesX();
  TH1D *TIQZE2pt1_1 = (TH1D*)gDirectory->Get("TIQZE2pt1_1");
  TIQZE[2][2]->FitSlicesX();
  TH1D *TIQZE2pt2_1 = (TH1D*)gDirectory->Get("TIQZE2pt2_1");
  TIQZE[3][0]->FitSlicesX();
  TH1D *TIQZE3pt0_1 = (TH1D*)gDirectory->Get("TIQZE3pt0_1");
  TIQZE[3][1]->FitSlicesX();
  TH1D *TIQZE3pt1_1 = (TH1D*)gDirectory->Get("TIQZE3pt1_1");
  TIQZE[3][2]->FitSlicesX();
  TH1D *TIQZE3pt2_1 = (TH1D*)gDirectory->Get("TIQZE3pt2_1");
  TIQZE[4][0]->FitSlicesX();
  TH1D *TIQZE4pt0_1 = (TH1D*)gDirectory->Get("TIQZE4pt0_1");
  TIQZE[4][1]->FitSlicesX();
  TH1D *TIQZE4pt1_1 = (TH1D*)gDirectory->Get("TIQZE4pt1_1");
  TIQZE[4][2]->FitSlicesX();
  TH1D *TIQZE4pt2_1 = (TH1D*)gDirectory->Get("TIQZE4pt2_1");

  TIRQZE[0][0]->FitSlicesX();
  TH1D *TIRQZE0pt0_1 = (TH1D*)gDirectory->Get("TIRQZE0pt0_1");
  TIRQZE[0][1]->FitSlicesX();
  TH1D *TIRQZE0pt1_1 = (TH1D*)gDirectory->Get("TIRQZE0pt1_1");
  TIRQZE[0][2]->FitSlicesX();
  TH1D *TIRQZE0pt2_1 = (TH1D*)gDirectory->Get("TIRQZE0pt2_1");
  TIRQZE[1][0]->FitSlicesX();
  TH1D *TIRQZE1pt0_1 = (TH1D*)gDirectory->Get("TIRQZE1pt0_1");
  TIRQZE[1][1]->FitSlicesX();
  TH1D *TIRQZE1pt1_1 = (TH1D*)gDirectory->Get("TIRQZE1pt1_1");
  TIRQZE[1][2]->FitSlicesX();
  TH1D *TIRQZE1pt2_1 = (TH1D*)gDirectory->Get("TIRQZE1pt2_1");
  TIRQZE[2][0]->FitSlicesX();
  TH1D *TIRQZE2pt0_1 = (TH1D*)gDirectory->Get("TIRQZE2pt0_1");
  TIRQZE[2][1]->FitSlicesX();
  TH1D *TIRQZE2pt1_1 = (TH1D*)gDirectory->Get("TIRQZE2pt1_1");
  TIRQZE[2][2]->FitSlicesX();
  TH1D *TIRQZE2pt2_1 = (TH1D*)gDirectory->Get("TIRQZE2pt2_1");
  TIRQZE[3][0]->FitSlicesX();
  TH1D *TIRQZE3pt0_1 = (TH1D*)gDirectory->Get("TIRQZE3pt0_1");
  TIRQZE[3][1]->FitSlicesX();
  TH1D *TIRQZE3pt1_1 = (TH1D*)gDirectory->Get("TIRQZE3pt1_1");
  TIRQZE[3][2]->FitSlicesX();
  TH1D *TIRQZE3pt2_1 = (TH1D*)gDirectory->Get("TIRQZE3pt2_1");
  TIRQZE[4][0]->FitSlicesX();
  TH1D *TIRQZE4pt0_1 = (TH1D*)gDirectory->Get("TIRQZE4pt0_1");
  TIRQZE[4][1]->FitSlicesX();
  TH1D *TIRQZE4pt1_1 = (TH1D*)gDirectory->Get("TIRQZE4pt1_1");
  TIRQZE[4][2]->FitSlicesX();
  TH1D *TIRQZE4pt2_1 = (TH1D*)gDirectory->Get("TIRQZE4pt2_1");

  hfile->Write();
  hfile->Close();


} // speed_calib::LoopVal

// Method to compute the t0 of the event from the tables
// =====================================================
double speed_calib::CalculateT0(int nvtx, double *T0Error) {
  
  // Define useful constants
  // =======================
  int _guess = 0; 

  double tzero = -999., etzero = -999.;

  // If we are requiring unbiased t0, make sure there's more than one track
  // ======================================================================
  if (!tof_ntt>1) return tzero;
  
  // First fill the data arrays
  // ==========================
  bool GoodT0 = true;
  Nmatches_sl = 0;
  _skip_sl = -1;
  for (int t = 0; t<tof_ntt; t++ ) 
    if (fabs(tof_ctz0bc[t]-zv_zv[nvtx])<2.0 /*|| zv_zv[nvtx]==-999.*/) {
      Npulses[Nmatches_sl] = tof_nttp[t];
      if (tof_nttp[t]!=tof_nttpTrk[t]) GoodT0 = false;
      for (int j = 0; j<tof_nttp[t]; j++) {
	
	int h = j/2;
	double z = tof_hbInZ[t][h];
	_s[Nmatches_sl][j] = tof_hb_s[t][h];
	_p[Nmatches_sl][j] = tof_pp[t];  
	int side = 0;
	if (j==1 || j==3) side = 1;
	int channel = tof_ttpBar[t][j] + side*216;  
	double adc  = tof_pulseQ[t][j];
// 	double fq = sqrt(1./adc);   
	_t[Nmatches_sl][j] = TimeCorrT0(tp_tRaw[1-side][tof_ttpBar[t][j]], z, adc, channel);
// 	_t[Nmatches_sl][j] = CorrectedTime(tp_tRaw[1-side][tof_ttpBar[t][j]], z, adc, channel);
// 	_t[Nmatches_sl][j] = -1.*tp_tRaw[1-side][tof_ttpBar[t][j]];
// 	_t[Nmatches_sl][j] += parc[channel][0] + parc[channel][1]*z + parc[channel][2]*fq  + 
// 	  parc[channel][3]*z*z + parc[channel][4]*fq*fq + parc[channel][5]*z*fq + 
// 	  parc[channel][6]*z*z*z;
	_et[Nmatches_sl][j] = t0parr[channel][0] + t0parr[channel][1]*z;

      }
      if (tof_pp[t]==_skip_p) _skip_sl = Nmatches_sl;
      Nmatches_sl++;
    }
  if (Nmatches_sl<4) return -999.; 
  

  // Do a single track fit to obtain a reasonable starting value
  // ===========================================================
  double p = 0.;
  int best = 0;
  for(int t = 0; t<Nmatches_sl; t++)
    if (_p[t][0]>p) {
      p = _p[t][0];
      best = t;
    }
  
  double tzs [4], etzs[4];
  for (int j = 0; j<Npulses[best]; j++) {
    tzs[j]  = _t[best][j] - getTofSL(_s[best][j], _p[best][j], _m_sl[_guess]);
    etzs[j] = _et[best][j];
  }
 
  getMeanSL(tzero, etzero, tof_nttp[best], tzs, etzs);
  
  // Then do the pid likelihood fit, without skipping any matches
  // ============================================================
//   _skip_sl = -1;
  static TMinuit minuit(1);
  double arglis[10];
  int ierr = 0;
//   instance = this;

  minuit.SetFCN ( pid_likelihood_fcn_sl );
  minuit.SetPrintLevel(-1);

  arglis[0] = 1;
  minuit.mnexcm("SET SETNOWARNINGS",arglis,1,ierr);
  minuit.mncler();
  minuit.mnparm(0, "t0_mean", tzero, 0.01, -10., 10., ierr);
  arglis[0] = 1000;
  minuit.mnexcm ("MIGRAD",arglis,1, ierr);
  minuit.GetParameter(0, tzero, etzero );

  // Now loop over the matches, calculate in unbiased t0 for each one,
  // and fill the match with its t0, t0error, and tof
  // ================================================================= 

  if (!GoodT0) {
    tzero = -999.;
    etzero = -999.;
  }
  
  *T0Error = etzero;
  
  return tzero;

} // end double JPsi_calib::CalculateT0()