//.L mppc.C++
//mppc("/afs/f9.ijs.si/scratch/rok/stm_timing0.dat.gz",0,-1,9,250,400)
//plothisto("h2dadctdc0",0,3,0,"","%s%02d",0,1500)
//plothisto("hadc0",9,11,1,"","%s%02d",0,1500)
// mppc(char *fname, int firsteve=0, int numeve=-1, int ch=-1, float xmin=1500, float xmax=2500 )
#include <stdlib.h>
#include <stdio.h>

#include "TH1F.h"
#include "TH2F.h"
#include "TH3F.h"
#include "TMath.h"
#include "TVector3.h"
#include "data.h"
#include <zlib.h>
#include <vector>
#include <fstream>
#include <iostream>

#define MAXCH 32
#define ADCRANGE 0x1000	 

#define RUNREC_ID 1
#define ENDREC_ID 2
#define POSREC_ID 3
#define EVTREC_ID 4
#define DACREC_ID 5 

DACREC dac;

class Hit {
public:
  TVector3 r;
  float energy;
  float time;
  Hit(){};
  Hit(TVector3 r0, float E, float t){r=r0;energy=E;time=t;};
  ~Hit(){};

};

std::vector <Hit> Hit_list[2];

std::vector <TVector3> Hit_Position;

class Coincidence {
public:
  Hit hit[2];
  Coincidence(Hit a,Hit b){hit[0]=a;hit[1]=b;}
};

std::vector <Coincidence> Coincidence_list;


TH3F * hvoxel;

TH1F * hadc[2][MAXCH];
TH1F * hadcmax[2][MAXCH];
TH2F * h2d[2][MAXCH];
TH2F * h2dadcdac[2][MAXCH];
TH2F * h2dmax[2][MAXCH];
TH2F * h2dadctdc[2][MAXCH];
TH2F * h2dadcadc[MAXCH];
TH2F * h2dposadcfrac[MAXCH];
TH2F * h2dposadc[MAXCH];

#define TDCRANGE 0x10000
TH1F * htdc[2][MAXCH];
TH1F * henergy[2][MAXCH];
TH2F * h2dtdc[2][MAXCH];

double aadc[2][MAXCH];
double valtdc[2][MAXCH];
double valadc[2][MAXCH];

int threshold[2]={0,0};
  
int SetThreshold(int a, int b){
  threshold[0]=a;
  threshold[1]=b;
  return 0;
}

typedef struct {
float dx[32];
float sx[32];
float n0[32];
} CALIBRATION;
CALIBRATION cal; 

float GetEnergy(float x, int ch){
	//return x;
if (ch>15) return x;
return -cal.sx[ch]*log(1-((x-cal.dx[ch])/cal.n0[ch]));
}

int readcalib(char *fname,int n,float *x,int nvar){

FILE *fp=fopen(fname,"r");
if (!fp) {
printf("No File %s\n",fname);
return -1;
} else {
printf("File %s\n",fname);
}
 int ndim=400;
  char *line = new char[ndim];

char fmt[200]="";
for (int i=0;i<nvar;i++) sprintf(fmt,"%s%s",fmt,"%" "*s");
sprintf(fmt,"%s%s",fmt,"%" "f");
printf("%s\n",fmt);
int nl=0;
while (fgets(line,ndim,fp)!=NULL) {

 if (nl<n) {
   int nb=sscanf(line,fmt,&x[nl]);
//   printf("%d %f-----%s\n",nb,x[nl],line);
   nl++;
 } else break;
}
fclose(fp);
return 0;
 
}

TVector3 ExtrapolateToZ(float z, TVector3 a, TVector3 b){

  float t= (a.z()-z)/(a.z()-b.z());
  return TVector3(a.x()-(a.x()-b.x())*t,a.y()-(a.y()-b.y())*t, z);
}


