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