/[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.22 - (show annotations) (download)
Thu Jun 19 20:05:44 2008 UTC (16 years, 5 months ago) by mocchiut
Branch: MAIN
Changes since 1.21: +891 -392 lines
New calorimeter calibration

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

  ViewVC Help
Powered by ViewVC 1.1.23