int hvoxel_fill(){
  int n= hvoxel->GetNbinsZ();
  for (int k=0;k<n;k++){
    float z=hvoxel->GetZaxis()->GetBinCenter(k);
    for (int i=0; i< Coincidence_list.size();i++){
      TVector3 r= ExtrapolateToZ(z,Coincidence_list[i].hit[0].r, Coincidence_list[i].hit[1].r) ;
      hvoxel->Fill(r[0],r[1],r[2]);
      // std::cout << i << " r="<< r[0]<<" " << r[1] << " " << r[2]  << std::endl;
    }
  }
  return 0;
}

int hfill(unsigned int *x, int len, int print, int ix, int iy){

  int recid,  reclen, i,j, edge, ch, rg, val;
  unsigned int data;
  
  for (int rg=0;rg<2;rg++) for (int ch=0;ch<MAXCH;ch++){
   valtdc[rg][ch]=-1111111;
   aadc[rg][ch]=0;
  }
  
  // trgigger timing
  int trgch=16;
  i=0;
  
  
  while (i<len){
    
    recid =x[i++];
    reclen=x[i++];
    if (print) printf("recid 0x%08x reclen=%d\n",recid, reclen);
    /*
    for (j=0;j<reclen;j++){
    
     if (i<len)    data=x[i++]; else break;
     if (recid==0x350000) && (((data>>27)&0x1F)==0)) {
	    edge  = (data >> 26 ) & 0x1;
	    ch    = (data >> 21 ) & 0x1F;
	    val   =   data & 0x1FFFFF; 
	    if (ch==trgch && edge==1){ ttrg=val; j= reclen;}
	 }
	} 
      
    */
      
      
    for (j=0;j<reclen;j++){
      if (i<len)    data=x[i++];     
      else break;
      
      switch (recid) {
      case 12:
	
	val =   data &0xFFFF;		       
        ch  =  (data  >> 17 ) &0x1F;               	  	
        edge = (data  >> 16 ) &0x1;
      //  htdc[edge][ch]->Fill(val);
	if (print) printf("%d [%d] %d  data=%d\n",j,ch,edge,val);
	break;
      case 0x350000:// 1290 TDC
   
	if (((data>>27)&0x1F)==0){
	    edge  = (data >> 26 ) & 0x1;
	    ch    = (data >> 21 ) & 0x1F;
	    val   =   data & 0x1FFFFF;
	    valtdc[edge][ch]=val;  
	   //  for (int krg=0;krg<1;krg++)  if (h2dadctdc[krg][ch]) h2dadctdc[krg][ch]->Fill(aadc[krg][ch],val/40.); 
	     //printf("%d [%d] %d  data=%d\n",j,ch,edge,val);
	}
	break;	
      case 0x340000:
      case 0x550000:
	
	edge=(data>>25)&0x3;
	if (edge==0) {
	  ch=(data>>17)&0xf;
	  if (recid==0x550000) ch+=16;
	  rg=(data>>16)&0x1;
	  val=data&0xfff;
	  if (print)  printf("%d [ch=%d] rg=%d  (ix=%d  iy=%d ) data=%d\n",j,ch,rg,ix, iy,val); 
        hadc[rg][ch]->Fill(val);   // Default action - if things go bad, go back to this option
        if (h2dadcdac[rg][ch]) h2dadcdac[rg][ch]->Fill(val,dac.dac);
          if (aadc[rg][ch]==0) aadc[rg][ch]=val;
	  
	}
	
	break;
      default:
	printf("Error in buffer\n");
	
      }
      
    } 
    
  }
  
 float vtdc,e1,e2;
  for (int rg=0;rg<2;rg++){
   for (int i=0;i<MAXCH;i++) {
          vtdc=valtdc[rg][i]-valtdc[rg][trgch];
          //vtdc=valtdc[rg][i];
          //vtdc/=40.;
          //printf("tdc:%f %d %d\n",vtdc, valtdc[rg][i],valtdc[rg][trgch]); 
          if (htdc[rg][i]) { 
            htdc[rg][i]->Fill(vtdc);// TDC   
            h2dadctdc[rg][i]-> Fill( aadc[rg][10],vtdc);
          }
    }  
  }
  return 0;
}  	 


