/[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.5 - (show annotations) (download)
Fri Jun 12 18:39:11 2009 UTC (15 years, 5 months ago) by pam-rm2
Branch: MAIN
CVS Tags: v1r0
Changes since 1.1: +1 -1 lines
- Introduced user-defined names of output files and random seeds number.
Users can do it use options of PamVMCApplication constructor:
PamVMCApplication(const char* name,  const char *title, const char*
filename="pamtest", Int_t seed=0).
The Random object that I use is TRandom3 object which has astronomical
large period (in case of default initialization 0). All random generators
in the code use this object by calling of gRandom singleton which keeps
it.

- Corrected TOF digitization routine. No problems with TDC hits due to
hadronic interactions anymore.

- Some small changes was done to compile code under Root 5.23. +
geant4_vmc v. 2.6 without any warnings

- Some classes of PamG4RunConfiguartion was changed for geant4_vmc v.
2.6.Some obsolete classes was deleted as soon as developers implemented
regions.

- Navigation was changed from "geomRootToGeant4" to "geomRoot", because on
VMC web page written that as soon as Geant4 has no option ONLY/MANY
translation of overlapped geometry to Geant4 through VGM could be wrong.
I'd like to stay with Root navigation:
http://root.cern.ch/root/vmc/Geant4VMC.html. This should be default
option.

- New Tracker digitization routine written by Sergio was implemented

- PamVMC again became compatible with geant4_vmc v.2.5 and ROOT 5.20.
 The problem was that ROOT developers introduced in TVirtualMC class a new
method SetMagField and new base class:TVirtualMagField from which
user-defined classes shoukd be derived

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 = frandom->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