cal/src/PamVMCCaloDig.cxx

#include "PamVMCCaloDig.h"
#include <TTree.h>
#include "CalibCalPedEvent.h"

using namespace pamela;

ClassImp(PamVMCCaloDig)

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


void PamVMCCaloDig::LoadCalib(){

  cout<<"Loading Tracker Calibrations..."<<endl; 


  fhc =(TClonesArray*)fhitscolmap.GetValue("CAST");
  if (!fhc) cout <<"!!!WARNING Hit Collection of CAL not found by PamVMCCaloDig!!!!"<<endl;

  Int_t GivenCaloCalib = 0; //@@@ should be given as input par @@@

  Int_t time=3;
  Int_t type=101;

  fdberr = fsql->Query_GL_PARAM(time,type); 

  if(fdberr<0){
    cout<<"No such record in DB for CAL: time="<<time<<" type="<<type<<endl;
    ThrowCalFileWarning("CAL");

  } else {

    fquery.str("");
    fquery << fsql->GetPAR()->PATH.Data() << "/";
    fquery << fsql->GetPAR()->NAME.Data();

    ThrowCalFileUsage("CAL",fquery.str().c_str());


    FILE *f;
    f =  fopen(fquery.str().c_str(),"rb");

    if(f==NULL) ThrowCalFileWarning("CAL");
    else{
      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
      
      UInt_t idcalib, calibno;
      UInt_t utime = 0;
      TString calname;
      for (UInt_t s=0; s<4; s++){
        
        // clear calo calib variables for section s
        
        ClearCaloCalib(s);
        
        if ( GivenCaloCalib ){
          
          // a time has been given as input on the command line 
          // so retrieve the calibration that preceed that time
          
          fsql->Query_GL_CALO_CALIB(GivenCaloCalib,utime,s);
            
          calibno = fsql->GetCaloCalib()->EV_ROOT;
          idcalib = fsql->GetCaloCalib()->ID_ROOT_L0;
          
          // determine path and name and entry of the calibration file
          
          printf("\n");
          printf(" ** SECTION %i **\n",s);
          
          fsql->Query_GL_ROOT(idcalib);
          
          fquery.str("");
          fquery << fsql->GetROOT()->PATH.Data() << "/";
          fquery << fsql->GetROOT()->NAME.Data();
          
          calname = (TString)fquery.str().c_str();
          
          printf("\n Section %i : using  file %s calibration at entry %i: \n",s,calname.Data(),calibno);
          
        } else {

          fdberr = 0;
          fdberr = fsql->Query_GL_PARAM(1,104); 

          fquery.str("");
          fquery << fsql->GetPAR()->PATH.Data() << "/";
          fquery << fsql->GetPAR()->NAME.Data(); 
          
          printf("\n Section %i : using default calorimeter calibration: \n %s \n",s,fquery.str().c_str());

          calname = (TString)fquery.str().c_str();
          calibno = s;
          
        };
        
        // load calibration variables in memory
        
        CaloLoadCalib(s,calname,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
      //

      DigitizeCALOCALIB();
    }
  }        
}


void PamVMCCaloDig::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 PamVMCCaloDig::CaloLoadCalib(Int_t s,TString calname, UInt_t calibno){
 
  UInt_t e = 0;
  switch(s){
  case 0: e = 0;
    break;
  case 1: e = 2;
    break;
  case 2: e = 3;
    break;
  case 3: e = 1;
    break;
  default:
    break;
  }
  fcfile.open(calname.Data());
  ThrowCalFileUsage("CAL",calname.Data());
  if ( !fcfile ){
    ThrowCalFileWarning("CAL");
    return(-107);
  };
  fcfile.close();
  
  fcrfile = new TFile(calname.Data());
  if ( !fcrfile ){
    ThrowCalFileWarning("CAL");
    return(-108);
  }
  TTree *tr = (TTree*)fcrfile->Get("CalibCalPed");
  if ( !tr ){
    cout<<"!!!WARNING Tree CalibCalPed in file "<<calname.Data()
        <<" was NOT found!!!"<<endl;
    return(-109);
  }
  TBranch *calo = tr->GetBranch("CalibCalPed");
  CalibCalPedEvent *ce = 0;
  tr->SetBranchAddress("CalibCalPed", &ce);
  Long64_t ncalibs = calo->GetEntries();
  if ( !ncalibs ){
    cout<<"!!!WARNING No entries in calibration file!!!"<<endl; 
    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 ");
    fcrfile->Close();
    return(-111);
  };
  fcrfile->Close();
  delete ce;
  return(0);
}

void PamVMCCaloDig::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 

  Short_t CRC;

  for(Int_t sec=0; sec < 4; sec++){
    //
    // sec =    0 -> XE      1 -> XO        2-> YE         3 -> YO
    //
    fData.clear();

    fData.push_back(fSecCalo[sec]);      
    fData.push_back(fSecCaloLength[sec]);
    
    if ( sec==1 ){
      FillCalPedBack(&fData,sec);
      FillCalThrBack(&fData,sec);
      FillCalRmsBack(&fData,sec);
      FillCalBaseVarBack(&fData,sec);
    } else {
      FillCalPedNorm(&fData,sec);
      FillCalThrNorm(&fData,sec);
      FillCalRmsNorm(&fData,sec);
      FillCalBaseVarNorm(&fData,sec);
    }
    //CRC
    
    fData.push_back(0x0000);
    CRC = 0;

  
    USBuffer::const_iterator p = fData.begin();
    while( p!= fData.end() ){
      CRC=crc(CRC,*p);
      p++;
    }
    fData.push_back(CRC);

    //Writing calibrations
    DigitizePSCU(fData.size()*2,0x18);
    fraw->WritePSCU(&fDataPSCU);
    fraw->CopyUShortToBuff(&fData);

    //No padding for now...
  }
}

void PamVMCCaloDig::FillCalPedNorm(USBuffer *buff, Int_t sec){

  Int_t chksum = 0;
  for (Int_t plane=0; plane < 11; plane++){
      for (Int_t strip=0; strip < 96; strip++){   
        chksum += (Int_t)fcalped[sec][plane][strip];
        buff->push_back((Int_t)fcalped[sec][plane][strip]);  
        buff->push_back((Int_t)fcalgood[sec][plane][strip]);   
      }
     
  }//3-2114
  
  buff->push_back(((UShort_t)chksum));
  buff->push_back(0x0000);
  buff->push_back(((UShort_t)((Int_t)(chksum >> 16))));
}

void PamVMCCaloDig::FillCalPedBack(USBuffer *buff, Int_t sec){

  Int_t chksum = 0;
  for (Int_t plane=0; plane < 11; plane++){
      for (Int_t strip=95; strip >= 0; strip--){   
        chksum += (Int_t)fcalped[sec][plane][strip];
        buff->push_back((Int_t)fcalped[sec][plane][strip]);  
        buff->push_back((Int_t)fcalgood[sec][plane][strip]);    
      }     
  }//3-2114

  buff->push_back(((UShort_t)chksum));
  buff->push_back(0x0000);
  buff->push_back(((UShort_t)((Int_t)(chksum >> 16))));
}


void PamVMCCaloDig::FillCalThrNorm(USBuffer *buff, Int_t sec){

  Int_t chksum = 0;
  for (Int_t plane=0; plane < 11; plane++){
      for (Int_t strip=0; strip < 6; strip++){   
        chksum += (Int_t)fcalthr[sec][plane][strip];
          buff->push_back(0x0000);  
          buff->push_back((Int_t)fcalthr[sec][plane][strip]);   
      }
     
  }//2118-2249

  buff->push_back(0x0000);
  buff->push_back(((UShort_t)chksum));
  buff->push_back(0x0000);
  buff->push_back(((UShort_t)((Int_t)(chksum >> 16))));
}

void PamVMCCaloDig::FillCalThrBack(USBuffer *buff, Int_t sec){

  Int_t chksum = 0;
  for (Int_t plane=0; plane < 11; plane++){
      for (Int_t strip=5; strip >= 0; strip--){   
        chksum += (Int_t)fcalthr[sec][plane][strip];
          buff->push_back(0x0000);  
          buff->push_back((Int_t)fcalthr[sec][plane][strip]);    
      }     
  }//2118-2249

  buff->push_back(0x0000);
  buff->push_back(((UShort_t)chksum));
  buff->push_back(0x0000);
  buff->push_back(((UShort_t)((Int_t)(chksum >> 16))));
}

void PamVMCCaloDig::FillCalRmsNorm(USBuffer *buff, Int_t sec){

  for (Int_t plane=0; plane < 11; plane++){
      for (Int_t strip=0; strip < 96; strip++){   
          buff->push_back(0x0000);  
          buff->push_back((Int_t)fcalrms[sec][plane][strip]);   
      }
  }//2254-4365
}

void PamVMCCaloDig::FillCalRmsBack(USBuffer *buff, Int_t sec){

  for (Int_t plane=0; plane < 11; plane++){
      for (Int_t strip=95; strip >= 0; strip--){   
          buff->push_back(0x0000);  
          buff->push_back((Int_t)fcalrms[sec][plane][strip]);    
      }     
  }//2254-4365
}

void PamVMCCaloDig::FillCalBaseVarNorm(USBuffer *buff, Int_t sec){

  for (Int_t plane=0; plane < 11; plane++){
      for (Int_t strip=0; strip < 6; strip++){   
          buff->push_back(0x0000);
          buff->push_back(((Int_t)fcalbase[sec][plane][strip]));   
          buff->push_back(0x0000);  
          buff->push_back(((Int_t)fcalvar[sec][plane][strip]));  
      }
  }//4366-4629
}

void PamVMCCaloDig::FillCalBaseVarBack(USBuffer *buff, Int_t sec){

  for (Int_t plane=0; plane < 11; plane++){
      for (Int_t strip=5; strip >= 0; strip--){   
        buff->push_back(0x0000);  
        buff->push_back((Int_t)fcalbase[sec][plane][strip]);   
        buff->push_back(0x0000);  
        buff->push_back((Int_t)fcalvar[sec][plane][strip]);     
      }     
  }//4366-4629
}

void PamVMCCaloDig::Digitize(){

  cout<<"Digitizing CALO..."<<endl;

  Int_t ModCalo = 1; // 0 is RAW, 1 is COMPRESS, 2 is FULL    
                     //####@@@@ should be given as input par @@@@####


  // call different routines depending on the acq mode you want to simulate
  
  switch ( ModCalo ){
  case 0:
    this->DigitizeCaloRaw();
    break;
  case 1:
    this->DigitizeCaloCompress();
    break;
  case 2:
    this->DigitizeCaloFull();
    break;
  };
}

void PamVMCCaloDig::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;
  Int_t fCALOlength;
  Double_t pedenoise;
  Float_t rms = 0.;
  Float_t pedestal = 0.;
  //
  // clean the data structure
  //

  UShort_t DataCALO[4264];
  fData.clear();

  memset(DataCALO,0,sizeof(UShort_t)*4264);

  static const Float_t CALOGeV2MIPratio = 0.0001059994;

  // 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
    //
    DataCALO[fCALOlength] = fSecCalo[sec];   //1
    fCALOlength++;
    DataCALO[fCALOlength] = fSecCaloLength[sec]; //2
    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++){
      DataCALO[fCALOlength] = 0x0000;
      fCALOlength++;
    }; //3-8
    //
    //
    // here comes data
    //
    //
    // Section XO is read in the opposite direction respect to the others
    //
    if ( sec == 1 ){      
      tstrip = 96*11 + fCALOlength; //tstrip = 1064
    } 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++){

        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->GetCaloErel(sec,plane,strip);       
        //
        // convert it into ADC channels
        //      
        adcsig = int(ens*fCalomip[l][lpl][strip]/CALOGeV2MIPratio);
        //
        // 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 = frandom->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 ){
          DataCALO[tstrip] = adc;
          tstrip++;
        } else {
          DataCALO[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,DataCALO[i+fSecPointer]);
    };
    DataCALO[fCALOlength] = (UShort_t)CRC;
    fCALOlength++;
  };

  for (Int_t i = 0; i<fCALOlength; i++) fData.push_back(DataCALO[i]);
  //reverse(fData.begin(),fData.end());
  //for (Int_t i = 0; i<fCALOlength; i++) fData.push_back(DataCALO[fCALOlength-i-1]);
}


