#include "Digitizer.h"

extern "C"{
  short crc(short, short);
};

void Digitizer::ClearCaloCalib(Int_t s){
  //
  fcstwerr[s] = 0;
  fcperror[s] = 0.;
  for ( Int_t d=0 ; d<11 ;d++  ){
    Int_t pre = -1;
    for ( Int_t j=0; j<96 ;j++){
      if ( j%16 == 0 ) pre++;
      fcalped[s][d][j] = 0.;
      fcstwerr[s] = 0.;
      fcperror[s] = 0.;
      fcalgood[s][d][j] = 0.;
      fcalthr[s][d][pre] = 0.;
      fcalrms[s][d][j] = 0.;
      fcalbase[s][d][pre] = 0.;
      fcalvar[s][d][pre] = 0.;
    };
  };
  return;
}

Int_t Digitizer::CaloLoadCalib(Int_t s,TString fcalname, UInt_t calibno){
  //
  //
  UInt_t e = 0;
  if ( s == 0 ) e = 0;
  if ( s == 1 ) e = 2;
  if ( s == 2 ) e = 3;
  if ( s == 3 ) e = 1;
  //
  ifstream myfile;
  myfile.open(fcalname.Data());
  if ( !myfile ){    
    return(-107);
  };
  myfile.close();
  //
  TFile *File = new TFile(fcalname.Data());
  if ( !File ) return(-108);
  TTree *tr = (TTree*)File->Get("CalibCalPed");
  if ( !tr ) return(-109);
  //
  TBranch *calo = tr->GetBranch("CalibCalPed");
  //
  pamela::CalibCalPedEvent *ce = 0;
  tr->SetBranchAddress("CalibCalPed", &ce);
  //
  Long64_t ncalibs = calo->GetEntries();
  //
  if ( !ncalibs ) return(-110); 
  //
  calo->GetEntry(calibno);
  //
  if (ce->cstwerr[s] != 0 && ce->cperror[s] == 0 ) {
    fcstwerr[s] = ce->cstwerr[s];
    fcperror[s] = ce->cperror[s];
    for ( Int_t d=0 ; d<11 ;d++  ){
      Int_t pre = -1;
      for ( Int_t j=0; j<96 ;j++){
	if ( j%16 == 0 ) pre++;
	fcalped[s][d][j] = ce->calped[e][d][j];
	fcalgood[s][d][j] = ce->calgood[e][d][j];
	fcalthr[s][d][pre] = ce->calthr[e][d][pre];
	fcalrms[s][d][j] = ce->calrms[e][d][j];
	fcalbase[s][d][pre] = ce->calbase[e][d][pre];
	fcalvar[s][d][pre] = ce->calvar[e][d][pre];
      };
    };
  } else {
    printf(" CALORIMETER - ERROR: problems finding a good calibration in this file! \n\n ");
    File->Close();
    return(-111);
  };
  File->Close();
  return(0);
}


