/[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.8 - (show annotations) (download)
Wed Oct 10 16:01:30 2007 UTC (17 years, 1 month ago) by mocchiut
Branch: MAIN
Changes since 1.7: +12 -0 lines
Bug in ToFPMT fixed, set zerofill option as defatult, DV exit code changed, set time 0 on each DB connection

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

  ViewVC Help
Powered by ViewVC 1.1.23