/*
read adc125 data with random trigger determine the pedestal and average
every 1000 events. Then change the DAC settings for each channel. Do
this iteratively till reaching convergens to the setvalue baseline.
The code starts out with the channels set to DAC value 0x9000. In the
second step it set the DAC values to 0x8a00. This serves as basis for
calibrating the DAC of each channel to the base line.
the base line is to be either 100 or 20 times the pedestal width, whichever
is gerater. 
*/

#include <iostream>
#include <fstream>
#include <iomanip>

using namespace std;

extern "C"{
#include <stdlib.h>
#include <string.h>
#include <time.h>
#include <jvme.h>
#include <tiLib.h>
#include <sdLib.h>
#include <fa125Lib.h>
#include "math.h"
#include "time.h"   
#include "sys/types.h"
#include <unistd.h>

  extern int nfa125;
  extern int fa125A32Base;
  extern int tiA32Base;
  
  DMA_MEM_ID vmeIN,vmeOUT;
  extern unsigned int *dma_dabufp;
  extern DMANODE *the_event;
  
}

//UINT32 F125DATA[40000];

#define FA125_ADDR (5<<19)
#define TI_ADDR (21<<19)

ofstream OUTF;
ofstream OUTFPed;
ofstream DACval;

void myReadoutCode();
void setDACvalues();
void mkpedestals(int);

char DETECTOR[128];

unsigned int START;
int HOSTROC = 0;
int NADCSAMPLES = 100;
int LATENCY = 300;
int THRESHOLD = 300;
int NADCS = 1;
int TIMETEST = 0;
int NEVENT = 0;
int RAWDATA = 0;
int DEBUG =0;
int DEBUGGER = 0;

int NInRange = 0;
int NConverged = 0;
int MINMEAN = 1000;
int MAXMEAN = 4000;
int DETECTOR_FACTOR = 1;
unsigned int timelaps = 0;

unsigned int Pedestal[20][72];
unsigned int PedestalOLD[2][20][72];
int ADCSamples[20][72][200];

float Ped[20][72];
int EventCounter = 0;

float DACGain[20][72];
int BASELINE = 20;
int EVCOUNTER = 20;

unsigned short int CurrentDACValue[20][72];
unsigned short int PreviousDACValue[20][72];
float PreviousMean[20][72];
float CurrentMean[20][72];

bool BaselineInRange[20][72];
bool BaselineConverged[20][72];
float PedMean[20][72];
float Q[20][72];
unsigned int NSamples[20][72];
float PedestalSigma[20][72];
float PedestalValues[20][72];

int BaseLine[20][72];

int ReferenceBaseLine = 100;

int PedestalCounter = 0;
int ERRORCounter = 0;
int ERRORLOG = 10000;

