#include <sstream>
#include <fstream>
#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include <ctype.h>
#include <time.h>
#include "Riostream.h"
#include "TFile.h"
#include "TDirectory.h"
#include "TTree.h"
#include "TLeafI.h"
#include "TH1.h"
#include "TH2.h"
#include "TF1.h"
#include "TMath.h"
#include "TRandom.h"
#include "TSQLServer.h"
#include "TSystem.h"
#include "CalibTrk1Event.h"
#include "CalibTrk2Event.h"
//
#include "Digitizer.h"
#include "CRC.h"
//
#include <PamelaRun.h>
#include <physics/calorimeter/CalorimeterEvent.h>
#include <CalibCalPedEvent.h>
#include "GLTables.h"


void Digitizer::DigitizeAC(TF1 *attenAC) {
  // 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
  
  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;

  const UInt_t nReg = 6;

  // FPGA Registers (dummy)
  for (UInt_t 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=Ievnt%65536;
  //Int_t evntMSB=(Int_t)(Ievnt/65536);
  Int_t evntLSB=(UShort_t)Ievnt;
  Int_t evntMSB=Ievnt >> 16;

  // singles counters are dummy
  for (UInt_t i=0; i<=15; i++){  //SO Oct '07: // for (UInt_t i=0; i<=16; i++){
    //     fDataAC[i+26] = 0x0000;   
    //     fDataAC[i+90] = 0x0000;
    fDataAC[i+26] = evntLSB;
    fDataAC[i+90] = evntLSB;
  };
  
  for (UInt_t 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 (UInt_t 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[5];
  Float_t SumEcas[5];
  Float_t SumEcard[5];
  for (Int_t k= 0;k<5;k++){
    SumEcat[k]=0.;
    SumEcas[k]=0.;
    SumEcard[k]=0.;
  };

  if (Nthcat>50 || Nthcas>50 || Nthcard>50)
    printf("*** ERROR AC! NthAC out of range!\n\n");

  // 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

  //TF1 *attenAC = new TF1("fAttAC",".825+.64*atan(9.8/x)",0.,45.);

  // 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
  
  // look in CAT
  //  for (UInt_t k= 0;k<50;k++){
  for (Int_t k= 0;k<Nthcat;k++){
    if (Erelcat[k] > 0)
      SumEcat[Icat[k]] += Erelcat[k]; 
  };

  // look in CAS
  for (Int_t k= 0;k<Nthcas;k++){
    if (Erelcas[k] >0) {
      dAC=sqrt(pow((Xincas[k]+Xoutcas[k])/2 - posCasPmt[Icas[k]-1][0],2) + pow((Yincas[k]+Youtcas[k])/2 - posCasPmt[Icas[k]-1][1],2) + pow((Zincas[k]+Zoutcas[k])/2 - posCasPmt[Icas[k]-1][2],2));
      SumEcas[Icas[k]] += Erelcas[k]*attenAC->Eval(dAC);
    }
   };
 // attenAC->Delete();	
  // look in CARD
  for (Int_t k= 0;k<Nthcard;k++){
    if (Erelcard[k] >0)
      SumEcard[Icard[k]] += 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.

  // Threshold: 0.8 MeV.

  Float_t thr = 8e-4;

  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;
  //   }
};