/[PAMELA software]/PamVMC/cal/src/PamVMCCaloDig.cxx
ViewVC logotype

Contents of /PamVMC/cal/src/PamVMCCaloDig.cxx

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1 - (show annotations) (download)
Thu Feb 19 16:00:14 2009 UTC (15 years, 9 months ago) by nikolas
Branch: MAIN
Cleaning files before release

1 #include "PamVMCCaloDig.h"
2 #include <TTree.h>
3 #include "CalibCalPedEvent.h"
4
5 using namespace pamela;
6
7 ClassImp(PamVMCCaloDig)
8
9 extern "C"{
10 short crc(short, short);
11 };
12
13
14 void PamVMCCaloDig::LoadCalib(){
15
16 cout<<"Loading Tracker Calibrations..."<<endl;
17
18
19
20 fhc =(TClonesArray*)fhitscolmap.GetValue("CAST");
21 if (!fhc) cout <<"!!!WARNING Hit Collection of CAL not found by PamVMCCaloDig!!!!"<<endl;
22
23 Int_t GivenCaloCalib = 0; //@@@ should be given as input par @@@
24
25 Int_t time=3;
26 Int_t type=101;
27
28 fdberr = fsql->Query_GL_PARAM(time,type);
29
30 if(fdberr<0){
31 cout<<"No such record in DB for CAL: time="<<time<<" type="<<type<<endl;
32 ThrowCalFileWarning("CAL");
33
34 } else {
35
36 fquery.str("");
37 fquery << fsql->GetPAR()->PATH.Data() << "/";
38 fquery << fsql->GetPAR()->NAME.Data();
39
40 ThrowCalFileUsage("CAL",fquery.str().c_str());
41
42
43 FILE *f;
44 f = fopen(fquery.str().c_str(),"rb");
45
46 if(f==NULL) ThrowCalFileWarning("CAL");
47 else{
48 memset(fCalomip,0,4224*sizeof(fCalomip[0][0][0]));
49
50 for (Int_t m = 0; m < 2 ; m++ ){
51 for (Int_t k = 0; k < 22; k++ ){
52 for (Int_t l = 0; l < 96; l++ ){
53 fread(&fCalomip[m][k][l],sizeof(fCalomip[m][k][l]),1,f);
54 };
55 };
56 };
57 fclose(f);
58
59
60 // determine which calibration has to be used
61 // and load it for each section
62
63 UInt_t idcalib, calibno;
64 UInt_t utime = 0;
65 TString calname;
66 for (UInt_t s=0; s<4; s++){
67
68 // clear calo calib variables for section s
69
70 ClearCaloCalib(s);
71
72 if ( GivenCaloCalib ){
73
74 // a time has been given as input on the command line
75 // so retrieve the calibration that preceed that time
76
77 fsql->Query_GL_CALO_CALIB(GivenCaloCalib,utime,s);
78
79 calibno = fsql->GetCaloCalib()->EV_ROOT;
80 idcalib = fsql->GetCaloCalib()->ID_ROOT_L0;
81
82 // determine path and name and entry of the calibration file
83
84 printf("\n");
85 printf(" ** SECTION %i **\n",s);
86
87 fsql->Query_GL_ROOT(idcalib);
88
89 fquery.str("");
90 fquery << fsql->GetROOT()->PATH.Data() << "/";
91 fquery << fsql->GetROOT()->NAME.Data();
92
93 calname = (TString)fquery.str().c_str();
94
95 printf("\n Section %i : using file %s calibration at entry %i: \n",s,calname.Data(),calibno);
96
97 } else {
98
99 fdberr = 0;
100 fdberr = fsql->Query_GL_PARAM(1,104);
101
102 fquery.str("");
103 fquery << fsql->GetPAR()->PATH.Data() << "/";
104 fquery << fsql->GetPAR()->NAME.Data();
105
106 printf("\n Section %i : using default calorimeter calibration: \n %s \n",s,fquery.str().c_str());
107
108 calname = (TString)fquery.str().c_str();
109 calibno = s;
110
111 };
112
113 // load calibration variables in memory
114
115 CaloLoadCalib(s,calname,calibno);
116
117 };
118 //
119 // at this point we have in memory the calorimeter calibration
120 // and we can save it to disk in the correct format and use it to digitize the data
121 //
122
123 DigitizeCALOCALIB();
124 }
125 }
126 }
127
128
129 void PamVMCCaloDig::ClearCaloCalib(Int_t s){
130
131 fcstwerr[s] = 0;
132 fcperror[s] = 0.;
133 for ( Int_t d=0 ; d<11 ;d++ ){
134 Int_t pre = -1;
135 for ( Int_t j=0; j<96 ;j++){
136 if ( j%16 == 0 ) pre++;
137 fcalped[s][d][j] = 0.;
138 fcstwerr[s] = 0.;
139 fcperror[s] = 0.;
140 fcalgood[s][d][j] = 0.;
141 fcalthr[s][d][pre] = 0.;
142 fcalrms[s][d][j] = 0.;
143 fcalbase[s][d][pre] = 0.;
144 fcalvar[s][d][pre] = 0.;
145 };
146 };
147 return;
148 }
149
150 Int_t PamVMCCaloDig::CaloLoadCalib(Int_t s,TString calname, UInt_t calibno){
151
152 UInt_t e = 0;
153 switch(s){
154 case 0: e = 0;
155 break;
156 case 1: e = 2;
157 break;
158 case 2: e = 3;
159 break;
160 case 3: e = 1;
161 break;
162 default:
163 break;
164 }
165 fcfile.open(calname.Data());
166 ThrowCalFileUsage("CAL",calname.Data());
167 if ( !fcfile ){
168 ThrowCalFileWarning("CAL");
169 return(-107);
170 };
171 fcfile.close();
172
173 fcrfile = new TFile(calname.Data());
174 if ( !fcrfile ){
175 ThrowCalFileWarning("CAL");
176 return(-108);
177 }
178 TTree *tr = (TTree*)fcrfile->Get("CalibCalPed");
179 if ( !tr ){
180 cout<<"!!!WARNING Tree CalibCalPed in file "<<calname.Data()
181 <<" was NOT found!!!"<<endl;
182 return(-109);
183 }
184 TBranch *calo = tr->GetBranch("CalibCalPed");
185 CalibCalPedEvent *ce = 0;
186 tr->SetBranchAddress("CalibCalPed", &ce);
187 Long64_t ncalibs = calo->GetEntries();
188 if ( !ncalibs ){
189 cout<<"!!!WARNING No entries in calibration file!!!"<<endl;
190 return(-110);
191 }
192
193 calo->GetEntry(calibno);
194
195 if (ce->cstwerr[s] != 0 && ce->cperror[s] == 0 ) {
196 fcstwerr[s] = ce->cstwerr[s];
197 fcperror[s] = ce->cperror[s];
198 for ( Int_t d=0 ; d<11 ;d++ ){
199 Int_t pre = -1;
200 for ( Int_t j=0; j<96 ;j++){
201 if ( j%16 == 0 ) pre++;
202 fcalped[s][d][j] = ce->calped[e][d][j];
203 fcalgood[s][d][j] = ce->calgood[e][d][j];
204 fcalthr[s][d][pre] = ce->calthr[e][d][pre];
205 fcalrms[s][d][j] = ce->calrms[e][d][j];
206 fcalbase[s][d][pre] = ce->calbase[e][d][pre];
207 fcalvar[s][d][pre] = ce->calvar[e][d][pre];
208 };
209 };
210 } else {
211 printf(" CALORIMETER - ERROR: problems finding a good calibration in this file! \n\n ");
212 fcrfile->Close();
213 return(-111);
214 };
215 fcrfile->Close();
216 delete ce;
217 return(0);
218 }
219
220 void PamVMCCaloDig::DigitizeCALOCALIB() {
221
222
223 // Header of the four sections
224
225 fSecCalo[0] = 0xAA00; // XE
226 fSecCalo[1] = 0xB100; // XO
227 fSecCalo[2] = 0xB600; // YE
228 fSecCalo[3] = 0xAD00; // YO
229
230 // length of the data is 0x1215
231
232 fSecCaloLength[0] = 0x1215; // XE
233 fSecCaloLength[1] = 0x1215; // XO
234 fSecCaloLength[2] = 0x1215; // YE
235 fSecCaloLength[3] = 0x1215; // YO
236
237 Short_t CRC;
238
239 for(Int_t sec=0; sec < 4; sec++){
240 //
241 // sec = 0 -> XE 1 -> XO 2-> YE 3 -> YO
242 //
243 fData.clear();
244
245 fData.push_back(fSecCalo[sec]);
246 fData.push_back(fSecCaloLength[sec]);
247
248 if ( sec==1 ){
249 FillCalPedBack(&fData,sec);
250 FillCalThrBack(&fData,sec);
251 FillCalRmsBack(&fData,sec);
252 FillCalBaseVarBack(&fData,sec);
253 } else {
254 FillCalPedNorm(&fData,sec);
255 FillCalThrNorm(&fData,sec);
256 FillCalRmsNorm(&fData,sec);
257 FillCalBaseVarNorm(&fData,sec);
258 }
259 //CRC
260
261 fData.push_back(0x0000);
262 CRC = 0;
263
264
265 USBuffer::const_iterator p = fData.begin();
266 while( p!= fData.end() ){
267 CRC=crc(CRC,*p);
268 p++;
269 }
270 fData.push_back(CRC);
271
272 //Writing calibrations
273 DigitizePSCU(fData.size()*2,0x18);
274 fraw->WritePSCU(&fDataPSCU);
275 fraw->CopyUShortToBuff(&fData);
276
277 //No padding for now...
278 }
279 }
280
281 void PamVMCCaloDig::FillCalPedNorm(USBuffer *buff, Int_t sec){
282
283 Int_t chksum = 0;
284 for (Int_t plane=0; plane < 11; plane++){
285 for (Int_t strip=0; strip < 96; strip++){
286 chksum += (Int_t)fcalped[sec][plane][strip];
287 buff->push_back((Int_t)fcalped[sec][plane][strip]);
288 buff->push_back((Int_t)fcalgood[sec][plane][strip]);
289 }
290
291 }//3-2114
292
293 buff->push_back(((UShort_t)chksum));
294 buff->push_back(0x0000);
295 buff->push_back(((UShort_t)((Int_t)(chksum >> 16))));
296 }
297
298 void PamVMCCaloDig::FillCalPedBack(USBuffer *buff, Int_t sec){
299
300 Int_t chksum = 0;
301 for (Int_t plane=0; plane < 11; plane++){
302 for (Int_t strip=95; strip >= 0; strip--){
303 chksum += (Int_t)fcalped[sec][plane][strip];
304 buff->push_back((Int_t)fcalped[sec][plane][strip]);
305 buff->push_back((Int_t)fcalgood[sec][plane][strip]);
306 }
307 }//3-2114
308
309 buff->push_back(((UShort_t)chksum));
310 buff->push_back(0x0000);
311 buff->push_back(((UShort_t)((Int_t)(chksum >> 16))));
312 }
313
314
315 void PamVMCCaloDig::FillCalThrNorm(USBuffer *buff, Int_t sec){
316
317 Int_t chksum = 0;
318 for (Int_t plane=0; plane < 11; plane++){
319 for (Int_t strip=0; strip < 6; strip++){
320 chksum += (Int_t)fcalthr[sec][plane][strip];
321 buff->push_back(0x0000);
322 buff->push_back((Int_t)fcalthr[sec][plane][strip]);
323 }
324
325 }//2118-2249
326
327 buff->push_back(0x0000);
328 buff->push_back(((UShort_t)chksum));
329 buff->push_back(0x0000);
330 buff->push_back(((UShort_t)((Int_t)(chksum >> 16))));
331 }
332
333 void PamVMCCaloDig::FillCalThrBack(USBuffer *buff, Int_t sec){
334
335 Int_t chksum = 0;
336 for (Int_t plane=0; plane < 11; plane++){
337 for (Int_t strip=5; strip >= 0; strip--){
338 chksum += (Int_t)fcalthr[sec][plane][strip];
339 buff->push_back(0x0000);
340 buff->push_back((Int_t)fcalthr[sec][plane][strip]);
341 }
342 }//2118-2249
343
344 buff->push_back(0x0000);
345 buff->push_back(((UShort_t)chksum));
346 buff->push_back(0x0000);
347 buff->push_back(((UShort_t)((Int_t)(chksum >> 16))));
348 }
349
350 void PamVMCCaloDig::FillCalRmsNorm(USBuffer *buff, Int_t sec){
351
352 for (Int_t plane=0; plane < 11; plane++){
353 for (Int_t strip=0; strip < 96; strip++){
354 buff->push_back(0x0000);
355 buff->push_back((Int_t)fcalrms[sec][plane][strip]);
356 }
357 }//2254-4365
358 }
359
360 void PamVMCCaloDig::FillCalRmsBack(USBuffer *buff, Int_t sec){
361
362 for (Int_t plane=0; plane < 11; plane++){
363 for (Int_t strip=95; strip >= 0; strip--){
364 buff->push_back(0x0000);
365 buff->push_back((Int_t)fcalrms[sec][plane][strip]);
366 }
367 }//2254-4365
368 }
369
370 void PamVMCCaloDig::FillCalBaseVarNorm(USBuffer *buff, Int_t sec){
371
372 for (Int_t plane=0; plane < 11; plane++){
373 for (Int_t strip=0; strip < 6; strip++){
374 buff->push_back(0x0000);
375 buff->push_back(((Int_t)fcalbase[sec][plane][strip]));
376 buff->push_back(0x0000);
377 buff->push_back(((Int_t)fcalvar[sec][plane][strip]));
378 }
379 }//4366-4629
380 }
381
382 void PamVMCCaloDig::FillCalBaseVarBack(USBuffer *buff, Int_t sec){
383
384 for (Int_t plane=0; plane < 11; plane++){
385 for (Int_t strip=5; strip >= 0; strip--){
386 buff->push_back(0x0000);
387 buff->push_back((Int_t)fcalbase[sec][plane][strip]);
388 buff->push_back(0x0000);
389 buff->push_back((Int_t)fcalvar[sec][plane][strip]);
390 }
391 }//4366-4629
392 }
393
394 void PamVMCCaloDig::Digitize(){
395
396 cout<<"Digitizing CALO..."<<endl;
397
398 Int_t ModCalo = 0; // 0 is RAW, 1 is COMPRESS, 2 is FULL
399 //####@@@@ should be given as input par @@@@####
400
401
402 // call different routines depending on the acq mode you want to simulate
403
404 switch ( ModCalo ){
405 case 0:
406 this->DigitizeCaloRaw();
407 break;
408 case 1:
409 this->DigitizeCaloCompress();
410 break;
411 case 2:
412 this->DigitizeCaloFull();
413 break;
414 };
415 }
416
417 void PamVMCCaloDig::DigitizeCaloRaw(){
418
419
420 // some variables
421 //
422 Float_t ens = 0.;
423 UInt_t adcsig = 0;
424 UInt_t adcbase = 0;
425 UInt_t adc = 0;
426 Int_t pre = 0;
427 UInt_t l = 0;
428 UInt_t lpl = 0;
429 UInt_t tstrip = 0;
430 UInt_t fSecPointer = 0;
431 Int_t fCALOlength;
432 Double_t pedenoise;
433 Float_t rms = 0.;
434 Float_t pedestal = 0.;
435 //
436 // clean the data structure
437 //
438
439 UShort_t DataCALO[4264];
440 fData.clear();
441
442 memset(DataCALO,0,sizeof(UShort_t)*4264);
443
444 static const Float_t CALOGeV2MIPratio = 0.0001059994;
445
446 // Header of the four sections
447
448 fSecCalo[0] = 0xEA08; // XE
449 fSecCalo[1] = 0xF108; // XO
450 fSecCalo[2] = 0xF608; // YE
451 fSecCalo[3] = 0xED08; // YO
452
453 // length of the data is 0x0428 in RAW mode
454
455 fSecCaloLength[0] = 0x0428; // XE
456 fSecCaloLength[1] = 0x0428; // XO
457 fSecCaloLength[2] = 0x0428; // YE
458 fSecCaloLength[3] = 0x0428; // YO
459 // let's start
460 //
461 fCALOlength = 0;
462 //
463 for (Int_t sec=0; sec < 4; sec++){
464 //
465 // sec = 0 -> XE 1 -> XO 2-> YE 3 -> YO
466 //
467 l = 0; // XE and XO are Y planes
468 if ( sec < 2 ) l = 1; // while YE and YO are X planes
469 //
470 fSecPointer = fCALOlength;
471 //
472 // First of all we have section header and packet length
473 //
474 DataCALO[fCALOlength] = fSecCalo[sec]; //1
475 fCALOlength++;
476 DataCALO[fCALOlength] = fSecCaloLength[sec]; //2
477 fCALOlength++;
478 //
479 // selftrigger coincidences - in the future we should add here some code
480 // to simulate timing response of pre-amplifiers
481 //
482 for (Int_t autoplane=0; autoplane < 7; autoplane++){
483 DataCALO[fCALOlength] = 0x0000;
484 fCALOlength++;
485 }; //3-8
486 //
487 //
488 // here comes data
489 //
490 //
491 // Section XO is read in the opposite direction respect to the others
492 //
493 if ( sec == 1 ){
494 tstrip = 96*11 + fCALOlength; //tstrip = 1064
495 } else {
496 tstrip = 0;
497 };
498 //
499 pre = -1;
500 //
501 for (Int_t strip=0; strip < 96; strip++){
502 //
503 // which is the pre for this strip?
504 //
505 if (strip%16 == 0) {
506 pre++;
507 };
508 //
509 if ( sec == 1 ) tstrip -= 11;
510 //
511 for (Int_t plane=0; plane < 11; plane++){
512
513 if ( sec == 0 || sec == 3 ) lpl = plane * 2;
514 if ( sec == 1 || sec == 2 ) lpl = (plane * 2) + 1;
515 //
516 // get the energy in GeV from the simulation for that strip
517 //
518 ens = this->GetCaloErel(sec,plane,strip);
519 //
520 // convert it into ADC channels
521 //
522 adcsig = int(ens*fCalomip[l][lpl][strip]/CALOGeV2MIPratio);
523 //
524 // sum baselines
525 //
526 adcbase = (UInt_t)fcalbase[sec][plane][pre];
527 //
528 // add noise and pedestals
529 //
530 pedestal = fcalped[sec][plane][strip];
531 rms = fcalrms[sec][plane][strip]/4.;
532 //
533 // Add random gaussian noise of RMS rms and Centered in the pedestal
534 //
535 pedenoise = gRandom->Gaus((Double_t)pedestal,(Double_t)rms);
536 //
537 // Sum all contribution
538 //
539 adc = adcsig + adcbase + (Int_t)round(pedenoise);
540 //
541 // Signal saturation
542 //
543 if ( adc > 0x7FFF ) adc = 0x7FFF;
544 //
545 // save value
546 //
547 if ( sec == 1 ){
548 DataCALO[tstrip] = adc;
549 tstrip++;
550 } else {
551 DataCALO[fCALOlength] = adc;
552 };
553 fCALOlength++;
554 //
555 };
556 //
557 if ( sec == 1 ) tstrip -= 11;
558 //
559 };
560 //
561 // here we calculate and save the CRC
562 //
563 Short_t CRC = 0;
564 for (UInt_t i=0; i<(fCALOlength-fSecPointer); i++){
565 CRC=crc(CRC,DataCALO[i+fSecPointer]);
566 };
567 DataCALO[fCALOlength] = (UShort_t)CRC;
568 fCALOlength++;
569 };
570
571 for (Int_t i = 0; i<fCALOlength; i++) fData.push_back(DataCALO[i]);
572 //reverse(fData.begin(),fData.end());
573 //for (Int_t i = 0; i<fCALOlength; i++) fData.push_back(DataCALO[fCALOlength-i-1]);
574 }
575
576
577 Float_t PamVMCCaloDig::GetCaloErel(Int_t sec, Int_t plane, Int_t strip){
578
579 // determine plane and strip
580
581 Int_t mplane = 0;
582
583 switch (sec){
584 case 0: //3
585 mplane = plane * 4 ; // it must be 0, 4, 8, ... (+1) from plane = 0, 11 YO
586 break;
587 case 1: //2
588 mplane = plane * 4 + 2; // it must be 2, 6, 10, ... (+1) from plane = 0, 11 YE
589 break;
590 case 2: //1
591 mplane = plane * 4 + 3 ; // it must be 3, 7, 11, ... (+1) from plane = 0, 11 XO
592 break;
593 case 3: //0
594 mplane = plane * 4 + 1 ; // it must be 1, 5, 9, ... (+1) from plane = 0, 11 XE
595 break;
596 }
597
598
599 Float_t Erel = 0.;
600
601 // search energy release in VMC output
602 if (fhc){
603 PamVMCDetectorHit * hit = (PamVMCDetectorHit*)fhc->At(mplane*96+strip);
604 if (hit) Erel=hit->GetEREL();
605 // if(Erel) cout<<"sec:"<<sec<<" plane:"<<plane<<" mplane:"<<mplane<<" AT:"<<mplane*96+strip<<" POS:"<<hit->GetPOS()<<endl;
606 } //else cout<<"CALO HIT Collection pointer not found"<<endl;
607
608 // if ((sec == 2) && (plane == 0)) Erel = 0.01; else Erel = 0.;
609
610 return Erel;
611 }
612
613
614
615 void PamVMCCaloDig::DigitizeCaloCompress(){
616
617 cout<<"!!!WARNING PamVMCCaloDig COMPRESS MODE STILL NOT IMPLEMENTED!!!"<<endl;
618
619 this->DigitizeCaloRaw();
620 return;
621
622 }
623
624 void PamVMCCaloDig::DigitizeCaloFull(){
625
626 cout<<"!!!WARNING PamVMCCaloDig FULL MODE STILL NOT IMPLEMENTED!!!"<<endl;
627
628 this->DigitizeCaloRaw();
629 return;
630
631 }

  ViewVC Help
Powered by ViewVC 1.1.23