#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.h>
#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 beta1    = par[3];
//   double beta2    = par[4];
//   double sigma    = par[5];
  double sigma    = par[3];
//   double alfaz1   = par[6];
//   double alfaz2   = par[7];
//   double betaz1   = par[8];
//   double betaz2   = par[9];
//   double gamma1   = par[10];
//   double gamma2   = par[11];
//   double gammz1   = par[12];
//   double gammz2   = par[13];
  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 z1    = z1_trk[indx];
        //double z2    = z2_trk[indx];
	double alfat1= alfa1;// + z1*alfaz1 + (beta1+z1*betaz1)/sqrt(q1);// + (gamma1+z1*gammz1)/q1;// + gamma/q1;
	double alfat2= alfa2;// + z2*alfaz2 + (beta2+z2*betaz2)/sqrt(q2);// + (gamma2+z2*gammz2)/q2;// + gamma/q2;
// 	double alfat1= alfa1 + z1*beta1 + (alfa1+z1*beta1)/sqrt(q1);// + gamma/q1;
// 	double alfat2= alfa2 + z2*beta2 + (alfa2+z2*beta2)/sqrt(q2);// + gamma/q1;
// 	double func  = offset  - alfat1/sqrt(q1)  + alfat2/sqrt(q2) + alfa2*(1./sqrt(q2)-1./sqrt(q1))*(1./sqrt(q2)-1./sqrt(q1))+ beta1*(1./sqrt(q2)-1./sqrt(q1))*(1./sqrt(q2)-1./sqrt(q1))*(1./sqrt(q2)-1./sqrt(q1));
	double func  = offset  - alfat1/sqrt(q1)  + alfat2/sqrt(q2); 
	double de    = dt - func;
// 	double arg   = exp(-de*de/(2*sigma*sigma))/(r2pi*sigma);
	double db    = dt - (offsetb - alfab/sqrt(q1)  + alfab/sqrt(q2));
	double arg   = frac*exp(-de*de/(2*sigma*sigma))/(r2pi*sigma) + (1.-frac)*exp(-db*db/(2*sigmab*sigmab))/(r2pi*sigmab);
// 	double arg   = exp(-de*de/(2*sigma*sigma))/(r2pi*sigma);
	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;

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

    if (jentry%100000 == 0) printf(" EventCounter:track %8d (run %d) \n",jentry,timewalk_run);
    //if (jentry%10000 == 0)cout <<  timewalk_d0wrt << ", "<<  timewalk_inBar << ", " 
    //                          << timewalk_run << ", " <<  timewalk_event << endl;

    if (timewalk_event != currentEvent) {
      currentEvent = timewalk_event;
      ++eventCounter;
    }
  } // 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<216; nb++) twevtbar[nb] = 0;
  double adc_comb_e, adc_comb_w, difft1, difft2, pl1, z1, z2;
  //double zadc_comb_e, zadc_comb_w;
  double resulte[NumBar][mparm], resultw[NumBar][mparm], resulsp[NumBar][2];
  double datae[mparm], dataw[mparm], datasp[9];

  //double erroffe[NumBar], erpae[NumBar], erpce[NumBar], erpde[NumBar];
  //double par_Ae[NumBar], par_Be[NumBar], par_Ce[NumBar], par_De[NumBar];
  //double erroffw[NumBar], erpaw[NumBar], erpcw[NumBar], erpdw[NumBar];
  //double par_Aw[NumBar], par_Bw[NumBar], par_Cw[NumBar], par_Dw[NumBar];

  //double erroffe2file[NumBar], erpae2file[NumBar], erpce2file[NumBar], erpde2file[NumBar];
  //double par_Ae2file[NumBar], par_Be2file[NumBar], par_Ce2file[NumBar], par_De2file[NumBar];
  //double erroffw2file[NumBar], erpaw2file[NumBar], erpcw2file[NumBar], erpdw2file[NumBar];
  //double par_Aw2file[NumBar], par_Bw2file[NumBar], par_Cw2file[NumBar], par_Dw2file[NumBar];

  //double cove[3][NumBar], covw[3][NumBar], csqe[NumBar], csqw[NumBar];

  int cont[NumBar];

  // Book histograms
  TH1F *para    = new TH1F("para",   "Par_A",  100, 50,120);
  TH1F *parb    = new TH1F("parb",   "Par_B",  100,-10, 10);
  TH1F *para1   = new TH1F("para1",  "Par_A1", 100, 50,120);
  TH1F *parb1   = new TH1F("parb1",  "Par_B1", 100,-10, 10);
  TH1F *parape  = new TH1F("parape", "parape", 100, 50,120);
  TH1F *parapw  = new TH1F("parapw", "parapw", 100, 50,120);
  TH1F *alphas_diff_e  = new TH1F("alphas_diff_e", "alphas_diff_e", 100, -20,20);
  TH1F *alphas_diff_w  = new TH1F("alphas_diff_w", "alphas_diff_w", 100, -20,20);
  TH1F *chime   = new TH1F("chime",  "chime",   40,  0, 10);
  TH1F *chimw   = new TH1F("chimw",  "chimw",   40,  0, 10);
  TH1F *residme = new TH1F("residme","residme",100, -3,  3);
  TH1F *residmw = new TH1F("residmw","residmw",100, -3,  3);

  TH1F *diff    = new TH1F("diff",   "diff",   120,-30, 30);
  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);