int  main(int argc, char *argv[])  {

  int Nevents = 1000;
  int iFlag;
  int stat;
  char ofnam[128];
  char ofnamped[128];
  char ofnamdac[128];
  char CrateName[128] = "One";
  char ifnam[128];

  char DataDir[128];

  printf("\nJLAB f125ADC pedestal Tests\n");
  printf("----------------------------\n");
  
  int k=0;
  int SLOT=3;
  if (argc<2) {
    cout<<argc<<endl;
    cout<<"usage: pulserTestNew [Options]"<<endl;
    cout<<"       -S  #        Slot number of first ADC to be used in test (default 3)"<<endl;
    cout<<"       -N  #        Number of triggers default is 1000"<<endl;
    cout<<"       -X  #        Number of samples to read (default 100)"<<endl;
    cout<<"       -A  #        All # ADCs from slot 3 to 3+# (default 1)"<<endl;
    cout<<"       -C  string   Crate Name (default \"One\")"<<endl;
    cout<<"       -B  #        set Base line to # (Default is 100)"<<endl;
    cout<<"       -L  #        set latency to # (Defaulat is 500)"<<endl;
    cout<<"       -TT          Timing Test" <<endl;
    cout<<"       -R           write out RAWDATA" <<endl;
    cout<<"       -D           write out DEBUG info" <<endl;
    cout<<"       -DD          write out DEBUGGGER info" <<endl;
    cout<<"       -h           print this help" <<endl;
    return 0;
  }

  for (k=0;k<argc;k++){
    
    if (!strcmp(argv[k],"-S")){
      SLOT = atoi(argv[k+1]);
      k++;
    }
    
    if (!strcmp(argv[k],"-B")){
      ReferenceBaseLine = atoi(argv[k+1]);
      k++;
    }
    if (!strcmp(argv[k],"-L")){
      LATENCY = atoi(argv[k+1]);
      k++;
    }
    
    if (!strcmp(argv[k],"-N")){
      Nevents = atoi(argv[k+1]);
      k++;
    }
    
    if (!strcmp(argv[k],"-X")){
      NADCSAMPLES = atoi(argv[k+1]);
      k++;
      if (NADCSAMPLES%2){
	NADCSAMPLES +=1;
      }
    }
    if (!strcmp(argv[k],"-A")){
      NADCS = atoi(argv[k+1]);
      k++;
    }
    if (!strcmp(argv[k],"-C")){
      sprintf(CrateName,"%s",argv[k+1]);
      k++;
    }
    if (!strcmp(argv[k],"-TT")){
      TIMETEST = 1;
    }
    if (!strcmp(argv[k],"-R")){
      RAWDATA = 1;
    }
    if (!strcmp(argv[k],"-D")){
      DEBUG = 1;
    }
    if (!strcmp(argv[k],"-DD")){
      DEBUGGER = 1;
    }
    if (!strcmp(argv[k],"-T")){
      ReferenceBaseLine = atoi(argv[k+1]);
      k++;
    }
    
    if (!strcmp(argv[k],"-h")){
      cout<<"usage: pulserTestNew [Options]"<<endl;
      cout<<"       -S  #        Slot number of first ADC to be used in test (default 3)"<<endl;
      cout<<"       -N  #        Number of triggers default is 1000"<<endl;
      cout<<"       -X  #        Number of samples to read (default 100)"<<endl;
      cout<<"       -A  #        All # ADCs from slot 3 to 3+# (default 1)"<<endl;
      cout<<"       -C  string   Crate Name (default \"One\")"<<endl;
      cout<<"       -B  #        set Base line to # (Default is 100)"<<endl;
      cout<<"       -L  #        set latency to # (Defaulat is 500)"<<endl;
      cout<<"       -TT          Timing Test" <<endl;
      cout<<"       -R           write out RAWDATA" <<endl;
      cout<<"       -D           write out DEBUG info" <<endl;
      cout<<"       -DD          write out DEBUGGGER info" <<endl;
      cout<<"       -h           print this help" <<endl;
      return 0;
    }
    
  }

  if(vmeOpenDefaultWindows() != OK) {
    cout<<"vmeOpenDefaultWindows failed. Terminate right here."<<endl;
    return 0;
  } 


  for (int n=0;n<20;n++){
    for (int k=0;k<72;k++){
	PedMean[n][k] = 0;
	Q[n][k] = 0;
	NSamples[n][k] = 0;
      PreviousDACValue[n][k] = -1;
      CurrentDACValue[n][k] = 0x9000;
      PreviousMean[n][k] = -1.0;
      CurrentMean[n][k] = -1.0;
      BaselineInRange[n][k] = false;
      BaselineConverged[n][k] = false;
    }
  }

  // determine username, host name, roc number and detector
  
  char username[128];
  char *val;
  //getlogin_r(username, 128);
  val = getenv("USER");
  strcpy(username,val);
  char hnam[128];
  gethostname(hnam,128);
  char s[128];
  strncpy(s,&hnam[6],128);
  HOSTROC = atoi(s);
  strncpy(DETECTOR,&hnam[3],3);
  DETECTOR[3] = '\0';
  cout<<"This code is run by user: "<<username<<endl;
  cout<<" HOST "<<s<<"   Detector: "<<DETECTOR<<endl;

  if(DETECTOR[0] == 'f')  DETECTOR_FACTOR = 2;

  sprintf(DataDir, "/gluonraid1/Users/%s/%sdata",username,DETECTOR,DETECTOR );  
  if (NADCS==1){
    sprintf(ofnam,"%s/adcdata_%s_slot%d.dat",DataDir,CrateName,SLOT); 
    sprintf(ofnamped,"%s/adcdata_%s_slot%d_pedestals.dat",DataDir,CrateName,SLOT); 
  } else {
    sprintf(ofnam,"%s/adcdata_%s_ALLslots.dat",DataDir,CrateName); 
    sprintf(ofnamped,"%s/adcdata_%s_ALLslots_pedestals.dat",DataDir,CrateName); 
    sprintf(ofnamdac,"%s/DAC/dacvalues_%s_ALLslots.dat",DataDir,CrateName); 
  }


  vmeDmaConfig(2,5,1);
  /* INIT dmaPList */
  dmaPFreeAll();
  vmeIN  = dmaPCreate("vmeIN",1024*4000,5,0);
  vmeOUT = dmaPCreate("vmeOUT",0,0,0);
  dmaPStatsAll();
  dmaPReInitAll();
    
  tiInit(TI_ADDR, TI_READOUT_EXT_POLL, 2);
  tiCheckAddresses();
  tiDisableDataReadout();  // do not read out TI
  tiDisableA32();          // do not read out TI   
  
  tiLoadTriggerTable(0);
  tiSetTriggerHoldoff(1,10,1);
  //tiSetTriggerHoldoff(2,40,1);

  tiSetPrescale(0);
  tiSetBlockLevel(1);

  tiSetTriggerSource(TI_TRIGGER_PULSER);
  //tiEnableTSInput(0x1);
  tiSetBlockLimit(0);

  tiSetBusySource(TI_BUSY_LOOPBACK, 1);
  tiSetBlockBufferLevel(1);
  tiSetFiberDelay(1,2);
  tiSetSyncDelayWidth(1, 0x3f, 1);
  
  // setup clock, etc for single-crate system only
  tiClockReset();
  taskDelay(1);
  tiTrigLinkReset();
  taskDelay(1);
  tiEnableVXSSignals();
  taskDelay(1);
  tiSyncReset(0);
  taskDelay(1);
  tiStatus(0);

  //tiSetOutputPort(1, 0, 0, 0);
  // tiClockResync();
  // tiClockReset();
  // tiSyncReset();
   
  extern unsigned int fa125AddrList[FA125_MAX_BOARDS];
  
  if (NADCS==1){
    fa125AddrList[0] = (SLOT<<19);  // adc in slot SLOT
  } else {
    int offset = 0;
    for (int k=0;k<NADCS;k++){
      fa125AddrList[k] = ((k+SLOT+offset)<<19);
      if (k+SLOT==10){
	offset = 2;
      }
    }
  }
  
  iFlag  = (1<<17); /* use fa125AddrList */
  /*
    bits 2-1:  defines Trigger source
    (0) 0 0  VXS (P0)
    (1) 0 1  Internal Timer
    (2) 1 0  Internal Multiplicity Sum
    (3) 1 1  P2 Connector (Backplane)
    bit    3:  NOT USED WITH THIS FIRMWARE VERSION
    bits 5-4:  defines Clock Source
    (0) 0 0  P2 Connector (Backplane)
    (1) 0 1  VXS (P0)
    (2) 1 0  Internal 125MHz Clock
  */
  iFlag |= (0<<1);  /* Trigger Source */
  iFlag |= (1<<4);  /* Clock Source */
  iFlag |= (1<<18); /* Skip firmware check*/

  stat = fa125Init(0, 0, NADCS, iFlag);
  if (stat != OK) {
    cout<<"ERROR: fa125Init failed, stop right here!!!"<<endl;
    return 0;
  }
  
  // the following needs to be done after initializing
  // only then fa125ScanMask() will work properly
  
  sdInit(0);
  sdSetActiveVmeSlots( fa125ScanMask());
  sdStatus();
  
  for (int k=0;k<NADCS;k++){
    int S = fa125Slot(k);
    fa125PowerOn(S);

    //fa125Reset(S, 0);
    //taskDelay(1);
 
   cout<<"Set DAC values for all channels in ADC slot "<<S<<endl;
    for(int ichan=0; ichan<72; ichan++) {
      fa125SetOffset(S,ichan,0x0000);
      CurrentDACValue[k][ichan] = 0x0000;
    }
    //set threshold to zero
    fa125SetCommonThreshold(S, 1);
    
    cout<<"Configure readout window size and latency to "<<NADCSAMPLES<<" and "<< LATENCY<<" resp."<<endl;
    //int pmode = 8;
    unsigned int PL = LATENCY;   // look back PL samples
    unsigned int PTW = NADCSAMPLES;  // window size in samples
    //unsigned int NSB = 3;    // sampels before threshold
    //unsigned int NSA = 30;   // sampels after threshold
    //unsigned int NP = 1;     // number of pulses in sample
    //fa125SetProcMode(S, pmode, PL, PTW, NSB, NSA, NP);

    if (!strcmp(DETECTOR,"cdc")){
      cout<<"Configure for CDC_long"<<endl; 
      fa125SetProcMode(S,"CDC_long", 500, NADCSAMPLES, NADCSAMPLES, 4, 1, 4, 4);
    }else {
      cout<<"Configure for FDC_long"<<endl; 
      fa125SetProcMode(S,"FDC_sum_long", 500, NADCSAMPLES, NADCSAMPLES, 4, 1, 4, 4);
    }

    // the following settings are no really need for this pedestal adjustment program
    fa125SetCommonTimingThreshold(S, 25, 100);
    fa125SetScaleFactors(S, 4, 3, 0);


    // enable all channels in module
    fa125SetChannelEnableMask(S,0xFFFFFF,0xFFFFFF,0xFFFFFF);
    // set blocklevel to 1 
    fa125SetBlocklevel(S, 1);

    fa125PrintTemps(S);
    //fa125Clear(S);  
    fa125Reset(S,0);
    taskDelay(1);
    fa125Enable(S);
    //taskDelay(1);
    //fa125Status(S,0);
  }
  

  tiSyncReset(0);

  //  tiIntConnect(TI_INT_VEC, myReadoutCode, tiIntArg);
  //  TI_INT_VEC =  0xec defined in tiLib.h
  //  tiIntArg = 0 defined in tiLib.c
  tiIntConnect(TI_INT_VEC, myReadoutCode, 0);
  vmeCheckMutexHealth(10);
  fa125GStatus(0);

  printf("Hit return to start Triggers...\n");
  getchar();
  
  /* Enable the TI and clear the trigger counter */
  tiIntEnable(0);
  taskDelay(1);
  tiStatus(0);
  //fa125Status(0,0);
  fa125GStatus(0);
  fa125ResetToken(0);
  //taskDelay(1);
 
  if (!TIMETEST){
    OUTF.open(ofnam,ios::out);
    OUTFPed.open(ofnamped,ios::out);
  }
  
  cout<<clock()<<"   START OF LOOP"<<endl; 
  START = clock();
  
  tiSetRandomTrigger(1,0x8);
  taskDelay(10);

  while (tiGetIntCount()<Nevents){
    
    //cout << "\r" <<"Event : "<<tiGetIntCount();
    if (NConverged == NADCS*72) break;
  }

  taskDelay(1);

    //cout << endl;	
  // end of loop overall events
  tiStatus(0);
  //fa125Status(0,0);

  tiIntDisable();
  tiIntDisconnect(); 
  
  for (int k=0;k<NADCS;k++){
    int S = fa125Slot(k);
    fa125PrintTemps(S);
    fa125PowerOff(S);
  }
  
  cout<<clock()<<"    END OF LOOP"<<endl;    
  unsigned int ENDL = clock();
  
  //tiIntDisable();
  cout<<"TOTAL TIME: "<<(ENDL-START)/CLOCKS_PER_SEC<<" seconds"<<endl;
    
  vmeCloseDefaultWindows();      

  if (!TIMETEST){
    OUTF.close();
    //mkpedestals(1);
    OUTFPed.close();
  }
  DACval.open(ofnamdac,ios::out);

  for (int k=0;k<NADCS;k++){
    int VMESlot = k+3;
    if (VMESlot>10){
      VMESlot +=2;
    }
    for (int j=0;j<72;j++){
	float locMean = PreviousMean[k][j];
	float locPedestalSigma = PedestalSigma[k][j];
	int locDACValue = CurrentDACValue[k][j];
	float locDACGain = DACGain[k][j];
	if(locDACGain == 0){
		printf("ERROR: DAC gain 0 (DAC=%i, mean=%f, sigma=%f, crate=%i, slot=%i, channel=%i) Setting sigma = 250\n", 
		       locDACValue, locMean, locPedestalSigma, HOSTROC, VMESlot, j);
                //locDACValue = 0;
                locPedestalSigma = 250.0;
                if (locMean = 0.0) locMean = 4095;
        }
	else if (locDACValue < 0){
		printf("ERROR: DAC Value Negative (DAC=%i, mean=%f, sigma=%f, crate=%i, slot=%i, channel=%i) Setting sigma = 250\n", 
		       locDACValue, locMean, locPedestalSigma, HOSTROC, VMESlot, j);
		//locDACValue = 0; 
		locPedestalSigma = 250.0;
		//if (locMean < 0.0) locMean = 0.0;
	}
	else if ((locPedestalSigma != locPedestalSigma) || locPedestalSigma == 0.0){
		printf("ERROR: Sigma nan or Equals 0 (DAC=%i, mean=%f, sigma=%f, crate=%i, slot=%i, channel=%i) Setting sigma = 250\n", 
		       locDACValue, locMean, locPedestalSigma, HOSTROC, VMESlot, j);
		//locDACValue = 0; 
		locPedestalSigma = 250.0;
		//if (locMean < 0.0) locMean = 0.0;
	}
	else if ((locMean > (ReferenceBaseLine + 100))|| (locMean < 20)){
		printf("ERROR: locMean not within [20,%i] (DAC=%i, mean=%f, sigma=%f, crate=%i, slot=%i, channel=%i) Setting sigma = 250\n", 
		       ReferenceBaseLine + 100, locDACValue, locMean, locPedestalSigma, HOSTROC, VMESlot, j);
		locPedestalSigma = 250.0;
	}
      DACval<<HOSTROC<<"  "<<VMESlot<<" "<<j<<" "<<DACGain[k][j]<<"  "
	    <<locDACValue<<" "<<locMean<<"  "
	    << locPedestalSigma << endl;     
    }
  }

  DACval.close();

  cout<<"Total time used: "<<timelaps<<" s"<<endl;


  return 0;
}

