#include <stdlib.h>
#include <stdio.h>
#include <iostream>
#include <sys/types.h>
#include <sys/stat.h>

#include "CLHEP/Vector/ThreeVector.h"
#include <math.h>

float rint(float a){
  if ( a>0 ) return int ( a + 0.5 ); 
  else       return int ( a - 0.5 ); 

}

// *********************************************************************
//                         koda za izpis daq1 fileov
// *********************************************************************
//  ce hoces MC tracke in hite iz denisove rutine 
//  nastavi v kumacu:  setenv WRITEMODE 1
//  sicer               setenv WRITEMODE 0
 
//  rezultati se pisejo v $HBSCRATCH/...

//  v usevnt: daq1_write();
//  v usstop: daq1_stop();

//  polnjenje RICHa: :
 int AddHit(int id);
 int  FillfromRICH( Hep3Vector r, 
 	      float txf,float tyf, 
  	      float momentum,int lund,float sig_tx=0, float sig_ty=0); 
// *********************************************************************
// *********************************************************************
#define OK 0
#define ERROR 1
#define MAXHIT 1024

#include "arte/arte.hh"
#include "farm/farm.hh"
#include "arte/RSEG_cmp.hh"
#include "rich/rich.h"
#include "rise/RICHEventMC.h"  
#define MAXFED 28

struct FED_record {
  unsigned short bc,fedid,recid,bxn;
  unsigned char data[128];
};
typedef struct FED_record;

FED_record fedr[MAXFED];                                 


struct LOG_HDR {
  unsigned long  bc,exp,run,evt,dtm;
}; 
typedef struct LOG_HDR LOG_HDR;

struct EVT_trailer
{
  unsigned long slt;
};
typedef struct EVT_trailer EVT_trailer;

struct MC_INFO {
  short lund;
  short mother;
};
typedef struct  MC_INFO   MC_INFO;  

struct ECAL_HIT {
  short tx;
  short ty;
  short sigma;
  unsigned short E;
  short sig_E;
  short dL;     
};
typedef struct ECAL_HIT ECAL_HIT;

struct OTR_TRACK {
  short x;     // x*100/cm at z=1000 cm
  short y;     // y*100/cm at z=1000 cm
  short sig_x; // sig_x*100/cm at z=1000 cm
  short sig_y; // sig_x*100/cm at z=1000 cm
  short tx;    // tx*100000 at z=1000 cm
  short ty;    // ty*100000 at z=1000 cm
  short sig_tx;   // sig_tx*100000 at z=1000 cm
  short sig_ty;   // sig_ty*100000 at z=1000 cm
  unsigned short pf;   // 1/p*10000 GeV     
  short sig_pf;        // sig(1/p)*100000 GeV 
};
typedef struct OTR_TRACK OTR_TRACK;  
// *********************************************************************


LOG_HDR  logh;
EVENT_HDR evthdr;
EVENT_TAIL evttail;
EVT_DATA_RECORD evt_rec;
SLB_HDR   ecal_SLB_HDR;
ECAL_HIT  ecal_hit[MAXHIT];
SLB_TAIL  ecal_SLB_TAIL;
OTR_TRACK track[MAXHIT];
SLB_HDR   otr_SLB_HDR;
SLB_TAIL  otr_SLB_TAIL;
SLB_HDR   rich_SLB_HDR;
ETB_fed *etb_fed =new ETB_fed[MAXFED];
SLB_TAIL   rich_SLB_TAIL;

MC_INFO  otr_mc[MAXHIT];
MC_INFO  ecal_mc[MAXHIT];
// *********************************************************************

RICHDAQ *daq;

static float zmiddle=1000; // middle RICH
int myeve=0;    // counter of output events
int myifile=1;    // counter of output files 
int allitra=0;
int allotrtra=0;
int allpc=0;
int alltc=0;

int itra_OTR=0;
int itra=0;

FILE *mf;
int numfedevt=0;
float Ecut=0;  // energy cut for cluster selection (GeV)
float Mcut=0;  // invariant mass cut for event selection (GeV)
float Ecutall=0;  // energy cut for cluster selection (GeV)
float Mcutall=0;  // invariant mass cut for event selection (GeV)
int Mode=0;       // mode za izpis 
      char myname[80];
// *********************************************************************

// ----------------------------SWAPING-----------------------------
#define swap_4(x) swap_long((long *) x)
#define swap_2(x) swap_short((short *) x)

int isLE=0;
int is_little_endian()
{
  int test=1;
  return (int) *((char *) &test);
}

int endianinit()
{
  if ((isLE=is_little_endian()))
    printf("Running on LittleEndian machine\n");
   else
    printf("Running on BigEndian machine\n");

  return 0;
}

void swap_short(short *shtosw)
{
  short stmp;
  char *b1 = (char *) shtosw, *b2 = (char *) &stmp;
  
  stmp=*shtosw;
  b1[0]=b2[1];
  b1[1]=b2[0];
  return;
}

