/[PAMELA software]/DarthVader/CalorimeterLevel2/src/CaloProcessing.cpp
ViewVC logotype

Annotation of /DarthVader/CalorimeterLevel2/src/CaloProcessing.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.13 - (hide annotations) (download)
Tue Feb 13 16:55:08 2007 UTC (17 years, 10 months ago) by mocchiut
Branch: MAIN
Changes since 1.12: +29 -20 lines
Calorimeter: --no-crosstalk option added

1 mocchiut 1.1 /**
2     * \file src/CaloProcessing.cpp
3     * \author Emiliano Mocchiutti
4     **/
5     //
6     // C/C++ headers
7     //
8     #include <sstream>
9     #include <fstream>
10     //
11     // ROOT headers
12     //
13     #include <TTree.h>
14     #include <TBranch.h>
15     #include <TFile.h>
16     #include <TObject.h>
17     //
18     // YODA headers
19     //
20     #include <PamelaRun.h>
21     #include <physics/calorimeter/CalorimeterEvent.h>
22     #include <CalibCalPedEvent.h>
23     //
24     //
25     //
26     #include <GLTables.h>
27     //
28     // this package headers
29     //
30     #include <delay.h>
31     #include <CaloProcessing.h>
32     //
33     //
34     // Declaration of the core fortran routines
35     //
36     #define calol2cm calol2cm_
37     extern "C" int calol2cm();
38     #define calol2tr calol2tr_
39     extern "C" int calol2tr();
40     //
41     using namespace std;
42     //
43     //
44     // Public methods
45     //
46    
47     CaloProcessing::~CaloProcessing(){
48 mocchiut 1.3 delete de;
49 mocchiut 1.1 delete this;
50     }
51    
52     CaloProcessing::CaloProcessing(){
53     //
54     extern struct FlCaLevel1 clevel1_;
55     extern struct FlCaLevel2 clevel2_;
56     clevel1 = &clevel1_;
57     clevel2 = &clevel2_;
58     //
59     trkseqno = 0;
60     ClearStructs();
61 mocchiut 1.3 //
62     memset(dexy, 0, 2*22*96*sizeof(Float_t));
63     memset(dexyc, 0, 2*22*96*sizeof(Float_t));
64     memset(mip, 0, 2*22*96*sizeof(Float_t));
65     memset(base, 0, 2*22*6*sizeof(Float_t));
66     memset(sbase, 0, 2*22*6*sizeof(Float_t));
67 mocchiut 1.1 calopar1 = true;
68     calopar2 = true;
69 mocchiut 1.6 calopar3 = true;
70 mocchiut 1.13 crosst = true;
71 mocchiut 1.4 ftcalopar1 = 0;
72     ttcalopar1 = 0;
73     ftcalopar2 = 0;
74     ttcalopar2 = 0;
75 mocchiut 1.6 ftcalopar3 = 0;
76     ttcalopar3 = 0;
77 mocchiut 1.1 }
78    
79 mocchiut 1.13 void CaloProcessing::SetCrossTalk(Bool_t ct){
80     crosst = ct;
81     };
82    
83 mocchiut 1.1 /**
84     * Initialize CaloProcessing object
85     **/
86 mocchiut 1.4 void CaloProcessing::ProcessingInit(TSQLServer *dbc, UInt_t hs, Int_t &sgnl, TTree *l0tree, Bool_t isdeb, Bool_t isverb){
87 mocchiut 1.1 //
88     debug = isdeb;
89     verbose = isverb;
90     //
91     l0tr=(TTree*)l0tree;
92     de = new pamela::calorimeter::CalorimeterEvent();
93     l0calo = (TBranch*)l0tr->GetBranch("Calorimeter");
94     l0tr->SetBranchAddress("Calorimeter", &de);
95     //
96     trkseqno = 0;
97     ClearStructs();
98     //
99     GL_CALO_CALIB *glcalo = new GL_CALO_CALIB();
100     //
101     sgnl = 0;
102 mocchiut 1.5 UInt_t uptime = 0;
103 mocchiut 1.1 //
104     for (Int_t s = 0; s < 4; s++){
105 mocchiut 1.4 idcalib[s] = 0;
106     fromtime[s] = 0;
107     totime[s] = 0;
108 mocchiut 1.1 calibno[s] = 0;
109     ClearCalibVals(s);
110 mocchiut 1.4 //
111 mocchiut 1.5 sgnl = glcalo->Query_GL_CALO_CALIB(hs,uptime,s,dbc);
112 mocchiut 1.1 if ( sgnl < 0 ){
113     if ( verbose ) printf(" CALORIMETER - ERROR: error from GLTables\n");
114     return;
115     };
116 mocchiut 1.4 //
117     idcalib[s] = glcalo->ID_ROOT_L0;
118 mocchiut 1.1 fromtime[s] = glcalo->FROM_TIME;
119 mocchiut 1.5 if ( glcalo->TO_TIME < hs ){ // calibration is corrupted and we are using the one that preceed the good one
120     totime[s] = uptime;
121     } else {
122     totime[s] = glcalo->TO_TIME;
123     };
124 mocchiut 1.4 calibno[s] = glcalo->EV_ROOT;
125     //
126     if ( totime[s] == 0 ){
127 mocchiut 1.1 if ( verbose ) printf(" CALORIMETER - WARNING: data with no associated calibration\n");
128     ClearCalibVals(s);
129     sgnl = 100;
130     };
131     };
132     //
133     // determine path and name and entry of the calibration file
134     //
135     GL_ROOT *glroot = new GL_ROOT();
136     if ( verbose ) printf("\n");
137     for (Int_t s = 0; s < 4; s++){
138     if ( verbose ) printf(" ** SECTION %i **\n",s);
139 mocchiut 1.4 if ( totime[s] > 0 ){
140     //
141 mocchiut 1.1 sgnl = glroot->Query_GL_ROOT(idcalib[s],dbc);
142     if ( sgnl < 0 ){
143     if ( verbose ) printf(" CALORIMETER - ERROR: error from GLTables\n");
144     return;
145     };
146 mocchiut 1.4 //
147 mocchiut 1.1 stringstream name;
148     name.str("");
149     name << glroot->PATH.Data() << "/";
150     name << glroot->NAME.Data();
151     //
152     fcalname[s] = (TString)name.str().c_str();
153 mocchiut 1.4 if ( verbose ) printf(" - runheader at time %u. From time %u to time %u \n use file %s \n calibration at entry %i \n\n",hs,fromtime[s],totime[s],fcalname[s].Data(),calibno[s]);
154 mocchiut 1.1 } else {
155 mocchiut 1.4 if ( verbose ) printf(" - runheader at time %u. NO CALIBRATION INCLUDE THE RUNHEADER! ",hs);
156 mocchiut 1.1 };
157     sgnl = LoadCalib(s);
158     if ( sgnl ) break;
159     };
160     //
161 mocchiut 1.2 delete glcalo;
162     delete glroot;
163     //
164 mocchiut 1.1 return;
165     //
166     }
167    
168 mocchiut 1.4 Int_t CaloProcessing::ChkCalib(TSQLServer *dbc, UInt_t atime){
169 mocchiut 1.1 Int_t sgnl = 0;
170     for ( Int_t s = 0; s < 4; s++){
171     if ( atime > totime[s] ){
172 mocchiut 1.7 if ( !dbc->IsConnected() ) throw -116;
173 mocchiut 1.1 sgnl = Update(dbc,atime,s);
174     if ( sgnl < 0 ) return(sgnl);
175     };
176     };
177     return(sgnl);
178     }
179    
180 mocchiut 1.4 Int_t CaloProcessing::ChkParam(TSQLServer *dbc, UInt_t runheader){
181 mocchiut 1.1 stringstream calfile;
182 mocchiut 1.6 stringstream bmfile;
183 mocchiut 1.1 stringstream aligfile;
184     Int_t error = 0;
185     FILE *f = 0;
186 mocchiut 1.6 ifstream badfile;
187 mocchiut 1.1 GL_PARAM *glparam = new GL_PARAM();
188     //
189     if ( calopar1 || ( ttcalopar1 != 0 && ttcalopar1 < runheader ) ){
190     calopar1 = false;
191     //
192     // determine where I can find calorimeter ADC to MIP conversion file
193     //
194     if ( verbose ) printf(" Querying DB for calorimeter parameters files...\n");
195     //
196     error = 0;
197 mocchiut 1.4 error = glparam->Query_GL_PARAM(runheader,101,dbc);
198 mocchiut 1.1 if ( error < 0 ) return(error);
199     //
200     calfile.str("");
201     calfile << glparam->PATH.Data() << "/";
202     calfile << glparam->NAME.Data();
203     ftcalopar1 = glparam->FROM_TIME;
204     ttcalopar1 = glparam->TO_TIME;
205     //
206     if ( verbose ) printf("\n Using ADC to MIP conversion file: \n %s \n",calfile.str().c_str());
207     f = fopen(calfile.str().c_str(),"rb");
208     if ( !f ){
209     if ( verbose ) printf(" CALORIMETER - ERROR: no ADC to MIP file!\n");
210     return(-105);
211     };
212     //
213     for (Int_t m = 0; m < 2 ; m++ ){
214     for (Int_t k = 0; k < 22; k++ ){
215     for (Int_t l = 0; l < 96; l++ ){
216     fread(&mip[m][k][l],sizeof(mip[m][k][l]),1,f);
217     };
218     };
219     };
220     fclose(f);
221     };
222     //
223     if ( calopar2 || ( ttcalopar2 != 0 && ttcalopar2 < runheader ) ){
224     calopar2 = false;
225     //
226     // determine where I can find calorimeter alignment file
227     //
228     //
229     error = 0;
230 mocchiut 1.4 error = glparam->Query_GL_PARAM(runheader,102,dbc);
231 mocchiut 1.1 if ( error < 0 ) return(error);
232     //
233     aligfile.str("");
234     aligfile << glparam->PATH.Data() << "/";
235     aligfile << glparam->NAME.Data();
236     ftcalopar2 = glparam->FROM_TIME;
237     ttcalopar2 = glparam->TO_TIME;
238     //
239 mocchiut 1.6 if ( verbose ) printf("\n Using alignment file: \n %s \n",aligfile.str().c_str());
240 mocchiut 1.1 f = fopen(aligfile.str().c_str(),"rb");
241     if ( !f ){
242     if ( verbose ) printf(" CALORIMETER - ERROR: no alignement file!\n");
243     return(-106);
244     };
245     //
246     fread(&clevel1->xalig,sizeof(clevel1->xalig),1,f);
247     if ( debug ) printf(" xalig = %f \n",clevel1->xalig);
248     fread(&clevel1->yalig,sizeof(clevel1->yalig),1,f);
249     if ( debug ) printf(" yalig = %f \n",clevel1->yalig);
250     fread(&clevel1->zalig,sizeof(clevel1->zalig),1,f);
251     if ( debug ) printf(" zalig = %f \n",clevel1->zalig);
252     fread(&clevel1->emin,sizeof(clevel1->emin),1,f);
253     if ( debug ) printf(" signal threshold = %f \n",clevel1->emin);
254     //
255     fclose(f);
256     };
257     //
258 mocchiut 1.6 // Load offline bad strip mask
259     //
260     if ( calopar3 || ( ttcalopar3 != 0 && ttcalopar3 < runheader ) ){
261     calopar3 = false;
262     //
263     // determine where I can find calorimeter alignment file
264     //
265     //
266     error = 0;
267     error = glparam->Query_GL_PARAM(runheader,103,dbc);
268     if ( error < 0 ) return(error);
269     //
270     bmfile.str("");
271     bmfile << glparam->PATH.Data() << "/";
272     bmfile << glparam->NAME.Data();
273     ftcalopar3 = glparam->FROM_TIME;
274     ttcalopar3 = glparam->TO_TIME;
275     //
276     if ( verbose ) printf("\n Using bad strip offline mask file: \n %s \n\n",bmfile.str().c_str());
277     badfile.open(bmfile.str().c_str());
278     if ( !badfile ){
279     if ( verbose ) printf(" CALORIMETER - ERROR: no bad strip offline mask file!\n");
280     return(-115);
281     };
282     //
283     Bool_t isdone = false;
284     Int_t bad = 0;
285     Int_t view = 1;
286     Int_t strip = 0;
287     Int_t plane = 21;
288     while ( !isdone ) {
289     badfile >> bad;
290     obadmask[view][plane][strip] = bad;
291     if ( debug && bad ) printf(" SETTING view %i plane %i strip %i BAD = %i \n",view,plane,strip,bad);
292     strip++;
293     if ( strip > 95 ){
294     strip = 0;
295     plane--;
296     if ( plane < 0 ){
297     plane = 21;
298     view--;
299     };
300     if ( view < 0 ) isdone = true;
301     };
302     };
303     //
304     badfile.close();
305     };
306     //
307 mocchiut 1.2 delete glparam;
308     //
309 mocchiut 1.1 return(0);
310     }
311    
312    
313    
314     void CaloProcessing::FindBaseRaw(Int_t l, Int_t m, Int_t pre){
315     Float_t minstrip = 100000.;
316     Float_t rms = 0.;
317     base[l][m][pre] = 0.;
318     for (Int_t e = pre*16; e < (pre+1)*16 ; e++){
319 mocchiut 1.6 if ( calgood[l][m][e] == 0. && obadmask[l][m][e] == 0 && dexy[l][m][e]-calped[l][m][e] < minstrip && dexy[l][m][e] > 0.) {
320 mocchiut 1.1 minstrip = dexy[l][m][e]-calped[l][m][e];
321     rms = calthr[l][m][pre];
322     };
323     };
324     if ( minstrip != 100000. ) {
325     Float_t strip6s = 0.;
326     for (Int_t e = pre*16; e < (pre+1)*16 ; e++){
327     if ( (dexy[l][m][e]-calped[l][m][e]) >= minstrip && (dexy[l][m][e]-calped[l][m][e]) <= (minstrip+rms) ) {
328     strip6s += 1.;
329     base[l][m][pre] += (dexy[l][m][e] - calped[l][m][e]);
330     };
331     //
332     // compression
333     //
334     if ( abs((int)(dexy[l][m][e]-calped[l][m][e])) <= (minstrip+rms) ) {
335     dexyc[l][m][e] = 0.;
336     } else {
337     dexyc[l][m][e] = dexy[l][m][e];
338     };
339     };
340     if ( strip6s >= 9. ){
341     Double_t arro = base[l][m][pre]/strip6s;
342     Float_t deci = 1000.*((float)arro - float(int(arro)));
343     if ( deci < 500. ) {
344     arro = double(int(arro));
345     } else {
346     arro = 1. + double(int(arro));
347     };
348     base[l][m][pre] = arro;
349     } else {
350     base[l][m][pre] = 31000.;
351     for (Int_t e = pre*16; e < (pre+1)*16 ; e++){
352     dexyc[l][m][e] = dexy[l][m][e];
353     };
354     };
355     } else {
356     base[l][m][pre] = 31000.;
357     };
358     }
359    
360     Int_t CaloProcessing::Calibrate(Int_t ei){
361     //
362     // get entry ei
363     //
364     l0calo->GetEntry(ei);
365     //
366     // if it was not a selftrigger event, could it ever been a selftrigger event? if so trigty = 3.
367     //
368     Int_t val = 0;
369     Int_t del = 1100;
370     if ( clevel2->trigty != 2. ){
371 mocchiut 1.10 Bool_t ck = false;
372 mocchiut 1.1 for (Int_t sec = 0; sec < 4; sec++){
373     val = (Int_t)de->calselftrig[sec][6];
374     del = delay(val);
375     if ( del < 1100 ){
376 mocchiut 1.10 clevel2->wartrig = 0.;
377 mocchiut 1.1 clevel2->trigty = 3.;
378 mocchiut 1.10 ck = true;
379 mocchiut 1.1 break;
380     };
381     };
382 mocchiut 1.10 if ( !ck ) clevel2->wartrig = 100.;
383     } else {
384     Bool_t ck = false;
385     for (Int_t sec = 0; sec < 4; sec++){
386     val = (Int_t)de->calselftrig[sec][6];
387     del = delay(val);
388     if ( del < 1100 ){
389     clevel2->wartrig = 0.;
390     ck = true;
391     };
392     };
393     if ( !ck ) clevel2->wartrig = 100.;
394 mocchiut 1.1 };
395     //
396     Int_t se = 5;
397     Int_t done = 0;
398     Int_t pre = -1;
399     Bool_t isCOMP = false;
400     Bool_t isFULL = false;
401     Bool_t isRAW = false;
402     Float_t ener;
403     Int_t doneb = 0;
404     Int_t donec = 0;
405     Int_t ck = 0;
406     Int_t ipre = 0;
407     Int_t ip[3] = {0};
408     Float_t base0, base1, base2;
409     base0 = 0.;
410     base1 = 0.;
411     base2 = 0.;
412     Float_t qpre[6] = {0.,0.,0.,0.,0.,0.};
413     Float_t ene[96];
414     Int_t chdone[4] = {0,0,0,0};
415     Int_t pe = 0;
416     //
417 mocchiut 1.6 Float_t ener0 = 0.;
418     Float_t cbase0 = 0.;
419     Bool_t pproblem = false;
420     //
421     Float_t tim = 0.;
422     Int_t plo = 0;
423     Int_t fbi = 0;
424     Int_t cle = 0;
425     //
426 mocchiut 1.1 // run over views and planes
427     //
428     for (Int_t l = 0; l < 2; l++){
429     for (Int_t m = 0; m < 22; m++){
430     //
431     // determine the section number
432     //
433     se = 5;
434     if (l == 0 && m%2 == 0) se = 3;
435     if (l == 0 && m%2 != 0) se = 2;
436     if (l == 1 && m%2 == 0) se = 1;
437     if (l == 1 && m%2 != 0) se = 0;
438     //
439     // determine what kind of event we are going to analyze
440     //
441     isCOMP = false;
442     isFULL = false;
443     isRAW = false;
444     if ( de->stwerr[se] & (1 << 16) ) isCOMP = true;
445     if ( de->stwerr[se] & (1 << 17) ) isFULL = true;
446     if ( de->stwerr[se] & (1 << 3) ) isRAW = true;
447     if ( !chdone[se] ){
448     //
449     // check for any error in the event
450     //
451     clevel2->crc[se] = 0;
452     if ( de->perror[se] == 132 ){
453     clevel2->crc[se] = 1;
454     pe++;
455     };
456     clevel2->perr[se] = 0;
457     if ( de->perror[se] != 0 ){
458     clevel2->perr[se] = 1;
459     pe++;
460     };
461     clevel2->swerr[se] = 0;
462     for (Int_t j = 0; j < 7 ; j++){
463     if ( (j != 3) && (de->stwerr[se] & (1 << j)) ){
464     clevel2->swerr[se] = 1;
465     pe++;
466     };
467     };
468     chdone[se] = 1;
469     };
470     if ( clevel2->crc[se] == 0 && (clevel1->good2 == 1 || clevel2->trigty >= 2) ){
471     pre = -1;
472     //
473     for (Int_t nn = 0; nn < 96; nn++){
474     ene[nn] = 0.;
475     dexy[l][m][nn] = de->dexy[l][m][nn] ;
476     dexyc[l][m][nn] = de->dexyc[l][m][nn] ;
477     };
478     //
479     // run over preamplifiers
480     //
481     pre = -1;
482 mocchiut 1.6 cbase0 = 0.;
483 mocchiut 1.1 for (Int_t i = 0; i < 3; i++){
484     for (Int_t j = 0; j < 2; j++){
485     pre = j + i*2;
486     //
487     // baseline check and calculation
488     //
489     if ( !isRAW ) {
490     base[l][m][pre] = de->base[l][m][pre] ;
491 mocchiut 1.6 cbase0 += base[l][m][pre];
492 mocchiut 1.1 } else {
493     //
494     // if it is a raw event and we haven't checked
495     // yet, calculate the baseline.
496     //
497     FindBaseRaw(l,m,pre);
498 mocchiut 1.6 cbase0 += base[l][m][pre];
499 mocchiut 1.1 };
500     };
501     };
502     //
503     // run over strips
504     //
505     pre = -1;
506 mocchiut 1.6 ener0 = 0.;
507 mocchiut 1.1 for (Int_t i = 0 ; i < 3 ; i++){
508     ip[i] = 0;
509     for (Int_t n = i*32 ; n < (i+1)*32 ; n++){
510     if (n%16 == 0) {
511     ck = 0;
512     done = 0;
513     doneb = 0;
514     donec = 0;
515     pre++;
516 mocchiut 1.3 qpre[pre] = 0.;
517 mocchiut 1.1 };
518     //
519     // baseline check and calculation
520     //
521     // no suitable new baseline, use old ones!
522     //
523     if ( !done ){
524     if ( (base[l][m][pre] == 31000. || base[l][m][pre] == 0.) ){
525     ck = 1;
526     if (pre%2 == 0) {
527     ip[i] = pre + 1;
528     } else {
529     ip[i] = pre - 1;
530     };
531     if ( (base[l][m][ip[i]] == 31000. || base[l][m][ip[i]] == 0.) ){
532     //
533     ck = 2;
534     if ( sbase[l][m][pre] == 31000. || sbase[l][m][pre] == 0. ) {
535     ck = 3;
536     };
537     };
538     done = 1;
539     };
540     };
541     //
542     // CALIBRATION ALGORITHM
543     //
544     if ( !doneb ){
545     switch (ck) {
546     case 0:
547     base0 = base[l][m][pre];
548     base2 = calbase[l][m][pre];
549     break;
550     case 1:
551     base0 = base[l][m][ip[i]];
552     base2 = calbase[l][m][ip[i]];
553     break;
554     case 2:
555     base0 = sbase[l][m][pre];
556     base2 = calbase[l][m][pre];
557     break;
558     case 3:
559     base0 = calbase[l][m][pre];
560     base2 = calbase[l][m][pre];
561     break;
562     };
563     base1 = calbase[l][m][pre];
564     doneb = 1;
565     };
566     ener = dexyc[l][m][n];
567 mocchiut 1.6 ener0 += ener;
568 mocchiut 1.1 clevel1->estrip[n][m][l] = 0.;
569     if ( base0>0 && base0 < 30000. ){
570     if ( !donec && (base0 - base1 + base2) != 0. ){
571     sbase[l][m][pre] = base0 - base1 + base2;
572     donec = 1;
573     };
574     if ( ener > 0. ){
575     clevel1->estrip[n][m][l] = (ener - calped[l][m][n] - base0 - base1 + base2)/mip[l][m][n] ;
576     //
577     // OK, now in estrip we have the energy deposit in MIP of all the strips for this event (at the end of loops of course)
578     //
579     qpre[pre] += clevel1->estrip[n][m][l];
580     };
581     };
582     };
583 mocchiut 1.13 if ( crosst ){
584     if (ck == 1){
585     if (ip[i]%2 == 0) {
586     ipre = ip[i] + 1;
587     } else {
588     ipre = ip[i] - 1;
589     };
590     for (Int_t j = ipre*16 ; j < (ipre+1)*16 ; j++){
591     clevel1->estrip[j][m][l] += (qpre[ipre] - qpre[ip[i]]) * 0.00478;
592     };
593 mocchiut 1.1 };
594 mocchiut 1.13 if (ck == 2){
595     for (Int_t j = i*32 ; j < (i+1)*32 ; j++){
596     ipre = j/16 + 1;
597     clevel1->estrip[j][m][l] += qpre[ipre] * 0.00478;
598     };
599 mocchiut 1.1 };
600     };
601 mocchiut 1.6 };
602     //
603 mocchiut 1.11 if ( ener0 == 0. && cbase0 == 0. && !pproblem && clevel2->perr[se] == 0){
604 mocchiut 1.12 if ( verbose ) printf(" L0 entry %i : calorimeter power problems! event marked as bad perr %f swerr %X view %i plane %i \n",ei,de->perror[se],de->stwerr[se],l,m);
605 mocchiut 1.6 pproblem = true;
606     pe++;
607     };
608 mocchiut 1.1 //
609     Int_t j4 = -4;
610     Int_t jjj = -3;
611     Int_t jj = -2;
612     for (Int_t j = 0 ; j < 100 ; j++){
613     jj++;
614     jjj++;
615     j4++;
616     if ( j < 96 ) ene[j] = clevel1->estrip[j][m][l];
617 mocchiut 1.13 if ( crosst ){
618     if ( jj >= 0 && jj < 96 ){
619     if ( jj != 0 && jj != 32 && jj != 64 ) ene[jj-1] += -clevel1->estrip[jj][m][l] * 0.01581;
620     if ( jj != 31 && jj != 63 && jj != 95 ) ene[jj+1] += -clevel1->estrip[jj][m][l] * 0.01581;
621     };
622     if ( jjj >= 0 && jjj < 96 ){
623     if ( jjj != 0 && jjj != 32 && jjj != 64 ) clevel1->estrip[jjj-1][m][l] += -ene[jjj] * 0.01581;
624     if ( jjj != 31 && jjj != 63 && jjj != 95 ) clevel1->estrip[jjj+1][m][l] += -ene[jjj] * 0.01581;
625     };
626 mocchiut 1.1 };
627     if ( j4 >= 0 && j4 < 96 ){
628     //
629 mocchiut 1.6 // NOTICE: THE FOLLOWING LINE EXCLUDE ALL STRIPS FOR WHICH THE RMS*4 IS GREATER THAN 26 !!! <=============== IMPORTANT! =================>
630 mocchiut 1.1 //
631 mocchiut 1.6 if ( obadmask[l][m][j4] == 1 || clevel1->estrip[j4][m][l] <= clevel1->emin || calrms[l][m][j4] > 26 ){
632 mocchiut 1.1 clevel1->estrip[j4][m][l] = 0.;
633     };
634 mocchiut 1.6 //
635     // code and save the energy for each strip in svstrip
636     //
637 mocchiut 1.1 if ( clevel1->estrip[j4][m][l] > clevel1->emin ){
638 mocchiut 1.6 //
639     tim = 100000.;
640     plo = m;
641     fbi = 0;
642 mocchiut 1.8 if ( clevel1->estrip[j4][m][l] > 0.99995 ){
643 mocchiut 1.6 tim = 10000.;
644     plo = m;
645     fbi = 1;
646     };
647 mocchiut 1.8 if ( clevel1->estrip[j4][m][l] > 9.9995 ){
648 mocchiut 1.6 tim = 1000.;
649     plo = 22 + m;
650     fbi = 1;
651     };
652 mocchiut 1.8 if ( clevel1->estrip[j4][m][l] > 99.995 ){
653 mocchiut 1.6 tim = 100.;
654     plo = 22 + m;
655     fbi = 0;
656     };
657 mocchiut 1.8 if ( clevel1->estrip[j4][m][l] > 999.95 ){
658 mocchiut 1.6 tim = 10.;
659     plo = 44 + m;
660     fbi = 0;
661     };
662 mocchiut 1.8 if ( clevel1->estrip[j4][m][l] > 9999.5 ){
663 mocchiut 1.6 tim = 1.;
664     plo = 66 + m;
665     fbi = 0;
666     };
667     //
668     cle = (Int_t)lroundf(tim*clevel1->estrip[j4][m][l]);
669     //
670 mocchiut 1.1 if ( l == 0 ){
671     //
672     // +-PPSSmmmm.mmmm
673     //
674 mocchiut 1.6 svstrip[istrip] = fbi*1000000000 + plo*10000000 + j4*100000 + cle;
675 mocchiut 1.1 } else {
676 mocchiut 1.6 svstrip[istrip] = -(fbi*1000000000 + plo*10000000 + j4*100000 + cle);
677 mocchiut 1.1 };
678 mocchiut 1.8 //
679     // if ( ei >= -770 ) printf(" j %i l %i m %i estrip %f \n",j4,l,m,clevel1->estrip[j4][m][l]);
680     // if ( ei >= -770 ) printf(" num lim %i fbi %i tim %f plo %i cle %i \n",numeric_limits<Int_t>::max(),fbi,tim,plo,cle);
681     // if ( ei >= -770 ) printf(" svstrip %i \n",svstrip[istrip]);
682 mocchiut 1.6 //
683 mocchiut 1.1 istrip++;
684     };
685     };
686     };
687     //
688     } else {
689     for (Int_t nn = 0; nn < 96; nn++){
690     clevel1->estrip[nn][m][l] = 0.;
691     };
692     };
693     };
694     };
695     if ( !pe ){
696     clevel2->good = 1;
697     } else {
698     clevel2->good = 0;
699     };
700     return(0);
701     }
702    
703     void CaloProcessing::GetTrkVar(){
704     calol2tr();
705     }
706    
707     void CaloProcessing::FillTrkVar(CaloLevel2 *ca, Int_t nutrk){
708     //
709     CaloTrkVar *t_ca = new CaloTrkVar();
710     //
711     t_ca->trkseqno = trkseqno;
712     t_ca->ncore = (Int_t)clevel2->ncore;
713     t_ca->qcore = clevel2->qcore;
714     t_ca->noint = (Int_t)clevel2->noint;
715     t_ca->ncyl = (Int_t)clevel2->ncyl;
716     t_ca->qcyl = clevel2->qcyl;
717     t_ca->qtrack = clevel2->qtrack;
718     t_ca->qtrackx = clevel2->qtrackx;
719     t_ca->qtracky = clevel2->qtracky;
720     t_ca->dxtrack = clevel2->dxtrack;
721     t_ca->dytrack = clevel2->dytrack;
722     t_ca->qlast = clevel2->qlast;
723     t_ca->nlast = (Int_t)clevel2->nlast;
724     t_ca->qpre = clevel2->qpre;
725     t_ca->npre = (Int_t)clevel2->npre;
726     t_ca->qpresh = clevel2->qpresh;
727     t_ca->npresh = (Int_t)clevel2->npresh;
728     t_ca->qtr = clevel2->qtr;
729     t_ca->ntr = (Int_t)clevel2->ntr;
730     t_ca->planetot = (Int_t)clevel2->planetot;
731     t_ca->qmean = clevel2->qmean;
732     t_ca->dX0l = clevel2->dX0l;
733     t_ca->qlow = clevel2->qlow;
734     t_ca->nlow = (Int_t)clevel2->nlow;
735     //
736     if ( trkseqno == -1 ){
737 mocchiut 1.11 // ca->impx = clevel2->impx;
738     // ca->impy = clevel2->impy;
739     ca->tanx[1] = clevel2->tanx;
740     ca->tany[1] = clevel2->tany;
741 mocchiut 1.1 ca->elen = clevel2->elen;
742     ca->selen = clevel2->selen;
743 mocchiut 1.11 // memcpy(ca->cibar,clevel2->cibar,sizeof(clevel2->cibar));
744     // memcpy(ca->cbar,clevel2->cbar,sizeof(clevel2->cbar));
745     memcpy(t_ca->tibar,clevel2->cibar,sizeof(clevel2->cibar));
746     memcpy(t_ca->tbar,clevel2->cbar,sizeof(clevel2->cbar));
747 mocchiut 1.1 memcpy(ca->planemax,clevel2->planemax,sizeof(clevel2->planemax));
748 mocchiut 1.11 ca->varcfit[2] = clevel2->varcfit[0];
749     ca->varcfit[3] = clevel2->varcfit[1];
750     ca->npcfit[2] = clevel2->npcfit[0];
751     ca->npcfit[3] = clevel2->npcfit[1];
752     // memcpy(ca->varcfit,clevel2->varcfit,sizeof(clevel2->varcfit));
753     // memcpy(ca->npcfit,clevel2->npcfit,sizeof(clevel2->npcfit));
754     } else {
755     memcpy(t_ca->tibar,clevel2->tibar,sizeof(clevel2->tibar));
756     memcpy(t_ca->tbar,clevel2->tbar,sizeof(clevel2->tbar));
757 mocchiut 1.1 };
758 mocchiut 1.11 //
759 mocchiut 1.1 //
760 mocchiut 1.11 if ( !(ca->CaloTrk) ) ca->CaloTrk = new TClonesArray("CaloTrkVar",1); //ELENA
761 mocchiut 1.1 TClonesArray &t = *ca->CaloTrk;
762     new(t[nutrk]) CaloTrkVar(*t_ca);
763     //
764 mocchiut 1.2 delete t_ca;
765     //
766 mocchiut 1.1 ClearTrkVar();
767     }
768    
769     void CaloProcessing::GetCommonVar(){
770     calol2cm();
771     }
772    
773 mocchiut 1.8 void CaloProcessing::FillCommonVar(CaloLevel1 *c1, CaloLevel2 *ca){
774 mocchiut 1.1 //
775     ca->good = clevel2->good;
776     if ( clevel2->trigty == 2. ){
777     ca->selftrigger = 1;
778     } else {
779     ca->selftrigger = 0;
780     };
781 mocchiut 1.10 //
782     ca->selftrigger += (Int_t)clevel2->wartrig;
783     //
784 mocchiut 1.1 memcpy(ca->perr,clevel2->perr,sizeof(clevel2->perr));
785     memcpy(ca->swerr,clevel2->swerr,sizeof(clevel2->swerr));
786     memcpy(ca->crc,clevel2->crc,sizeof(clevel2->crc));
787     ca->nstrip = (Int_t)clevel2->nstrip;
788     ca->qtot = clevel2->qtot;
789 mocchiut 1.11 // ca->impx = clevel2->impx;
790     // ca->impy = clevel2->impy;
791     ca->tanx[0] = clevel2->tanx;
792     ca->tany[0] = clevel2->tany;
793 mocchiut 1.1 ca->nx22 = (Int_t)clevel2->nx22;
794     ca->qx22 = clevel2->qx22;
795     ca->qmax = clevel2->qmax;
796     ca->elen = clevel2->elen;
797     ca->selen = clevel2->selen;
798     memcpy(ca->qq,clevel2->qq,sizeof(clevel2->qq));
799     memcpy(ca->planemax,clevel2->planemax,sizeof(clevel2->planemax));
800 mocchiut 1.11 ca->varcfit[0] = clevel2->varcfit[0];
801     ca->varcfit[1] = clevel2->varcfit[1];
802     ca->npcfit[0] = clevel2->npcfit[0];
803     ca->npcfit[1] = clevel2->npcfit[1];
804     ca->fitmode[0] = clevel2->fmode[0];
805     ca->fitmode[1] = clevel2->fmode[1];
806     // memcpy(ca->varcfit,clevel2->varcfit,sizeof(clevel2->varcfit));
807     // memcpy(ca->npcfit,clevel2->npcfit,sizeof(clevel2->npcfit));
808 mocchiut 1.1 memcpy(ca->cibar,clevel2->cibar,sizeof(clevel2->cibar));
809     memcpy(ca->cbar,clevel2->cbar,sizeof(clevel2->cbar));
810     //
811 mocchiut 1.8 if ( c1 ){
812     c1->istrip = istrip;
813     c1->estrip = TArrayI(istrip,svstrip);
814     };
815     //
816 mocchiut 1.1 }
817    
818     void CaloProcessing::ClearStructs(){
819     ClearTrkVar();
820     ClearCommonVar();
821     }
822    
823     void CaloProcessing::RunClose(){
824     l0tr->Delete();
825     ClearStructs();
826 mocchiut 1.3 //
827     memset(dexy, 0, 2*22*96*sizeof(Float_t));
828     memset(dexyc, 0, 2*22*96*sizeof(Float_t));
829     memset(base, 0, 2*22*6*sizeof(Float_t));
830     memset(sbase, 0, 2*22*6*sizeof(Float_t));
831     //
832 mocchiut 1.1 }
833    
834     //
835     // Private methods
836     //
837    
838     void CaloProcessing::ClearTrkVar(){
839     clevel2->ncore = 0;
840     clevel2->qcore = 0.;
841     clevel2->noint = 0.;
842     clevel2->ncyl = 0.;
843     clevel2->qcyl = 0.;
844     clevel2->qtrack = 0.;
845     clevel2->qtrackx = 0.;
846     clevel2->qtracky = 0.;
847     clevel2->dxtrack = 0.;
848     clevel2->dytrack = 0.;
849     clevel2->qlast = 0.;
850     clevel2->nlast = 0.;
851     clevel2->qpre = 0.;
852     clevel2->npre = 0.;
853     clevel2->qpresh = 0.;
854     clevel2->npresh = 0.;
855     clevel2->qlow = 0.;
856     clevel2->nlow = 0.;
857     clevel2->qtr = 0.;
858     clevel2->ntr = 0.;
859     clevel2->planetot = 0.;
860     clevel2->qmean = 0.;
861     clevel2->dX0l = 0.;
862     clevel2->elen = 0.;
863     clevel2->selen = 0.;
864 mocchiut 1.3 memset(clevel1->al_p, 0, 5*2*sizeof(Double_t));
865     memset(clevel2->tibar, 0, 2*22*sizeof(Int_t));
866     memset(clevel2->tbar, 0, 2*22*sizeof(Float_t));
867 mocchiut 1.1 }
868    
869     void CaloProcessing::ClearCommonVar(){
870     istrip = 0;
871     clevel2->trigty = -1.;
872 mocchiut 1.10 clevel2->wartrig = 0.;
873 mocchiut 1.1 clevel2->good = 0;
874     clevel2->nstrip = 0.;
875     clevel2->qtot = 0.;
876 mocchiut 1.11 // clevel2->impx = 0.;
877     // clevel2->impy = 0.;
878     clevel2->tanx = 0.; // this is correct since it refers to the fortran structure
879     clevel2->tany = 0.; // this is correct since it refers to the fortran structure
880 mocchiut 1.1 clevel2->qmax = 0.;
881     clevel2->nx22 = 0.;
882     clevel2->qx22 = 0.;
883 mocchiut 1.3 memset(clevel2->perr, 0, 4*sizeof(Int_t));
884     memset(clevel2->swerr, 0, 4*sizeof(Int_t));
885     memset(clevel2->crc, 0, 4*sizeof(Int_t));
886     memset(clevel2->qq, 0, 4*sizeof(Int_t));
887 mocchiut 1.11 memset(clevel2->varcfit, 0, 4*sizeof(Float_t));
888     memset(clevel2->npcfit, 0, 4*sizeof(Int_t));
889 mocchiut 1.3 memset(clevel2->planemax, 0, 2*sizeof(Int_t));
890 mocchiut 1.11 memset(clevel2->fmode, 0, 2*sizeof(Int_t));
891 mocchiut 1.3 memset(clevel2->cibar, 0, 2*22*sizeof(Int_t));
892     memset(clevel2->cbar, 0, 2*22*sizeof(Float_t));
893 mocchiut 1.1 }
894    
895     void CaloProcessing::ClearCalibVals(Int_t s){
896     //
897     for ( Int_t d=0 ; d<11 ;d++ ){
898     Int_t pre = -1;
899     for ( Int_t j=0; j<96 ;j++){
900     if ( j%16 == 0 ) pre++;
901     if ( s == 2 ){
902     calped[0][2*d+1][j] = 0.;
903     cstwerr[3] = 0.;
904     cperror[3] = 0.;
905     calgood[0][2*d+1][j] = 0.;
906     calthr[0][2*d+1][pre] = 0.;
907     calrms[0][2*d+1][j] = 0.;
908     calbase[0][2*d+1][pre] = 0.;
909     calvar[0][2*d+1][pre] = 0.;
910     };
911     if ( s == 3 ){
912     calped[0][2*d][j] = 0.;
913     cstwerr[1] = 0.;
914     cperror[1] = 0.;
915     calgood[0][2*d][j] = 0.;
916     calthr[0][2*d][pre] = 0.;
917     calrms[0][2*d][j] = 0.;
918     calbase[0][2*d][pre] = 0.;
919     calvar[0][2*d][pre] = 0.;
920     };
921     if ( s == 0 ){
922     calped[1][2*d][j] = 0.;
923     cstwerr[0] = 0.;
924     cperror[0] = 0.;
925     calgood[1][2*d][j] = 0.;
926     calthr[1][2*d][pre] = 0.;
927     calrms[1][2*d][j] = 0.;
928     calbase[1][2*d][pre] = 0.;
929     calvar[1][2*d][pre] = 0.;
930     };
931     if ( s == 1 ){
932     calped[1][2*d+1][j] = 0.;
933     cstwerr[2] = 0.;
934     cperror[2] = 0.;
935     calgood[1][2*d+1][j] = 0.;
936     calthr[1][2*d+1][pre] = 0.;
937     calrms[1][2*d+1][j] = 0.;
938     calbase[1][2*d+1][pre] = 0.;
939     calvar[1][2*d+1][pre] = 0.;
940     };
941     };
942     };
943     return;
944     }
945    
946 mocchiut 1.4 Int_t CaloProcessing::Update(TSQLServer *dbc, UInt_t atime, Int_t s){
947 mocchiut 1.1 //
948     Int_t sgnl = 0;
949     //
950     GL_CALO_CALIB *glcalo = new GL_CALO_CALIB();
951     //
952     sgnl = 0;
953     //
954 mocchiut 1.4 idcalib[s] = 0;
955     fromtime[s] = 0;
956     totime[s] = 0;
957 mocchiut 1.1 calibno[s] = 0;
958     ClearCalibVals(s);
959     //
960 mocchiut 1.5 UInt_t uptime = 0;
961     //
962     sgnl = glcalo->Query_GL_CALO_CALIB(atime,uptime,s,dbc);
963 mocchiut 1.1 if ( sgnl < 0 ){
964     if ( verbose ) printf(" CALORIMETER - ERROR: error from GLTables\n");
965     return(sgnl);
966     };
967     //
968 mocchiut 1.4 idcalib[s] = glcalo->ID_ROOT_L0;
969 mocchiut 1.1 fromtime[s] = glcalo->FROM_TIME;
970 mocchiut 1.5 if ( glcalo->TO_TIME < atime ){ // calibration is corrupted and we are using the one that preceed the good one
971     totime[s] = uptime;
972     } else {
973     totime[s] = glcalo->TO_TIME;
974     };
975     // totime[s] = glcalo->TO_TIME;
976 mocchiut 1.4 calibno[s] = glcalo->EV_ROOT;
977 mocchiut 1.1 //
978 mocchiut 1.4 if ( totime[s] == 0 ){
979 mocchiut 1.1 if ( verbose ) printf(" CALORIMETER - WARNING: data with no associated calibration\n");
980     ClearCalibVals(s);
981     sgnl = 100;
982     };
983     //
984     // determine path and name and entry of the calibration file
985     //
986     GL_ROOT *glroot = new GL_ROOT();
987     if ( verbose ) printf("\n");
988     if ( verbose ) printf(" ** SECTION %i **\n",s);
989     //
990     sgnl = glroot->Query_GL_ROOT(idcalib[s],dbc);
991     if ( sgnl < 0 ){
992     if ( verbose ) printf(" CALORIMETER - ERROR: error from GLTables\n");
993     return(sgnl);
994     };
995     //
996     stringstream name;
997     name.str("");
998     name << glroot->PATH.Data() << "/";
999     name << glroot->NAME.Data();
1000     //
1001     fcalname[s] = (TString)name.str().c_str();
1002 mocchiut 1.4 if ( verbose ) printf(" - event at time %u. From time %u to time %u \n use file %s \n calibration at entry %i \n\n",atime,fromtime[s],totime[s],fcalname[s].Data(),calibno[s]);
1003 mocchiut 1.1 //
1004     sgnl = LoadCalib(s);
1005     //
1006     if ( sgnl != 0 ) return(sgnl);
1007 mocchiut 1.2 delete glcalo;
1008     delete glroot;
1009 mocchiut 1.1 //
1010     return(0);
1011     //
1012     }
1013    
1014     Int_t CaloProcessing::LoadCalib(Int_t s){
1015     //
1016     ifstream myfile;
1017     myfile.open(fcalname[s].Data());
1018     if ( !myfile ){
1019     return(-107);
1020     };
1021     myfile.close();
1022     //
1023     TFile *File = new TFile(fcalname[s].Data());
1024     if ( !File ) return(-108);
1025     TTree *tr = (TTree*)File->Get("CalibCalPed");
1026     if ( !tr ) return(-109);
1027     //
1028     TBranch *calo = tr->GetBranch("CalibCalPed");
1029     //
1030     pamela::CalibCalPedEvent *ce = 0;
1031     tr->SetBranchAddress("CalibCalPed", &ce);
1032     //
1033 mocchiut 1.4 Long64_t ncalibs = calo->GetEntries();
1034 mocchiut 1.1 //
1035     if ( !ncalibs ) return(-110);
1036     //
1037 mocchiut 1.4 calo->GetEntry(calibno[s]);
1038 mocchiut 1.1 //
1039     if (ce->cstwerr[s] != 0 && ce->cperror[s] == 0 ) {
1040     for ( Int_t d=0 ; d<11 ;d++ ){
1041     Int_t pre = -1;
1042     for ( Int_t j=0; j<96 ;j++){
1043     if ( j%16 == 0 ) pre++;
1044     if ( s == 2 ){
1045     calped[0][2*d+1][j] = ce->calped[3][d][j];
1046     cstwerr[3] = ce->cstwerr[3];
1047     cperror[3] = ce->cperror[3];
1048     calgood[0][2*d+1][j] = ce->calgood[3][d][j];
1049     calthr[0][2*d+1][pre] = ce->calthr[3][d][pre];
1050     calrms[0][2*d+1][j] = ce->calrms[3][d][j];
1051     calbase[0][2*d+1][pre] = ce->calbase[3][d][pre];
1052     calvar[0][2*d+1][pre] = ce->calvar[3][d][pre];
1053     };
1054     if ( s == 3 ){
1055     calped[0][2*d][j] = ce->calped[1][d][j];
1056     cstwerr[1] = ce->cstwerr[1];
1057     cperror[1] = ce->cperror[1];
1058     calgood[0][2*d][j] = ce->calgood[1][d][j];
1059     calthr[0][2*d][pre] = ce->calthr[1][d][pre];
1060     calrms[0][2*d][j] = ce->calrms[1][d][j];
1061     calbase[0][2*d][pre] = ce->calbase[1][d][pre];
1062     calvar[0][2*d][pre] = ce->calvar[1][d][pre];
1063     };
1064     if ( s == 0 ){
1065     calped[1][2*d][j] = ce->calped[0][d][j];
1066     cstwerr[0] = ce->cstwerr[0];
1067     cperror[0] = ce->cperror[0];
1068     calgood[1][2*d][j] = ce->calgood[0][d][j];
1069     calthr[1][2*d][pre] = ce->calthr[0][d][pre];
1070     calrms[1][2*d][j] = ce->calrms[0][d][j];
1071     calbase[1][2*d][pre] = ce->calbase[0][d][pre];
1072     calvar[1][2*d][pre] = ce->calvar[0][d][pre];
1073     };
1074     if ( s == 1 ){
1075     calped[1][2*d+1][j] = ce->calped[2][d][j];
1076     cstwerr[2] = ce->cstwerr[2];
1077     cperror[2] = ce->cperror[2];
1078     calgood[1][2*d+1][j] = ce->calgood[2][d][j];
1079     calthr[1][2*d+1][pre] = ce->calthr[2][d][pre];
1080     calrms[1][2*d+1][j] = ce->calrms[2][d][j];
1081     calbase[1][2*d+1][pre] = ce->calbase[2][d][pre];
1082     calvar[1][2*d+1][pre] = ce->calvar[2][d][pre];
1083     };
1084     };
1085     };
1086     } else {
1087     if ( verbose ) printf(" CALORIMETER - ERROR: problems finding a good calibration in this file! \n\n ");
1088     return(-111);
1089     };
1090     File->Close();
1091     return(0);
1092     }

  ViewVC Help
Powered by ViewVC 1.1.23