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

Contents of /PamelaDigitizer/DigitizeCalo.cxx

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.2 - (show annotations) (download)
Wed Oct 15 14:03:16 2008 UTC (16 years, 1 month ago) by pamelats
Branch: MAIN
CVS Tags: v3r04, v3r05, v3r03
Changes since 1.1: +0 -28 lines
Cambiamenti principali: TOF, AC; cambiamenti di struttura (tutti *h esterni si trovano in Digitizer.h)

1 #include "Digitizer.h"
2
3 extern "C"{
4 short crc(short, short);
5 };
6
7 void Digitizer::ClearCaloCalib(Int_t s){
8 //
9 fcstwerr[s] = 0;
10 fcperror[s] = 0.;
11 for ( Int_t d=0 ; d<11 ;d++ ){
12 Int_t pre = -1;
13 for ( Int_t j=0; j<96 ;j++){
14 if ( j%16 == 0 ) pre++;
15 fcalped[s][d][j] = 0.;
16 fcstwerr[s] = 0.;
17 fcperror[s] = 0.;
18 fcalgood[s][d][j] = 0.;
19 fcalthr[s][d][pre] = 0.;
20 fcalrms[s][d][j] = 0.;
21 fcalbase[s][d][pre] = 0.;
22 fcalvar[s][d][pre] = 0.;
23 };
24 };
25 return;
26 }
27
28 Int_t Digitizer::CaloLoadCalib(Int_t s,TString fcalname, UInt_t calibno){
29 //
30 //
31 UInt_t e = 0;
32 if ( s == 0 ) e = 0;
33 if ( s == 1 ) e = 2;
34 if ( s == 2 ) e = 3;
35 if ( s == 3 ) e = 1;
36 //
37 ifstream myfile;
38 myfile.open(fcalname.Data());
39 if ( !myfile ){
40 return(-107);
41 };
42 myfile.close();
43 //
44 TFile *File = new TFile(fcalname.Data());
45 if ( !File ) return(-108);
46 TTree *tr = (TTree*)File->Get("CalibCalPed");
47 if ( !tr ) return(-109);
48 //
49 TBranch *calo = tr->GetBranch("CalibCalPed");
50 //
51 pamela::CalibCalPedEvent *ce = 0;
52 tr->SetBranchAddress("CalibCalPed", &ce);
53 //
54 Long64_t ncalibs = calo->GetEntries();
55 //
56 if ( !ncalibs ) return(-110);
57 //
58 calo->GetEntry(calibno);
59 //
60 if (ce->cstwerr[s] != 0 && ce->cperror[s] == 0 ) {
61 fcstwerr[s] = ce->cstwerr[s];
62 fcperror[s] = ce->cperror[s];
63 for ( Int_t d=0 ; d<11 ;d++ ){
64 Int_t pre = -1;
65 for ( Int_t j=0; j<96 ;j++){
66 if ( j%16 == 0 ) pre++;
67 fcalped[s][d][j] = ce->calped[e][d][j];
68 fcalgood[s][d][j] = ce->calgood[e][d][j];
69 fcalthr[s][d][pre] = ce->calthr[e][d][pre];
70 fcalrms[s][d][j] = ce->calrms[e][d][j];
71 fcalbase[s][d][pre] = ce->calbase[e][d][pre];
72 fcalvar[s][d][pre] = ce->calvar[e][d][pre];
73 };
74 };
75 } else {
76 printf(" CALORIMETER - ERROR: problems finding a good calibration in this file! \n\n ");
77 File->Close();
78 return(-111);
79 };
80 File->Close();
81 return(0);
82 }
83
84
85 void Digitizer::DigitizeCALOCALIB() {
86 //
87 // Header of the four sections
88 //
89 fSecCalo[0] = 0xAA00; // XE
90 fSecCalo[1] = 0xB100; // XO
91 fSecCalo[2] = 0xB600; // YE
92 fSecCalo[3] = 0xAD00; // YO
93 //
94 // length of the data is 0x1215
95 //
96 fSecCALOLength[0] = 0x1215; // XE
97 fSecCALOLength[1] = 0x1215; // XO
98 fSecCALOLength[2] = 0x1215; // YE
99 fSecCALOLength[3] = 0x1215; // YO
100 //
101 Int_t chksum = 0;
102 UInt_t tstrip = 0;
103 UInt_t fSecPointer = 0;
104 //
105 for (Int_t sec=0; sec < 4; sec++){
106 //
107 // sec = 0 -> XE 1 -> XO 2-> YE 3 -> YO
108 //
109 fCALOlength = 0;
110 memset(fDataCALO,0,sizeof(UShort_t)*fCALObuffer);
111 fSecPointer = fCALOlength;
112 //
113 // First of all we have section header and packet length
114 //
115 fDataCALO[fCALOlength] = fSecCalo[sec];
116 fCALOlength++;
117 fDataCALO[fCALOlength] = fSecCALOLength[sec];
118 fCALOlength++;
119 //
120 // Section XO is read in the opposite direction respect to the others
121 //
122 chksum = 0;
123 //
124 for (Int_t plane=0; plane < 11; plane++){
125 //
126 if ( sec == 1 ) tstrip = fCALOlength + 96*2;
127 //
128 for (Int_t strip=0; strip < 96; strip++){
129 //
130 chksum += (Int_t)fcalped[sec][plane][strip];
131 //
132 // save value
133 //
134 if ( sec == 1 ){
135 tstrip -= 2;
136 fDataCALO[tstrip] = (Int_t)fcalped[sec][plane][strip];
137 fDataCALO[tstrip+1] = (Int_t)fcalgood[sec][plane][strip];
138 } else {
139 fDataCALO[fCALOlength] = (Int_t)fcalped[sec][plane][strip];
140 fDataCALO[fCALOlength+1] = (Int_t)fcalgood[sec][plane][strip];
141 };
142 fCALOlength +=2;
143 };
144 //
145 };
146 //
147 fDataCALO[fCALOlength] = (UShort_t)chksum;
148 fCALOlength++;
149 fDataCALO[fCALOlength] = 0;
150 fCALOlength++;
151 fDataCALO[fCALOlength] = (UShort_t)((Int_t)(chksum >> 16));
152 fCALOlength++;
153 //
154 // Section XO is read in the opposite direction respect to the others
155 //
156 chksum = 0;
157 //
158 for (Int_t plane=0; plane < 11; plane++){
159 //
160 if ( sec == 1 ) tstrip = fCALOlength+6*2;
161 //
162 for (Int_t strip=0; strip < 6; strip++){
163 //
164 chksum += (Int_t)fcalthr[sec][plane][strip];
165 //
166 // save value
167 //
168 if ( sec == 1 ){
169 tstrip -= 2;
170 fDataCALO[tstrip] = 0;
171 fDataCALO[tstrip+1] = (Int_t)fcalthr[sec][plane][strip];
172 } else {
173 fDataCALO[fCALOlength] = 0;
174 fDataCALO[fCALOlength+1] = (Int_t)fcalthr[sec][plane][strip];
175 };
176 fCALOlength +=2;
177 };
178 //
179 };
180 //
181 fDataCALO[fCALOlength] = 0;
182 fCALOlength++;
183 fDataCALO[fCALOlength] = (UShort_t)chksum;
184 fCALOlength++;
185 fDataCALO[fCALOlength] = 0;
186 fCALOlength++;
187 fDataCALO[fCALOlength] = (UShort_t)((Int_t)(chksum >> 16));
188 fCALOlength++;
189 //
190 // Section XO is read in the opposite direction respect to the others
191 //
192 for (Int_t plane=0; plane < 11; plane++){
193 //
194 if ( sec == 1 ) tstrip = fCALOlength+96*2;
195 //
196 for (Int_t strip=0; strip < 96; strip++){
197 //
198 // save value
199 //
200 if ( sec == 1 ){
201 tstrip -= 2;
202 fDataCALO[tstrip] = 0;
203 fDataCALO[tstrip+1] = (Int_t)fcalrms[sec][plane][strip];
204 } else {
205 fDataCALO[fCALOlength] = 0;
206 fDataCALO[fCALOlength+1] = (Int_t)fcalrms[sec][plane][strip];
207 };
208 fCALOlength += 2;
209 };
210 //
211 };
212 //
213 // Section XO is read in the opposite direction respect to the others
214 //
215 for (Int_t plane=0; plane < 11; plane++){
216 //
217 if ( sec == 1 ) tstrip = fCALOlength+6*4;
218 //
219 for (Int_t strip=0; strip < 6; strip++){
220 //
221 // save value
222 //
223 if ( sec == 1 ){
224 tstrip -= 4;
225 fDataCALO[tstrip] = 0;
226 fDataCALO[tstrip+1] = (Int_t)fcalbase[sec][plane][strip];
227 fDataCALO[tstrip+2] = 0;
228 fDataCALO[tstrip+3] = (Int_t)fcalvar[sec][plane][strip];
229 } else {
230 fDataCALO[fCALOlength] = 0;
231 fDataCALO[fCALOlength+1] = (Int_t)fcalbase[sec][plane][strip];
232 fDataCALO[fCALOlength+2] = 0;
233 fDataCALO[fCALOlength+3] = (Int_t)fcalvar[sec][plane][strip];
234 };
235 fCALOlength +=4;
236 };
237 //
238 };
239 //
240 //
241 // here we calculate and save the CRC
242 //
243 fDataCALO[fCALOlength] = 0;
244 fCALOlength++;
245 Short_t CRC = 0;
246 for (UInt_t i=0; i<(fCALOlength-fSecPointer); i++){
247 CRC=crc(CRC,fDataCALO[i+fSecPointer]);
248 };
249 fDataCALO[fCALOlength] = (UShort_t)CRC;
250 fCALOlength++;
251 //
252 UInt_t length=fCALOlength*2;
253 DigitizePSCU(length,0x18,fDataPSCU);
254 //
255 // Add padding to 64 bits
256 //
257 AddPadding();
258 //
259 fOutputfile.write(reinterpret_cast<char*>(fDataPSCU),sizeof(UShort_t)*fPSCUbuffer);
260 UShort_t temp[1000000];
261 memset(temp,0,sizeof(UShort_t)*1000000);
262 swab(fDataCALO,temp,sizeof(UShort_t)*fCALOlength); // WE MUST SWAP THE BYTES!!!
263 fOutputfile.write(reinterpret_cast<char*>(temp),sizeof(UShort_t)*fCALOlength);
264 //
265 // padding to 64 bytes
266 //
267 if ( fPadding ){
268 fOutputfile.write(reinterpret_cast<char*>(fDataPadding),sizeof(UChar_t)*fPadding);
269 };
270 //
271 //
272 };
273 //
274 };
275
276 void Digitizer::CaloLoadCalib() {
277 //
278 fGivenCaloCalib = 0; // ####@@@@ should be given as input par @@@@####
279 //
280 // first of all load the MIP to ADC conversion values
281 //
282 stringstream calfile;
283 Int_t error = 0;
284 GL_PARAM *glparam = new GL_PARAM();
285 //
286 // determine where I can find calorimeter ADC to MIP conversion file
287 //
288 error = 0;
289 error = glparam->Query_GL_PARAM(3,101,fDbc);
290 //
291 calfile.str("");
292 calfile << glparam->PATH.Data() << "/";
293 calfile << glparam->NAME.Data();
294 //
295 printf("\n Using Calorimeter ADC to MIP conversion file: \n %s \n",calfile.str().c_str());
296 FILE *f;
297 f = fopen(calfile.str().c_str(),"rb");
298 //
299 memset(fCalomip,0,4224*sizeof(fCalomip[0][0][0]));
300 //
301 for (Int_t m = 0; m < 2 ; m++ ){
302 for (Int_t k = 0; k < 22; k++ ){
303 for (Int_t l = 0; l < 96; l++ ){
304 fread(&fCalomip[m][k][l],sizeof(fCalomip[m][k][l]),1,f);
305 };
306 };
307 };
308 fclose(f);
309 //
310 // determine which calibration has to be used and load it for each section
311 //
312 GL_CALO_CALIB *glcalo = new GL_CALO_CALIB();
313 GL_ROOT *glroot = new GL_ROOT();
314 TString fcalname;
315 UInt_t idcalib;
316 UInt_t calibno;
317 UInt_t utime = 0;
318 //
319 for (UInt_t s=0; s<4; s++){
320 //
321 // clear calo calib variables for section s
322 //
323 ClearCaloCalib(s);
324 //
325 if ( fGivenCaloCalib ){
326 //
327 // a time has been given as input on the command line so retrieve the calibration that preceed that time
328 //
329 glcalo->Query_GL_CALO_CALIB(fGivenCaloCalib,utime,s,fDbc);
330 //
331 calibno = glcalo->EV_ROOT;
332 idcalib = glcalo->ID_ROOT_L0;
333 //
334 // determine path and name and entry of the calibration file
335 //
336 printf("\n");
337 printf(" ** SECTION %i **\n",s);
338 //
339 glroot->Query_GL_ROOT(idcalib,fDbc);
340 //
341 stringstream name;
342 name.str("");
343 name << glroot->PATH.Data() << "/";
344 name << glroot->NAME.Data();
345 //
346 fcalname = (TString)name.str().c_str();
347 //
348 printf("\n Section %i : using file %s calibration at entry %i: \n",s,fcalname.Data(),calibno);
349 //
350 } else {
351 error = 0;
352 error = glparam->Query_GL_PARAM(1,104,fDbc);
353 //
354 calfile.str("");
355 calfile << glparam->PATH.Data() << "/";
356 calfile << glparam->NAME.Data();
357 //
358 printf("\n Section %i : using default calorimeter calibration: \n %s \n",s,calfile.str().c_str());
359 //
360 fcalname = (TString)calfile.str().c_str();
361 calibno = s;
362 //
363 };
364 //
365 // load calibration variables in memory
366 //
367 CaloLoadCalib(s,fcalname,calibno);
368 //
369 };
370 //
371 // 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
372 //
373 delete glparam;
374 delete glcalo;
375 delete glroot;
376 };
377
378 void Digitizer::DigitizeCALO() {
379 //
380 fModCalo = 0; // 0 is RAW, 1 is COMPRESS, 2 is FULL ####@@@@ should be given as input par @@@@####
381 //
382 //
383 //
384 fCALOlength = 0; // reset total dimension of calo data
385 //
386 // gpamela variables to be used
387 //
388 // fhBookTree->SetBranchStatus("Nthcali",1);//modified by E.Vannuccini 03/08
389 // fhBookTree->SetBranchStatus("Icaplane",1);
390 // fhBookTree->SetBranchStatus("Icastrip",1);
391 // fhBookTree->SetBranchStatus("Icamod",1);
392 // fhBookTree->SetBranchStatus("Enestrip",1);
393 //
394 // call different routines depending on the acq mode you want to simulate
395 //
396 switch ( fModCalo ){
397 case 0:
398 this->DigitizeCALORAW();
399 break;
400 case 1:
401 this->DigitizeCALOCOMPRESS();
402 break;
403 case 2:
404 this->DigitizeCALOFULL();
405 break;
406 };
407 //
408 };
409
410 Float_t Digitizer::GetCALOen(Int_t sec, Int_t plane, Int_t strip){
411 //
412 // determine plane and strip
413 //
414 Int_t mplane = 0;
415 //
416 // wrong!
417 //
418 // if ( sec == 0 || sec == 3 ) mplane = (plane * 4) + sec + 1;
419 // if ( sec == 1 ) mplane = (plane * 4) + 2 + 1;
420 // if ( sec == 2 ) mplane = (plane * 4) + 1 + 1;
421 //
422 if ( sec == 0 ) mplane = plane * 4 + 1; // it must be 0, 4, 8, ... (+1) from plane = 0, 11
423 if ( sec == 1 ) mplane = plane * 4 + 2 + 1; // it must be 2, 6, 10, ... (+1) from plane = 0, 11
424 if ( sec == 2 ) mplane = plane * 4 + 3 + 1; // it must be 3, 7, 11, ... (+1) from plane = 0, 11
425 if ( sec == 3 ) mplane = plane * 4 + 1 + 1; // it must be 1, 5, 9, ... (+1) from plane = 0, 11
426 //
427 Int_t mstrip = strip + 1;
428 //
429 // search energy release in gpamela output
430 //
431 for (Int_t i=0; i<Nthcali;i++){
432 if ( Icaplane[i] == mplane && Icastrip[i] == mstrip ){
433 return (Enestrip[i]);
434 };
435 };
436 //
437 // if not found it means no energy release so return 0.
438 //
439 return(0.);
440 };
441
442 void Digitizer::DigitizeCALORAW() {
443 //
444 // some variables
445 //
446 Float_t ens = 0.;
447 UInt_t adcsig = 0;
448 UInt_t adcbase = 0;
449 UInt_t adc = 0;
450 Int_t pre = 0;
451 UInt_t l = 0;
452 UInt_t lpl = 0;
453 UInt_t tstrip = 0;
454 UInt_t fSecPointer = 0;
455 Double_t pedenoise;
456 Float_t rms = 0.;
457 Float_t pedestal = 0.;
458 //
459 // clean the data structure
460 //
461 memset(fDataCALO,0,sizeof(UShort_t)*fCALObuffer);
462 //
463 // Header of the four sections
464 //
465 fSecCalo[0] = 0xEA08; // XE
466 fSecCalo[1] = 0xF108; // XO
467 fSecCalo[2] = 0xF608; // YE
468 fSecCalo[3] = 0xED08; // YO
469 //
470 // length of the data is 0x0428 in RAW mode
471 //
472 fSecCALOLength[0] = 0x0428; // XE
473 fSecCALOLength[1] = 0x0428; // XO
474 fSecCALOLength[2] = 0x0428; // YE
475 fSecCALOLength[3] = 0x0428; // YO
476 //
477 // let's start
478 //
479 fCALOlength = 0;
480 //
481 for (Int_t sec=0; sec < 4; sec++){
482 //
483 // sec = 0 -> XE 1 -> XO 2-> YE 3 -> YO
484 //
485 l = 0; // XE and XO are Y planes
486 if ( sec < 2 ) l = 1; // while YE and YO are X planes
487 //
488 fSecPointer = fCALOlength;
489 //
490 // First of all we have section header and packet length
491 //
492 fDataCALO[fCALOlength] = fSecCalo[sec];
493 fCALOlength++;
494 fDataCALO[fCALOlength] = fSecCALOLength[sec];
495 fCALOlength++;
496 //
497 // selftrigger coincidences - in the future we should add here some code to simulate timing response of pre-amplifiers
498 //
499 for (Int_t autoplane=0; autoplane < 7; autoplane++){
500 fDataCALO[fCALOlength] = 0x0000;
501 fCALOlength++;
502 };
503 //
504 //
505 // here comes data
506 //
507 //
508 // Section XO is read in the opposite direction respect to the others
509 //
510 if ( sec == 1 ){
511 tstrip = 96*11 + fCALOlength;
512 } else {
513 tstrip = 0;
514 };
515 //
516 pre = -1;
517 //
518 for (Int_t strip=0; strip < 96; strip++){
519 //
520 // which is the pre for this strip?
521 //
522 if (strip%16 == 0) {
523 pre++;
524 };
525 //
526 if ( sec == 1 ) tstrip -= 11;
527 //
528 for (Int_t plane=0; plane < 11; plane++){
529 //
530 // here is wrong!!!!
531 //
532 //
533 // if ( plane%2 == 0 && sec%2 != 0){
534 // lpl = plane*2;
535 // } else {
536 // lpl = (plane*2) + 1;
537 // };
538 //
539 if ( sec == 0 || sec == 3 ) lpl = plane * 2;
540 if ( sec == 1 || sec == 2 ) lpl = (plane * 2) + 1;
541 //
542 // get the energy in GeV from the simulation for that strip
543 //
544 ens = this->GetCALOen(sec,plane,strip);
545 //
546 // convert it into ADC channels
547 //
548 adcsig = int(ens*fCalomip[l][lpl][strip]/fCALOGeV2MIPratio);
549 //
550 // sum baselines
551 //
552 adcbase = (UInt_t)fcalbase[sec][plane][pre];
553 //
554 // add noise and pedestals
555 //
556 pedestal = fcalped[sec][plane][strip];
557 rms = fcalrms[sec][plane][strip]/4.;
558 //
559 // Add random gaussian noise of RMS rms and Centered in the pedestal
560 //
561 pedenoise = gRandom->Gaus((Double_t)pedestal,(Double_t)rms);
562 //
563 // Sum all contribution
564 //
565 adc = adcsig + adcbase + (Int_t)round(pedenoise);
566 //
567 // Signal saturation
568 //
569 if ( adc > 0x7FFF ) adc = 0x7FFF;
570 //
571 // save value
572 //
573 if ( sec == 1 ){
574 fDataCALO[tstrip] = adc;
575 tstrip++;
576 } else {
577 fDataCALO[fCALOlength] = adc;
578 };
579 fCALOlength++;
580 //
581 };
582 //
583 if ( sec == 1 ) tstrip -= 11;
584 //
585 };
586 //
587 // here we calculate and save the CRC
588 //
589 Short_t CRC = 0;
590 for (UInt_t i=0; i<(fCALOlength-fSecPointer); i++){
591 CRC=crc(CRC,fDataCALO[i+fSecPointer]);
592 };
593 fDataCALO[fCALOlength] = (UShort_t)CRC;
594 fCALOlength++;
595 //
596 };
597 //
598 // for (Int_t i=0; i<fCALOlength; i++){
599 // printf(" WORD %i DIGIT %0x \n",i,fDataCALO[i]);
600 // };
601 //
602 };
603
604 void Digitizer::DigitizeCALOCOMPRESS() {
605 //
606 printf(" COMPRESS MODE STILL NOT IMPLEMENTED! \n");
607 //
608 this->DigitizeCALORAW();
609 return;
610 //
611 //
612 //
613 fSecCalo[0] = 0xEA00;
614 fSecCalo[1] = 0xF100;
615 fSecCalo[2] = 0xF600;
616 fSecCalo[3] = 0xED00;
617 //
618 // length of the data in DSP mode must be calculated on fly during digitization
619 //
620 memset(fSecCALOLength,0x0,4*sizeof(UShort_t));
621 //
622 // here comes raw data
623 //
624 Int_t en = 0;
625 //
626 for (Int_t sec=0; sec < 4; sec++){
627 fDataCALO[en] = fSecCalo[sec];
628 en++;
629 fDataCALO[en] = fSecCALOLength[sec];
630 en++;
631 for (Int_t plane=0; plane < 11; plane++){
632 for (Int_t strip=0; strip < 11; strip++){
633 fDataCALO[en] = 0x0;
634 en++;
635 };
636 };
637 };
638 //
639 };
640
641 void Digitizer::DigitizeCALOFULL() {
642 //
643 printf(" FULL MODE STILL NOT IMPLEMENTED! \n");
644 //
645 this->DigitizeCALORAW();
646 return;
647 //
648 fSecCalo[0] = 0xEA00;
649 fSecCalo[1] = 0xF100;
650 fSecCalo[2] = 0xF600;
651 fSecCalo[3] = 0xED00;
652 //
653 // length of the data in DSP mode must be calculated on fly during digitization
654 //
655 memset(fSecCALOLength,0x0,4*sizeof(UShort_t));
656 //
657 // here comes raw data
658 //
659 Int_t en = 0;
660 //
661 for (Int_t sec=0; sec < 4; sec++){
662 fDataCALO[en] = fSecCalo[sec];
663 en++;
664 fDataCALO[en] = fSecCALOLength[sec];
665 en++;
666 for (Int_t plane=0; plane < 11; plane++){
667 for (Int_t strip=0; strip < 11; strip++){
668 fDataCALO[en] = 0x0;
669 en++;
670 };
671 };
672 };
673 //
674 };

  ViewVC Help
Powered by ViewVC 1.1.23