//
//********************************************************************
//!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
//********************************************************************
//

void myReadoutCode(){
  //cout << "Entering myReadoutCode " << endl;
  DMANODE *outEvent;

  int iread = 0;
  int timeout = 100000;
  int dCnt = 0;
  int ReadyCode = 0;

  EventCounter++;
  GETEVENT(vmeIN,EventCounter);

  //cout << "Getting event " << EventCounter << endl;
  for(iread=0; iread<timeout; iread++) {
    if(fa125GBready()==fa125ScanMask())
      break;
  }

  if(iread==timeout) {
    ERRORCounter++;
    if (ERRORCounter<100){
      cout<<"ADC read TIMEOUT event: "<<EventCounter<<"     Return Mask is "
          << hex << ReadyCode<<"  should be "<< fa125ScanMask()<<dec<<endl;
    } else if (ERRORCounter%ERRORLOG){
      cout<<"ADC read TIMEOUT event: "<<EventCounter<<"     Return Mask is "
          << hex<<ReadyCode<<"  should be "<< fa125ScanMask()<<dec<<endl;
      ERRORLOG *= 10;
    }

  }  else {
    dCnt = fa125ReadBlock(0,(volatile UINT32 *)dma_dabufp, 0xF000, 0x2);
    if(dCnt<=0) {
      fa125ResetToken(0);
      if (ERRORCounter<100){
        printf("No fa125 data or error.  dCnt = %d\n",dCnt);
      } else if (ERRORCounter%ERRORLOG){
        ERRORLOG *=10;
        printf("No fa125 data or error.  dCnt = %d\n",dCnt);
      }
    } else {
      dma_dabufp += dCnt;
    }
  }

  fa125ResetToken(0);

  //cout<<"event all read out and ready to analyze"<<endl;
  //OUTF<<EventCounter<< "    Event Size: "<<dCnt<<endl;
  PUTEVENT(vmeOUT);
  outEvent = dmaPGetItem(vmeOUT);
  int eventLen = outEvent->length;

  int SLOTNUM = 0;
  int EVENT = -1;
  int CHANNEL = -1;
  int WSize = 0;
  int idx = 0;
  int pedcnt = 0;
  int adc1, adc2;
  int DATAReady=0;
  int EventsInBlock=0;
  int TriggerTime_1=0;
  int TriggerTime_2=0;
  int TrueSize=0;

  if (DEBUGGER){
    OUTF<<"EVENT "<<EventCounter<<"           Event Length: "<<eventLen<<endl; 
  }

  int DataSize = int(((float)NADCSAMPLES/2+1.7)*72.*(float)NADCS);

  DEBUG = 0;
  if ((eventLen<(DataSize-40)) || (eventLen>(DataSize+40))){
    //if ((eventLen<(NADCS*3866-40)) || (eventLen>(NADCS*3866+40))){
    //if ((eventLen<7730) || (eventLen>7800)){
    //if ((eventLen<950) || (eventLen>970)){
    OUTF<<"EVENT "<<EventCounter<<"           Event Length: "<<eventLen<<endl; 
    DEBUG = 1 ;
  } 

  int Indx = 0;

  //if ((EventCounter==56)|| (EventCounter==1)|| (EventCounter==111)|| 
  //    (EventCounter==167)){
  //  DEBUG = 1;
  //} else {
  //  DEBUG = 0;
  //}
 
  //cout<<"EVENTSIZE = "<<eventLen<<endl;
  for (Indx=0;Indx<eventLen;Indx++){
    unsigned int data = (unsigned int)LSWAP(outEvent->data[Indx]);

    //cout<<"loop index "<<Indx<<endl;

    if (DEBUG)
      OUTF<<Indx<<" "<<hex<<data<<" "<<dec<<"  :";


    if (data & 0x80000000) { // Block Header
      if (((data & 0xf8000000) >>27) == 0x10) { // Block Header
	SLOTNUM = ((data& 0x7C00000)>>22);
	if (SLOTNUM<11){
	  SLOTNUM -=3;
	} else {
	  SLOTNUM -= 5;
	}
	EventsInBlock = (data & 0x000FF);
	if (DEBUG)
	  OUTF<<"Block Header"<<endl;
      } else if (((data & 0xf8000000)>>27) == 0x12) { // Event Header
	EVENT = data & 0x003FFFFF;
	if (DEBUG)
	  OUTF<<"Event Header"<<endl;
      } else if (((data & 0xf8000000)>>27) == 0x13) { // Trigger Time
	TriggerTime_1 = (data & 0x00FFFFFF);
	data = (unsigned int)LSWAP(outEvent->data[Indx+1]);
	//OUTF<<Indx<<" "<<hex<<data<<" "<<dec<<endl;
	TriggerTime_2 = (data & 0x00FFFFFF);
	if (DEBUG)
	  OUTF<<"Trigger Time"<<endl;
      } else if (((data & 0xf8000000)>>27) == 0x14) { // Window Raw Data header
	CHANNEL = ((data & 0x7F00000)>>20); // flash is couning channels from 1 to 72 need 0 to 71
	WSize =  (data & 0xFFF);
	
	if (DEBUGGER){
	  if (DATAReady!=0){ 
	    OUTF<<"EVENT "<<EventCounter<<"           DATAREADY: "<<DATAReady<<endl; 
	    dmaPFreeItem(outEvent);
	    return;
	  }
	  if (WSize!=10){ 
	    OUTF<<"EVENT "<<EventCounter<<"           WINDOWSIZE: "<<WSize<<endl;
	    if (!DEBUG){
	      dmaPFreeItem(outEvent);
	      return;
	    }
	  }
	}
	//cout<<"Window Size: "<<WSize<<endl;
	if (NADCSAMPLES != WSize ) cout << "NADCSamples not equal WSize!" << NADCSAMPLES << " " << WSize << endl;
	DATAReady = WSize/2;
	idx = 0;
	Ped[SLOTNUM][CHANNEL] = 0.;
	pedcnt = 0;
	if (DEBUG)
	  OUTF<<"Window Raw Data Header, WSize is "<<WSize<<endl;
      } else if (((data & 0xf8000000)>>27) == 0x1d){// End of event
	if (DEBUG)
	  OUTF<<"End of EVENT"<<endl;
      } else if (((data & 0xf8000000)>>27) == 0x1f){
	if (DEBUG)
	  OUTF<<"Filler Word"<<endl;
      } else if (((data & 0xf8000000)>>27) == 0x11){
	if (DEBUG)
	  OUTF<<"Block Trailer"<<endl;
      } else if (((data & 0xf8000000)>>27) == 0x18){
	if (DEBUG)
	  OUTF<<"0x18 Unknown type"<<endl;
      } else if (((data & 0xf8000000)>>27) == 0x19){
	if (DEBUG)
	  OUTF<<"Streaming RAW DATA!! WHY???"<<endl;
      } else if (((data & 0xf8000000)>>27) == 0x16){
	if (DEBUG)
	  OUTF<<"Puls RAW DATA!! WHY???"<<endl;
      } else if (((data & 0xf8000000)>>27) == 0x15){
	if (DEBUG)
	  OUTF<<"Window SUM!! WHY???"<<endl;
      } else if (((data & 0xf8000000)>>27) == 0x1e){
	if (DEBUG)
	  OUTF<<"DATA NOT VALID!! WHY???"<<endl;
      } else if (((data & 0xf8000000)>>27) == 0x17){
	if (DEBUG)
	  OUTF<<"PULSE INTEGRAL!! WHY???"<<endl;
      } else if (((data & 0xf8000000)>>27) == 0x16){
	if (DEBUG)
	  OUTF<<"Puls RAW DATA!! WHY???"<<endl;
      } else if (((data & 0xf8000000)>>27) == 0x15){
	if (DEBUG)
	  OUTF<<"Window SUM!! WHY???"<<endl;
      } else if (((data & 0xf8000000)>>27) == 0x1a){
	if (DEBUG)
	  OUTF<<"UNKNOWN!! WHY???"<<endl;
      } else if (((data & 0xf8000000)>>27) == 0x1b){
	if (DEBUG)
	  OUTF<<"UNDEFINED TYPE!! WHY???"<<endl;
      } else if (((data & 0xf8000000)>>27) == 0x1c){
	if (DEBUG)
	  OUTF<<"SCALER HEADER!! WHY???"<<endl;
      }
    } else if (DATAReady>0) { // Window Raw Data values
      if (DEBUG){
	OUTF<<"The actual DATA!"<<endl;
      }
      DATAReady--; 
      ADCSamples[SLOTNUM][CHANNEL][idx++] =  (data & 0x1FFF0000) >> 16;
      ADCSamples[SLOTNUM][CHANNEL][idx++] =  (data & 0x1FFF);
      adc1 = ((data & 0x1FFF0000) >> 16);
      adc2 = (data & 0x1FFF);
      //cout<<SLOTNUM<<" "<<CHANNEL<<" "<<idx<<" "<<adc1<<"  "<<adc2<<endl;

	// This is a method to calculate the mean and sigma by updating the two values iteratively with each new measurement

	NSamples[SLOTNUM][CHANNEL] += 1;
	//cout << Q[SLOTNUM][CHANNEL] << " "<< NSamples[SLOTNUM][CHANNEL] << " "<< adc1  << " " << PedMean[SLOTNUM][CHANNEL] <<endl;
	Q[SLOTNUM][CHANNEL] = Q[SLOTNUM][CHANNEL] + ( ((float)NSamples[SLOTNUM][CHANNEL] - 1) / (float)NSamples[SLOTNUM][CHANNEL])*(adc1 - PedMean[SLOTNUM][CHANNEL])*(adc1 - PedMean[SLOTNUM][CHANNEL]);
      	PedMean[SLOTNUM][CHANNEL] = PedMean[SLOTNUM][CHANNEL] + (adc1 - PedMean[SLOTNUM][CHANNEL])/NSamples[SLOTNUM][CHANNEL];

	NSamples[SLOTNUM][CHANNEL] += 1;
	Q[SLOTNUM][CHANNEL] = Q[SLOTNUM][CHANNEL] + (((float)NSamples[SLOTNUM][CHANNEL] - 1)/(float)NSamples[SLOTNUM][CHANNEL])*(adc2-PedMean[SLOTNUM][CHANNEL])*(adc2-PedMean[SLOTNUM][CHANNEL]);
      	PedMean[SLOTNUM][CHANNEL] = PedMean[SLOTNUM][CHANNEL] + (adc2 - PedMean[SLOTNUM][CHANNEL])/NSamples[SLOTNUM][CHANNEL];

/*
      if (pedcnt<20){
	//cout << "adding "<< (adc1+adc2) << " to pedestal" <<endl;
	Ped[SLOTNUM][CHANNEL] += (adc1+adc2); 
	pedcnt += 2;
	if (pedcnt==20){
	  Ped[SLOTNUM][CHANNEL] /= 20.;
	  PedSum[SLOTNUM][CHANNEL] += Ped[SLOTNUM][CHANNEL];
	  PedSum2[SLOTNUM][CHANNEL] += (Ped[SLOTNUM][CHANNEL]*Ped[SLOTNUM][CHANNEL]);
	  PedSumCntr[SLOTNUM][CHANNEL] += 1.;
	  pedcnt = 0;
       	  //cout << "cntr = " << PedSumCntr[SLOTNUM][CHANNEL] << " pedSum = " << PedSum[SLOTNUM][CHANNEL] << endl;
	}
      }
*/
    } else {
      if (DEBUG)
	OUTF<<"Other stuff"<<endl;
    }
    
  } // end of loop over all data

  dmaPFreeItem(outEvent);
  if (DEBUG){
    OUTF<<"Event "<<EventCounter<<"            WSIZE =  "<<WSize<<endl;
  }
  
  if (RAWDATA){
    for (int j=0; j<NADCS;j++) {
      for (int i=0; i<72;i++) {
	OUTF<<EventCounter<<" "<<j<<" "<< i <<"  ";
	for (int k=0; k<WSize;k++){
	  //for (int k=0; k<8;k++){
	  OUTF<< ADCSamples[j][i][k]<< " "; 
	OUTF<<endl;
        }
      }
    }  
  }

  if (!(EventCounter%EVCOUNTER)){
    cout << "Calling Pedestal code" << endl;
    unsigned int time = (clock()-START)/CLOCKS_PER_SEC;
    cout<<time<<" s  Entering Pedestal Calculation :"<<EVENT<<" "<<EventCounter<<endl;
    mkpedestals(0);
    taskDelay(1);
    START = clock();
    timelaps += time;
  }

  //cout<<"END OF SUBROUTINE"<<endl;
}
//
//
//
void mkpedestals( int last){

  PedestalCounter++;
    
  for (int n=0;n<NADCS;n++){

    int VMESlot = n + 3;
    if (VMESlot > 10){
      VMESlot +=2;
    }

    for (int k=0;k<72;k++){

      float mu=0,r=0,sig=0;
/*
      if (PedSumCntr[n][k]>1){
	mu  = ((double)PedSum[n][k]) / PedSumCntr[n][k];
        sig = sqrt( (PedSumCntr[n][k] * (double)PedSum2[n][k] - (double)PedSum[n][k] * PedSum[n][k]) / (PedSumCntr[n][k]*(PedSumCntr[n][k]-1)) );
      }
*/
	//cout << Q[n][k] << "<=Q[n][k] "<< mu << " NSamples[n][k]=>"<<NSamples[n][k] << endl;
        mu = PedMean[n][k];
	sig = sqrt(Q[n][k] / (NSamples[n][k] - 1));
        PedestalSigma[n][k] = sig;
        CurrentMean[n][k] = mu;

	PedMean[n][k] = 0;
	Q[n][k] = 0;
	NSamples[n][k] = 0;


     // OUTFPed<< std::setprecision(3) << std::fixed <<PedestalCounter<<"\t"<<VMESlot<<"\t"<<k<<"\t"
	//     <<mu<<"\t\t"<<sig<<"\t0x"<<hex<<DACValues[n][k]<<dec
	//     <<"\t\t"<<BaseLine[n][k]<<"\t\t"<<DACGain[n][k]<< endl;

      if (!last && !BaselineConverged[n][k]){

	//if(!BaselineInRange[n][k]){
	if (!BaselineInRange[n][k] || mu == 0.0 ){
	  //cout << " mu = " << mu << endl;
	  if( mu > MINMEAN && mu < MAXMEAN){
	    //cout << "Baseline in range!" << endl;
	        NInRange++;
                BaselineInRange[n][k] = true;
                PreviousDACValue[n][k] = CurrentDACValue[n][k];
                PreviousMean[n][k] = CurrentMean[n][k];
                CurrentDACValue[n][k] -= 1000 / DETECTOR_FACTOR;
          }
	  else {
		PreviousDACValue[n][k] = CurrentDACValue[n][k];
                CurrentDACValue[n][k] += 2000 / DETECTOR_FACTOR;
		PreviousMean[n][k] = CurrentMean[n][k];
	  }

	  BaseLine[n][k] = ReferenceBaseLine;
	  int S = fa125Slot(n);
	  fa125SetOffset(S,k,CurrentDACValue[n][k]);	
	  continue;
	} 

	else {
	  float dped = float(PreviousMean[n][k] - CurrentMean[n][k]);
          float ddac = float(PreviousDACValue[n][k] - CurrentDACValue[n][k]);

          //if ( (fabs(dped) < 0.01 &&  fabs(mu - BaseLine[n][k]) < 50) || fabs(mu - BaseLine[n][k]) < 1.0) {
	  if ( (fabs(dped) < 0.4 &&  fabs(mu - BaseLine[n][k]) < 50)){
	  //if ( (fabs(dped) < 0.5 &&  mu > 50)) {
	    if (( mu > (ReferenceBaseLine + 100))|| (mu < 20)) cout << "Pedestal Converged far from target: dped " << dped << " ddac " << ddac << endl;
	    BaselineConverged[n][k]=true; NConverged++; 
	    PreviousMean[n][k] = CurrentMean[n][k];
	    PreviousDACValue[n][k] = CurrentDACValue[n][k];
	    continue;
	  }

	//cout << ddaq << "<=ddaq dped=>" << dped << endl;
	  float dacGain;
	  if (dped!= 0.0){
	    dacGain = ddac/dped;
	  } 
	  else{
	    dacGain = 1.;
	  }
	  if(fabs(dacGain - DACGain[n][k]) > 10 && DACGain[n][k] != 0.0) dacGain = DACGain[n][k];
          else if(ddac == 1000 / DETECTOR_FACTOR ){
                DACGain[n][k] = dacGain;
          }

	  if (dacGain!=0.0 && mu != 4095){
	    float newval = ((mu - (float)BaseLine[n][k]) * dacGain) ;
	    //cout<<n<<"    "<<k<<"   0x"<<hex<<DACValues[n][k] <<"    "<<dec<<newval
	    //	<<"      "<< (int)(mu+0.5)<<"    "<<BaseLine[n][k]<<"   "<<DACGain[n][k]<<endl;
	    int newDAC = CurrentDACValue[n][k] - int(newval);

	    PreviousDACValue[n][k] = CurrentDACValue[n][k];
	    CurrentDACValue[n][k] = newDAC;

	    int S = fa125Slot(n);
	    fa125SetOffset(S,k,CurrentDACValue[n][k]);
	    PreviousMean[n][k] = CurrentMean[n][k];
	  } 
	  else{
            PreviousMean[n][k] = CurrentMean[n][k];
            PreviousDACValue[n][k] = CurrentDACValue[n][k];
	    CurrentDACValue[n][k] -= 1500 / DETECTOR_FACTOR;
	    int S = fa125Slot(n);
            fa125SetOffset(S,k,CurrentDACValue[n][k]);
	  }
	} 
      }
     
    }
  }
  cout << NInRange   << "/" << NADCS*72 << " channels in range" << endl;
  cout << NConverged << "/" << NADCS*72 << " channels converged" << endl;
}