/[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.30 - (show annotations) (download)
Fri Jan 29 05:49:24 2010 UTC (14 years, 10 months ago) by mocchiut
Branch: MAIN
CVS Tags: v9r00, v9r01
Changes since 1.29: +4 -4 lines
Check I/O problems while getting entry

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

  ViewVC Help
Powered by ViewVC 1.1.23