Float_t PamVMCCaloDig::GetCaloErel(Int_t sec, Int_t plane, Int_t strip){
  
  // determine plane and strip
  
  Int_t mplane = 0;

  switch (sec){
  case 0: //3
    mplane =  plane * 4 ; // it must be 0, 4, 8, ... (+1)  from plane = 0, 11  YO
    break;
  case 1: //2
    mplane =  plane * 4 + 2; // it must be 2, 6, 10, ... (+1)  from plane = 0, 11 YE
    break;
  case 2: //1
    mplane =  plane * 4 + 3 ; // it must be 3, 7, 11, ... (+1)  from plane = 0, 11 XO
    break;
  case 3: //0
    mplane =  plane * 4 + 1 ; // it must be 1, 5, 9, ... (+1)  from plane = 0, 11  XE
    break;
  }
  
 
  Float_t Erel = 0.;

  // search energy release in VMC output
  if (fhc){
    PamVMCDetectorHit * hit = (PamVMCDetectorHit*)fhc->At(mplane*96+strip);
    if (hit) Erel=hit->GetEREL();
    // if(Erel) cout<<"sec:"<<sec<<" plane:"<<plane<<" mplane:"<<mplane<<" AT:"<<mplane*96+strip<<" POS:"<<hit->GetPOS()<<endl;
  } //else cout<<"CALO HIT Collection pointer not found"<<endl;

  // if ((sec == 2) && (plane == 0)) Erel = 0.01; else Erel = 0.;
 
  return Erel;
}