void swap_long(long *lotosw)
{
  long ltmp;
  char *b1 = (char *) lotosw, *b2 = (char *) &ltmp;
  
  ltmp=*lotosw;
  b1[0]=b2[3];
  b1[1]=b2[2];
  b1[2]=b2[1];
  b1[3]=b2[0];
  return;
}

void swap( unsigned char *data) {
  for (int j=0;j<64;j++) {
    char tmp=data[2*j];
    data[2*j]=data[2*j+1];
    data[2*j+1]=tmp;      
  }
  return;
}

void swap(EVENT_TAIL  &a) { 
  swap_4(&a.sltevt);
  return;
}

void swap(SLB_TAIL &a ){
  swap_2(&a.recid);
  swap_2(&a.fltevt);
  return ;
}

void swap(MC_INFO &a ){
  swap_2(&a.lund);
  swap_2(&a.mother);
  return ;
}
void swap(OTR_TRACK &otrdata) {
  swap_2(&otrdata.x);
  swap_2(&otrdata.y);
  swap_2(&otrdata.sig_x);
  swap_2(&otrdata.sig_y);
  swap_2(&otrdata.tx);
  swap_2(&otrdata.ty);
  swap_2(&otrdata.sig_tx);
  swap_2(&otrdata.sig_ty);
  swap_2(&otrdata.pf);
  swap_2(&otrdata.sig_pf);
  return;
}
void swap(ECAL_HIT &ecaldata) {
  swap_2(&ecaldata.tx);
  swap_2(&ecaldata.ty);
  swap_2(&ecaldata.sigma);
  swap_2(&ecaldata.E);
  swap_2(&ecaldata.sig_E);
  swap_2(&ecaldata.dL);
  return;
}

void swap(LOG_HDR    &a) { 
  swap_4(&a.bc);
  swap_4(&a.exp);
  swap_4(&a.run);
  swap_4(&a.evt);
  swap_4(&a.dtm);
  return;
}
void swap(EVENT_HDR  &evh) {
  swap_4(&evh.length);
  swap_2(&evh.recid);
  swap_2(&evh.format);
  swap_4(&evh.sltevt);
  swap_4(&evh.compmask);
  swap_4(&evh.evtmask);
  return;
}


void swap(EVT_DATA  &evdata) {
  swap_4(&evdata.run_number);
  swap_4(&evdata.event_time);
  swap_4(&evdata.event_status);
  swap_4(&evdata.flt_status);
  swap_4(&evdata.slp_node);
  swap_4(&evdata.slt_status);
  swap_4(&evdata.exp_number);
  swap_4(&evdata.farm_node);
  swap_4(&evdata.BX_long);
  swap_2(&evdata.BX_upper);
  swap_2(&evdata.FCS_tcode);
  swap_2(&evdata.physBX);
  swap_2(&evdata.FLT_BX);
  swap_2(&evdata.fltevt);
  return;
}


void swap(SLB_HDR &slbh ){
  swap_4(&slbh.length);
  swap_2(&slbh.recid);
  swap_2(&slbh.fltevt);
  swap_2(&slbh.compid);
  swap_2(&slbh.slbid);
  swap_4(&slbh.nodeid);
  return ;
}

// ------------------------building fed record rich --------------

inline Hep3Vector toHep3Vector(const Vec3 &x) {
  return Hep3Vector(x.x(), x.y(), x.z());
}

int  Filltracks_MIMP(){
  // DENIS's function
  RICHrings *rich=new RICHrings;
  // loop over RICH MC particles (= at least one photon is radiated)
  int nhits=0;
  for ( rich->begin(); rich->valid(); rich->next() ) {
    float tx=rich->tx();
    float ty=rich->ty();
//     float x=rich->x().x();
//      float y=rich->x().y();
//      float z=rich->x().z();
//      Vec3 rr(x,y,z);
    Vec3 entr=rich->x();
    //    cout << "tracks " << entr<< endl;
    if (rich->vertex()->z > 850 ) continue;
    int lund=rich->topp()->lund;
    //if (lund !=11) continue; // select only electrons
    if (lund !=211) continue; // select only pions
    float momentum = rich->momentum(); 
    if (momentum<15) continue; // cut  15 Gev
    FillfromRICH( 0.1*toHep3Vector(entr), tx, ty, momentum,lund,0,0); 

    RICHphotons *ph=rich->GetPhotons();
      
    for (ph->begin(); ph->valid(); ph->next()) { 
      if (rich->mtra()==ph->mtra()) {
	//ArteRelationIterator<HITB>::iterator ihitb= ph->mimp()->hitb.begin();
	if (!ph->detected()) continue;   
        int id=ph->id();
	if (id>0) {
	  AddHit(id);
	  nhits++;
	}
      } 
    }
    
    break; 
  }
  cout <<    itra_OTR << " MIMP  nhits => " <<  nhits <<endl;
  return itra;
}



