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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.13 - (show 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 /**
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 delete de;
49 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 //
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 calopar1 = true;
68 calopar2 = true;
69 calopar3 = true;
70 crosst = true;
71 ftcalopar1 = 0;
72 ttcalopar1 = 0;
73 ftcalopar2 = 0;
74 ttcalopar2 = 0;
75 ftcalopar3 = 0;
76 ttcalopar3 = 0;
77 }
78
79 void CaloProcessing::SetCrossTalk(Bool_t ct){
80 crosst = ct;
81 };
82
83 /**
84 * Initialize CaloProcessing object
85 **/
86 void CaloProcessing::ProcessingInit(TSQLServer *dbc, UInt_t hs, Int_t &sgnl, TTree *l0tree, Bool_t isdeb, Bool_t isverb){
87 //
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 UInt_t uptime = 0;
103 //
104 for (Int_t s = 0; s < 4; s++){
105 idcalib[s] = 0;
106 fromtime[s] = 0;
107 totime[s] = 0;
108 calibno[s] = 0;
109 ClearCalibVals(s);
110 //
111 sgnl = glcalo->Query_GL_CALO_CALIB(hs,uptime,s,dbc);
112 if ( sgnl < 0 ){
113 if ( verbose ) printf(" CALORIMETER - ERROR: error from GLTables\n");
114 return;
115 };
116 //
117 idcalib[s] = glcalo->ID_ROOT_L0;
118 fromtime[s] = glcalo->FROM_TIME;
119 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 calibno[s] = glcalo->EV_ROOT;
125 //
126 if ( totime[s] == 0 ){
127 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 if ( totime[s] > 0 ){
140 //
141 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 //
147 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 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 } else {
155 if ( verbose ) printf(" - runheader at time %u. NO CALIBRATION INCLUDE THE RUNHEADER! ",hs);
156 };
157 sgnl = LoadCalib(s);
158 if ( sgnl ) break;
159 };
160 //
161 delete glcalo;
162 delete glroot;
163 //
164 return;
165 //
166 }
167
168 Int_t CaloProcessing::ChkCalib(TSQLServer *dbc, UInt_t atime){
169 Int_t sgnl = 0;
170 for ( Int_t s = 0; s < 4; s++){
171 if ( atime > totime[s] ){
172 if ( !dbc->IsConnected() ) throw -116;
173 sgnl = Update(dbc,atime,s);
174 if ( sgnl < 0 ) return(sgnl);
175 };
176 };
177 return(sgnl);
178 }
179
180 Int_t CaloProcessing::ChkParam(TSQLServer *dbc, UInt_t runheader){
181 stringstream calfile;
182 stringstream bmfile;
183 stringstream aligfile;
184 Int_t error = 0;
185 FILE *f = 0;
186 ifstream badfile;
187 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 error = glparam->Query_GL_PARAM(runheader,101,dbc);
198 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 error = glparam->Query_GL_PARAM(runheader,102,dbc);
231 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 if ( verbose ) printf("\n Using alignment file: \n %s \n",aligfile.str().c_str());
240 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 // 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 delete glparam;
308 //
309 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 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 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 Bool_t ck = false;
372 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 clevel2->wartrig = 0.;
377 clevel2->trigty = 3.;
378 ck = true;
379 break;
380 };
381 };
382 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 };
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 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 // 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 cbase0 = 0.;
483 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 cbase0 += base[l][m][pre];
492 } 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 cbase0 += base[l][m][pre];
499 };
500 };
501 };
502 //
503 // run over strips
504 //
505 pre = -1;
506 ener0 = 0.;
507 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 qpre[pre] = 0.;
517 };
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 ener0 += ener;
568 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 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 };
594 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 };
600 };
601 };
602 //
603 if ( ener0 == 0. && cbase0 == 0. && !pproblem && clevel2->perr[se] == 0){
604 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 pproblem = true;
606 pe++;
607 };
608 //
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 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 };
627 if ( j4 >= 0 && j4 < 96 ){
628 //
629 // NOTICE: THE FOLLOWING LINE EXCLUDE ALL STRIPS FOR WHICH THE RMS*4 IS GREATER THAN 26 !!! <=============== IMPORTANT! =================>
630 //
631 if ( obadmask[l][m][j4] == 1 || clevel1->estrip[j4][m][l] <= clevel1->emin || calrms[l][m][j4] > 26 ){
632 clevel1->estrip[j4][m][l] = 0.;
633 };
634 //
635 // code and save the energy for each strip in svstrip
636 //
637 if ( clevel1->estrip[j4][m][l] > clevel1->emin ){
638 //
639 tim = 100000.;
640 plo = m;
641 fbi = 0;
642 if ( clevel1->estrip[j4][m][l] > 0.99995 ){
643 tim = 10000.;
644 plo = m;
645 fbi = 1;
646 };
647 if ( clevel1->estrip[j4][m][l] > 9.9995 ){
648 tim = 1000.;
649 plo = 22 + m;
650 fbi = 1;
651 };
652 if ( clevel1->estrip[j4][m][l] > 99.995 ){
653 tim = 100.;
654 plo = 22 + m;
655 fbi = 0;
656 };
657 if ( clevel1->estrip[j4][m][l] > 999.95 ){
658 tim = 10.;
659 plo = 44 + m;
660 fbi = 0;
661 };
662 if ( clevel1->estrip[j4][m][l] > 9999.5 ){
663 tim = 1.;
664 plo = 66 + m;
665 fbi = 0;
666 };
667 //
668 cle = (Int_t)lroundf(tim*clevel1->estrip[j4][m][l]);
669 //
670 if ( l == 0 ){
671 //
672 // +-PPSSmmmm.mmmm
673 //
674 svstrip[istrip] = fbi*1000000000 + plo*10000000 + j4*100000 + cle;
675 } else {
676 svstrip[istrip] = -(fbi*1000000000 + plo*10000000 + j4*100000 + cle);
677 };
678 //
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 //
683 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 // ca->impx = clevel2->impx;
738 // ca->impy = clevel2->impy;
739 ca->tanx[1] = clevel2->tanx;
740 ca->tany[1] = clevel2->tany;
741 ca->elen = clevel2->elen;
742 ca->selen = clevel2->selen;
743 // 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 memcpy(ca->planemax,clevel2->planemax,sizeof(clevel2->planemax));
748 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 };
758 //
759 //
760 if ( !(ca->CaloTrk) ) ca->CaloTrk = new TClonesArray("CaloTrkVar",1); //ELENA
761 TClonesArray &t = *ca->CaloTrk;
762 new(t[nutrk]) CaloTrkVar(*t_ca);
763 //
764 delete t_ca;
765 //
766 ClearTrkVar();
767 }
768
769 void CaloProcessing::GetCommonVar(){
770 calol2cm();
771 }
772
773 void CaloProcessing::FillCommonVar(CaloLevel1 *c1, CaloLevel2 *ca){
774 //
775 ca->good = clevel2->good;
776 if ( clevel2->trigty == 2. ){
777 ca->selftrigger = 1;
778 } else {
779 ca->selftrigger = 0;
780 };
781 //
782 ca->selftrigger += (Int_t)clevel2->wartrig;
783 //
784 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 // ca->impx = clevel2->impx;
790 // ca->impy = clevel2->impy;
791 ca->tanx[0] = clevel2->tanx;
792 ca->tany[0] = clevel2->tany;
793 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 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 memcpy(ca->cibar,clevel2->cibar,sizeof(clevel2->cibar));
809 memcpy(ca->cbar,clevel2->cbar,sizeof(clevel2->cbar));
810 //
811 if ( c1 ){
812 c1->istrip = istrip;
813 c1->estrip = TArrayI(istrip,svstrip);
814 };
815 //
816 }
817
818 void CaloProcessing::ClearStructs(){
819 ClearTrkVar();
820 ClearCommonVar();
821 }
822
823 void CaloProcessing::RunClose(){
824 l0tr->Delete();
825 ClearStructs();
826 //
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 }
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 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 }
868
869 void CaloProcessing::ClearCommonVar(){
870 istrip = 0;
871 clevel2->trigty = -1.;
872 clevel2->wartrig = 0.;
873 clevel2->good = 0;
874 clevel2->nstrip = 0.;
875 clevel2->qtot = 0.;
876 // 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 clevel2->qmax = 0.;
881 clevel2->nx22 = 0.;
882 clevel2->qx22 = 0.;
883 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 memset(clevel2->varcfit, 0, 4*sizeof(Float_t));
888 memset(clevel2->npcfit, 0, 4*sizeof(Int_t));
889 memset(clevel2->planemax, 0, 2*sizeof(Int_t));
890 memset(clevel2->fmode, 0, 2*sizeof(Int_t));
891 memset(clevel2->cibar, 0, 2*22*sizeof(Int_t));
892 memset(clevel2->cbar, 0, 2*22*sizeof(Float_t));
893 }
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 Int_t CaloProcessing::Update(TSQLServer *dbc, UInt_t atime, Int_t s){
947 //
948 Int_t sgnl = 0;
949 //
950 GL_CALO_CALIB *glcalo = new GL_CALO_CALIB();
951 //
952 sgnl = 0;
953 //
954 idcalib[s] = 0;
955 fromtime[s] = 0;
956 totime[s] = 0;
957 calibno[s] = 0;
958 ClearCalibVals(s);
959 //
960 UInt_t uptime = 0;
961 //
962 sgnl = glcalo->Query_GL_CALO_CALIB(atime,uptime,s,dbc);
963 if ( sgnl < 0 ){
964 if ( verbose ) printf(" CALORIMETER - ERROR: error from GLTables\n");
965 return(sgnl);
966 };
967 //
968 idcalib[s] = glcalo->ID_ROOT_L0;
969 fromtime[s] = glcalo->FROM_TIME;
970 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 calibno[s] = glcalo->EV_ROOT;
977 //
978 if ( totime[s] == 0 ){
979 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 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 //
1004 sgnl = LoadCalib(s);
1005 //
1006 if ( sgnl != 0 ) return(sgnl);
1007 delete glcalo;
1008 delete glroot;
1009 //
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 Long64_t ncalibs = calo->GetEntries();
1034 //
1035 if ( !ncalibs ) return(-110);
1036 //
1037 calo->GetEntry(calibno[s]);
1038 //
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