/[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.12 - (show annotations) (download)
Tue Jan 23 11:52:25 2007 UTC (17 years, 10 months ago) by mocchiut
Branch: MAIN
CVS Tags: v3r00
Changes since 1.11: +1 -1 lines
OrbitalInfo updated, small changes in calo code

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

  ViewVC Help
Powered by ViewVC 1.1.23