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

Annotation of /PamelaDigitizer/DigitizeCalo.cxx

Parent Directory Parent Directory | Revision Log Revision Log


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