/[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.37 - (show annotations) (download)
Tue Oct 14 14:07:10 2014 UTC (10 years, 4 months ago) by mocchiut
Branch: MAIN
CVS Tags: v10REDr01, v10RED, HEAD
Changes since 1.36: +15 -10 lines
10RED: lost sync bug fixed

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 // extern struct FlEventi eventi_;
60 // extern struct FlGruppo gruppo_;
61 // extern struct FlGruppo2 gruppo2_;
62 // extern struct FlGruppo4 gruppo4_;
63 // extern struct FlTaglioen taglioen_;
64 // extern struct FlAngolo angolo_;
65 // extern struct FlWhere where_;
66 // extern struct FlGeneral general_;
67 // extern struct FlCh ch_;
68 // extern struct FlCalofit calofit_;
69 // extern struct FlPawcd pawcd_;
70 // extern struct FlQuestd questd_;
71 // eventi = &eventi_;
72 // gruppo = &gruppo_;
73 // gruppo2 = &gruppo2_;
74 // gruppo4 = &gruppo4_;
75 // taglioen = &taglioen_;
76 // angolo = &angolo_;
77 // where = &where_;
78 // general = &general_;
79 // ch = &ch_;
80 // calofit = &calofit_;
81 // pawcd = &pawcd_;
82 // questd = &questd_;
83 //
84 trkseqno = 0;
85 ClearStructs();
86 //
87 memset(dexy, 0, 2*22*96*sizeof(Float_t));
88 memset(dexyc, 0, 2*22*96*sizeof(Float_t));
89 memset(mip, 0, 2*22*96*sizeof(Float_t));
90 memset(base, 0, 2*22*6*sizeof(Float_t));
91 memset(sbase, 0, 2*22*6*sizeof(Float_t));
92 memset(obadmask, 0, 2*22*96*sizeof(Int_t));
93 memset(obadpulsemask, 0, 2*22*6*sizeof(Int_t));
94 memset(ctprecor, 0, 2*22*6*sizeof(Float_t));
95 memset(ctsicor, 0, 2*22*9*sizeof(Float_t));
96 memset(ctneigcor, 0, 2*22*6*sizeof(Float_t));
97 calopar1 = true;
98 calopar2 = true;
99 calopar3 = true;
100 calopar4 = true;
101 calopar5 = true;
102 crosst = true;
103 mask18 = false;
104 ftcalopar1 = 0;
105 ttcalopar1 = 0;
106 ftcalopar2 = 0;
107 ttcalopar2 = 0;
108 ftcalopar3 = 0;
109 ttcalopar3 = 0;
110 ftcalopar4 = 0;
111 ttcalopar4 = 0;
112 ftcalopar5 = 0;
113 ttcalopar5 = 0;
114 }
115
116 void CaloLevel0::SetCrossTalk(Bool_t ct){
117 crosst = ct;
118 }
119
120 void CaloLevel0::SetCrossTalkType(Bool_t ct){
121 ctground = ct;
122 }
123
124 void CaloLevel0::SetCrossTalkType(Int_t ct){
125 if ( ct == 0 ) ctground = true;
126 if ( ct == 1 ){
127 ctground = false;
128 noselfct = false;
129 };
130 if ( ct == 2 ){
131 ctground = false;
132 noselfct = true;
133 };
134 }
135
136 void CaloLevel0::SetVerbose(Bool_t ct){
137 verbose = ct;
138 }
139
140 /**
141 * Initialize CaloLevel0 object
142 **/
143 void CaloLevel0::ProcessingInit(TSQLServer *dbc, UInt_t hs, Int_t &sgnl, TTree *l0tree, Bool_t isdeb, Bool_t isverb){
144 if ( !dbc->IsConnected() ) throw -116;
145 this->InitDo(dbc,hs,sgnl,l0tree,isdeb,isverb);
146 }
147
148 /**
149 * Initialize CaloLevel0 object
150 **/
151 void CaloLevel0::ProcessingInit(GL_TABLES *glt, UInt_t hs, Int_t &sgnl, TTree *l0tree, Bool_t isdeb, Bool_t isverb){
152 //
153 const TString host = glt->CGetHost();
154 const TString user = glt->CGetUser();
155 const TString psw = glt->CGetPsw();
156 TSQLServer *dbc = TSQLServer::Connect(host.Data(),user.Data(),psw.Data());
157 if ( !dbc->IsConnected() ) throw -116;
158 this->InitDo(dbc,hs,sgnl,l0tree,isdeb,isverb);
159 dbc->Close();
160 delete dbc;
161 dbc = 0;
162 }
163
164
165 void CaloLevel0::InitDo(TSQLServer *dbc, UInt_t hs, Int_t &sgnl, TTree *l0tree, Bool_t isdeb, Bool_t isverb){
166 stringstream myquery;
167 myquery.str("");
168 myquery << "SET time_zone='+0:00';";
169 delete dbc->Query(myquery.str().c_str());
170 delete dbc->Query("SET sql_mode = 'NO_UNSIGNED_SUBTRACTION';");
171 //
172 debug = isdeb;
173 verbose = isverb;
174 //
175 l0tr=(TTree*)l0tree;
176 de = new pamela::calorimeter::CalorimeterEvent();
177 l0calo = (TBranch*)l0tr->GetBranch("Calorimeter");
178 l0tr->SetBranchAddress("Calorimeter", &de);
179 //
180 trkseqno = 0;
181 ClearStructs();
182 //
183 GL_CALO_CALIB *glcalo = new GL_CALO_CALIB();
184 //
185 sgnl = 0;
186 UInt_t uptime = 0;
187 //
188 for (Int_t s = 0; s < 4; s++){
189 idcalib[s] = 0;
190 fromtime[s] = 0;
191 totime[s] = 0;
192 calibno[s] = 0;
193 ClearCalibVals(s);
194 //
195 sgnl = glcalo->Query_GL_CALO_CALIB(hs,uptime,s,dbc);
196 if ( sgnl < 0 ){
197 if ( verbose ) printf(" CALORIMETER - ERROR: error from GLTables\n");
198 return;
199 };
200 //
201 idcalib[s] = glcalo->ID_ROOT_L0;
202 fromtime[s] = glcalo->FROM_TIME;
203 if ( glcalo->TO_TIME < hs ){ // calibration is corrupted and we are using the one that preceed the good one
204 totime[s] = uptime;
205 } else {
206 totime[s] = glcalo->TO_TIME;
207 };
208 calibno[s] = glcalo->EV_ROOT;
209 //
210 if ( totime[s] == 0 ){
211 if ( verbose ) printf(" CALORIMETER - WARNING: data with no associated calibration\n");
212 ClearCalibVals(s);
213 sgnl = 100;
214 };
215 };
216 //
217 // determine path and name and entry of the calibration file
218 //
219 GL_ROOT *glroot = new GL_ROOT();
220 if ( verbose ) printf("\n");
221 for (Int_t s = 0; s < 4; s++){
222 if ( verbose ) printf(" ** SECTION %i **\n",s);
223 if ( totime[s] > 0 ){
224 //
225 sgnl = glroot->Query_GL_ROOT(idcalib[s],dbc);
226 if ( sgnl < 0 ){
227 if ( verbose ) printf(" CALORIMETER - ERROR: error from GLTables\n");
228 return;
229 };
230 //
231 stringstream name;
232 name.str("");
233 name << glroot->PATH.Data() << "/";
234 name << glroot->NAME.Data();
235 //
236 fcalname[s] = (TString)name.str().c_str();
237 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]);
238 } else {
239 if ( verbose ) printf(" - runheader at time %u. NO CALIBRATION INCLUDE THE RUNHEADER! ",hs);
240 };
241 sgnl = LoadCalib(s);
242 if ( sgnl ) break;
243 };
244 //
245 delete glcalo;
246 delete glroot;
247 //
248 return;
249 //
250 }
251
252 Int_t CaloLevel0::ChkCalib(GL_TABLES *glt, UInt_t atime){
253 Int_t sgnl = 0;
254 for ( Int_t s = 0; s < 4; s++){
255 if ( atime > totime[s] ){
256 sgnl = Update(glt,atime,s);
257 if ( sgnl < 0 ) return(sgnl);
258 };
259 };
260 return(sgnl);
261 }
262
263 Int_t CaloLevel0::ChkParam(TSQLServer *dbc, UInt_t runheader, Bool_t mechal){
264 Int_t sig = this->ChkParamDo(dbc,runheader,mechal);
265 return(sig);
266 }
267
268 Int_t CaloLevel0::ChkParam(GL_TABLES *glt, UInt_t runheader, Bool_t mechal){
269 const TString host = glt->CGetHost();
270 const TString user = glt->CGetUser();
271 const TString psw = glt->CGetPsw();
272 TSQLServer *dbc = TSQLServer::Connect(host.Data(),user.Data(),psw.Data());
273 if ( !dbc->IsConnected() ) throw -116;
274 stringstream myquery;
275 myquery.str("");
276 myquery << "SET time_zone='+0:00';";
277 delete dbc->Query(myquery.str().c_str());
278 delete dbc->Query("SET sql_mode = 'NO_UNSIGNED_SUBTRACTION';");
279 //
280 Int_t sig = this->ChkParamDo(dbc,runheader,mechal);
281 dbc->Close();
282 delete dbc;
283 dbc = 0;
284 return(sig);
285 }
286
287 Int_t CaloLevel0::ChkParamDo(TSQLServer *dbc, UInt_t runheader, Bool_t mechal){
288 //
289 stringstream calfile;
290 stringstream bmfile;
291 stringstream aligfile;
292 Int_t error = 0;
293 FILE *f = 0;
294 ifstream badfile;
295 GL_PARAM *glparam = new GL_PARAM();
296 //
297 if ( calopar1 || ( ttcalopar1 != 0 && ttcalopar1 < runheader ) ){
298 //
299 if ( debug ) printf(" calopar1 %i ftcalopar1 %u ttcalopar1 %u runheader %u \n",calopar1,ftcalopar1,ttcalopar1,runheader);
300 //
301 if ( calopar1 ){
302 //
303 // determine where I can find calorimeter ADC to MIP conversion file
304 //
305 if ( verbose ) printf(" Querying DB for calorimeter parameters files...\n");
306 //
307 error = 0;
308 error = glparam->Query_GL_PARAM(runheader,101,dbc);
309 if ( error < 0 ) return(error);
310 //
311 calfile.str("");
312 calfile << glparam->PATH.Data() << "/";
313 calfile << glparam->NAME.Data();
314 //
315 if ( verbose ) printf("\n Using ADC to MIP conversion file: \n %s \n",calfile.str().c_str());
316 f = fopen(calfile.str().c_str(),"rb");
317 if ( !f ){
318 if ( verbose ) printf(" CALORIMETER - ERROR: no ADC to MIP file!\n");
319 return(-105);
320 };
321 //
322 for (Int_t m = 0; m < 2 ; m++ ){
323 for (Int_t k = 0; k < 22; k++ ){
324 for (Int_t l = 0; l < 96; l++ ){
325 fread(&mip[m][k][l],sizeof(mip[m][k][l]),1,f);
326 if ( debug ) printf(" %f \n",mip[m][k][l]);
327 };
328 };
329 };
330 fclose(f);
331 };
332 //
333 calopar1 = false;
334 //
335 // flight extra corrections:
336 //
337 if ( verbose ) printf(" Querying DB for calorimeter flight ADC to MIP files...\n");
338 //
339 error = 0;
340 error = glparam->Query_GL_PARAM(runheader,110,dbc);
341 if ( error < 0 ) return(error);
342 //
343 calfile.str("");
344 calfile << glparam->PATH.Data() << "/";
345 calfile << glparam->NAME.Data();
346 ftcalopar1 = glparam->FROM_TIME;
347 ttcalopar1 = glparam->TO_TIME;
348 //
349 if ( verbose ) printf("\n Using ADC to MIP special conversion file: \n %s \n",calfile.str().c_str());
350 ifstream spfile;
351 spfile.open(calfile.str().c_str());
352 if ( !spfile ){
353 if ( verbose ) printf(" CALORIMETER - ERROR: no special calibration file!\n");
354 return(-123);
355 };
356 //
357 Int_t vview = 0;
358 Int_t vplane = 0;
359 Int_t vstrip = 0;
360 Float_t vval = 0.;
361 while ( spfile >> vview && spfile >> vplane && spfile >> vstrip && spfile >> vval){
362 if ( debug ) printf(" Setting ADC to MIP conversion factor: view %i plane %i strip %i mip %f \n",vview,vplane,vstrip,vval);
363 mip[vview][vplane][vstrip] = vval;
364 };
365 //
366 };
367 //
368 //
369 if ( calopar2 || ( ttcalopar2 != 0 && ttcalopar2 < runheader ) ){
370 //
371 if ( debug ) printf(" calopar2 %i ftcalopar2 %u ttcalopar2 %u runheader %u \n",calopar2,ftcalopar2,ttcalopar2,runheader);
372 calopar2 = false;
373 //
374 // determine where I can find calorimeter alignment file
375 //
376 //
377 error = 0;
378 error = glparam->Query_GL_PARAM(runheader,102,dbc);
379 if ( error < 0 ) return(error);
380 //
381 aligfile.str("");
382 aligfile << glparam->PATH.Data() << "/";
383 aligfile << glparam->NAME.Data();
384 ftcalopar2 = glparam->FROM_TIME;
385 ttcalopar2 = glparam->TO_TIME;
386 //
387 if ( verbose ) printf("\n Using parameter file: \n %s \n",aligfile.str().c_str());
388 f = fopen(aligfile.str().c_str(),"rb");
389 if ( !f ){
390 if ( verbose ) printf(" CALORIMETER - ERROR: no parameter file!\n");
391 return(-106);
392 };
393 //
394 if ( !mechal ){
395 //
396 fread(&clevel1->xalig,sizeof(clevel1->xalig),1,f);
397 if ( debug ) printf(" xalig = %f \n",clevel1->xalig);
398 fread(&clevel1->yalig,sizeof(clevel1->yalig),1,f);
399 if ( debug ) printf(" yalig = %f \n",clevel1->yalig);
400 fread(&clevel1->zalig,sizeof(clevel1->zalig),1,f);
401 if ( debug ) printf(" zalig = %f \n",clevel1->zalig);
402 } else {
403 if ( verbose ) printf("\n Using MECHANICAL alignement parameters \n");
404 //
405 CaloStrip cs = CaloStrip();
406 cs.UseMechanicalAlig();
407 clevel1->xalig = cs.GetXalig();
408 if ( debug ) printf(" xalig = %f \n",clevel1->xalig);
409 clevel1->yalig = cs.GetYalig();
410 if ( debug ) printf(" yalig = %f \n",clevel1->yalig);
411 clevel1->zalig = cs.GetZalig();
412 if ( debug ) printf(" zalig = %f \n",clevel1->zalig);
413 //
414 Float_t tmp = 0;
415 fread(&tmp,sizeof(clevel1->xalig),1,f);
416 fread(&tmp,sizeof(clevel1->yalig),1,f);
417 fread(&tmp,sizeof(clevel1->zalig),1,f);
418 // clevel1->zalig = -265.82;
419 //
420 };
421 fread(&clevel1->emin,sizeof(clevel1->emin),1,f);
422 if ( debug ) printf(" signal threshold = %f \n",clevel1->emin);
423 //
424 fclose(f);
425 };
426 //
427 // Load offline bad strip mask
428 //
429 if ( calopar3 || ( ttcalopar3 != 0 && ttcalopar3 < runheader ) ){
430 if ( debug ) printf(" calopar3 %i ftcalopar3 %u ttcalopar3 %u runheader %u \n",calopar3,ftcalopar3,ttcalopar3,runheader);
431 calopar3 = false;
432 //
433 // determine where I can find calorimeter alignment file
434 //
435 //
436 error = 0;
437 error = glparam->Query_GL_PARAM(runheader,103,dbc);
438 if ( error < 0 ) return(error);
439 //
440 bmfile.str("");
441 bmfile << glparam->PATH.Data() << "/";
442 bmfile << glparam->NAME.Data();
443 ftcalopar3 = glparam->FROM_TIME;
444 ttcalopar3 = glparam->TO_TIME;
445 //
446 if ( verbose ) printf("\n Using bad strip offline mask file: \n %s \n\n",bmfile.str().c_str());
447 badfile.open(bmfile.str().c_str());
448 if ( !badfile ){
449 if ( verbose ) printf(" CALORIMETER - ERROR: no bad strip offline mask file!\n");
450 return(-115);
451 };
452 //
453 Bool_t isdone = false;
454 Int_t bad = 0;
455 Int_t view = 1;
456 Int_t strip = 0;
457 Int_t plane = 21;
458 while ( !isdone ) {
459 badfile >> bad;
460 obadmask[view][plane][strip] = bad;
461 if ( debug && bad ) printf(" SETTING view %i plane %i strip %i BAD = %i \n",view,plane,strip,bad);
462 strip++;
463 if ( strip > 95 ){
464 strip = 0;
465 plane--;
466 if ( plane < 0 ){
467 plane = 21;
468 view--;
469 };
470 if ( view < 0 ) isdone = true;
471 };
472 };
473 //
474 badfile.close();
475 };
476 //
477 // calopar4
478 //
479 if ( calopar4 || ( ttcalopar4 != 0 && ttcalopar4 < runheader ) ){
480 //
481 if ( debug ) printf(" calopar4 %i ftcalopar4 %u ttcalopar4 %u runheader %u \n",calopar4,ftcalopar4,ttcalopar4,runheader);
482 //
483 calopar4 = false;
484 //
485 // flight extra corrections:
486 //
487 if ( verbose ) printf(" Querying DB for calorimeter max rms file...\n");
488 //
489 error = 0;
490 error = glparam->Query_GL_PARAM(runheader,109,dbc);
491 if ( error < 0 ) return(error);
492 //
493 calfile.str("");
494 calfile << glparam->PATH.Data() << "/";
495 calfile << glparam->NAME.Data();
496 ftcalopar4 = glparam->FROM_TIME;
497 ttcalopar4 = glparam->TO_TIME;
498 //
499 if ( verbose ) printf("\n Using calorimeter max rms file: \n %s \n",calfile.str().c_str());
500 ifstream spfile;
501 spfile.open(calfile.str().c_str());
502 if ( !spfile ){
503 if ( verbose ) printf(" CALORIMETER - ERROR: no max rms file!\n");
504 return(-124);
505 };
506 //
507 Int_t vview = 0;
508 Int_t vplane = 0;
509 Int_t vval = 0;
510 for (Int_t l=0; l<2; l++){
511 for (Int_t m=0; m<22; m++){
512 maxrms[l][m] = 26;
513 };
514 };
515 while ( spfile >> vview && spfile >> vplane && spfile >> vval){
516 if ( debug ) printf(" Setting view %i plane %i max rms %i \n",vview,vplane,vval);
517 maxrms[vview][vplane] = vval;
518 };
519 spfile.close();
520 //
521 };
522 //
523 // calopar5
524 //
525 if ( calopar5 || ( ttcalopar5 != 0 && ttcalopar5 < runheader ) ){
526 //
527 if ( debug ) printf(" calopar5 %i ftcalopar5 %u ttcalopar5 %u runheader %u \n",calopar5,ftcalopar5,ttcalopar5,runheader);
528 //
529 calopar5 = false;
530 //
531 // flight extra corrections:
532 //
533 if ( verbose ) printf(" Querying DB for calorimeter noise to signal threshold file...\n");
534 //
535 error = 0;
536 error = glparam->Query_GL_PARAM(runheader,111,dbc);
537 if ( error < 0 ) return(error);
538 //
539 calfile.str("");
540 calfile << glparam->PATH.Data() << "/";
541 calfile << glparam->NAME.Data();
542 ftcalopar5 = glparam->FROM_TIME;
543 ttcalopar5 = glparam->TO_TIME;
544 //
545 if ( verbose ) printf("\n Using calorimeter noise to signal threshold file: \n %s \n",calfile.str().c_str());
546 ifstream spfile;
547 spfile.open(calfile.str().c_str());
548 if ( !spfile ){
549 if ( verbose ) printf(" CALORIMETER - ERROR: no noise to signal threshold file!\n");
550 return(-125);
551 };
552 //
553 Int_t vview = 0;
554 Int_t vplane = 0;
555 Int_t vstrip = 0;
556 Float_t vval = 0.;
557 for (Int_t l=0; l<2; l++){
558 for (Int_t m=0; m<22; m++){
559 for (Int_t n=0; n<96; n++){
560 memin[l][m][n] = 0.7;
561 };
562 };
563 };
564 while ( spfile >> vview && spfile >> vplane && spfile >> vstrip && spfile >> vval){
565 if ( vstrip == -1 ){
566 for (Int_t ll=0; ll<96; ll++){
567 if ( debug ) printf(" Setting view %i plane %i strip %i noise to signal ratio %f \n",vview,vplane,ll,vval);
568 memin[vview][vplane][ll] = vval;
569 };
570 } else {
571 if ( debug ) printf(" Setting view %i plane %i strip %i noise to signal ratio %f \n",vview,vplane,vstrip,vval);
572 memin[vview][vplane][vstrip] = vval;
573 };
574 };
575 spfile.close();
576 //
577 };
578 //
579 //
580 delete glparam;
581 //
582 return(0);
583 }
584
585 Int_t CaloLevel0::CalcCrossTalkCorr(TSQLServer *dbc, UInt_t runheader, Bool_t ctusetable){
586 Int_t sig = CalcCrossTalkCorrDo(dbc,runheader,ctusetable);
587 return(sig);
588 };
589
590 Int_t CaloLevel0::CalcCrossTalkCorr(TSQLServer *dbc, UInt_t runheader){
591 Int_t sig = CalcCrossTalkCorrDo(dbc,runheader,true);
592 return(sig);
593 }
594
595 Int_t CaloLevel0::CalcCrossTalkCorr(GL_TABLES *glt, UInt_t runheader, Bool_t usetable){
596 const TString host = glt->CGetHost();
597 const TString user = glt->CGetUser();
598 const TString psw = glt->CGetPsw();
599 TSQLServer *dbc = TSQLServer::Connect(host.Data(),user.Data(),psw.Data());
600 if ( !dbc->IsConnected() ) throw -116;
601 stringstream myquery;
602 myquery.str("");
603 myquery << "SET time_zone='+0:00';";
604 delete dbc->Query(myquery.str().c_str());
605 delete dbc->Query("SET sql_mode = 'NO_UNSIGNED_SUBTRACTION';");
606 //
607 Int_t sig = CalcCrossTalkCorrDo(dbc,runheader,usetable);
608 dbc->Close();
609 delete dbc;
610 dbc = 0;
611 //
612 return(sig);
613 //
614 };
615
616 Int_t CaloLevel0::CalcCrossTalkCorr(GL_TABLES *glt, UInt_t runheader){
617 const TString host = glt->CGetHost();
618 const TString user = glt->CGetUser();
619 const TString psw = glt->CGetPsw();
620 TSQLServer *dbc = TSQLServer::Connect(host.Data(),user.Data(),psw.Data());
621 if ( !dbc->IsConnected() ) throw -116;
622 stringstream myquery;
623 myquery.str("");
624 myquery << "SET time_zone='+0:00';";
625 delete dbc->Query(myquery.str().c_str());
626 delete dbc->Query("SET sql_mode = 'NO_UNSIGNED_SUBTRACTION';");
627 //
628 Int_t sig = CalcCrossTalkCorrDo(dbc,runheader,true);
629 dbc->Close();
630 delete dbc;
631 dbc = 0;
632 //
633 return(sig);
634 //
635 }
636
637 Int_t CaloLevel0::CalcCrossTalkCorrDo(TSQLServer *dbc, UInt_t runheader, Bool_t usetable){
638 //
639 if ( ctground ) return(0);
640 //
641 Int_t error = 0;
642 GL_PARAM *glparam = new GL_PARAM();
643 //
644 // determine where I can find file with offline bad pulser mask
645 //
646 stringstream bmfile;
647 error = 0;
648 error = glparam->Query_GL_PARAM(runheader,105,dbc);
649 if ( error < 0 ) return(error);
650 //
651 bmfile.str("");
652 bmfile << glparam->PATH.Data() << "/";
653 bmfile << glparam->NAME.Data();
654 //
655 ifstream badfile;
656 if ( verbose ) printf("\n Using bad pulser offline mask file: \n %s \n\n",bmfile.str().c_str());
657 badfile.open(bmfile.str().c_str());
658 if ( !badfile ){
659 if ( verbose ) printf(" CALORIMETER - ERROR: no bad pulser offline mask file!\n");
660 return(-115);
661 };
662 //
663 Bool_t isdone = false;
664 Int_t bad = 0;
665 Int_t view = 1;
666 Int_t pre = 0;
667 Int_t plane = 21;
668 while ( !isdone ) {
669 badfile >> bad;
670 obadpulsemask[view][plane][pre] = bad;
671 if ( debug && bad ) printf(" SETTING view %i plane %i pre %i BAD = %i \n",view,plane,pre,bad);
672 pre++;
673 if ( pre > 5 ){
674 pre = 0;
675 plane--;
676 if ( plane < 0 ){
677 plane = 21;
678 view--;
679 };
680 if ( view < 0 ) isdone = true;
681 };
682 };
683 //
684 badfile.close();
685 if ( !usetable ){
686 //
687 // Let's start with cross-talk correction calculation
688 //
689 GL_CALOPULSE_CALIB *glp = new GL_CALOPULSE_CALIB();
690 Float_t adcp[2][22][96];
691 Float_t adcpcal[2][22][96];
692 memset(adcp , 0, 2*22*96*sizeof(Float_t));
693 memset(adcpcal , 0, 2*22*96*sizeof(Float_t));
694 //
695 UInt_t pampli = 0;
696 for (Int_t s=0; s<4; s++){
697 //
698 // Save into matrix adcp the values of the highest pulse calibration (pulse amplitude = 2)
699 //
700 pampli = 2;
701 error = 0;
702 error = glp->Query_GL_CALOPULSE_CALIB(runheader,s,pampli,dbc);
703 if ( error < 0 ){
704 if ( verbose ) printf(" CALORIMETER - ERROR: error from GLTables\n");
705 return(error);
706 };
707 //
708 UInt_t idcalib = glp->ID_ROOT_L0;
709 UInt_t fromtime = glp->FROM_TIME;
710 UInt_t calibno = glp->EV_ROOT;
711 //
712 // determine path and name and entry of the calibration file
713 //
714 GL_ROOT *glroot = new GL_ROOT();
715 if ( verbose ) printf("\n");
716 if ( verbose ) printf(" ** SECTION %i **\n",s);
717 //
718 error = 0;
719 error = glroot->Query_GL_ROOT(idcalib,dbc);
720 if ( error < 0 ){
721 if ( verbose ) printf(" CALORIMETER - ERROR: error from GLTables\n");
722 return(error);
723 };
724 //
725 stringstream name;
726 name.str("");
727 name << glroot->PATH.Data() << "/";
728 name << glroot->NAME.Data();
729 //
730 TString fcalname = (TString)name.str().c_str();
731 ifstream myfile;
732 myfile.open(fcalname.Data());
733 if ( !myfile ){
734 return(-107);
735 };
736 myfile.close();
737 //
738 TFile *File = new TFile(fcalname.Data());
739 if ( !File ) return(-108);
740 TTree *tr = (TTree*)File->Get("CalibCalPulse2");
741 if ( !tr ) return(-119);
742 //
743 TBranch *calo = tr->GetBranch("CalibCalPulse2");
744 //
745 pamela::CalibCalPulse2Event *ce = 0;
746 tr->SetBranchAddress("CalibCalPulse2", &ce);
747 //
748 Long64_t ncalibs = calo->GetEntries();
749 //
750 if ( !ncalibs ) return(-110);
751 //
752 if ( calo->GetEntry(calibno) <= 0) throw -36;
753 if ( verbose ) printf(" PULSE2 using entry %u from file %s",calibno,fcalname.Data());
754 //
755 // retrieve calibration table
756 //
757 if ( ce->pstwerr[s] && ce->pperror[s] == 0 && ce->unpackError == 0 ){
758 for ( Int_t d=0 ; d<11 ;d++ ){
759 for ( Int_t j=0; j<96 ;j++){
760 if ( s == 2 ){
761 adcp[0][2*d+1][j] = ce->calpuls[3][d][j];
762 };
763 if ( s == 3 ){
764 adcp[0][2*d][j] = ce->calpuls[1][d][j];
765 };
766 if ( s == 0 ){
767 adcp[1][2*d][j] = ce->calpuls[0][d][j];
768 };
769 if ( s == 1 ){
770 adcp[1][2*d+1][j] = ce->calpuls[2][d][j];
771 };
772 };
773 };
774 } else {
775 if ( verbose ) printf(" CALORIMETER - ERROR: problems finding a good calibration in this file! \n\n ");
776 return(-111);
777 };
778 //
779 File->Close();
780 delete glroot;
781 //
782 // Save into matrix adcpcal the calibrated values of the pulse calibration (subtraction of pulse amplitude = 0 relative to the pulse2 calibration used)
783 //
784 pampli = 0;
785 error = 0;
786 error = glp->Query_GL_CALOPULSE_CALIB(fromtime,s,pampli,dbc);
787 if ( error < 0 ){
788 if ( verbose ) printf(" CALORIMETER - ERROR: error from GLTables\n");
789 return(error);
790 };
791 //
792 idcalib = glp->ID_ROOT_L0;
793 calibno = glp->EV_ROOT;
794 //
795 // determine path and name and entry of the calibration file
796 //
797 glroot = new GL_ROOT();
798 if ( verbose ) printf("\n");
799 if ( verbose ) printf(" ** SECTION %i **\n",s);
800 //
801 error = 0;
802 error = glroot->Query_GL_ROOT(idcalib,dbc);
803 if ( error < 0 ){
804 if ( verbose ) printf(" CALORIMETER - ERROR: error from GLTables\n");
805 return(error);
806 };
807 //
808 name.str("");
809 name << glroot->PATH.Data() << "/";
810 name << glroot->NAME.Data();
811 //
812 fcalname = (TString)name.str().c_str();
813 myfile.open(fcalname.Data());
814 if ( !myfile ){
815 return(-107);
816 };
817 myfile.close();
818 //
819 TFile *File1 = new TFile(fcalname.Data());
820 if ( !File1 ) return(-108);
821 TTree *tr1 = (TTree*)File1->Get("CalibCalPulse1");
822 if ( !tr1 ) return(-120);
823 //
824 TBranch *calo1 = tr1->GetBranch("CalibCalPulse1");
825 //
826 pamela::CalibCalPulse1Event *ce1 = 0;
827 tr1->SetBranchAddress("CalibCalPulse1", &ce1);
828 //
829 ncalibs = calo1->GetEntries();
830 //
831 if ( !ncalibs ) return(-110);
832 //
833 if ( calo1->GetEntry(calibno) <= 0 ) throw -36;
834 if ( verbose ) printf(" PULSE1 using entry %u from file %s",calibno,fcalname.Data());
835 //
836 // retrieve calibration table
837 //
838 if ( ce1->pstwerr[s] && ce1->pperror[s] == 0 && ce1->unpackError == 0 ){
839 for ( Int_t d=0 ; d<11 ;d++ ){
840 for ( Int_t j=0; j<96 ;j++){
841 if ( s == 2 ){
842 adcpcal[0][2*d+1][j] = adcp[0][2*d+1][j] - ce1->calpuls[3][d][j];
843 };
844 if ( s == 3 ){
845 adcpcal[0][2*d][j] = adcp[0][2*d][j] - ce1->calpuls[1][d][j];
846 };
847 if ( s == 0 ){
848 adcpcal[1][2*d][j] = adcp[1][2*d][j] - ce1->calpuls[0][d][j];
849 };
850 if ( s == 1 ){
851 adcpcal[1][2*d+1][j] = adcp[1][2*d+1][j] - ce1->calpuls[2][d][j];
852 };
853 };
854 };
855 } else {
856 if ( verbose ) printf(" CALORIMETER - ERROR: problems finding a good calibration in this file! \n\n ");
857 return(-111);
858 };
859 //
860 File1->Close();
861 //
862 delete glroot;
863 //
864 };// loop on the four sections
865 //
866 //
867 delete glp;
868 //
869 // Ok, now we can try to calculate the cross-talk correction for each pre-amplifier
870 //
871 for ( Int_t v=0; v<2; v++){
872 if ( debug ) printf(" \n\n NEW VIEW \n");
873 for ( Int_t p=0; p<22; p++){
874 for ( Int_t npre=0; npre<6; npre++){
875 ctprecor[v][p][npre] = 1000.;
876 ctneigcor[v][p][npre] = 1000.;
877 Int_t str0=npre*16;
878 Int_t str16= -1 + (1+npre)*16;
879 //
880 UInt_t neigc = 0;
881 UInt_t farc = 0;
882 UInt_t pulsc = 0;
883 Float_t sigpulsed = 0.;
884 Float_t neigbase = 0.;
885 Float_t farbase = 0.;
886 //
887 // Loop over the strip of the pre and sum all signal far away from pulsed strip, signal in the neighbour(s) strip(s) and save the pulsed signal
888 // moreover count the number of strips in each case
889 //
890 for (Int_t s=str0; s<=str16; s++){
891 if ( adcpcal[v][p][s] > 10000.){
892 sigpulsed = adcpcal[v][p][s];
893 pulsc++;
894 if ( s > str0 ){
895 neigbase += adcpcal[v][p][s-1];
896 neigc++;
897 farbase -= adcpcal[v][p][s-1];
898 farc--;
899 };
900 if ( s < str16 ){
901 neigbase += adcpcal[v][p][s+1];
902 neigc++;
903 farbase -= adcpcal[v][p][s+1];
904 farc--;
905 };
906 } else {
907 farc++;
908 farbase += adcpcal[v][p][s];
909 };
910 };
911 //
912 // Now calculate the corrections
913 //
914 Float_t avefarbase = 0.;
915 if ( farc ) avefarbase = farbase/(Float_t)farc;
916 Float_t aveneigbase = 0.;
917 if ( neigc ) aveneigbase = neigbase/(Float_t)neigc;
918 //
919 if ( pulsc == 1 && farc && neigc ){
920 ctprecor[v][p][npre] = -avefarbase/(sigpulsed+fabs(avefarbase));
921 ctneigcor[v][p][npre] = fabs(aveneigbase-avefarbase)/(sigpulsed+fabs(avefarbase));
922 if ( debug ) printf(" Cross-talk correction View %i Plane %i Pre %i : pre-correction: %f neighbour strips correction %f \n",v,p,npre,ctprecor[v][p][npre],ctneigcor[v][p][npre]);
923 } else {
924 //
925 // did not find the pulsed strip or more than one pulsed strip found!
926 //
927 if ( debug ) printf(" Problems finding the cross-talk corrections: \n View %i Plane %i Pre %i number of pulsed strip %i \n Average faraway baseline %f number of strips %i Average neighbour baseline %f number of neighbour strips %i \n",v,p,npre,pulsc,avefarbase,farc,aveneigbase,neigc);
928 //
929 };
930 };
931 if ( debug ) printf(" \n ==================== \n");
932 };
933 };
934 } else {
935 //
936 // use pre-amply table
937 //
938 //
939 // determine where I can find file with offline neighbour correction table
940 //
941 stringstream bmfile2;
942 error = 0;
943 error = glparam->Query_GL_PARAM(runheader,106,dbc);
944 if ( error < 0 ) return(error);
945 //
946 bmfile2.str("");
947 bmfile2 << glparam->PATH.Data() << "/";
948 bmfile2 << glparam->NAME.Data();
949 //
950 ifstream badfile2;
951 if ( verbose ) printf("\n Using pre-amply neighbour crosstalk table file: \n %s \n\n",bmfile2.str().c_str());
952 badfile2.open(bmfile2.str().c_str());
953 if ( !badfile2 ){
954 if ( verbose ) printf(" CALORIMETER - ERROR: no pre-amply neighbour crosstalk table file!\n");
955 return(-121);
956 };
957 //
958 Int_t vview = 0;
959 Int_t vplane = 0;
960 Int_t vpre = 0;
961 Float_t vcorr = 0.;
962 while ( badfile2 >> vview && badfile2 >> vplane && badfile2 >> vpre && badfile2 >> vcorr){
963 if ( debug ) printf(" Pre-amply neighbour correction: view %i plane %i pre %i correction %f \n",vview,vplane,vpre,vcorr);
964 ctneigcor[vview][vplane][vpre] = vcorr;
965 };
966 //
967 // determine where I can find file with offline SECOND neighbour correction table
968 //
969 stringstream bmfile3;
970 error = 0;
971 error = glparam->Query_GL_PARAM(runheader,107,dbc);
972 if ( error < 0 ) return(error);
973 //
974 bmfile3.str("");
975 bmfile3 << glparam->PATH.Data() << "/";
976 bmfile3 << glparam->NAME.Data();
977 //
978 ifstream badfile3;
979 if ( verbose ) printf("\n Using pre-amply second neighbour crosstalk table file: \n %s \n\n",bmfile3.str().c_str());
980 badfile3.open(bmfile3.str().c_str());
981 if ( !badfile3 ){
982 if ( verbose ) printf(" CALORIMETER - ERROR: no pre-amply second neighbour crosstalk table file!\n");
983 return(-122);
984 };
985 //
986 Int_t pview = 0;
987 Int_t pplane = 0;
988 Int_t ppre = 0;
989 Float_t pcorr = 0.;
990 while ( badfile3 >> pview && badfile3 >> pplane && badfile3 >> ppre && badfile3 >> pcorr){
991 if ( debug ) printf(" Pre-amply second neighbour correction: view %i plane %i pre %i correction %f \n",pview,pplane,ppre,-pcorr);
992 ctprecor[pview][pplane][ppre] = -pcorr; // data are saved as negatives in the file
993 };
994 //
995 // determine where to find the file containing the Silicon crosstalk correction table
996 //
997 stringstream bmfile4;
998 error = 0;
999 error = glparam->Query_GL_PARAM(runheader,108,dbc);
1000 if ( error < 0 ) return(error);
1001 //
1002 bmfile4.str("");
1003 bmfile4 << glparam->PATH.Data() << "/";
1004 bmfile4 << glparam->NAME.Data();
1005 //
1006 ifstream badfile4;
1007 if ( verbose ) printf("\n Using Silicon crosstalk table file: \n %s \n\n",bmfile4.str().c_str());
1008 badfile4.open(bmfile4.str().c_str());
1009 if ( !badfile4 ){
1010 if ( verbose ) printf(" CALORIMETER - ERROR: no Silicon crosstalk table file!\n");
1011 return(-125);
1012 };
1013 //
1014 Int_t spview = 0;
1015 Int_t spplane = 0;
1016 Int_t psil = 0;
1017 Float_t spcorr = 0.;
1018 memset(ctsicor, 0, 2*22*9*sizeof(Float_t));
1019 while ( badfile4 >> spview && badfile4 >> spplane && badfile4 >> psil && badfile4 >> spcorr){
1020 if ( debug ) printf(" Silicon correction: view %i plane %i silicon %i correction %f \n",spview,spplane,psil,-spcorr);
1021 ctsicor[spview][spplane][psil] = -spcorr; // data are saved as negatives in the file
1022 };
1023 //
1024 };
1025 //
1026 delete glparam;
1027 //
1028 // Check the calculated corrections
1029 //
1030 Int_t opre=0;
1031 Int_t ppre=0;
1032 Bool_t found = false;
1033 for ( Int_t v=0; v<2; v++){
1034 for ( Int_t p=0; p<22; p++){
1035 for ( Int_t npre=0; npre<6; npre++){
1036 if ( ctprecor[v][p][npre] == 1000. || ctneigcor[v][p][npre] == 1000. || obadpulsemask[v][p][npre] != 0 ){
1037 if ( debug ) printf(" Cross-talk correction CHANGED for view %i Plane %i Pre %i\n BEFORE: pre-correction: %f neighbour strips correction %f \n",v,p,npre,ctprecor[v][p][npre],ctneigcor[v][p][npre]);
1038 if ( npre%2 ){
1039 opre = npre-1;
1040 } else {
1041 opre = npre+1;
1042 };
1043 if ( ctprecor[v][p][opre] == 1000. || ctneigcor[v][p][opre] == 1000. || obadpulsemask[v][p][opre] != 0 ){
1044 ppre=0;
1045 found = false;
1046 while ( ppre < 6 ){
1047 if ( ctprecor[v][p][ppre] != 1000. && ctneigcor[v][p][ppre] != 1000. && !obadpulsemask[v][p][ppre] ){
1048 found = true;
1049 ctprecor[v][p][npre] = ctprecor[v][p][ppre];
1050 ctneigcor[v][p][npre] = ctneigcor[v][p][ppre];
1051 break;
1052 };
1053 ppre++;
1054 };
1055 if ( !found ){
1056 if ( verbose ) printf(" WARNING: cannot find a good cross-talk correction for view %i plane %i pre %i \n Setting to default values 0.002 0.002\n",v,p,npre);
1057 ctprecor[v][p][npre] = 0.002;
1058 ctneigcor[v][p][npre] = 0.002;
1059 };
1060 } else {
1061 ctprecor[v][p][npre] = ctprecor[v][p][opre];
1062 ctneigcor[v][p][npre] = ctneigcor[v][p][opre];
1063 };
1064 if ( debug ) printf(" AFTER: pre-correction: %f neighbour strips correction %f \n",ctprecor[v][p][npre],ctneigcor[v][p][npre]);
1065 };
1066 };
1067 };
1068 };
1069 //
1070 return(0);
1071 }
1072
1073 void CaloLevel0::FindBaseCompress(Int_t l, Int_t m, Int_t pre){
1074 Int_t n = 0;
1075 Float_t q = 0;
1076 this->FindBaseCompress(l,m,pre,n,q);
1077 }
1078
1079 void CaloLevel0::FindBaseCompress(Int_t l, Int_t m, Int_t pre, Int_t &nst, Float_t &qp){
1080 for (Int_t e = pre*16; e < (pre+1)*16 ; e++){
1081 dexy[l][m][e] = dexyc[l][m][e];
1082 };
1083 this->FindBaseRaw(l,m,pre,nst,qp);
1084 }
1085
1086 void CaloLevel0::FindBaseRaw(Int_t l, Int_t m, Int_t pre){
1087 Int_t n = 0;
1088 Float_t q = 0;
1089 this->FindBaseRaw(l,m,pre,n,q);
1090 }
1091
1092 void CaloLevel0::FindBaseRaw(Int_t l, Int_t m, Int_t pre, Int_t &nst, Float_t &qp){
1093 //
1094 Float_t minstrip = 100000.;
1095 Float_t rms = 0.;
1096 Int_t process = 0;
1097 Int_t onlmask[16];
1098 memset(onlmask, 0, 16*sizeof(Int_t));
1099 //
1100 while ( process < 2 ){
1101 //
1102 minstrip = 100000.;
1103 rms = 0.;
1104 base[l][m][pre] = 0.;
1105 qp = 0.;
1106 //
1107 Int_t spos = -1;
1108 Int_t ee = 0;
1109 for (Int_t e = pre*16; e < (pre+1)*16 ; e++){
1110 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. && onlmask[ee] == 0 ) {
1111 minstrip = dexy[l][m][e]-calped[l][m][e];
1112 rms = calthr[l][m][pre];
1113 spos = ee;
1114 };
1115 ee++;
1116 qp += (dexy[l][m][e]-calped[l][m][e]-sbase[l][m][e]);
1117 };
1118 //
1119 if ( debug && l==0 ){
1120 printf("\n BASELINE CALCULATION for view %i pl %i pre %i: \n => minstrip %f rms %f \n => qp = %f \n",l,m,pre,minstrip,rms,qp);
1121 };
1122 if ( minstrip != 100000. ) {
1123 Float_t strip6s = 0.;
1124 for (Int_t e = pre*16; e < (pre+1)*16 ; e++){
1125 if ( (dexy[l][m][e]-calped[l][m][e]) >= minstrip && (dexy[l][m][e]-calped[l][m][e]) <= (minstrip+rms) ) {
1126 strip6s += 1.;
1127 base[l][m][pre] += (dexy[l][m][e] - calped[l][m][e]);
1128 };
1129 //
1130 // compression
1131 //
1132 // if ( abs((int)(dexy[l][m][e]-calped[l][m][e])) <= (minstrip+rms) ) {
1133 // dexyc[l][m][e] = 0.;
1134 // } else {
1135 dexyc[l][m][e] = dexy[l][m][e];
1136 // };
1137 };
1138 //
1139 if ( strip6s == 1. && process < 1 ){
1140 onlmask[spos] = 1;
1141 process++;
1142 if ( debug ) printf(" Warning, only one strip to calculate baseline: minstrip %f rms %f spos %i l %i m %i pre %i \n",minstrip,rms,spos,l,m,pre);
1143 continue;
1144 };
1145 process += 2;
1146 nst = (Int_t)strip6s;
1147 //
1148 if ( debug ){
1149 printf(" strip6s %f \n",strip6s);
1150 };
1151 // if ( strip6s >= 9. ){
1152 if ( (strip6s >= 2. && process == 2) || (strip6s >= 9. && process > 2) ){
1153 //if ( (strip6s >= 4. && process == 2) || (strip6s >= 9. && process > 2) ){
1154 Double_t arro = base[l][m][pre]/strip6s;
1155 Float_t deci = 1000.*((float)arro - float(int(arro)));
1156 if ( deci < 500. ) {
1157 arro = double(int(arro));
1158 } else {
1159 arro = 1. + double(int(arro));
1160 };
1161 base[l][m][pre] = arro;
1162 //
1163 // if too few strips were used to determine the baseline check if it is comparable with the previous event, if not mark it as bad
1164 //
1165 if ( debug && process > 2 ) printf(" AGH low strip value was discarded process %i strip6s %f minstrip %f rms %f spos %i\n",process,strip6s,minstrip,rms,spos);
1166 if ( debug ) printf(" Calculated baseline: base %f sbase-0.02*qp %f \n",base[l][m][pre],(-qp*0.02+sbase[l][m][pre]));
1167 //
1168 if ( strip6s < 4 && base[l][m][pre] > (-0.015*qp+sbase[l][m][pre]) && sbase[l][m][pre] > 0. ){
1169 if ( debug ) printf(" Suspicious calculated baseline: base %f sbase-0.02*qp %f strip6s %i \n",base[l][m][pre],(-qp*0.02+sbase[l][m][pre]),(Int_t)strip6s);
1170 base[l][m][pre] = 31000.;
1171 nst = 0; // 9RED BUG
1172 qp = 0.; // 9RED BUG
1173 for (Int_t e = pre*16; e < (pre+1)*16 ; e++){
1174 dexyc[l][m][e] = dexy[l][m][e];
1175 };
1176 };
1177 } else {
1178 if ( debug ) printf(" reset baseline here if ! ( (strip6s >=2 && process == 2) || (strip6s >= 9 and process > 2) ) \n");
1179 base[l][m][pre] = 31000.;
1180 nst = 0; // 9RED BUG
1181 qp = 0.; // 9RED BUG
1182 for (Int_t e = pre*16; e < (pre+1)*16 ; e++){
1183 dexyc[l][m][e] = dexy[l][m][e];
1184 };
1185 };
1186 } else {
1187 if ( debug ) printf(" reset baseline here if no minimum find\n");
1188 nst = 0; // 9RED BUG
1189 qp = 0.; // 9RED BUG
1190 process += 2;
1191 base[l][m][pre] = 31000.;
1192 for (Int_t e = pre*16; e < (pre+1)*16 ; e++){
1193 dexyc[l][m][e] = dexy[l][m][e];
1194 };
1195 };
1196 };
1197 if ( debug ) printf(" Baseline calculation: baseline for view %i plane %i pre %i is %f nst %i qp %f \n",l,m,pre,base[l][m][pre],nst,qp);
1198 }
1199
1200 Int_t CaloLevel0::Calibrate(Int_t ei){
1201 //
1202 // get entry ei
1203 //
1204 if ( l0calo->GetEntry(ei) <= 0 ) throw -36;
1205 //
1206 // if it was not a selftrigger event, could it ever been a selftrigger event? if so trigty = 3.
1207 //
1208 clevel2->nsatstrip = 0.;
1209 Int_t val = 0;
1210 Int_t del = 1000;
1211 for (Int_t sec = 0; sec < 4; sec++){
1212 for (Int_t dsec = 0; dsec < 7; dsec++){
1213 val = (Int_t)de->calselftrig[sec][dsec];
1214 del = delay(val);
1215 clevel2->selfdelay[sec][dsec] = del;
1216 };
1217 };
1218 val = 0;
1219 del = 1000;
1220 if ( clevel2->trigty < 2. ){
1221 // Bool_t ck = false;
1222 for (Int_t sec = 0; sec < 4; sec++){
1223 val = (Int_t)de->calselftrig[sec][6];
1224 del = delay(val);
1225 if ( del < 1000 ){
1226 clevel2->wartrig = 0.;
1227 clevel2->trigty = 3.;
1228 // ck = true;
1229 break;
1230 };
1231 };
1232 // if ( !ck ) clevel2->wartrig = 100.;
1233 } else {
1234 Bool_t ck = false;
1235 for (Int_t sec = 0; sec < 4; sec++){
1236 val = (Int_t)de->calselftrig[sec][6];
1237 del = delay(val);
1238 if ( del < 1000 ){
1239 clevel2->wartrig = 0.;
1240 ck = true;
1241 };
1242 };
1243 if ( !ck ) clevel2->wartrig = 100.;
1244 };
1245 //
1246 Int_t se = 5;
1247 Int_t done = 0;
1248 Int_t pre = -1;
1249 // Bool_t isCOMP = false;
1250 // Bool_t isFULL = false;
1251 Bool_t isRAW = false;
1252 Float_t ener;
1253 Int_t doneb = 0;
1254 Int_t donec = 0;
1255 Int_t ck[2][22][6];
1256 memset(ck, 0, 2*22*6*sizeof(Int_t));
1257 Int_t ipre = 0;
1258 // Int_t ip[3] = {0};
1259 Int_t ip[3] = {0,0,0};
1260 Int_t ipp = 0;
1261 Float_t base0, base1, base2;
1262 base0 = 0.;
1263 base1 = 0.;
1264 base2 = 0.;
1265 Float_t qpre[2][22][6];
1266 memset(qpre, 0, 2*22*6*sizeof(Float_t));
1267 Float_t ene[96];
1268 Int_t chdone[4] = {0,0,0,0};
1269 Int_t pe = 0;
1270 //
1271 Float_t ener0 = 0.;
1272 Float_t cbase0 = 0.;
1273 Float_t totbase = 0.;
1274 Float_t totped = 0.;
1275 Bool_t pproblem = false;
1276 Bool_t negbase = false;
1277 //
1278 Float_t tim = 0.;
1279 Int_t plo = 0;
1280 Int_t fbi = 0;
1281 Int_t cle = 0;
1282 //
1283 // run over views and planes
1284 //
1285 for (Int_t l = 0; l < 2; l++){
1286 for (Int_t m = 0; m < 22; m++){
1287 //
1288 // determine the section number
1289 //
1290 negbase = false;
1291 se = 5;
1292 if (l == 0 && m%2 == 0) se = 3;
1293 if (l == 0 && m%2 != 0) se = 2;
1294 if (l == 1 && m%2 != 0) se = 1;
1295 if (l == 1 && m%2 == 0) se = 0;
1296 //
1297 // determine what kind of event we are going to analyze
1298 //
1299 // isCOMP = false;
1300 // isFULL = false;
1301 isRAW = false;
1302 // if ( de->stwerr[se] & (1 << 16) ) isCOMP = true;
1303 // if ( de->stwerr[se] & (1 << 17) ) isFULL = true;
1304 if ( de->stwerr[se] & (1 << 3) ) isRAW = true;
1305 if ( !chdone[se] ){
1306 //
1307 // check for any error in the event
1308 //
1309 clevel2->crc[se] = 0;
1310 if ( de->perror[se] == 132 ){
1311 clevel2->crc[se] = 1;
1312 pe++;
1313 };
1314 clevel2->perr[se] = 0;
1315 if ( de->perror[se] != 0 ){
1316 clevel2->perr[se] = (Int_t)de->perror[se];
1317 pe++;
1318 };
1319 clevel2->swerr[se] = 0;
1320 for (Int_t j = 0; j < 7 ; j++){
1321 if ( (j != 3) && (de->stwerr[se] & (1 << j)) ){
1322 clevel2->swerr[se] = 1;
1323 pe++;
1324 };
1325 };
1326 chdone[se] = 1;
1327 };
1328 if ( clevel2->crc[se] == 0 && (clevel1->good2 == 1 || clevel2->trigty >= 2) ){
1329 pre = -1;
1330 //
1331 for (Int_t nn = 0; nn < 96; nn++){
1332 // ene[nn] = 0.;
1333 dexy[l][m][nn] = de->dexy[l][m][nn] ;
1334 dexyc[l][m][nn] = de->dexyc[l][m][nn] ;
1335 };
1336 //
1337 // run over preamplifiers
1338 //
1339 pre = -1;
1340 cbase0 = 0.;
1341 Int_t nstt[2];
1342 Float_t rqp[2];
1343 for (Int_t i = 0; i < 3; i++){
1344 nstt[0] = 1000;
1345 nstt[1] = 1000;
1346 rqp[0] = 0.;
1347 rqp[1] = 0.;
1348 for (Int_t j = 0; j < 2; j++){
1349 pre = j + i*2;
1350 //
1351 // baseline check and calculation
1352 //
1353 if ( !isRAW ){
1354 //
1355 // if it is a compress event with fully transmitted pre try to calculate the baseline
1356 //
1357 if ( de->base[l][m][pre] != 0. && de->base[l][m][pre]<31000. ) {
1358 base[l][m][pre] = de->base[l][m][pre] ;
1359 } else {
1360 FindBaseCompress(l,m,pre,nstt[j],rqp[j]);
1361 };
1362 cbase0 += base[l][m][pre];
1363 } else {
1364 //
1365 // if it is a raw event calculate the baseline.
1366 //
1367 FindBaseRaw(l,m,pre,nstt[j],rqp[j]);
1368 cbase0 += base[l][m][pre];
1369 };
1370 };
1371 //
1372 // if we are able to calculate the baseline with more than 3 strips on one pre and not in the other one choose the pre with more calculated strips
1373 //
1374 if ( nstt[0] < 4 && nstt[1] >= 4 && nstt[0] != 1000 && nstt[1] != 1000 ) base[l][m][pre-1] = 31000.;
1375 if ( nstt[0] >= 4 && nstt[1] < 4 && nstt[0] != 1000 && nstt[1] != 1000 ) base[l][m][pre] = 31000.;
1376 // //
1377 // // if we are NOT able to calculate the baseline with more than 3 strips on both pres take the baseline (if any) of the one which has less energy
1378 // //
1379 // if ( nstt[0] < 4 && nstt[1] < 4 ){
1380 // if ( rqp[0] >= rqp[1] ) base[l][m][pre-1] = 31000.;
1381 // if ( rqp[0] < rqp[1] ) base[l][m][pre] = 31000.;
1382 // };
1383 };
1384 //
1385 // run over strips
1386 //
1387 pre = -1;
1388 ener0 = 0.;
1389 totbase = 0.;
1390 totped = 0.;
1391 for (Int_t i = 0 ; i < 3 ; i++){
1392 ip[i] = 0;
1393 for (Int_t n = i*32 ; n < (i+1)*32 ; n++){
1394 if (n%16 == 0) {
1395 done = 0;
1396 doneb = 0;
1397 donec = 0;
1398 pre++;
1399 ck[l][m][pre] = 0;
1400 qpre[l][m][pre] = 0.;
1401 };
1402 //
1403 // baseline check and calculation
1404 //
1405 // no suitable new baseline, use old ones!
1406 //
1407 if ( !done ){
1408 if ( debug ) printf(" l %i m %i pre %i ip[i] %i base %f base ip[i] %f sbase %f \n",l,m,pre,ip[i],base[l][m][pre],base[l][m][ip[i]],sbase[l][m][pre]);
1409 if ( (base[l][m][pre] == 31000. || base[l][m][pre] == 0.) ){
1410 ck[l][m][pre] = 1;
1411 if (pre%2 == 0) {
1412 ip[i] = pre + 1;
1413 } else {
1414 ip[i] = pre - 1;
1415 };
1416 if ( (base[l][m][ip[i]] == 31000. || base[l][m][ip[i]] == 0. || !crosst ) ){
1417 //
1418 ck[l][m][pre] = 2;
1419 if ( sbase[l][m][pre] == 31000. || sbase[l][m][pre] == 0. ) {
1420 ck[l][m][pre] = 3;
1421 };
1422 };
1423 };
1424 done = 1;
1425 };
1426 //
1427 // CALIBRATION ALGORITHM
1428 //
1429 if ( !doneb ){
1430 if ( debug ) printf(" ck[l][m][pre] is %i \n",ck[l][m][pre]);
1431 switch (ck[l][m][pre]) {
1432 case 0:
1433 base0 = base[l][m][pre];
1434 base2 = calbase[l][m][pre];
1435 if ( debug ) printf(" base0 = base l%i m%i pre%i = %f base2 = calbase l m pre = %f \n",l,m,pre,base[l][m][pre],calbase[l][m][pre]);
1436 break;
1437 case 1:
1438 base0 = base[l][m][ip[i]];
1439 base2 = calbase[l][m][ip[i]];
1440 if ( debug ) printf(" base0 = base l%i m%i ip(i)%i = %f base2 = calbase l m ip(i) = %f \n",l,m,ip[i],base[l][m][ip[i]],calbase[l][m][ip[i]]);
1441 break;
1442 case 2:
1443 base0 = sbase[l][m][pre];
1444 base2 = calbase[l][m][pre];
1445 if ( debug ) printf(" base0 = sbase l%i m%i pre%i = %f base2 = calbase l m pre = %f \n",l,m,pre,sbase[l][m][pre],calbase[l][m][pre]);
1446 break;
1447 case 3:
1448 base0 = calbase[l][m][pre];
1449 base2 = calbase[l][m][pre];
1450 if ( debug ) printf(" base0 = calbase l%i m%i pre%i = %f base2 = calbase l m pre = %f \n",l,m,pre,calbase[l][m][pre],calbase[l][m][pre]);
1451 break;
1452 };
1453 base1 = calbase[l][m][pre];
1454 doneb = 1;
1455 };
1456 ener = dexyc[l][m][n];
1457 ener0 += ener;
1458 clevel1->estrip[n][m][l] = 0.;
1459 totbase += de->base[l][m][pre]/96.;
1460 totped += fabs(calped[l][m][n]);
1461 if ( de->base[l][m][pre] < 0 ) negbase = true;
1462 if ( base0>0 && base0 < 30000. ){
1463 //
1464 // save the baseline only if the energy release is "small"
1465 //
1466 if ( !donec && (base0 + base1 - base2) != 0. && (n+1)%16==0 ){
1467 if ( qpre[l][m][pre] < 200. ) sbase[l][m][pre] = base0 + base1 - base2;
1468 donec = 1;
1469 };
1470 if ( ener > 0. ){
1471 clevel1->estrip[n][m][l] = (ener - calped[l][m][n] - base0 - base1 + base2)/mip[l][m][n] ;
1472 //
1473 // 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)
1474 //
1475 if ( debug && l==0 && (m==17 || m==18) ) printf(" view %i plane %i strip %i ener %f calped %f base0 %f base1 %f base2 %f mip %f ENERGIA %f \n",l,m,n,ener,calped[l][m][n],base0,base1,base2,mip[l][m][n],clevel1->estrip[n][m][l]);
1476 if ( clevel1->estrip[n][m][l] > 0. ) qpre[l][m][pre] += clevel1->estrip[n][m][l];
1477 //
1478 //
1479 };
1480 };
1481 };
1482 };
1483 //
1484 // check if there were problems with 5.7 or glitches in the power supply
1485 //
1486 // if ( ((ener0 == 0. && cbase0 == 0.) || negbase || totbase > 196600. || totped < 1. ) && !pproblem && clevel2->perr[se] == 0){ // check pedestal and baseline values for one plane, if all zeros calibration is not valid (calorimeter power problems) [8th data reduction bug, fixed on 25/11/2009 by E.M.]
1487 if ( ((ener0 == 0. && cbase0 == 0.) || negbase || totbase > 32700. || totped < 1. ) && !pproblem && clevel2->perr[se] == 0){ // check pedestal and baseline values for one plane, if all zeros calibration is not valid (calorimeter power problems) [8th data reduction bug, fixed on 25/11/2009 by E.M.]
1488 if ( verbose ) printf(" L0 entry %i : calorimeter power problems! event marked as bad perr %f swerr %X view %i plane %i negbase %i totbase %f totped %f\n",ei,de->perror[se],de->stwerr[se],l,m, negbase, totbase, totped);
1489 pproblem = true;
1490 pe++;
1491 };
1492 //
1493 } else {
1494 for (Int_t nn = 0; nn < 96; nn++){
1495 clevel1->estrip[nn][m][l] = 0.;
1496 };
1497 };
1498 };
1499 };
1500 //
1501 // run over views and planes to apply crosstalk corrections
1502 //
1503 for (Int_t l = 0; l < 2; l++){
1504 for (Int_t m = 0; m < 22; m++){
1505 //
1506 // determine the section number
1507 //
1508 se = 5;
1509 if (l == 0 && m%2 == 0) se = 3;
1510 if (l == 0 && m%2 != 0) se = 2;
1511 if (l == 1 && m%2 != 0) se = 1;
1512 if (l == 1 && m%2 == 0) se = 0;
1513 //
1514 // check for any error in the event
1515 //
1516 if ( clevel2->crc[se] == 0 && (clevel1->good2 == 1 || clevel2->trigty >= 2) ){
1517 //
1518 // Cross-talk corrections
1519 //
1520 if ( crosst ){
1521 //
1522 // energy on silicon ladders
1523 //
1524 Float_t qsi[3];
1525 qsi[0] = qpre[l][m][0]+qpre[l][m][1];
1526 qsi[1] = qpre[l][m][2]+qpre[l][m][3];
1527 qsi[2] = qpre[l][m][4]+qpre[l][m][5];
1528 //
1529 for ( pre = 1; pre < 6; pre += 2 ){
1530 Int_t ladder = (pre - 1)/2;
1531 //
1532 // If the noselfct flag is set the strip doesn't suffer the self crosstalk due to electronics so we must subtract some energy
1533 //
1534 if ( noselfct ){
1535 for (Int_t j = ladder*32 ; j < (ladder+1)*32 ; j++){
1536 ipre = j/16 ;
1537 if ( debug ) printf(" CT STEP1 %i %i %i estrip %f ctprecor %f \n",j,m,l,clevel1->estrip[j][m][l],ctprecor[l][m][ipre]);
1538 if ( clevel1->estrip[j][m][l] != 0. ) clevel1->estrip[j][m][l] -= clevel1->estrip[j][m][l] * ctprecor[l][m][ipre];
1539 if ( debug ) printf(" CT STEP2 %i %i %i estrip %f ctprecor %f \n",j,m,l,clevel1->estrip[j][m][l],ctprecor[l][m][ipre]);
1540 };
1541 };
1542 //
1543 // Using the neighbour pre baseline
1544 //
1545 if (ck[l][m][pre] == 1 || ck[l][m][pre-1] == 1){
1546 //
1547 // pre-amplifier effect on baseline when using the neighbour pre (ck=1)
1548 //
1549 if (ck[l][m][pre] == 1){
1550 ipre = pre;
1551 ipp = pre - 1;
1552 } else {
1553 ipre = pre - 1;
1554 ipp = pre;
1555 };
1556 Int_t it = 0;
1557 Float_t nqpre = 0.;
1558 //
1559 if ( debug ) printf(" CK1 Limit for while: 0.07 \n");
1560 for (Int_t j = ipre*16 ; j < (ipre+1)*16 ; j++){
1561 if ( debug ) printf(" CT STEP3 %i %i %i estrip %f ctprecor %f \n",j,m,l,clevel1->estrip[j][m][l],ctprecor[l][m][ipp]);
1562 if ( !ctground ){
1563 if ( clevel1->estrip[j][m][l] != 0. ) clevel1->estrip[j][m][l] += - qpre[l][m][ipp] * ctprecor[l][m][ipp];
1564 } else {
1565 if ( clevel1->estrip[j][m][l] != 0. ) clevel1->estrip[j][m][l] += - qpre[l][m][ipp] * 0.00478;
1566 };
1567 if ( clevel1->estrip[j][m][l] > 0. ) nqpre += clevel1->estrip[j][m][l] ;
1568 if ( debug ) printf(" CT STEP4 %i %i %i estrip %f ctprecor %f \n",j,m,l,clevel1->estrip[j][m][l],ctprecor[l][m][ipp]);
1569 };
1570 qpre[l][m][ipre] = nqpre;
1571 nqpre = 0.;
1572 Float_t deltaqpre = qpre[l][m][ipre];
1573 //
1574 // these values are empirically determined, usually the routine converge due to deltaqsi and the latest applied correction is based on less than 1 mip
1575 //
1576 while ( it < 10 && deltaqpre > 0.07 ){
1577 nqpre = 0.;
1578 for (Int_t j = ipre*16 ; j < (ipre+1)*16 ; j++){
1579 if ( debug ) printf(" CT STEP5 %i %i %i estrip %f ctprecor %f \n",j,m,l,clevel1->estrip[j][m][l],ctprecor[l][m][ipre]);
1580 if ( !ctground ){
1581 if ( debug ) printf(" CK1 pre correction: iteration %i deltaqpre %f ctprecor %f TOTAL CORRECTION %f \n",it,deltaqpre,ctprecor[l][m][ipre],deltaqpre * ctprecor[l][m][ipre]);
1582 if ( clevel1->estrip[j][m][l] != 0. ) clevel1->estrip[j][m][l] += deltaqpre * ctprecor[l][m][ipre];
1583 } else {
1584 if ( clevel1->estrip[j][m][l] != 0. ) clevel1->estrip[j][m][l] += deltaqpre * 0.00478;
1585 };
1586 if ( clevel1->estrip[j][m][l] > 0. ) nqpre += clevel1->estrip[j][m][l] ;
1587 if ( debug ) printf(" CT STEP6 %i %i %i estrip %f ctprecor %f \n",j,m,l,clevel1->estrip[j][m][l],ctprecor[l][m][ipre]);
1588 };
1589 if ( ctground ) it = 100;
1590 it++;
1591 deltaqpre = nqpre - qpre[l][m][ipre];
1592 if ( debug ) printf(" CK1 BEFORE: qpre %f \n",qpre[l][m][ipre]);
1593 qpre[l][m][ipre] = nqpre;
1594 if ( debug ) printf(" CK1 AFTER: qpre %f \n",qpre[l][m][ipre]);
1595 };
1596 //
1597 };
1598 //
1599 // No baseline calculation due to high energy release
1600 //
1601 if (ck[l][m][pre] == 2 && ck[l][m][pre-1] == 2){
1602 //
1603 // y^
1604 // |
1605 // | 6 7 8
1606 // | 3 4 5
1607 // | 0 1 2
1608 // | --------------------------------------> x
1609 //
1610 Int_t si1 = 0;
1611 Int_t si2 = 0;
1612 Int_t si3 = 0;
1613 if ( l == 0 ){
1614 if ( ladder == 0 ){
1615 si1 = 0;
1616 si2 = 3;
1617 si3 = 6;
1618 };
1619 if ( ladder == 1 ){
1620 si1 = 1;
1621 si2 = 4;
1622 si3 = 7;
1623 };
1624 if ( ladder == 2 ){
1625 si1 = 2;
1626 si2 = 5;
1627 si3 = 8;
1628 };
1629 } else {
1630 if ( ladder == 0 ){
1631 si1 = 0;
1632 si2 = 1;
1633 si3 = 2;
1634 };
1635 if ( ladder == 1 ){
1636 si1 = 3;
1637 si2 = 4;
1638 si3 = 5;
1639 };
1640 if ( ladder == 2 ){
1641 si1 = 6;
1642 si2 = 7;
1643 si3 = 8;
1644 };
1645 };
1646 //
1647 // Find the energy distribution along the considered plane looking at the two sandwiching plane of the other view.
1648 //
1649 Float_t sied[3] = {0.,0.,0.};
1650 Int_t othv = !l;
1651 Int_t othpl1 = m - 1;
1652 Int_t othpl2 = m + 1;
1653 Float_t oprof[3] = {0.,0.,0.};
1654 for(Int_t s=0; s<3; s++){
1655 for(Int_t t=(s*32); t<32*(s + 1); t++){
1656 if ( othpl1 > -1 ) {
1657 oprof[s] += clevel1->estrip[othv][othpl1][t];
1658 };
1659 if ( othpl2 < 22 ) {
1660 oprof[s] += clevel1->estrip[othv][othpl2][t];
1661 };
1662 };
1663 };
1664 Float_t otote = fabs(oprof[0]) + fabs(oprof[1]) + fabs(oprof[2]);
1665 for(Int_t g=0; g<3; g++){
1666 if ( otote > 0. ){
1667 sied[g] = fabs(oprof[g])/otote;
1668 } else {
1669 sied[g] = 1./3.;
1670 };
1671 };
1672 //
1673 //
1674 //
1675 Int_t it = 0;
1676 Int_t jpre = 0;
1677 Float_t nqsi = 0.;
1678 Float_t snqsi = qsi[ladder];
1679 Float_t nqpre[2] = {0.,0.};
1680 Float_t deltaqsi = qsi[ladder];
1681 Float_t deltaqpre[2];
1682 deltaqpre[0] = qpre[l][m][pre-1];
1683 deltaqpre[1] = qpre[l][m][pre];
1684 //
1685 if ( debug ) printf(" Limit for while: 0.07 it < 10 \n");
1686 //
1687 // these values are empirically determined, usually the routine converge due to deltaqsi and the latest applied correction is based on less than 1 mip
1688 //
1689 while ( it < 10 && (deltaqsi > 0.07 || deltaqpre[0] > 0.07 || deltaqpre[1] > 0.07) ){
1690 nqsi = 0.;
1691 nqpre[0] = 0.;
1692 nqpre[1] = 0.;
1693 for (Int_t j = ladder*32 ; j < (ladder+1)*32 ; j++){
1694 if ( debug ) printf(" CT STEP6 %i %i %i estrip %f ctprecor %f \n",j,m,l,clevel1->estrip[j][m][l],ctsicor[l][m][si2]);
1695 ipre = 0;
1696 if ( j > (ladder*32)+15 ) ipre = 1;
1697 jpre = j/16 ;
1698 //
1699 // Silicon effect on the baseline when using the same pre previous baseline (ck = 2) + pre-amply effect
1700 //
1701 if ( !ctground ){
1702 if ( debug ) printf(" silicon correction: iteration %i deltaqsi[%i] %f ctsicor %f %f %f sied %f %f %f si %i %i %i TOTAL CORRECTION %f \n",it,ladder,deltaqsi,ctsicor[l][m][si1],ctsicor[l][m][si2],ctsicor[l][m][si3],sied[0],sied[1],sied[2],si1,si2,si3,deltaqsi * (ctsicor[l][m][si1] * sied[0] + ctsicor[l][m][si2] * sied[1] + ctsicor[l][m][si3] * sied[2]));
1703 if ( debug ) printf(" pre correction: iteration %i deltaqpre[0] %f deltaqpre[1] %f ctprecor %f TOTAL CORRECTION %f \n",it,deltaqpre[0],deltaqpre[1],ctprecor[l][m][jpre],deltaqpre[ipre] * ctprecor[l][m][jpre]);
1704 if ( clevel1->estrip[j][m][l] != 0. ) clevel1->estrip[j][m][l] += (deltaqsi * (ctsicor[l][m][si1] * sied[0] + ctsicor[l][m][si2] * sied[1] + ctsicor[l][m][si3] * sied[2])/mip[l][m][j]) + deltaqpre[ipre] * ctprecor[l][m][jpre];
1705 } else {
1706 if ( clevel1->estrip[j][m][l] != 0. ) clevel1->estrip[j][m][l] += 0. + qpre[l][m][jpre] * 0.00478; // no correction
1707 };
1708 if ( clevel1->estrip[j][m][l] > 0. ) nqsi += clevel1->estrip[j][m][l] ;
1709 if ( clevel1->estrip[j][m][l] > 0. ) nqpre[ipre] += clevel1->estrip[j][m][l] ;
1710 if ( debug ) printf(" CT STEP7 %i %i %i estrip %f ctprecor %f \n",j,m,l,clevel1->estrip[j][m][l],ctsicor[l][m][si2]);
1711 };
1712 if ( ctground ) it = 100;
1713 deltaqsi = nqsi-snqsi;
1714 deltaqpre[0] = nqpre[0] - qpre[l][m][pre-1];
1715 deltaqpre[1] = nqpre[1] - qpre[l][m][pre];
1716 //
1717 // Check for divergence and stop if it happens! [9RED bug noticed with plane 18X]
1718 //
1719 if ( deltaqpre[0] > qpre[l][m][pre-1] || deltaqpre[1] > qpre[l][m][pre] || deltaqsi >snqsi ){
1720 if ( debug ) printf(" WARNING!! DIVERGING CORRECTION EXIT IMMEDIATLY FROM THE LOOP!! dqpre0 %f qpre0 %f // dqpre1 %f qpre1 %f // dqsi %f qsi %f \n",deltaqpre[0],qpre[l][m][pre-1],deltaqpre[1],qpre[l][m][pre],deltaqsi,snqsi);
1721 it = 1000;
1722 };
1723 //
1724 snqsi = nqsi;
1725 it++;
1726 if ( debug ) printf(" BEFORE: qpre 0 %f qpre 1 %f \n",qpre[l][m][pre-1],qpre[l][m][pre]);
1727 qpre[l][m][pre-1] = nqpre[0];
1728 qpre[l][m][pre] = nqpre[1];
1729 if ( debug ) printf(" AFTER: qpre 0 %f qpre 1 %f \n",qpre[l][m][pre-1],qpre[l][m][pre]);
1730 };
1731 //
1732 //
1733 //
1734 // for (Int_t j = ladder*32 ; j < (ladder+1)*32 ; j++){
1735 // ipre = j/16 ;
1736 // //
1737 // // pre-amplifier effect on baseline when using the same pre previous event baseline (ck=2)
1738 // //
1739 // if ( !ctground ){
1740 // if ( clevel1->estrip[j][m][l] != 0. ) clevel1->estrip[j][m][l] += qpre[l][m][ipre] * ctprecor[l][m][ipre];
1741 // } else {
1742 // if ( clevel1->estrip[j][m][l] != 0. ) clevel1->estrip[j][m][l] += qpre[l][m][ipre] * 0.00478;
1743 // };
1744 // };
1745 };
1746 };
1747 };
1748 };
1749 //
1750 Int_t j4 = -4;
1751 Int_t jjj = -3;
1752 Int_t jj = -2;
1753 Int_t jjpre = -1;
1754 Int_t jjjpre = -1;
1755 memset(ene, 0, 96*sizeof(Float_t));
1756 for (Int_t j = 0 ; j < 100 ; j++){
1757 jj++;
1758 jjj++;
1759 j4++;
1760 if ( j < 96 ) ene[j] = clevel1->estrip[j][m][l];
1761 if ( crosst ){
1762 //
1763 // "Real" crosstalk effect on the neighbour strips respect to the one which have seen the energy deposit
1764 //
1765 if ( jj >= 0 && jj < 96 ){
1766 if ( !ctground ){
1767 if ( jj%16 == 0 ) jjpre++;
1768 if ( jj != 0 && jj != 32 && jj != 64 && ene[jj-1] != 0. ) ene[jj-1] += -clevel1->estrip[jj][m][l] * ctneigcor[l][m][jjpre];
1769 if ( jj != 31 && jj != 63 && jj != 95 && ene[jj+1] != 0. ) ene[jj+1] += -clevel1->estrip[jj][m][l] * ctneigcor[l][m][jjpre];
1770 } else {
1771 if ( jj != 0 && jj != 32 && jj != 64 && ene[jj-1] != 0. ) ene[jj-1] += -clevel1->estrip[jj][m][l] * 0.01581;
1772 if ( jj != 31 && jj != 63 && jj != 95 && ene[jj+1] != 0. ) ene[jj+1] += -clevel1->estrip[jj][m][l] * 0.01581;
1773 };
1774 };
1775 if ( jjj >= 0 && jjj < 96 ){
1776 if ( !ctground ){
1777 if ( jjj%16 == 0 ) jjjpre++;
1778 if ( jjj != 0 && jjj != 32 && jjj != 64 && clevel1->estrip[jjj-1][m][l] != 0. ) clevel1->estrip[jjj-1][m][l] += -ene[jjj] * ctneigcor[l][m][jjjpre];
1779 if ( jjj != 31 && jjj != 63 && jjj != 95 && clevel1->estrip[jjj+1][m][l] !=0. ) clevel1->estrip[jjj+1][m][l] += -ene[jjj] * ctneigcor[l][m][jjjpre];
1780 } else {
1781 if ( jjj != 0 && jjj != 32 && jjj != 64 && clevel1->estrip[jjj-1][m][l] != 0. ) clevel1->estrip[jjj-1][m][l] += -ene[jjj] * 0.01581;
1782 if ( jjj != 31 && jjj != 63 && jjj != 95 && clevel1->estrip[jjj+1][m][l] != 0. ) clevel1->estrip[jjj+1][m][l] += -ene[jjj] * 0.01581;
1783 };
1784 };
1785 };
1786 if ( j4 >= 0 && j4 < 96 ){
1787 //
1788 // CALOLEVEL1 CODING AND FILLING
1789 //
1790 //
1791 // NOTICE: THE FOLLOWING LINE EXCLUDE ALL STRIPS FOR WHICH THE RMS*4 IS GREATER THAN 26 !!! <=============== IMPORTANT! =================> // not true anymore, now it trust parameter files
1792 //
1793 if ( obadmask[l][m][j4] == 1 || clevel1->estrip[j4][m][l] <= clevel1->emin || clevel1->estrip[j4][m][l] <= memin[l][m][j4] || calrms[l][m][j4] > maxrms[l][m] || (l==0 && m == 18 && mask18 ) ){
1794 clevel1->estrip[j4][m][l] = 0.;
1795 };
1796 //
1797 if ( debug ) printf(" STRIP: view %i plane %i strip %i energy: %f \n",l,m,j4,clevel1->estrip[j4][m][l]);
1798 //
1799 // code and save the energy for each strip in svstrip
1800 //
1801 if ( clevel1->estrip[j4][m][l] > clevel1->emin ){
1802 //
1803 Float_t savel1 = clevel1->estrip[j4][m][l];
1804 //
1805 if ( m == 18 && l == 0 ){
1806 if ( debug ) printf(" Resetting plane 18X for variable calculation: view %i plane %i strip %i \n",l,m,j4);
1807 clevel1->estrip[j4][m][l] = 0.; // SAVE STRIPS VALUE FOR PLANE 18 X but DO NOT USE IT FOR VARIABLE CALCULATION
1808 };
1809 if ( debug ) printf(" HIT STRIP: view %i plane %i strip %i energy: %f \n",l,m,j4,clevel1->estrip[j4][m][l]);
1810 // if ( dexyc[l][m][j4] == 32767. ){
1811 if ( dexyc[l][m][j4] > 32000. || savel1 > 5000.){ // CaloLevel1 bug with plane 18X [9RED 14/04/2010]
1812 if ( savel1 > 5000 ){
1813 if ( debug ) printf(" Absurd plane 18X energy... resetting value to 1100 MIP \n");
1814 savel1 = 1100.; // CaloLevel1 bug with plane 18x [9RED 14/04/2010]
1815 };
1816 savel1 += 5000.;
1817 clevel2->nsatstrip += 1.;
1818 };
1819 //
1820 tim = 100000.;
1821 plo = m;
1822 fbi = 0;
1823 if ( savel1 > 0.99995 ){
1824 tim = 10000.;
1825 plo = m;
1826 fbi = 1;
1827 };
1828 if ( savel1 > 9.9995 ){
1829 tim = 1000.;
1830 plo = 22 + m;
1831 fbi = 1;
1832 };
1833 if ( savel1 > 99.995 ){
1834 tim = 100.;
1835 plo = 22 + m;
1836 fbi = 0;
1837 };
1838 if ( savel1 > 999.95 ){
1839 tim = 10.;
1840 plo = 44 + m;
1841 fbi = 0;
1842 };
1843 if ( savel1 > 9999.5 ){
1844 tim = 1.;
1845 plo = 66 + m;
1846 fbi = 0;
1847 };
1848 //
1849 cle = (Int_t)lroundf(tim*savel1);
1850 //
1851 if ( l == 0 ){
1852 //
1853 // +-PPSSmmmm.mmmm
1854 //
1855 svstrip[istrip] = fbi*1000000000 + plo*10000000 + j4*100000 + cle;
1856 } else {
1857 svstrip[istrip] = -(fbi*1000000000 + plo*10000000 + j4*100000 + cle);
1858 };
1859 if ( debug ) printf(" svstrip[%i] = %i fbi %i plo %i j4 %i cle %i \n",istrip,svstrip[istrip],fbi,plo,j4,cle);
1860 //
1861 istrip++;
1862 };
1863 };
1864 };
1865 //
1866 };
1867 };
1868 //
1869 // store goodness flag
1870 //
1871 if ( !pe ){
1872 clevel2->good = 1;
1873 } else {
1874 clevel2->good = 0;
1875 };
1876 //
1877 // done
1878 //
1879 return(0);
1880 }
1881
1882 void CaloLevel0::GetTrkVar(){
1883 calol2tr();
1884 }
1885
1886 void CaloLevel0::FillTrkVar(CaloLevel2 *ca, Int_t nutrk){
1887 //
1888 CaloTrkVar *t_ca = new CaloTrkVar();
1889 //
1890 t_ca->trkseqno = trkseqno;
1891 t_ca->ncore = (Int_t)clevel2->ncore;
1892 t_ca->qcore = clevel2->qcore;
1893 t_ca->noint = (Int_t)clevel2->noint;
1894 t_ca->ncyl = (Int_t)clevel2->ncyl;
1895 t_ca->qcyl = clevel2->qcyl;
1896 t_ca->qtrack = clevel2->qtrack;
1897 t_ca->qtrackx = clevel2->qtrackx;
1898 t_ca->qtracky = clevel2->qtracky;
1899 t_ca->dxtrack = clevel2->dxtrack;
1900 t_ca->dytrack = clevel2->dytrack;
1901 t_ca->qlast = clevel2->qlast;
1902 t_ca->nlast = (Int_t)clevel2->nlast;
1903 t_ca->qpre = clevel2->qpre;
1904 t_ca->npre = (Int_t)clevel2->npre;
1905 t_ca->qpresh = clevel2->qpresh;
1906 t_ca->npresh = (Int_t)clevel2->npresh;
1907 t_ca->qtr = clevel2->qtr;
1908 t_ca->ntr = (Int_t)clevel2->ntr;
1909 t_ca->planetot = (Int_t)clevel2->planetot;
1910 t_ca->qmean = clevel2->qmean;
1911 t_ca->dX0l = clevel2->dX0l;
1912 t_ca->qlow = clevel2->qlow;
1913 t_ca->nlow = (Int_t)clevel2->nlow;
1914 //
1915 if ( trkseqno == -1 ){
1916 // ca->impx = clevel2->impx;
1917 // ca->impy = clevel2->impy;
1918 ca->tanx[1] = clevel2->tanx;
1919 ca->tany[1] = clevel2->tany;
1920 ca->elen = clevel2->elen;
1921 ca->selen = clevel2->selen;
1922 // memcpy(ca->cibar,clevel2->cibar,sizeof(clevel2->cibar));
1923 // memcpy(ca->cbar,clevel2->cbar,sizeof(clevel2->cbar));
1924 memcpy(t_ca->tibar,clevel2->cibar,sizeof(clevel2->cibar));
1925 memcpy(t_ca->tbar,clevel2->cbar,sizeof(clevel2->cbar));
1926 memcpy(ca->planemax,clevel2->planemax,sizeof(clevel2->planemax));
1927 memcpy(ca->selfdelay,clevel2->selfdelay,sizeof(clevel2->selfdelay));
1928 ca->varcfit[2] = clevel2->varcfit[0];
1929 ca->varcfit[3] = clevel2->varcfit[1];
1930 ca->npcfit[2] = clevel2->npcfit[0];
1931 ca->npcfit[3] = clevel2->npcfit[1];
1932 // memcpy(ca->varcfit,clevel2->varcfit,sizeof(clevel2->varcfit));
1933 // memcpy(ca->npcfit,clevel2->npcfit,sizeof(clevel2->npcfit));
1934 } else {
1935 memcpy(t_ca->tibar,clevel2->tibar,sizeof(clevel2->tibar));
1936 memcpy(t_ca->tbar,clevel2->tbar,sizeof(clevel2->tbar));
1937 };
1938 //
1939 //
1940 if ( !(ca->CaloTrk) ) ca->CaloTrk = new TClonesArray("CaloTrkVar",1); //ELENA
1941 TClonesArray &t = *ca->CaloTrk;
1942 new(t[nutrk]) CaloTrkVar(*t_ca);
1943 //
1944 delete t_ca;
1945 //
1946 ClearTrkVar();
1947 }
1948
1949 void CaloLevel0::FillTrkVar(TClonesArray *tcl, Int_t nutrk){
1950 //
1951 CaloTrkVar *t_ca = new CaloTrkVar();
1952 //
1953 t_ca->trkseqno = trkseqno;
1954 t_ca->ncore = (Int_t)clevel2->ncore;
1955 t_ca->qcore = clevel2->qcore;
1956 t_ca->noint = (Int_t)clevel2->noint;
1957 t_ca->ncyl = (Int_t)clevel2->ncyl;
1958 t_ca->qcyl = clevel2->qcyl;
1959 t_ca->qtrack = clevel2->qtrack;
1960 t_ca->qtrackx = clevel2->qtrackx;
1961 t_ca->qtracky = clevel2->qtracky;
1962 t_ca->dxtrack = clevel2->dxtrack;
1963 t_ca->dytrack = clevel2->dytrack;
1964 t_ca->qlast = clevel2->qlast;
1965 t_ca->nlast = (Int_t)clevel2->nlast;
1966 t_ca->qpre = clevel2->qpre;
1967 t_ca->npre = (Int_t)clevel2->npre;
1968 t_ca->qpresh = clevel2->qpresh;
1969 t_ca->npresh = (Int_t)clevel2->npresh;
1970 t_ca->qtr = clevel2->qtr;
1971 t_ca->ntr = (Int_t)clevel2->ntr;
1972 t_ca->planetot = (Int_t)clevel2->planetot;
1973 t_ca->qmean = clevel2->qmean;
1974 t_ca->dX0l = clevel2->dX0l;
1975 t_ca->qlow = clevel2->qlow;
1976 t_ca->nlow = (Int_t)clevel2->nlow;
1977 //
1978 memcpy(t_ca->tibar,clevel2->tibar,sizeof(clevel2->tibar));
1979 memcpy(t_ca->tbar,clevel2->tbar,sizeof(clevel2->tbar));
1980 //
1981 //
1982 TClonesArray &t = *tcl;
1983 new(t[nutrk]) CaloTrkVar(*t_ca);
1984 //
1985 delete t_ca;
1986 //
1987 ClearTrkVar();
1988 }
1989
1990 void CaloLevel0::GetCommonVar(){
1991 calol2cm();
1992 }
1993
1994 void CaloLevel0::FillCommonVar(CaloLevel1 *c1, CaloLevel2 *ca){
1995 //
1996 ca->good = clevel2->good;
1997 // if ( clevel2->trigty == 2. ){
1998 // ca->selftrigger = 1;
1999 // } else {
2000 // ca->selftrigger = 0;
2001 // };
2002 //
2003 ca->selftrigger = (Int_t)clevel2->trigty + (Int_t)clevel2->wartrig;
2004 //
2005 memcpy(ca->perr,clevel2->perr,sizeof(clevel2->perr));
2006 memcpy(ca->swerr,clevel2->swerr,sizeof(clevel2->swerr));
2007 memcpy(ca->crc,clevel2->crc,sizeof(clevel2->crc));
2008 ca->nstrip = (Int_t)clevel2->nstrip;
2009 ca->nsatstrip = (Int_t)clevel2->nsatstrip;
2010 ca->qtot = clevel2->qtot;
2011 // ca->impx = clevel2->impx;
2012 // ca->impy = clevel2->impy;
2013 ca->tanx[0] = clevel2->tanx;
2014 ca->tany[0] = clevel2->tany;
2015 ca->nx22 = (Int_t)clevel2->nx22;
2016 ca->qx22 = clevel2->qx22;
2017 ca->qmax = clevel2->qmax;
2018 ca->elen = clevel2->elen;
2019 ca->selen = clevel2->selen;
2020 memcpy(ca->qq,clevel2->qq,sizeof(clevel2->qq));
2021 memcpy(ca->planemax,clevel2->planemax,sizeof(clevel2->planemax));
2022 memcpy(ca->selfdelay,clevel2->selfdelay,sizeof(clevel2->selfdelay));
2023 ca->varcfit[0] = clevel2->varcfit[0];
2024 ca->varcfit[1] = clevel2->varcfit[1];
2025 ca->npcfit[0] = clevel2->npcfit[0];
2026 ca->npcfit[1] = clevel2->npcfit[1];
2027 ca->fitmode[0] = clevel2->fmode[0];
2028 ca->fitmode[1] = clevel2->fmode[1];
2029 // memcpy(ca->varcfit,clevel2->varcfit,sizeof(clevel2->varcfit));
2030 // memcpy(ca->npcfit,clevel2->npcfit,sizeof(clevel2->npcfit));
2031 memcpy(ca->cibar,clevel2->cibar,sizeof(clevel2->cibar));
2032 memcpy(ca->cbar,clevel2->cbar,sizeof(clevel2->cbar));
2033 //
2034 if ( c1 ){
2035 c1->istrip = istrip;
2036 c1->estrip = TArrayI(istrip,svstrip);
2037 };
2038 //
2039 }
2040 void CaloLevel0::FillCommonVar(CaloLevel1 *c1){
2041 if ( c1 ){
2042 c1->istrip = istrip;
2043 c1->estrip = TArrayI(istrip,svstrip);
2044 };
2045 //
2046 }
2047
2048 void CaloLevel0::ClearStructs(){
2049 ClearTrkVar();
2050 ClearCommonVar();
2051 }
2052
2053 void CaloLevel0::Delete(Option_t *t){
2054 if ( de ) delete de;
2055 delete this;
2056 }
2057
2058
2059 void CaloLevel0::RunClose(){
2060 l0tr->Delete();
2061 ClearStructs();
2062 //
2063 memset(dexy, 0, 2*22*96*sizeof(Float_t));
2064 memset(dexyc, 0, 2*22*96*sizeof(Float_t));
2065 memset(base, 0, 2*22*6*sizeof(Float_t));
2066 memset(sbase, 0, 2*22*6*sizeof(Float_t));
2067 memset(ctprecor, 0, 2*22*6*sizeof(Float_t));
2068 memset(ctsicor, 0, 2*22*9*sizeof(Float_t));
2069 memset(ctneigcor, 0, 2*22*6*sizeof(Float_t));
2070 //
2071 }
2072
2073 //
2074 // Private methods
2075 //
2076
2077 void CaloLevel0::ClearTrkVar(){
2078 clevel2->ncore = 0;
2079 clevel2->qcore = 0.;
2080 clevel2->noint = 0.;
2081 clevel2->ncyl = 0.;
2082 clevel2->qcyl = 0.;
2083 clevel2->qtrack = 0.;
2084 clevel2->qtrackx = 0.;
2085 clevel2->qtracky = 0.;
2086 clevel2->dxtrack = 0.;
2087 clevel2->dytrack = 0.;
2088 clevel2->qlast = 0.;
2089 clevel2->nlast = 0.;
2090 clevel2->qpre = 0.;
2091 clevel2->npre = 0.;
2092 clevel2->qpresh = 0.;
2093 clevel2->npresh = 0.;
2094 clevel2->qlow = 0.;
2095 clevel2->nlow = 0.;
2096 clevel2->qtr = 0.;
2097 clevel2->ntr = 0.;
2098 clevel2->planetot = 0.;
2099 clevel2->qmean = 0.;
2100 clevel2->dX0l = 0.;
2101 clevel2->elen = 0.;
2102 clevel2->selen = 0.;
2103 memset(clevel1->al_p, 0, 5*2*sizeof(Double_t));
2104 memset(clevel2->tibar, 0, 2*22*sizeof(Int_t));
2105 memset(clevel2->tbar, 0, 2*22*sizeof(Float_t));
2106 }
2107
2108 void CaloLevel0::ClearCommonVar(){
2109 istrip = 0;
2110 clevel2->trigty = -1.;
2111 clevel2->wartrig = 0.;
2112 clevel2->good = 0;
2113 clevel2->nstrip = 0.;
2114 clevel2->nsatstrip = 0.;
2115 clevel2->qtot = 0.;
2116 // clevel2->impx = 0.;
2117 // clevel2->impy = 0.;
2118 clevel2->tanx = 0.; // this is correct since it refers to the fortran structure
2119 clevel2->tany = 0.; // this is correct since it refers to the fortran structure
2120 clevel2->qmax = 0.;
2121 clevel2->nx22 = 0.;
2122 clevel2->qx22 = 0.;
2123 memset(clevel2->perr, 0, 4*sizeof(Int_t));
2124 memset(clevel2->swerr, 0, 4*sizeof(Int_t));
2125 memset(clevel2->crc, 0, 4*sizeof(Int_t));
2126 memset(clevel2->qq, 0, 4*sizeof(Int_t));
2127 memset(clevel2->varcfit, 0, 4*sizeof(Float_t));
2128 memset(clevel2->npcfit, 0, 4*sizeof(Int_t));
2129 memset(clevel2->planemax, 0, 2*sizeof(Int_t));
2130 memset(clevel2->selfdelay, 0, 4*7*sizeof(Int_t));
2131 memset(clevel2->fmode, 0, 2*sizeof(Int_t));
2132 memset(clevel2->cibar, 0, 2*22*sizeof(Int_t));
2133 memset(clevel2->cbar, 0, 2*22*sizeof(Float_t));
2134 }
2135
2136 void CaloLevel0::ClearCalibVals(Int_t s){
2137 //
2138 for ( Int_t d=0 ; d<11 ;d++ ){
2139 Int_t pre = -1;
2140 for ( Int_t j=0; j<96 ;j++){
2141 if ( j%16 == 0 ) pre++;
2142 if ( s == 2 ){
2143 calped[0][2*d+1][j] = 0.;
2144 cstwerr[3] = 0.;
2145 cperror[3] = 0.;
2146 calgood[0][2*d+1][j] = 0.;
2147 calthr[0][2*d+1][pre] = 0.;
2148 calrms[0][2*d+1][j] = 0.;
2149 calbase[0][2*d+1][pre] = 0.;
2150 calvar[0][2*d+1][pre] = 0.;
2151 };
2152 if ( s == 3 ){
2153 calped[0][2*d][j] = 0.;
2154 cstwerr[1] = 0.;
2155 cperror[1] = 0.;
2156 calgood[0][2*d][j] = 0.;
2157 calthr[0][2*d][pre] = 0.;
2158 calrms[0][2*d][j] = 0.;
2159 calbase[0][2*d][pre] = 0.;
2160 calvar[0][2*d][pre] = 0.;
2161 };
2162 if ( s == 0 ){
2163 calped[1][2*d][j] = 0.;
2164 cstwerr[0] = 0.;
2165 cperror[0] = 0.;
2166 calgood[1][2*d][j] = 0.;
2167 calthr[1][2*d][pre] = 0.;
2168 calrms[1][2*d][j] = 0.;
2169 calbase[1][2*d][pre] = 0.;
2170 calvar[1][2*d][pre] = 0.;
2171 };
2172 if ( s == 1 ){
2173 calped[1][2*d+1][j] = 0.;
2174 cstwerr[2] = 0.;
2175 cperror[2] = 0.;
2176 calgood[1][2*d+1][j] = 0.;
2177 calthr[1][2*d+1][pre] = 0.;
2178 calrms[1][2*d+1][j] = 0.;
2179 calbase[1][2*d+1][pre] = 0.;
2180 calvar[1][2*d+1][pre] = 0.;
2181 };
2182 };
2183 };
2184 return;
2185 }
2186
2187 Int_t CaloLevel0::Update(GL_TABLES *glt, UInt_t atime, Int_t s){
2188 //
2189 const TString host = glt->CGetHost();
2190 const TString user = glt->CGetUser();
2191 const TString psw = glt->CGetPsw();
2192 TSQLServer *dbc = TSQLServer::Connect(host.Data(),user.Data(),psw.Data());
2193 if ( !dbc->IsConnected() ) throw -116;
2194 stringstream myquery;
2195 myquery.str("");
2196 myquery << "SET time_zone='+0:00';";
2197 delete dbc->Query(myquery.str().c_str());
2198 delete dbc->Query("SET sql_mode = 'NO_UNSIGNED_SUBTRACTION';");
2199 Int_t sgnl = 0;
2200 //
2201 GL_CALO_CALIB *glcalo = new GL_CALO_CALIB();
2202 //
2203 sgnl = 0;
2204 //
2205 idcalib[s] = 0;
2206 fromtime[s] = 0;
2207 totime[s] = 0;
2208 calibno[s] = 0;
2209 ClearCalibVals(s);
2210 //
2211 UInt_t uptime = 0;
2212 //
2213 sgnl = glcalo->Query_GL_CALO_CALIB(atime,uptime,s,dbc);
2214 if ( sgnl < 0 ){
2215 if ( verbose ) printf(" CALORIMETER - ERROR: error from GLTables\n");
2216 return(sgnl);
2217 };
2218 //
2219 idcalib[s] = glcalo->ID_ROOT_L0;
2220 fromtime[s] = glcalo->FROM_TIME;
2221 if ( glcalo->TO_TIME < atime ){ // calibration is corrupted and we are using the one that preceed the good one
2222 totime[s] = uptime;
2223 } else {
2224 totime[s] = glcalo->TO_TIME;
2225 };
2226 // totime[s] = glcalo->TO_TIME;
2227 calibno[s] = glcalo->EV_ROOT;
2228 //
2229 if ( totime[s] == 0 ){
2230 if ( verbose ) printf(" CALORIMETER - WARNING: data with no associated calibration\n");
2231 ClearCalibVals(s);
2232 sgnl = 100;
2233 };
2234 //
2235 // determine path and name and entry of the calibration file
2236 //
2237 GL_ROOT *glroot = new GL_ROOT();
2238 if ( verbose ) printf("\n");
2239 if ( verbose ) printf(" ** SECTION %i **\n",s);
2240 //
2241 sgnl = glroot->Query_GL_ROOT(idcalib[s],dbc);
2242 if ( sgnl < 0 ){
2243 if ( verbose ) printf(" CALORIMETER - ERROR: error from GLTables\n");
2244 return(sgnl);
2245 };
2246 //
2247 stringstream name;
2248 name.str("");
2249 name << glroot->PATH.Data() << "/";
2250 name << glroot->NAME.Data();
2251 //
2252 fcalname[s] = (TString)name.str().c_str();
2253 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]);
2254 //
2255 sgnl = LoadCalib(s);
2256 //
2257 if ( sgnl != 0 ) return(sgnl);
2258 delete glcalo;
2259 delete glroot;
2260 //
2261 return(0);
2262 //
2263 }
2264
2265 Int_t CaloLevel0::LoadCalib(Int_t s){
2266 //
2267 ifstream myfile;
2268 myfile.open(fcalname[s].Data());
2269 if ( !myfile ){
2270 return(-107);
2271 };
2272 myfile.close();
2273 //
2274 TFile *File = new TFile(fcalname[s].Data());
2275 if ( !File ) return(-108);
2276 TTree *tr = (TTree*)File->Get("CalibCalPed");
2277 if ( !tr ) return(-109);
2278 //
2279 TBranch *calo = tr->GetBranch("CalibCalPed");
2280 //
2281 pamela::CalibCalPedEvent *ce = 0;
2282 tr->SetBranchAddress("CalibCalPed", &ce);
2283 //
2284 Long64_t ncalibs = calo->GetEntries();
2285 //
2286 if ( !ncalibs ) return(-110);
2287 //
2288 if ( calo->GetEntry(calibno[s]) <= 0 ) throw -36;
2289 //
2290 if (ce->cstwerr[s] != 0 && ce->cperror[s] == 0 ) {
2291 for ( Int_t d=0 ; d<11 ;d++ ){
2292 Int_t pre = -1;
2293 for ( Int_t j=0; j<96 ;j++){
2294 if ( j%16 == 0 ) pre++;
2295 if ( s == 2 ){
2296 calped[0][2*d+1][j] = ce->calped[3][d][j];
2297 cstwerr[3] = ce->cstwerr[3];
2298 cperror[3] = ce->cperror[3];
2299 calgood[0][2*d+1][j] = ce->calgood[3][d][j];
2300 calthr[0][2*d+1][pre] = ce->calthr[3][d][pre];
2301 calrms[0][2*d+1][j] = ce->calrms[3][d][j];
2302 calbase[0][2*d+1][pre] = ce->calbase[3][d][pre];
2303 calvar[0][2*d+1][pre] = ce->calvar[3][d][pre];
2304 };
2305 if ( s == 3 ){
2306 calped[0][2*d][j] = ce->calped[1][d][j];
2307 cstwerr[1] = ce->cstwerr[1];
2308 cperror[1] = ce->cperror[1];
2309 calgood[0][2*d][j] = ce->calgood[1][d][j];
2310 calthr[0][2*d][pre] = ce->calthr[1][d][pre];
2311 calrms[0][2*d][j] = ce->calrms[1][d][j];
2312 calbase[0][2*d][pre] = ce->calbase[1][d][pre];
2313 calvar[0][2*d][pre] = ce->calvar[1][d][pre];
2314 };
2315 if ( s == 0 ){
2316 calped[1][2*d][j] = ce->calped[0][d][j];
2317 cstwerr[0] = ce->cstwerr[0];
2318 cperror[0] = ce->cperror[0];
2319 calgood[1][2*d][j] = ce->calgood[0][d][j];
2320 calthr[1][2*d][pre] = ce->calthr[0][d][pre];
2321 calrms[1][2*d][j] = ce->calrms[0][d][j];
2322 calbase[1][2*d][pre] = ce->calbase[0][d][pre];
2323 calvar[1][2*d][pre] = ce->calvar[0][d][pre];
2324 };
2325 if ( s == 1 ){
2326 calped[1][2*d+1][j] = ce->calped[2][d][j];
2327 cstwerr[2] = ce->cstwerr[2];
2328 cperror[2] = ce->cperror[2];
2329 calgood[1][2*d+1][j] = ce->calgood[2][d][j];
2330 calthr[1][2*d+1][pre] = ce->calthr[2][d][pre];
2331 calrms[1][2*d+1][j] = ce->calrms[2][d][j];
2332 calbase[1][2*d+1][pre] = ce->calbase[2][d][pre];
2333 calvar[1][2*d+1][pre] = ce->calvar[2][d][pre];
2334 };
2335 };
2336 };
2337 } else {
2338 if ( verbose ) printf(" CALORIMETER - ERROR: problems finding a good calibration in this file! \n\n ");
2339 return(-111);
2340 };
2341 File->Close();
2342 return(0);
2343 }

  ViewVC Help
Powered by ViewVC 1.1.23