// this is the new analysis code to read and analyze the data from 
// 5 f1TDC in slots 4 to 8 and 2 FADCs in slot 13 and 14 and one 
// CAEN ADC in slot 19
//
// f1TDC slot 1: IU-chambers x bottom input 1, 2, 3, 4
//       slot 2: IU-chambers y bottom input 1, 2, 3,       IU-chamber x top 4 
//       slot 3: IU-chambers x top 1, 2,                   IU-chambers y top 3, 4  
//       slot 4: IU-chambers y top 1, 
//       slot 5: FADC-chamber 1,                           Trigger Scintillators 3, 4
//


// include files */
#include <iostream>
#include <iomanip>
#include <stdio.h>
#include <stdlib.h>
#include <fstream>
using namespace std;

// root include files
#include "TFile.h"
#include "TH1.h"
#include "TH1F.h"
#include "TH2F.h"
#include "TH3F.h"
#include "TF1.h"
#include "TMath.h"
#include "TTree.h"
#include "TGraph.h"
#include "TCanvas.h"



#include <string.h>

#include "channelnumbers.h"

Int_t DEBUG1 = 0;


Float_t PI = 3.1415926;
Float_t PIh = 1.5707963;
Float_t strip_pitch = 5.0;  // distance between strips is 5mm (4mm width!)
Float_t strip_angle = PI/180.*15.0;
Float_t HWindow = 30.*strip_pitch/sin(strip_angle); // horizontal width of the cathode frame
Float_t YOffset = HWindow*tan(PI-strip_angle) ; // vertical width of the cathode frame
Float_t sd = strip_pitch;       // distance between two strips projected on vertical axis (wire axis)
Float_t CellRadius = 12.5 ; //units in mm

Float_t XCathOffSet = -47.5; // in units of strips
Float_t YCathOffSet = 0.5;   // in units of strips

Float_t getCathodeXpos(Float_t xpos, Float_t ypos){

 // return (xpos*sin(strip_angle)+ypos*cos(strip_angle));
 // return (xpos-ypos)/2./tan(strip_angle);
  Float_t b1, b2,x0;
  //b1 = YOffset  + (xpos-15.5)*sd;
  //b2 = (ypos-15.5)*sd;
  //x0 = (b1-b2)/2./tan(strip_angle);
  //b1 = xpos*(sin(strip_angle)+cos(strip_angle)*cos(strip_angle)/sin(strip_angle))*sd;
  //b2 = ypos/tan(strip_angle)*sd;
  //x0 = (b1-b2)/2./tan(PI-strip_angle);

  //b1 = xpos/cos(strip_angle)*sd;
  //b2 = ypos/cos(strip_angle)*sd;
  //x0 = (b2-b1)/2./tan(strip_angle);

  x0 = (xpos-ypos)*strip_pitch/2./sin(strip_angle);


  return x0;
}
Float_t getCathodeYpos(Float_t xpos, Float_t ypos){

  //return (xpos*cos(strip_angle)-ypos*sin(strip_angle));
  //return (xpos-ypos)/2.;
  Float_t b1, b2 ,y0;
  //b1 = YOffset - sd/2. + (xpos-15.)*sd;
  //b2 = -sd/2. + (ypos-15.)*sd;
  //y0 = (b1-b2)/2.+b2;
  //b1 = xpos*(sin(strip_angle)+cos(strip_angle)*cos(strip_angle)/sin(strip_angle))*sd;
  //b2 = ypos/tan(strip_angle)*sd;
  //y0 = (b1-b2)/2.+b2 - 200.;

  //b1 = xpos/cos(strip_angle)*sd;
  //b2 = ypos/cos(strip_angle)*sd;
  //y0 = (b1+b2)/2.;

  y0 = (xpos+ypos)*strip_pitch/2./cos(strip_angle);

  return y0;
}

Int_t GetMaxBin(Float_t *p,Int_t k){

  Float_t max = 0;
  Int_t ind = -9999;

  for (Int_t i=0;i<k;i++){
    if (p[i]>max){
      max = p[i];
      ind = i;
    }
  }
    
  return ind;

}

Float_t T2D_FDC[2][42];
void ReadT2Ddata_FDC(char* fnam){

  ifstream INF;
  INF.open(fnam,ios::in);
  if (!INF.is_open()){
    cout<<"Error open file "<<fnam<<endl;
    return;
  }

  Int_t j =0;
  while(!INF.eof()){
    INF >> T2D_FDC[0][j] >> T2D_FDC[1][j];
    //cout<< T2D_FDC[0][j] <<"  "<<  T2D_FDC[1][j]<<endl;
    j++;
  }

  INF.close(); 
}



Float_t R_Cell_FDC = 5.;


// time to distance transformation using analytic function from 
// fit to FDC data
Float_t time2distFDC(Float_t time){
  
  // time to distance calculation for FDC prototype Ar/CO2 mixture 60/40
  // used channel number 268 
  // F(x) = sqrt(t)*p0 + p1 + p2*t

  //P0 = 0.09425
  //P1 = 0.0251859
  //P2 = -0.000389234

  //P0 = 0.0996152
  //P1 = 0.210598
  //P2 = -0.000841775

  // from run 1336 Ar/CO2 is 40/60
  //P0 = 0.0732676
  //P1 = -0.284804
  //P2 = 0.000538546


  Float_t t = time;
  Float_t r = sqrt(t)*0.0732676 - 0.284804 + 0.000538546 * t;

  return r*R_Cell_FDC;

}


// time to distance transformation using interpolation with data from
// FDC prototype.
Float_t time2distFDC_Inter(Float_t time){

  if (time>T2D_FDC[0][42])
    return CellRadius;
  
  if (time<T2D_FDC[0][0])
    return 99999.;

  Int_t k = 21;
  Int_t j = 21;
  Int_t notfound = 1;
  while (notfound){
    if (time>T2D_FDC[0][j]){
      k = k/2;
      j = j + k;
    } else{
      j = j - k;
    }
    if (k<1)
      notfound = 0;
  }

  //cout<<"FDC t2d: "<<j<<"  "<<time<<"  "<<T2D_FDC[0][j]<<"  "<<T2D_FDC[1][j]<<endl;
  
  Float_t x = (time-T2D_FDC[0][j])/(T2D_FDC[0][j+1]-T2D_FDC[0][j]);
  Float_t y = T2D_FDC[1][j] + (T2D_FDC[1][j+1] - T2D_FDC[1][j])*x;
  Float_t dist = y*CellRadius;
    
  //cout<<x<<"  "<<y<<"  "<<dist<<endl;

  return dist;

}



// code for drift time to distance conversion of IU aluminum wire chambers 



Float_t time2distIU(Float_t time){

  // time to distance for IU chamberse 90/10 Ar/CO2 at 2200V
  // used wire channel 20 for calibration of 7th order polynomial fit
  // with time zero at 412.45. These numbers are based on run 138
  // P0 = 2516.82
  // P1 = -32.2102
  // P2 = 0.175796
  // P3 = -0.000530658
  // P4 = 9.57064e-07
  // P5 = -1.03134e-09
  // P6 = 6.14844e-13
  // P7 = -1.5643e-16

  const Float_t Toffset = 412.45;

  Float_t t = time + Toffset;
  
  if (t<Toffset)
    return 88888.;
    
  if (t>701.)
    return 77777.;

  Float_t dist = 2516.82 -32.2102*t + 0.175796*t*t - 0.000530658*t*t*t + 9.57064e-07*t*t*t*t
    - 1.03134e-09*t*t*t*t*t + 6.14844e-13*t*t*t*t*t*t - 1.5643e-16*t*t*t*t*t*t*t;

  dist *= CellRadius;

  return dist;
}


Int_t ReadDriftTimeZeros(char * fnam, Double_t * DTZ){
  
  ifstream INF;
  INF.open(fnam, ios::in);
  if (!INF.is_open()){
    cout<<"Error open DriftTime zero file: "<<fnam<<endl;
  }
  Int_t k;
  Int_t cnt=0;
  while (!INF.eof()){
    INF>>k;
    INF>>DTZ[k];
    cnt++;
  }
  INF.close();
  return cnt;
}