void Digitizer::DigitizeCALOCALIB() {
  // 
  // Header of the four sections
  //
  fSecCalo[0] = 0xAA00; // XE
  fSecCalo[1] = 0xB100; // XO
  fSecCalo[2] = 0xB600; // YE
  fSecCalo[3] = 0xAD00; // YO
  //
  // length of the data is 0x1215 
  // 
  fSecCALOLength[0] = 0x1215; // XE
  fSecCALOLength[1] = 0x1215; // XO
  fSecCALOLength[2] = 0x1215; // YE
  fSecCALOLength[3] = 0x1215; // YO
  //
  Int_t chksum = 0;
  UInt_t tstrip = 0;
  UInt_t fSecPointer = 0;
  // 
  for (Int_t sec=0; sec < 4; sec++){
    //
    // sec =    0 -> XE      1 -> XO        2-> YE         3 -> YO
    //
    fCALOlength = 0;
    memset(fDataCALO,0,sizeof(UShort_t)*fCALObuffer);
    fSecPointer = fCALOlength;
    //
    // First of all we have section header and packet length
    //
    fDataCALO[fCALOlength] = fSecCalo[sec];
    fCALOlength++;
    fDataCALO[fCALOlength] = fSecCALOLength[sec];
    fCALOlength++;
    //
    // Section XO is read in the opposite direction respect to the others
    //
    chksum = 0;
    //
    for (Int_t plane=0; plane < 11; plane++){
      //
      if ( sec == 1 ) tstrip = fCALOlength + 96*2;
      //
      for (Int_t strip=0; strip < 96; strip++){   
	//
	chksum += (Int_t)fcalped[sec][plane][strip];
        //
        // save value
        // 
        if ( sec == 1 ){
          tstrip -= 2;
          fDataCALO[tstrip] = (Int_t)fcalped[sec][plane][strip];
          fDataCALO[tstrip+1] = (Int_t)fcalgood[sec][plane][strip];
        } else {
          fDataCALO[fCALOlength] = (Int_t)fcalped[sec][plane][strip];
          fDataCALO[fCALOlength+1] = (Int_t)fcalgood[sec][plane][strip];
        }; 
        fCALOlength +=2;
      };
      //
    };
    //
    fDataCALO[fCALOlength] = (UShort_t)chksum;
    fCALOlength++;
    fDataCALO[fCALOlength] = 0;
    fCALOlength++;
    fDataCALO[fCALOlength] = (UShort_t)((Int_t)(chksum >> 16));
    fCALOlength++;
    //
    // Section XO is read in the opposite direction respect to the others
    //
    chksum = 0;
    //
    for (Int_t plane=0; plane < 11; plane++){
      //
      if ( sec == 1 ) tstrip = fCALOlength+6*2;
      //
      for (Int_t strip=0; strip < 6; strip++){
        //
        chksum += (Int_t)fcalthr[sec][plane][strip];
        //
        // save value
        //
        if ( sec == 1 ){
          tstrip -= 2;
          fDataCALO[tstrip] = 0;
          fDataCALO[tstrip+1] = (Int_t)fcalthr[sec][plane][strip];
        } else {
          fDataCALO[fCALOlength] = 0;
          fDataCALO[fCALOlength+1] = (Int_t)fcalthr[sec][plane][strip];
        };
        fCALOlength +=2;
      };
      //
    };
    //
    fDataCALO[fCALOlength] = 0;
    fCALOlength++;
    fDataCALO[fCALOlength] = (UShort_t)chksum;
    fCALOlength++;
    fDataCALO[fCALOlength] = 0;
    fCALOlength++;
    fDataCALO[fCALOlength] = (UShort_t)((Int_t)(chksum >> 16));
    fCALOlength++;
    //
    // Section XO is read in the opposite direction respect to the others
    //
    for (Int_t plane=0; plane < 11; plane++){
      //
      if ( sec == 1 ) tstrip = fCALOlength+96*2;
      //
      for (Int_t strip=0; strip < 96; strip++){
        //
        // save value
        //
        if ( sec == 1 ){
          tstrip -= 2;
          fDataCALO[tstrip] = 0;
          fDataCALO[tstrip+1] = (Int_t)fcalrms[sec][plane][strip];
        } else {
          fDataCALO[fCALOlength] = 0;
          fDataCALO[fCALOlength+1] = (Int_t)fcalrms[sec][plane][strip];
        };
        fCALOlength += 2;
      };
      //
    };     
    //
    // Section XO is read in the opposite direction respect to the others
    //
    for (Int_t plane=0; plane < 11; plane++){
      //
      if ( sec == 1 ) tstrip = fCALOlength+6*4;
      //
      for (Int_t strip=0; strip < 6; strip++){
	//
        // save value
        //
        if ( sec == 1 ){
	  tstrip -= 4;
          fDataCALO[tstrip] = 0;
          fDataCALO[tstrip+1] = (Int_t)fcalbase[sec][plane][strip];
          fDataCALO[tstrip+2] = 0;
          fDataCALO[tstrip+3] = (Int_t)fcalvar[sec][plane][strip];
        } else { 
          fDataCALO[fCALOlength] = 0;
          fDataCALO[fCALOlength+1] = (Int_t)fcalbase[sec][plane][strip];
          fDataCALO[fCALOlength+2] = 0;
          fDataCALO[fCALOlength+3] = (Int_t)fcalvar[sec][plane][strip];
        };
        fCALOlength +=4;
      };
      //
    };      
    //
    //
    // here we calculate and save the CRC
    //
    fDataCALO[fCALOlength] = 0;
    fCALOlength++;
    Short_t CRC = 0;
    for (UInt_t i=0; i<(fCALOlength-fSecPointer); i++){
      CRC=crc(CRC,fDataCALO[i+fSecPointer]);
    };
    fDataCALO[fCALOlength] = (UShort_t)CRC;
    fCALOlength++;
    //
    UInt_t length=fCALOlength*2;
    DigitizePSCU(length,0x18,fDataPSCU);
    //
    // Add padding to 64 bits
    //
    AddPadding();
    //
    fOutputfile.write(reinterpret_cast<char*>(fDataPSCU),sizeof(UShort_t)*fPSCUbuffer);
    UShort_t temp[1000000];
    memset(temp,0,sizeof(UShort_t)*1000000);
    swab(fDataCALO,temp,sizeof(UShort_t)*fCALOlength);  // WE MUST SWAP THE BYTES!!!
    fOutputfile.write(reinterpret_cast<char*>(temp),sizeof(UShort_t)*fCALOlength);
    //
    // padding to 64 bytes
    //
    if ( fPadding ){
      fOutputfile.write(reinterpret_cast<char*>(fDataPadding),sizeof(UChar_t)*fPadding);
    };
    //
    //    
  };
  //
};

