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

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

Parent Directory Parent Directory | Revision Log Revision Log

Revision 1.5 - (hide 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 nikolas 1.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 pam-rm2 1.5 pedenoise = frandom->Gaus((Double_t)pedestal,(Double_t)rms);
536 nikolas 1.1 //
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