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

  ViewVC Help
Powered by ViewVC 1.1.23