//static char vcid[] = "$Id: main.cpp 254 2011-04-28 16:52:15Z xiangming Sun $";
//static const char __attribute__ ((used )) *Get_vcid(){return vcid;}


#include "log/def.h"

#include <TROOT.h>
#include <TStyle.h>
#include <TFile.h>
#include <TMath.h>
#include <TGraph.h>
#include <TGraphErrors.h>
#include <TCanvas.h>
#include <TH2D.h>
#include <TGraph2D.h>
#include <TH1D.h>
#include <TF1.h>
#include <TNtuple.h>
#include <TPaveStats.h>
#include <TLegend.h>





#include <stdio.h>
#include <stdlib.h>
#include <iostream>
#include <fstream>
#include <vector>
#include <bitset>

using namespace std;
#define  SEQUENCE_OFFSET 5



#ifdef PHASE1_H

#define NSENSOR 4
#define NROW 640
#define NCOL 640

#else

#define NSENSOR 3
#define NCOL 960
#define NROW 928

#endif

#define MAX_CLUSTER 1000



int help=0;
int gargc;char** gargv;int para(char *a);int para(char *a){if(help==1)cout<<a<<" ";int i;for(i=0;i<gargc;i++){if(strcmp(a,gargv[i])==0) return i;}return -1;}


int get_phase1_Event(ifstream & datafile, vector <int>  *x,vector <int>  *y);
int get_ultimate_Event(ifstream & datafile, vector <int>  *x,vector <int>  *y);


unsigned int shtoi(string strg){
unsigned int ival;
sscanf(strg.data(), "%x", &ival);
return ival;

};
unsigned int mid(int start,int end, unsigned int input){
unsigned int buf;
buf=input<<(32-end);
buf=buf>>(32-end);
//cout<<"mid "<<input<<endl;
return buf>>start;
};
int element_getbit(int data,int position){

int sd=data>>position;
return sd%2;

}

int element_setbit(unsigned int & data,int position){
if(position>31) return -1;

unsigned int mask=0;
mask=1<<position;
data=data|mask;
return 1;

}


/////////for ultimate^
int read_depth=1;
int nevent;
int nline;
int debug;
char buf[20];
unsigned int length[NSENSOR];
unsigned int header[5];
string sval;
int frame_number;
int sensor_data_length;
int frame_counter;
int number;
int data_length;
int row;
int column;
int OVF;
int coding;
int new_state;
int trigger_row;
int sensor;
int image;
int state=0;
int read_read=1;  
/////////for ultimate v





class pxl_cluster{
public:
int nhits;
vector<int> x; 
vector<int> y; 
float xc;
float yc;
void centering(){
float xsum=0;
float ysum=0;

for(int i=0;i<nhits;i++){
xsum=xsum+x[i];
ysum=ysum+y[i];
}
xc=xsum/float(nhits);
yc=ysum/float(nhits);
}

void clear(){
nhits=0;
x.clear();
y.clear();
}
};


class pxl_cluster_finder{
public:
  vector<pxl_cluster> clu[NSENSOR];
  bitset<NCOL> bsbit[NSENSOR][NROW];
  pxl_cluster tc; 
  // int cluster_hist[NSENSOR][MAX_CLUSTER];
 

 void clear_bsbit(){

 for(int k=0;k<NSENSOR;k++){
  for(int i=0;i<NROW;i++){
  bsbit[k][i].reset();
  } 
}
}; 

 void load_bsbit(vector <int>  *x,vector <int>  *y){
	for(int k=0;k<NSENSOR;k++){
		for(int i=0;i<x[k].size();i++){
		bsbit[k][y[k][i]].set(x[k][i]);
		}
	}
}; 
 
 void clear_cluster(){

 for(int k=0;k<NSENSOR;k++){
  // for(int i=0;i<MAX_CLUSTER;i++){
  // cluster_hist[k][i]=0;
  // } 
 for(int i=0;i<clu[k].size();i++) clu[k][i].clear();
 clu[k].clear();
 }
 }; 
 void clustering(vector <int>  *x,vector <int>  *y){
 
  for(int k=0;k<NSENSOR;k++){

  for (int j=0; j<x[k].size(); j++) {
  tc.clear();
  find_cluster(k,x[k][j],y[k][j],1);
  if(tc.nhits!=0) 
  {
    tc.centering();
  clu[k].push_back(tc);
    if(tc.nhits>MAX_CLUSTER) return;
  }
  }
  }
 
 };
 