void PamVMCCaloDig::DigitizeCaloCompress(){
  //
  // CompressMode implemented by C.Pizzolotto october 2009 
  //
  // some variables
  //
  Float_t ens = 0.;
  UInt_t adcsig = 0;
  UInt_t adcbase = 0;
  UInt_t adc[16];
  Float_t rms = 0.;
  Double_t pedenoise=0.;
  Float_t  pedestal = 0.;
  UInt_t   pedround[16];
  Float_t  thres[16];
  Float_t  goodflag[16];
  UInt_t min_adc = 0x7FFF;
  UInt_t min_adc_ch = 0;
  UInt_t l = 0;
  UInt_t lpl = 0;
  Int_t  plane = 0;
  Int_t pre;  
  Int_t  npre = 0; // number of pre between 0-5
  UInt_t strip = 0;
  UInt_t remainder;
  Float_t basesum=0.;
  Float_t basenof=0.;
  UInt_t  baseline=0;
  UInt_t fSecPointer = 0;
  UInt_t fNofTStripsPointer = 0;
  UInt_t NofTransmittedStrips = 0 ;
  Int_t fCALOlength;
  UShort_t DataCALO[9040]; //TOO LONG? 4264ma non e' vero che e' cosi' lungo...... CECI CECI CECI 
  fData.clear();
  static const Float_t CALOGeV2MIPratio = 0.0001059994;
  //
  // clean the data structure
  //
  memset(adc, 0,sizeof(adc));
  memset(pedround, 0,sizeof(pedround));
  memset(thres, 0,sizeof(thres));
  memset(goodflag, 0,sizeof(goodflag));
  //
  memset(DataCALO,0,sizeof(UShort_t)*9040);
  //
  // Header of the four sections
  //
  fSecCalo[0] = 0xEA00; // XE
  fSecCalo[1] = 0xF100; // XO
  fSecCalo[2] = 0xF600; // YE 
  fSecCalo[3] = 0xED00; // YO
  //
  // here comes raw data
  //
  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
    //
    DataCALO[fCALOlength] = fSecCalo[sec];
    fCALOlength++;
    DataCALO[fCALOlength] = 0; // Unknown: length must be calculated on fly
    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++){
      DataCALO[fCALOlength] = 0x0000;
      fCALOlength++;
    };
    //
    // second level trigger 
    //
    DataCALO[fCALOlength] = 0x0000;
    fCALOlength++;
    //
    // Nof strips transmetted: must be calculated on fly
    //
    fNofTStripsPointer = fCALOlength;
    DataCALO[fCALOlength] = 0x0000;
    fCALOlength++;
    NofTransmittedStrips=0;
    //
    // Identifier of calo data
    //
    DataCALO[fCALOlength] = 0xCA50;
    fCALOlength++;
    DataCALO[fCALOlength] = 0xCA50;
    fCALOlength++;
    DataCALO[fCALOlength] = 0xFFFF; // compresso
    fCALOlength++;
    //
    // Pedestal threashold table checksum
    //
    DataCALO[fCALOlength] = 0x0000; 
    fCALOlength++;
    //
    // Calorimeter event counter
    //
    DataCALO[fCALOlength] = Getevtcalo() ;
    fCALOlength++;
    //cout<<" evtcalo?"<<Getevtcalo()<<endl;
    //
    // Start here with data 
    //
    plane=-1;
    npre =-1;
    for (Int_t ipre=0; ipre< 66; ipre++){ // (11 planes*6 preampl)  
      //
      // which plane
      if ( (ipre % 6) == 0) {
        plane++;
      } 
      //
      pre=ipre;
      //
      // Adjust counter for plane X0
      if (sec==1) // conto invertito
       {
           remainder = pre % 6 ;
           pre = ((plane+1)*6) - remainder ;
       }  
      //
      if ( sec == 0 || sec == 3 ) lpl = plane * 2;
      if ( sec == 1 || sec == 2 ) lpl = (plane * 2) + 1;
      //
      // initialize min_adc
      min_adc =  0x7FFF;
      for (Int_t ch=0; ch <16; ch++){ // 16 channels each pre
        //
        // strip number
        //
        strip=((pre-(6*plane))*16)+ch;
        if(sec==1) strip = ((pre-(6*plane))*16)+(15-ch)-16;
        //
        // calculate npre a number between 0-5
        //
        if( sec==1) {
          if ( ((95-strip) % 16) == 0) {
            npre++;
            if(npre>5) npre=0;
          } 
        } else {
          if ( (strip % 16) == 0) {
            npre++;
            if(npre>5) npre=0;
          } 
        }
        //
        ens = this->GetCaloErel(sec,plane,strip);       
        //
        // convert it into ADC channels
        //      
        adcsig = int(ens*fCalomip[l][lpl][strip]/CALOGeV2MIPratio);
        //
        // sum baselines
        //
        adcbase = (UInt_t)fcalbase[sec][plane][npre]; 
        //
        // 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 = frandom->Gaus((Double_t)pedestal,(Double_t)rms); 
        //
        // Sum all contribution
        //
        adc[ch] = adcsig + adcbase + (Int_t)round(pedenoise);
        //
        // Signal saturation
        //
        if ( adc[ch] > 0x7FFF ) adc[ch] = 0x7FFF;
        //
        // save infos
        //
        pedround[ch] = (Int_t)round(pedestal) ;
        thres[ch]    = ( fcalthr[sec][plane][npre] );
        goodflag[ch] = ( fcalgood[sec][plane][strip] ); // if bad should be 255
        //
        // Find minimum adc in this preamp
        //
        if (  goodflag[ch]==0 && (adc[ch]-pedround[ch])<min_adc ) 
          {
            min_adc = ( adc[ch]-pedround[ch] ) ;
            min_adc_ch = ch ;
          }
      }; // close channel loop ch
      //
      // Find how many channels are below threshold 
      //
      Int_t nof_chs_below = 0;
      for (Int_t ch=0; ch <16; ch++){ // 16 channels each pre
        if ( goodflag[ch]==0 &&  ((adc[ch]-pedround[ch]) < (min_adc+thres[min_adc_ch])) ) 
          nof_chs_below++;
      };
      //
      // Transmit data: CASE nof_chs_below<9
      // 
      if(nof_chs_below<9)
        {
          if(sec==1) {
            DataCALO[fCALOlength] = 0x1000 + ipre ;
          } else {
            DataCALO[fCALOlength] = 0x1000 + pre ;
          }
          fCALOlength++;
          for (Int_t ch=0; ch <16; ch++)
            { 
              DataCALO[fCALOlength] = adc[ch]; 
              fCALOlength++;
              NofTransmittedStrips++;
            };
        }
      else 
        //
        // Transmit data: CASE nof_chs_below>=9 
        //
        { 
          if(sec==1) {
            DataCALO[fCALOlength] = 0x0800 + ipre ;
          } else {
            DataCALO[fCALOlength] = 0x0800 + pre; 
          }
          fCALOlength++;
          //
          // calculate baseline and save it
          //
          basenof=0;
          baseline=0;
          basesum=0;
          for (Int_t ch=0; ch <16; ch++){ 
            if(  goodflag[ch]==0 && ( (adc[ch]-pedround[ch])<(min_adc+thres[ch]) ) )
              {
                basesum = basesum +  (adc[ch]-pedround[ch]) ;
                basenof++;
              }
          }; 
          baseline = (Int_t)round( basesum / basenof );
          DataCALO[fCALOlength] = baseline;
          fCALOlength++;
          //
          // Transmit only channels > (min_adc+thres[ch])
          //
          for (Int_t ch=0; ch <16; ch++){ 
            if ( (adc[ch]-pedround[ch] )>(min_adc+thres[ch]) )
              { 
                DataCALO[fCALOlength] = ch;
                fCALOlength++;
                DataCALO[fCALOlength] = adc[ch];
                fCALOlength++; 
                NofTransmittedStrips++;
              }
          }; 
        } // close if nof_chs_below
    }; // close preampl loop
    //
    // Write the length 
    //
    DataCALO[fSecPointer+1] = (fCALOlength-fSecPointer+1)-2 ;
    // total length of the packet: -2: because the words with status and length are not included
    DataCALO[fNofTStripsPointer] = NofTransmittedStrips ;
    //
    // here we calculate and save the CRC
    //
    Short_t CRC = 0;
    for (UInt_t i=0; i<(fCALOlength-fSecPointer); i++){
      CRC=crc(CRC,DataCALO[i+fSecPointer]);
    };
    DataCALO[fCALOlength] = (UShort_t)CRC;
    fCALOlength++;
     //
  }; // close sec loop

  Incrementevtcalo();

  for (Int_t i = 0; i<fCALOlength; i++) fData.push_back(DataCALO[i]);

  return;


}

void PamVMCCaloDig::DigitizeCaloFull(){

  cout<<"!!!WARNING PamVMCCaloDig FULL MODE STILL NOT IMPLEMENTED!!!"<<endl;  
  
  this->DigitizeCaloRaw();
  return;

}