// 
struct Hit{
  Int_t plane;
  Int_t wire;
  Float_t drifttime;
  Float_t driftdist;
  Float_t x0;
  Float_t x;
  Float_t z;
  Float_t dx;
};

struct Hit HitListX[4][20]; 
struct Hit HitListY[4][20]; 
struct Hit FDC_X[16];

struct fitres{
  Float_t val;
  Float_t dval;
};



void iu_fit(struct Hit *list, Int_t l, fitres *a, fitres *b, Float_t *chi2){
  
  unsigned int i;
  Float_t wt,t,sxoss,sx=0.0,sy=0.0,st2=0.0,ss=0.0;
  
  //b->val=0.0;
  // Accumulate sums
  //for (i=0;i<l;i++){
  //  wt=1./list[i].dx/list[i].dx;
  //  ss+=wt;
  //  sx+=list[i].z*wt;
  //  sy+=list[i].x*wt;
  // }
  //sxoss=sx/ss;
  //for (i=0;i<l;i++){
  //  t=(list[i].z-sxoss)/list[i].dx;
  // st2+=t*t;
  // b->val+=t*list[i].x/list[i].dx;
  //}
  //b->val/=st2;
  // Solve for a,b,siga,sigb
  //a->val=(sy-sx*(b->val))/ss;
  //a->dval=sqrt((1.+sx*sx/(ss*st2))/ss);
  //b->dval=sqrt(1.0/st2);
  
  // Compute the chi2 for the fit
  //*chi2=0.0;
  //for (i=0;i<l;i++){
  //  
  //  *chi2+=((list[i].x-(a->val)-(b->val)*list[i].z)/list[i].dx) * ((list[i].x-(a->val)-(b->val)*
  //								    list[i].z)/list[i].dx);
  //}
  
  Float_t Sw = 0.;
  Float_t Sx = 0.;
  Float_t Sz = 0.;
  Float_t Sxx = 0.;
  Float_t Sxz = 0.;
  
  for (i=0;i<l;i++){
    Float_t w = 1./list[i].dx/list[i].dx;
    Sw += w;
    Sx += list[i].z*w;
    Sz += list[i].x*w;
    Sxx += list[i].z*list[i].z*w;
    Sxz += list[i].z*list[i].x*w;    
  }
  
  Float_t delta = Sw*Sxx-Sx*Sx;
  a->val = (Sxx*Sz-Sx*Sxz)/delta;
  b->val = (Sw*Sxz-Sx*Sz)/delta;
  
  *chi2 = 0;
  for (i=0;i<l;i++){
    
    *chi2+=( (list[i].x-(a->val)-(b->val)*list[i].z) / list[i].dx ) * 
      ( (list[i].x-(a->val)-(b->val)*list[i].z) / list[i].dx );
  }


}

Float_t fit_trackX(Int_t Npoints, Int_t* FDC_N, fitres *interseptx, fitres *slopex){

  struct Hit HLX[10];
  Float_t oldchi2 = 99999.;
  Float_t chi2 = 99999;
  fitres a,b;

  for (Int_t n1=-1;n1<2;n1+=2){
    HLX[0].x = HitListX[0][0].x + ((Float_t)n1)*HitListX[0][0].driftdist;
    HLX[0].z = HitListX[0][0].z;
    HLX[0].dx =  HitListX[0][0].dx;
    
    for (Int_t n2=-1;n2<2;n2+=2){
      HLX[1].x =   HitListX[1][0].x + ((Float_t)n2)*HitListX[1][0].driftdist;
      HLX[1].z =   HitListX[1][0].z;
      HLX[1].dx =  HitListX[1][0].dx;
      
      for (Int_t n3=-1;n3<2;n3+=2){
	HLX[2].x =   HitListX[2][0].x + ((Float_t)n3)*HitListX[2][0].driftdist;
	HLX[2].z =   HitListX[2][0].z;
	HLX[2].dx =  HitListX[2][0].dx;
	
	for (Int_t n4=-1;n4<2;n4+=2){
	  HLX[3].x =   HitListX[3][0].x + ((Float_t)n4)*HitListX[3][0].driftdist;
	  HLX[3].z =   HitListX[3][0].z;
	  HLX[3].dx =  HitListX[3][0].dx;
	  
	  if (Npoints==5){

	    Int_t cntx = 0; // get the wire in the second plane
	    while (FDC_X[cntx].wire>271){
	      cntx++;
	    }
	    
	    for  (Int_t n5=-1;n5<2;n5+=2){

	      HLX[4].x  = FDC_X[cntx].x  + ((Float_t)n5)*FDC_X[cntx].driftdist;
	      HLX[4].z  = 684.;
	      HLX[4].dx = 0.5;
	      iu_fit(HLX, 5, &a, &b, &chi2);
	      
	      if (chi2<oldchi2){
		oldchi2 = chi2;
		interseptx->val = a.val;
		interseptx->dval = a.dval;
		slopex->val = b.val;
		slopex->dval = b.dval;
		*FDC_N = n5;
	      }
	    }
	  } else {

	    iu_fit(HLX, 4, &a, &b, &chi2);
	      
	    if (chi2<oldchi2){
	      oldchi2 = chi2;
	      interseptx->val = a.val;
	      interseptx->dval = a.dval;
	      slopex->val = b.val;
	      slopex->dval = b.dval;
	    }
	    
	  }
	  
	}
	
      }
      
    }
  }   
  
  return oldchi2;

}

Float_t fit_trackY(Int_t Npoints, Float_t ypos, fitres *intersepty, fitres *slopey){

  struct Hit HLY[10];
  Float_t oldchi2 = 99999.;
  Float_t chi2 = 99999;
  fitres a,b;
  
  for (Int_t n1=-1;n1<2;n1+=2){
    HLY[0].x = HitListY[0][0].x + ((Float_t)n1)*HitListY[0][0].driftdist;
    HLY[0].z = HitListY[0][0].z;
    HLY[0].dx =  HitListY[0][0].dx;
    
    for (Int_t n2=-1;n2<2;n2+=2){
      HLY[1].x =   HitListY[1][0].x + ((Float_t)n2)*HitListY[1][0].driftdist;
      HLY[1].z =   HitListY[1][0].z;
      HLY[1].dx =  HitListY[1][0].dx;
      
      for (Int_t n3=-1;n3<2;n3+=2){
	HLY[2].x =   HitListY[2][0].x + ((Float_t)n3)*HitListY[2][0].driftdist;
	HLY[2].z =   HitListY[2][0].z;
	HLY[2].dx =  HitListY[2][0].dx;
	
	for (Int_t n4=-1;n4<2;n4+=2){
	  HLY[3].x =   HitListY[3][0].x + ((Float_t)n4)*HitListY[3][0].driftdist;
	  HLY[3].z =   HitListY[3][0].z;
	  HLY[3].dx =  HitListY[3][0].dx;
	  
	  if (Npoints==5){

	    Int_t cntx = 0; // get the wire in the second plane
	    while (FDC_X[cntx].wire>271){
	      cntx++;
	    }
	    
	    HLY[4].x = ypos + 206.8;
	    HLY[4].z = 684;
	    HLY[4].dx = 0.5;
	    iu_fit(HLY, 5, &a, &b, &chi2);
	  } else {
	    iu_fit(HLY, 4, &a, &b, &chi2);
	  }
	  if (chi2<oldchi2){
	    oldchi2 = chi2;
	    intersepty->val = a.val;
	    intersepty->dval = a.dval;
	    slopey->val = b.val;
	    slopey->dval = b.dval;
	  }
	  
	}
	
      }
      
    }
    
  }
  
  return oldchi2;
  
}


