/[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.2 by orsi, Fri Sep 28 10:46:23 2007 UTC revision 1.5 by silvio, Wed Nov 28 18:54:31 2007 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 <sstream>
8  #include <fstream>  #include <fstream>
9  #include <stdlib.h>  #include <stdlib.h>
10  #include <stdio.h>  #include <stdio.h>
11  #include <string.h>  #include <string.h>
12  #include <ctype.h>  #include <ctype.h>
13  #include <time.h>  #include <time.h>
14  #include "Riostream.h"  #include "Riostream.h"
15  #include "TFile.h"  #include "TFile.h"
16  #include "TDirectory.h"  #include "TDirectory.h"
17  #include "TTree.h"  #include "TTree.h"
18  #include "TLeafI.h"  #include "TLeafI.h"
19  #include "TH1.h"  #include "TH1.h"
20  #include "TH2.h"  #include "TH2.h"
21  #include "TMath.h"  #include "TMath.h"
22  #include "TRandom.h"  #include "TRandom.h"
23  #include "TSQLServer.h"  #include "TSQLServer.h"
24  #include "TSystem.h"  #include "TSystem.h"
25  //  //
26  #include "Digitizer.h"  #include "Digitizer.h"
27  #include "CRC.h"  #include "CRC.h"
28  //  //
29  #include <PamelaRun.h>  #include <PamelaRun.h>
30  #include <physics/calorimeter/CalorimeterEvent.h>  #include <physics/calorimeter/CalorimeterEvent.h>
31  #include <CalibCalPedEvent.h>  #include <CalibCalPedEvent.h>
32  #include "GLTables.h"  #include "GLTables.h"
33  //  //
34  extern "C"{  extern "C"{
35    short crc(short, short);    short crc(short, short);
36  };  };
37  //  //
38    
39  Digitizer::Digitizer(TTree* tree, char* &file_raw){  Digitizer::Digitizer(TTree* tree, char* &file_raw){
40    fhBookTree = tree;    fhBookTree = tree;
41    fFilename =  file_raw;    fFilename =  file_raw;
42    fCounter = 0;    fCounter = 0;
43    fOBT = 0;    fOBT = 0;
44    
45    //    //
46    // DB connections    // DB connections
47    //    //
48    TString host = "mysql://localhost/pamelaprod";    TString host = "mysql://localhost/pamelaprod";
49    TString user = "anonymous";    TString user = "anonymous";
50    TString psw = "";    TString psw = "";
51    //    //
52    const char *pamdbhost=gSystem->Getenv("PAM_DBHOST");    const char *pamdbhost=gSystem->Getenv("PAM_DBHOST");
53    const char *pamdbuser=gSystem->Getenv("PAM_DBUSER");    const char *pamdbuser=gSystem->Getenv("PAM_DBUSER");
54    const char *pamdbpsw=gSystem->Getenv("PAM_DBPSW");    const char *pamdbpsw=gSystem->Getenv("PAM_DBPSW");
55    if ( !pamdbhost ) pamdbhost = "";    if ( !pamdbhost ) pamdbhost = "";
56    if ( !pamdbuser ) pamdbuser = "";    if ( !pamdbuser ) pamdbuser = "";
57    if ( !pamdbpsw ) pamdbpsw = "";    if ( !pamdbpsw ) pamdbpsw = "";
58    if ( strcmp(pamdbhost,"") ) host = pamdbhost;    if ( strcmp(pamdbhost,"") ) host = pamdbhost;
59    if ( strcmp(pamdbuser,"") ) user = pamdbuser;    if ( strcmp(pamdbuser,"") ) user = pamdbuser;
60    if ( strcmp(pamdbpsw,"") ) psw = pamdbpsw;    if ( strcmp(pamdbpsw,"") ) psw = pamdbpsw;
61    fDbc = TSQLServer::Connect(host.Data(),user.Data(),psw.Data());    fDbc = TSQLServer::Connect(host.Data(),user.Data(),psw.Data());
62    //    //
63    GL_TABLES *glt = new GL_TABLES(host,user,psw);    GL_TABLES *glt = new GL_TABLES(host,user,psw);
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());    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    //    //
66    // Use UTC in the DB and make timeout bigger    // Use UTC in the DB and make timeout bigger
67    //    //
68    stringstream myquery;    stringstream myquery;
69    myquery.str("");    myquery.str("");
70    myquery << "SET time_zone='+0:00'";    myquery << "SET time_zone='+0:00'";
71    fDbc->Query(myquery.str().c_str());    fDbc->Query(myquery.str().c_str());
72    myquery.str("");    myquery.str("");
73    myquery << "SET wait_timeout=173000;";    myquery << "SET wait_timeout=173000;";
74    fDbc->Query(myquery.str().c_str());    fDbc->Query(myquery.str().c_str());
75    //    //
76        
77    std:: cout << "preparing tree" << endl;    std:: cout << "preparing tree" << endl;
78    
79    // prepare tree    // prepare tree
80    fhBookTree->SetBranchAddress("Irun",&Irun);    fhBookTree->SetBranchAddress("Irun",&Irun);
81    fhBookTree->SetBranchAddress("Ievnt",&Ievnt);    fhBookTree->SetBranchAddress("Ievnt",&Ievnt);
82    fhBookTree->SetBranchAddress("Ipa",&Ipa);    fhBookTree->SetBranchAddress("Ipa",&Ipa);
83    fhBookTree->SetBranchAddress("X0",&X0);    fhBookTree->SetBranchAddress("X0",&X0);
84    fhBookTree->SetBranchAddress("Y0",&Y0);    fhBookTree->SetBranchAddress("Y0",&Y0);
85    fhBookTree->SetBranchAddress("Z0",&Z0);    fhBookTree->SetBranchAddress("Z0",&Z0);
86    fhBookTree->SetBranchAddress("Theta",&Theta);    fhBookTree->SetBranchAddress("Theta",&Theta);
87    fhBookTree->SetBranchAddress("Phi",&Phi);    fhBookTree->SetBranchAddress("Phi",&Phi);
88    fhBookTree->SetBranchAddress("P0",&P0);    fhBookTree->SetBranchAddress("P0",&P0);
89    fhBookTree->SetBranchAddress("Nthtof",&Nthtof);    fhBookTree->SetBranchAddress("Nthtof",&Nthtof);
90    fhBookTree->SetBranchAddress("Ipltof",Ipltof);    fhBookTree->SetBranchAddress("Ipltof",Ipltof);
91    fhBookTree->SetBranchAddress("Ipaddle",Ipaddle);    fhBookTree->SetBranchAddress("Ipaddle",Ipaddle);
92    fhBookTree->SetBranchAddress("Ipartof",Ipartof);    fhBookTree->SetBranchAddress("Ipartof",Ipartof);
93    fhBookTree->SetBranchAddress("Xintof",Xintof);    fhBookTree->SetBranchAddress("Xintof",Xintof);
94    fhBookTree->SetBranchAddress("Yintof",Yintof);    fhBookTree->SetBranchAddress("Yintof",Yintof);
95    fhBookTree->SetBranchAddress("Zintof",Zintof);    fhBookTree->SetBranchAddress("Zintof",Zintof);
96    fhBookTree->SetBranchAddress("Xouttof",Xouttof);    fhBookTree->SetBranchAddress("Xouttof",Xouttof);
97    fhBookTree->SetBranchAddress("Youttof",Youttof);    fhBookTree->SetBranchAddress("Youttof",Youttof);
98    fhBookTree->SetBranchAddress("Zouttof",Zouttof);    fhBookTree->SetBranchAddress("Zouttof",Zouttof);
99    fhBookTree->SetBranchAddress("Ereltof",Ereltof);    fhBookTree->SetBranchAddress("Ereltof",Ereltof);
100    fhBookTree->SetBranchAddress("Timetof",Timetof);    fhBookTree->SetBranchAddress("Timetof",Timetof);
101    fhBookTree->SetBranchAddress("Pathtof",Pathtof);    fhBookTree->SetBranchAddress("Pathtof",Pathtof);
102    fhBookTree->SetBranchAddress("P0tof",P0tof);    fhBookTree->SetBranchAddress("P0tof",P0tof);
103    fhBookTree->SetBranchAddress("Nthcat",&Nthcat);    fhBookTree->SetBranchAddress("Nthcat",&Nthcat);
104    fhBookTree->SetBranchAddress("Iparcat",Iparcat);    fhBookTree->SetBranchAddress("Iparcat",Iparcat);
105    fhBookTree->SetBranchAddress("Icat",Icat);    fhBookTree->SetBranchAddress("Icat",Icat);
106    fhBookTree->SetBranchAddress("Xincat",Xincat);    fhBookTree->SetBranchAddress("Xincat",Xincat);
107    fhBookTree->SetBranchAddress("Yincat",Yincat);    fhBookTree->SetBranchAddress("Yincat",Yincat);
108    fhBookTree->SetBranchAddress("Zincat",Zincat);    fhBookTree->SetBranchAddress("Zincat",Zincat);
109    fhBookTree->SetBranchAddress("Xoutcat",Xoutcat);    fhBookTree->SetBranchAddress("Xoutcat",Xoutcat);
110    fhBookTree->SetBranchAddress("Youtcat",Youtcat);    fhBookTree->SetBranchAddress("Youtcat",Youtcat);
111    fhBookTree->SetBranchAddress("Zoutcat",Zoutcat);    fhBookTree->SetBranchAddress("Zoutcat",Zoutcat);
112    fhBookTree->SetBranchAddress("Erelcat",Erelcat);    fhBookTree->SetBranchAddress("Erelcat",Erelcat);
113    fhBookTree->SetBranchAddress("Timecat",Timecat);    fhBookTree->SetBranchAddress("Timecat",Timecat);
114    fhBookTree->SetBranchAddress("Pathcat",Pathcat);    fhBookTree->SetBranchAddress("Pathcat",Pathcat);
115    fhBookTree->SetBranchAddress("P0cat",P0cat);    fhBookTree->SetBranchAddress("P0cat",P0cat);
116    fhBookTree->SetBranchAddress("Nthcas",&Nthcas);    fhBookTree->SetBranchAddress("Nthcas",&Nthcas);
117    fhBookTree->SetBranchAddress("Iparcas",Iparcas);    fhBookTree->SetBranchAddress("Iparcas",Iparcas);
118    fhBookTree->SetBranchAddress("Icas",Icas);    fhBookTree->SetBranchAddress("Icas",Icas);
119    fhBookTree->SetBranchAddress("Xincas",Xincas);    fhBookTree->SetBranchAddress("Xincas",Xincas);
120    fhBookTree->SetBranchAddress("Yincas",Yincas);    fhBookTree->SetBranchAddress("Yincas",Yincas);
121    fhBookTree->SetBranchAddress("Zincas",Zincas);    fhBookTree->SetBranchAddress("Zincas",Zincas);
122    fhBookTree->SetBranchAddress("Xoutcas",Xoutcas);    fhBookTree->SetBranchAddress("Xoutcas",Xoutcas);
123    fhBookTree->SetBranchAddress("Youtcas",Youtcas);    fhBookTree->SetBranchAddress("Youtcas",Youtcas);
124    fhBookTree->SetBranchAddress("Zoutcas",Zoutcas);    fhBookTree->SetBranchAddress("Zoutcas",Zoutcas);
125    fhBookTree->SetBranchAddress("Erelcas",Erelcas);    fhBookTree->SetBranchAddress("Erelcas",Erelcas);
126    fhBookTree->SetBranchAddress("Timecas",Timecas);    fhBookTree->SetBranchAddress("Timecas",Timecas);
127    fhBookTree->SetBranchAddress("Pathcas",Pathcas);    fhBookTree->SetBranchAddress("Pathcas",Pathcas);
128    fhBookTree->SetBranchAddress("P0cas",P0cas);    fhBookTree->SetBranchAddress("P0cas",P0cas);
129    fhBookTree->SetBranchAddress("Nthspe",&Nthspe);    fhBookTree->SetBranchAddress("Nthspe",&Nthspe);
130    fhBookTree->SetBranchAddress("Iparspe",Iparspe);    fhBookTree->SetBranchAddress("Iparspe",Iparspe);
131    fhBookTree->SetBranchAddress("Itrpb",Itrpb);    fhBookTree->SetBranchAddress("Itrpb",Itrpb);
132    fhBookTree->SetBranchAddress("Itrsl",Itrsl);    fhBookTree->SetBranchAddress("Itrsl",Itrsl);
133    fhBookTree->SetBranchAddress("Itspa",Itspa);    fhBookTree->SetBranchAddress("Itspa",Itspa);
134    fhBookTree->SetBranchAddress("Xinspe",Xinspe);    fhBookTree->SetBranchAddress("Xinspe",Xinspe);
135    fhBookTree->SetBranchAddress("Yinspe",Yinspe);    fhBookTree->SetBranchAddress("Yinspe",Yinspe);
136    fhBookTree->SetBranchAddress("Zinspe",Zinspe);    fhBookTree->SetBranchAddress("Zinspe",Zinspe);
137    fhBookTree->SetBranchAddress("Xoutspe",Xoutspe);    fhBookTree->SetBranchAddress("Xoutspe",Xoutspe);
138    fhBookTree->SetBranchAddress("Youtspe",Youtspe);    fhBookTree->SetBranchAddress("Youtspe",Youtspe);
139    fhBookTree->SetBranchAddress("Zoutspe",Zoutspe);    fhBookTree->SetBranchAddress("Zoutspe",Zoutspe);
140    fhBookTree->SetBranchAddress("Xavspe",Xavspe);    fhBookTree->SetBranchAddress("Xavspe",Xavspe);
141    fhBookTree->SetBranchAddress("Yavspe",Yavspe);    fhBookTree->SetBranchAddress("Yavspe",Yavspe);
142    fhBookTree->SetBranchAddress("Zavspe",Zavspe);    fhBookTree->SetBranchAddress("Zavspe",Zavspe);
143    fhBookTree->SetBranchAddress("Erelspe",Erelspe);    fhBookTree->SetBranchAddress("Erelspe",Erelspe);
144    fhBookTree->SetBranchAddress("Pathspe",Pathspe);    fhBookTree->SetBranchAddress("Pathspe",Pathspe);
145    fhBookTree->SetBranchAddress("P0spe",P0spe);    fhBookTree->SetBranchAddress("P0spe",P0spe);
146    fhBookTree->SetBranchAddress("Nxmult",Nxmult);    fhBookTree->SetBranchAddress("Nxmult",Nxmult);
147    fhBookTree->SetBranchAddress("Nymult",Nymult);    fhBookTree->SetBranchAddress("Nymult",Nymult);
148    fhBookTree->SetBranchAddress("Nstrpx",&Nstrpx);    fhBookTree->SetBranchAddress("Nstrpx",&Nstrpx);
149    fhBookTree->SetBranchAddress("Npstripx",Npstripx);    fhBookTree->SetBranchAddress("Npstripx",Npstripx);
150    fhBookTree->SetBranchAddress("Ntstripx",Ntstripx);    fhBookTree->SetBranchAddress("Ntstripx",Ntstripx);
151    fhBookTree->SetBranchAddress("Istripx",Istripx);    fhBookTree->SetBranchAddress("Istripx",Istripx);
152    fhBookTree->SetBranchAddress("Qstripx",Qstripx);    fhBookTree->SetBranchAddress("Qstripx",Qstripx);
153    fhBookTree->SetBranchAddress("Xstripx",Xstripx);    fhBookTree->SetBranchAddress("Xstripx",Xstripx);
154    fhBookTree->SetBranchAddress("Nstrpy",&Nstrpy);    fhBookTree->SetBranchAddress("Nstrpy",&Nstrpy);
155    fhBookTree->SetBranchAddress("Npstripy",Npstripy);    fhBookTree->SetBranchAddress("Npstripy",Npstripy);
156    fhBookTree->SetBranchAddress("Ntstripy",Ntstripy);    fhBookTree->SetBranchAddress("Ntstripy",Ntstripy);
157    fhBookTree->SetBranchAddress("Istripy",Istripy);    fhBookTree->SetBranchAddress("Istripy",Istripy);
158    fhBookTree->SetBranchAddress("Qstripy",Qstripy);    fhBookTree->SetBranchAddress("Qstripy",Qstripy);
159    fhBookTree->SetBranchAddress("Ystripy",Ystripy);    fhBookTree->SetBranchAddress("Ystripy",Ystripy);
160    fhBookTree->SetBranchAddress("Nthcali",&Nthcali);    fhBookTree->SetBranchAddress("Nthcali",&Nthcali);
161    fhBookTree->SetBranchAddress("Icaplane",Icaplane);    fhBookTree->SetBranchAddress("Icaplane",Icaplane);
162    fhBookTree->SetBranchAddress("Icastrip",Icastrip);    fhBookTree->SetBranchAddress("Icastrip",Icastrip);
163    fhBookTree->SetBranchAddress("Icamod",Icamod);    fhBookTree->SetBranchAddress("Icamod",Icamod);
164    fhBookTree->SetBranchAddress("Enestrip",Enestrip);    fhBookTree->SetBranchAddress("Enestrip",Enestrip);
165    fhBookTree->SetBranchAddress("Nthcal",&Nthcal);    fhBookTree->SetBranchAddress("Nthcal",&Nthcal);
166    fhBookTree->SetBranchAddress("Icapl",Icapl);    fhBookTree->SetBranchAddress("Icapl",Icapl);
167    fhBookTree->SetBranchAddress("Icasi",Icasi);    fhBookTree->SetBranchAddress("Icasi",Icasi);
168    fhBookTree->SetBranchAddress("Icast",Icast);    fhBookTree->SetBranchAddress("Icast",Icast);
169    fhBookTree->SetBranchAddress("Xincal",Xincal);    fhBookTree->SetBranchAddress("Xincal",Xincal);
170    fhBookTree->SetBranchAddress("Yincal",Yincal);    fhBookTree->SetBranchAddress("Yincal",Yincal);
171    fhBookTree->SetBranchAddress("Zincal",Zincal);    fhBookTree->SetBranchAddress("Zincal",Zincal);
172    fhBookTree->SetBranchAddress("Erelcal",Erelcal);    fhBookTree->SetBranchAddress("Erelcal",Erelcal);
173    fhBookTree->SetBranchAddress("Nthnd",&Nthnd);    fhBookTree->SetBranchAddress("Nthnd",&Nthnd);
174    fhBookTree->SetBranchAddress("Itubend",Itubend);    fhBookTree->SetBranchAddress("Itubend",Itubend);
175    fhBookTree->SetBranchAddress("Iparnd",Iparnd);    fhBookTree->SetBranchAddress("Iparnd",Iparnd);
176    fhBookTree->SetBranchAddress("Xinnd",Xinnd);    fhBookTree->SetBranchAddress("Xinnd",Xinnd);
177    fhBookTree->SetBranchAddress("Yinnd",Yinnd);    fhBookTree->SetBranchAddress("Yinnd",Yinnd);
178    fhBookTree->SetBranchAddress("Zinnd",Zinnd);    fhBookTree->SetBranchAddress("Zinnd",Zinnd);
179    fhBookTree->SetBranchAddress("Xoutnd",Xoutnd);    fhBookTree->SetBranchAddress("Xoutnd",Xoutnd);
180    fhBookTree->SetBranchAddress("Youtnd",Youtnd);    fhBookTree->SetBranchAddress("Youtnd",Youtnd);
181    fhBookTree->SetBranchAddress("Zoutnd",Zoutnd);    fhBookTree->SetBranchAddress("Zoutnd",Zoutnd);
182    fhBookTree->SetBranchAddress("Erelnd",Erelnd);    fhBookTree->SetBranchAddress("Erelnd",Erelnd);
183    fhBookTree->SetBranchAddress("Timend",Timend);    fhBookTree->SetBranchAddress("Timend",Timend);
184    fhBookTree->SetBranchAddress("Pathnd",Pathnd);    fhBookTree->SetBranchAddress("Pathnd",Pathnd);
185    fhBookTree->SetBranchAddress("P0nd",P0nd);    fhBookTree->SetBranchAddress("P0nd",P0nd);
186    fhBookTree->SetBranchAddress("Nthcard",&Nthcard);    fhBookTree->SetBranchAddress("Nthcard",&Nthcard);
187    fhBookTree->SetBranchAddress("Iparcard",Iparcard);    fhBookTree->SetBranchAddress("Iparcard",Iparcard);
188    fhBookTree->SetBranchAddress("Icard",Icard);    fhBookTree->SetBranchAddress("Icard",Icard);
189    fhBookTree->SetBranchAddress("Xincard",Xincard);    fhBookTree->SetBranchAddress("Xincard",Xincard);
190    fhBookTree->SetBranchAddress("Yincard",Yincard);    fhBookTree->SetBranchAddress("Yincard",Yincard);
191    fhBookTree->SetBranchAddress("Zincard",Zincard);    fhBookTree->SetBranchAddress("Zincard",Zincard);
192    fhBookTree->SetBranchAddress("Xoutcard",Xoutcard);    fhBookTree->SetBranchAddress("Xoutcard",Xoutcard);
193    fhBookTree->SetBranchAddress("Youtcard",Youtcard);    fhBookTree->SetBranchAddress("Youtcard",Youtcard);
194    fhBookTree->SetBranchAddress("Zoutcard",Zoutcard);    fhBookTree->SetBranchAddress("Zoutcard",Zoutcard);
195    fhBookTree->SetBranchAddress("Erelcard",Erelcard);    fhBookTree->SetBranchAddress("Erelcard",Erelcard);
196    fhBookTree->SetBranchAddress("Timecard",Timecard);    fhBookTree->SetBranchAddress("Timecard",Timecard);
197    fhBookTree->SetBranchAddress("Pathcard",Pathcard);    fhBookTree->SetBranchAddress("Pathcard",Pathcard);
198    fhBookTree->SetBranchAddress("P0card",P0card);    fhBookTree->SetBranchAddress("P0card",P0card);
199    
200    fhBookTree->SetBranchStatus("*",0);    fhBookTree->SetBranchStatus("*",0);
201    
202  };  };
203    
204    
205    
206  void Digitizer::Close(){  void Digitizer::Close(){
207    
208    delete fhBookTree;    delete fhBookTree;
209    
210  };  };
211    
212    
213    
214    
215  void Digitizer::Loop() {  void Digitizer::Loop() {
216    //    //
217    // opens the raw output file and loops over the events    // opens the raw output file and loops over the events
218    //    //
219    fOutputfile.open(fFilename, ios::out | ios::binary);    fOutputfile.open(fFilename, ios::out | ios::binary);
220    //fOutputfile.open(Form("Output%s",fFilename), ios::out | ios::binary);    //fOutputfile.open(Form("Output%s",fFilename), ios::out | ios::binary);
221    //    //
222    // Load in memory and save at the beginning of file the calorimeter calibration    // Load in memory and save at the beginning of file the calorimeter calibration
223    //    //
224    CaloLoadCalib();    CaloLoadCalib();
225    DigitizeCALOCALIB();    DigitizeCALOCALIB();
226    
227    //  load, digitize and write tracker calibration    //  load, digitize and write tracker calibration
228    LoadTrackCalib();    LoadTrackCalib();
229        
230    DigitizeTrackCalib(1);    DigitizeTrackCalib(1);
231    UInt_t length=fTracklength*2;    UInt_t length=fTracklength*2;
232    DigitizePSCU(length,0x12);    DigitizePSCU(length,0x12);
233    AddPadding();    AddPadding();
234    WriteTrackCalib();    WriteTrackCalib();
235        
236    DigitizeTrackCalib(2);    DigitizeTrackCalib(2);
237    length=fTracklength*2;    length=fTracklength*2;
238    DigitizePSCU(length,0x13);    DigitizePSCU(length,0x13);
239    AddPadding();    AddPadding();
240    WriteTrackCalib();    WriteTrackCalib();
241        
242    LoadMipCor();  // some initialization of parameters -not used now-    LoadMipCor();  // some initialization of parameters -not used now-
243    //  end loading, digitizing and writing tracker calibration    //  end loading, digitizing and writing tracker calibration
244    
245    //    //
246    // loops over the events    // loops over the events
247    //    //
248        
249    Int_t nentries = fhBookTree->GetEntriesFast();    Int_t nentries = fhBookTree->GetEntriesFast();
250    Long64_t nbytes = 0;    Long64_t nbytes = 0;
251    for (Int_t i=0; i<nentries;i++) {    for (Int_t i=0; i<nentries;i++) {
252      //      //
253      nbytes += fhBookTree->GetEntry(i);      nbytes += fhBookTree->GetEntry(i);
254      // read detectors sequentially:      // read detectors sequentially:
255      // http://www.ts.infn.it/fileadmin/documents/physics/experiments/wizard/cpu/gen_arch/RM_Acquisition.pdf      // http://www.ts.infn.it/fileadmin/documents/physics/experiments/wizard/cpu/gen_arch/RM_Acquisition.pdf
256      // on pamelatov:      // on pamelatov:
257      // /cvs/yoda/techmodel/physics/NeutronDetectorReader.cpp      // /cvs/yoda/techmodel/physics/NeutronDetectorReader.cpp
258      DigitizeTRIGGER();      DigitizeTRIGGER();
259      DigitizeTOF();      DigitizeTOF();
260      DigitizeAC();      DigitizeAC();
261      DigitizeCALO();      DigitizeCALO();
262      DigitizeTrack();      DigitizeTrack();
263      //DigitizeS4();      DigitizeS4();
264      DigitizeND();      DigitizeND();
265        //        //
266      // Add padding to 64 bits      // Add padding to 64 bits
267      //      //
268      AddPadding();      AddPadding();
269  //  //
270      // Create CPU header, we need packet type (0x10 = physics data) and packet length.      // Create CPU header, we need packet type (0x10 = physics data) and packet length.
271      //      //
272      //UInt_t length=2*(fCALOlength+fACbuffer+fTracklength+fNDbuffer+fS4buffer)+fPadding+fTOFbuffer+fTRIGGERbuffer;      UInt_t length=2*(fCALOlength+fACbuffer+fTracklength+fNDbuffer+fS4buffer)+fPadding+fTOFbuffer+fTRIGGERbuffer;
273      UInt_t length=2*(fCALOlength+fACbuffer+fTracklength+fNDbuffer)+fPadding+fTOFbuffer+fTRIGGERbuffer;      //UInt_t length=2*(fCALOlength+fACbuffer+fTracklength+fNDbuffer)+fPadding+fTOFbuffer+fTRIGGERbuffer;
274      DigitizePSCU(length,0x10);      DigitizePSCU(length,0x10);
275      if ( !i%100 ) std::cout << "writing event " << i << endl;      if ( !i%100 ) std::cout << "writing event " << i << endl;
276      WriteData();      WriteData();
277    };    };
278    
279    fOutputfile.close();    fOutputfile.close();
280    std::cout << "files closed" << endl << flush;    std::cout << "files closed" << endl << flush;
281    
282  };  };
283    
284  void Digitizer::AddPadding(){  void Digitizer::AddPadding(){
285    //    //
286    Float_t pd0 = (fLen+16)/64.;    Float_t pd0 = (fLen+16)/64.;
287    Float_t pd1 = pd0 - (Float_t)int(pd0);    Float_t pd1 = pd0 - (Float_t)int(pd0);
288    Float_t padfrac = 64. - pd1 * 64.;    Float_t padfrac = 64. - pd1 * 64.;
289    //    //
290    UInt_t padbytes = (UInt_t)padfrac;    UInt_t padbytes = (UInt_t)padfrac;
291    if ( padbytes > 0 && padbytes < 64 ){    if ( padbytes > 0 && padbytes < 64 ){
292      //      //
293      // here the padding length      // here the padding length
294      //      //
295      fPadding = padbytes+64;      fPadding = padbytes+64;
296      //      //
297      // random padding bytes      // random padding bytes
298      //      //
299      for (Int_t ur=0; ur<32; ur++){      for (Int_t ur=0; ur<32; ur++){
300        fDataPadding[ur] = (UShort_t)rand();        fDataPadding[ur] = (UShort_t)rand();
301      };      };
302    };    };
303  };  };
304    
305    
306  void Digitizer::DigitizePSCU(UInt_t length, UChar_t type) {  void Digitizer::DigitizePSCU(UInt_t length, UChar_t type) {
307    //    //
308    UChar_t buff[16];    UChar_t buff[16];
309    //    //
310    // CPU signature    // CPU signature
311    //      //  
312    buff[0] = 0xFA;    buff[0] = 0xFA;
313    buff[1] = 0xFE;    buff[1] = 0xFE;
314    buff[2] = 0xDE;    buff[2] = 0xDE;
315    //    //
316    // packet type (twice)    // packet type (twice)
317    //    //
318    buff[3] = type;    buff[3] = type;
319    buff[4] = type;    buff[4] = type;
320    //    //
321    // counter    // counter
322    //    //
323    fCounter++;    fCounter++;
324    while ( fCounter > 16777215 ){    while ( fCounter > 16777215 ){
325      fCounter -= 16777215;      fCounter -= 16777215;
326    };    };
327    //    //
328    buff[5] = (UChar_t)(fCounter >> 16);    buff[5] = (UChar_t)(fCounter >> 16);
329    buff[6] = (UChar_t)(fCounter >> 8);    buff[6] = (UChar_t)(fCounter >> 8);
330    buff[7] = (UChar_t)fCounter;    buff[7] = (UChar_t)fCounter;
331    //    //
332    // on board time    // on board time
333    //    //
334    ULong64_t obt = fOBT + 30LL;    ULong64_t obt = fOBT + 30LL;
335    //    //
336    while ( obt > 4294967295LL ){    while ( obt > 4294967295LL ){
337      obt -= 4294967295LL;      obt -= 4294967295LL;
338    };    };
339    fOBT = (UInt_t)obt;    fOBT = (UInt_t)obt;
340    //    //
341    buff[8] = (UChar_t)(fOBT >> 24);    buff[8] = (UChar_t)(fOBT >> 24);
342    buff[9] = (UChar_t)(fOBT >> 16);    buff[9] = (UChar_t)(fOBT >> 16);
343    buff[10] = (UChar_t)(fOBT >> 8);    buff[10] = (UChar_t)(fOBT >> 8);
344    buff[11] = (UChar_t)fOBT;    buff[11] = (UChar_t)fOBT;
345    //    //
346    // Packet length    // Packet length
347    //    //
348    fLen = length;    fLen = length;
349    //    //
350    buff[12] = (UChar_t)(fLen >> 16);    buff[12] = (UChar_t)(fLen >> 16);
351    buff[13] = (UChar_t)(fLen >> 8);    buff[13] = (UChar_t)(fLen >> 8);
352    buff[14] = (UChar_t)fLen;    buff[14] = (UChar_t)fLen;
353    //    //
354    // CPU header CRC    // CPU header CRC
355    //    //
356    buff[15] = (BYTE)CM_Compute_CRC16((UINT16)0, (BYTE*)&buff, (UINT32)15);    buff[15] = (BYTE)CM_Compute_CRC16((UINT16)0, (BYTE*)&buff, (UINT32)15);
357    //    //
358    memcpy(fDataPSCU,buff,16*sizeof(UChar_t));    memcpy(fDataPSCU,buff,16*sizeof(UChar_t));
359    //    //
360  };  };
361    
362  void Digitizer::ClearCaloCalib(Int_t s){  void Digitizer::ClearCaloCalib(Int_t s){
363    //    //
364    fcstwerr[s] = 0;    fcstwerr[s] = 0;
365    fcperror[s] = 0.;    fcperror[s] = 0.;
366    for ( Int_t d=0 ; d<11 ;d++  ){    for ( Int_t d=0 ; d<11 ;d++  ){
367      Int_t pre = -1;      Int_t pre = -1;
368      for ( Int_t j=0; j<96 ;j++){      for ( Int_t j=0; j<96 ;j++){
369        if ( j%16 == 0 ) pre++;        if ( j%16 == 0 ) pre++;
370        fcalped[s][d][j] = 0.;        fcalped[s][d][j] = 0.;
371        fcstwerr[s] = 0.;        fcstwerr[s] = 0.;
372        fcperror[s] = 0.;        fcperror[s] = 0.;
373        fcalgood[s][d][j] = 0.;        fcalgood[s][d][j] = 0.;
374        fcalthr[s][d][pre] = 0.;        fcalthr[s][d][pre] = 0.;
375        fcalrms[s][d][j] = 0.;        fcalrms[s][d][j] = 0.;
376        fcalbase[s][d][pre] = 0.;        fcalbase[s][d][pre] = 0.;
377        fcalvar[s][d][pre] = 0.;        fcalvar[s][d][pre] = 0.;
378      };      };
379    };    };
380    return;    return;
381  }  }
382    
383  Int_t Digitizer::CaloLoadCalib(Int_t s,TString fcalname, UInt_t calibno){  Int_t Digitizer::CaloLoadCalib(Int_t s,TString fcalname, UInt_t calibno){
384    //    //
385    //    //
386    UInt_t e = 0;    UInt_t e = 0;
387    if ( s == 0 ) e = 0;    if ( s == 0 ) e = 0;
388    if ( s == 1 ) e = 2;    if ( s == 1 ) e = 2;
389    if ( s == 2 ) e = 3;    if ( s == 2 ) e = 3;
390    if ( s == 3 ) e = 1;    if ( s == 3 ) e = 1;
391    //    //
392    ifstream myfile;    ifstream myfile;
393    myfile.open(fcalname.Data());    myfile.open(fcalname.Data());
394    if ( !myfile ){        if ( !myfile ){    
395      return(-107);      return(-107);
396    };    };
397    myfile.close();    myfile.close();
398    //    //
399    TFile *File = new TFile(fcalname.Data());    TFile *File = new TFile(fcalname.Data());
400    if ( !File ) return(-108);    if ( !File ) return(-108);
401    TTree *tr = (TTree*)File->Get("CalibCalPed");    TTree *tr = (TTree*)File->Get("CalibCalPed");
402    if ( !tr ) return(-109);    if ( !tr ) return(-109);
403    //    //
404    TBranch *calo = tr->GetBranch("CalibCalPed");    TBranch *calo = tr->GetBranch("CalibCalPed");
405    //    //
406    pamela::CalibCalPedEvent *ce = 0;    pamela::CalibCalPedEvent *ce = 0;
407    tr->SetBranchAddress("CalibCalPed", &ce);    tr->SetBranchAddress("CalibCalPed", &ce);
408    //    //
409    Long64_t ncalibs = calo->GetEntries();    Long64_t ncalibs = calo->GetEntries();
410    //    //
411    if ( !ncalibs ) return(-110);    if ( !ncalibs ) return(-110);
412    //    //
413    calo->GetEntry(calibno);    calo->GetEntry(calibno);
414    //    //
415    if (ce->cstwerr[s] != 0 && ce->cperror[s] == 0 ) {    if (ce->cstwerr[s] != 0 && ce->cperror[s] == 0 ) {
416      fcstwerr[s] = ce->cstwerr[s];      fcstwerr[s] = ce->cstwerr[s];
417      fcperror[s] = ce->cperror[s];      fcperror[s] = ce->cperror[s];
418      for ( Int_t d=0 ; d<11 ;d++  ){      for ( Int_t d=0 ; d<11 ;d++  ){
419        Int_t pre = -1;        Int_t pre = -1;
420        for ( Int_t j=0; j<96 ;j++){        for ( Int_t j=0; j<96 ;j++){
421          if ( j%16 == 0 ) pre++;          if ( j%16 == 0 ) pre++;
422          fcalped[s][d][j] = ce->calped[e][d][j];          fcalped[s][d][j] = ce->calped[e][d][j];
423          fcalgood[s][d][j] = ce->calgood[e][d][j];          fcalgood[s][d][j] = ce->calgood[e][d][j];
424          fcalthr[s][d][pre] = ce->calthr[e][d][pre];          fcalthr[s][d][pre] = ce->calthr[e][d][pre];
425          fcalrms[s][d][j] = ce->calrms[e][d][j];          fcalrms[s][d][j] = ce->calrms[e][d][j];
426          fcalbase[s][d][pre] = ce->calbase[e][d][pre];          fcalbase[s][d][pre] = ce->calbase[e][d][pre];
427          fcalvar[s][d][pre] = ce->calvar[e][d][pre];          fcalvar[s][d][pre] = ce->calvar[e][d][pre];
428        };        };
429      };      };
430    } else {    } else {
431      printf(" CALORIMETER - ERROR: problems finding a good calibration in this file! \n\n ");      printf(" CALORIMETER - ERROR: problems finding a good calibration in this file! \n\n ");
432      File->Close();      File->Close();
433      return(-111);      return(-111);
434    };    };
435    File->Close();    File->Close();
436    return(0);    return(0);
437  }  }
438    
439    
440  void Digitizer::DigitizeCALOCALIB() {  void Digitizer::DigitizeCALOCALIB() {
441    //    //
442    // Header of the four sections    // Header of the four sections
443    //    //
444    fSecCalo[0] = 0xAA00; // XE    fSecCalo[0] = 0xAA00; // XE
445    fSecCalo[1] = 0xB100; // XO    fSecCalo[1] = 0xB100; // XO
446    fSecCalo[2] = 0xB600; // YE    fSecCalo[2] = 0xB600; // YE
447    fSecCalo[3] = 0xAD00; // YO    fSecCalo[3] = 0xAD00; // YO
448    //    //
449    // length of the data is 0x1215    // length of the data is 0x1215
450    //    //
451    fSecCALOLength[0] = 0x1215; // XE    fSecCALOLength[0] = 0x1215; // XE
452    fSecCALOLength[1] = 0x1215; // XO    fSecCALOLength[1] = 0x1215; // XO
453    fSecCALOLength[2] = 0x1215; // YE    fSecCALOLength[2] = 0x1215; // YE
454    fSecCALOLength[3] = 0x1215; // YO    fSecCALOLength[3] = 0x1215; // YO
455    //    //
456    Int_t chksum = 0;    Int_t chksum = 0;
457    UInt_t tstrip = 0;    UInt_t tstrip = 0;
458    UInt_t fSecPointer = 0;    UInt_t fSecPointer = 0;
459    //    //
460    for (Int_t sec=0; sec < 4; sec++){    for (Int_t sec=0; sec < 4; sec++){
461      //      //
462      // sec =    0 -> XE      1 -> XO        2-> YE         3 -> YO      // sec =    0 -> XE      1 -> XO        2-> YE         3 -> YO
463      //      //
464      fCALOlength = 0;      fCALOlength = 0;
465      memset(fDataCALO,0,sizeof(UShort_t)*fCALObuffer);      memset(fDataCALO,0,sizeof(UShort_t)*fCALObuffer);
466      fSecPointer = fCALOlength;      fSecPointer = fCALOlength;
467      //      //
468      // First of all we have section header and packet length      // First of all we have section header and packet length
469      //      //
470      fDataCALO[fCALOlength] = fSecCalo[sec];      fDataCALO[fCALOlength] = fSecCalo[sec];
471      fCALOlength++;      fCALOlength++;
472      fDataCALO[fCALOlength] = fSecCALOLength[sec];      fDataCALO[fCALOlength] = fSecCALOLength[sec];
473      fCALOlength++;      fCALOlength++;
474      //      //
475      // Section XO is read in the opposite direction respect to the others      // Section XO is read in the opposite direction respect to the others
476      //      //
477      chksum = 0;      chksum = 0;
478      //      //
479      for (Int_t plane=0; plane < 11; plane++){      for (Int_t plane=0; plane < 11; plane++){
480        //        //
481        if ( sec == 1 ) tstrip = fCALOlength + 96*2;        if ( sec == 1 ) tstrip = fCALOlength + 96*2;
482        //        //
483        for (Int_t strip=0; strip < 96; strip++){          for (Int_t strip=0; strip < 96; strip++){  
484          //          //
485          chksum += (Int_t)fcalped[sec][plane][strip];          chksum += (Int_t)fcalped[sec][plane][strip];
486          //          //
487          // save value          // save value
488          //          //
489          if ( sec == 1 ){          if ( sec == 1 ){
490            tstrip -= 2;            tstrip -= 2;
491            fDataCALO[tstrip] = (Int_t)fcalped[sec][plane][strip];            fDataCALO[tstrip] = (Int_t)fcalped[sec][plane][strip];
492            fDataCALO[tstrip+1] = (Int_t)fcalgood[sec][plane][strip];            fDataCALO[tstrip+1] = (Int_t)fcalgood[sec][plane][strip];
493          } else {          } else {
494            fDataCALO[fCALOlength] = (Int_t)fcalped[sec][plane][strip];            fDataCALO[fCALOlength] = (Int_t)fcalped[sec][plane][strip];
495            fDataCALO[fCALOlength+1] = (Int_t)fcalgood[sec][plane][strip];            fDataCALO[fCALOlength+1] = (Int_t)fcalgood[sec][plane][strip];
496          };          };
497          fCALOlength +=2;          fCALOlength +=2;
498        };        };
499        //        //
500      };      };
501      //      //
502      fDataCALO[fCALOlength] = (UShort_t)chksum;      fDataCALO[fCALOlength] = (UShort_t)chksum;
503      fCALOlength++;      fCALOlength++;
504      fDataCALO[fCALOlength] = 0;      fDataCALO[fCALOlength] = 0;
505      fCALOlength++;      fCALOlength++;
506      fDataCALO[fCALOlength] = (UShort_t)((Int_t)(chksum >> 16));      fDataCALO[fCALOlength] = (UShort_t)((Int_t)(chksum >> 16));
507      fCALOlength++;      fCALOlength++;
508      //      //
509      // Section XO is read in the opposite direction respect to the others      // Section XO is read in the opposite direction respect to the others
510      //      //
511      chksum = 0;      chksum = 0;
512      //      //
513      for (Int_t plane=0; plane < 11; plane++){      for (Int_t plane=0; plane < 11; plane++){
514        //        //
515        if ( sec == 1 ) tstrip = fCALOlength+6*2;        if ( sec == 1 ) tstrip = fCALOlength+6*2;
516        //        //
517        for (Int_t strip=0; strip < 6; strip++){        for (Int_t strip=0; strip < 6; strip++){
518          //          //
519          chksum += (Int_t)fcalthr[sec][plane][strip];          chksum += (Int_t)fcalthr[sec][plane][strip];
520          //          //
521          // save value          // save value
522          //          //
523          if ( sec == 1 ){          if ( sec == 1 ){
524            tstrip -= 2;            tstrip -= 2;
525            fDataCALO[tstrip] = 0;            fDataCALO[tstrip] = 0;
526            fDataCALO[tstrip+1] = (Int_t)fcalthr[sec][plane][strip];            fDataCALO[tstrip+1] = (Int_t)fcalthr[sec][plane][strip];
527          } else {          } else {
528            fDataCALO[fCALOlength] = 0;            fDataCALO[fCALOlength] = 0;
529            fDataCALO[fCALOlength+1] = (Int_t)fcalthr[sec][plane][strip];            fDataCALO[fCALOlength+1] = (Int_t)fcalthr[sec][plane][strip];
530          };          };
531          fCALOlength +=2;          fCALOlength +=2;
532        };        };
533        //        //
534      };      };
535      //      //
536      fDataCALO[fCALOlength] = 0;      fDataCALO[fCALOlength] = 0;
537      fCALOlength++;      fCALOlength++;
538      fDataCALO[fCALOlength] = (UShort_t)chksum;      fDataCALO[fCALOlength] = (UShort_t)chksum;
539      fCALOlength++;      fCALOlength++;
540      fDataCALO[fCALOlength] = 0;      fDataCALO[fCALOlength] = 0;
541      fCALOlength++;      fCALOlength++;
542      fDataCALO[fCALOlength] = (UShort_t)((Int_t)(chksum >> 16));      fDataCALO[fCALOlength] = (UShort_t)((Int_t)(chksum >> 16));
543      fCALOlength++;      fCALOlength++;
544      //      //
545      // Section XO is read in the opposite direction respect to the others      // Section XO is read in the opposite direction respect to the others
546      //      //
547      for (Int_t plane=0; plane < 11; plane++){      for (Int_t plane=0; plane < 11; plane++){
548        //        //
549        if ( sec == 1 ) tstrip = fCALOlength+96*2;        if ( sec == 1 ) tstrip = fCALOlength+96*2;
550        //        //
551        for (Int_t strip=0; strip < 96; strip++){        for (Int_t strip=0; strip < 96; strip++){
552          //          //
553          // save value          // save value
554          //          //
555          if ( sec == 1 ){          if ( sec == 1 ){
556            tstrip -= 2;            tstrip -= 2;
557            fDataCALO[tstrip] = 0;            fDataCALO[tstrip] = 0;
558            fDataCALO[tstrip+1] = (Int_t)fcalrms[sec][plane][strip];            fDataCALO[tstrip+1] = (Int_t)fcalrms[sec][plane][strip];
559          } else {          } else {
560            fDataCALO[fCALOlength] = 0;            fDataCALO[fCALOlength] = 0;
561            fDataCALO[fCALOlength+1] = (Int_t)fcalrms[sec][plane][strip];            fDataCALO[fCALOlength+1] = (Int_t)fcalrms[sec][plane][strip];
562          };          };
563          fCALOlength += 2;          fCALOlength += 2;
564        };        };
565        //        //
566      };          };    
567      //      //
568      // Section XO is read in the opposite direction respect to the others      // Section XO is read in the opposite direction respect to the others
569      //      //
570      for (Int_t plane=0; plane < 11; plane++){      for (Int_t plane=0; plane < 11; plane++){
571        //        //
572        if ( sec == 1 ) tstrip = fCALOlength+6*4;        if ( sec == 1 ) tstrip = fCALOlength+6*4;
573        //        //
574        for (Int_t strip=0; strip < 6; strip++){        for (Int_t strip=0; strip < 6; strip++){
575          //          //
576          // save value          // save value
577          //          //
578          if ( sec == 1 ){          if ( sec == 1 ){
579            tstrip -= 4;            tstrip -= 4;
580            fDataCALO[tstrip] = 0;            fDataCALO[tstrip] = 0;
581            fDataCALO[tstrip+1] = (Int_t)fcalbase[sec][plane][strip];            fDataCALO[tstrip+1] = (Int_t)fcalbase[sec][plane][strip];
582            fDataCALO[tstrip+2] = 0;            fDataCALO[tstrip+2] = 0;
583            fDataCALO[tstrip+3] = (Int_t)fcalvar[sec][plane][strip];            fDataCALO[tstrip+3] = (Int_t)fcalvar[sec][plane][strip];
584          } else {          } else {
585            fDataCALO[fCALOlength] = 0;            fDataCALO[fCALOlength] = 0;
586            fDataCALO[fCALOlength+1] = (Int_t)fcalbase[sec][plane][strip];            fDataCALO[fCALOlength+1] = (Int_t)fcalbase[sec][plane][strip];
587            fDataCALO[fCALOlength+2] = 0;            fDataCALO[fCALOlength+2] = 0;
588            fDataCALO[fCALOlength+3] = (Int_t)fcalvar[sec][plane][strip];            fDataCALO[fCALOlength+3] = (Int_t)fcalvar[sec][plane][strip];
589          };          };
590          fCALOlength +=4;          fCALOlength +=4;
591        };        };
592        //        //
593      };            };      
594      //      //
595      //      //
596      // here we calculate and save the CRC      // here we calculate and save the CRC
597      //      //
598      fDataCALO[fCALOlength] = 0;      fDataCALO[fCALOlength] = 0;
599      fCALOlength++;      fCALOlength++;
600      Short_t CRC = 0;      Short_t CRC = 0;
601      for (UInt_t i=0; i<(fCALOlength-fSecPointer); i++){      for (UInt_t i=0; i<(fCALOlength-fSecPointer); i++){
602        CRC=crc(CRC,fDataCALO[i+fSecPointer]);        CRC=crc(CRC,fDataCALO[i+fSecPointer]);
603      };      };
604      fDataCALO[fCALOlength] = (UShort_t)CRC;      fDataCALO[fCALOlength] = (UShort_t)CRC;
605      fCALOlength++;      fCALOlength++;
606      //      //
607      UInt_t length=fCALOlength*2;      UInt_t length=fCALOlength*2;
608      DigitizePSCU(length,0x18);      DigitizePSCU(length,0x18);
609      //      //
610      // Add padding to 64 bits      // Add padding to 64 bits
611      //      //
612      AddPadding();      AddPadding();
613      //      //
614      fOutputfile.write(reinterpret_cast<char*>(fDataPSCU),sizeof(UShort_t)*fPSCUbuffer);      fOutputfile.write(reinterpret_cast<char*>(fDataPSCU),sizeof(UShort_t)*fPSCUbuffer);
615      UShort_t temp[1000000];      UShort_t temp[1000000];
616      memset(temp,0,sizeof(UShort_t)*1000000);      memset(temp,0,sizeof(UShort_t)*1000000);
617      swab(fDataCALO,temp,sizeof(UShort_t)*fCALOlength);  // WE MUST SWAP THE BYTES!!!      swab(fDataCALO,temp,sizeof(UShort_t)*fCALOlength);  // WE MUST SWAP THE BYTES!!!
618      fOutputfile.write(reinterpret_cast<char*>(temp),sizeof(UShort_t)*fCALOlength);      fOutputfile.write(reinterpret_cast<char*>(temp),sizeof(UShort_t)*fCALOlength);
619      //      //
620      // padding to 64 bytes      // padding to 64 bytes
621      //      //
622      if ( fPadding ){      if ( fPadding ){
623        fOutputfile.write(reinterpret_cast<char*>(fDataPadding),sizeof(UChar_t)*fPadding);        fOutputfile.write(reinterpret_cast<char*>(fDataPadding),sizeof(UChar_t)*fPadding);
624      };      };
625      //      //
626      //          //    
627    };    };
628    //    //
629  };  };
630    
631  void Digitizer::CaloLoadCalib() {  void Digitizer::CaloLoadCalib() {
632    //    //
633    fGivenCaloCalib = 0; //                                  ####@@@@ should be given as input par @@@@####    fGivenCaloCalib = 0; //                                  ####@@@@ should be given as input par @@@@####
634    //    //
635    // first of all load the MIP to ADC conversion values    // first of all load the MIP to ADC conversion values
636    //    //
637    stringstream calfile;    stringstream calfile;
638    Int_t error = 0;    Int_t error = 0;
639    GL_PARAM *glparam = new GL_PARAM();    GL_PARAM *glparam = new GL_PARAM();
640    //    //
641    // determine where I can find calorimeter ADC to MIP conversion file      // determine where I can find calorimeter ADC to MIP conversion file  
642    //    //
643    error = 0;    error = 0;
644    error = glparam->Query_GL_PARAM(3,101,fDbc);    error = glparam->Query_GL_PARAM(3,101,fDbc);
645    //    //
646    calfile.str("");    calfile.str("");
647    calfile << glparam->PATH.Data() << "/";    calfile << glparam->PATH.Data() << "/";
648    calfile << glparam->NAME.Data();    calfile << glparam->NAME.Data();
649    //    //
650    printf("\n Using Calorimeter ADC to MIP conversion file: \n %s \n",calfile.str().c_str());    printf("\n Using Calorimeter ADC to MIP conversion file: \n %s \n",calfile.str().c_str());
651    FILE *f;    FILE *f;
652    f = fopen(calfile.str().c_str(),"rb");    f = fopen(calfile.str().c_str(),"rb");
653    //    //
654    memset(fCalomip,0,4224*sizeof(fCalomip[0][0][0]));    memset(fCalomip,0,4224*sizeof(fCalomip[0][0][0]));
655    //    //
656    for (Int_t m = 0; m < 2 ; m++ ){    for (Int_t m = 0; m < 2 ; m++ ){
657      for (Int_t k = 0; k < 22; k++ ){      for (Int_t k = 0; k < 22; k++ ){
658        for (Int_t l = 0; l < 96; l++ ){        for (Int_t l = 0; l < 96; l++ ){
659          fread(&fCalomip[m][k][l],sizeof(fCalomip[m][k][l]),1,f);          fread(&fCalomip[m][k][l],sizeof(fCalomip[m][k][l]),1,f);
660        };        };
661      };      };
662    };    };
663    fclose(f);    fclose(f);
664    //    //
665    // determine which calibration has to be used and load it for each section    // determine which calibration has to be used and load it for each section
666    //      //  
667    GL_CALO_CALIB *glcalo = new GL_CALO_CALIB();    GL_CALO_CALIB *glcalo = new GL_CALO_CALIB();
668    GL_ROOT *glroot = new GL_ROOT();      GL_ROOT *glroot = new GL_ROOT();  
669    TString fcalname;    TString fcalname;
670    UInt_t idcalib;    UInt_t idcalib;
671    UInt_t calibno;    UInt_t calibno;
672    UInt_t utime = 0;    UInt_t utime = 0;
673    //    //
674    for (UInt_t s=0; s<4; s++){    for (UInt_t s=0; s<4; s++){
675      //      //
676      // clear calo calib variables for section s      // clear calo calib variables for section s
677      //      //
678      ClearCaloCalib(s);      ClearCaloCalib(s);
679      //      //
680      if ( fGivenCaloCalib ){      if ( fGivenCaloCalib ){
681        //        //
682        // a time has been given as input on the command line so retrieve the calibration that preceed that time        // a time has been given as input on the command line so retrieve the calibration that preceed that time
683        //        //
684        glcalo->Query_GL_CALO_CALIB(fGivenCaloCalib,utime,s,fDbc);        glcalo->Query_GL_CALO_CALIB(fGivenCaloCalib,utime,s,fDbc);
685        //          //  
686        calibno = glcalo->EV_ROOT;        calibno = glcalo->EV_ROOT;
687        idcalib = glcalo->ID_ROOT_L0;        idcalib = glcalo->ID_ROOT_L0;
688        //        //
689        // determine path and name and entry of the calibration file        // determine path and name and entry of the calibration file
690        //        //
691        printf("\n");        printf("\n");
692        printf(" ** SECTION %i **\n",s);        printf(" ** SECTION %i **\n",s);
693        //        //
694        glroot->Query_GL_ROOT(idcalib,fDbc);        glroot->Query_GL_ROOT(idcalib,fDbc);
695        //        //
696        stringstream name;        stringstream name;
697        name.str("");        name.str("");
698        name << glroot->PATH.Data() << "/";        name << glroot->PATH.Data() << "/";
699        name << glroot->NAME.Data();        name << glroot->NAME.Data();
700        //        //
701        fcalname = (TString)name.str().c_str();        fcalname = (TString)name.str().c_str();
702        //        //
703        printf("\n Section %i : using  file %s calibration at entry %i: \n",s,fcalname.Data(),calibno);        printf("\n Section %i : using  file %s calibration at entry %i: \n",s,fcalname.Data(),calibno);
704        //        //
705      } else {      } else {
706        error = 0;        error = 0;
707        error = glparam->Query_GL_PARAM(1,104,fDbc);        error = glparam->Query_GL_PARAM(1,104,fDbc);
708        //        //
709        calfile.str("");        calfile.str("");
710        calfile << glparam->PATH.Data() << "/";        calfile << glparam->PATH.Data() << "/";
711        calfile << glparam->NAME.Data();        calfile << glparam->NAME.Data();
712        //        //
713        printf("\n Section %i : using default calorimeter calibration: \n %s \n",s,calfile.str().c_str());        printf("\n Section %i : using default calorimeter calibration: \n %s \n",s,calfile.str().c_str());
714        //        //
715        fcalname = (TString)calfile.str().c_str();        fcalname = (TString)calfile.str().c_str();
716        calibno = s;        calibno = s;
717        //        //
718      };      };
719      //      //
720      // load calibration variables in memory      // load calibration variables in memory
721      //      //
722      CaloLoadCalib(s,fcalname,calibno);      CaloLoadCalib(s,fcalname,calibno);
723      //      //
724    };    };
725    //    //
726    // 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    // 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
727    //    //
728    delete glparam;    delete glparam;
729    delete glcalo;    delete glcalo;
730    delete glroot;    delete glroot;
731  };  };
732    
733  void Digitizer::DigitizeCALO() {  void Digitizer::DigitizeCALO() {
734    //    //
735    fModCalo = 0; // 0 is RAW, 1 is COMPRESS, 2 is FULL     ####@@@@ should be given as input par @@@@####    fModCalo = 0; // 0 is RAW, 1 is COMPRESS, 2 is FULL     ####@@@@ should be given as input par @@@@####
736    //    //
737    //    //
738    //    //
739    fCALOlength = 0;  // reset total dimension of calo data    fCALOlength = 0;  // reset total dimension of calo data
740    //    //
741    // gpamela variables to be used    // gpamela variables to be used
742    //    //
743    fhBookTree->SetBranchStatus("Nthcali",1);    fhBookTree->SetBranchStatus("Nthcali",1);
744    fhBookTree->SetBranchStatus("Icaplane",1);    fhBookTree->SetBranchStatus("Icaplane",1);
745    fhBookTree->SetBranchStatus("Icastrip",1);    fhBookTree->SetBranchStatus("Icastrip",1);
746    fhBookTree->SetBranchStatus("Icamod",1);    fhBookTree->SetBranchStatus("Icamod",1);
747    fhBookTree->SetBranchStatus("Enestrip",1);    fhBookTree->SetBranchStatus("Enestrip",1);
748    //    //
749    // call different routines depending on the acq mode you want to simulate    // call different routines depending on the acq mode you want to simulate
750    //    //
751    switch ( fModCalo ){    switch ( fModCalo ){
752    case 0:    case 0:
753      this->DigitizeCALORAW();      this->DigitizeCALORAW();
754      break;      break;
755    case 1:    case 1:
756      this->DigitizeCALOCOMPRESS();      this->DigitizeCALOCOMPRESS();
757      break;      break;
758    case 2:    case 2:
759      this->DigitizeCALOFULL();      this->DigitizeCALOFULL();
760      break;      break;
761    };    };
762    //    //
763  };  };
764    
765  Float_t Digitizer::GetCALOen(Int_t sec, Int_t plane, Int_t strip){  Float_t Digitizer::GetCALOen(Int_t sec, Int_t plane, Int_t strip){
766    //    //
767    // determine plane and strip    // determine plane and strip
768    //    //
769    Int_t mplane = 0;    Int_t mplane = 0;
770    //    //
771    // wrong!    // wrong!
772    //    //
773    //  if ( sec == 0 || sec == 3 ) mplane = (plane * 4) + sec + 1;      //  if ( sec == 0 || sec == 3 ) mplane = (plane * 4) + sec + 1;  
774    //  if ( sec == 1 ) mplane = (plane * 4) + 2 + 1;      //  if ( sec == 1 ) mplane = (plane * 4) + 2 + 1;  
775    //  if ( sec == 2 ) mplane = (plane * 4) + 1 + 1;      //  if ( sec == 2 ) mplane = (plane * 4) + 1 + 1;  
776    //    //
777    if ( sec == 0 ) mplane = plane * 4 + 1; // it must be 0, 4, 8, ... (+1)  from plane = 0, 11    if ( sec == 0 ) mplane = plane * 4 + 1; // it must be 0, 4, 8, ... (+1)  from plane = 0, 11
778    if ( sec == 1 ) mplane = plane * 4 + 2 + 1; // it must be 2, 6, 10, ... (+1)  from plane = 0, 11    if ( sec == 1 ) mplane = plane * 4 + 2 + 1; // it must be 2, 6, 10, ... (+1)  from plane = 0, 11
779    if ( sec == 2 ) mplane = plane * 4 + 3 + 1; // it must be 3, 7, 11, ... (+1)  from plane = 0, 11    if ( sec == 2 ) mplane = plane * 4 + 3 + 1; // it must be 3, 7, 11, ... (+1)  from plane = 0, 11
780    if ( sec == 3 ) mplane = plane * 4 + 1 + 1; // it must be 1, 5, 9, ... (+1)  from plane = 0, 11    if ( sec == 3 ) mplane = plane * 4 + 1 + 1; // it must be 1, 5, 9, ... (+1)  from plane = 0, 11
781    //    //
782    Int_t mstrip = strip + 1;    Int_t mstrip = strip + 1;
783    //    //
784    // search energy release in gpamela output    // search energy release in gpamela output
785    //    //
786    for (Int_t i=0; i<Nthcali;i++){    for (Int_t i=0; i<Nthcali;i++){
787      if ( Icaplane[i] == mplane && Icastrip[i] == mstrip ){      if ( Icaplane[i] == mplane && Icastrip[i] == mstrip ){
788        return (Enestrip[i]);        return (Enestrip[i]);
789      };      };
790    };    };
791    //    //
792    // if not found it means no energy release so return 0.    // if not found it means no energy release so return 0.
793    //    //
794    return(0.);    return(0.);
795  };  };
796    
797  void Digitizer::DigitizeCALORAW() {  void Digitizer::DigitizeCALORAW() {
798    //    //
799    // some variables    // some variables
800    //    //
801    Float_t ens = 0.;    Float_t ens = 0.;
802    UInt_t adcsig = 0;    UInt_t adcsig = 0;
803    UInt_t adcbase = 0;    UInt_t adcbase = 0;
804    UInt_t adc = 0;    UInt_t adc = 0;
805    Int_t pre = 0;    Int_t pre = 0;
806    UInt_t l = 0;    UInt_t l = 0;
807    UInt_t lpl = 0;    UInt_t lpl = 0;
808    UInt_t tstrip = 0;    UInt_t tstrip = 0;
809    UInt_t fSecPointer = 0;    UInt_t fSecPointer = 0;
810    Double_t pedenoise;    Double_t pedenoise;
811    Float_t rms = 0.;    Float_t rms = 0.;
812    Float_t pedestal = 0.;    Float_t pedestal = 0.;
813    //    //
814    // clean the data structure    // clean the data structure
815    //    //
816    memset(fDataCALO,0,sizeof(UShort_t)*fCALObuffer);    memset(fDataCALO,0,sizeof(UShort_t)*fCALObuffer);
817    //    //
818    // Header of the four sections    // Header of the four sections
819    //    //
820    fSecCalo[0] = 0xEA08; // XE    fSecCalo[0] = 0xEA08; // XE
821    fSecCalo[1] = 0xF108; // XO    fSecCalo[1] = 0xF108; // XO
822    fSecCalo[2] = 0xF608; // YE    fSecCalo[2] = 0xF608; // YE
823    fSecCalo[3] = 0xED08; // YO    fSecCalo[3] = 0xED08; // YO
824    //    //
825    // length of the data is 0x0428 in RAW mode    // length of the data is 0x0428 in RAW mode
826    //    //
827    fSecCALOLength[0] = 0x0428; // XE    fSecCALOLength[0] = 0x0428; // XE
828    fSecCALOLength[1] = 0x0428; // XO    fSecCALOLength[1] = 0x0428; // XO
829    fSecCALOLength[2] = 0x0428; // YE    fSecCALOLength[2] = 0x0428; // YE
830    fSecCALOLength[3] = 0x0428; // YO    fSecCALOLength[3] = 0x0428; // YO
831    //    //
832    // let's start    // let's start
833    //    //
834    fCALOlength = 0;    fCALOlength = 0;
835    //    //
836    for (Int_t sec=0; sec < 4; sec++){    for (Int_t sec=0; sec < 4; sec++){
837      //      //
838      // sec =    0 -> XE      1 -> XO        2-> YE         3 -> YO      // sec =    0 -> XE      1 -> XO        2-> YE         3 -> YO
839      //      //
840      l = 0;                 // XE and XO are Y planes      l = 0;                 // XE and XO are Y planes
841      if ( sec < 2 ) l = 1;  // while YE and  YO are X planes      if ( sec < 2 ) l = 1;  // while YE and  YO are X planes
842      //      //
843      fSecPointer = fCALOlength;      fSecPointer = fCALOlength;
844      //      //
845      // First of all we have section header and packet length      // First of all we have section header and packet length
846      //      //
847      fDataCALO[fCALOlength] = fSecCalo[sec];      fDataCALO[fCALOlength] = fSecCalo[sec];
848      fCALOlength++;      fCALOlength++;
849      fDataCALO[fCALOlength] = fSecCALOLength[sec];      fDataCALO[fCALOlength] = fSecCALOLength[sec];
850      fCALOlength++;      fCALOlength++;
851      //      //
852      // selftrigger coincidences - in the future we should add here some code to simulate timing response of pre-amplifiers      // selftrigger coincidences - in the future we should add here some code to simulate timing response of pre-amplifiers
853      //      //
854      for (Int_t autoplane=0; autoplane < 7; autoplane++){      for (Int_t autoplane=0; autoplane < 7; autoplane++){
855        fDataCALO[fCALOlength] = 0x0000;        fDataCALO[fCALOlength] = 0x0000;
856        fCALOlength++;        fCALOlength++;
857      };      };
858      //      //
859      //      //
860      // here comes data      // here comes data
861      //      //
862      //      //
863      // Section XO is read in the opposite direction respect to the others      // Section XO is read in the opposite direction respect to the others
864      //      //
865      if ( sec == 1 ){            if ( sec == 1 ){      
866        tstrip = 96*11 + fCALOlength;        tstrip = 96*11 + fCALOlength;
867      } else {      } else {
868        tstrip = 0;        tstrip = 0;
869      };      };
870      //      //
871      pre = -1;      pre = -1;
872      //      //
873      for (Int_t strip=0; strip < 96; strip++){        for (Int_t strip=0; strip < 96; strip++){  
874        //        //
875        // which is the pre for this strip?        // which is the pre for this strip?
876        //        //
877        if (strip%16 == 0) {        if (strip%16 == 0) {
878          pre++;          pre++;
879        };        };
880        //        //
881        if ( sec == 1 ) tstrip -= 11;        if ( sec == 1 ) tstrip -= 11;
882        //        //
883        for (Int_t plane=0; plane < 11; plane++){        for (Int_t plane=0; plane < 11; plane++){
884          //          //
885          // here is wrong!!!!          // here is wrong!!!!
886          //          //
887          //          //
888          //      if ( plane%2 == 0 && sec%2 != 0){          //      if ( plane%2 == 0 && sec%2 != 0){
889          //        lpl = plane*2;          //        lpl = plane*2;
890          //      } else {          //      } else {
891          //        lpl = (plane*2) + 1;          //        lpl = (plane*2) + 1;
892          //      };          //      };
893          //          //
894          if ( sec == 0 || sec == 3 ) lpl = plane * 2;          if ( sec == 0 || sec == 3 ) lpl = plane * 2;
895          if ( sec == 1 || sec == 2 ) lpl = (plane * 2) + 1;          if ( sec == 1 || sec == 2 ) lpl = (plane * 2) + 1;
896          //          //
897          // get the energy in GeV from the simulation for that strip          // get the energy in GeV from the simulation for that strip
898          //          //
899          ens = this->GetCALOen(sec,plane,strip);          ens = this->GetCALOen(sec,plane,strip);
900          //          //
901          // convert it into ADC channels          // convert it into ADC channels
902          //                //      
903          adcsig = int(ens*fCalomip[l][lpl][strip]/fCALOGeV2MIPratio);          adcsig = int(ens*fCalomip[l][lpl][strip]/fCALOGeV2MIPratio);
904          //          //
905          // sum baselines          // sum baselines
906          //          //
907          adcbase = (UInt_t)fcalbase[sec][plane][pre];          adcbase = (UInt_t)fcalbase[sec][plane][pre];
908          //          //
909          // add noise and pedestals          // add noise and pedestals
910          //                //      
911          pedestal = fcalped[sec][plane][strip];          pedestal = fcalped[sec][plane][strip];
912          rms = fcalrms[sec][plane][strip]/4.;          rms = fcalrms[sec][plane][strip]/4.;
913          //          //
914          // Add random gaussian noise of RMS rms and Centered in the pedestal          // Add random gaussian noise of RMS rms and Centered in the pedestal
915          //          //
916          pedenoise = gRandom->Gaus((Double_t)pedestal,(Double_t)rms);          pedenoise = gRandom->Gaus((Double_t)pedestal,(Double_t)rms);
917          //          //
918          // Sum all contribution          // Sum all contribution
919          //          //
920          adc = adcsig + adcbase + (Int_t)round(pedenoise);          adc = adcsig + adcbase + (Int_t)round(pedenoise);
921          //          //
922          // Signal saturation          // Signal saturation
923          //          //
924          if ( adc > 0x7FFF ) adc = 0x7FFF;          if ( adc > 0x7FFF ) adc = 0x7FFF;
925          //          //
926          // save value          // save value
927          //          //
928          if ( sec == 1 ){          if ( sec == 1 ){
929            fDataCALO[tstrip] = adc;            fDataCALO[tstrip] = adc;
930            tstrip++;            tstrip++;
931          } else {          } else {
932            fDataCALO[fCALOlength] = adc;            fDataCALO[fCALOlength] = adc;
933          };          };
934          fCALOlength++;          fCALOlength++;
935          //          //
936        };        };
937        //        //
938        if ( sec == 1 ) tstrip -= 11;        if ( sec == 1 ) tstrip -= 11;
939        //        //
940      };      };
941      //      //
942      // here we calculate and save the CRC      // here we calculate and save the CRC
943      //      //
944      Short_t CRC = 0;      Short_t CRC = 0;
945      for (UInt_t i=0; i<(fCALOlength-fSecPointer); i++){      for (UInt_t i=0; i<(fCALOlength-fSecPointer); i++){
946        CRC=crc(CRC,fDataCALO[i+fSecPointer]);        CRC=crc(CRC,fDataCALO[i+fSecPointer]);
947      };      };
948      fDataCALO[fCALOlength] = (UShort_t)CRC;      fDataCALO[fCALOlength] = (UShort_t)CRC;
949      fCALOlength++;      fCALOlength++;
950      //      //
951    };    };
952    //    //
953    //   for (Int_t i=0; i<fCALOlength; i++){    //   for (Int_t i=0; i<fCALOlength; i++){
954    //     printf(" WORD %i       DIGIT %0x   \n",i,fDataCALO[i]);    //     printf(" WORD %i       DIGIT %0x   \n",i,fDataCALO[i]);
955    //   };    //   };
956    //    //
957  };  };
958    
959  void Digitizer::DigitizeCALOCOMPRESS() {  void Digitizer::DigitizeCALOCOMPRESS() {
960    //    //
961    printf(" COMPRESS MODE STILL NOT IMPLEMENTED! \n");      printf(" COMPRESS MODE STILL NOT IMPLEMENTED! \n");  
962    //    //
963    this->DigitizeCALORAW();    this->DigitizeCALORAW();
964    return;    return;
965    //    //
966    //    //
967    //    //
968    fSecCalo[0] = 0xEA00;    fSecCalo[0] = 0xEA00;
969    fSecCalo[1] = 0xF100;    fSecCalo[1] = 0xF100;
970    fSecCalo[2] = 0xF600;    fSecCalo[2] = 0xF600;
971    fSecCalo[3] = 0xED00;    fSecCalo[3] = 0xED00;
972    //    //
973    // length of the data in DSP mode must be calculated on fly during digitization    // length of the data in DSP mode must be calculated on fly during digitization
974    //    //
975    memset(fSecCALOLength,0x0,4*sizeof(UShort_t));    memset(fSecCALOLength,0x0,4*sizeof(UShort_t));
976    //    //
977    // here comes raw data    // here comes raw data
978    //    //
979    Int_t en = 0;    Int_t en = 0;
980    //    //
981    for (Int_t sec=0; sec < 4; sec++){    for (Int_t sec=0; sec < 4; sec++){
982      fDataCALO[en] = fSecCalo[sec];      fDataCALO[en] = fSecCalo[sec];
983      en++;      en++;
984      fDataCALO[en] = fSecCALOLength[sec];      fDataCALO[en] = fSecCALOLength[sec];
985      en++;      en++;
986      for (Int_t plane=0; plane < 11; plane++){      for (Int_t plane=0; plane < 11; plane++){
987        for (Int_t strip=0; strip < 11; strip++){        for (Int_t strip=0; strip < 11; strip++){
988          fDataCALO[en] = 0x0;          fDataCALO[en] = 0x0;
989          en++;          en++;
990        };        };
991      };      };
992    };    };
993    //    //
994  };  };
995    
996  void Digitizer::DigitizeCALOFULL() {  void Digitizer::DigitizeCALOFULL() {
997    //    //
998    printf(" FULL MODE STILL NOT IMPLEMENTED! \n");      printf(" FULL MODE STILL NOT IMPLEMENTED! \n");  
999    //    //
1000    this->DigitizeCALORAW();    this->DigitizeCALORAW();
1001    return;    return;
1002    //    //
1003    fSecCalo[0] = 0xEA00;    fSecCalo[0] = 0xEA00;
1004    fSecCalo[1] = 0xF100;    fSecCalo[1] = 0xF100;
1005    fSecCalo[2] = 0xF600;    fSecCalo[2] = 0xF600;
1006    fSecCalo[3] = 0xED00;    fSecCalo[3] = 0xED00;
1007    //    //
1008    // length of the data in DSP mode must be calculated on fly during digitization    // length of the data in DSP mode must be calculated on fly during digitization
1009    //    //
1010    memset(fSecCALOLength,0x0,4*sizeof(UShort_t));    memset(fSecCALOLength,0x0,4*sizeof(UShort_t));
1011    //    //
1012    // here comes raw data    // here comes raw data
1013    //    //
1014    Int_t  en = 0;    Int_t  en = 0;
1015    //    //
1016    for (Int_t sec=0; sec < 4; sec++){    for (Int_t sec=0; sec < 4; sec++){
1017      fDataCALO[en] = fSecCalo[sec];      fDataCALO[en] = fSecCalo[sec];
1018      en++;      en++;
1019      fDataCALO[en] = fSecCALOLength[sec];      fDataCALO[en] = fSecCALOLength[sec];
1020      en++;      en++;
1021      for (Int_t plane=0; plane < 11; plane++){      for (Int_t plane=0; plane < 11; plane++){
1022        for (Int_t strip=0; strip < 11; strip++){        for (Int_t strip=0; strip < 11; strip++){
1023          fDataCALO[en] = 0x0;          fDataCALO[en] = 0x0;
1024          en++;          en++;
1025        };        };
1026      };      };
1027    };    };
1028    //    //
1029  };  };
1030    
1031  void Digitizer::DigitizeTRIGGER() {  void Digitizer::DigitizeTRIGGER() {
1032    //fDataTrigger: 153 bytes    //fDataTrigger: 153 bytes
1033    for (Int_t j=0; j < 153; j++)    for (Int_t j=0; j < 153; j++)
1034      fDataTrigger[0]=0x00;      fDataTrigger[j]=0x00;
1035  };  };
1036    
1037  Int_t Digitizer::DigitizeTOF() {  Int_t Digitizer::DigitizeTOF() {
1038    //fDataTof: 12 x 23 bytes (=276 bytes)    //fDataTof: 12 x 23 bytes (=276 bytes)
1039    UChar_t *pTof=fDataTof;    UChar_t *pTof=fDataTof;
1040      Bool_t DEBUG=false;
1041    // --- activate branches:  
1042    fhBookTree->SetBranchStatus("Nthtof",1);    // --- activate branches:
1043    fhBookTree->SetBranchStatus("Ipltof",1);    fhBookTree->SetBranchStatus("Nthtof",1);
1044    fhBookTree->SetBranchStatus("Ipaddle",1);    fhBookTree->SetBranchStatus("Ipltof",1);
1045    fhBookTree->SetBranchStatus("Xintof",1);    fhBookTree->SetBranchStatus("Ipaddle",1);
1046    fhBookTree->SetBranchStatus("Yintof",1);    fhBookTree->SetBranchStatus("Xintof",1);
1047    fhBookTree->SetBranchStatus("Xouttof",1);    fhBookTree->SetBranchStatus("Yintof",1);
1048    fhBookTree->SetBranchStatus("Youttof",1);    fhBookTree->SetBranchStatus("Xouttof",1);
1049    fhBookTree->SetBranchStatus("Ereltof",1);    fhBookTree->SetBranchStatus("Youttof",1);
1050    fhBookTree->SetBranchStatus("Timetof",1);    fhBookTree->SetBranchStatus("Ereltof",1);
1051    // not yet used: Zintof, Xouttof, Youttof, Zouttof    fhBookTree->SetBranchStatus("Timetof",1);
1052      // not yet used: Zintof, Xouttof, Youttof, Zouttof
1053    // ------ evaluate energy in each pmt: ------  
1054    // strip geometry (lenght/width)    // ------ evaluate energy in each pmt: ------
1055    Float_t dimel[6] = {33.0, 40.8 ,18.0, 15.0, 15.0, 18.0};    // strip geometry (lenght/width)
1056    //Float_t dimes[6] = {5.1, 5.5, 7.5, 9.0, 6.0, 5.0};    Float_t dimel[6] = {33.0, 40.8 ,18.0, 15.0, 15.0, 18.0};
1057        //Float_t dimes[6] = {5.1, 5.5, 7.5, 9.0, 6.0, 5.0};
1058    //  S11 8 paddles  33.0 x 5.1 cm    
1059    //  S12 6 paddles  40.8 x 5.5 cm    //  S11 8 paddles  33.0 x 5.1 cm
1060    //  S21 2 paddles  18.0 x 7.5 cm    //  S12 6 paddles  40.8 x 5.5 cm
1061    //  S22 2 paddles  15.0 x 9.0 cm    //  S21 2 paddles  18.0 x 7.5 cm
1062    //  S31 3 paddles  15.0 x 6.0 cm    //  S22 2 paddles  15.0 x 9.0 cm
1063    //  S32 3 paddles  18.0 x 5.0 cm    //  S31 3 paddles  15.0 x 6.0 cm
1064      //  S32 3 paddles  18.0 x 5.0 cm
1065    // distance from the interaction point to the pmts (right,left)  
1066    Float_t xpath[2]={0., 0.}; /*path(cm) in X per S12,S21,S32 verso il pmt DX o SX*/    Float_t FGeo[2]={0., 0.}; /* geometrical factor */
1067    Float_t ypath[2]={0., 0.}; /*path(cm) in Y per S11,S22,S31 verso il pmt DX o SX*/  
1068    Float_t FGeo[2]={0., 0.}; /* fattore geometrico */    const Float_t Pho_keV = 10.;     // photons per keV in scintillator
1069      const Float_t echarge = 1.6e-19; // electron charge
1070    const Float_t Pho_keV = 10.;     // photons per keV in scintillator    Float_t Npho=0.;
1071    const Float_t echarge = 1.6e-19; // carica dell'elettrone    Float_t QevePmt_pC[48];
1072    Float_t Npho=0.;    Float_t QhitPad_pC[2]={0., 0.};
1073    Float_t QevePmt_pC[48];    Float_t QhitPmt_pC[2]={0., 0.};
1074    Float_t QhitPad_pC[2]={0., 0.};    Float_t pmGain = 3.5e6;  /* PMT Gain: the same for all PMTs */
1075    Float_t QhitPmt_pC[2]={0., 0.};    Float_t effi=0.21;       /* Efficienza di fotocatodo */
1076    Float_t pmGain = 3.5e6;  /* Gain: per il momento uguale per tutti */  
1077    Float_t effi=0.21;       /* Efficienza di fotocatodo */    //   Float_t ADC_pC0=-58.1;      //  ADC/pC conversion coefficient 0
1078    Float_t ADC_pC=1.666667; // ADC_ch/pC conversion = 0.6 pC/channel (+30 di offset)    //   Float_t ADC_pC1=1.728;      //  ADC/pC conversion coefficient 1
1079    Float_t ADCoffset=30.;    //   Float_t ADC_pC2=-4.063e-05; //  ADC/pC conversion coefficient 2
1080    Int_t ADClast=4095;      // no signal --> ADC ch=4095    //   Float_t ADC_pC3=-5.763e-08; //  ADC/pC conversion coefficient 3
1081    Int_t ADCtof[48];    
1082    //Float_t ADCsat=3100;  ci pensiamo in futuro !    // pC < 800
1083    //Float_t pCsat=2500;      Float_t ADC_pC0A =   -4.437616e+01   ;
1084    for(Int_t i=0; i<48; i++){    Float_t ADC_pC1A =   1.573329e+00    ;
1085      QevePmt_pC[i] = 0;    Float_t ADC_pC2A =   2.780518e-04  ;
1086      ADCtof[i]=0;    Float_t ADC_pC3A =  -2.302160e-07 ;
1087    }    
1088        // pC   > 800:
1089    // ------ read calibration file (get A1, A2, lambda1, lambda2)    Float_t ADC_pC0B =    -2.245756e+02  ;
1090    ifstream fileTriggerCalib;    Float_t ADC_pC1B =     2.184156e+00  ;
1091    TString ftrigname="TrigCalibParam.txt";    Float_t ADC_pC2B =     -4.171825e-04  ;
1092    fileTriggerCalib.open(ftrigname.Data());    Float_t ADC_pC3B =     3.789715e-08  ;
1093    if ( !fileTriggerCalib ) {  
1094      printf("debug: no trigger calib file!\n");    Float_t pCthres=40.;     // threshold in charge
1095      return(-117); //check output!    Int_t ADClast=4095;      // no signal --> ADC ch=4095
1096    };    Int_t ADCsat=3100;       // saturation value for the ADCs
1097    Float_t atte1[48],atte2[48],lambda1[48],lambda2[48];    Int_t ADCtof[48];
1098    Int_t temp=0;  
1099    for(Int_t i=0; i<48; i++){  
1100      fileTriggerCalib >> temp;    // ---- introduce scale factors to tune simul ADC to real data   24-oct DC
1101      fileTriggerCalib >> atte1[i];    //   Float_t ScaleFact[48]={0.18,0.22,0.35,0.26,0.47,0.35,0.31,0.37,
1102      fileTriggerCalib >> atte2[i];    //                     0.44,0.23,0.38,0.60,0.39,0.29,0.40,0.23,
1103      fileTriggerCalib >> lambda1[i];    //                     0.30,0.66,0.22,1.53,0.17,0.55,
1104      fileTriggerCalib >> lambda2[i];    //                     0.84,0.19,0.21,1.64,0.62,0.13,
1105      fileTriggerCalib >> temp;    //                     0.18,0.15,0.10,0.14,0.14,0.14,0.14,0.12,
1106    }    //                     0.26,0.18,0.25,0.23,0.20,0.40,0.19,0.23,0.25,0.23,0.25,0.20};
1107    fileTriggerCalib.close();    
1108      // new scale factors: WM 30-Oct-07
1109    //  Read from file the 48*4 values of the attenuation fit function    //   Float_t ScaleFact[48]={0.35,0.41,0.32,0.34,0.58,0.47,0.42,0.44,
1110    //  NB: lambda<0; x,y defined in gpamela (=0 in the centre of the cavity)    //                     0.50,0.34,0.50,0.50,0.51,0.42,0.46,0.25,
1111    //    Qhitpmt_pC =  atte1 * exp(x/lambda1) + atte2 * exp(x/lambda2)    //                     0.20,0.38,0.29,0.49,0.24,0.68,
1112      //                     0.30,0.26,0.28,0.79,0.31,0.12,
1113    //    fine lettura dal file */    //                     0.25,0.21,0.14,0.20,
1114      //                     0.16,0.17,0.19,0.18,
1115    //const Int_t nmax=??; = Nthtof    //                     0.34,0.27,0.34,0.31,0.25,0.57,
1116    Int_t ip, ipad;    //                     0.24,0.34,0.34,0.32,0.31,0.30};
1117    //Int_t ipmt;    
1118    Int_t pmtleft=0, pmtright=0;    Float_t ScaleFact[48]={0.39, 0.49, 0.38, 0.40, 0.65, 0.51, 0.43,
1119    Int_t *pl, *pr;                           0.49, 0.58, 0.38, 0.53, 0.57, 0.53, 0.45, 0.49, 0.16,
1120    pl = &pmtleft;                           0.15, 0.44, 0.28, 0.57, 0.26, 0.72, 0.37, 0.29, 0.30, 0.89,
1121    pr = &pmtright;                           0.37, 0.08,  0.27, 0.23, 0.12, 0.22, 0.15, 0.16, 0.21,
1122                             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 };
1123    // TDC variables:  
1124    Int_t TDClast=4095;      // no signal --> ADC ch=4095    for(Int_t i=0; i<48; i++){
1125    Int_t TDCint[48];      QevePmt_pC[i] = 0;
1126    Float_t  tdc[48],tdc1[48],tdcpmt[48];      ADCtof[i]=0;
1127    for(Int_t i=0; i<48; i++)    }
1128      tdcpmt[i] = 1000.;    
1129    Float_t thresh=1.; // to be defined better... (Wolfgang)    // ------ read calibration file (get A1, A2, lambda1, lambda2)
1130      ifstream fileTriggerCalib;
1131      // === TDC: simulate timing for each paddle    TString ftrigname="TrigCalibParam.txt";
1132      Float_t dt1 = 285.e-12 ;   // single PMT resolution    fileTriggerCalib.open(ftrigname.Data());
1133      Float_t tdcres[50],c1_S[50],c2_S[50],c3_S[50];    if ( !fileTriggerCalib ) {
1134      for(Int_t j=0;j<48;j++)  tdcres[j] = 50.E-12;   // TDC resolution 50 picosec      printf("debug: no trigger calib file!\n");
1135      for(Int_t j=0;j<48;j++)  c1_S[j] = 500.;  // cable length in channels      return(-117); //check output!
1136      for(Int_t j=0;j<48;j++)  c2_S[j] = 0.;    };
1137      for(Int_t j=0;j<48;j++)  c3_S[j] = 1000.;    Float_t atte1[48],atte2[48],lambda1[48],lambda2[48];
1138      for(Int_t j=0;j<48;j++)  c1_S[j] = c1_S[j]*tdcres[j];   // cable length in sec    Int_t temp=0;
1139      for(Int_t j=0;j<48;j++)  c2_S[j] = c2_S[j]*tdcres[j];    // correct readout WM Oct '07
1140       // ih = 0 + i1;  // not used?? (Silvio)    for(Int_t i=0; i<48; i++){
1141        fileTriggerCalib >> temp;
1142    /* **********************************  inizio loop sugli hit */      fileTriggerCalib >> atte1[i];
1143          fileTriggerCalib >> lambda1[i];
1144    for(Int_t nh=0; nh<Nthtof; nh++){      fileTriggerCalib >> atte2[i];
1145            fileTriggerCalib >> lambda2[i];
1146      for(Int_t j=0; j<2; j++) { // already done!! remove???      fileTriggerCalib >> temp;
1147        xpath[j]=0.;    }
1148        ypath[j]=0.;    fileTriggerCalib.close();
1149        FGeo[j]=0.;  
1150      }    Int_t ip, ipad;
1151      //Int_t ipmt;
1152      Float_t s_l_g[6] = {8.0, 8.0, 20.9, 22.0, 9.8, 8.3 };  // length of the lightguide    Int_t pmtleft=0, pmtright=0;
1153      Float_t t1,t2,veff,veff1,veff0 ;    Int_t *pl, *pr;
1154      veff0 = 100.*1.0e8 ; // light velocity in the scintillator in m/sec    pl = &pmtleft;
1155      veff1 = 100.*1.5e8; // light velocity in the lightguide in m/sec    pr = &pmtright;
1156      veff=veff0;         // signal velocity in the paddle  
1157      // TDC variables:
1158      t1 = Timetof[nh] ;  // Start    Int_t TDClast=4095;      // no signal --> TDC ch=4095
1159      t2 = Timetof[nh] ;    Int_t TDCint[48];
1160      Float_t  tdc[48],tdc1[48],tdcpmt[48];
1161      // Donatella        for(Int_t i=0; i<48; i++) {
1162      // ridefiniz. piano e pad per i vettori in C      tdcpmt[i] = 1000.;
1163      ip = Ipltof[nh]-1;      tdc[i]  = 0.;            // 18-oct WM
1164      ipad = Ipaddle[nh]-1;      tdc1[i] = 0.;            // 18-oct WM
1165      pmtleft=0;    }
1166      pmtright=0;  
1167        Float_t thresh=10.; // to be defined better... (Wolfgang)
1168      //Paddle2Pmt((Int_t)ip, (Int_t) ipad, (Int_t*) &pmtleft, (Int_t*) &pmtright);  
1169      Paddle2Pmt(ip, ipad, &pmtleft, &pmtright);    // === TDC: simulate timing for each paddle
1170      //Paddle2Pmt(ip, ipad, pl, pr);    //Float_t dt1 = 285.e-12 ;   // single PMT resolution
1171          Float_t dt1 = 425.e-12 ;   // single PMT resolution (WM, Nov'07)
1172      // per avere anche la corrispondenza pmt --> half board e canale    Float_t tdcres[50],c1_S[50],c2_S[50],c3_S[50];
1173      // metodo GetPMTIndex(Int_t ipmt, Int_t &hb, Int_t &ch) // non lo usiamo x ora    for(Int_t j=0;j<48;j++)  tdcres[j] = 50.E-12;   // TDC resolution 50 picosec
1174        for(Int_t j=0;j<48;j++)  c1_S[j] = 500.;  // cable length in channels
1175      /*calcola la pos media e il path all'interno della paddle */    for(Int_t j=0;j<48;j++)  c2_S[j] = 0.;
1176        for(Int_t j=0;j<48;j++)  c3_S[j] = 1000.;
1177      Float_t tpos=0.;    for(Int_t j=0;j<48;j++)  c1_S[j] = c1_S[j]*tdcres[j];   // cable length in sec
1178      Float_t path[2] = {0., 0.};    for(Int_t j=0;j<48;j++)  c2_S[j] = c2_S[j]*tdcres[j];
1179      //--- Strip in Y = S11,S22,S31 ------    // ih = 0 + i1;  // not used?? (Silvio)
1180      if(ip==0 || ip==3 || ip==4)  
1181        tpos = (Yintof[nh]+Youttof[nh])/2.;    /* **********************************  start loop over hits */
1182      else    
1183        if(ip==1 || ip==2 || ip==5)   //--- Strip in X per S12,S21,S32    for(Int_t nh=0; nh<Nthtof; nh++){
1184          tpos = (Xintof[nh]+Xouttof[nh])/2.;      
1185        else if (ip!=6)      Float_t s_l_g[6] = {8.0, 8.0, 20.9, 22.0, 9.8, 8.3 };  // length of the lightguide
1186          printf("*** Warning: this option should never occur! (ip=%2i, nh=%2i)\n",ip,nh);      Float_t t1,t2,veff,veff1,veff0 ;
1187      path[0]= tpos + dimel[ip]/2.;      veff0 = 100.*1.0e8 ; // light velocity in the scintillator in m/sec
1188      path[1]= dimel[ip]/2.- tpos;      veff1 = 100.*1.5e8; // light velocity in the lightguide in m/sec
1189        veff=veff0;         // signal velocity in the paddle
1190      //  cout <<"Strip N. ="<< ipaddle <<" piano n.= "<< iplane <<" POSIZ = "<< tpos <<"\n";  
1191        t1 = Timetof[nh] ;  // Start
1192      /* per il momento metto un fattore geometrico costante*/        t2 = Timetof[nh] ;
1193      FGeo[0] =0.5;  
1194      FGeo[1] =0.5;      // Donatella: redefinition plane and pad for vectors in C
1195      //  FGeo[1] = atan(path[1]/dimes[ip])/6.28318; // frazione fotoni verso SX      ip = Ipltof[nh]-1;
1196      //  FGeo[2] = atan(path[2]/dimes[ip])/6.28318; // e verso DX      ipad = Ipaddle[nh]-1;
1197              pmtleft=0;
1198      /*  rimando la fluttuazione poissoniana sui fotoni prodotti      pmtright=0;
1199          sto studiando come funziona la funzione:  
1200          long int i = sto.Poisson(double x);  */      // WM: S12 paddles are "reversed" (Nov'07)
1201      //  Npho = Poisson(ERELTOF[nh])*Pho_keV*1e6   Eloss in GeV ?      if (ip==2)
1202      Npho = Ereltof[nh]*Pho_keV*10.0e6;  // Eloss in GeV ?        if (ipad==0)
1203              ipad=1;
1204      Float_t knorm[2]={0., 0.}; // Donatella        else
1205      Float_t Atten[2]={0., 0.}; // Donatella          ipad=0;
1206      for(Int_t j=0; j<2; j++){  
1207        QhitPad_pC[j]= Npho*FGeo[j]*effi*pmGain*echarge;     if (ip<6) {
1208        knorm[j]=QhitPad_pC[j]/(atte1[pmtleft+j]*exp((dimel[ip]/2.*pow(-1,j+1))/lambda1[pmtleft+j]) +        Paddle2Pmt(ip, ipad, &pmtleft, &pmtright);
1209                                atte2[pmtleft+j]*exp((dimel[ip]/2.*pow(-1,j+1))/lambda2[pmtleft+j]));  
1210          // DC: evaluates mean position and path inside the paddle
1211        Atten[j]=knorm[j]*(atte1[pmtleft+j]*exp(tpos/lambda1[pmtleft+j]) +    
1212                           atte2[pmtleft+j]*exp(tpos/lambda2[pmtleft+j]));        Float_t tpos=0.;
1213                Float_t path[2] = {0., 0.};
1214        QhitPmt_pC[j]= QhitPad_pC[j]*Atten[j];        //--- Strip in Y = S11,S22,S31 ------
1215      }        if(ip==0 || ip==3 || ip==4)
1216                tpos = (Yintof[nh]+Youttof[nh])/2.;
1217      QevePmt_pC[pmtleft]  += QhitPmt_pC[0];        else
1218      QevePmt_pC[pmtright] += QhitPmt_pC[1];          if(ip==1 || ip==2 || ip==5)   //--- Strip in X for S12,S21,S32
1219                  tpos = (Xintof[nh]+Xouttof[nh])/2.;
1220      // TDC          else //if (ip!=6)
1221      t2 = t2 + fabs(path[0]/veff) + s_l_g[ip]/veff1 ;  // Signal reaches PMT            printf("*** WARNING TOF: this option should never occur! (ip=%2i, nh=%2i)\n",ip,nh);
1222      t1 = t1 + fabs(path[1]/veff) + s_l_g[ip]/veff1;  
1223              path[0]= tpos + dimel[ip]/2.;   // path to left PMT
1224      TRandom r;        path[1]= dimel[ip]/2.- tpos;    // path to right PMT
1225      t1 = r.Gaus(t1,dt1);  //apply gaussian error  dt  
1226      t2 = r.Gaus(t2,dt1);  //apply gaussian error  dt        //  cout <<"Strip N. ="<< ipaddle <<" piano n.= "<< iplane <<" POSIZ = "<< tpos <<"\n";
1227              
1228      t1 = t1 + c1_S[pmtleft] ;  // Signal reaches Discriminator ,TDC starts  to run        if (DEBUG) {
1229      t2 = t2 + c1_S[pmtright] ;          cout <<" plane "<<ip<<" strip # ="<< ipad <<" tpos  "<< tpos <<"\n";
1230                cout <<"pmtleft, pmtright "<<pmtleft<<" "<<pmtright<<endl;
1231      // check if signal is above threshold        }
1232      // then check if tdcpmt is already filled by another hit...        
1233      // only re-fill if time is smaller        // constant geometric factor, the rest is handled by the scaling factor
1234              FGeo[0] =0.5;
1235      if (QhitPmt_pC[0] > thresh)        FGeo[1] =0.5;
1236        if (tdcpmt[pmtleft] < 1000.) // is already filled!        //  FGeo[1] = atan(path[1]/dimes[ip])/6.28318; // fraction of photons toward left
1237          if (t1 <  tdcpmt[pmtleft]) {        //  FGeo[2] = atan(path[2]/dimes[ip])/6.28318; // toward right
1238            tdcpmt[pmtleft] = t1;        
1239            t1 = t1 + c2_S[pmtleft] ;  // Signal reaches Coincidence        
1240            tdc[pmtleft] = t1;        //  Npho = Poisson(ERELTOF[nh])*Pho_keV*1e6  Poissonian fluctuations to be inserted-DC
1241          }        Npho = Ereltof[nh]*Pho_keV*1.0e6;  // Eloss in GeV
1242              
1243      if (QhitPmt_pC[1] > thresh)        Float_t knorm[2]={0., 0.}; // Donatella
1244          if (tdcpmt[pmtright] < 1000.)  // is already filled!        Float_t Atten[2]={0., 0.}; // Donatella
1245            if (t2 <  tdcpmt[pmtright]) {        for(Int_t j=0; j<2; j++){
1246              tdcpmt[pmtright] = t2;          QhitPad_pC[j]= Npho*FGeo[j]*effi*pmGain*echarge*1.E12*ScaleFact[pmtleft+j];
1247              t2 = t2 + c2_S[pmtright] ;          // WM
1248              tdc[pmtright] = t2;          knorm[j]=atte1[pmtleft+j]*exp(lambda1[pmtleft+j]*dimel[ip]/2.*pow(-1,j+1)) +
1249            }                  atte2[pmtleft+j]*exp(lambda2[pmtleft+j]*dimel[ip]/2.*pow(-1,j+1));
1250                Atten[j]=atte1[pmtleft+j]*exp(tpos*lambda1[pmtleft+j]) +
1251    } // ****************************************       end loop over hits            atte2[pmtleft+j]*exp(tpos*lambda2[pmtleft+j]) ;
1252            QhitPmt_pC[j]= QhitPad_pC[j]*Atten[j]/knorm[j];
1253    // ======  ADC ======          //              QhitPmt_pC[j]= QhitPad_pC[j]; //no attenuation
1254    for(Int_t i=0; i<48; i++){  
1255      if(QevePmt_pC[i] != 0.){  
1256        ADCtof[i]= (Int_t)(ADC_pC*QevePmt_pC[i] + ADCoffset);          if (DEBUG) {
1257        if(ADCtof[i]> ADClast) ADCtof[i]=ADClast;            cout<<"pmtleft "<<pmtleft<<" j "<<j<<endl;
1258      } else            cout<<" atte1 "<<atte1[pmtleft+j]<<"lambda1 "<<lambda1[pmtleft+j]<<" atte2 "<<atte2[pmtleft+j]<<"lambda2 "<<lambda2[pmtleft+j] <<endl;    
1259        ADCtof[i]= ADClast;            cout<<j<<" tpos "<<tpos<<" knorm "<<knorm[j]<<" "<<Atten[j]<<" "<<"QhitPmt_pC "<<QhitPmt_pC[j]<<endl;    
1260    };          }
1261          }
1262              
1263    // ======  build  TDC coincidence  ======        if (DEBUG)
1264            cout<<"Npho "<<Npho<<" QhitPmt_pC "<<QhitPmt_pC[0]<<" "<<QhitPmt_pC[1]<<endl;  
1265    Float_t t_coinc = 0;        
1266    Int_t ilast = 100;        QevePmt_pC[pmtleft]  += QhitPmt_pC[0];
1267    for (Int_t ii=0; ii<48;ii++)        QevePmt_pC[pmtright] += QhitPmt_pC[1];
1268      if (tdc[ii] > t_coinc) {        
1269        t_coinc = tdc[ii];        // TDC
1270        ilast = ii;        // WM right and left <->
1271      }        //      t2 = t2 + fabs(path[0]/veff) + s_l_g[ip]/veff1 ;  // Signal reaches PMT
1272            //      t1 = t1 + fabs(path[1]/veff) + s_l_g[ip]/veff1;
1273    //     cout<<ilast<<" "<<t_coinc<<endl;  
1274    //     At t_coinc  trigger condition is fulfilled        t1 = t1 + fabs(path[0]/veff) + s_l_g[ip]/veff1;
1275            t2 = t2 + fabs(path[1]/veff) + s_l_g[ip]/veff1 ;  // Signal reaches PMT
1276    for (Int_t ii=0; ii<48;ii++){        
1277      //      if (tdc[ii] != 0) tdc1[ii] = t_coinc - tdc[ii];   // test 1        t1 = gRandom->Gaus(t1,dt1); //apply gaussian error  dt
1278      if (tdc[ii] != 0) tdc1[ii] = t_coinc - tdcpmt[ii];  // test 2        t2 = gRandom->Gaus(t2,dt1); //apply gaussian error  dt
1279      tdc1[ii] = tdc1[ii]/tdcres[ii];                     // divide by TDC resolution  
1280      if (tdc[ii] != 0) tdc1[ii] = tdc1[ii] + c3_S[ii];  // add cable length c3        t1 = t1 + c1_S[pmtleft] ;  // Signal reaches Discriminator ,TDC starts  to run
1281          t2 = t2 + c1_S[pmtright] ;
1282    } // missing parenthesis inserted! (Silvio)        
1283          // check if signal is above threshold
1284    for(Int_t i=0; i<48; i++){        // then check if tdcpmt is already filled by another hit...
1285      if(tdc1[i] != 0.){        // only re-fill if time is smaller
1286        TDCint[i]=(Int_t)tdc1[i];        
1287        //ADC[i]= ADC_pC * QevePmt_pC[i] + ADCoffset;        if (QhitPmt_pC[0] > thresh) {
1288        //if(ADC[i]> ADClast) ADC[i]=ADClast;          if (tdcpmt[pmtleft] == 1000.) {  // fill for the first time
1289      } else            tdcpmt[pmtleft] = t1;
1290        TDCint[i]= TDClast;            tdc[pmtleft] = t1 + c2_S[pmtleft] ;  // Signal reaches Coincidence
1291    }          }
1292            if (tdcpmt[pmtleft] < 1000.) // is already filled!
1293  // ======  write fDataTof  =======            if (t1 <  tdcpmt[pmtleft]) {
1294    UChar_t tofBin;              tdcpmt[pmtleft] = t1;
1295    for (Int_t j=0; j < 12; j++){              t1 = t1 + c2_S[pmtleft] ;  // Signal reaches Coincidence
1296      Int_t j12=j*12;              tdc[pmtleft] = t1;
1297      fDataTof[j12+0]=0x00;   // TDC_ID            }
1298      fDataTof[j12+1]=0x00;   // EV_COUNT        }      
1299      fDataTof[j12+2]=0x00;   // TDC_MASK (1)        if (QhitPmt_pC[1] > thresh) {
1300      fDataTof[j12+3]=0x00;   // TDC_MASK (2)          if (tdcpmt[pmtright] == 1000.) {  // fill for the first time
1301      for (Int_t k=0; k < 4; k++){            tdcpmt[pmtright] = t2;
1302        Int_t jk12=j12+k;            tdc[pmtright] = t2 + c2_S[pmtright] ;  // Signal reaches Coincidence
1303        tofBin=(UChar_t)(ADCtof[k+4*j]/256);   // ADC# (msb) (#=1+k+4*j)          }
1304        fDataTof[jk12+4] = Bin2GrayTof(tofBin,fDataTof[jk12+4]);          if (tdcpmt[pmtright] < 1000.)  // is already filled!
1305        tofBin=(UChar_t)(ADCtof[k+4*j]%256);   // ADC# (lsb)            if (t2 <  tdcpmt[pmtright]) {
1306        fDataTof[jk12+5] = Bin2GrayTof(tofBin,fDataTof[jk12+5]);              tdcpmt[pmtright] = t2;
1307        tofBin=(UChar_t)(TDCint[k+4*j]/256);   // TDC# (msb)              t2 = t2 + c2_S[pmtright] ;
1308        fDataTof[jk12+6]=Bin2GrayTof(tofBin,fDataTof[jk12+6]);              tdc[pmtright] = t2;
1309        tofBin=(UChar_t)(TDCint[k+4*j]%256);   // TDC# (lsb)            }      
1310        fDataTof[jk12+7]=Bin2GrayTof(tofBin,fDataTof[jk12+7]);        }
1311      };  
1312      fDataTof[j12+20]=0x00;   // TEMP1        if (DEBUG)
1313      fDataTof[j12+21]=0x00;   // TEMP2          cout<<nh<<" "<<Timetof[nh]<<" "<<t1<<" "<<t2<<endl;
1314      fDataTof[j12+22]= EvaluateCrcTof(pTof);   // CRC  
1315      pTof+=23;      } // ip < 6
1316    };  
1317    return(0);    }; // ****************************************       end loop over hits
1318  };    
1319      // ======  ADC ======
1320  UChar_t Digitizer::Bin2GrayTof(UChar_t binaTOF,UChar_t grayTOF){  
1321    union graytof_data {    for(Int_t i=0; i<48; i++){
1322      UChar_t word;      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));
1323      struct bit_field {      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));
1324        unsigned b0:1;      if (QevePmt_pC[i] > 2485.) ADCtof[i]= (Int_t)(1758. + 0.54*QevePmt_pC[i]);  //assuming a fictional 0.54 ch/pC above ADCsat
1325        unsigned b1:1;      if (ADCtof[i]>ADCsat) ADCtof[i]=ADCsat;
1326        unsigned b2:1;      if (QevePmt_pC[i] < pCthres)  ADCtof[i]= ADClast;
1327        unsigned b3:1;      if (ADCtof[i] < 0) ADCtof[i]=ADClast;    
1328        unsigned b4:1;      if (ADCtof[i] > ADClast) ADCtof[i]=ADClast;
1329        unsigned b5:1;    }
1330        unsigned b6:1;  
1331        unsigned b7:1;  //   for(Int_t i=0; i<48; i++){
1332      } bit;  //     if(QevePmt_pC[i] >= pCthres){
1333    } bi,gr;  //       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));
1334    //  //   } else
1335    bi.word = binaTOF;  //       ADCtof[i]= ADClast;  
1336    gr.word = grayTOF;  // }
1337    //    
1338    gr.bit.b0 = bi.bit.b1 ^ bi.bit.b0;  // // ---- introduce scale factors to tune simul ADC to real data   24-oct DC
1339    gr.bit.b1 = bi.bit.b2 ^ bi.bit.b1;  
1340    gr.bit.b2 = bi.bit.b3 ^ bi.bit.b2;  //   for(Int_t i=0; i<48; i++){
1341    gr.bit.b3 = bi.bit.b3;  //     if(ADCtof[i] != ADClast){
1342    //  //                              //      printf("%3d, %4d, %4.2f\n",i, ADCtof[i],ScaleFact[i]);
1343    /* bin to gray conversion 4 bit per time*/  //       ADCtof[i]= Int_t (ADCtof[i]*ScaleFact[i]);
1344    //  //                              //      printf("%3d, %4d,\n",i, ADCtof[i]);
1345    gr.bit.b4 = bi.bit.b5 ^ bi.bit.b4;  //     }
1346    gr.bit.b5 = bi.bit.b6 ^ bi.bit.b5;  //   }
1347    gr.bit.b6 = bi.bit.b7 ^ bi.bit.b6;  
1348    gr.bit.b7 = bi.bit.b7;  //   for(Int_t i=0; i<48; i++){
1349    //  //     if(ADCtof[i] != ADClast){
1350    return(gr.word);  //       if(ADCtof[i]> ADCsat) ADCtof[i]=ADCsat;
1351  }  //       else if(ADCtof[i]< 0) ADCtof[i]=ADClast;    
1352    //   }
1353  UChar_t Digitizer::EvaluateCrcTof(UChar_t *pTof) {  //   }    
1354    // UChar_t crcTof=0x00;  
1355    //   for (Int_t jp=0; jp < 23; jp++){  
1356    //     crcTof = crc8(...)  // ======  build  TDC coincidence  ======
1357    //   }  
1358    return(0x00);    Float_t t_coinc = 0;
1359  };    Int_t ilast = 100;
1360      for (Int_t ii=0; ii<48;ii++)
1361  //void Digitizer::Paddle2Pmt(Int_t plane, Int_t paddle, Int_t* &pmtleft, Int_t* &pmtright){      if (tdc[ii] > t_coinc) {
1362  void Digitizer::Paddle2Pmt(Int_t plane, Int_t paddle, Int_t *pl, Int_t *pr){        t_coinc = tdc[ii];
1363    //* @param plane    (0 - 5)        ilast = ii;
1364    //* @param paddle   (plane=0, paddle = 0,...5)      }
1365    //* @param padid    (0 - 23)      
1366    //    //     cout<<ilast<<" "<<t_coinc<<endl;
1367    Int_t padid=-1;    //     At t_coinc  trigger condition is fulfilled
1368    Int_t pads[6]={8,6,2,2,3,3};    
1369    //    for (Int_t ii=0; ii<48;ii++){
1370    Int_t somma=0;      //      if (tdc[ii] != 0) tdc1[ii] = t_coinc - tdc[ii];   // test 1
1371    Int_t np=plane;      if (tdc[ii] != 0) tdc1[ii] = t_coinc - tdcpmt[ii];  // test 2
1372    for(Int_t j=0; j<np; j++)      tdc1[ii] = tdc1[ii]/tdcres[ii];                     // divide by TDC resolution
1373      somma+=pads[j];      if (tdc[ii] != 0) tdc1[ii] = tdc1[ii] + c3_S[ii];  // add cable length c3
1374    padid=paddle+somma;    } // missing parenthesis inserted! (Silvio)
1375    *pl = padid*2;  
1376    *pr = *pr + 1;    for(Int_t i=0; i<48; i++){
1377  };      if(tdc1[i] != 0.){
1378          TDCint[i]=(Int_t)tdc1[i];
1379  void Digitizer::DigitizeAC() {        if (TDCint[i]>4093) TDCint[i]=TDClast;  // 18-oct WM
1380    // created:  J. Conrad, KTH        if (DEBUG)
1381    // modified: S. Orsi, INFN Roma2          cout<<i<<" "<<TDCint[i]<<endl;
1382          //ADC[i]= ADC_pC * QevePmt_pC[i] + ADCoffset;
1383    fDataAC[0] = 0xACAC;        //if(ADC[i]> ADClast) ADC[i]=ADClast;
1384    fDataAC[64]= 0xACAC;      } else
1385    fDataAC[1] = 0xAC11;   // main card        TDCint[i]= TDClast;
1386    fDataAC[65] = 0xAC22;   // extra card    }
1387    
1388    // the third word is a status word (dummy)    if (DEBUG)
1389    fDataAC[2] = 0xFFFF; //FFEF?      cout<<"-----------"<<endl;
1390    fDataAC[66] = 0xFFFF;  
1391    
1392    const UInt_t nReg = 6;  //------ use channelmap  18-oct WM
1393    
1394    // Registers (dummy)    Int_t  channelmap[] =  {3,21,11,29,19,45,27,37,36,28,44,20,5,12,13,4,
1395    for (UInt_t i=0; i<=nReg; i++){                          6,47,14,39,22,31,30,23,38,15,46,7,0,33,16,24,
1396      fDataAC[i+4] = 0xFFFF;                          8,41,32,40,25,17,34,9,42,1,2,10,18,26,35,43};
1397      fDataAC[i+68] = 0xFFFF;  
1398    }    Int_t ADChelp[48];
1399      Int_t TDChelp[48];
1400    // the last word is a CRC  
1401    // Dummy for the time being, but it might need to be calculated in the end    for(Int_t i=0; i<48; i++){
1402    fDataAC[63] = 0xABCD;      Int_t ii=channelmap[i];
1403    fDataAC[127] = 0xABCD;      ADChelp[ii]= ADCtof[i];
1404        TDChelp[ii]= TDCint[i];
1405    // shift registers, which one is with respect to PMT, where in    }
1406    // shift registers is a question of time relative trigger  
1407    // In level2: hitmap, hitmap-status (synchronised with a trigger),    for(Int_t i=0; i<48; i++){
1408    // status      ADCtof[i]= ADChelp[i];
1409                  TDCint[i]= TDChelp[i];
1410    for (UInt_t i=0; i<=15; i++){    }
1411      fDataAC[i+11] = 0x0000;    
1412      fDataAC[i+75] = 0x0000;  
1413    }  // ======  write fDataTof  =======
1414    
1415    // singles counters are dummy  
1416    //  UChar_t tdcadd[8]={1,0,3,2,5,4,7,6};  (coded in 3 bit)
1417    for (UInt_t i=0; i<=16; i++){    UChar_t Ctrl3bit[8]={32,0,96,64,160,128,224,192};  // DC (msb in 8 bit word )
1418      fDataAC[i+26] = 0x0000;      
1419      fDataAC[i+90] = 0x0000;    UChar_t tofBin;
1420    }    for (Int_t j=0; j < 12; j++){   // loop on TDC #12
1421          Int_t j12=j*23;               // for each TDC 23 bytes (8 bits)
1422    // coincidences are dummy      fDataTof[j12+0]=0x00;   // TDC_ID
1423        fDataTof[j12+1]=0x00;   // EV_COUNT
1424    for (UInt_t i=0; i<=7; i++){      fDataTof[j12+2]=0x00;   // TDC_MASK (1)
1425      fDataAC[i+42] = 0x0000;      fDataTof[j12+3]=0x00;   // TDC_MASK (2)
1426      fDataAC[i+106] = 0x0000;      for (Int_t k=0; k < 4; k++){   // for each TDC 4 channels (ADC+TDC)
1427    }            
1428          Int_t jk12=j12+4*k;         // ADC,TDC channel (0-47)
1429    // increments for every trigger might be needed at some point.  
1430    // dummy for now        tofBin =(UChar_t)(ADCtof[k+4*j]/256);   // ADC# (msb)
1431    fDataAC[50]  = 0x0000;        fDataTof[jk12+4] = Bin2GrayTof(tofBin,fDataTof[jk12+4]);
1432    fDataAC[114] = 0x0000;        /* control bits inserted here, after the bin to gray conv - DC*/
1433          fDataTof[jk12+4] = Ctrl3bit[2*k] | fDataTof[jk12+4];
1434    // dummy FPGA clock        tofBin=(UChar_t)(ADCtof[k+4*j]%256);   // ADC# (lsb)
1435          fDataTof[jk12+5] = Bin2GrayTof(tofBin,fDataTof[jk12+5]);
1436    fDataAC[51] = 0x006C;        tofBin=(UChar_t)(TDCint[k+4*j]/256);   // TDC# (msb)
1437    fDataAC[52] = 0x6C6C;        fDataTof[jk12+6]=Bin2GrayTof(tofBin,fDataTof[jk12+6]);
1438    fDataAC[115] = 0x006C;        /* control bits inserted here, after the bin to gray conv - DC*/
1439    fDataAC[116] = 0x6C6C;        fDataTof[jk12+6] = Ctrl3bit[2*k+1] | fDataTof[jk12+6];
1440          tofBin=(UChar_t)(TDCint[k+4*j]%256);   // TDC# (lsb)
1441          fDataTof[jk12+7]=Bin2GrayTof(tofBin,fDataTof[jk12+7]);
1442    // dummy temperatures      };
1443    fDataAC[53] = 0x0000;      fDataTof[j12+20]=0x00;   // TEMP1
1444    fDataAC[54] = 0x0000;      fDataTof[j12+21]=0x00;   // TEMP2
1445    fDataAC[117] = 0x0000;      fDataTof[j12+22]= EvaluateCrcTof(pTof);   // CRC
1446    fDataAC[118] = 0x0000;      pTof+=23;
1447      };
1448      return(0);
1449    // dummy DAC thresholds  };
1450    for (UInt_t i=0; i<=7; i++){  
1451      fDataAC[i+55] = 0x1A13;    
1452      fDataAC[i+119] = 0x1A13;  UChar_t Digitizer::Bin2GrayTof(UChar_t binaTOF,UChar_t grayTOF){
1453    }    union graytof_data {
1454          UChar_t word;
1455    // We activate all branches. Once the digitization algorithm      struct bit_field {
1456    // is determined only the branches need to activated which involve needed        unsigned b0:1;
1457    // information        unsigned b1:1;
1458            unsigned b2:1;
1459    fhBookTree->SetBranchStatus("Nthcat",1);        unsigned b3:1;
1460    fhBookTree->SetBranchStatus("Iparcat",1);        unsigned b4:1;
1461    fhBookTree->SetBranchStatus("Icat",1);        unsigned b5:1;
1462    fhBookTree->SetBranchStatus("Xincat",1);        unsigned b6:1;
1463    fhBookTree->SetBranchStatus("Yincat",1);        unsigned b7:1;
1464    fhBookTree->SetBranchStatus("Zincat",1);      } bit;
1465    fhBookTree->SetBranchStatus("Xoutcat",1);    } bi,gr;
1466    fhBookTree->SetBranchStatus("Youtcat",1);    //
1467    fhBookTree->SetBranchStatus("Zoutcat",1);    bi.word = binaTOF;
1468    fhBookTree->SetBranchStatus("Erelcat",1);    gr.word = grayTOF;
1469    fhBookTree->SetBranchStatus("Timecat",1);    //
1470    fhBookTree->SetBranchStatus("Pathcat",1);    gr.bit.b0 = bi.bit.b1 ^ bi.bit.b0;
1471    fhBookTree->SetBranchStatus("P0cat",1);    gr.bit.b1 = bi.bit.b2 ^ bi.bit.b1;
1472    fhBookTree->SetBranchStatus("Nthcas",1);    gr.bit.b2 = bi.bit.b3 ^ bi.bit.b2;
1473    fhBookTree->SetBranchStatus("Iparcas",1);    gr.bit.b3 = bi.bit.b3;
1474    fhBookTree->SetBranchStatus("Icas",1);    //
1475    fhBookTree->SetBranchStatus("Xincas",1);    /* bin to gray conversion 4 bit per time*/
1476    fhBookTree->SetBranchStatus("Yincas",1);    //
1477    fhBookTree->SetBranchStatus("Zincas",1);    gr.bit.b4 = bi.bit.b5 ^ bi.bit.b4;
1478    fhBookTree->SetBranchStatus("Xoutcas",1);    gr.bit.b5 = bi.bit.b6 ^ bi.bit.b5;
1479    fhBookTree->SetBranchStatus("Youtcas",1);    gr.bit.b6 = bi.bit.b7 ^ bi.bit.b6;
1480    fhBookTree->SetBranchStatus("Zoutcas",1);    gr.bit.b7 = bi.bit.b7;
1481    fhBookTree->SetBranchStatus("Erelcas",1);    //
1482    fhBookTree->SetBranchStatus("Timecas",1);    return(gr.word);
1483    fhBookTree->SetBranchStatus("Pathcas",1);  }
1484    fhBookTree->SetBranchStatus("P0cas",1);  
1485    fhBookTree->SetBranchStatus("Nthcard",1);  UChar_t Digitizer::EvaluateCrcTof(UChar_t *pTof) {
1486    fhBookTree->SetBranchStatus("Iparcard",1);    Bool_t DEBUG=false;
1487    fhBookTree->SetBranchStatus("Icard",1);    if (DEBUG)
1488    fhBookTree->SetBranchStatus("Xincard",1);      return(0x00);
1489    fhBookTree->SetBranchStatus("Yincard",1);  
1490    fhBookTree->SetBranchStatus("Zincard",1);    UChar_t crcTof=0x00;
1491    fhBookTree->SetBranchStatus("Xoutcard",1);    UChar_t *pc=&crcTof, *pc2;
1492    fhBookTree->SetBranchStatus("Youtcard",1);    pc2=pTof;
1493    fhBookTree->SetBranchStatus("Zoutcard",1);    for (Int_t jp=0; jp < 23; jp++){
1494    fhBookTree->SetBranchStatus("Erelcard",1);      //crcTof = crc8(...)
1495    fhBookTree->SetBranchStatus("Timecard",1);      Crc8Tof(pc2++,pc);
1496    fhBookTree->SetBranchStatus("Pathcard",1);      //    printf("%2i --- %x\n",jp,crcTof);
1497    fhBookTree->SetBranchStatus("P0card",1);    }
1498      return(crcTof);
1499    // In this simpliefied approach we will assume that once  }
1500    // a particle releases > 0.5 mip in one of the 12 AC detectors it  
1501    // will fire. We will furthermore assume that both cards read out  void Digitizer::Crc8Tof(UChar_t *oldCRC, UChar_t *crcTof){
1502    // identical data.    union crctof_data {
1503        UChar_t word;
1504    // If you develop you digitization algorithm, you should start by      struct bit_field {
1505    // identifying the information present in level2 (post-darth-vader)        unsigned b0:1;
1506    // data.        unsigned b1:1;
1507          unsigned b2:1;
1508    Float_t SumEcat[5];        unsigned b3:1;
1509    Float_t SumEcas[5];        unsigned b4:1;
1510    Float_t SumEcard[5];        unsigned b5:1;
1511    for (Int_t k= 0;k<5;k++){        unsigned b6:1;
1512      SumEcat[k]=0.;        unsigned b7:1;
1513      SumEcas[k]=0.;      } bit;
1514      SumEcard[k]=0.;    } c,d,r;
1515    };  
1516      c.word = *oldCRC;
1517    if (Nthcat>50 || Nthcas>50 || Nthcard>50)    //d.word = *newCRC;
1518      printf("Error! NthAC out of range!\n\n");    d.word = *crcTof;
1519      r.word = 0;
1520    // look in CAT  
1521    //  for (UInt_t k= 0;k<50;k++){    r.bit.b0 = c.bit.b7 ^ c.bit.b6 ^ c.bit.b0 ^
1522    for (Int_t k= 0;k<Nthcat;k++){               d.bit.b0 ^ d.bit.b6 ^ d.bit.b7;
1523      if (Erelcat[k] > 0)  
1524        SumEcat[Icat[k]] += Erelcat[k];    r.bit.b1 = c.bit.b6 ^ c.bit.b1 ^ c.bit.b0 ^
1525    };               d.bit.b0 ^ d.bit.b1 ^ d.bit.b6;
1526    
1527    // look in CAS    r.bit.b2 = c.bit.b6 ^ c.bit.b2 ^ c.bit.b1 ^ c.bit.b0 ^
1528    for (Int_t k= 0;k<Nthcas;k++){               d.bit.b0 ^ d.bit.b1 ^ d.bit.b2 ^ d.bit.b6;
1529      if (Erelcas[k] >0)  
1530        SumEcas[Icas[k]] += Erelcas[k];    r.bit.b3 = c.bit.b7 ^ c.bit.b3 ^ c.bit.b2 ^ c.bit.b1 ^
1531    };               d.bit.b1 ^ d.bit.b2 ^ d.bit.b3 ^ d.bit.b7;
1532    
1533    // look in CARD    r.bit.b4 = c.bit.b4 ^ c.bit.b3 ^ c.bit.b2 ^
1534    for (Int_t k= 0;k<Nthcard;k++){               d.bit.b2 ^ d.bit.b3 ^ d.bit.b4;
1535      if (Erelcard[k] >0)  
1536        SumEcard[Icard[k]] += Erelcard[k];    r.bit.b5 = c.bit.b5 ^ c.bit.b4 ^ c.bit.b3 ^
1537    };               d.bit.b3 ^ d.bit.b4 ^ d.bit.b5;
1538    
1539    // channel mapping              Hit Map    r.bit.b6 = c.bit.b6 ^ c.bit.b5 ^ c.bit.b4 ^
1540    // 1 CARD4                      0          LSB               d.bit.b4 ^ d.bit.b5 ^ d.bit.b6;
1541    // 2 CAT2                       0  
1542    // 3 CAS1                       0    r.bit.b7 = c.bit.b7 ^ c.bit.b6 ^ c.bit.b5 ^
1543    // 4 NC                         0               d.bit.b5 ^ d.bit.b6 ^ d.bit.b7 ;
1544    // 5 CARD2                      0  
1545    // 6 CAT4                       1    *crcTof=r.word;
1546    // 7 CAS4                       0      //return r.word;
1547    // 8 NC                         0  };
1548    // 9 CARD3                      0  
1549    // 10 CAT3                      0  //void Digitizer::Paddle2Pmt(Int_t plane, Int_t paddle, Int_t* &pmtleft, Int_t* &pmtright){
1550    // 11 CAS3                      0  void Digitizer::Paddle2Pmt(Int_t plane, Int_t paddle, Int_t *pl, Int_t *pr){
1551    // 12 NC                        0    //* @param plane    (0 - 5)
1552    // 13 CARD1                     0    //* @param paddle   (plane=0, paddle = 0,...5)
1553    // 14 CAT1                      0    //* @param padid    (0 - 23)  
1554    // 15 CAS2                      0    //
1555    // 16 NC                        0          MSB    Int_t padid=-1;
1556      Int_t pads[6]={8,6,2,2,3,3};
1557    // In the first version only the hit-map is filled, not the SR.    //
1558      Int_t somma=0;
1559    // Threshold: 0.8 MeV.    Int_t np=plane;
1560      for(Int_t j=0; j<np; j++)
1561    Float_t thr = 8e-4;      somma+=pads[j];
1562      padid=paddle+somma;
1563    fDataAC[3] = 0x0000;    *pl = padid*2;
1564      //  *pr = *pr + 1;
1565    if (SumEcas[0] > thr)  fDataAC[3]  = 0x0004;    *pr = *pl + 1; // WM
1566    if (SumEcas[1] > thr)  fDataAC[3] += 0x4000;  };
1567    if (SumEcas[2] > thr)  fDataAC[3] += 0x0400;  
1568    if (SumEcas[3] > thr)  fDataAC[3] += 0x0040;    void Digitizer::DigitizeAC() {
1569      // created:  J. Conrad, KTH
1570    if (SumEcat[0] > thr)  fDataAC[3] += 0x2000;    // modified: S. Orsi, INFN Roma2
1571    if (SumEcat[1] > thr)  fDataAC[3] += 0x0002;    // fDataAC[0-63]:   main AC board
1572    if (SumEcat[2] > thr)  fDataAC[3] += 0x0200;    // fDataAC[64-127]: extra AC board
1573    if (SumEcat[3] > thr)  fDataAC[3] += 0x0020;  
1574        fDataAC[0] = 0xACAC;
1575    if (SumEcard[0] > thr)  fDataAC[3] += 0x1000;    fDataAC[64]= 0xACAC;
1576    if (SumEcard[1] > thr)  fDataAC[3] += 0x0010;    fDataAC[1] = 0xAC11;  
1577    if (SumEcard[2] > thr)  fDataAC[3] += 0x0100;    fDataAC[65] = 0xAC22;  
1578    if (SumEcard[3] > thr)  fDataAC[3] += 0x0001;  
1579      // the third word is a status word (dummy: "no errors are present in the AC boards")
1580    fDataAC[67] = fDataAC[3];    fDataAC[2] = 0xFFFF; //FFEF?
1581      fDataAC[66] = 0xFFFF;
1582    //    for (Int_t i=0; i<fACbuffer; i++){  
1583    //    printf("%0x  ",fDataAC[i]);      const UInt_t nReg = 6;
1584    //    if ((i+1)%8 ==0) cout << endl;  
1585    //   }    // FPGA Registers (dummy)
1586  };    for (UInt_t i=0; i<=nReg; i++){
1587        fDataAC[i+4] = 0xFFFF;
1588        fDataAC[i+68] = 0xFFFF;
1589      }
1590  void Digitizer::DigitizeS4(){  
1591    // creato:  S. Borisov, INFN Roma2 e MEPHI, Sett 2007    // the last word is a CRC
1592    TString ciao,modo="ns";    // Dummy for the time being, but it might need to be calculated in the end
1593    Int_t i,j,t,NdF,pmt,NdFT,S4,S4v=0,S4p=32;    fDataAC[63] = 0xABCD;
1594    Float_t E0,E1=1e-6,Ert,X,Y,Z,x,y,z,V[3],Xs[2],Ys[2],Zs[2],Yp[6],q,w,p=0.1,l,l0=500.;    fDataAC[127] = 0xABCD;
1595    Xs[0]=-24.1;  
1596    Xs[1]=24.1;    // shift registers (moved to the end of the routine)
1597    Ys[0]=-24.1;  
1598    Ys[1]=24.1;    Int_t evntLSB=Ievnt%65536;
1599    Zs[0]=-0.5;    Int_t evntMSB=(Int_t)(Ievnt/65536);
1600    Zs[1]=0.5;  
1601    Yp[0]=-20.;    // singles counters are dummy
1602    Yp[2]=-1.;    for (UInt_t i=0; i<=15; i++){  //SO Oct '07: // for (UInt_t i=0; i<=16; i++){
1603    Yp[4]=17.;      //     fDataAC[i+26] = 0x0000;  
1604    for(i=0;i<3;i++)      //     fDataAC[i+90] = 0x0000;
1605      Yp[2*i+1]=Yp[2*i]+3;      fDataAC[i+26] = evntLSB;
1606    srand(time(NULL));      fDataAC[i+90] = evntLSB;
1607    // --- activate branches:    };
1608    fhBookTree->SetBranchStatus("Nthtof",1);    
1609    fhBookTree->SetBranchStatus("Ipltof",1);    // coincidences are dummy (increment by 1 at each event)
1610    fhBookTree->SetBranchStatus("Ipaddle",1);    // for (UInt_t i=0; i<=7; i++){
1611        //    fDataAC[i+42] = 0x0000;
1612    fhBookTree->SetBranchStatus("Xintof",1);    //    fDataAC[i+106] = 0x0000;
1613    fhBookTree->SetBranchStatus("Yintof",1);    //   }
1614    fhBookTree->SetBranchStatus("Xouttof",1);    for (UInt_t i=0; i<=7; i++){
1615    fhBookTree->SetBranchStatus("Youttof",1);      fDataAC[i+42] = evntLSB;
1616          fDataAC[i+106] = evntLSB;
1617    fhBookTree->SetBranchStatus("Ereltof",1);    };
1618    fhBookTree->SetBranchStatus("Timetof",1);  
1619    NdFT=0;    // increments for every trigger might be needed at some point.
1620    Ert=0;    // dummy for now
1621    for(i=0;i<Nthtof;i++){    fDataAC[50]  = 0x0000;
1622      if(Ipltof[i]!=6) continue;    fDataAC[114] = 0x0000;
1623      Ert+=Ereltof[i];  
1624            // dummy FPGA clock (increment by 1 at each event)
1625          /*    
1626      if(modo=="ns") continue;      fDataAC[51] = 0x006C;
1627      NdF=Int_t(Ereltof[i]/E1);      fDataAC[52] = 0x6C6C;
1628      NdFT=0;      fDataAC[115] = 0x006C;
1629      X=Xintof[i];      fDataAC[116] = 0x6C6C;
1630      Y=Yintof[i];    */
1631      Z=((Float_t)random()/(Float_t)0x7fffffff)-0.5;    if (Ievnt<=0xFFFF) {
1632      //cout<<"XYZ "<<X<<"  "<<Y<<"  "<<Z<<endl;      fDataAC[51] = 0x0000;
1633      for(j=0;j<NdF;j++){      fDataAC[52] = Ievnt;
1634        q=(Float_t)random()/(Float_t)0x7fffffff;      fDataAC[115] = 0x0000;
1635        w=(Float_t)random()/(Float_t)0x7fffffff;      fDataAC[116] = Ievnt;
1636        // cout<<"qw "<<q<<" "<<w<<endl;    } else {
1637        V[0]=p*cos(6.28318*q);      fDataAC[51] = evntMSB;
1638        V[1]=p*sin(6.28318*q);      fDataAC[52] = evntLSB;
1639        V[2]=p*(2.*w-1.);      fDataAC[115] = fDataAC[51];
1640        pmt=0;      fDataAC[116] = fDataAC[52];
1641        x=X;    }
1642        y=Y;  
1643        z=Z;    // dummy temperatures
1644        while(pmt==0 && (x>Xs[0] && x<Xs[1])&&(y>Ys[0] && y<Ys[1])&&(z>Zs[0] && z<Zs[1])){    fDataAC[53] = 0x0000;
1645          l=0;    fDataAC[54] = 0x0000;
1646          while(pmt==0 && (x>Xs[0] && x<Xs[1])&&(y>Ys[0] && y<Ys[1])&&(z>Zs[0] && z<Zs[1])){    fDataAC[117] = 0x0000;
1647            x+=V[0];    fDataAC[118] = 0x0000;
1648            y+=V[1];  
1649            z+=V[2];  
1650            l+=p;    // dummy DAC thresholds
1651            //cout<<x<<"  "<<y<<"  "<<z<<"  "<<l<<endl;    for (UInt_t i=0; i<=7; i++){
1652            //cin>>ciao;      fDataAC[i+55] = 0x1A13;  
1653          }      fDataAC[i+119] = 0x1A13;
1654          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)){    }
1655            for(t=0;t<3;t++){    
1656              if(y>=Yp[2*t] && y<Yp[2*t+1]){    // We activate all branches. Once the digitization algorithm is determined
1657                if(pmt==0)NdFT++;    // only the branches that involve needed information will be activated
1658                pmt=1;    
1659                //cout<<NdFT<<endl;    fhBookTree->SetBranchAddress("Ievnt",&Ievnt);
1660                break;    fhBookTree->SetBranchStatus("Nthcat",1);
1661              }    fhBookTree->SetBranchStatus("Iparcat",1);
1662            }    fhBookTree->SetBranchStatus("Icat",1);
1663            if(pmt==1)break;    fhBookTree->SetBranchStatus("Xincat",1);
1664            V[0]=-V[0];    fhBookTree->SetBranchStatus("Yincat",1);
1665          }    fhBookTree->SetBranchStatus("Zincat",1);
1666          q=(Float_t)random()/(Float_t)0x7fffffff;    fhBookTree->SetBranchStatus("Xoutcat",1);
1667          w=1-exp(-l/l0);    fhBookTree->SetBranchStatus("Youtcat",1);
1668          if(q<w)break;    fhBookTree->SetBranchStatus("Zoutcat",1);
1669          q=(Float_t)random()/(Float_t)0x7fffffff;    fhBookTree->SetBranchStatus("Erelcat",1);
1670          w=0.5;    fhBookTree->SetBranchStatus("Timecat",1);
1671          if(q<w)break;    fhBookTree->SetBranchStatus("Pathcat",1);
1672          if((x>Xs[0]+p && x<Xs[1]-p)&&(y<Ys[0]+p || y>Ys[1]-p)&&(z>Zs[0]+p && z<Zs[1]-p))V[1]=-V[1];    fhBookTree->SetBranchStatus("P0cat",1);
1673          if((x>Xs[0]+p && x<Xs[1]-p)&&(y>Ys[0]+p && y<Ys[1]-p)&&(z<Zs[0]+p || z>Zs[1]-p))V[2]=-V[2];    fhBookTree->SetBranchStatus("Nthcas",1);
1674          x+=V[0];    fhBookTree->SetBranchStatus("Iparcas",1);
1675          y+=V[1];    fhBookTree->SetBranchStatus("Icas",1);
1676          z+=V[2];    fhBookTree->SetBranchStatus("Xincas",1);
1677          l=0;    fhBookTree->SetBranchStatus("Yincas",1);
1678          //cout<<x<<"  "<<y<<"  "<<z<<"  "<<l<<endl;    fhBookTree->SetBranchStatus("Zincas",1);
1679                  //cin>>ciao;    fhBookTree->SetBranchStatus("Xoutcas",1);
1680        }    fhBookTree->SetBranchStatus("Youtcas",1);
1681      }    fhBookTree->SetBranchStatus("Zoutcas",1);
1682    }    fhBookTree->SetBranchStatus("Erelcas",1);
1683    Ert=Ert/0.002;    fhBookTree->SetBranchStatus("Timecas",1);
1684    q=(Float_t)(random())/(Float_t)(0x7fffffff);    fhBookTree->SetBranchStatus("Pathcas",1);
1685    w=0.7;    fhBookTree->SetBranchStatus("P0cas",1);
1686    //E0=Float_t(4064)/7;    fhBookTree->SetBranchStatus("Nthcard",1);
1687    E0=4064./7.;    fhBookTree->SetBranchStatus("Iparcard",1);
1688    S4=(Int_t)(4064.*(1.-exp(-int(Ert)/E0)));    fhBookTree->SetBranchStatus("Icard",1);
1689    //S4=Ert*7;    fhBookTree->SetBranchStatus("Xincard",1);
1690    i=S4/4;    fhBookTree->SetBranchStatus("Yincard",1);
1691    if(S4%4==0)    fhBookTree->SetBranchStatus("Zincard",1);
1692      S4v=S4+S4p;    fhBookTree->SetBranchStatus("Xoutcard",1);
1693    else if(S4%4==1) {    fhBookTree->SetBranchStatus("Youtcard",1);
1694      if(q<w) S4v=S4-1+S4p;    fhBookTree->SetBranchStatus("Zoutcard",1);
1695      else S4v=S4+1+S4p;    fhBookTree->SetBranchStatus("Erelcard",1);
1696    } else if(S4%4==2)    fhBookTree->SetBranchStatus("Timecard",1);
1697      S4v=S4+S4p;    fhBookTree->SetBranchStatus("Pathcard",1);
1698    else if(S4%4==3){    fhBookTree->SetBranchStatus("P0card",1);
1699      if(q<w) S4v=S4+1+S4p;  
1700      else S4v=S4-1+S4p;    // In this simpliefied approach we will assume that once
1701    }    // a particle releases > 0.5 mip in one of the 12 AC detectors it
1702      // will fire. We will furthermore assume that both cards read out
1703    cout << "Ert= " <<Ert<<"; S4v= "<<S4v<<"; S4= "<<S4<<endl;    // identical data.
1704    fDataS4[0]=S4v;//0xf028;  
1705    fDataS4[1]=0xd800;    // If you develop your digitization algorithm, you should start by
1706    fDataS4[2]=0x0300;    // identifying the information present in level2 (post-darth-vader)
1707    //  cout<<"  PMT  "<<NdFT<<"  "<<NdF<<endl;    // data.
1708    //cin>>ciao;  
1709  }    Float_t SumEcat[5];
1710      Float_t SumEcas[5];
1711      Float_t SumEcard[5];
1712      for (Int_t k= 0;k<5;k++){
1713  void Digitizer::DigitizeND(){      SumEcat[k]=0.;
1714    // creato:  S. Borisov, INFN Roma2 e MEPHI, Sett 2007      SumEcas[k]=0.;
1715    Int_t i=0;      SumEcard[k]=0.;
1716    UShort_t NdN=0;    };
1717    fhBookTree->SetBranchStatus("Nthnd",1);  
1718    fhBookTree->SetBranchStatus("Itubend",1);    if (Nthcat>50 || Nthcas>50 || Nthcard>50)
1719    fhBookTree->SetBranchStatus("Iparnd",1);        printf("*** ERROR AC! NthAC out of range!\n\n");
1720    fhBookTree->SetBranchStatus("Xinnd",1);  
1721    fhBookTree->SetBranchStatus("Yinnd",1);    // energy dependence on position (see file AcFitOutputDistancePmt.C by S.Orsi)
1722    fhBookTree->SetBranchStatus("Zinnd",1);    // based on J.Lundquist's calculations (PhD thesis, page 94)
1723    fhBookTree->SetBranchStatus("Xoutnd",1);    // function: [0]+[1]*atan([2]/(x+1)), where the 3 parameters are:
1724    fhBookTree->SetBranchStatus("Youtnd",1);    //   8.25470e-01   +- 1.79489e-02
1725    fhBookTree->SetBranchStatus("Zoutnd",1);    //   6.41609e-01   +- 2.65846e-02
1726    fhBookTree->SetBranchStatus("Erelnd",1);    //   9.81177e+00   +- 1.21284e+00
1727    fhBookTree->SetBranchStatus("Timend",1);    // hp: 1 minimum ionising particle at 35cm from the PMT releases 1mip
1728    fhBookTree->SetBranchStatus("Pathnd",1);    //
1729    fhBookTree->SetBranchStatus("P0nd",1);    // NB: the PMT positions are needed!
1730    //cout<<"n="<<Nthnd<<"  "<<NdN<<"\n";  
1731    for(i=0;i<Nthnd;i++){    // look in CAT
1732      if(Iparnd[i]==13){    //  for (UInt_t k= 0;k<50;k++){
1733        NdN++;    for (Int_t k= 0;k<Nthcat;k++){
1734      }      if (Erelcat[k] > 0)
1735    }        SumEcat[Icat[k]] += Erelcat[k];
1736    NdN=100;    };
1737    for(i=0;i<3;i++){  
1738      fDataND[2*i]=0x0000;    // look in CAS
1739      fDataND[2*i+1]=0x010F;    for (Int_t k= 0;k<Nthcas;k++){
1740    }      if (Erelcas[k] >0)
1741    fDataND[0]=0xFF00 & (256*NdN);        SumEcas[Icas[k]] += Erelcas[k];
1742  }    };
1743    
1744      // look in CARD
1745  void Digitizer::DigitizeDummy() {    for (Int_t k= 0;k<Nthcard;k++){
1746        if (Erelcard[k] >0)
1747    fhBookTree->SetBranchStatus("Enestrip",1);        SumEcard[Icard[k]] += Erelcard[k];
1748      };
1749    // dumy header  
1750    fDataDummy[0] = 0xCAAA;    // channel mapping              Hit Map
1751      // 1 CARD4                      0          LSB
1752    for (Int_t i=1; i<fDummybuffer; i++){    // 2 CAT2                       0
1753      fDataDummy[i] = 0xFFFF;    // 3 CAS1                       0
1754      //   printf("%0x  ",fDataDummy[i]);      // 4 NC                         0
1755      //if ((i+1)%8 ==0) cout << endl;    // 5 CARD2                      0
1756    }    // 6 CAT4                       1
1757  };    // 7 CAS4                       0  
1758      // 8 NC                         0
1759      // 9 CARD3                      0
1760  void Digitizer::WriteData(){    // 10 CAT3                      0
1761      // 11 CAS3                      0
1762    // Routine that writes the data to a binary file    // 12 NC                        0
1763    // PSCU data are already swapped    // 13 CARD1                     0
1764    fOutputfile.write(reinterpret_cast<char*>(fDataPSCU),sizeof(UShort_t)*fPSCUbuffer);    // 14 CAT1                      0
1765    // TRG    // 15 CAS2                      0
1766    fOutputfile.write(reinterpret_cast<char*>(fDataTrigger),sizeof(UChar_t)*153);    // 16 NC                        0          MSB
1767    // TOF  
1768    fOutputfile.write(reinterpret_cast<char*>(fDataTof),sizeof(UChar_t)*276);    // In the first version only the hit-map is filled, not the SR.
1769    // AC  
1770    UShort_t temp[1000000];    // Threshold: 0.8 MeV.
1771    memset(temp,0,sizeof(UShort_t)*1000000);  
1772    swab(fDataAC,temp,sizeof(UShort_t)*fACbuffer);  // WE MUST SWAP THE BYTES!!!    Float_t thr = 8e-4;
1773    fOutputfile.write(reinterpret_cast<char*>(temp),sizeof(UShort_t)*fACbuffer);  
1774    // CALO    fDataAC[3] = 0x0000;
1775    memset(temp,0,sizeof(UShort_t)*1000000);  
1776    swab(fDataCALO,temp,sizeof(UShort_t)*fCALOlength); // WE MUST SWAP THE BYTES!!!    if (SumEcas[0] > thr)  fDataAC[3]  = 0x0004;
1777    fOutputfile.write(reinterpret_cast<char*>(temp),sizeof(UShort_t)*fCALOlength);    if (SumEcas[1] > thr)  fDataAC[3] += 0x4000;
1778    // TRK    if (SumEcas[2] > thr)  fDataAC[3] += 0x0400;
1779    memset(temp,0,sizeof(UShort_t)*1000000);    if (SumEcas[3] > thr)  fDataAC[3] += 0x0040;  
1780    swab(fDataTrack,temp,sizeof(UShort_t)*fTracklength);  // WE MUST SWAP THE BYTES!!!  
1781    fOutputfile.write(reinterpret_cast<char*>(temp),sizeof(UShort_t)*fTracklength);    if (SumEcat[0] > thr)  fDataAC[3] += 0x2000;
1782    fTracklength=0;    if (SumEcat[1] > thr)  fDataAC[3] += 0x0002;
1783     // padding to 64 bytes    if (SumEcat[2] > thr)  fDataAC[3] += 0x0200;
1784    //    if (SumEcat[3] > thr)  fDataAC[3] += 0x0020;
1785    if ( fPadding ){    
1786      fOutputfile.write(reinterpret_cast<char*>(fDataPadding),sizeof(UChar_t)*fPadding);    if (SumEcard[0] > thr)  fDataAC[3] += 0x1000;
1787    };    if (SumEcard[1] > thr)  fDataAC[3] += 0x0010;
1788    // S4    if (SumEcard[2] > thr)  fDataAC[3] += 0x0100;
1789    memset(temp,0,sizeof(UShort_t)*1000000);    if (SumEcard[3] > thr)  fDataAC[3] += 0x0001;
1790    swab(fDataS4,temp,sizeof(UShort_t)*fS4buffer);  // WE MUST SWAP THE BYTES!!!  
1791    fOutputfile.write(reinterpret_cast<char*>(temp),sizeof(UShort_t)*fS4buffer);    fDataAC[67] = fDataAC[3];
1792    // ND  
1793    memset(temp,0,sizeof(UShort_t)*1000000);    // shift registers
1794    swab(fDataND,temp,sizeof(UShort_t)*fNDbuffer);  // WE MUST SWAP THE BYTES!!!    // the central bin is equal to the hitmap, all other bins in the shift register are 0
1795    fOutputfile.write(reinterpret_cast<char*>(temp),sizeof(UShort_t)*fNDbuffer);    for (UInt_t i=0; i<=15; i++){
1796  };      fDataAC[i+11] = 0x0000;  
1797        fDataAC[i+75] = 0x0000;
1798      }
1799  void Digitizer::ReadData(){    fDataAC[18] = fDataAC[3];
1800      fDataAC[82] = fDataAC[3];
1801    UShort_t InData[64];  
1802      //    for (Int_t i=0; i<fACbuffer; i++){
1803    // for debuggigng purposes only, write your own routine if you like (many    //    printf("%0x  ",fDataAC[i]);  
1804    // hardwired things.    //    if ((i+1)%8 ==0) cout << endl;
1805        //   }
1806    ifstream InputFile;  };
1807    
1808    // if (!InputFile) {  
1809    void Digitizer::DigitizeS4(){
1810    //     std::cout << "ERROR" << endl;    Int_t DEBUG=0;
1811    //     // An error occurred!    // creato:  S. Borisov, INFN Roma2 e MEPHI, Sett 2007
1812    //     // myFile.gcount() returns the number of bytes read.    TString ciao,modo="ns";
1813    //     // calling myFile.clear() will reset the stream state    Int_t i,j,t,NdF,pmt,NdFT,S4,S4v=0,S4p=32;
1814    //     // so it is usable again.    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;
1815    //   };    Xs[0]=-24.1;
1816      Xs[1]=24.1;
1817      Ys[0]=-24.1;
1818        Ys[1]=24.1;
1819    //InputFile.seekg(0);    Zs[0]=-0.5;
1820      Zs[1]=0.5;
1821    InputFile.open(fFilename, ios::in | ios::binary);    Yp[0]=-20.;
1822    //    fOutputfile.seekg(0);    Yp[2]=-1.;
1823    if (!InputFile.is_open()) std::cout << "ERROR" << endl;    Yp[4]=17.;
1824      for(i=0;i<3;i++)
1825    InputFile.seekg(0);      Yp[2*i+1]=Yp[2*i]+3;
1826        srand(time(NULL));
1827    for (Int_t k=0; k<=1000; k++){    // --- activate branches:
1828      InputFile.read(reinterpret_cast<char*>(InData),384*sizeof(UShort_t));    fhBookTree->SetBranchStatus("Nthtof",1);
1829      fhBookTree->SetBranchStatus("Ipltof",1);
1830      std::cout << "Read back: " << endl << endl;    fhBookTree->SetBranchStatus("Ipaddle",1);
1831      
1832      for (Int_t i=0; i<=384; i++){    fhBookTree->SetBranchStatus("Xintof",1);
1833        printf("%4x ", InData[i]);      fhBookTree->SetBranchStatus("Yintof",1);
1834        if ((i+1)%8 ==0) cout << endl;    fhBookTree->SetBranchStatus("Xouttof",1);
1835      }    fhBookTree->SetBranchStatus("Youttof",1);
1836      
1837    }    fhBookTree->SetBranchStatus("Ereltof",1);
1838    cout << endl;    fhBookTree->SetBranchStatus("Timetof",1);
1839    InputFile.close();    NdFT=0;
1840      Ert=0;
1841  };    for(i=0;i<Nthtof;i++){
1842        if(Ipltof[i]!=6) continue;
1843        Ert+=Ereltof[i];
1844              
1845  void Digitizer::DigitizeTrack() {      if(modo=="ns") continue;
1846  //std:: cout << "Entering DigitizeTrack " << endl;      NdF=Int_t(Ereltof[i]/E1);
1847  Float_t  AdcTrack[fNviews][fNstrips_view];  //  Vector of strips to be compressed      NdFT=0;
1848        X=Xintof[i];
1849  Int_t Iview;      Y=Yintof[i];
1850  Int_t Nstrip;      Z=(Float_t)(random())/(Float_t)(0x7fffffff)-0.5;
1851        //cout<<"XYZ "<<X<<"  "<<Y<<"  "<<Z<<endl;
1852    for (Int_t j=0; j<fNviews;j++) {      for(j=0;j<NdF;j++){
1853          q=(Float_t)random()/(Float_t)0x7fffffff;
1854      for (Int_t i=0; i<fNladder;i++) {        w=(Float_t)random()/(Float_t)0x7fffffff;
1855          // cout<<"qw "<<q<<" "<<w<<endl;
1856        Float_t commonN1=gRandom->Gaus(0.,fSigmaCommon);        V[0]=p*cos(6.28318*q);
1857        Float_t commonN2=gRandom->Gaus(0.,fSigmaCommon);        V[1]=p*sin(6.28318*q);
1858        for (Int_t k=0; k<fNstrips_ladder;k++) {        V[2]=p*(2.*w-1.);
1859        Nstrip=i*fNstrips_ladder+k;        pmt=0;
1860        AdcTrack[j][Nstrip]=gRandom->Gaus(fPedeTrack[j][Nstrip],fSigmaTrack[j][Nstrip]);        x=X;
1861        if(k<4*128) {AdcTrack[j][Nstrip] += commonN1;}  // full correlation of 4 VA1 Com. Noise        y=Y;
1862        else {AdcTrack[j][Nstrip] += commonN2;}   // full correlation of 4 VA1 Com. Noise        z=Z;
1863          while(pmt==0 && (x>Xs[0] && x<Xs[1])&&(y>Ys[0] && y<Ys[1])&&(z>Zs[0] && z<Zs[1])){
1864        };          l=0;
1865                  while(pmt==0 && (x>Xs[0] && x<Xs[1])&&(y>Ys[0] && y<Ys[1])&&(z>Zs[0] && z<Zs[1])){
1866                      x+=V[0];
1867      };            y+=V[1];
1868              z+=V[2];
1869              l+=p;
1870    };            //cout<<x<<"  "<<y<<"  "<<z<<"  "<<l<<endl;
1871              //cin>>ciao;
1872            }
1873    fhBookTree->SetBranchStatus("Nstrpx",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)){
1874    fhBookTree->SetBranchStatus("Npstripx",1);            for(t=0;t<3;t++){
1875    fhBookTree->SetBranchStatus("Ntstripx",1);              if(y>=Yp[2*t] && y<Yp[2*t+1]){
1876    fhBookTree->SetBranchStatus("Istripx",1);                if(pmt==0)NdFT++;
1877    fhBookTree->SetBranchStatus("Qstripx",1);                pmt=1;
1878    fhBookTree->SetBranchStatus("Xstripx",1);                //cout<<NdFT<<endl;
1879    fhBookTree->SetBranchStatus("Nstrpy",1);                break;
1880    fhBookTree->SetBranchStatus("Npstripy",1);              }
1881    fhBookTree->SetBranchStatus("Ntstripy",1);            }
1882    fhBookTree->SetBranchStatus("Istripy",1);            if(pmt==1)break;
1883    fhBookTree->SetBranchStatus("Qstripy",1);            V[0]=-V[0];
1884    fhBookTree->SetBranchStatus("Ystripy",1);          }
1885            q=(Float_t)random()/(Float_t)0x7fffffff;
1886            w=1-exp(-l/l0);
1887            if(q<w)break;
1888    Float_t ADCfull;          q=(Float_t)random()/(Float_t)0x7fffffff;
1889    Int_t iladd=0;          w=0.5;
1890    for (Int_t ix=0; ix<Nstrpx;ix++) {          if(q<w)break;
1891    Iview=Npstripx[ix]*2-1;          if((x>Xs[0]+p && x<Xs[1]-p)&&(y<Ys[0]+p || y>Ys[1]-p)&&(z>Zs[0]+p && z<Zs[1]-p))V[1]=-V[1];
1892    Nstrip=(Int_t)Istripx[ix]-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];
1893    if(Nstrip<fNstrips_ladder) iladd=0;          x+=V[0];
1894    if((Nstrip>=fNstrips_ladder)&&(Nstrip<2*fNstrips_ladder)) iladd=1;          y+=V[1];
1895    if((Nstrip>=2*fNstrips_ladder)&&(Nstrip<3*fNstrips_ladder)) iladd=2;          z+=V[2];
1896    ADCfull=AdcTrack[Iview][Nstrip] += Qstripx[ix]*fMipCor[iladd][Iview];          l=0;
1897    AdcTrack[Iview][Nstrip] *= SaturationTrack(ADCfull);          //cout<<x<<"  "<<y<<"  "<<z<<"  "<<l<<endl;
1898                    //cin>>ciao;
1899    };        }
1900        }
1901      }
1902    for (Int_t iy=0; iy<Nstrpy;iy++) {    Ert=Ert/0.002;
1903    Iview=Npstripy[iy]*2-2;    q=(Float_t)(random())/(Float_t)0x7fffffff;
1904    Nstrip=(Int_t)Istripy[iy]-1;    w=0.7;
1905    if(Nstrip<fNstrips_ladder) iladd=0;    //E0=(Float_t)(4064./7.);
1906    if((Nstrip>=fNstrips_ladder)&&(Nstrip<2*fNstrips_ladder)) iladd=1;    E0=4064./7.;
1907    if((Nstrip>=2*fNstrips_ladder)&&(Nstrip<3*fNstrips_ladder)) iladd=2;    if(Ert<1) S4=0;
1908    ADCfull=AdcTrack[Iview][Nstrip] -= Qstripy[iy]*fMipCor[iladd][Iview];    else S4=(Int_t)(4064.*(1.-exp(-(Ert-1.)/E0)));
1909    AdcTrack[Iview][Nstrip] *= SaturationTrack(ADCfull);    i=S4/4;
1910      if(S4%4==0)
1911    };        S4v=S4+S4p;
1912            else if(S4%4==1){
1913  CompressTrackData(AdcTrack);  // Compress and Digitize data of one Ladder  in turn for all ladders      if(q<w) S4v=S4-1+S4p;
1914        else S4v=S4+1+S4p;
1915  };    } else if(S4%4==2) S4v=S4+S4p;
1916      else if(S4%4==3){
1917        if(q<w) S4v=S4+1+S4p;
1918        else S4v=S4-1+S4p;
1919  void Digitizer::DigitizeTrackCalib(Int_t ii) {    }
1920      if (DEBUG)
1921  std:: cout << "Entering DigitizeTrackCalib " << ii << endl;      cout<<"Ert_S4 = " << Ert << " --- S4v = " << S4v << endl;
1922  if( (ii!=1)&&(ii!=2) ) {    fDataS4[0]=S4v;//0xf028;
1923   std:: cout << "error wrong DigitizeTrackCalib argument" << endl;    fDataS4[1]=0xd800;
1924   return;    fDataS4[2]=0x0300;
1925  };    //cout<<"  PMT  "<<NdFT<<"  "<<NdF<<endl;
1926      //cin>>ciao;
1927  memset(fDataTrack,0,sizeof(UShort_t)*fTRACKbuffer);  }
1928  fTracklength=0;  
1929    
1930  UShort_t Dato;  
1931    void Digitizer::DigitizeND(){
1932  Float_t dato1;    // creato:  S. Borisov, INFN Roma2 e MEPHI, Sett 2007
1933  Float_t dato2;    Int_t i=0;
1934  Float_t dato3;    UShort_t NdN=0;
1935  Float_t dato4;    fhBookTree->SetBranchStatus("Nthnd",1);
1936      fhBookTree->SetBranchStatus("Itubend",1);
1937  UShort_t DatoDec;    fhBookTree->SetBranchStatus("Iparnd",1);  
1938  UShort_t DatoDec1;    fhBookTree->SetBranchStatus("Xinnd",1);
1939  UShort_t DatoDec2;    fhBookTree->SetBranchStatus("Yinnd",1);
1940  UShort_t DatoDec3;    fhBookTree->SetBranchStatus("Zinnd",1);
1941  UShort_t DatoDec4;    fhBookTree->SetBranchStatus("Xoutnd",1);
1942      fhBookTree->SetBranchStatus("Youtnd",1);
1943  UShort_t EVENT_CAL;    fhBookTree->SetBranchStatus("Zoutnd",1);
1944  UShort_t PED_L1;    fhBookTree->SetBranchStatus("Erelnd",1);
1945  UShort_t ReLength;    fhBookTree->SetBranchStatus("Timend",1);
1946  UShort_t OveCheckCode;    fhBookTree->SetBranchStatus("Pathnd",1);
1947  //UShort_t PED_L2;    fhBookTree->SetBranchStatus("P0nd",1);
1948  //UShort_t PED_L3HI;    //cout<<"n="<<Nthnd<<"  "<<NdN<<"\n";
1949  //UShort_t PED_L3LO;    for(i=0;i<Nthnd;i++){
1950  //UShort_t SIG_L1HI;      if(Iparnd[i]==13){
1951  //UShort_t SIG_L1LO;        NdN++;
1952  //UShort_t SIG_L2HI;      }
1953  //UShort_t SIG_L2LO;    }
1954  //UShort_t SIG_L3;    //NdN=100; //only for debug
1955  //UShort_t BAD_L1;  
1956  //UShort_t BAD_L2LO;    for(i=0;i<3;i++){
1957  //UShort_t BAD_L3HI;      fDataND[2*i]=0x0000;
1958  //UShort_t BAD_L3LO;      fDataND[2*i+1]=0x010F;
1959  //UShort_t FLAG;    }
1960      fDataND[0]=0xFF00 & (256*NdN);
1961    }
1962    Int_t DSPpos;  
1963    for (Int_t j=ii-1; j<fNviews;j+=2) {  
1964    UShort_t CkSum=0;  void Digitizer::DigitizeDummy() {
1965    // here skip the dsp header and his trailer , to be written later  
1966    DSPpos=fTracklength;    fhBookTree->SetBranchStatus("Enestrip",1);
1967    fTracklength=fTracklength+13+3;  
1968      // dumy header
1969      fDataDummy[0] = 0xCAAA;
1970      for (Int_t i=0; i<fNladder;i++) {  
1971        for (Int_t k=0; k<fNstrips_ladder;k++) {    for (Int_t i=1; i<fDummybuffer; i++){
1972        // write in buffer the current LADDER      fDataDummy[i] = 0xFFFF;
1973        Dato=(UShort_t)fPedeTrack[j][i*fNstrips_ladder+k];      //   printf("%0x  ",fDataDummy[i]);  
1974        dato1=fPedeTrack[j][i*fNstrips_ladder+k]-Dato;      //if ((i+1)%8 ==0) cout << endl;
1975      }
1976        DatoDec1=(UShort_t)(dato1*2);  };
1977        dato2=dato1*2-DatoDec1;  
1978    
1979        DatoDec2=(UShort_t)(dato2*2);  void Digitizer::WriteData(){
1980        dato3=dato2*2-DatoDec2;  
1981            // Routine that writes the data to a binary file
1982        DatoDec3=(UShort_t)(dato3*2);    // PSCU data are already swapped
1983        dato4=dato3*2-DatoDec3;    fOutputfile.write(reinterpret_cast<char*>(fDataPSCU),sizeof(UShort_t)*fPSCUbuffer);
1984            // TRG
1985        DatoDec4=(UShort_t)(dato4*2);    fOutputfile.write(reinterpret_cast<char*>(fDataTrigger),sizeof(UChar_t)*153);
1986            // TOF
1987        DatoDec=DatoDec1*0x0008+DatoDec2*0x0004+DatoDec3*0x0002+DatoDec4*0x0001;    fOutputfile.write(reinterpret_cast<char*>(fDataTof),sizeof(UChar_t)*276);
1988        fDataTrack[fTracklength]=( (Dato << 4) | (DatoDec & 0x000F) );    // AC
1989        CkSum=CkSum^fDataTrack[fTracklength];    UShort_t temp[1000000];
1990        fTracklength++;    memset(temp,0,sizeof(UShort_t)*1000000);
1991        };    swab(fDataAC,temp,sizeof(UShort_t)*fACbuffer);  // WE MUST SWAP THE BYTES!!!
1992      fOutputfile.write(reinterpret_cast<char*>(temp),sizeof(UShort_t)*fACbuffer);
1993        for (Int_t k=0; k<fNstrips_ladder;k++) {    // CALO
1994        // write in buffer the current LADDER    memset(temp,0,sizeof(UShort_t)*1000000);
1995        Dato=(UShort_t)fSigmaTrack[j][i*fNstrips_ladder+k];    swab(fDataCALO,temp,sizeof(UShort_t)*fCALOlength); // WE MUST SWAP THE BYTES!!!
1996        dato1=fSigmaTrack[j][i*fNstrips_ladder+k]-Dato;    fOutputfile.write(reinterpret_cast<char*>(temp),sizeof(UShort_t)*fCALOlength);
1997      // TRK
1998        DatoDec1=(UShort_t)(dato1*2);    memset(temp,0,sizeof(UShort_t)*1000000);
1999        dato2=dato1*2-DatoDec1;    swab(fDataTrack,temp,sizeof(UShort_t)*fTracklength);  // WE MUST SWAP THE BYTES!!!
2000      fOutputfile.write(reinterpret_cast<char*>(temp),sizeof(UShort_t)*fTracklength);
2001        DatoDec2=(UShort_t)(dato2*2);    fTracklength=0;
2002        dato3=dato2*2-DatoDec2;     // padding to 64 bytes
2003            //
2004        DatoDec3=(UShort_t)(dato3*2);    if ( fPadding ){
2005        dato4=dato3*2-DatoDec3;      fOutputfile.write(reinterpret_cast<char*>(fDataPadding),sizeof(UChar_t)*fPadding);
2006            };
2007        DatoDec4=(UShort_t)(dato4*2);    // S4
2008            memset(temp,0,sizeof(UShort_t)*1000000);
2009        DatoDec=DatoDec1*0x0008+DatoDec2*0x0004+DatoDec3*0x0002+DatoDec4*0x0001;    swab(fDataS4,temp,sizeof(UShort_t)*fS4buffer);  // WE MUST SWAP THE BYTES!!!
2010            fOutputfile.write(reinterpret_cast<char*>(temp),sizeof(UShort_t)*fS4buffer);
2011        fDataTrack[fTracklength]=( (Dato << 4) | (DatoDec & 0x000F) );    // ND
2012        CkSum=CkSum^fDataTrack[fTracklength];    memset(temp,0,sizeof(UShort_t)*1000000);
2013        fTracklength++;    swab(fDataND,temp,sizeof(UShort_t)*fNDbuffer);  // WE MUST SWAP THE BYTES!!!
2014        };    fOutputfile.write(reinterpret_cast<char*>(temp),sizeof(UShort_t)*fNDbuffer);
2015          };
2016        for (Int_t k=0; k<64;k++) {  
2017        fDataTrack[fTracklength]=0x0000;  
2018        CkSum=CkSum^fDataTrack[fTracklength];  void Digitizer::ReadData(){
2019        fTracklength++;  
2020      UShort_t InData[64];
2021        };  
2022        // end ladder    // for debuggigng purposes only, write your own routine if you like (many
2023      // hardwired things.
2024      // write in buffer the end ladder word    
2025      if(i==0) fDataTrack[fTracklength]=0x1807;    ifstream InputFile;
2026      if(i==1) fDataTrack[fTracklength]=0x1808;  
2027      if(i==2) fDataTrack[fTracklength]=0x1809;    // if (!InputFile) {
2028      CkSum=CkSum^fDataTrack[fTracklength];  
2029      fTracklength++;    //     std::cout << "ERROR" << endl;
2030      //     // An error occurred!
2031      // write in buffer the TRAILER    //     // myFile.gcount() returns the number of bytes read.
2032      ReLength=(UShort_t)((fNstrips_ladder*2+64+1)*2+3);    //     // calling myFile.clear() will reset the stream state
2033      OveCheckCode=0x0000;    //     // so it is usable again.
2034      //   };
2035      fDataTrack[fTracklength]=0x0000;  
2036      fTracklength++;  
2037        
2038      fDataTrack[fTracklength]=(ReLength >> 8);    //InputFile.seekg(0);
2039      fTracklength++;  
2040        InputFile.open(fFilename, ios::in | ios::binary);
2041      fDataTrack[fTracklength]=( (ReLength << 8) | (OveCheckCode & 0x00FF) );    //    fOutputfile.seekg(0);
2042      fTracklength++;    if (!InputFile.is_open()) std::cout << "ERROR" << endl;
2043    
2044      // end TRAILER            InputFile.seekg(0);
2045      };    
2046      for (Int_t k=0; k<=1000; k++){
2047    // write in buffer the DSP header      InputFile.read(reinterpret_cast<char*>(InData),384*sizeof(UShort_t));
2048    
2049    fDataTrack[DSPpos]=(0xE800 | ( ((j+1) << 3) | 0x0005) );      std::cout << "Read back: " << endl << endl;
2050      
2051    fDataTrack[DSPpos+1]=0x01A9;      for (Int_t i=0; i<=384; i++){
2052          printf("%4x ", InData[i]);  
2053    fDataTrack[DSPpos+2]=0x8740;        if ((i+1)%8 ==0) cout << endl;
2054          }
2055    EVENT_CAL=0;  
2056    fDataTrack[DSPpos+3]=(0x1A00 | ( (0x03FF & EVENT_CAL)>> 1) );    }
2057        cout << endl;
2058    PED_L1=0;    InputFile.close();
2059    fDataTrack[DSPpos+4]=( ((EVENT_CAL << 15) | 0x5002 ) | ((0x03FF & PED_L1) << 2) );  
2060      };
2061    // FROM HERE WE WRITE AS ALL VARIABLE apart CkSum are =0    
2062    
2063    fDataTrack[DSPpos+5]=0x8014;  
2064      void Digitizer::DigitizeTrack() {
2065    fDataTrack[DSPpos+6]=0x00A0;  //std:: cout << "Entering DigitizeTrack " << endl;
2066      Float_t  AdcTrack[fNviews][fNstrips_view];  //  Vector of strips to be compressed
2067    fDataTrack[DSPpos+7]=0x0500;  
2068      Int_t Iview;
2069    fDataTrack[DSPpos+8]=0x2801;  Int_t Nstrip;
2070      
2071    fDataTrack[DSPpos+9]=0x400A;    for (Int_t j=0; j<fNviews;j++) {
2072      
2073    fDataTrack[DSPpos+10]=0x0050;      for (Int_t i=0; i<fNladder;i++) {
2074      
2075    CkSum=(CkSum >> 8)^(CkSum&0x00FF);        Float_t commonN1=gRandom->Gaus(0.,fSigmaCommon);
2076    fDataTrack[DSPpos+11]=(0x0280 | (CkSum >> 3));        Float_t commonN2=gRandom->Gaus(0.,fSigmaCommon);
2077            for (Int_t k=0; k<fNstrips_ladder;k++) {
2078    fDataTrack[DSPpos+12]=(0x1FFF | (CkSum << 13) );        Nstrip=i*fNstrips_ladder+k;
2079          AdcTrack[j][Nstrip]=gRandom->Gaus(fPedeTrack[j][Nstrip],fSigmaTrack[j][Nstrip]);
2080    // end dsp header        if(k<4*128) {AdcTrack[j][Nstrip] += commonN1;}  // full correlation of 4 VA1 Com. Noise
2081          else {AdcTrack[j][Nstrip] += commonN2;}   // full correlation of 4 VA1 Com. Noise
2082    // write in buffer the TRAILER  
2083            };
2084    ReLength=(UShort_t)((13*2)+3);        
2085    OveCheckCode=0x0000;          
2086    fDataTrack[DSPpos+13]=0x0000;      };
2087    
2088    fDataTrack[DSPpos+14]=(ReLength >> 8);  
2089        };
2090    fDataTrack[DSPpos+15]=( (ReLength << 8) | (OveCheckCode & 0x00FF) );  
2091      
2092    // end TRAILER    fhBookTree->SetBranchStatus("Nstrpx",1);
2093      fhBookTree->SetBranchStatus("Npstripx",1);
2094      fhBookTree->SetBranchStatus("Ntstripx",1);
2095      fhBookTree->SetBranchStatus("Istripx",1);
2096        fhBookTree->SetBranchStatus("Qstripx",1);
2097    // end DSP        fhBookTree->SetBranchStatus("Xstripx",1);
2098    };        fhBookTree->SetBranchStatus("Nstrpy",1);
2099      fhBookTree->SetBranchStatus("Npstripy",1);
2100      fhBookTree->SetBranchStatus("Ntstripy",1);
2101      fhBookTree->SetBranchStatus("Istripy",1);
2102  };    fhBookTree->SetBranchStatus("Qstripy",1);
2103      fhBookTree->SetBranchStatus("Ystripy",1);
2104  void Digitizer::WriteTrackCalib() {  
2105    
2106    
2107  std:: cout << " Entering WriteTrackCalib " << endl;    Float_t ADCfull;
2108      Int_t iladd=0;
2109  fOutputfile.write(reinterpret_cast<char*>(fDataPSCU),sizeof(UShort_t)*fPSCUbuffer);    for (Int_t ix=0; ix<Nstrpx;ix++) {
2110      Iview=Npstripx[ix]*2-1;
2111  UShort_t temp[1000000];    Nstrip=(Int_t)Istripx[ix]-1;
2112  memset(temp,0,sizeof(UShort_t)*1000000);    if(Nstrip<fNstrips_ladder) iladd=0;
2113  swab(fDataTrack,temp,sizeof(UShort_t)*fTracklength);  // WE MUST SWAP THE BYTES!!!    if((Nstrip>=fNstrips_ladder)&&(Nstrip<2*fNstrips_ladder)) iladd=1;
2114  fOutputfile.write(reinterpret_cast<char*>(temp),sizeof(UShort_t)*fTracklength);    if((Nstrip>=2*fNstrips_ladder)&&(Nstrip<3*fNstrips_ladder)) iladd=2;
2115  fTracklength=0;    ADCfull=AdcTrack[Iview][Nstrip] += Qstripx[ix]*fMipCor[iladd][Iview];
2116  if ( fPadding ){    AdcTrack[Iview][Nstrip] *= SaturationTrack(ADCfull);
2117        fOutputfile.write(reinterpret_cast<char*>(fDataPadding),sizeof(UChar_t)*fPadding);  
2118  };    };
2119    
2120  };  
2121      for (Int_t iy=0; iy<Nstrpy;iy++) {
2122      Iview=Npstripy[iy]*2-2;
2123  void Digitizer::ClearTrackCalib() {    Nstrip=(Int_t)Istripy[iy]-1;
2124      if(Nstrip<fNstrips_ladder) iladd=0;
2125  std:: cout << "Entering ClearTrackCalib " << endl;    if((Nstrip>=fNstrips_ladder)&&(Nstrip<2*fNstrips_ladder)) iladd=1;
2126          if((Nstrip>=2*fNstrips_ladder)&&(Nstrip<3*fNstrips_ladder)) iladd=2;
2127        ADCfull=AdcTrack[Iview][Nstrip] -= Qstripy[iy]*fMipCor[iladd][Iview];
2128  };    AdcTrack[Iview][Nstrip] *= SaturationTrack(ADCfull);
2129    
2130      };  
2131  void Digitizer::LoadTrackCalib() {        
2132  std:: cout << "Entering LoadTrackCalib " << endl;  CompressTrackData(AdcTrack);  // Compress and Digitize data of one Ladder  in turn for all ladders
2133    
2134  // Generate the pedestals and sigmas according to parametrization  };
2135    for (Int_t j=0; j<fNviews;j++) {  
2136      for (Int_t i=0; i<fNstrips_view;i++) {  
2137        
2138      if((j+1)%2==0) {  void Digitizer::DigitizeTrackCalib(Int_t ii) {
2139      fPedeTrack[j][i]=gRandom->Gaus(fAvePedex,fSigmaPedex);  
2140      fSigmaTrack[j][i]=gRandom->Gaus(fAveSigmax,fSigmaSigmax);  std:: cout << "Entering DigitizeTrackCalib " << ii << endl;
2141      };  if( (ii!=1)&&(ii!=2) ) {
2142      if((j+1)%2==1) {   std:: cout << "error wrong DigitizeTrackCalib argument" << endl;
2143      fPedeTrack[j][i]=gRandom->Gaus(fAvePedey,fSigmaPedey);   return;
2144      fSigmaTrack[j][i]=gRandom->Gaus(fAveSigmay,fSigmaSigmay);  };
2145      };  
2146        memset(fDataTrack,0,sizeof(UShort_t)*fTRACKbuffer);
2147      };  fTracklength=0;
2148    };      
2149    UShort_t Dato;
2150      
2151      Float_t dato1;
2152  };  Float_t dato2;
2153    Float_t dato3;
2154  void Digitizer::LoadMipCor() {  Float_t dato4;
2155  std:: cout << "Entering LoadMipCor" << endl;  
2156    Float_t xfactor=1./151.6*1.04;  UShort_t DatoDec;
2157    Float_t yfactor=1./152.1;  UShort_t DatoDec1;
2158    UShort_t DatoDec2;
2159    fMipCor[0][0]=140.02*yfactor;  UShort_t DatoDec3;
2160    fMipCor[0][1]=140.99*xfactor;  UShort_t DatoDec4;
2161    fMipCor[0][2]=134.48*yfactor;  
2162    fMipCor[0][3]=144.41*xfactor;  UShort_t EVENT_CAL;
2163    fMipCor[0][4]=140.74*yfactor;  UShort_t PED_L1;
2164    fMipCor[0][5]=142.28*xfactor;  UShort_t ReLength;
2165    fMipCor[0][6]=134.53*yfactor;  UShort_t OveCheckCode;
2166    fMipCor[0][7]=140.63*xfactor;  //UShort_t PED_L2;
2167    fMipCor[0][8]=135.55*yfactor;  //UShort_t PED_L3HI;
2168    fMipCor[0][9]=138.00*xfactor;  //UShort_t PED_L3LO;
2169    fMipCor[0][10]=154.95*yfactor;  //UShort_t SIG_L1HI;
2170    fMipCor[0][11]=158.44*xfactor;  //UShort_t SIG_L1LO;
2171      //UShort_t SIG_L2HI;
2172      //UShort_t SIG_L2LO;
2173    fMipCor[1][0]=136.07*yfactor;  //UShort_t SIG_L3;
2174    fMipCor[1][1]=135.59*xfactor;  //UShort_t BAD_L1;
2175    fMipCor[1][2]=142.69*yfactor;  //UShort_t BAD_L2LO;
2176    fMipCor[1][3]=138.19*xfactor;  //UShort_t BAD_L3HI;
2177    fMipCor[1][4]=137.35*yfactor;  //UShort_t BAD_L3LO;
2178    fMipCor[1][5]=140.23*xfactor;  //UShort_t FLAG;
2179    fMipCor[1][6]=153.15*yfactor;  
2180    fMipCor[1][7]=151.42*xfactor;  
2181    fMipCor[1][8]=129.76*yfactor;    Int_t DSPpos;
2182    fMipCor[1][9]=140.63*xfactor;    for (Int_t j=ii-1; j<fNviews;j+=2) {
2183    fMipCor[1][10]=157.87*yfactor;    UShort_t CkSum=0;
2184    fMipCor[1][11]=153.64*xfactor;    // here skip the dsp header and his trailer , to be written later
2185      DSPpos=fTracklength;
2186    fMipCor[2][0]=134.98*yfactor;    fTracklength=fTracklength+13+3;
2187    fMipCor[2][1]=143.95*xfactor;  
2188    fMipCor[2][2]=140.23*yfactor;  
2189    fMipCor[2][3]=138.88*xfactor;      for (Int_t i=0; i<fNladder;i++) {
2190    fMipCor[2][4]=137.95*yfactor;        for (Int_t k=0; k<fNstrips_ladder;k++) {
2191    fMipCor[2][5]=134.87*xfactor;        // write in buffer the current LADDER
2192    fMipCor[2][6]=157.56*yfactor;        Dato=(UShort_t)fPedeTrack[j][i*fNstrips_ladder+k];
2193    fMipCor[2][7]=157.31*xfactor;        dato1=fPedeTrack[j][i*fNstrips_ladder+k]-Dato;
2194    fMipCor[2][8]=141.37*yfactor;  
2195    fMipCor[2][9]=143.39*xfactor;        DatoDec1=(UShort_t)(dato1*2);
2196    fMipCor[2][10]=156.15*yfactor;        dato2=dato1*2-DatoDec1;
2197    fMipCor[2][11]=158.79*xfactor;  
2198          DatoDec2=(UShort_t)(dato2*2);
2199  /*        dato3=dato2*2-DatoDec2;
2200    for (Int_t j=0; j<fNviews;j++) {        
2201      for (Int_t i=0; i<fNstrips_view;i++) {        DatoDec3=(UShort_t)(dato3*2);
2202      fMipCor[j][i]=1.;        dato4=dato3*2-DatoDec3;
2203      };        
2204    };            DatoDec4=(UShort_t)(dato4*2);
2205          
2206              DatoDec=DatoDec1*0x0008+DatoDec2*0x0004+DatoDec3*0x0002+DatoDec4*0x0001;
2207  */          fDataTrack[fTracklength]=( (Dato << 4) | (DatoDec & 0x000F) );
2208  };        CkSum=CkSum^fDataTrack[fTracklength];
2209          fTracklength++;
2210  void Digitizer::CompressTrackData(Float_t AdcTrack[fNviews][fNstrips_view]) {        };
2211  // copy of the corresponding compression fortran routine + new digitization  
2212  // std:: cout << "Entering CompressTrackData " << endl;        for (Int_t k=0; k<fNstrips_ladder;k++) {
2213  Int_t oldval=0;        // write in buffer the current LADDER
2214  Int_t newval=0;        Dato=(UShort_t)fSigmaTrack[j][i*fNstrips_ladder+k];
2215  Int_t trasmesso=0;        dato1=fSigmaTrack[j][i*fNstrips_ladder+k]-Dato;
2216  Int_t ntrastot=0;  
2217  Float_t real;        DatoDec1=(UShort_t)(dato1*2);
2218  Float_t inte;        dato2=dato1*2-DatoDec1;
2219  Int_t cercacluster=0;  
2220  Int_t kt=0;        DatoDec2=(UShort_t)(dato2*2);
2221  static const int DSPbufferSize = 4000; // 13 bit buffer to be rearranged in 16 bit Track buffer        dato3=dato2*2-DatoDec2;
2222  UShort_t DataDSP[DSPbufferSize];  // 13 bit  buffer to be rearranged in 16 bit Track buffer        
2223  UShort_t DSPlength;  // 13 bit buffer to be rearranged in 16 bit Track buffer        DatoDec3=(UShort_t)(dato3*2);
2224          dato4=dato3*2-DatoDec3;
2225  memset(fDataTrack,0,sizeof(UShort_t)*fTRACKbuffer); // probably not necessary becouse already done ?        
2226  fTracklength=0;        DatoDec4=(UShort_t)(dato4*2);
2227          
2228    for (Int_t iv=0; iv<fNviews;iv++) {        DatoDec=DatoDec1*0x0008+DatoDec2*0x0004+DatoDec3*0x0002+DatoDec4*0x0001;
2229    memset(DataDSP,0,sizeof(UShort_t)*DSPbufferSize);        
2230    DSPlength=16; // skip the header, to be written later        fDataTrack[fTracklength]=( (Dato << 4) | (DatoDec & 0x000F) );
2231    UShort_t CheckSum=0;        CkSum=CkSum^fDataTrack[fTracklength];
2232  // write dsp header on buffer        fTracklength++;
2233          };
2234  //    fDataTrack[fTracklength]=0xE805;        
2235  //    fTracklength++;        for (Int_t k=0; k<64;k++) {
2236          fDataTrack[fTracklength]=0x0000;
2237  //    fDataTrack[fTracklength]=0x01A9;        CkSum=CkSum^fDataTrack[fTracklength];
2238  //    fTracklength++;        fTracklength++;
2239    
2240  // end dsp header        };
2241          // end ladder
2242     //  
2243     // INIZIO VISTA IV - TAKE PROPER ACTION      // write in buffer the end ladder word
2244     //      if(i==0) fDataTrack[fTracklength]=0x1807;
2245        if(i==1) fDataTrack[fTracklength]=0x1808;
2246        if(i==2) fDataTrack[fTracklength]=0x1809;
2247        CkSum=CkSum^fDataTrack[fTracklength];
2248      for (Int_t ladder=0; ladder<fNladder;ladder++) {      fTracklength++;
2249        Int_t k=0;  
2250        while (k<fNstrips_ladder) {      // write in buffer the TRAILER
2251          // compress write in buffer the current LADDER      ReLength=(UShort_t)((fNstrips_ladder*2+64+1)*2+3);
2252          if ( k == 0)  {      OveCheckCode=0x0000;
2253            real=modff(AdcTrack[iv][ladder*fNstrips_ladder+k],&inte);  
2254            if (real > 0.5) inte=inte+1;      fDataTrack[fTracklength]=0x0000;
2255            newval=(Int_t)inte -(Int_t)fPedeTrack[iv][ladder*fNstrips_ladder+k];      fTracklength++;
2256            // first strip of ladder is transmitted    
2257            // DC_TOT first " << AdcTrack[iv][ladder*fNstrips_ladder+k] << endl;      fDataTrack[fTracklength]=(ReLength >> 8);
2258            DataDSP[DSPlength]=( ((UShort_t)inte) & 0x0FFF);      fTracklength++;
2259            DSPlength++;      
2260            ntrastot++;      fDataTrack[fTracklength]=( (ReLength << 8) | (OveCheckCode & 0x00FF) );
2261            trasmesso=1;      fTracklength++;
2262            oldval=newval;  
2263            kt=k;      // end TRAILER        
2264            k++;      };
2265            continue;  
2266          };    // write in buffer the DSP header
2267          real=modff(AdcTrack[iv][ladder*fNstrips_ladder+k],&inte);  
2268          if (real > 0.5) inte=inte+1;    fDataTrack[DSPpos]=(0xE800 | ( ((j+1) << 3) | 0x0005) );
2269          newval=(Int_t)inte -(Int_t)(fPedeTrack[iv][ladder*fNstrips_ladder+k]);    
2270          cercacluster=1; // ?????????    fDataTrack[DSPpos+1]=0x01A9;
2271          if (cercacluster==1) {  
2272              fDataTrack[DSPpos+2]=0x8740;
2273   // controlla l'ordine di tutti queste strip ladder e DSP !!!!!!!          
2274            Int_t diff=0;    EVENT_CAL=0;
2275      fDataTrack[DSPpos+3]=(0x1A00 | ( (0x03FF & EVENT_CAL)>> 1) );
2276                
2277            switch ((iv+1)%2) {    PED_L1=0;
2278            case 0: diff=newval-oldval;    fDataTrack[DSPpos+4]=( ((EVENT_CAL << 15) | 0x5002 ) | ((0x03FF & PED_L1) << 2) );
2279            break;    
2280            case 1: diff=oldval-newval;        // FROM HERE WE WRITE AS ALL VARIABLE apart CkSum are =0  
2281            break;  
2282            };    fDataTrack[DSPpos+5]=0x8014;
2283      
2284            if (diff>fCutclu*(Int_t)fSigmaTrack[iv][ladder*fNstrips_ladder+k]) {    fDataTrack[DSPpos+6]=0x00A0;
2285              Int_t clval=newval;    
2286              Int_t klp=k;        // go on to search for maximum    fDataTrack[DSPpos+7]=0x0500;
2287              klp++;    
2288      fDataTrack[DSPpos+8]=0x2801;
2289              while(klp<fNstrips_ladder) {    
2290                real=modff(AdcTrack[iv][ladder*fNstrips_ladder+klp],&inte);    fDataTrack[DSPpos+9]=0x400A;
2291                if (real > 0.5) inte=inte+1;    
2292                Int_t clvalp=(Int_t)inte -(Int_t)fPedeTrack[iv][ladder*fNstrips_ladder+klp];    fDataTrack[DSPpos+10]=0x0050;
2293                if((iv+1)%2==0) {    
2294      CkSum=(CkSum >> 8)^(CkSum&0x00FF);
2295                  if(clvalp>clval) {    fDataTrack[DSPpos+11]=(0x0280 | (CkSum >> 3));
2296                     clval=clvalp;    
2297                     k=klp;}    fDataTrack[DSPpos+12]=(0x1FFF | (CkSum << 13) );
2298                  else break; // max of cluster found  
2299      // end dsp header
2300                } else {    
2301      // write in buffer the TRAILER
2302                  if(clvalp<clval) {    
2303                     clval=clvalp;    ReLength=(UShort_t)((13*2)+3);
2304                     k=klp;}    OveCheckCode=0x0000;
2305                  else break; // max of cluster found    fDataTrack[DSPpos+13]=0x0000;
2306    
2307                };    fDataTrack[DSPpos+14]=(ReLength >> 8);
2308                    
2309                klp++;    fDataTrack[DSPpos+15]=( (ReLength << 8) | (OveCheckCode & 0x00FF) );
2310              };    
2311      // end TRAILER
2312              Int_t kl1=k-fNclst; // max of cluster (or end of ladder ?)  
2313              trasmesso=0;  
2314              if(kl1<0)  kl1=0;  
2315      
2316              if(kt>=kl1) kl1=kt+1;    // end DSP    
2317              if( (kt+1)==kl1 ) trasmesso=1;    };    
2318    
2319    
2320                
2321              Int_t kl2=k+fNclst;  };
2322              if(kl2>=fNstrips_ladder) kl2=fNstrips_ladder-1;  
2323    void Digitizer::WriteTrackCalib() {
2324              for(Int_t klt=kl1 ; klt<=kl2 ; klt++) {  
2325                if(trasmesso==0) {  
2326                //  std:: cout << "STRIP " << klt << endl;  std:: cout << " Entering WriteTrackCalib " << endl;
2327                //  std:: cout << "ADC_TOT " <<AdcTrack[iv][ladder*fNstrips_ladder+klt] << endl;  
2328    fOutputfile.write(reinterpret_cast<char*>(fDataPSCU),sizeof(UShort_t)*fPSCUbuffer);
2329                  DataDSP[DSPlength]=( ((UShort_t)klt) | 0x1000);  
2330                  DSPlength++;    UShort_t temp[1000000];
2331                  ntrastot++;  memset(temp,0,sizeof(UShort_t)*1000000);
2332                swab(fDataTrack,temp,sizeof(UShort_t)*fTracklength);  // WE MUST SWAP THE BYTES!!!
2333    fOutputfile.write(reinterpret_cast<char*>(temp),sizeof(UShort_t)*fTracklength);
2334                  real=modff(AdcTrack[iv][ladder*fNstrips_ladder+klt],&inte);  fTracklength=0;
2335                  if (real > 0.5) inte=inte+1;  if ( fPadding ){
2336                  DataDSP[DSPlength]=( ((UShort_t)inte) & 0x0FFF);        fOutputfile.write(reinterpret_cast<char*>(fDataPadding),sizeof(UChar_t)*fPadding);
2337                  DSPlength++;  };
2338                  ntrastot++;  
2339                };
2340                 }  
2341                 else {  
2342                 //   std:: cout << "ADC_TOT " <<AdcTrack[iv][ladder*fNstrips_ladder+klt] << endl;  void Digitizer::ClearTrackCalib() {
2343                  real=modff(AdcTrack[iv][ladder*fNstrips_ladder+klt],&inte);  
2344                  if (real > 0.5) inte=inte+1;  std:: cout << "Entering ClearTrackCalib " << endl;
2345                  DataDSP[DSPlength]=( ((UShort_t)inte) & 0x0FFF);      
2346                  DSPlength++;                  
2347                  ntrastot++;  };
2348                  };  
2349                  trasmesso=1;                          
2350              };  // end trasmission  void Digitizer::LoadTrackCalib() {
2351              kt=kl2;  std:: cout << "Entering LoadTrackCalib " << endl;
2352              k=kl2;  
2353              real=modff(AdcTrack[iv][ladder*fNstrips_ladder+kt],&inte);  // Generate the pedestals and sigmas according to parametrization
2354              if (real > 0.5) inte=inte+1;    for (Int_t j=0; j<fNviews;j++) {
2355              oldval=(Int_t)inte -(Int_t)fPedeTrack[iv][ladder*fNstrips_ladder+kt];      for (Int_t i=0; i<fNstrips_view;i++) {
2356              k++;      
2357              continue;      if((j+1)%2==0) {
2358        fPedeTrack[j][i]=gRandom->Gaus(fAvePedex,fSigmaPedex);
2359        fSigmaTrack[j][i]=gRandom->Gaus(fAveSigmax,fSigmaSigmax);
2360            }; // end cercacluster      };
2361          }; // end cercacluster      if((j+1)%2==1) {
2362              fPedeTrack[j][i]=gRandom->Gaus(fAvePedey,fSigmaPedey);
2363  // start ZOP check for strips no      fSigmaTrack[j][i]=gRandom->Gaus(fAveSigmay,fSigmaSigmay);
2364              };
2365        if(abs(newval-oldval)>=fCutzop*(Int_t)fSigmaTrack[iv][ladder*fNstrips_ladder+k]) {      
2366        };
2367         if(trasmesso==0) {    };    
2368         // std:: cout << "STRIP " << k << endl;  
2369         // std:: cout << "ADC_TOT " << AdcTrack[iv][ladder*fNstrips_ladder+k] << endl;    
2370      
2371           DataDSP[DSPlength]=( ((UShort_t)k) | 0x1000);  };
2372           DSPlength++;    
2373           ntrastot++;  void Digitizer::LoadMipCor() {
2374                std:: cout << "Entering LoadMipCor" << endl;
2375      Float_t xfactor=1./151.6*1.04;
2376           real=modff(AdcTrack[iv][ladder*fNstrips_ladder+k],&inte);    Float_t yfactor=1./152.1;
2377           if (real > 0.5) inte=inte+1;  
2378           DataDSP[DSPlength]=( ((UShort_t)inte) & 0x0FFF);    fMipCor[0][0]=140.02*yfactor;
2379           DSPlength++;    fMipCor[0][1]=140.99*xfactor;
2380           ntrastot++;    fMipCor[0][2]=134.48*yfactor;
2381      fMipCor[0][3]=144.41*xfactor;
2382         }    fMipCor[0][4]=140.74*yfactor;
2383         else {    fMipCor[0][5]=142.28*xfactor;
2384         //  std:: cout << "ADC_TOT " << AdcTrack[iv][ladder*fNstrips_ladder+k] << endl;    fMipCor[0][6]=134.53*yfactor;
2385           real=modff(AdcTrack[iv][ladder*fNstrips_ladder+k],&inte);    fMipCor[0][7]=140.63*xfactor;
2386           if (real > 0.5) inte=inte+1;    fMipCor[0][8]=135.55*yfactor;
2387           DataDSP[DSPlength]=(  ((UShort_t)inte) & 0x0FFF);    fMipCor[0][9]=138.00*xfactor;
2388           DSPlength++;    fMipCor[0][10]=154.95*yfactor;
2389           ntrastot++;          fMipCor[0][11]=158.44*xfactor;
2390         };    
2391         trasmesso=1;    
2392         oldval=newval;    fMipCor[1][0]=136.07*yfactor;
2393         kt=k;    fMipCor[1][1]=135.59*xfactor;
2394      fMipCor[1][2]=142.69*yfactor;
2395        }      fMipCor[1][3]=138.19*xfactor;
2396        else trasmesso=0;    fMipCor[1][4]=137.35*yfactor;
2397        // end zop    fMipCor[1][5]=140.23*xfactor;
2398            fMipCor[1][6]=153.15*yfactor;
2399        k++;          fMipCor[1][7]=151.42*xfactor;
2400        };  // end cycle inside ladder    fMipCor[1][8]=129.76*yfactor;
2401  // write here the end ladder bytes    fMipCor[1][9]=140.63*xfactor;
2402  //            std:: cout << "FINE LADDER " << ladder+1 << endl;    fMipCor[1][10]=157.87*yfactor;
2403      fMipCor[1][11]=153.64*xfactor;
2404        DataDSP[DSPlength]=( ((UShort_t)(ladder+1))  | 0x1800);  
2405        DSPlength++;    fMipCor[2][0]=134.98*yfactor;
2406        ntrastot++;    fMipCor[2][1]=143.95*xfactor;
2407        trasmesso=0;    fMipCor[2][2]=140.23*yfactor;
2408      fMipCor[2][3]=138.88*xfactor;
2409      };  //end cycle inside dsp    fMipCor[2][4]=137.95*yfactor;
2410  //  std:: cout << "FINE DSP " << iv+1 << endl;    fMipCor[2][5]=134.87*xfactor;
2411  // here put DSP header    fMipCor[2][6]=157.56*yfactor;
2412      DataDSP[0]=(0x1CA0 | ((UShort_t)(iv+1)) );    fMipCor[2][7]=157.31*xfactor;
2413      UShort_t Nword=(DSPlength*13)/16;    fMipCor[2][8]=141.37*yfactor;
2414      if( ((DSPlength*13)%16)!=0) Nword++;    fMipCor[2][9]=143.39*xfactor;
2415      DataDSP[1]=(0x1400 | ( Nword >> 10));    fMipCor[2][10]=156.15*yfactor;
2416      DataDSP[2]=(0x1400 | ( Nword & 0x03FF) );    fMipCor[2][11]=158.79*xfactor;
2417      DataDSP[3]=(0x1400 | (( (UShort_t)(fCounter >> 10) ) & 0x03FF) );  
2418      DataDSP[4]=(0x1400 | (( (UShort_t)(fCounter) )  & 0x03FF) );  /*
2419      DataDSP[5]=(0x1400 | ( (UShort_t)(fNclst << 7) ) | ( (UShort_t)(fCutzop << 4) )    for (Int_t j=0; j<fNviews;j++) {
2420       | ( (UShort_t)fCutzop  ) );      for (Int_t i=0; i<fNstrips_view;i++) {
2421      DataDSP[6]=0x1400;      fMipCor[j][i]=1.;
2422      DataDSP[7]=0x1400;      };
2423      DataDSP[8]=0x1400;    };    
2424      DataDSP[9]=0x1400;  
2425      DataDSP[10]=0x1400;      
2426      DataDSP[11]=0x1400;  */  
2427      DataDSP[12]=0x1400;  };
2428      DataDSP[13]=0x1400;  
2429      DataDSP[14]=(0x1400 | (CheckSum & 0x00FF) );  void Digitizer::CompressTrackData(Float_t AdcTrack[fNviews][fNstrips_view]) {
2430      DataDSP[15]=0x1C00;  // copy of the corresponding compression fortran routine + new digitization
2431  // end DSP header      // std:: cout << "Entering CompressTrackData " << endl;
2432    Int_t oldval=0;
2433    Int_t newval=0;
2434  // write 13 bit DataDSP bufer inside 16 bit fDataTrack buffer  Int_t trasmesso=0;
2435      Int_t Bit16free=16;  Int_t ntrastot=0;
2436      UShort_t Dato;  Float_t real;
2437      for (Int_t NDSP=0; NDSP<DSPlength;NDSP++) {  Float_t inte;
2438        Int_t Bit13ToWrite=13;  Int_t cercacluster=0;
2439        while(Bit13ToWrite>0) {  Int_t kt=0;
2440          if(Bit13ToWrite<=Bit16free) {  static const int DSPbufferSize = 4000; // 13 bit buffer to be rearranged in 16 bit Track buffer
2441            Dato=((DataDSP[NDSP]&(0xFFFF >> (16-Bit13ToWrite)))<<(Bit16free-Bit13ToWrite));  UShort_t DataDSP[DSPbufferSize];  // 13 bit  buffer to be rearranged in 16 bit Track buffer
2442            fDataTrack[fTracklength]=fDataTrack[fTracklength] | Dato ;  UShort_t DSPlength;  // 13 bit buffer to be rearranged in 16 bit Track buffer
2443            Bit16free=Bit16free-Bit13ToWrite;  
2444            Bit13ToWrite=0;  memset(fDataTrack,0,sizeof(UShort_t)*fTRACKbuffer); // probably not necessary becouse already done ?
2445            if(Bit16free==0) {  fTracklength=0;
2446              if(NDSP>15) CheckSum=CheckSum^fDataTrack[fTracklength];  
2447              fTracklength++;    for (Int_t iv=0; iv<fNviews;iv++) {
2448              Bit16free=16;    memset(DataDSP,0,sizeof(UShort_t)*DSPbufferSize);
2449            };              DSPlength=16; // skip the header, to be written later
2450          }    UShort_t CheckSum=0;
2451          else if(Bit13ToWrite>Bit16free) {  // write dsp header on buffer
2452            Dato=( (DataDSP[NDSP]&(0xFFFF >> (16-Bit13ToWrite) ) ) >> (Bit13ToWrite-Bit16free) );  
2453            fDataTrack[fTracklength]=fDataTrack[fTracklength] | Dato ;  //    fDataTrack[fTracklength]=0xE805;
2454            if(NDSP>15) CheckSum=CheckSum^fDataTrack[fTracklength];  //    fTracklength++;
2455            fTracklength++;  
2456            Bit13ToWrite=Bit13ToWrite-Bit16free;  //    fDataTrack[fTracklength]=0x01A9;
2457            Bit16free=16;    //    fTracklength++;
2458          };  
2459            // end dsp header
2460        }; // end cycle while(Bit13ToWrite>0)  
2461                 //
2462      }; // end cycle DataDSP     // INIZIO VISTA IV - TAKE PROPER ACTION
2463      if(Bit16free!=16) { fTracklength++; CheckSum=CheckSum^fDataTrack[fTracklength]; };     //
2464      CheckSum=(CheckSum >> 8)^(CheckSum&0x00FF);  
2465      fDataTrack[fTracklength-Nword+11]=(0x0280 | (CheckSum >> 3));  
2466      fDataTrack[fTracklength-Nword+12]=(0x1C00 | (CheckSum << 13) );  
2467        for (Int_t ladder=0; ladder<fNladder;ladder++) {
2468  // end write 13 bit DataDSP bufer inside 16 bit fDataTrack buffer        Int_t k=0;
2469          while (k<fNstrips_ladder) {
2470  //write trailer on buffer          // compress write in buffer the current LADDER
2471      UShort_t ReLength=(UShort_t)((Nword+13)*2+3);          if ( k == 0)  {
2472      UShort_t OveCheckCode=0x0000;            real=modff(AdcTrack[iv][ladder*fNstrips_ladder+k],&inte);
2473              if (real > 0.5) inte=inte+1;
2474      fDataTrack[fTracklength]=0x0000;            newval=(Int_t)inte -(Int_t)fPedeTrack[iv][ladder*fNstrips_ladder+k];
2475      fTracklength++;            // first strip of ladder is transmitted
2476                // DC_TOT first " << AdcTrack[iv][ladder*fNstrips_ladder+k] << endl;
2477      fDataTrack[fTracklength]=(ReLength >> 8);            DataDSP[DSPlength]=( ((UShort_t)inte) & 0x0FFF);
2478      fTracklength++;            DSPlength++;  
2479                ntrastot++;
2480      fDataTrack[fTracklength]=( (ReLength << 8) | (OveCheckCode & 0x00FF) );            trasmesso=1;
2481      fTracklength++;              oldval=newval;
2482  // end trailer            kt=k;
2483  //    std:: cout  << "DSPlength  " <<DSPlength << endl;            k++;
2484  //    std:: cout << "Nword " << Nword  << endl;            continue;
2485  //    std:: cout <<  "ReLength " << ReLength << endl;          };
2486    };              real=modff(AdcTrack[iv][ladder*fNstrips_ladder+k],&inte);
2487  //    std:: cout << "ntrastot " << ntrastot << endl;              if (real > 0.5) inte=inte+1;
2488            newval=(Int_t)inte -(Int_t)(fPedeTrack[iv][ladder*fNstrips_ladder+k]);
2489  };          cercacluster=1; // ?????????
2490            if (cercacluster==1) {
2491            
2492  Float_t Digitizer::SaturationTrack(Float_t ADC) {   // controlla l'ordine di tutti queste strip ladder e DSP !!!!!!!      
2493    Float_t SatFact=1.;            Int_t diff=0;
2494    if(ADC<70.) { SatFact=80./ADC; };  
2495    if(ADC>3000.) { SatFact=3000./ADC; };            
2496    return SatFact;            switch ((iv+1)%2) {
2497  };            case 0: diff=newval-oldval;
2498              break;
2499              case 1: diff=oldval-newval;    
2500              break;
2501              };
2502    
2503              if (diff>fCutclu*(Int_t)fSigmaTrack[iv][ladder*fNstrips_ladder+k]) {
2504                Int_t clval=newval;
2505                Int_t klp=k;        // go on to search for maximum
2506                klp++;
2507    
2508                while(klp<fNstrips_ladder) {
2509                  real=modff(AdcTrack[iv][ladder*fNstrips_ladder+klp],&inte);
2510                  if (real > 0.5) inte=inte+1;
2511                  Int_t clvalp=(Int_t)inte -(Int_t)fPedeTrack[iv][ladder*fNstrips_ladder+klp];
2512                  if((iv+1)%2==0) {
2513    
2514                    if(clvalp>clval) {
2515                       clval=clvalp;
2516                       k=klp;}
2517                    else break; // max of cluster found
2518    
2519                  } else {  
2520    
2521                    if(clvalp<clval) {
2522                       clval=clvalp;
2523                       k=klp;}
2524                    else break; // max of cluster found
2525    
2526                  };
2527                  
2528                  klp++;
2529                };
2530    
2531                Int_t kl1=k-fNclst; // max of cluster (or end of ladder ?)
2532                trasmesso=0;
2533                if(kl1<0)  kl1=0;
2534    
2535                if(kt>=kl1) kl1=kt+1;
2536                if( (kt+1)==kl1 ) trasmesso=1;
2537    
2538    
2539                
2540                Int_t kl2=k+fNclst;
2541                if(kl2>=fNstrips_ladder) kl2=fNstrips_ladder-1;
2542    
2543                for(Int_t klt=kl1 ; klt<=kl2 ; klt++) {
2544                  if(trasmesso==0) {
2545                  //  std:: cout << "STRIP " << klt << endl;
2546                  //  std:: cout << "ADC_TOT " <<AdcTrack[iv][ladder*fNstrips_ladder+klt] << endl;
2547    
2548                    DataDSP[DSPlength]=( ((UShort_t)klt) | 0x1000);
2549                    DSPlength++;  
2550                    ntrastot++;
2551                
2552    
2553                    real=modff(AdcTrack[iv][ladder*fNstrips_ladder+klt],&inte);
2554                    if (real > 0.5) inte=inte+1;
2555                    DataDSP[DSPlength]=( ((UShort_t)inte) & 0x0FFF);
2556                    DSPlength++;
2557                    ntrastot++;
2558                
2559                   }
2560                   else {
2561                   //   std:: cout << "ADC_TOT " <<AdcTrack[iv][ladder*fNstrips_ladder+klt] << endl;
2562                    real=modff(AdcTrack[iv][ladder*fNstrips_ladder+klt],&inte);
2563                    if (real > 0.5) inte=inte+1;
2564                    DataDSP[DSPlength]=( ((UShort_t)inte) & 0x0FFF);
2565                    DSPlength++;              
2566                    ntrastot++;
2567                    };
2568                    trasmesso=1;                        
2569                };  // end trasmission
2570                kt=kl2;
2571                k=kl2;
2572                real=modff(AdcTrack[iv][ladder*fNstrips_ladder+kt],&inte);
2573                if (real > 0.5) inte=inte+1;
2574                oldval=(Int_t)inte -(Int_t)fPedeTrack[iv][ladder*fNstrips_ladder+kt];
2575                k++;
2576                continue;
2577    
2578    
2579              }; // end cercacluster
2580            }; // end cercacluster
2581          
2582    // start ZOP check for strips no
2583          
2584          if(abs(newval-oldval)>=fCutzop*(Int_t)fSigmaTrack[iv][ladder*fNstrips_ladder+k]) {
2585    
2586           if(trasmesso==0) {
2587           // std:: cout << "STRIP " << k << endl;
2588           // std:: cout << "ADC_TOT " << AdcTrack[iv][ladder*fNstrips_ladder+k] << endl;
2589    
2590             DataDSP[DSPlength]=( ((UShort_t)k) | 0x1000);
2591             DSPlength++;  
2592             ntrastot++;
2593                
2594    
2595             real=modff(AdcTrack[iv][ladder*fNstrips_ladder+k],&inte);
2596             if (real > 0.5) inte=inte+1;
2597             DataDSP[DSPlength]=( ((UShort_t)inte) & 0x0FFF);
2598             DSPlength++;
2599             ntrastot++;
2600    
2601           }
2602           else {
2603           //  std:: cout << "ADC_TOT " << AdcTrack[iv][ladder*fNstrips_ladder+k] << endl;
2604             real=modff(AdcTrack[iv][ladder*fNstrips_ladder+k],&inte);
2605             if (real > 0.5) inte=inte+1;
2606             DataDSP[DSPlength]=(  ((UShort_t)inte) & 0x0FFF);
2607             DSPlength++;
2608             ntrastot++;      
2609           };
2610           trasmesso=1;
2611           oldval=newval;
2612           kt=k;
2613    
2614          }  
2615          else trasmesso=0;
2616          // end zop
2617          
2618          k++;      
2619          };  // end cycle inside ladder
2620    // write here the end ladder bytes
2621    //            std:: cout << "FINE LADDER " << ladder+1 << endl;
2622    
2623          DataDSP[DSPlength]=( ((UShort_t)(ladder+1))  | 0x1800);
2624          DSPlength++;
2625          ntrastot++;
2626          trasmesso=0;
2627    
2628        };  //end cycle inside dsp
2629    //  std:: cout << "FINE DSP " << iv+1 << endl;
2630    // here put DSP header
2631        DataDSP[0]=(0x1CA0 | ((UShort_t)(iv+1)) );
2632        UShort_t Nword=(DSPlength*13)/16;
2633        if( ((DSPlength*13)%16)!=0) Nword++;
2634        DataDSP[1]=(0x1400 | ( Nword >> 10));
2635        DataDSP[2]=(0x1400 | ( Nword & 0x03FF) );
2636        DataDSP[3]=(0x1400 | (( (UShort_t)(fCounter >> 10) ) & 0x03FF) );
2637        DataDSP[4]=(0x1400 | (( (UShort_t)(fCounter) )  & 0x03FF) );
2638        DataDSP[5]=(0x1400 | ( (UShort_t)(fNclst << 7) ) | ( (UShort_t)(fCutzop << 4) )
2639         | ( (UShort_t)fCutzop  ) );
2640        DataDSP[6]=0x1400;
2641        DataDSP[7]=0x1400;
2642        DataDSP[8]=0x1400;
2643        DataDSP[9]=0x1400;
2644        DataDSP[10]=0x1400;
2645        DataDSP[11]=0x1400;
2646        DataDSP[12]=0x1400;
2647        DataDSP[13]=0x1400;
2648        DataDSP[14]=(0x1400 | (CheckSum & 0x00FF) );
2649        DataDSP[15]=0x1C00;
2650    // end DSP header    
2651    
2652    
2653    // write 13 bit DataDSP bufer inside 16 bit fDataTrack buffer
2654        Int_t Bit16free=16;
2655        UShort_t Dato;
2656        for (Int_t NDSP=0; NDSP<DSPlength;NDSP++) {
2657          Int_t Bit13ToWrite=13;
2658          while(Bit13ToWrite>0) {
2659            if(Bit13ToWrite<=Bit16free) {
2660              Dato=((DataDSP[NDSP]&(0xFFFF >> (16-Bit13ToWrite)))<<(Bit16free-Bit13ToWrite));
2661              fDataTrack[fTracklength]=fDataTrack[fTracklength] | Dato ;
2662              Bit16free=Bit16free-Bit13ToWrite;
2663              Bit13ToWrite=0;
2664              if(Bit16free==0) {
2665                if(NDSP>15) CheckSum=CheckSum^fDataTrack[fTracklength];
2666                fTracklength++;
2667                Bit16free=16;
2668              };          
2669            }
2670            else if(Bit13ToWrite>Bit16free) {
2671              Dato=( (DataDSP[NDSP]&(0xFFFF >> (16-Bit13ToWrite) ) ) >> (Bit13ToWrite-Bit16free) );
2672              fDataTrack[fTracklength]=fDataTrack[fTracklength] | Dato ;
2673              if(NDSP>15) CheckSum=CheckSum^fDataTrack[fTracklength];
2674              fTracklength++;
2675              Bit13ToWrite=Bit13ToWrite-Bit16free;
2676              Bit16free=16;  
2677            };
2678            
2679          }; // end cycle while(Bit13ToWrite>0)
2680              
2681        }; // end cycle DataDSP
2682        if(Bit16free!=16) { fTracklength++; CheckSum=CheckSum^fDataTrack[fTracklength]; };
2683        CheckSum=(CheckSum >> 8)^(CheckSum&0x00FF);
2684        fDataTrack[fTracklength-Nword+11]=(0x0280 | (CheckSum >> 3));
2685        fDataTrack[fTracklength-Nword+12]=(0x1C00 | (CheckSum << 13) );
2686    
2687    // end write 13 bit DataDSP bufer inside 16 bit fDataTrack buffer
2688    
2689    //write trailer on buffer
2690        UShort_t ReLength=(UShort_t)((Nword+13)*2+3);
2691        UShort_t OveCheckCode=0x0000;
2692    
2693        fDataTrack[fTracklength]=0x0000;
2694        fTracklength++;
2695      
2696        fDataTrack[fTracklength]=(ReLength >> 8);
2697        fTracklength++;
2698      
2699        fDataTrack[fTracklength]=( (ReLength << 8) | (OveCheckCode & 0x00FF) );
2700        fTracklength++;  
2701    // end trailer
2702    //    std:: cout  << "DSPlength  " <<DSPlength << endl;
2703    //    std:: cout << "Nword " << Nword  << endl;
2704    //    std:: cout <<  "ReLength " << ReLength << endl;
2705      };    
2706    //    std:: cout << "ntrastot " << ntrastot << endl;    
2707    
2708    };
2709    
2710    
2711    Float_t Digitizer::SaturationTrack(Float_t ADC) {
2712      Float_t SatFact=1.;
2713      if(ADC<70.) { SatFact=80./ADC; };
2714      if(ADC>3000.) { SatFact=3000./ADC; };
2715      return SatFact;
2716    };
2717    
2718    
2719    
2720    
2721    
2722    

Legend:
Removed from v.1.2  
changed lines
  Added in v.1.5

  ViewVC Help
Powered by ViewVC 1.1.23