#include "rootspy_hist.hh"
#include <fcntl.h>
using namespace evio;
using namespace std;

static Hist *hist ;
static int event_cnt;
static FILE *RunFile;

extern "C" {
  //-----------------------------------------------------
  int   close_hist() {
    //printf("INFO::: cMsg Receive Stop  domainId=%d...\n",0);
    /*
    cMsgReceiveStop(domainId);
    cMsgUnSubscribe(domainId, subHandle1);
    cMsgUnSubscribe(domainId, subHandle2);
    cMsgDisconnect(&domainId);
    */
    printf("INFO::: Delete HIST class p=%p ...\n",hist);
    delete hist;  hist=NULL;
  }  
  //-----------------------------------------------------
  int   init_hist() {
    printf("INFO:::: Create HIST class:  ptr: before=%p ",hist);
    hist = new Hist;
    printf(" after=%p  \n",hist);
  }
  //-----------------------------------------------------
  int   reset_hist(int flag) {
    printf("INFO:::: Reset all histos  ...\n");
    hist->Reset(flag);
  }
  //-----------------------------------------------------
  int   fill_hist() {
    //printf("Hptr2 = %p  \n",hist);
    hist->test_fill();    
  }
  //-----------------------------------------------------
  int   fill_ped(int ped) {
    hist->ped_fill(ped);    
  }
  //-----------------------------------------------------
  int   fill_fadc(int slot, int chan, int adc ) {
    //printf("fill_fadc:: slot=%d   chan=%d,  adc=%d \n",slot, chan, adc);
    hist->Fill_FADC(slot,chan,adc);
  }
  //-----------------------------------------------------
  int   fill_bank(unsigned int  *bank, int size) {
    scan_bank(bank, size);
  }
}  // ---- end extern C  ---

//-----------------------------------------------------------------------------------------
//           WRITE BUFFER
//-----------------------------------------------------------------------------------------
void write_bank(unsigned int  *bank, int size) {

  //evioFileChannel *chan = (evioFileChannel*) bank;

  // create event tree and analyze
  //evioDOMTree eventTree(chan);

  // get list of all banks in event
  //evioDOMNodeListP bankList = eventTree->getNodeList();
  char host_name[80];
  gethostname(host_name, 80);    
  event_cnt++;
  char name[128];
  
  BankHeader Bank(bank); Bank.print();
  RocHeader ROC_Header(bank);    ROC_Header.print();

  sprintf(name,"/home/furletov/event_tof_%s_%d.evio",host_name,event_cnt);
  printf("open file %s\n",name);
  int rc=0;
  if((RunFile = fopen(name,"w")) == NULL) {
    printf("Can not open file %s\n",name); 
  } else {
    rc=fwrite(bank,4,size,RunFile);
  }
  if (RunFile)  fclose(RunFile);  RunFile=NULL;
  printf("close file  rc=%d\n",rc);
  //int fd=open(name,O_WRONLY);
  //int rc = write(fd, bank, size*4);
  //close(fd);
}