int SelectionFunc(unsigned int *x, int len, int ich,float xmin, float xmax){
if (ich<0) return 1;
int recid,  reclen, i,j, edge, ch, rg, val;
 unsigned int data;
 int istdc=0;
 int which=ich/1000;
 ich %= 1000;
 
 
 i=0;
 while (i<len){
   
   recid =x[i++];
   reclen=x[i++];
   
   for (j=0;j<reclen;j++){
     if (i<len)    data=x[i++];     
     else break;
     
     switch (recid) {
     case 12:
       istdc=1; 
       val =   data &0xFFFF;		       
       ch  =  (data  >> 17 ) &0x1F;               	  	
       edge = (data  >> 16 ) &0x1;
       break;
      case 0x350000:// 1290 TDC
	break;	
     
     case 0x340000:
     case 0x550000:
       istdc=0;
       edge=(data>>25)&0x3;
       if (edge==0) {
	 ch=(data>>17)&0xf;
	 if (recid==0x550000) ch+=16;
	 rg=(data>>16)&0x1;
	 val=data&0xfff;
	 edge=rg;
	 valadc[rg][ch] = val;
       }
       
       break;
     default:
       printf("Error in buffer\n");
       
     }
     
     if (ich==ch+100*edge && which==istdc && val>=xmin && val<=xmax ) return 1;
     
     
   } 
   
 }
 
 
 return 0;
}  	 



int histo_init(){
int rg, ch;
char hn[128];
 hvoxel= new TH3F("hvoxel","hvoxel",100,-25,25,100,-25,25,100,0,200);

 std::ifstream myfile( "channels.map", ios::in  );
 float x,y,z; 
 
 
 if (myfile.is_open())
   {
     while (! myfile.eof() )
       {
	 myfile >> x >> y >> z;
	 Hit_Position.push_back(TVector3(x,y,z));
	 std::cout << x <<" " <<y <<" " << z<< std::endl;
       }
     myfile.close();
   }
 
 




 for (rg=0;rg<2;rg++){
    for (ch=0;ch<MAXCH;ch++){
       sprintf(hn,"hadc%d%02d",rg,ch);
       hadc[rg][ch]= new TH1F(hn,hn,ADCRANGE,-0.5,ADCRANGE-0.5);   

       sprintf(hn,"henergy%d%02d",rg,ch);
       henergy[rg][ch]= new TH1F(hn,hn,300,0,1300);   

       sprintf(hn,"hadcmax%d%02d",rg,ch);
       hadcmax[rg][ch]= new TH1F(hn,hn,ADCRANGE,-0.5,ADCRANGE-0.5);   

       sprintf(hn,"htdc%d%02d",rg,ch);
       htdc[rg][ch]= new TH1F(hn,hn,TDCRANGE,-0.5-TDCRANGE/2,TDCRANGE/2-0.5); 
       printf("%s\n",hn);
       
       sprintf(hn,"h2dadctdc%d%02d",rg,ch);
       int adcrange=3000;
       int tdcrange=2000;
       h2dadctdc[rg][ch]= new TH2F(hn,hn,adcrange,-0.5,adcrange-0.5, tdcrange,-tdcrange/2,tdcrange/2); 
       printf("%s\n",hn);
       
       if (!rg && ch<16){
	     sprintf(hn,"h2dadcadc%d%02d",rg,ch);
	     adcrange=1000;
	 
	     h2dadcadc[ch]= new TH2F(hn,hn,100,0,adcrange, 100,0,1); 
         printf("%s\n",hn);
       }
       
    }
  }
  return 0;
 }  
    	 
  	 
