/[PAMELA software]/PamelaDigitizer/Digitizer.cxx

Diff of /PamelaDigitizer/Digitizer.cxx

Parent Directory | Revision Log | View Patch Patch

-revision 1.2 by orsi,
Fri Sep 28 10:46:23 2007 UTC
+revision 1.4 by orsi,
Wed Oct 31 18:17:58 2007 UTC
 Line 1
  //               ------ 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();
+     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+fS4buffer)+fPadding+fTOFbuffer+fTRIGGERbuffer;
-     UInt_t length=2*(fCALOlength+fACbuffer+fTracklength+fNDbuffer)+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[0]=0x00;
+     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);
+   // --- activate branches:
-   fhBookTree->SetBranchStatus("Ipltof",1);
+   fhBookTree->SetBranchStatus("Nthtof",1);
-   fhBookTree->SetBranchStatus("Ipaddle",1);
+   fhBookTree->SetBranchStatus("Ipltof",1);
-   fhBookTree->SetBranchStatus("Xintof",1);
+   fhBookTree->SetBranchStatus("Ipaddle",1);
-   fhBookTree->SetBranchStatus("Yintof",1);
+   fhBookTree->SetBranchStatus("Xintof",1);
-   fhBookTree->SetBranchStatus("Xouttof",1);
+   fhBookTree->SetBranchStatus("Yintof",1);
-   fhBookTree->SetBranchStatus("Youttof",1);
+   fhBookTree->SetBranchStatus("Xouttof",1);
-   fhBookTree->SetBranchStatus("Ereltof",1);
+   fhBookTree->SetBranchStatus("Youttof",1);
-   fhBookTree->SetBranchStatus("Timetof",1);
+   fhBookTree->SetBranchStatus("Ereltof",1);
-   // not yet used: Zintof, Xouttof, Youttof, Zouttof
+   fhBookTree->SetBranchStatus("Timetof",1);
+   // not yet used: Zintof, Xouttof, Youttof, Zouttof
-   // ------ evaluate energy in each pmt: ------
-   // strip geometry (lenght/width)
+   // ------ evaluate energy in each pmt: ------
-   Float_t dimel[6] = {33.0, 40.8 ,18.0, 15.0, 15.0, 18.0};
+   // strip geometry (lenght/width)
-   //Float_t dimes[6] = {5.1, 5.5, 7.5, 9.0, 6.0, 5.0};
+   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
+   //  S11 8 paddles  33.0 x 5.1 cm
-   //  S21 2 paddles  18.0 x 7.5 cm
+   //  S12 6 paddles  40.8 x 5.5 cm
-   //  S22 2 paddles  15.0 x 9.0 cm
+   //  S21 2 paddles  18.0 x 7.5 cm
-   //  S31 3 paddles  15.0 x 6.0 cm
+   //  S22 2 paddles  15.0 x 9.0 cm
-   //  S32 3 paddles  18.0 x 5.0 cm
+   //  S31 3 paddles  15.0 x 6.0 cm
+   //  S32 3 paddles  18.0 x 5.0 cm
-   // distance from the interaction point to the pmts (right,left)
-   Float_t xpath[2]={0., 0.}; /*path(cm) in X per S12,S21,S32 verso il pmt DX o SX*/
+   Float_t FGeo[2]={0., 0.}; /* geometrical factor */
-   Float_t ypath[2]={0., 0.}; /*path(cm) in Y per S11,S22,S31 verso il pmt DX o SX*/
-   Float_t FGeo[2]={0., 0.}; /* fattore geometrico */
+   const Float_t Pho_keV = 10.;     // photons per keV in scintillator
+   const Float_t echarge = 1.6e-19; // electron charge
-   const Float_t Pho_keV = 10.;     // photons per keV in scintillator
+   Float_t Npho=0.;
-   const Float_t echarge = 1.6e-19; // carica dell'elettrone
+   Float_t QevePmt_pC[48];
-   Float_t Npho=0.;
+   Float_t QhitPad_pC[2]={0., 0.};
-   Float_t QevePmt_pC[48];
+   Float_t QhitPmt_pC[2]={0., 0.};
-   Float_t QhitPad_pC[2]={0., 0.};
+   Float_t pmGain = 3.5e6;  /* PMT Gain: the same for all PMTs */
-   Float_t QhitPmt_pC[2]={0., 0.};
+   Float_t effi=0.21;       /* Efficienza di fotocatodo */
-   Float_t pmGain = 3.5e6;  /* Gain: per il momento uguale per tutti */
-   Float_t effi=0.21;       /* Efficienza di fotocatodo */
+   Float_t ADC_pC0=-58.1;      //  ADC/pC conversion coefficient 0
-   Float_t ADC_pC=1.666667; // ADC_ch/pC conversion = 0.6 pC/channel (+30 di offset)
+   Float_t ADC_pC1=1.728;      //  ADC/pC conversion coefficient 1
-   Float_t ADCoffset=30.;
+   Float_t ADC_pC2=-4.063e-05; //  ADC/pC conversion coefficient 2
-   Int_t ADClast=4095;      // no signal --> ADC ch=4095
+   Float_t ADC_pC3=-5.763e-08; //  ADC/pC conversion coefficient 3
-   Int_t ADCtof[48];
-   //Float_t ADCsat=3100;  ci pensiamo in futuro !
+   Float_t pCthres=40.;     // threshold in charge
-   //Float_t pCsat=2500;
+   Int_t ADClast=4095;      // no signal --> ADC ch=4095
-   for(Int_t i=0; i<48; i++){
+   Int_t ADCsat=3100;       // saturation value for the ADCs
-     QevePmt_pC[i] = 0;
+   Int_t ADCtof[48];
-     ADCtof[i]=0;
-   }
+ // ---- introduce scale factors to tune simul ADC to real data   24-oct DC
-   // ------ read calibration file (get A1, A2, lambda1, lambda2)
+   Float_t ScaleFact[48]={0.18,0.22,0.35,0.26,0.47,0.35,0.31,0.37,
-   ifstream fileTriggerCalib;
+.44,0.23,0.38,0.60,0.39,0.29,0.40,0.23,
-   TString ftrigname="TrigCalibParam.txt";
+.30,0.66,0.22,1.53,0.17,0.55,
-   fileTriggerCalib.open(ftrigname.Data());
+.84,0.19,0.21,1.64,0.62,0.13,
-   if ( !fileTriggerCalib ) {
+.18,0.15,0.10,0.14,0.14,0.14,0.14,0.12,
-     printf("debug: no trigger calib file!\n");
+.26,0.18,0.25,0.23,0.20,0.40,
-     return(-117); //check output!
+.19,0.23,0.25,0.23,0.25,0.20};
-   };
-   Float_t atte1[48],atte2[48],lambda1[48],lambda2[48];
+   for(Int_t i=0; i<48; i++){
-   Int_t temp=0;
+     QevePmt_pC[i] = 0;
-   for(Int_t i=0; i<48; i++){
+     ADCtof[i]=0;
-     fileTriggerCalib >> temp;
+   }
-     fileTriggerCalib >> atte1[i];
-     fileTriggerCalib >> atte2[i];
+   // ------ read calibration file (get A1, A2, lambda1, lambda2)
-     fileTriggerCalib >> lambda1[i];
+   ifstream fileTriggerCalib;
-     fileTriggerCalib >> lambda2[i];
+   TString ftrigname="TrigCalibParam.txt";
-     fileTriggerCalib >> temp;
+   fileTriggerCalib.open(ftrigname.Data());
-   }
+   if ( !fileTriggerCalib ) {
-   fileTriggerCalib.close();
+     printf("debug: no trigger calib file!\n");
+     return(-117); //check output!
-   //  Read from file the 48*4 values of the attenuation fit function
+   };
-   //  NB: lambda<0; x,y defined in gpamela (=0 in the centre of the cavity)
+   Float_t atte1[48],atte2[48],lambda1[48],lambda2[48];
-   //    Qhitpmt_pC =  atte1 * exp(x/lambda1) + atte2 * exp(x/lambda2)
+   Int_t temp=0;
+   // correct readout WM Oct '07
-   //    fine lettura dal file */
+   for(Int_t i=0; i<48; i++){
+     fileTriggerCalib >> temp;
-   //const Int_t nmax=??; = Nthtof
+     fileTriggerCalib >> atte1[i];
-   Int_t ip, ipad;
+     fileTriggerCalib >> lambda1[i];
-   //Int_t ipmt;
+     fileTriggerCalib >> atte2[i];
-   Int_t pmtleft=0, pmtright=0;
+     fileTriggerCalib >> lambda2[i];
-   Int_t *pl, *pr;
+     fileTriggerCalib >> temp;
-   pl = &pmtleft;
+   }
-   pr = &pmtright;
+   fileTriggerCalib.close();
-   // TDC variables:
+   Int_t ip, ipad;
-   Int_t TDClast=4095;      // no signal --> ADC ch=4095
+   //Int_t ipmt;
-   Int_t TDCint[48];
+   Int_t pmtleft=0, pmtright=0;
-   Float_t  tdc[48],tdc1[48],tdcpmt[48];
+   Int_t *pl, *pr;
-   for(Int_t i=0; i<48; i++)
+   pl = &pmtleft;
-     tdcpmt[i] = 1000.;
+   pr = &pmtright;
-   Float_t thresh=1.; // to be defined better... (Wolfgang)
+   // TDC variables:
-     // === TDC: simulate timing for each paddle
+   Int_t TDClast=4095;      // no signal --> TDC ch=4095
-     Float_t dt1 = 285.e-12 ;   // single PMT resolution
+   Int_t TDCint[48];
-     Float_t tdcres[50],c1_S[50],c2_S[50],c3_S[50];
+   Float_t  tdc[48],tdc1[48],tdcpmt[48];
-     for(Int_t j=0;j<48;j++)  tdcres[j] = 50.E-12;   // TDC resolution 50 picosec
+   for(Int_t i=0; i<48; i++) {
-     for(Int_t j=0;j<48;j++)  c1_S[j] = 500.;  // cable length in channels
+     tdcpmt[i] = 1000.;
-     for(Int_t j=0;j<48;j++)  c2_S[j] = 0.;
+     tdc[i]  = 0.;            // 18-oct WM
-     for(Int_t j=0;j<48;j++)  c3_S[j] = 1000.;
+     tdc1[i] = 0.;            // 18-oct WM
-     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)
+   Float_t thresh=10.; // to be defined better... (Wolfgang)
-   /* **********************************  inizio loop sugli hit */
+   // === TDC: simulate timing for each paddle
+   Float_t dt1 = 285.e-12 ;   // single PMT resolution
-   for(Int_t nh=0; nh<Nthtof; nh++){
+   //    Float_t dt1 = 10.e-12 ;   // TEST
+   Float_t tdcres[50],c1_S[50],c2_S[50],c3_S[50];
-     for(Int_t j=0; j<2; j++) { // already done!! remove???
+   for(Int_t j=0;j<48;j++)  tdcres[j] = 50.E-12;   // TDC resolution 50 picosec
-       xpath[j]=0.;
+   for(Int_t j=0;j<48;j++)  c1_S[j] = 500.;  // cable length in channels
-       ypath[j]=0.;
+   for(Int_t j=0;j<48;j++)  c2_S[j] = 0.;
-       FGeo[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];
-     Float_t s_l_g[6] = {8.0, 8.0, 20.9, 22.0, 9.8, 8.3 };  // length of the lightguide
+   // ih = 0 + i1;  // not used?? (Silvio)
-     Float_t t1,t2,veff,veff1,veff0 ;
-     veff0 = 100.*1.0e8 ; // light velocity in the scintillator in m/sec
+   /* **********************************  start loop over hits */
-     veff1 = 100.*1.5e8; // light velocity in the lightguide in m/sec
-     veff=veff0;         // signal velocity in the paddle
+   for(Int_t nh=0; nh<Nthtof; nh++){
-     t1 = Timetof[nh] ;  // Start
+     Float_t s_l_g[6] = {8.0, 8.0, 20.9, 22.0, 9.8, 8.3 };  // length of the lightguide
-     t2 = Timetof[nh] ;
+     Float_t t1,t2,veff,veff1,veff0 ;
+     veff0 = 100.*1.0e8 ; // light velocity in the scintillator in m/sec
-     // Donatella
+     veff1 = 100.*1.5e8; // light velocity in the lightguide in m/sec
-     // ridefiniz. piano e pad per i vettori in C
+     veff=veff0;         // signal velocity in the paddle
-     ip = Ipltof[nh]-1;
-     ipad = Ipaddle[nh]-1;
+     t1 = Timetof[nh] ;  // Start
-     pmtleft=0;
+     t2 = Timetof[nh] ;
-     pmtright=0;
+     // Donatella: redefinition plane and pad for vectors in C
-     //Paddle2Pmt((Int_t)ip, (Int_t) ipad, (Int_t*) &pmtleft, (Int_t*) &pmtright);
+     ip = Ipltof[nh]-1;
-     Paddle2Pmt(ip, ipad, &pmtleft, &pmtright);
+     ipad = Ipaddle[nh]-1;
-     //Paddle2Pmt(ip, ipad, pl, pr);
+     pmtleft=0;
+     pmtright=0;
-     // per avere anche la corrispondenza pmt --> half board e canale
-     // metodo GetPMTIndex(Int_t ipmt, Int_t &hb, Int_t &ch) // non lo usiamo x ora
+     if (ip<6) {
+       Paddle2Pmt(ip, ipad, &pmtleft, &pmtright);
-     /*calcola la pos media e il path all'interno della paddle */
+       // DC: evaluates mean position and path inside the paddle
-     Float_t tpos=0.;
-     Float_t path[2] = {0., 0.};
+       Float_t tpos=0.;
-     //--- Strip in Y = S11,S22,S31 ------
+       Float_t path[2] = {0., 0.};
-     if(ip==0 || ip==3 || ip==4)
+       //--- Strip in Y = S11,S22,S31 ------
-       tpos = (Yintof[nh]+Youttof[nh])/2.;
+       if(ip==0 || ip==3 || ip==4)
-     else
+         tpos = (Yintof[nh]+Youttof[nh])/2.;
-       if(ip==1 || ip==2 || ip==5)   //--- Strip in X per S12,S21,S32
+       else
-         tpos = (Xintof[nh]+Xouttof[nh])/2.;
+         if(ip==1 || ip==2 || ip==5)   //--- Strip in X per S12,S21,S32
-       else if (ip!=6)
+           tpos = (Xintof[nh]+Xouttof[nh])/2.;
-         printf("*** Warning: this option should never occur! (ip=%2i, nh=%2i)\n",ip,nh);
+         else //if (ip!=6)
-     path[0]= tpos + dimel[ip]/2.;
+           printf("*** WARNING TOF: this option should never occur! (ip=%2i, nh=%2i)\n",ip,nh);
-     path[1]= dimel[ip]/2.- tpos;
+       path[0]= tpos + dimel[ip]/2.;   // path to left PMT
-     //  cout <<"Strip N. ="<< ipaddle <<" piano n.= "<< iplane <<" POSIZ = "<< tpos <<"\n";
+       path[1]= dimel[ip]/2.- tpos;    // path to right PMT
-     /* per il momento metto un fattore geometrico costante*/
+       //  cout <<"Strip N. ="<< ipaddle <<" piano n.= "<< iplane <<" POSIZ = "<< tpos <<"\n";
-     FGeo[0] =0.5;
-     FGeo[1] =0.5;
+       if (DEBUG) {
-     //  FGeo[1] = atan(path[1]/dimes[ip])/6.28318; // frazione fotoni verso SX
+         cout <<" plane "<<ip<<" strip # ="<< ipad <<" tpos  "<< tpos <<"\n";
-     //  FGeo[2] = atan(path[2]/dimes[ip])/6.28318; // e verso DX
+         cout <<"pmtleft, pmtright "<<pmtleft<<" "<<pmtright<<endl;
+       }
-     /*  rimando la fluttuazione poissoniana sui fotoni prodotti
-         sto studiando come funziona la funzione:
+       // constant geometric factor, for the moment
-         long int i = sto.Poisson(double x);  */
+       FGeo[0] =0.5;
-     //  Npho = Poisson(ERELTOF[nh])*Pho_keV*1e6   Eloss in GeV ?
+       FGeo[1] =0.5;
-     Npho = Ereltof[nh]*Pho_keV*10.0e6;  // Eloss in GeV ?
+       //  FGeo[1] = atan(path[1]/dimes[ip])/6.28318; // fraction of photons toward SX
+       //  FGeo[2] = atan(path[2]/dimes[ip])/6.28318; // toward DX
-     Float_t knorm[2]={0., 0.}; // Donatella
-     Float_t Atten[2]={0., 0.}; // Donatella
-     for(Int_t j=0; j<2; j++){
+       //  Npho = Poisson(ERELTOF[nh])*Pho_keV*1e6  Poissonian fluctuations to be inserted-DC
-       QhitPad_pC[j]= Npho*FGeo[j]*effi*pmGain*echarge;
+       Npho = Ereltof[nh]*Pho_keV*1.0e6;  // Eloss in GeV
-       knorm[j]=QhitPad_pC[j]/(atte1[pmtleft+j]*exp((dimel[ip]/2.*pow(-1,j+1))/lambda1[pmtleft+j]) +
-                               atte2[pmtleft+j]*exp((dimel[ip]/2.*pow(-1,j+1))/lambda2[pmtleft+j]));
+       Float_t knorm[2]={0., 0.}; // Donatella
+       Float_t Atten[2]={0., 0.}; // Donatella
-       Atten[j]=knorm[j]*(atte1[pmtleft+j]*exp(tpos/lambda1[pmtleft+j]) +
+       for(Int_t j=0; j<2; j++){
-                          atte2[pmtleft+j]*exp(tpos/lambda2[pmtleft+j]));
+         QhitPad_pC[j]= Npho*FGeo[j]*effi*pmGain*echarge*1.E12; // corrected WM
+         // WM
-       QhitPmt_pC[j]= QhitPad_pC[j]*Atten[j];
+         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]) +
-     QevePmt_pC[pmtleft]  += QhitPmt_pC[0];
+           atte2[pmtleft+j]*exp(tpos*lambda2[pmtleft+j]) ;
-     QevePmt_pC[pmtright] += QhitPmt_pC[1];
+         QhitPmt_pC[j]= QhitPad_pC[j]*Atten[j]/knorm[j];
+         //              QhitPmt_pC[j]= QhitPad_pC[j]; //no attenuation
-     // TDC
-     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;
+         if (DEBUG) {
+           cout<<"pmtleft "<<pmtleft<<" j "<<j<<endl;
-     TRandom r;
+           cout<<" atte1 "<<atte1[pmtleft+j]<<"lambda1 "<<lambda1[pmtleft+j]<<" atte2 "<<atte2[pmtleft+j]<<"lambda2 "<<lambda2[pmtleft+j] <<endl;
-     t1 = r.Gaus(t1,dt1);  //apply gaussian error  dt
+           cout<<j<<" tpos "<<tpos<<" knorm "<<knorm[j]<<" "<<Atten[j]<<" "<<"QhitPmt_pC "<<QhitPmt_pC[j]<<endl;
-     t2 = r.Gaus(t2,dt1);  //apply gaussian error  dt
+         }
+       }
-     t1 = t1 + c1_S[pmtleft] ;  // Signal reaches Discriminator ,TDC starts  to run
-     t2 = t2 + c1_S[pmtright] ;
+       if (DEBUG)
+         cout<<"Npho "<<Npho<<" QhitPmt_pC "<<QhitPmt_pC[0]<<" "<<QhitPmt_pC[1]<<endl;
-     // check if signal is above threshold
-     // then check if tdcpmt is already filled by another hit...
+       QevePmt_pC[pmtleft]  += QhitPmt_pC[0];
-     // only re-fill if time is smaller
+       QevePmt_pC[pmtright] += QhitPmt_pC[1];
-     if (QhitPmt_pC[0] > thresh)
+       // TDC
-       if (tdcpmt[pmtleft] < 1000.) // is already filled!
+ // WM right and left <->
-         if (t1 <  tdcpmt[pmtleft]) {
+ //      t2 = t2 + fabs(path[0]/veff) + s_l_g[ip]/veff1 ;  // Signal reaches PMT
-           tdcpmt[pmtleft] = t1;
+ //      t1 = t1 + fabs(path[1]/veff) + s_l_g[ip]/veff1;
-           t1 = t1 + c2_S[pmtleft] ;  // Signal reaches Coincidence
-           tdc[pmtleft] = t1;
+       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
-     if (QhitPmt_pC[1] > thresh)
+       Float_t t1save = t1;
-         if (tdcpmt[pmtright] < 1000.)  // is already filled!
+       Float_t t2save = t2;
-           if (t2 <  tdcpmt[pmtright]) {
-             tdcpmt[pmtright] = t2;
+ /*
-             t2 = t2 + c2_S[pmtright] ;
+       TRandom r;
-             tdc[pmtright] = t2;
+ // This does not work... WM -  but works in my simulation code ??
-           }
+ //      t1 = r.Gaus(t1,dt1);  //apply gaussian error  dt
+ //      t2 = r.Gaus(t2,dt1);  //apply gaussian error  dt
-   } // ****************************************       end loop over hits
+ */
+         t1 = gRandom->Gaus(t1,dt1); //apply gaussian error  dt
-   // ======  ADC ======
+         t2 = gRandom->Gaus(t2,dt1); //apply gaussian error  dt
-   for(Int_t i=0; i<48; i++){
-     if(QevePmt_pC[i] != 0.){
+ //      cout<<1E12*(t1save-t1)<<" "<<1E12*(t2save-t2)<<endl;
-       ADCtof[i]= (Int_t)(ADC_pC*QevePmt_pC[i] + ADCoffset);
-       if(ADCtof[i]> ADClast) ADCtof[i]=ADClast;
+       t1 = t1 + c1_S[pmtleft] ;  // Signal reaches Discriminator ,TDC starts  to run
-     } else
+       t2 = t2 + c1_S[pmtright] ;
-       ADCtof[i]= ADClast;
-   };
+       // check if signal is above threshold
+       // then check if tdcpmt is already filled by another hit...
+       // only re-fill if time is smaller
-   // ======  build  TDC coincidence  ======
+       if (QhitPmt_pC[0] > thresh) {
-   Float_t t_coinc = 0;
+         if (tdcpmt[pmtleft] == 1000.) {  // fill for the first time
-   Int_t ilast = 100;
+           tdcpmt[pmtleft] = t1;
-   for (Int_t ii=0; ii<48;ii++)
+           tdc[pmtleft] = t1 + c2_S[pmtleft] ;  // Signal reaches Coincidence
-     if (tdc[ii] > t_coinc) {
+         }
-       t_coinc = tdc[ii];
+         if (tdcpmt[pmtleft] < 1000.) // is already filled!
-       ilast = ii;
+           if (t1 <  tdcpmt[pmtleft]) {
-     }
+             tdcpmt[pmtleft] = t1;
+             t1 = t1 + c2_S[pmtleft] ;  // Signal reaches Coincidence
-   //     cout<<ilast<<" "<<t_coinc<<endl;
+             tdc[pmtleft] = t1;
-   //     At t_coinc  trigger condition is fulfilled
+           }
+       }
-   for (Int_t ii=0; ii<48;ii++){
+       if (QhitPmt_pC[1] > thresh) {
-     //      if (tdc[ii] != 0) tdc1[ii] = t_coinc - tdc[ii];   // test 1
+         if (tdcpmt[pmtright] == 1000.) {  // fill for the first time
-     if (tdc[ii] != 0) tdc1[ii] = t_coinc - tdcpmt[ii];  // test 2
+           tdcpmt[pmtright] = t2;
-     tdc1[ii] = tdc1[ii]/tdcres[ii];                     // divide by TDC resolution
+           tdc[pmtright] = t2 + c2_S[pmtright] ;  // Signal reaches Coincidence
-     if (tdc[ii] != 0) tdc1[ii] = tdc1[ii] + c3_S[ii];  // add cable length c3
+         }
+         if (tdcpmt[pmtright] < 1000.)  // is already filled!
-   } // missing parenthesis inserted! (Silvio)
+           if (t2 <  tdcpmt[pmtright]) {
+             tdcpmt[pmtright] = t2;
-   for(Int_t i=0; i<48; i++){
+             t2 = t2 + c2_S[pmtright] ;
-     if(tdc1[i] != 0.){
+             tdc[pmtright] = t2;
-       TDCint[i]=(Int_t)tdc1[i];
+           }
-       //ADC[i]= ADC_pC * QevePmt_pC[i] + ADCoffset;
+       }
-       //if(ADC[i]> ADClast) ADC[i]=ADClast;
-     } else
+       if (DEBUG)
-       TDCint[i]= TDClast;
+         cout<<nh<<" "<<Timetof[nh]<<" "<<t1<<" "<<t2<<endl;
-   }
+     } // ip < 6
- // ======  write fDataTof  =======
-   UChar_t tofBin;
+   }; // ****************************************       end loop over hits
-   for (Int_t j=0; j < 12; j++){
-     Int_t j12=j*12;
+   // ======  ADC ======
-     fDataTof[j12+0]=0x00;   // TDC_ID
-     fDataTof[j12+1]=0x00;   // EV_COUNT
-     fDataTof[j12+2]=0x00;   // TDC_MASK (1)
+   for(Int_t i=0; i<48; i++){
-     fDataTof[j12+3]=0x00;   // TDC_MASK (2)
+     if(QevePmt_pC[i] >= pCthres){
-     for (Int_t k=0; k < 4; k++){
+       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));
-       Int_t jk12=j12+k;
+   } else
-       tofBin=(UChar_t)(ADCtof[k+4*j]/256);   // ADC# (msb) (#=1+k+4*j)
+       ADCtof[i]= ADClast;
-       fDataTof[jk12+4] = Bin2GrayTof(tofBin,fDataTof[jk12+4]);
+ }
-       tofBin=(UChar_t)(ADCtof[k+4*j]%256);   // ADC# (lsb)
-       fDataTof[jk12+5] = Bin2GrayTof(tofBin,fDataTof[jk12+5]);
+ // ---- introduce scale factors to tune simul ADC to real data   24-oct DC
-       tofBin=(UChar_t)(TDCint[k+4*j]/256);   // TDC# (msb)
-       fDataTof[jk12+6]=Bin2GrayTof(tofBin,fDataTof[jk12+6]);
+   for(Int_t i=0; i<48; i++){
-       tofBin=(UChar_t)(TDCint[k+4*j]%256);   // TDC# (lsb)
+     if(ADCtof[i] != ADClast){
-       fDataTof[jk12+7]=Bin2GrayTof(tofBin,fDataTof[jk12+7]);
+                              //      printf("%3d, %4d, %4.2f\n",i, ADCtof[i],ScaleFact[i]);
-     };
+       ADCtof[i]= Int_t (ADCtof[i]*ScaleFact[i]);
-     fDataTof[j12+20]=0x00;   // TEMP1
+                              //      printf("%3d, %4d,\n",i, ADCtof[i]);
-     fDataTof[j12+21]=0x00;   // TEMP2
+     }
-     fDataTof[j12+22]= EvaluateCrcTof(pTof);   // CRC
+   }
-     pTof+=23;
-   };
+   for(Int_t i=0; i<48; i++){
-   return(0);
+     if(ADCtof[i] != ADClast){
- };
+       if(ADCtof[i]> ADCsat) ADCtof[i]=ADCsat;
+       else if(ADCtof[i]< 0) ADCtof[i]=ADClast;
- UChar_t Digitizer::Bin2GrayTof(UChar_t binaTOF,UChar_t grayTOF){
+   }
-   union graytof_data {
+   }
-     UChar_t word;
+   // ======  build  TDC coincidence  ======
-     struct bit_field {
-       unsigned b0:1;
+   Float_t t_coinc = 0;
-       unsigned b1:1;
+   Int_t ilast = 100;
-       unsigned b2:1;
+   for (Int_t ii=0; ii<48;ii++)
-       unsigned b3:1;
+     if (tdc[ii] > t_coinc) {
-       unsigned b4:1;
+       t_coinc = tdc[ii];
-       unsigned b5:1;
+       ilast = ii;
-       unsigned b6:1;
+     }
-       unsigned b7:1;
-     } bit;
+   //     cout<<ilast<<" "<<t_coinc<<endl;
-   } bi,gr;
+   //     At t_coinc  trigger condition is fulfilled
-   //
-   bi.word = binaTOF;
+   for (Int_t ii=0; ii<48;ii++){
-   gr.word = grayTOF;
+     //      if (tdc[ii] != 0) tdc1[ii] = t_coinc - tdc[ii];   // test 1
-   //
+     if (tdc[ii] != 0) tdc1[ii] = t_coinc - tdcpmt[ii];  // test 2
-   gr.bit.b0 = bi.bit.b1 ^ bi.bit.b0;
+     tdc1[ii] = tdc1[ii]/tdcres[ii];                     // divide by TDC resolution
-   gr.bit.b1 = bi.bit.b2 ^ bi.bit.b1;
+     if (tdc[ii] != 0) tdc1[ii] = tdc1[ii] + c3_S[ii];  // add cable length c3
-   gr.bit.b2 = bi.bit.b3 ^ bi.bit.b2;
-   gr.bit.b3 = bi.bit.b3;
+   } // missing parenthesis inserted! (Silvio)
-   //
-   /* bin to gray conversion 4 bit per time*/
+   for(Int_t i=0; i<48; i++){
-   //
+     if(tdc1[i] != 0.){
-   gr.bit.b4 = bi.bit.b5 ^ bi.bit.b4;
+       TDCint[i]=(Int_t)tdc1[i];
-   gr.bit.b5 = bi.bit.b6 ^ bi.bit.b5;
+       if (TDCint[i]>4093) TDCint[i]=TDClast;  // 18-oct WM
-   gr.bit.b6 = bi.bit.b7 ^ bi.bit.b6;
+       if (DEBUG)
-   gr.bit.b7 = bi.bit.b7;
+         cout<<i<<" "<<TDCint[i]<<endl;
-   //
+       //ADC[i]= ADC_pC * QevePmt_pC[i] + ADCoffset;
-   return(gr.word);
+       //if(ADC[i]> ADClast) ADC[i]=ADClast;
- }
+     } else
+       TDCint[i]= TDClast;
- UChar_t Digitizer::EvaluateCrcTof(UChar_t *pTof) {
+   }
-   // UChar_t crcTof=0x00;
-   //   for (Int_t jp=0; jp < 23; jp++){
+   if (DEBUG)
-   //     crcTof = crc8(...)
+     cout<<"-----------"<<endl;
-   //   }
-   return(0x00);
- };
+ //------ use channelmap  18-oct WM
- //void Digitizer::Paddle2Pmt(Int_t plane, Int_t paddle, Int_t* &pmtleft, Int_t* &pmtright){
+ Int_t  channelmap[] =  {3,21,11,29,19,45,27,37,36,28,44,20,5,12,13,4,
- void Digitizer::Paddle2Pmt(Int_t plane, Int_t paddle, Int_t *pl, Int_t *pr){
+,47,14,39,22,31,30,23,38,15,46,7,0,33,16,24,
-   //* @param plane    (0 - 5)
+,41,32,40,25,17,34,9,42,1,2,10,18,26,35,43};
-   //* @param paddle   (plane=0, paddle = 0,...5)
-   //* @param padid    (0 - 23)
+   Int_t ADChelp[48];
-   //
+   Int_t TDChelp[48];
-   Int_t padid=-1;
-   Int_t pads[6]={8,6,2,2,3,3};
+   for(Int_t i=0; i<48; i++){
-   //
+       Int_t ii=channelmap[i];
-   Int_t somma=0;
+       ADChelp[ii]= ADCtof[i];
-   Int_t np=plane;
+       TDChelp[ii]= TDCint[i];
-   for(Int_t j=0; j<np; j++)
+                           }
-     somma+=pads[j];
-   padid=paddle+somma;
+   for(Int_t i=0; i<48; i++){
-   *pl = padid*2;
+       ADCtof[i]= ADChelp[i];
-   *pr = *pr + 1;
+       TDCint[i]= TDChelp[i];
- };
+                           }
- void Digitizer::DigitizeAC() {
-   // created:  J. Conrad, KTH
+ /*
-   // modified: S. Orsi, INFN Roma2
+ //--- fake data ------------------------
+   for(Int_t i=0; i<48; i++){
-   fDataAC[0] = 0xACAC;
+       ADCtof[i]= 100 + 10*i;
-   fDataAC[64]= 0xACAC;
+       TDCint[i]= 800 + 10*i;
-   fDataAC[1] = 0xAC11;   // main card
+ //      cout<<i<<" "<<ADCtof[i]<<" "<<TDCint[i]<<endl;
-   fDataAC[65] = 0xAC22;   // extra card
+                           }
+ */
-   // the third word is a status word (dummy)
-   fDataAC[2] = 0xFFFF; //FFEF?
+ /*
-   fDataAC[66] = 0xFFFF;
+   for(Int_t i=0; i<48; i++){
+    if (((ADCtof[i]>0)&&(ADCtof[i]<4095)) || ((TDCint[i]>0)&&(TDCint[i]<4095)))  cout<<i<<" "<<ADCtof[i]<<" "<<TDCint[i]<<endl;
-   const UInt_t nReg = 6;
+                           }
+ */
-   // Registers (dummy)
-   for (UInt_t i=0; i<=nReg; i++){
-     fDataAC[i+4] = 0xFFFF;
+ // ======  write fDataTof  =======
-     fDataAC[i+68] = 0xFFFF;
-   }
+ //  UChar_t tdcadd[8]={1,0,3,2,5,4,7,6};  (coded in 3 bit)
-   // the last word is a CRC
+   UChar_t Ctrl3bit[8]={32,0,96,64,160,128,224,192};  // DC (msb in 8 bit word )
-   // Dummy for the time being, but it might need to be calculated in the end
-   fDataAC[63] = 0xABCD;
+   UChar_t tofBin;
-   fDataAC[127] = 0xABCD;
+   for (Int_t j=0; j < 12; j++){   // loop on TDC #12
+     Int_t j12=j*23;               // for each TDC 23 bytes (8 bits)
-   // shift registers, which one is with respect to PMT, where in
+     fDataTof[j12+0]=0x00;   // TDC_ID
-   // shift registers is a question of time relative trigger
+     fDataTof[j12+1]=0x00;   // EV_COUNT
-   // In level2: hitmap, hitmap-status (synchronised with a trigger),
+     fDataTof[j12+2]=0x00;   // TDC_MASK (1)
-   // status
+     fDataTof[j12+3]=0x00;   // TDC_MASK (2)
+     for (Int_t k=0; k < 4; k++){   // for each TDC 4 channels (ADC+TDC)
-   for (UInt_t i=0; i<=15; i++){
-     fDataAC[i+11] = 0x0000;
+       Int_t jk12=j12+4*k;         // ADC,TDC channel (0-47)
-     fDataAC[i+75] = 0x0000;
-   }
+       tofBin =(UChar_t)(ADCtof[k+4*j]/256);   // ADC# (msb)
+       fDataTof[jk12+4] = Bin2GrayTof(tofBin,fDataTof[jk12+4]);
-   // singles counters are dummy
+       /* control bits inserted here, after the bin to gray conv - DC*/
+       fDataTof[jk12+4] = Ctrl3bit[2*k] | fDataTof[jk12+4];
-   for (UInt_t i=0; i<=16; i++){
+       tofBin=(UChar_t)(ADCtof[k+4*j]%256);   // ADC# (lsb)
-     fDataAC[i+26] = 0x0000;
+       fDataTof[jk12+5] = Bin2GrayTof(tofBin,fDataTof[jk12+5]);
-     fDataAC[i+90] = 0x0000;
+       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*/
-   // coincidences are dummy
+       fDataTof[jk12+6] = Ctrl3bit[2*k+1] | fDataTof[jk12+6];
+       tofBin=(UChar_t)(TDCint[k+4*j]%256);   // TDC# (lsb)
-   for (UInt_t i=0; i<=7; i++){
+       fDataTof[jk12+7]=Bin2GrayTof(tofBin,fDataTof[jk12+7]);
-     fDataAC[i+42] = 0x0000;
+     };
-     fDataAC[i+106] = 0x0000;
+     fDataTof[j12+20]=0x00;   // TEMP1
-   }
+     fDataTof[j12+21]=0x00;   // TEMP2
+     fDataTof[j12+22]= EvaluateCrcTof(pTof);   // CRC
-   // increments for every trigger might be needed at some point.
+     pTof+=23;
-   // dummy for now
+   };
-   fDataAC[50]  = 0x0000;
+   return(0);
-   fDataAC[114] = 0x0000;
+ };
-   // dummy FPGA clock
+ UChar_t Digitizer::Bin2GrayTof(UChar_t binaTOF,UChar_t grayTOF){
-   fDataAC[51] = 0x006C;
+   union graytof_data {
-   fDataAC[52] = 0x6C6C;
+     UChar_t word;
-   fDataAC[115] = 0x006C;
+     struct bit_field {
-   fDataAC[116] = 0x6C6C;
+       unsigned b0:1;
+       unsigned b1:1;
+       unsigned b2:1;
-   // dummy temperatures
+       unsigned b3:1;
-   fDataAC[53] = 0x0000;
+       unsigned b4:1;
-   fDataAC[54] = 0x0000;
+       unsigned b5:1;
-   fDataAC[117] = 0x0000;
+       unsigned b6:1;
-   fDataAC[118] = 0x0000;
+       unsigned b7:1;
+     } bit;
+   } bi,gr;
-   // dummy DAC thresholds
+   //
-   for (UInt_t i=0; i<=7; i++){
+   bi.word = binaTOF;
-     fDataAC[i+55] = 0x1A13;
+   gr.word = grayTOF;
-     fDataAC[i+119] = 0x1A13;
+   //
-   }
+   gr.bit.b0 = bi.bit.b1 ^ bi.bit.b0;
+   gr.bit.b1 = bi.bit.b2 ^ bi.bit.b1;
-   // We activate all branches. Once the digitization algorithm
+   gr.bit.b2 = bi.bit.b3 ^ bi.bit.b2;
-   // is determined only the branches need to activated which involve needed
+   gr.bit.b3 = bi.bit.b3;
-   // information
+   //
+   /* bin to gray conversion 4 bit per time*/
-   fhBookTree->SetBranchStatus("Nthcat",1);
+   //
-   fhBookTree->SetBranchStatus("Iparcat",1);
+   gr.bit.b4 = bi.bit.b5 ^ bi.bit.b4;
-   fhBookTree->SetBranchStatus("Icat",1);
+   gr.bit.b5 = bi.bit.b6 ^ bi.bit.b5;
-   fhBookTree->SetBranchStatus("Xincat",1);
+   gr.bit.b6 = bi.bit.b7 ^ bi.bit.b6;
-   fhBookTree->SetBranchStatus("Yincat",1);
+   gr.bit.b7 = bi.bit.b7;
-   fhBookTree->SetBranchStatus("Zincat",1);
+   //
-   fhBookTree->SetBranchStatus("Xoutcat",1);
+   return(gr.word);
-   fhBookTree->SetBranchStatus("Youtcat",1);
+ }
-   fhBookTree->SetBranchStatus("Zoutcat",1);
-   fhBookTree->SetBranchStatus("Erelcat",1);
+ UChar_t Digitizer::EvaluateCrcTof(UChar_t *pTof) {
-   fhBookTree->SetBranchStatus("Timecat",1);
+   Bool_t DEBUG=false;
-   fhBookTree->SetBranchStatus("Pathcat",1);
+   if (DEBUG)
-   fhBookTree->SetBranchStatus("P0cat",1);
+     return(0x00);
-   fhBookTree->SetBranchStatus("Nthcas",1);
-   fhBookTree->SetBranchStatus("Iparcas",1);
+   UChar_t crcTof=0x00;
-   fhBookTree->SetBranchStatus("Icas",1);
+   UChar_t *pc=&crcTof, *pc2;
-   fhBookTree->SetBranchStatus("Xincas",1);
+   pc2=pTof;
-   fhBookTree->SetBranchStatus("Yincas",1);
+   for (Int_t jp=0; jp < 23; jp++){
-   fhBookTree->SetBranchStatus("Zincas",1);
+     //crcTof = crc8(...)
-   fhBookTree->SetBranchStatus("Xoutcas",1);
+     Crc8Tof(pc2++,pc);
-   fhBookTree->SetBranchStatus("Youtcas",1);
+     //    printf("%2i --- %x\n",jp,crcTof);
-   fhBookTree->SetBranchStatus("Zoutcas",1);
+   }
-   fhBookTree->SetBranchStatus("Erelcas",1);
+   return(crcTof);
-   fhBookTree->SetBranchStatus("Timecas",1);
+ }
-   fhBookTree->SetBranchStatus("Pathcas",1);
-   fhBookTree->SetBranchStatus("P0cas",1);
+ void Digitizer::Crc8Tof(UChar_t *oldCRC, UChar_t *crcTof){
-   fhBookTree->SetBranchStatus("Nthcard",1);
+   union crctof_data {
-   fhBookTree->SetBranchStatus("Iparcard",1);
+     UChar_t word;
-   fhBookTree->SetBranchStatus("Icard",1);
+     struct bit_field {
-   fhBookTree->SetBranchStatus("Xincard",1);
+       unsigned b0:1;
-   fhBookTree->SetBranchStatus("Yincard",1);
+       unsigned b1:1;
-   fhBookTree->SetBranchStatus("Zincard",1);
+       unsigned b2:1;
-   fhBookTree->SetBranchStatus("Xoutcard",1);
+       unsigned b3:1;
-   fhBookTree->SetBranchStatus("Youtcard",1);
+       unsigned b4:1;
-   fhBookTree->SetBranchStatus("Zoutcard",1);
+       unsigned b5:1;
-   fhBookTree->SetBranchStatus("Erelcard",1);
+       unsigned b6:1;
-   fhBookTree->SetBranchStatus("Timecard",1);
+       unsigned b7:1;
-   fhBookTree->SetBranchStatus("Pathcard",1);
+     } bit;
-   fhBookTree->SetBranchStatus("P0card",1);
+   } c,d,r;
-   // In this simpliefied approach we will assume that once
+   c.word = *oldCRC;
-   // a particle releases > 0.5 mip in one of the 12 AC detectors it
+   //d.word = *newCRC;
-   // will fire. We will furthermore assume that both cards read out
+   d.word = *crcTof;
-   // identical data.
+   r.word = 0;
-   // If you develop you digitization algorithm, you should start by
+   r.bit.b0 = c.bit.b7 ^ c.bit.b6 ^ c.bit.b0 ^
-   // identifying the information present in level2 (post-darth-vader)
+              d.bit.b0 ^ d.bit.b6 ^ d.bit.b7;
-   // data.
+   r.bit.b1 = c.bit.b6 ^ c.bit.b1 ^ c.bit.b0 ^
-   Float_t SumEcat[5];
+              d.bit.b0 ^ d.bit.b1 ^ d.bit.b6;
-   Float_t SumEcas[5];
-   Float_t SumEcard[5];
+   r.bit.b2 = c.bit.b6 ^ c.bit.b2 ^ c.bit.b1 ^ c.bit.b0 ^
-   for (Int_t k= 0;k<5;k++){
+              d.bit.b0 ^ d.bit.b1 ^ d.bit.b2 ^ d.bit.b6;
-     SumEcat[k]=0.;
-     SumEcas[k]=0.;
+   r.bit.b3 = c.bit.b7 ^ c.bit.b3 ^ c.bit.b2 ^ c.bit.b1 ^
-     SumEcard[k]=0.;
+              d.bit.b1 ^ d.bit.b2 ^ d.bit.b3 ^ d.bit.b7;
-   };
+   r.bit.b4 = c.bit.b4 ^ c.bit.b3 ^ c.bit.b2 ^
-   if (Nthcat>50 || Nthcas>50 || Nthcard>50)
+              d.bit.b2 ^ d.bit.b3 ^ d.bit.b4;
-     printf("Error! NthAC out of range!\n\n");
+   r.bit.b5 = c.bit.b5 ^ c.bit.b4 ^ c.bit.b3 ^
-   // look in CAT
+              d.bit.b3 ^ d.bit.b4 ^ d.bit.b5;
-   //  for (UInt_t k= 0;k<50;k++){
-   for (Int_t k= 0;k<Nthcat;k++){
+   r.bit.b6 = c.bit.b6 ^ c.bit.b5 ^ c.bit.b4 ^
-     if (Erelcat[k] > 0)
+              d.bit.b4 ^ d.bit.b5 ^ d.bit.b6;
-       SumEcat[Icat[k]] += Erelcat[k];
-   };
+   r.bit.b7 = c.bit.b7 ^ c.bit.b6 ^ c.bit.b5 ^
+              d.bit.b5 ^ d.bit.b6 ^ d.bit.b7 ;
-   // look in CAS
-   for (Int_t k= 0;k<Nthcas;k++){
+   *crcTof=r.word;
-     if (Erelcas[k] >0)
+   //return r.word;
-       SumEcas[Icas[k]] += Erelcas[k];
+ };
-   };
+ //void Digitizer::Paddle2Pmt(Int_t plane, Int_t paddle, Int_t* &pmtleft, Int_t* &pmtright){
-   // look in CARD
+ void Digitizer::Paddle2Pmt(Int_t plane, Int_t paddle, Int_t *pl, Int_t *pr){
-   for (Int_t k= 0;k<Nthcard;k++){
+   //* @param plane    (0 - 5)
-     if (Erelcard[k] >0)
+   //* @param paddle   (plane=0, paddle = 0,...5)
-       SumEcard[Icard[k]] += Erelcard[k];
+   //* @param padid    (0 - 23)
-   };
+   //
+   Int_t padid=-1;
-   // channel mapping              Hit Map
+   Int_t pads[6]={8,6,2,2,3,3};
-   // 1 CARD4                      0          LSB
+   //
-   // 2 CAT2                       0
+   Int_t somma=0;
-   // 3 CAS1                       0
+   Int_t np=plane;
-   // 4 NC                         0
+   for(Int_t j=0; j<np; j++)
-   // 5 CARD2                      0
+     somma+=pads[j];
-   // 6 CAT4                       1
+   padid=paddle+somma;
-   // 7 CAS4                       0
+   *pl = padid*2;
-   // 8 NC                         0
+   //  *pr = *pr + 1;
-   // 9 CARD3                      0
+   *pr = *pl + 1; // WM
-   // 10 CAT3                      0
+ };
-   // 11 CAS3                      0
-   // 12 NC                        0
+ void Digitizer::DigitizeAC() {
-   // 13 CARD1                     0
+   // created:  J. Conrad, KTH
-   // 14 CAT1                      0
+   // modified: S. Orsi, INFN Roma2
-   // 15 CAS2                      0
+   // fDataAC[0-63]:   main AC board
-   // 16 NC                        0          MSB
+   // fDataAC[64-127]: extra AC board
-   // In the first version only the hit-map is filled, not the SR.
+   fDataAC[0] = 0xACAC;
+   fDataAC[64]= 0xACAC;
-   // Threshold: 0.8 MeV.
+   fDataAC[1] = 0xAC11;
+   fDataAC[65] = 0xAC22;
-   Float_t thr = 8e-4;
+   // the third word is a status word (dummy: "no errors are present in the AC boards")
-   fDataAC[3] = 0x0000;
+   fDataAC[2] = 0xFFFF; //FFEF?
+   fDataAC[66] = 0xFFFF;
-   if (SumEcas[0] > thr)  fDataAC[3]  = 0x0004;
-   if (SumEcas[1] > thr)  fDataAC[3] += 0x4000;
+   const UInt_t nReg = 6;
-   if (SumEcas[2] > thr)  fDataAC[3] += 0x0400;
-   if (SumEcas[3] > thr)  fDataAC[3] += 0x0040;
+   // FPGA Registers (dummy)
+   for (UInt_t i=0; i<=nReg; i++){
-   if (SumEcat[0] > thr)  fDataAC[3] += 0x2000;
+     fDataAC[i+4] = 0xFFFF;
-   if (SumEcat[1] > thr)  fDataAC[3] += 0x0002;
+     fDataAC[i+68] = 0xFFFF;
-   if (SumEcat[2] > thr)  fDataAC[3] += 0x0200;
+   }
-   if (SumEcat[3] > thr)  fDataAC[3] += 0x0020;
+   // the last word is a CRC
-   if (SumEcard[0] > thr)  fDataAC[3] += 0x1000;
+   // Dummy for the time being, but it might need to be calculated in the end
-   if (SumEcard[1] > thr)  fDataAC[3] += 0x0010;
+   fDataAC[63] = 0xABCD;
-   if (SumEcard[2] > thr)  fDataAC[3] += 0x0100;
+   fDataAC[127] = 0xABCD;
-   if (SumEcard[3] > thr)  fDataAC[3] += 0x0001;
+   // shift registers (moved to the end of the routine)
-   fDataAC[67] = fDataAC[3];
+   Int_t evntLSB=Ievnt%65536;
-   //    for (Int_t i=0; i<fACbuffer; i++){
+   Int_t evntMSB=(Int_t)(Ievnt/65536);
-   //    printf("%0x  ",fDataAC[i]);
-   //    if ((i+1)%8 ==0) cout << endl;
+   // 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;
- void Digitizer::DigitizeS4(){
+   };
-   // creato:  S. Borisov, INFN Roma2 e MEPHI, Sett 2007
-   TString ciao,modo="ns";
+   // coincidences are dummy (increment by 1 at each event)
-   Int_t i,j,t,NdF,pmt,NdFT,S4,S4v=0,S4p=32;
+   // for (UInt_t i=0; i<=7; i++){
-   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.;
+   //    fDataAC[i+42] = 0x0000;
-   Xs[0]=-24.1;
+   //    fDataAC[i+106] = 0x0000;
-   Xs[1]=24.1;
+   //   }
-   Ys[0]=-24.1;
+   for (UInt_t i=0; i<=7; i++){
-   Ys[1]=24.1;
+     fDataAC[i+42] = evntLSB;
-   Zs[0]=-0.5;
+     fDataAC[i+106] = evntLSB;
-   Zs[1]=0.5;
+   };
-   Yp[0]=-20.;
-   Yp[2]=-1.;
+   // increments for every trigger might be needed at some point.
-   Yp[4]=17.;
+   // dummy for now
-   for(i=0;i<3;i++)
+   fDataAC[50]  = 0x0000;
-     Yp[2*i+1]=Yp[2*i]+3;
+   fDataAC[114] = 0x0000;
-   srand(time(NULL));
-   // --- activate branches:
+   // dummy FPGA clock (increment by 1 at each event)
-   fhBookTree->SetBranchStatus("Nthtof",1);
+   /*
-   fhBookTree->SetBranchStatus("Ipltof",1);
+     fDataAC[51] = 0x006C;
-   fhBookTree->SetBranchStatus("Ipaddle",1);
+     fDataAC[52] = 0x6C6C;
+     fDataAC[115] = 0x006C;
-   fhBookTree->SetBranchStatus("Xintof",1);
+     fDataAC[116] = 0x6C6C;
-   fhBookTree->SetBranchStatus("Yintof",1);
+   */
-   fhBookTree->SetBranchStatus("Xouttof",1);
+   if (Ievnt<=0xFFFF) {
-   fhBookTree->SetBranchStatus("Youttof",1);
+     fDataAC[51] = 0x0000;
+     fDataAC[52] = Ievnt;
-   fhBookTree->SetBranchStatus("Ereltof",1);
+     fDataAC[115] = 0x0000;
-   fhBookTree->SetBranchStatus("Timetof",1);
+     fDataAC[116] = Ievnt;
-   NdFT=0;
+   } else {
-   Ert=0;
+     fDataAC[51] = evntMSB;
-   for(i=0;i<Nthtof;i++){
+     fDataAC[52] = evntLSB;
-     if(Ipltof[i]!=6) continue;
+     fDataAC[115] = fDataAC[51];
-     Ert+=Ereltof[i];
+     fDataAC[116] = fDataAC[52];
+   }
-     if(modo=="ns") continue;
+   // dummy temperatures
-     NdF=Int_t(Ereltof[i]/E1);
+   fDataAC[53] = 0x0000;
-     NdFT=0;
+   fDataAC[54] = 0x0000;
-     X=Xintof[i];
+   fDataAC[117] = 0x0000;
-     Y=Yintof[i];
+   fDataAC[118] = 0x0000;
-     Z=((Float_t)random()/(Float_t)0x7fffffff)-0.5;
-     //cout<<"XYZ "<<X<<"  "<<Y<<"  "<<Z<<endl;
-     for(j=0;j<NdF;j++){
+   // dummy DAC thresholds
-       q=(Float_t)random()/(Float_t)0x7fffffff;
+   for (UInt_t i=0; i<=7; i++){
-       w=(Float_t)random()/(Float_t)0x7fffffff;
+     fDataAC[i+55] = 0x1A13;
-       // cout<<"qw "<<q<<" "<<w<<endl;
+     fDataAC[i+119] = 0x1A13;
-       V[0]=p*cos(6.28318*q);
+   }
-       V[1]=p*sin(6.28318*q);
-       V[2]=p*(2.*w-1.);
+   // We activate all branches. Once the digitization algorithm is determined
-       pmt=0;
+   // only the branches that involve needed information will be activated
-       x=X;
-       y=Y;
+   fhBookTree->SetBranchAddress("Ievnt",&Ievnt);
-       z=Z;
+   fhBookTree->SetBranchStatus("Nthcat",1);
-       while(pmt==0 && (x>Xs[0] && x<Xs[1])&&(y>Ys[0] && y<Ys[1])&&(z>Zs[0] && z<Zs[1])){
+   fhBookTree->SetBranchStatus("Iparcat",1);
-         l=0;
+   fhBookTree->SetBranchStatus("Icat",1);
-         while(pmt==0 && (x>Xs[0] && x<Xs[1])&&(y>Ys[0] && y<Ys[1])&&(z>Zs[0] && z<Zs[1])){
+   fhBookTree->SetBranchStatus("Xincat",1);
-           x+=V[0];
+   fhBookTree->SetBranchStatus("Yincat",1);
-           y+=V[1];
+   fhBookTree->SetBranchStatus("Zincat",1);
-           z+=V[2];
+   fhBookTree->SetBranchStatus("Xoutcat",1);
-           l+=p;
+   fhBookTree->SetBranchStatus("Youtcat",1);
-           //cout<<x<<"  "<<y<<"  "<<z<<"  "<<l<<endl;
+   fhBookTree->SetBranchStatus("Zoutcat",1);
-           //cin>>ciao;
+   fhBookTree->SetBranchStatus("Erelcat",1);
-         }
+   fhBookTree->SetBranchStatus("Timecat",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)){
+   fhBookTree->SetBranchStatus("Pathcat",1);
-           for(t=0;t<3;t++){
+   fhBookTree->SetBranchStatus("P0cat",1);
-             if(y>=Yp[2*t] && y<Yp[2*t+1]){
+   fhBookTree->SetBranchStatus("Nthcas",1);
-               if(pmt==0)NdFT++;
+   fhBookTree->SetBranchStatus("Iparcas",1);
-               pmt=1;
+   fhBookTree->SetBranchStatus("Icas",1);
-               //cout<<NdFT<<endl;
+   fhBookTree->SetBranchStatus("Xincas",1);
-               break;
+   fhBookTree->SetBranchStatus("Yincas",1);
-             }
+   fhBookTree->SetBranchStatus("Zincas",1);
-           }
+   fhBookTree->SetBranchStatus("Xoutcas",1);
-           if(pmt==1)break;
+   fhBookTree->SetBranchStatus("Youtcas",1);
-           V[0]=-V[0];
+   fhBookTree->SetBranchStatus("Zoutcas",1);
-         }
+   fhBookTree->SetBranchStatus("Erelcas",1);
-         q=(Float_t)random()/(Float_t)0x7fffffff;
+   fhBookTree->SetBranchStatus("Timecas",1);
-         w=1-exp(-l/l0);
+   fhBookTree->SetBranchStatus("Pathcas",1);
-         if(q<w)break;
+   fhBookTree->SetBranchStatus("P0cas",1);
-         q=(Float_t)random()/(Float_t)0x7fffffff;
+   fhBookTree->SetBranchStatus("Nthcard",1);
-         w=0.5;
+   fhBookTree->SetBranchStatus("Iparcard",1);
-         if(q<w)break;
+   fhBookTree->SetBranchStatus("Icard",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[1]=-V[1];
+   fhBookTree->SetBranchStatus("Xincard",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];
+   fhBookTree->SetBranchStatus("Yincard",1);
-         x+=V[0];
+   fhBookTree->SetBranchStatus("Zincard",1);
-         y+=V[1];
+   fhBookTree->SetBranchStatus("Xoutcard",1);
-         z+=V[2];
+   fhBookTree->SetBranchStatus("Youtcard",1);
-         l=0;
+   fhBookTree->SetBranchStatus("Zoutcard",1);
-         //cout<<x<<"  "<<y<<"  "<<z<<"  "<<l<<endl;
+   fhBookTree->SetBranchStatus("Erelcard",1);
-                 //cin>>ciao;
+   fhBookTree->SetBranchStatus("Timecard",1);
-       }
+   fhBookTree->SetBranchStatus("Pathcard",1);
-     }
+   fhBookTree->SetBranchStatus("P0card",1);
-   }
-   Ert=Ert/0.002;
+   // In this simpliefied approach we will assume that once
-   q=(Float_t)(random())/(Float_t)(0x7fffffff);
+   // a particle releases > 0.5 mip in one of the 12 AC detectors it
-   w=0.7;
+   // will fire. We will furthermore assume that both cards read out
-   //E0=Float_t(4064)/7;
+   // identical data.
-   E0=4064./7.;
-   S4=(Int_t)(4064.*(1.-exp(-int(Ert)/E0)));
+   // If you develop your digitization algorithm, you should start by
-   //S4=Ert*7;
+   // identifying the information present in level2 (post-darth-vader)
-   i=S4/4;
+   // data.
-   if(S4%4==0)
-     S4v=S4+S4p;
+   Float_t SumEcat[5];
-   else if(S4%4==1) {
+   Float_t SumEcas[5];
-     if(q<w) S4v=S4-1+S4p;
+   Float_t SumEcard[5];
-     else S4v=S4+1+S4p;
+   for (Int_t k= 0;k<5;k++){
-   } else if(S4%4==2)
+     SumEcat[k]=0.;
-     S4v=S4+S4p;
+     SumEcas[k]=0.;
-   else if(S4%4==3){
+     SumEcard[k]=0.;
-     if(q<w) S4v=S4+1+S4p;
+   };
-     else S4v=S4-1+S4p;
-   }
+   if (Nthcat>50 || Nthcas>50 || Nthcard>50)
+     printf("*** ERROR AC! NthAC out of range!\n\n");
-   cout << "Ert= " <<Ert<<"; S4v= "<<S4v<<"; S4= "<<S4<<endl;
-   fDataS4[0]=S4v;//0xf028;
+   // energy dependence on position (see file AcFitOutputDistancePmt.C by S.Orsi)
-   fDataS4[1]=0xd800;
+   // based on J.Lundquist's calculations (PhD thesis, page 94)
-   fDataS4[2]=0x0300;
+   // function: [0]+[1]*atan([2]/(x+1)), where the 3 parameters are:
-   //  cout<<"  PMT  "<<NdFT<<"  "<<NdF<<endl;
+   //   8.25470e-01   +- 1.79489e-02
-   //cin>>ciao;
+   //   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!
- void Digitizer::DigitizeND(){
-   // creato:  S. Borisov, INFN Roma2 e MEPHI, Sett 2007
+   // look in CAT
-   Int_t i=0;
+   //  for (UInt_t k= 0;k<50;k++){
-   UShort_t NdN=0;
+   for (Int_t k= 0;k<Nthcat;k++){
-   fhBookTree->SetBranchStatus("Nthnd",1);
+     if (Erelcat[k] > 0)
-   fhBookTree->SetBranchStatus("Itubend",1);
+       SumEcat[Icat[k]] += Erelcat[k];
-   fhBookTree->SetBranchStatus("Iparnd",1);
+   };
-   fhBookTree->SetBranchStatus("Xinnd",1);
-   fhBookTree->SetBranchStatus("Yinnd",1);
+   // look in CAS
-   fhBookTree->SetBranchStatus("Zinnd",1);
+   for (Int_t k= 0;k<Nthcas;k++){
-   fhBookTree->SetBranchStatus("Xoutnd",1);
+     if (Erelcas[k] >0)
-   fhBookTree->SetBranchStatus("Youtnd",1);
+       SumEcas[Icas[k]] += Erelcas[k];
-   fhBookTree->SetBranchStatus("Zoutnd",1);
+   };
-   fhBookTree->SetBranchStatus("Erelnd",1);
-   fhBookTree->SetBranchStatus("Timend",1);
+   // look in CARD
-   fhBookTree->SetBranchStatus("Pathnd",1);
+   for (Int_t k= 0;k<Nthcard;k++){
-   fhBookTree->SetBranchStatus("P0nd",1);
+     if (Erelcard[k] >0)
-   //cout<<"n="<<Nthnd<<"  "<<NdN<<"\n";
+       SumEcard[Icard[k]] += Erelcard[k];
-   for(i=0;i<Nthnd;i++){
+   };
-     if(Iparnd[i]==13){
-       NdN++;
+   // channel mapping              Hit Map
-     }
+   // 1 CARD4                      0          LSB
-   }
+   // 2 CAT2                       0
-   NdN=100;
+   // 3 CAS1                       0
-   for(i=0;i<3;i++){
+   // 4 NC                         0
-     fDataND[2*i]=0x0000;
+   // 5 CARD2                      0
-     fDataND[2*i+1]=0x010F;
+   // 6 CAT4                       1
-   }
+   // 7 CAS4                       0
-   fDataND[0]=0xFF00 & (256*NdN);
+   // 8 NC                         0
- }
+   // 9 CARD3                      0
+   // 10 CAT3                      0
+   // 11 CAS3                      0
- void Digitizer::DigitizeDummy() {
+   // 12 NC                        0
+   // 13 CARD1                     0
-   fhBookTree->SetBranchStatus("Enestrip",1);
+   // 14 CAT1                      0
+   // 15 CAS2                      0
-   // dumy header
+   // 16 NC                        0          MSB
-   fDataDummy[0] = 0xCAAA;
+   // In the first version only the hit-map is filled, not the SR.
-   for (Int_t i=1; i<fDummybuffer; i++){
-     fDataDummy[i] = 0xFFFF;
+   // Threshold: 0.8 MeV.
-     //   printf("%0x  ",fDataDummy[i]);
-     //if ((i+1)%8 ==0) cout << endl;
+   Float_t thr = 8e-4;
-   }
- };
+   fDataAC[3] = 0x0000;
+   if (SumEcas[0] > thr)  fDataAC[3]  = 0x0004;
- void Digitizer::WriteData(){
+   if (SumEcas[1] > thr)  fDataAC[3] += 0x4000;
+   if (SumEcas[2] > thr)  fDataAC[3] += 0x0400;
-   // Routine that writes the data to a binary file
+   if (SumEcas[3] > thr)  fDataAC[3] += 0x0040;
-   // PSCU data are already swapped
-   fOutputfile.write(reinterpret_cast<char*>(fDataPSCU),sizeof(UShort_t)*fPSCUbuffer);
+   if (SumEcat[0] > thr)  fDataAC[3] += 0x2000;
-   // TRG
+   if (SumEcat[1] > thr)  fDataAC[3] += 0x0002;
-   fOutputfile.write(reinterpret_cast<char*>(fDataTrigger),sizeof(UChar_t)*153);
+   if (SumEcat[2] > thr)  fDataAC[3] += 0x0200;
-   // TOF
+   if (SumEcat[3] > thr)  fDataAC[3] += 0x0020;
-   fOutputfile.write(reinterpret_cast<char*>(fDataTof),sizeof(UChar_t)*276);
-   // AC
+   if (SumEcard[0] > thr)  fDataAC[3] += 0x1000;
-   UShort_t temp[1000000];
+   if (SumEcard[1] > thr)  fDataAC[3] += 0x0010;
-   memset(temp,0,sizeof(UShort_t)*1000000);
+   if (SumEcard[2] > thr)  fDataAC[3] += 0x0100;
-   swab(fDataAC,temp,sizeof(UShort_t)*fACbuffer);  // WE MUST SWAP THE BYTES!!!
+   if (SumEcard[3] > thr)  fDataAC[3] += 0x0001;
-   fOutputfile.write(reinterpret_cast<char*>(temp),sizeof(UShort_t)*fACbuffer);
-   // CALO
+   fDataAC[67] = fDataAC[3];
-   memset(temp,0,sizeof(UShort_t)*1000000);
-   swab(fDataCALO,temp,sizeof(UShort_t)*fCALOlength); // WE MUST SWAP THE BYTES!!!
+   // shift registers
-   fOutputfile.write(reinterpret_cast<char*>(temp),sizeof(UShort_t)*fCALOlength);
+   // the central bin is equal to the hitmap, all other bins in the shift register are 0
-   // TRK
+   for (UInt_t i=0; i<=15; i++){
-   memset(temp,0,sizeof(UShort_t)*1000000);
+     fDataAC[i+11] = 0x0000;
-   swab(fDataTrack,temp,sizeof(UShort_t)*fTracklength);  // WE MUST SWAP THE BYTES!!!
+     fDataAC[i+75] = 0x0000;
-   fOutputfile.write(reinterpret_cast<char*>(temp),sizeof(UShort_t)*fTracklength);
+   }
-   fTracklength=0;
+   fDataAC[18] = fDataAC[3];
-    // padding to 64 bytes
+   fDataAC[82] = fDataAC[3];
-   //
-   if ( fPadding ){
+   //    for (Int_t i=0; i<fACbuffer; i++){
-     fOutputfile.write(reinterpret_cast<char*>(fDataPadding),sizeof(UChar_t)*fPadding);
+   //    printf("%0x  ",fDataAC[i]);
-   };
+   //    if ((i+1)%8 ==0) cout << endl;
-   // 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
+ void Digitizer::DigitizeS4(){
-   memset(temp,0,sizeof(UShort_t)*1000000);
+   Int_t DEBUG=0;
-   swab(fDataND,temp,sizeof(UShort_t)*fNDbuffer);  // WE MUST SWAP THE BYTES!!!
+   // creato:  S. Borisov, INFN Roma2 e MEPHI, Sett 2007
-   fOutputfile.write(reinterpret_cast<char*>(temp),sizeof(UShort_t)*fNDbuffer);
+   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;
- void Digitizer::ReadData(){
+   Xs[1]=24.1;
+   Ys[0]=-24.1;
-   UShort_t InData[64];
+   Ys[1]=24.1;
+   Zs[0]=-0.5;
-   // for debuggigng purposes only, write your own routine if you like (many
+   Zs[1]=0.5;
-   // hardwired things.
+   Yp[0]=-20.;
+   Yp[2]=-1.;
-   ifstream InputFile;
+   Yp[4]=17.;
+   for(i=0;i<3;i++)
-   // if (!InputFile) {
+     Yp[2*i+1]=Yp[2*i]+3;
+   srand(time(NULL));
-   //     std::cout << "ERROR" << endl;
+   // --- activate branches:
-   //     // An error occurred!
+   fhBookTree->SetBranchStatus("Nthtof",1);
-   //     // myFile.gcount() returns the number of bytes read.
+   fhBookTree->SetBranchStatus("Ipltof",1);
-   //     // calling myFile.clear() will reset the stream state
+   fhBookTree->SetBranchStatus("Ipaddle",1);
-   //     // so it is usable again.
-   //   };
+   fhBookTree->SetBranchStatus("Xintof",1);
+   fhBookTree->SetBranchStatus("Yintof",1);
+   fhBookTree->SetBranchStatus("Xouttof",1);
+   fhBookTree->SetBranchStatus("Youttof",1);
-   //InputFile.seekg(0);
+   fhBookTree->SetBranchStatus("Ereltof",1);
-   InputFile.open(fFilename, ios::in | ios::binary);
+   fhBookTree->SetBranchStatus("Timetof",1);
-   //    fOutputfile.seekg(0);
+   NdFT=0;
-   if (!InputFile.is_open()) std::cout << "ERROR" << endl;
+   Ert=0;
+   for(i=0;i<Nthtof;i++){
-   InputFile.seekg(0);
+     if(Ipltof[i]!=6) continue;
+     Ert+=Ereltof[i];
-   for (Int_t k=0; k<=1000; k++){
-     InputFile.read(reinterpret_cast<char*>(InData),384*sizeof(UShort_t));
+     if(modo=="ns") continue;
+     NdF=Int_t(Ereltof[i]/E1);
-     std::cout << "Read back: " << endl << endl;
+     NdFT=0;
+     X=Xintof[i];
-     for (Int_t i=0; i<=384; i++){
+     Y=Yintof[i];
-       printf("%4x ", InData[i]);
+     Z=(Float_t)(random())/(Float_t)(0x7fffffff)-0.5;
-       if ((i+1)%8 ==0) cout << endl;
+     //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 << endl;
+       // cout<<"qw "<<q<<" "<<w<<endl;
-   InputFile.close();
+       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;
- void Digitizer::DigitizeTrack() {
+       z=Z;
- //std:: cout << "Entering DigitizeTrack " << endl;
+       while(pmt==0 && (x>Xs[0] && x<Xs[1])&&(y>Ys[0] && y<Ys[1])&&(z>Zs[0] && z<Zs[1])){
- Float_t  AdcTrack[fNviews][fNstrips_view];  //  Vector of strips to be compressed
+         l=0;
+         while(pmt==0 && (x>Xs[0] && x<Xs[1])&&(y>Ys[0] && y<Ys[1])&&(z>Zs[0] && z<Zs[1])){
- Int_t Iview;
+           x+=V[0];
- Int_t Nstrip;
+           y+=V[1];
+           z+=V[2];
-   for (Int_t j=0; j<fNviews;j++) {
+           l+=p;
+           //cout<<x<<"  "<<y<<"  "<<z<<"  "<<l<<endl;
-     for (Int_t i=0; i<fNladder;i++) {
+           //cin>>ciao;
+         }
-       Float_t commonN1=gRandom->Gaus(0.,fSigmaCommon);
+         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)){
-       Float_t commonN2=gRandom->Gaus(0.,fSigmaCommon);
+           for(t=0;t<3;t++){
-       for (Int_t k=0; k<fNstrips_ladder;k++) {
+             if(y>=Yp[2*t] && y<Yp[2*t+1]){
-       Nstrip=i*fNstrips_ladder+k;
+               if(pmt==0)NdFT++;
-       AdcTrack[j][Nstrip]=gRandom->Gaus(fPedeTrack[j][Nstrip],fSigmaTrack[j][Nstrip]);
+               pmt=1;
-       if(k<4*128) {AdcTrack[j][Nstrip] += commonN1;}  // full correlation of 4 VA1 Com. Noise
+               //cout<<NdFT<<endl;
-       else {AdcTrack[j][Nstrip] += commonN2;}   // full correlation of 4 VA1 Com. Noise
+               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;
-   fhBookTree->SetBranchStatus("Nstrpx",1);
+         if(q<w)break;
-   fhBookTree->SetBranchStatus("Npstripx",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[1]=-V[1];
-   fhBookTree->SetBranchStatus("Ntstripx",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];
-   fhBookTree->SetBranchStatus("Istripx",1);
+         x+=V[0];
-   fhBookTree->SetBranchStatus("Qstripx",1);
+         y+=V[1];
-   fhBookTree->SetBranchStatus("Xstripx",1);
+         z+=V[2];
-   fhBookTree->SetBranchStatus("Nstrpy",1);
+         l=0;
-   fhBookTree->SetBranchStatus("Npstripy",1);
+         //cout<<x<<"  "<<y<<"  "<<z<<"  "<<l<<endl;
-   fhBookTree->SetBranchStatus("Ntstripy",1);
+                 //cin>>ciao;
-   fhBookTree->SetBranchStatus("Istripy",1);
+       }
-   fhBookTree->SetBranchStatus("Qstripy",1);
+     }
-   fhBookTree->SetBranchStatus("Ystripy",1);
+   }
+   Ert=Ert/0.002;
+   q=(Float_t)(random())/(Float_t)0x7fffffff;
+   w=0.7;
-   Float_t ADCfull;
+   //E0=(Float_t)(4064./7.);
-   Int_t iladd=0;
+   E0=4064./7.;
-   for (Int_t ix=0; ix<Nstrpx;ix++) {
+   if(Ert<1) S4=0;
-   Iview=Npstripx[ix]*2-1;
+   else S4=(Int_t)(4064.*(1.-exp(-(Ert-1.)/E0)));
-   Nstrip=(Int_t)Istripx[ix]-1;
+   i=S4/4;
-   if(Nstrip<fNstrips_ladder) iladd=0;
+   if(S4%4==0)
-   if((Nstrip>=fNstrips_ladder)&&(Nstrip<2*fNstrips_ladder)) iladd=1;
+     S4v=S4+S4p;
-   if((Nstrip>=2*fNstrips_ladder)&&(Nstrip<3*fNstrips_ladder)) iladd=2;
+   else if(S4%4==1){
-   ADCfull=AdcTrack[Iview][Nstrip] += Qstripx[ix]*fMipCor[iladd][Iview];
+     if(q<w) S4v=S4-1+S4p;
-   AdcTrack[Iview][Nstrip] *= SaturationTrack(ADCfull);
+     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;
-   for (Int_t iy=0; iy<Nstrpy;iy++) {
+   }
-   Iview=Npstripy[iy]*2-2;
+   if (DEBUG)
-   Nstrip=(Int_t)Istripy[iy]-1;
+     cout<<"Ert_S4 = " << Ert << " --- S4v = " << S4v << endl;
-   if(Nstrip<fNstrips_ladder) iladd=0;
+   fDataS4[0]=S4v;//0xf028;
-   if((Nstrip>=fNstrips_ladder)&&(Nstrip<2*fNstrips_ladder)) iladd=1;
+   fDataS4[1]=0xd800;
-   if((Nstrip>=2*fNstrips_ladder)&&(Nstrip<3*fNstrips_ladder)) iladd=2;
+   fDataS4[2]=0x0300;
-   ADCfull=AdcTrack[Iview][Nstrip] -= Qstripy[iy]*fMipCor[iladd][Iview];
+   //cout<<"  PMT  "<<NdFT<<"  "<<NdF<<endl;
-   AdcTrack[Iview][Nstrip] *= SaturationTrack(ADCfull);
+   //cin>>ciao;
+ }
-   };
- CompressTrackData(AdcTrack);  // Compress and Digitize data of one Ladder  in turn for all ladders
+ void Digitizer::DigitizeND(){
- };
+   // creato:  S. Borisov, INFN Roma2 e MEPHI, Sett 2007
+   Int_t i=0;
+   UShort_t NdN=0;
+   fhBookTree->SetBranchStatus("Nthnd",1);
- void Digitizer::DigitizeTrackCalib(Int_t ii) {
+   fhBookTree->SetBranchStatus("Itubend",1);
+   fhBookTree->SetBranchStatus("Iparnd",1);
- std:: cout << "Entering DigitizeTrackCalib " << ii << endl;
+   fhBookTree->SetBranchStatus("Xinnd",1);
- if( (ii!=1)&&(ii!=2) ) {
+   fhBookTree->SetBranchStatus("Yinnd",1);
-  std:: cout << "error wrong DigitizeTrackCalib argument" << endl;
+   fhBookTree->SetBranchStatus("Zinnd",1);
-  return;
+   fhBookTree->SetBranchStatus("Xoutnd",1);
- };
+   fhBookTree->SetBranchStatus("Youtnd",1);
+   fhBookTree->SetBranchStatus("Zoutnd",1);
- memset(fDataTrack,0,sizeof(UShort_t)*fTRACKbuffer);
+   fhBookTree->SetBranchStatus("Erelnd",1);
- fTracklength=0;
+   fhBookTree->SetBranchStatus("Timend",1);
+   fhBookTree->SetBranchStatus("Pathnd",1);
- UShort_t Dato;
+   fhBookTree->SetBranchStatus("P0nd",1);
+   //cout<<"n="<<Nthnd<<"  "<<NdN<<"\n";
- Float_t dato1;
+   for(i=0;i<Nthnd;i++){
- Float_t dato2;
+     if(Iparnd[i]==13){
- Float_t dato3;
+       NdN++;
- Float_t dato4;
+     }
+   }
- UShort_t DatoDec;
+   //NdN=100; //only for debug
- UShort_t DatoDec1;
- UShort_t DatoDec2;
+   for(i=0;i<3;i++){
- UShort_t DatoDec3;
+     fDataND[2*i]=0x0000;
- UShort_t DatoDec4;
+     fDataND[2*i+1]=0x010F;
+   }
- UShort_t EVENT_CAL;
+   fDataND[0]=0xFF00 & (256*NdN);
- UShort_t PED_L1;
+ }
- UShort_t ReLength;
- UShort_t OveCheckCode;
- //UShort_t PED_L2;
+ void Digitizer::DigitizeDummy() {
- //UShort_t PED_L3HI;
- //UShort_t PED_L3LO;
+   fhBookTree->SetBranchStatus("Enestrip",1);
- //UShort_t SIG_L1HI;
- //UShort_t SIG_L1LO;
+   // dumy header
- //UShort_t SIG_L2HI;
+   fDataDummy[0] = 0xCAAA;
- //UShort_t SIG_L2LO;
- //UShort_t SIG_L3;
+   for (Int_t i=1; i<fDummybuffer; i++){
- //UShort_t BAD_L1;
+     fDataDummy[i] = 0xFFFF;
- //UShort_t BAD_L2LO;
+     //   printf("%0x  ",fDataDummy[i]);
- //UShort_t BAD_L3HI;
+     //if ((i+1)%8 ==0) cout << endl;
- //UShort_t BAD_L3LO;
+   }
- //UShort_t FLAG;
+ };
-   Int_t DSPpos;
+ void Digitizer::WriteData(){
-   for (Int_t j=ii-1; j<fNviews;j+=2) {
-   UShort_t CkSum=0;
+   // Routine that writes the data to a binary file
-   // here skip the dsp header and his trailer , to be written later
+   // PSCU data are already swapped
-   DSPpos=fTracklength;
+   fOutputfile.write(reinterpret_cast<char*>(fDataPSCU),sizeof(UShort_t)*fPSCUbuffer);
-   fTracklength=fTracklength+13+3;
+   // TRG
+   fOutputfile.write(reinterpret_cast<char*>(fDataTrigger),sizeof(UChar_t)*153);
+   // TOF
-     for (Int_t i=0; i<fNladder;i++) {
+   fOutputfile.write(reinterpret_cast<char*>(fDataTof),sizeof(UChar_t)*276);
-       for (Int_t k=0; k<fNstrips_ladder;k++) {
+   // AC
-       // write in buffer the current LADDER
+   UShort_t temp[1000000];
-       Dato=(UShort_t)fPedeTrack[j][i*fNstrips_ladder+k];
+   memset(temp,0,sizeof(UShort_t)*1000000);
-       dato1=fPedeTrack[j][i*fNstrips_ladder+k]-Dato;
+   swab(fDataAC,temp,sizeof(UShort_t)*fACbuffer);  // WE MUST SWAP THE BYTES!!!
+   fOutputfile.write(reinterpret_cast<char*>(temp),sizeof(UShort_t)*fACbuffer);
-       DatoDec1=(UShort_t)(dato1*2);
+   // CALO
-       dato2=dato1*2-DatoDec1;
+   memset(temp,0,sizeof(UShort_t)*1000000);
+   swab(fDataCALO,temp,sizeof(UShort_t)*fCALOlength); // WE MUST SWAP THE BYTES!!!
-       DatoDec2=(UShort_t)(dato2*2);
+   fOutputfile.write(reinterpret_cast<char*>(temp),sizeof(UShort_t)*fCALOlength);
-       dato3=dato2*2-DatoDec2;
+   // TRK
+   memset(temp,0,sizeof(UShort_t)*1000000);
-       DatoDec3=(UShort_t)(dato3*2);
+   swab(fDataTrack,temp,sizeof(UShort_t)*fTracklength);  // WE MUST SWAP THE BYTES!!!
-       dato4=dato3*2-DatoDec3;
+   fOutputfile.write(reinterpret_cast<char*>(temp),sizeof(UShort_t)*fTracklength);
+   fTracklength=0;
-       DatoDec4=(UShort_t)(dato4*2);
+    // padding to 64 bytes
+   //
-       DatoDec=DatoDec1*0x0008+DatoDec2*0x0004+DatoDec3*0x0002+DatoDec4*0x0001;
+   if ( fPadding ){
-       fDataTrack[fTracklength]=( (Dato << 4) | (DatoDec & 0x000F) );
+     fOutputfile.write(reinterpret_cast<char*>(fDataPadding),sizeof(UChar_t)*fPadding);
-       CkSum=CkSum^fDataTrack[fTracklength];
+   };
-       fTracklength++;
+   // S4
-       };
+   memset(temp,0,sizeof(UShort_t)*1000000);
+   swab(fDataS4,temp,sizeof(UShort_t)*fS4buffer);  // WE MUST SWAP THE BYTES!!!
-       for (Int_t k=0; k<fNstrips_ladder;k++) {
+   fOutputfile.write(reinterpret_cast<char*>(temp),sizeof(UShort_t)*fS4buffer);
-       // write in buffer the current LADDER
+   // ND
-       Dato=(UShort_t)fSigmaTrack[j][i*fNstrips_ladder+k];
+   memset(temp,0,sizeof(UShort_t)*1000000);
-       dato1=fSigmaTrack[j][i*fNstrips_ladder+k]-Dato;
+   swab(fDataND,temp,sizeof(UShort_t)*fNDbuffer);  // WE MUST SWAP THE BYTES!!!
+   fOutputfile.write(reinterpret_cast<char*>(temp),sizeof(UShort_t)*fNDbuffer);
-       DatoDec1=(UShort_t)(dato1*2);
+ };
-       dato2=dato1*2-DatoDec1;
-       DatoDec2=(UShort_t)(dato2*2);
+ void Digitizer::ReadData(){
-       dato3=dato2*2-DatoDec2;
+   UShort_t InData[64];
-       DatoDec3=(UShort_t)(dato3*2);
-       dato4=dato3*2-DatoDec3;
+   // for debuggigng purposes only, write your own routine if you like (many
+   // hardwired things.
-       DatoDec4=(UShort_t)(dato4*2);
+   ifstream InputFile;
-       DatoDec=DatoDec1*0x0008+DatoDec2*0x0004+DatoDec3*0x0002+DatoDec4*0x0001;
+   // if (!InputFile) {
-       fDataTrack[fTracklength]=( (Dato << 4) | (DatoDec & 0x000F) );
-       CkSum=CkSum^fDataTrack[fTracklength];
+   //     std::cout << "ERROR" << endl;
-       fTracklength++;
+   //     // An error occurred!
-       };
+   //     // myFile.gcount() returns the number of bytes read.
+   //     // calling myFile.clear() will reset the stream state
-       for (Int_t k=0; k<64;k++) {
+   //     // so it is usable again.
-       fDataTrack[fTracklength]=0x0000;
+   //   };
-       CkSum=CkSum^fDataTrack[fTracklength];
-       fTracklength++;
-       };
+   //InputFile.seekg(0);
-       // end ladder
+   InputFile.open(fFilename, ios::in | ios::binary);
-     // write in buffer the end ladder word
+   //    fOutputfile.seekg(0);
-     if(i==0) fDataTrack[fTracklength]=0x1807;
+   if (!InputFile.is_open()) std::cout << "ERROR" << endl;
-     if(i==1) fDataTrack[fTracklength]=0x1808;
-     if(i==2) fDataTrack[fTracklength]=0x1809;
+   InputFile.seekg(0);
-     CkSum=CkSum^fDataTrack[fTracklength];
-     fTracklength++;
+   for (Int_t k=0; k<=1000; k++){
+     InputFile.read(reinterpret_cast<char*>(InData),384*sizeof(UShort_t));
-     // write in buffer the TRAILER
-     ReLength=(UShort_t)((fNstrips_ladder*2+64+1)*2+3);
+     std::cout << "Read back: " << endl << endl;
-     OveCheckCode=0x0000;
+     for (Int_t i=0; i<=384; i++){
-     fDataTrack[fTracklength]=0x0000;
+       printf("%4x ", InData[i]);
-     fTracklength++;
+       if ((i+1)%8 ==0) cout << endl;
+     }
-     fDataTrack[fTracklength]=(ReLength >> 8);
-     fTracklength++;
+   }
+   cout << endl;
-     fDataTrack[fTracklength]=( (ReLength << 8) | (OveCheckCode & 0x00FF) );
+   InputFile.close();
-     fTracklength++;
+ };
-     // end TRAILER
-     };
-   // write in buffer the DSP header
+ void Digitizer::DigitizeTrack() {
+ //std:: cout << "Entering DigitizeTrack " << endl;
-   fDataTrack[DSPpos]=(0xE800 | ( ((j+1) << 3) | 0x0005) );
+ Float_t  AdcTrack[fNviews][fNstrips_view];  //  Vector of strips to be compressed
-   fDataTrack[DSPpos+1]=0x01A9;
+ Int_t Iview;
+ Int_t Nstrip;
-   fDataTrack[DSPpos+2]=0x8740;
+   for (Int_t j=0; j<fNviews;j++) {
-   EVENT_CAL=0;
-   fDataTrack[DSPpos+3]=(0x1A00 | ( (0x03FF & EVENT_CAL)>> 1) );
+     for (Int_t i=0; i<fNladder;i++) {
-   PED_L1=0;
+       Float_t commonN1=gRandom->Gaus(0.,fSigmaCommon);
-   fDataTrack[DSPpos+4]=( ((EVENT_CAL << 15) | 0x5002 ) | ((0x03FF & PED_L1) << 2) );
+       Float_t commonN2=gRandom->Gaus(0.,fSigmaCommon);
+       for (Int_t k=0; k<fNstrips_ladder;k++) {
-   // FROM HERE WE WRITE AS ALL VARIABLE apart CkSum are =0
+       Nstrip=i*fNstrips_ladder+k;
+       AdcTrack[j][Nstrip]=gRandom->Gaus(fPedeTrack[j][Nstrip],fSigmaTrack[j][Nstrip]);
-   fDataTrack[DSPpos+5]=0x8014;
+       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
-   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);
+   fhBookTree->SetBranchStatus("Nstrpx",1);
-   fDataTrack[DSPpos+11]=(0x0280 | (CkSum >> 3));
+   fhBookTree->SetBranchStatus("Npstripx",1);
+   fhBookTree->SetBranchStatus("Ntstripx",1);
-   fDataTrack[DSPpos+12]=(0x1FFF | (CkSum << 13) );
+   fhBookTree->SetBranchStatus("Istripx",1);
+   fhBookTree->SetBranchStatus("Qstripx",1);
-   // end dsp header
+   fhBookTree->SetBranchStatus("Xstripx",1);
+   fhBookTree->SetBranchStatus("Nstrpy",1);
-   // write in buffer the TRAILER
+   fhBookTree->SetBranchStatus("Npstripy",1);
+   fhBookTree->SetBranchStatus("Ntstripy",1);
-   ReLength=(UShort_t)((13*2)+3);
+   fhBookTree->SetBranchStatus("Istripy",1);
-   OveCheckCode=0x0000;
+   fhBookTree->SetBranchStatus("Qstripy",1);
-   fDataTrack[DSPpos+13]=0x0000;
+   fhBookTree->SetBranchStatus("Ystripy",1);
-   fDataTrack[DSPpos+14]=(ReLength >> 8);
-   fDataTrack[DSPpos+15]=( (ReLength << 8) | (OveCheckCode & 0x00FF) );
+   Float_t ADCfull;
+   Int_t iladd=0;
-   // end TRAILER
+   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;
-   // end DSP
+   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);
+   };
- };
- void Digitizer::WriteTrackCalib() {
+   for (Int_t iy=0; iy<Nstrpy;iy++) {
+   Iview=Npstripy[iy]*2-2;
+   Nstrip=(Int_t)Istripy[iy]-1;
- std:: cout << " Entering WriteTrackCalib " << endl;
+   if(Nstrip<fNstrips_ladder) iladd=0;
+   if((Nstrip>=fNstrips_ladder)&&(Nstrip<2*fNstrips_ladder)) iladd=1;
- fOutputfile.write(reinterpret_cast<char*>(fDataPSCU),sizeof(UShort_t)*fPSCUbuffer);
+   if((Nstrip>=2*fNstrips_ladder)&&(Nstrip<3*fNstrips_ladder)) iladd=2;
+   ADCfull=AdcTrack[Iview][Nstrip] -= Qstripy[iy]*fMipCor[iladd][Iview];
- UShort_t temp[1000000];
+   AdcTrack[Iview][Nstrip] *= SaturationTrack(ADCfull);
- 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;
+ CompressTrackData(AdcTrack);  // Compress and Digitize data of one Ladder  in turn for all ladders
- if ( fPadding ){
-       fOutputfile.write(reinterpret_cast<char*>(fDataPadding),sizeof(UChar_t)*fPadding);
+ };
- };
- };
+ void Digitizer::DigitizeTrackCalib(Int_t ii) {
- void Digitizer::ClearTrackCalib() {
+ std:: cout << "Entering DigitizeTrackCalib " << ii << endl;
+ if( (ii!=1)&&(ii!=2) ) {
- std:: cout << "Entering ClearTrackCalib " << endl;
+  std:: cout << "error wrong DigitizeTrackCalib argument" << endl;
+  return;
+ };
- };
+ memset(fDataTrack,0,sizeof(UShort_t)*fTRACKbuffer);
+ fTracklength=0;
- void Digitizer::LoadTrackCalib() {
- std:: cout << "Entering LoadTrackCalib " << endl;
+ UShort_t Dato;
- // Generate the pedestals and sigmas according to parametrization
+ Float_t dato1;
-   for (Int_t j=0; j<fNviews;j++) {
+ Float_t dato2;
-     for (Int_t i=0; i<fNstrips_view;i++) {
+ Float_t dato3;
+ Float_t dato4;
-     if((j+1)%2==0) {
-     fPedeTrack[j][i]=gRandom->Gaus(fAvePedex,fSigmaPedex);
+ UShort_t DatoDec;
-     fSigmaTrack[j][i]=gRandom->Gaus(fAveSigmax,fSigmaSigmax);
+ UShort_t DatoDec1;
-     };
+ UShort_t DatoDec2;
-     if((j+1)%2==1) {
+ UShort_t DatoDec3;
-     fPedeTrack[j][i]=gRandom->Gaus(fAvePedey,fSigmaPedey);
+ UShort_t DatoDec4;
-     fSigmaTrack[j][i]=gRandom->Gaus(fAveSigmay,fSigmaSigmay);
-     };
+ 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;
- void Digitizer::LoadMipCor() {
+ //UShort_t SIG_L2HI;
- std:: cout << "Entering LoadMipCor" << endl;
+ //UShort_t SIG_L2LO;
-   Float_t xfactor=1./151.6*1.04;
+ //UShort_t SIG_L3;
-   Float_t yfactor=1./152.1;
+ //UShort_t BAD_L1;
+ //UShort_t BAD_L2LO;
-   fMipCor[0][0]=140.02*yfactor;
+ //UShort_t BAD_L3HI;
-   fMipCor[0][1]=140.99*xfactor;
+ //UShort_t BAD_L3LO;
-   fMipCor[0][2]=134.48*yfactor;
+ //UShort_t FLAG;
-   fMipCor[0][3]=144.41*xfactor;
-   fMipCor[0][4]=140.74*yfactor;
-   fMipCor[0][5]=142.28*xfactor;
+   Int_t DSPpos;
-   fMipCor[0][6]=134.53*yfactor;
+   for (Int_t j=ii-1; j<fNviews;j+=2) {
-   fMipCor[0][7]=140.63*xfactor;
+   UShort_t CkSum=0;
-   fMipCor[0][8]=135.55*yfactor;
+   // here skip the dsp header and his trailer , to be written later
-   fMipCor[0][9]=138.00*xfactor;
+   DSPpos=fTracklength;
-   fMipCor[0][10]=154.95*yfactor;
+   fTracklength=fTracklength+13+3;
-   fMipCor[0][11]=158.44*xfactor;
+     for (Int_t i=0; i<fNladder;i++) {
-   fMipCor[1][0]=136.07*yfactor;
+       for (Int_t k=0; k<fNstrips_ladder;k++) {
-   fMipCor[1][1]=135.59*xfactor;
+       // write in buffer the current LADDER
-   fMipCor[1][2]=142.69*yfactor;
+       Dato=(UShort_t)fPedeTrack[j][i*fNstrips_ladder+k];
-   fMipCor[1][3]=138.19*xfactor;
+       dato1=fPedeTrack[j][i*fNstrips_ladder+k]-Dato;
-   fMipCor[1][4]=137.35*yfactor;
-   fMipCor[1][5]=140.23*xfactor;
+       DatoDec1=(UShort_t)(dato1*2);
-   fMipCor[1][6]=153.15*yfactor;
+       dato2=dato1*2-DatoDec1;
-   fMipCor[1][7]=151.42*xfactor;
-   fMipCor[1][8]=129.76*yfactor;
+       DatoDec2=(UShort_t)(dato2*2);
-   fMipCor[1][9]=140.63*xfactor;
+       dato3=dato2*2-DatoDec2;
-   fMipCor[1][10]=157.87*yfactor;
-   fMipCor[1][11]=153.64*xfactor;
+       DatoDec3=(UShort_t)(dato3*2);
+       dato4=dato3*2-DatoDec3;
-   fMipCor[2][0]=134.98*yfactor;
-   fMipCor[2][1]=143.95*xfactor;
+       DatoDec4=(UShort_t)(dato4*2);
-   fMipCor[2][2]=140.23*yfactor;
-   fMipCor[2][3]=138.88*xfactor;
+       DatoDec=DatoDec1*0x0008+DatoDec2*0x0004+DatoDec3*0x0002+DatoDec4*0x0001;
-   fMipCor[2][4]=137.95*yfactor;
+       fDataTrack[fTracklength]=( (Dato << 4) | (DatoDec & 0x000F) );
-   fMipCor[2][5]=134.87*xfactor;
+       CkSum=CkSum^fDataTrack[fTracklength];
-   fMipCor[2][6]=157.56*yfactor;
+       fTracklength++;
-   fMipCor[2][7]=157.31*xfactor;
+       };
-   fMipCor[2][8]=141.37*yfactor;
-   fMipCor[2][9]=143.39*xfactor;
+       for (Int_t k=0; k<fNstrips_ladder;k++) {
-   fMipCor[2][10]=156.15*yfactor;
+       // write in buffer the current LADDER
-   fMipCor[2][11]=158.79*xfactor;
+       Dato=(UShort_t)fSigmaTrack[j][i*fNstrips_ladder+k];
+       dato1=fSigmaTrack[j][i*fNstrips_ladder+k]-Dato;
- /*
-   for (Int_t j=0; j<fNviews;j++) {
+       DatoDec1=(UShort_t)(dato1*2);
-     for (Int_t i=0; i<fNstrips_view;i++) {
+       dato2=dato1*2-DatoDec1;
-     fMipCor[j][i]=1.;
-     };
+       DatoDec2=(UShort_t)(dato2*2);
-   };
+       dato3=dato2*2-DatoDec2;
+       DatoDec3=(UShort_t)(dato3*2);
- */
+       dato4=dato3*2-DatoDec3;
- };
+       DatoDec4=(UShort_t)(dato4*2);
- void Digitizer::CompressTrackData(Float_t AdcTrack[fNviews][fNstrips_view]) {
- // copy of the corresponding compression fortran routine + new digitization
+       DatoDec=DatoDec1*0x0008+DatoDec2*0x0004+DatoDec3*0x0002+DatoDec4*0x0001;
- // std:: cout << "Entering CompressTrackData " << endl;
- Int_t oldval=0;
+       fDataTrack[fTracklength]=( (Dato << 4) | (DatoDec & 0x000F) );
- Int_t newval=0;
+       CkSum=CkSum^fDataTrack[fTracklength];
- Int_t trasmesso=0;
+       fTracklength++;
- Int_t ntrastot=0;
+       };
- Float_t real;
- Float_t inte;
+       for (Int_t k=0; k<64;k++) {
- Int_t cercacluster=0;
+       fDataTrack[fTracklength]=0x0000;
- Int_t kt=0;
+       CkSum=CkSum^fDataTrack[fTracklength];
- static const int DSPbufferSize = 4000; // 13 bit buffer to be rearranged in 16 bit Track buffer
+       fTracklength++;
- 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
+       };
+       // end ladder
- memset(fDataTrack,0,sizeof(UShort_t)*fTRACKbuffer); // probably not necessary becouse already done ?
- fTracklength=0;
+     // write in buffer the end ladder word
+     if(i==0) fDataTrack[fTracklength]=0x1807;
-   for (Int_t iv=0; iv<fNviews;iv++) {
+     if(i==1) fDataTrack[fTracklength]=0x1808;
-   memset(DataDSP,0,sizeof(UShort_t)*DSPbufferSize);
+     if(i==2) fDataTrack[fTracklength]=0x1809;
-   DSPlength=16; // skip the header, to be written later
+     CkSum=CkSum^fDataTrack[fTracklength];
-   UShort_t CheckSum=0;
+     fTracklength++;
- // write dsp header on buffer
+     // write in buffer the TRAILER
- //    fDataTrack[fTracklength]=0xE805;
+     ReLength=(UShort_t)((fNstrips_ladder*2+64+1)*2+3);
- //    fTracklength++;
+     OveCheckCode=0x0000;
- //    fDataTrack[fTracklength]=0x01A9;
+     fDataTrack[fTracklength]=0x0000;
- //    fTracklength++;
+     fTracklength++;
- // end dsp header
+     fDataTrack[fTracklength]=(ReLength >> 8);
+     fTracklength++;
-    //
-    // INIZIO VISTA IV - TAKE PROPER ACTION
+     fDataTrack[fTracklength]=( (ReLength << 8) | (OveCheckCode & 0x00FF) );
-    //
+     fTracklength++;
+     // end TRAILER
+     };
-     for (Int_t ladder=0; ladder<fNladder;ladder++) {
-       Int_t k=0;
+   // write in buffer the DSP header
-       while (k<fNstrips_ladder) {
-         // compress write in buffer the current LADDER
+   fDataTrack[DSPpos]=(0xE800 | ( ((j+1) << 3) | 0x0005) );
-         if ( k == 0)  {
-           real=modff(AdcTrack[iv][ladder*fNstrips_ladder+k],&inte);
+   fDataTrack[DSPpos+1]=0x01A9;
-           if (real > 0.5) inte=inte+1;
-           newval=(Int_t)inte -(Int_t)fPedeTrack[iv][ladder*fNstrips_ladder+k];
+   fDataTrack[DSPpos+2]=0x8740;
-           // first strip of ladder is transmitted
-           // DC_TOT first " << AdcTrack[iv][ladder*fNstrips_ladder+k] << endl;
+   EVENT_CAL=0;
-           DataDSP[DSPlength]=( ((UShort_t)inte) & 0x0FFF);
+   fDataTrack[DSPpos+3]=(0x1A00 | ( (0x03FF & EVENT_CAL)>> 1) );
-           DSPlength++;
-           ntrastot++;
+   PED_L1=0;
-           trasmesso=1;
+   fDataTrack[DSPpos+4]=( ((EVENT_CAL << 15) | 0x5002 ) | ((0x03FF & PED_L1) << 2) );
-           oldval=newval;
-           kt=k;
+   // FROM HERE WE WRITE AS ALL VARIABLE apart CkSum are =0
-           k++;
-           continue;
+   fDataTrack[DSPpos+5]=0x8014;
-         };
-         real=modff(AdcTrack[iv][ladder*fNstrips_ladder+k],&inte);
+   fDataTrack[DSPpos+6]=0x00A0;
-         if (real > 0.5) inte=inte+1;
-         newval=(Int_t)inte -(Int_t)(fPedeTrack[iv][ladder*fNstrips_ladder+k]);
+   fDataTrack[DSPpos+7]=0x0500;
-         cercacluster=1; // ?????????
-         if (cercacluster==1) {
+   fDataTrack[DSPpos+8]=0x2801;
-  // controlla l'ordine di tutti queste strip ladder e DSP !!!!!!!
+   fDataTrack[DSPpos+9]=0x400A;
-           Int_t diff=0;
+   fDataTrack[DSPpos+10]=0x0050;
-           switch ((iv+1)%2) {
+   CkSum=(CkSum >> 8)^(CkSum&0x00FF);
-           case 0: diff=newval-oldval;
+   fDataTrack[DSPpos+11]=(0x0280 | (CkSum >> 3));
-           break;
-           case 1: diff=oldval-newval;
+   fDataTrack[DSPpos+12]=(0x1FFF | (CkSum << 13) );
-           break;
-           };
+   // end dsp header
-           if (diff>fCutclu*(Int_t)fSigmaTrack[iv][ladder*fNstrips_ladder+k]) {
+   // write in buffer the TRAILER
-             Int_t clval=newval;
-             Int_t klp=k;        // go on to search for maximum
+   ReLength=(UShort_t)((13*2)+3);
-             klp++;
+   OveCheckCode=0x0000;
+   fDataTrack[DSPpos+13]=0x0000;
-             while(klp<fNstrips_ladder) {
-               real=modff(AdcTrack[iv][ladder*fNstrips_ladder+klp],&inte);
+   fDataTrack[DSPpos+14]=(ReLength >> 8);
-               if (real > 0.5) inte=inte+1;
-               Int_t clvalp=(Int_t)inte -(Int_t)fPedeTrack[iv][ladder*fNstrips_ladder+klp];
+   fDataTrack[DSPpos+15]=( (ReLength << 8) | (OveCheckCode & 0x00FF) );
-               if((iv+1)%2==0) {
+   // end TRAILER
-                 if(clvalp>clval) {
-                    clval=clvalp;
-                    k=klp;}
-                 else break; // max of cluster found
+   // end DSP
-               } else {
+   };
-                 if(clvalp<clval) {
-                    clval=clvalp;
-                    k=klp;}
+ };
-                 else break; // max of cluster found
+ void Digitizer::WriteTrackCalib() {
-               };
-               klp++;
+ std:: cout << " Entering WriteTrackCalib " << endl;
-             };
+ fOutputfile.write(reinterpret_cast<char*>(fDataPSCU),sizeof(UShort_t)*fPSCUbuffer);
-             Int_t kl1=k-fNclst; // max of cluster (or end of ladder ?)
-             trasmesso=0;
+ UShort_t temp[1000000];
-             if(kl1<0)  kl1=0;
+ memset(temp,0,sizeof(UShort_t)*1000000);
+ swab(fDataTrack,temp,sizeof(UShort_t)*fTracklength);  // WE MUST SWAP THE BYTES!!!
-             if(kt>=kl1) kl1=kt+1;
+ fOutputfile.write(reinterpret_cast<char*>(temp),sizeof(UShort_t)*fTracklength);
-             if( (kt+1)==kl1 ) trasmesso=1;
+ fTracklength=0;
+ if ( fPadding ){
+       fOutputfile.write(reinterpret_cast<char*>(fDataPadding),sizeof(UChar_t)*fPadding);
+ };
-             Int_t kl2=k+fNclst;
-             if(kl2>=fNstrips_ladder) kl2=fNstrips_ladder-1;
+ };
-             for(Int_t klt=kl1 ; klt<=kl2 ; klt++) {
-               if(trasmesso==0) {
+ void Digitizer::ClearTrackCalib() {
-               //  std:: cout << "STRIP " << klt << endl;
-               //  std:: cout << "ADC_TOT " <<AdcTrack[iv][ladder*fNstrips_ladder+klt] << endl;
+ std:: cout << "Entering ClearTrackCalib " << endl;
-                 DataDSP[DSPlength]=( ((UShort_t)klt) | 0x1000);
-                 DSPlength++;
+ };
-                 ntrastot++;
+ void Digitizer::LoadTrackCalib() {
-                 real=modff(AdcTrack[iv][ladder*fNstrips_ladder+klt],&inte);
+ std:: cout << "Entering LoadTrackCalib " << endl;
-                 if (real > 0.5) inte=inte+1;
-                 DataDSP[DSPlength]=( ((UShort_t)inte) & 0x0FFF);
+ // Generate the pedestals and sigmas according to parametrization
-                 DSPlength++;
+   for (Int_t j=0; j<fNviews;j++) {
-                 ntrastot++;
+     for (Int_t i=0; i<fNstrips_view;i++) {
-                }
+     if((j+1)%2==0) {
-                else {
+     fPedeTrack[j][i]=gRandom->Gaus(fAvePedex,fSigmaPedex);
-                //   std:: cout << "ADC_TOT " <<AdcTrack[iv][ladder*fNstrips_ladder+klt] << endl;
+     fSigmaTrack[j][i]=gRandom->Gaus(fAveSigmax,fSigmaSigmax);
-                 real=modff(AdcTrack[iv][ladder*fNstrips_ladder+klt],&inte);
+     };
-                 if (real > 0.5) inte=inte+1;
+     if((j+1)%2==1) {
-                 DataDSP[DSPlength]=( ((UShort_t)inte) & 0x0FFF);
+     fPedeTrack[j][i]=gRandom->Gaus(fAvePedey,fSigmaPedey);
-                 DSPlength++;
+     fSigmaTrack[j][i]=gRandom->Gaus(fAveSigmay,fSigmaSigmay);
-                 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++;
+ void Digitizer::LoadMipCor() {
-             continue;
+ std:: cout << "Entering LoadMipCor" << endl;
+   Float_t xfactor=1./151.6*1.04;
+   Float_t yfactor=1./152.1;
-           }; // end cercacluster
-         }; // end cercacluster
+   fMipCor[0][0]=140.02*yfactor;
+   fMipCor[0][1]=140.99*xfactor;
- // start ZOP check for strips no
+   fMipCor[0][2]=134.48*yfactor;
+   fMipCor[0][3]=144.41*xfactor;
-       if(abs(newval-oldval)>=fCutzop*(Int_t)fSigmaTrack[iv][ladder*fNstrips_ladder+k]) {
+   fMipCor[0][4]=140.74*yfactor;
+   fMipCor[0][5]=142.28*xfactor;
-        if(trasmesso==0) {
+   fMipCor[0][6]=134.53*yfactor;
-        // std:: cout << "STRIP " << k << endl;
+   fMipCor[0][7]=140.63*xfactor;
-        // std:: cout << "ADC_TOT " << AdcTrack[iv][ladder*fNstrips_ladder+k] << endl;
+   fMipCor[0][8]=135.55*yfactor;
+   fMipCor[0][9]=138.00*xfactor;
-          DataDSP[DSPlength]=( ((UShort_t)k) | 0x1000);
+   fMipCor[0][10]=154.95*yfactor;
-          DSPlength++;
+   fMipCor[0][11]=158.44*xfactor;
-          ntrastot++;
+   fMipCor[1][0]=136.07*yfactor;
-          real=modff(AdcTrack[iv][ladder*fNstrips_ladder+k],&inte);
+   fMipCor[1][1]=135.59*xfactor;
-          if (real > 0.5) inte=inte+1;
+   fMipCor[1][2]=142.69*yfactor;
-          DataDSP[DSPlength]=( ((UShort_t)inte) & 0x0FFF);
+   fMipCor[1][3]=138.19*xfactor;
-          DSPlength++;
+   fMipCor[1][4]=137.35*yfactor;
-          ntrastot++;
+   fMipCor[1][5]=140.23*xfactor;
+   fMipCor[1][6]=153.15*yfactor;
-        }
+   fMipCor[1][7]=151.42*xfactor;
-        else {
+   fMipCor[1][8]=129.76*yfactor;
-        //  std:: cout << "ADC_TOT " << AdcTrack[iv][ladder*fNstrips_ladder+k] << endl;
+   fMipCor[1][9]=140.63*xfactor;
-          real=modff(AdcTrack[iv][ladder*fNstrips_ladder+k],&inte);
+   fMipCor[1][10]=157.87*yfactor;
-          if (real > 0.5) inte=inte+1;
+   fMipCor[1][11]=153.64*xfactor;
-          DataDSP[DSPlength]=(  ((UShort_t)inte) & 0x0FFF);
-          DSPlength++;
+   fMipCor[2][0]=134.98*yfactor;
-          ntrastot++;
+   fMipCor[2][1]=143.95*xfactor;
-        };
+   fMipCor[2][2]=140.23*yfactor;
-        trasmesso=1;
+   fMipCor[2][3]=138.88*xfactor;
-        oldval=newval;
+   fMipCor[2][4]=137.95*yfactor;
-        kt=k;
+   fMipCor[2][5]=134.87*xfactor;
+   fMipCor[2][6]=157.56*yfactor;
-       }
+   fMipCor[2][7]=157.31*xfactor;
-       else trasmesso=0;
+   fMipCor[2][8]=141.37*yfactor;
-       // end zop
+   fMipCor[2][9]=143.39*xfactor;
+   fMipCor[2][10]=156.15*yfactor;
-       k++;
+   fMipCor[2][11]=158.79*xfactor;
-       };  // end cycle inside ladder
- // write here the end ladder bytes
+ /*
- //            std:: cout << "FINE LADDER " << ladder+1 << endl;
+   for (Int_t j=0; j<fNviews;j++) {
+     for (Int_t i=0; i<fNstrips_view;i++) {
-       DataDSP[DSPlength]=( ((UShort_t)(ladder+1))  | 0x1800);
+     fMipCor[j][i]=1.;
-       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)) );
+ void Digitizer::CompressTrackData(Float_t AdcTrack[fNviews][fNstrips_view]) {
-     UShort_t Nword=(DSPlength*13)/16;
+ // copy of the corresponding compression fortran routine + new digitization
-     if( ((DSPlength*13)%16)!=0) Nword++;
+ // std:: cout << "Entering CompressTrackData " << endl;
-     DataDSP[1]=(0x1400 | ( Nword >> 10));
+ Int_t oldval=0;
-     DataDSP[2]=(0x1400 | ( Nword & 0x03FF) );
+ Int_t newval=0;
-     DataDSP[3]=(0x1400 | (( (UShort_t)(fCounter >> 10) ) & 0x03FF) );
+ Int_t trasmesso=0;
-     DataDSP[4]=(0x1400 | (( (UShort_t)(fCounter) )  & 0x03FF) );
+ Int_t ntrastot=0;
-     DataDSP[5]=(0x1400 | ( (UShort_t)(fNclst << 7) ) | ( (UShort_t)(fCutzop << 4) )
+ Float_t real;
-      | ( (UShort_t)fCutzop  ) );
+ Float_t inte;
-     DataDSP[6]=0x1400;
+ Int_t cercacluster=0;
-     DataDSP[7]=0x1400;
+ Int_t kt=0;
-     DataDSP[8]=0x1400;
+ static const int DSPbufferSize = 4000; // 13 bit buffer to be rearranged in 16 bit Track buffer
-     DataDSP[9]=0x1400;
+ UShort_t DataDSP[DSPbufferSize];  // 13 bit  buffer to be rearranged in 16 bit Track buffer
-     DataDSP[10]=0x1400;
+ UShort_t DSPlength;  // 13 bit buffer to be rearranged in 16 bit Track buffer
-     DataDSP[11]=0x1400;
-     DataDSP[12]=0x1400;
+ memset(fDataTrack,0,sizeof(UShort_t)*fTRACKbuffer); // probably not necessary becouse already done ?
-     DataDSP[13]=0x1400;
+ fTracklength=0;
-     DataDSP[14]=(0x1400 | (CheckSum & 0x00FF) );
-     DataDSP[15]=0x1C00;
+   for (Int_t iv=0; iv<fNviews;iv++) {
- // end DSP header
+   memset(DataDSP,0,sizeof(UShort_t)*DSPbufferSize);
+   DSPlength=16; // skip the header, to be written later
+   UShort_t CheckSum=0;
- // write 13 bit DataDSP bufer inside 16 bit fDataTrack buffer
+ // write dsp header on buffer
-     Int_t Bit16free=16;
-     UShort_t Dato;
+ //    fDataTrack[fTracklength]=0xE805;
-     for (Int_t NDSP=0; NDSP<DSPlength;NDSP++) {
+ //    fTracklength++;
-       Int_t Bit13ToWrite=13;
-       while(Bit13ToWrite>0) {
+ //    fDataTrack[fTracklength]=0x01A9;
-         if(Bit13ToWrite<=Bit16free) {
+ //    fTracklength++;
-           Dato=((DataDSP[NDSP]&(0xFFFF >> (16-Bit13ToWrite)))<<(Bit16free-Bit13ToWrite));
-           fDataTrack[fTracklength]=fDataTrack[fTracklength] | Dato ;
+ // end dsp header
-           Bit16free=Bit16free-Bit13ToWrite;
-           Bit13ToWrite=0;
+    //
-           if(Bit16free==0) {
+    // INIZIO VISTA IV - TAKE PROPER ACTION
-             if(NDSP>15) CheckSum=CheckSum^fDataTrack[fTracklength];
+    //
-             fTracklength++;
-             Bit16free=16;
-           };
-         }
+     for (Int_t ladder=0; ladder<fNladder;ladder++) {
-         else if(Bit13ToWrite>Bit16free) {
+       Int_t k=0;
-           Dato=( (DataDSP[NDSP]&(0xFFFF >> (16-Bit13ToWrite) ) ) >> (Bit13ToWrite-Bit16free) );
+       while (k<fNstrips_ladder) {
-           fDataTrack[fTracklength]=fDataTrack[fTracklength] | Dato ;
+         // compress write in buffer the current LADDER
-           if(NDSP>15) CheckSum=CheckSum^fDataTrack[fTracklength];
+         if ( k == 0)  {
-           fTracklength++;
+           real=modff(AdcTrack[iv][ladder*fNstrips_ladder+k],&inte);
-           Bit13ToWrite=Bit13ToWrite-Bit16free;
+           if (real > 0.5) inte=inte+1;
-           Bit16free=16;
+           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;
-       }; // end cycle while(Bit13ToWrite>0)
+           DataDSP[DSPlength]=( ((UShort_t)inte) & 0x0FFF);
+           DSPlength++;
-     }; // end cycle DataDSP
+           ntrastot++;
-     if(Bit16free!=16) { fTracklength++; CheckSum=CheckSum^fDataTrack[fTracklength]; };
+           trasmesso=1;
-     CheckSum=(CheckSum >> 8)^(CheckSum&0x00FF);
+           oldval=newval;
-     fDataTrack[fTracklength-Nword+11]=(0x0280 | (CheckSum >> 3));
+           kt=k;
-     fDataTrack[fTracklength-Nword+12]=(0x1C00 | (CheckSum << 13) );
+           k++;
+           continue;
- // end write 13 bit DataDSP bufer inside 16 bit fDataTrack buffer
+         };
+         real=modff(AdcTrack[iv][ladder*fNstrips_ladder+k],&inte);
- //write trailer on buffer
+         if (real > 0.5) inte=inte+1;
-     UShort_t ReLength=(UShort_t)((Nword+13)*2+3);
+         newval=(Int_t)inte -(Int_t)(fPedeTrack[iv][ladder*fNstrips_ladder+k]);
-     UShort_t OveCheckCode=0x0000;
+         cercacluster=1; // ?????????
+         if (cercacluster==1) {
-     fDataTrack[fTracklength]=0x0000;
-     fTracklength++;
+  // controlla l'ordine di tutti queste strip ladder e DSP !!!!!!!
+           Int_t diff=0;
-     fDataTrack[fTracklength]=(ReLength >> 8);
-     fTracklength++;
+           switch ((iv+1)%2) {
-     fDataTrack[fTracklength]=( (ReLength << 8) | (OveCheckCode & 0x00FF) );
+           case 0: diff=newval-oldval;
-     fTracklength++;
+           break;
- // end trailer
+           case 1: diff=oldval-newval;
- //    std:: cout  << "DSPlength  " <<DSPlength << endl;
+           break;
- //    std:: cout << "Nword " << Nword  << endl;
+           };
- //    std:: cout <<  "ReLength " << ReLength << endl;
-   };
+           if (diff>fCutclu*(Int_t)fSigmaTrack[iv][ladder*fNstrips_ladder+k]) {
- //    std:: cout << "ntrastot " << ntrastot << endl;
+             Int_t clval=newval;
+             Int_t klp=k;        // go on to search for maximum
- };
+             klp++;
+             while(klp<fNstrips_ladder) {
- Float_t Digitizer::SaturationTrack(Float_t ADC) {
+               real=modff(AdcTrack[iv][ladder*fNstrips_ladder+klp],&inte);
-   Float_t SatFact=1.;
+               if (real > 0.5) inte=inte+1;
-   if(ADC<70.) { SatFact=80./ADC; };
+               Int_t clvalp=(Int_t)inte -(Int_t)fPedeTrack[iv][ladder*fNstrips_ladder+klp];
-   if(ADC>3000.) { SatFact=3000./ADC; };
+               if((iv+1)%2==0) {
-   return SatFact;
- };
+                 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;
+ };

 Legend:



Removed from v.1.2
 


changed lines


 
Added in v.1.4
 Legend:



Removed from v.1.2
 


changed lines


 
Added in v.1.4
-Removed from v.1.2
+Added in v.1.4

	ViewVC Help
Powered by ViewVC 1.1.23