//-----------------------------------------------------------------------------------------
//       Banks decoding
//-----------------------------------------------------------------------------------------
int scan_bank(unsigned int  *buffer, int size) {
  
  //write_bank(buffer, size);

  BankHeader Bank(buffer); Bank.print();
  RocHeader ROC_Header(buffer);    ROC_Header.print();
  buffer+=2;


  unsigned i = 0;
  while (i < (size-2)) {
    printf("Loop:: i=%d size=%d\n",i,size-2);
    BankHeader Bank(buffer); Bank.print();
    if ((Bank.Tag&0xFF00)==0xFF00) { //---  Trigger INFO
      TriggerHeader Trigger_Header(buffer);    Trigger_Header.print();
      printf("       .... Tag=0x%X, Skip %d words \n",Bank.Tag,Bank.Size+1);
    }
    else if (Bank.Tag==6) {   //-- FADC250 event
      printf("   FADC250 event:: i=%d Tag=0x%X, size=%d (0x%X) words \n",i,Bank.Tag,Bank.Size,Bank.Size);
      analyzeBank_FA250(buffer);    
    }
    else if (Bank.Tag==20) {   //-- CAEN TDC  event
      printf("   CAEN TDC event:: i=%d Tag=0x%X, size=%d (0x%X) words \n",i,Bank.Tag,Bank.Size,Bank.Size);
      analyzeBank_CEAN_TDC(buffer);    
    }
    else if (Bank.Tag==26) {   //-- F1 TDC  event
      printf("   F1 TDC event:: i=%d Tag=0x%X, size=%d (0x%X) words \n",i,Bank.Tag,Bank.Size,Bank.Size);
      analyzeBank_F1_TDC(buffer);    
    }
    buffer+=(Bank.Size+1);
    i+=(Bank.Size+1);
  } //-- end while
}
//-----------------------------------------------------------------------------------------
int analyzeBank_F1_TDC(unsigned int  *buffer) {

#define F1_BLOCK_HEAD    16
#define F1_BLOCK_TRAIL   17
#define F1_EVENT_HEAD    18
#define F1_TRIG_TIME     19

#define F1_MEAS          23
#define F1_CHIP_HEAD     24
#define F1_NOT_VALID     30
#define F1_FILLER        31

#define F1_SLOT   0x7C00000
#define F1_ID     0x3C0000
#define F1_BLOCK  0x3FF00
#define F1_NEVT_BLOCK 0xFF

#define F1_TRIG_NUMB 0x3FFFFF

#define F1_TDE  0xFFFF
#define F1_CHIP 0x380000
#define F1_CHAN 0x70000
#define F1_TIME 0xFFFF


  int debug_tdc = 1;
  
  BankHeader Bank(buffer); Bank.print();
  
  int bank_size=Bank.Size-1;
  unsigned int  *bank = (buffer+2); 

  printf("F1_TDC::  bank ptr = %p size=%d\n",bank, bank_size);
  if (bank_size<=0) return 0;
  printf("F1_TDC::  H0=0x%08X   H1=0x%08X   \n",*bank,*(bank+1));
  fflush(stdout);


  debug_tdc = 1;
  
  int block_head  =   0;
  
  int slot        =  -1;
  int mod_id      =  -1;
  int block_numb  =  -1;
  int nevt_block  =  -1;
  
  int trig_numb   =  -1;
  int trig_time   =  -1;
  
  int chip_stat     =  -1;
  int hit_fifo_over =  -1 ;
  int out_fifo_over =  -1;
  
  int  chip_numb  =  -1;
  int  chip_chan  =  -1;
  int  channel    =  -1;
  int  time       =  -1;
  
  int  ch_tmp  = -1;
  
  
  for (unsigned i = 0; i < bank_size; i++){
    
     bswap(&bank[i]);
    
    if( (bank[i] >> 27) == F1_BLOCK_HEAD){
      if(block_head != 0){
	cout << "  FATAL:  Missing global trailer " << endl;
        
	//            exit(EXIT_FAILURE); 
      }
      
      slot        =  (bank[i] & F1_SLOT)   >>  22;
      mod_id      =  (bank[i] & F1_ID)     >>  18;
      block_numb  =  (bank[i] & F1_BLOCK)  >>  8;
      nevt_block  =  (bank[i] & F1_NEVT_BLOCK);
      
      if(debug_tdc)
	cout << " BLOCK HEADER " <<  " SLOT = " <<  slot  <<  "  ID = " << mod_id << 
	  "  BLOCK NUMBER = " << block_numb <<  "  NEVT_BLOCK = " << nevt_block << endl;
      
      block_head = 1;
      
      
    }
    
    if( (bank[i] >> 27) == F1_EVENT_HEAD){
      
      slot        =  (bank[i] & F1_SLOT)   >>  22;
      trig_numb = (bank[i] & F1_TRIG_NUMB);
      
      if(debug_tdc)
	cout << "  EVENT HEADER " <<  " SLOT = " <<  slot  <<  "  Trigger number = " << trig_numb << endl;
      
    }
    
    // Trigger Time
    if( (bank[i] >> 27) == F1_TRIG_TIME){
      
      trig_time      =  (bank[i] & F1_TDE);
      
      if(debug_tdc)
	cout << "  TRIGGER TIME DE  " <<  trig_time << endl;          
    }
    
    if( (bank[i] >> 27) == 0){
      trig_time      =  (bank[i] & F1_TDE);
      
      if(debug_tdc)
	cout << "  TRIGGER TIME AB  " <<  trig_time << endl;          
    }
    
    
    
    if( (bank[i] >> 27) == F1_MEAS){
      
      chip_stat     =  ((bank[i] >> 26) & 0x1);
      hit_fifo_over =  ((bank[i] >> 25) & 0x1);
      out_fifo_over =  ((bank[i] >> 24) & 0x1);
      
      chip_numb     =   (bank[i] & F1_CHIP)  >> 19;
      chip_chan     =   (bank[i] & F1_CHAN)  >> 16;         
      time          =   (bank[i] & F1_TIME);
      
      ch_tmp = -1;
      
      if( (chip_chan == 0) || (chip_chan == 1)) ch_tmp = 0;
      if( (chip_chan == 2) || (chip_chan == 3)) ch_tmp = 1;
      if( (chip_chan == 4) || (chip_chan == 5)) ch_tmp = 2;
      if( (chip_chan == 6) || (chip_chan == 7)) ch_tmp = 3;
      
      if( (ch_tmp >= 0) && (ch_tmp < 4) ){
	channel  = chip_numb*4 + ch_tmp;
      }
      
      cout << "  HITS:  " <<   " Chip  status = "  << chip_stat << 
	"  Hit FIFO over = " <<  hit_fifo_over <<  "  Out FIFO over " << out_fifo_over << endl;
      
      cout <<  "  Chip = " << chip_numb <<  "   Chip Chan  = " << chip_chan <<  
	"   Channel = "  <<  channel <<  "   Time = " << time << endl;
      
      if ( 0<= slot) { //-- fill hist ---
	hist->Fill_F1TDC(slot,channel,time);
      }
      
    }
    
    
    if( (bank[i] >> 27) == F1_BLOCK_TRAIL){
      if(block_head != 1){
	cout << "  FATAL:  Missing global header " << endl;
	//      exit(EXIT_FAILURE); 
      }
      
      block_head = 0;
      
    }
    
    
  } // Loop over TDC data words
  
  


}
//-----------------------------------------------------------------------------------------
int analyzeBank_CEAN_TDC(unsigned int  *buffer) {

#define CAEN_GLOB_HEAD   8
#define  EVT_COUNT  0x7FFFFD0
#define  GEO  0x1F
#define CAEN_GLOB_TRAIL 16

  BankHeader Bank(buffer); Bank.print();
  
  int bank_size=Bank.Size-1;
  unsigned int  *bank = (buffer+2); 

  printf("CAEN_TDC::  bank ptr = %p size=%d\n",bank, bank_size);
  if (bank_size<=0) return 0;
  printf("CAEN_TDC::  H0=0x%08X   H1=0x%08X   \n",*bank,*(bank+1));
  fflush(stdout);

  int debug = 1;
  
  //const vector<uint32_t> *vec = bankPtr->getVector<uint32_t>();
  //if(vec == NULL) {cerr << "?unable to get tdc bank vector" << endl; return;}
  //  if(debug) cout << " CAEN1290 TDC BANK FOUND " <<  vec->size() << " Event number  " << (*vec)[0] << endl;

  int global_head  = 0;     
  int event_count  = -1;
  int caen_slot    = -1;
  
  for (unsigned i = 0; i < bank_size; i++){

    bswap(&bank[i]);

    if( (bank[i] >> 27) == CAEN_GLOB_HEAD){
      if(global_head != 0){
	cout << "CEAN_TDC::  FATAL:  Missing global trailer " << endl;
	return 1;
      }
      event_count = (bank[i] & EVT_COUNT) >> 5;
      cout << "CEAN_TDC:: EVENT COUNT = " << event_count << endl;
      caen_slot = (bank[i] & GEO);
      cout << "CEAN_TDC:: SLOT = " << caen_slot << endl;        
      global_head = 1;      
    }
    
    if( (bank[i] >> 27) == CAEN_GLOB_TRAIL){
      if(global_head != 1){
	cout << "  FATAL:  Missing global header " << endl;
	return 1;
      }      
      global_head = 0;      
    }
    
    if( (bank[i] >> 27) == 1){
      unsigned int tdc = (bank[i] & 0x3000000) >> 23;
      //  if(debug)  cout << " TDC header found  " <<  " TDC = " << tdc <<  "  i = " << i << endl;
    }
    // TDC measurement
    if( (bank[i] >> 27) == 0){
      unsigned int tdc_channel = (bank[i] & 0x3E00000) >> 21;
      unsigned int tdc_time  =  (bank[i] & 0x1FFFFF);
      if(debug) cout << " TDC channel   " <<   tdc_channel << " TDC time =  " << tdc_time  <<  endl;
      if (caen_slot>=0 ) { //-- fill hist ---
	hist->Fill_CAEN(caen_slot,tdc_channel,tdc_time);
      }
    }
  }  //--  end for ---
}
//-----------------------------------------------------------------------------------------
int analyzeBank_FA250(unsigned int  *buffer) {

  BankHeader Bank(buffer); Bank.print();
  
  //int bank_size=size;
  int bank_size=Bank.Size-1;
  unsigned int  *bank = (buffer+2); 

  printf("FA250::  bank ptr = %p size=%d\n",bank, bank_size);
  if (bank_size<=0) return 0;
  printf("FA250::  H0=0x%08X   H1=0x%08X   \n",*bank,*(bank+1));
  fflush(stdout);

  int debug_fadc = 0;
  int SlotID  =  -1;
  const int max_hits = 3000;
  
  int    nhit;
  int    hcrate[max_hits];
  int    hslot[max_hits];
  int    hch[max_hits];
  int    nsamp[max_hits];
  int    samp[max_hits][100];
  
  // Loop over data      
  unsigned int data_new_type  =   0;
  unsigned int data_type      =   0;
  int channel = -1 ;
  unsigned int fadc_raw_wind_width = 0;
  
  SlotID = -1;
  unsigned int time_stamp     =   0;
  
  unsigned int pulse_number   =  -1;
  unsigned int first_sample   =  -1;
  
  
  //      cout << " Vector Size = " << vec->size() << endl;
  

  unsigned char* word;
  
  for (unsigned i = 0; i < bank_size; i++){
    
    
    //if (i==0) printf("Bank:1: H0=0x%08X ",bank[0]);
    
    word=(unsigned char*)&bank[i];
    unsigned char byte;   // swap 
    byte=word[0];  word[0]=word[3];  word[3]=byte;
    byte=word[1];  word[1]=word[2];  word[2]=byte;
    
    //if (i==0) printf("Bank:2: H0=0x%08X ",bank[0]);
    
    // Bank header
    // bit  31       =  1
    // bits 27- 30   =  0
    // bits 22 - 26  =  SlotID 
    // bits 11 - 21  =  Number of events in block;  bits 0 - 10 event block number
    if( bank[i] >> 27 == 0x10)         /* Bank header 1000 0 */
      {
	SlotID =  (bank[i] & 0x7c00000) >> 22;
	if(debug_fadc)
	  cout << " FADC found in slot:  " << SlotID << endl;
      }
    
    // Event header
    // bit  31       =  1
    // bits 27- 30   =  0x1001
    // bits 0 - 26   =  Trigger number 
    if( bank[i] >> 27 == 0x12)         /* Event header 1001 0 */
      {
	//            trigger_numb[nboard] =  (bank[i] & 0x7FFFFFF);
	
	if(debug_fadc)
	  cout << " FADC found in slot:  " << SlotID << " Trigger time " << (bank[i] & 0x7FFFFFF) << endl;
      }
    
    
    
    if( bank[i] & 0x80000000 )         /* data type defining word */
      {
	data_type = (bank[i] & 0x78000000) >> 27;
	// data_type == 4  - window raw data
	// data_type == 6  - pulse raw data
	data_new_type = 1;
      }
    else
      {
	data_new_type = 0;
      };
    
    
    //	cout << " Data type = " << data_type << endl;
    
    
#if 1
    // WINDOW RAW DATA	
    if(data_type == 4 ){
      
      int raw_debug = 0;
      
      if(data_new_type == 1){
	fadc_raw_wind_width = (bank[i] & 0xFFF);
	channel    = (bank[i] & 0x7800000) >> 23;
	//  cout << " Data Type = " << data_type << " Fadc_width  = " << fadc_width << " Channel = " << channel << endl;
	
	
	time_stamp = 1;
	nhit ++;
	//	    cout << " NEW CHANNEL = " << channel << " CRATE = " << crate << " SLOT = " << SlotID <<  " NHIT = " << nhit << endl;
	
      } else{
	
	
	unsigned int valid_1 = 1;
	unsigned int valid_2 = 1;
	unsigned int adc_1 = (bank[i] & 0x1FFF0000) >> 16;
	
	if(bank[i] & 0x20000000 )
	  valid_1 = 0;
	
	unsigned int adc_2 = (bank[i] & 0x1FFF);
	if( bank[i] & 0x2000 )
	  valid_2 = 0;
	
	
	
	if(raw_debug)
	  cout << " Time stamp = " << time_stamp <<  " SlotID = " << SlotID << " Channel = " << channel << 
	    "  ADC1 = " <<  float(adc_1) <<  
	    " " << float(adc_2) <<  endl;
	
	
	if((SlotID>=0) && (channel >= 0) && (channel < 16)){
	  
	  hist->Fill_FADC(SlotID,channel,adc_1);
	  hist->Fill_FADC(SlotID,channel,adc_2);
	  
	} else {cout << " Wrong slot/channel number  slot=" <<  SlotID << "  channel =" << channel << endl;}
	
	time_stamp++;
	
      }
    }   // Data type 4, raw window data
#endif	
    
    int pulse_debug = 0;
    
    
    
    // Pulse Integral data
    if(data_type == 7 ){
      
      channel    = (bank[i] & 0x7800000) >> 23;
      unsigned int pulse_numb = (bank[i] & 0x200000) >> 21;
      unsigned int pulse_int  = (bank[i] & 0x3FFFF);
      
    }
    
  }
}
//-----------------------------------------------------------------------------------------
//           FILL ADC
//-----------------------------------------------------------------------------------------
void Hist::Fill_FADC(int slot, int chan, int adc) {
  if (!FADC[slot][chan]) {
    sprintf(hist_name,"%s_Slot_%02d_FADC_%02d",host_name,slot,chan);
    printf("Book hist: %s \n",hist_name);
    FADC[slot][chan] = new TH1D(hist_name, "ADC value", 500, -0.5, 499.5); 
    ALL_1D_HIST->Add(FADC[slot][chan]);
    //main->cd();
  }
  //printf("Fill::: Slot=%02d FADC=%02d  ptr=%p \n",slot,chan,FADC[slot][chan]);
  FADC[slot][chan]->Fill((double)adc);
  ped_prof->Fill(float(slot*MAXCHAN+chan),float(adc));
}