int mppc(char *fname, int firsteve=0, int numeve=-1, int ch=-1, float xmin=1500, float xmax=2500 ){ 	 



RUNREC *runrec;
ENDREC *endrec;
POSREC *posrec;
POSREC pos;
EVTREC *evtrec;
DACREC *dacrec;


#define BSIZE 10000
  unsigned int data[BSIZE];
  unsigned int *pdata;
  gzFile fp = gzopen(fname,"rb");	
  if (!fp) return 0;
  unsigned int hdr[4]={1};//recid
  int ntotal=0;
  int nb=0;
  int neve=0;
  histo_init();
 
  readcalib("linearity/calibration_constants.txt",16,cal.dx,1);
  readcalib("linearity/calibration_constants.txt",16,cal.n0,2);
  readcalib("linearity/calibration_constants.txt",16,cal.sx,3);
    while (!gzeof(fp)){
    nb=gzread(fp,hdr,   sizeof(int)*2 );
      if (gzeof(fp)) {
         printf("EOF %s nb=%d\n",fname,nb);  
	     break;
      }   
      nb = hdr[1]-2*sizeof(int);
      neve++;

      int recid=hdr[0];
      int len = hdr[1];
      if ((neve % 10000) == 0 || (neve<10)) {
	printf("Event %d %d  nb=%d %d\n",neve,(int) hdr[1],nb,nb/4);
	printf("recid=%d len=%d\n",recid,len);
      }
      ntotal = gzread(fp,data, nb);
      if (neve<firsteve) continue;
      if (numeve>0 && neve>(firsteve+numeve)) break;
      
      if (ntotal==nb){


	switch (recid){
	  
	case  RUNREC_ID:
	  runrec = (RUNREC *)&data[-2];
          printf("nev %d nch %d ped %d xy %d\n",(int) runrec->nev, (int) runrec->nch, (int) runrec->ped, (int) runrec->xy);
	  printf("nx= %d x0=%d dx=%d\n", (int) runrec->nx,(int) runrec->x0,(int) runrec->dx);
	  printf("ny= %d y0=%d dy=%d\n", (int ) runrec->ny, (int) runrec->y0, (int) runrec->dy);
	  for (int rg=0;rg<2;rg++){
	    for (int ch=0;ch<MAXCH;ch++){
              char hn[25];
              int nx=runrec->nx;
              int ny=runrec->ny;
              float minx =  runrec->x0-0.5*runrec->dx;
              float maxx =  runrec->x0+(nx-0.5)*runrec->dx;
              float miny =  runrec->y0-0.5*runrec->dy;
              float maxy =  runrec->y0+(ny-0.5)*runrec->dy;
              
	      sprintf(hn,"h2d%d%02d",rg,ch);
	      h2d[rg][ch] = new TH2F(hn,hn, nx,minx, maxx,ny,miny, maxy );
	      sprintf(hn,"h2dmax%d%02d",rg,ch);
	      h2dmax[rg][ch] = new TH2F(hn,hn, nx,minx, maxx,ny,miny, maxy  );



	      sprintf(hn,"h2dtdc%d%02d",rg,ch);
	      h2dtdc[rg][ch] = new TH2F(hn,hn, 3000,-3000, 0,ny,miny, maxy  );
	      
          sprintf(hn,"h2dposadcfrac%d%02d",rg,ch);
	      if (!rg && ch<16) h2dposadcfrac[ch] = new TH2F(hn,hn, nx,minx, maxx,50,0,1  );

		  sprintf(hn,"h2dposadc%d%02d",rg,ch);
	      if (!rg && ch<16) h2dposadc[ch] = new TH2F(hn,hn, nx,minx, maxx,300,0,1300  );


	    }
	  }
	  break;
	case  ENDREC_ID:
	  endrec = (ENDREC *) &data[-2];
	  break;
	case  POSREC_ID:
	  posrec = (POSREC *) &data[-2];
	  printf("POSREC x=%d y=%d\n",posrec->ix,posrec->iy);
	  pos.ix=posrec->ix;
	  pos.iy=posrec->iy;
	  pos.x=posrec->x;
	  pos.y=posrec->y;
	  pos.xset=posrec->xset;
	  pos.yset=posrec->yset;

	  break;
	case  EVTREC_ID:
	  evtrec = (EVTREC *) &data[-2];
	  // fprintf(stderr,"EVTREC len=%d \n",evtrec->len);

	  pdata =  &data[2];
	  {
	  
	    int selected = SelectionFunc(pdata,nb/sizeof(int),ch, xmin, xmax);
	    
	    if (selected) hfill(pdata,nb/sizeof(int),0, pos.xset, pos.yset);
	  }
	  break;
	case  DACREC_ID: 
	  {
	    dacrec = (DACREC *) &data[-2];
		dac.dac=dacrec->dac;
		dac.dac0=dacrec->dac0;
		dac.dac1=dacrec->dac1;
		dac.ddac=dacrec->ddac;
		
		

	    printf("DACREC len=%d dac=%f\n",dacrec->len, dacrec->dac);
	    
	    for (int rg=0;rg<2;rg++){
	      for (int ch=0;ch<MAXCH;ch++){
		char hn[25];
		
		int ny=int( (dac.dac1-dac.dac0)/dac.ddac ) + 1;
		float miny =  dac.dac0-dac.ddac/2.;
		float maxy =  dac.dac1+dac.ddac/2.;
	      
		
		sprintf(hn,"h2dadcdac%d%02d",rg,ch);
		if (!h2dadcdac[rg][ch] && ny>0)
		  h2dadcdac[rg][ch] = new TH2F(hn,hn,ADCRANGE,-0.5,ADCRANGE-0.5,ny,miny, maxy );
	      }
	    }
	  }
	  break;
	default:
          printf("Unknown recid!\n");
	  break;
	}
      } else {
        printf("Event %d nread=%d nrequested=%d\n",neve,ntotal,nb);    
        break;
      }  
  } 
 
  printf("Number of Events: %d\n",neve);
  

  gzclose(fp);  
 
 
  return 0;
}