//     sprintf(file,"%s/files/%s_timew-fitres_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;
//     sprintf(file,"%s/files/%s_timew-fitres_w.txt",directory,file_name);
    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++) {
    //par_Ae[j] = 0;
    //par_Aw[j] = 0;
    //par_Be[j] = 0;
    //par_Bw[j] = 0;
    //par_Ce[j] = 0;
    //par_Cw[j] = 0;
    //par_De[j] = 0;
    //par_Dw[j] = 0;
 
    //par_Ae2file[j] = 0;
    //par_Aw2file[j] = 0;
    //par_Be2file[j] = 0;
    //par_Bw2file[j] = 0;
    //par_Ce2file[j] = 0;
    //par_Cw2file[j] = 0;
    //par_De2file[j] = 0;
    //par_Dw2file[j] = 0;

    //erroffe2file[j] = 0;
    //erpae2file  [j] = 0;
    //erpce2file  [j] = 0;
    //erpde2file  [j] = 0;
    //erroffw2file[j] = 0;
    //erpaw2file  [j] = 0;
    //erpcw2file  [j] = 0;
    //erpdw2file  [j] = 0;
/*
    for(int i=0; i<9; i++) {
      datae [i] = 0;
      dataw [i] = 0;
      datasp[i] = 0;
      twe_result >> datae[i];
      tww_result >> dataw[i];
    }

    resulte[j][0] = datae[1];
    resulte[j][1] = datae[3];
    resultw[j][0] = dataw[1];
    resultw[j][1] = dataw[3];
*/
// 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();

  // 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 (jentry%100000 == 0) printf(" LoopTW:track %8d (run %d) \n",jentry,timewalk_run);

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

    // Track selection
    if (   //(fabs(timewalk_d0wrt)  <  1.5)  
         (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)
       ) {

      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);

//       if (timewalk_qe2<400. || timewalk_qe1<400. || timewalk_qw2<400. || timewalk_qw1<400.)  continue;
//       if (timewalk_qe2>1500. || timewalk_qe1>1500. || timewalk_qw2>1500. || timewalk_qw1>1500.)  continue;
//       if (fabs(timewalk_z1)>20. || fabs(timewalk_z2)>20.) continue;

      int bar1 = int(timewalk_bar1);
      int bar2 = int(timewalk_bar2);
      ++cont[bar1];

      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;

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

      if (timewalk_qe1 > 0 && timewalk_qe2 > 0) {
        adc_comb_e = (sqrt(timewalk_qe1)-sqrt(timewalk_qe2))/
                     (sqrt(timewalk_qe1)*sqrt(timewalk_qe2));
        //zadc_comb_e= (sqrt(timewalk_qe1)*timewalk_z2-sqrt(timewalk_qe2)*timewalk_z1)/
        //             (sqrt(timewalk_qe1)*sqrt(timewalk_qe2));
      }

      if (timewalk_qw1 > 0 && timewalk_qw2 > 0) {
        adc_comb_w = (sqrt(timewalk_qw1)-sqrt(timewalk_qw2))/
                     (sqrt(timewalk_qw1)*sqrt(timewalk_qw2));
        //zadc_comb_w= (sqrt(timewalk_qw1)*timewalk_z2-sqrt(timewalk_qw2)*timewalk_z1)/
        //             (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);
      if (twevtbar[bar1]<maxevtnum) {
	z1_bar[twevtbar[bar1]][bar1]  = timewalk_z1;
	z2_bar[twevtbar[bar1]][bar1]  = timewalk_z2;
	qe1_bar[twevtbar[bar1]][bar1] = timewalk_qe1;
	qe2_bar[twevtbar[bar1]][bar1] = timewalk_qe2;
	qw1_bar[twevtbar[bar1]][bar1] = timewalk_qw1;
	qw2_bar[twevtbar[bar1]][bar1] = timewalk_qw2;
	dte_bar[twevtbar[bar1]][bar1] = difft1;
	dtw_bar[twevtbar[bar1]][bar1] = difft2;
	twevtbar[bar1]++;
      }
   
      double e0 = 0., e1 = 0., w0 = 0., w1 = 0.;
      if (flag==1) {
	//       double e0 = (resulte[bar1][0]+resulte[bar1][2]/sqrt(timewalk_qe2)+resulte[bar1][3]/timewalk_qe2)/sqrt(timewalk_qe2) - 
	// 	          (resulte[bar1][0]+resulte[bar1][2]/sqrt(timewalk_qe1)+resulte[bar1][3]/timewalk_qe1)/sqrt(timewalk_qe1);
	//double e0 = (resulte[bar1][0])*adc_comb_e;
// 	e0  = (resulte[bar1][4]+timewalk_z*resulte[bar1][14])/sqrt(timewalk_qe2) - (resulte[bar1][2]+timewalk_z1*resulte[bar1][12])/sqrt(timewalk_qe1);
// 	e0 += (resulte[bar1][8]+timewalk_z*resulte[bar1][18])/timewalk_qe2	 - (resulte[bar1][6]+timewalk_z1*resulte[bar1][16])/timewalk_qe1;
	e0  = resulte[bar1][4]/sqrt(timewalk_qe2) - resulte[bar1][2]/sqrt(timewalk_qe1);
// 	e0 += (resulte[bar1][22]+timewalk_z2*resulte[bar1][26])/(sqrt(timewalk_qe2)*timewalk_qe2)	
// 	  -   (resulte[bar1][20]+timewalk_z1*resulte[bar1][24])/(sqrt(timewalk_qe1)*timewalk_qe1);
// 	e0  = (resulte[bar1][4])/sqrt(timewalk_qe2) - (resulte[bar1][2])/sqrt(timewalk_qe1);
// 	e0 += (resulte[bar1][8])/timewalk_qe2       - (resulte[bar1][6])/timewalk_qe1;
// 	e0  = (resulte[bar1][2])/sqrt(timewalk_qe2) - (resulte[bar1][2])/sqrt(timewalk_qe1) + resulte[bar1][4]*adc_comb_e*adc_comb_e;
	//e0 += 0.*resulte[bar1][2]*zadc_comb_e;
	e1 = resulte[bar1][0];
	//       double w0 = (resultw[bar1][0]+resultw[bar1][2]/sqrt(timewalk_qw2)+resultw[bar1][3]/timewalk_qw2)/sqrt(timewalk_qw2) - 
	// 	          (resultw[bar1][0]+resultw[bar1][2]/sqrt(timewalk_qw1)+resultw[bar1][3]/timewalk_qw1)/sqrt(timewalk_qw1);
	//double w0 = (resultw[bar1][0])*adc_comb_w;
// 	w0  = (resultw[bar1][4]+timewalk_z2*resultw[bar1][14])/sqrt(timewalk_qw2) - (resultw[bar1][2]+timewalk_z1*resultw[bar1][12])/sqrt(timewalk_qw1);
// 	w0 += (resultw[bar1][8]+timewalk_z2*resultw[bar1][18])/timewalk_qw2	  - (resultw[bar1][6]+timewalk_z1*resultw[bar1][16])/timewalk_qw1;
// 	w0  = (resultw[bar1][4])/sqrt(timewalk_qw2) - (resultw[bar1][2])/sqrt(timewalk_qw1);
// 	w0 += (resultw[bar1][8])/timewalk_qw2       - (resultw[bar1][6])/timewalk_qw1;
// 	w0  = (resultw[bar1][2])/sqrt(timewalk_qw2) - (resultw[bar1][2])/sqrt(timewalk_qw1) + resultw[bar1][4]*adc_comb_w*adc_comb_w;
	w0  = resultw[bar1][4]/sqrt(timewalk_qw2) - resultw[bar1][2]/sqrt(timewalk_qw1);
	//w0 += 0.*resultw[bar1][2]*zadc_comb_w;
	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));
//       residgq->Fill(timewalk_qw1, difft2-(w0+w1));
      
      residgz->Fill(timewalk_z1, difft1-(e0+e1));

//       residgx->Fill(adc_comb_e, difft1-(e0+e1));
//       residgx->Fill(adc_comb_w, difft2-(w0+w1));
   
      int zbin = (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 = (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 = (timewalk_qe1/800.);
	    if (qbin2>=0 && qbin2<=4) residgz_q[qbin2]->Fill(timewalk_z1, difft1-(e0+e1));
	    int xbin2 = (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));
	}
      }
 
    } // track selection
  } // loop over tracks

