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

Contents of /PamelaDigitizer/Digitizer.cxx

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1.1.1 - (show annotations) (download) (vendor branch)
Thu Sep 13 11:00:53 2007 UTC (17 years, 2 months ago) by silvio
Branch: PamelaDigitizer/
CVS Tags: v1r00, start
Changes since 1.1: +0 -0 lines
First release

1 //
2 #include <sstream>
3 #include <fstream>
4 #include <stdlib.h>
5 #include <string.h>
6 #include <ctype.h>
7 #include "Riostream.h"
8 #include "TFile.h"
9 #include "TDirectory.h"
10 #include "TTree.h"
11 #include "TLeafI.h"
12 #include "TH1.h"
13 #include "TH2.h"
14 #include "TMath.h"
15 #include "TRandom.h"
16 #include "TSQLServer.h"
17 #include "TSystem.h"
18 //
19 #include "Digitizer.h"
20 #include "CRC.h"
21 //
22 #include <PamelaRun.h>
23 #include <physics/calorimeter/CalorimeterEvent.h>
24 #include <CalibCalPedEvent.h>
25 #include "GLTables.h"
26 //
27 extern "C"{
28 short crc(short, short);
29 };
30 //
31
32 Digitizer::Digitizer(TTree* tree, char* &file_raw){
33 fhBookTree = tree;
34 fFilename = file_raw;
35 fCounter = 0;
36 fOBT = 0;
37
38 //
39 // DB connections
40 //
41 TString host = "mysql://localhost/pamelaprod";
42 TString user = "anonymous";
43 TString psw = "";
44 //
45 const char *pamdbhost=gSystem->Getenv("PAM_DBHOST");
46 const char *pamdbuser=gSystem->Getenv("PAM_DBUSER");
47 const char *pamdbpsw=gSystem->Getenv("PAM_DBPSW");
48 if ( !pamdbhost ) pamdbhost = "";
49 if ( !pamdbuser ) pamdbuser = "";
50 if ( !pamdbpsw ) pamdbpsw = "";
51 if ( strcmp(pamdbhost,"") ) host = pamdbhost;
52 if ( strcmp(pamdbuser,"") ) user = pamdbuser;
53 if ( strcmp(pamdbpsw,"") ) psw = pamdbpsw;
54 fDbc = TSQLServer::Connect(host.Data(),user.Data(),psw.Data());
55 //
56 GL_TABLES *glt = new GL_TABLES(host,user,psw);
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());
58 //
59 // Use UTC in the DB and make timeout bigger
60 //
61 stringstream myquery;
62 myquery.str("");
63 myquery << "SET time_zone='+0:00'";
64 fDbc->Query(myquery.str().c_str());
65 myquery.str("");
66 myquery << "SET wait_timeout=173000;";
67 fDbc->Query(myquery.str().c_str());
68 //
69
70 std:: cout << "preparing tree" << endl;
71
72 // prepare tree
73 fhBookTree->SetBranchAddress("Irun",&Irun);
74 fhBookTree->SetBranchAddress("Ievnt",&Ievnt);
75 fhBookTree->SetBranchAddress("Ipa",&Ipa);
76 fhBookTree->SetBranchAddress("X0",&X0);
77 fhBookTree->SetBranchAddress("Y0",&Y0);
78 fhBookTree->SetBranchAddress("Z0",&Z0);
79 fhBookTree->SetBranchAddress("Theta",&Theta);
80 fhBookTree->SetBranchAddress("Phi",&Phi);
81 fhBookTree->SetBranchAddress("P0",&P0);
82 fhBookTree->SetBranchAddress("Nthtof",&Nthtof);
83 fhBookTree->SetBranchAddress("Ipltof",Ipltof);
84 fhBookTree->SetBranchAddress("Ipaddle",Ipaddle);
85 fhBookTree->SetBranchAddress("Ipartof",Ipartof);
86 fhBookTree->SetBranchAddress("Xintof",Xintof);
87 fhBookTree->SetBranchAddress("Yintof",Yintof);
88 fhBookTree->SetBranchAddress("Zintof",Zintof);
89 fhBookTree->SetBranchAddress("Xouttof",Xouttof);
90 fhBookTree->SetBranchAddress("Youttof",Youttof);
91 fhBookTree->SetBranchAddress("Zouttof",Zouttof);
92 fhBookTree->SetBranchAddress("Ereltof",Ereltof);
93 fhBookTree->SetBranchAddress("Timetof",Timetof);
94 fhBookTree->SetBranchAddress("Pathtof",Pathtof);
95 fhBookTree->SetBranchAddress("P0tof",P0tof);
96 fhBookTree->SetBranchAddress("Nthcat",&Nthcat);
97 fhBookTree->SetBranchAddress("Iparcat",Iparcat);
98 fhBookTree->SetBranchAddress("Icat",Icat);
99 fhBookTree->SetBranchAddress("Xincat",Xincat);
100 fhBookTree->SetBranchAddress("Yincat",Yincat);
101 fhBookTree->SetBranchAddress("Zincat",Zincat);
102 fhBookTree->SetBranchAddress("Xoutcat",Xoutcat);
103 fhBookTree->SetBranchAddress("Youtcat",Youtcat);
104 fhBookTree->SetBranchAddress("Zoutcat",Zoutcat);
105 fhBookTree->SetBranchAddress("Erelcat",Erelcat);
106 fhBookTree->SetBranchAddress("Timecat",Timecat);
107 fhBookTree->SetBranchAddress("Pathcat",Pathcat);
108 fhBookTree->SetBranchAddress("P0cat",P0cat);
109 fhBookTree->SetBranchAddress("Nthcas",&Nthcas);
110 fhBookTree->SetBranchAddress("Iparcas",Iparcas);
111 fhBookTree->SetBranchAddress("Icas",Icas);
112 fhBookTree->SetBranchAddress("Xincas",Xincas);
113 fhBookTree->SetBranchAddress("Yincas",Yincas);
114 fhBookTree->SetBranchAddress("Zincas",Zincas);
115 fhBookTree->SetBranchAddress("Xoutcas",Xoutcas);
116 fhBookTree->SetBranchAddress("Youtcas",Youtcas);
117 fhBookTree->SetBranchAddress("Zoutcas",Zoutcas);
118 fhBookTree->SetBranchAddress("Erelcas",Erelcas);
119 fhBookTree->SetBranchAddress("Timecas",Timecas);
120 fhBookTree->SetBranchAddress("Pathcas",Pathcas);
121 fhBookTree->SetBranchAddress("P0cas",P0cas);
122 fhBookTree->SetBranchAddress("Nthspe",&Nthspe);
123 fhBookTree->SetBranchAddress("Iparspe",Iparspe);
124 fhBookTree->SetBranchAddress("Itrpb",Itrpb);
125 fhBookTree->SetBranchAddress("Itrsl",Itrsl);
126 fhBookTree->SetBranchAddress("Itspa",Itspa);
127 fhBookTree->SetBranchAddress("Xinspe",Xinspe);
128 fhBookTree->SetBranchAddress("Yinspe",Yinspe);
129 fhBookTree->SetBranchAddress("Zinspe",Zinspe);
130 fhBookTree->SetBranchAddress("Xoutspe",Xoutspe);
131 fhBookTree->SetBranchAddress("Youtspe",Youtspe);
132 fhBookTree->SetBranchAddress("Zoutspe",Zoutspe);
133 fhBookTree->SetBranchAddress("Xavspe",Xavspe);
134 fhBookTree->SetBranchAddress("Yavspe",Yavspe);
135 fhBookTree->SetBranchAddress("Zavspe",Zavspe);
136 fhBookTree->SetBranchAddress("Erelspe",Erelspe);
137 fhBookTree->SetBranchAddress("Pathspe",Pathspe);
138 fhBookTree->SetBranchAddress("P0spe",P0spe);
139 fhBookTree->SetBranchAddress("Nxmult",Nxmult);
140 fhBookTree->SetBranchAddress("Nymult",Nymult);
141 fhBookTree->SetBranchAddress("Nstrpx",&Nstrpx);
142 fhBookTree->SetBranchAddress("Npstripx",Npstripx);
143 fhBookTree->SetBranchAddress("Ntstripx",Ntstripx);
144 fhBookTree->SetBranchAddress("Istripx",Istripx);
145 fhBookTree->SetBranchAddress("Qstripx",Qstripx);
146 fhBookTree->SetBranchAddress("Xstripx",Xstripx);
147 fhBookTree->SetBranchAddress("Nstrpy",&Nstrpy);
148 fhBookTree->SetBranchAddress("Npstripy",Npstripy);
149 fhBookTree->SetBranchAddress("Ntstripy",Ntstripy);
150 fhBookTree->SetBranchAddress("Istripy",Istripy);
151 fhBookTree->SetBranchAddress("Qstripy",Qstripy);
152 fhBookTree->SetBranchAddress("Ystripy",Ystripy);
153 fhBookTree->SetBranchAddress("Nthcali",&Nthcali);
154 fhBookTree->SetBranchAddress("Icaplane",Icaplane);
155 fhBookTree->SetBranchAddress("Icastrip",Icastrip);
156 fhBookTree->SetBranchAddress("Icamod",Icamod);
157 fhBookTree->SetBranchAddress("Enestrip",Enestrip);
158 fhBookTree->SetBranchAddress("Nthcal",&Nthcal);
159 fhBookTree->SetBranchAddress("Icapl",Icapl);
160 fhBookTree->SetBranchAddress("Icasi",Icasi);
161 fhBookTree->SetBranchAddress("Icast",Icast);
162 fhBookTree->SetBranchAddress("Xincal",Xincal);
163 fhBookTree->SetBranchAddress("Yincal",Yincal);
164 fhBookTree->SetBranchAddress("Zincal",Zincal);
165 fhBookTree->SetBranchAddress("Erelcal",Erelcal);
166 fhBookTree->SetBranchAddress("Nthnd",&Nthnd);
167 fhBookTree->SetBranchAddress("Itubend",Itubend);
168 fhBookTree->SetBranchAddress("Iparnd",Iparnd);
169 fhBookTree->SetBranchAddress("Xinnd",Xinnd);
170 fhBookTree->SetBranchAddress("Yinnd",Yinnd);
171 fhBookTree->SetBranchAddress("Zinnd",Zinnd);
172 fhBookTree->SetBranchAddress("Xoutnd",Xoutnd);
173 fhBookTree->SetBranchAddress("Youtnd",Youtnd);
174 fhBookTree->SetBranchAddress("Zoutnd",Zoutnd);
175 fhBookTree->SetBranchAddress("Erelnd",Erelnd);
176 fhBookTree->SetBranchAddress("Timend",Timend);
177 fhBookTree->SetBranchAddress("Pathnd",Pathnd);
178 fhBookTree->SetBranchAddress("P0nd",P0nd);
179 fhBookTree->SetBranchAddress("Nthcard",&Nthcard);
180 fhBookTree->SetBranchAddress("Iparcard",Iparcard);
181 fhBookTree->SetBranchAddress("Icard",Icard);
182 fhBookTree->SetBranchAddress("Xincard",Xincard);
183 fhBookTree->SetBranchAddress("Yincard",Yincard);
184 fhBookTree->SetBranchAddress("Zincard",Zincard);
185 fhBookTree->SetBranchAddress("Xoutcard",Xoutcard);
186 fhBookTree->SetBranchAddress("Youtcard",Youtcard);
187 fhBookTree->SetBranchAddress("Zoutcard",Zoutcard);
188 fhBookTree->SetBranchAddress("Erelcard",Erelcard);
189 fhBookTree->SetBranchAddress("Timecard",Timecard);
190 fhBookTree->SetBranchAddress("Pathcard",Pathcard);
191 fhBookTree->SetBranchAddress("P0card",P0card);
192
193 fhBookTree->SetBranchStatus("*",0);
194
195 };
196
197
198
199 void Digitizer::Close(){
200
201 delete fhBookTree;
202
203 };
204
205
206
207
208 void Digitizer::Loop() {
209 //
210 // opens the raw output file and loops over the events
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();
218 DigitizeCALOCALIB();
219
220 // load, digitize and write tracker calibration
221 LoadTrackCalib();
222
223 DigitizeTrackCalib(1);
224 UInt_t length=fTracklength*2;
225 DigitizePSCU(length,0x12);
226 AddPadding();
227 WriteTrackCalib();
228
229 DigitizeTrackCalib(2);
230 length=fTracklength*2;
231 DigitizePSCU(length,0x13);
232 AddPadding();
233 WriteTrackCalib();
234
235 LoadMipCor(); // some initialization of parameters -not used now-
236 // end loading, digitizing and writing tracker calibration
237
238 //
239 // loops over the events
240 //
241
242 Int_t nentries = fhBookTree->GetEntriesFast();
243 Long64_t nbytes = 0;
244 for (Int_t i=0; i<nentries;i++) {
245 //
246 nbytes += fhBookTree->GetEntry(i);
247 // read detectors sequentially:
248 // http://www.ts.infn.it/fileadmin/documents/physics/experiments/wizard/cpu/gen_arch/RM_Acquisition.pdf
249 // on pamelatov:
250 // /cvs/yoda/techmodel/physics/NeutronDetectorReader.cpp
251 DigitizeTRIGGER();
252 DigitizeTOF();
253 DigitizeAC();
254 DigitizeCALO();
255 DigitizeTrack();
256 //DigitizeS4();
257 DigitizeND();
258 //
259 // Create CPU header, we need packet type (0x10 = physics data) and packet length.
260 //
261 UInt_t length = (fCALOlength + fACbuffer + fTracklength)*2;
262 DigitizePSCU(length,0x10);
263 //
264 // Add padding to 64 bits
265 //
266 AddPadding();
267 //
268 if ( !i%100 ) std::cout << "writing event " << i << endl;
269 WriteData();
270 };
271
272 fOutputfile.close();
273 std::cout << "files closed" << endl << flush;
274
275 };
276
277 void Digitizer::AddPadding(){
278 //
279 Float_t pd0 = (fLen+16)/64.;
280 Float_t pd1 = pd0 - (Float_t)int(pd0);
281 Float_t padfrac = 64. - pd1 * 64.;
282 //
283 UInt_t padbytes = (UInt_t)padfrac;
284 if ( padbytes > 0 && padbytes < 64 ){
285 //
286 // here the padding length
287 //
288 fPadding = padbytes+64;
289 //
290 // random padding bytes
291 //
292 for (Int_t ur=0; ur<32; ur++){
293 fDataPadding[ur] = (UShort_t)rand();
294 };
295 };
296 };
297
298
299 void Digitizer::DigitizePSCU(UInt_t length, UChar_t type) {
300 //
301 UChar_t buff[16];
302 //
303 // CPU signature
304 //
305 buff[0] = 0xFA;
306 buff[1] = 0xFE;
307 buff[2] = 0xDE;
308 //
309 // packet type (twice)
310 //
311 buff[3] = type;
312 buff[4] = type;
313 //
314 // counter
315 //
316 fCounter++;
317 while ( fCounter > 16777215 ){
318 fCounter -= 16777215;
319 };
320 //
321 buff[5] = (UChar_t)(fCounter >> 16);
322 buff[6] = (UChar_t)(fCounter >> 8);
323 buff[7] = (UChar_t)fCounter;
324 //
325 // on board time
326 //
327 ULong64_t obt = fOBT + 30LL;
328 //
329 while ( obt > 4294967295LL ){
330 obt -= 4294967295LL;
331 };
332 fOBT = (UInt_t)obt;
333 //
334 buff[8] = (UChar_t)(fOBT >> 24);
335 buff[9] = (UChar_t)(fOBT >> 16);
336 buff[10] = (UChar_t)(fOBT >> 8);
337 buff[11] = (UChar_t)fOBT;
338 //
339 // Packet length
340 //
341 fLen = length;
342 //
343 buff[12] = (UChar_t)(fLen >> 16);
344 buff[13] = (UChar_t)(fLen >> 8);
345 buff[14] = (UChar_t)fLen;
346 //
347 // CPU header CRC
348 //
349 buff[15] = (BYTE)CM_Compute_CRC16((UINT16)0, (BYTE*)&buff, (UINT32)15);
350 //
351 memcpy(fDataPSCU,buff,16*sizeof(UChar_t));
352 //
353 };
354
355 void Digitizer::ClearCaloCalib(Int_t s){
356 //
357 fcstwerr[s] = 0;
358 fcperror[s] = 0.;
359 for ( Int_t d=0 ; d<11 ;d++ ){
360 Int_t pre = -1;
361 for ( Int_t j=0; j<96 ;j++){
362 if ( j%16 == 0 ) pre++;
363 fcalped[s][d][j] = 0.;
364 fcstwerr[s] = 0.;
365 fcperror[s] = 0.;
366 fcalgood[s][d][j] = 0.;
367 fcalthr[s][d][pre] = 0.;
368 fcalrms[s][d][j] = 0.;
369 fcalbase[s][d][pre] = 0.;
370 fcalvar[s][d][pre] = 0.;
371 };
372 };
373 return;
374 }
375
376 Int_t Digitizer::CaloLoadCalib(Int_t s,TString fcalname, UInt_t calibno){
377 //
378 //
379 UInt_t e = 0;
380 if ( s == 0 ) e = 0;
381 if ( s == 1 ) e = 2;
382 if ( s == 2 ) e = 3;
383 if ( s == 3 ) e = 1;
384 //
385 ifstream myfile;
386 myfile.open(fcalname.Data());
387 if ( !myfile ){
388 return(-107);
389 };
390 myfile.close();
391 //
392 TFile *File = new TFile(fcalname.Data());
393 if ( !File ) return(-108);
394 TTree *tr = (TTree*)File->Get("CalibCalPed");
395 if ( !tr ) return(-109);
396 //
397 TBranch *calo = tr->GetBranch("CalibCalPed");
398 //
399 pamela::CalibCalPedEvent *ce = 0;
400 tr->SetBranchAddress("CalibCalPed", &ce);
401 //
402 Long64_t ncalibs = calo->GetEntries();
403 //
404 if ( !ncalibs ) return(-110);
405 //
406 calo->GetEntry(calibno);
407 //
408 if (ce->cstwerr[s] != 0 && ce->cperror[s] == 0 ) {
409 fcstwerr[s] = ce->cstwerr[s];
410 fcperror[s] = ce->cperror[s];
411 for ( Int_t d=0 ; d<11 ;d++ ){
412 Int_t pre = -1;
413 for ( Int_t j=0; j<96 ;j++){
414 if ( j%16 == 0 ) pre++;
415 fcalped[s][d][j] = ce->calped[e][d][j];
416 fcalgood[s][d][j] = ce->calgood[e][d][j];
417 fcalthr[s][d][pre] = ce->calthr[e][d][pre];
418 fcalrms[s][d][j] = ce->calrms[e][d][j];
419 fcalbase[s][d][pre] = ce->calbase[e][d][pre];
420 fcalvar[s][d][pre] = ce->calvar[e][d][pre];
421 };
422 };
423 } else {
424 printf(" CALORIMETER - ERROR: problems finding a good calibration in this file! \n\n ");
425 File->Close();
426 return(-111);
427 };
428 File->Close();
429 return(0);
430 }
431
432
433 void Digitizer::DigitizeCALOCALIB() {
434 //
435 // Header of the four sections
436 //
437 fSecCalo[0] = 0xAA00; // XE
438 fSecCalo[1] = 0xB100; // XO
439 fSecCalo[2] = 0xB600; // YE
440 fSecCalo[3] = 0xAD00; // YO
441 //
442 // length of the data is 0x1215
443 //
444 fSecCALOLength[0] = 0x1215; // XE
445 fSecCALOLength[1] = 0x1215; // XO
446 fSecCALOLength[2] = 0x1215; // YE
447 fSecCALOLength[3] = 0x1215; // YO
448 //
449 Int_t chksum = 0;
450 UInt_t tstrip = 0;
451 UInt_t fSecPointer = 0;
452 //
453 for (Int_t sec=0; sec < 4; sec++){
454 //
455 // sec = 0 -> XE 1 -> XO 2-> YE 3 -> YO
456 //
457 fCALOlength = 0;
458 memset(fDataCALO,0,sizeof(UShort_t)*fCALObuffer);
459 fSecPointer = fCALOlength;
460 //
461 // First of all we have section header and packet length
462 //
463 fDataCALO[fCALOlength] = fSecCalo[sec];
464 fCALOlength++;
465 fDataCALO[fCALOlength] = fSecCALOLength[sec];
466 fCALOlength++;
467 //
468 // Section XO is read in the opposite direction respect to the others
469 //
470 chksum = 0;
471 //
472 for (Int_t plane=0; plane < 11; plane++){
473 //
474 if ( sec == 1 ) tstrip = fCALOlength + 96*2;
475 //
476 for (Int_t strip=0; strip < 96; strip++){
477 //
478 chksum += (Int_t)fcalped[sec][plane][strip];
479 //
480 // save value
481 //
482 if ( sec == 1 ){
483 tstrip -= 2;
484 fDataCALO[tstrip] = (Int_t)fcalped[sec][plane][strip];
485 fDataCALO[tstrip+1] = (Int_t)fcalgood[sec][plane][strip];
486 } else {
487 fDataCALO[fCALOlength] = (Int_t)fcalped[sec][plane][strip];
488 fDataCALO[fCALOlength+1] = (Int_t)fcalgood[sec][plane][strip];
489 };
490 fCALOlength +=2;
491 };
492 //
493 };
494 //
495 fDataCALO[fCALOlength] = (UShort_t)chksum;
496 fCALOlength++;
497 fDataCALO[fCALOlength] = 0;
498 fCALOlength++;
499 fDataCALO[fCALOlength] = (UShort_t)((Int_t)(chksum >> 16));
500 fCALOlength++;
501 //
502 // Section XO is read in the opposite direction respect to the others
503 //
504 chksum = 0;
505 //
506 for (Int_t plane=0; plane < 11; plane++){
507 //
508 if ( sec == 1 ) tstrip = fCALOlength+6*2;
509 //
510 for (Int_t strip=0; strip < 6; strip++){
511 //
512 chksum += (Int_t)fcalthr[sec][plane][strip];
513 //
514 // save value
515 //
516 if ( sec == 1 ){
517 tstrip -= 2;
518 fDataCALO[tstrip] = 0;
519 fDataCALO[tstrip+1] = (Int_t)fcalthr[sec][plane][strip];
520 } else {
521 fDataCALO[fCALOlength] = 0;
522 fDataCALO[fCALOlength+1] = (Int_t)fcalthr[sec][plane][strip];
523 };
524 fCALOlength +=2;
525 };
526 //
527 };
528 //
529 fDataCALO[fCALOlength] = 0;
530 fCALOlength++;
531 fDataCALO[fCALOlength] = (UShort_t)chksum;
532 fCALOlength++;
533 fDataCALO[fCALOlength] = 0;
534 fCALOlength++;
535 fDataCALO[fCALOlength] = (UShort_t)((Int_t)(chksum >> 16));
536 fCALOlength++;
537 //
538 // Section XO is read in the opposite direction respect to the others
539 //
540 for (Int_t plane=0; plane < 11; plane++){
541 //
542 if ( sec == 1 ) tstrip = fCALOlength+96*2;
543 //
544 for (Int_t strip=0; strip < 96; strip++){
545 //
546 // save value
547 //
548 if ( sec == 1 ){
549 tstrip -= 2;
550 fDataCALO[tstrip] = 0;
551 fDataCALO[tstrip+1] = (Int_t)fcalrms[sec][plane][strip];
552 } else {
553 fDataCALO[fCALOlength] = 0;
554 fDataCALO[fCALOlength+1] = (Int_t)fcalrms[sec][plane][strip];
555 };
556 fCALOlength += 2;
557 };
558 //
559 };
560 //
561 // Section XO is read in the opposite direction respect to the others
562 //
563 for (Int_t plane=0; plane < 11; plane++){
564 //
565 if ( sec == 1 ) tstrip = fCALOlength+6*4;
566 //
567 for (Int_t strip=0; strip < 6; strip++){
568 //
569 // save value
570 //
571 if ( sec == 1 ){
572 tstrip -= 4;
573 fDataCALO[tstrip] = 0;
574 fDataCALO[tstrip+1] = (Int_t)fcalbase[sec][plane][strip];
575 fDataCALO[tstrip+2] = 0;
576 fDataCALO[tstrip+3] = (Int_t)fcalvar[sec][plane][strip];
577 } else {
578 fDataCALO[fCALOlength] = 0;
579 fDataCALO[fCALOlength+1] = (Int_t)fcalbase[sec][plane][strip];
580 fDataCALO[fCALOlength+2] = 0;
581 fDataCALO[fCALOlength+3] = (Int_t)fcalvar[sec][plane][strip];
582 };
583 fCALOlength +=4;
584 };
585 //
586 };
587 //
588 //
589 // here we calculate and save the CRC
590 //
591 fDataCALO[fCALOlength] = 0;
592 fCALOlength++;
593 Short_t CRC = 0;
594 for (UInt_t i=0; i<(fCALOlength-fSecPointer); i++){
595 CRC=crc(CRC,fDataCALO[i+fSecPointer]);
596 };
597 fDataCALO[fCALOlength] = (UShort_t)CRC;
598 fCALOlength++;
599 //
600 UInt_t length=fCALOlength*2;
601 DigitizePSCU(length,0x18);
602 //
603 // Add padding to 64 bits
604 //
605 AddPadding();
606 //
607 fOutputfile.write(reinterpret_cast<char*>(fDataPSCU),sizeof(UShort_t)*fPSCUbuffer);
608 UShort_t temp[1000000];
609 memset(temp,0,sizeof(UShort_t)*1000000);
610 swab(fDataCALO,temp,sizeof(UShort_t)*fCALOlength); // WE MUST SWAP THE BYTES!!!
611 fOutputfile.write(reinterpret_cast<char*>(temp),sizeof(UShort_t)*fCALOlength);
612 //
613 // padding to 64 bytes
614 //
615 if ( fPadding ){
616 fOutputfile.write(reinterpret_cast<char*>(fDataPadding),sizeof(UChar_t)*fPadding);
617 };
618 //
619 //
620 };
621 //
622 };
623
624 void Digitizer::CaloLoadCalib() {
625 //
626 fGivenCaloCalib = 0; // ####@@@@ should be given as input par @@@@####
627 //
628 // first of all load the MIP to ADC conversion values
629 //
630 stringstream calfile;
631 Int_t error = 0;
632 GL_PARAM *glparam = new GL_PARAM();
633 //
634 // determine where I can find calorimeter ADC to MIP conversion file
635 //
636 error = 0;
637 error = glparam->Query_GL_PARAM(3,101,fDbc);
638 //
639 calfile.str("");
640 calfile << glparam->PATH.Data() << "/";
641 calfile << glparam->NAME.Data();
642 //
643 printf("\n Using Calorimeter ADC to MIP conversion file: \n %s \n",calfile.str().c_str());
644 FILE *f;
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++ ){
650 for (Int_t k = 0; k < 22; k++ ){
651 for (Int_t l = 0; l < 96; l++ ){
652 fread(&fCalomip[m][k][l],sizeof(fCalomip[m][k][l]),1,f);
653 };
654 };
655 };
656 fclose(f);
657 //
658 // determine which calibration has to be used and load it for each section
659 //
660 GL_CALO_CALIB *glcalo = new GL_CALO_CALIB();
661 GL_ROOT *glroot = new GL_ROOT();
662 TString fcalname;
663 UInt_t idcalib;
664 UInt_t calibno;
665 UInt_t utime = 0;
666 //
667 for (UInt_t s=0; s<4; s++){
668 //
669 // clear calo calib variables for section s
670 //
671 ClearCaloCalib(s);
672 //
673 if ( fGivenCaloCalib ){
674 //
675 // a time has been given as input on the command line so retrieve the calibration that preceed that time
676 //
677 glcalo->Query_GL_CALO_CALIB(fGivenCaloCalib,utime,s,fDbc);
678 //
679 calibno = glcalo->EV_ROOT;
680 idcalib = glcalo->ID_ROOT_L0;
681 //
682 // determine path and name and entry of the calibration file
683 //
684 printf("\n");
685 printf(" ** SECTION %i **\n",s);
686 //
687 glroot->Query_GL_ROOT(idcalib,fDbc);
688 //
689 stringstream name;
690 name.str("");
691 name << glroot->PATH.Data() << "/";
692 name << glroot->NAME.Data();
693 //
694 fcalname = (TString)name.str().c_str();
695 //
696 printf("\n Section %i : using file %s calibration at entry %i: \n",s,fcalname.Data(),calibno);
697 //
698 } else {
699 error = 0;
700 error = glparam->Query_GL_PARAM(1,104,fDbc);
701 //
702 calfile.str("");
703 calfile << glparam->PATH.Data() << "/";
704 calfile << glparam->NAME.Data();
705 //
706 printf("\n Section %i : using default calorimeter calibration: \n %s \n",s,calfile.str().c_str());
707 //
708 fcalname = (TString)calfile.str().c_str();
709 calibno = s;
710 //
711 };
712 //
713 // load calibration variables in memory
714 //
715 CaloLoadCalib(s,fcalname,calibno);
716 //
717 };
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
720 //
721 delete glparam;
722 delete glcalo;
723 delete glroot;
724 };
725
726 void Digitizer::DigitizeCALO() {
727 //
728 fModCalo = 0; // 0 is RAW, 1 is COMPRESS, 2 is FULL ####@@@@ should be given as input par @@@@####
729 //
730 //
731 //
732 fCALOlength = 0; // reset total dimension of calo data
733 //
734 // gpamela variables to be used
735 //
736 fhBookTree->SetBranchStatus("Nthcali",1);
737 fhBookTree->SetBranchStatus("Icaplane",1);
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 //
744 switch ( fModCalo ){
745 case 0:
746 this->DigitizeCALORAW();
747 break;
748 case 1:
749 this->DigitizeCALOCOMPRESS();
750 break;
751 case 2:
752 this->DigitizeCALOFULL();
753 break;
754 };
755 //
756 };
757
758 Float_t Digitizer::GetCALOen(Int_t sec, Int_t plane, Int_t strip){
759 //
760 // determine plane and strip
761 //
762 Int_t mplane = 0;
763 //
764 // wrong!
765 //
766 // if ( sec == 0 || sec == 3 ) mplane = (plane * 4) + sec + 1;
767 // if ( sec == 1 ) mplane = (plane * 4) + 2 + 1;
768 // if ( sec == 2 ) mplane = (plane * 4) + 1 + 1;
769 //
770 if ( sec == 0 ) mplane = plane * 4 + 1; // it must be 0, 4, 8, ... (+1) from plane = 0, 11
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 //
779 for (Int_t i=0; i<Nthcali;i++){
780 if ( Icaplane[i] == mplane && Icastrip[i] == mstrip ){
781 return (Enestrip[i]);
782 };
783 };
784 //
785 // if not found it means no energy release so return 0.
786 //
787 return(0.);
788 };
789
790 void Digitizer::DigitizeCALORAW() {
791 //
792 // some variables
793 //
794 Float_t ens = 0.;
795 UInt_t adcsig = 0;
796 UInt_t adcbase = 0;
797 UInt_t adc = 0;
798 Int_t pre = 0;
799 UInt_t l = 0;
800 UInt_t lpl = 0;
801 UInt_t tstrip = 0;
802 UInt_t fSecPointer = 0;
803 Double_t pedenoise;
804 Float_t rms = 0.;
805 Float_t pedestal = 0.;
806 //
807 // clean the data structure
808 //
809 memset(fDataCALO,0,sizeof(UShort_t)*fCALObuffer);
810 //
811 // Header of the four sections
812 //
813 fSecCalo[0] = 0xEA08; // XE
814 fSecCalo[1] = 0xF108; // XO
815 fSecCalo[2] = 0xF608; // YE
816 fSecCalo[3] = 0xED08; // YO
817 //
818 // length of the data is 0x0428 in RAW mode
819 //
820 fSecCALOLength[0] = 0x0428; // XE
821 fSecCALOLength[1] = 0x0428; // XO
822 fSecCALOLength[2] = 0x0428; // YE
823 fSecCALOLength[3] = 0x0428; // YO
824 //
825 // let's start
826 //
827 fCALOlength = 0;
828 //
829 for (Int_t sec=0; sec < 4; sec++){
830 //
831 // sec = 0 -> XE 1 -> XO 2-> YE 3 -> YO
832 //
833 l = 0; // XE and XO are Y planes
834 if ( sec < 2 ) l = 1; // while YE and YO are X planes
835 //
836 fSecPointer = fCALOlength;
837 //
838 // First of all we have section header and packet length
839 //
840 fDataCALO[fCALOlength] = fSecCalo[sec];
841 fCALOlength++;
842 fDataCALO[fCALOlength] = fSecCALOLength[sec];
843 fCALOlength++;
844 //
845 // selftrigger coincidences - in the future we should add here some code to simulate timing response of pre-amplifiers
846 //
847 for (Int_t autoplane=0; autoplane < 7; autoplane++){
848 fDataCALO[fCALOlength] = 0x0000;
849 fCALOlength++;
850 };
851 //
852 //
853 // here comes data
854 //
855 //
856 // Section XO is read in the opposite direction respect to the others
857 //
858 if ( sec == 1 ){
859 tstrip = 96*11 + fCALOlength;
860 } else {
861 tstrip = 0;
862 };
863 //
864 pre = -1;
865 //
866 for (Int_t strip=0; strip < 96; strip++){
867 //
868 // which is the pre for this strip?
869 //
870 if (strip%16 == 0) {
871 pre++;
872 };
873 //
874 if ( sec == 1 ) tstrip -= 11;
875 //
876 for (Int_t plane=0; plane < 11; plane++){
877 //
878 // here is wrong!!!!
879 //
880 //
881 // if ( plane%2 == 0 && sec%2 != 0){
882 // lpl = plane*2;
883 // } else {
884 // lpl = (plane*2) + 1;
885 // };
886 //
887 if ( sec == 0 || sec == 3 ) lpl = plane * 2;
888 if ( sec == 1 || sec == 2 ) lpl = (plane * 2) + 1;
889 //
890 // get the energy in GeV from the simulation for that strip
891 //
892 ens = this->GetCALOen(sec,plane,strip);
893 //
894 // convert it into ADC channels
895 //
896 adcsig = int(ens*fCalomip[l][lpl][strip]/fCALOGeV2MIPratio);
897 //
898 // sum baselines
899 //
900 adcbase = (UInt_t)fcalbase[sec][plane][pre];
901 //
902 // add noise and pedestals
903 //
904 pedestal = fcalped[sec][plane][strip];
905 rms = fcalrms[sec][plane][strip]/4.;
906 //
907 // Add random gaussian noise of RMS rms and Centered in the pedestal
908 //
909 pedenoise = gRandom->Gaus((Double_t)pedestal,(Double_t)rms);
910 //
911 // Sum all contribution
912 //
913 adc = adcsig + adcbase + (Int_t)round(pedenoise);
914 //
915 // Signal saturation
916 //
917 if ( adc > 0x7FFF ) adc = 0x7FFF;
918 //
919 // save value
920 //
921 if ( sec == 1 ){
922 fDataCALO[tstrip] = adc;
923 tstrip++;
924 } else {
925 fDataCALO[fCALOlength] = adc;
926 };
927 fCALOlength++;
928 //
929 };
930 //
931 if ( sec == 1 ) tstrip -= 11;
932 //
933 };
934 //
935 // here we calculate and save the CRC
936 //
937 Short_t CRC = 0;
938 for (UInt_t i=0; i<(fCALOlength-fSecPointer); i++){
939 CRC=crc(CRC,fDataCALO[i+fSecPointer]);
940 };
941 fDataCALO[fCALOlength] = (UShort_t)CRC;
942 fCALOlength++;
943 //
944 };
945 //
946 // for (Int_t i=0; i<fCALOlength; i++){
947 // printf(" WORD %i DIGIT %0x \n",i,fDataCALO[i]);
948 // };
949 //
950 };
951
952 void Digitizer::DigitizeCALOCOMPRESS() {
953 //
954 printf(" COMPRESS MODE STILL NOT IMPLEMENTED! \n");
955 //
956 this->DigitizeCALORAW();
957 return;
958 //
959 //
960 //
961 fSecCalo[0] = 0xEA00;
962 fSecCalo[1] = 0xF100;
963 fSecCalo[2] = 0xF600;
964 fSecCalo[3] = 0xED00;
965 //
966 // length of the data in DSP mode must be calculated on fly during digitization
967 //
968 memset(fSecCALOLength,0x0,4*sizeof(UShort_t));
969 //
970 // here comes raw data
971 //
972 Int_t en = 0;
973 //
974 for (Int_t sec=0; sec < 4; sec++){
975 fDataCALO[en] = fSecCalo[sec];
976 en++;
977 fDataCALO[en] = fSecCALOLength[sec];
978 en++;
979 for (Int_t plane=0; plane < 11; plane++){
980 for (Int_t strip=0; strip < 11; strip++){
981 fDataCALO[en] = 0x0;
982 en++;
983 };
984 };
985 };
986 //
987 };
988
989 void Digitizer::DigitizeCALOFULL() {
990 //
991 printf(" FULL MODE STILL NOT IMPLEMENTED! \n");
992 //
993 this->DigitizeCALORAW();
994 return;
995 //
996 fSecCalo[0] = 0xEA00;
997 fSecCalo[1] = 0xF100;
998 fSecCalo[2] = 0xF600;
999 fSecCalo[3] = 0xED00;
1000 //
1001 // length of the data in DSP mode must be calculated on fly during digitization
1002 //
1003 memset(fSecCALOLength,0x0,4*sizeof(UShort_t));
1004 //
1005 // here comes raw data
1006 //
1007 Int_t en = 0;
1008 //
1009 for (Int_t sec=0; sec < 4; sec++){
1010 fDataCALO[en] = fSecCalo[sec];
1011 en++;
1012 fDataCALO[en] = fSecCALOLength[sec];
1013 en++;
1014 for (Int_t plane=0; plane < 11; plane++){
1015 for (Int_t strip=0; strip < 11; strip++){
1016 fDataCALO[en] = 0x0;
1017 en++;
1018 };
1019 };
1020 };
1021 //
1022 };
1023
1024 void Digitizer::DigitizeTRIGGER() {
1025 //fDataTrigger: 153 bytes
1026 for (Int_t j=0; j < 153; j++)
1027 fDataTrigger[0]=0x00;
1028 };
1029
1030 Int_t Digitizer::DigitizeTOF() {
1031 //fDataTof: 12 x 23 bytes (=276 bytes)
1032 UChar_t *pTof=fDataTof;
1033
1034 // --- activate branches:
1035 fhBookTree->SetBranchStatus("Nthtof",1);
1036 fhBookTree->SetBranchStatus("Ipltof",1);
1037 fhBookTree->SetBranchStatus("Ipaddle",1);
1038 fhBookTree->SetBranchStatus("Xintof",1);
1039 fhBookTree->SetBranchStatus("Yintof",1);
1040 fhBookTree->SetBranchStatus("Xouttof",1);
1041 fhBookTree->SetBranchStatus("Youttof",1);
1042 fhBookTree->SetBranchStatus("Ereltof",1);
1043 fhBookTree->SetBranchStatus("Timetof",1);
1044 // not yet used: Zintof, Xouttof, Youttof, Zouttof
1045
1046 // ------ evaluate energy in each pmt: ------
1047 // strip geometry (lenght/width)
1048 Float_t dimel[6] = {33.0, 40.8 ,18.0, 15.0, 15.0, 18.0};
1049 Float_t dimes[6] = {5.1, 5.5, 7.5, 9.0, 6.0, 5.0};
1050
1051 // S11 8 paddles 33.0 x 5.1 cm
1052 // S12 6 paddles 40.8 x 5.5 cm
1053 // S21 2 paddles 18.0 x 7.5 cm
1054 // S22 2 paddles 15.0 x 9.0 cm
1055 // S31 3 paddles 15.0 x 6.0 cm
1056 // S32 3 paddles 18.0 x 5.0 cm
1057
1058 // distance from the interaction point to the pmts (right,left)
1059 Float_t xpath[2]={0., 0.}; /*path(cm) in X per S12,S21,S32 verso il pmt DX o SX*/
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 */
1062
1063 const Float_t Pho_keV = 10.; // photons per keV in scintillator
1064 const Float_t echarge = 1.6e-19; // carica dell'elettrone
1065 Float_t Npho=0.;
1066 Float_t QevePmt_pC[48];
1067 Float_t QhitPad_pC[2]={0., 0.};
1068 Float_t QhitPmt_pC[2]={0., 0.};
1069 Float_t pmGain = 3.5e6; /* Gain: per il momento uguale per tutti */
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)
1072 Float_t ADCoffset=30.;
1073 Int_t ADClast=4095; // no signal --> ADC ch=4095
1074 Int_t ADCtof[48];
1075 //Float_t ADCsat=3100; ci pensiamo in futuro !
1076 //Float_t pCsat=2500;
1077 for(Int_t i=0; i<48; i++){
1078 QevePmt_pC[i] = 0;
1079 ADCtof[i]=0;
1080 }
1081
1082 // ------ read calibration file (get A1, A2, lambda1, lambda2)
1083 ifstream fileTriggerCalib;
1084 TString ftrigname="TrigCalibParam.txt";
1085 fileTriggerCalib.open(ftrigname.Data());
1086 if ( !fileTriggerCalib ) {
1087 printf("debug: no trigger calib file!\n");
1088 return(-117); //check output!
1089 };
1090 Float_t atte1[48],atte2[48],lambda1[48],lambda2[48];
1091 Int_t temp=0;
1092 for(Int_t i=0; i<48; i++){
1093 fileTriggerCalib >> temp;
1094 fileTriggerCalib >> atte1[i];
1095 fileTriggerCalib >> atte2[i];
1096 fileTriggerCalib >> lambda1[i];
1097 fileTriggerCalib >> lambda2[i];
1098 fileTriggerCalib >> temp;
1099 }
1100 fileTriggerCalib.close();
1101
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 */
1107
1108 //const Int_t nmax=??; = Nthtof
1109 Int_t nh, ip, ipad, ipmt;
1110 Int_t pmtleft=0, pmtright=0;
1111 Int_t *pl, *pr;
1112 pl = &pmtleft;
1113 pr = &pmtright;
1114
1115 /* ********************************** inizio loop sugli hit */
1116
1117 for(Int_t nh=0; nh<Nthtof; nh++){
1118
1119 for(Int_t j=0; j<2; j++) { // already done!! remove???
1120 xpath[j]=0.;
1121 ypath[j]=0.;
1122 FGeo[j]=0.;
1123 }
1124
1125 // ridefiniz. piano e pad per i vettori in C
1126 ip = Ipltof[nh]-1;
1127 ipad = Ipaddle[nh]-1;
1128 pmtleft=0;
1129 pmtright=0;
1130
1131 //Paddle2Pmt((Int_t)ip, (Int_t) ipad, (Int_t*) &pmtleft, (Int_t*) &pmtright);
1132 Paddle2Pmt(ip, ipad, &pmtleft, &pmtright);
1133 //Paddle2Pmt(ip, ipad, pl, pr);
1134
1135 // per avere anche la corrispondenza pmt --> half board e canale
1136 // metodo GetPMTIndex(Int_t ipmt, Int_t &hb, Int_t &ch) // non lo usiamo x ora
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)
1144 tpos = (Yintof[nh]+Youttof[nh])/2.;
1145 else
1146 if(ip==1 || ip==2 || ip==5) //--- Strip in X per S12,S21,S32
1147 tpos = (Xintof[nh]+Xouttof[nh])/2.;
1148 else if (ip!=6)
1149 printf("*** Warning: this option should never occur! (ip=%2i, nh=%2i)\n",ip,nh);
1150 path[0]= tpos + dimel[ip]/2.;
1151 path[1]= dimel[ip]/2.- tpos;
1152
1153 // cout <<"Strip N. ="<< ipaddle <<" piano n.= "<< iplane <<" POSIZ = "<< tpos <<"\n";
1154
1155 /* per il momento metto un fattore geometrico costante*/
1156 FGeo[0] =0.5;
1157 FGeo[1] =0.5;
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
1160
1161 /* rimando la fluttuazione poissoniana sui fotoni prodotti
1162 sto studiando come funziona la funzione:
1163 long int i = sto.Poisson(double x); */
1164 // Npho = Poisson(ERELTOF[nh])*Pho_keV*1e6 Eloss in GeV ?
1165 Npho = Ereltof[nh]*Pho_keV*10.0e6; // Eloss in GeV ?
1166
1167 Float_t knorm[2]={0., 0.}; // Donatella
1168 Float_t Atten[2]={0., 0.}; // Donatella
1169 for(Int_t j=0; j<2; j++){
1170 QhitPad_pC[j]= Npho*FGeo[j]*effi*pmGain*echarge;
1171 knorm[j]=QhitPad_pC[j]/(atte1[pmtleft+j]*exp((dimel[ip]/2.*pow(-1,j+1))/lambda1[pmtleft+j]) +
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 }
1179
1180 QevePmt_pC[pmtleft] += QhitPmt_pC[0];
1181 QevePmt_pC[pmtright] += QhitPmt_pC[1];
1182
1183 } // **************************************** fine loop sugli hit
1184
1185 for(Int_t i=0; i<48; i++){
1186 if(QevePmt_pC[i] != 0.){
1187 ADCtof[i]= (Int_t)(ADC_pC*QevePmt_pC[i] + ADCoffset);
1188 if(ADCtof[i]> ADClast) ADCtof[i]=ADClast;
1189 } else
1190 ADCtof[i]= ADClast;
1191 };
1192
1193 UChar_t tofBin;
1194 // --- write fDataTof:
1195 for (Int_t j=0; j < 12; j++){
1196 Int_t j12=j*12;
1197 fDataTof[j12+0]=0x00; // TDC_ID
1198 fDataTof[j12+1]=0x00; // EV_COUNT
1199 fDataTof[j12+2]=0x00; // TDC_MASK (1)
1200 fDataTof[j12+3]=0x00; // TDC_MASK (2)
1201 for (Int_t k=0; k < 4; k++){
1202 Int_t jk12=j12+k;
1203 tofBin=(UChar_t)(ADCtof[k+4*j]/256); // ADC# (msb) (#=1+k+4*j)
1204 fDataTof[jk12+4] = Bin2GrayTof(tofBin,fDataTof[jk12+4]);
1205 tofBin=(UChar_t)(ADCtof[k+4*j]%256); // ADC# (lsb)
1206 fDataTof[jk12+5] = Bin2GrayTof(tofBin,fDataTof[jk12+5]);
1207 fDataTof[jk12+6]=0x00; // TDC# (msb) -- Wolfgang
1208 fDataTof[jk12+7]=0x00; // TDC# (lsb) -- Wolfgang
1209 };
1210 fDataTof[j12+20]=0x00; // TEMP1
1211 fDataTof[j12+21]=0x00; // TEMP2
1212 fDataTof[j12+22]= EvaluateCrcTof(pTof); // CRC
1213 pTof+=23;
1214 };
1215 return(0);
1216 };
1217
1218 UChar_t Digitizer::Bin2GrayTof(UChar_t binaTOF,UChar_t grayTOF){
1219 union graytof_data {
1220 UChar_t word;
1221 struct bit_field {
1222 unsigned b0:1;
1223 unsigned b1:1;
1224 unsigned b2:1;
1225 unsigned b3:1;
1226 unsigned b4:1;
1227 unsigned b5:1;
1228 unsigned b6:1;
1229 unsigned b7:1;
1230 } bit;
1231 } bi,gr;
1232 //
1233 bi.word = binaTOF;
1234 gr.word = grayTOF;
1235 //
1236 gr.bit.b0 = bi.bit.b1 ^ bi.bit.b0;
1237 gr.bit.b1 = bi.bit.b2 ^ bi.bit.b1;
1238 gr.bit.b2 = bi.bit.b3 ^ bi.bit.b2;
1239 gr.bit.b3 = bi.bit.b3;
1240 //
1241 /* bin to gray conversion 4 bit per time*/
1242 //
1243 gr.bit.b4 = bi.bit.b5 ^ bi.bit.b4;
1244 gr.bit.b5 = bi.bit.b6 ^ bi.bit.b5;
1245 gr.bit.b6 = bi.bit.b7 ^ bi.bit.b6;
1246 gr.bit.b7 = bi.bit.b7;
1247 //
1248 return(gr.word);
1249 }
1250
1251 UChar_t Digitizer::EvaluateCrcTof(UChar_t *pTof) {
1252 // UChar_t crcTof=0x00;
1253 // for (Int_t jp=0; jp < 23; jp++){
1254 // crcTof = crc8(...)
1255 // }
1256 return(0x00);
1257 };
1258
1259 //void Digitizer::Paddle2Pmt(Int_t plane, Int_t paddle, Int_t* &pmtleft, Int_t* &pmtright){
1260 void Digitizer::Paddle2Pmt(Int_t plane, Int_t paddle, Int_t *pl, Int_t *pr){
1261 //* @param plane (0 - 5)
1262 //* @param paddle (plane=0, paddle = 0,...5)
1263 //* @param padid (0 - 23)
1264 //
1265 Int_t padid=-1;
1266 Int_t pads[6]={8,6,2,2,3,3};
1267 //
1268 Int_t somma=0;
1269 Int_t np=plane;
1270 for(Int_t j=0; j<np; j++)
1271 somma+=pads[j];
1272 padid=paddle+somma;
1273 *pl = padid*2;
1274 *pr = *pr + 1;
1275 };
1276
1277 void Digitizer::DigitizeAC() {
1278 // created: J. Conrad, KTH
1279 // modified: S. Orsi, INFN Roma2
1280
1281 fDataAC[0] = 0xACAC;
1282 fDataAC[64]= 0xACAC;
1283 fDataAC[1] = 0xAC11; // main card
1284 fDataAC[65] = 0xAC22; // extra card
1285
1286 // the third word is a status word (dummy)
1287 fDataAC[2] = 0xFFFF; //FFEF?
1288 fDataAC[66] = 0xFFFF;
1289
1290 const UInt_t nReg = 6;
1291
1292 // Registers (dummy)
1293 for (UInt_t i=0; i<=nReg; i++){
1294 fDataAC[i+4] = 0xFFFF;
1295 fDataAC[i+68] = 0xFFFF;
1296 }
1297
1298 // the last word is a CRC
1299 // Dummy for the time being, but it might need to be calculated in the end
1300 fDataAC[63] = 0xABCD;
1301 fDataAC[127] = 0xABCD;
1302
1303 // shift registers, which one is with respect to PMT, where in
1304 // shift registers is a question of time relative trigger
1305 // In level2: hitmap, hitmap-status (synchronised with a trigger),
1306 // status
1307
1308 for (UInt_t i=0; i<=15; i++){
1309 fDataAC[i+11] = 0x0000;
1310 fDataAC[i+75] = 0x0000;
1311 }
1312
1313 // singles counters are dummy
1314
1315 for (UInt_t i=0; i<=16; i++){
1316 fDataAC[i+26] = 0x0000;
1317 fDataAC[i+90] = 0x0000;
1318 }
1319
1320 // coincidences are dummy
1321
1322 for (UInt_t i=0; i<=7; i++){
1323 fDataAC[i+42] = 0x0000;
1324 fDataAC[i+106] = 0x0000;
1325 }
1326
1327 // increments for every trigger might be needed at some point.
1328 // dummy for now
1329 fDataAC[50] = 0x0000;
1330 fDataAC[114] = 0x0000;
1331
1332 // dummy FPGA clock
1333
1334 fDataAC[51] = 0x006C;
1335 fDataAC[52] = 0x6C6C;
1336 fDataAC[115] = 0x006C;
1337 fDataAC[116] = 0x6C6C;
1338
1339
1340 // dummy temperatures
1341 fDataAC[53] = 0x0000;
1342 fDataAC[54] = 0x0000;
1343 fDataAC[117] = 0x0000;
1344 fDataAC[118] = 0x0000;
1345
1346
1347 // dummy DAC thresholds
1348 for (UInt_t i=0; i<=7; i++){
1349 fDataAC[i+55] = 0x1A13;
1350 fDataAC[i+119] = 0x1A13;
1351 }
1352
1353 // We activate all branches. Once the digitization algorithm
1354 // is determined only the branches need to activated which involve needed
1355 // information
1356
1357 fhBookTree->SetBranchStatus("Nthcat",1);
1358 fhBookTree->SetBranchStatus("Iparcat",1);
1359 fhBookTree->SetBranchStatus("Icat",1);
1360 fhBookTree->SetBranchStatus("Xincat",1);
1361 fhBookTree->SetBranchStatus("Yincat",1);
1362 fhBookTree->SetBranchStatus("Zincat",1);
1363 fhBookTree->SetBranchStatus("Xoutcat",1);
1364 fhBookTree->SetBranchStatus("Youtcat",1);
1365 fhBookTree->SetBranchStatus("Zoutcat",1);
1366 fhBookTree->SetBranchStatus("Erelcat",1);
1367 fhBookTree->SetBranchStatus("Timecat",1);
1368 fhBookTree->SetBranchStatus("Pathcat",1);
1369 fhBookTree->SetBranchStatus("P0cat",1);
1370 fhBookTree->SetBranchStatus("Nthcas",1);
1371 fhBookTree->SetBranchStatus("Iparcas",1);
1372 fhBookTree->SetBranchStatus("Icas",1);
1373 fhBookTree->SetBranchStatus("Xincas",1);
1374 fhBookTree->SetBranchStatus("Yincas",1);
1375 fhBookTree->SetBranchStatus("Zincas",1);
1376 fhBookTree->SetBranchStatus("Xoutcas",1);
1377 fhBookTree->SetBranchStatus("Youtcas",1);
1378 fhBookTree->SetBranchStatus("Zoutcas",1);
1379 fhBookTree->SetBranchStatus("Erelcas",1);
1380 fhBookTree->SetBranchStatus("Timecas",1);
1381 fhBookTree->SetBranchStatus("Pathcas",1);
1382 fhBookTree->SetBranchStatus("P0cas",1);
1383 fhBookTree->SetBranchStatus("Nthcard",1);
1384 fhBookTree->SetBranchStatus("Iparcard",1);
1385 fhBookTree->SetBranchStatus("Icard",1);
1386 fhBookTree->SetBranchStatus("Xincard",1);
1387 fhBookTree->SetBranchStatus("Yincard",1);
1388 fhBookTree->SetBranchStatus("Zincard",1);
1389 fhBookTree->SetBranchStatus("Xoutcard",1);
1390 fhBookTree->SetBranchStatus("Youtcard",1);
1391 fhBookTree->SetBranchStatus("Zoutcard",1);
1392 fhBookTree->SetBranchStatus("Erelcard",1);
1393 fhBookTree->SetBranchStatus("Timecard",1);
1394 fhBookTree->SetBranchStatus("Pathcard",1);
1395 fhBookTree->SetBranchStatus("P0card",1);
1396
1397 // In this simpliefied approach we will assume that once
1398 // a particle releases > 0.5 mip in one of the 12 AC detectors it
1399 // will fire. We will furthermore assume that both cards read out
1400 // identical data.
1401
1402 // If you develop you digitization algorithm, you should start by
1403 // identifying the information present in level2 (post-darth-vader)
1404 // data.
1405
1406 Float_t SumEcat[5];
1407 Float_t SumEcas[5];
1408 Float_t SumEcard[5];
1409 for (Int_t k= 0;k<5;k++){
1410 SumEcat[k]=0.;
1411 SumEcas[k]=0.;
1412 SumEcard[k]=0.;
1413 };
1414
1415 if (Nthcat>50 || Nthcas>50 || Nthcard>50)
1416 printf("Error! NthAC out of range!\n\n");
1417
1418 // look in CAT
1419 // for (UInt_t k= 0;k<50;k++){
1420 for (Int_t k= 0;k<Nthcat;k++){
1421 if (Erelcat[k] > 0)
1422 SumEcat[Icat[k]] += Erelcat[k];
1423 };
1424
1425 // look in CAS
1426 for (Int_t k= 0;k<Nthcas;k++){
1427 if (Erelcas[k] >0)
1428 SumEcas[Icas[k]] += Erelcas[k];
1429 };
1430
1431 // look in CARD
1432 for (Int_t k= 0;k<Nthcard;k++){
1433 if (Erelcard[k] >0)
1434 SumEcard[Icard[k]] += Erelcard[k];
1435 };
1436
1437 // channel mapping Hit Map
1438 // 1 CARD4 0 LSB
1439 // 2 CAT2 0
1440 // 3 CAS1 0
1441 // 4 NC 0
1442 // 5 CARD2 0
1443 // 6 CAT4 1
1444 // 7 CAS4 0
1445 // 8 NC 0
1446 // 9 CARD3 0
1447 // 10 CAT3 0
1448 // 11 CAS3 0
1449 // 12 NC 0
1450 // 13 CARD1 0
1451 // 14 CAT1 0
1452 // 15 CAS2 0
1453 // 16 NC 0 MSB
1454
1455 // In the first version only the hit-map is filled, not the SR.
1456
1457 // Threshold: 0.8 MeV.
1458
1459 Float_t thr = 8e-4;
1460
1461 fDataAC[3] = 0x0000;
1462
1463 if (SumEcas[0] > thr) fDataAC[3] = 0x0004;
1464 if (SumEcas[1] > thr) fDataAC[3] += 0x4000;
1465 if (SumEcas[2] > thr) fDataAC[3] += 0x0400;
1466 if (SumEcas[3] > thr) fDataAC[3] += 0x0040;
1467
1468 if (SumEcat[0] > thr) fDataAC[3] += 0x2000;
1469 if (SumEcat[1] > thr) fDataAC[3] += 0x0002;
1470 if (SumEcat[2] > thr) fDataAC[3] += 0x0200;
1471 if (SumEcat[3] > thr) fDataAC[3] += 0x0020;
1472
1473 if (SumEcard[0] > thr) fDataAC[3] += 0x1000;
1474 if (SumEcard[1] > thr) fDataAC[3] += 0x0010;
1475 if (SumEcard[2] > thr) fDataAC[3] += 0x0100;
1476 if (SumEcard[3] > thr) fDataAC[3] += 0x0001;
1477
1478 fDataAC[67] = fDataAC[3];
1479
1480 // for (Int_t i=0; i<fACbuffer; i++){
1481 // printf("%0x ",fDataAC[i]);
1482 // if ((i+1)%8 ==0) cout << endl;
1483 // }
1484 };
1485
1486
1487 void Digitizer::DigitizeND(){
1488 // creato: S. Borisov, INFN Roma2 e MEPHI, Sept 2007
1489 // 4 bytes: 16bit header, 8bit trigPhysics, 16bit up&low background
1490
1491 // ND header
1492 fDataND[0] = 0x0000;
1493 fDataND[1] = 0x000F;
1494
1495 fhBookTree->SetBranchStatus("Nthnd",1);
1496 fhBookTree->SetBranchStatus("Itubend",1);
1497 fhBookTree->SetBranchStatus("Iparnd",1);
1498 fhBookTree->SetBranchStatus("Xinnd",1);
1499 fhBookTree->SetBranchStatus("Yinnd",1);
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);
1505 fhBookTree->SetBranchStatus("Timend",1);
1506 fhBookTree->SetBranchStatus("Pathnd",1);
1507 fhBookTree->SetBranchStatus("P0nd",1);
1508
1509 UShort_t NdN=0;
1510 for(Int_t i=0;i<Nthnd;i++)
1511 if(Iparnd[i]==13)
1512 NdN++;
1513
1514 NdN=10; // test!
1515 fDataND[2]=0x0F00 & (NdN*256);
1516 //fDataND[2]=0xFFFF; //test
1517 fDataND[2]=0x0000; //background neutrons
1518 }
1519
1520
1521 void Digitizer::DigitizeDummy() {
1522
1523 fhBookTree->SetBranchStatus("Enestrip",1);
1524
1525 // dumy header
1526 fDataDummy[0] = 0xCAAA;
1527
1528 for (Int_t i=1; i<fDummybuffer; i++){
1529 fDataDummy[i] = 0xFFFF;
1530 // printf("%0x ",fDataDummy[i]);
1531 //if ((i+1)%8 ==0) cout << endl;
1532 }
1533
1534
1535
1536 };
1537
1538
1539 void Digitizer::WriteData(){
1540
1541 // Routine that writes the data to a binary file
1542 // PSCU data are already swapped
1543 fOutputfile.write(reinterpret_cast<char*>(fDataPSCU),sizeof(UShort_t)*fPSCUbuffer);
1544 // TRG
1545 fOutputfile.write(reinterpret_cast<char*>(fDataTrigger),sizeof(UChar_t)*153);
1546 // TOF
1547 fOutputfile.write(reinterpret_cast<char*>(fDataTof),sizeof(UChar_t)*276);
1548 // AC
1549 UShort_t temp[1000000];
1550 memset(temp,0,sizeof(UShort_t)*1000000);
1551 swab(fDataAC,temp,sizeof(UShort_t)*fACbuffer); // WE MUST SWAP THE BYTES!!!
1552 fOutputfile.write(reinterpret_cast<char*>(temp),sizeof(UShort_t)*fACbuffer);
1553 // CALO
1554 memset(temp,0,sizeof(UShort_t)*1000000);
1555 swab(fDataCALO,temp,sizeof(UShort_t)*fCALOlength); // WE MUST SWAP THE BYTES!!!
1556 fOutputfile.write(reinterpret_cast<char*>(temp),sizeof(UShort_t)*fCALOlength);
1557 // TRK
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;
1562 // S4
1563 // ...to be done...
1564 // ND
1565 memset(temp,0,sizeof(UShort_t)*1000000);
1566 swab(fDataND,temp,sizeof(UShort_t)*4); // WE MUST SWAP THE BYTES!!!
1567 fOutputfile.write(reinterpret_cast<char*>(temp),sizeof(UShort_t)*4);
1568
1569 //
1570 // fOutputfile.write(reinterpret_cast<char*>(fDataDummy),sizeof(UShort_t)*fDummybuffer);
1571 //
1572 // padding to 64 bytes
1573 //
1574 if ( fPadding ){
1575 fOutputfile.write(reinterpret_cast<char*>(fDataPadding),sizeof(UChar_t)*fPadding);
1576 };
1577 //
1578 };
1579
1580
1581 void Digitizer::ReadData(){
1582
1583 UShort_t InData[64];
1584
1585 // for debuggigng purposes only, write your own routine if you like (many
1586 // hardwired things.
1587
1588 ifstream InputFile;
1589
1590 // if (!InputFile) {
1591
1592 // std::cout << "ERROR" << endl;
1593 // // An error occurred!
1594 // // myFile.gcount() returns the number of bytes read.
1595 // // calling myFile.clear() will reset the stream state
1596 // // so it is usable again.
1597 // };
1598
1599
1600
1601 //InputFile.seekg(0);
1602
1603 InputFile.open(fFilename, ios::in | ios::binary);
1604 // fOutputfile.seekg(0);
1605 if (!InputFile.is_open()) std::cout << "ERROR" << endl;
1606
1607 InputFile.seekg(0);
1608
1609 for (Int_t k=0; k<=1000; k++){
1610 InputFile.read(reinterpret_cast<char*>(InData),384*sizeof(UShort_t));
1611
1612 std::cout << "Read back: " << endl << endl;
1613
1614 for (Int_t i=0; i<=384; i++){
1615 printf("%4x ", InData[i]);
1616 if ((i+1)%8 ==0) cout << endl;
1617 }
1618
1619 }
1620 cout << endl;
1621 InputFile.close();
1622
1623 };
1624
1625
1626
1627 void Digitizer::DigitizeTrack() {
1628 //std:: cout << "Entering DigitizeTrack " << endl;
1629 Float_t AdcTrack[fNviews][fNstrips_view]; // Vector of strips to be compressed
1630
1631 Int_t Iview;
1632 Int_t Nstrip;
1633
1634 for (Int_t j=0; j<fNviews;j++) {
1635
1636 for (Int_t i=0; i<fNladder;i++) {
1637
1638 Float_t commonN1=gRandom->Gaus(0.,fSigmaCommon);
1639 Float_t commonN2=gRandom->Gaus(0.,fSigmaCommon);
1640 for (Int_t k=0; k<fNstrips_ladder;k++) {
1641 Nstrip=i*fNstrips_ladder+k;
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
1644 else {AdcTrack[j][Nstrip] += commonN2;} // full correlation of 4 VA1 Com. Noise
1645
1646 };
1647
1648
1649 };
1650
1651
1652 };
1653
1654
1655 fhBookTree->SetBranchStatus("Nstrpx",1);
1656 fhBookTree->SetBranchStatus("Npstripx",1);
1657 fhBookTree->SetBranchStatus("Ntstripx",1);
1658 fhBookTree->SetBranchStatus("Istripx",1);
1659 fhBookTree->SetBranchStatus("Qstripx",1);
1660 fhBookTree->SetBranchStatus("Xstripx",1);
1661 fhBookTree->SetBranchStatus("Nstrpy",1);
1662 fhBookTree->SetBranchStatus("Npstripy",1);
1663 fhBookTree->SetBranchStatus("Ntstripy",1);
1664 fhBookTree->SetBranchStatus("Istripy",1);
1665 fhBookTree->SetBranchStatus("Qstripy",1);
1666 fhBookTree->SetBranchStatus("Ystripy",1);
1667
1668
1669
1670
1671 Float_t ADCfull;
1672 for (Int_t ix=0; ix<Nstrpx;ix++) {
1673 Iview=Npstripx[ix]*2-1;
1674 Nstrip=(Int_t)Istripx[ix]-1;
1675 ADCfull=AdcTrack[Iview][Nstrip] += Qstripx[ix]*fMipCor;
1676 AdcTrack[Iview][Nstrip] *= SaturationTrack(ADCfull);
1677
1678 };
1679
1680
1681 for (Int_t iy=0; iy<Nstrpy;iy++) {
1682 Iview=Npstripy[iy]*2-2;
1683 Nstrip=(Int_t)Istripy[iy]-1;
1684 ADCfull=AdcTrack[Iview][Nstrip] -= Qstripy[iy]*fMipCor;
1685 AdcTrack[Iview][Nstrip] *= SaturationTrack(ADCfull);
1686
1687 };
1688
1689 CompressTrackData(AdcTrack); // Compress and Digitize data of one Ladder in turn for all ladders
1690
1691 };
1692
1693
1694
1695 void Digitizer::DigitizeTrackCalib(Int_t ii) {
1696
1697 std:: cout << "Entering DigitizeTrackCalib " << ii << endl;
1698 if( (ii!=1)&&(ii!=2) ) {
1699 std:: cout << "error wrong DigitizeTrackCalib argument" << endl;
1700 return;
1701 };
1702
1703 memset(fDataTrack,0,sizeof(UShort_t)*fTRACKbuffer);
1704 fTracklength=0;
1705
1706 UShort_t Dato;
1707
1708 Float_t dato1;
1709 Float_t dato2;
1710 Float_t dato3;
1711 Float_t dato4;
1712
1713 UShort_t DatoDec;
1714 UShort_t DatoDec1;
1715 UShort_t DatoDec2;
1716 UShort_t DatoDec3;
1717 UShort_t DatoDec4;
1718
1719 UShort_t EVENT_CAL;
1720 UShort_t PED_L1;
1721 UShort_t ReLength;
1722 UShort_t OveCheckCode;
1723 //UShort_t PED_L2;
1724 //UShort_t PED_L3HI;
1725 //UShort_t PED_L3LO;
1726 //UShort_t SIG_L1HI;
1727 //UShort_t SIG_L1LO;
1728 //UShort_t SIG_L2HI;
1729 //UShort_t SIG_L2LO;
1730 //UShort_t SIG_L3;
1731 //UShort_t BAD_L1;
1732 //UShort_t BAD_L2LO;
1733 //UShort_t BAD_L3HI;
1734 //UShort_t BAD_L3LO;
1735 //UShort_t FLAG;
1736
1737
1738 Int_t DSPpos;
1739 for (Int_t j=ii-1; j<fNviews;j+=2) {
1740 UShort_t CkSum=0;
1741 // here skip the dsp header and his trailer , to be written later
1742 DSPpos=fTracklength;
1743 fTracklength=fTracklength+13+3;
1744
1745
1746 for (Int_t i=0; i<fNladder;i++) {
1747 for (Int_t k=0; k<fNstrips_ladder;k++) {
1748 // write in buffer the current LADDER
1749 Dato=(UShort_t)fPedeTrack[j][i*fNstrips_ladder+k];
1750 dato1=fPedeTrack[j][i*fNstrips_ladder+k]-Dato;
1751
1752 DatoDec1=(UShort_t)(dato1*2);
1753 dato2=dato1*2-DatoDec1;
1754
1755 DatoDec2=(UShort_t)(dato2*2);
1756 dato3=dato2*2-DatoDec2;
1757
1758 DatoDec3=(UShort_t)(dato3*2);
1759 dato4=dato3*2-DatoDec3;
1760
1761 DatoDec4=(UShort_t)(dato4*2);
1762
1763 DatoDec=DatoDec1*0x0008+DatoDec2*0x0004+DatoDec3*0x0002+DatoDec4*0x0001;
1764 fDataTrack[fTracklength]=( (Dato << 4) | (DatoDec & 0x000F) );
1765 CkSum=CkSum^fDataTrack[fTracklength];
1766 fTracklength++;
1767 };
1768
1769 for (Int_t k=0; k<fNstrips_ladder;k++) {
1770 // write in buffer the current LADDER
1771 Dato=(UShort_t)fSigmaTrack[j][i*fNstrips_ladder+k];
1772 dato1=fSigmaTrack[j][i*fNstrips_ladder+k]-Dato;
1773
1774 DatoDec1=(UShort_t)(dato1*2);
1775 dato2=dato1*2-DatoDec1;
1776
1777 DatoDec2=(UShort_t)(dato2*2);
1778 dato3=dato2*2-DatoDec2;
1779
1780 DatoDec3=(UShort_t)(dato3*2);
1781 dato4=dato3*2-DatoDec3;
1782
1783 DatoDec4=(UShort_t)(dato4*2);
1784
1785 DatoDec=DatoDec1*0x0008+DatoDec2*0x0004+DatoDec3*0x0002+DatoDec4*0x0001;
1786
1787 fDataTrack[fTracklength]=( (Dato << 4) | (DatoDec & 0x000F) );
1788 CkSum=CkSum^fDataTrack[fTracklength];
1789 fTracklength++;
1790 };
1791
1792 for (Int_t k=0; k<64;k++) {
1793 fDataTrack[fTracklength]=0x0000;
1794 CkSum=CkSum^fDataTrack[fTracklength];
1795 fTracklength++;
1796
1797 };
1798 // end ladder
1799
1800 // write in buffer the end ladder word
1801 if(i==0) fDataTrack[fTracklength]=0x1807;
1802 if(i==1) fDataTrack[fTracklength]=0x1808;
1803 if(i==2) fDataTrack[fTracklength]=0x1809;
1804 CkSum=CkSum^fDataTrack[fTracklength];
1805 fTracklength++;
1806
1807 // write in buffer the TRAILER
1808 ReLength=(UShort_t)((fNstrips_ladder*2+64+1)*2+3);
1809 OveCheckCode=0x0000;
1810
1811 fDataTrack[fTracklength]=0x0000;
1812 fTracklength++;
1813
1814 fDataTrack[fTracklength]=(ReLength >> 8);
1815 fTracklength++;
1816
1817 fDataTrack[fTracklength]=( (ReLength << 8) | (OveCheckCode & 0x00FF) );
1818 fTracklength++;
1819
1820 // end TRAILER
1821 };
1822
1823 // write in buffer the DSP header
1824
1825 fDataTrack[DSPpos]=(0xE800 | ( ((j+1) << 3) | 0x0005) );
1826
1827 fDataTrack[DSPpos+1]=0x01A9;
1828
1829 fDataTrack[DSPpos+2]=0x8740;
1830
1831 EVENT_CAL=0;
1832 fDataTrack[DSPpos+3]=(0x1A00 | ( (0x03FF & EVENT_CAL)>> 1) );
1833
1834 PED_L1=0;
1835 fDataTrack[DSPpos+4]=( ((EVENT_CAL << 15) | 0x5002 ) | ((0x03FF & PED_L1) << 2) );
1836
1837 // FROM HERE WE WRITE AS ALL VARIABLE apart CkSum are =0
1838
1839 fDataTrack[DSPpos+5]=0x8014;
1840
1841 fDataTrack[DSPpos+6]=0x00A0;
1842
1843 fDataTrack[DSPpos+7]=0x0500;
1844
1845 fDataTrack[DSPpos+8]=0x2801;
1846
1847 fDataTrack[DSPpos+9]=0x400A;
1848
1849 fDataTrack[DSPpos+10]=0x0050;
1850
1851 CkSum=(CkSum >> 8)^(CkSum&0x00FF);
1852 fDataTrack[DSPpos+11]=(0x0280 | (CkSum >> 3));
1853
1854 fDataTrack[DSPpos+12]=(0x1FFF | (CkSum << 13) );
1855
1856 // end dsp header
1857
1858 // write in buffer the TRAILER
1859
1860 ReLength=(UShort_t)((13*2)+3);
1861 OveCheckCode=0x0000;
1862 fDataTrack[DSPpos+13]=0x0000;
1863
1864 fDataTrack[DSPpos+14]=(ReLength >> 8);
1865
1866 fDataTrack[DSPpos+15]=( (ReLength << 8) | (OveCheckCode & 0x00FF) );
1867
1868 // end TRAILER
1869
1870
1871
1872
1873 // end DSP
1874 };
1875
1876
1877
1878 };
1879
1880 void Digitizer::WriteTrackCalib() {
1881
1882
1883 std:: cout << " Entering WriteTrackCalib " << endl;
1884
1885 fOutputfile.write(reinterpret_cast<char*>(fDataPSCU),sizeof(UShort_t)*fPSCUbuffer);
1886
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!!!
1890 fOutputfile.write(reinterpret_cast<char*>(temp),sizeof(UShort_t)*fTracklength);
1891 fTracklength=0;
1892 if ( fPadding ){
1893 fOutputfile.write(reinterpret_cast<char*>(fDataPadding),sizeof(UChar_t)*fPadding);
1894 };
1895
1896 };
1897
1898
1899 void Digitizer::ClearTrackCalib() {
1900
1901 std:: cout << "Entering ClearTrackCalib " << endl;
1902
1903
1904 };
1905
1906
1907 void Digitizer::LoadTrackCalib() {
1908 std:: cout << "Entering LoadTrackCalib " << endl;
1909
1910 // Generate the pedestals and sigmas according to parametrization
1911 for (Int_t j=0; j<fNviews;j++) {
1912 for (Int_t i=0; i<fNstrips_view;i++) {
1913
1914 if((j+1)%2==0) {
1915 fPedeTrack[j][i]=gRandom->Gaus(fAvePedex,fSigmaPedex);
1916 fSigmaTrack[j][i]=gRandom->Gaus(fAveSigmax,fSigmaSigmax);
1917 };
1918 if((j+1)%2==1) {
1919 fPedeTrack[j][i]=gRandom->Gaus(fAvePedey,fSigmaPedey);
1920 fSigmaTrack[j][i]=gRandom->Gaus(fAveSigmay,fSigmaSigmay);
1921 };
1922
1923 };
1924 };
1925
1926
1927
1928 };
1929
1930 void Digitizer::LoadMipCor() {
1931 std:: cout << "Entering LoadMipCor" << endl;
1932 /*
1933 for (Int_t j=0; j<fNviews;j++) {
1934 for (Int_t i=0; i<fNstrips_view;i++) {
1935 fMipCor[j][i]=1.;
1936 };
1937 };
1938
1939
1940 */
1941 };
1942
1943 void Digitizer::CompressTrackData(Float_t AdcTrack[fNviews][fNstrips_view]) {
1944 // copy of the corresponding compression fortran routine + new digitization
1945 // std:: cout << "Entering CompressTrackData " << endl;
1946 Int_t oldval=0;
1947 Int_t newval=0;
1948 Int_t trasmesso=0;
1949 Int_t ntrastot=0;
1950 Float_t real;
1951 Float_t inte;
1952 Int_t cercacluster=0;
1953 Int_t kt=0;
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
1956 UShort_t DSPlength; // 13 bit buffer to be rearranged in 16 bit Track buffer
1957
1958 memset(fDataTrack,0,sizeof(UShort_t)*fTRACKbuffer); // probably not necessary becouse already done ?
1959 fTracklength=0;
1960
1961 for (Int_t iv=0; iv<fNviews;iv++) {
1962 memset(DataDSP,0,sizeof(UShort_t)*DSPbufferSize);
1963 DSPlength=16; // skip the header, to be written later
1964 UShort_t CheckSum=0;
1965 // write dsp header on buffer
1966
1967 // fDataTrack[fTracklength]=0xE805;
1968 // fTracklength++;
1969
1970 // fDataTrack[fTracklength]=0x01A9;
1971 // fTracklength++;
1972
1973 // end dsp header
1974
1975 //
1976 // INIZIO VISTA IV - TAKE PROPER ACTION
1977 //
1978
1979
1980
1981 for (Int_t ladder=0; ladder<fNladder;ladder++) {
1982 Int_t k=0;
1983 while (k<fNstrips_ladder) {
1984 // compress write in buffer the current LADDER
1985 if ( k == 0) {
1986 real=modff(AdcTrack[iv][ladder*fNstrips_ladder+k],&inte);
1987 if (real > 0.5) inte=inte+1;
1988 newval=(Int_t)inte -(Int_t)fPedeTrack[iv][ladder*fNstrips_ladder+k];
1989 // first strip of ladder is transmitted
1990 // DC_TOT first " << AdcTrack[iv][ladder*fNstrips_ladder+k] << endl;
1991 DataDSP[DSPlength]=( ((UShort_t)inte) & 0x0FFF);
1992 DSPlength++;
1993 ntrastot++;
1994 trasmesso=1;
1995 oldval=newval;
1996 kt=k;
1997 k++;
1998 continue;
1999 };
2000 real=modff(AdcTrack[iv][ladder*fNstrips_ladder+k],&inte);
2001 if (real > 0.5) inte=inte+1;
2002 newval=(Int_t)inte -(Int_t)(fPedeTrack[iv][ladder*fNstrips_ladder+k]);
2003 cercacluster=1; // ?????????
2004 if (cercacluster==1) {
2005
2006 // controlla l'ordine di tutti queste strip ladder e DSP !!!!!!!
2007 Int_t diff=0;
2008
2009
2010 switch ((iv+1)%2) {
2011 case 0: diff=newval-oldval;
2012 break;
2013 case 1: diff=oldval-newval;
2014 break;
2015 };
2016
2017 if (diff>fCutclu*(Int_t)fSigmaTrack[iv][ladder*fNstrips_ladder+k]) {
2018 Int_t clval=newval;
2019 Int_t klp=k; // go on to search for maximum
2020 klp++;
2021
2022 while(klp<fNstrips_ladder) {
2023 real=modff(AdcTrack[iv][ladder*fNstrips_ladder+klp],&inte);
2024 if (real > 0.5) inte=inte+1;
2025 Int_t clvalp=(Int_t)inte -(Int_t)fPedeTrack[iv][ladder*fNstrips_ladder+klp];
2026 if((iv+1)%2==0) {
2027
2028 if(clvalp>clval) {
2029 clval=clvalp;
2030 k=klp;}
2031 else break; // max of cluster found
2032
2033 } else {
2034
2035 if(clvalp<clval) {
2036 clval=clvalp;
2037 k=klp;}
2038 else break; // max of cluster found
2039
2040 };
2041
2042 klp++;
2043 };
2044
2045 Int_t kl1=k-fNclst; // max of cluster (or end of ladder ?)
2046 trasmesso=0;
2047 if(kl1<0) kl1=0;
2048
2049 if(kt>=kl1) kl1=kt+1;
2050 if( (kt+1)==kl1 ) trasmesso=1;
2051
2052
2053
2054 Int_t kl2=k+fNclst;
2055 if(kl2>=fNstrips_ladder) kl2=fNstrips_ladder-1;
2056
2057 for(Int_t klt=kl1 ; klt<=kl2 ; klt++) {
2058 if(trasmesso==0) {
2059 // std:: cout << "STRIP " << klt << endl;
2060 // std:: cout << "ADC_TOT " <<AdcTrack[iv][ladder*fNstrips_ladder+klt] << endl;
2061
2062 DataDSP[DSPlength]=( ((UShort_t)klt) | 0x1000);
2063 DSPlength++;
2064 ntrastot++;
2065
2066
2067 real=modff(AdcTrack[iv][ladder*fNstrips_ladder+klt],&inte);
2068 if (real > 0.5) inte=inte+1;
2069 DataDSP[DSPlength]=( ((UShort_t)inte) & 0x0FFF);
2070 DSPlength++;
2071 ntrastot++;
2072
2073 }
2074 else {
2075 // std:: cout << "ADC_TOT " <<AdcTrack[iv][ladder*fNstrips_ladder+klt] << endl;
2076 real=modff(AdcTrack[iv][ladder*fNstrips_ladder+klt],&inte);
2077 if (real > 0.5) inte=inte+1;
2078 DataDSP[DSPlength]=( ((UShort_t)inte) & 0x0FFF);
2079 DSPlength++;
2080 ntrastot++;
2081 };
2082 trasmesso=1;
2083 }; // end trasmission
2084 kt=kl2;
2085 k=kl2;
2086 real=modff(AdcTrack[iv][ladder*fNstrips_ladder+kt],&inte);
2087 if (real > 0.5) inte=inte+1;
2088 oldval=(Int_t)inte -(Int_t)fPedeTrack[iv][ladder*fNstrips_ladder+kt];
2089 k++;
2090 continue;
2091
2092
2093 }; // end cercacluster
2094 }; // end cercacluster
2095
2096 // start ZOP check for strips no
2097
2098 if(abs(newval-oldval)>=fCutzop*(Int_t)fSigmaTrack[iv][ladder*fNstrips_ladder+k]) {
2099
2100 if(trasmesso==0) {
2101 // std:: cout << "STRIP " << k << endl;
2102 // std:: cout << "ADC_TOT " << AdcTrack[iv][ladder*fNstrips_ladder+k] << endl;
2103
2104 DataDSP[DSPlength]=( ((UShort_t)k) | 0x1000);
2105 DSPlength++;
2106 ntrastot++;
2107
2108
2109 real=modff(AdcTrack[iv][ladder*fNstrips_ladder+k],&inte);
2110 if (real > 0.5) inte=inte+1;
2111 DataDSP[DSPlength]=( ((UShort_t)inte) & 0x0FFF);
2112 DSPlength++;
2113 ntrastot++;
2114
2115 }
2116 else {
2117 // std:: cout << "ADC_TOT " << AdcTrack[iv][ladder*fNstrips_ladder+k] << endl;
2118 real=modff(AdcTrack[iv][ladder*fNstrips_ladder+k],&inte);
2119 if (real > 0.5) inte=inte+1;
2120 DataDSP[DSPlength]=( ((UShort_t)inte) & 0x0FFF);
2121 DSPlength++;
2122 ntrastot++;
2123 };
2124 trasmesso=1;
2125 oldval=newval;
2126 kt=k;
2127
2128 }
2129 else trasmesso=0;
2130 // end zop
2131
2132 k++;
2133 }; // end cycle inside ladder
2134 // write here the end ladder bytes
2135 // std:: cout << "FINE LADDER " << ladder+1 << endl;
2136
2137 DataDSP[DSPlength]=( ((UShort_t)(ladder+1)) | 0x1800);
2138 DSPlength++;
2139 ntrastot++;
2140 trasmesso=0;
2141
2142 }; //end cycle inside dsp
2143 // std:: cout << "FINE DSP " << iv+1 << endl;
2144 // here put DSP header
2145 DataDSP[0]=(0x1CA0 | ((UShort_t)(iv+1)) );
2146 UShort_t Nword=(DSPlength*13)/16;
2147 if( ((DSPlength*13)%16)!=0) Nword++;
2148 DataDSP[1]=(0x1400 | ( Nword >> 10));
2149 DataDSP[2]=(0x1400 | ( Nword & 0x03FF) );
2150 DataDSP[3]=(0x1400 | (( (UShort_t)(fCounter >> 10) ) & 0x03FF) );
2151 DataDSP[4]=(0x1400 | (( (UShort_t)(fCounter) ) & 0x03FF) );
2152 DataDSP[5]=(0x1400 | ( (UShort_t)(fNclst << 7) ) | ( (UShort_t)(fCutzop << 4) )
2153 | ( (UShort_t)fCutzop ) );
2154 DataDSP[6]=0x1400;
2155 DataDSP[7]=0x1400;
2156 DataDSP[8]=0x1400;
2157 DataDSP[9]=0x1400;
2158 DataDSP[10]=0x1400;
2159 DataDSP[11]=0x1400;
2160 DataDSP[12]=0x1400;
2161 DataDSP[13]=0x1400;
2162 DataDSP[14]=(0x1400 | (CheckSum & 0x00FF) );
2163 DataDSP[15]=0x1C00;
2164 // end DSP header
2165
2166
2167 // write 13 bit DataDSP bufer inside 16 bit fDataTrack buffer
2168 Int_t Bit16free=16;
2169 UShort_t Dato;
2170 for (Int_t NDSP=0; NDSP<DSPlength;NDSP++) {
2171 Int_t Bit13ToWrite=13;
2172 while(Bit13ToWrite>0) {
2173 if(Bit13ToWrite<=Bit16free) {
2174 Dato=((DataDSP[NDSP]&(0xFFFF >> (16-Bit13ToWrite)))<<(Bit16free-Bit13ToWrite));
2175 fDataTrack[fTracklength]=fDataTrack[fTracklength] | Dato ;
2176 Bit16free=Bit16free-Bit13ToWrite;
2177 Bit13ToWrite=0;
2178 if(Bit16free==0) {
2179 if(NDSP>15) CheckSum=CheckSum^fDataTrack[fTracklength];
2180 fTracklength++;
2181 Bit16free=16;
2182 };
2183 }
2184 else if(Bit13ToWrite>Bit16free) {
2185 Dato=( (DataDSP[NDSP]&(0xFFFF >> (16-Bit13ToWrite) ) ) >> (Bit13ToWrite-Bit16free) );
2186 fDataTrack[fTracklength]=fDataTrack[fTracklength] | Dato ;
2187 if(NDSP>15) CheckSum=CheckSum^fDataTrack[fTracklength];
2188 fTracklength++;
2189 Bit13ToWrite=Bit13ToWrite-Bit16free;
2190 Bit16free=16;
2191 };
2192
2193 }; // end cycle while(Bit13ToWrite>0)
2194
2195 }; // end cycle DataDSP
2196 if(Bit16free!=16) { fTracklength++; CheckSum=CheckSum^fDataTrack[fTracklength]; };
2197 CheckSum=(CheckSum >> 8)^(CheckSum&0x00FF);
2198 fDataTrack[fTracklength-Nword+11]=(0x0280 | (CheckSum >> 3));
2199 fDataTrack[fTracklength-Nword+12]=(0x1C00 | (CheckSum << 13) );
2200
2201 // end write 13 bit DataDSP bufer inside 16 bit fDataTrack buffer
2202
2203 //write trailer on buffer
2204 UShort_t ReLength=(UShort_t)((Nword+13)*2+3);
2205 UShort_t OveCheckCode=0x0000;
2206
2207 fDataTrack[fTracklength]=0x0000;
2208 fTracklength++;
2209
2210 fDataTrack[fTracklength]=(ReLength >> 8);
2211 fTracklength++;
2212
2213 fDataTrack[fTracklength]=( (ReLength << 8) | (OveCheckCode & 0x00FF) );
2214 fTracklength++;
2215 // end trailer
2216 // std:: cout << "DSPlength " <<DSPlength << endl;
2217 // std:: cout << "Nword " << Nword << endl;
2218 // std:: cout << "ReLength " << ReLength << endl;
2219 };
2220 // std:: cout << "ntrastot " << ntrastot << endl;
2221
2222 };
2223
2224 Float_t Digitizer::SaturationTrack(Float_t ADC) {
2225
2226 Float_t SatFact=1.;
2227 return SatFact;
2228 };
2229
2230
2231
2232
2233
2234

  ViewVC Help
Powered by ViewVC 1.1.23