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

Contents of /PamelaDigitizer/DigitizeCalo.cxx

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1 - (show annotations) (download)
Wed May 21 09:50:39 2008 UTC (16 years, 6 months ago) by pamelats
Branch: MAIN
CVS Tags: v3r00, v3r01, v3r02
*** empty log message ***

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

  ViewVC Help
Powered by ViewVC 1.1.23