#define timew_calib_cxx
#include "timew_calib.h"

#include "TH2.h"
#include "TStyle.h"
#include "TCanvas.h"
#include "TCanvas.h"
#include "TGraph.h"
#include "TH2.h"
#include "TMath.h"
#include "TStyle.h"

#include <fstream>
#include <iomanip.h>      
#include <iostream.h>
#include <string.h>
#include <TMinuit.h>

const int    NumBar = 216;
const double L      = 279.5; 

char file[50], title[50];

double gx[5000], gy[5000];
const int maxevtnum = 40000;
int mynevts, twevts;


// 
const int nparm = 8;
const int mparm = nparm*2+1;
double q1_trk[maxevtnum], q2_trk[maxevtnum], z1_trk[maxevtnum], z2_trk[maxevtnum], dt_trk[maxevtnum];
double qe1_bar[maxevtnum][NumBar], qe2_bar[maxevtnum][NumBar], qw1_bar[maxevtnum][NumBar], qw2_bar[maxevtnum][NumBar];
double z1_bar[maxevtnum][NumBar], z2_bar[maxevtnum][NumBar], dte_bar[maxevtnum][NumBar], dtw_bar[maxevtnum][NumBar];
void fcntw(int    &npar, double *gin, double &f, double *par, int     iflag) {

  const double r2pi = 2.50662827479464932329;
  double offset   = par[0];
  double alfa1    = par[1];
  double alfa2    = par[2];
  double sigma    = par[3];
  double frac     = par[4];
  double offsetb  = par[5];
  double alfab    = par[6];
  double sigmab   = par[7];
  
  double logl = 0.;
  for ( int indx=0; indx<twevts; indx++ ) {
 
        double dt    = dt_trk[indx];
        double q1    = q1_trk[indx];
        double q2    = q2_trk[indx];
	double q1i( 1/sqrt(q1) );
	double q2i( 1/sqrt(q2) );

	double alfat1= alfa1;
	double alfat2= alfa2;
	double func  = offset  - alfat1*q1i  + alfat2*q2i; 
	double de    = dt - func;
	double db    = dt - (offsetb - alfab*q1i  + alfab*q2i);
	double arg   = frac*exp(-de*de/(2*sigma*sigma))/(r2pi*sigma) + 
		  (1.-frac)*exp(-db*db/(2*sigmab*sigmab))/(r2pi*sigmab);
	if(arg>0)logl += log(arg);
  }
  f = -2*logl;
}

void MyFittingFunction(int &npar, double *gin, double &f, double *par, int iflag) {
  const double r2pi = 2.50662827479464932329;
  double arg;
  double logl = 0;
  for (int i=0; i<mynevts; i++) {
     double de = gy[i] - (par[0] + par[1]*gx[i] + par[4]*gx[i]*gx[i]*gx[i]) ;
     double pe = par[2] * exp(-de*de/(2*par[3]*par[3]))/(r2pi*par[3]);
   //double db = par[4]; // background
   //double pb = (1-par[2])*exp(-db*db/(2*par[5]*par[5]))/(r2pi*par[5]);
   //arg = (pe+pb); // with background
     arg = pe; // without background
     if(arg > 0.0) logl += log(arg);
  }
  f = -2*logl;
}

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

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

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

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

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

    if (currentEvent == timewalk_event) newEvent = 0;
    else
    {
      if (currentEvent >= 0) newEvent = 1;
      currentEvent = timewalk_event;
    }

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

    if (kentry%100000 == 0) printf(" EventCounter:track %8d (run %d) \n",kentry,timewalk_run);

    if ( newEvent ) ++eventCounter;

    } // trkPrescale
  } // for..jentry

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

}

//--------------------------------------------------------------------------
// Time walk
//--------------------------------------------------------------------------

