#include "Digitizer.h"

void Digitizer::DigitizeAC() {
  // created:  J. Conrad, KTH
  // modified: S. Orsi, INFN Roma2
  // fDataAC[0-63]:   main AC board
  // fDataAC[64-127]: extra AC board (identical to main board, for now)
  // We activate all branches. Once the digitization algorithm is determined 
  // only the branches that involve needed information will be activated

  // Threshold: thr=0.8 MeV.
  Float_t thr = 8e-4;
  Int_t nReg=6,i;
  fDataAC[0] = 0xACAC;
  fDataAC[64]= 0xACAC;
  fDataAC[1] = 0xAC11;   
  fDataAC[65] = 0xAC22;  
  
  // the third word is a status word (dummy: "no errors are present in the AC boards")
  fDataAC[2] = 0xFFFF; //FFEF?
  fDataAC[66] = 0xFFFF;
  // FPGA Registers (dummy)
  for (i=0; i<=nReg; i++){
    fDataAC[i+4] = 0xFFFF;
    fDataAC[i+68] = 0xFFFF;
  }
  // the last word is a CRC
  // Dummy for the time being, but it might need to be calculated in the end
  fDataAC[63] = 0xABCD;
  fDataAC[127] = 0xABCD;

  // shift registers (moved to the end of the routine)

  Int_t evntLSB=(UShort_t)Ievnt;
  Int_t evntMSB=Ievnt >> 16;

  // singles counters are dummy
  for(i=0; i<=15; i++){  //SO Oct '07:
    fDataAC[i+26] = evntLSB;
    fDataAC[i+90] = evntLSB;
  }
  for(i=0; i<=7; i++){
    fDataAC[i+42] = evntLSB;
    fDataAC[i+106] = evntLSB;
  }
  // increments for every trigger might be needed at some point.
  // dummy for now
  fDataAC[50]  = 0x0000;
  fDataAC[114] = 0x0000;

  // dummy FPGA clock (increment by 1 at each event)
  if(Ievnt<=0xFFFF){
    fDataAC[51] = 0x0000;
    fDataAC[52] = Ievnt;
    fDataAC[115] = 0x0000;
    fDataAC[116] = Ievnt;
  }
  else{
    fDataAC[51] = evntMSB;
    fDataAC[52] = evntLSB;
    fDataAC[115] = fDataAC[51];
    fDataAC[116] = fDataAC[52];
  }
  // dummy temperatures
  fDataAC[53] = 0x0000;
  fDataAC[54] = 0x0000;
  fDataAC[117] = 0x0000;
  fDataAC[118] = 0x0000;
  // dummy DAC thresholds
  for(i=0; i<=7; i++){
    fDataAC[i+55] = 0x1A13;  
    fDataAC[i+119] = 0x1A13;
  }
  // In this simpliefied approach we will assume that once a particle releases > 0.5 mip in one of the 12 AC detectors it will fire. We will furthermore assume that both cards read out identical data.

  // If you develop your digitization algorithm, you should start by identifying the information present in level2 (post-darth-vader) data. 

  Float_t SumEcat[4],SumEcas[4],SumEcard[4];
  for(i=0;i<4;i++){
    SumEcat[i]=0.;
    SumEcas[i]=0.;
    SumEcard[i]=0.;
  }
  // energy dependence on position (see file AcFitOutputDistancePmt.C by S.Orsi)
  // based on J.Lundquist's calculations (PhD thesis, page 94)
  // function: [0]+[1]*atan([2]/(x+1)), where the 3 parameters are:
  //   8.25470e-01   +- 1.79489e-02
  //   6.41609e-01   +- 2.65846e-02
  //   9.81177e+00   +- 1.21284e+00
  // hp: 1 minimum ionising particle at 35cm from the PMT releases 1mip

  // PMT positions: x,y,z: (average position of the 2 PMTs)
  Float_t posCasPmt[4][3]={{28.308, -17.168, 63.644},    // 1 - CAS CPU: x,y,z
			   {18.893, 24.913, 63.644},     // 2 - CAS DCDC 
			   {-24.307, 17.162, 63.644},    // 3 - CAS VME  
			   {-17.765, -28.300, 63.644}};  // 4 - CAS IPM  
  Float_t dAC=0.; // distance from PMT

  if(Nthcat>*ncat){
    cout<<"*** ERROR AC! Nthcat=  "<<Nthcat<<" out of range! "<<endl;
    for(i=0;i<4;i++)SumEcat[i]=2*thr;
  }
  else{
    for(i=0;i<Nthcat;i++){
      if(Icat[i]>0 && Icat[i]<5)SumEcat[Icat[i]-1]+=Erelcat[i]; 
    }
  }
  if(Nthcas>*ncas){
    cout<<"*** ERROR AC! Nthcas=  "<<Nthcas<<" out of range!"<<endl;
    for(i=0;i<4;i++)SumEcas[i]=2*thr;
  }
  else{
    for (i=0;i<Nthcas;i++){
      if(Icas[i]>0 && Icas[i]<5){
	dAC=sqrt(pow((Xincas[i]+Xoutcas[i])/2 - posCasPmt[Icas[i]-1][0],2) + pow((Yincas[i]+Youtcas[i])/2 - posCasPmt[Icas[i]-1][1],2) + pow((Zincas[i]+Zoutcas[i])/2 - posCasPmt[Icas[i]-1][2],2));
	SumEcas[Icas[i]-1] += Erelcas[i]*attenAC->Eval(dAC);
      }
    }
  }
  if(Nthcard>*ncar){
    cout<<"*** ERROR AC! Nthcard= "<<Nthcard<<" out of range!"<<endl;
    for(i=0;i<4;i++)SumEcard[i]=2*thr;
  }
  else{
    for(Int_t k= 0;k<Nthcard;k++){
      if(Icard[i]>0 && Icard[i]<5)SumEcard[Icard[k]-1] += Erelcard[k];
    }
  }
 
  // channel mapping              Hit Map
  // 1 CARD4                      0          LSB
  // 2 CAT2                       0
  // 3 CAS1                       0
  // 4 NC                         0
  // 5 CARD2                      0
  // 6 CAT4                       1 
  // 7 CAS4                       0  
  // 8 NC                         0
  // 9 CARD3                      0
  // 10 CAT3                      0
  // 11 CAS3                      0 
  // 12 NC                        0
  // 13 CARD1                     0
  // 14 CAT1                      0
  // 15 CAS2                      0
  // 16 NC                        0          MSB

  // In the first version only the hit-map is filled, not the SR.

  fDataAC[3] = 0x0000;

  if (SumEcas[0] > thr)  fDataAC[3]  = 0x0004;
  if (SumEcas[1] > thr)  fDataAC[3] += 0x4000;
  if (SumEcas[2] > thr)  fDataAC[3] += 0x0400;
  if (SumEcas[3] > thr)  fDataAC[3] += 0x0040;   

  if (SumEcat[0] > thr)  fDataAC[3] += 0x2000;
  if (SumEcat[1] > thr)  fDataAC[3] += 0x0002;
  if (SumEcat[2] > thr)  fDataAC[3] += 0x0200;
  if (SumEcat[3] > thr)  fDataAC[3] += 0x0020; 
  
  if (SumEcard[0] > thr)  fDataAC[3] += 0x1000;
  if (SumEcard[1] > thr)  fDataAC[3] += 0x0010;
  if (SumEcard[2] > thr)  fDataAC[3] += 0x0100;
  if (SumEcard[3] > thr)  fDataAC[3] += 0x0001; 

  fDataAC[67] = fDataAC[3];

  // shift registers
  // the central bin is equal to the hitmap, all other bins in the shift register are 0
  for (UInt_t i=0; i<=15; i++){
    fDataAC[i+11] = 0x0000;   
    fDataAC[i+75] = 0x0000;
  }
  fDataAC[18] = fDataAC[3];
  fDataAC[82] = fDataAC[3];

  //    for (Int_t i=0; i<fACbuffer; i++){
  //    printf("%0x  ",fDataAC[i]);   
  //    if ((i+1)%8 ==0) cout << endl;
  //   }
};