/[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.6 - (show annotations) (download)
Tue Nov 14 14:08:50 2006 UTC (18 years, 1 month ago) by mocchiut
Branch: MAIN
Changes since 1.5: +121 -9 lines
Major calorimeter release, some news in DarthVader main

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

  ViewVC Help
Powered by ViewVC 1.1.23