#include "TROOT.h"
#include "TStyle.h"
#include "TF1.h"
#include "TDirectory.h"
#include "TCanvas.h"
#include "TH2F.h"

int setplain(){
gROOT->SetStyle("Plain");

gStyle->SetPalette(1, 0);
gStyle->SetOptStat(0);
gStyle->SetOptFit(111);


gStyle->SetPaperSize(TStyle::kA4) ;
gStyle->SetStatBorderSize(1);
gStyle->SetFrameBorderMode(0);
gStyle->SetCanvasBorderMode(0);
gStyle->SetPadBorderMode(0);
gStyle->SetPadColor(0);
gStyle->SetCanvasColor(0);
gStyle->SetStatColor(0);
return 0;
}

const int ind[16]={0,1,2,3,7,6,5,4,8,9,10,11,15,14,13,12};
int plothisto(char *name,int min, int max,int logy=1, char *opt="", char *fmt="%s%d", float xmin=-1, float xmax=-1){
setplain();

TCanvas *c2 = (TCanvas*)gROOT->GetListOfCanvases()->FindObject("cres");
if(!c2)  {
 c2= new TCanvas("cx","cx",900,900);
 c2->Connect("ProcessedEvent(Int_t,Int_t,Int_t, TObject *)", 0, 0 ,"ExtractHisto(Int_t, Int_t, Int_t,TObject *)");
}
int npads=max-min;
if (!npads) return 0;

int nrows= int(TMath::Sqrt(npads));
int ncol = int(ceil(float(npads)/nrows));

if (npads>0) c2->Divide(nrows,ncol);
TDirectory * dir=gDirectory;
printf("npads=%d ncols=%d nrows=%d %s\n",npads,nrows,ncol, dir->GetName());

char mname[128];


for (int i=0;i<npads;i++){
 
  if (npads==16){
   c2->cd(ind[i]+1)->SetLogy(logy);
  } else {
   c2->cd(i+1)->SetLogy(logy);
  }
  
 sprintf(mname,fmt,name,min+i);
 TH1F *h1= 0;
 //gDirectory->GetObject(mname,h1);
 h1 = ((TH1F *) dir->Get(mname));
 //h1->SetMinimum(-1);
 if (xmin!=xmax) h1->GetXaxis()->SetRange(xmin,xmax);
 if (h1!=NULL) h1->Draw(opt);
 else printf("not found ");
 
 
 printf("%s\n",mname);
}
c2->Modified();
c2->Update();

/*TCanvas *c3 = new TCanvas("c3","",900,900);
c3->cd();
TH1F *htotal = new TH1F("h1","",ADCRANGE,-0.5,ADCRANGE-0.5);

sprintf(mname,fmt,name,1);htotal->Add( ((TH1F *) dir->Get(mname)) );
sprintf(mname,fmt,name,4);htotal->Add( ((TH1F *) dir->Get(mname)) );
sprintf(mname,fmt,name,5);htotal->Add( ((TH1F *) dir->Get(mname)) );
sprintf(mname,fmt,name,6);htotal->Add( ((TH1F *) dir->Get(mname)) );
sprintf(mname,fmt,name,9);htotal->Add( ((TH1F *) dir->Get(mname)) );
htotal->Draw();*/

return 0;
}