void Digitizer::CaloLoadCalib() {
  //
  fGivenCaloCalib = 0; //                                  ####@@@@ should be given as input par @@@@####
  //
  // first of all load the MIP to ADC conversion values
  //
  stringstream calfile;
  Int_t error = 0;
  GL_PARAM *glparam = new GL_PARAM();
  //
  // determine where I can find calorimeter ADC to MIP conversion file  
  //
  error = 0;
  error = glparam->Query_GL_PARAM(3,101,fDbc);
  //
  calfile.str("");
  calfile << glparam->PATH.Data() << "/";
  calfile << glparam->NAME.Data();
  //
  printf("\n Using Calorimeter ADC to MIP conversion file: \n %s \n",calfile.str().c_str());
  FILE *f;
  f = fopen(calfile.str().c_str(),"rb");
  //
  memset(fCalomip,0,4224*sizeof(fCalomip[0][0][0]));
  //
  for (Int_t m = 0; m < 2 ; m++ ){
    for (Int_t k = 0; k < 22; k++ ){
      for (Int_t l = 0; l < 96; l++ ){
	fread(&fCalomip[m][k][l],sizeof(fCalomip[m][k][l]),1,f);
      };
    };
  };
  fclose(f);
  //
  // determine which calibration has to be used and load it for each section
  //  
  GL_CALO_CALIB *glcalo = new GL_CALO_CALIB();
  GL_ROOT *glroot = new GL_ROOT();  
  TString fcalname;
  UInt_t idcalib;
  UInt_t calibno;
  UInt_t utime = 0;
  //
  for (UInt_t s=0; s<4; s++){
    //
    // clear calo calib variables for section s
    //
    ClearCaloCalib(s);
    //
    if ( fGivenCaloCalib ){
      //
      // a time has been given as input on the command line so retrieve the calibration that preceed that time
      //
      glcalo->Query_GL_CALO_CALIB(fGivenCaloCalib,utime,s,fDbc);
      //  
      calibno = glcalo->EV_ROOT;
      idcalib = glcalo->ID_ROOT_L0;
      //
      // determine path and name and entry of the calibration file
      //
      printf("\n");
      printf(" ** SECTION %i **\n",s);
      //
      glroot->Query_GL_ROOT(idcalib,fDbc);
      //
      stringstream name;
      name.str("");
      name << glroot->PATH.Data() << "/";
      name << glroot->NAME.Data();
      //
      fcalname = (TString)name.str().c_str();
      //
      printf("\n Section %i : using  file %s calibration at entry %i: \n",s,fcalname.Data(),calibno);
      //
    } else {
      error = 0;
      error = glparam->Query_GL_PARAM(1,104,fDbc);
      //
      calfile.str("");
      calfile << glparam->PATH.Data() << "/";
      calfile << glparam->NAME.Data();
      //
      printf("\n Section %i : using default calorimeter calibration: \n %s \n",s,calfile.str().c_str());
      //
      fcalname = (TString)calfile.str().c_str();
      calibno = s;
      //
    };
    //
    // load calibration variables in memory
    //
    CaloLoadCalib(s,fcalname,calibno);
    //
  };
  //
  // at this point we have in memory the calorimeter calibration and we can save it to disk in the correct format and use it to digitize the data
  //
  delete glparam;
  delete glcalo;
  delete glroot;
};

