/[PAMELA software]/trieste/pamVMC/cal/src/PamVMCCaloDig.cxx_my
ViewVC logotype

Contents of /trieste/pamVMC/cal/src/PamVMCCaloDig.cxx_my

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.1.1.1 - (show annotations) (download) (vendor branch)
Wed Mar 4 12:51:26 2009 UTC (15 years, 11 months ago) by pamelats
Branch: MAIN, pamVMC
CVS Tags: start, v0r00, HEAD
Changes since 1.1: +0 -0 lines
Test pamVMC

1 #include "PamVMCCaloDig.h"
2 #include <TTree.h>
3 #include "CalibCalPedEvent.h"
4
5 ClassImp(PamVMCCaloDig)
6
7 extern "C"{
8 short crc(short, short);
9 };
10
11
12 void PamVMCCaloDig::LoadCalib(){
13
14 cout<<"Loading Tracker Calibrations..."<<endl;
15
16
17
18 fhc =(TClonesArray*)fhitscolmap.GetValue("CAST");
19 if (!fhc) cout <<"!!!WARNING Hit Collection of CAL not found by PamVMCCaloDig!!!!"<<endl;
20
21 Int_t GivenCaloCalib = 0; //@@@ should be given as input par @@@
22
23 Int_t time=3;
24 Int_t type=101;
25
26 fdberr = fsql->Query_GL_PARAM(time,type);
27
28 if(fdberr<0){
29 cout<<"No such record in DB for CAL: time="<<time<<" type="<<type<<endl;
30 ThrowCalFileWarning("CAL");
31
32 } else {
33
34 fquery.str("");
35 fquery << fsql->GetPAR()->PATH.Data() << "/";
36 fquery << fsql->GetPAR()->NAME.Data();
37
38 ThrowCalFileUsage("CAL",fquery.str().c_str());
39
40 FILE *f;
41 f = fopen(fquery.str().c_str(),"rb");
42
43 if(f==NULL) ThrowCalFileWarning("CAL");
44 else{
45 memset(fCalomip,0,4224*sizeof(fCalomip[0][0][0]));
46
47 for (Int_t m = 0; m < 2 ; m++ ){
48 for (Int_t k = 0; k < 22; k++ ){
49 for (Int_t l = 0; l < 96; l++ ){
50 fread(&fCalomip[m][k][l],sizeof(fCalomip[m][k][l]),1,f);
51 };
52 };
53 };
54 fclose(f);
55
56 // determine which calibration has to be used
57 // and load it for each section
58
59
60 UInt_t idcalib, calibno;
61 UInt_t utime = 0;
62 TString calname;
63 for (UInt_t s=0; s<4; s++){
64
65 // clear calo calib variables for section s
66
67 ClearCaloCalib(s);
68
69 if ( GivenCaloCalib ){
70
71 // a time has been given as input on the command line
72 // so retrieve the calibration that preceed that time
73
74 fsql->Query_GL_CALO_CALIB(GivenCaloCalib,utime,s);
75
76 calibno = fsql->GetCaloCalib()->EV_ROOT;
77 idcalib = fsql->GetCaloCalib()->ID_ROOT_L0;
78
79 // determine path and name and entry of the calibration file
80
81 printf("\n");
82 printf(" ** SECTION %i **\n",s);
83
84 fsql->Query_GL_ROOT(idcalib);
85
86 fquery.str("");
87 fquery << fsql->GetROOT()->PATH.Data() << "/";
88 fquery << fsql->GetROOT()->NAME.Data();
89
90 calname = (TString)fquery.str().c_str();
91
92 printf("\n Section %i : using file %s calibration at entry %i: \n",s,calname.Data(),calibno);
93
94 } else {
95
96 fdberr = 0;
97 fdberr = fsql->Query_GL_PARAM(1,104);
98
99 fquery.str("");
100 fquery << fsql->GetPAR()->PATH.Data() << "/";
101 fquery << fsql->GetPAR()->NAME.Data();
102
103 printf("\n Section %i : using default calorimeter calibration: \n %s \n",s,fquery.str().c_str());
104
105 calname = (TString)fquery.str().c_str();
106 calibno = s;
107
108 };
109
110 // load calibration variables in memory
111
112 CaloLoadCalib(s,calname,calibno);
113
114 };
115 //
116 // at this point we have in memory the calorimeter calibration
117 // and we can save it to disk in the correct format and use it to digitize the data
118 //
119
120 DigitizeCALOCALIB();
121 }
122 }
123 }
124
125
126 void PamVMCCaloDig::ClearCaloCalib(Int_t s){
127
128 fcstwerr[s] = 0;
129 fcperror[s] = 0.;
130 for ( Int_t d=0 ; d<11 ;d++ ){
131 Int_t pre = -1;
132 for ( Int_t j=0; j<96 ;j++){
133 if ( j%16 == 0 ) pre++;
134 fcalped[s][d][j] = 0.;
135 fcstwerr[s] = 0.;
136 fcperror[s] = 0.;
137 fcalgood[s][d][j] = 0.;
138 fcalthr[s][d][pre] = 0.;
139 fcalrms[s][d][j] = 0.;
140 fcalbase[s][d][pre] = 0.;
141 fcalvar[s][d][pre] = 0.;
142 };
143 };
144 return;
145 }
146
147 Int_t PamVMCCaloDig::CaloLoadCalib(Int_t s,TString calname, UInt_t calibno){
148
149 UInt_t e = 0;
150 switch(s){
151 case 0: e = 0;
152 break;
153 case 1: e = 2;
154 break;
155 case 2: e = 3;
156 break;
157 case 3: e = 1;
158 break;
159 default:
160 break;
161 }
162 fcfile.open(calname.Data());
163 ThrowCalFileUsage("CAL",calname.Data());
164 if ( !fcfile ){
165 ThrowCalFileWarning("CAL");
166 return(-107);
167 };
168 fcfile.close();
169
170 fcrfile = new TFile(calname.Data());
171 if ( !fcrfile ){
172 ThrowCalFileWarning("CAL");
173 return(-108);
174 }
175 TTree *tr = (TTree*)fcrfile->Get("CalibCalPed");
176 if ( !tr ){
177 cout<<"!!!WARNING Tree CalibCalPed in file "<<calname.Data()
178 <<" was NOT found!!!"<<endl;
179 return(-109);
180 }
181
182 TBranch *calo = tr->GetBranch("CalibCalPed");
183
184 pamela::CalibCalPedEvent *ce = 0;
185 tr->SetBranchAddress("CalibCalPed", &ce);
186
187 Long64_t ncalibs = calo->GetEntries();
188
189 if ( !ncalibs ){
190 cout<<"!!!WARNING No entries in calibration file!!!"<<endl;
191 return(-110);
192 }
193
194 calo->GetEntry(calibno);
195
196 if (ce->cstwerr[s] != 0 && ce->cperror[s] == 0 ) {
197 fcstwerr[s] = ce->cstwerr[s];
198 fcperror[s] = ce->cperror[s];
199 for ( Int_t d=0 ; d<11 ;d++ ){
200 Int_t pre = -1;
201 for ( Int_t j=0; j<96 ;j++){
202 if ( j%16 == 0 ) pre++;
203 fcalped[s][d][j] = ce->calped[e][d][j];
204 fcalgood[s][d][j] = ce->calgood[e][d][j];
205 fcalthr[s][d][pre] = ce->calthr[e][d][pre];
206 fcalrms[s][d][j] = ce->calrms[e][d][j];
207 fcalbase[s][d][pre] = ce->calbase[e][d][pre];
208 fcalvar[s][d][pre] = ce->calvar[e][d][pre];
209 };
210 };
211 } else {
212 printf(" CALORIMETER - ERROR: problems finding a good calibration in this file! \n\n ");
213 fcrfile->Close();
214 return(-111);
215 };
216 fcrfile->Close();
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 Double_t pedenoise;
430 Float_t rms = 0.;
431 Float_t pedestal = 0.;
432 UInt_t tmpadc[11];
433 USBuffer tempdig;
434
435 static const Float_t CALOGeV2MIPratio = 0.0001059994;
436
437 // Header of the four sections
438
439 fSecCalo[0] = 0xEA08; // XE
440 fSecCalo[1] = 0xF108; // XO
441 fSecCalo[2] = 0xF608; // YE
442 fSecCalo[3] = 0xED08; // YO
443
444 // length of the data is 0x0428 in RAW mode
445
446 fSecCaloLength[0] = 0x0428; // XE
447 fSecCaloLength[1] = 0x0428; // XO
448 fSecCaloLength[2] = 0x0428; // YE
449 fSecCaloLength[3] = 0x0428; // YO
450
451 fData.clear();
452
453 for (Int_t sec=0; sec < 4; sec++){
454
455 // sec = 0 -> XE 1 -> XO 2-> YE 3 -> YO
456
457 tempdig.clear();
458
459 // First of all we have section header and packet length
460 if (sec !=1){
461 tempdig.push_back(fSecCalo[sec]);//1
462 tempdig.push_back(fSecCaloLength[sec]);//2
463 }
464
465 // selftrigger coincidences - in the future we should add here some code
466 // to simulate timing response of pre-amplifiers
467 // For now from application point of view we have 1 hit = 1 strip (nHits <=4224)
468 // It's possible to introduce 1 hit = 1 particle like in other detectors... if need
469
470 for (Int_t autoplane=0; autoplane < 7; autoplane++) tempdig.push_back(0x0000); //3-8
471
472 // here comes data
473 // Section XO is read in the opposite direction respect to the others
474
475 pre = -1;
476
477 //l=(sec<2)?1:0;
478 l = 0; // XE and XO are Y planes
479 if ( sec < 2 ) l = 1; // while YE and YO are X planes
480
481 for (Int_t strip=0; strip < 96; strip++){
482
483 // which is the pre for this strip?
484
485 if (strip%16 == 0) pre++;
486
487 for (Int_t plane=0; plane < 11; plane++){
488
489 switch (sec){
490 case 0:
491 case 3:
492 lpl = plane * 2;
493 break;
494 case 1:
495 case 2:
496 lpl = (plane * 2) + 1;
497 break;
498 default:
499 break;
500 }
501
502 // get the energy in GeV from the simulation for that strip
503
504 ens = GetCaloErel(sec,plane,strip);
505
506 // convert it into ADC channels
507
508 adcsig = Int_t(ens*fCalomip[l][lpl][strip]/CALOGeV2MIPratio);
509
510 // sum baselines
511
512 adcbase = (UInt_t)fcalbase[sec][plane][pre];
513
514 // add noise and pedestals
515
516 pedestal = fcalped[sec][plane][strip];
517 rms = fcalrms[sec][plane][strip]/4.;
518
519 // Add random gaussian noise of RMS rms and Centered in the pedestal
520
521 pedenoise = gRandom->Gaus((Double_t)pedestal,(Double_t)rms);
522
523 // Sum all contribution
524
525 adc = adcsig + adcbase + (Int_t)round(pedenoise);
526
527 // Signal saturation
528
529 if ( adc > 0x7FFF ) adc = 0x7FFF;
530
531 // save value
532 if (adcsig>0){
533 cout <<"STRIP HIT"<<endl
534 <<"SEC: "<<sec<<" PLANE: "<<plane<<" STRIP:"<<strip<<endl
535 <<"adcsig: "<< adcsig<<" adc: "<<adc<<endl;
536 }
537 if ( sec != 1 ) tempdig.push_back(adc);
538 else tmpadc[plane]= adc;
539
540 }
541
542 if ( sec == 1 ) for (Int_t i=10; i>=0; i--) tempdig.push_back(tmpadc[i]);
543 }
544
545 if ( sec == 1 ){
546 tempdig.push_back(fSecCaloLength[1]);//2
547 tempdig.push_back(fSecCalo[1]);//1
548 reverse(tempdig.begin(), tempdig.end());
549 }
550 // here we calculate and save the CRC
551 //++++++ COPY TO DIGITIZER'S BUFFER +++//
552 Short_t CRC = 0;
553
554 USBuffer::const_iterator p = tempdig.begin();
555 while( p!= tempdig.end() ){
556 CRC=crc(CRC,*p);
557 fData.push_back(*p);
558 p++;
559 }
560 fData.push_back(CRC);
561 }
562 }
563
564
565 Float_t PamVMCCaloDig::GetCaloErel(Int_t sec, Int_t plane, Int_t strip){
566
567 // determine plane and strip
568
569 Int_t mplane = 0;
570
571 switch (sec){
572 case 0:
573 mplane = plane * 4 + 1; // it must be 0, 4, 8, ... (+1) from plane = 0, 11
574 break;
575 case 1:
576 mplane = plane * 4 + 2 + 1; // it must be 2, 6, 10, ... (+1) from plane = 0, 11
577 break;
578 case 2:
579 mplane = plane * 4 + 3 + 1; // it must be 3, 7, 11, ... (+1) from plane = 0, 11
580 break;
581 case 3:
582 mplane = plane * 4 + 1 + 1; // it must be 1, 5, 9, ... (+1) from plane = 0, 11
583 break;
584 }
585
586
587 Float_t Erel = 0.;
588
589
590 // search energy release in gpamela output
591 if (fhc){
592 PamVMCDetectorHit *hit = (PamVMCDetectorHit*)fhc->At((mplane-1)*96+strip);
593 if (hit) Erel=hit->GetEREL();
594 } else cout<<"CALO HIT Collection pointer not found"<<endl;
595
596 //if (sec == 3) Erel = 0.01; else Erel = 0.;
597
598 return Erel;
599 }
600
601
602
603 void PamVMCCaloDig::DigitizeCaloCompress(){
604
605 cout<<"!!!WARNING PamVMCCaloDig COMPRESS MODE STILL NOT IMPLEMENTED!!!"<<endl;
606
607 this->DigitizeCaloRaw();
608 return;
609
610 }
611
612 void PamVMCCaloDig::DigitizeCaloFull(){
613
614 cout<<"!!!WARNING PamVMCCaloDig FULL MODE STILL NOT IMPLEMENTED!!!"<<endl;
615
616 this->DigitizeCaloRaw();
617 return;
618
619 }

  ViewVC Help
Powered by ViewVC 1.1.23