  void clustering_clearbit(vector <int>  *x,vector <int>  *y){
  
  for(int k=0;k<NSENSOR;k++){

  for (int j=0; j<x[k].size(); j++) {
  tc.clear();
  find_cluster(k,x[k][j],y[k][j],0);
  if(tc.nhits!=0) 
  {
    tc.centering();
  clu[k].push_back(tc);
    if(tc.nhits>MAX_CLUSTER) return;
  }
  }
  }
 
 };
 
  void find_cluster(int sen,int x,int y,int keep){
  if(bsbit[sen][y][x]==0) return;
  bsbit[sen][y][x]=0;
  if(x+1<NCOL) find_cluster(sen,x+1,y,keep);
  if(x-1>=0) find_cluster(sen,x-1,y,keep);
  if(y+1<NROW) find_cluster(sen,x,y+1,keep);
  if(y-1>=0) find_cluster(sen,x,y-1,keep);
  tc.nhits++;
  tc.x.push_back(x);
  tc.y.push_back(y);
  bsbit[sen][y][x]=keep;
  };
};








//------------------------------------------------------
int correlation(int type, char *fn)
{



#ifdef PHASE1_H
int offsetx[4]={0,0,-4,-14};
int offsety[4]={0,-22,-49,-70};

cout<<"phase1"<<endl;

#ifdef PDSF_H
char filename[2000]="/eliza17/rnc/sunxm/cern_beam_test/phase1/phase1.txt";
char path[2000]="/eliza17/rnc/sunxm/cern_beam_test/phase1/";
char result_path[2000]="/eliza17/rnc/sunxm/cern_beam_test/phase1_result/";
#else
char filename[2000]="/data1/CERN_beam_test/phase1/phase1.txt";
char path[2000]="/data1/CERN_beam_test/phase1/";
char result_path[2000]="/data1/CERN_beam_test/phase1_result/";
#endif


#else
int offsetx[4]={0,15,31,0};
int offsety[4]={0,-38,-60,0};


cout<<"ultimate"<<endl;


#ifdef PDSF_H
char filename[2000]="/eliza17/rnc/sunxm/cern_beam_test/ultimate/ultimate.txt";
char path[2000]="/eliza17/rnc/sunxm/cern_beam_test/ultimate/";
char result_path[2000]="/eliza17/rnc/sunxm/cern_beam_test/ultimate_result/";
#else
char filename[2000]="/data1/CERN_beam_test/ultimate/ultimate.txt";
char path[2000]="/data1/CERN_beam_test/ultimate/";
char result_path[2000]="/data1/CERN_beam_test/ultimate_result/";
#endif



#endif



int nfile=0;
char buffer[2000];
char data_file[2000];
char root_file[2000];


pxl_cluster_finder cf;





#ifdef PDSF_H
sprintf(buffer,"%s",fn);
if(strstr(buffer,".txt")==NULL)return -1;

#else
ifstream in;
in.open(filename);
do{
in>>buffer;
if(!in.good()) break;
if(in.eof()) break;
if(strstr(buffer,".txt")==NULL)continue;

#endif



// cout<<data_file<<endl;
// cout<<root_file<<endl;


#ifdef PDSF_H

#else
nfile++;
 if(nfile!=type) continue;
// if(nfile<88) continue;
char command[2000];
sprintf(command,"./ultimate2 %s",buffer);
system(command);
cout<<buffer<<" converted"<<endl;


#endif

sprintf(data_file,"%s%s",path,buffer);
sprintf(root_file,"%s%s.root",result_path,buffer);


vector <int>  x[NSENSOR];
vector <int>  y[NSENSOR];
ifstream datafile;

TFile * correlation_root_file=new TFile(root_file,"RECREATE");

TH1F *hx[NSENSOR];
TH1F *hy[NSENSOR];
TH1F *hx_align[NSENSOR];
TH1F *hy_align[NSENSOR];
TH1F *cluster_nhits[NSENSOR];
for(int i=0;i<NSENSOR;i++){
char n1[200];
char n2[200];
sprintf(n1,"correlation-x_%d",i);
sprintf(n2,"correlation-x_%d",i);
hx[i]=new TH1F(n1,n2,2000,-1000,1000);
sprintf(n1,"correlation-y_%d",i);
sprintf(n2,"correlation-y_%d",i);
hy[i]=new TH1F(n1,n2,2000,-1000,1000);

sprintf(n1,"correlation_align-x_%d",i);
sprintf(n2,"correlation_align-x_%d",i);
hx_align[i]=new TH1F(n1,n2,2000,-1000,1000);
sprintf(n1,"correlation_align-y_%d",i);
sprintf(n2,"correlation_align-y_%d",i);
hy_align[i]=new TH1F(n1,n2,2000,-1000,1000);

sprintf(n1,"cluster_%d",i);
sprintf(n2,"nhits_%d",i);
cluster_nhits[i]=new TH1F(n1,n2,100,0,100);
}


TH1F *candidate=new TH1F("candidate","candidate_ref",100,0,100);
TH1F *candidate_match=new TH1F("candidate_match","candidate_match_ref",100,0,100);


datafile.open(data_file);

if(datafile.good()) cout<<"good"<<endl;

cf.clear_cluster();

nevent=0;
#ifdef PHASE1_H

while(get_phase1_Event(datafile,x,y)==0){
nevent++;
#else

read_depth=1;
read_read=1; 
state=0;

while(get_ultimate_Event(datafile,x,y)==0){
#endif
if(nevent%1000==0) cout<<nevent<<endl;
// if(nevent>500) break;

cf.load_bsbit(x,y);
cf.clustering_clearbit(x,y);


for(int i=0;i<x[0].size();i++){
for(int k=1;k<NSENSOR;k++){
	for(int j=0;j<x[k].size();j++){
		hx[k]->Fill((x[k][j])-(x[0][i]));
		hy[k]->Fill((y[k][j])-(y[0][i]));
	}
}
}




for(int i=0;i<cf.clu[0].size();i++){
for(int k=1;k<NSENSOR;k++){
	for(int j=0;j<cf.clu[k].size();j++){
		hx_align[k]->Fill((cf.clu[k][j].xc+offsetx[k])-(cf.clu[0][i].xc+offsetx[0]));
		hy_align[k]->Fill((cf.clu[k][j].yc+offsety[k])-(cf.clu[0][i].yc+offsety[0]));
	}
}
}

for(int k=0;k<NSENSOR;k++){
for(int i=0;i<cf.clu[k].size();i++){
cluster_nhits[k]->Fill(cf.clu[k][i].nhits);
}
}


int ref1=0;
int ref2=NSENSOR-1;
int meas=1;
float window=8;
int cand=0;
int cand_match=0;
for(int i=0;i<cf.clu[ref1].size();i++){
	if(cf.clu[ref1][i].nhits==0) continue;
	float ref1x=cf.clu[ref1][i].xc+offsetx[ref1];
	float ref1y=cf.clu[ref1][i].yc+offsety[ref1];
for(int h=0;h<cf.clu[ref2].size();h++){
	if(cf.clu[ref2][h].nhits==0) continue;
	float ref2x=cf.clu[ref2][h].xc+offsetx[ref2];
	float ref2y=cf.clu[ref2][h].yc+offsety[ref2];
	
if(fabs(ref1x-ref2x)>window) continue;
if(fabs(ref1y-ref2y)>window) continue;		
	
	cand++;
for(int j=0;j<cf.clu[meas].size();j++){
	if(cf.clu[meas][j].nhits==0) continue;
	float measx=cf.clu[meas][j].xc+offsetx[meas];
	float measy=cf.clu[meas][j].yc+offsety[meas];
	
if(fabs(ref1x-measx)>window) continue;
if(fabs(ref1y-measy)>window) continue;	
cand_match++;
cf.clu[meas][j].nhits=0;
cf.clu[ref2][h].nhits=0;
cf.clu[ref1][i].nhits=0;
break;
	}
}
}

candidate->Fill(cand);
candidate_match->Fill(cand_match);

cf.clear_cluster();


for(int i=0;i<NSENSOR;i++){
x[i].clear();
y[i].clear();
}

}
datafile.close();
datafile.clear();
correlation_root_file->Write();
correlation_root_file->Close();


#ifdef PDSF_H

#else
}while(!in.eof());
#endif




return 1;




// int wait;cout<<"press any key + enter to go!"<<endl;
// scanf("%d",&wait);
  return 1;
}

/*
int DisplayPixelDialog::get_phase1_Event()
{
  unsigned int value;
  unsigned int sequence;
  bool bSensor1and2 = false;
  int sensor;
  int row;
  int col;
  int remainder;
  int numEventFlags;

  numEventFlags = 0;
  datafile >> hex;

#ifdef LOCALDEBUG
  cout << hex;
#endif

  while (datafile.good()) {
    datafile >> value;
    if (!datafile.good()) break;

#ifdef LOCALDEBUG
    cout << "0x";
    cout.width(8); cout.fill('0');
    cout  << value;
#endif

    int dataType = int( (value & 0xffff0000)>> 16 );
    switch (dataType) {

    case 0xaaaa: // header
#ifdef LOCALDEBUG
      if ( (value&0xaaaa) == 0xaaaa)
	cout << "\tStart of  event\n";
      else
	cout << "\t1st header word\n";
#endif
      break;

    case 0xbbbb: // header
#ifdef LOCALDEBUG
      cout << "\t2nd header word\n";
#endif
      break;

    case 0xccc0: // header
    case 0xccc1: // header
    case 0xcccc: // header
#ifdef LOCALDEBUG
      cout << "\t3rd header word\n";
#endif
      break;

    case 0xdddd: // header
#ifdef LOCALDEBUG
      cout << "\t4th header word\n";
#endif
      break;

    case 0xeeee: // end-of-event
#ifdef LOCALDEBUG
      cout << "\tend of event\n";
#endif
      if (numEventFlags == 2) return 0;
      break;

    case 0xf777:
    case 0x7777: // event flag
      numEventFlags++;
      bSensor1and2 = ((value & 0xf) == 0);
#ifdef LOCALDEBUG
      cout << "\tevent flag: ";
      if ( bSensor1and2 )
	cout << "sensor 1 and 2\n";
      else
	cout << "sensor 3 and 4\n";
#endif
      break;

    default: // data
      sequence = ( (value & 0xfffff000)>>12 );
#ifdef LOCALDEBUG
      cout << "\tsequence 0x";
      cout.width(5); cout.fill('0');
      cout << sequence;
      cout.fill(' ');
#endif
      
      sequence += SEQUENCE_OFFSET;
      row = ((int)(sequence / 160))%640;
      remainder = sequence % 160;
      col = ((remainder % 4)*40 + (remainder/4))%160;

#ifdef LOCALDEBUG
      cout << dec;
#endif

      if (bSensor1and2) { // Sensors 1 and 2
	if ( (value&0xf) != 0) { // Sensor 1
	  sensor = 1;
	  if ( (value&0x1) == 0x1) {
#ifdef LOCALDEBUG
	    _PS(sensor, row, col);
#endif
	    y[0].push_back(row);x[0].push_back(col);
	  }
	  if ( (value&0x2) == 0x2) {
#ifdef LOCALDEBUG
	    _PS(sensor, row, col+160);
#endif
	    y[0].push_back(row);x[0].push_back(col+160);
	  }
	  if ( (value&0x4) == 0x4) {
#ifdef LOCALDEBUG
	    _PS(sensor, row, col+320);
#endif
	    y[0].push_back(row);x[0].push_back(col+320);
	  }
	  if ( (value&0x8) == 0x8) {
#ifdef LOCALDEBUG
	    _PS(sensor, row, col+480);
#endif
	    y[0].push_back(row);x[0].push_back(col+480);
	  }
	}
	if ( (value&0xf0) != 0 ) { // Sensor 2
	  sensor = 2;
	  if ( (value&0x10) == 0x10) {
#ifdef LOCALDEBUG
	    _PS(sensor, row, col);
#endif
	    y[1].push_back(row);x[1].push_back(col);
	  }
	  if ( (value&0x20) == 0x20) {
#ifdef LOCALDEBUG
	    _PS(sensor, row, col+160);
#endif
	    y[1].push_back(row);x[1].push_back(col+160);
	  }
	  if ( (value&0x40) == 0x40) {
#ifdef LOCALDEBUG
	    _PS(sensor, row, col+320);
#endif
	    y[1].push_back(row);x[1].push_back(col+320);
	  }
	  if ( (value&0x80) == 0x80) {
#ifdef LOCALDEBUG
	    _PS(sensor, row, col+480);
#endif
	    y[1].push_back(row);x[1].push_back(col+480);
	  }
	}
      }

      else { // Sensors 3 and 4
	if ( (value&0xf) != 0 ) { // Sensor 3
	  sensor = 3;
	  if ( (value&0x1) == 0x1) {
#ifdef LOCALDEBUG
	    _PS(sensor, row, col);
#endif
	    y[2].push_back(row);x[2].push_back(col);
	  }
	  if ( (value&0x2) == 0x2) {
#ifdef LOCALDEBUG
	    _PS(sensor, row, col+160);
#endif
	    y[2].push_back(row);x[2].push_back(col+160);
	  }
	  if ( (value&0x4) == 0x4) {
#ifdef LOCALDEBUG
	    _PS(sensor, row, col+320);
#endif
	    y[2].push_back(row);x[2].push_back(col+320);
	  }
	  if ( (value&0x8) == 0x8) {
#ifdef LOCALDEBUG
	    _PS(sensor, row, col+480);
#endif
	    y[2].push_back(row);x[2].push_back(col+480);
	  }
	}
	if ( (value&0xf0) != 0 ) { // Sensor 4
	  sensor = 4;
	  if ( (value&0x10) == 0x10) {
#ifdef LOCALDEBUG
	    _PS(sensor, row, col);
#endif
	    y[3].push_back(row);x[3].push_back(col);
	  }
	  if ( (value&0x20) == 0x20) {
#ifdef LOCALDEBUG
	    _PS(sensor, row, col+160);
#endif
	    y[3].push_back(row);x[3].push_back(col+160);
	  }
	  if ( (value&0x40) == 0x40) {
#ifdef LOCALDEBUG
	    _PS(sensor, row, col+320);
#endif
	    y[3].push_back(row);x[3].push_back(col+320);
	  }
	  if ( (value&0x80) == 0x80) {
#ifdef LOCALDEBUG
	    _PS(sensor, row, col+480);
#endif
	    y[3].push_back(row);x[3].push_back(col+480);
	  }
	}
      }
      
      // Marker
      if ( (value&0x100) == 0x100) {
#ifdef LOCALDEBUG
	cout << " Marker ";
#endif
      }

#ifdef LOCALDEBUG
      cout << endl << hex;
#endif
      break;
    }


  }   // end while loop 


  // got out of the "while" loop, so must have bad file status
  if (datafile.eof()) {
#ifdef LOCALDEBUG
    cout << "End of File\n";
#endif
    posLabel->setText("End of File");
  }
  else {
#ifdef LOCALDEBUG
    cout << "Bad datafile status\n";
#endif
    posLabel->setText("Bad file status on datafile");
  }
  return 1;

}
*/



int get_phase1_Event(ifstream & datafile, vector <int>  *x,vector <int>  *y)
{



  
  unsigned int value;
  unsigned int sequence;
  bool bSensor1and2 = false;
  int sensor;
  int row;
  int col;
  int remainder;
  int numEventFlags;

  numEventFlags = 0;
  datafile >> hex;

#ifdef LOCALDEBUG
  cout << hex;
#endif

  while (datafile.good()) {
    datafile >> value;
    if (!datafile.good()) break;

#ifdef LOCALDEBUG
    cout << "0x";
    cout.width(8); cout.fill('0');
    cout  << value;
#endif

    int dataType = int( (value & 0xffff0000)>> 16 );
    switch (dataType) {

    case 0xaaaa: // header
#ifdef LOCALDEBUG
      if ( (value&0xaaaa) == 0xaaaa)
	cout << "\tStart of  event\n";
      else
	cout << "\t1st header word\n";
#endif
      break;

    case 0xbbbb: // header
#ifdef LOCALDEBUG
      cout << "\t2nd header word\n";
#endif
      break;

    case 0xccc0: // header
    case 0xccc1: // header
    case 0xcccc: // header
#ifdef LOCALDEBUG
      cout << "\t3rd header word\n";
#endif
      break;

    case 0xdddd: // header
#ifdef LOCALDEBUG
      cout << "\t4th header word\n";
#endif
      break;

    case 0xeeee: // end-of-event
#ifdef LOCALDEBUG
      cout << "\tend of event\n";
#endif
      if (numEventFlags == 2) return 0;
      break;

    case 0xf777:
    case 0x7777: // event flag
      numEventFlags++;
      bSensor1and2 = ((value & 0xf) == 0);
#ifdef LOCALDEBUG
      cout << "\tevent flag: ";
      if ( bSensor1and2 )
	cout << "sensor 1 and 2\n";
      else
	cout << "sensor 3 and 4\n";
#endif
      break;

    default: // data
      sequence = ( (value & 0xfffff000)>>12 );
#ifdef LOCALDEBUG
      cout << "\tsequence 0x";
      cout.width(5); cout.fill('0');
      cout << sequence;
      cout.fill(' ');
#endif
      
      sequence += SEQUENCE_OFFSET;
      row = ((int)(sequence / 160))%640;
      remainder = sequence % 160;
      col = ((remainder % 4)*40 + (remainder/4))%160;

#ifdef LOCALDEBUG
      cout << dec;
#endif

      if (bSensor1and2) { // Sensors 1 and 2
	if ( (value&0xf) != 0) { // Sensor 1
	  sensor = 1;
	  if ( (value&0x1) == 0x1) {
#ifdef LOCALDEBUG
	    _PS(sensor, row, col);
#endif
	    y[0].push_back(row);x[0].push_back(col);
	  }
	  if ( (value&0x2) == 0x2) {
#ifdef LOCALDEBUG
	    _PS(sensor, row, col+160);
#endif
	    y[0].push_back(row);x[0].push_back(col+160);
	  }
	  if ( (value&0x4) == 0x4) {
#ifdef LOCALDEBUG
	    _PS(sensor, row, col+320);
#endif
	    y[0].push_back(row);x[0].push_back(col+320);
	  }
	  if ( (value&0x8) == 0x8) {
#ifdef LOCALDEBUG
	    _PS(sensor, row, col+480);
#endif
	    y[0].push_back(row);x[0].push_back(col+480);
	  }
	}
	if ( (value&0xf0) != 0 ) { // Sensor 2
	  sensor = 2;
	  if ( (value&0x10) == 0x10) {
#ifdef LOCALDEBUG
	    _PS(sensor, row, col);
#endif
	    y[1].push_back(row);x[1].push_back(col);
	  }
	  if ( (value&0x20) == 0x20) {
#ifdef LOCALDEBUG
	    _PS(sensor, row, col+160);
#endif
	    y[1].push_back(row);x[1].push_back(col+160);
	  }
	  if ( (value&0x40) == 0x40) {
#ifdef LOCALDEBUG
	    _PS(sensor, row, col+320);
#endif
	    y[1].push_back(row);x[1].push_back(col+320);
	  }
	  if ( (value&0x80) == 0x80) {
#ifdef LOCALDEBUG
	    _PS(sensor, row, col+480);
#endif
	    y[1].push_back(row);x[1].push_back(col+480);
	  }
	}
      }

      else { // Sensors 3 and 4
	if ( (value&0xf) != 0 ) { // Sensor 3
	  sensor = 3;
	  if ( (value&0x1) == 0x1) {
#ifdef LOCALDEBUG
	    _PS(sensor, row, col);
#endif
	    y[2].push_back(row);x[2].push_back(col);
	  }
	  if ( (value&0x2) == 0x2) {
#ifdef LOCALDEBUG
	    _PS(sensor, row, col+160);
#endif
	    y[2].push_back(row);x[2].push_back(col+160);
	  }
	  if ( (value&0x4) == 0x4) {
#ifdef LOCALDEBUG
	    _PS(sensor, row, col+320);
#endif
	    y[2].push_back(row);x[2].push_back(col+320);
	  }
	  if ( (value&0x8) == 0x8) {
#ifdef LOCALDEBUG
	    _PS(sensor, row, col+480);
#endif
	    y[2].push_back(row);x[2].push_back(col+480);
	  }
	}
	if ( (value&0xf0) != 0 ) { // Sensor 4
	  sensor = 4;
	  if ( (value&0x10) == 0x10) {
#ifdef LOCALDEBUG
	    _PS(sensor, row, col);
#endif
	    y[3].push_back(row);x[3].push_back(col);
	  }
	  if ( (value&0x20) == 0x20) {
#ifdef LOCALDEBUG
	    _PS(sensor, row, col+160);
#endif
	    y[3].push_back(row);x[3].push_back(col+160);
	  }
	  if ( (value&0x40) == 0x40) {
#ifdef LOCALDEBUG
	    _PS(sensor, row, col+320);
#endif
	    y[3].push_back(row);x[3].push_back(col+320);
	  }
	  if ( (value&0x80) == 0x80) {
#ifdef LOCALDEBUG
	    _PS(sensor, row, col+480);
#endif
	    y[3].push_back(row);x[3].push_back(col+480);
	  }
	}
      }
      
      // Marker
      if ( (value&0x100) == 0x100) {
#ifdef LOCALDEBUG
	cout << " Marker ";
#endif
      }

#ifdef LOCALDEBUG
      cout << endl << hex;
#endif
      break;
    }


  }   // end while loop 


  // got out of the "while" loop, so must have bad file status
  if (datafile.eof()) {
#ifdef LOCALDEBUG
    cout << "End of File\n";
#endif
     cout << "End of File"<<endl;
	 return 1;
  }
  else {
#ifdef LOCALDEBUG
    cout << "Bad datafile status\n";
#endif
     cout << "Bad file status on datafile"<<endl;


  }
  return 1;

};



int get_ultimate_Event(ifstream & datafile, vector <int>  *x,vector <int>  *y)
{






  //char bs[4*NROW*NCOL];

 // int selectedSensor;
 // int selectedQuarter;
 // bool useFiltered;

 //bitset<NCOL> bsbit[NSENSOR][NROW];
  //bitset<NCOL> bsbit_rotate[NROW];
 // bitset<NCOL> bsbit_filter[NSENSOR][NROW];






int nEvents=0;


do{



	if(datafile.eof()) break;
	if(!datafile.good()) break;

	if(read_read==1){
	datafile.getline(buf,12);nline++;
	sval=buf;
	}
//////////////////////////	
	if(state==0||state==3){
	
		
	


		if(sval.substr(0,10)=="0xaaaaaaaa"){
		state=1;// find the header worder
			
		if(debug>3)cout<<"0xaaaaaaaa found"<<endl;
		
//		fprintf(out,"0xaaaaaaaa\n");	


			nevent++;
	//		if(do_cluster==1){
 //if(acumulate_event==0||nevent==0){ 
 // clear_cluster();
 //} 
 // clustering();
 // // print_cluster();
 // }
		}else{
			if(!datafile.eof())
				//cout<<"no 0xaaaaaaaa found "<<sval.data()<<" nline "<<nline<<endl;
		break;
		}

		if(read_read==0) read_read=1;
	
	}

//////////////////////////
	if(state==1){
		for(int i=0;i<5;i++){
		if(datafile.eof()) break;
		datafile.getline(buf,12);nline++;
		sval=buf;
		header[i]=shtoi(sval);	
		}
		trigger_row=(header[1]&(0xfffff))>>5;

		if(sval.substr(0,6)=="0xdddd"){
		state=4;// find the last word in header;
		//write header;
		for(int i=0;i<5;i++){
//		fprintf(out,"0x%08lx\n",header[i]);	
		}
		frame_number=0;

		continue;
		}else{
			cout<<"no 0xdddd found"<<endl;
		break;
		}	
	}

////////////////////////
	if(state==2||state==4){

		if(sval.substr(0,10)=="0xeeeeeeee"){
			frame_number++;
			if(frame_number==6){
				state=3; //one event is end
				nEvents++;
				if(debug>1)cout<<"frame found"<<endl;
			}else state=4;// one sensor is end, the start of the new sensor
				if(debug>1)cout<<"from eeeeeeeeee"<<endl;
			continue;
		}else if(sval.substr(0,10)=="0xffffffff"&&state==2){
			if(debug>0)cout<<"trucated frame!!!"<<" nline "<<nline<<endl;
			frame_number++;
			if(frame_number>=6){
				state=3; //one event is end
				nEvents++;
				cout<<"unexpected ffffffff"<<" nline "<<nline<<endl;
				break;
			}else state=4;// one sensor is end
			read_read=0;
			continue;
		}else if(sval.substr(0,10)=="0xaaaaaaaa"){
			if(frame_number!=5){
				// cout<<"unexpected aaaaaaaa"<<" nline "<<nline<<endl;
				cout<<"?";
				read_read=0;
				state=3;
				continue;			
			}
			state=3;
			nEvents++;
			read_read=0;
			continue;
		}else{

			int val=shtoi(sval); 
			int val1=val>>16;
			int val0=val&(0xffff);

			if(state==4){
				if(debug>2)cout<<"start of sensor "<<sval.data()<<endl;
				state=2;
				if(debug>2)cout<<"previous length "<<sensor_data_length<<endl;		
				sensor_data_length=0;
				new_state=0;
			}
			if(sensor_data_length==1){
				if(debug>2)cout<<"frame_counter "<<sval.data()<<endl;
				frame_counter=val;
			}
			if(sensor_data_length==2){
				int dl=val;
				int data_length1=dl>>16;
				int data_length0=dl&(0xffff);
				if(data_length1!=data_length0)
				{
					//cout<<"data_length1!=data_length0 "<<" nline "<<nline<<" sval "<<sval.data()<<endl;
					state=0;	
					break;
				}
				data_length=data_length0;

			}			
//int frame_number=0;
//int sensor_data_length=0;
//int frame_counter=0;
//int data_length;
//int number=0;
//int row;
//int column;
//int OVF;
//int coding;
//int new_state;
sensor=frame_number%3;
image=frame_number/3;

			if(sensor_data_length>2){
					
				if(new_state==0){
				
					number=val0&(0xf);
					row=(val0>>4)&(0x3ff);
					OVF=element_getbit(val0,15);
					new_state++;
					coding=val1&3;
					column=(val1>>2)&(0x3ff); 				
					new_state++;
					if((row>trigger_row&&image==0)||(row<=trigger_row&&image==1)){
						if(debug>0) cout<<"hit founded 1"<<" nline "<<nline<<endl;
						//hit_number[sensor][column/240]=//hit_number[sensor][column/240]+coding+1;
						for(int c=0;c<coding+1;c++)
						if(row<NROW&&column<NCOL)  {x[sensor].push_back(column+c);y[sensor].push_back(row);}
						// if(bs[4*(row*NCOL+column)+0]<239)
						// bs[4*(row*NCOL+column)+0]+=80;
						// bs[4*(row*NCOL+column)+sensor+1]=0;

					}
					if(new_state==number+1){
						new_state=0;
					}
					//cout<<"val0 "<<val0<<" val1 "<<val1<<" val "<<val<<endl;
					//cout<<"number"<<number<<" sval.data() "<<sval.data()<<endl;
				}else{
					coding=val0&3;
					column=(val0>>2)&(0x3ff);  new_state++;	
					if((row>trigger_row&&image==0)||(row<=trigger_row&&image==1)){
						if(debug>0) cout<<"hit founded 2"<<" nline "<<nline<<endl;
						//hit_number[sensor][column/240]=//hit_number[sensor][column/240]+coding+1;
						for(int c=0;c<coding+1;c++)
						if(row<NROW&&column<NCOL)  {x[sensor].push_back(column+c);y[sensor].push_back(row);}

						// if(bs[4*(row*NCOL+column)+0]<239)
						// bs[4*(row*NCOL+column)+0]+=80;
						// bs[4*(row*NCOL+column)+sensor+1]=0;

					}

					if(new_state==number+1){
						new_state=0;
						
					number=val1&(0xf);
					row=(val1>>4)&(0x3ff);
					OVF=element_getbit(val1,15);
					new_state++;
					if(new_state==number+1){
						new_state=0;
					}						
						
					}else{
					coding=val1&3;
					column=(val1>>2)&(0x3ff);  new_state++;					
					if((row>trigger_row&&image==0)||(row<=trigger_row&&image==1)){
						if(debug>0) cout<<"hit founded 3"<<" nline "<<nline<<endl;
						//hit_number[sensor][column/240]=//hit_number[sensor][column/240]+coding+1;
						for(int c=0;c<coding+1;c++)
						if(row<NROW&&column<NCOL)  {x[sensor].push_back(column+c);y[sensor].push_back(row);}
						
						// if(bs[4*(row*NCOL+column)+0]<239)
						// bs[4*(row*NCOL+column)+0]+=80;
						// bs[4*(row*NCOL+column)+sensor+1]=0;

					}
						
						if(new_state==number+1){
						new_state=0;
					}					
					}				
				
				}

				

			}
			sensor_data_length++;
			if(read_read==0) read_read=1;
			
		}


	}

	if(read_depth>0){
	if(nEvents==read_depth){
	if(state!=0&&state!=3) break;
	}
	}
}while(!datafile.eof());


//fclose(out);

if(datafile.eof()) return 1;

  return 0;

}


int batch_correlation(int start, int stop){
for(int i=start;i<stop;i++){
cout<<"********************* "<<i<<"*********************"<<endl;
correlation(i,"");
}
return 0;
}