//-----------------------------------------------------------------------------------------
//           FILL F1TDC TDC
//-----------------------------------------------------------------------------------------
void Hist::Fill_F1TDC(int slot, int chan,unsigned int tdc) {
  if (!F1TDC[slot][chan]) {
    sprintf(hist_name,"%s_Slot_%02d_F1_TDC_%02d",host_name,slot,chan);
    printf("Book hist: %s \n",hist_name);
    F1TDC[slot][chan] = new TH1D(hist_name, "F1 TDC value", 500, 0., 70000.0); 
    ALL_1D_HIST->Add(F1TDC[slot][chan]);
    //main->cd();
  }
  //printf("Fill::: Slot=%02d FADC=%02d  ptr=%p \n",slot,chan,FADC[slot][chan]);
  F1TDC[slot][chan]->Fill((double)tdc);
  tdc_prof_f1->Fill(float(slot*MAXCHAN_F1+chan),float(tdc));
}

//-----------------------------------------------------------------------------------------
//           FILL CAEN TDC
//-----------------------------------------------------------------------------------------
void Hist::Fill_CAEN(int slot, int chan,unsigned int tdc) {
  if (!CTDC[slot][chan]) {
    sprintf(hist_name,"%s_Slot_%02d_C_TDC_%02d",host_name,slot,chan);
    printf("Book hist: %s \n",hist_name);
    CTDC[slot][chan] = new TH1D(hist_name, "CAEN TDC value", 500, 0., 10000.0); 
    ALL_1D_HIST->Add(CTDC[slot][chan]);
    //main->cd();
  }
  //printf("Fill::: Slot=%02d FADC=%02d  ptr=%p \n",slot,chan,FADC[slot][chan]);
  CTDC[slot][chan]->Fill((double)tdc);
  tdc_prof_caen->Fill(float(slot*MAXCHAN_CAEN+chan),float(tdc));
}