void Digitizer::DigitizeCALO() {
  //
  fModCalo = 0; // 0 is RAW, 1 is COMPRESS, 2 is FULL     ####@@@@ should be given as input par @@@@####
  //
  //
  // 
  fCALOlength = 0;  // reset total dimension of calo data
  //
  // gpamela variables to be used
  //
//   fhBookTree->SetBranchStatus("Nthcali",1);//modified by E.Vannuccini 03/08
//   fhBookTree->SetBranchStatus("Icaplane",1);
//   fhBookTree->SetBranchStatus("Icastrip",1);
//   fhBookTree->SetBranchStatus("Icamod",1);
//   fhBookTree->SetBranchStatus("Enestrip",1);
  //
  // call different routines depending on the acq mode you want to simulate
  //
  switch ( fModCalo ){
  case 0:
    this->DigitizeCALORAW();
    break;
  case 1:
    this->DigitizeCALOCOMPRESS();
    break;
  case 2:
    this->DigitizeCALOFULL();
    break;
  };
  //
};

Float_t Digitizer::GetCALOen(Int_t sec, Int_t plane, Int_t strip){
  //
  // determine plane and strip
  //
  Int_t mplane = 0;
  //
  // wrong!
  // 
  //  if ( sec == 0 || sec == 3 ) mplane = (plane * 4) + sec + 1;  
  //  if ( sec == 1 ) mplane = (plane * 4) + 2 + 1;  
  //  if ( sec == 2 ) mplane = (plane * 4) + 1 + 1;  
  //
  if ( sec == 0 ) mplane = plane * 4 + 1; // it must be 0, 4, 8, ... (+1)  from plane = 0, 11
  if ( sec == 1 ) mplane = plane * 4 + 2 + 1; // it must be 2, 6, 10, ... (+1)  from plane = 0, 11
  if ( sec == 2 ) mplane = plane * 4 + 3 + 1; // it must be 3, 7, 11, ... (+1)  from plane = 0, 11
  if ( sec == 3 ) mplane = plane * 4 + 1 + 1; // it must be 1, 5, 9, ... (+1)  from plane = 0, 11
  //
  Int_t mstrip = strip + 1;
  //
  // search energy release in gpamela output
  //
  for (Int_t i=0; i<Nthcali;i++){
    if ( Icaplane[i] == mplane && Icastrip[i] == mstrip ){
      return (Enestrip[i]);
    };
  };
  //
  // if not found it means no energy release so return 0.
  //
  return(0.);
};

void Digitizer::DigitizeCALORAW() {
  //
  // some variables
  //
  Float_t ens = 0.;
  UInt_t adcsig = 0;
  UInt_t adcbase = 0;
  UInt_t adc = 0;
  Int_t pre = 0;
  UInt_t l = 0;
  UInt_t lpl = 0;
  UInt_t tstrip = 0;
  UInt_t fSecPointer = 0;
  Double_t pedenoise;
  Float_t rms = 0.;
  Float_t pedestal = 0.;
  //
  // clean the data structure
  //
  memset(fDataCALO,0,sizeof(UShort_t)*fCALObuffer);
  //
  // Header of the four sections
  //
  fSecCalo[0] = 0xEA08; // XE
  fSecCalo[1] = 0xF108; // XO
  fSecCalo[2] = 0xF608; // YE 
  fSecCalo[3] = 0xED08; // YO
  //
  // length of the data is 0x0428 in RAW mode
  //
  fSecCALOLength[0] = 0x0428; // XE
  fSecCALOLength[1] = 0x0428; // XO
  fSecCALOLength[2] = 0x0428; // YE
  fSecCALOLength[3] = 0x0428; // YO
  //
  // let's start
  //
  fCALOlength = 0;
  //
  for (Int_t sec=0; sec < 4; sec++){
    //
    // sec =    0 -> XE      1 -> XO        2-> YE         3 -> YO
    //
    l = 0;                 // XE and XO are Y planes 
    if ( sec < 2 ) l = 1;  // while YE and  YO are X planes
    //
    fSecPointer = fCALOlength;
    //
    // First of all we have section header and packet length
    //
    fDataCALO[fCALOlength] = fSecCalo[sec];
    fCALOlength++;
    fDataCALO[fCALOlength] = fSecCALOLength[sec];
    fCALOlength++;
    //
    // selftrigger coincidences - in the future we should add here some code to simulate timing response of pre-amplifiers
    //
    for (Int_t autoplane=0; autoplane < 7; autoplane++){
      fDataCALO[fCALOlength] = 0x0000;
      fCALOlength++;
    };
    //
    //
    // here comes data
    //
    //
    // Section XO is read in the opposite direction respect to the others
    //
    if ( sec == 1 ){	  
      tstrip = 96*11 + fCALOlength;
    } else {
      tstrip = 0;
    };
    //
    pre = -1;
    //
    for (Int_t strip=0; strip < 96; strip++){   
      //
      // which is the pre for this strip?
      //
      if (strip%16 == 0) {
	pre++;
      };
      //
      if ( sec == 1 ) tstrip -= 11;
      //
      for (Int_t plane=0; plane < 11; plane++){
	//
	// here is wrong!!!!
	//
	//
	//	if ( plane%2 == 0 && sec%2 != 0){
	//	  lpl = plane*2;
	//	} else {
	//	  lpl = (plane*2) + 1;
	//	};
	//
	if ( sec == 0 || sec == 3 ) lpl = plane * 2;
	if ( sec == 1 || sec == 2 ) lpl = (plane * 2) + 1;
	//
	// get the energy in GeV from the simulation for that strip
	//
	ens = this->GetCALOen(sec,plane,strip);	
	//
	// convert it into ADC channels
	//	
	adcsig = int(ens*fCalomip[l][lpl][strip]/fCALOGeV2MIPratio);
	//
	// sum baselines
	//
	adcbase = (UInt_t)fcalbase[sec][plane][pre]; 
	//
	// add noise and pedestals
	//	
	pedestal = fcalped[sec][plane][strip];
	rms = fcalrms[sec][plane][strip]/4.;
	//
	// Add random gaussian noise of RMS rms and Centered in the pedestal
	// 
	pedenoise = gRandom->Gaus((Double_t)pedestal,(Double_t)rms); 
	//
	// Sum all contribution
	//
	adc = adcsig + adcbase + (Int_t)round(pedenoise);
	//
	// Signal saturation
	//
	if ( adc > 0x7FFF ) adc = 0x7FFF;
	//
	// save value
	//
	if ( sec == 1 ){
	  fDataCALO[tstrip] = adc;
	  tstrip++;
	} else {
	  fDataCALO[fCALOlength] = adc;
	};
	fCALOlength++;
	//
      };
      //
      if ( sec == 1 ) tstrip -= 11;
      //
    };
    //
    // here we calculate and save the CRC
    //
    Short_t CRC = 0;
    for (UInt_t i=0; i<(fCALOlength-fSecPointer); i++){
      CRC=crc(CRC,fDataCALO[i+fSecPointer]);
    };
    fDataCALO[fCALOlength] = (UShort_t)CRC;
    fCALOlength++;
    //
  };
  //
  //   for (Int_t i=0; i<fCALOlength; i++){
  //     printf(" WORD %i       DIGIT %0x   \n",i,fDataCALO[i]);
  //   };
  //
};