int PrintInfo(int id){ 
  RICHDetector &rd=*RICH_Detector;
  rd.ChannelExist(id);
  int ci=rd.ChannelIndex();
  RICH_DAQ->ChannelExist(ci);
  cout << id << " fedid()        "<< RICH_DAQ->fedid() << endl; 
  cout << id << " connector()    "<< RICH_DAQ->connector() << endl; 
  cout << id << " daughter()     "<< RICH_DAQ->daughter() << endl; 
  cout << id << " channel()      "<< RICH_DAQ->channel() << endl; 
  cout << id << " bit()          "<< RICH_DAQ->bit() << endl; 
  cout << id << " status()       "<< RICH_DAQ->status() << endl;
  cout << id << " sm()           "<< RICH_DAQ->sm()  << endl;
  cout << id << " row()          "<< RICH_DAQ->row() << endl;
  cout << id << " col()          "<< RICH_DAQ->col()  << endl; 
  cout << id << " pm_channel()   "<< RICH_DAQ->pm_channel() << endl; 
  cout << id << " pm_numchannel()"<< RICH_DAQ->pm_numchannel() << endl;
  cout << id << " ci " << ci << endl;
  cout << id << "--------------------------------------"<<endl;
  return 0;
}

int write_fed_rich(){
  
  for (int nr=0;nr< MAXFED ; nr++) {
    fedr[nr].bc    = 136;
    fedr[nr].fedid = nr;
    fedr[nr].recid = 0x1002;
    fedr[nr].bxn   = 0;
    if (isLE) {      
      swap_2(&fedr[nr].bc);
      swap_2(&fedr[nr].fedid);
      swap_2(&fedr[nr].recid);
      swap_2(&fedr[nr].bxn);
    } else {
      //      swap(fedr[nr].data);
    }  
    if ( fwrite((void *) &fedr[nr], 1, 8, mf)!= 8 ) return 0;
    if ( fwrite((void *) fedr[nr].data, 1, 128, mf) != 128 ) return 0;
    // reset buffer
    for (int i=0;i<128;i++ )fedr[nr].data[i]=0;
  }
  return 0;
}
//-----------------------RPNT--------------------------------------
int FillfromRPNT( ArtePointer<RTRA> rtra)
{
  float z= 1000; // middle RICH
  register float zmin= 850; // begining of the rich
  //look for RPNT closer to the begin of RICH
  if (!(itra_OTR<MAXHIT)) return 0;

  ArtePointer<RPNT> closerpnt;
  for(ArteRelationIterator<RPNT> rpntit=rtra->rpnt.begin();
      rpntit!=rtra->rpnt.end();rpntit++)
    {
      if( ((*rpntit)->zf > closerpnt->zf) && 
          ((*rpntit)->zf < zmin) )	closerpnt=*rpntit;
    }
  if ( closerpnt->zf< 450 ) return 0;
  ArtePointer<MTRA> imtra(rtra->mtra.begin());
  ArtePointer<MTRA> moth_mtra=imtra;
  int eva = moth_mtra->topp->lund ;
  int bin=moth_mtra->topp->lund/100;
      bin=abs(bin%10);

  while (  (moth_mtra->topp->lund > 10 )&& ( bin!=5 )   ) {
    eva = moth_mtra->topp->lund ;
    bin=moth_mtra->topp->lund/100;
    bin=abs(bin%10);
    //cout << i++ << " topp " << moth_mtra->topp->lund << endl;
    moth_mtra= moth_mtra->btra;
  }
  //cout << i << " eva " << eva << endl;
  float x=(100*(closerpnt->xf + closerpnt->txf*(z- closerpnt->zf)) );
  if (fabs(x)>32000) return 0;
  track[itra_OTR].x = (short) x;
  
  float y=(100*(closerpnt->yf + closerpnt->tyf*(z- closerpnt->zf)));
  if (fabs(y)>32000) return 0;
  track[itra_OTR].y =(short) y;
  
  track[itra_OTR].sig_x    =(short)( 100*sqrt(closerpnt->cvf[0]));      
  track[itra_OTR].sig_y    =(short)( 100*sqrt(closerpnt->cvf[2])); 
  
  float tx=(  100000*closerpnt->txf);
  if (fabs(tx)>32000) return 0;
  track[itra_OTR].tx = (short)tx;
  
  float ty=(  100000*closerpnt->tyf);  
  if (fabs(ty)>32000) return 0;
  track[itra_OTR].ty = (short) ty;
 
  
  track[itra_OTR].sig_tx   =(short)( 100000*sqrt(closerpnt->cvf[5]));
  track[itra_OTR].sig_ty   =(short)( 100000*sqrt(closerpnt->cvf[9])); 
  
  float pf = rint(10000.*fabs(closerpnt->pf));
  track[itra_OTR].pf = (short unsigned) pf; 
  track[itra_OTR].sig_pf    =(short) rint( 100000*sqrt(closerpnt->cvf[14]));     

  otr_mc[itra_OTR].lund=imtra->topp->lund;
  otr_mc[itra_OTR].mother =moth_mtra->topp->lund  ;
  
  itra_OTR++;
  return 1;
}