//-----------------------------------------------------------------------------------------
//          RESET ALL HISTS 
//-----------------------------------------------------------------------------------------
int Hist::Reset(int flag) {
  printf(" INFO::: HIST Reset()::  \n");
  ALL_1D_HIST->Print();
  TObject *obj;
  TH1D *h1d;
  TIter next(ALL_1D_HIST);
  while ((obj = next())) { 
    obj->Print();
    h1d=(TH1D *) obj;
    h1d->Reset(" ");
  }  
}
//-----------------------------------------------------------------------------------------
//      BOOK HIST
//-----------------------------------------------------------------------------------------
Hist::Hist () {
    cout << " Constructor Hist !!! " << "\n";
    gethostname(host_name, 80);    
    rootspy = new DRootSpy();
    // Random number generator
    rndm = new TRandom();
    ALL_1D_HIST=new TList();
    FADC250_HIST=new TList();
    FADC125_HIST=new TList();
    F1TDC_HIST=new TList();
    CAENTDC_HIST=new TList();

    MAXCHAN=16;
    MAXSLOT=25;
    MAXBIN=MAXSLOT*MAXCHAN;
    
    MAXCHAN_CAEN=32;
    MAXSLOT_CAEN=25;
    MAXBIN_CAEN=MAXSLOT_CAEN*MAXCHAN_CAEN;

    
    MAXCHAN_F1=32;
    MAXSLOT_F1=25;
    MAXBIN_F1=MAXSLOT_F1*MAXCHAN_F1;

    //    signal(SIGINT, sigHandler);
    DONE = false;
    
    // Define some histograms to file
    
    main = gDirectory;
    gDirectory->mkdir("components")->cd();
    
    for (int is=0; is<MAXSLOT; is++) {
      for (int ic=0; ic<MAXCHAN; ic++) {
	FADC[is][ic] = NULL;
      }
    }
    
    for (int is=0; is<MAXSLOT_CAEN; is++) {
      for (int ic=0; ic<MAXCHAN_CAEN; ic++) {
	CTDC[is][ic] = NULL;
      }
    }
    
    for (int is=0; is<MAXSLOT_F1; is++) {
      for (int ic=0; ic<MAXCHAN_F1; ic++) {
	F1TDC[is][ic] = NULL;
      }
    }
    
     sprintf(hist_name,"%s_ped_prof",host_name);
    printf("Book hist: %s \n",hist_name);
    ped_prof = new TProfile(hist_name,"Pedestals profile",MAXBIN,-0.5,MAXBIN-0.5,-100.,1000.); 
    ALL_1D_HIST->Add(ped_prof);

    sprintf(hist_name,"%s_tdc_prof_caen",host_name);
    printf("Book hist: %s \n",hist_name);
    tdc_prof_caen = new TProfile(hist_name,"TDC profile Caen ",MAXBIN_CAEN,-0.5,MAXBIN_CAEN-0.5,-10.,10000.); 
    ALL_1D_HIST->Add(tdc_prof_caen);
    
    sprintf(hist_name,"%s_tdc_prof_f1",host_name);
    printf("Book hist: %s \n",hist_name);
    tdc_prof_f1 = new TProfile(hist_name,"F1TDC profile  ",MAXBIN_F1,-0.5,MAXBIN_F1-0.5,-10.,70000.); 
    ALL_1D_HIST->Add(tdc_prof_f1);
    
    /*
      char host_name[80];
      gethostname(host_name, 80);    
      char hist_name[80];
      for (int is=0; is<MAXSLOT; is++) {
      for (int ic=0; ic<MAXCHAN; ic++) {
      sprintf(hist_name,"%s_Slot_%02d_FADC_%02d",host_name,is,ic);
      FADC[is][ic] = new TH1D(hist_name, "ADC value", 500, 0.0, 2000.0); 
      }
      }
      main->cd();
    */
}
  