// 
// 
// Main analysis code looping over events
//
Int_t analyze1(Int_t runnumber, char* dfil, TH1F** fadc, TH1F** fadcNP, 
	       TH1F** histT, TH2F **hist2d, TH1F **qadc, 
	       TH1F ** DriftTimes, Int_t Flag1, Int_t FLAG2) {
  
  //////////////////////////////////////////////////////////
  //   This file has been automatically generated 
  //     (Wed Jul 22 17:53:42 2009 by ROOT version5.18/00)
  //   from TTree fdctest/FDC response
  //   found on file: fdcTest_000012.root
  //////////////////////////////////////////////////////////
  
  
  //Reset ROOT and connect tree file
  
  //TFile *f = (TFile*)gROOT->GetListOfFiles()->FindObject(dfil);
  //if (!f) {
  TFile *f = new TFile(dfil);
  //}
  TTree *fdctest = (TTree*)gDirectory->Get("fdctest");
  
  //Declaration of leaves types
  Int_t           nevent;
  Int_t           ntdc;
  Float_t         tdc[500];
  Int_t           tdc_chan[500];
  Int_t           nFADC;
  Float_t         FADC_SUM[32];
  Float_t         FADC_SUMSUB[32];
  Int_t           nCAENADC;
  Float_t         CAEN_ADC[64];
  Float_t         CAEN_OVFL[64];
  
  Float_t         QDC_ped[2][64];
  Float_t         QDC_sig[2][64];
  Float_t         QDC_lim[64];
  Float_t         FADC_ped[2][32];
  Float_t         FADC_sig[2][32];
  Float_t         FADC_lim[32];

  Float_t         SigCut = 3.5;  // cut above pedestal for strip signals

  Float_t         Chi2X;
  Float_t         Chi2Y;

  
  // Set branch addresses.
  fdctest->SetBranchAddress("nevent",&nevent);
  fdctest->SetBranchAddress("ntdc",&ntdc);
  fdctest->SetBranchAddress("tdc",tdc);
  fdctest->SetBranchAddress("tdc_chan",tdc_chan);
  fdctest->SetBranchAddress("nFADC",&nFADC);
  fdctest->SetBranchAddress("FADC_SUM",FADC_SUM);
  fdctest->SetBranchAddress("FADC_SUMSUB",FADC_SUMSUB);
  fdctest->SetBranchAddress("nCAENADC",&nCAENADC);
  fdctest->SetBranchAddress("CAEN_ADC",&CAEN_ADC);
  fdctest->SetBranchAddress("CAEN_OVFL",&CAEN_OVFL);
  
  //     This is the loop skeleton
  //       To read only selected branches, Insert statements like:
  // fdctest->SetBranchStatus("*",0);  // disable all branches
  // TTreePlayer->SetBranchStatus("branchname",1);  // activate branchname
  
  Long64_t nentries = fdctest->GetEntries();

  cout<<"Number of Entries: "<<nentries<<endl;

  Float_t Timing[2][16];
  Float_t deltaT[16];
  Float_t sumT[16];
  Float_t Distance[16];
  Float_t Time2Dist = 16.0; //Time to distance 16cm/ns
  Float_t timeOfFlight[16]; 
  Float_t TDC_Resolution = 0.1283; // TDC resolution in ns => bin size
  Float_t T_Ref = 0; // TDC time reference take always the last channel
  
  Float_t wireTime[300];
  Int_t wireChan[300];
  Int_t wireNum;
  Int_t wireHit[300];
  Int_t FDCHit[80];
  Int_t FDCHits;
  Int_t FDC_Hits[2];
  Int_t FHits;
  Int_t FDCHits_Chan[300];
  Float_t FDCHits_Time[300];
  Float_t Planes[2][48];
  Float_t Xaxis[48];  
  Int_t PlaneCnt[2];

  Float_t WTOffset = 500.;
  
  Int_t IUHits = 0;

  Int_t FDC_Hit_Counter[2][16];


  for (Int_t k=0;k<48;k++)
    Xaxis[k] = k;

  // read pedestal data for QADC and FADC

  ifstream INF; 
  char fn1[128];
  sprintf(fn1,"qadc_pedestals_run%d.dat",runnumber);

  //if (!access(fn1,F_OK)){ // test for file existance
  //  cout << "Pedestal file for curernt run not found try previous run"<<endl;
  //  Int_t nr = runnumber-1;
  //  sprintf(fn1,"qadc_pedestals_run%d.dat",nr);
  // }

  INF.open(fn1,ios::in );
  if (INF.is_open()){
    Int_t l;
    Float_t dummy;
    Int_t Nchan =0;
    INF>>Nchan;
    for (Int_t k=0;k<Nchan;k++){
      INF>>l>>QDC_ped[0][k]>>QDC_ped[1][k]>>QDC_sig[0][k]>>QDC_sig[1][k]>>dummy;
      //cout<<l<<" "<<QDC_ped[0][k]<<" "<<QDC_ped[1][k]<<" "<<
      // QDC_sig[0][k]<<" "<<QDC_sig[1][k]<<endl;
      QDC_lim[k] = QDC_sig[0][k]*SigCut;
      //cout<<k<<"  "<<QDC_lim[k]<<endl;
    }
    INF.close();
    cout<<"QDC pedestal file "<<fn1<<" found and read in!"<<endl;
    Flag1=1;
  } else{
    cout<<"No QDC pedestal file "<<fn1<<" found!"<<endl;
    cout<<"-> set Flag1 to zero"<<endl;
    Flag1 = 0;
  }
 
 
  //sprintf(fn1,"fadc_pedestals_run%d.dat",runnumber);
  //INF.open(fn1,ios::in );
  //if (INF.is_open()){
  //  for (Int_t k=0;k<32;k++){
  //    INF>>l>>FADC_ped[0][k]>>FADC_ped[1][k]>>FADC_sig[0][k]>>FADC_sig[1][k]>>dummy;
  //    FADC_lim[k] = FADC_ped[0][k] - FADC_sig[0][k]*SigCut;
  //  }
  //  INF.close();
  //} else {
  //  cout<<"No FADC pedestal file "<<fn1<<" found!"<<endl;
  //}

  // Read data for time to distance conversion


  //read gain calibration data for cathode strips

  ifstream INF1;
  sprintf(fn1,"cathode_gain_calib.txt");
  INF1.open(fn1,ios::in );

  Float_t Gains[64];
  Float_t Pedestals[64];
  Int_t chan1=0;
  Double_t ped;
  Int_t k1 = 0;
  if (INF1.is_open()){
    while (k1<64){
      INF1 >> chan1 >> Pedestals[k1] >> Gains[k1];
      //cout<<chan1<<" "<<Pedestals[k1]<<" "<<Gains[k1]<<endl;
      k1++;
    }
    INF1.close();
    cout<<"Gain calibration file "<<fn1<<" found and read in!"<<endl;
  } else {


    memset(Pedestals,0,256);
    //memset(Gains,0,256);
      
    for (k1=0;k1<64;k1++){
      //Pedestals[k1] = 0.0;
      Gains[k1] = 1.0;
    }
    Gains[53]=0.1;
    Gains[60]=0.1;
    Gains[15]=0.3;
      cout<<"NO Gain calibration file "<<fn1<<" found! All gains set to 1! " <<endl;
  }


  // Read data for Drift time zero
  Double_t DriftTimeZero[320];
  char fnamTZ[128] = "drifttimezero.txt";
  Int_t nZeros = ReadDriftTimeZeros(fnamTZ, DriftTimeZero);
  if (nZeros){
    cout<<"Drift time Zeros read in"<<endl;
  }
 
  // Read time 2 distance file for FDC
  // fitting the data with an analytic function does not a good job
  // so we do interpolation.
  char t2d_fdc_file[128] = "time_to_dist_FDCprototype.dat";
  ReadT2Ddata_FDC(t2d_fdc_file);
  cout<<"Time 2 distance file for FDC read in"<<endl;

 
  // initialize FDC_Hit_Counter 

  memset(&FDC_Hit_Counter[0][0],0,128);
  //for (Int_t k=0;k<16;k++){
  //  FDC_Hit_Counter[0][k] = 0;
  //  FDC_Hit_Counter[1][k] = 0;
  // }



  // calculate the linear functions that represent 
  // the FDC wires in the IU-chamber coordinate system

  Double_t slopeW  = 1.7316;    // slope of the FDC wire in the IU chamber coordinate system
  Double_t B3      = -39.3; // constant of 3rd FDC wire for f(x) = slopeW*x + B3
  Double_t alphaFDC= atan(slopeW);
  
  Double_t B[16];
  
  B[2] = B3;
  B[1] = B[2] - 10./cos(alphaFDC);
  B[0] = B[1] - 10./cos(alphaFDC);
  
  for (Int_t i=3;i<16;i++){
    
    B[i] = B[i-1] + 10/cos(alphaFDC);
    
  }

  
  TCanvas *c1 = new TCanvas("c1","A Simple Graph Example",200,10,700,500);
  
  c1->SetFillColor(42);
  c1->SetGrid();

  Int_t imod = 10000;

  Float_t TrackHitCounter[6] = {0., 0., 0., 0., 0., 0.};

  //                             !!!THE BIG LOOP!!!
  // loop over data
  // 
  Long64_t nbytes = 0;
  for (Long64_t i=0; i<nentries;i++) {
  //for (Long64_t i=0; i<5000000;i++) {
    
    //if (i<16961)
    //  continue;

    nbytes += fdctest->GetEntry(i);

    
    //cout<<"Event "<<i<<"    number of tdc hits" <<ntdc<<endl;

    if(i%imod==0){
      cout << i  << " events done."<<endl;
      //	   << " events" << endl;
      //      if (evtCount==100000) imod=100000;
    }
    
    
    if (ntdc<100){ 
   
      //memset(&Timing[0][0],0,128);
   
      for(Int_t j=0;j<16;j++){
	Timing[0][j] = 0.;
	Timing[1][j] = 9999999.;
	Distance[j]=9999999.;
	sumT[j]=999999.;
      }

      wireNum = 0;

      memset(wireHit,0,1200);
      memset(wireChan,0,1200);
      for(Int_t j=0;j<300;j++){
	//wireHit[j] =  0;
	//wireChan[j] = 0;
	wireTime[j] = 99999.;
      }
      //cout<<"Event: "<<i<<"   TDC hits: "<<ntdc<<endl;
      IUHits = 0;
      FDCHits = 0;
      FHits = 0;
      FDC_Hits[0] = 0;
      FDC_Hits[1] = 0;

      memset(FDCHit,0,320);
      //for (Int_t k=0;k<80;k++)
      //FDCHit[k] = 0;

      for (Int_t j=0;j<ntdc;j++){  // loop over all TDC hits
	
	Int_t lr = 0;
	Int_t chan = tdc_chan[j];
	Int_t channel = 0;

	// TDC channels from 288 are trigger scintillators
	if (chan>287){ 
	  channel = chan-288;
	  if (channel>15) {
	    lr = 1;
	    channel -=16;
	  }
	  
	  Timing[lr][channel] = (Float_t)tdc[j]*TDC_Resolution;
	  
	  //cout << j<<"  "<< tdc_chan[j]<<"  :"<<tdc[j]<<" "<<lr<<" "<<chan<<endl;
	} else { // all the rest are the wires chambers FDC and IU

	  // Sort wire channel numbers such that the IU chambers are
          // first and the FDC prototype last.
	  // IU chambers run from 0 up to 175(at most);
	  // FDC wire chamber runs from 176 to 195(at most)
	  // Also only take the first TDC hit per channel and ignore any secondary hit

	  if ((wireHit[chan]==0)){
	    wireChan[wireNum] = chan;
	    wireTime[wireNum++] = (Float_t)tdc[j]*TDC_Resolution;
	    wireHit[chan] += 1;
	    if (chan>207){
	      if (chan>255){
		Int_t wnum = chan-256;
		FDCHit[wnum] = 1;
	      }
	    } else
	      IUHits++;
	  }
	}
      }
      
      // Reference time is the trigger signal fed to channel 16 of the second
      // Camac discriminator  
      T_Ref = Timing[1][15];
      
     
      // calculate TDC time difference and position of trigger scintillators
      UInt_t TOP;
      UInt_t BOT;
      TOP = 0;
      BOT = 0;
      Int_t Tcnt = 0;
      Int_t Bcnt = 0;
      Int_t Top, Bot;

      Float_t DistT = 0.0;
      Float_t DistB = 0.0;

      for(Int_t j=0;j<8;j++){
	deltaT[j] = Timing[0][j] - Timing[1][j] ;
	Distance[j] = deltaT[j]*Time2Dist/2.;
	sumT[j] = Timing[0][j] + Timing[1][j] ;
	if (TMath::Abs(Distance[j])<120){
	  BOT = BOT | (0x01<<j);
	  Bot = j;
	  Bcnt++;
	  DistB = Distance[j];

	}
	histT[j]->Fill(deltaT[j],1);  
	histT[j+16]->Fill(Distance[j],1);     

	hist2d[15]->Fill(Distance[j],(Float_t)j);

      }
      for(Int_t j=8;j<16;j++){
	deltaT[j] = Timing[0][j] - Timing[1][j] ;
	Distance[j] = deltaT[j]*Time2Dist/2.;
	sumT[j] = Timing[0][j] + Timing[1][j] ;
	if (TMath::Abs(Distance[j])<120){
	  //cout<<"intop!"<<(0x01<<(j-8))<<endl;
	  TOP = TOP | (0x01<<(j-8));
	  Top = j;
	  Tcnt++;	
	  DistT = Distance[j];
	}
	histT[j]->Fill(deltaT[j],1);  
	histT[j+16]->Fill(Distance[j],1);     

	hist2d[16]->Fill(Distance[j],(Float_t)(j-8));
      }
      
      // only look at events with one hit per scintillator plane
      if ((Tcnt>1)||(Bcnt>1)){
      	continue;
      }
      if ((Tcnt==0)||(Bcnt==0)){
      	continue;
      }

      if (((Top-8<2)&&(Top-8>3)) ||
      	  ((Bot<2)&&(Bot>3))){
	continue;
      }
      
      if (TMath::Abs(DistB)>50.)
	continue;

      if (TMath::Abs(DistT)>50.)
	continue;


      histT[32]->Fill((sumT[Top]-sumT[Bot])/2.,1);
      //cout<<"TOP AND BOT:   "<<TOP<<"   "<<BOT<<endl;
      hist2d[0]->Fill((Float_t)Bot+0.5, (Float_t)Top-8.0+0.5, 1);
      
      Int_t Flag2 = FLAG2;
      Flag2 = 1;
      Int_t HitCounterX[4];
      Int_t HitCounterY[4];
 
      for (Int_t j=0;j<4;j++){
	HitCounterX[j] = 0;
	HitCounterY[j] = 0;
      }

      // loop over all wire chamber hits
 

      Float_t ChamberOff = 0.;

      
      for (Int_t k=0;k<wireNum;k++){
	if (T_Ref>wireTime[k]){
	  //cout<<i<<"   T_Ref is bigger than wire time!!!\n";
	  wireTime[k] -= T_Ref;
	  wireTime[k] += WTOffset;
	  //wireTime[k] = 0.; 
	} else {
	  wireTime[k] -= T_Ref;
	  wireTime[k] += WTOffset;
	}
	//cout<<"wire number: "<<wireChan[k]<<"    Drift time= "<<wireTime[k]<<endl;
	//flush(cout);
	DriftTimes[wireChan[k]]->Fill(wireTime[k],1);
	hist2d[3]->Fill(wireTime[k],(Float_t)wireChan[k]-0.5,1);

	//if (wireChan[k]>255){
	//  FDCHits_Chan[FHits] = wireChan[k];
	//  FDCHits_Time[FHits] = wireTime[k];
	// FHits++;
	//}
	
	Int_t Chan = wireChan[k];
	
	if ( (Flag2)){
	  
	  if (Chan<208) { // these are the IU ALuminum chambers
	    Float_t DriftDist = 99999.;
	    Float_t dtime = wireTime[k]-DriftTimeZero[wireChan[k]];
	    if (dtime>=0.)
	      DriftDist = time2distIU(dtime);
	    
	    Int_t channel = IU_Channels[Chan];
	    Int_t plane   = IU_Planes[Chan]-1;
	    ChamberOff = (Float_t)IU_ChamberOff[Chan];
	    
	    Int_t xplane = 0;
	    Int_t yplane = 0;

	    if ( (dtime>=0.) && 
		 (DriftDist<=CellRadius) ){

	      if ((Chan<64)||
		  ((Chan>111)&&(Chan<160))){
		xplane = IU_Planes[Chan];
	      } else{
		yplane = IU_Planes[Chan];
	      }
	      

	      // The ChamberOff is taking care of the fact that the chamber aluminum frame
	      // is about 1.5mm thick resulting in an additional 3mm gap between the last
	      // cell of the previous chamber and the first cell of the next chamber.

	      // Distance between the top of the IU-x chamber below and above is 1540.5mm
              // This is also the distance between plane 1y and plane 3y or plane 2y and plane 4y
              // or plane 1x and plane 3x or plane 2x and plane 4x

	      //          
	      //          --------------------- plane 4X
	      //          --------------------- plane 3X
	      //          --------------------- plane 4Y   -> set z = 0
	      //          --------------------- plane 3Y
	      //          +++++++++++++++++++++ TOP scintillators
	      //
	      //
	      //
	      //
	      //          --------------------- plane 2X
	      //          --------------------- plane 1X
	      //          --------------------- plane 2Y   -> z = 1540.5 = DeltaDist
	      //          --------------------- plane 1Y
              //          +++++++++++++++++++++ BOTTOM scintillators
              //
              //        _________________________________________
	      //   _____| 14    12    10    8     6     4     2 |____            -> IU-Chamber has 15 channels
	      //   |__15_____13____11____9_____7_____5_____3_____1__|            -> view from the readout side
              //


	      // the z direction is vertial and z=0 is defined to be the 4th y plane

	      Float_t IU_pos_err = 1.5;     // IU-chamber resolution
	      Float_t DeltaDist = 1540.5;   
	      Float_t XOffset = 103.0;      // Top x-chambers (3,4) shifted w.r.t. bottom chambers
	      Float_t ChamberWall = 2.5;    // gap between two IU-chambers

	      if (xplane){
		//cout<<"xplane: "<<xplane<<endl;
		HitListX[plane][HitCounterX[plane]].plane = xplane;
		HitListX[plane][HitCounterX[plane]].wire = channel;
		HitListX[plane][HitCounterX[plane]].drifttime = dtime;
		HitListX[plane][HitCounterX[plane]].driftdist = DriftDist;
		Float_t x = ((Float_t)channel) * 2.* CellRadius + CellRadius + ChamberOff*ChamberWall;
		Float_t z = DeltaDist - 2.* CellRadius;
		if (xplane==2){
		  x -= CellRadius;
		  z -= 2.*CellRadius ;
		}
		if (xplane==3) {
		  x = x + XOffset - CellRadius;
		  z = -2. * CellRadius;
		}
		if (xplane==4){
		  x += XOffset; 
		  z = -4. *  CellRadius;
		}
		HitListX[plane][HitCounterX[plane]].x = x;
		HitListX[plane][HitCounterX[plane]].z = z;
		HitListX[plane][HitCounterX[plane]].dx = IU_pos_err;
		if (HitCounterX[plane]<19){
		  HitCounterX[plane]++;
		}
	      }

	      Float_t YOffset = 67.0; // bottom yplanes are shifted by 67mm w.r.t. top yplanes

	      if (yplane){
		//cout <<"Event "<<i;
		//cout<<"   yplane: "<<yplane<<"  "<<channel<<endl;
		HitListY[plane][HitCounterY[plane]].plane = yplane;
		HitListY[plane][HitCounterY[plane]].wire = channel;
		HitListY[plane][HitCounterY[plane]].drifttime = dtime;
		HitListY[plane][HitCounterY[plane]].driftdist = DriftDist;
		Float_t x = (Float_t)channel * 2.* CellRadius + YOffset + ChamberOff*ChamberWall;
		Float_t z = DeltaDist + 2.* CellRadius;


		if (yplane==2){
		  x += CellRadius ; 
		  z -= 2.0 * CellRadius;
		}
		if (yplane==3) {
		  x -= YOffset ;
		  z = 2.0 * CellRadius;
		}
		if (yplane==4){
		  x = x - YOffset + CellRadius ;
		  z = 0.;
		}
		HitListY[plane][HitCounterY[plane]].x = x;
		HitListY[plane][HitCounterY[plane]].z = z;
		HitListY[plane][HitCounterY[plane]].dx = IU_pos_err;
		if (HitCounterY[plane]<19){
		  HitCounterY[plane]++;
		}
	      }
	    }

	  } else { 

	    // now the FDC wires
	    
	    if (Chan>255){
	      Int_t wnum = Chan-256; 
	      Float_t drifttime = wireTime[k]-DriftTimeZero[wireChan[k]];
	      FDC_X[FDCHits].plane = 0;
	      FDC_X[FDCHits].wire = wnum; 
	      if ((drifttime<=164.5)&& (drifttime>=0.)){
		FDC_X[FDCHits].drifttime = drifttime;
		FDC_X[FDCHits].driftdist = time2distFDC_Inter(drifttime);
		FDC_X[FDCHits].x = -1.*(wnum*10.- 130.) + 200.;
		//cout << FDCHits<<": drifttime/driftdist   " << drifttime << "  " << FDC_X[FDCHits].driftdist << endl;
	      }else{
		FDC_X[FDCHits].drifttime = 99999.;
		FDC_X[FDCHits].driftdist = 99999.;
	      }
	      FDCHits_Chan[FDCHits] = wireChan[k];
	      FDCHits_Time[FDCHits] = wireTime[k];
	      FDCHits++;

	      if (Chan<272){
		FDC_Hits[0]++;
	      } else {
		FDC_Hits[1]++;
	      }

	    }
	    
	    
	  }
	  
	}// Flag2
	
      }// end loop over wire hits
      

      if (FDC_Hits[0]==1){

	//fcout<<"two hits"<<endl;
	Int_t cntx = 0;
	while (( FDCHits_Chan[cntx]>271) && (cntx<FDCHits)){
	  cntx++;
	}

	DriftTimes[FDCHits_Chan[cntx]+44]->Fill(FDCHits_Time[cntx],1);
	hist2d[3]->Fill(FDCHits_Time[cntx],(Float_t)FDCHits_Chan[cntx]+44.-0.5,1);

      }

      // do hit accounting to test efficiency

      if (FDC_Hits[1]==1){

	//fcout<<"two hits"<<endl;
	Int_t cntx = 0;
	while (( FDCHits_Chan[cntx]<272) && (cntx<FDCHits)){
	  cntx++;
	}

	DriftTimes[FDCHits_Chan[cntx]+44]->Fill(FDCHits_Time[cntx],1);
	hist2d[3]->Fill(FDCHits_Time[cntx],(Float_t)FDCHits_Chan[cntx]+44.-0.5,1);
      }

      if (FDC_Hits[0]==1){
	Int_t cntx = 0;
	while (( FDCHits_Chan[cntx]>271) && (cntx<FDCHits)){
	  cntx++;
	}
	Int_t wire1 = FDCHits_Chan[cntx] - 256;
	FDC_Hit_Counter[0][wire1] +=1;
	for (Int_t j=0;j<FDCHits;j++){
	  if ( FDCHits_Chan[j]>271){
	    Int_t wire0 = FDCHits_Chan[j] - 272;
	    if ((wire0==wire1) || (wire0-1==wire1) || (wire0+1==wire1))
	      FDC_Hit_Counter[1][wire1] +=1;
	  }
	}	
      }
      

      Float_t FDCDeltaT = 999.;

      if  ((FDC_Hits[0]==1) &&
	   (FDC_Hits[1]==1)){

	//cout<<"TWO HITS"<<endl;

	Float_t ch1 = (Float_t)(FDCHits_Chan[0]-256)+0.5;
	Float_t ch2 = (Float_t)(FDCHits_Chan[1]-272)+0.5;

	Float_t t1 = FDCHits_Time[0];
	Float_t t2 = FDCHits_Time[1];


	if (FDCHits_Chan[1]<FDCHits_Chan[0]){
	  ch1 = (Float_t)(FDCHits_Chan[1]-256)+0.5;
	  ch2 = (Float_t)(FDCHits_Chan[0]-272)+0.5;	  

	  Float_t t1 = FDCHits_Time[1];
	  Float_t t2 = FDCHits_Time[0];	  
	}

	hist2d[35]->Fill(ch2, ch1);
	
	FDCDeltaT = t1-t2;


      }

      struct fitres a, b, interseptx, slopex, intersepty, slopey;
      Float_t chi2 = 99999.;
      Int_t mulx = HitCounterX[0]*HitCounterX[1]*HitCounterX[2]*HitCounterX[3];

      if (DEBUG1) {
	cout<<"Event "<<i<<endl;
	cout<<"X: ";
	for (Int_t kkk=0;kkk<4;kkk++)
	  cout<<  HitCounterX[kkk]<<"  ";
	cout<<endl;
	cout<<"Y: ";
	for (Int_t kkk=0;kkk<4;kkk++)
	  cout<<  HitCounterY[kkk]<<"  ";
	cout<<endl;
      }

      Float_t oldchi2 = 99999.;
      struct Hit HLX[10];
      struct Hit HLY[10];
      Int_t HLXcnt;
      Int_t HLYcnt;
      HLXcnt = HLYcnt = 0;
      
      //mulx *= FDCHits;
      Int_t FDC_N = 0;

      if (mulx==1) {

	oldchi2 = fit_trackX(4, &FDC_N, &interseptx, &slopex); // fit without FDC wire to be unbiased.

	//if  ((FDC_Hits[0]==1) &&
	//     (FDC_Hits[1]==1)){
	//
	//  oldchi2 = fit_trackX(5, &FDC_N, &interseptx, &slopex); // fit with FDC second wire
	//}

	//cout<<"f(x) = "<<slopex.val<<"x + "<<interseptx.val<<endl;
	histT[40]->Fill(oldchi2);
	histT[42]->Fill(slopex.val);
	if (TMath::Abs(slopex.val)<0.1)
	  histT[44]->Fill(oldchi2);
	
	Chi2X = oldchi2;
      }
      

      oldchi2 = 99999.;
      Int_t muly = HitCounterY[0]*HitCounterY[1]*HitCounterY[2]*HitCounterY[3];
      Float_t dummy =0.;
      if (muly==1){

	oldchi2 =  fit_trackY(4, dummy , &intersepty, &slopey);
	
	//cout<<"f(y) = "<<slopey.val<<"y + "<<intersepty.val<<endl;
	histT[41]->Fill(oldchi2);
	histT[43]->Fill(slopey.val);
	if (TMath::Abs(slopey.val)<0.1)
	  histT[45]->Fill(oldchi2);
	Chi2Y = oldchi2;
      }

      Float_t IU_x = 0.;
      Float_t IU_y = 0.;
      

      // FDC wire plane z-position between the IU chambers
      //Float_t WirePlanePos = 771.5;
      Float_t WirePlanePos = 781.5;
      if ((muly==1)&&(mulx==1)){

	//cout <<"Event "<<i<<endl;

	IU_x = slopex.val*WirePlanePos + interseptx.val;
	IU_y = slopey.val*WirePlanePos + intersepty.val;
	//IU_x = slopex.val*1178.925 + interseptx.val;
	//IU_y = slopey.val*1178.925 + intersepty.val;
	
	hist2d[5]->Fill(IU_x, IU_y,1.);         
	hist2d[18]->Fill(slopey.val, slopex.val);


	if ((TMath::Abs(IU_x-170.)<30.) && 
	    (TMath::Abs(IU_y-410.)<30.)){
	  
	  TrackHitCounter[0] += 1.;
	  if (FDC_Hits[0]>0){
	    TrackHitCounter[1] += 1.;	
	  }
	  
	}
	if ((TMath::Abs(IU_x-200.)<30.) && 
	    (TMath::Abs(IU_y-470.)<20.)){
	  
	  TrackHitCounter[2] += 1.;
	  if (FDC_Hits[0]>0){
	    TrackHitCounter[3] += 1.;	
	  }
	  
	}
	
	
	if ((TMath::Abs(IU_x*slopeW - 26.844 - IU_y)<11.5) &&
	    (TMath::Abs(IU_x-225)<50.)){
	  TrackHitCounter[4] += 1.;
	  if (FDC_Hits[0]>0){
	    TrackHitCounter[5] += 1.;	
	  }
	}


	if (FDC_Hits[0]==1){
	  
	  Int_t cnt1 = 0;
	  while ( cnt1<FDCHits){

	    if (FDCHits_Chan[cnt1]<272){
	      
	      Int_t wir =  FDCHits_Chan[cnt1] - 256 ;
	      Int_t c = wir + 19;
	      Float_t WPos = WirePlanePos;
	      
	      IU_x = slopex.val*WPos + interseptx.val;
	      IU_y = slopey.val*WPos + intersepty.val;
	      
	      //cout<<c<<"   "<<WPos<<"    "<< FDCHits_Chan[0]<<endl;
	      
	      hist2d[c]->Fill(IU_x, IU_y, 1.);

	      Double_t yfdc = slopeW*IU_x+B[wir];
	      Double_t d = (yfdc-IU_y)*cos(alphaFDC);
	      
	      c = wir + 55;
	      histT[c]->Fill(d); 
	      Double_t tim =  FDCHits_Time[cnt1] - DriftTimeZero[FDCHits_Chan[cnt1]];;

	      //cout<<"FDC time "<<tim<<endl;
	      if ((tim<168.)&&(tim>=0.)) {
		Double_t dFDC = time2distFDC(tim);
		Double_t d_diff = dFDC-d;

		
		// unit vector along track
		Float_t x1 = slopex.val;
		Float_t y1 = slopey.val;
		Float_t z1 = 1.;
		
		// vector normal to the direction of the FDC wire 
		// the norm of this vector is =1.15477541 
		Float_t x0 = -slopeW;
		Float_t y0 = 1.;
		Float_t z0 = 0.;
		
		// angle between plane and track
		Float_t alpha = (x1*x0 + y1*y0) / sqrt(x1*x1+y1*y1+z1*z1)/ sqrt(1.+slopeW*slopeW) ;
		
		alpha = acos(alpha)-PIh;
		
		//cout<<"angle "<<alpha<<endl;
		
		hist2d[36]->Fill(alpha, d_diff);
		
		// vector in the direction of the FDC wire 
		// the norm of this vector is =1.15477541 
		Float_t x2 = 1.;
		Float_t y2 = slopeW;
		Float_t z2 = 0.;
		alpha = (x2*x1 + y2*y1) / sqrt(x2*x2+y2*y2+z2*z2)/ sqrt(x1*x1+y1*y1+z1*z1);
		alpha = acos(alpha) - PIh;
		//cout<<"angle "<<alpha<<endl;
		hist2d[37]->Fill(alpha, d_diff);
		
		c = wir + 71;
		if ((TMath::Abs(slopex.val)<0.05) && (TMath::Abs(slopey.val)<0.05)){
		  histT[c]->Fill(d_diff);
		}
		//cout<<"FDC time "<<tim<<endl;
	      }

	    }
	    cnt1++;
	  }


	  // calculate angle between track projection in the x-y plane and FDC wire orientation
          // then add 90 degree
	  // angle between FDC wires and y axis is 0.53356 rad or 30.571 degree
	  if (FDC_Hits[1]==1){

	    // unit vector along track
	    Float_t x1 = slopex.val;
	    Float_t y1 = slopey.val;
	    Float_t z1 = 1.;

	    // vector normal to the direction of the FDC wire 
	    // the norm of this vector is 344.120168
	    Float_t x0 = 296.3;
	    Float_t y0 = -175.0;
	    Float_t z0 = 0.;

	    // angle between plane and track
	    Float_t alpha = (x1*x0 + y1*y0) / sqrt(x1*x1+y1*y1+z1*z1)/ 344.120168;

	    alpha = asin(alpha);

	    //cout<<"angle "<<alpha<<endl;

	    hist2d[36]->Fill(alpha, FDCDeltaT);

	    // also calculate angle between track and vertical axis (0,0,1)

	    alpha = 1./ sqrt(x1*x1+y1*y1+z1*z1);
	    alpha = acos(alpha);
	    //cout<<"angle "<<alpha<<endl;
	    hist2d[37]->Fill(alpha, FDCDeltaT);
	    


	  }

	}
	
	// fill the wire hits from the second wire plane 
	if (FDC_Hits[1]==1){ 
	  
	  Int_t cnt1 = 0;
	  while ( cnt1<FDCHits){
	    
	    if (FDCHits_Chan[cnt1]>271){
	      
	      Int_t c = FDCHits_Chan[cnt1] - 272 + 38;
	      
	      Float_t WPos = WirePlanePos+20.0;
	      
	      IU_x = slopex.val*WPos + interseptx.val;
	      IU_y = slopey.val*WPos + intersepty.val;
	      
	      //cout<<c<<"   "<<WPos<<"    "<< FDCHits_Chan[0]<<endl;
	      
	      hist2d[c]->Fill(IU_x, IU_y,1.);
	      
	      
	    }
	    cnt1++;
	  }
	}


      }

      // Now the FDC cathode data analysis
      // initialize planes
      
      
      memset(&Planes[0][0],0,384);
      
      //for (Int_t k=0;k<32;k++){
      //  Planes[0][k] = 0.;
      //  Planes[1][k] = 0.;
      //}
      PlaneCnt[0] = PlaneCnt[1] = 0;
      
      // fill raw ADC spectra for V792 CAEN ADC
      
      for (Int_t k=0; k<nCAENADC; k++){
	qadc[k]->Fill(CAEN_ADC[k]);
	hist2d[8]->Fill(CAEN_ADC[k],(Float_t)k);
	//hist2d[9]->Fill((CAEN_ADC[k]-Pedestals[k])*Gains[k]+200.,(Float_t)k);
	//cout << QDC_ped[0][k] << "  " << QDC_sig[0][k] << "  "<< Gains[k]<<endl;
	
      }
      
	
      if (Flag1){
	
	for (Int_t k=0; k<nCAENADC; k++){
	  
	  CAEN_ADC[k] = ( CAEN_ADC[k] - QDC_ped[0][k] ) * Gains[k];
	  
	}

	CAEN_ADC[15] = 0.;
	
	// try to do a common mode subtraction in groups of 8

	Float_t CM0[8];
	Float_t CM1[8];
	for (Int_t k = 0;k<8;k++){
	  CM0[k] = 0.;
	  CM1[k] = 0.;
	}

	for (Int_t k = 0;k<8;k++){
	  for (Int_t j = 0;j<8;j++){
	    Int_t ch = k*8 + j;
	    CM0[k] += CAEN_ADC[ch];
	  }
	  CM0[k] /= 8.;
	  CM0[k] += 100.;
	}
	for (Int_t k = 0;k<8;k++){
	  Float_t chcnt = 0.;
	  for (Int_t j = 0;j<8;j++){
	    Int_t ch = k*8 + j;
	    if (CAEN_ADC[ch]<CM0[k]){
	      CM1[k] += CAEN_ADC[ch];
	      chcnt +=1.0;
	    }
	  }
	  CM1[k] /= chcnt;
	}


	//for (Int_t k=0; k<8;k++){
	//  cout<< CM1[k]<< "  ";
	//}
	//cout<<endl;

	Int_t plane = 0;
	Int_t strip = 0;

	hist2d[54]->Fill(CAEN_ADC[2] + 200. ,CAEN_ADC[1] + 200.);
	hist2d[55]->Fill(CAEN_ADC[3] + 200. ,CAEN_ADC[1] + 200.);
	hist2d[56]->Fill(CAEN_ADC[4] + 200. ,CAEN_ADC[1] + 200.);
	hist2d[57]->Fill(CAEN_ADC[3] + 200. ,CAEN_ADC[2] + 200.);
	hist2d[58]->Fill(CAEN_ADC[4] + 200. ,CAEN_ADC[2] + 200.);
	hist2d[59]->Fill(CAEN_ADC[5] + 200. ,CAEN_ADC[2] + 200.);

	hist2d[60]->Fill(CAEN_ADC[34] + 200. ,CAEN_ADC[33] + 200.);
	hist2d[61]->Fill(CAEN_ADC[35] + 200. ,CAEN_ADC[34] + 200.);
	hist2d[62]->Fill(CAEN_ADC[36] + 200. ,CAEN_ADC[35] + 200.);
	hist2d[63]->Fill(CAEN_ADC[37] + 200. ,CAEN_ADC[36] + 200.);
	hist2d[64]->Fill(CAEN_ADC[38] + 200. ,CAEN_ADC[37] + 200.);
	hist2d[65]->Fill(CAEN_ADC[39] + 200. ,CAEN_ADC[38] + 200.);

	for (Int_t k=0; k<nCAENADC; k++){

	  Int_t ch =  k/8;

	  Float_t QADC =  CAEN_ADC[k] + 200. ;//- CM1[ch] ;

	  hist2d[9]->Fill( QADC,(Float_t)k);
	  qadc[k+64]->Fill(QADC);
	  
	  Int_t kk = k;
	  plane = 0;
	  if (k>31){
	    plane = 1;
	  }
	  if (CAEN_ADC[k]>QDC_lim[k]){
	    //cout<<i<<"   "<<k<<"  "<<CAEN_ADC[k]<<endl;

	    strip = FDC_QDC_Chan[k];

	    Planes[plane][strip] = CAEN_ADC[k];
	    PlaneCnt[plane]++;
	  }
	}
	//cout<<i<<"  " <<  PlaneCnt[0] << "  " << PlaneCnt[1] <<endl;	
      }
      

      // Fill FADC histograms
      //plane = 1;
      //for (Int_t k=0;k<nFADC;k++){
      
      //cout<<"FADC-SUM: "<<FADC_SUM[k]<<"           SUM-PED:"<<FADC_SUMSUB[k]<<endl;
      
      //fadc[k]->Fill(FADC_SUM[k],1);
      //FADC_SUMSUB[k]*= -1.;
      //fadcNP[k]->Fill(FADC_SUMSUB[k],1);
      
      //strip = k;
      //Float_t rnorm = 1.;
      //if (k<16)
      //strip = 15-k;
      //strip = k;
      //if (k>15)
      // strip = 47-k;
      //
      //if ((FADC_SUM[k])<FADC_lim[k]){
      //  if (k>15)
      //    rnorm = 0.2;
      //  Planes[plane][strip] = FADC_SUMSUB[k]*rnorm;
      //  PlaneCnt[plane]++;
      //}
      //
      //}
      
      
      if (Flag1){
	
	// determine Cathode positions
	Float_t xpos = 99999;
	Float_t ypos = 99999;
	Float_t wpos = 99999;
	Float_t lpos = 99999;
	Int_t inp;
	Int_t FitBoth = 0;
	if (PlaneCnt[0]>2){

	  // check if at least 3 hits are next neighours
	  Int_t plane = 0;
	  Int_t bool1=0;
	  for (Int_t k=0;k<46;k++){
	   
	    if ((Planes[plane][k]>0)&&(Planes[plane][k+1]>0)&&(Planes[plane][k+2]>0))
	      bool1++;
	  }
	  
	  if (bool1){

	    //cout<<"fit FDC Plane 0"<<endl;
	    TGraph* tgr = new TGraph(48,Xaxis,Planes[0]);
	    Double_t max = (Double_t)GetMaxBin(Planes[0],48);
	    tgr->Fit("gaus","Q","R",max-3.,max+3.);
	    TF1* gf = tgr->GetFunction("gaus");
	    //gf->SetParameter(1,max);
	    //gf->SetParameter(2,1.);
	    //tgr->Fit("gaus","Q","R",max-2.,max+3.);
	    xpos = gf->GetParameter(1) + XCathOffSet;
	    //gf->SetLineColor(2);
	    //tgr->Draw("AP*C");
	    //c1->Update();
	    //cout<<"INP: ";
	    //cin>>inp;
	    FitBoth++;
	  }
	}
	if (PlaneCnt[1]>2){

	  // check if at least 3 hits are next neighours
	  Int_t plane = 1;
	  Int_t bool1=0;
	  for (Int_t k=0;k<46;k++){
	   
	    if ((Planes[plane][k]>0)&&(Planes[plane][k+1]>0)&&(Planes[plane][k+2]>0))
	      bool1++;
	  }
	  
	  if (bool1){
	    
	    //cout<<"fit FDC Plane 1"<<endl;
	    TGraph* tgr = new TGraph(48,Xaxis,Planes[1]);
	    Double_t max = (Double_t)GetMaxBin(Planes[1],48);
	    //cout<<"MaxBin= "<<max<<endl;
	    tgr->Fit("gaus","Q","R",max-3.,max+3.);
	    TF1* gf = tgr->GetFunction("gaus");
	    //gf->SetParameter(1,max);
	    //gf->SetParameter(2,1.);
	    //tgr->Fit("gaus","Q","R",max-2.,max+3.);
	    ypos = gf->GetParameter(1) + YCathOffSet;
	    //gf->SetLineColor(4);
	    //tgr->Draw("AP*C");
	    //c1->Update();
	    //cout<<"INP: ";
	    //cin>>inp;
	    //c1->SaveAs("test2.jpg");
	    FitBoth++;
	  }
	}
	if (FitBoth==2){
	  wpos =  getCathodeXpos(xpos, ypos);
	  lpos =  getCathodeYpos(xpos, ypos);
	  
	  
	  //cout<<xpos<<"   "<<ypos<<"       ";
	  //cout<<wpos<<"   "<<lpos<<endl;
	  hist2d[1]->Fill(xpos,ypos,1);
	  hist2d[2]->Fill(wpos,lpos,1);
	  histT[33]->Fill(wpos,1);
	  histT[34]->Fill(lpos,1);
	  histT[35]->Fill(xpos,1);
	  histT[36]->Fill(ypos,1);
	  

	  if (muly==1){
	    //oldchi2 =  fit_trackY(5, lpos , &intersepty, &slopey);
	  
	  //cout<<"f(y) = "<<slopey.val<<"y + "<<intersepty.val<<endl;
	  //  histT[41]->Fill(oldchi2);
	  //  histT[43]->Fill(slopey.val);
	  //  if (TMath::Abs(slopey.val)<0.1)
	  //    histT[45]->Fill(oldchi2);
	  }
	  

	  if (FDCHits==1){
	    histT[38]->Fill(wpos,1);
	    hist2d[7]->Fill(wpos,lpos,1);

	    Int_t cntfdchits=0;
	    for (Int_t n=0;n<16;n++){
	      if (FDCHit[n]==1){
		cntfdchits++;
		hist2d[4]->Fill(wpos,(Float_t)n, 1);
		hist2d[6]->Fill(lpos,(Float_t)n, 1);

		//cout<<HitCounterX[0]<<" "<<HitCounterX[1]<<" "<<HitCounterX[2]<<" "<<HitCounterX[3]<<endl;
		//cout<<HitCounterX[0]<<" "<<HitCounterX[1]<<" "<<HitCounterX[2]<<" "<<HitCounterX[3]<<endl;

		//cout<<"wire: "<<n<<"   wpos="<<wpos<<endl;
		if ((muly==1)&&(mulx==1)){
		  hist2d[10]->Fill(IU_x,(Float_t)n, 1);
		  hist2d[13]->Fill(IU_y,(Float_t)n, 1);
		}
	      }
	    }
	    // if ((muly==1)&&(mulx==1))
	    //  cout<<"Event "<<i<<" total FDCHits is: "<<cntfdchits<<"  "<<IU_x<<"  "<<IU_y<<endl;
 
	    //if ((lpos>30.)&&(lpos<70.)){
	      
	      histT[37]->Fill(wpos,1);
	      
	      //}	
	  }

	  if ((muly==1)&&(mulx==1)){
	    
	    //if (TMath::Abs(slopex.val+0.02)<0.05)
	    //hist2d[11]->Fill(IU_x,wpos, 1);
	    
	    if ( (TMath::Abs(slopey.val)<0.45) && 
		 (TMath::Abs(slopex.val)<0.45) )  {
	      hist2d[12]->Fill(IU_y,lpos, 1);
	      histT[52]->Fill(lpos-IU_y,1);

	      // project lpol to f(x)=ax+b with a=-1 and b=0
	      Float_t x0 = (-IU_y+lpos)/2.;

	      if (FDCHits==1){
		Float_t fdcx = FDC_X[0].x + ((Float_t)FDC_N) *  FDC_X[0].driftdist;
		//Float_t fdcx = FDC_X[0].x + FDC_X[0].driftdist;
		hist2d[17]->Fill(IU_x, fdcx, 1);
		histT[53]->Fill(fdcx-IU_x,1);
		if (IU_x<280.){
		  histT[54]->Fill(fdcx-IU_x,1);
		}
	      }

	      histT[39]->Fill(TMath::Sqrt(2.)*TMath::Abs(x0));
	      if (IU_y<360){
		
		hist2d[11]->Fill(IU_x,wpos, 1);
		
		hist2d[14]->Fill(TMath::Sqrt(2.)*TMath::Abs(x0),slopey.val);
		
	      }

	      if (TMath::Abs(IU_y-260.)<20.)
		  histT[46]->Fill(TMath::Sqrt(2.)*TMath::Abs(x0));
	      if (TMath::Abs(IU_y-300.)<20.)
		  histT[47]->Fill(TMath::Sqrt(2.)*TMath::Abs(x0));
	      if (TMath::Abs(IU_y-340.)<20.)
		  histT[48]->Fill(TMath::Sqrt(2.)*TMath::Abs(x0));
	      if (TMath::Abs(IU_y-380.)<20.)
		  histT[49]->Fill(TMath::Sqrt(2.)*TMath::Abs(x0));
	      if (TMath::Abs(IU_y-320.)<20.)
		  histT[50]->Fill(TMath::Sqrt(2.)*TMath::Abs(x0));
	      if (TMath::Abs(IU_y-320.)<40.)
		  histT[51]->Fill(TMath::Sqrt(2.)*TMath::Abs(x0));


	    }	  
	  }
	}
      }
    }
  }
  
  fdctest->Delete();
  f->Close();


  cout<<"FDC WIRE COUNTERS:"<<endl;
  cout<<"Plane1   Plane2"<<endl;
  for (Int_t j=0;j<16;j++){
    cout << FDC_Hit_Counter[0][j] << "     " << FDC_Hit_Counter[1][j] << endl; 
  }

  Float_t effi[6] = {0., 0., 0., 0., 0., 0.};

  effi[0] = TrackHitCounter[1]/TrackHitCounter[0];
  effi[1] = effi[0]*sqrt(1./TrackHitCounter[1] + 1./TrackHitCounter[0]);
  effi[2] = TrackHitCounter[3]/TrackHitCounter[2];
  effi[3] = effi[2]*sqrt(1./TrackHitCounter[3] + 1./TrackHitCounter[2]);
  effi[4] = TrackHitCounter[5]/TrackHitCounter[4];
  effi[5] = effi[4]*sqrt(1./TrackHitCounter[5] + 1./TrackHitCounter[4]);


  cout<<"FDC plane 0 test efficiencies "<<endl;
  char out1[128];
  for (Int_t i = 0;i<3;i++){
    sprintf(out1,"%5.3f +/- %5.3f",effi[i*2],effi[i*2+1]);
    cout<<out1<<endl;
  }

  return 0;
}