//   cout << "Performing linear fits for time walk on profiles..." << endl;

//   double err [NumBar], errsig [NumBar], sigma [NumBar], csq [NumBar];
//   double err2[NumBar], errsig2[NumBar], sigma2[NumBar], csq2[NumBar];

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

//     par_Ae[i] = ajuste.fitlin3(ADC_TDCE[i], erpae[i], par_Be[i], erroffe[i],
//                                csqe[i], cove[0][i], cove[1][i], cove[2][i]);

//     ajuste.gaus1(reside[i],-3.5,3.5,err[i],sigma[i],errsig[i],csq[i]);
//     residme->Fill(err[i]); 
//     if(par_Ae[i] != 0) {
//       para->Fill(par_Ae[i]);
//       chime->Fill(csqe[i]);
//     }

//     par_Aw[i] = ajuste.fitlin3(ADC_TDCW[i], erpaw[i], par_Bw[i], erroffw[i],
// 			       csqw[i], covw[0][i], covw[1][i], covw[2][i]);
			         
//     ajuste.gaus1(residw[i],-3.5,3.5,err2[i],sigma2[i],errsig2[i],csq2[i]);
//     residmw->Fill(err2[i]); 
//     if(par_Aw[i] !=0) {
//       para1->Fill(par_Aw[i]);
//       chimw->Fill(csqw[i]);
//     }
//   } // loop over bars

 if (flag==0) {

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

//    sprintf(file,"%s/files/%s_timew-fitres_e.txt",directory,file_name);
//    twe_fit.open(file);
   
//    sprintf(file,"%s/files/%s_timew-fitres_w.txt",directory,file_name);
//    tww_fit.open(file);
   
   sprintf(file,"%s/files/%s_timew-alphas_diff_e.txt",directory,file_name);
   twe_alphas_diff.open(file); 
   
   sprintf(file,"%s/files/%s_timew-alphas_diff_w.txt",directory,file_name);
   tww_alphas_diff.open(file);

   double parme[216][nparm], parmw[216][nparm], eparme[216][nparm], eparmw[216][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++) {
//    for (int j=0 ;j<1; j++) {
     
     int ntwevts = twevtbar[j];

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

     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++;
     }

     if (maxevtnum<twevts) cout << "Attention: number of tracks in bar " << j << "!" << endl;
     
     gMinuit->mncler();
     gMinuit->SetFCN(fcntw);
     gMinuit->mnparm(0,"A",eoff,0.1,-5,5,ierr);
     gMinuit->mnparm(1,"B",eslope,1,0,150,ierr);
     gMinuit->mnparm(2,"C",eslope,1,0,150,ierr);
//      gMinuit->mnparm(2,"C",250.,1,-1500,1500,ierr);
//      gMinuit->mnparm(3,"D",5000.,10,-15000,15000,ierr);
//      gMinuit->mnparm(3,"D",200.,1,0,1000,ierr);
//      gMinuit->mnparm(4,"E",200.,1,0,1000,ierr);
//      gMinuit->mnparm(5,"F",0.5,0.1,0,4,ierr);
     gMinuit->mnparm(3,"F",0.5,0.1,0,4,ierr);
//      gMinuit->mnparm(6,"G",0.,0.01,-1,1,ierr);
//      gMinuit->mnparm(7,"H",0.,0.01,-1,1,ierr);
//      gMinuit->mnparm(8,"I",0.,0.2,-20,20,ierr);
//      gMinuit->mnparm(9,"J",0.,0.2,-20,20,ierr);
//      gMinuit->mnparm(10,"K",0.,10,-5000,5000,ierr);
//      gMinuit->mnparm(11,"L",0.,10,-5000,5000,ierr);
//      gMinuit->mnparm(12,"M",0.,10,-5000,5000,ierr);
//      gMinuit->mnparm(13,"N",0.,10,-5000,5000,ierr);
     gMinuit->mnparm(4,"K",0.9,0.01,0.5,1,ierr);
     gMinuit->mnparm(5,"L",0.,0.1,-5,5,ierr);
     gMinuit->mnparm(6,"M",80.,0.1,40,150,ierr);
     gMinuit->mnparm(7,"N",2.0,0.1,0,9,ierr);

     arglis[0] = 1000;
     arglis[1] = 1;
//      cout << "Fitting channel " << j << endl;
     gMinuit->mnexcm("MIGRAD",arglis,1,ierr);

     if (ierr != 0 ) {

       cout << " east side bar " << j << " time walk corrections: this failed because ierr = " << ierr << 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 {
      
       for (int p = 0; p<nparm; p++)  
	 gMinuit->GetParameter(p, parme[j][p], eparme[j][p]);
       //      gMinuit->GetParameter(0,par_Be[j],erroffe[j]);
       //      gMinuit->GetParameter(1,par_Ae[j],erpae[j]  );
       //      gMinuit->GetParameter(2,csqe[j],  cove[0][j]);
       //      gMinuit->GetParameter(3,par_Ce[j],erpce[j]  );
       //      gMinuit->GetParameter(8,par_De[j],erpde[j]  );
//        parme[j][0] = 1.05780e-01; 
//        parme[j][1] = 8.47824e+01;
//        parme[j][2] = -1.11383e+02;

       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);

     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++;
     }

     if (maxevtnum<twevts) cout << "Attention: number of tracks in bar " << j+216 << "!" << endl;
     
     gMinuit->mncler();
     gMinuit->SetFCN(fcntw);
     gMinuit->mnparm(0,"A",eoff,0.1,-5,5,ierr);
     gMinuit->mnparm(1,"B",eslope,1,0,150,ierr);
     gMinuit->mnparm(2,"C",eslope,1,0,150,ierr);
//      gMinuit->mnparm(2,"C",250.,1,-1500,1500,ierr);
//      gMinuit->mnparm(3,"D",5000.,10,-15000,15000,ierr);
//      gMinuit->mnparm(3,"D",200.,1,0,1000,ierr);
//      gMinuit->mnparm(4,"E",200.,1,0,1000,ierr);
//      gMinuit->mnparm(5,"F",0.5,0.1,0,4,ierr);
     gMinuit->mnparm(3,"F",0.5,0.1,0,4,ierr);
//      gMinuit->mnparm(6,"G",0.,0.01,-1,1,ierr);
//      gMinuit->mnparm(7,"H",0.,0.01,-1,1,ierr);
//      gMinuit->mnparm(8,"I",0.,0.2,-20,20,ierr);
//      gMinuit->mnparm(9,"J",0.,0.2,-20,20,ierr);
//      gMinuit->mnparm(10,"K",0.,10,-5000,5000,ierr);
//      gMinuit->mnparm(11,"L",0.,10,-5000,5000,ierr);
//      gMinuit->mnparm(12,"M",0.,10,-5000,5000,ierr);
//      gMinuit->mnparm(13,"N",0.,10,-5000,5000,ierr);
     gMinuit->mnparm(4,"K",0.9,0.01,0.5,1,ierr);
     gMinuit->mnparm(5,"L",0.,0.1,-5,5,ierr);
     gMinuit->mnparm(6,"M",80.,0.1,40,150,ierr);
     gMinuit->mnparm(7,"N",2.0,0.1,0,9,ierr);


     arglis[0] = 1000;
     arglis[1] = 1;
//      cout << "Fitting channel " << j << endl;
     gMinuit->mnexcm("MIGRAD",arglis,1,ierr);

     if (ierr != 0) {

       cout << " west side bar " << j << " time walk corrections: this failed because ierr = " << ierr << 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 {
      
       for (int p = 0; p<nparm; p++)  
	 gMinuit->GetParameter(p, parmw[j][p], eparmw[j][p]);
       //      gMinuit->GetParameter(0,par_Be[j],erroffe[j]);
       //      gMinuit->GetParameter(1,par_Ae[j],erpae[j]  );
       //      gMinuit->GetParameter(2,csqe[j],  cove[0][j]);
       //      gMinuit->GetParameter(3,par_Ce[j],erpce[j]  );
       //      gMinuit->GetParameter(8,par_De[j],erpde[j]  );
//        parme[j][0] = 1.05780e-01; 
//        parme[j][1] = 8.47824e+01;
//        parme[j][2] = -1.11383e+02;

       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);
   
   sprintf(file,"%s/files/%s_timew-calib_w.txt",directory,file_name);
   tww_calib.open(file); 

   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();
   
//    sprintf(file,"%s/files/%s_timew-alphas_diff_e.txt",directory,file_name);
//    twe_alphas_diff.open(file);
   
//    sprintf(file,"%s/files/%s_timew-alphas_diff_w.txt",directory,file_name);
//    tww_alphas_diff.open(file);
   
//    // Average over pairs of adjacent bars
//    // We start with the par 215-0 and 0-1 for the bar 0
   
//    double min_alpha =  52;
//    double max_alpha = 103;
//    // bar 215
//    if (par_Ae[214] > min_alpha && par_Ae[215] > min_alpha &&
//        par_Ae[214] < max_alpha && par_Ae[215] < max_alpha) {
//      par_Ae2file [215] = (par_Ae[215]+par_Ae[0]) / 2.0;   
//      par_Be2file [215] = (par_Be[215]+par_Be[0]) / 2.0;   
//      par_Ce2file [215] = (par_Ce[215]+par_Ce[0]) / 2.0;   
//      par_De2file [215] = (par_De[215]+par_De[0]) / 2.0;
//      erpae2file  [215] = sqrt(pow(erpae[215],2) + pow(erpae[214],2)) / 2.0;
//      erroffe2file[215] = sqrt(pow(erroffe[215],2) + pow(erroffe[214],2)) / 2.0;
//      erpce2file  [215] = sqrt(pow(erpce[215],2) + pow(erpce[214],2)) / 2.0;
//      erpde2file  [215] = sqrt(pow(erpde[215],2) + pow(erpde[214],2)) / 2.0;
//    }else if (par_Ae[214] <= min_alpha || par_Ae[214] >= max_alpha){
//      printf("Warning! Alpha PMT's 214-215 (east) = %f\n",par_Ae[0]);
//      par_Ae2file [215] = par_Ae[215];   
//      par_Be2file [215] = par_Be[215];   
//      par_Ce2file [215] = par_Ce[215];   
//      par_De2file [215] = par_De[215];
//      erpae2file  [215] = erpae[215];
//      erroffe2file[215] = erroffe[215];
//      erpce2file  [215] = erpce[215];
//      erpde2file  [215] = erpde[215];
//    }else if (par_Ae[215] <= min_alpha || par_Ae[215] >= max_alpha){
//      printf("Warning! Alpha PMT's 215-0 (east) = %f\n",par_Ae[215]);
//      par_Ae2file [215] = par_Ae[214];   
//      par_Be2file [215] = par_Be[214];   
//      par_Ce2file [215] = par_Ce[214];  
//      par_De2file [215] = par_De[214];
//      erpae2file  [215] = erpae[214];
//      erroffe2file[215] = erroffe[214];
//      erpce2file  [215] = erpce[214];
//      erpde2file  [215] = erpde[214];
//    }
   
//    if (par_Aw[214] > min_alpha && par_Aw[215] > min_alpha &&
//        par_Aw[214] < max_alpha && par_Aw[215] < max_alpha) {
//      par_Aw2file [214] = (par_Aw[215]+par_Aw[214]) / 2.0;
//      par_Bw2file [214] = (par_Bw[215]+par_Bw[214]) / 2.0;
//      par_Cw2file [214] = (par_Cw[215]+par_Cw[214]) / 2.0;
//      par_Dw2file [214] = (par_Dw[215]+par_Dw[214]) / 2.0;
//      erpaw2file  [214] = sqrt(pow(erpaw[215],2) + pow(erpaw[214],2)) / 2.0;
//      erroffw2file[214] = sqrt(pow(erroffw[215],2) + pow(erroffw[214],2)) / 2.0;
//      erpcw2file  [214] = sqrt(pow(erpcw[215],2) + pow(erpcw[214],2)) / 2.0;
//      erpdw2file  [214] = sqrt(pow(erpdw[215],2) + pow(erpdw[214],2)) / 2.0;
//    }else if (par_Aw[214] < min_alpha || par_Aw[214] > max_alpha){
//      printf("Warning! Alpha PMT's 214-215 (west) = %f\n",par_Aw[0]);
//      par_Aw2file [214] = par_Aw[215];
//      par_Bw2file [214] = par_Bw[215];
//      par_Cw2file [214] = par_Cw[215];
//      par_Dw2file [214] = par_Dw[215];
//      erpaw2file  [214] = erpaw[215];
//      erroffw2file[214] = erroffw[215];
//      erpcw2file  [214] = erpcw[215];
//      erpdw2file  [214] = erpdw[215];
//    }else if (par_Aw[215] < min_alpha || par_Aw[214] > max_alpha){
//      printf("Warning! Alpha PMT's 215-0 (west) = %f\n",par_Aw[215]);
//      par_Aw2file [214] = par_Aw[214];
//      par_Bw2file [214] = par_Bw[214];
//      par_Cw2file [214] = par_Cw[214];
//      par_Dw2file [214] = par_Dw[214];
//      erpaw2file  [214] = erpaw[214];
//      erroffw2file[214] = erroffw[214];
//      erpcw2file  [214] = erpcw[214];
//      erpdw2file  [214] = erpdw[214];
//    }
   
//    // bar 0
//    if (par_Ae[215] > min_alpha && par_Ae[0] > min_alpha &&
//        par_Ae[215] < max_alpha && par_Ae[0] < max_alpha) {
//      par_Ae2file [0] = (par_Ae[215]+par_Ae[0]) / 2.0;   
//      par_Be2file [0] = (par_Be[215]+par_Be[0]) / 2.0;  
//      par_Ce2file [0] = (par_Ce[215]+par_Ce[0]) / 2.0; 
//      par_De2file [0] = (par_De[215]+par_De[0]) / 2.0;
//      erpae2file  [0] = sqrt(pow(erpae[215],2) + pow(erpae[0],2)) / 2.0;
//      erroffe2file[0] = sqrt(pow(erroffe[215],2) + pow(erroffe[0],2)) / 2.0;
//      erpce2file  [0] = sqrt(pow(erpce[215],2) + pow(erpce[0],2)) / 2.0;
//      erpde2file  [0] = sqrt(pow(erpde[215],2) + pow(erpde[0],2)) / 2.0;
//    }else if (par_Ae[0] <= min_alpha || par_Ae[0] >= max_alpha){
//      printf("Warning! Alpha PMT's 0-1 (east) = %f\n",par_Ae[0]);
//      par_Ae2file [0] = par_Ae[215];   
//      par_Be2file [0] = par_Be[215];   
//      par_Ce2file [0] = par_Ce[215];   
//      par_De2file [0] = par_De[215];
//      erpae2file  [0] = erpae[215];
//      erroffe2file[0] = erroffe[215];
//      erpce2file  [0] = erpce[215];
//      erpde2file  [0] = erpde[215];
//    }else if (par_Ae[215] <= min_alpha || par_Ae[215] >= max_alpha){
//      printf("Warning! Alpha PMT's 215-0 (east) = %f\n",par_Ae[215]);
//      par_Ae2file [0] = par_Ae[0];   
//      par_Be2file [0] = par_Be[0];  
//      par_Ce2file [0] = par_Ce[0];  
//      par_De2file [0] = par_De[0];
//      erpae2file  [0] = erpae[0];
//      erroffe2file[0] = erroffe[0];
//      erpce2file  [0] = erpce[0];
//      erpde2file  [0] = erpde[0];
//    }
   
//    if (par_Aw[215] > min_alpha && par_Aw[0] > min_alpha &&
//        par_Aw[215] < max_alpha && par_Aw[0] < max_alpha) {
//      par_Aw2file [0] = (par_Aw[215]+par_Aw[0]) / 2.0;
//      par_Bw2file [0] = (par_Bw[215]+par_Bw[0]) / 2.0;
//      par_Cw2file [0] = (par_Cw[215]+par_Cw[0]) / 2.0;
//      par_Dw2file [0] = (par_Dw[215]+par_Dw[0]) / 2.0;
//      erpaw2file  [0] = sqrt(pow(erpaw[215],2) + pow(erpaw[0],2)) / 2.0;
//      erroffw2file[0] = sqrt(pow(erroffw[215],2) + pow(erroffw[0],2)) / 2.0;
//      erpcw2file  [0] = sqrt(pow(erpcw[215],2) + pow(erpcw[0],2)) / 2.0;
//      erpdw2file  [0] = sqrt(pow(erpdw[215],2) + pow(erpdw[0],2)) / 2.0;
//    }else if (par_Aw[0] < min_alpha || par_Aw[0] > max_alpha){
//      printf("Warning! Alpha PMT's 0-1 (west) = %f\n",par_Aw[0]);
//      par_Aw2file [0] = par_Aw[215];
//      par_Bw2file [0] = par_Bw[215];
//      par_Cw2file [0] = par_Cw[215];
//      par_Dw2file [0] = par_Dw[215];
//      erpaw2file  [0] = erpaw[215];
//      erroffw2file[0] = erroffw[215];
//      erpcw2file  [0] = erpcw[215];
//      erpdw2file  [0] = erpdw[215];
//    }else if (par_Aw[215] < min_alpha || par_Aw[215] > max_alpha){
//      printf("Warning! Alpha PMT's 215-0 (west) = %f\n",par_Aw[215]);
//      par_Aw2file [0] = par_Aw[0];
//      par_Bw2file [0] = par_Bw[0];
//      par_Cw2file [0] = par_Cw[0];
//      par_Dw2file [0] = par_Dw[0];
//      erpaw2file  [0] = erpaw[0];
//      erroffw2file[0] = erroffw[0];
//      erpcw2file  [0] = erpcw[0];
//      erpdw2file  [0] = erpdw[0];
//    }
   
//    twe_alphas_diff << setw(13) << 0
// 		   << setw(13) << par_Ae[215]
// 		   << setw(13) << erpae[215]
// 		   << setw(13) << par_Ae[0]
// 		   << setw(13) << erpae[0]
// 		   << setw(13) << par_Ae[215] - par_Ae[0]
// 		   << setw(13) << erpae[215] - erpae[0] 
// 		   << setw(13) << par_Be[0]
// 		   << setw(13) << erroffe[0]
// 		   << setw(13) << par_Ce[0]
// 		   << setw(13) << erpce[0]
// 		   << setw(13) << par_De[0]
// 		   << setw(13) << erpde[0] << endl;
   
//    tww_alphas_diff << setw(13) << 216
// 		   << setw(13) << par_Aw[215]
// 		   << setw(13) << erpaw[215]
// 		   << setw(13) << par_Aw[0]
// 		   << setw(13) << erpaw[0]
// 		   << setw(13) << par_Aw[215] - par_Aw[0]
// 		   << setw(13) << erpaw[215] - erpaw[0]
// 		   << setw(13) << par_Bw[0]
// 		   << setw(13) << erroffw[0]
// 		   << setw(13) << par_Cw[0]
// 		   << setw(13) << erpcw[0]
// 		   << setw(13) << par_Dw[0]
// 		   << setw(13) << erpdw[0] << endl;
   
//    for(int i=0; i<NumBar-1; i++) {
     
//      twe_alphas_diff << setw(13) << i+1
// 		     << setw(13) << par_Ae[i  ] << setw(13) << erpae[i  ]
// 		     << setw(13) << par_Ae[i+1] << setw(13) << erpae[i+1]
// 		     << setw(13) << par_Ae[i] - par_Ae[i+1]
// 		     << setw(13) << sqrt(pow(erpae[i],2) + pow(erpae[i+1],2)) / 2.0
// 		     << setw(13) << par_Be[i+1] << setw(13) << erroffe[i+1]
// 		     << setw(13) << par_Ce[i+1] << setw(13) << erpce[i+1]
// 		     << setw(13) << par_De[i+1] << setw(13) << erpde[i+1]
// 		     << endl;
     
//      tww_alphas_diff << setw(13) << i+1+216
// 		     << setw(13) << par_Aw[i  ] << setw(13) << erpaw[i  ]
// 		     << setw(13) << par_Aw[i+1] << setw(13) << erpaw[i+1]
// 		     << setw(13) << par_Aw[i] - par_Aw[i+1]
// 		     << setw(13) << sqrt(pow(erroffw[i],2) + pow(erroffw[i+1],2)) / 2.0
// 		     << setw(13) << par_Bw[i+1] << setw(13) << erroffw[i+1]
// 		     << setw(13) << par_Cw[i+1] << setw(13) << erpcw[i+1]
// 		     << setw(13) << par_Dw[i+1] << setw(13) << erpdw[i+1]
// 		     << endl;
     
//      parape->Fill((par_Ae[i]+par_Ae[i+1])/2.0);    
//      parapw->Fill((par_Aw[i]+par_Aw[i+1])/2.0); 
//      alphas_diff_e->Fill(par_Ae[i]-par_Ae[i+1]);    
//      alphas_diff_w->Fill(par_Aw[i]-par_Aw[i+1]);    
     
//      if (par_Ae[i] > min_alpha && par_Ae[i+1] > min_alpha && 
// 	 par_Ae[i] < max_alpha && par_Ae[i+1] < max_alpha) {
//        par_Ae2file [i+1] = (par_Ae[i]+par_Ae[i+1]) / 2.0;   
//        par_Be2file [i+1] = (par_Be[i]+par_Be[i+1]) / 2.0;   
//        par_Ce2file [i+1] = (par_Ce[i]+par_Ce[i+1]) / 2.0;   
//        par_De2file [i+1] = (par_De[i]+par_De[i+1]) / 2.0;
//        erpae2file  [i+1] = sqrt(pow(erpae[i],2) + pow(erpae[i+1],2)) / 2.0;
//        erroffe2file[i+1] = sqrt(pow(erroffe[i],2) + pow(erroffe[i+1],2)) / 2.0;
//        erpce2file  [i+1] = sqrt(pow(erpce[i],2) + pow(erpce[i+1],2)) / 2.0;
//        erpde2file  [i+1] = sqrt(pow(erpde[i],2) + pow(erpde[i+1],2)) / 2.0;
//      }else if(par_Ae[i] <= min_alpha || par_Ae[i] >= max_alpha){ 
//        //printf("Warning! Alpha PMT's %d-%d (east) = %f\n",i,i+1,par_Ae[i]);
//        par_Ae2file [i+1] = par_Ae[i+1];   
//        par_Be2file [i+1] = par_Be[i+1];  
//        par_Ce2file [i+1] = par_Ce[i+1];  
//        par_De2file [i+1] = par_De[i+1];
//        erpae2file  [i+1] = erpae[i+1];
//        erroffe2file[i+1] = erroffe[i+1];
//        erpce2file  [i+1] = erpce[i+1];
//        erpde2file  [i+1] = erpde[i+1];
//      }else if(par_Ae[i+1] <= min_alpha || par_Ae[i+1] >= max_alpha) {
//        printf("Warning 2! Alpha PMT's %d-%d (east) = %f\n",i+1,i+2,par_Ae[i+1]);
//        par_Ae2file [i+1] = par_Ae[i];   
//        par_Be2file [i+1] = par_Be[i];  
//        par_Ce2file [i+1] = par_Ce[i];
//        par_De2file [i+1] = par_De[i];
//        erpae2file  [i+1] = erpae[i];
//        erroffe2file[i+1] = erroffe[i];
//        erpce2file  [i+1] = erpce[i];
//        erpde2file  [i+1] = erpde[i];
//     }
//      if (par_Aw[i] > min_alpha && par_Aw[i+1] > min_alpha && 
// 	 par_Aw[i] < max_alpha && par_Aw[i+1] < max_alpha) {
//        par_Aw2file [i+1] = (par_Aw[i]+par_Aw[i+1]) / 2.0;
//        par_Bw2file [i+1] = (par_Bw[i]+par_Bw[i+1]) / 2.0;
//        par_Cw2file [i+1] = (par_Cw[i]+par_Cw[i+1]) / 2.0;
//        par_Dw2file [i+1] = (par_Dw[i]+par_Dw[i+1]) / 2.0;
//        erpaw2file  [i+1] = sqrt(pow(erpaw[i],2) + pow(erpaw[i+1],2)) / 2.0;
//        erroffw2file[i+1] = sqrt(pow(erroffw[i],2) + pow(erroffw[i+1],2)) / 2.0;
//        erpcw2file  [i+1] = sqrt(pow(erpcw[i],2) + pow(erpcw[i+1],2)) / 2.0;
//        erpdw2file  [i+1] = sqrt(pow(erpdw[i],2) + pow(erpdw[i+1],2)) / 2.0;
//      }else if(par_Aw[i] <= min_alpha || par_Aw[i] >= max_alpha){
//        printf("Warning 1! Alpha PMT's %d-%d (west) = %f\n",i,i+1,par_Aw[i]);
//        par_Aw2file [i+1] = par_Aw[i+1];
//        par_Bw2file [i+1] = par_Bw[i+1];
//        par_Cw2file [i+1] = par_Cw[i+1];
//        par_Dw2file [i+1] = par_Dw[i+1];
//        erpaw2file  [i+1] = erpaw[i+1];
//        erroffw2file[i+1] = erroffw[i+1];
//        erpcw2file  [i+1] = erpcw[i+1];
//        erpdw2file  [i+1] = erpdw[i+1];
//      }else if(par_Aw[i+1] <= min_alpha || par_Aw[i+1] >= max_alpha){
//        printf("Warning 2! Alpha PMT's %d-%d (west) = %f\n",i+1,i+2,par_Aw[i+1]);
//        par_Aw2file [i+1] = par_Aw[i];
//        par_Bw2file [i+1] = par_Bw[i];
//        par_Cw2file [i+1] = par_Cw[i];
//        par_Dw2file [i+1] = par_Dw[i];
//        erpaw2file  [i+1] = erpaw[i];
//        erroffw2file[i+1] = erroffw[i];
//        erpcw2file  [i+1] = erpcw[i];
//        erpdw2file  [i+1] = erpdw[i];
//      }

//      for(int k=0; k<3; k++) {
//        cove[k][i] = (cove[k][i]+cove[k][i+1])/2.0;
//        covw[k][i] = (covw[k][i]+covw[k][i+1])/2.0;
//      }
     
//      twe_calib << setw(13) << i
// 	       << setw(13) << par_Ae2file[i] << setw(13) << erpae2file[i]
// 	       << setw(13) << par_Be2file[i] << setw(13) << erroffe2file[i]
// 	       << setw(13) << par_Ce2file[i] << setw(13) << erpce2file[i]
// 	       << setw(13) << par_De2file[i] << setw(13) << erpde2file[i]
// 	       << setw(13) << csqe[i]
// 	       << setw(13) << cove[0][i]
// 	       << setw(13) << cove[1][i]
// 	       << setw(13) << cove[2][i] << endl;
     
//      tww_calib << setw(13) << i+216 
// 	       << setw(13) << par_Aw2file[i] << setw(13) << erpaw2file[i]
// 	       << setw(13) << par_Bw2file[i] << setw(13) << erroffw2file[i]
// 	       << setw(13) << par_Cw2file[i] << setw(13) << erpcw2file[i]
// 	       << setw(13) << par_Dw2file[i] << setw(13) << erpdw2file[i]
// 	       << setw(13) << csqw[i]
// 	       << setw(13) << covw[0][i]
// 	       << setw(13) << covw[1][i]
// 	       << setw(13) << covw[2][i] << endl;
     
//    } // loop over bars-1
   
//    twe_calib << setw(13) << 215
// 	     << setw(13) << par_Ae2file[215] << setw(13) << erpae2file[215]
// 	     << setw(13) << par_Be2file[215] << setw(13) << erroffe2file[215]
// 	     << setw(13) << par_Ce2file[215] << setw(13) << erpce2file[215]
// 	     << setw(13) << par_De2file[215] << setw(13) << erpde2file[215]
// 	     << setw(13) << csqe[215]
// 	     << setw(13) << cove[0][215]
// 	     << setw(13) << cove[1][215]
// 	     << setw(13) << cove[2][215] << endl;
   
//    tww_calib << setw(13) << 431 
// 	     << setw(13) << par_Aw2file[215] << setw(13) << erpaw2file[215]
// 	     << setw(13) << par_Bw2file[215] << setw(13) << erroffw2file[215]
// 	     << setw(13) << par_Cw2file[215] << setw(13) << erpcw2file[215]
// 	     << setw(13) << par_Dw2file[215] << setw(13) << erpdw2file[215]
// 	     << setw(13) << csqw[215]
// 	     << setw(13) << covw[0][215]
// 	     << setw(13) << covw[1][215]
// 	     << setw(13) << covw[2][215] << endl;
   
//    twe_calib.close(); 
//    tww_calib.close();
//    twe_alphas_diff.close(); 
//    tww_alphas_diff.close(); 
   
 } // if (flag==0)

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