--- PamelaDigitizer/Digitizer.cxx	2007/11/28 18:54:31	1.5
+++ PamelaDigitizer/Digitizer.cxx	2009/10/16 09:15:50	1.14
@@ -1,2722 +1,423 @@
-//               ------ PAMELA Digitizer ------
-// 
-// Date, release and how-to: see file Pamelagp2Digits.cxx 
-// 
-// NB: Check length physics packet [packet type (0x10 = physics data)]
-//
-#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 "TMath.h"
-#include "TRandom.h"
-#include "TSQLServer.h"
-#include "TSystem.h"
-//
-#include "Digitizer.h"
-#include "CRC.h"
-//
-#include <PamelaRun.h>
-#include <physics/calorimeter/CalorimeterEvent.h>
-#include <CalibCalPedEvent.h>
-#include "GLTables.h"
-//
-extern "C"{
-  short crc(short, short);
-};
-//
-
-Digitizer::Digitizer(TTree* tree, char* &file_raw){
-  fhBookTree = tree;
-  fFilename =  file_raw;
-  fCounter = 0;
-  fOBT = 0;
-
-  //
-  // DB connections
-  //
-  TString host = "mysql://localhost/pamelaprod";
-  TString user = "anonymous";
-  TString psw = "";
-  //
-  const char *pamdbhost=gSystem->Getenv("PAM_DBHOST");
-  const char *pamdbuser=gSystem->Getenv("PAM_DBUSER");
-  const char *pamdbpsw=gSystem->Getenv("PAM_DBPSW");
-  if ( !pamdbhost ) pamdbhost = "";
-  if ( !pamdbuser ) pamdbuser = "";
-  if ( !pamdbpsw ) pamdbpsw = "";
-  if ( strcmp(pamdbhost,"") ) host = pamdbhost;
-  if ( strcmp(pamdbuser,"") ) user = pamdbuser;
-  if ( strcmp(pamdbpsw,"") ) psw = pamdbpsw;
-  fDbc = TSQLServer::Connect(host.Data(),user.Data(),psw.Data());
-  //
-  GL_TABLES *glt = new GL_TABLES(host,user,psw);
-  if ( glt->IsConnected(fDbc) ) printf("\n DB INFORMATION:\n SQL: %s Version: %s Host %s Port %i \n\n",fDbc->GetDBMS(),fDbc->ServerInfo(),fDbc->GetHost(),fDbc->GetPort());
-  //
-  // Use UTC in the DB and make timeout bigger
-  //
-  stringstream myquery;
-  myquery.str("");
-  myquery << "SET time_zone='+0:00'";
-  fDbc->Query(myquery.str().c_str());
-  myquery.str("");
-  myquery << "SET wait_timeout=173000;";
-  fDbc->Query(myquery.str().c_str());
-  //
-  
-  std:: cout << "preparing tree" << endl;
-
-  // prepare tree
-  fhBookTree->SetBranchAddress("Irun",&Irun);
-  fhBookTree->SetBranchAddress("Ievnt",&Ievnt);
-  fhBookTree->SetBranchAddress("Ipa",&Ipa);
-  fhBookTree->SetBranchAddress("X0",&X0);
-  fhBookTree->SetBranchAddress("Y0",&Y0);
-  fhBookTree->SetBranchAddress("Z0",&Z0);
-  fhBookTree->SetBranchAddress("Theta",&Theta);
-  fhBookTree->SetBranchAddress("Phi",&Phi);
-  fhBookTree->SetBranchAddress("P0",&P0);
-  fhBookTree->SetBranchAddress("Nthtof",&Nthtof);
-  fhBookTree->SetBranchAddress("Ipltof",Ipltof);
-  fhBookTree->SetBranchAddress("Ipaddle",Ipaddle);
-  fhBookTree->SetBranchAddress("Ipartof",Ipartof);
-  fhBookTree->SetBranchAddress("Xintof",Xintof);
-  fhBookTree->SetBranchAddress("Yintof",Yintof);
-  fhBookTree->SetBranchAddress("Zintof",Zintof);
-  fhBookTree->SetBranchAddress("Xouttof",Xouttof);
-  fhBookTree->SetBranchAddress("Youttof",Youttof);
-  fhBookTree->SetBranchAddress("Zouttof",Zouttof);
-  fhBookTree->SetBranchAddress("Ereltof",Ereltof);
-  fhBookTree->SetBranchAddress("Timetof",Timetof);
-  fhBookTree->SetBranchAddress("Pathtof",Pathtof);
-  fhBookTree->SetBranchAddress("P0tof",P0tof);
-  fhBookTree->SetBranchAddress("Nthcat",&Nthcat);
-  fhBookTree->SetBranchAddress("Iparcat",Iparcat);
-  fhBookTree->SetBranchAddress("Icat",Icat);
-  fhBookTree->SetBranchAddress("Xincat",Xincat);
-  fhBookTree->SetBranchAddress("Yincat",Yincat);
-  fhBookTree->SetBranchAddress("Zincat",Zincat);
-  fhBookTree->SetBranchAddress("Xoutcat",Xoutcat);
-  fhBookTree->SetBranchAddress("Youtcat",Youtcat);
-  fhBookTree->SetBranchAddress("Zoutcat",Zoutcat);
-  fhBookTree->SetBranchAddress("Erelcat",Erelcat);
-  fhBookTree->SetBranchAddress("Timecat",Timecat);
-  fhBookTree->SetBranchAddress("Pathcat",Pathcat);
-  fhBookTree->SetBranchAddress("P0cat",P0cat);
-  fhBookTree->SetBranchAddress("Nthcas",&Nthcas);
-  fhBookTree->SetBranchAddress("Iparcas",Iparcas);
-  fhBookTree->SetBranchAddress("Icas",Icas);
-  fhBookTree->SetBranchAddress("Xincas",Xincas);
-  fhBookTree->SetBranchAddress("Yincas",Yincas);
-  fhBookTree->SetBranchAddress("Zincas",Zincas);
-  fhBookTree->SetBranchAddress("Xoutcas",Xoutcas);
-  fhBookTree->SetBranchAddress("Youtcas",Youtcas);
-  fhBookTree->SetBranchAddress("Zoutcas",Zoutcas);
-  fhBookTree->SetBranchAddress("Erelcas",Erelcas);
-  fhBookTree->SetBranchAddress("Timecas",Timecas);
-  fhBookTree->SetBranchAddress("Pathcas",Pathcas);
-  fhBookTree->SetBranchAddress("P0cas",P0cas);
-  fhBookTree->SetBranchAddress("Nthspe",&Nthspe);
-  fhBookTree->SetBranchAddress("Iparspe",Iparspe);
-  fhBookTree->SetBranchAddress("Itrpb",Itrpb);
-  fhBookTree->SetBranchAddress("Itrsl",Itrsl);
-  fhBookTree->SetBranchAddress("Itspa",Itspa);
-  fhBookTree->SetBranchAddress("Xinspe",Xinspe);
-  fhBookTree->SetBranchAddress("Yinspe",Yinspe);
-  fhBookTree->SetBranchAddress("Zinspe",Zinspe);
-  fhBookTree->SetBranchAddress("Xoutspe",Xoutspe);
-  fhBookTree->SetBranchAddress("Youtspe",Youtspe);
-  fhBookTree->SetBranchAddress("Zoutspe",Zoutspe);
-  fhBookTree->SetBranchAddress("Xavspe",Xavspe);
-  fhBookTree->SetBranchAddress("Yavspe",Yavspe);
-  fhBookTree->SetBranchAddress("Zavspe",Zavspe);
-  fhBookTree->SetBranchAddress("Erelspe",Erelspe);
-  fhBookTree->SetBranchAddress("Pathspe",Pathspe);
-  fhBookTree->SetBranchAddress("P0spe",P0spe);
-  fhBookTree->SetBranchAddress("Nxmult",Nxmult);
-  fhBookTree->SetBranchAddress("Nymult",Nymult);
-  fhBookTree->SetBranchAddress("Nstrpx",&Nstrpx);
-  fhBookTree->SetBranchAddress("Npstripx",Npstripx);
-  fhBookTree->SetBranchAddress("Ntstripx",Ntstripx);
-  fhBookTree->SetBranchAddress("Istripx",Istripx);
-  fhBookTree->SetBranchAddress("Qstripx",Qstripx);
-  fhBookTree->SetBranchAddress("Xstripx",Xstripx);
-  fhBookTree->SetBranchAddress("Nstrpy",&Nstrpy);
-  fhBookTree->SetBranchAddress("Npstripy",Npstripy);
-  fhBookTree->SetBranchAddress("Ntstripy",Ntstripy);
-  fhBookTree->SetBranchAddress("Istripy",Istripy);
-  fhBookTree->SetBranchAddress("Qstripy",Qstripy);
-  fhBookTree->SetBranchAddress("Ystripy",Ystripy);
-  fhBookTree->SetBranchAddress("Nthcali",&Nthcali);
-  fhBookTree->SetBranchAddress("Icaplane",Icaplane);
-  fhBookTree->SetBranchAddress("Icastrip",Icastrip);
-  fhBookTree->SetBranchAddress("Icamod",Icamod);
-  fhBookTree->SetBranchAddress("Enestrip",Enestrip);
-  fhBookTree->SetBranchAddress("Nthcal",&Nthcal);
-  fhBookTree->SetBranchAddress("Icapl",Icapl);
-  fhBookTree->SetBranchAddress("Icasi",Icasi);
-  fhBookTree->SetBranchAddress("Icast",Icast);
-  fhBookTree->SetBranchAddress("Xincal",Xincal);
-  fhBookTree->SetBranchAddress("Yincal",Yincal);
-  fhBookTree->SetBranchAddress("Zincal",Zincal);
-  fhBookTree->SetBranchAddress("Erelcal",Erelcal);
-  fhBookTree->SetBranchAddress("Nthnd",&Nthnd);
-  fhBookTree->SetBranchAddress("Itubend",Itubend);
-  fhBookTree->SetBranchAddress("Iparnd",Iparnd);
-  fhBookTree->SetBranchAddress("Xinnd",Xinnd);
-  fhBookTree->SetBranchAddress("Yinnd",Yinnd);
-  fhBookTree->SetBranchAddress("Zinnd",Zinnd);
-  fhBookTree->SetBranchAddress("Xoutnd",Xoutnd);
-  fhBookTree->SetBranchAddress("Youtnd",Youtnd);
-  fhBookTree->SetBranchAddress("Zoutnd",Zoutnd);
-  fhBookTree->SetBranchAddress("Erelnd",Erelnd);
-  fhBookTree->SetBranchAddress("Timend",Timend);
-  fhBookTree->SetBranchAddress("Pathnd",Pathnd);
-  fhBookTree->SetBranchAddress("P0nd",P0nd);
-  fhBookTree->SetBranchAddress("Nthcard",&Nthcard);
-  fhBookTree->SetBranchAddress("Iparcard",Iparcard);
-  fhBookTree->SetBranchAddress("Icard",Icard);
-  fhBookTree->SetBranchAddress("Xincard",Xincard);
-  fhBookTree->SetBranchAddress("Yincard",Yincard);
-  fhBookTree->SetBranchAddress("Zincard",Zincard);
-  fhBookTree->SetBranchAddress("Xoutcard",Xoutcard);
-  fhBookTree->SetBranchAddress("Youtcard",Youtcard);
-  fhBookTree->SetBranchAddress("Zoutcard",Zoutcard);
-  fhBookTree->SetBranchAddress("Erelcard",Erelcard);
-  fhBookTree->SetBranchAddress("Timecard",Timecard);
-  fhBookTree->SetBranchAddress("Pathcard",Pathcard);
-  fhBookTree->SetBranchAddress("P0card",P0card);
-
-  fhBookTree->SetBranchStatus("*",0); 
-
-};
-
-
-
-void Digitizer::Close(){
-
-  delete fhBookTree;
-
-};
-
-
-
-
-void Digitizer::Loop() {
-  //
-  // opens the raw output file and loops over the events
-  //
-  fOutputfile.open(fFilename, ios::out | ios::binary);
-  //fOutputfile.open(Form("Output%s",fFilename), ios::out | ios::binary);
-  //
-  // Load in memory and save at the beginning of file the calorimeter calibration
-  //
-  CaloLoadCalib();
-  DigitizeCALOCALIB();
-
-  //  load, digitize and write tracker calibration
-  LoadTrackCalib();
-  
-  DigitizeTrackCalib(1);
-  UInt_t length=fTracklength*2;
-  DigitizePSCU(length,0x12);
-  AddPadding();
-  WriteTrackCalib(); 
-  
-  DigitizeTrackCalib(2);
-  length=fTracklength*2;
-  DigitizePSCU(length,0x13);
-  AddPadding();
-  WriteTrackCalib();
-  
-  LoadMipCor();  // some initialization of parameters -not used now-
-  //  end loading, digitizing and writing tracker calibration
-
-  //
-  // loops over the events
-  //
-  
-  Int_t nentries = fhBookTree->GetEntriesFast();
-  Long64_t nbytes = 0;
-  for (Int_t i=0; i<nentries;i++) {
-    //
-    nbytes += fhBookTree->GetEntry(i);
-    // read detectors sequentially:
-    // http://www.ts.infn.it/fileadmin/documents/physics/experiments/wizard/cpu/gen_arch/RM_Acquisition.pdf
-    // on pamelatov:
-    // /cvs/yoda/techmodel/physics/NeutronDetectorReader.cpp
-    DigitizeTRIGGER();
-    DigitizeTOF();
-    DigitizeAC();
-    DigitizeCALO();
-    DigitizeTrack();
-    DigitizeS4();
-    DigitizeND();
-      //
-    // Add padding to 64 bits
-    //
-    AddPadding();
-//
-    // Create CPU header, we need packet type (0x10 = physics data) and packet length.
-    //
-    UInt_t length=2*(fCALOlength+fACbuffer+fTracklength+fNDbuffer+fS4buffer)+fPadding+fTOFbuffer+fTRIGGERbuffer;
-    //UInt_t length=2*(fCALOlength+fACbuffer+fTracklength+fNDbuffer)+fPadding+fTOFbuffer+fTRIGGERbuffer;
-    DigitizePSCU(length,0x10);
-    if ( !i%100 ) std::cout << "writing event " << i << endl;
-    WriteData();
-  };
-
-  fOutputfile.close();
-  std::cout << "files closed" << endl << flush;
-
-};
-
-void Digitizer::AddPadding(){
-  //
-  Float_t pd0 = (fLen+16)/64.;
-  Float_t pd1 = pd0 - (Float_t)int(pd0);
-  Float_t padfrac = 64. - pd1 * 64.;
-  //
-  UInt_t padbytes = (UInt_t)padfrac;
-  if ( padbytes > 0 && padbytes < 64 ){
-    //
-    // here the padding length
-    //
-    fPadding = padbytes+64;
-    //
-    // random padding bytes
-    //
-    for (Int_t ur=0; ur<32; ur++){
-      fDataPadding[ur] = (UShort_t)rand();
-    };
-  };
-};
-
-
-void Digitizer::DigitizePSCU(UInt_t length, UChar_t type) {
-  //
-  UChar_t buff[16];
-  //
-  // CPU signature
-  //  
-  buff[0] = 0xFA;
-  buff[1] = 0xFE;
-  buff[2] = 0xDE;
-  //
-  // packet type (twice)
-  //
-  buff[3] = type;
-  buff[4] = type;
-  //
-  // counter
-  //
-  fCounter++;
-  while ( fCounter > 16777215 ){
-    fCounter -= 16777215;
-  };
-  //
-  buff[5] = (UChar_t)(fCounter >> 16);
-  buff[6] = (UChar_t)(fCounter >> 8);
-  buff[7] = (UChar_t)fCounter;
-  //
-  // on board time
-  //
-  ULong64_t obt = fOBT + 30LL;
-  //
-  while ( obt > 4294967295LL ){
-    obt -= 4294967295LL;
-  };
-  fOBT = (UInt_t)obt;
-  //
-  buff[8] = (UChar_t)(fOBT >> 24);
-  buff[9] = (UChar_t)(fOBT >> 16);
-  buff[10] = (UChar_t)(fOBT >> 8);
-  buff[11] = (UChar_t)fOBT;
-  //
-  // Packet length
-  //
-  fLen = length;
-  //
-  buff[12] = (UChar_t)(fLen >> 16);
-  buff[13] = (UChar_t)(fLen >> 8);
-  buff[14] = (UChar_t)fLen;
-  //
-  // CPU header CRC
-  //
-  buff[15] = (BYTE)CM_Compute_CRC16((UINT16)0, (BYTE*)&buff, (UINT32)15);
-  //
-  memcpy(fDataPSCU,buff,16*sizeof(UChar_t));
-  //
-};
-
-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);
-    //
-    // 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);
-  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++;
-      };
-    };
-  };
-  //
-};
-
-void Digitizer::DigitizeTRIGGER() {
-  //fDataTrigger: 153 bytes
-  for (Int_t j=0; j < 153; j++)
-    fDataTrigger[j]=0x00;
-};
-
-Int_t Digitizer::DigitizeTOF() {
-  //fDataTof: 12 x 23 bytes (=276 bytes)
-  UChar_t *pTof=fDataTof;
-  Bool_t DEBUG=false;
-
-  // --- activate branches:
-  fhBookTree->SetBranchStatus("Nthtof",1);
-  fhBookTree->SetBranchStatus("Ipltof",1);
-  fhBookTree->SetBranchStatus("Ipaddle",1);
-  fhBookTree->SetBranchStatus("Xintof",1);
-  fhBookTree->SetBranchStatus("Yintof",1);
-  fhBookTree->SetBranchStatus("Xouttof",1);
-  fhBookTree->SetBranchStatus("Youttof",1);
-  fhBookTree->SetBranchStatus("Ereltof",1);
-  fhBookTree->SetBranchStatus("Timetof",1);
-  // not yet used: Zintof, Xouttof, Youttof, Zouttof
-
-  // ------ evaluate energy in each pmt: ------
-  // strip geometry (lenght/width)
-  Float_t dimel[6] = {33.0, 40.8 ,18.0, 15.0, 15.0, 18.0};
-  //Float_t dimes[6] = {5.1, 5.5, 7.5, 9.0, 6.0, 5.0};
-  
-  //  S11 8 paddles  33.0 x 5.1 cm
-  //  S12 6 paddles  40.8 x 5.5 cm
-  //  S21 2 paddles  18.0 x 7.5 cm
-  //  S22 2 paddles  15.0 x 9.0 cm
-  //  S31 3 paddles  15.0 x 6.0 cm
-  //  S32 3 paddles  18.0 x 5.0 cm
-
-  Float_t FGeo[2]={0., 0.}; /* geometrical factor */
- 
-  const Float_t Pho_keV = 10.;     // photons per keV in scintillator
-  const Float_t echarge = 1.6e-19; // electron charge
-  Float_t Npho=0.;
-  Float_t QevePmt_pC[48];
-  Float_t QhitPad_pC[2]={0., 0.};
-  Float_t QhitPmt_pC[2]={0., 0.};
-  Float_t pmGain = 3.5e6;  /* PMT Gain: the same for all PMTs */
-  Float_t effi=0.21;       /* Efficienza di fotocatodo */
-
-  //   Float_t ADC_pC0=-58.1;      //  ADC/pC conversion coefficient 0
-  //   Float_t ADC_pC1=1.728;      //  ADC/pC conversion coefficient 1 
-  //   Float_t ADC_pC2=-4.063e-05; //  ADC/pC conversion coefficient 2
-  //   Float_t ADC_pC3=-5.763e-08; //  ADC/pC conversion coefficient 3
-  
-  // pC < 800 
-  Float_t ADC_pC0A =   -4.437616e+01   ;
-  Float_t ADC_pC1A =   1.573329e+00    ;
-  Float_t ADC_pC2A =   2.780518e-04  ;
-  Float_t ADC_pC3A =  -2.302160e-07 ;
-  
-  // pC   > 800:
-  Float_t ADC_pC0B =    -2.245756e+02  ;
-  Float_t ADC_pC1B =     2.184156e+00  ;
-  Float_t ADC_pC2B =     -4.171825e-04  ;
-  Float_t ADC_pC3B =     3.789715e-08  ;
-
-  Float_t pCthres=40.;     // threshold in charge
-  Int_t ADClast=4095;      // no signal --> ADC ch=4095
-  Int_t ADCsat=3100;       // saturation value for the ADCs
-  Int_t ADCtof[48];
-
-
-  // ---- introduce scale factors to tune simul ADC to real data   24-oct DC
-  //   Float_t ScaleFact[48]={0.18,0.22,0.35,0.26,0.47,0.35,0.31,0.37,
-  // 			 0.44,0.23,0.38,0.60,0.39,0.29,0.40,0.23,
-  // 			 0.30,0.66,0.22,1.53,0.17,0.55,
-  // 			 0.84,0.19,0.21,1.64,0.62,0.13,
-  // 			 0.18,0.15,0.10,0.14,0.14,0.14,0.14,0.12,
-  // 			 0.26,0.18,0.25,0.23,0.20,0.40,0.19,0.23,0.25,0.23,0.25,0.20};
-  
-  // new scale factors: WM 30-Oct-07
-  //   Float_t ScaleFact[48]={0.35,0.41,0.32,0.34,0.58,0.47,0.42,0.44,
-  // 			 0.50,0.34,0.50,0.50,0.51,0.42,0.46,0.25,
-  // 			 0.20,0.38,0.29,0.49,0.24,0.68,
-  // 			 0.30,0.26,0.28,0.79,0.31,0.12,
-  // 			 0.25,0.21,0.14,0.20,
-  // 			 0.16,0.17,0.19,0.18,
-  // 			 0.34,0.27,0.34,0.31,0.25,0.57,
-  // 			 0.24,0.34,0.34,0.32,0.31,0.30};
-  
-  Float_t ScaleFact[48]={0.39, 0.49, 0.38, 0.40, 0.65, 0.51, 0.43,
-			 0.49, 0.58, 0.38, 0.53, 0.57, 0.53, 0.45, 0.49, 0.16,
-			 0.15, 0.44, 0.28, 0.57, 0.26, 0.72, 0.37, 0.29, 0.30, 0.89,
-			 0.37, 0.08,  0.27, 0.23, 0.12, 0.22, 0.15, 0.16, 0.21,
-			 0.19, 0.41, 0.32, 0.39, 0.38, 0.28, 0.66, 0.28, 0.40, 0.39, 0.40, 0.37, 0.35 };
-
-  for(Int_t i=0; i<48; i++){
-    QevePmt_pC[i] = 0;
-    ADCtof[i]=0;
-  }
-  
-  // ------ read calibration file (get A1, A2, lambda1, lambda2)
-  ifstream fileTriggerCalib;
-  TString ftrigname="TrigCalibParam.txt";
-  fileTriggerCalib.open(ftrigname.Data());
-  if ( !fileTriggerCalib ) {
-    printf("debug: no trigger calib file!\n");
-    return(-117); //check output!
-  };
-  Float_t atte1[48],atte2[48],lambda1[48],lambda2[48];
-  Int_t temp=0;
-  // correct readout WM Oct '07
-  for(Int_t i=0; i<48; i++){
-    fileTriggerCalib >> temp;
-    fileTriggerCalib >> atte1[i];
-    fileTriggerCalib >> lambda1[i];
-    fileTriggerCalib >> atte2[i];
-    fileTriggerCalib >> lambda2[i];
-    fileTriggerCalib >> temp;
-  }
-  fileTriggerCalib.close();
-
-  Int_t ip, ipad;
-  //Int_t ipmt;
-  Int_t pmtleft=0, pmtright=0;
-  Int_t *pl, *pr;
-  pl = &pmtleft;
-  pr = &pmtright;
-
-  // TDC variables:
-  Int_t TDClast=4095;      // no signal --> TDC ch=4095
-  Int_t TDCint[48];
-  Float_t  tdc[48],tdc1[48],tdcpmt[48];
-  for(Int_t i=0; i<48; i++) {
-    tdcpmt[i] = 1000.; 
-    tdc[i]  = 0.;            // 18-oct WM
-    tdc1[i] = 0.;            // 18-oct WM
-  } 
-
-  Float_t thresh=10.; // to be defined better... (Wolfgang)
-
-  // === TDC: simulate timing for each paddle
-  //Float_t dt1 = 285.e-12 ;   // single PMT resolution
-  Float_t dt1 = 425.e-12 ;   // single PMT resolution (WM, Nov'07)
-  Float_t tdcres[50],c1_S[50],c2_S[50],c3_S[50]; 
-  for(Int_t j=0;j<48;j++)  tdcres[j] = 50.E-12;   // TDC resolution 50 picosec
-  for(Int_t j=0;j<48;j++)  c1_S[j] = 500.;  // cable length in channels
-  for(Int_t j=0;j<48;j++)  c2_S[j] = 0.;
-  for(Int_t j=0;j<48;j++)  c3_S[j] = 1000.;
-  for(Int_t j=0;j<48;j++)  c1_S[j] = c1_S[j]*tdcres[j];   // cable length in sec 
-  for(Int_t j=0;j<48;j++)  c2_S[j] = c2_S[j]*tdcres[j];
-  // ih = 0 + i1;  // not used?? (Silvio)
-
-  /* **********************************  start loop over hits */
-  
-  for(Int_t nh=0; nh<Nthtof; nh++){
-    
-    Float_t s_l_g[6] = {8.0, 8.0, 20.9, 22.0, 9.8, 8.3 };  // length of the lightguide
-    Float_t t1,t2,veff,veff1,veff0 ;
-    veff0 = 100.*1.0e8 ; // light velocity in the scintillator in m/sec
-    veff1 = 100.*1.5e8; // light velocity in the lightguide in m/sec
-    veff=veff0;         // signal velocity in the paddle
-
-    t1 = Timetof[nh] ;  // Start
-    t2 = Timetof[nh] ;
-
-    // Donatella: redefinition plane and pad for vectors in C
-    ip = Ipltof[nh]-1;
-    ipad = Ipaddle[nh]-1;
-    pmtleft=0; 
-    pmtright=0;
-
-    // WM: S12 paddles are "reversed" (Nov'07)
-    if (ip==2)
-      if (ipad==0)
-	ipad=1;
-      else
-	ipad=0;
-
-   if (ip<6) {
-      Paddle2Pmt(ip, ipad, &pmtleft, &pmtright);
- 
-      // DC: evaluates mean position and path inside the paddle
-  
-      Float_t tpos=0.;
-      Float_t path[2] = {0., 0.};
-      //--- Strip in Y = S11,S22,S31 ------
-      if(ip==0 || ip==3 || ip==4)
-	tpos = (Yintof[nh]+Youttof[nh])/2.;
-      else 
-	if(ip==1 || ip==2 || ip==5)   //--- Strip in X for S12,S21,S32 
-	  tpos = (Xintof[nh]+Xouttof[nh])/2.;
-	else //if (ip!=6)
-	  printf("*** WARNING TOF: this option should never occur! (ip=%2i, nh=%2i)\n",ip,nh);
-
-      path[0]= tpos + dimel[ip]/2.;   // path to left PMT
-      path[1]= dimel[ip]/2.- tpos;    // path to right PMT
-
-      //  cout <<"Strip N. ="<< ipaddle <<" piano n.= "<< iplane <<" POSIZ = "<< tpos <<"\n";
-      
-      if (DEBUG) {
-	cout <<" plane "<<ip<<" strip # ="<< ipad <<" tpos  "<< tpos <<"\n";
-	cout <<"pmtleft, pmtright "<<pmtleft<<" "<<pmtright<<endl;
-      }
-      
-      // constant geometric factor, the rest is handled by the scaling factor
-      FGeo[0] =0.5;
-      FGeo[1] =0.5; 
-      //  FGeo[1] = atan(path[1]/dimes[ip])/6.28318; // fraction of photons toward left
-      //  FGeo[2] = atan(path[2]/dimes[ip])/6.28318; // toward right
-      
-      
-      //  Npho = Poisson(ERELTOF[nh])*Pho_keV*1e6  Poissonian fluctuations to be inserted-DC
-      Npho = Ereltof[nh]*Pho_keV*1.0e6;  // Eloss in GeV 
-      
-      Float_t knorm[2]={0., 0.}; // Donatella
-      Float_t Atten[2]={0., 0.}; // Donatella
-      for(Int_t j=0; j<2; j++){
-	QhitPad_pC[j]= Npho*FGeo[j]*effi*pmGain*echarge*1.E12*ScaleFact[pmtleft+j];
-	// WM
-	knorm[j]=atte1[pmtleft+j]*exp(lambda1[pmtleft+j]*dimel[ip]/2.*pow(-1,j+1)) + 
-	  atte2[pmtleft+j]*exp(lambda2[pmtleft+j]*dimel[ip]/2.*pow(-1,j+1));
-	Atten[j]=atte1[pmtleft+j]*exp(tpos*lambda1[pmtleft+j]) + 
-	  atte2[pmtleft+j]*exp(tpos*lambda2[pmtleft+j]) ;
-	QhitPmt_pC[j]= QhitPad_pC[j]*Atten[j]/knorm[j];
-	//       	QhitPmt_pC[j]= QhitPad_pC[j]; //no attenuation
-
-
-	if (DEBUG) {
-	  cout<<"pmtleft "<<pmtleft<<" j "<<j<<endl;
-	  cout<<" atte1 "<<atte1[pmtleft+j]<<"lambda1 "<<lambda1[pmtleft+j]<<" atte2 "<<atte2[pmtleft+j]<<"lambda2 "<<lambda2[pmtleft+j] <<endl;    
-	  cout<<j<<" tpos "<<tpos<<" knorm "<<knorm[j]<<" "<<Atten[j]<<" "<<"QhitPmt_pC "<<QhitPmt_pC[j]<<endl;    
-	}
-      }
-      
-      if (DEBUG)
-	cout<<"Npho "<<Npho<<" QhitPmt_pC "<<QhitPmt_pC[0]<<" "<<QhitPmt_pC[1]<<endl;   
-      
-      QevePmt_pC[pmtleft]  += QhitPmt_pC[0];
-      QevePmt_pC[pmtright] += QhitPmt_pC[1];
-      
-      // TDC
-      // WM right and left <->
-      //      t2 = t2 + fabs(path[0]/veff) + s_l_g[ip]/veff1 ;  // Signal reaches PMT
-      //      t1 = t1 + fabs(path[1]/veff) + s_l_g[ip]/veff1;
-
-      t1 = t1 + fabs(path[0]/veff) + s_l_g[ip]/veff1;
-      t2 = t2 + fabs(path[1]/veff) + s_l_g[ip]/veff1 ;  // Signal reaches PMT
-      
-      t1 = gRandom->Gaus(t1,dt1); //apply gaussian error  dt
-      t2 = gRandom->Gaus(t2,dt1); //apply gaussian error  dt
-
-      t1 = t1 + c1_S[pmtleft] ;  // Signal reaches Discriminator ,TDC starts  to run
-      t2 = t2 + c1_S[pmtright] ;
-      
-      // check if signal is above threshold
-      // then check if tdcpmt is already filled by another hit...
-      // only re-fill if time is smaller
-      
-      if (QhitPmt_pC[0] > thresh) {
-	if (tdcpmt[pmtleft] == 1000.) {  // fill for the first time
-	  tdcpmt[pmtleft] = t1;
-	  tdc[pmtleft] = t1 + c2_S[pmtleft] ;  // Signal reaches Coincidence
-	}
-	if (tdcpmt[pmtleft] < 1000.) // is already filled!
-	  if (t1 <  tdcpmt[pmtleft]) {
-	    tdcpmt[pmtleft] = t1;
-	    t1 = t1 + c2_S[pmtleft] ;  // Signal reaches Coincidence
-	    tdc[pmtleft] = t1;
-	  }
-      }      
-      if (QhitPmt_pC[1] > thresh) {
-	if (tdcpmt[pmtright] == 1000.) {  // fill for the first time
-	  tdcpmt[pmtright] = t2;
-	  tdc[pmtright] = t2 + c2_S[pmtright] ;  // Signal reaches Coincidence
-	}
-	if (tdcpmt[pmtright] < 1000.)  // is already filled!
-	  if (t2 <  tdcpmt[pmtright]) {
-	    tdcpmt[pmtright] = t2;
-	    t2 = t2 + c2_S[pmtright] ;
-	    tdc[pmtright] = t2;
-	  }      
-      }
-
-      if (DEBUG)
-	cout<<nh<<" "<<Timetof[nh]<<" "<<t1<<" "<<t2<<endl;
-
-    } // ip < 6
-
-  }; // ****************************************       end loop over hits
-  
-  // ======  ADC ======
-
-  for(Int_t i=0; i<48; i++){
-    if (QevePmt_pC[i] < 800.)  ADCtof[i]= (Int_t)(ADC_pC0A + ADC_pC1A*QevePmt_pC[i] + ADC_pC2A*pow(QevePmt_pC[i],2) + ADC_pC3A*pow(QevePmt_pC[i],3));
-    if (QevePmt_pC[i] > 800.)  ADCtof[i]= (Int_t)(ADC_pC0B + ADC_pC1B*QevePmt_pC[i] + ADC_pC2B*pow(QevePmt_pC[i],2) + ADC_pC3B*pow(QevePmt_pC[i],3));
-    if (QevePmt_pC[i] > 2485.) ADCtof[i]= (Int_t)(1758. + 0.54*QevePmt_pC[i]);  //assuming a fictional 0.54 ch/pC above ADCsat
-    if (ADCtof[i]>ADCsat) ADCtof[i]=ADCsat;
-    if (QevePmt_pC[i] < pCthres)  ADCtof[i]= ADClast;
-    if (ADCtof[i] < 0) ADCtof[i]=ADClast;    
-    if (ADCtof[i] > ADClast) ADCtof[i]=ADClast;
-  }
-
-//   for(Int_t i=0; i<48; i++){
-//     if(QevePmt_pC[i] >= pCthres){
-//       ADCtof[i]= (Int_t)(ADC_pC0 + ADC_pC1*QevePmt_pC[i] + ADC_pC2*pow(QevePmt_pC[i],2) + ADC_pC3*pow(QevePmt_pC[i],3));
-//   } else 
-//       ADCtof[i]= ADClast;   
-// }
-  
-// // ---- introduce scale factors to tune simul ADC to real data   24-oct DC
-
-//   for(Int_t i=0; i<48; i++){
-//     if(ADCtof[i] != ADClast){
-//                              //      printf("%3d, %4d, %4.2f\n",i, ADCtof[i],ScaleFact[i]);
-//       ADCtof[i]= Int_t (ADCtof[i]*ScaleFact[i]);
-//                              //      printf("%3d, %4d,\n",i, ADCtof[i]);
-//     }
-//   }
-
-//   for(Int_t i=0; i<48; i++){
-//     if(ADCtof[i] != ADClast){
-//       if(ADCtof[i]> ADCsat) ADCtof[i]=ADCsat;
-//       else if(ADCtof[i]< 0) ADCtof[i]=ADClast;    
-//   }
-//   }    
-
-
-// ======  build  TDC coincidence  ======
-
-  Float_t t_coinc = 0;
-  Int_t ilast = 100;
-  for (Int_t ii=0; ii<48;ii++)
-    if (tdc[ii] > t_coinc) {
-      t_coinc = tdc[ii];
-      ilast = ii;
-    }
-  
-  //     cout<<ilast<<" "<<t_coinc<<endl;
-  //     At t_coinc  trigger condition is fulfilled
-  
-  for (Int_t ii=0; ii<48;ii++){
-    //      if (tdc[ii] != 0) tdc1[ii] = t_coinc - tdc[ii];   // test 1
-    if (tdc[ii] != 0) tdc1[ii] = t_coinc - tdcpmt[ii];  // test 2
-    tdc1[ii] = tdc1[ii]/tdcres[ii];                     // divide by TDC resolution 
-    if (tdc[ii] != 0) tdc1[ii] = tdc1[ii] + c3_S[ii];  // add cable length c3
-  } // missing parenthesis inserted! (Silvio)
-
-  for(Int_t i=0; i<48; i++){
-    if(tdc1[i] != 0.){
-      TDCint[i]=(Int_t)tdc1[i];
-      if (TDCint[i]>4093) TDCint[i]=TDClast;  // 18-oct WM
-      if (DEBUG)
-	cout<<i<<" "<<TDCint[i]<<endl;
-      //ADC[i]= ADC_pC * QevePmt_pC[i] + ADCoffset;
-      //if(ADC[i]> ADClast) ADC[i]=ADClast;
-    } else
-      TDCint[i]= TDClast;
-  }
-
-  if (DEBUG)
-    cout<<"-----------"<<endl;
- 
-
-//------ use channelmap  18-oct WM
-
-  Int_t  channelmap[] =  {3,21,11,29,19,45,27,37,36,28,44,20,5,12,13,4,
-			6,47,14,39,22,31,30,23,38,15,46,7,0,33,16,24,
-			8,41,32,40,25,17,34,9,42,1,2,10,18,26,35,43};
-
-  Int_t ADChelp[48];
-  Int_t TDChelp[48];
-
-  for(Int_t i=0; i<48; i++){
-    Int_t ii=channelmap[i];
-    ADChelp[ii]= ADCtof[i];
-    TDChelp[ii]= TDCint[i];
-  }
-
-  for(Int_t i=0; i<48; i++){
-    ADCtof[i]= ADChelp[i];
-    TDCint[i]= TDChelp[i];
-  }
-
-
-// ======  write fDataTof  =======
-
-
-//  UChar_t tdcadd[8]={1,0,3,2,5,4,7,6};  (coded in 3 bit)
-  UChar_t Ctrl3bit[8]={32,0,96,64,160,128,224,192};  // DC (msb in 8 bit word )
-   
-  UChar_t tofBin;
-  for (Int_t j=0; j < 12; j++){   // loop on TDC #12
-    Int_t j12=j*23;               // for each TDC 23 bytes (8 bits)
-    fDataTof[j12+0]=0x00;   // TDC_ID
-    fDataTof[j12+1]=0x00;   // EV_COUNT
-    fDataTof[j12+2]=0x00;   // TDC_MASK (1)
-    fDataTof[j12+3]=0x00;   // TDC_MASK (2)
-    for (Int_t k=0; k < 4; k++){   // for each TDC 4 channels (ADC+TDC)
-          
-      Int_t jk12=j12+4*k;         // ADC,TDC channel (0-47)
-
-      tofBin =(UChar_t)(ADCtof[k+4*j]/256);   // ADC# (msb) 
-      fDataTof[jk12+4] = Bin2GrayTof(tofBin,fDataTof[jk12+4]);
-      /* control bits inserted here, after the bin to gray conv - DC*/ 
-      fDataTof[jk12+4] = Ctrl3bit[2*k] | fDataTof[jk12+4];
-      tofBin=(UChar_t)(ADCtof[k+4*j]%256);   // ADC# (lsb)
-      fDataTof[jk12+5] = Bin2GrayTof(tofBin,fDataTof[jk12+5]);
-      tofBin=(UChar_t)(TDCint[k+4*j]/256);   // TDC# (msb)
-      fDataTof[jk12+6]=Bin2GrayTof(tofBin,fDataTof[jk12+6]);
-      /* control bits inserted here, after the bin to gray conv - DC*/ 
-      fDataTof[jk12+6] = Ctrl3bit[2*k+1] | fDataTof[jk12+6];
-      tofBin=(UChar_t)(TDCint[k+4*j]%256);   // TDC# (lsb)
-      fDataTof[jk12+7]=Bin2GrayTof(tofBin,fDataTof[jk12+7]);
-    };
-    fDataTof[j12+20]=0x00;   // TEMP1
-    fDataTof[j12+21]=0x00;   // TEMP2
-    fDataTof[j12+22]= EvaluateCrcTof(pTof);   // CRC
-    pTof+=23;
-  };
-  return(0);
-};
-
-
-UChar_t Digitizer::Bin2GrayTof(UChar_t binaTOF,UChar_t grayTOF){
-  union graytof_data {
-    UChar_t word;
-    struct bit_field {
-      unsigned b0:1;
-      unsigned b1:1;
-      unsigned b2:1;
-      unsigned b3:1;
-      unsigned b4:1;
-      unsigned b5:1;
-      unsigned b6:1;
-      unsigned b7:1;
-    } bit;
-  } bi,gr;
-  //
-  bi.word = binaTOF;
-  gr.word = grayTOF;
-  // 
-  gr.bit.b0 = bi.bit.b1 ^ bi.bit.b0;
-  gr.bit.b1 = bi.bit.b2 ^ bi.bit.b1;
-  gr.bit.b2 = bi.bit.b3 ^ bi.bit.b2;
-  gr.bit.b3 = bi.bit.b3;
-  //
-  /* bin to gray conversion 4 bit per time*/
-  //
-  gr.bit.b4 = bi.bit.b5 ^ bi.bit.b4;
-  gr.bit.b5 = bi.bit.b6 ^ bi.bit.b5;
-  gr.bit.b6 = bi.bit.b7 ^ bi.bit.b6;
-  gr.bit.b7 = bi.bit.b7;
-  //
-  return(gr.word); 
-}
-
-UChar_t Digitizer::EvaluateCrcTof(UChar_t *pTof) {
-  Bool_t DEBUG=false;
-  if (DEBUG)
-    return(0x00);
-
-  UChar_t crcTof=0x00;
-  UChar_t *pc=&crcTof, *pc2;
-  pc2=pTof;
-  for (Int_t jp=0; jp < 23; jp++){
-    //crcTof = crc8(...)
-    Crc8Tof(pc2++,pc);
-    //    printf("%2i --- %x\n",jp,crcTof);
-  }
-  return(crcTof);
-}
-
-void Digitizer::Crc8Tof(UChar_t *oldCRC, UChar_t *crcTof){
-  union crctof_data {
-    UChar_t word;
-    struct bit_field {
-      unsigned b0:1;
-      unsigned b1:1;
-      unsigned b2:1;
-      unsigned b3:1;
-      unsigned b4:1;
-      unsigned b5:1;
-      unsigned b6:1;
-      unsigned b7:1;
-    } bit;
-  } c,d,r;
-
-  c.word = *oldCRC;
-  //d.word = *newCRC;
-  d.word = *crcTof;
-  r.word = 0;
-
-  r.bit.b0 = c.bit.b7 ^ c.bit.b6 ^ c.bit.b0 ^ 
-             d.bit.b0 ^ d.bit.b6 ^ d.bit.b7;
-
-  r.bit.b1 = c.bit.b6 ^ c.bit.b1 ^ c.bit.b0 ^ 
-             d.bit.b0 ^ d.bit.b1 ^ d.bit.b6;
-
-  r.bit.b2 = c.bit.b6 ^ c.bit.b2 ^ c.bit.b1 ^ c.bit.b0 ^
-             d.bit.b0 ^ d.bit.b1 ^ d.bit.b2 ^ d.bit.b6;
-
-  r.bit.b3 = c.bit.b7 ^ c.bit.b3 ^ c.bit.b2 ^ c.bit.b1 ^ 
-             d.bit.b1 ^ d.bit.b2 ^ d.bit.b3 ^ d.bit.b7;
-
-  r.bit.b4 = c.bit.b4 ^ c.bit.b3 ^ c.bit.b2 ^
-             d.bit.b2 ^ d.bit.b3 ^ d.bit.b4;
-
-  r.bit.b5 = c.bit.b5 ^ c.bit.b4 ^ c.bit.b3 ^
-             d.bit.b3 ^ d.bit.b4 ^ d.bit.b5;
-
-  r.bit.b6 = c.bit.b6 ^ c.bit.b5 ^ c.bit.b4 ^ 
-             d.bit.b4 ^ d.bit.b5 ^ d.bit.b6;
-
-  r.bit.b7 = c.bit.b7 ^ c.bit.b6 ^ c.bit.b5 ^ 
-             d.bit.b5 ^ d.bit.b6 ^ d.bit.b7 ;
-
-  *crcTof=r.word;
-  //return r.word; 
-};
-
-//void Digitizer::Paddle2Pmt(Int_t plane, Int_t paddle, Int_t* &pmtleft, Int_t* &pmtright){
-void Digitizer::Paddle2Pmt(Int_t plane, Int_t paddle, Int_t *pl, Int_t *pr){
-  //* @param plane    (0 - 5)
-  //* @param paddle   (plane=0, paddle = 0,...5)
-  //* @param padid    (0 - 23)  
-  //
-  Int_t padid=-1;
-  Int_t pads[6]={8,6,2,2,3,3};
-  //
-  Int_t somma=0;
-  Int_t np=plane;
-  for(Int_t j=0; j<np; j++)
-    somma+=pads[j];
-  padid=paddle+somma;
-  *pl = padid*2;
-  //  *pr = *pr + 1; 
-  *pr = *pl + 1; // WM
-};
-
-void Digitizer::DigitizeAC() {
-  // created:  J. Conrad, KTH
-  // modified: S. Orsi, INFN Roma2
-  // fDataAC[0-63]:   main AC board
-  // fDataAC[64-127]: extra AC board
-
-  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);
-
-  // 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;
-  };
-  
-  // coincidences are dummy (increment by 1 at each event)
-  // for (UInt_t i=0; i<=7; i++){
-  //    fDataAC[i+42] = 0x0000;
-  //    fDataAC[i+106] = 0x0000;
-  //   }
-  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)
-  /*     
-    fDataAC[51] = 0x006C;
-    fDataAC[52] = 0x6C6C;
-    fDataAC[115] = 0x006C;
-    fDataAC[116] = 0x6C6C;
-  */
-  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;
-  }
-  
-  // We activate all branches. Once the digitization algorithm is determined 
-  // only the branches that involve needed information will be activated
-  
-  fhBookTree->SetBranchAddress("Ievnt",&Ievnt);
-  fhBookTree->SetBranchStatus("Nthcat",1);
-  fhBookTree->SetBranchStatus("Iparcat",1);
-  fhBookTree->SetBranchStatus("Icat",1);
-  fhBookTree->SetBranchStatus("Xincat",1);
-  fhBookTree->SetBranchStatus("Yincat",1);
-  fhBookTree->SetBranchStatus("Zincat",1);
-  fhBookTree->SetBranchStatus("Xoutcat",1);
-  fhBookTree->SetBranchStatus("Youtcat",1);
-  fhBookTree->SetBranchStatus("Zoutcat",1);
-  fhBookTree->SetBranchStatus("Erelcat",1);
-  fhBookTree->SetBranchStatus("Timecat",1);
-  fhBookTree->SetBranchStatus("Pathcat",1);
-  fhBookTree->SetBranchStatus("P0cat",1);
-  fhBookTree->SetBranchStatus("Nthcas",1);
-  fhBookTree->SetBranchStatus("Iparcas",1);
-  fhBookTree->SetBranchStatus("Icas",1);
-  fhBookTree->SetBranchStatus("Xincas",1);
-  fhBookTree->SetBranchStatus("Yincas",1);
-  fhBookTree->SetBranchStatus("Zincas",1);
-  fhBookTree->SetBranchStatus("Xoutcas",1);
-  fhBookTree->SetBranchStatus("Youtcas",1);
-  fhBookTree->SetBranchStatus("Zoutcas",1);
-  fhBookTree->SetBranchStatus("Erelcas",1);
-  fhBookTree->SetBranchStatus("Timecas",1);
-  fhBookTree->SetBranchStatus("Pathcas",1);
-  fhBookTree->SetBranchStatus("P0cas",1);
-  fhBookTree->SetBranchStatus("Nthcard",1);
-  fhBookTree->SetBranchStatus("Iparcard",1);
-  fhBookTree->SetBranchStatus("Icard",1);
-  fhBookTree->SetBranchStatus("Xincard",1);
-  fhBookTree->SetBranchStatus("Yincard",1);
-  fhBookTree->SetBranchStatus("Zincard",1);
-  fhBookTree->SetBranchStatus("Xoutcard",1);
-  fhBookTree->SetBranchStatus("Youtcard",1);
-  fhBookTree->SetBranchStatus("Zoutcard",1);
-  fhBookTree->SetBranchStatus("Erelcard",1);
-  fhBookTree->SetBranchStatus("Timecard",1);
-  fhBookTree->SetBranchStatus("Pathcard",1);
-  fhBookTree->SetBranchStatus("P0card",1);
-
-  // 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
-  //
-  // NB: the PMT positions are needed!
-
-  // 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) 
-      SumEcas[Icas[k]] += Erelcas[k];
-  };
-
-  // 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;
-  //   }
-};
-
-
-void Digitizer::DigitizeS4(){
-  Int_t DEBUG=0;
-  // creato:  S. Borisov, INFN Roma2 e MEPHI, Sett 2007
-  TString ciao,modo="ns";
-  Int_t i,j,t,NdF,pmt,NdFT,S4,S4v=0,S4p=32;
-  Float_t E0,E1=1e-6,Ert,X,Y,Z,x,y,z,V[3],Xs[2],Ys[2],Zs[2],Yp[6],q,w,p=0.1,l,l0=500;
-  Xs[0]=-24.1;
-  Xs[1]=24.1;
-  Ys[0]=-24.1;
-  Ys[1]=24.1;
-  Zs[0]=-0.5;
-  Zs[1]=0.5;
-  Yp[0]=-20.;
-  Yp[2]=-1.;
-  Yp[4]=17.;
-  for(i=0;i<3;i++)
-    Yp[2*i+1]=Yp[2*i]+3;
-  srand(time(NULL));
-  // --- activate branches:
-  fhBookTree->SetBranchStatus("Nthtof",1);
-  fhBookTree->SetBranchStatus("Ipltof",1);
-  fhBookTree->SetBranchStatus("Ipaddle",1);
-  
-  fhBookTree->SetBranchStatus("Xintof",1);
-  fhBookTree->SetBranchStatus("Yintof",1);
-  fhBookTree->SetBranchStatus("Xouttof",1);
-  fhBookTree->SetBranchStatus("Youttof",1);
-  
-  fhBookTree->SetBranchStatus("Ereltof",1);
-  fhBookTree->SetBranchStatus("Timetof",1);
-  NdFT=0;
-  Ert=0;
-  for(i=0;i<Nthtof;i++){
-    if(Ipltof[i]!=6) continue;
-    Ert+=Ereltof[i];
-           
-    if(modo=="ns") continue;
-    NdF=Int_t(Ereltof[i]/E1);
-    NdFT=0;
-    X=Xintof[i];
-    Y=Yintof[i];
-    Z=(Float_t)(random())/(Float_t)(0x7fffffff)-0.5;
-    //cout<<"XYZ "<<X<<"  "<<Y<<"  "<<Z<<endl;
-    for(j=0;j<NdF;j++){
-      q=(Float_t)random()/(Float_t)0x7fffffff;
-      w=(Float_t)random()/(Float_t)0x7fffffff;
-      // cout<<"qw "<<q<<" "<<w<<endl;
-      V[0]=p*cos(6.28318*q);
-      V[1]=p*sin(6.28318*q);
-      V[2]=p*(2.*w-1.);
-      pmt=0;
-      x=X;
-      y=Y;
-      z=Z;
-      while(pmt==0 && (x>Xs[0] && x<Xs[1])&&(y>Ys[0] && y<Ys[1])&&(z>Zs[0] && z<Zs[1])){
-	l=0;
-	while(pmt==0 && (x>Xs[0] && x<Xs[1])&&(y>Ys[0] && y<Ys[1])&&(z>Zs[0] && z<Zs[1])){
-	  x+=V[0];
-	  y+=V[1];
-	  z+=V[2];
-	  l+=p;
-	  //cout<<x<<"  "<<y<<"  "<<z<<"  "<<l<<endl;
-	  //cin>>ciao;
-	}
-	if((x<Xs[0]+p || x>Xs[1]-p)&&(y>Ys[0]+p && y<Ys[1]-p)&&(z>Zs[0]+p && z<Zs[1]-p)){
-	  for(t=0;t<3;t++){
-	    if(y>=Yp[2*t] && y<Yp[2*t+1]){
-	      if(pmt==0)NdFT++;
-	      pmt=1;
-	      //cout<<NdFT<<endl;
-	      break;
-	    }
-	  }
-	  if(pmt==1)break;
-	  V[0]=-V[0];
-	}
-	q=(Float_t)random()/(Float_t)0x7fffffff;
-	w=1-exp(-l/l0);
-	if(q<w)break;
-	q=(Float_t)random()/(Float_t)0x7fffffff;
-	w=0.5;
-	if(q<w)break;
-	if((x>Xs[0]+p && x<Xs[1]-p)&&(y<Ys[0]+p || y>Ys[1]-p)&&(z>Zs[0]+p && z<Zs[1]-p))V[1]=-V[1];
-	if((x>Xs[0]+p && x<Xs[1]-p)&&(y>Ys[0]+p && y<Ys[1]-p)&&(z<Zs[0]+p || z>Zs[1]-p))V[2]=-V[2];
-	x+=V[0];
-	y+=V[1];
-	z+=V[2];
-	l=0;
-	//cout<<x<<"  "<<y<<"  "<<z<<"  "<<l<<endl;
-		//cin>>ciao;
-      }
-    }
-  }
-  Ert=Ert/0.002;
-  q=(Float_t)(random())/(Float_t)0x7fffffff;
-  w=0.7;
-  //E0=(Float_t)(4064./7.);
-  E0=4064./7.;
-  if(Ert<1) S4=0;
-  else S4=(Int_t)(4064.*(1.-exp(-(Ert-1.)/E0)));
-  i=S4/4;
-  if(S4%4==0) 
-    S4v=S4+S4p;
-  else if(S4%4==1){
-    if(q<w) S4v=S4-1+S4p;
-    else S4v=S4+1+S4p;
-  } else if(S4%4==2) S4v=S4+S4p;
-  else if(S4%4==3){
-    if(q<w) S4v=S4+1+S4p;
-    else S4v=S4-1+S4p;
-  }
-  if (DEBUG)
-    cout<<"Ert_S4 = " << Ert << " --- S4v = " << S4v << endl;
-  fDataS4[0]=S4v;//0xf028;
-  fDataS4[1]=0xd800;
-  fDataS4[2]=0x0300;
-  //cout<<"  PMT  "<<NdFT<<"  "<<NdF<<endl;
-  //cin>>ciao;
-}
-
-
-
-void Digitizer::DigitizeND(){
-  // creato:  S. Borisov, INFN Roma2 e MEPHI, Sett 2007
-  Int_t i=0;
-  UShort_t NdN=0;
-  fhBookTree->SetBranchStatus("Nthnd",1);
-  fhBookTree->SetBranchStatus("Itubend",1);
-  fhBookTree->SetBranchStatus("Iparnd",1);  
-  fhBookTree->SetBranchStatus("Xinnd",1);
-  fhBookTree->SetBranchStatus("Yinnd",1);
-  fhBookTree->SetBranchStatus("Zinnd",1);
-  fhBookTree->SetBranchStatus("Xoutnd",1);
-  fhBookTree->SetBranchStatus("Youtnd",1);
-  fhBookTree->SetBranchStatus("Zoutnd",1);
-  fhBookTree->SetBranchStatus("Erelnd",1);
-  fhBookTree->SetBranchStatus("Timend",1);
-  fhBookTree->SetBranchStatus("Pathnd",1);
-  fhBookTree->SetBranchStatus("P0nd",1);
-  //cout<<"n="<<Nthnd<<"  "<<NdN<<"\n";
-  for(i=0;i<Nthnd;i++){
-    if(Iparnd[i]==13){
-      NdN++;
-    }
-  }
-  //NdN=100; //only for debug
-
-  for(i=0;i<3;i++){
-    fDataND[2*i]=0x0000;
-    fDataND[2*i+1]=0x010F;
-  }
-  fDataND[0]=0xFF00 & (256*NdN);
-}
-
-
-void Digitizer::DigitizeDummy() {
-
-  fhBookTree->SetBranchStatus("Enestrip",1);
-
-  // dumy header
-  fDataDummy[0] = 0xCAAA;
-
-  for (Int_t i=1; i<fDummybuffer; i++){
-    fDataDummy[i] = 0xFFFF;
-    //   printf("%0x  ",fDataDummy[i]);   
-    //if ((i+1)%8 ==0) cout << endl;
-  }
-};
-
-
-void Digitizer::WriteData(){
-
-  // Routine that writes the data to a binary file
-  // PSCU data are already swapped
-  fOutputfile.write(reinterpret_cast<char*>(fDataPSCU),sizeof(UShort_t)*fPSCUbuffer);
-  // TRG
-  fOutputfile.write(reinterpret_cast<char*>(fDataTrigger),sizeof(UChar_t)*153);
-  // TOF
-  fOutputfile.write(reinterpret_cast<char*>(fDataTof),sizeof(UChar_t)*276);
-  // AC
-  UShort_t temp[1000000];
-  memset(temp,0,sizeof(UShort_t)*1000000);
-  swab(fDataAC,temp,sizeof(UShort_t)*fACbuffer);  // WE MUST SWAP THE BYTES!!!
-  fOutputfile.write(reinterpret_cast<char*>(temp),sizeof(UShort_t)*fACbuffer);
-  // CALO
-  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);
-  // TRK
-  memset(temp,0,sizeof(UShort_t)*1000000);
-  swab(fDataTrack,temp,sizeof(UShort_t)*fTracklength);  // WE MUST SWAP THE BYTES!!!
-  fOutputfile.write(reinterpret_cast<char*>(temp),sizeof(UShort_t)*fTracklength); 
-  fTracklength=0; 
-   // padding to 64 bytes
-  //
-  if ( fPadding ){
-    fOutputfile.write(reinterpret_cast<char*>(fDataPadding),sizeof(UChar_t)*fPadding);
-  };
-  // S4
-  memset(temp,0,sizeof(UShort_t)*1000000);
-  swab(fDataS4,temp,sizeof(UShort_t)*fS4buffer);  // WE MUST SWAP THE BYTES!!!
-  fOutputfile.write(reinterpret_cast<char*>(temp),sizeof(UShort_t)*fS4buffer);
-  // ND
-  memset(temp,0,sizeof(UShort_t)*1000000);
-  swab(fDataND,temp,sizeof(UShort_t)*fNDbuffer);  // WE MUST SWAP THE BYTES!!!
-  fOutputfile.write(reinterpret_cast<char*>(temp),sizeof(UShort_t)*fNDbuffer);
-};
-
-
-void Digitizer::ReadData(){
-
-  UShort_t InData[64];
-
-  // for debuggigng purposes only, write your own routine if you like (many
-  // hardwired things.
-  
-  ifstream InputFile;
-
-  // if (!InputFile) {
-
-  //     std::cout << "ERROR" << endl;
-  //     // An error occurred!
-  //     // myFile.gcount() returns the number of bytes read.
-  //     // calling myFile.clear() will reset the stream state
-  //     // so it is usable again.
-  //   };
-
-
-  
-  //InputFile.seekg(0);
-
-  InputFile.open(fFilename, ios::in | ios::binary);
-  //    fOutputfile.seekg(0);
-  if (!InputFile.is_open()) std::cout << "ERROR" << endl;
-
-  InputFile.seekg(0);
-  
-  for (Int_t k=0; k<=1000; k++){
-    InputFile.read(reinterpret_cast<char*>(InData),384*sizeof(UShort_t));
-
-    std::cout << "Read back: " << endl << endl;
-
-    for (Int_t i=0; i<=384; i++){
-      printf("%4x ", InData[i]);   
-      if ((i+1)%8 ==0) cout << endl;
-    }
-
-  }
-  cout << endl;
-  InputFile.close();
-
-};
-
-
-
-void Digitizer::DigitizeTrack() {
-//std:: cout << "Entering DigitizeTrack " << endl;
-Float_t  AdcTrack[fNviews][fNstrips_view];  //  Vector of strips to be compressed
-
-Int_t Iview;
-Int_t Nstrip;
-
-  for (Int_t j=0; j<fNviews;j++) {
-
-    for (Int_t i=0; i<fNladder;i++) {
-
-      Float_t commonN1=gRandom->Gaus(0.,fSigmaCommon);
-      Float_t commonN2=gRandom->Gaus(0.,fSigmaCommon);
-      for (Int_t k=0; k<fNstrips_ladder;k++) {
-      Nstrip=i*fNstrips_ladder+k;
-      AdcTrack[j][Nstrip]=gRandom->Gaus(fPedeTrack[j][Nstrip],fSigmaTrack[j][Nstrip]);
-      if(k<4*128) {AdcTrack[j][Nstrip] += commonN1;}  // full correlation of 4 VA1 Com. Noise
-      else {AdcTrack[j][Nstrip] += commonN2;}   // full correlation of 4 VA1 Com. Noise
-
-      };
-      
-        
-    };
-
-
-  }; 
-
-
-  fhBookTree->SetBranchStatus("Nstrpx",1);
-  fhBookTree->SetBranchStatus("Npstripx",1);
-  fhBookTree->SetBranchStatus("Ntstripx",1);
-  fhBookTree->SetBranchStatus("Istripx",1);
-  fhBookTree->SetBranchStatus("Qstripx",1);
-  fhBookTree->SetBranchStatus("Xstripx",1);
-  fhBookTree->SetBranchStatus("Nstrpy",1);
-  fhBookTree->SetBranchStatus("Npstripy",1);
-  fhBookTree->SetBranchStatus("Ntstripy",1);
-  fhBookTree->SetBranchStatus("Istripy",1);
-  fhBookTree->SetBranchStatus("Qstripy",1);
-  fhBookTree->SetBranchStatus("Ystripy",1);
-
-
-
-  Float_t ADCfull;
-  Int_t iladd=0;
-  for (Int_t ix=0; ix<Nstrpx;ix++) {
-  Iview=Npstripx[ix]*2-1;
-  Nstrip=(Int_t)Istripx[ix]-1;
-  if(Nstrip<fNstrips_ladder) iladd=0;
-  if((Nstrip>=fNstrips_ladder)&&(Nstrip<2*fNstrips_ladder)) iladd=1;
-  if((Nstrip>=2*fNstrips_ladder)&&(Nstrip<3*fNstrips_ladder)) iladd=2;
-  ADCfull=AdcTrack[Iview][Nstrip] += Qstripx[ix]*fMipCor[iladd][Iview];
-  AdcTrack[Iview][Nstrip] *= SaturationTrack(ADCfull);
-
-  };
-
-
-  for (Int_t iy=0; iy<Nstrpy;iy++) {
-  Iview=Npstripy[iy]*2-2;
-  Nstrip=(Int_t)Istripy[iy]-1;
-  if(Nstrip<fNstrips_ladder) iladd=0;
-  if((Nstrip>=fNstrips_ladder)&&(Nstrip<2*fNstrips_ladder)) iladd=1;
-  if((Nstrip>=2*fNstrips_ladder)&&(Nstrip<3*fNstrips_ladder)) iladd=2;
-  ADCfull=AdcTrack[Iview][Nstrip] -= Qstripy[iy]*fMipCor[iladd][Iview];
-  AdcTrack[Iview][Nstrip] *= SaturationTrack(ADCfull);
-
-  };  
-      
-CompressTrackData(AdcTrack);  // Compress and Digitize data of one Ladder  in turn for all ladders
-
-};
-
-
-
-void Digitizer::DigitizeTrackCalib(Int_t ii) {
-
-std:: cout << "Entering DigitizeTrackCalib " << ii << endl;
-if( (ii!=1)&&(ii!=2) ) {
- std:: cout << "error wrong DigitizeTrackCalib argument" << endl;
- return;
-}; 
-
-memset(fDataTrack,0,sizeof(UShort_t)*fTRACKbuffer);
-fTracklength=0;
-
-UShort_t Dato;
-
-Float_t dato1;
-Float_t dato2;
-Float_t dato3;
-Float_t dato4;
-
-UShort_t DatoDec;
-UShort_t DatoDec1;
-UShort_t DatoDec2;
-UShort_t DatoDec3;
-UShort_t DatoDec4;
-
-UShort_t EVENT_CAL;
-UShort_t PED_L1;
-UShort_t ReLength;
-UShort_t OveCheckCode;
-//UShort_t PED_L2;
-//UShort_t PED_L3HI;
-//UShort_t PED_L3LO;
-//UShort_t SIG_L1HI;
-//UShort_t SIG_L1LO;
-//UShort_t SIG_L2HI;
-//UShort_t SIG_L2LO;
-//UShort_t SIG_L3;
-//UShort_t BAD_L1;
-//UShort_t BAD_L2LO;
-//UShort_t BAD_L3HI;
-//UShort_t BAD_L3LO;
-//UShort_t FLAG;
-
-
-  Int_t DSPpos;
-  for (Int_t j=ii-1; j<fNviews;j+=2) {
-  UShort_t CkSum=0;
-  // here skip the dsp header and his trailer , to be written later
-  DSPpos=fTracklength;
-  fTracklength=fTracklength+13+3;
-
-
-    for (Int_t i=0; i<fNladder;i++) {
-      for (Int_t k=0; k<fNstrips_ladder;k++) {
-      // write in buffer the current LADDER
-      Dato=(UShort_t)fPedeTrack[j][i*fNstrips_ladder+k];
-      dato1=fPedeTrack[j][i*fNstrips_ladder+k]-Dato; 
-
-      DatoDec1=(UShort_t)(dato1*2);
-      dato2=dato1*2-DatoDec1;
-
-      DatoDec2=(UShort_t)(dato2*2);
-      dato3=dato2*2-DatoDec2;
-      
-      DatoDec3=(UShort_t)(dato3*2);
-      dato4=dato3*2-DatoDec3;
-      
-      DatoDec4=(UShort_t)(dato4*2);
-      
-      DatoDec=DatoDec1*0x0008+DatoDec2*0x0004+DatoDec3*0x0002+DatoDec4*0x0001;
-      fDataTrack[fTracklength]=( (Dato << 4) | (DatoDec & 0x000F) );
-      CkSum=CkSum^fDataTrack[fTracklength];
-      fTracklength++;
-      };
-
-      for (Int_t k=0; k<fNstrips_ladder;k++) {
-      // write in buffer the current LADDER
-      Dato=(UShort_t)fSigmaTrack[j][i*fNstrips_ladder+k]; 
-      dato1=fSigmaTrack[j][i*fNstrips_ladder+k]-Dato; 
-
-      DatoDec1=(UShort_t)(dato1*2);
-      dato2=dato1*2-DatoDec1;
-
-      DatoDec2=(UShort_t)(dato2*2);
-      dato3=dato2*2-DatoDec2;
-      
-      DatoDec3=(UShort_t)(dato3*2);
-      dato4=dato3*2-DatoDec3;
-      
-      DatoDec4=(UShort_t)(dato4*2);
-      
-      DatoDec=DatoDec1*0x0008+DatoDec2*0x0004+DatoDec3*0x0002+DatoDec4*0x0001;
-      
-      fDataTrack[fTracklength]=( (Dato << 4) | (DatoDec & 0x000F) );
-      CkSum=CkSum^fDataTrack[fTracklength];
-      fTracklength++;
-      };
-      
-      for (Int_t k=0; k<64;k++) {
-      fDataTrack[fTracklength]=0x0000;
-      CkSum=CkSum^fDataTrack[fTracklength];
-      fTracklength++;
-
-      };
-      // end ladder
-
-    // write in buffer the end ladder word
-    if(i==0) fDataTrack[fTracklength]=0x1807;
-    if(i==1) fDataTrack[fTracklength]=0x1808;
-    if(i==2) fDataTrack[fTracklength]=0x1809;
-    CkSum=CkSum^fDataTrack[fTracklength];
-    fTracklength++;
-
-    // write in buffer the TRAILER
-    ReLength=(UShort_t)((fNstrips_ladder*2+64+1)*2+3);
-    OveCheckCode=0x0000;
-
-    fDataTrack[fTracklength]=0x0000;
-    fTracklength++;
-  
-    fDataTrack[fTracklength]=(ReLength >> 8);
-    fTracklength++;
-  
-    fDataTrack[fTracklength]=( (ReLength << 8) | (OveCheckCode & 0x00FF) );
-    fTracklength++;
-
-    // end TRAILER        
-    };
-
-  // write in buffer the DSP header
-
-  fDataTrack[DSPpos]=(0xE800 | ( ((j+1) << 3) | 0x0005) );
-  
-  fDataTrack[DSPpos+1]=0x01A9;
-
-  fDataTrack[DSPpos+2]=0x8740;
-   
-  EVENT_CAL=0;
-  fDataTrack[DSPpos+3]=(0x1A00 | ( (0x03FF & EVENT_CAL)>> 1) );
-   
-  PED_L1=0;
-  fDataTrack[DSPpos+4]=( ((EVENT_CAL << 15) | 0x5002 ) | ((0x03FF & PED_L1) << 2) );
-  
-  // FROM HERE WE WRITE AS ALL VARIABLE apart CkSum are =0   
-
-  fDataTrack[DSPpos+5]=0x8014;
-  
-  fDataTrack[DSPpos+6]=0x00A0;
-  
-  fDataTrack[DSPpos+7]=0x0500;
-  
-  fDataTrack[DSPpos+8]=0x2801;
-  
-  fDataTrack[DSPpos+9]=0x400A;
-   
-  fDataTrack[DSPpos+10]=0x0050;
-  
-  CkSum=(CkSum >> 8)^(CkSum&0x00FF);
-  fDataTrack[DSPpos+11]=(0x0280 | (CkSum >> 3));
-  
-  fDataTrack[DSPpos+12]=(0x1FFF | (CkSum << 13) );
-
-  // end dsp header
-
-  // write in buffer the TRAILER
-  
-  ReLength=(UShort_t)((13*2)+3);
-  OveCheckCode=0x0000;
-  fDataTrack[DSPpos+13]=0x0000;
-
-  fDataTrack[DSPpos+14]=(ReLength >> 8);
-   
-  fDataTrack[DSPpos+15]=( (ReLength << 8) | (OveCheckCode & 0x00FF) );
-  
-  // end TRAILER
-
-
-
-   
-  // end DSP    
-  };    
-
-
-
-};
-
-void Digitizer::WriteTrackCalib() {
-
-
-std:: cout << " Entering WriteTrackCalib " << endl;
-
-fOutputfile.write(reinterpret_cast<char*>(fDataPSCU),sizeof(UShort_t)*fPSCUbuffer);
-
-UShort_t temp[1000000];
-memset(temp,0,sizeof(UShort_t)*1000000);
-swab(fDataTrack,temp,sizeof(UShort_t)*fTracklength);  // WE MUST SWAP THE BYTES!!!
-fOutputfile.write(reinterpret_cast<char*>(temp),sizeof(UShort_t)*fTracklength); 
-fTracklength=0; 
-if ( fPadding ){
-      fOutputfile.write(reinterpret_cast<char*>(fDataPadding),sizeof(UChar_t)*fPadding);
-};
-
-};
-
-
-void Digitizer::ClearTrackCalib() {
-
-std:: cout << "Entering ClearTrackCalib " << endl;
-     
-  
-};
-
-
-void Digitizer::LoadTrackCalib() {
-std:: cout << "Entering LoadTrackCalib " << endl;
-
-// Generate the pedestals and sigmas according to parametrization
-  for (Int_t j=0; j<fNviews;j++) {
-    for (Int_t i=0; i<fNstrips_view;i++) {
-    
-    if((j+1)%2==0) {
-    fPedeTrack[j][i]=gRandom->Gaus(fAvePedex,fSigmaPedex);
-    fSigmaTrack[j][i]=gRandom->Gaus(fAveSigmax,fSigmaSigmax);
-    };
-    if((j+1)%2==1) {
-    fPedeTrack[j][i]=gRandom->Gaus(fAvePedey,fSigmaPedey);
-    fSigmaTrack[j][i]=gRandom->Gaus(fAveSigmay,fSigmaSigmay);
-    };
-    
-    };
-  };    
-
-   
-  
-};
-
-void Digitizer::LoadMipCor() {
-std:: cout << "Entering LoadMipCor" << endl;
-  Float_t xfactor=1./151.6*1.04;
-  Float_t yfactor=1./152.1;
-
-  fMipCor[0][0]=140.02*yfactor;
-  fMipCor[0][1]=140.99*xfactor;
-  fMipCor[0][2]=134.48*yfactor;
-  fMipCor[0][3]=144.41*xfactor;
-  fMipCor[0][4]=140.74*yfactor;
-  fMipCor[0][5]=142.28*xfactor;
-  fMipCor[0][6]=134.53*yfactor;
-  fMipCor[0][7]=140.63*xfactor;
-  fMipCor[0][8]=135.55*yfactor;
-  fMipCor[0][9]=138.00*xfactor;
-  fMipCor[0][10]=154.95*yfactor;
-  fMipCor[0][11]=158.44*xfactor;
-  
-  
-  fMipCor[1][0]=136.07*yfactor;
-  fMipCor[1][1]=135.59*xfactor;
-  fMipCor[1][2]=142.69*yfactor;
-  fMipCor[1][3]=138.19*xfactor;
-  fMipCor[1][4]=137.35*yfactor;
-  fMipCor[1][5]=140.23*xfactor;
-  fMipCor[1][6]=153.15*yfactor;
-  fMipCor[1][7]=151.42*xfactor;
-  fMipCor[1][8]=129.76*yfactor;
-  fMipCor[1][9]=140.63*xfactor;
-  fMipCor[1][10]=157.87*yfactor;
-  fMipCor[1][11]=153.64*xfactor;
-
-  fMipCor[2][0]=134.98*yfactor;
-  fMipCor[2][1]=143.95*xfactor;
-  fMipCor[2][2]=140.23*yfactor;
-  fMipCor[2][3]=138.88*xfactor;
-  fMipCor[2][4]=137.95*yfactor;
-  fMipCor[2][5]=134.87*xfactor;
-  fMipCor[2][6]=157.56*yfactor;
-  fMipCor[2][7]=157.31*xfactor;
-  fMipCor[2][8]=141.37*yfactor;
-  fMipCor[2][9]=143.39*xfactor;
-  fMipCor[2][10]=156.15*yfactor;
-  fMipCor[2][11]=158.79*xfactor;
-
-/*
-  for (Int_t j=0; j<fNviews;j++) {
-    for (Int_t i=0; i<fNstrips_view;i++) {
-    fMipCor[j][i]=1.;
-    };
-  };    
-
-     
-*/  
-};
-
-void Digitizer::CompressTrackData(Float_t AdcTrack[fNviews][fNstrips_view]) {
-// copy of the corresponding compression fortran routine + new digitization
-// std:: cout << "Entering CompressTrackData " << endl;
-Int_t oldval=0;
-Int_t newval=0;
-Int_t trasmesso=0;
-Int_t ntrastot=0;
-Float_t real;
-Float_t inte;
-Int_t cercacluster=0;
-Int_t kt=0;
-static const int DSPbufferSize = 4000; // 13 bit buffer to be rearranged in 16 bit Track buffer
-UShort_t DataDSP[DSPbufferSize];  // 13 bit  buffer to be rearranged in 16 bit Track buffer
-UShort_t DSPlength;  // 13 bit buffer to be rearranged in 16 bit Track buffer
-
-memset(fDataTrack,0,sizeof(UShort_t)*fTRACKbuffer); // probably not necessary becouse already done ?
-fTracklength=0;
-
-  for (Int_t iv=0; iv<fNviews;iv++) {
-  memset(DataDSP,0,sizeof(UShort_t)*DSPbufferSize);
-  DSPlength=16; // skip the header, to be written later
-  UShort_t CheckSum=0;
-// write dsp header on buffer
-
-//    fDataTrack[fTracklength]=0xE805;
-//    fTracklength++;
-
-//    fDataTrack[fTracklength]=0x01A9;
-//    fTracklength++;
-
-// end dsp header
-
-   //
-   // INIZIO VISTA IV - TAKE PROPER ACTION
-   //
-
-
-
-    for (Int_t ladder=0; ladder<fNladder;ladder++) {
-      Int_t k=0;
-      while (k<fNstrips_ladder) {
-        // compress write in buffer the current LADDER
-        if ( k == 0)  {
-	  real=modff(AdcTrack[iv][ladder*fNstrips_ladder+k],&inte);
-          if (real > 0.5) inte=inte+1;
-          newval=(Int_t)inte -(Int_t)fPedeTrack[iv][ladder*fNstrips_ladder+k];
-          // first strip of ladder is transmitted
-          // DC_TOT first " << AdcTrack[iv][ladder*fNstrips_ladder+k] << endl;
-	  DataDSP[DSPlength]=( ((UShort_t)inte) & 0x0FFF);
-	  DSPlength++;  
-          ntrastot++;
-          trasmesso=1;
-          oldval=newval;
-	  kt=k;
-	  k++;
-	  continue; 
-        };
-	real=modff(AdcTrack[iv][ladder*fNstrips_ladder+k],&inte);
-        if (real > 0.5) inte=inte+1;
-        newval=(Int_t)inte -(Int_t)(fPedeTrack[iv][ladder*fNstrips_ladder+k]);
-        cercacluster=1; // ????????? 
-        if (cercacluster==1) {
-	
- // controlla l'ordine di tutti queste strip ladder e DSP !!!!!!!	
-          Int_t diff=0;
-
-	  
-	  switch ((iv+1)%2) {
-	  case 0: diff=newval-oldval;
-	  break;
-	  case 1: diff=oldval-newval;	  
-          break;
-	  };
-
-	  if (diff>fCutclu*(Int_t)fSigmaTrack[iv][ladder*fNstrips_ladder+k]) {
-            Int_t clval=newval;
-            Int_t klp=k;	// go on to search for maximum 
-            klp++;
-
-            while(klp<fNstrips_ladder) {
-	      real=modff(AdcTrack[iv][ladder*fNstrips_ladder+klp],&inte);
-              if (real > 0.5) inte=inte+1;
-              Int_t clvalp=(Int_t)inte -(Int_t)fPedeTrack[iv][ladder*fNstrips_ladder+klp];
-	      if((iv+1)%2==0) {
-
- 	        if(clvalp>clval) {
-	           clval=clvalp;
-	  	   k=klp;}
-		else break; // max of cluster found
-
-	      }	else {  
-
-	        if(clvalp<clval) {
-	           clval=clvalp;
-		   k=klp;}
-		else break; // max of cluster found
-
-	      };
-	      
-	      klp++;
-	    };
-
-            Int_t kl1=k-fNclst; // max of cluster (or end of ladder ?)
-	    trasmesso=0;
-	    if(kl1<0)  kl1=0;
-
-	    if(kt>=kl1) kl1=kt+1;
-            if( (kt+1)==kl1 ) trasmesso=1;
-
-
-	    
-	    Int_t kl2=k+fNclst;
-	    if(kl2>=fNstrips_ladder) kl2=fNstrips_ladder-1;
-
-	    for(Int_t klt=kl1 ; klt<=kl2 ; klt++) {
-              if(trasmesso==0) { 
-	      //  std:: cout << "STRIP " << klt << endl;
-              //  std:: cout << "ADC_TOT " <<AdcTrack[iv][ladder*fNstrips_ladder+klt] << endl;
-
-	        DataDSP[DSPlength]=( ((UShort_t)klt) | 0x1000);
-	        DSPlength++;  
-                ntrastot++;
-	     
-
-	        real=modff(AdcTrack[iv][ladder*fNstrips_ladder+klt],&inte);
-                if (real > 0.5) inte=inte+1;
-	        DataDSP[DSPlength]=( ((UShort_t)inte) & 0x0FFF);
-	        DSPlength++; 
-	        ntrastot++;
-	     
-               }
-	       else { 
-               //   std:: cout << "ADC_TOT " <<AdcTrack[iv][ladder*fNstrips_ladder+klt] << endl;
- 	        real=modff(AdcTrack[iv][ladder*fNstrips_ladder+klt],&inte);
-                if (real > 0.5) inte=inte+1;
-	        DataDSP[DSPlength]=( ((UShort_t)inte) & 0x0FFF);
-	        DSPlength++;               
-		ntrastot++;
-                };
-                trasmesso=1;	      	   	    
-	    };  // end trasmission
-	    kt=kl2;
-	    k=kl2;
-	    real=modff(AdcTrack[iv][ladder*fNstrips_ladder+kt],&inte);
-            if (real > 0.5) inte=inte+1;
- 	    oldval=(Int_t)inte -(Int_t)fPedeTrack[iv][ladder*fNstrips_ladder+kt];
-            k++;
-	    continue;
-
-
-	  }; // end cercacluster
-	}; // end cercacluster
-      
-// start ZOP check for strips no
-      
-      if(abs(newval-oldval)>=fCutzop*(Int_t)fSigmaTrack[iv][ladder*fNstrips_ladder+k]) { 
-
-       if(trasmesso==0) {
-       // std:: cout << "STRIP " << k << endl;
-       // std:: cout << "ADC_TOT " << AdcTrack[iv][ladder*fNstrips_ladder+k] << endl;
-
-	 DataDSP[DSPlength]=( ((UShort_t)k) | 0x1000);
-	 DSPlength++;  
-         ntrastot++;
-	     
-
-	 real=modff(AdcTrack[iv][ladder*fNstrips_ladder+k],&inte);
-         if (real > 0.5) inte=inte+1;
-	 DataDSP[DSPlength]=( ((UShort_t)inte) & 0x0FFF);
-	 DSPlength++; 
-	 ntrastot++;
-
-       }
-       else {
-       //  std:: cout << "ADC_TOT " << AdcTrack[iv][ladder*fNstrips_ladder+k] << endl;
-	 real=modff(AdcTrack[iv][ladder*fNstrips_ladder+k],&inte);
-         if (real > 0.5) inte=inte+1;
-	 DataDSP[DSPlength]=(  ((UShort_t)inte) & 0x0FFF);
-	 DSPlength++; 
-	 ntrastot++;      
-       };
-       trasmesso=1;
-       oldval=newval;
-       kt=k;
-
-      }  
-      else trasmesso=0;
-      // end zop 
-      
-      k++;      
-      };  // end cycle inside ladder
-// write here the end ladder bytes
-//            std:: cout << "FINE LADDER " << ladder+1 << endl; 
-
-      DataDSP[DSPlength]=( ((UShort_t)(ladder+1))  | 0x1800);
-      DSPlength++; 
-      ntrastot++; 
-      trasmesso=0;
-
-    };  //end cycle inside dsp
-//  std:: cout << "FINE DSP " << iv+1 << endl;
-// here put DSP header
-    DataDSP[0]=(0x1CA0 | ((UShort_t)(iv+1)) );
-    UShort_t Nword=(DSPlength*13)/16; 
-    if( ((DSPlength*13)%16)!=0) Nword++;
-    DataDSP[1]=(0x1400 | ( Nword >> 10));
-    DataDSP[2]=(0x1400 | ( Nword & 0x03FF) );
-    DataDSP[3]=(0x1400 | (( (UShort_t)(fCounter >> 10) ) & 0x03FF) ); 
-    DataDSP[4]=(0x1400 | (( (UShort_t)(fCounter) )  & 0x03FF) );
-    DataDSP[5]=(0x1400 | ( (UShort_t)(fNclst << 7) ) | ( (UShort_t)(fCutzop << 4) )
-     | ( (UShort_t)fCutzop  ) ); 
-    DataDSP[6]=0x1400;
-    DataDSP[7]=0x1400;
-    DataDSP[8]=0x1400;
-    DataDSP[9]=0x1400;
-    DataDSP[10]=0x1400;
-    DataDSP[11]=0x1400;
-    DataDSP[12]=0x1400;
-    DataDSP[13]=0x1400;
-    DataDSP[14]=(0x1400 | (CheckSum & 0x00FF) ); 
-    DataDSP[15]=0x1C00;
-// end DSP header    
-
-
-// write 13 bit DataDSP bufer inside 16 bit fDataTrack buffer
-    Int_t Bit16free=16;
-    UShort_t Dato;
-    for (Int_t NDSP=0; NDSP<DSPlength;NDSP++) {
-      Int_t Bit13ToWrite=13;
-      while(Bit13ToWrite>0) { 
-        if(Bit13ToWrite<=Bit16free) {
-          Dato=((DataDSP[NDSP]&(0xFFFF >> (16-Bit13ToWrite)))<<(Bit16free-Bit13ToWrite));
-	  fDataTrack[fTracklength]=fDataTrack[fTracklength] | Dato ;
-	  Bit16free=Bit16free-Bit13ToWrite;
-	  Bit13ToWrite=0;
-	  if(Bit16free==0) {
-	    if(NDSP>15) CheckSum=CheckSum^fDataTrack[fTracklength];
-	    fTracklength++;
-	    Bit16free=16;
-	  };          
-        }
-	else if(Bit13ToWrite>Bit16free) {
-	  Dato=( (DataDSP[NDSP]&(0xFFFF >> (16-Bit13ToWrite) ) ) >> (Bit13ToWrite-Bit16free) );
-	  fDataTrack[fTracklength]=fDataTrack[fTracklength] | Dato ;
-	  if(NDSP>15) CheckSum=CheckSum^fDataTrack[fTracklength];
-	  fTracklength++;
-	  Bit13ToWrite=Bit13ToWrite-Bit16free;
-	  Bit16free=16;	 
-	};
-	
-      }; // end cycle while(Bit13ToWrite>0)
-          
-    }; // end cycle DataDSP
-    if(Bit16free!=16) { fTracklength++; CheckSum=CheckSum^fDataTrack[fTracklength]; }; 
-    CheckSum=(CheckSum >> 8)^(CheckSum&0x00FF);
-    fDataTrack[fTracklength-Nword+11]=(0x0280 | (CheckSum >> 3));
-    fDataTrack[fTracklength-Nword+12]=(0x1C00 | (CheckSum << 13) );
-
-// end write 13 bit DataDSP bufer inside 16 bit fDataTrack buffer
-
-//write trailer on buffer
-    UShort_t ReLength=(UShort_t)((Nword+13)*2+3);
-    UShort_t OveCheckCode=0x0000;
-
-    fDataTrack[fTracklength]=0x0000;
-    fTracklength++;
-  
-    fDataTrack[fTracklength]=(ReLength >> 8);
-    fTracklength++;
-  
-    fDataTrack[fTracklength]=( (ReLength << 8) | (OveCheckCode & 0x00FF) );
-    fTracklength++;  
-// end trailer 
-//    std:: cout  << "DSPlength  " <<DSPlength << endl;
-//    std:: cout << "Nword " << Nword  << endl;
-//    std:: cout <<  "ReLength " << ReLength << endl;
-  };    
-//    std:: cout << "ntrastot " << ntrastot << endl;    
-
-};
-
-
-Float_t Digitizer::SaturationTrack(Float_t ADC) {
-  Float_t SatFact=1.;
-  if(ADC<70.) { SatFact=80./ADC; }; 
-  if(ADC>3000.) { SatFact=3000./ADC; }; 
-  return SatFact;
-};
-
-
-
-
-
-
+//               ------ PAMELA Digitizer ------
+// 
+// Date, release and how-to: see file  Pamelagp2Digits.cxx 
+// 
+// NB: Check length  physics  packet [packet type (0x10 = physics data)]
+//
+#include "Digitizer.h"
+
+extern "C"{
+  short crc(short, short);
+};
+//
+
+Digitizer::Digitizer(TTree* tree, char* &file_raw,
+		     int nspe1=200,int ntof1=200,int ncat1=50,
+		     int ncas1=50,int ncar1=100,int ncal1=1000,
+		     int nnd1=200,int nstr1=1000, int comprcalomod=0){
+  
+  nspe=new int[1];
+  ntof=new int[1];
+  ncat=new int[1];
+  ncas=new int[1];
+  ncar=new int[1];
+  ncal=new int[1];
+  nnd=new int[1];
+  nstr=new int[1];
+
+  *nspe=nspe1;
+  *ntof=ntof1;
+  *ncat=ncat1;
+  *ncas=ncas1;
+  *ncar=ncar1;
+  *ncal=ncal1;
+  *nnd=nnd1;
+  *nstr=nstr1;
+
+  fhBookTree = tree;
+  fFilename =  file_raw;
+  fCounter = 0;
+  fCounterPhys = 0; // SO 5/12/'07
+  fOBT = 0;
+  fModCalo = comprcalomod ;
+
+  //
+  // DB  connections
+  //
+  TString host = "mysql://localhost/pamelaprod";
+  TString user = "anonymous";
+  TString psw = "";
+  //
+  const char *pamdbhost=gSystem->Getenv("PAM_DBHOST");
+  const char *pamdbuser=gSystem->Getenv("PAM_DBUSER");
+  const char *pamdbpsw=gSystem->Getenv("PAM_DBPSW");
+  if ( !pamdbhost ) pamdbhost = "";
+  if ( !pamdbuser ) pamdbuser = "";
+  if ( !pamdbpsw ) pamdbpsw = "";
+  if ( strcmp(pamdbhost,"") ) host = pamdbhost;
+  if ( strcmp(pamdbuser,"") ) user = pamdbuser;
+  if ( strcmp(pamdbpsw,"") ) psw = pamdbpsw;
+  fDbc = TSQLServer::Connect(host.Data(),user.Data(),psw.Data());
+  //
+  GL_TABLES *glt = new GL_TABLES(host,user,psw);
+  if ( glt->IsConnected(fDbc) ) printf("\n DB INFORMATION:\n SQL: %s Version: %s Host %s Port %i \n\n",fDbc->GetDBMS(),fDbc->ServerInfo(),fDbc->GetHost(),fDbc->GetPort());
+  //
+  // Use UTC in the DB and make timeout bigger
+  //
+  stringstream myquery;
+  myquery.str("");
+  myquery << "SET time_zone='+0:00'";
+  fDbc->Query(myquery.str().c_str());
+  myquery.str("");
+  myquery << "SET wait_timeout=173000;";
+  fDbc->Query(myquery.str().c_str());
+  //
+  
+  std:: cout << "preparing tree" << endl;
+
+  Ipltof=(UChar_t*)malloc(*ntof *sizeof(UChar_t));
+  Ipaddle=(UChar_t*)malloc(*ntof *sizeof(UChar_t));
+  Ipartof=(UShort_t*)malloc(*ntof *sizeof(UShort_t));
+  //  Ipartof=(UChar_t*)malloc(*ntof *sizeof(UChar_t));//DPMJET
+  Xintof=(Float_t*)malloc(*ntof *sizeof(Float_t));
+  Yintof=(Float_t*)malloc(*ntof *sizeof(Float_t));
+  Zintof=(Float_t*)malloc(*ntof *sizeof(Float_t));
+  Xouttof=(Float_t*)malloc(*ntof *sizeof(Float_t));
+  Youttof=(Float_t*)malloc(*ntof *sizeof(Float_t));
+  Zouttof=(Float_t*)malloc(*ntof *sizeof(Float_t));
+  Ereltof=(Float_t*)malloc(*ntof *sizeof(Float_t));
+  Timetof=(Float_t*)malloc(*ntof *sizeof(Float_t));
+  Pathtof=(Float_t*)malloc(*ntof *sizeof(Float_t));
+  P0tof=(Float_t*)malloc(*ntof *sizeof(Float_t));
+  Iparcat=(UChar_t*)malloc(*ncat *sizeof(UChar_t));
+  Icat=(UChar_t*)malloc(*ncat *sizeof(UChar_t));
+  Xincat=(Float_t*)malloc(*ncat *sizeof(Float_t));
+  Yincat=(Float_t*)malloc(*ncat *sizeof(Float_t));
+  Zincat=(Float_t*)malloc(*ncat *sizeof(Float_t));
+  Xoutcat=(Float_t*)malloc(*ncat *sizeof(Float_t));
+  Youtcat=(Float_t*)malloc(*ncat *sizeof(Float_t));
+  Zoutcat=(Float_t*)malloc(*ncat *sizeof(Float_t));
+  Erelcat=(Float_t*)malloc(*ncat *sizeof(Float_t));
+  Timecat=(Float_t*)malloc(*ncat *sizeof(Float_t));
+  Pathcat=(Float_t*)malloc(*ncat *sizeof(Float_t));
+  P0cat=(Float_t*)malloc(*ncat *sizeof(Float_t));
+  Iparcas=(UChar_t*)malloc(*ncas *sizeof(UChar_t));
+  Icas=(UChar_t*)malloc(*ncas *sizeof(UChar_t));
+  Xincas=(Float_t*)malloc(*ncas *sizeof(Float_t));
+  Yincas=(Float_t*)malloc(*ncas *sizeof(Float_t));
+  Zincas=(Float_t*)malloc(*ncas *sizeof(Float_t));
+  Xoutcas=(Float_t*)malloc(*ncas *sizeof(Float_t));
+  Youtcas=(Float_t*)malloc(*ncas *sizeof(Float_t));
+  Zoutcas=(Float_t*)malloc(*ncas *sizeof(Float_t));
+  Erelcas=(Float_t*)malloc(*ncas *sizeof(Float_t));
+  Timecas=(Float_t*)malloc(*ncas *sizeof(Float_t));
+  Pathcas=(Float_t*)malloc(*ncas *sizeof(Float_t));
+  P0cas=(Float_t*)malloc(*ncas *sizeof(Float_t));
+  //  Iparspe=(UShort_t*)malloc(*nspe *sizeof(UShort_t));
+  //  Iparspe=(UChar_t*)malloc(*nspe *sizeof(UChar_t));
+  Itrpb=(UChar_t*)malloc(*nspe *sizeof(UChar_t));
+  Itrsl=(UChar_t*)malloc(*nspe *sizeof(UChar_t));
+  Itspa=(UChar_t*)malloc(*nspe *sizeof(UChar_t));
+  Xinspe=(Float_t*)malloc(*nspe *sizeof(Float_t));
+  Yinspe=(Float_t*)malloc(*nspe *sizeof(Float_t));
+  Zinspe=(Float_t*)malloc(*nspe *sizeof(Float_t));
+  Xoutspe=(Float_t*)malloc(*nspe *sizeof(Float_t));
+  Youtspe=(Float_t*)malloc(*nspe *sizeof(Float_t));
+  Zoutspe=(Float_t*)malloc(*nspe *sizeof(Float_t));
+  Xavspe=(Float_t*)malloc(*nspe *sizeof(Float_t));
+  Yavspe=(Float_t*)malloc(*nspe *sizeof(Float_t));
+  Zavspe=(Float_t*)malloc(*nspe *sizeof(Float_t));
+  Erelspe=(Float_t*)malloc(*nspe *sizeof(Float_t));
+  Pathspe=(Float_t*)malloc(*nspe *sizeof(Float_t));
+  P0spe=(Float_t*)malloc(*nspe *sizeof(Float_t));;
+  Nxmult=(UChar_t*)malloc(*nspe *sizeof(UChar_t));
+  Nymult=(UChar_t*)malloc(*nspe *sizeof(UChar_t));
+  Istripx=(UShort_t*)malloc(*nstr *sizeof(UShort_t));
+  Qstripx=(Float_t*)malloc(*nstr *sizeof(Float_t));
+  Xstripx=(Float_t*)malloc(*nstr *sizeof(Float_t));
+  Npstripx=(UChar_t*)malloc(*nstr *sizeof(UChar_t));
+  Ntstripx=(UChar_t*)malloc(*nstr *sizeof(UChar_t));
+  Npstripy=(UChar_t*)malloc(*nstr *sizeof(UChar_t));
+  Ntstripy=(UChar_t*)malloc(*nstr *sizeof(UChar_t));
+  Istripy=(UShort_t*)malloc(*nstr *sizeof(UShort_t));
+  Qstripy=(Float_t*)malloc(*nstr *sizeof(Float_t));
+  Ystripy=(Float_t*)malloc(*nstr *sizeof(Float_t));
+  Icapl=(UChar_t*)malloc(*ncal *sizeof(UChar_t));
+  Icasi=(UChar_t*)malloc(*ncal *sizeof(UChar_t));
+  Icast=(UChar_t*)malloc(*ncal *sizeof(UChar_t));
+  Xincal=(Float_t*)malloc(*ncal *sizeof(Float_t));
+  Yincal=(Float_t*)malloc(*ncal *sizeof(Float_t));
+  Zincal=(Float_t*)malloc(*ncal *sizeof(Float_t));
+  Erelcal=(Float_t*)malloc(*ncal *sizeof(Float_t));
+  Itubend=(UChar_t*)malloc(*nnd *sizeof(UChar_t));
+  Iparnd=(UChar_t*)malloc(*nnd *sizeof(UChar_t));
+  Xinnd=(Float_t*)malloc(*nnd *sizeof(Float_t));
+  Yinnd=(Float_t*)malloc(*nnd *sizeof(Float_t));
+  Zinnd=(Float_t*)malloc(*nnd *sizeof(Float_t));
+  Xoutnd=(Float_t*)malloc(*nnd *sizeof(Float_t));
+  Youtnd=(Float_t*)malloc(*nnd *sizeof(Float_t));
+  Zoutnd=(Float_t*)malloc(*nnd *sizeof(Float_t));
+  Erelnd=(Float_t*)malloc(*nnd *sizeof(Float_t));
+  Timend=(Float_t*)malloc(*nnd *sizeof(Float_t));
+  Pathnd=(Float_t*)malloc(*nnd *sizeof(Float_t));
+  P0nd=(Float_t*)malloc(*nnd *sizeof(Float_t));
+  Iparcard=(UChar_t*)malloc(*ncar *sizeof(UChar_t));
+  Icard=(UChar_t*)malloc(*ncar *sizeof(UChar_t));
+  Xincard=(Float_t*)malloc(*ncar *sizeof(Float_t));
+  Yincard=(Float_t*)malloc(*ncar *sizeof(Float_t));
+  Zincard=(Float_t*)malloc(*ncar *sizeof(Float_t));
+  Xoutcard=(Float_t*)malloc(*ncar *sizeof(Float_t));
+  Youtcard=(Float_t*)malloc(*ncar *sizeof(Float_t));
+  Zoutcard=(Float_t*)malloc(*ncar *sizeof(Float_t));
+  Erelcard=(Float_t*)malloc(*ncar *sizeof(Float_t));
+  Timecard=(Float_t*)malloc(*ncar *sizeof(Float_t));
+  Pathcard=(Float_t*)malloc(*ncar *sizeof(Float_t));
+  P0card=(Float_t*)malloc(*ncar *sizeof(Float_t));
+
+
+
+  // prepare tree//modified by E.Vannuccini 03/08
+  if(fhBookTree->GetBranch("Irun"))fhBookTree->SetBranchAddress("Irun",&Irun);
+  if(fhBookTree->GetBranch("Ievnt"))fhBookTree->SetBranchAddress("Ievnt",&Ievnt);
+  if(fhBookTree->GetBranch("Ipa"))fhBookTree->SetBranchAddress("Ipa",&Ipa);
+  if(fhBookTree->GetBranch("X0"))fhBookTree->SetBranchAddress("X0",&X0);
+  if(fhBookTree->GetBranch("Y0"))fhBookTree->SetBranchAddress("Y0",&Y0);
+  if(fhBookTree->GetBranch("Z0"))fhBookTree->SetBranchAddress("Z0",&Z0);
+  if(fhBookTree->GetBranch("Theta"))fhBookTree->SetBranchAddress("Theta",&Theta);
+  if(fhBookTree->GetBranch("Phi"))fhBookTree->SetBranchAddress("Phi",&Phi);
+  if(fhBookTree->GetBranch("P0"))fhBookTree->SetBranchAddress("P0",&P0);
+  if(fhBookTree->GetBranch("Nthtof"))fhBookTree->SetBranchAddress("Nthtof",&Nthtof);
+  if(fhBookTree->GetBranch("Ipltof"))fhBookTree->SetBranchAddress("Ipltof",Ipltof);///////////////////////////
+  if(fhBookTree->GetBranch("Ipaddle"))fhBookTree->SetBranchAddress("Ipaddle",Ipaddle);
+  if(fhBookTree->GetBranch("Ipartof"))fhBookTree->SetBranchAddress("Ipartof",Ipartof);
+  if(fhBookTree->GetBranch("Xintof"))fhBookTree->SetBranchAddress("Xintof",Xintof);
+  if(fhBookTree->GetBranch("Yintof"))fhBookTree->SetBranchAddress("Yintof",Yintof);
+  if(fhBookTree->GetBranch("Zintof"))fhBookTree->SetBranchAddress("Zintof",Zintof);
+  if(fhBookTree->GetBranch("Xouttof"))fhBookTree->SetBranchAddress("Xouttof",Xouttof);
+  if(fhBookTree->GetBranch("Youttof"))fhBookTree->SetBranchAddress("Youttof",Youttof);
+  if(fhBookTree->GetBranch("Zouttof"))fhBookTree->SetBranchAddress("Zouttof",Zouttof);
+  if(fhBookTree->GetBranch("Ereltof"))fhBookTree->SetBranchAddress("Ereltof",Ereltof);
+  if(fhBookTree->GetBranch("Timetof"))fhBookTree->SetBranchAddress("Timetof",Timetof);
+  if(fhBookTree->GetBranch("Pathtof"))fhBookTree->SetBranchAddress("Pathtof",Pathtof);
+  if(fhBookTree->GetBranch("P0tof"))fhBookTree->SetBranchAddress("P0tof",P0tof);
+  if(fhBookTree->GetBranch("Nthcat"))fhBookTree->SetBranchAddress("Nthcat",&Nthcat);
+  if(fhBookTree->GetBranch("Iparcat"))fhBookTree->SetBranchAddress("Iparcat",Iparcat);
+  if(fhBookTree->GetBranch("Icat"))fhBookTree->SetBranchAddress("Icat",Icat);
+  if(fhBookTree->GetBranch("Xincat"))fhBookTree->SetBranchAddress("Xincat",Xincat);
+  if(fhBookTree->GetBranch("Yincat"))fhBookTree->SetBranchAddress("Yincat",Yincat);
+  if(fhBookTree->GetBranch("Zincat"))fhBookTree->SetBranchAddress("Zincat",Zincat);
+  if(fhBookTree->GetBranch("Xoutcat"))fhBookTree->SetBranchAddress("Xoutcat",Xoutcat);
+  if(fhBookTree->GetBranch("Youtcat"))fhBookTree->SetBranchAddress("Youtcat",Youtcat);
+  if(fhBookTree->GetBranch("Zoutcat"))fhBookTree->SetBranchAddress("Zoutcat",Zoutcat);
+  if(fhBookTree->GetBranch("Erelcat"))fhBookTree->SetBranchAddress("Erelcat",Erelcat);
+  if(fhBookTree->GetBranch("Timecat"))fhBookTree->SetBranchAddress("Timecat",Timecat);
+  if(fhBookTree->GetBranch("Pathcat"))fhBookTree->SetBranchAddress("Pathcat",Pathcat);
+  if(fhBookTree->GetBranch("P0cat"))fhBookTree->SetBranchAddress("P0cat",P0cat);
+  if(fhBookTree->GetBranch("Nthcas"))fhBookTree->SetBranchAddress("Nthcas",&Nthcas);
+  if(fhBookTree->GetBranch("Iparcas"))fhBookTree->SetBranchAddress("Iparcas",Iparcas);
+  if(fhBookTree->GetBranch("Icas"))fhBookTree->SetBranchAddress("Icas",Icas);///////////////////////////////
+  if(fhBookTree->GetBranch("Xincas"))fhBookTree->SetBranchAddress("Xincas",Xincas);
+  if(fhBookTree->GetBranch("Yincas"))fhBookTree->SetBranchAddress("Yincas",Yincas);
+  if(fhBookTree->GetBranch("Zincas"))fhBookTree->SetBranchAddress("Zincas",Zincas);
+  if(fhBookTree->GetBranch("Xoutcas"))fhBookTree->SetBranchAddress("Xoutcas",Xoutcas);
+  if(fhBookTree->GetBranch("Youtcas"))fhBookTree->SetBranchAddress("Youtcas",Youtcas);
+  if(fhBookTree->GetBranch("Zoutcas"))fhBookTree->SetBranchAddress("Zoutcas",Zoutcas);
+  if(fhBookTree->GetBranch("Erelcas"))fhBookTree->SetBranchAddress("Erelcas",Erelcas);
+  if(fhBookTree->GetBranch("Timecas"))fhBookTree->SetBranchAddress("Timecas",Timecas);
+  if(fhBookTree->GetBranch("Pathcas"))fhBookTree->SetBranchAddress("Pathcas",Pathcas);
+  if(fhBookTree->GetBranch("P0cas"))fhBookTree->SetBranchAddress("P0cas",P0cas);
+  if(fhBookTree->GetBranch("Nthspe"))fhBookTree->SetBranchAddress("Nthspe",&Nthspe);
+  //  if(fhBookTree->GetBranch("Iparspe"))fhBookTree->SetBranchAddress("Iparspe",Iparspe);
+  if(fhBookTree->GetBranch("Itrpb"))fhBookTree->SetBranchAddress("Itrpb",Itrpb);
+  if(fhBookTree->GetBranch("Itrsl"))fhBookTree->SetBranchAddress("Itrsl",Itrsl);
+  if(fhBookTree->GetBranch("Itspa"))fhBookTree->SetBranchAddress("Itspa",Itspa);
+  if(fhBookTree->GetBranch("Xinspe"))fhBookTree->SetBranchAddress("Xinspe",Xinspe);
+  if(fhBookTree->GetBranch("Yinspe"))fhBookTree->SetBranchAddress("Yinspe",Yinspe);
+  if(fhBookTree->GetBranch("Zinspe"))fhBookTree->SetBranchAddress("Zinspe",Zinspe);
+  if(fhBookTree->GetBranch("Xoutspe"))fhBookTree->SetBranchAddress("Xoutspe",Xoutspe);
+  if(fhBookTree->GetBranch("Youtspe"))fhBookTree->SetBranchAddress("Youtspe",Youtspe);
+  if(fhBookTree->GetBranch("Zoutspe"))fhBookTree->SetBranchAddress("Zoutspe",Zoutspe);
+  if(fhBookTree->GetBranch("Xavspe"))fhBookTree->SetBranchAddress("Xavspe",Xavspe);
+  if(fhBookTree->GetBranch("Yavspe"))fhBookTree->SetBranchAddress("Yavspe",Yavspe);
+  if(fhBookTree->GetBranch("Zavspe"))fhBookTree->SetBranchAddress("Zavspe",Zavspe);
+  if(fhBookTree->GetBranch("Erelspe"))fhBookTree->SetBranchAddress("Erelspe",Erelspe);
+  if(fhBookTree->GetBranch("Pathspe"))fhBookTree->SetBranchAddress("Pathspe",Pathspe);
+  if(fhBookTree->GetBranch("P0spe"))fhBookTree->SetBranchAddress("P0spe",P0spe);
+  if(fhBookTree->GetBranch("Nxmult"))fhBookTree->SetBranchAddress("Nxmult",Nxmult);
+  if(fhBookTree->GetBranch("Nymult"))fhBookTree->SetBranchAddress("Nymult",Nymult);
+  if(fhBookTree->GetBranch("Nstrpx"))fhBookTree->SetBranchAddress("Nstrpx",&Nstrpx);
+  if(fhBookTree->GetBranch("Npstripx"))fhBookTree->SetBranchAddress("Npstripx",Npstripx);
+  if(fhBookTree->GetBranch("Ntstripx"))fhBookTree->SetBranchAddress("Ntstripx",Ntstripx);
+  if(fhBookTree->GetBranch("Istripx"))fhBookTree->SetBranchAddress("Istripx",Istripx);
+  if(fhBookTree->GetBranch("Qstripx"))fhBookTree->SetBranchAddress("Qstripx",Qstripx);
+  if(fhBookTree->GetBranch("Xstripx"))fhBookTree->SetBranchAddress("Xstripx",Xstripx);
+  if(fhBookTree->GetBranch("Nstrpy"))fhBookTree->SetBranchAddress("Nstrpy",&Nstrpy);
+  if(fhBookTree->GetBranch("Npstripy"))fhBookTree->SetBranchAddress("Npstripy",Npstripy);
+  if(fhBookTree->GetBranch("Ntstripy"))fhBookTree->SetBranchAddress("Ntstripy",Ntstripy);
+  if(fhBookTree->GetBranch("Istripy"))fhBookTree->SetBranchAddress("Istripy",Istripy);
+  if(fhBookTree->GetBranch("Qstripy"))fhBookTree->SetBranchAddress("Qstripy",Qstripy);///////////////////////
+  if(fhBookTree->GetBranch("Ystripy"))fhBookTree->SetBranchAddress("Ystripy",Ystripy);
+  if(fhBookTree->GetBranch("Nthcali"))fhBookTree->SetBranchAddress("Nthcali",&Nthcali);
+  if(fhBookTree->GetBranch("Icaplane"))fhBookTree->SetBranchAddress("Icaplane",Icaplane);
+  if(fhBookTree->GetBranch("Icastrip"))fhBookTree->SetBranchAddress("Icastrip",Icastrip);
+  if(fhBookTree->GetBranch("Icamod"))fhBookTree->SetBranchAddress("Icamod",Icamod);
+  if(fhBookTree->GetBranch("Enestrip"))fhBookTree->SetBranchAddress("Enestrip",Enestrip);
+  if(fhBookTree->GetBranch("Nthcal"))fhBookTree->SetBranchAddress("Nthcal",&Nthcal);
+  if(fhBookTree->GetBranch("Icapl"))fhBookTree->SetBranchAddress("Icapl",Icapl);
+  if(fhBookTree->GetBranch("Icasi"))fhBookTree->SetBranchAddress("Icasi",Icasi);
+  if(fhBookTree->GetBranch("Icast"))fhBookTree->SetBranchAddress("Icast",Icast);
+  if(fhBookTree->GetBranch("Xincal"))fhBookTree->SetBranchAddress("Xincal",Xincal);
+  if(fhBookTree->GetBranch("Yincal"))fhBookTree->SetBranchAddress("Yincal",Yincal);
+  if(fhBookTree->GetBranch("Zincal"))fhBookTree->SetBranchAddress("Zincal",Zincal);
+  if(fhBookTree->GetBranch("Erelcal"))fhBookTree->SetBranchAddress("Erelcal",Erelcal);
+  if(fhBookTree->GetBranch("Nthnd"))fhBookTree->SetBranchAddress("Nthnd",&Nthnd);
+  if(fhBookTree->GetBranch("Itubend"))fhBookTree->SetBranchAddress("Itubend",Itubend);
+  if(fhBookTree->GetBranch("Iparnd"))fhBookTree->SetBranchAddress("Iparnd",Iparnd);
+  if(fhBookTree->GetBranch("Xinnd"))fhBookTree->SetBranchAddress("Xinnd",Xinnd);/////////////////////////
+  if(fhBookTree->GetBranch("Yinnd"))fhBookTree->SetBranchAddress("Yinnd",Yinnd);
+  if(fhBookTree->GetBranch("Zinnd"))fhBookTree->SetBranchAddress("Zinnd",Zinnd);
+  if(fhBookTree->GetBranch("Xoutnd"))fhBookTree->SetBranchAddress("Xoutnd",Xoutnd);
+  if(fhBookTree->GetBranch("Youtnd"))fhBookTree->SetBranchAddress("Youtnd",Youtnd);
+  if(fhBookTree->GetBranch("Zoutnd"))fhBookTree->SetBranchAddress("Zoutnd",Zoutnd);
+  if(fhBookTree->GetBranch("Erelnd"))fhBookTree->SetBranchAddress("Erelnd",Erelnd);
+  if(fhBookTree->GetBranch("Timend"))fhBookTree->SetBranchAddress("Timend",Timend);
+  if(fhBookTree->GetBranch("Pathnd"))fhBookTree->SetBranchAddress("Pathnd",Pathnd);
+  if(fhBookTree->GetBranch("P0nd"))fhBookTree->SetBranchAddress("P0nd",P0nd);
+  if(fhBookTree->GetBranch("Nthcard"))fhBookTree->SetBranchAddress("Nthcard",&Nthcard);/////////////////////
+  if(fhBookTree->GetBranch("Iparcard"))fhBookTree->SetBranchAddress("Iparcard",Iparcard);
+  if(fhBookTree->GetBranch("Icard"))fhBookTree->SetBranchAddress("Icard",Icard);
+  if(fhBookTree->GetBranch("Xincard"))fhBookTree->SetBranchAddress("Xincard",Xincard);
+  if(fhBookTree->GetBranch("Yincard"))fhBookTree->SetBranchAddress("Yincard",Yincard);
+  if(fhBookTree->GetBranch("Zincard"))fhBookTree->SetBranchAddress("Zincard",Zincard);
+  if(fhBookTree->GetBranch("Xoutcard"))fhBookTree->SetBranchAddress("Xoutcard",Xoutcard);
+  if(fhBookTree->GetBranch("Youtcard"))fhBookTree->SetBranchAddress("Youtcard",Youtcard);/////////////////
+  if(fhBookTree->GetBranch("Zoutcard"))fhBookTree->SetBranchAddress("Zoutcard",Zoutcard);
+  if(fhBookTree->GetBranch("Erelcard"))fhBookTree->SetBranchAddress("Erelcard",Erelcard);
+  if(fhBookTree->GetBranch("Timecard"))fhBookTree->SetBranchAddress("Timecard",Timecard);
+  if(fhBookTree->GetBranch("Pathcard"))fhBookTree->SetBranchAddress("Pathcard",Pathcard);
+  //  if(fhBookTree->GetBranch("P0card"))fhBookTree->SetBranchAddress("P0card",P0card);
+//    fhBookTree->SetBranchStatus("*",0); //modified by E.Vannuccini 03/08
+}
+
+void Digitizer::Close(){
+  delete fhBookTree;
+}
+
+void Digitizer::Loop() {
+  //
+  // opens the raw output file and loops over the events
+  //
+  fOutputfile.open(fFilename, ios::out | ios::binary);
+  //fOutputfile.open(Form("Output%s",fFilename), ios::out | ios::binary);
+  //
+  // Load in memory and save at the beginning of file the calorimeter calibration
+  //
+  CaloLoadCalib();
+  DigitizeCALOCALIB();
+
+  //  load, digitize and write tracker calibration
+  LoadTrackCalib();
+  
+  DigitizeTrackCalib(1);
+  UInt_t length=fTracklength*2;
+  DigitizePSCU(length,0x12,fDataPSCU);
+  AddPadding();
+  WriteTrackCalib(); 
+  
+  DigitizeTrackCalib(2);
+  length=fTracklength*2;
+  DigitizePSCU(length,0x13,fDataPSCU);
+  AddPadding();
+  WriteTrackCalib();
+  LoadMipCor();  // some initialization of parameters -not used now-
+  //  end loading, digitizing and writing tracker calibration
+  // TOF ------ read calibration file (get A1, A2, lambda1, lambda2)
+  const int np=48;
+  float *atte1,*atte2,*lambda1,*lambda2;
+  atte1=(float *)malloc(np *sizeof(float));
+  atte2=(float *)malloc(np *sizeof(float));
+  lambda1=(float *)malloc(np *sizeof(float));
+  lambda2=(float *)malloc(np *sizeof(float));
+  LoadTOFCalib(np,atte1,atte2,lambda1,lambda2);
+  attenAC = new TF1("fAttAC",".825+.64*atan(9.8/x)",0.,45.);
+  //end tof calib
+  //
+  // loops over the events
+  //
+  
+  Int_t nentries = fhBookTree->GetEntriesFast();
+  Long64_t nbytes = 0;
+  for (Int_t i=0; i<nentries;i++) {
+     nbytes += fhBookTree->GetEntry(i);
+     fEvent=i; //  cecilia for calo compress mode
+      // read detectors sequentially:
+      // http://www.ts.infn.it/fileadmin/documents/physics/experiments/wizard/cpu/gen_arch/RM_Acquisition.pdf
+      // on pamelatov: /cvs/yoda/techmodel/physics/NeutronDetectorReader.cpp
+     DigitizeTOF(np,atte1,atte2,lambda1,lambda2);
+      DigitizeAC();
+      DigitizeCALO();
+      DigitizeTrack();
+      DigitizeS4();
+      DigitizeND();
+      //
+      // Add padding to 64 bits
+      //
+      AddPadding();
+      //
+      // Create CPU header, we need packet type (0x10 = physics data) and packet length.
+      //
+      UInt_t length=2*(fCALOlength+fACbuffer+fTracklength+fNDbuffer+fS4buffer)+fPadding+fTOFbuffer+fTRIGGERbuffer;
+      //UInt_t length=2*(fCALOlength+fACbuffer+fTracklength+fNDbuffer)+fPadding+fTOFbuffer+fTRIGGERbuffer;
+      DigitizePSCU(length,0x10,fDataPSCU);
+      if ((i%1000)==0)cout << "writing event " << i << endl;
+      WriteData();
+  }
+  
+  fOutputfile.close();
+  cout << "files closed" << endl;
+};
+
+
+void Digitizer::ReadData(){
+
+  UShort_t InData[64];
+
+  // for debuggigng purposes only, write your own routine if you like (many
+  // hardwired things.
+  
+  ifstream InputFile;
+
+  // if (!InputFile) {
+
+  //     std::cout << "ERROR" << endl;
+  //     // An error occurred!
+  //     // myFile.gcount() returns the number of bytes read.
+  //     // calling myFile.clear() will reset the stream state
+  //     // so it is usable again.
+  //   };
+
+
+  
+  //InputFile.seekg(0);
+
+  InputFile.open(fFilename, ios::in | ios::binary);
+  //    fOutputfile.seekg(0);
+  if (!InputFile.is_open()) std::cout << "ERROR" << endl;
+
+  InputFile.seekg(0);
+  
+  for (Int_t k=0; k<=1000; k++){
+    InputFile.read(reinterpret_cast<char*>(InData),384*sizeof(UShort_t));
+
+    std::cout << "Read back: " << endl << endl;
+
+    for (Int_t i=0; i<=384; i++){
+      printf("%4x ", InData[i]);   
+      if ((i+1)%8 ==0) cout << endl;
+    }
+
+  }
+  cout << endl;
+  InputFile.close();
+
+};