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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.6 - (show annotations) (download)
Mon Sep 3 08:42:13 2007 UTC (17 years, 3 months ago) by mocchiut
Branch: MAIN
Changes since 1.5: +12 -0 lines
Calorimeter: added selfdelay variable

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

  ViewVC Help
Powered by ViewVC 1.1.23