/[PAMELA software]/PamelaDigitizer/Digitizer.cxx
ViewVC logotype

Diff of /PamelaDigitizer/Digitizer.cxx

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.5 by silvio, Wed Nov 28 18:54:31 2007 UTC revision 1.14 by pizzolot, Fri Oct 16 09:15:50 2009 UTC
# Line 1  Line 1 
1  //               ------ PAMELA Digitizer ------  //               ------ PAMELA Digitizer ------
2  //  //
3  // Date, release and how-to: see file Pamelagp2Digits.cxx  // Date, release and how-to: see file  Pamelagp2Digits.cxx
4  //  //
5  // NB: Check length physics packet [packet type (0x10 = physics data)]  // NB: Check length  physics  packet [packet type (0x10 = physics data)]
6  //  //
7  #include <sstream>  #include "Digitizer.h"
8  #include <fstream>  
9  #include <stdlib.h>  extern "C"{
10  #include <stdio.h>    short crc(short, short);
11  #include <string.h>  };
12  #include <ctype.h>  //
13  #include <time.h>  
14  #include "Riostream.h"  Digitizer::Digitizer(TTree* tree, char* &file_raw,
15  #include "TFile.h"                       int nspe1=200,int ntof1=200,int ncat1=50,
16  #include "TDirectory.h"                       int ncas1=50,int ncar1=100,int ncal1=1000,
17  #include "TTree.h"                       int nnd1=200,int nstr1=1000, int comprcalomod=0){
18  #include "TLeafI.h"    
19  #include "TH1.h"    nspe=new int[1];
20  #include "TH2.h"    ntof=new int[1];
21  #include "TMath.h"    ncat=new int[1];
22  #include "TRandom.h"    ncas=new int[1];
23  #include "TSQLServer.h"    ncar=new int[1];
24  #include "TSystem.h"    ncal=new int[1];
25  //    nnd=new int[1];
26  #include "Digitizer.h"    nstr=new int[1];
27  #include "CRC.h"  
28  //    *nspe=nspe1;
29  #include <PamelaRun.h>    *ntof=ntof1;
30  #include <physics/calorimeter/CalorimeterEvent.h>    *ncat=ncat1;
31  #include <CalibCalPedEvent.h>    *ncas=ncas1;
32  #include "GLTables.h"    *ncar=ncar1;
33  //    *ncal=ncal1;
34  extern "C"{    *nnd=nnd1;
35    short crc(short, short);    *nstr=nstr1;
36  };  
37  //    fhBookTree = tree;
38      fFilename =  file_raw;
39  Digitizer::Digitizer(TTree* tree, char* &file_raw){    fCounter = 0;
40    fhBookTree = tree;    fCounterPhys = 0; // SO 5/12/'07
41    fFilename =  file_raw;    fOBT = 0;
42    fCounter = 0;    fModCalo = comprcalomod ;
43    fOBT = 0;  
44      //
45    //    // DB  connections
46    // DB connections    //
47    //    TString host = "mysql://localhost/pamelaprod";
48    TString host = "mysql://localhost/pamelaprod";    TString user = "anonymous";
49    TString user = "anonymous";    TString psw = "";
50    TString psw = "";    //
51    //    const char *pamdbhost=gSystem->Getenv("PAM_DBHOST");
52    const char *pamdbhost=gSystem->Getenv("PAM_DBHOST");    const char *pamdbuser=gSystem->Getenv("PAM_DBUSER");
53    const char *pamdbuser=gSystem->Getenv("PAM_DBUSER");    const char *pamdbpsw=gSystem->Getenv("PAM_DBPSW");
54    const char *pamdbpsw=gSystem->Getenv("PAM_DBPSW");    if ( !pamdbhost ) pamdbhost = "";
55    if ( !pamdbhost ) pamdbhost = "";    if ( !pamdbuser ) pamdbuser = "";
56    if ( !pamdbuser ) pamdbuser = "";    if ( !pamdbpsw ) pamdbpsw = "";
57    if ( !pamdbpsw ) pamdbpsw = "";    if ( strcmp(pamdbhost,"") ) host = pamdbhost;
58    if ( strcmp(pamdbhost,"") ) host = pamdbhost;    if ( strcmp(pamdbuser,"") ) user = pamdbuser;
59    if ( strcmp(pamdbuser,"") ) user = pamdbuser;    if ( strcmp(pamdbpsw,"") ) psw = pamdbpsw;
60    if ( strcmp(pamdbpsw,"") ) psw = pamdbpsw;    fDbc = TSQLServer::Connect(host.Data(),user.Data(),psw.Data());
61    fDbc = TSQLServer::Connect(host.Data(),user.Data(),psw.Data());    //
62    //    GL_TABLES *glt = new GL_TABLES(host,user,psw);
63    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());
64    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());    //
65    //    // Use UTC in the DB and make timeout bigger
66    // Use UTC in the DB and make timeout bigger    //
67    //    stringstream myquery;
68    stringstream myquery;    myquery.str("");
69    myquery.str("");    myquery << "SET time_zone='+0:00'";
70    myquery << "SET time_zone='+0:00'";    fDbc->Query(myquery.str().c_str());
71    fDbc->Query(myquery.str().c_str());    myquery.str("");
72    myquery.str("");    myquery << "SET wait_timeout=173000;";
73    myquery << "SET wait_timeout=173000;";    fDbc->Query(myquery.str().c_str());
74    fDbc->Query(myquery.str().c_str());    //
75    //    
76        std:: cout << "preparing tree" << endl;
77    std:: cout << "preparing tree" << endl;  
78      Ipltof=(UChar_t*)malloc(*ntof *sizeof(UChar_t));
79    // prepare tree    Ipaddle=(UChar_t*)malloc(*ntof *sizeof(UChar_t));
80    fhBookTree->SetBranchAddress("Irun",&Irun);    Ipartof=(UShort_t*)malloc(*ntof *sizeof(UShort_t));
81    fhBookTree->SetBranchAddress("Ievnt",&Ievnt);    //  Ipartof=(UChar_t*)malloc(*ntof *sizeof(UChar_t));//DPMJET
82    fhBookTree->SetBranchAddress("Ipa",&Ipa);    Xintof=(Float_t*)malloc(*ntof *sizeof(Float_t));
83    fhBookTree->SetBranchAddress("X0",&X0);    Yintof=(Float_t*)malloc(*ntof *sizeof(Float_t));
84    fhBookTree->SetBranchAddress("Y0",&Y0);    Zintof=(Float_t*)malloc(*ntof *sizeof(Float_t));
85    fhBookTree->SetBranchAddress("Z0",&Z0);    Xouttof=(Float_t*)malloc(*ntof *sizeof(Float_t));
86    fhBookTree->SetBranchAddress("Theta",&Theta);    Youttof=(Float_t*)malloc(*ntof *sizeof(Float_t));
87    fhBookTree->SetBranchAddress("Phi",&Phi);    Zouttof=(Float_t*)malloc(*ntof *sizeof(Float_t));
88    fhBookTree->SetBranchAddress("P0",&P0);    Ereltof=(Float_t*)malloc(*ntof *sizeof(Float_t));
89    fhBookTree->SetBranchAddress("Nthtof",&Nthtof);    Timetof=(Float_t*)malloc(*ntof *sizeof(Float_t));
90    fhBookTree->SetBranchAddress("Ipltof",Ipltof);    Pathtof=(Float_t*)malloc(*ntof *sizeof(Float_t));
91    fhBookTree->SetBranchAddress("Ipaddle",Ipaddle);    P0tof=(Float_t*)malloc(*ntof *sizeof(Float_t));
92    fhBookTree->SetBranchAddress("Ipartof",Ipartof);    Iparcat=(UChar_t*)malloc(*ncat *sizeof(UChar_t));
93    fhBookTree->SetBranchAddress("Xintof",Xintof);    Icat=(UChar_t*)malloc(*ncat *sizeof(UChar_t));
94    fhBookTree->SetBranchAddress("Yintof",Yintof);    Xincat=(Float_t*)malloc(*ncat *sizeof(Float_t));
95    fhBookTree->SetBranchAddress("Zintof",Zintof);    Yincat=(Float_t*)malloc(*ncat *sizeof(Float_t));
96    fhBookTree->SetBranchAddress("Xouttof",Xouttof);    Zincat=(Float_t*)malloc(*ncat *sizeof(Float_t));
97    fhBookTree->SetBranchAddress("Youttof",Youttof);    Xoutcat=(Float_t*)malloc(*ncat *sizeof(Float_t));
98    fhBookTree->SetBranchAddress("Zouttof",Zouttof);    Youtcat=(Float_t*)malloc(*ncat *sizeof(Float_t));
99    fhBookTree->SetBranchAddress("Ereltof",Ereltof);    Zoutcat=(Float_t*)malloc(*ncat *sizeof(Float_t));
100    fhBookTree->SetBranchAddress("Timetof",Timetof);    Erelcat=(Float_t*)malloc(*ncat *sizeof(Float_t));
101    fhBookTree->SetBranchAddress("Pathtof",Pathtof);    Timecat=(Float_t*)malloc(*ncat *sizeof(Float_t));
102    fhBookTree->SetBranchAddress("P0tof",P0tof);    Pathcat=(Float_t*)malloc(*ncat *sizeof(Float_t));
103    fhBookTree->SetBranchAddress("Nthcat",&Nthcat);    P0cat=(Float_t*)malloc(*ncat *sizeof(Float_t));
104    fhBookTree->SetBranchAddress("Iparcat",Iparcat);    Iparcas=(UChar_t*)malloc(*ncas *sizeof(UChar_t));
105    fhBookTree->SetBranchAddress("Icat",Icat);    Icas=(UChar_t*)malloc(*ncas *sizeof(UChar_t));
106    fhBookTree->SetBranchAddress("Xincat",Xincat);    Xincas=(Float_t*)malloc(*ncas *sizeof(Float_t));
107    fhBookTree->SetBranchAddress("Yincat",Yincat);    Yincas=(Float_t*)malloc(*ncas *sizeof(Float_t));
108    fhBookTree->SetBranchAddress("Zincat",Zincat);    Zincas=(Float_t*)malloc(*ncas *sizeof(Float_t));
109    fhBookTree->SetBranchAddress("Xoutcat",Xoutcat);    Xoutcas=(Float_t*)malloc(*ncas *sizeof(Float_t));
110    fhBookTree->SetBranchAddress("Youtcat",Youtcat);    Youtcas=(Float_t*)malloc(*ncas *sizeof(Float_t));
111    fhBookTree->SetBranchAddress("Zoutcat",Zoutcat);    Zoutcas=(Float_t*)malloc(*ncas *sizeof(Float_t));
112    fhBookTree->SetBranchAddress("Erelcat",Erelcat);    Erelcas=(Float_t*)malloc(*ncas *sizeof(Float_t));
113    fhBookTree->SetBranchAddress("Timecat",Timecat);    Timecas=(Float_t*)malloc(*ncas *sizeof(Float_t));
114    fhBookTree->SetBranchAddress("Pathcat",Pathcat);    Pathcas=(Float_t*)malloc(*ncas *sizeof(Float_t));
115    fhBookTree->SetBranchAddress("P0cat",P0cat);    P0cas=(Float_t*)malloc(*ncas *sizeof(Float_t));
116    fhBookTree->SetBranchAddress("Nthcas",&Nthcas);    //  Iparspe=(UShort_t*)malloc(*nspe *sizeof(UShort_t));
117    fhBookTree->SetBranchAddress("Iparcas",Iparcas);    //  Iparspe=(UChar_t*)malloc(*nspe *sizeof(UChar_t));
118    fhBookTree->SetBranchAddress("Icas",Icas);    Itrpb=(UChar_t*)malloc(*nspe *sizeof(UChar_t));
119    fhBookTree->SetBranchAddress("Xincas",Xincas);    Itrsl=(UChar_t*)malloc(*nspe *sizeof(UChar_t));
120    fhBookTree->SetBranchAddress("Yincas",Yincas);    Itspa=(UChar_t*)malloc(*nspe *sizeof(UChar_t));
121    fhBookTree->SetBranchAddress("Zincas",Zincas);    Xinspe=(Float_t*)malloc(*nspe *sizeof(Float_t));
122    fhBookTree->SetBranchAddress("Xoutcas",Xoutcas);    Yinspe=(Float_t*)malloc(*nspe *sizeof(Float_t));
123    fhBookTree->SetBranchAddress("Youtcas",Youtcas);    Zinspe=(Float_t*)malloc(*nspe *sizeof(Float_t));
124    fhBookTree->SetBranchAddress("Zoutcas",Zoutcas);    Xoutspe=(Float_t*)malloc(*nspe *sizeof(Float_t));
125    fhBookTree->SetBranchAddress("Erelcas",Erelcas);    Youtspe=(Float_t*)malloc(*nspe *sizeof(Float_t));
126    fhBookTree->SetBranchAddress("Timecas",Timecas);    Zoutspe=(Float_t*)malloc(*nspe *sizeof(Float_t));
127    fhBookTree->SetBranchAddress("Pathcas",Pathcas);    Xavspe=(Float_t*)malloc(*nspe *sizeof(Float_t));
128    fhBookTree->SetBranchAddress("P0cas",P0cas);    Yavspe=(Float_t*)malloc(*nspe *sizeof(Float_t));
129    fhBookTree->SetBranchAddress("Nthspe",&Nthspe);    Zavspe=(Float_t*)malloc(*nspe *sizeof(Float_t));
130    fhBookTree->SetBranchAddress("Iparspe",Iparspe);    Erelspe=(Float_t*)malloc(*nspe *sizeof(Float_t));
131    fhBookTree->SetBranchAddress("Itrpb",Itrpb);    Pathspe=(Float_t*)malloc(*nspe *sizeof(Float_t));
132    fhBookTree->SetBranchAddress("Itrsl",Itrsl);    P0spe=(Float_t*)malloc(*nspe *sizeof(Float_t));;
133    fhBookTree->SetBranchAddress("Itspa",Itspa);    Nxmult=(UChar_t*)malloc(*nspe *sizeof(UChar_t));
134    fhBookTree->SetBranchAddress("Xinspe",Xinspe);    Nymult=(UChar_t*)malloc(*nspe *sizeof(UChar_t));
135    fhBookTree->SetBranchAddress("Yinspe",Yinspe);    Istripx=(UShort_t*)malloc(*nstr *sizeof(UShort_t));
136    fhBookTree->SetBranchAddress("Zinspe",Zinspe);    Qstripx=(Float_t*)malloc(*nstr *sizeof(Float_t));
137    fhBookTree->SetBranchAddress("Xoutspe",Xoutspe);    Xstripx=(Float_t*)malloc(*nstr *sizeof(Float_t));
138    fhBookTree->SetBranchAddress("Youtspe",Youtspe);    Npstripx=(UChar_t*)malloc(*nstr *sizeof(UChar_t));
139    fhBookTree->SetBranchAddress("Zoutspe",Zoutspe);    Ntstripx=(UChar_t*)malloc(*nstr *sizeof(UChar_t));
140    fhBookTree->SetBranchAddress("Xavspe",Xavspe);    Npstripy=(UChar_t*)malloc(*nstr *sizeof(UChar_t));
141    fhBookTree->SetBranchAddress("Yavspe",Yavspe);    Ntstripy=(UChar_t*)malloc(*nstr *sizeof(UChar_t));
142    fhBookTree->SetBranchAddress("Zavspe",Zavspe);    Istripy=(UShort_t*)malloc(*nstr *sizeof(UShort_t));
143    fhBookTree->SetBranchAddress("Erelspe",Erelspe);    Qstripy=(Float_t*)malloc(*nstr *sizeof(Float_t));
144    fhBookTree->SetBranchAddress("Pathspe",Pathspe);    Ystripy=(Float_t*)malloc(*nstr *sizeof(Float_t));
145    fhBookTree->SetBranchAddress("P0spe",P0spe);    Icapl=(UChar_t*)malloc(*ncal *sizeof(UChar_t));
146    fhBookTree->SetBranchAddress("Nxmult",Nxmult);    Icasi=(UChar_t*)malloc(*ncal *sizeof(UChar_t));
147    fhBookTree->SetBranchAddress("Nymult",Nymult);    Icast=(UChar_t*)malloc(*ncal *sizeof(UChar_t));
148    fhBookTree->SetBranchAddress("Nstrpx",&Nstrpx);    Xincal=(Float_t*)malloc(*ncal *sizeof(Float_t));
149    fhBookTree->SetBranchAddress("Npstripx",Npstripx);    Yincal=(Float_t*)malloc(*ncal *sizeof(Float_t));
150    fhBookTree->SetBranchAddress("Ntstripx",Ntstripx);    Zincal=(Float_t*)malloc(*ncal *sizeof(Float_t));
151    fhBookTree->SetBranchAddress("Istripx",Istripx);    Erelcal=(Float_t*)malloc(*ncal *sizeof(Float_t));
152    fhBookTree->SetBranchAddress("Qstripx",Qstripx);    Itubend=(UChar_t*)malloc(*nnd *sizeof(UChar_t));
153    fhBookTree->SetBranchAddress("Xstripx",Xstripx);    Iparnd=(UChar_t*)malloc(*nnd *sizeof(UChar_t));
154    fhBookTree->SetBranchAddress("Nstrpy",&Nstrpy);    Xinnd=(Float_t*)malloc(*nnd *sizeof(Float_t));
155    fhBookTree->SetBranchAddress("Npstripy",Npstripy);    Yinnd=(Float_t*)malloc(*nnd *sizeof(Float_t));
156    fhBookTree->SetBranchAddress("Ntstripy",Ntstripy);    Zinnd=(Float_t*)malloc(*nnd *sizeof(Float_t));
157    fhBookTree->SetBranchAddress("Istripy",Istripy);    Xoutnd=(Float_t*)malloc(*nnd *sizeof(Float_t));
158    fhBookTree->SetBranchAddress("Qstripy",Qstripy);    Youtnd=(Float_t*)malloc(*nnd *sizeof(Float_t));
159    fhBookTree->SetBranchAddress("Ystripy",Ystripy);    Zoutnd=(Float_t*)malloc(*nnd *sizeof(Float_t));
160    fhBookTree->SetBranchAddress("Nthcali",&Nthcali);    Erelnd=(Float_t*)malloc(*nnd *sizeof(Float_t));
161    fhBookTree->SetBranchAddress("Icaplane",Icaplane);    Timend=(Float_t*)malloc(*nnd *sizeof(Float_t));
162    fhBookTree->SetBranchAddress("Icastrip",Icastrip);    Pathnd=(Float_t*)malloc(*nnd *sizeof(Float_t));
163    fhBookTree->SetBranchAddress("Icamod",Icamod);    P0nd=(Float_t*)malloc(*nnd *sizeof(Float_t));
164    fhBookTree->SetBranchAddress("Enestrip",Enestrip);    Iparcard=(UChar_t*)malloc(*ncar *sizeof(UChar_t));
165    fhBookTree->SetBranchAddress("Nthcal",&Nthcal);    Icard=(UChar_t*)malloc(*ncar *sizeof(UChar_t));
166    fhBookTree->SetBranchAddress("Icapl",Icapl);    Xincard=(Float_t*)malloc(*ncar *sizeof(Float_t));
167    fhBookTree->SetBranchAddress("Icasi",Icasi);    Yincard=(Float_t*)malloc(*ncar *sizeof(Float_t));
168    fhBookTree->SetBranchAddress("Icast",Icast);    Zincard=(Float_t*)malloc(*ncar *sizeof(Float_t));
169    fhBookTree->SetBranchAddress("Xincal",Xincal);    Xoutcard=(Float_t*)malloc(*ncar *sizeof(Float_t));
170    fhBookTree->SetBranchAddress("Yincal",Yincal);    Youtcard=(Float_t*)malloc(*ncar *sizeof(Float_t));
171    fhBookTree->SetBranchAddress("Zincal",Zincal);    Zoutcard=(Float_t*)malloc(*ncar *sizeof(Float_t));
172    fhBookTree->SetBranchAddress("Erelcal",Erelcal);    Erelcard=(Float_t*)malloc(*ncar *sizeof(Float_t));
173    fhBookTree->SetBranchAddress("Nthnd",&Nthnd);    Timecard=(Float_t*)malloc(*ncar *sizeof(Float_t));
174    fhBookTree->SetBranchAddress("Itubend",Itubend);    Pathcard=(Float_t*)malloc(*ncar *sizeof(Float_t));
175    fhBookTree->SetBranchAddress("Iparnd",Iparnd);    P0card=(Float_t*)malloc(*ncar *sizeof(Float_t));
176    fhBookTree->SetBranchAddress("Xinnd",Xinnd);  
177    fhBookTree->SetBranchAddress("Yinnd",Yinnd);  
178    fhBookTree->SetBranchAddress("Zinnd",Zinnd);  
179    fhBookTree->SetBranchAddress("Xoutnd",Xoutnd);    // prepare tree//modified by E.Vannuccini 03/08
180    fhBookTree->SetBranchAddress("Youtnd",Youtnd);    if(fhBookTree->GetBranch("Irun"))fhBookTree->SetBranchAddress("Irun",&Irun);
181    fhBookTree->SetBranchAddress("Zoutnd",Zoutnd);    if(fhBookTree->GetBranch("Ievnt"))fhBookTree->SetBranchAddress("Ievnt",&Ievnt);
182    fhBookTree->SetBranchAddress("Erelnd",Erelnd);    if(fhBookTree->GetBranch("Ipa"))fhBookTree->SetBranchAddress("Ipa",&Ipa);
183    fhBookTree->SetBranchAddress("Timend",Timend);    if(fhBookTree->GetBranch("X0"))fhBookTree->SetBranchAddress("X0",&X0);
184    fhBookTree->SetBranchAddress("Pathnd",Pathnd);    if(fhBookTree->GetBranch("Y0"))fhBookTree->SetBranchAddress("Y0",&Y0);
185    fhBookTree->SetBranchAddress("P0nd",P0nd);    if(fhBookTree->GetBranch("Z0"))fhBookTree->SetBranchAddress("Z0",&Z0);
186    fhBookTree->SetBranchAddress("Nthcard",&Nthcard);    if(fhBookTree->GetBranch("Theta"))fhBookTree->SetBranchAddress("Theta",&Theta);
187    fhBookTree->SetBranchAddress("Iparcard",Iparcard);    if(fhBookTree->GetBranch("Phi"))fhBookTree->SetBranchAddress("Phi",&Phi);
188    fhBookTree->SetBranchAddress("Icard",Icard);    if(fhBookTree->GetBranch("P0"))fhBookTree->SetBranchAddress("P0",&P0);
189    fhBookTree->SetBranchAddress("Xincard",Xincard);    if(fhBookTree->GetBranch("Nthtof"))fhBookTree->SetBranchAddress("Nthtof",&Nthtof);
190    fhBookTree->SetBranchAddress("Yincard",Yincard);    if(fhBookTree->GetBranch("Ipltof"))fhBookTree->SetBranchAddress("Ipltof",Ipltof);///////////////////////////
191    fhBookTree->SetBranchAddress("Zincard",Zincard);    if(fhBookTree->GetBranch("Ipaddle"))fhBookTree->SetBranchAddress("Ipaddle",Ipaddle);
192    fhBookTree->SetBranchAddress("Xoutcard",Xoutcard);    if(fhBookTree->GetBranch("Ipartof"))fhBookTree->SetBranchAddress("Ipartof",Ipartof);
193    fhBookTree->SetBranchAddress("Youtcard",Youtcard);    if(fhBookTree->GetBranch("Xintof"))fhBookTree->SetBranchAddress("Xintof",Xintof);
194    fhBookTree->SetBranchAddress("Zoutcard",Zoutcard);    if(fhBookTree->GetBranch("Yintof"))fhBookTree->SetBranchAddress("Yintof",Yintof);
195    fhBookTree->SetBranchAddress("Erelcard",Erelcard);    if(fhBookTree->GetBranch("Zintof"))fhBookTree->SetBranchAddress("Zintof",Zintof);
196    fhBookTree->SetBranchAddress("Timecard",Timecard);    if(fhBookTree->GetBranch("Xouttof"))fhBookTree->SetBranchAddress("Xouttof",Xouttof);
197    fhBookTree->SetBranchAddress("Pathcard",Pathcard);    if(fhBookTree->GetBranch("Youttof"))fhBookTree->SetBranchAddress("Youttof",Youttof);
198    fhBookTree->SetBranchAddress("P0card",P0card);    if(fhBookTree->GetBranch("Zouttof"))fhBookTree->SetBranchAddress("Zouttof",Zouttof);
199      if(fhBookTree->GetBranch("Ereltof"))fhBookTree->SetBranchAddress("Ereltof",Ereltof);
200    fhBookTree->SetBranchStatus("*",0);    if(fhBookTree->GetBranch("Timetof"))fhBookTree->SetBranchAddress("Timetof",Timetof);
201      if(fhBookTree->GetBranch("Pathtof"))fhBookTree->SetBranchAddress("Pathtof",Pathtof);
202  };    if(fhBookTree->GetBranch("P0tof"))fhBookTree->SetBranchAddress("P0tof",P0tof);
203      if(fhBookTree->GetBranch("Nthcat"))fhBookTree->SetBranchAddress("Nthcat",&Nthcat);
204      if(fhBookTree->GetBranch("Iparcat"))fhBookTree->SetBranchAddress("Iparcat",Iparcat);
205      if(fhBookTree->GetBranch("Icat"))fhBookTree->SetBranchAddress("Icat",Icat);
206  void Digitizer::Close(){    if(fhBookTree->GetBranch("Xincat"))fhBookTree->SetBranchAddress("Xincat",Xincat);
207      if(fhBookTree->GetBranch("Yincat"))fhBookTree->SetBranchAddress("Yincat",Yincat);
208    delete fhBookTree;    if(fhBookTree->GetBranch("Zincat"))fhBookTree->SetBranchAddress("Zincat",Zincat);
209      if(fhBookTree->GetBranch("Xoutcat"))fhBookTree->SetBranchAddress("Xoutcat",Xoutcat);
210  };    if(fhBookTree->GetBranch("Youtcat"))fhBookTree->SetBranchAddress("Youtcat",Youtcat);
211      if(fhBookTree->GetBranch("Zoutcat"))fhBookTree->SetBranchAddress("Zoutcat",Zoutcat);
212      if(fhBookTree->GetBranch("Erelcat"))fhBookTree->SetBranchAddress("Erelcat",Erelcat);
213      if(fhBookTree->GetBranch("Timecat"))fhBookTree->SetBranchAddress("Timecat",Timecat);
214      if(fhBookTree->GetBranch("Pathcat"))fhBookTree->SetBranchAddress("Pathcat",Pathcat);
215  void Digitizer::Loop() {    if(fhBookTree->GetBranch("P0cat"))fhBookTree->SetBranchAddress("P0cat",P0cat);
216    //    if(fhBookTree->GetBranch("Nthcas"))fhBookTree->SetBranchAddress("Nthcas",&Nthcas);
217    // opens the raw output file and loops over the events    if(fhBookTree->GetBranch("Iparcas"))fhBookTree->SetBranchAddress("Iparcas",Iparcas);
218    //    if(fhBookTree->GetBranch("Icas"))fhBookTree->SetBranchAddress("Icas",Icas);///////////////////////////////
219    fOutputfile.open(fFilename, ios::out | ios::binary);    if(fhBookTree->GetBranch("Xincas"))fhBookTree->SetBranchAddress("Xincas",Xincas);
220    //fOutputfile.open(Form("Output%s",fFilename), ios::out | ios::binary);    if(fhBookTree->GetBranch("Yincas"))fhBookTree->SetBranchAddress("Yincas",Yincas);
221    //    if(fhBookTree->GetBranch("Zincas"))fhBookTree->SetBranchAddress("Zincas",Zincas);
222    // Load in memory and save at the beginning of file the calorimeter calibration    if(fhBookTree->GetBranch("Xoutcas"))fhBookTree->SetBranchAddress("Xoutcas",Xoutcas);
223    //    if(fhBookTree->GetBranch("Youtcas"))fhBookTree->SetBranchAddress("Youtcas",Youtcas);
224    CaloLoadCalib();    if(fhBookTree->GetBranch("Zoutcas"))fhBookTree->SetBranchAddress("Zoutcas",Zoutcas);
225    DigitizeCALOCALIB();    if(fhBookTree->GetBranch("Erelcas"))fhBookTree->SetBranchAddress("Erelcas",Erelcas);
226      if(fhBookTree->GetBranch("Timecas"))fhBookTree->SetBranchAddress("Timecas",Timecas);
227    //  load, digitize and write tracker calibration    if(fhBookTree->GetBranch("Pathcas"))fhBookTree->SetBranchAddress("Pathcas",Pathcas);
228    LoadTrackCalib();    if(fhBookTree->GetBranch("P0cas"))fhBookTree->SetBranchAddress("P0cas",P0cas);
229        if(fhBookTree->GetBranch("Nthspe"))fhBookTree->SetBranchAddress("Nthspe",&Nthspe);
230    DigitizeTrackCalib(1);    //  if(fhBookTree->GetBranch("Iparspe"))fhBookTree->SetBranchAddress("Iparspe",Iparspe);
231    UInt_t length=fTracklength*2;    if(fhBookTree->GetBranch("Itrpb"))fhBookTree->SetBranchAddress("Itrpb",Itrpb);
232    DigitizePSCU(length,0x12);    if(fhBookTree->GetBranch("Itrsl"))fhBookTree->SetBranchAddress("Itrsl",Itrsl);
233    AddPadding();    if(fhBookTree->GetBranch("Itspa"))fhBookTree->SetBranchAddress("Itspa",Itspa);
234    WriteTrackCalib();    if(fhBookTree->GetBranch("Xinspe"))fhBookTree->SetBranchAddress("Xinspe",Xinspe);
235        if(fhBookTree->GetBranch("Yinspe"))fhBookTree->SetBranchAddress("Yinspe",Yinspe);
236    DigitizeTrackCalib(2);    if(fhBookTree->GetBranch("Zinspe"))fhBookTree->SetBranchAddress("Zinspe",Zinspe);
237    length=fTracklength*2;    if(fhBookTree->GetBranch("Xoutspe"))fhBookTree->SetBranchAddress("Xoutspe",Xoutspe);
238    DigitizePSCU(length,0x13);    if(fhBookTree->GetBranch("Youtspe"))fhBookTree->SetBranchAddress("Youtspe",Youtspe);
239    AddPadding();    if(fhBookTree->GetBranch("Zoutspe"))fhBookTree->SetBranchAddress("Zoutspe",Zoutspe);
240    WriteTrackCalib();    if(fhBookTree->GetBranch("Xavspe"))fhBookTree->SetBranchAddress("Xavspe",Xavspe);
241        if(fhBookTree->GetBranch("Yavspe"))fhBookTree->SetBranchAddress("Yavspe",Yavspe);
242    LoadMipCor();  // some initialization of parameters -not used now-    if(fhBookTree->GetBranch("Zavspe"))fhBookTree->SetBranchAddress("Zavspe",Zavspe);
243    //  end loading, digitizing and writing tracker calibration    if(fhBookTree->GetBranch("Erelspe"))fhBookTree->SetBranchAddress("Erelspe",Erelspe);
244      if(fhBookTree->GetBranch("Pathspe"))fhBookTree->SetBranchAddress("Pathspe",Pathspe);
245    //    if(fhBookTree->GetBranch("P0spe"))fhBookTree->SetBranchAddress("P0spe",P0spe);
246    // loops over the events    if(fhBookTree->GetBranch("Nxmult"))fhBookTree->SetBranchAddress("Nxmult",Nxmult);
247    //    if(fhBookTree->GetBranch("Nymult"))fhBookTree->SetBranchAddress("Nymult",Nymult);
248        if(fhBookTree->GetBranch("Nstrpx"))fhBookTree->SetBranchAddress("Nstrpx",&Nstrpx);
249    Int_t nentries = fhBookTree->GetEntriesFast();    if(fhBookTree->GetBranch("Npstripx"))fhBookTree->SetBranchAddress("Npstripx",Npstripx);
250    Long64_t nbytes = 0;    if(fhBookTree->GetBranch("Ntstripx"))fhBookTree->SetBranchAddress("Ntstripx",Ntstripx);
251    for (Int_t i=0; i<nentries;i++) {    if(fhBookTree->GetBranch("Istripx"))fhBookTree->SetBranchAddress("Istripx",Istripx);
252      //    if(fhBookTree->GetBranch("Qstripx"))fhBookTree->SetBranchAddress("Qstripx",Qstripx);
253      nbytes += fhBookTree->GetEntry(i);    if(fhBookTree->GetBranch("Xstripx"))fhBookTree->SetBranchAddress("Xstripx",Xstripx);
254      // read detectors sequentially:    if(fhBookTree->GetBranch("Nstrpy"))fhBookTree->SetBranchAddress("Nstrpy",&Nstrpy);
255      // http://www.ts.infn.it/fileadmin/documents/physics/experiments/wizard/cpu/gen_arch/RM_Acquisition.pdf    if(fhBookTree->GetBranch("Npstripy"))fhBookTree->SetBranchAddress("Npstripy",Npstripy);
256      // on pamelatov:    if(fhBookTree->GetBranch("Ntstripy"))fhBookTree->SetBranchAddress("Ntstripy",Ntstripy);
257      // /cvs/yoda/techmodel/physics/NeutronDetectorReader.cpp    if(fhBookTree->GetBranch("Istripy"))fhBookTree->SetBranchAddress("Istripy",Istripy);
258      DigitizeTRIGGER();    if(fhBookTree->GetBranch("Qstripy"))fhBookTree->SetBranchAddress("Qstripy",Qstripy);///////////////////////
259      DigitizeTOF();    if(fhBookTree->GetBranch("Ystripy"))fhBookTree->SetBranchAddress("Ystripy",Ystripy);
260      DigitizeAC();    if(fhBookTree->GetBranch("Nthcali"))fhBookTree->SetBranchAddress("Nthcali",&Nthcali);
261      DigitizeCALO();    if(fhBookTree->GetBranch("Icaplane"))fhBookTree->SetBranchAddress("Icaplane",Icaplane);
262      DigitizeTrack();    if(fhBookTree->GetBranch("Icastrip"))fhBookTree->SetBranchAddress("Icastrip",Icastrip);
263      DigitizeS4();    if(fhBookTree->GetBranch("Icamod"))fhBookTree->SetBranchAddress("Icamod",Icamod);
264      DigitizeND();    if(fhBookTree->GetBranch("Enestrip"))fhBookTree->SetBranchAddress("Enestrip",Enestrip);
265        //    if(fhBookTree->GetBranch("Nthcal"))fhBookTree->SetBranchAddress("Nthcal",&Nthcal);
266      // Add padding to 64 bits    if(fhBookTree->GetBranch("Icapl"))fhBookTree->SetBranchAddress("Icapl",Icapl);
267      //    if(fhBookTree->GetBranch("Icasi"))fhBookTree->SetBranchAddress("Icasi",Icasi);
268      AddPadding();    if(fhBookTree->GetBranch("Icast"))fhBookTree->SetBranchAddress("Icast",Icast);
269  //    if(fhBookTree->GetBranch("Xincal"))fhBookTree->SetBranchAddress("Xincal",Xincal);
270      // Create CPU header, we need packet type (0x10 = physics data) and packet length.    if(fhBookTree->GetBranch("Yincal"))fhBookTree->SetBranchAddress("Yincal",Yincal);
271      //    if(fhBookTree->GetBranch("Zincal"))fhBookTree->SetBranchAddress("Zincal",Zincal);
272      UInt_t length=2*(fCALOlength+fACbuffer+fTracklength+fNDbuffer+fS4buffer)+fPadding+fTOFbuffer+fTRIGGERbuffer;    if(fhBookTree->GetBranch("Erelcal"))fhBookTree->SetBranchAddress("Erelcal",Erelcal);
273      //UInt_t length=2*(fCALOlength+fACbuffer+fTracklength+fNDbuffer)+fPadding+fTOFbuffer+fTRIGGERbuffer;    if(fhBookTree->GetBranch("Nthnd"))fhBookTree->SetBranchAddress("Nthnd",&Nthnd);
274      DigitizePSCU(length,0x10);    if(fhBookTree->GetBranch("Itubend"))fhBookTree->SetBranchAddress("Itubend",Itubend);
275      if ( !i%100 ) std::cout << "writing event " << i << endl;    if(fhBookTree->GetBranch("Iparnd"))fhBookTree->SetBranchAddress("Iparnd",Iparnd);
276      WriteData();    if(fhBookTree->GetBranch("Xinnd"))fhBookTree->SetBranchAddress("Xinnd",Xinnd);/////////////////////////
277    };    if(fhBookTree->GetBranch("Yinnd"))fhBookTree->SetBranchAddress("Yinnd",Yinnd);
278      if(fhBookTree->GetBranch("Zinnd"))fhBookTree->SetBranchAddress("Zinnd",Zinnd);
279    fOutputfile.close();    if(fhBookTree->GetBranch("Xoutnd"))fhBookTree->SetBranchAddress("Xoutnd",Xoutnd);
280    std::cout << "files closed" << endl << flush;    if(fhBookTree->GetBranch("Youtnd"))fhBookTree->SetBranchAddress("Youtnd",Youtnd);
281      if(fhBookTree->GetBranch("Zoutnd"))fhBookTree->SetBranchAddress("Zoutnd",Zoutnd);
282  };    if(fhBookTree->GetBranch("Erelnd"))fhBookTree->SetBranchAddress("Erelnd",Erelnd);
283      if(fhBookTree->GetBranch("Timend"))fhBookTree->SetBranchAddress("Timend",Timend);
284  void Digitizer::AddPadding(){    if(fhBookTree->GetBranch("Pathnd"))fhBookTree->SetBranchAddress("Pathnd",Pathnd);
285    //    if(fhBookTree->GetBranch("P0nd"))fhBookTree->SetBranchAddress("P0nd",P0nd);
286    Float_t pd0 = (fLen+16)/64.;    if(fhBookTree->GetBranch("Nthcard"))fhBookTree->SetBranchAddress("Nthcard",&Nthcard);/////////////////////
287    Float_t pd1 = pd0 - (Float_t)int(pd0);    if(fhBookTree->GetBranch("Iparcard"))fhBookTree->SetBranchAddress("Iparcard",Iparcard);
288    Float_t padfrac = 64. - pd1 * 64.;    if(fhBookTree->GetBranch("Icard"))fhBookTree->SetBranchAddress("Icard",Icard);
289    //    if(fhBookTree->GetBranch("Xincard"))fhBookTree->SetBranchAddress("Xincard",Xincard);
290    UInt_t padbytes = (UInt_t)padfrac;    if(fhBookTree->GetBranch("Yincard"))fhBookTree->SetBranchAddress("Yincard",Yincard);
291    if ( padbytes > 0 && padbytes < 64 ){    if(fhBookTree->GetBranch("Zincard"))fhBookTree->SetBranchAddress("Zincard",Zincard);
292      //    if(fhBookTree->GetBranch("Xoutcard"))fhBookTree->SetBranchAddress("Xoutcard",Xoutcard);
293      // here the padding length    if(fhBookTree->GetBranch("Youtcard"))fhBookTree->SetBranchAddress("Youtcard",Youtcard);/////////////////
294      //    if(fhBookTree->GetBranch("Zoutcard"))fhBookTree->SetBranchAddress("Zoutcard",Zoutcard);
295      fPadding = padbytes+64;    if(fhBookTree->GetBranch("Erelcard"))fhBookTree->SetBranchAddress("Erelcard",Erelcard);
296      //    if(fhBookTree->GetBranch("Timecard"))fhBookTree->SetBranchAddress("Timecard",Timecard);
297      // random padding bytes    if(fhBookTree->GetBranch("Pathcard"))fhBookTree->SetBranchAddress("Pathcard",Pathcard);
298      //    //  if(fhBookTree->GetBranch("P0card"))fhBookTree->SetBranchAddress("P0card",P0card);
299      for (Int_t ur=0; ur<32; ur++){  //    fhBookTree->SetBranchStatus("*",0); //modified by E.Vannuccini 03/08
300        fDataPadding[ur] = (UShort_t)rand();  }
301      };  
302    };  void Digitizer::Close(){
303  };    delete fhBookTree;
304    }
305    
306  void Digitizer::DigitizePSCU(UInt_t length, UChar_t type) {  void Digitizer::Loop() {
307    //    //
308    UChar_t buff[16];    // opens the raw output file and loops over the events
309    //    //
310    // CPU signature    fOutputfile.open(fFilename, ios::out | ios::binary);
311    //      //fOutputfile.open(Form("Output%s",fFilename), ios::out | ios::binary);
312    buff[0] = 0xFA;    //
313    buff[1] = 0xFE;    // Load in memory and save at the beginning of file the calorimeter calibration
314    buff[2] = 0xDE;    //
315    //    CaloLoadCalib();
316    // packet type (twice)    DigitizeCALOCALIB();
317    //  
318    buff[3] = type;    //  load, digitize and write tracker calibration
319    buff[4] = type;    LoadTrackCalib();
320    //    
321    // counter    DigitizeTrackCalib(1);
322    //    UInt_t length=fTracklength*2;
323    fCounter++;    DigitizePSCU(length,0x12,fDataPSCU);
324    while ( fCounter > 16777215 ){    AddPadding();
325      fCounter -= 16777215;    WriteTrackCalib();
326    };    
327    //    DigitizeTrackCalib(2);
328    buff[5] = (UChar_t)(fCounter >> 16);    length=fTracklength*2;
329    buff[6] = (UChar_t)(fCounter >> 8);    DigitizePSCU(length,0x13,fDataPSCU);
330    buff[7] = (UChar_t)fCounter;    AddPadding();
331    //    WriteTrackCalib();
332    // on board time    LoadMipCor();  // some initialization of parameters -not used now-
333    //    //  end loading, digitizing and writing tracker calibration
334    ULong64_t obt = fOBT + 30LL;    // TOF ------ read calibration file (get A1, A2, lambda1, lambda2)
335    //    const int np=48;
336    while ( obt > 4294967295LL ){    float *atte1,*atte2,*lambda1,*lambda2;
337      obt -= 4294967295LL;    atte1=(float *)malloc(np *sizeof(float));
338    };    atte2=(float *)malloc(np *sizeof(float));
339    fOBT = (UInt_t)obt;    lambda1=(float *)malloc(np *sizeof(float));
340    //    lambda2=(float *)malloc(np *sizeof(float));
341    buff[8] = (UChar_t)(fOBT >> 24);    LoadTOFCalib(np,atte1,atte2,lambda1,lambda2);
342    buff[9] = (UChar_t)(fOBT >> 16);    attenAC = new TF1("fAttAC",".825+.64*atan(9.8/x)",0.,45.);
343    buff[10] = (UChar_t)(fOBT >> 8);    //end tof calib
344    buff[11] = (UChar_t)fOBT;    //
345    //    // loops over the events
346    // Packet length    //
347    //    
348    fLen = length;    Int_t nentries = fhBookTree->GetEntriesFast();
349    //    Long64_t nbytes = 0;
350    buff[12] = (UChar_t)(fLen >> 16);    for (Int_t i=0; i<nentries;i++) {
351    buff[13] = (UChar_t)(fLen >> 8);       nbytes += fhBookTree->GetEntry(i);
352    buff[14] = (UChar_t)fLen;       fEvent=i; //  cecilia for calo compress mode
353    //        // read detectors sequentially:
354    // CPU header CRC        // http://www.ts.infn.it/fileadmin/documents/physics/experiments/wizard/cpu/gen_arch/RM_Acquisition.pdf
355    //        // on pamelatov: /cvs/yoda/techmodel/physics/NeutronDetectorReader.cpp
356    buff[15] = (BYTE)CM_Compute_CRC16((UINT16)0, (BYTE*)&buff, (UINT32)15);       DigitizeTOF(np,atte1,atte2,lambda1,lambda2);
357    //        DigitizeAC();
358    memcpy(fDataPSCU,buff,16*sizeof(UChar_t));        DigitizeCALO();
359    //        DigitizeTrack();
360  };        DigitizeS4();
361          DigitizeND();
362  void Digitizer::ClearCaloCalib(Int_t s){        //
363    //        // Add padding to 64 bits
364    fcstwerr[s] = 0;        //
365    fcperror[s] = 0.;        AddPadding();
366    for ( Int_t d=0 ; d<11 ;d++  ){        //
367      Int_t pre = -1;        // Create CPU header, we need packet type (0x10 = physics data) and packet length.
368      for ( Int_t j=0; j<96 ;j++){        //
369        if ( j%16 == 0 ) pre++;        UInt_t length=2*(fCALOlength+fACbuffer+fTracklength+fNDbuffer+fS4buffer)+fPadding+fTOFbuffer+fTRIGGERbuffer;
370        fcalped[s][d][j] = 0.;        //UInt_t length=2*(fCALOlength+fACbuffer+fTracklength+fNDbuffer)+fPadding+fTOFbuffer+fTRIGGERbuffer;
371        fcstwerr[s] = 0.;        DigitizePSCU(length,0x10,fDataPSCU);
372        fcperror[s] = 0.;        if ((i%1000)==0)cout << "writing event " << i << endl;
373        fcalgood[s][d][j] = 0.;        WriteData();
374        fcalthr[s][d][pre] = 0.;    }
375        fcalrms[s][d][j] = 0.;    
376        fcalbase[s][d][pre] = 0.;    fOutputfile.close();
377        fcalvar[s][d][pre] = 0.;    cout << "files closed" << endl;
378      };  };
379    };  
380    return;  
381  }  void Digitizer::ReadData(){
382    
383  Int_t Digitizer::CaloLoadCalib(Int_t s,TString fcalname, UInt_t calibno){    UShort_t InData[64];
384    //  
385    //    // for debuggigng purposes only, write your own routine if you like (many
386    UInt_t e = 0;    // hardwired things.
387    if ( s == 0 ) e = 0;    
388    if ( s == 1 ) e = 2;    ifstream InputFile;
389    if ( s == 2 ) e = 3;  
390    if ( s == 3 ) e = 1;    // if (!InputFile) {
391    //  
392    ifstream myfile;    //     std::cout << "ERROR" << endl;
393    myfile.open(fcalname.Data());    //     // An error occurred!
394    if ( !myfile ){        //     // myFile.gcount() returns the number of bytes read.
395      return(-107);    //     // calling myFile.clear() will reset the stream state
396    };    //     // so it is usable again.
397    myfile.close();    //   };
398    //  
399    TFile *File = new TFile(fcalname.Data());  
400    if ( !File ) return(-108);    
401    TTree *tr = (TTree*)File->Get("CalibCalPed");    //InputFile.seekg(0);
402    if ( !tr ) return(-109);  
403    //    InputFile.open(fFilename, ios::in | ios::binary);
404    TBranch *calo = tr->GetBranch("CalibCalPed");    //    fOutputfile.seekg(0);
405    //    if (!InputFile.is_open()) std::cout << "ERROR" << endl;
406    pamela::CalibCalPedEvent *ce = 0;  
407    tr->SetBranchAddress("CalibCalPed", &ce);    InputFile.seekg(0);
408    //    
409    Long64_t ncalibs = calo->GetEntries();    for (Int_t k=0; k<=1000; k++){
410    //      InputFile.read(reinterpret_cast<char*>(InData),384*sizeof(UShort_t));
411    if ( !ncalibs ) return(-110);  
412    //      std::cout << "Read back: " << endl << endl;
413    calo->GetEntry(calibno);  
414    //      for (Int_t i=0; i<=384; i++){
415    if (ce->cstwerr[s] != 0 && ce->cperror[s] == 0 ) {        printf("%4x ", InData[i]);  
416      fcstwerr[s] = ce->cstwerr[s];        if ((i+1)%8 ==0) cout << endl;
417      fcperror[s] = ce->cperror[s];      }
418      for ( Int_t d=0 ; d<11 ;d++  ){  
419        Int_t pre = -1;    }
420        for ( Int_t j=0; j<96 ;j++){    cout << endl;
421          if ( j%16 == 0 ) pre++;    InputFile.close();
422          fcalped[s][d][j] = ce->calped[e][d][j];  
423          fcalgood[s][d][j] = ce->calgood[e][d][j];  };
         fcalthr[s][d][pre] = ce->calthr[e][d][pre];  
         fcalrms[s][d][j] = ce->calrms[e][d][j];  
         fcalbase[s][d][pre] = ce->calbase[e][d][pre];  
         fcalvar[s][d][pre] = ce->calvar[e][d][pre];  
       };  
     };  
   } else {  
     printf(" CALORIMETER - ERROR: problems finding a good calibration in this file! \n\n ");  
     File->Close();  
     return(-111);  
   };  
   File->Close();  
   return(0);  
 }  
   
   
 void Digitizer::DigitizeCALOCALIB() {  
   //  
   // Header of the four sections  
   //  
   fSecCalo[0] = 0xAA00; // XE  
   fSecCalo[1] = 0xB100; // XO  
   fSecCalo[2] = 0xB600; // YE  
   fSecCalo[3] = 0xAD00; // YO  
   //  
   // length of the data is 0x1215  
   //  
   fSecCALOLength[0] = 0x1215; // XE  
   fSecCALOLength[1] = 0x1215; // XO  
   fSecCALOLength[2] = 0x1215; // YE  
   fSecCALOLength[3] = 0x1215; // YO  
   //  
   Int_t chksum = 0;  
   UInt_t tstrip = 0;  
   UInt_t fSecPointer = 0;  
   //  
   for (Int_t sec=0; sec < 4; sec++){  
     //  
     // sec =    0 -> XE      1 -> XO        2-> YE         3 -> YO  
     //  
     fCALOlength = 0;  
     memset(fDataCALO,0,sizeof(UShort_t)*fCALObuffer);  
     fSecPointer = fCALOlength;  
     //  
     // First of all we have section header and packet length  
     //  
     fDataCALO[fCALOlength] = fSecCalo[sec];  
     fCALOlength++;  
     fDataCALO[fCALOlength] = fSecCALOLength[sec];  
     fCALOlength++;  
     //  
     // Section XO is read in the opposite direction respect to the others  
     //  
     chksum = 0;  
     //  
     for (Int_t plane=0; plane < 11; plane++){  
       //  
       if ( sec == 1 ) tstrip = fCALOlength + 96*2;  
       //  
       for (Int_t strip=0; strip < 96; strip++){    
         //  
         chksum += (Int_t)fcalped[sec][plane][strip];  
         //  
         // save value  
         //  
         if ( sec == 1 ){  
           tstrip -= 2;  
           fDataCALO[tstrip] = (Int_t)fcalped[sec][plane][strip];  
           fDataCALO[tstrip+1] = (Int_t)fcalgood[sec][plane][strip];  
         } else {  
           fDataCALO[fCALOlength] = (Int_t)fcalped[sec][plane][strip];  
           fDataCALO[fCALOlength+1] = (Int_t)fcalgood[sec][plane][strip];  
         };  
         fCALOlength +=2;  
       };  
       //  
     };  
     //  
     fDataCALO[fCALOlength] = (UShort_t)chksum;  
     fCALOlength++;  
     fDataCALO[fCALOlength] = 0;  
     fCALOlength++;  
     fDataCALO[fCALOlength] = (UShort_t)((Int_t)(chksum >> 16));  
     fCALOlength++;  
     //  
     // Section XO is read in the opposite direction respect to the others  
     //  
     chksum = 0;  
     //  
     for (Int_t plane=0; plane < 11; plane++){  
       //  
       if ( sec == 1 ) tstrip = fCALOlength+6*2;  
       //  
       for (Int_t strip=0; strip < 6; strip++){  
         //  
         chksum += (Int_t)fcalthr[sec][plane][strip];  
         //  
         // save value  
         //  
         if ( sec == 1 ){  
           tstrip -= 2;  
           fDataCALO[tstrip] = 0;  
           fDataCALO[tstrip+1] = (Int_t)fcalthr[sec][plane][strip];  
         } else {  
           fDataCALO[fCALOlength] = 0;  
           fDataCALO[fCALOlength+1] = (Int_t)fcalthr[sec][plane][strip];  
         };  
         fCALOlength +=2;  
       };  
       //  
     };  
     //  
     fDataCALO[fCALOlength] = 0;  
     fCALOlength++;  
     fDataCALO[fCALOlength] = (UShort_t)chksum;  
     fCALOlength++;  
     fDataCALO[fCALOlength] = 0;  
     fCALOlength++;  
     fDataCALO[fCALOlength] = (UShort_t)((Int_t)(chksum >> 16));  
     fCALOlength++;  
     //  
     // Section XO is read in the opposite direction respect to the others  
     //  
     for (Int_t plane=0; plane < 11; plane++){  
       //  
       if ( sec == 1 ) tstrip = fCALOlength+96*2;  
       //  
       for (Int_t strip=0; strip < 96; strip++){  
         //  
         // save value  
         //  
         if ( sec == 1 ){  
           tstrip -= 2;  
           fDataCALO[tstrip] = 0;  
           fDataCALO[tstrip+1] = (Int_t)fcalrms[sec][plane][strip];  
         } else {  
           fDataCALO[fCALOlength] = 0;  
           fDataCALO[fCALOlength+1] = (Int_t)fcalrms[sec][plane][strip];  
         };  
         fCALOlength += 2;  
       };  
       //  
     };      
     //  
     // Section XO is read in the opposite direction respect to the others  
     //  
     for (Int_t plane=0; plane < 11; plane++){  
       //  
       if ( sec == 1 ) tstrip = fCALOlength+6*4;  
       //  
       for (Int_t strip=0; strip < 6; strip++){  
         //  
         // save value  
         //  
         if ( sec == 1 ){  
           tstrip -= 4;  
           fDataCALO[tstrip] = 0;  
           fDataCALO[tstrip+1] = (Int_t)fcalbase[sec][plane][strip];  
           fDataCALO[tstrip+2] = 0;  
           fDataCALO[tstrip+3] = (Int_t)fcalvar[sec][plane][strip];  
         } else {  
           fDataCALO[fCALOlength] = 0;  
           fDataCALO[fCALOlength+1] = (Int_t)fcalbase[sec][plane][strip];  
           fDataCALO[fCALOlength+2] = 0;  
           fDataCALO[fCALOlength+3] = (Int_t)fcalvar[sec][plane][strip];  
         };  
         fCALOlength +=4;  
       };  
       //  
     };        
     //  
     //  
     // here we calculate and save the CRC  
     //  
     fDataCALO[fCALOlength] = 0;  
     fCALOlength++;  
     Short_t CRC = 0;  
     for (UInt_t i=0; i<(fCALOlength-fSecPointer); i++){  
       CRC=crc(CRC,fDataCALO[i+fSecPointer]);  
     };  
     fDataCALO[fCALOlength] = (UShort_t)CRC;  
     fCALOlength++;  
     //  
     UInt_t length=fCALOlength*2;  
     DigitizePSCU(length,0x18);  
     //  
     // Add padding to 64 bits  
     //  
     AddPadding();  
     //  
     fOutputfile.write(reinterpret_cast<char*>(fDataPSCU),sizeof(UShort_t)*fPSCUbuffer);  
     UShort_t temp[1000000];  
     memset(temp,0,sizeof(UShort_t)*1000000);  
     swab(fDataCALO,temp,sizeof(UShort_t)*fCALOlength);  // WE MUST SWAP THE BYTES!!!  
     fOutputfile.write(reinterpret_cast<char*>(temp),sizeof(UShort_t)*fCALOlength);  
     //  
     // padding to 64 bytes  
     //  
     if ( fPadding ){  
       fOutputfile.write(reinterpret_cast<char*>(fDataPadding),sizeof(UChar_t)*fPadding);  
     };  
     //  
     //      
   };  
   //  
 };  
   
 void Digitizer::CaloLoadCalib() {  
   //  
   fGivenCaloCalib = 0; //                                  ####@@@@ should be given as input par @@@@####  
   //  
   // first of all load the MIP to ADC conversion values  
   //  
   stringstream calfile;  
   Int_t error = 0;  
   GL_PARAM *glparam = new GL_PARAM();  
   //  
   // determine where I can find calorimeter ADC to MIP conversion file    
   //  
   error = 0;  
   error = glparam->Query_GL_PARAM(3,101,fDbc);  
   //  
   calfile.str("");  
   calfile << glparam->PATH.Data() << "/";  
   calfile << glparam->NAME.Data();  
   //  
   printf("\n Using Calorimeter ADC to MIP conversion file: \n %s \n",calfile.str().c_str());  
   FILE *f;  
   f = fopen(calfile.str().c_str(),"rb");  
   //  
   memset(fCalomip,0,4224*sizeof(fCalomip[0][0][0]));  
   //  
   for (Int_t m = 0; m < 2 ; m++ ){  
     for (Int_t k = 0; k < 22; k++ ){  
       for (Int_t l = 0; l < 96; l++ ){  
         fread(&fCalomip[m][k][l],sizeof(fCalomip[m][k][l]),1,f);  
       };  
     };  
   };  
   fclose(f);  
   //  
   // determine which calibration has to be used and load it for each section  
   //    
   GL_CALO_CALIB *glcalo = new GL_CALO_CALIB();  
   GL_ROOT *glroot = new GL_ROOT();    
   TString fcalname;  
   UInt_t idcalib;  
   UInt_t calibno;  
   UInt_t utime = 0;  
   //  
   for (UInt_t s=0; s<4; s++){  
     //  
     // clear calo calib variables for section s  
     //  
     ClearCaloCalib(s);  
     //  
     if ( fGivenCaloCalib ){  
       //  
       // a time has been given as input on the command line so retrieve the calibration that preceed that time  
       //  
       glcalo->Query_GL_CALO_CALIB(fGivenCaloCalib,utime,s,fDbc);  
       //    
       calibno = glcalo->EV_ROOT;  
       idcalib = glcalo->ID_ROOT_L0;  
       //  
       // determine path and name and entry of the calibration file  
       //  
       printf("\n");  
       printf(" ** SECTION %i **\n",s);  
       //  
       glroot->Query_GL_ROOT(idcalib,fDbc);  
       //  
       stringstream name;  
       name.str("");  
       name << glroot->PATH.Data() << "/";  
       name << glroot->NAME.Data();  
       //  
       fcalname = (TString)name.str().c_str();  
       //  
       printf("\n Section %i : using  file %s calibration at entry %i: \n",s,fcalname.Data(),calibno);  
       //  
     } else {  
       error = 0;  
       error = glparam->Query_GL_PARAM(1,104,fDbc);  
       //  
       calfile.str("");  
       calfile << glparam->PATH.Data() << "/";  
       calfile << glparam->NAME.Data();  
       //  
       printf("\n Section %i : using default calorimeter calibration: \n %s \n",s,calfile.str().c_str());  
       //  
       fcalname = (TString)calfile.str().c_str();  
       calibno = s;  
       //  
     };  
     //  
     // load calibration variables in memory  
     //  
     CaloLoadCalib(s,fcalname,calibno);  
     //  
   };  
   //  
   // at this point we have in memory the calorimeter calibration and we can save it to disk in the correct format and use it to digitize the data  
   //  
   delete glparam;  
   delete glcalo;  
   delete glroot;  
 };  
   
 void Digitizer::DigitizeCALO() {  
   //  
   fModCalo = 0; // 0 is RAW, 1 is COMPRESS, 2 is FULL     ####@@@@ should be given as input par @@@@####  
   //  
   //  
   //  
   fCALOlength = 0;  // reset total dimension of calo data  
   //  
   // gpamela variables to be used  
   //  
   fhBookTree->SetBranchStatus("Nthcali",1);  
   fhBookTree->SetBranchStatus("Icaplane",1);  
   fhBookTree->SetBranchStatus("Icastrip",1);  
   fhBookTree->SetBranchStatus("Icamod",1);  
   fhBookTree->SetBranchStatus("Enestrip",1);  
   //  
   // call different routines depending on the acq mode you want to simulate  
   //  
   switch ( fModCalo ){  
   case 0:  
     this->DigitizeCALORAW();  
     break;  
   case 1:  
     this->DigitizeCALOCOMPRESS();  
     break;  
   case 2:  
     this->DigitizeCALOFULL();  
     break;  
   };  
   //  
 };  
   
 Float_t Digitizer::GetCALOen(Int_t sec, Int_t plane, Int_t strip){  
   //  
   // determine plane and strip  
   //  
   Int_t mplane = 0;  
   //  
   // wrong!  
   //  
   //  if ( sec == 0 || sec == 3 ) mplane = (plane * 4) + sec + 1;    
   //  if ( sec == 1 ) mplane = (plane * 4) + 2 + 1;    
   //  if ( sec == 2 ) mplane = (plane * 4) + 1 + 1;    
   //  
   if ( sec == 0 ) mplane = plane * 4 + 1; // it must be 0, 4, 8, ... (+1)  from plane = 0, 11  
   if ( sec == 1 ) mplane = plane * 4 + 2 + 1; // it must be 2, 6, 10, ... (+1)  from plane = 0, 11  
   if ( sec == 2 ) mplane = plane * 4 + 3 + 1; // it must be 3, 7, 11, ... (+1)  from plane = 0, 11  
   if ( sec == 3 ) mplane = plane * 4 + 1 + 1; // it must be 1, 5, 9, ... (+1)  from plane = 0, 11  
   //  
   Int_t mstrip = strip + 1;  
   //  
   // search energy release in gpamela output  
   //  
   for (Int_t i=0; i<Nthcali;i++){  
     if ( Icaplane[i] == mplane && Icastrip[i] == mstrip ){  
       return (Enestrip[i]);  
     };  
   };  
   //  
   // if not found it means no energy release so return 0.  
   //  
   return(0.);  
 };  
   
 void Digitizer::DigitizeCALORAW() {  
   //  
   // some variables  
   //  
   Float_t ens = 0.;  
   UInt_t adcsig = 0;  
   UInt_t adcbase = 0;  
   UInt_t adc = 0;  
   Int_t pre = 0;  
   UInt_t l = 0;  
   UInt_t lpl = 0;  
   UInt_t tstrip = 0;  
   UInt_t fSecPointer = 0;  
   Double_t pedenoise;  
   Float_t rms = 0.;  
   Float_t pedestal = 0.;  
   //  
   // clean the data structure  
   //  
   memset(fDataCALO,0,sizeof(UShort_t)*fCALObuffer);  
   //  
   // Header of the four sections  
   //  
   fSecCalo[0] = 0xEA08; // XE  
   fSecCalo[1] = 0xF108; // XO  
   fSecCalo[2] = 0xF608; // YE  
   fSecCalo[3] = 0xED08; // YO  
   //  
   // length of the data is 0x0428 in RAW mode  
   //  
   fSecCALOLength[0] = 0x0428; // XE  
   fSecCALOLength[1] = 0x0428; // XO  
   fSecCALOLength[2] = 0x0428; // YE  
   fSecCALOLength[3] = 0x0428; // YO  
   //  
   // let's start  
   //  
   fCALOlength = 0;  
   //  
   for (Int_t sec=0; sec < 4; sec++){  
     //  
     // sec =    0 -> XE      1 -> XO        2-> YE         3 -> YO  
     //  
     l = 0;                 // XE and XO are Y planes  
     if ( sec < 2 ) l = 1;  // while YE and  YO are X planes  
     //  
     fSecPointer = fCALOlength;  
     //  
     // First of all we have section header and packet length  
     //  
     fDataCALO[fCALOlength] = fSecCalo[sec];  
     fCALOlength++;  
     fDataCALO[fCALOlength] = fSecCALOLength[sec];  
     fCALOlength++;  
     //  
     // selftrigger coincidences - in the future we should add here some code to simulate timing response of pre-amplifiers  
     //  
     for (Int_t autoplane=0; autoplane < 7; autoplane++){  
       fDataCALO[fCALOlength] = 0x0000;  
       fCALOlength++;  
     };  
     //  
     //  
     // here comes data  
     //  
     //  
     // Section XO is read in the opposite direction respect to the others  
     //  
     if ( sec == 1 ){        
       tstrip = 96*11 + fCALOlength;  
     } else {  
       tstrip = 0;  
     };  
     //  
     pre = -1;  
     //  
     for (Int_t strip=0; strip < 96; strip++){    
       //  
       // which is the pre for this strip?  
       //  
       if (strip%16 == 0) {  
         pre++;  
       };  
       //  
       if ( sec == 1 ) tstrip -= 11;  
       //  
       for (Int_t plane=0; plane < 11; plane++){  
         //  
         // here is wrong!!!!  
         //  
         //  
         //      if ( plane%2 == 0 && sec%2 != 0){  
         //        lpl = plane*2;  
         //      } else {  
         //        lpl = (plane*2) + 1;  
         //      };  
         //  
         if ( sec == 0 || sec == 3 ) lpl = plane * 2;  
         if ( sec == 1 || sec == 2 ) lpl = (plane * 2) + 1;  
         //  
         // get the energy in GeV from the simulation for that strip  
         //  
         ens = this->GetCALOen(sec,plane,strip);  
         //  
         // convert it into ADC channels  
         //        
         adcsig = int(ens*fCalomip[l][lpl][strip]/fCALOGeV2MIPratio);  
         //  
         // sum baselines  
         //  
         adcbase = (UInt_t)fcalbase[sec][plane][pre];  
         //  
         // add noise and pedestals  
         //        
         pedestal = fcalped[sec][plane][strip];  
         rms = fcalrms[sec][plane][strip]/4.;  
         //  
         // Add random gaussian noise of RMS rms and Centered in the pedestal  
         //  
         pedenoise = gRandom->Gaus((Double_t)pedestal,(Double_t)rms);  
         //  
         // Sum all contribution  
         //  
         adc = adcsig + adcbase + (Int_t)round(pedenoise);  
         //  
         // Signal saturation  
         //  
         if ( adc > 0x7FFF ) adc = 0x7FFF;  
         //  
         // save value  
         //  
         if ( sec == 1 ){  
           fDataCALO[tstrip] = adc;  
           tstrip++;  
         } else {  
           fDataCALO[fCALOlength] = adc;  
         };  
         fCALOlength++;  
         //  
       };  
       //  
       if ( sec == 1 ) tstrip -= 11;  
       //  
     };  
     //  
     // here we calculate and save the CRC  
     //  
     Short_t CRC = 0;  
     for (UInt_t i=0; i<(fCALOlength-fSecPointer); i++){  
       CRC=crc(CRC,fDataCALO[i+fSecPointer]);  
     };  
     fDataCALO[fCALOlength] = (UShort_t)CRC;  
     fCALOlength++;  
     //  
   };  
   //  
   //   for (Int_t i=0; i<fCALOlength; i++){  
   //     printf(" WORD %i       DIGIT %0x   \n",i,fDataCALO[i]);  
   //   };  
   //  
 };  
   
 void Digitizer::DigitizeCALOCOMPRESS() {  
   //  
   printf(" COMPRESS MODE STILL NOT IMPLEMENTED! \n");    
   //  
   this->DigitizeCALORAW();  
   return;  
   //  
   //  
   //  
   fSecCalo[0] = 0xEA00;  
   fSecCalo[1] = 0xF100;  
   fSecCalo[2] = 0xF600;  
   fSecCalo[3] = 0xED00;  
   //  
   // length of the data in DSP mode must be calculated on fly during digitization  
   //  
   memset(fSecCALOLength,0x0,4*sizeof(UShort_t));  
   //  
   // here comes raw data  
   //  
   Int_t en = 0;  
   //  
   for (Int_t sec=0; sec < 4; sec++){  
     fDataCALO[en] = fSecCalo[sec];  
     en++;  
     fDataCALO[en] = fSecCALOLength[sec];  
     en++;  
     for (Int_t plane=0; plane < 11; plane++){  
       for (Int_t strip=0; strip < 11; strip++){  
         fDataCALO[en] = 0x0;  
         en++;  
       };  
     };  
   };  
   //  
 };  
   
 void Digitizer::DigitizeCALOFULL() {  
   //  
   printf(" FULL MODE STILL NOT IMPLEMENTED! \n");    
   //  
   this->DigitizeCALORAW();  
   return;  
   //  
   fSecCalo[0] = 0xEA00;  
   fSecCalo[1] = 0xF100;  
   fSecCalo[2] = 0xF600;  
   fSecCalo[3] = 0xED00;  
   //  
   // length of the data in DSP mode must be calculated on fly during digitization  
   //  
   memset(fSecCALOLength,0x0,4*sizeof(UShort_t));  
   //  
   // here comes raw data  
   //  
   Int_t  en = 0;  
   //  
   for (Int_t sec=0; sec < 4; sec++){  
     fDataCALO[en] = fSecCalo[sec];  
     en++;  
     fDataCALO[en] = fSecCALOLength[sec];  
     en++;  
     for (Int_t plane=0; plane < 11; plane++){  
       for (Int_t strip=0; strip < 11; strip++){  
         fDataCALO[en] = 0x0;  
         en++;  
       };  
     };  
   };  
   //  
 };  
   
 void Digitizer::DigitizeTRIGGER() {  
   //fDataTrigger: 153 bytes  
   for (Int_t j=0; j < 153; j++)  
     fDataTrigger[j]=0x00;  
 };  
   
 Int_t Digitizer::DigitizeTOF() {  
   //fDataTof: 12 x 23 bytes (=276 bytes)  
   UChar_t *pTof=fDataTof;  
   Bool_t DEBUG=false;  
   
   // --- activate branches:  
   fhBookTree->SetBranchStatus("Nthtof",1);  
   fhBookTree->SetBranchStatus("Ipltof",1);  
   fhBookTree->SetBranchStatus("Ipaddle",1);  
   fhBookTree->SetBranchStatus("Xintof",1);  
   fhBookTree->SetBranchStatus("Yintof",1);  
   fhBookTree->SetBranchStatus("Xouttof",1);  
   fhBookTree->SetBranchStatus("Youttof",1);  
   fhBookTree->SetBranchStatus("Ereltof",1);  
   fhBookTree->SetBranchStatus("Timetof",1);  
   // not yet used: Zintof, Xouttof, Youttof, Zouttof  
   
   // ------ evaluate energy in each pmt: ------  
   // strip geometry (lenght/width)  
   Float_t dimel[6] = {33.0, 40.8 ,18.0, 15.0, 15.0, 18.0};  
   //Float_t dimes[6] = {5.1, 5.5, 7.5, 9.0, 6.0, 5.0};  
     
   //  S11 8 paddles  33.0 x 5.1 cm  
   //  S12 6 paddles  40.8 x 5.5 cm  
   //  S21 2 paddles  18.0 x 7.5 cm  
   //  S22 2 paddles  15.0 x 9.0 cm  
   //  S31 3 paddles  15.0 x 6.0 cm  
   //  S32 3 paddles  18.0 x 5.0 cm  
   
   Float_t FGeo[2]={0., 0.}; /* geometrical factor */  
   
   const Float_t Pho_keV = 10.;     // photons per keV in scintillator  
   const Float_t echarge = 1.6e-19; // electron charge  
   Float_t Npho=0.;  
   Float_t QevePmt_pC[48];  
   Float_t QhitPad_pC[2]={0., 0.};  
   Float_t QhitPmt_pC[2]={0., 0.};  
   Float_t pmGain = 3.5e6;  /* PMT Gain: the same for all PMTs */  
   Float_t effi=0.21;       /* Efficienza di fotocatodo */  
   
   //   Float_t ADC_pC0=-58.1;      //  ADC/pC conversion coefficient 0  
   //   Float_t ADC_pC1=1.728;      //  ADC/pC conversion coefficient 1  
   //   Float_t ADC_pC2=-4.063e-05; //  ADC/pC conversion coefficient 2  
   //   Float_t ADC_pC3=-5.763e-08; //  ADC/pC conversion coefficient 3  
     
   // pC < 800  
   Float_t ADC_pC0A =   -4.437616e+01   ;  
   Float_t ADC_pC1A =   1.573329e+00    ;  
   Float_t ADC_pC2A =   2.780518e-04  ;  
   Float_t ADC_pC3A =  -2.302160e-07 ;  
     
   // pC   > 800:  
   Float_t ADC_pC0B =    -2.245756e+02  ;  
   Float_t ADC_pC1B =     2.184156e+00  ;  
   Float_t ADC_pC2B =     -4.171825e-04  ;  
   Float_t ADC_pC3B =     3.789715e-08  ;  
   
   Float_t pCthres=40.;     // threshold in charge  
   Int_t ADClast=4095;      // no signal --> ADC ch=4095  
   Int_t ADCsat=3100;       // saturation value for the ADCs  
   Int_t ADCtof[48];  
   
   
   // ---- introduce scale factors to tune simul ADC to real data   24-oct DC  
   //   Float_t ScaleFact[48]={0.18,0.22,0.35,0.26,0.47,0.35,0.31,0.37,  
   //                     0.44,0.23,0.38,0.60,0.39,0.29,0.40,0.23,  
   //                     0.30,0.66,0.22,1.53,0.17,0.55,  
   //                     0.84,0.19,0.21,1.64,0.62,0.13,  
   //                     0.18,0.15,0.10,0.14,0.14,0.14,0.14,0.12,  
   //                     0.26,0.18,0.25,0.23,0.20,0.40,0.19,0.23,0.25,0.23,0.25,0.20};  
     
   // new scale factors: WM 30-Oct-07  
   //   Float_t ScaleFact[48]={0.35,0.41,0.32,0.34,0.58,0.47,0.42,0.44,  
   //                     0.50,0.34,0.50,0.50,0.51,0.42,0.46,0.25,  
   //                     0.20,0.38,0.29,0.49,0.24,0.68,  
   //                     0.30,0.26,0.28,0.79,0.31,0.12,  
   //                     0.25,0.21,0.14,0.20,  
   //                     0.16,0.17,0.19,0.18,  
   //                     0.34,0.27,0.34,0.31,0.25,0.57,  
   //                     0.24,0.34,0.34,0.32,0.31,0.30};  
     
   Float_t ScaleFact[48]={0.39, 0.49, 0.38, 0.40, 0.65, 0.51, 0.43,  
                          0.49, 0.58, 0.38, 0.53, 0.57, 0.53, 0.45, 0.49, 0.16,  
                          0.15, 0.44, 0.28, 0.57, 0.26, 0.72, 0.37, 0.29, 0.30, 0.89,  
                          0.37, 0.08,  0.27, 0.23, 0.12, 0.22, 0.15, 0.16, 0.21,  
                          0.19, 0.41, 0.32, 0.39, 0.38, 0.28, 0.66, 0.28, 0.40, 0.39, 0.40, 0.37, 0.35 };  
   
   for(Int_t i=0; i<48; i++){  
     QevePmt_pC[i] = 0;  
     ADCtof[i]=0;  
   }  
     
   // ------ read calibration file (get A1, A2, lambda1, lambda2)  
   ifstream fileTriggerCalib;  
   TString ftrigname="TrigCalibParam.txt";  
   fileTriggerCalib.open(ftrigname.Data());  
   if ( !fileTriggerCalib ) {  
     printf("debug: no trigger calib file!\n");  
     return(-117); //check output!  
   };  
   Float_t atte1[48],atte2[48],lambda1[48],lambda2[48];  
   Int_t temp=0;  
   // correct readout WM Oct '07  
   for(Int_t i=0; i<48; i++){  
     fileTriggerCalib >> temp;  
     fileTriggerCalib >> atte1[i];  
     fileTriggerCalib >> lambda1[i];  
     fileTriggerCalib >> atte2[i];  
     fileTriggerCalib >> lambda2[i];  
     fileTriggerCalib >> temp;  
   }  
   fileTriggerCalib.close();  
   
   Int_t ip, ipad;  
   //Int_t ipmt;  
   Int_t pmtleft=0, pmtright=0;  
   Int_t *pl, *pr;  
   pl = &pmtleft;  
   pr = &pmtright;  
   
   // TDC variables:  
   Int_t TDClast=4095;      // no signal --> TDC ch=4095  
   Int_t TDCint[48];  
   Float_t  tdc[48],tdc1[48],tdcpmt[48];  
   for(Int_t i=0; i<48; i++) {  
     tdcpmt[i] = 1000.;  
     tdc[i]  = 0.;            // 18-oct WM  
     tdc1[i] = 0.;            // 18-oct WM  
   }  
   
   Float_t thresh=10.; // to be defined better... (Wolfgang)  
   
   // === TDC: simulate timing for each paddle  
   //Float_t dt1 = 285.e-12 ;   // single PMT resolution  
   Float_t dt1 = 425.e-12 ;   // single PMT resolution (WM, Nov'07)  
   Float_t tdcres[50],c1_S[50],c2_S[50],c3_S[50];  
   for(Int_t j=0;j<48;j++)  tdcres[j] = 50.E-12;   // TDC resolution 50 picosec  
   for(Int_t j=0;j<48;j++)  c1_S[j] = 500.;  // cable length in channels  
   for(Int_t j=0;j<48;j++)  c2_S[j] = 0.;  
   for(Int_t j=0;j<48;j++)  c3_S[j] = 1000.;  
   for(Int_t j=0;j<48;j++)  c1_S[j] = c1_S[j]*tdcres[j];   // cable length in sec  
   for(Int_t j=0;j<48;j++)  c2_S[j] = c2_S[j]*tdcres[j];  
   // ih = 0 + i1;  // not used?? (Silvio)  
   
   /* **********************************  start loop over hits */  
     
   for(Int_t nh=0; nh<Nthtof; nh++){  
       
     Float_t s_l_g[6] = {8.0, 8.0, 20.9, 22.0, 9.8, 8.3 };  // length of the lightguide  
     Float_t t1,t2,veff,veff1,veff0 ;  
     veff0 = 100.*1.0e8 ; // light velocity in the scintillator in m/sec  
     veff1 = 100.*1.5e8; // light velocity in the lightguide in m/sec  
     veff=veff0;         // signal velocity in the paddle  
   
     t1 = Timetof[nh] ;  // Start  
     t2 = Timetof[nh] ;  
   
     // Donatella: redefinition plane and pad for vectors in C  
     ip = Ipltof[nh]-1;  
     ipad = Ipaddle[nh]-1;  
     pmtleft=0;  
     pmtright=0;  
   
     // WM: S12 paddles are "reversed" (Nov'07)  
     if (ip==2)  
       if (ipad==0)  
         ipad=1;  
       else  
         ipad=0;  
   
    if (ip<6) {  
       Paddle2Pmt(ip, ipad, &pmtleft, &pmtright);  
   
       // DC: evaluates mean position and path inside the paddle  
     
       Float_t tpos=0.;  
       Float_t path[2] = {0., 0.};  
       //--- Strip in Y = S11,S22,S31 ------  
       if(ip==0 || ip==3 || ip==4)  
         tpos = (Yintof[nh]+Youttof[nh])/2.;  
       else  
         if(ip==1 || ip==2 || ip==5)   //--- Strip in X for S12,S21,S32  
           tpos = (Xintof[nh]+Xouttof[nh])/2.;  
         else //if (ip!=6)  
           printf("*** WARNING TOF: this option should never occur! (ip=%2i, nh=%2i)\n",ip,nh);  
   
       path[0]= tpos + dimel[ip]/2.;   // path to left PMT  
       path[1]= dimel[ip]/2.- tpos;    // path to right PMT  
   
       //  cout <<"Strip N. ="<< ipaddle <<" piano n.= "<< iplane <<" POSIZ = "<< tpos <<"\n";  
         
       if (DEBUG) {  
         cout <<" plane "<<ip<<" strip # ="<< ipad <<" tpos  "<< tpos <<"\n";  
         cout <<"pmtleft, pmtright "<<pmtleft<<" "<<pmtright<<endl;  
       }  
         
       // constant geometric factor, the rest is handled by the scaling factor  
       FGeo[0] =0.5;  
       FGeo[1] =0.5;  
       //  FGeo[1] = atan(path[1]/dimes[ip])/6.28318; // fraction of photons toward left  
       //  FGeo[2] = atan(path[2]/dimes[ip])/6.28318; // toward right  
         
         
       //  Npho = Poisson(ERELTOF[nh])*Pho_keV*1e6  Poissonian fluctuations to be inserted-DC  
       Npho = Ereltof[nh]*Pho_keV*1.0e6;  // Eloss in GeV  
         
       Float_t knorm[2]={0., 0.}; // Donatella  
       Float_t Atten[2]={0., 0.}; // Donatella  
       for(Int_t j=0; j<2; j++){  
         QhitPad_pC[j]= Npho*FGeo[j]*effi*pmGain*echarge*1.E12*ScaleFact[pmtleft+j];  
         // WM  
         knorm[j]=atte1[pmtleft+j]*exp(lambda1[pmtleft+j]*dimel[ip]/2.*pow(-1,j+1)) +  
           atte2[pmtleft+j]*exp(lambda2[pmtleft+j]*dimel[ip]/2.*pow(-1,j+1));  
         Atten[j]=atte1[pmtleft+j]*exp(tpos*lambda1[pmtleft+j]) +  
           atte2[pmtleft+j]*exp(tpos*lambda2[pmtleft+j]) ;  
         QhitPmt_pC[j]= QhitPad_pC[j]*Atten[j]/knorm[j];  
         //              QhitPmt_pC[j]= QhitPad_pC[j]; //no attenuation  
   
   
         if (DEBUG) {  
           cout<<"pmtleft "<<pmtleft<<" j "<<j<<endl;  
           cout<<" atte1 "<<atte1[pmtleft+j]<<"lambda1 "<<lambda1[pmtleft+j]<<" atte2 "<<atte2[pmtleft+j]<<"lambda2 "<<lambda2[pmtleft+j] <<endl;      
           cout<<j<<" tpos "<<tpos<<" knorm "<<knorm[j]<<" "<<Atten[j]<<" "<<"QhitPmt_pC "<<QhitPmt_pC[j]<<endl;      
         }  
       }  
         
       if (DEBUG)  
         cout<<"Npho "<<Npho<<" QhitPmt_pC "<<QhitPmt_pC[0]<<" "<<QhitPmt_pC[1]<<endl;    
         
       QevePmt_pC[pmtleft]  += QhitPmt_pC[0];  
       QevePmt_pC[pmtright] += QhitPmt_pC[1];  
         
       // TDC  
       // WM right and left <->  
       //      t2 = t2 + fabs(path[0]/veff) + s_l_g[ip]/veff1 ;  // Signal reaches PMT  
       //      t1 = t1 + fabs(path[1]/veff) + s_l_g[ip]/veff1;  
   
       t1 = t1 + fabs(path[0]/veff) + s_l_g[ip]/veff1;  
       t2 = t2 + fabs(path[1]/veff) + s_l_g[ip]/veff1 ;  // Signal reaches PMT  
         
       t1 = gRandom->Gaus(t1,dt1); //apply gaussian error  dt  
       t2 = gRandom->Gaus(t2,dt1); //apply gaussian error  dt  
   
       t1 = t1 + c1_S[pmtleft] ;  // Signal reaches Discriminator ,TDC starts  to run  
       t2 = t2 + c1_S[pmtright] ;  
         
       // check if signal is above threshold  
       // then check if tdcpmt is already filled by another hit...  
       // only re-fill if time is smaller  
         
       if (QhitPmt_pC[0] > thresh) {  
         if (tdcpmt[pmtleft] == 1000.) {  // fill for the first time  
           tdcpmt[pmtleft] = t1;  
           tdc[pmtleft] = t1 + c2_S[pmtleft] ;  // Signal reaches Coincidence  
         }  
         if (tdcpmt[pmtleft] < 1000.) // is already filled!  
           if (t1 <  tdcpmt[pmtleft]) {  
             tdcpmt[pmtleft] = t1;  
             t1 = t1 + c2_S[pmtleft] ;  // Signal reaches Coincidence  
             tdc[pmtleft] = t1;  
           }  
       }        
       if (QhitPmt_pC[1] > thresh) {  
         if (tdcpmt[pmtright] == 1000.) {  // fill for the first time  
           tdcpmt[pmtright] = t2;  
           tdc[pmtright] = t2 + c2_S[pmtright] ;  // Signal reaches Coincidence  
         }  
         if (tdcpmt[pmtright] < 1000.)  // is already filled!  
           if (t2 <  tdcpmt[pmtright]) {  
             tdcpmt[pmtright] = t2;  
             t2 = t2 + c2_S[pmtright] ;  
             tdc[pmtright] = t2;  
           }        
       }  
   
       if (DEBUG)  
         cout<<nh<<" "<<Timetof[nh]<<" "<<t1<<" "<<t2<<endl;  
   
     } // ip < 6  
   
   }; // ****************************************       end loop over hits  
     
   // ======  ADC ======  
   
   for(Int_t i=0; i<48; i++){  
     if (QevePmt_pC[i] < 800.)  ADCtof[i]= (Int_t)(ADC_pC0A + ADC_pC1A*QevePmt_pC[i] + ADC_pC2A*pow(QevePmt_pC[i],2) + ADC_pC3A*pow(QevePmt_pC[i],3));  
     if (QevePmt_pC[i] > 800.)  ADCtof[i]= (Int_t)(ADC_pC0B + ADC_pC1B*QevePmt_pC[i] + ADC_pC2B*pow(QevePmt_pC[i],2) + ADC_pC3B*pow(QevePmt_pC[i],3));  
     if (QevePmt_pC[i] > 2485.) ADCtof[i]= (Int_t)(1758. + 0.54*QevePmt_pC[i]);  //assuming a fictional 0.54 ch/pC above ADCsat  
     if (ADCtof[i]>ADCsat) ADCtof[i]=ADCsat;  
     if (QevePmt_pC[i] < pCthres)  ADCtof[i]= ADClast;  
     if (ADCtof[i] < 0) ADCtof[i]=ADClast;      
     if (ADCtof[i] > ADClast) ADCtof[i]=ADClast;  
   }  
   
 //   for(Int_t i=0; i<48; i++){  
 //     if(QevePmt_pC[i] >= pCthres){  
 //       ADCtof[i]= (Int_t)(ADC_pC0 + ADC_pC1*QevePmt_pC[i] + ADC_pC2*pow(QevePmt_pC[i],2) + ADC_pC3*pow(QevePmt_pC[i],3));  
 //   } else  
 //       ADCtof[i]= ADClast;    
 // }  
     
 // // ---- introduce scale factors to tune simul ADC to real data   24-oct DC  
   
 //   for(Int_t i=0; i<48; i++){  
 //     if(ADCtof[i] != ADClast){  
 //                              //      printf("%3d, %4d, %4.2f\n",i, ADCtof[i],ScaleFact[i]);  
 //       ADCtof[i]= Int_t (ADCtof[i]*ScaleFact[i]);  
 //                              //      printf("%3d, %4d,\n",i, ADCtof[i]);  
 //     }  
 //   }  
   
 //   for(Int_t i=0; i<48; i++){  
 //     if(ADCtof[i] != ADClast){  
 //       if(ADCtof[i]> ADCsat) ADCtof[i]=ADCsat;  
 //       else if(ADCtof[i]< 0) ADCtof[i]=ADClast;      
 //   }  
 //   }      
   
   
 // ======  build  TDC coincidence  ======  
   
   Float_t t_coinc = 0;  
   Int_t ilast = 100;  
   for (Int_t ii=0; ii<48;ii++)  
     if (tdc[ii] > t_coinc) {  
       t_coinc = tdc[ii];  
       ilast = ii;  
     }  
     
   //     cout<<ilast<<" "<<t_coinc<<endl;  
   //     At t_coinc  trigger condition is fulfilled  
     
   for (Int_t ii=0; ii<48;ii++){  
     //      if (tdc[ii] != 0) tdc1[ii] = t_coinc - tdc[ii];   // test 1  
     if (tdc[ii] != 0) tdc1[ii] = t_coinc - tdcpmt[ii];  // test 2  
     tdc1[ii] = tdc1[ii]/tdcres[ii];                     // divide by TDC resolution  
     if (tdc[ii] != 0) tdc1[ii] = tdc1[ii] + c3_S[ii];  // add cable length c3  
   } // missing parenthesis inserted! (Silvio)  
   
   for(Int_t i=0; i<48; i++){  
     if(tdc1[i] != 0.){  
       TDCint[i]=(Int_t)tdc1[i];  
       if (TDCint[i]>4093) TDCint[i]=TDClast;  // 18-oct WM  
       if (DEBUG)  
         cout<<i<<" "<<TDCint[i]<<endl;  
       //ADC[i]= ADC_pC * QevePmt_pC[i] + ADCoffset;  
       //if(ADC[i]> ADClast) ADC[i]=ADClast;  
     } else  
       TDCint[i]= TDClast;  
   }  
   
   if (DEBUG)  
     cout<<"-----------"<<endl;  
   
   
 //------ use channelmap  18-oct WM  
   
   Int_t  channelmap[] =  {3,21,11,29,19,45,27,37,36,28,44,20,5,12,13,4,  
                         6,47,14,39,22,31,30,23,38,15,46,7,0,33,16,24,  
                         8,41,32,40,25,17,34,9,42,1,2,10,18,26,35,43};  
   
   Int_t ADChelp[48];  
   Int_t TDChelp[48];  
   
   for(Int_t i=0; i<48; i++){  
     Int_t ii=channelmap[i];  
     ADChelp[ii]= ADCtof[i];  
     TDChelp[ii]= TDCint[i];  
   }  
   
   for(Int_t i=0; i<48; i++){  
     ADCtof[i]= ADChelp[i];  
     TDCint[i]= TDChelp[i];  
   }  
   
   
 // ======  write fDataTof  =======  
   
   
 //  UChar_t tdcadd[8]={1,0,3,2,5,4,7,6};  (coded in 3 bit)  
   UChar_t Ctrl3bit[8]={32,0,96,64,160,128,224,192};  // DC (msb in 8 bit word )  
     
   UChar_t tofBin;  
   for (Int_t j=0; j < 12; j++){   // loop on TDC #12  
     Int_t j12=j*23;               // for each TDC 23 bytes (8 bits)  
     fDataTof[j12+0]=0x00;   // TDC_ID  
     fDataTof[j12+1]=0x00;   // EV_COUNT  
     fDataTof[j12+2]=0x00;   // TDC_MASK (1)  
     fDataTof[j12+3]=0x00;   // TDC_MASK (2)  
     for (Int_t k=0; k < 4; k++){   // for each TDC 4 channels (ADC+TDC)  
             
       Int_t jk12=j12+4*k;         // ADC,TDC channel (0-47)  
   
       tofBin =(UChar_t)(ADCtof[k+4*j]/256);   // ADC# (msb)  
       fDataTof[jk12+4] = Bin2GrayTof(tofBin,fDataTof[jk12+4]);  
       /* control bits inserted here, after the bin to gray conv - DC*/  
       fDataTof[jk12+4] = Ctrl3bit[2*k] | fDataTof[jk12+4];  
       tofBin=(UChar_t)(ADCtof[k+4*j]%256);   // ADC# (lsb)  
       fDataTof[jk12+5] = Bin2GrayTof(tofBin,fDataTof[jk12+5]);  
       tofBin=(UChar_t)(TDCint[k+4*j]/256);   // TDC# (msb)  
       fDataTof[jk12+6]=Bin2GrayTof(tofBin,fDataTof[jk12+6]);  
       /* control bits inserted here, after the bin to gray conv - DC*/  
       fDataTof[jk12+6] = Ctrl3bit[2*k+1] | fDataTof[jk12+6];  
       tofBin=(UChar_t)(TDCint[k+4*j]%256);   // TDC# (lsb)  
       fDataTof[jk12+7]=Bin2GrayTof(tofBin,fDataTof[jk12+7]);  
     };  
     fDataTof[j12+20]=0x00;   // TEMP1  
     fDataTof[j12+21]=0x00;   // TEMP2  
     fDataTof[j12+22]= EvaluateCrcTof(pTof);   // CRC  
     pTof+=23;  
   };  
   return(0);  
 };  
   
   
 UChar_t Digitizer::Bin2GrayTof(UChar_t binaTOF,UChar_t grayTOF){  
   union graytof_data {  
     UChar_t word;  
     struct bit_field {  
       unsigned b0:1;  
       unsigned b1:1;  
       unsigned b2:1;  
       unsigned b3:1;  
       unsigned b4:1;  
       unsigned b5:1;  
       unsigned b6:1;  
       unsigned b7:1;  
     } bit;  
   } bi,gr;  
   //  
   bi.word = binaTOF;  
   gr.word = grayTOF;  
   //  
   gr.bit.b0 = bi.bit.b1 ^ bi.bit.b0;  
   gr.bit.b1 = bi.bit.b2 ^ bi.bit.b1;  
   gr.bit.b2 = bi.bit.b3 ^ bi.bit.b2;  
   gr.bit.b3 = bi.bit.b3;  
   //  
   /* bin to gray conversion 4 bit per time*/  
   //  
   gr.bit.b4 = bi.bit.b5 ^ bi.bit.b4;  
   gr.bit.b5 = bi.bit.b6 ^ bi.bit.b5;  
   gr.bit.b6 = bi.bit.b7 ^ bi.bit.b6;  
   gr.bit.b7 = bi.bit.b7;  
   //  
   return(gr.word);  
 }  
   
 UChar_t Digitizer::EvaluateCrcTof(UChar_t *pTof) {  
   Bool_t DEBUG=false;  
   if (DEBUG)  
     return(0x00);  
   
   UChar_t crcTof=0x00;  
   UChar_t *pc=&crcTof, *pc2;  
   pc2=pTof;  
   for (Int_t jp=0; jp < 23; jp++){  
     //crcTof = crc8(...)  
     Crc8Tof(pc2++,pc);  
     //    printf("%2i --- %x\n",jp,crcTof);  
   }  
   return(crcTof);  
 }  
   
 void Digitizer::Crc8Tof(UChar_t *oldCRC, UChar_t *crcTof){  
   union crctof_data {  
     UChar_t word;  
     struct bit_field {  
       unsigned b0:1;  
       unsigned b1:1;  
       unsigned b2:1;  
       unsigned b3:1;  
       unsigned b4:1;  
       unsigned b5:1;  
       unsigned b6:1;  
       unsigned b7:1;  
     } bit;  
   } c,d,r;  
   
   c.word = *oldCRC;  
   //d.word = *newCRC;  
   d.word = *crcTof;  
   r.word = 0;  
   
   r.bit.b0 = c.bit.b7 ^ c.bit.b6 ^ c.bit.b0 ^  
              d.bit.b0 ^ d.bit.b6 ^ d.bit.b7;  
   
   r.bit.b1 = c.bit.b6 ^ c.bit.b1 ^ c.bit.b0 ^  
              d.bit.b0 ^ d.bit.b1 ^ d.bit.b6;  
   
   r.bit.b2 = c.bit.b6 ^ c.bit.b2 ^ c.bit.b1 ^ c.bit.b0 ^  
              d.bit.b0 ^ d.bit.b1 ^ d.bit.b2 ^ d.bit.b6;  
   
   r.bit.b3 = c.bit.b7 ^ c.bit.b3 ^ c.bit.b2 ^ c.bit.b1 ^  
              d.bit.b1 ^ d.bit.b2 ^ d.bit.b3 ^ d.bit.b7;  
   
   r.bit.b4 = c.bit.b4 ^ c.bit.b3 ^ c.bit.b2 ^  
              d.bit.b2 ^ d.bit.b3 ^ d.bit.b4;  
   
   r.bit.b5 = c.bit.b5 ^ c.bit.b4 ^ c.bit.b3 ^  
              d.bit.b3 ^ d.bit.b4 ^ d.bit.b5;  
   
   r.bit.b6 = c.bit.b6 ^ c.bit.b5 ^ c.bit.b4 ^  
              d.bit.b4 ^ d.bit.b5 ^ d.bit.b6;  
   
   r.bit.b7 = c.bit.b7 ^ c.bit.b6 ^ c.bit.b5 ^  
              d.bit.b5 ^ d.bit.b6 ^ d.bit.b7 ;  
   
   *crcTof=r.word;  
   //return r.word;  
 };  
   
 //void Digitizer::Paddle2Pmt(Int_t plane, Int_t paddle, Int_t* &pmtleft, Int_t* &pmtright){  
 void Digitizer::Paddle2Pmt(Int_t plane, Int_t paddle, Int_t *pl, Int_t *pr){  
   //* @param plane    (0 - 5)  
   //* @param paddle   (plane=0, paddle = 0,...5)  
   //* @param padid    (0 - 23)    
   //  
   Int_t padid=-1;  
   Int_t pads[6]={8,6,2,2,3,3};  
   //  
   Int_t somma=0;  
   Int_t np=plane;  
   for(Int_t j=0; j<np; j++)  
     somma+=pads[j];  
   padid=paddle+somma;  
   *pl = padid*2;  
   //  *pr = *pr + 1;  
   *pr = *pl + 1; // WM  
 };  
   
 void Digitizer::DigitizeAC() {  
   // created:  J. Conrad, KTH  
   // modified: S. Orsi, INFN Roma2  
   // fDataAC[0-63]:   main AC board  
   // fDataAC[64-127]: extra AC board  
   
   fDataAC[0] = 0xACAC;  
   fDataAC[64]= 0xACAC;  
   fDataAC[1] = 0xAC11;    
   fDataAC[65] = 0xAC22;    
   
   // the third word is a status word (dummy: "no errors are present in the AC boards")  
   fDataAC[2] = 0xFFFF; //FFEF?  
   fDataAC[66] = 0xFFFF;  
   
   const UInt_t nReg = 6;  
   
   // FPGA Registers (dummy)  
   for (UInt_t i=0; i<=nReg; i++){  
     fDataAC[i+4] = 0xFFFF;  
     fDataAC[i+68] = 0xFFFF;  
   }  
   
   // the last word is a CRC  
   // Dummy for the time being, but it might need to be calculated in the end  
   fDataAC[63] = 0xABCD;  
   fDataAC[127] = 0xABCD;  
   
   // shift registers (moved to the end of the routine)  
   
   Int_t evntLSB=Ievnt%65536;  
   Int_t evntMSB=(Int_t)(Ievnt/65536);  
   
   // singles counters are dummy  
   for (UInt_t i=0; i<=15; i++){  //SO Oct '07: // for (UInt_t i=0; i<=16; i++){  
     //     fDataAC[i+26] = 0x0000;    
     //     fDataAC[i+90] = 0x0000;  
     fDataAC[i+26] = evntLSB;  
     fDataAC[i+90] = evntLSB;  
   };  
     
   // coincidences are dummy (increment by 1 at each event)  
   // for (UInt_t i=0; i<=7; i++){  
   //    fDataAC[i+42] = 0x0000;  
   //    fDataAC[i+106] = 0x0000;  
   //   }  
   for (UInt_t i=0; i<=7; i++){  
     fDataAC[i+42] = evntLSB;  
     fDataAC[i+106] = evntLSB;  
   };  
   
   // increments for every trigger might be needed at some point.  
   // dummy for now  
   fDataAC[50]  = 0x0000;  
   fDataAC[114] = 0x0000;  
   
   // dummy FPGA clock (increment by 1 at each event)  
   /*      
     fDataAC[51] = 0x006C;  
     fDataAC[52] = 0x6C6C;  
     fDataAC[115] = 0x006C;  
     fDataAC[116] = 0x6C6C;  
   */  
   if (Ievnt<=0xFFFF) {  
     fDataAC[51] = 0x0000;  
     fDataAC[52] = Ievnt;  
     fDataAC[115] = 0x0000;  
     fDataAC[116] = Ievnt;  
   } else {  
     fDataAC[51] = evntMSB;  
     fDataAC[52] = evntLSB;  
     fDataAC[115] = fDataAC[51];  
     fDataAC[116] = fDataAC[52];  
   }  
   
   // dummy temperatures  
   fDataAC[53] = 0x0000;  
   fDataAC[54] = 0x0000;  
   fDataAC[117] = 0x0000;  
   fDataAC[118] = 0x0000;  
   
   
   // dummy DAC thresholds  
   for (UInt_t i=0; i<=7; i++){  
     fDataAC[i+55] = 0x1A13;    
     fDataAC[i+119] = 0x1A13;  
   }  
     
   // We activate all branches. Once the digitization algorithm is determined  
   // only the branches that involve needed information will be activated  
     
   fhBookTree->SetBranchAddress("Ievnt",&Ievnt);  
   fhBookTree->SetBranchStatus("Nthcat",1);  
   fhBookTree->SetBranchStatus("Iparcat",1);  
   fhBookTree->SetBranchStatus("Icat",1);  
   fhBookTree->SetBranchStatus("Xincat",1);  
   fhBookTree->SetBranchStatus("Yincat",1);  
   fhBookTree->SetBranchStatus("Zincat",1);  
   fhBookTree->SetBranchStatus("Xoutcat",1);  
   fhBookTree->SetBranchStatus("Youtcat",1);  
   fhBookTree->SetBranchStatus("Zoutcat",1);  
   fhBookTree->SetBranchStatus("Erelcat",1);  
   fhBookTree->SetBranchStatus("Timecat",1);  
   fhBookTree->SetBranchStatus("Pathcat",1);  
   fhBookTree->SetBranchStatus("P0cat",1);  
   fhBookTree->SetBranchStatus("Nthcas",1);  
   fhBookTree->SetBranchStatus("Iparcas",1);  
   fhBookTree->SetBranchStatus("Icas",1);  
   fhBookTree->SetBranchStatus("Xincas",1);  
   fhBookTree->SetBranchStatus("Yincas",1);  
   fhBookTree->SetBranchStatus("Zincas",1);  
   fhBookTree->SetBranchStatus("Xoutcas",1);  
   fhBookTree->SetBranchStatus("Youtcas",1);  
   fhBookTree->SetBranchStatus("Zoutcas",1);  
   fhBookTree->SetBranchStatus("Erelcas",1);  
   fhBookTree->SetBranchStatus("Timecas",1);  
   fhBookTree->SetBranchStatus("Pathcas",1);  
   fhBookTree->SetBranchStatus("P0cas",1);  
   fhBookTree->SetBranchStatus("Nthcard",1);  
   fhBookTree->SetBranchStatus("Iparcard",1);  
   fhBookTree->SetBranchStatus("Icard",1);  
   fhBookTree->SetBranchStatus("Xincard",1);  
   fhBookTree->SetBranchStatus("Yincard",1);  
   fhBookTree->SetBranchStatus("Zincard",1);  
   fhBookTree->SetBranchStatus("Xoutcard",1);  
   fhBookTree->SetBranchStatus("Youtcard",1);  
   fhBookTree->SetBranchStatus("Zoutcard",1);  
   fhBookTree->SetBranchStatus("Erelcard",1);  
   fhBookTree->SetBranchStatus("Timecard",1);  
   fhBookTree->SetBranchStatus("Pathcard",1);  
   fhBookTree->SetBranchStatus("P0card",1);  
   
   // In this simpliefied approach we will assume that once  
   // a particle releases > 0.5 mip in one of the 12 AC detectors it  
   // will fire. We will furthermore assume that both cards read out  
   // identical data.  
   
   // If you develop your digitization algorithm, you should start by  
   // identifying the information present in level2 (post-darth-vader)  
   // data.  
   
   Float_t SumEcat[5];  
   Float_t SumEcas[5];  
   Float_t SumEcard[5];  
   for (Int_t k= 0;k<5;k++){  
     SumEcat[k]=0.;  
     SumEcas[k]=0.;  
     SumEcard[k]=0.;  
   };  
   
   if (Nthcat>50 || Nthcas>50 || Nthcard>50)  
     printf("*** ERROR AC! NthAC out of range!\n\n");  
   
   // energy dependence on position (see file AcFitOutputDistancePmt.C by S.Orsi)  
   // based on J.Lundquist's calculations (PhD thesis, page 94)  
   // function: [0]+[1]*atan([2]/(x+1)), where the 3 parameters are:  
   //   8.25470e-01   +- 1.79489e-02  
   //   6.41609e-01   +- 2.65846e-02  
   //   9.81177e+00   +- 1.21284e+00  
   // hp: 1 minimum ionising particle at 35cm from the PMT releases 1mip  
   //  
   // NB: the PMT positions are needed!  
   
   // look in CAT  
   //  for (UInt_t k= 0;k<50;k++){  
   for (Int_t k= 0;k<Nthcat;k++){  
     if (Erelcat[k] > 0)  
       SumEcat[Icat[k]] += Erelcat[k];  
   };  
   
   // look in CAS  
   for (Int_t k= 0;k<Nthcas;k++){  
     if (Erelcas[k] >0)  
       SumEcas[Icas[k]] += Erelcas[k];  
   };  
   
   // look in CARD  
   for (Int_t k= 0;k<Nthcard;k++){  
     if (Erelcard[k] >0)  
       SumEcard[Icard[k]] += Erelcard[k];  
   };  
   
   // channel mapping              Hit Map  
   // 1 CARD4                      0          LSB  
   // 2 CAT2                       0  
   // 3 CAS1                       0  
   // 4 NC                         0  
   // 5 CARD2                      0  
   // 6 CAT4                       1  
   // 7 CAS4                       0    
   // 8 NC                         0  
   // 9 CARD3                      0  
   // 10 CAT3                      0  
   // 11 CAS3                      0  
   // 12 NC                        0  
   // 13 CARD1                     0  
   // 14 CAT1                      0  
   // 15 CAS2                      0  
   // 16 NC                        0          MSB  
   
   // In the first version only the hit-map is filled, not the SR.  
   
   // Threshold: 0.8 MeV.  
   
   Float_t thr = 8e-4;  
   
   fDataAC[3] = 0x0000;  
   
   if (SumEcas[0] > thr)  fDataAC[3]  = 0x0004;  
   if (SumEcas[1] > thr)  fDataAC[3] += 0x4000;  
   if (SumEcas[2] > thr)  fDataAC[3] += 0x0400;  
   if (SumEcas[3] > thr)  fDataAC[3] += 0x0040;    
   
   if (SumEcat[0] > thr)  fDataAC[3] += 0x2000;  
   if (SumEcat[1] > thr)  fDataAC[3] += 0x0002;  
   if (SumEcat[2] > thr)  fDataAC[3] += 0x0200;  
   if (SumEcat[3] > thr)  fDataAC[3] += 0x0020;  
     
   if (SumEcard[0] > thr)  fDataAC[3] += 0x1000;  
   if (SumEcard[1] > thr)  fDataAC[3] += 0x0010;  
   if (SumEcard[2] > thr)  fDataAC[3] += 0x0100;  
   if (SumEcard[3] > thr)  fDataAC[3] += 0x0001;  
   
   fDataAC[67] = fDataAC[3];  
   
   // shift registers  
   // the central bin is equal to the hitmap, all other bins in the shift register are 0  
   for (UInt_t i=0; i<=15; i++){  
     fDataAC[i+11] = 0x0000;    
     fDataAC[i+75] = 0x0000;  
   }  
   fDataAC[18] = fDataAC[3];  
   fDataAC[82] = fDataAC[3];  
   
   //    for (Int_t i=0; i<fACbuffer; i++){  
   //    printf("%0x  ",fDataAC[i]);    
   //    if ((i+1)%8 ==0) cout << endl;  
   //   }  
 };  
   
   
 void Digitizer::DigitizeS4(){  
   Int_t DEBUG=0;  
   // creato:  S. Borisov, INFN Roma2 e MEPHI, Sett 2007  
   TString ciao,modo="ns";  
   Int_t i,j,t,NdF,pmt,NdFT,S4,S4v=0,S4p=32;  
   Float_t E0,E1=1e-6,Ert,X,Y,Z,x,y,z,V[3],Xs[2],Ys[2],Zs[2],Yp[6],q,w,p=0.1,l,l0=500;  
   Xs[0]=-24.1;  
   Xs[1]=24.1;  
   Ys[0]=-24.1;  
   Ys[1]=24.1;  
   Zs[0]=-0.5;  
   Zs[1]=0.5;  
   Yp[0]=-20.;  
   Yp[2]=-1.;  
   Yp[4]=17.;  
   for(i=0;i<3;i++)  
     Yp[2*i+1]=Yp[2*i]+3;  
   srand(time(NULL));  
   // --- activate branches:  
   fhBookTree->SetBranchStatus("Nthtof",1);  
   fhBookTree->SetBranchStatus("Ipltof",1);  
   fhBookTree->SetBranchStatus("Ipaddle",1);  
     
   fhBookTree->SetBranchStatus("Xintof",1);  
   fhBookTree->SetBranchStatus("Yintof",1);  
   fhBookTree->SetBranchStatus("Xouttof",1);  
   fhBookTree->SetBranchStatus("Youttof",1);  
     
   fhBookTree->SetBranchStatus("Ereltof",1);  
   fhBookTree->SetBranchStatus("Timetof",1);  
   NdFT=0;  
   Ert=0;  
   for(i=0;i<Nthtof;i++){  
     if(Ipltof[i]!=6) continue;  
     Ert+=Ereltof[i];  
             
     if(modo=="ns") continue;  
     NdF=Int_t(Ereltof[i]/E1);  
     NdFT=0;  
     X=Xintof[i];  
     Y=Yintof[i];  
     Z=(Float_t)(random())/(Float_t)(0x7fffffff)-0.5;  
     //cout<<"XYZ "<<X<<"  "<<Y<<"  "<<Z<<endl;  
     for(j=0;j<NdF;j++){  
       q=(Float_t)random()/(Float_t)0x7fffffff;  
       w=(Float_t)random()/(Float_t)0x7fffffff;  
       // cout<<"qw "<<q<<" "<<w<<endl;  
       V[0]=p*cos(6.28318*q);  
       V[1]=p*sin(6.28318*q);  
       V[2]=p*(2.*w-1.);  
       pmt=0;  
       x=X;  
       y=Y;  
       z=Z;  
       while(pmt==0 && (x>Xs[0] && x<Xs[1])&&(y>Ys[0] && y<Ys[1])&&(z>Zs[0] && z<Zs[1])){  
         l=0;  
         while(pmt==0 && (x>Xs[0] && x<Xs[1])&&(y>Ys[0] && y<Ys[1])&&(z>Zs[0] && z<Zs[1])){  
           x+=V[0];  
           y+=V[1];  
           z+=V[2];  
           l+=p;  
           //cout<<x<<"  "<<y<<"  "<<z<<"  "<<l<<endl;  
           //cin>>ciao;  
         }  
         if((x<Xs[0]+p || x>Xs[1]-p)&&(y>Ys[0]+p && y<Ys[1]-p)&&(z>Zs[0]+p && z<Zs[1]-p)){  
           for(t=0;t<3;t++){  
             if(y>=Yp[2*t] && y<Yp[2*t+1]){  
               if(pmt==0)NdFT++;  
               pmt=1;  
               //cout<<NdFT<<endl;  
               break;  
             }  
           }  
           if(pmt==1)break;  
           V[0]=-V[0];  
         }  
         q=(Float_t)random()/(Float_t)0x7fffffff;  
         w=1-exp(-l/l0);  
         if(q<w)break;  
         q=(Float_t)random()/(Float_t)0x7fffffff;  
         w=0.5;  
         if(q<w)break;  
         if((x>Xs[0]+p && x<Xs[1]-p)&&(y<Ys[0]+p || y>Ys[1]-p)&&(z>Zs[0]+p && z<Zs[1]-p))V[1]=-V[1];  
         if((x>Xs[0]+p && x<Xs[1]-p)&&(y>Ys[0]+p && y<Ys[1]-p)&&(z<Zs[0]+p || z>Zs[1]-p))V[2]=-V[2];  
         x+=V[0];  
         y+=V[1];  
         z+=V[2];  
         l=0;  
         //cout<<x<<"  "<<y<<"  "<<z<<"  "<<l<<endl;  
                 //cin>>ciao;  
       }  
     }  
   }  
   Ert=Ert/0.002;  
   q=(Float_t)(random())/(Float_t)0x7fffffff;  
   w=0.7;  
   //E0=(Float_t)(4064./7.);  
   E0=4064./7.;  
   if(Ert<1) S4=0;  
   else S4=(Int_t)(4064.*(1.-exp(-(Ert-1.)/E0)));  
   i=S4/4;  
   if(S4%4==0)  
     S4v=S4+S4p;  
   else if(S4%4==1){  
     if(q<w) S4v=S4-1+S4p;  
     else S4v=S4+1+S4p;  
   } else if(S4%4==2) S4v=S4+S4p;  
   else if(S4%4==3){  
     if(q<w) S4v=S4+1+S4p;  
     else S4v=S4-1+S4p;  
   }  
   if (DEBUG)  
     cout<<"Ert_S4 = " << Ert << " --- S4v = " << S4v << endl;  
   fDataS4[0]=S4v;//0xf028;  
   fDataS4[1]=0xd800;  
   fDataS4[2]=0x0300;  
   //cout<<"  PMT  "<<NdFT<<"  "<<NdF<<endl;  
   //cin>>ciao;  
 }  
   
   
   
 void Digitizer::DigitizeND(){  
   // creato:  S. Borisov, INFN Roma2 e MEPHI, Sett 2007  
   Int_t i=0;  
   UShort_t NdN=0;  
   fhBookTree->SetBranchStatus("Nthnd",1);  
   fhBookTree->SetBranchStatus("Itubend",1);  
   fhBookTree->SetBranchStatus("Iparnd",1);    
   fhBookTree->SetBranchStatus("Xinnd",1);  
   fhBookTree->SetBranchStatus("Yinnd",1);  
   fhBookTree->SetBranchStatus("Zinnd",1);  
   fhBookTree->SetBranchStatus("Xoutnd",1);  
   fhBookTree->SetBranchStatus("Youtnd",1);  
   fhBookTree->SetBranchStatus("Zoutnd",1);  
   fhBookTree->SetBranchStatus("Erelnd",1);  
   fhBookTree->SetBranchStatus("Timend",1);  
   fhBookTree->SetBranchStatus("Pathnd",1);  
   fhBookTree->SetBranchStatus("P0nd",1);  
   //cout<<"n="<<Nthnd<<"  "<<NdN<<"\n";  
   for(i=0;i<Nthnd;i++){  
     if(Iparnd[i]==13){  
       NdN++;  
     }  
   }  
   //NdN=100; //only for debug  
   
   for(i=0;i<3;i++){  
     fDataND[2*i]=0x0000;  
     fDataND[2*i+1]=0x010F;  
   }  
   fDataND[0]=0xFF00 & (256*NdN);  
 }  
   
   
 void Digitizer::DigitizeDummy() {  
   
   fhBookTree->SetBranchStatus("Enestrip",1);  
   
   // dumy header  
   fDataDummy[0] = 0xCAAA;  
   
   for (Int_t i=1; i<fDummybuffer; i++){  
     fDataDummy[i] = 0xFFFF;  
     //   printf("%0x  ",fDataDummy[i]);    
     //if ((i+1)%8 ==0) cout << endl;  
   }  
 };  
   
   
 void Digitizer::WriteData(){  
   
   // Routine that writes the data to a binary file  
   // PSCU data are already swapped  
   fOutputfile.write(reinterpret_cast<char*>(fDataPSCU),sizeof(UShort_t)*fPSCUbuffer);  
   // TRG  
   fOutputfile.write(reinterpret_cast<char*>(fDataTrigger),sizeof(UChar_t)*153);  
   // TOF  
   fOutputfile.write(reinterpret_cast<char*>(fDataTof),sizeof(UChar_t)*276);  
   // AC  
   UShort_t temp[1000000];  
   memset(temp,0,sizeof(UShort_t)*1000000);  
   swab(fDataAC,temp,sizeof(UShort_t)*fACbuffer);  // WE MUST SWAP THE BYTES!!!  
   fOutputfile.write(reinterpret_cast<char*>(temp),sizeof(UShort_t)*fACbuffer);  
   // CALO  
   memset(temp,0,sizeof(UShort_t)*1000000);  
   swab(fDataCALO,temp,sizeof(UShort_t)*fCALOlength); // WE MUST SWAP THE BYTES!!!  
   fOutputfile.write(reinterpret_cast<char*>(temp),sizeof(UShort_t)*fCALOlength);  
   // TRK  
   memset(temp,0,sizeof(UShort_t)*1000000);  
   swab(fDataTrack,temp,sizeof(UShort_t)*fTracklength);  // WE MUST SWAP THE BYTES!!!  
   fOutputfile.write(reinterpret_cast<char*>(temp),sizeof(UShort_t)*fTracklength);  
   fTracklength=0;  
    // padding to 64 bytes  
   //  
   if ( fPadding ){  
     fOutputfile.write(reinterpret_cast<char*>(fDataPadding),sizeof(UChar_t)*fPadding);  
   };  
   // S4  
   memset(temp,0,sizeof(UShort_t)*1000000);  
   swab(fDataS4,temp,sizeof(UShort_t)*fS4buffer);  // WE MUST SWAP THE BYTES!!!  
   fOutputfile.write(reinterpret_cast<char*>(temp),sizeof(UShort_t)*fS4buffer);  
   // ND  
   memset(temp,0,sizeof(UShort_t)*1000000);  
   swab(fDataND,temp,sizeof(UShort_t)*fNDbuffer);  // WE MUST SWAP THE BYTES!!!  
   fOutputfile.write(reinterpret_cast<char*>(temp),sizeof(UShort_t)*fNDbuffer);  
 };  
   
   
 void Digitizer::ReadData(){  
   
   UShort_t InData[64];  
   
   // for debuggigng purposes only, write your own routine if you like (many  
   // hardwired things.  
     
   ifstream InputFile;  
   
   // if (!InputFile) {  
   
   //     std::cout << "ERROR" << endl;  
   //     // An error occurred!  
   //     // myFile.gcount() returns the number of bytes read.  
   //     // calling myFile.clear() will reset the stream state  
   //     // so it is usable again.  
   //   };  
   
   
     
   //InputFile.seekg(0);  
   
   InputFile.open(fFilename, ios::in | ios::binary);  
   //    fOutputfile.seekg(0);  
   if (!InputFile.is_open()) std::cout << "ERROR" << endl;  
   
   InputFile.seekg(0);  
     
   for (Int_t k=0; k<=1000; k++){  
     InputFile.read(reinterpret_cast<char*>(InData),384*sizeof(UShort_t));  
   
     std::cout << "Read back: " << endl << endl;  
   
     for (Int_t i=0; i<=384; i++){  
       printf("%4x ", InData[i]);    
       if ((i+1)%8 ==0) cout << endl;  
     }  
   
   }  
   cout << endl;  
   InputFile.close();  
   
 };  
   
   
   
 void Digitizer::DigitizeTrack() {  
 //std:: cout << "Entering DigitizeTrack " << endl;  
 Float_t  AdcTrack[fNviews][fNstrips_view];  //  Vector of strips to be compressed  
   
 Int_t Iview;  
 Int_t Nstrip;  
   
   for (Int_t j=0; j<fNviews;j++) {  
   
     for (Int_t i=0; i<fNladder;i++) {  
   
       Float_t commonN1=gRandom->Gaus(0.,fSigmaCommon);  
       Float_t commonN2=gRandom->Gaus(0.,fSigmaCommon);  
       for (Int_t k=0; k<fNstrips_ladder;k++) {  
       Nstrip=i*fNstrips_ladder+k;  
       AdcTrack[j][Nstrip]=gRandom->Gaus(fPedeTrack[j][Nstrip],fSigmaTrack[j][Nstrip]);  
       if(k<4*128) {AdcTrack[j][Nstrip] += commonN1;}  // full correlation of 4 VA1 Com. Noise  
       else {AdcTrack[j][Nstrip] += commonN2;}   // full correlation of 4 VA1 Com. Noise  
   
       };  
         
           
     };  
   
   
   };  
   
   
   fhBookTree->SetBranchStatus("Nstrpx",1);  
   fhBookTree->SetBranchStatus("Npstripx",1);  
   fhBookTree->SetBranchStatus("Ntstripx",1);  
   fhBookTree->SetBranchStatus("Istripx",1);  
   fhBookTree->SetBranchStatus("Qstripx",1);  
   fhBookTree->SetBranchStatus("Xstripx",1);  
   fhBookTree->SetBranchStatus("Nstrpy",1);  
   fhBookTree->SetBranchStatus("Npstripy",1);  
   fhBookTree->SetBranchStatus("Ntstripy",1);  
   fhBookTree->SetBranchStatus("Istripy",1);  
   fhBookTree->SetBranchStatus("Qstripy",1);  
   fhBookTree->SetBranchStatus("Ystripy",1);  
   
   
   
   Float_t ADCfull;  
   Int_t iladd=0;  
   for (Int_t ix=0; ix<Nstrpx;ix++) {  
   Iview=Npstripx[ix]*2-1;  
   Nstrip=(Int_t)Istripx[ix]-1;  
   if(Nstrip<fNstrips_ladder) iladd=0;  
   if((Nstrip>=fNstrips_ladder)&&(Nstrip<2*fNstrips_ladder)) iladd=1;  
   if((Nstrip>=2*fNstrips_ladder)&&(Nstrip<3*fNstrips_ladder)) iladd=2;  
   ADCfull=AdcTrack[Iview][Nstrip] += Qstripx[ix]*fMipCor[iladd][Iview];  
   AdcTrack[Iview][Nstrip] *= SaturationTrack(ADCfull);  
   
   };  
   
   
   for (Int_t iy=0; iy<Nstrpy;iy++) {  
   Iview=Npstripy[iy]*2-2;  
   Nstrip=(Int_t)Istripy[iy]-1;  
   if(Nstrip<fNstrips_ladder) iladd=0;  
   if((Nstrip>=fNstrips_ladder)&&(Nstrip<2*fNstrips_ladder)) iladd=1;  
   if((Nstrip>=2*fNstrips_ladder)&&(Nstrip<3*fNstrips_ladder)) iladd=2;  
   ADCfull=AdcTrack[Iview][Nstrip] -= Qstripy[iy]*fMipCor[iladd][Iview];  
   AdcTrack[Iview][Nstrip] *= SaturationTrack(ADCfull);  
   
   };    
         
 CompressTrackData(AdcTrack);  // Compress and Digitize data of one Ladder  in turn for all ladders  
   
 };  
   
   
   
 void Digitizer::DigitizeTrackCalib(Int_t ii) {  
   
 std:: cout << "Entering DigitizeTrackCalib " << ii << endl;  
 if( (ii!=1)&&(ii!=2) ) {  
  std:: cout << "error wrong DigitizeTrackCalib argument" << endl;  
  return;  
 };  
   
 memset(fDataTrack,0,sizeof(UShort_t)*fTRACKbuffer);  
 fTracklength=0;  
   
 UShort_t Dato;  
   
 Float_t dato1;  
 Float_t dato2;  
 Float_t dato3;  
 Float_t dato4;  
   
 UShort_t DatoDec;  
 UShort_t DatoDec1;  
 UShort_t DatoDec2;  
 UShort_t DatoDec3;  
 UShort_t DatoDec4;  
   
 UShort_t EVENT_CAL;  
 UShort_t PED_L1;  
 UShort_t ReLength;  
 UShort_t OveCheckCode;  
 //UShort_t PED_L2;  
 //UShort_t PED_L3HI;  
 //UShort_t PED_L3LO;  
 //UShort_t SIG_L1HI;  
 //UShort_t SIG_L1LO;  
 //UShort_t SIG_L2HI;  
 //UShort_t SIG_L2LO;  
 //UShort_t SIG_L3;  
 //UShort_t BAD_L1;  
 //UShort_t BAD_L2LO;  
 //UShort_t BAD_L3HI;  
 //UShort_t BAD_L3LO;  
 //UShort_t FLAG;  
   
   
   Int_t DSPpos;  
   for (Int_t j=ii-1; j<fNviews;j+=2) {  
   UShort_t CkSum=0;  
   // here skip the dsp header and his trailer , to be written later  
   DSPpos=fTracklength;  
   fTracklength=fTracklength+13+3;  
   
   
     for (Int_t i=0; i<fNladder;i++) {  
       for (Int_t k=0; k<fNstrips_ladder;k++) {  
       // write in buffer the current LADDER  
       Dato=(UShort_t)fPedeTrack[j][i*fNstrips_ladder+k];  
       dato1=fPedeTrack[j][i*fNstrips_ladder+k]-Dato;  
   
       DatoDec1=(UShort_t)(dato1*2);  
       dato2=dato1*2-DatoDec1;  
   
       DatoDec2=(UShort_t)(dato2*2);  
       dato3=dato2*2-DatoDec2;  
         
       DatoDec3=(UShort_t)(dato3*2);  
       dato4=dato3*2-DatoDec3;  
         
       DatoDec4=(UShort_t)(dato4*2);  
         
       DatoDec=DatoDec1*0x0008+DatoDec2*0x0004+DatoDec3*0x0002+DatoDec4*0x0001;  
       fDataTrack[fTracklength]=( (Dato << 4) | (DatoDec & 0x000F) );  
       CkSum=CkSum^fDataTrack[fTracklength];  
       fTracklength++;  
       };  
   
       for (Int_t k=0; k<fNstrips_ladder;k++) {  
       // write in buffer the current LADDER  
       Dato=(UShort_t)fSigmaTrack[j][i*fNstrips_ladder+k];  
       dato1=fSigmaTrack[j][i*fNstrips_ladder+k]-Dato;  
   
       DatoDec1=(UShort_t)(dato1*2);  
       dato2=dato1*2-DatoDec1;  
   
       DatoDec2=(UShort_t)(dato2*2);  
       dato3=dato2*2-DatoDec2;  
         
       DatoDec3=(UShort_t)(dato3*2);  
       dato4=dato3*2-DatoDec3;  
         
       DatoDec4=(UShort_t)(dato4*2);  
         
       DatoDec=DatoDec1*0x0008+DatoDec2*0x0004+DatoDec3*0x0002+DatoDec4*0x0001;  
         
       fDataTrack[fTracklength]=( (Dato << 4) | (DatoDec & 0x000F) );  
       CkSum=CkSum^fDataTrack[fTracklength];  
       fTracklength++;  
       };  
         
       for (Int_t k=0; k<64;k++) {  
       fDataTrack[fTracklength]=0x0000;  
       CkSum=CkSum^fDataTrack[fTracklength];  
       fTracklength++;  
   
       };  
       // end ladder  
   
     // write in buffer the end ladder word  
     if(i==0) fDataTrack[fTracklength]=0x1807;  
     if(i==1) fDataTrack[fTracklength]=0x1808;  
     if(i==2) fDataTrack[fTracklength]=0x1809;  
     CkSum=CkSum^fDataTrack[fTracklength];  
     fTracklength++;  
   
     // write in buffer the TRAILER  
     ReLength=(UShort_t)((fNstrips_ladder*2+64+1)*2+3);  
     OveCheckCode=0x0000;  
   
     fDataTrack[fTracklength]=0x0000;  
     fTracklength++;  
     
     fDataTrack[fTracklength]=(ReLength >> 8);  
     fTracklength++;  
     
     fDataTrack[fTracklength]=( (ReLength << 8) | (OveCheckCode & 0x00FF) );  
     fTracklength++;  
   
     // end TRAILER          
     };  
   
   // write in buffer the DSP header  
   
   fDataTrack[DSPpos]=(0xE800 | ( ((j+1) << 3) | 0x0005) );  
     
   fDataTrack[DSPpos+1]=0x01A9;  
   
   fDataTrack[DSPpos+2]=0x8740;  
     
   EVENT_CAL=0;  
   fDataTrack[DSPpos+3]=(0x1A00 | ( (0x03FF & EVENT_CAL)>> 1) );  
     
   PED_L1=0;  
   fDataTrack[DSPpos+4]=( ((EVENT_CAL << 15) | 0x5002 ) | ((0x03FF & PED_L1) << 2) );  
     
   // FROM HERE WE WRITE AS ALL VARIABLE apart CkSum are =0    
   
   fDataTrack[DSPpos+5]=0x8014;  
     
   fDataTrack[DSPpos+6]=0x00A0;  
     
   fDataTrack[DSPpos+7]=0x0500;  
     
   fDataTrack[DSPpos+8]=0x2801;  
     
   fDataTrack[DSPpos+9]=0x400A;  
     
   fDataTrack[DSPpos+10]=0x0050;  
     
   CkSum=(CkSum >> 8)^(CkSum&0x00FF);  
   fDataTrack[DSPpos+11]=(0x0280 | (CkSum >> 3));  
     
   fDataTrack[DSPpos+12]=(0x1FFF | (CkSum << 13) );  
   
   // end dsp header  
   
   // write in buffer the TRAILER  
     
   ReLength=(UShort_t)((13*2)+3);  
   OveCheckCode=0x0000;  
   fDataTrack[DSPpos+13]=0x0000;  
   
   fDataTrack[DSPpos+14]=(ReLength >> 8);  
     
   fDataTrack[DSPpos+15]=( (ReLength << 8) | (OveCheckCode & 0x00FF) );  
     
   // end TRAILER  
   
   
   
     
   // end DSP      
   };      
   
   
   
 };  
   
 void Digitizer::WriteTrackCalib() {  
   
   
 std:: cout << " Entering WriteTrackCalib " << endl;  
   
 fOutputfile.write(reinterpret_cast<char*>(fDataPSCU),sizeof(UShort_t)*fPSCUbuffer);  
   
 UShort_t temp[1000000];  
 memset(temp,0,sizeof(UShort_t)*1000000);  
 swab(fDataTrack,temp,sizeof(UShort_t)*fTracklength);  // WE MUST SWAP THE BYTES!!!  
 fOutputfile.write(reinterpret_cast<char*>(temp),sizeof(UShort_t)*fTracklength);  
 fTracklength=0;  
 if ( fPadding ){  
       fOutputfile.write(reinterpret_cast<char*>(fDataPadding),sizeof(UChar_t)*fPadding);  
 };  
   
 };  
   
   
 void Digitizer::ClearTrackCalib() {  
   
 std:: cout << "Entering ClearTrackCalib " << endl;  
       
     
 };  
   
   
 void Digitizer::LoadTrackCalib() {  
 std:: cout << "Entering LoadTrackCalib " << endl;  
   
 // Generate the pedestals and sigmas according to parametrization  
   for (Int_t j=0; j<fNviews;j++) {  
     for (Int_t i=0; i<fNstrips_view;i++) {  
       
     if((j+1)%2==0) {  
     fPedeTrack[j][i]=gRandom->Gaus(fAvePedex,fSigmaPedex);  
     fSigmaTrack[j][i]=gRandom->Gaus(fAveSigmax,fSigmaSigmax);  
     };  
     if((j+1)%2==1) {  
     fPedeTrack[j][i]=gRandom->Gaus(fAvePedey,fSigmaPedey);  
     fSigmaTrack[j][i]=gRandom->Gaus(fAveSigmay,fSigmaSigmay);  
     };  
       
     };  
   };      
   
     
     
 };  
   
 void Digitizer::LoadMipCor() {  
 std:: cout << "Entering LoadMipCor" << endl;  
   Float_t xfactor=1./151.6*1.04;  
   Float_t yfactor=1./152.1;  
   
   fMipCor[0][0]=140.02*yfactor;  
   fMipCor[0][1]=140.99*xfactor;  
   fMipCor[0][2]=134.48*yfactor;  
   fMipCor[0][3]=144.41*xfactor;  
   fMipCor[0][4]=140.74*yfactor;  
   fMipCor[0][5]=142.28*xfactor;  
   fMipCor[0][6]=134.53*yfactor;  
   fMipCor[0][7]=140.63*xfactor;  
   fMipCor[0][8]=135.55*yfactor;  
   fMipCor[0][9]=138.00*xfactor;  
   fMipCor[0][10]=154.95*yfactor;  
   fMipCor[0][11]=158.44*xfactor;  
     
     
   fMipCor[1][0]=136.07*yfactor;  
   fMipCor[1][1]=135.59*xfactor;  
   fMipCor[1][2]=142.69*yfactor;  
   fMipCor[1][3]=138.19*xfactor;  
   fMipCor[1][4]=137.35*yfactor;  
   fMipCor[1][5]=140.23*xfactor;  
   fMipCor[1][6]=153.15*yfactor;  
   fMipCor[1][7]=151.42*xfactor;  
   fMipCor[1][8]=129.76*yfactor;  
   fMipCor[1][9]=140.63*xfactor;  
   fMipCor[1][10]=157.87*yfactor;  
   fMipCor[1][11]=153.64*xfactor;  
   
   fMipCor[2][0]=134.98*yfactor;  
   fMipCor[2][1]=143.95*xfactor;  
   fMipCor[2][2]=140.23*yfactor;  
   fMipCor[2][3]=138.88*xfactor;  
   fMipCor[2][4]=137.95*yfactor;  
   fMipCor[2][5]=134.87*xfactor;  
   fMipCor[2][6]=157.56*yfactor;  
   fMipCor[2][7]=157.31*xfactor;  
   fMipCor[2][8]=141.37*yfactor;  
   fMipCor[2][9]=143.39*xfactor;  
   fMipCor[2][10]=156.15*yfactor;  
   fMipCor[2][11]=158.79*xfactor;  
   
 /*  
   for (Int_t j=0; j<fNviews;j++) {  
     for (Int_t i=0; i<fNstrips_view;i++) {  
     fMipCor[j][i]=1.;  
     };  
   };      
   
       
 */    
 };  
   
 void Digitizer::CompressTrackData(Float_t AdcTrack[fNviews][fNstrips_view]) {  
 // copy of the corresponding compression fortran routine + new digitization  
 // std:: cout << "Entering CompressTrackData " << endl;  
 Int_t oldval=0;  
 Int_t newval=0;  
 Int_t trasmesso=0;  
 Int_t ntrastot=0;  
 Float_t real;  
 Float_t inte;  
 Int_t cercacluster=0;  
 Int_t kt=0;  
 static const int DSPbufferSize = 4000; // 13 bit buffer to be rearranged in 16 bit Track buffer  
 UShort_t DataDSP[DSPbufferSize];  // 13 bit  buffer to be rearranged in 16 bit Track buffer  
 UShort_t DSPlength;  // 13 bit buffer to be rearranged in 16 bit Track buffer  
   
 memset(fDataTrack,0,sizeof(UShort_t)*fTRACKbuffer); // probably not necessary becouse already done ?  
 fTracklength=0;  
   
   for (Int_t iv=0; iv<fNviews;iv++) {  
   memset(DataDSP,0,sizeof(UShort_t)*DSPbufferSize);  
   DSPlength=16; // skip the header, to be written later  
   UShort_t CheckSum=0;  
 // write dsp header on buffer  
   
 //    fDataTrack[fTracklength]=0xE805;  
 //    fTracklength++;  
   
 //    fDataTrack[fTracklength]=0x01A9;  
 //    fTracklength++;  
   
 // end dsp header  
   
    //  
    // INIZIO VISTA IV - TAKE PROPER ACTION  
    //  
   
   
   
     for (Int_t ladder=0; ladder<fNladder;ladder++) {  
       Int_t k=0;  
       while (k<fNstrips_ladder) {  
         // compress write in buffer the current LADDER  
         if ( k == 0)  {  
           real=modff(AdcTrack[iv][ladder*fNstrips_ladder+k],&inte);  
           if (real > 0.5) inte=inte+1;  
           newval=(Int_t)inte -(Int_t)fPedeTrack[iv][ladder*fNstrips_ladder+k];  
           // first strip of ladder is transmitted  
           // DC_TOT first " << AdcTrack[iv][ladder*fNstrips_ladder+k] << endl;  
           DataDSP[DSPlength]=( ((UShort_t)inte) & 0x0FFF);  
           DSPlength++;    
           ntrastot++;  
           trasmesso=1;  
           oldval=newval;  
           kt=k;  
           k++;  
           continue;  
         };  
         real=modff(AdcTrack[iv][ladder*fNstrips_ladder+k],&inte);  
         if (real > 0.5) inte=inte+1;  
         newval=(Int_t)inte -(Int_t)(fPedeTrack[iv][ladder*fNstrips_ladder+k]);  
         cercacluster=1; // ?????????  
         if (cercacluster==1) {  
           
  // controlla l'ordine di tutti queste strip ladder e DSP !!!!!!!        
           Int_t diff=0;  
   
             
           switch ((iv+1)%2) {  
           case 0: diff=newval-oldval;  
           break;  
           case 1: diff=oldval-newval;      
           break;  
           };  
   
           if (diff>fCutclu*(Int_t)fSigmaTrack[iv][ladder*fNstrips_ladder+k]) {  
             Int_t clval=newval;  
             Int_t klp=k;        // go on to search for maximum  
             klp++;  
   
             while(klp<fNstrips_ladder) {  
               real=modff(AdcTrack[iv][ladder*fNstrips_ladder+klp],&inte);  
               if (real > 0.5) inte=inte+1;  
               Int_t clvalp=(Int_t)inte -(Int_t)fPedeTrack[iv][ladder*fNstrips_ladder+klp];  
               if((iv+1)%2==0) {  
   
                 if(clvalp>clval) {  
                    clval=clvalp;  
                    k=klp;}  
                 else break; // max of cluster found  
   
               } else {    
   
                 if(clvalp<clval) {  
                    clval=clvalp;  
                    k=klp;}  
                 else break; // max of cluster found  
   
               };  
                 
               klp++;  
             };  
   
             Int_t kl1=k-fNclst; // max of cluster (or end of ladder ?)  
             trasmesso=0;  
             if(kl1<0)  kl1=0;  
   
             if(kt>=kl1) kl1=kt+1;  
             if( (kt+1)==kl1 ) trasmesso=1;  
   
   
               
             Int_t kl2=k+fNclst;  
             if(kl2>=fNstrips_ladder) kl2=fNstrips_ladder-1;  
   
             for(Int_t klt=kl1 ; klt<=kl2 ; klt++) {  
               if(trasmesso==0) {  
               //  std:: cout << "STRIP " << klt << endl;  
               //  std:: cout << "ADC_TOT " <<AdcTrack[iv][ladder*fNstrips_ladder+klt] << endl;  
   
                 DataDSP[DSPlength]=( ((UShort_t)klt) | 0x1000);  
                 DSPlength++;    
                 ntrastot++;  
               
   
                 real=modff(AdcTrack[iv][ladder*fNstrips_ladder+klt],&inte);  
                 if (real > 0.5) inte=inte+1;  
                 DataDSP[DSPlength]=( ((UShort_t)inte) & 0x0FFF);  
                 DSPlength++;  
                 ntrastot++;  
               
                }  
                else {  
                //   std:: cout << "ADC_TOT " <<AdcTrack[iv][ladder*fNstrips_ladder+klt] << endl;  
                 real=modff(AdcTrack[iv][ladder*fNstrips_ladder+klt],&inte);  
                 if (real > 0.5) inte=inte+1;  
                 DataDSP[DSPlength]=( ((UShort_t)inte) & 0x0FFF);  
                 DSPlength++;                
                 ntrastot++;  
                 };  
                 trasmesso=1;                          
             };  // end trasmission  
             kt=kl2;  
             k=kl2;  
             real=modff(AdcTrack[iv][ladder*fNstrips_ladder+kt],&inte);  
             if (real > 0.5) inte=inte+1;  
             oldval=(Int_t)inte -(Int_t)fPedeTrack[iv][ladder*fNstrips_ladder+kt];  
             k++;  
             continue;  
   
   
           }; // end cercacluster  
         }; // end cercacluster  
         
 // start ZOP check for strips no  
         
       if(abs(newval-oldval)>=fCutzop*(Int_t)fSigmaTrack[iv][ladder*fNstrips_ladder+k]) {  
   
        if(trasmesso==0) {  
        // std:: cout << "STRIP " << k << endl;  
        // std:: cout << "ADC_TOT " << AdcTrack[iv][ladder*fNstrips_ladder+k] << endl;  
   
          DataDSP[DSPlength]=( ((UShort_t)k) | 0x1000);  
          DSPlength++;    
          ntrastot++;  
               
   
          real=modff(AdcTrack[iv][ladder*fNstrips_ladder+k],&inte);  
          if (real > 0.5) inte=inte+1;  
          DataDSP[DSPlength]=( ((UShort_t)inte) & 0x0FFF);  
          DSPlength++;  
          ntrastot++;  
   
        }  
        else {  
        //  std:: cout << "ADC_TOT " << AdcTrack[iv][ladder*fNstrips_ladder+k] << endl;  
          real=modff(AdcTrack[iv][ladder*fNstrips_ladder+k],&inte);  
          if (real > 0.5) inte=inte+1;  
          DataDSP[DSPlength]=(  ((UShort_t)inte) & 0x0FFF);  
          DSPlength++;  
          ntrastot++;        
        };  
        trasmesso=1;  
        oldval=newval;  
        kt=k;  
   
       }    
       else trasmesso=0;  
       // end zop  
         
       k++;        
       };  // end cycle inside ladder  
 // write here the end ladder bytes  
 //            std:: cout << "FINE LADDER " << ladder+1 << endl;  
   
       DataDSP[DSPlength]=( ((UShort_t)(ladder+1))  | 0x1800);  
       DSPlength++;  
       ntrastot++;  
       trasmesso=0;  
   
     };  //end cycle inside dsp  
 //  std:: cout << "FINE DSP " << iv+1 << endl;  
 // here put DSP header  
     DataDSP[0]=(0x1CA0 | ((UShort_t)(iv+1)) );  
     UShort_t Nword=(DSPlength*13)/16;  
     if( ((DSPlength*13)%16)!=0) Nword++;  
     DataDSP[1]=(0x1400 | ( Nword >> 10));  
     DataDSP[2]=(0x1400 | ( Nword & 0x03FF) );  
     DataDSP[3]=(0x1400 | (( (UShort_t)(fCounter >> 10) ) & 0x03FF) );  
     DataDSP[4]=(0x1400 | (( (UShort_t)(fCounter) )  & 0x03FF) );  
     DataDSP[5]=(0x1400 | ( (UShort_t)(fNclst << 7) ) | ( (UShort_t)(fCutzop << 4) )  
      | ( (UShort_t)fCutzop  ) );  
     DataDSP[6]=0x1400;  
     DataDSP[7]=0x1400;  
     DataDSP[8]=0x1400;  
     DataDSP[9]=0x1400;  
     DataDSP[10]=0x1400;  
     DataDSP[11]=0x1400;  
     DataDSP[12]=0x1400;  
     DataDSP[13]=0x1400;  
     DataDSP[14]=(0x1400 | (CheckSum & 0x00FF) );  
     DataDSP[15]=0x1C00;  
 // end DSP header      
   
   
 // write 13 bit DataDSP bufer inside 16 bit fDataTrack buffer  
     Int_t Bit16free=16;  
     UShort_t Dato;  
     for (Int_t NDSP=0; NDSP<DSPlength;NDSP++) {  
       Int_t Bit13ToWrite=13;  
       while(Bit13ToWrite>0) {  
         if(Bit13ToWrite<=Bit16free) {  
           Dato=((DataDSP[NDSP]&(0xFFFF >> (16-Bit13ToWrite)))<<(Bit16free-Bit13ToWrite));  
           fDataTrack[fTracklength]=fDataTrack[fTracklength] | Dato ;  
           Bit16free=Bit16free-Bit13ToWrite;  
           Bit13ToWrite=0;  
           if(Bit16free==0) {  
             if(NDSP>15) CheckSum=CheckSum^fDataTrack[fTracklength];  
             fTracklength++;  
             Bit16free=16;  
           };            
         }  
         else if(Bit13ToWrite>Bit16free) {  
           Dato=( (DataDSP[NDSP]&(0xFFFF >> (16-Bit13ToWrite) ) ) >> (Bit13ToWrite-Bit16free) );  
           fDataTrack[fTracklength]=fDataTrack[fTracklength] | Dato ;  
           if(NDSP>15) CheckSum=CheckSum^fDataTrack[fTracklength];  
           fTracklength++;  
           Bit13ToWrite=Bit13ToWrite-Bit16free;  
           Bit16free=16;    
         };  
           
       }; // end cycle while(Bit13ToWrite>0)  
             
     }; // end cycle DataDSP  
     if(Bit16free!=16) { fTracklength++; CheckSum=CheckSum^fDataTrack[fTracklength]; };  
     CheckSum=(CheckSum >> 8)^(CheckSum&0x00FF);  
     fDataTrack[fTracklength-Nword+11]=(0x0280 | (CheckSum >> 3));  
     fDataTrack[fTracklength-Nword+12]=(0x1C00 | (CheckSum << 13) );  
   
 // end write 13 bit DataDSP bufer inside 16 bit fDataTrack buffer  
   
 //write trailer on buffer  
     UShort_t ReLength=(UShort_t)((Nword+13)*2+3);  
     UShort_t OveCheckCode=0x0000;  
   
     fDataTrack[fTracklength]=0x0000;  
     fTracklength++;  
     
     fDataTrack[fTracklength]=(ReLength >> 8);  
     fTracklength++;  
     
     fDataTrack[fTracklength]=( (ReLength << 8) | (OveCheckCode & 0x00FF) );  
     fTracklength++;    
 // end trailer  
 //    std:: cout  << "DSPlength  " <<DSPlength << endl;  
 //    std:: cout << "Nword " << Nword  << endl;  
 //    std:: cout <<  "ReLength " << ReLength << endl;  
   };      
 //    std:: cout << "ntrastot " << ntrastot << endl;      
   
 };  
   
   
 Float_t Digitizer::SaturationTrack(Float_t ADC) {  
   Float_t SatFact=1.;  
   if(ADC<70.) { SatFact=80./ADC; };  
   if(ADC>3000.) { SatFact=3000./ADC; };  
   return SatFact;  
 };  
   
   
   
   
   
   
Legend:
Removed from v.1.5  
changed lines
  Added in v.1.14
  ViewVC Help
Powered by ViewVC 1.1.23