#include "PamVMCAcDig.h"
#include <TF1.h>
#include <TMath.h>

using TMath::Power;
using TMath::Sqrt;

ClassImp(PamVMCAcDig)

void PamVMCAcDig::DigitizeAC(Int_t EventNo){

  // cout<<"Digitizing AC... "<<endl;

  fDEBUG = kTRUE;
  
  for (Int_t k= 0;k<5;k++) fSumEcat[k]=fSumEcas[k]=fSumEcard[k]=0.;
  //Zeroing arrays holds energy releases in AC

  UShort_t DataAC[128];
  for(Int_t i=0; i<128; i++) DataAC[i]=0;
    

  // created:  J. Conrad, KTH
  // modified: S. Orsi, INFN Roma2
  // DataAC[0-63]:   main AC board
  // DataAC[64-127]: extra AC board (identical to main board, for now)

  DataAC[0] = 0xACAC;
  DataAC[64]= 0xACAC;
  DataAC[1] = 0xAC11;   
  DataAC[65] = 0xAC22;  

  // the third word is a status word 
  //(dummy: "no errors are present in the AC boards")
  DataAC[2] = 0xFFFF; //FFEF?
  DataAC[66] = 0xFFFF;

  const UInt_t nReg = 6;

  // FPGA Registers (dummy)
  for (UInt_t i=0; i<=nReg; i++){
    DataAC[i+4] = 0xFFFF;
    DataAC[i+68] = 0xFFFF;
  }

  // the last word is a CRC
  // Dummy for the time being, but it might need 
  // to be calculated in the end
  DataAC[63] = 0xABCD;
  DataAC[127] = 0xABCD;

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

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

  // singles counters are dummy
  for (UInt_t i=0; i<=15; i++){  
    DataAC[i+26] = evntLSB;
    DataAC[i+90] = evntLSB;
  };
  
  for (UInt_t i=0; i<=7; i++){
    DataAC[i+42] = evntLSB;
    DataAC[i+106] = evntLSB;
  };


  // increments for every trigger might be needed at some point.
  // dummy for now
  DataAC[50]  = 0x0000;
  DataAC[114] = 0x0000;

  // dummy FPGA clock (increment by 1 at each event)
  if (EventNo<=0xFFFF) {
    DataAC[51] = 0x0000;
    DataAC[52] = EventNo;
    DataAC[115] = 0x0000;
    DataAC[116] = EventNo;
  } else {
    DataAC[51] = evntMSB;
    DataAC[52] = evntLSB;
    DataAC[115] = DataAC[51];
    DataAC[116] = DataAC[52];
  }

  // dummy temperatures
  DataAC[53] = 0x0000;
  DataAC[54] = 0x0000;
  DataAC[117] = 0x0000;
  DataAC[118] = 0x0000;


  // dummy DAC thresholds
  for (UInt_t i=0; i<=7; i++){
    DataAC[i+55] = 0x1A13;  
    DataAC[i+119] = 0x1A13;
  }

  FillErel();
  
  // 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. 

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

  const Float_t thr = 8e-4; //(GeV)

  //Important Note: In simulation model we have
  //CAT as solid plasic plate. So, I set values:

  fSumEcat[1]=fSumEcat[2]=fSumEcat[3]=fSumEcat[0];
  
  DataAC[3] = 0x0000;

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

  if (fSumEcat[0] > thr)  DataAC[3] += 0x2000; //We have solid CAT scintillator
  if (fSumEcat[1] > thr)  DataAC[3] += 0x0002; //See note above
  if (fSumEcat[2] > thr)  DataAC[3] += 0x0200; //
  if (fSumEcat[3] > thr)  DataAC[3] += 0x0020; //
  
  if (fSumEcard[0] > thr) DataAC[3] += 0x1000; //C1D1(1)
  if (fSumEcard[1] > thr) DataAC[3] += 0x0010; //C1D1(2)
  if (fSumEcard[2] > thr) DataAC[3] += 0x0100; //C2D1(1)
  if (fSumEcard[3] > thr) DataAC[3] += 0x0001; //C2D1(2)

  DataAC[67] = DataAC[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++){
    DataAC[i+11] = 0x0000;   
    DataAC[i+75] = 0x0000;
  }
  DataAC[18] = DataAC[3];
  DataAC[82] = DataAC[3];

  //++++++ WRITE EVERYTHING TO DIGITIZER'S BUFFER +++//
  fData.clear(); //clearing data buffer

  for(Int_t i= 0; i<128; i++) fData.push_back(DataAC[i]);


  if (fDEBUG){
    cout<<"====== AC hit info ======"<<endl;
    cout<<"CAT:   Erel="<<fSumEcat[0]*1.e6<<" (keV)"<<endl
	<<"CARD1: Erel="<<fSumEcard[0]*1.e6<<" (keV)"<<endl
	<<"CARD2: Erel="<<fSumEcard[1]*1.e6<<" (keV)"<<endl
	<<"CARD3: Erel="<<fSumEcard[2]*1.e6<<" (keV)"<<endl
	<<"CARD4: Erel="<<fSumEcard[3]*1.e6<<" (keV)"<<endl
	<<"CAS1:  Erel="<<fSumEcas[0]*1.e6<<" (keV)"<<endl
	<<"CAS2:  Erel="<<fSumEcas[1]*1.e6<<" (keV)"<<endl
	<<"CAS3:  Erel="<<fSumEcas[2]*1.e6<<" (keV)"<<endl
	<<"CAS4:  Erel="<<fSumEcas[3]*1.e6<<" (keV)"<<endl;
    cout<<"========================="<<endl;

  }
}