void Digitizer::DigitizeCALOCOMPRESS() {
  //
  printf(" COMPRESS MODE STILL NOT IMPLEMENTED! \n");  
  //
  this->DigitizeCALORAW();
  return;
  //
  //
  //
  fSecCalo[0] = 0xEA00;
  fSecCalo[1] = 0xF100;
  fSecCalo[2] = 0xF600;
  fSecCalo[3] = 0xED00;
  //
  // length of the data in DSP mode must be calculated on fly during digitization
  //
  memset(fSecCALOLength,0x0,4*sizeof(UShort_t));
  //
  // here comes raw data
  //
  Int_t en = 0;
  //
  for (Int_t sec=0; sec < 4; sec++){
    fDataCALO[en] = fSecCalo[sec];
    en++;
    fDataCALO[en] = fSecCALOLength[sec];
    en++;
    for (Int_t plane=0; plane < 11; plane++){
      for (Int_t strip=0; strip < 11; strip++){
	fDataCALO[en] = 0x0;
	en++;
      };
    };
  };
  //
};

void Digitizer::DigitizeCALOFULL() {
  //
  printf(" FULL MODE STILL NOT IMPLEMENTED! \n");  
  //
  this->DigitizeCALORAW();
  return;
  //
  fSecCalo[0] = 0xEA00;
  fSecCalo[1] = 0xF100;
  fSecCalo[2] = 0xF600;
  fSecCalo[3] = 0xED00;
  //
  // length of the data in DSP mode must be calculated on fly during digitization
  //
  memset(fSecCALOLength,0x0,4*sizeof(UShort_t));
  //
  // here comes raw data
  //
  Int_t  en = 0;
  //
  for (Int_t sec=0; sec < 4; sec++){
    fDataCALO[en] = fSecCalo[sec];
    en++;
    fDataCALO[en] = fSecCALOLength[sec];
    en++;
    for (Int_t plane=0; plane < 11; plane++){
      for (Int_t strip=0; strip < 11; strip++){
	fDataCALO[en] = 0x0;
	en++;
      };
    };
  };
  //
};