//-----------------------------------------------------------------------------------------
//      TEST FILL
//-----------------------------------------------------------------------------------------
int   Hist::ped_fill(int ped) {	
  double rp = rndm->Gaus(ped,10);

  for (int is=0; is<MAXSLOT; is++) {
    for (int ic=0; ic<MAXCHAN; ic++) {
      hist->Fill_FADC(is,ic,rp+is-ic);
    }
  }
}  
//-----------------------------------------------------------------------------------------
int   Hist::test_fill() {	
  
  // Loop forever while filling the hists
  //cout<<endl<<"Filling histograms ..."<<endl;
  unsigned int Nevents=0;
  
  
  // Randomly choose values to put into an x and y coordinate on some plane 6 meters away from
  double x = (rndm->Rndm()-0.5)*2.0;
  double y = (rndm->Rndm()-0.5)*2.0;
  double z = 6.0;
  double p = rndm->Landau(0.5, 0.075); // total momentum
  double m = 0.134;
  double Etot = sqrt(p*p + m*m) + rndm->Gaus(0.0, 0.005);
  
  px = p*x/z;
  py = p*y/z;
  pz = sqrt(p*p*(1.0 - (x*x+y*y)/(z*z)));
  E = Etot;
  Mass = sqrt(Etot*Etot - p*p);

  if (!h_px) {  //-- first call, Book hist --
    h_px = new TH1D("px", "Momentum X-component", 500, 0.0, 10.0);
    h_py = new TH1D("py", "Momentum Y-component", 500, 0.0, 10.0);
    h_pz = new TH1D("pz", "Momentum Z-component", 500, 0.0, 10.0);

    h_E = new TH1D("E", "Energy", 500, 0.0, 10.0);
    h_Mass = new TH1D("Mass", "Mass", 1000, 0.0, 2.0);

    for(int i=0; i< 100; i++) data[i] = (float)i;
    // Define some trees to save
    
    T = new TTree("T", "Event Info");
    T->Branch("px", &px, "px/D");
    T->Branch("py", &py, "py/D");
    T->Branch("pz", &pz, "pz/D");
    T->Branch("E", &E, "E/D");
    T->Branch("Mass", &Mass, "Mass/D");
    T->Branch("data", &data, "data[100]/F");
    
    // Set nice labels for X-axes
    h_px->SetXTitle("p_{x} (GeV/c)");
    h_py->SetXTitle("p_{y} (GeV/c)");
    h_pz->SetXTitle("p_{z} (GeV/c)");
    h_E->SetXTitle("Energy (GeV)");
    h_Mass->SetXTitle("mass (GeV/c^2)");
  }
  
  h_px->Fill(px);
  h_py->Fill(py);
  h_pz->Fill(pz);
  h_E->Fill(E);
  h_Mass->Fill(Mass);
  
  
  // Limit how large the TTree can get
  if(Nevents<MAX_TREE_EVENTS){
    T->Fill();
  }else if(Nevents == MAX_TREE_EVENTS){
    cout << MAX_TREE_EVENTS/1000 << "k events in tree. Only histograms will be filled from here on." << endl;
  }
  
  
  if(((++Nevents) % 100)==0){
    //gDirectory->ls();
    
    cout <<" "<<Nevents<<" events simulated "<<endl;
    
    //cout<<" "<<Nevents<<" events simulated "<<"\r";
    //cout.flush();
  }
  
  // sleep a bit in order to limit how fast the histos are filled
  //usleep(5000);
  
  // cout << endl;
  // cout << "Ending" << endl;
  
  return 0;
}
//-----------------------------------------------------------------------------------------
unsigned int aswap(unsigned int *bank, int len) {
  for (int i=0; i<len; i++) {
    unsigned char* bytes=(unsigned char*)&bank[i];
    unsigned char b;   // swap 
    //printf("i=%d  Before=0x%08X  ...",i, bank[i]);
    b=bytes[0];  bytes[0]=bytes[3];  bytes[3]=b;
    b=bytes[1];  bytes[1]=bytes[2];  bytes[2]=b;
    //printf("After=0x%08X \n",bank[i]);
  }
}
unsigned int bswap(unsigned int *bank) {
    unsigned char* bytes=(unsigned char*)bank;
    unsigned char b;   // swap 
    b=bytes[0];  bytes[0]=bytes[3];  bytes[3]=b;
    b=bytes[1];  bytes[1]=bytes[2];  bytes[2]=b;
}