//------------------------------------------------------------

int FillfromRICH( Hep3Vector r, float txf,float tyf, float momentum,int lund,float sig_tx, float sig_ty)
{
  // ARTE r units [cm]
  float z= 1000; // middle RICH
  if (!(itra_OTR<MAXHIT)) return 0;  

  float x=100*(r.x() + txf*(z- r.z() ) );
  if (fabs(x)>32000) return 0;
  track[itra_OTR].x = (short) x;
  
  float y=100*(r.y() + tyf*(z- r.z() ) );
  if (fabs(y)>32000) return 0;
  track[itra_OTR].y =(short) y;
  
  track[itra_OTR].sig_x    =0;
  track[itra_OTR].sig_y    =0;
  
  float tx=(  100000*txf);

  if (fabs(tx)>32000) return 0;
  track[itra_OTR].tx = (short)tx;
  
  float ty=(  100000*tyf);  
  

  if (fabs(ty)>32000) return 0;
  track[itra_OTR].ty = (short) ty;
 
  
  track[itra_OTR].sig_tx   =(short)( 100000*sqrt(sig_tx));
  track[itra_OTR].sig_ty   =(short)( 100000*sqrt(sig_ty)); 
  
  float pf=0;
  if (momentum !=0) pf =1/momentum;
  track[itra_OTR].pf= (short unsigned) rint(10000*fabs(pf)); 
  track[itra_OTR].sig_pf    =0;     

  otr_mc[itra_OTR].lund=lund;
  otr_mc[itra_OTR].mother =0  ;
  cout << "TRACK " << r << r-(r.z()-860)*Hep3Vector(txf,tyf,1) << endl;
  printf(" > TRACK :x %d y %d tx %d ty %d \n",
			     track[itra_OTR].x,
			     track[itra_OTR].y,
			     track[itra_OTR].tx,
			     track[itra_OTR].ty);  
			     
  itra_OTR++;
  return 1;
}

//------------------------------------------------------------
int clear_rich(){
  for (int i=0; i<MAXFED ;i++)
    for (int byte=0; byte<128 ;byte++)
      fedr[i].data[byte]=0;

return 0;
}



// add hit
int AddHit(int id){
  if (id<0) return 0;
  //      PrintInfo(id);
  if (!RICH_DAQ)      
  {
    //    printf("RICH_DAQ\n");
    RICH_DAQ=new RICHDAQ();
  }
  if (!RICH_Detector) {
    //  printf("RICH_Detector\n");
    RICH_Detector=new RICHDetector();
  }
  RICHDetector &rd=*RICH_Detector;
  RICHDAQ      &daq= *RICH_DAQ;  
  //  if (!rd.ChannelExist(id)) return 0;
  rd.ChannelExist(id);
  int ci=RICH_Detector->ChannelIndex();
  daq.ChannelExist(ci);
  int fedid = daq.fedid();
  int ch    = daq.channel();
  int con   = daq.connector();
  int dau   = daq.daughter() ;
  int byte  = ch/8 + 2* (con + dau * 16 );
  int bit   = ch%8;
  fedr[fedid].data[byte]=(fedr[fedid].data[byte] | (1<<bit) );
  return 0;
}


//------------------------------------------------------------

