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

Annotation of /PamelaDigitizer/DigitizeCalo.cxx

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.2 - (hide 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 pamelats 1.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