TObject *FindClass(TPad *c, char * classname){
         printf("Searching for Object of class name %s in %s!\n",classname, c->GetName());
        TObject *obj=0;
        TIter    next(c->GetListOfPrimitives());
        while ((obj = next())) {
   	      printf("%s\n",obj->GetName());
   	      if (obj->InheritsFrom("TPad")) {
   	         printf("SubPad\n");
   	         TObject *o = FindClass( (TPad *) obj, classname );
   	         if (o) return o; 
   	      }
	      if (obj->InheritsFrom(classname))   return obj;	   	  
        }
        printf("Object of class name %s not found in %s!\n",classname, c->GetName());
        return 0;
}

TPad *FindPad(TPad *c, TObject *x){
        TObject *obj=0;
        TIter    next(c->GetListOfPrimitives());
        while ((obj = next())) {
          if (obj==x)   return c;	  
   	      printf("%s\n",obj->GetName());
   	      if (obj->InheritsFrom("TPad")) {
   	         printf("SubPad\n");
   	         TPad *o = FindPad( (TPad *) obj, x );
   	         if (o) return o; 
   	      }
        }
        return 0;
}


void  ExtractHisto(Int_t event, Int_t ex, Int_t ey,TObject *selected){
    const int print=0;
    TCanvas *c = (TCanvas *) gTQSender; 
        
    if (print) printf("Canvas %s: event=%d, x=%d, y=%d, selected=%s from TPad=%d\n", c->GetName(),
          event, ex, ey, selected->IsA()->GetName(),selected->InheritsFrom("TPad"));

    if (event != kButton1Down)      return ;
    if(!selected) return ;   
  
    TPad *mpad = FindPad( (TPad *) c , selected); 
    if (!mpad) return;
    
    if (print) printf("Found Pad %s\n",mpad->GetName());
    TObject *obj = FindClass( mpad ,"TH1");
    
    if (!obj){
	   printf("%s\n",obj->GetName());
	   c->SetUniqueID(0);
	   return ;
    }
    if (print) printf("Found Object %s\n",obj->GetName());
    gPad->GetCanvas()->FeedbackMode(kTRUE);
 
    TCanvas *c2 = (TCanvas*)gROOT->GetListOfCanvases()->FindObject("cres");
    if (!c2) c2= new TCanvas("cres","cres",500,400);
    
    if (obj) {
    
      TH1 *h = (TH1*)obj;
      c2->cd();
      h->Draw();
      c2->Modified();
      c2->Update();
    }
    return;
}



#ifdef MAIN
int main(int argc, char ** argv){

ana(argv[1]);

return 0;
}

#endif