int mc2daq1(){
  for(ArteTable<HITB>::iterator hitbit=ArteTable<HITB>::begin();
      hitbit!=ArteTable<HITB>::end();++hitbit) {
    ArtePointer<HITB> hitb(hitbit);
    ArtePointer<GEDE> gede(hitb->gede); 
    if(gede->cmp==5){        // 5 = RICH 
      int id=(*hitb->digb.begin())->id;
      //      PrintInfo(id);
      AddHit(id);
    }
  }    

  return 0;
}
//-----------------------------------------------------------
int get_rccl() {
  // get ECAL clusters above Ecut
  itra=0;
  for(ArteTable<RCCL>::iterator rcclit=ArteTable<RCCL>::begin();
      rcclit!=ArteTable<RCCL>::end();++rcclit) {
    ArtePointer<RCCL> rccl(rcclit);
    if (rccl->esum <Ecut) continue;
    ecal_hit[itra].tx=(short) rint((rccl->x/rccl->z)*100000);
    ecal_hit[itra].ty=(short) rint((rccl->y/rccl->z)*100000);
    ecal_hit[itra].sigma=
      (short) rint(sqrt((rccl->cve[0]+rccl->cve[2])/2)/rccl->z*100000);
    ecal_hit[itra].E= (short) rint(100*rccl->esum);
    ecal_hit[itra].sig_E=(short) rint(100*sqrt(rccl->cve[14]));
    float dl= rccl->lemp-rccl->lhad;
    if(dl>300) {dl=300;}
    if(dl<-300) {dl=-300;}
    ecal_hit[itra].dL=(short) rint(100*dl);
    ArtePointer<MTRA> imtra(rccl->mtra.begin());
    ArtePointer<MTRA> moth_mtra(imtra->btra );
    ecal_mc[itra].lund= imtra->topp->lund;


    ecal_mc[itra].mother=moth_mtra->topp->lund;

    itra++;
    if(itra==MAXHIT) break;
  }
  return 0;
}
//-------------------------------------------------------------------
int findTC( ArtePointer<RSEG>  rseg ){
  int nseg=0;
  for(ArteRelationIterator<RSEG> crsegIt = rseg->rend.begin() ;
      crsegIt != rseg->rend.end() ; ++crsegIt) {
    if( !((*crsegIt)->cmp & Rsegc::bitsign) &&
	(*crsegIt)->cmp & Rsegc::trit ) {
      nseg++;
    }
  }
  return (nseg>0);
}
//-------------------------------------------------------------------
int findPC( ArtePointer<RSEG>  rseg ){
  int nseg=0;
  for(ArteRelationIterator<RSEG> crsegIt = rseg->rbeg.begin() ;
      crsegIt != rseg->rbeg.end() ; ++crsegIt) {
    if( !((*crsegIt)->cmp & Rsegc::bitsign) &&
	(*crsegIt)->cmp & Rsegc::patt ) {
      nseg++;
    }
  }
  return (nseg>0);
}
//-------------------------------------------------------------------
int get_otr(){
   // fill track structure **************************************
  itra_OTR=0;
  /*
  for(ArteTable<RTRA>::iterator rtra1it=ArteTable<RTRA>::begin();
      rtra1it!=ArteTable<RTRA>::end();++rtra1it) {
    ArtePointer<RTRA> rtra(rtra1it);
    FillfromRPNT(rtra);
  }
  printf( "get_otr()\n");
  return 0;
  */
  
  for(ArteTable<RSEG>::iterator rtra1it=ArteTable<RSEG>::begin();
      rtra1it!=ArteTable<RSEG>::end();++rtra1it) {
    
    ArtePointer<RSEG> rtra(rtra1it);
    // ************************************************************
    // track selection is done here
    // ************************************************************
    // before RICH
    if (( rtra->cmp & Rsegc::bitsign ) !=0) continue; 

    int pc= (rtra->cmp & ((Rsegc::wch )|(Rsegc::patt))) ;
    if (pc) allpc++;
    // between RICH and ECAL
    int tc= (rtra->cmp & ((Rsegc::wch )|(Rsegc::trit))) ;
    int istc = 0, ispc=0;

    if (pc) istc = findTC( rtra );
    if (tc) ispc = findPC( rtra );

    if (tc) alltc++;

    if (pc||tc) {
      float x=(100*(rtra->xf + rtra->txf*(zmiddle- rtra->zf)) );
      if (fabs(x)>32000) continue;
      track[itra_OTR].x = (short) x;
      
      float y=(100*(rtra->yf + rtra->tyf*(zmiddle- rtra->zf)));
      if (fabs(y)>32000) continue;
      track[itra_OTR].y =(short) y;
      
      track[itra_OTR].sig_x    =(short)( 100*sqrt(rtra->cvf[0]));      
      track[itra_OTR].sig_y    =(short)( 100*sqrt(rtra->cvf[2])); 
      
      
      float tx=(  100000*rtra->txf);
      if (fabs(tx)>32000) continue;
      track[itra_OTR].tx = (short)tx;
      
      float ty=(  100000*rtra->tyf);  
      if (fabs(ty)>32000) continue;
      track[itra_OTR].ty = (short) ty;
      
      
      track[itra_OTR].sig_tx   =(short)( 100000*sqrt(rtra->cvf[5]));
      track[itra_OTR].sig_ty   =(short)( 100000*sqrt(rtra->cvf[9])); 
      
      
      if (fabs(rtra->pf)>0) 
	track[itra_OTR].pf=(int) rint(10000*fabs(rtra->pf)); 
      else  
	track[itra_OTR].pf=0;
      //track[itra_OTR].sig_pf    =rint( 100000*sqrt(rtra->cvf[14]));           
      track[itra_OTR].sig_pf    =short (tc + (ispc << 1 )+ 
					(pc << 2 )+ (istc << 3) + rtra->nhit << 4);     
      
      ArtePointer<MTRA>   imtra(rtra->mtra.begin());
      ArtePointer<MTRA> moth_mtra(imtra->btra );

      otr_mc[itra_OTR].lund=imtra->topp->lund;
      otr_mc[itra_OTR].mother=moth_mtra->topp->lund;

      itra_OTR++;
      if(itra_OTR==MAXHIT) break;
      
    }
  }
  return 0;
}