void timew_calib::LoopTW(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 LoopTW: Number of input tracks = " << nentries << " \n" << endl;

  if (flag==0) 
    sprintf(file,"%s/root/%s_timew-calib.root",directory,file_name);
  else 
    sprintf(file,"%s/root/%s_timew-calib_res.root",directory,file_name);
  TFile *hfile = new TFile(file,"recreate");

  gStyle->SetOptFit(111);

  // Variables
  int    twevtbar[NumBar]; for (int nb = 0; nb<NumBar; nb++) twevtbar[nb] = 0;
  double adc_comb_e, adc_comb_w, difft1, difft2, pl1, z1, z2;
  double resulte[NumBar][mparm], resultw[NumBar][mparm], resulsp[NumBar][2];
  double datae[mparm], dataw[mparm], datasp[9];

  double tqe1_bar[NumBar], tqe2_bar[NumBar], tqw1_bar[NumBar], tqw2_bar[NumBar];
  double  tz1_bar[NumBar],  tz2_bar[NumBar], tdte_bar[NumBar], tdtw_bar[NumBar];
  int cont[NumBar], cont1[NumBar];

  // Book histograms
  TH1F *barst   = new TH1F("barst",  "barst",  216,  0,216);
  TH1F *timec_e = new TH1F("timec_e","timec_e",216,  0,216);
  TH1F *timec_w = new TH1F("timec_w","timec_w",216,  0,216);

  TH1F *ADC      [NumBar];
  TH1F *reside   [NumBar];
  TH1F *residw   [NumBar];

  TH2F *residet  [NumBar];  
  TH2F *residez  [NumBar];  
  TH2F *residwt  [NumBar];
  TH2F *residwz  [NumBar];  

  TH2F *ADC_TDCE         [NumBar];
  TH2F *ADC_TDCE_bef_cuts[NumBar];

  TH2F *ADC_TDCW [NumBar];

  TH2F *ADC_TDCE_FIT[NumBar];
  TH2F *ADC_TDCW_FIT[NumBar];
                         
  TH2F *timeza   [NumBar];
  TH2F *timezb   [NumBar];

  TH2F *twqp;
  TH2F *twqn;

  TH2F *twgx;

  TH2F *resideq[NumBar];
  TH2F *residwq[NumBar];
  TH2F *residex[NumBar];
  TH2F *residwx[NumBar];
  TH2F *residgq;
  TH2F *residgz;

  TH2F *resideq_z[NumBar][5];
  TH2F *residwq_z[NumBar][5];
  TH2F *residex_z[NumBar][5];
  TH2F *residwx_z[NumBar][5];

  TH2F *residgq_z[5];
  TH2F *residgx_z[5];

  TH2F *residegq_z[5];
  TH2F *residegx_z[5];
  TH2F *residwgq_z[5];
  TH2F *residwgx_z[5];

  TH2F *residegx_z_d[5][6];
  TH2F *residwgx_z_d[5][6];

  TH2F *resideq_p[NumBar][5];
  TH2F *residwq_p[NumBar][5];
  TH2F *residex_p[NumBar][5];
  TH2F *residwx_p[NumBar][5];

  TH2F *residgq_p[5];
  TH2F *residgx_p[5];

  TH2F *residgq_q[5];
  TH2F *residgx_q[5];

  TH2F *residgq_qz[5];
  TH2F *residgx_qz[5];

  TH2F *residgz_q[5];

  TH2F *residgz_x[5];

  TH2F *residgq_q2[5];

  TH2F *QRE;
  TH2F *QRW;

  TH2F *QQDIS;
  TH1F *Q1DIS;
  TH1F *Q2DIS;

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

    sprintf(title,"ADC%d",i);
    ADC[i] = new TH1F(title,title,1000,0.,10000.); 
    sprintf(title,"bar_par_%d_%d",i,i+1);
    ADC_TDCE[i] = new TH2F(title,title,25,-0.05,0.05,100,-10,10);
    sprintf(title,"bar_par_bef_cuts_%d_%d",i,i+1);
    ADC_TDCE_bef_cuts[i] = new TH2F(title,title,25,-0.05,0.05,100,-10,10);
    sprintf(title,"bar_par2_%d_%d",i,i+1);
    ADC_TDCW[i] = new TH2F(title,title,25,-0.05,0.05,100,-10,10);
    sprintf(title,"resid_e_%d",i);
    reside[i] = new TH1F(title,title,100,-10,10);

    sprintf(title,"reside_vs_dt_%d",i);
    residet[i] = new TH2F(title,title,100,-10,10,100,-10,10);
    sprintf(title,"reside_vs_z_%d",i);
    residez[i] = new TH2F(title,title,100,-150,150,100,-10,10);
    sprintf(title,"resid_w_%d",i);
    residw[i] = new TH1F(title,title,150,-10,10);

    sprintf(title,"residw_vs_dt_%d",i);
    residwt[i] = new TH2F(title,title,100,-10,10,100,-10,10);
    sprintf(title,"residw_vs_z_%d",i);
    residwz[i] = new TH2F(title,title,100,-10,10,100,-10,10);
    sprintf(title,"dt_vs_za_%d",i);
    timeza[i] = new TH2F(title,title,300,-150,150,100,-10,10);
    sprintf(title,"dt_vs_zb_%d",i);
    timezb[i] = new TH2F(title,title,300,-150,150,100,-10,10); 

    sprintf(title,"bar_par_fit_%d_%d",i,i+1);
    ADC_TDCE_FIT[i] = new TH2F(title,title,25,-0.05,0.05,100,-10,10);
    sprintf(title,"bar_par2_fit_%d_%d",i,i+1);
    ADC_TDCW_FIT[i] = new TH2F(title,title,25,-0.05,0.05,100,-10,10);

    sprintf(title,"reside_vs_q_%d",i);
    resideq[i] = new TH2F(title,title,40,0,4000,100,-10,10);
    sprintf(title,"residw_vs_q_%d",i);
    residwq[i] = new TH2F(title,title,40,0,4000,100,-10,10);

    sprintf(title,"reside_vs_x_%d",i);
    residex[i] = new TH2F(title,title,25,-0.05,0.05,100,-10,10);
    sprintf(title,"residw_vs_x_%d",i);
    residwx[i] = new TH2F(title,title,25,-0.05,0.05,100,-10,10);

    for (int b = 0; b<5; b++) {
      sprintf(title,"reside_vs_q_z_%d_%d",i,b);
      resideq_z[i][b] = new TH2F(title,title,40,0,4000,100,-10,10);
      sprintf(title,"residw_vs_q_z_%d_%d",i,b);
      residwq_z[i][b] = new TH2F(title,title,40,0,4000,100,-10,10);
    }

    for (int b = 0; b<5; b++) {
      sprintf(title,"reside_vs_x_z_%d_%d",i,b);
      residex_z[i][b] = new TH2F(title,title,25,-0.05,0.05,100,-10,10);
      sprintf(title,"residw_vs_x_z_%d_%d",i,b);
      residwx_z[i][b] = new TH2F(title,title,25,-0.05,0.05,100,-10,10);
    }

    for (int b = 0; b<5; b++) {
      sprintf(title,"reside_vs_q_p_%d_%d",i,b);
      resideq_p[i][b] = new TH2F(title,title,40,0,4000,100,-10,10);
      sprintf(title,"residw_vs_q_p_%d_%d",i,b);
      residwq_p[i][b] = new TH2F(title,title,40,0,4000,100,-10,10);
    }

    for (int b = 0; b<5; b++) {
      sprintf(title,"reside_vs_x_p_%d_%d",i,b);
      residex_p[i][b] = new TH2F(title,title,25,-0.05,0.05,100,-10,10);
      sprintf(title,"residw_vs_x_p_%d_%d",i,b);
      residwx_p[i][b] = new TH2F(title,title,25,-0.05,0.05,100,-10,10);
    }

  }

  sprintf(title,"tw_pos");
  twqp = new TH2F(title,title,25,-0.05,0.05,100,-10,10);

  sprintf(title,"tw_neg");
  twqn = new TH2F(title,title,25,-0.05,0.05,100,-10,10);

  sprintf(title,"tw_vs_x");
  twgx = new TH2F(title,title,25,-0.05,0.05,100,-10,10);

  sprintf(title,"resid_vs_q");
  residgq = new TH2F(title,title,40,0,4000,100,-10,10);

  sprintf(title,"resid_vs_z");
  residgz = new TH2F(title,title,30,-150,150,100,-10,10);

  for (int b = 0; b<5; b++) {

    sprintf(title,"resid_vs_q_z_%d",b);
    residgq_z[b] = new TH2F(title,title,40,0,4000,100,-10,10);
    
    sprintf(title,"resid_vs_x_z_%d",b);
    residgx_z[b] = new TH2F(title,title,25,-0.05,0.05,100,-10,10);

    sprintf(title,"reside_vs_q_z_%d",b);
    residegq_z[b] = new TH2F(title,title,40,0,4000,100,-10,10);
    
    sprintf(title,"reside_vs_x_z_%d",b);
    residegx_z[b] = new TH2F(title,title,25,-0.05,0.05,100,-10,10);
    for (int y = 0; y<6; y++) {
      sprintf(title,"reside_vs_x_z_%d_d_%d",b,y);
      residegx_z_d[b][y] = new TH2F(title,title,25,-0.05,0.05,100,-10,10);
    }

    sprintf(title,"residw_vs_q_z_%d",b);
    residwgq_z[b] = new TH2F(title,title,40,0,4000,100,-10,10);
    for (int y = 0; y<6; y++) {
      sprintf(title,"residw_vs_x_z_%d_d_%d",b,y);
      residwgx_z_d[b][y] = new TH2F(title,title,25,-0.05,0.05,100,-10,10);
    }
    
    sprintf(title,"residw_vs_x_z_%d",b);
    residwgx_z[b] = new TH2F(title,title,25,-0.05,0.05,100,-10,10);
    
    sprintf(title,"resid_vs_q_p_%d",b);
    residgq_p[b] = new TH2F(title,title,40,0,4000,100,-10,10);
    
    sprintf(title,"resid_vs_x_p_%d",b);
    residgx_p[b] = new TH2F(title,title,25,-0.05,0.05,100,-10,10);
    
    sprintf(title,"resid_vs_q_q_%d",b);
    residgq_q[b] = new TH2F(title,title,40,0,4000,100,-10,10);
    
    sprintf(title,"resid_vs_x_q_%d",b);
    residgx_q[b] = new TH2F(title,title,25,-0.05,0.05,100,-10,10);
    
    sprintf(title,"resid_vs_q_qz_%d",b);
    residgq_qz[b] = new TH2F(title,title,40,0,4000,100,-10,10);
    
    sprintf(title,"resid_vs_x_qz_%d",b);
    residgx_qz[b] = new TH2F(title,title,25,-0.05,0.05,100,-10,10);
    
    sprintf(title,"resid_vs_z_q_%d",b);
    residgz_q[b] = new TH2F(title,title,30,-150,150,100,-10,10);
    
    sprintf(title,"resid_vs_z_x_%d",b);
    residgz_x[b] = new TH2F(title,title,30,-150,150,100,-10,10);
    
    sprintf(title,"resid_vs_q_q2_%d",b);
    residgq_q2[b] = new TH2F(title,title,40,0,4000,100,-10,10);
  }

  sprintf(title,"QQDIS");
  QQDIS = new TH2F(title,title,40,0.,4000.,40,0.,4000.);
  sprintf(title,"Q1DIS");
  Q1DIS = new TH1F(title,title,40,0.,4000.);
  sprintf(title,"Q2DIS");
  Q2DIS = new TH1F(title,title,40,0.,4000.);
  sprintf(title,"QRE");
  QRE = new TH2F(title,title,30,-150.,150., 50,0.,5);
  sprintf(title,"QRW");
  QRW = new TH2F(title,title,30,-150.,150., 50,0.,5);
  
  ifstream twe_result, tww_result, t0e_result, t0w_result, speed_result; 
  ifstream twe_slopes, tww_slopes;
  ofstream twe_calib, tww_calib, ofe_calib, ofw_calib;
  ofstream twe_alphas_diff, tww_alphas_diff;
  ofstream twe_fit, tww_fit;
  
  //-------------------------------------------------------------------------
  // Open the speed of light and time walk calibration files
  // flag <= 0      --> No input file
  // 0 < flag < 100 --> Speed of light + TW input files
  // flag > 100     --> Speed of light after TW + TW input files
  //-------------------------------------------------------------------------

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

    sprintf(file,"%s/files/%s_timew-alphas_diff_e.txt",directory,file_name);
    twe_slopes.open(file,ios::in);
    sprintf(file,"%s/files/%s_timew-alphas_diff_w.txt",directory,file_name); cout << file << endl;
    tww_slopes.open(file,ios::in);
  }

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

  for(int j=0; j<NumBar; j++) {
      
    cont[j]     = 0;
    cont1[j]    = 0;
    twevtbar[j] = 0;

    // new mine...
    for(int i=0; i<mparm; i++) {
      datae [i] = 0;
      dataw [i] = 0;
      twe_slopes >> datae[i];
      tww_slopes >> dataw[i];
    }

    for(int i=0; i<mparm-1; i++) {
      resulte[j][i] = datae[i+1];
      resultw[j][i] = dataw[i+1];
    }

    // new mine...
    for(int i=0; i<9; i++) datasp[i] = 0;
    for(int i=0; i<9; i++) speed_result >> datasp[i];

    resulsp[j][0] = datasp[1];
    resulsp[j][1] = datasp[3];
  } // for..int j<NumBar

  twe_result.close();
  tww_result.close();
  twe_slopes.close();
  tww_slopes.close();
  speed_result.close();

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

    if (currentEvent == timewalk_event) newEvent = 0;
    else
    {
      if (currentEvent >= 0) newEvent = 1;
      currentEvent = timewalk_event;
      ++rawEvents;
    }

    if ( (rawEvents%trkPrescale)<10 )
    {
    // Prescaled event count: kentry
    ++kentry;
    
    //if (kentry%100000 == 0) printf(" LoopTW:track %8d (run %d) \n",kentry,timewalk_run);

    if ( newEvent )		// Reset for next event
    {
      for (int j=0; j<NumBar; ++j)
      {
	if ( (cont[j]==1) && (cont1[j]==1) && (twevtbar[j]<maxevtnum) ) 
	{
	  z1_bar[twevtbar[j]][j]  = tz1_bar[j];
	  z2_bar[twevtbar[j]][j]  = tz2_bar[j];
	  qe1_bar[twevtbar[j]][j] = tqe1_bar[j];
	  qe2_bar[twevtbar[j]][j] = tqe2_bar[j];
	  qw1_bar[twevtbar[j]][j] = tqw1_bar[j];
	  qw2_bar[twevtbar[j]][j] = tqw2_bar[j];
	  dte_bar[twevtbar[j]][j] = tdte_bar[j];
	  dtw_bar[twevtbar[j]][j] = tdtw_bar[j];
	  twevtbar[j]++;
	}
	cont[j]=cont1[j]=0;
      }
    }

    int bar1 = int(timewalk_bar1);
    int bar2 = int(timewalk_bar2);

    if ( (bar1<0) || (bar1>216)) continue;
    if ( (bar2<0) || (bar2>216)) continue;

    // Track selection
    if (   (timewalk_track_type  !=    2)         
        && (timewalk_pt           >  0.4) // different from lightspeed	 
        && (timewalk_SiPhiHits    >    2)         
        && (timewalk_cotAxHits    >   10)  
        && (timewalk_cotStHits    >   10)  
        && (timewalk_sl8Hits      >    4)
        && (timewalk_inBar        ==   1) // Is this 1 track per event ??
       )
      {

      ++cont1[bar1];

      QQDIS->Fill(timewalk_qe1); QQDIS->Fill(timewalk_qw1, timewalk_qw2);
      Q1DIS->Fill(timewalk_qe1); Q1DIS->Fill(timewalk_qw1);
      Q2DIS->Fill(timewalk_qe2); Q2DIS->Fill(timewalk_qw2);

      QRE->Fill( (timewalk_z1+timewalk_z2)/2., timewalk_qe1/timewalk_qe2);
      QRW->Fill( (timewalk_z1+timewalk_z2)/2., timewalk_qw1/timewalk_qw2);

      ADC[bar1]->Fill(timewalk_qe1);
      barst->Fill(timewalk_bar1);

      if (fabs((timewalk_te1-timewalk_tw1)-(2*timewalk_z1/resulsp[bar1][0]+resulsp[bar1][1])) > cut) continue;
      if (fabs((timewalk_te2-timewalk_tw2)-(2*timewalk_z2/resulsp[bar2][0]+resulsp[bar2][1])) > cut) continue;

      // Single track per bar per event

      if (cont1[bar1] != 1) continue;

      if ( !finite(timewalk_qe1) ) continue;
      if ( !finite(timewalk_qe2) ) continue;
      if (timewalk_qe1<=1e-64 || timewalk_qe2<1e-64) continue;

      if ( !finite(timewalk_qw1) ) continue;
      if ( !finite(timewalk_qw2) ) continue;
      if (timewalk_qw1<=1e-64 || timewalk_qw2<1e-64) continue;

      if ( !finite(timewalk_te1) ) continue;
      if ( !finite(timewalk_te2) ) continue;

      if ( !finite(timewalk_tw1) ) continue;
      if ( !finite(timewalk_tw2) ) continue;

      if ( !finite(timewalk_z1)  ) continue;
      if ( !finite(timewalk_z2)  ) continue;

      ++cont[bar1];

      adc_comb_e = (sqrt(timewalk_qe1)-sqrt(timewalk_qe2))/
                   (sqrt(timewalk_qe1)*sqrt(timewalk_qe2));

      adc_comb_w = (sqrt(timewalk_qw1)-sqrt(timewalk_qw2))/
                   (sqrt(timewalk_qw1)*sqrt(timewalk_qw2));

      timec_e->Fill(resulte[bar1][0]/sqrt(timewalk_qe1));
      timec_w->Fill(resultw[bar1][0]/sqrt(timewalk_qw1));
      timec_e->Fill(resulte[bar1][0]/sqrt(timewalk_qe2));
      timec_w->Fill(resultw[bar1][0]/sqrt(timewalk_qw2));

      if (timewalk_curv > 0) {
        z1  = timewalk_z1;
        z2  = timewalk_z2;
        pl1 = timewalk_path1;
      } else {
        z1  = timewalk_z2;
        z2  = timewalk_z1;
        pl1 = timewalk_path2; 
      }

      difft1 = timewalk_te1 - timewalk_te2;
      difft2 = timewalk_tw1 - timewalk_tw2;


      if (resulsp[bar1][0] > 0 && resulsp[bar2][0] > 0) {

        difft1 -= (L/2)*(1.0/resulsp[bar2][0] - 1.0/resulsp[bar1][0])
                  - (z1/resulsp[bar1][0] - z2/resulsp[bar2][0])
		  - pl1/timewalk_beta;

        difft2 -= (L/2)*(1.0/resulsp[bar2][0] - 1.0/resulsp[bar1][0]) 
                  + (z1/resulsp[bar1][0] - z2/resulsp[bar2][0])
                  - pl1/timewalk_beta;
      }

      ADC_TDCE_bef_cuts[bar1]->Fill(adc_comb_e,difft1);

      ADC_TDCE[bar1]->Fill(adc_comb_e,difft1);
      ADC_TDCW[bar1]->Fill(adc_comb_w,difft2); 

      if (timewalk_curv > 0) {
	twqp->Fill(adc_comb_e,difft1);
	twqp->Fill(adc_comb_w,difft2);
      } else { 
	twqn->Fill(adc_comb_e,difft1);
	twqn->Fill(adc_comb_w,difft2);
      }
 
      twgx->Fill(adc_comb_e,difft1);
      twgx->Fill(adc_comb_w,difft2);

      // Data used in fits
      tz1_bar[bar1]  = timewalk_z1;
      tz2_bar[bar1]  = timewalk_z2;
      tqe1_bar[bar1] = timewalk_qe1;
      tqe2_bar[bar1] = timewalk_qe2;
      tqw1_bar[bar1] = timewalk_qw1;
      tqw2_bar[bar1] = timewalk_qw2;
      tdte_bar[bar1] = difft1;
      tdtw_bar[bar1] = difft2;

   
      double e0 = 0., e1 = 0., w0 = 0., w1 = 0.;
      if (flag==1) {
	e0  = resulte[bar1][4]/sqrt(timewalk_qe2) - resulte[bar1][2]/sqrt(timewalk_qe1);
	e1 = resulte[bar1][0];
	w0  = resultw[bar1][4]/sqrt(timewalk_qw2) - resultw[bar1][2]/sqrt(timewalk_qw1);
	w1 = resultw[bar1][0];
      }	

      timeza [bar1]->Fill(timewalk_z1,timewalk_te1 - timewalk_tw1);
      timezb [bar1]->Fill(timewalk_z2,timewalk_te2 - timewalk_tw2); 
 
      reside [bar1]->Fill(difft1-(e0+e1));
      residw [bar1]->Fill(difft2-(w0+w1));  

      residet[bar1]->Fill(difft1,difft1-(e0+e1));
      residwt[bar1]->Fill(difft2,difft2-(w0+w1));

      residez[bar1]->Fill(timewalk_z1,difft1-(e0+e1));
      residwz[bar1]->Fill(timewalk_z1,difft2-(w0+w1));

      resideq[bar1]->Fill(timewalk_qe1, difft1-(e0+e1));
      residex[bar1]->Fill(adc_comb_e, difft1-(e0+e1));

      residwq[bar1]->Fill(timewalk_qw1, difft2-(w0+w1));
      residwx[bar1]->Fill(adc_comb_w, difft2-(w0+w1));
      
      residgq->Fill(timewalk_qe1, difft1-(e0+e1));
      
      residgz->Fill(timewalk_z1, difft1-(e0+e1));

   
      int zbin = int ( (timewalk_z1+L/2)/60. );
      double dang = asin((z2-z1)/pl1); 
      int dbin = -1; 
      if (dang<-3.1428/6.) dbin = 0;
      else if (dang<-3.1428/6.) dbin = 1;
      else if (dang<0.) dbin = 2;  
      else if (dang<3.1428/6.) dbin = 3;
      else if (dang<3.1428/12.) dbin = 4;
      else dbin = 5;  
   
      if (zbin>=0 && zbin<=4) {
	
	resideq_z[bar1][zbin]->Fill(timewalk_qe1, difft1-(e0+e1));
	residex_z[bar1][zbin]->Fill(adc_comb_e, difft1-(e0+e1));
	
	residwq_z[bar1][zbin]->Fill(timewalk_qw1, difft2-(w0+w1));
	residwx_z[bar1][zbin]->Fill(adc_comb_w, difft2-(w0+w1));
	
	residgq_z[zbin]->Fill(timewalk_qe1, difft1-(e0+e1));
	residgq_z[zbin]->Fill(timewalk_qw1, difft2-(w0+w1));
	
	residgx_z[zbin]->Fill(adc_comb_e, difft1-(e0+e1));
	residgx_z[zbin]->Fill(adc_comb_w, difft2-(w0+w1));
	
	residegq_z[zbin]->Fill(timewalk_qe1, difft1-(e0+e1));
	residwgq_z[zbin]->Fill(timewalk_qw1, difft2-(w0+w1));
	
	residegx_z[zbin]->Fill(adc_comb_e, difft1-(e0+e1));
	residwgx_z[zbin]->Fill(adc_comb_w, difft2-(w0+w1));
	
	if (dbin>=0 && dbin<=5) {
	  residegx_z_d[zbin][dbin]->Fill(adc_comb_e, difft1-(e0+e1));
	  residwgx_z_d[zbin][dbin]->Fill(adc_comb_w, difft2-(w0+w1));
	}

      }

      int pbin = int ( (pl1)/1.6 );
      if (pbin>=0 && pbin<=4) {
	
	resideq_p[bar1][pbin]->Fill(timewalk_qe1, difft1-(e0+e1));
	residex_p[bar1][pbin]->Fill(adc_comb_e, difft1-(e0+e1));
	
	residwq_p[bar1][pbin]->Fill(timewalk_qw1, difft2-(w0+w1));
	residwx_p[bar1][pbin]->Fill(adc_comb_w, difft2-(w0+w1));
	
	residgq_p[pbin]->Fill(timewalk_qe1, difft1-(e0+e1));
	residgq_p[pbin]->Fill(timewalk_qw1, difft2-(w0+w1));
	
	residgx_p[pbin]->Fill(adc_comb_e, difft1-(e0+e1));
	residgx_p[pbin]->Fill(adc_comb_w, difft2-(w0+w1));
	
      }
 
      for (int qbin = 0; qbin<5; qbin++) {
	if (timewalk_qe2>3000+qbin*250) {
	  residgq_q[qbin]->Fill(timewalk_qe1, difft1-(e0+e1));
	  residgx_q[qbin]->Fill(adc_comb_e, difft1-(e0+e1));
	  if (qbin==2) {
	    int binz = -1; 
	    if (fabs(timewalk_z1+100.)<10.) binz = 0;
	    if (fabs(timewalk_z1+50.)<10.) binz = 1;
	    if (fabs(timewalk_z1)<10.) binz = 2;
	    if (fabs(timewalk_z1-50.)<10.) binz = 3;
	    if (fabs(timewalk_z1-100.)<10.) binz = 4;
	    if (binz!=-1) {
	      residgq_qz[binz]->Fill(timewalk_qe1, difft1-(e0+e1));
	      residgx_qz[binz]->Fill(adc_comb_e, difft1-(e0+e1));
	    }
	  }
	  if (qbin>=2) {
	    int qbin2 = int (timewalk_qe1/800.);
	    if (qbin2>=0 && qbin2<=4) residgz_q[qbin2]->Fill(timewalk_z1, difft1-(e0+e1));
	    int xbin2 = int ((adc_comb_e+0.05)/0.02);
	    if (xbin2>=0 && xbin2<=4) residgz_x[xbin2]->Fill(timewalk_z1, difft1-(e0+e1));
	  }
	}
      }
 
      for (int qbin = 0; qbin<5; qbin++) {
	if (timewalk_qw2>3000+qbin*250) {
	  residgq_q[qbin]->Fill(timewalk_qw1, difft2-(w0+w1));
	  residgx_q[qbin]->Fill(adc_comb_w, difft2-(w0+w1));
	  if (fabs(timewalk_z1)<20.) {
	    residgq_qz[qbin]->Fill(timewalk_qw1, difft2-(w0+w1));
	    residgx_qz[qbin]->Fill(adc_comb_w, difft2-(w0+w1));
	  }
	}
      }
 
      for (int qbin = 0; qbin<5; qbin++) {
	if (timewalk_qe1>3000+qbin*250) {
	  residgq_q2[qbin]->Fill(timewalk_qe2, difft1-(e0+e1));
	}
      }
 
      for (int qbin = 0; qbin<5; qbin++) {
	if (timewalk_qw1>3000+qbin*250) {
	  residgq_q2[qbin]->Fill(timewalk_qw2, difft2-(w0+w1));
	}
      }
 
    } // End track selection


    } // trkPrescale
  } // loop over tracks

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


 if (flag==0) {

   cout << " Performing likelihood fits for time walk" << endl;

   sprintf(file,"%s/files/%s_timew-alphas_diff_e.txt",directory,file_name);
   twe_alphas_diff.open(file,ios::out); 
   
   sprintf(file,"%s/files/%s_timew-alphas_diff_w.txt",directory,file_name);
   tww_alphas_diff.open(file,ios::out);

   double parme[NumBar][nparm], eparme[NumBar][nparm];
   double parmw[NumBar][nparm], eparmw[NumBar][nparm];

   TMinuit *gMinuit = new TMinuit(nparm);
   int    ierr     =    0;
   double arglis[10];
   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);
   for (int j=0 ;j<NumBar; j++) {
     
     int ntwevts = twevtbar[j];

     if (maxevtnum<=ntwevts) cout << "Attention: MAX number of tracks in bar " << j << "!" << endl;

     double p1, p2, p3, p4, p5, p6, eoff;
     double eslope = ajuste.fitlin3(ADC_TDCE[j], p1, eoff, p2, p3, p4, p5, p6);

     const double slpmin(   0 );
     const double slpmax( 150 );
     const double offmin(  -5 );
     const double offmax(   5 );
     
     if (eslope < slpmin || eslope > slpmax ) eslope = 0.5*(slpmin+slpmax);
     if (  eoff < offmin ||   eoff > offmax )   eoff = 0.5*(offmin+offmax);

     twevts = 0;
     for (int k=0; k<ntwevts; k++) {
       if (fabs(dte_bar[k][j]-eoff-eslope*(1./sqrt(qe2_bar[k][j])-1./sqrt(qe1_bar[k][j])))>2.) continue; 
       z1_trk[twevts] = z1_bar[k][j];
       z2_trk[twevts] = z2_bar[k][j];
       q1_trk[twevts] = qe1_bar[k][j];
       q2_trk[twevts] = qe2_bar[k][j];
       dt_trk[twevts] = dte_bar[k][j];
       ADC_TDCE_FIT[j]->Fill((sqrt(q1_trk[twevts])-sqrt(q2_trk[twevts]))/sqrt(q1_trk[twevts]*q2_trk[twevts]),dt_trk[twevts]);
       twevts++;
     }

     // Fit is unstable presumably because signal and background
     // shape functional forms are very similar. Start them far apart. 
      
     gMinuit->mncler();
     gMinuit->SetFCN(fcntw);
     gMinuit->mnparm(0,"A",eoff,0.1,offmin,offmax,ierr);
     gMinuit->mnparm(1,"B",eslope,1,slpmin,slpmax,ierr);
     gMinuit->mnparm(2,"C",eslope,1,slpmin,slpmax,ierr);
     gMinuit->mnparm(3,"F",0.3,0.1,0.1,4,ierr);
     gMinuit->mnparm(4,"K",0.9,0.01,0.6,1,ierr);
     gMinuit->mnparm(5,"L",0,0.1,offmin,offmax,ierr);
     gMinuit->mnparm(6,"M",100,0.1,40,150,ierr);
     gMinuit->mnparm(7,"N",2.0,0.1,0.1,9,ierr);

     arglis[0] = 1000;
     arglis[1] = 1;
     gMinuit->mnexcm("MIGRAD",arglis,1,ierr);

     for (int p = 0; p<nparm; p++)
	  gMinuit->GetParameter(p, parme[j][p], eparme[j][p]);
     
     if (ierr != 0 ) {

       cout << " East bar " << j << " tw corr ierr " << ierr << ": ";
       for (int p = 0; p<nparm; p++) cout << parme[j][p] << " ";
       cout << endl;

       twe_alphas_diff << j;
       for (int p = 0; p<nparm; p++) 
	 twe_alphas_diff  << setw(13) << -999. << setw(13) << -999.;
       twe_alphas_diff << endl;

       for (int p = 0; p<nparm; p++) parme[j][p] = -999.;		 

     } else {
      
       int pnan( 0 );			// nan protection
       for (int p = 0; p<nparm; p++)
       {  
	 if (!finite(eparme[j][p])) { pnan=1; eparme[j][p]=999; }
	 if (!finite( parme[j][p])) { pnan=1; parme[j][p]=-999; eparme[j][p]=999; }
       }
       if (pnan) cout << " east side bar " << j 
		      << " time walk corr error: nan replaced (999)" << endl;

       if (parme[j][1]<40. || parme[j][1]>110) parme[j][1] = -999.; 
       if (parme[j][2]<40. || parme[j][2]>110) parme[j][2] = -999.; 
       
       twe_alphas_diff << j;
       for (int p = 0; p<nparm; p++) 
	 twe_alphas_diff  << setw(13) << parme[j][p] << setw(13) << eparme[j][p];
       twe_alphas_diff << endl;
       
     }

     double woff;
     double wslope = ajuste.fitlin3(ADC_TDCW[j], p1, woff, p2, p3, p4, p5, p6);
     
     if (wslope < slpmin || wslope > slpmax ) wslope = 0.5*(slpmin+slpmax);
     if (  woff < offmin ||   woff > offmax )   woff = 0.5*(offmin+offmax);

     twevts = 0;
     for (int k=0; k<ntwevts; k++) {
       if (fabs(dtw_bar[k][j]-woff-wslope*(1./sqrt(qw2_bar[k][j])-1./sqrt(qw1_bar[k][j])))>2.) continue; 
       z1_trk[twevts] = z1_bar[k][j];
       z2_trk[twevts] = z2_bar[k][j];
       q1_trk[twevts] = qw1_bar[k][j];
       q2_trk[twevts] = qw2_bar[k][j];
       dt_trk[twevts] = dtw_bar[k][j];
       ADC_TDCW_FIT[j]->Fill((sqrt(q1_trk[twevts])-sqrt(q2_trk[twevts]))/sqrt(q1_trk[twevts]*q2_trk[twevts]),dt_trk[twevts]);
       twevts++;
     }

     
     gMinuit->mncler();
     gMinuit->SetFCN(fcntw);
     gMinuit->mnparm(0,"A",woff,0.1,offmin,offmax,ierr);
     gMinuit->mnparm(1,"B",wslope,1,slpmin,slpmax,ierr);
     gMinuit->mnparm(2,"C",wslope,1,slpmin,slpmax,ierr);
     gMinuit->mnparm(3,"F",0.3,0.2,0.1,4,ierr);
     gMinuit->mnparm(4,"K",0.9,0.01,0.6,1,ierr);
     gMinuit->mnparm(5,"L",0,0.1,offmin,offmax,ierr);
     gMinuit->mnparm(6,"M",100,0.1,40,150,ierr);
     gMinuit->mnparm(7,"N",2.0,0.1,0.1,9,ierr);


     arglis[0] = 1000;
     arglis[1] = 1;
     gMinuit->mnexcm("MIGRAD",arglis,1,ierr);

     for (int p = 0; p<nparm; p++)
	  gMinuit->GetParameter(p, parmw[j][p], eparmw[j][p]);

     if (ierr != 0) {

       cout << " West bar " << j << " tw corr ierr " << ierr << ": ";
       for (int p = 0; p<nparm; p++) cout << parmw[j][p] << " ";
       cout << endl;
       
       tww_alphas_diff << j;
       for (int p = 0; p<nparm; p++) 
	 tww_alphas_diff  << setw(13) << -999. << setw(13) << -999.;
       tww_alphas_diff << endl;

       for (int p = 0; p<nparm; p++) parmw[j][p] = -999.;	

     } else {
      
       int pnan( 0 );			// nan protection
       for (int p = 0; p<nparm; p++)
       {
	 if (!finite(eparmw[j][p])) { pnan=1; eparmw[j][p]=999; }
	 if (!finite( parmw[j][p])) { pnan=1; parmw[j][p]=-999; eparmw[j][p]=999; }
       }
       if (pnan) cout << " west side bar " << j
		      << " time walk corr error: nan replaced (999)" << endl;

       if (parmw[j][1]<40. || parmw[j][1]>110) parmw[j][1] = -999.; 
       if (parmw[j][2]<40. || parmw[j][2]>110) parmw[j][2] = -999.; 
       
       tww_alphas_diff << j;
       for (int p = 0; p<nparm; p++) 
	 tww_alphas_diff  << setw(13) << parmw[j][p] << setw(13) << eparmw[j][p];
       tww_alphas_diff << endl;
       
     }

   }
   
   twe_alphas_diff.close();
   tww_alphas_diff.close();
     
   
   sprintf(file,"%s/files/%s_timew-calib_e.txt",directory,file_name);
   twe_calib.open(file,ios::out);
   
   sprintf(file,"%s/files/%s_timew-calib_w.txt",directory,file_name);
   tww_calib.open(file,ios::out); 

   for (int j = 0; j<216; j++) {

     int bar1 = j;
     int bar2 = j - 1; if (j==0) bar2 = 215;

     // east
     twe_calib << bar1;
     int p1[4] = {1, 3, 6, 8};
     int p2[4] = {2, 4, 7, 9};
     for (int p = 0; p<1; p++) {
       double mean, emean;
       if (parme[bar1][p1[p]]!=-999 && parme[bar2][p2[p]]!=-999) {
	 mean = (parme[bar1][p1[p]]+parme[bar2][p2[p]])/2.;
	 emean = sqrt(pow(eparme[bar1][p1[p]],2) + pow(eparme[bar2][p2[p]],2))/2.;
       } else if (parme[bar1][p1[p]]!=-999 && parme[bar2][p2[p]]==-999) {
	 mean = parme[bar1][p1[p]];
	 emean = eparme[bar1][p1[p]];
       } else if (parme[bar1][p1[p]]==-999 && parme[bar2][p2[p]]!=-999) {
	 mean = parme[bar2][p2[p]];
	 emean = eparme[bar2][p2[p]];
       } else {
	 mean = -999;
	 emean = -999;
       }
       twe_calib << setw(13) << mean << setw(13) << emean;
     }
     twe_calib << endl;
     
     // west
     tww_calib << bar1;
     for (int p = 0; p<1; p++) {
       double mean, emean;
       if (parmw[bar1][p1[p]]!=-999 && parmw[bar2][p2[p]]!=-999) {
	 mean = (parmw[bar1][p1[p]]+parmw[bar2][p2[p]])/2.;
	 emean = sqrt(pow(eparmw[bar1][p1[p]],2) + pow(eparmw[bar2][p2[p]],2))/2.;
       } else if (parmw[bar1][p1[p]]!=-999 && parmw[bar2][p2[p]]==-999) {
	 mean = parmw[bar1][p1[p]];
	 emean = eparmw[bar1][p1[p]];
       } else if (parmw[bar1][p1[p]]==-999 && parmw[bar2][p2[p]]!=-999) {
	 mean = parmw[bar2][p2[p]];
	 emean = eparmw[bar2][p2[p]];
       } else {
	 mean = -999;
	 emean = -999;
       }
       tww_calib << setw(13) << mean << setw(13) << emean;
     }
     tww_calib << endl;

   }
   
   twe_calib.close();
   tww_calib.close();
   
 } // if (flag==0)

 hfile->Write();
 hfile->Close();
   
} // timew_calib::LoopTW