// HiLumUtils
// 08 November 2009
// 08 Nov 2009 Augmented for Nov 2009 run
// 23 Apr 2010 Add Cenerkov event finder and reader
// and update to newer header format on Cerenkov files
// 27 Mar 2013 Update unix time stamp
// 19 Apr 2013 Accommodate various bunch structures
//
// Utility routines to be used by HiLum analysis programs
//
// ***************************************************
// Convert unix time stamp into Yr/Mo/Day/Hr/Min/Sec integers
// for a clock in Protvino on Standard Time, not Daylight Savings Time

void TimeStamp(int TS, unsigned int* pntr) {
  int dyr=0; // days in the year
  int Mth[]={31,28,31,30,31,30,31,31,30,31,30,31};
  int UTSref=1199134800; // 2008/01/01 00:00:00 in Protvino
  int df=TS-UTSref;  // seconds since UTSref time
  int dy=df/24/3600; // days since UTSref time
  int Year=2008;
  int Month=1;
  int leap;

  while(1) {
    leap=((Year%4==0)&&(Year%100!=0))?1:0;
    dyr=365+leap;
    if(dy<dyr) {break;}
    dy-=dyr; df-=dyr*24*3600; Year++;
  }
  *(pntr++)=Year;
  Mth[1]+=leap;
  for(int i=0;i<12;i++) {
    if(dy<Mth[i]) {break;}
    dy-=Mth[i]; df-=Mth[i]*24*3600; Month++; 
  }
  *(pntr++)=Month;
  int Day=df/24/3600;
  *(pntr++)=1+Day;
  df-=Day*24*3600;
  int Hour=df/3600;
  *(pntr++)=Hour;
  df-=Hour*3600;
  int Min=df/60;
  *(pntr++)=Min;
  df-=Min*60;
  int Sec=df;
  *(pntr++)=Sec;
  return;
}
// ***************************************************
// Read in Kozelov Cerenkov data, find pointers
// to event flags, and summarize history

unsigned int nw16=0; // Number of Cerenkov ADC words in file 
// Cherenkov file header data
short	voltm_type;	// 11
int	voltm_length;
double	voltm_data[8];
  //	0	unused
  //	1	Cerenkov HV
  //	2	S1 HV
  //	3	S2 HV
  //	4	S3 HV
  //	5	S4 HV
  //	6	S5 HV
  //	7	S6 HV
short	time_type;	// 12
int	time_length;
int	time_data;
unsigned int ymdhms[6];
short	adcbuf_type;	// 15
const int BufSize=1024*1024*16; 
unsigned short Cbuf[BufSize]; // Cerenkov ADC data 
unsigned  char Cflg[BufSize],Covf[BufSize]; // Begin-event flag, Overflow flag 
int nCovf=0; // Number of ADC overflows
const int nSp=5; // Number of spills (in order starting from the beginning
                 // of the run) with C data
char* pfn[nSp];
// Data relative to the event
  const int nEvt=1000;    // Max number of events in spill
  int EvtPnt[nEvt];       // Point to beginning of bunches with event flag set
  int TnPntr[nEvt];       // Point to beginning of first turn in event
  //int PkPntr[nEvt];       // Point to first isolated peak in event
  int PdPntr[nEvt];       // Pointer to pedestal in event
  int PedC[nEvt];         // Pedestal for Cerenkov in this event
  int NoCflg[nEvt];       // Number of readings with event flag set
  int iii[nEvt];          // Cbuf line number of EvtPnt (before offset)
  double CHist1[nEvt],CHist2[nEvt];
  int JBeg[nEvt],JEnd[nEvt];
  int nsum[nEvt];         // Number of bunches in history
  double AvePed,AvePedBeg,AvePedEnd,AvePed1,AvePedSqr,AvePedBegRMS; // Average value of pedestals
  int nPed,nPedBeg,nPedEnd; // Number of bunches in average Ped
  int kevents;           // Number of events in spill seen in C data
  const int nBnTrn=30;   // Number of bunches in one turn (or one lap) around U70
  int nTurns;            // Number of complete turns
  double HiBin[nBnTrn];  // Histogram of Cph minus ped versus bunch within turn
  int nModFlg;           // Number flagged as isolated filld bunches in a turn
  int mModFlg;           // Number flagged as filled bunches in a turn
  int ModFlg[nBnTrn];    // Flag isolated bunches
  int iModFlg[15];       // Bunch number in a turn which is flagged as isolated
  int jModFlg[30];       // Bunch number in a turn which is flagged as filled