// ----------------------------- writing data ----------------
int write_raw_data( int mode=0) {
  ArtePointer<EVHD> evhd( 1 );
  int mc= (evhd->run < 2);

  //  printf("mode=%d MC=%d\n",mode, mc);
  if (!mode) {
    get_rccl();
    get_otr();
    if (mc) mc2daq1(); else   ARTE_get_feds_DRIC(MAXFED,etb_fed,&numfedevt);
  } else {
    Filltracks_MIMP();
  }


  if (mc) {    numfedevt= MAXFED; }
  // fill Event_data_record
  allitra+=itra;
  allotrtra+=itra_OTR;
  // ************************************************************

  evt_rec.hdr.length = sizeof(evt_rec);
  evt_rec.hdr.recid  = 0x2000;
  evt_rec.hdr.fltevt = 0;   // where to get this?
  evt_rec.hdr.compid = 0x000;
  evt_rec.hdr.slbid  = 0;   // where to get this?
  evt_rec.hdr.nodeid = 0;   // where to get this?
  evt_rec.data.run_number   = faDDAQrun_number();   // Run number
  evt_rec.data.event_time   = faDDAQevent_time();   // Time event was built
  evt_rec.data.event_status = faDDAQevent_status(); // event status/error code
  evt_rec.data.flt_status   = faDDAQflt_status();   // FLT trigger number
  evt_rec.data.slp_node     = faDDAQslp_node();     // SLP processing node
  evt_rec.data.slt_status   = faDDAQslt_status();   // SLT trigger number/status
  evt_rec.data.exp_number   = faDDAQexp_number();   // Experiment Number
  evt_rec.data.farm_node    = faDDAQfarm_node();    // Farm node
  evt_rec.data.BX_long      = faDDAQBX_long();      // bits 0 to 31 of the BX tag
  evt_rec.data.BX_upper     = faDDAQBX_upper();     // bits 32 to 47 of the BX tag
  evt_rec.data.FCS_tcode    = faDDAQFCS_tcode();    // FCS trigger code for event
  evt_rec.data.physBX       = faDDAQphysBX();       // physical BX == hera bucket ID of event
  evt_rec.data.FLT_BX       = faDDAQFLT_BX();       // pipeline cell triggered on by FLT
  evt_rec.data.fltevt       = faDDAQfltevt();       // flt event number as delivered by FCS
  evt_rec.tail.recid  = evt_rec.hdr.recid;
  evt_rec.tail.fltevt = evt_rec.hdr.fltevt;
  
  // fill ECAL hits hdr, tail
  ecal_SLB_HDR.length = sizeof(SLB_HDR) + 
    sizeof(SLB_TAIL) + 
    (sizeof(ECAL_HIT)+sizeof(MC_INFO)*mc) * itra;
  ecal_SLB_HDR.recid = 0x7500+ mc;
  ecal_SLB_HDR.fltevt = (sizeof(ECAL_HIT)+sizeof(MC_INFO)*mc);  
  ecal_SLB_HDR.compid = 0x0700;
  ecal_SLB_HDR.slbid = (short) rint(Ecut*100.);   // cut on cluster energy
  ecal_SLB_HDR.nodeid = (short) rint(Mcut*100.);   // cut on invariant mass
  ecal_SLB_TAIL.recid = ecal_SLB_HDR.recid;
  ecal_SLB_TAIL.fltevt = ecal_SLB_HDR.fltevt;

  // fill OTR hits hdr, tail
  otr_SLB_HDR.length = sizeof(SLB_HDR) + sizeof(SLB_TAIL) 
    + ( sizeof(OTR_TRACK)+sizeof(MC_INFO)*mc ) * itra_OTR;
  otr_SLB_HDR.recid = 0x4500 +mc;
  otr_SLB_HDR.fltevt =( sizeof(OTR_TRACK)+sizeof(MC_INFO)*mc );   // where to get this?
  otr_SLB_HDR.compid = 0x0400;
  otr_SLB_HDR.slbid =   (short) rint(zmiddle*10.);   //ce z v [cm]  
  otr_SLB_HDR.nodeid = 0;
  otr_SLB_TAIL.recid = otr_SLB_HDR.recid;
  otr_SLB_TAIL.fltevt = otr_SLB_HDR.fltevt;
  
  // fill RICH data hdr, tail
  int siz = sizeof(etb_fed[0].hdr) + 128; // + sizeof(etb_fed[0].data);
  rich_SLB_HDR.length = sizeof(SLB_HDR) + sizeof(SLB_TAIL) + siz * numfedevt;
  rich_SLB_HDR.recid = 0x2000;
  rich_SLB_HDR.fltevt = 0;   // where to get this?
  rich_SLB_HDR.compid = 0x0500;
  rich_SLB_HDR.slbid = 0;   // where to get this?
  rich_SLB_HDR.nodeid = numfedevt;   // where to get this?
  rich_SLB_TAIL.recid = rich_SLB_HDR.recid;
  rich_SLB_TAIL.fltevt = rich_SLB_HDR.fltevt;
  
  // fill event hdr, tail
  evthdr.length = sizeof(evthdr) + sizeof(evttail) 
    + evt_rec.hdr.length + 
    ecal_SLB_HDR.length + 
    rich_SLB_HDR.length +
    otr_SLB_HDR.length;
  evthdr.recid = 0x3000;
  evthdr.format = 0;   // where to get this?
  ArtePointer<DIHD> dihd( 1 );
  evthdr.sltevt=dihd->slt;
  evthdr.compmask = 0;   // where to get this?
  evthdr.evtmask = 0;   // where to get this?
  evttail.sltevt = evthdr.sltevt;
  
  // fill log header
  logh.bc= sizeof(logh) + evthdr.length;
  logh.exp= evhd->exp;
  logh.run= evhd->run;
  logh.evt= evhd->evt;
  logh.dtm= evhd->date;
  // Endian inversion 
  if (isLE) {
    swap(logh);
    swap(evthdr);
    swap(evt_rec.hdr); 
    swap(evt_rec.data); 

    swap(ecal_SLB_HDR);
    for(int i=0; i<itra; i++)  { 
      swap(ecal_hit[i]);   
      if (mc) swap(ecal_mc[i]);
    }
    swap(ecal_SLB_TAIL);

    swap(otr_SLB_HDR);
    for(int i=0; i<itra_OTR; i++)  {
      swap(track[i]);
      if (mc) swap(otr_mc[i]);
    }   
    swap(otr_SLB_TAIL);

    swap(rich_SLB_HDR);
    swap(rich_SLB_TAIL);

    swap(evttail);
  }  
   
  // write to file
  // Event **********************************************
  int ib=fwrite((void *) &logh,sizeof(LOG_HDR),1,mf);
  ib=fwrite((void *) &evthdr,sizeof(EVENT_HDR),1,mf);
  ib=fwrite((void *) &evt_rec,sizeof(evt_rec),1,mf);
  
  // ECAL **********************************************
  ib=fwrite((void *) &ecal_SLB_HDR,sizeof(SLB_HDR),1,mf);
  for(int i=0; i<itra; i++) {
    ib=fwrite((void *) &ecal_hit[i],sizeof(ECAL_HIT),1,mf);
  if (mc)  ib=fwrite((void *) &ecal_mc[i],sizeof(MC_INFO),1,mf);
  }
  ib=fwrite((void *) &ecal_SLB_TAIL,sizeof(SLB_TAIL),1,mf);
  
  // OTR *************************************************
  ib=fwrite((void *) &otr_SLB_HDR,sizeof(SLB_HDR),1,mf);
  for(int i=0; i<itra_OTR; i++) {
    ib=fwrite((void *) &track[i],sizeof(OTR_TRACK),1,mf);
    if (mc)  ib=fwrite((void *) &otr_mc[i],sizeof(MC_INFO),1,mf);
  }
  ib=fwrite((void *) &otr_SLB_TAIL,sizeof(SLB_TAIL),1,mf);
  
  // RICH   ***********************************************
  ib=fwrite((void *) &rich_SLB_HDR,sizeof(SLB_HDR),1,mf);
  if (mc) {
    write_fed_rich();
  } else { 
    for (int i=0;i<numfedevt;i++) {
      ib=fwrite((void *) &(etb_fed[i].hdr) ,sizeof(etb_fed[i].hdr),1,mf);
      ib=fwrite((void *) etb_fed[i].data ,128,1,mf); // sizeof(data)
    }
  }
  ib=fwrite((void *) &rich_SLB_TAIL,sizeof(SLB_TAIL),1,mf);
  // Event **********************************************
  ib=fwrite((void *) &evttail,sizeof(EVENT_TAIL),1,mf);
  
  if (myeve%1000==0) {
    struct stat *buf = new struct stat;
    stat(myname,buf);
    printf("File size:%d EVE %d>number of ECAL CLUSTERS:%d\n",
	   buf->st_size,myeve,itra);
    if (buf->st_size > 200000000) {
      fclose(mf);
      myifile++;
      char* daq1name= getenv("RECODAQ1");
      if (daq1name)
	sprintf(myname,"%s_file%.3d.daq1",daq1name,myifile); 
      else {
	
	char* path=getenv("HBSCRATCH"); 
	ArtePointer<EVHD> evhd(1);
	if (evhd->run>1) 
	  sprintf(myname,"%s/run%.5d_file%.3d.daq1",path,evhd->run,myifile);
	else
	  sprintf(myname,"%s/run_mc_file%.3d.daq1",path,myifile); 
      }
      mf=fopen(myname,"w");
      if (mf==NULL) {
	cout << "ERROR fopen" << endl;
	exit(1);
      }
      cout << " FILE fopen :" << myname << endl;
    }
  }                 
  itra=0;
  itra_OTR=0;
  clear_rich();
  return 0;
}