void PamVMCAcDig::FillErel(){
  /*
    Detector     KEY          Number of paddles

    CAT          TOP1                  1 (solid volume with id=1)
    CAS          SID1                  4 (1-4)
    CARD         C1D1 & C2D1           2 & 2 (1-2)
  */
  TClonesArray* hc=0;
  const char* keyac[4] = {"TOP1","SID1","C1D1","C2D1"};
  for(Int_t i=0; i<4; i++){
    hc = (TClonesArray*)fhitscolmap.GetValue(keyac[i]);
    if(hc) AddErel(hc, i);
    hc = 0;
  }
}

void PamVMCAcDig::AddErel(TClonesArray* HitColl, Int_t detcode){


  // 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.);

  const 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 dx,dy,dz,dAC; // distance from PMT
  Int_t hitsNo=HitColl->GetEntriesFast();
  PamVMCDetectorHit * hit = 0;

  if (hitsNo>50) cout<<"WARNING: PamVMCAcDig::AddErel... "
                     <<"Number of hits in AC is out of range!!!"<<endl;

  switch(detcode){
  case 0: //CAT
    for(Int_t k=0; k<hitsNo; k++){
      hit = (PamVMCDetectorHit*)HitColl->At(k);
      fSumEcat[(hit->GetPOS())-1] += hit->GetEREL();
    }
    break;

  case 1: //CAS
    for(Int_t k=0; k<hitsNo; k++){
      hit = (PamVMCDetectorHit*)HitColl->At(k);
      dx = (hit->GetXIN()+hit->GetXOUT())/2.-posCasPmt[hit->GetPOS()-1][0];
      dy = (hit->GetYIN()+hit->GetYOUT())/2.-posCasPmt[hit->GetPOS()-1][1];
      dz = (hit->GetZIN()+hit->GetZOUT())/2.-posCasPmt[hit->GetPOS()-1][2];
      dAC = Sqrt(dx*dx + dy*dy + dz*dz);
      //fSumEcas[(hit->GetPOS())-1] += (hit->GetEREL())*attenAC->Eval(dAC);
      if ((hit->GetPOS()==3)) fSumEcas[3] += (hit->GetEREL())*fattenAC->Eval(dAC);
      if ((hit->GetPOS()==4)) fSumEcas[1] += (hit->GetEREL())*fattenAC->Eval(dAC);
      if ((hit->GetPOS()==2)) fSumEcas[0] += (hit->GetEREL())*fattenAC->Eval(dAC);
      if ((hit->GetPOS()==1)) fSumEcas[2] += (hit->GetEREL())*fattenAC->Eval(dAC);
    }
    break;

  case 2: //C1D1
    for(Int_t k=0; k<hitsNo; k++){
      hit = (PamVMCDetectorHit*)HitColl->At(k);
      if ((hit->GetPOS()==1)) fSumEcard[3] += hit->GetEREL();
      if ((hit->GetPOS()==2)) fSumEcard[0] += hit->GetEREL(); 
      //fSumEcard[(hit->GetPOS())-1] += hit->GetEREL();
    }
    break;
  case 3: //C2D1
    for(Int_t k=0; k<hitsNo; k++){
      hit = (PamVMCDetectorHit*)HitColl->At(k);
      if ((hit->GetPOS()==1)) fSumEcard[2] += hit->GetEREL();
      if ((hit->GetPOS()==2)) fSumEcard[1] += hit->GetEREL(); 
      //fSumEcard[(hit->GetPOS())+1] += hit->GetEREL();
    }
    break;

  default:
    break;
  }


  

}