//
unsigned char Bite1,Bite2; 
unsigned short word1,word2;

void Cerenkov (int SpillNo, char* GateName) {
  // Read in Cerenkov data - Note that Cerenkov DAQ time is on Day-Light savings
  // time while Sepp's DAQ is on standard time.
  //char cerenDir [] = "/scratch/raid/HiLum/Apr10data0/Apr10Ceren/"; // directory 
  //char cerenFile00[] = "spill-2010-04-17-13-27-41.dat"; pfn[ 0]=cerenFile00; // goes with run 833
  char cerenDir [] = "/scratch/raid/HiLum/Mar13data0/Mar13Ceren/"; // directory 
  //char cerenFile00[] = ""; pfn[ 0]=cerenFile00; //
  //char cerenFile00[] = "spsim-2013-05-31-07-02-24.dat"; pfn[ 0]=cerenFile00; //
  //char cerenFile00[] = "spsim-2013-05-30-13-26-59.dat"; pfn[ 0]=cerenFile00; //
  //char cerenFile00[] = "spsim-2013-05-10-11-19-03.dat"; pfn[ 0]=cerenFile00; //
  //char cerenFile00[] = "spill-2013-03-24-00-54-37.dat"; pfn[ 0]=cerenFile00; // 
  //char cerenFile00[] = "spill-2013-03-24-03-09-53.dat"; pfn[ 0]=cerenFile00; // 
  //char cerenFile00[] = "spill-2013-03-24-03-10-03.dat"; pfn[ 0]=cerenFile00; // 
  //char cerenFile00[] = "spill-2013-03-24-03-10-13.dat"; pfn[ 0]=cerenFile00; // 
  //char cerenFile00[] = "spill-2013-03-24-03-10-23.dat"; pfn[ 0]=cerenFile00; // 
  //char cerenFile00[] = "spill-2013-03-24-00-54-37.dat"; pfn[ 0]=cerenFile00; // 
  //char cerenFile00[] = "spill-2013-03-25-16-31-05.dat"; pfn[ 0]=cerenFile00; // One bunch
  //char cerenFile00[] = "spill-2010-04-17-06-40-11.dat"; pfn[ 0]=cerenFile00; // goes with run 826
  //char cerenFile00[] = "spill-2010-04-16-21-16-51.dat"; pfn[ 0]=cerenFile00; // goes with run 801
  //char cerenFile00[] = "spill-2013-03-24-16-17-53.dat"; pfn[ 0]=cerenFile00; // goes with run 1102
  //char cerenFile00[] = "integ-2012-04-10-18-22-14.binB"; pfn[ 0]=cerenFile00; // Integrated a la Alisha
  char cerenFile00[] = "spill-2013-03-24-09-47-41.dat"; pfn[ 0]=cerenFile00; // goes with run 1092
  char cerenFile01[] = "spill-2013-03-25-17-06-38.dat"; pfn[ 1]=cerenFile01; // goes with run 1141
  char cerenFile02[] = "spill-2013-03-25-17-06-57.dat"; pfn[ 2]=cerenFile02;
  char cerenFile03[] = "spill-2013-03-25-17-07-07.dat"; pfn[ 3]=cerenFile03;
  char cerenFile04[] = "spill-2013-03-25-17-07-17.dat"; pfn[ 4]=cerenFile04;
  //char cerenFile01[] = "spill-2010-04-16-21-17-01.dat"; pfn[ 1]=cerenFile01;
  //char cerenFile02[] = "spill-2010-04-16-21-17-10.dat"; pfn[ 2]=cerenFile02;
  //char cerenFile03[] = "spill-2010-04-16-21-17-19.dat"; pfn[ 3]=cerenFile03;
  //char cerenFile04[] = "spill-2010-04-16-21-17-28.dat"; pfn[ 4]=cerenFile04;
  char cerenFull[120];
  strcpy(GateName,pfn[SpillNo-1]);
 
  if(SpillNo>nSp) {
    cout <<"No Cerenkov file for this spill"<<endl;
    return;
  }
  strcpy(cerenFull,cerenDir); strcat(cerenFull,pfn[SpillNo-1]); // <----**********
  cout << cerenFull<<endl; 
  std::ifstream inC; 
  inC.open(cerenFull,ios::binary); 
  if (!inC.good()) {cout <<"Cannot read Cerenkov file"<<endl; return;} 
  Bite1=inC.get(); Bite2=inC.get(); voltm_type=256*Bite2+Bite1;
  Bite1=inC.get(); Bite2=inC.get(); word1=256*Bite2+Bite1; 
  Bite1=inC.get(); Bite2=inC.get(); word2=256*Bite2+Bite1; 
  voltm_length=65536*word2+word1; 
  union mix_t {double D; unsigned short sc[4]; unsigned int ic[2];
    unsigned char c[8];} mix; 
  for(int j=0;j<8;j++) {
    //for(int i=0;i<4;i++) {mix.c[2*i+1]=inC.get(); mix.c[2*i]=inC.get();} 
    for(int i=0;i<4;i++) {mix.c[2*i]=inC.get(); mix.c[2*i+1]=inC.get();} 
    voltm_data[j]=mix.D; 
  }
  printf(" voltm_length=%3d\n",voltm_length);
  for(int j=0;j<8;j++) {
    printf("%1d %17.12f\n",j,voltm_data[j]);
  }
  Bite1=inC.get(); Bite2=inC.get(); time_type=256*Bite2+Bite1;
  Bite1=inC.get(); Bite2=inC.get(); word1=256*Bite2+Bite1; 
  Bite1=inC.get(); Bite2=inC.get(); word2=256*Bite2+Bite1; 
  time_length=65536*word2+word1; 
  Bite1=inC.get(); Bite2=inC.get(); word1=256*Bite2+Bite1; 
  Bite1=inC.get(); Bite2=inC.get(); word2=256*Bite2+Bite1; 
  time_data=65536*word2+word1; 
  TimeStamp(time_data,ymdhms);  
  printf("time_length/data = %12d %12d\n",time_length,time_data);
  printf("%4d.%02d.%02d %02d:%02d:%02d\n",
    ymdhms[0],ymdhms[1],ymdhms[2],ymdhms[3],ymdhms[4],ymdhms[5]);
//****************** Test
  UInt_t YMDHMS[6];
  time_t t=time_data;
  TTimeStamp *tm = new TTimeStamp(t,0);
  tm->GetDate(kTRUE,0,&YMDHMS[0],&YMDHMS[1],&YMDHMS[2]);
  tm->GetTime(kTRUE,0,&YMDHMS[3],&YMDHMS[4],&YMDHMS[5]);
  printf("%4d.%02d.%02d %02d:%02d:%02d\n",
    YMDHMS[0],YMDHMS[1],YMDHMS[2],YMDHMS[3],YMDHMS[4],YMDHMS[5]);
//****************** End test
  Bite1=inC.get(); Bite2=inC.get(); adcbuf_type=256*Bite2+Bite1;
  Bite1=inC.get(); Bite2=inC.get(); word1=256*Bite2+Bite1; 
  Bite1=inC.get(); Bite2=inC.get(); word2=256*Bite2+Bite1; 
  nw16=65536*word2+word1; // adcbuf_length
  nTurns=nw16/nBnTrn; // Number of complete turns
  cout <<"Getting Cerenkov Data, nw16 = "<<nw16<<endl; 
  int k=0, kk, kMax, jbeg, jend;
  int i,Ediff;
  int Offset=25; // Bunches to go backwards to align with MPI samples
  for(unsigned int ii=0;ii<nw16;ii++) { 
    i=ii;
    Bite1=inC.get(); Bite2=inC.get();
    Cbuf[i]=((256*Bite2+Bite1)&0x0FFF); 
    Cflg[i]=(Bite2/16)&0x0001; Covf[i]=(Bite2/32)&0x0001;
    if(Covf[i]>0) {
      Cbuf[i]=4095; // *************************************
      nCovf++;
    }
    //if(ii<32) printf("ii Cbuf[ii]: %8d %4d\n",ii,Cbuf[ii]);
    if(((i==0)&&(Cflg[0]==1))||((i>0)&&(Cflg[i]==1)&&(Cflg[i-1]==0))) {
      EvtPnt[k]=i; // Mark event
      iii[k]=i;    // Mark line number
      // Go backwards to bunches corresponding to MPI samples
      if(EvtPnt[k]>Offset) {EvtPnt[k]-=Offset;}
      k++;
    }
  }
  printf("End of %s with %8d ADC overflows out of %10d bunches\n",pfn[SpillNo-1],nCovf,nw16);
  inC.close();
  kevents=k;

  AvePed=0.; AvePedBeg=0.; AvePedEnd=0.; AvePedSqr=0; nPed=0; nPedBeg=0; nPedEnd=0;
  // First way to get pedestals
  int pedevent=kevents-1;
  if(kevents>6) pedevent=6;
  jend=EvtPnt[pedevent]; jbeg=jend-240000; if(jbeg<0) jbeg=0; // 4 events (40 ms) wide region
  for(int j=jbeg;j<jend;j++) {
    AvePedBeg+=double(Cbuf[j]);
    AvePedSqr+=pow(double(Cbuf[j]),2);
    nPedBeg++;
  }
  AvePedBeg/=double(nPedBeg); AvePedSqr/=double(nPedBeg);
  AvePedBegRMS=sqrt(AvePedSqr-AvePedBeg*AvePedBeg);

  // Find isolated bunches
  for(int j=0;j<nBnTrn;j++) {HiBin[j]=0.0; ModFlg[j]=0;}
  for(unsigned int ii=0;ii<nw16;ii++) {HiBin[ii%nBnTrn]+=double(Cbuf[ii])-AvePedBeg;}
  int ModMxB=-1;
  double HiMax=-1.0E14;
  for(int j=0;j<nBnTrn;j++) {if(HiBin[j]>HiMax) {HiMax=HiBin[j]; ModMxB=j;}}
  for(int j=0;j<nBnTrn;j++) {if(HiBin[j]>0.25*HiMax) {ModFlg[j]=2;}}
  int jl,jj,jh;
  for(int j=nBnTrn;j<2*nBnTrn;j++) {
    jl=(j-1)%nBnTrn; jj=j%nBnTrn; jh=(j+1)%nBnTrn;
    if(ModFlg[jl]==0&&ModFlg[jj]==2&&ModFlg[jh]==0) ModFlg[jj]=1;
  }
  nModFlg=0; mModFlg=0;
  for(int j=0;j<nBnTrn;j++) {
    if(ModFlg[j]==1) {
      iModFlg[nModFlg]=j;
      nModFlg++;
    }
    if(ModFlg[j]>0) {
      jModFlg[mModFlg]=j;
      mModFlg++;
    }
  }
  printf("\n Isolated filled bunches\n nModFlg=%2d\n i iModFlg\n",nModFlg);
  for(int j=0;j<nModFlg;j++) {printf("%2d   %2d\n",j,iModFlg[j]);}
  printf("\n All filled bunches\n mModFlg=%2d\n j jModFlg\n",mModFlg);
  for(int j=0;j<mModFlg;j++) {printf("%2d   %2d\n",j,jModFlg[j]);}
  printf("\n  i       HiBin ModFlg\n");
  for(int j=0;j<nBnTrn;j++) {
    printf("%3d %12.0f %2d\n",j,HiBin[j],ModFlg[j]);
  }

  // Find event-related pointers
  for(int ki=0;ki<kevents;ki++) {
    k=ki;
    // point to beginning of first full turn after event flag
    // The beginning of a turn corresponds to the first bunch (filled or not)
    // in the gated Cerenkov file and to every bunch stepping by 30 bunches
    kMax=(EvtPnt[k]/nBnTrn)*nBnTrn;  if(kMax<EvtPnt[k]) {kMax+=nBnTrn;}
    TnPntr[k]=kMax; 
    //PdPntr[k]= (kMax>EvtPnt[k]) ? kMax-1: kMax+2; PedC[k]=Cbuf[PdPntr[k]];
    PdPntr[k]=0; PedC[k]=0;
    // Point to first isolated filled bunch in turn after event flag
    // and then take the empty bunch just before it as candidate for pedestal
    if(nModFlg>0) {
      PdPntr[k]=TnPntr[k]+iModFlg[0]-1; PedC[k]=Cbuf[PdPntr[k]];
    }
    // Calculate pedestal values
    // Second way to get pedestals
    AvePed+=double(PedC[k]); nPed++;
  }
  printf("\nAvePedBeg integration region, jbeg, jend %8d %8d %8d\n\n",jbeg,jend,nPedBeg);
  // Third way to get pedestals
  jbeg=EvtPnt[kevents-2]; jend=jbeg+500;
  for(int j=jbeg;j<jend;j++) {AvePedEnd+=double(Cbuf[j]); nPedEnd++;}  
  AvePed/=double(nPed);
  AvePedEnd/=double(nPedEnd);
  AvePed1=(AvePedBeg+AvePedEnd)/2.;
  printf("Pedestals calculated different ways:\n");
  printf(" %8.1f At beginning of spill\n",AvePedBeg);
  printf(" %8.1f RMS of ped at beginning of spill\n",AvePedBegRMS);
  printf(" %8.1f Average ped in events\n",AvePed);
  printf(" %8.1f At end of spill\n",AvePedEnd);
  printf(" %8.1f Ave of Beginning and End\n\n",AvePed1);

  cout <<" Found "<<kevents<<" Events.   Ave Ped Beg End 1 RMS = " 
       <<AvePed<<" "<<AvePedBeg<<" "<<AvePedEnd<<" "<<AvePed1<<" "<<AvePedBegRMS<<endl; 

  printf("   k   EvtPnt  Ediff    jbeg jend Cbuf[isolat] PedC  nsum   CHist1\n"); 
  // int kback=8000*nBnTrn; // For ionization history, go back 40 ms 
  int kback=1600*nBnTrn; // For ionization history, go back 8 ms 
  for(int ki=0;ki<kevents;ki++) {
    k=ki;
    // Construct ionization history before EvtPnt[k]
    // Go back 33 ms for V=1200 across 2 mm gap
    jend=TnPntr[k]; if(jend<0) {jend=0;}
    jbeg=jend-kback;
    if(jbeg<0) {jbeg=0;}
    nsum[k]=0; CHist1[k]=0.;
    JBeg[k]=jbeg; JEnd[k]=jend;
    if(jend>=nBnTrn) {
      for(int j=jbeg;j<=jend;j++) {
        nsum[k]++; CHist1[k]+=double(Cbuf[j])-AvePedBeg;
      }
    }
  }
  for(int ki=0;ki<kevents;ki++) {
    k=ki;
    kk=(k>0) ? k-1:k;
    Ediff=EvtPnt[k]-EvtPnt[kk];
    if(nsum[k]>0) {
      CHist1[k]/=double(nsum[k]*mModFlg)/double(nBnTrn);
    }
    printf(" %3d %8d %6d %7d %7d %6d %6d %6d %8.0f\n",
           k,EvtPnt[k],Ediff,JBeg[k],JEnd[k],Cbuf[TnPntr[k]+iModFlg[0]],PedC[k],nsum[k],CHist1[k]);
  }  
  return;
}
// ***************************************************
int GetCData(int burst, int EvtNo) {
  int EP;
  EP=-1; if((burst>0)&&(burst<=nSp)) {EP=EvtPnt[EvtNo];}
  return EP;
}
// ***************************************************