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

Contents of /PamelaDigitizer/Digitizer.cxx

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.5 - (show annotations) (download)
Wed Nov 28 18:54:31 2007 UTC (17 years, 3 months ago) by silvio
Branch: MAIN
Changes since 1.4: +102 -85 lines
Various ToF improvements

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

  ViewVC Help
Powered by ViewVC 1.1.23