int  evnt_select(){
  float mcut2=Mcut*Mcut;
  for(ArteTable<RCCL>::iterator rcclit=ArteTable<RCCL>::begin();
      rcclit!=ArteTable<RCCL>::end();++rcclit) {
    ArtePointer<RCCL> rccl(rcclit);
    if (rccl->esum <Ecut) continue;
    for(ArteTable<RCCL>::iterator rcclit1=rcclit;
	rcclit1!=ArteTable<RCCL>::end();++rcclit1) {
      ArtePointer<RCCL> rccl1(rcclit1);
      if (rccl1->esum<Ecut) continue;
      float invmass=(rccl->esum+rccl1->esum)*(rccl->esum+rccl1->esum)-
      (rccl->esum* Hep3Vector(rccl->x ,rccl->y ,rccl->z ).unit()+
       rccl1->esum*Hep3Vector(rccl1->x,rccl1->y,rccl1->z).unit() ).mag2();
      if (invmass > mcut2) return 1;
    }
  }
  
  return 0;
}

int daq1_init() {
char* mode= getenv("WRITEMODE");
if (mode) Mode=atoi(mode); else Mode=0;
char* mcut= getenv("MCUT");
if (mcut) Mcutall=atof(mcut); else Mcutall=0;
char* ecut= getenv("ECUT");
if (ecut) Ecutall=atof(ecut); else Ecutall=0;
endianinit();
clear_rich();
return 0;
}

int daq1_stop() {
if (mf) fclose(mf);
 if (myeve){
   cout << "Average number of ECAL Clusters (cuts applied)" << (float) allitra/myeve;
   cout << "\nAverage number of OTR tracks " << (float) allotrtra/myeve;
   cout << "\nAverage number of PC tracks " << (float) allpc/myeve;
   cout << "\nAverage number of TC tracks " << (float) alltc/myeve;
 }
 return 0;
}

void initrich_();

int  daq1_write(){
  // ecut energy cut for cluster selection (GeV)
  // mcut invariant mass cut for event selection (GeV)

  static int ib=0;
  ArtePointer<EVHD> evhd(1);
  if (!myeve) {
    daq1_init();
    //    initrich_();
    char* daq1name= getenv("RECODAQ1");
    if (daq1name)
      sprintf(myname,"%s_file%.3d.daq1",daq1name,myifile); 
    else {
      char* path=getenv("HBSCRATCH");
      if (evhd->run>1) 
	sprintf(myname,"%s/run%.5d_file%.3d.daq1",path,evhd->run,myifile); 
      else
	sprintf(myname,"%s/run_mc_file%.3d.daq1",path,myifile); 
    }
    mf=fopen(myname,"w");
    if (mf==NULL) { 
      cout << "ERROR fopen" << endl;
      exit(1);
    }
    cout << " FILE fopen  :" << myname << endl;
    //    hrccl  =new   HBookHistogram("ECAL clusters",        100,0,100); 

  }
  myeve++;
  if (myeve==5000) { Ecut=Ecutall; Mcut=Mcutall;}
  if (evnt_select()||(Mcut==0)) write_raw_data( Mode );
  return OK;
}






void printevent(FILE *fp){ 
  int   NumChannelPerFED=1024;
  cout << "PRINT EVENT" << numfedevt << endl; 
  for (int i=0;i<numfedevt;i++) {
    SLB_SUBHDR* slb_header=&etb_fed[i].hdr;
    char* data=(char*)etb_fed[i].data;
    if (slb_header==0) {
      fprintf(fp,"\n No data in this event for fed %08X.\n",i);
      return;
      }
    fprintf(fp,"\n");
    fprintf(fp,"%d .fed\n",i);
    fprintf(fp,"    Length ................: %04X %6d\n",
          slb_header->length,slb_header->length);
    fprintf(fp,"    MotherID (BoardID) ....: %04X\n",slb_header->fedid);
    fprintf(fp,"    RecordID ..............: %04X\n",slb_header->recid);
    fprintf(fp,"    BX-ID .................: %04X\n",slb_header->bxid);
    fflush(fp);
    int i;
    for (i=0;i< NumChannelPerFED/8;i++) {
      fprintf(fp,"%02X",data[i]);
      if (i%32==31) fprintf(fp,"\n");
      else if (i%2==1) fprintf(fp," ");
      }
    fflush(fp);
    }
    return;
}