/[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.21 - (show annotations) (download)
Mon May 12 14:36:08 2008 UTC (16 years, 6 months ago) by mocchiut
Branch: MAIN
Changes since 1.20: +130 -64 lines
New method in CaloStrip + cross-talk bugs fixed

1 /**
2 * \file src/CaloLevel0.cpp
3 * \author Emiliano Mocchiutti
4 **/
5 //
6 // C/C++ headers
7 //
8 #include <sstream>
9 #include <fstream>
10 //
11 // ROOT headers
12 //
13 #include <TTree.h>
14 #include <TBranch.h>
15 #include <TFile.h>
16 #include <TObject.h>
17 //
18 // YODA headers
19 //
20 #include <PamelaRun.h>
21 #include <physics/calorimeter/CalorimeterEvent.h>
22 #include <CalibCalPedEvent.h>
23 //
24 //
25 //
26 #include <GLTables.h>
27 //
28 // this package headers
29 //
30 #include <delay.h>
31 #include <CaloLevel0.h>
32 //
33 //
34 // Declaration of the core fortran routines
35 //
36 #define calol2cm calol2cm_
37 extern "C" int calol2cm();
38 #define calol2tr calol2tr_
39 extern "C" int calol2tr();
40 //
41 using namespace std;
42 //
43 //
44 // Public methods
45 //
46
47 CaloLevel0::~CaloLevel0(){
48 if ( de ) delete de;
49 delete this;
50 }
51
52 CaloLevel0::CaloLevel0(){
53 //
54 extern struct FlCaLevel1 clevel1_;
55 extern struct FlCaLevel2 clevel2_;
56 clevel1 = &clevel1_;
57 clevel2 = &clevel2_;
58 //
59 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(ctneigcor, 0, 2*22*6*sizeof(Float_t));
71 calopar1 = true;
72 calopar2 = true;
73 calopar3 = true;
74 crosst = true;
75 ftcalopar1 = 0;
76 ttcalopar1 = 0;
77 ftcalopar2 = 0;
78 ttcalopar2 = 0;
79 ftcalopar3 = 0;
80 ttcalopar3 = 0;
81 }
82
83 void CaloLevel0::SetCrossTalk(Bool_t ct){
84 crosst = ct;
85 }
86
87 void CaloLevel0::SetCrossTalkType(Bool_t ct){
88 ctground = ct;
89 }
90
91 void CaloLevel0::SetVerbose(Bool_t ct){
92 verbose = ct;
93 }
94
95 /**
96 * Initialize CaloLevel0 object
97 **/
98 void CaloLevel0::ProcessingInit(TSQLServer *dbc, UInt_t hs, Int_t &sgnl, TTree *l0tree, Bool_t isdeb, Bool_t isverb){
99 if ( !dbc->IsConnected() ) throw -116;
100 this->InitDo(dbc,hs,sgnl,l0tree,isdeb,isverb);
101 }
102
103 /**
104 * Initialize CaloLevel0 object
105 **/
106 void CaloLevel0::ProcessingInit(GL_TABLES *glt, UInt_t hs, Int_t &sgnl, TTree *l0tree, Bool_t isdeb, Bool_t isverb){
107 //
108 const TString host = glt->CGetHost();
109 const TString user = glt->CGetUser();
110 const TString psw = glt->CGetPsw();
111 TSQLServer *dbc = TSQLServer::Connect(host.Data(),user.Data(),psw.Data());
112 if ( !dbc->IsConnected() ) throw -116;
113 this->InitDo(dbc,hs,sgnl,l0tree,isdeb,isverb);
114 dbc->Close();
115 delete dbc;
116 }
117
118
119 void CaloLevel0::InitDo(TSQLServer *dbc, UInt_t hs, Int_t &sgnl, TTree *l0tree, Bool_t isdeb, Bool_t isverb){
120 stringstream myquery;
121 myquery.str("");
122 myquery << "SET time_zone='+0:00'";
123 dbc->Query(myquery.str().c_str());
124 //
125 debug = isdeb;
126 verbose = isverb;
127 //
128 l0tr=(TTree*)l0tree;
129 de = new pamela::calorimeter::CalorimeterEvent();
130 l0calo = (TBranch*)l0tr->GetBranch("Calorimeter");
131 l0tr->SetBranchAddress("Calorimeter", &de);
132 //
133 trkseqno = 0;
134 ClearStructs();
135 //
136 GL_CALO_CALIB *glcalo = new GL_CALO_CALIB();
137 //
138 sgnl = 0;
139 UInt_t uptime = 0;
140 //
141 for (Int_t s = 0; s < 4; s++){
142 idcalib[s] = 0;
143 fromtime[s] = 0;
144 totime[s] = 0;
145 calibno[s] = 0;
146 ClearCalibVals(s);
147 //
148 sgnl = glcalo->Query_GL_CALO_CALIB(hs,uptime,s,dbc);
149 if ( sgnl < 0 ){
150 if ( verbose ) printf(" CALORIMETER - ERROR: error from GLTables\n");
151 return;
152 };
153 //
154 idcalib[s] = glcalo->ID_ROOT_L0;
155 fromtime[s] = glcalo->FROM_TIME;
156 if ( glcalo->TO_TIME < hs ){ // calibration is corrupted and we are using the one that preceed the good one
157 totime[s] = uptime;
158 } else {
159 totime[s] = glcalo->TO_TIME;
160 };
161 calibno[s] = glcalo->EV_ROOT;
162 //
163 if ( totime[s] == 0 ){
164 if ( verbose ) printf(" CALORIMETER - WARNING: data with no associated calibration\n");
165 ClearCalibVals(s);
166 sgnl = 100;
167 };
168 };
169 //
170 // determine path and name and entry of the calibration file
171 //
172 GL_ROOT *glroot = new GL_ROOT();
173 if ( verbose ) printf("\n");
174 for (Int_t s = 0; s < 4; s++){
175 if ( verbose ) printf(" ** SECTION %i **\n",s);
176 if ( totime[s] > 0 ){
177 //
178 sgnl = glroot->Query_GL_ROOT(idcalib[s],dbc);
179 if ( sgnl < 0 ){
180 if ( verbose ) printf(" CALORIMETER - ERROR: error from GLTables\n");
181 return;
182 };
183 //
184 stringstream name;
185 name.str("");
186 name << glroot->PATH.Data() << "/";
187 name << glroot->NAME.Data();
188 //
189 fcalname[s] = (TString)name.str().c_str();
190 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]);
191 } else {
192 if ( verbose ) printf(" - runheader at time %u. NO CALIBRATION INCLUDE THE RUNHEADER! ",hs);
193 };
194 sgnl = LoadCalib(s);
195 if ( sgnl ) break;
196 };
197 //
198 delete glcalo;
199 delete glroot;
200 //
201 return;
202 //
203 }
204
205 Int_t CaloLevel0::ChkCalib(GL_TABLES *glt, UInt_t atime){
206 Int_t sgnl = 0;
207 for ( Int_t s = 0; s < 4; s++){
208 if ( atime > totime[s] ){
209 sgnl = Update(glt,atime,s);
210 if ( sgnl < 0 ) return(sgnl);
211 };
212 };
213 return(sgnl);
214 }
215
216 Int_t CaloLevel0::ChkParam(TSQLServer *dbc, UInt_t runheader, Bool_t mechal){
217 Int_t sig = this->ChkParamDo(dbc,runheader,mechal);
218 return(sig);
219 }
220
221 Int_t CaloLevel0::ChkParam(GL_TABLES *glt, UInt_t runheader, Bool_t mechal){
222 const TString host = glt->CGetHost();
223 const TString user = glt->CGetUser();
224 const TString psw = glt->CGetPsw();
225 TSQLServer *dbc = TSQLServer::Connect(host.Data(),user.Data(),psw.Data());
226 if ( !dbc->IsConnected() ) throw -116;
227 stringstream myquery;
228 myquery.str("");
229 myquery << "SET time_zone='+0:00'";
230 dbc->Query(myquery.str().c_str());
231 //
232 Int_t sig = this->ChkParamDo(dbc,runheader,mechal);
233 dbc->Close();
234 delete dbc;
235 return(sig);
236 }
237
238 Int_t CaloLevel0::ChkParamDo(TSQLServer *dbc, UInt_t runheader, Bool_t mechal){
239 //
240 stringstream calfile;
241 stringstream bmfile;
242 stringstream aligfile;
243 Int_t error = 0;
244 FILE *f = 0;
245 ifstream badfile;
246 GL_PARAM *glparam = new GL_PARAM();
247 //
248 if ( calopar1 || ( ttcalopar1 != 0 && ttcalopar1 < runheader ) ){
249 //
250 if ( debug ) printf(" calopar1 %i ftcalopar1 %u ttcalopar1 %u runheader %u \n",calopar1,ftcalopar1,ttcalopar1,runheader);
251 //
252 calopar1 = false;
253 //
254 // determine where I can find calorimeter ADC to MIP conversion file
255 //
256 if ( verbose ) printf(" Querying DB for calorimeter parameters files...\n");
257 //
258 error = 0;
259 error = glparam->Query_GL_PARAM(runheader,101,dbc);
260 if ( error < 0 ) return(error);
261 //
262 calfile.str("");
263 calfile << glparam->PATH.Data() << "/";
264 calfile << glparam->NAME.Data();
265 ftcalopar1 = glparam->FROM_TIME;
266 ttcalopar1 = glparam->TO_TIME;
267 //
268 if ( verbose ) printf("\n Using ADC to MIP conversion file: \n %s \n",calfile.str().c_str());
269 f = fopen(calfile.str().c_str(),"rb");
270 if ( !f ){
271 if ( verbose ) printf(" CALORIMETER - ERROR: no ADC to MIP file!\n");
272 return(-105);
273 };
274 //
275 for (Int_t m = 0; m < 2 ; m++ ){
276 for (Int_t k = 0; k < 22; k++ ){
277 for (Int_t l = 0; l < 96; l++ ){
278 fread(&mip[m][k][l],sizeof(mip[m][k][l]),1,f);
279 if ( debug ) printf(" %f \n",mip[m][k][l]);
280 };
281 };
282 };
283 fclose(f);
284 };
285 //
286 if ( calopar2 || ( ttcalopar2 != 0 && ttcalopar2 < runheader ) ){
287 //
288 if ( debug ) printf(" calopar2 %i ftcalopar2 %u ttcalopar2 %u runheader %u \n",calopar2,ftcalopar2,ttcalopar2,runheader);
289 calopar2 = false;
290 //
291 // determine where I can find calorimeter alignment file
292 //
293 //
294 error = 0;
295 error = glparam->Query_GL_PARAM(runheader,102,dbc);
296 if ( error < 0 ) return(error);
297 //
298 aligfile.str("");
299 aligfile << glparam->PATH.Data() << "/";
300 aligfile << glparam->NAME.Data();
301 ftcalopar2 = glparam->FROM_TIME;
302 ttcalopar2 = glparam->TO_TIME;
303 //
304 if ( verbose ) printf("\n Using parameter file: \n %s \n",aligfile.str().c_str());
305 f = fopen(aligfile.str().c_str(),"rb");
306 if ( !f ){
307 if ( verbose ) printf(" CALORIMETER - ERROR: no parameter file!\n");
308 return(-106);
309 };
310 //
311 if ( !mechal ){
312 //
313 fread(&clevel1->xalig,sizeof(clevel1->xalig),1,f);
314 if ( debug ) printf(" xalig = %f \n",clevel1->xalig);
315 fread(&clevel1->yalig,sizeof(clevel1->yalig),1,f);
316 if ( debug ) printf(" yalig = %f \n",clevel1->yalig);
317 fread(&clevel1->zalig,sizeof(clevel1->zalig),1,f);
318 if ( debug ) printf(" zalig = %f \n",clevel1->zalig);
319 } else {
320 if ( verbose ) printf("\n Using MECHANICAL alignement parameters \n");
321 //
322 CaloStrip cs = CaloStrip();
323 cs.UseMechanicalAlig();
324 clevel1->xalig = cs.GetXalig();
325 if ( debug ) printf(" xalig = %f \n",clevel1->xalig);
326 clevel1->yalig = cs.GetYalig();
327 if ( debug ) printf(" yalig = %f \n",clevel1->yalig);
328 clevel1->zalig = cs.GetZalig();
329 if ( debug ) printf(" zalig = %f \n",clevel1->zalig);
330 //
331 Float_t tmp = 0;
332 fread(&tmp,sizeof(clevel1->xalig),1,f);
333 fread(&tmp,sizeof(clevel1->yalig),1,f);
334 fread(&tmp,sizeof(clevel1->zalig),1,f);
335 // clevel1->zalig = -265.82;
336 //
337 };
338 fread(&clevel1->emin,sizeof(clevel1->emin),1,f);
339 if ( debug ) printf(" signal threshold = %f \n",clevel1->emin);
340 //
341 fclose(f);
342 };
343 //
344 // Load offline bad strip mask
345 //
346 if ( calopar3 || ( ttcalopar3 != 0 && ttcalopar3 < runheader ) ){
347 if ( debug ) printf(" calopar3 %i ftcalopar3 %u ttcalopar3 %u runheader %u \n",calopar3,ftcalopar3,ttcalopar3,runheader);
348 calopar3 = false;
349 //
350 // determine where I can find calorimeter alignment file
351 //
352 //
353 error = 0;
354 error = glparam->Query_GL_PARAM(runheader,103,dbc);
355 if ( error < 0 ) return(error);
356 //
357 bmfile.str("");
358 bmfile << glparam->PATH.Data() << "/";
359 bmfile << glparam->NAME.Data();
360 ftcalopar3 = glparam->FROM_TIME;
361 ttcalopar3 = glparam->TO_TIME;
362 //
363 if ( verbose ) printf("\n Using bad strip offline mask file: \n %s \n\n",bmfile.str().c_str());
364 badfile.open(bmfile.str().c_str());
365 if ( !badfile ){
366 if ( verbose ) printf(" CALORIMETER - ERROR: no bad strip offline mask file!\n");
367 return(-115);
368 };
369 //
370 Bool_t isdone = false;
371 Int_t bad = 0;
372 Int_t view = 1;
373 Int_t strip = 0;
374 Int_t plane = 21;
375 while ( !isdone ) {
376 badfile >> bad;
377 obadmask[view][plane][strip] = bad;
378 if ( debug && bad ) printf(" SETTING view %i plane %i strip %i BAD = %i \n",view,plane,strip,bad);
379 strip++;
380 if ( strip > 95 ){
381 strip = 0;
382 plane--;
383 if ( plane < 0 ){
384 plane = 21;
385 view--;
386 };
387 if ( view < 0 ) isdone = true;
388 };
389 };
390 //
391 badfile.close();
392 };
393 //
394 delete glparam;
395 //
396 return(0);
397 }
398
399 Int_t CaloLevel0::CalcCrossTalkCorr(TSQLServer *dbc, UInt_t runheader){
400 Int_t sig = CalcCrossTalkCorrDo(dbc,runheader);
401 return(sig);
402 }
403
404 Int_t CaloLevel0::CalcCrossTalkCorr(GL_TABLES *glt, UInt_t runheader){
405 const TString host = glt->CGetHost();
406 const TString user = glt->CGetUser();
407 const TString psw = glt->CGetPsw();
408 TSQLServer *dbc = TSQLServer::Connect(host.Data(),user.Data(),psw.Data());
409 if ( !dbc->IsConnected() ) throw -116;
410 stringstream myquery;
411 myquery.str("");
412 myquery << "SET time_zone='+0:00'";
413 dbc->Query(myquery.str().c_str());
414 //
415 Int_t sig = CalcCrossTalkCorrDo(dbc,runheader);
416 dbc->Close();
417 delete dbc;
418 //
419 return(sig);
420 //
421 }
422
423 Int_t CaloLevel0::CalcCrossTalkCorrDo(TSQLServer *dbc, UInt_t runheader){
424 //
425 if ( ctground ) return(0);
426 //
427 stringstream bmfile;
428 Int_t error = 0;
429 ifstream badfile;
430 GL_PARAM *glparam = new GL_PARAM();
431 //
432 // determine where I can find file with offline bad pulser mask
433 //
434 error = 0;
435 error = glparam->Query_GL_PARAM(runheader,105,dbc);
436 if ( error < 0 ) return(error);
437 //
438 bmfile.str("");
439 bmfile << glparam->PATH.Data() << "/";
440 bmfile << glparam->NAME.Data();
441 //
442 if ( verbose ) printf("\n Using bad pulser 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 pulser 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 pre = 0;
453 Int_t plane = 21;
454 while ( !isdone ) {
455 badfile >> bad;
456 obadpulsemask[view][plane][pre] = bad;
457 if ( debug && bad ) printf(" SETTING view %i plane %i pre %i BAD = %i \n",view,plane,pre,bad);
458 pre++;
459 if ( pre > 5 ){
460 pre = 0;
461 plane--;
462 if ( plane < 0 ){
463 plane = 21;
464 view--;
465 };
466 if ( view < 0 ) isdone = true;
467 };
468 };
469 //
470 delete glparam;
471 badfile.close();
472 //
473 // Let's start with cross-talk correction calculation
474 //
475 GL_CALOPULSE_CALIB *glp = new GL_CALOPULSE_CALIB();
476 Float_t adcp[2][22][96];
477 Float_t adcpcal[2][22][96];
478 memset(adcp , 0, 2*22*96*sizeof(Float_t));
479 memset(adcpcal , 0, 2*22*96*sizeof(Float_t));
480 //
481 UInt_t pampli = 0;
482 for (Int_t s=0; s<4; s++){
483 //
484 // Save into matrix adcp the values of the highest pulse calibration (pulse amplitude = 2)
485 //
486 pampli = 2;
487 error = 0;
488 error = glp->Query_GL_CALOPULSE_CALIB(runheader,s,pampli,dbc);
489 if ( error < 0 ){
490 if ( verbose ) printf(" CALORIMETER - ERROR: error from GLTables\n");
491 return(error);
492 };
493 //
494 UInt_t idcalib = glp->ID_ROOT_L0;
495 UInt_t fromtime = glp->FROM_TIME;
496 UInt_t calibno = glp->EV_ROOT;
497 //
498 // determine path and name and entry of the calibration file
499 //
500 GL_ROOT *glroot = new GL_ROOT();
501 if ( verbose ) printf("\n");
502 if ( verbose ) printf(" ** SECTION %i **\n",s);
503 //
504 error = 0;
505 error = glroot->Query_GL_ROOT(idcalib,dbc);
506 if ( error < 0 ){
507 if ( verbose ) printf(" CALORIMETER - ERROR: error from GLTables\n");
508 return(error);
509 };
510 //
511 stringstream name;
512 name.str("");
513 name << glroot->PATH.Data() << "/";
514 name << glroot->NAME.Data();
515 //
516 TString fcalname = (TString)name.str().c_str();
517 ifstream myfile;
518 myfile.open(fcalname.Data());
519 if ( !myfile ){
520 return(-107);
521 };
522 myfile.close();
523 //
524 TFile *File = new TFile(fcalname.Data());
525 if ( !File ) return(-108);
526 TTree *tr = (TTree*)File->Get("CalibCalPulse2");
527 if ( !tr ) return(-119);
528 //
529 TBranch *calo = tr->GetBranch("CalibCalPulse2");
530 //
531 pamela::CalibCalPulse2Event *ce = 0;
532 tr->SetBranchAddress("CalibCalPulse2", &ce);
533 //
534 Long64_t ncalibs = calo->GetEntries();
535 //
536 if ( !ncalibs ) return(-110);
537 //
538 calo->GetEntry(calibno);
539 if ( verbose ) printf(" PULSE2 using entry %u from file %s",calibno,fcalname.Data());
540 //
541 // retrieve calibration table
542 //
543 if ( ce->pstwerr[s] && ce->pperror[s] == 0 && ce->unpackError == 0 ){
544 for ( Int_t d=0 ; d<11 ;d++ ){
545 for ( Int_t j=0; j<96 ;j++){
546 if ( s == 2 ){
547 adcp[0][2*d+1][j] = ce->calpuls[3][d][j];
548 };
549 if ( s == 3 ){
550 adcp[0][2*d][j] = ce->calpuls[1][d][j];
551 };
552 if ( s == 0 ){
553 adcp[1][2*d][j] = ce->calpuls[0][d][j];
554 };
555 if ( s == 1 ){
556 adcp[1][2*d+1][j] = ce->calpuls[2][d][j];
557 };
558 };
559 };
560 } else {
561 if ( verbose ) printf(" CALORIMETER - ERROR: problems finding a good calibration in this file! \n\n ");
562 return(-111);
563 };
564 //
565 File->Close();
566 delete glroot;
567 //
568 // Save into matrix adcpcal the calibrated values of the pulse calibration (subtraction of pulse amplitude = 0 relative to the pulse2 calibration used)
569 //
570 pampli = 0;
571 error = 0;
572 error = glp->Query_GL_CALOPULSE_CALIB(fromtime,s,pampli,dbc);
573 if ( error < 0 ){
574 if ( verbose ) printf(" CALORIMETER - ERROR: error from GLTables\n");
575 return(error);
576 };
577 //
578 idcalib = glp->ID_ROOT_L0;
579 calibno = glp->EV_ROOT;
580 //
581 // determine path and name and entry of the calibration file
582 //
583 glroot = new GL_ROOT();
584 if ( verbose ) printf("\n");
585 if ( verbose ) printf(" ** SECTION %i **\n",s);
586 //
587 error = 0;
588 error = glroot->Query_GL_ROOT(idcalib,dbc);
589 if ( error < 0 ){
590 if ( verbose ) printf(" CALORIMETER - ERROR: error from GLTables\n");
591 return(error);
592 };
593 //
594 name.str("");
595 name << glroot->PATH.Data() << "/";
596 name << glroot->NAME.Data();
597 //
598 fcalname = (TString)name.str().c_str();
599 myfile.open(fcalname.Data());
600 if ( !myfile ){
601 return(-107);
602 };
603 myfile.close();
604 //
605 TFile *File1 = new TFile(fcalname.Data());
606 if ( !File1 ) return(-108);
607 TTree *tr1 = (TTree*)File1->Get("CalibCalPulse1");
608 if ( !tr1 ) return(-120);
609 //
610 TBranch *calo1 = tr1->GetBranch("CalibCalPulse1");
611 //
612 pamela::CalibCalPulse1Event *ce1 = 0;
613 tr1->SetBranchAddress("CalibCalPulse1", &ce1);
614 //
615 ncalibs = calo1->GetEntries();
616 //
617 if ( !ncalibs ) return(-110);
618 //
619 calo1->GetEntry(calibno);
620 if ( verbose ) printf(" PULSE1 using entry %u from file %s",calibno,fcalname.Data());
621 //
622 // retrieve calibration table
623 //
624 if ( ce1->pstwerr[s] && ce1->pperror[s] == 0 && ce1->unpackError == 0 ){
625 for ( Int_t d=0 ; d<11 ;d++ ){
626 for ( Int_t j=0; j<96 ;j++){
627 if ( s == 2 ){
628 adcpcal[0][2*d+1][j] = adcp[0][2*d+1][j] - ce1->calpuls[3][d][j];
629 };
630 if ( s == 3 ){
631 adcpcal[0][2*d][j] = adcp[0][2*d][j] - ce1->calpuls[1][d][j];
632 };
633 if ( s == 0 ){
634 adcpcal[1][2*d][j] = adcp[1][2*d][j] - ce1->calpuls[0][d][j];
635 };
636 if ( s == 1 ){
637 adcpcal[1][2*d+1][j] = adcp[1][2*d+1][j] - ce1->calpuls[2][d][j];
638 };
639 };
640 };
641 } else {
642 if ( verbose ) printf(" CALORIMETER - ERROR: problems finding a good calibration in this file! \n\n ");
643 return(-111);
644 };
645 //
646 File1->Close();
647 //
648 delete glroot;
649 //
650 };// loop on the four sections
651 //
652 //
653 delete glp;
654 //
655 // Ok, now we can try to calculate the cross-talk correction for each pre-amplifier
656 //
657 for ( Int_t v=0; v<2; v++){
658 if ( debug ) printf(" \n\n NEW VIEW \n");
659 for ( Int_t p=0; p<22; p++){
660 for ( Int_t npre=0; npre<6; npre++){
661 ctprecor[v][p][npre] = 1000.;
662 ctneigcor[v][p][npre] = 1000.;
663 Int_t str0=npre*16;
664 Int_t str16= -1 + (1+npre)*16;
665 //
666 UInt_t neigc = 0;
667 UInt_t farc = 0;
668 UInt_t pulsc = 0;
669 Float_t sigpulsed = 0.;
670 Float_t neigbase = 0.;
671 Float_t farbase = 0.;
672 //
673 // 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
674 // moreover count the number of strips in each case
675 //
676 for (Int_t s=str0; s<=str16; s++){
677 if ( adcpcal[v][p][s] > 10000.){
678 sigpulsed = adcpcal[v][p][s];
679 pulsc++;
680 if ( s > str0 ){
681 neigbase += adcpcal[v][p][s-1];
682 neigc++;
683 farbase -= adcpcal[v][p][s-1];
684 farc--;
685 };
686 if ( s < str16 ){
687 neigbase += adcpcal[v][p][s+1];
688 neigc++;
689 farbase -= adcpcal[v][p][s+1];
690 farc--;
691 };
692 } else {
693 farc++;
694 farbase += adcpcal[v][p][s];
695 };
696 };
697 //
698 // Now calculate the corrections
699 //
700 Float_t avefarbase = 0.;
701 if ( farc ) avefarbase = farbase/(Float_t)farc;
702 Float_t aveneigbase = 0.;
703 if ( neigc ) aveneigbase = neigbase/(Float_t)neigc;
704 //
705 if ( pulsc == 1 && farc && neigc ){
706 ctprecor[v][p][npre] = -avefarbase/(sigpulsed+fabs(avefarbase));
707 ctneigcor[v][p][npre] = fabs(aveneigbase-avefarbase)/(sigpulsed+fabs(avefarbase));
708 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]);
709 } else {
710 //
711 // did not find the pulsed strip or more than one pulsed strip found!
712 //
713 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);
714 //
715 };
716 };
717 if ( debug ) printf(" \n ==================== \n");
718 };
719 };
720 //
721 // Check the calculated corrections
722 //
723 Int_t opre=0;
724 Int_t ppre=0;
725 Bool_t found = false;
726 for ( Int_t v=0; v<2; v++){
727 for ( Int_t p=0; p<22; p++){
728 for ( Int_t npre=0; npre<6; npre++){
729 if ( ctprecor[v][p][npre] == 1000. || ctneigcor[v][p][npre] == 1000. || obadpulsemask[v][p][npre] != 0 ){
730 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]);
731 if ( npre%2 ){
732 opre = npre-1;
733 } else {
734 opre = npre+1;
735 };
736 if ( ctprecor[v][p][opre] == 1000. || ctneigcor[v][p][opre] == 1000. || obadpulsemask[v][p][opre] != 0 ){
737 ppre=0;
738 found = false;
739 while ( ppre < 6 ){
740 if ( ctprecor[v][p][ppre] != 1000. && ctneigcor[v][p][ppre] != 1000. && !obadpulsemask[v][p][ppre] ){
741 found = true;
742 ctprecor[v][p][npre] = ctprecor[v][p][ppre];
743 ctneigcor[v][p][npre] = ctneigcor[v][p][ppre];
744 break;
745 };
746 ppre++;
747 };
748 if ( !found ){
749 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);
750 ctprecor[v][p][npre] = 0.002;
751 ctneigcor[v][p][npre] = 0.002;
752 };
753 } else {
754 ctprecor[v][p][npre] = ctprecor[v][p][opre];
755 ctneigcor[v][p][npre] = ctneigcor[v][p][opre];
756 };
757 if ( debug ) printf(" AFTER: pre-correction: %f neighbour strips correction %f \n",ctprecor[v][p][npre],ctneigcor[v][p][npre]);
758 };
759 };
760 };
761 };
762 //
763 return(0);
764 }
765
766 void CaloLevel0::FindBaseCompress(Int_t l, Int_t m, Int_t pre){
767 Int_t n = 0;
768 Float_t q = 0;
769 this->FindBaseCompress(l,m,pre,n,q);
770 }
771
772 void CaloLevel0::FindBaseCompress(Int_t l, Int_t m, Int_t pre, Int_t &nst, Float_t &qp){
773 for (Int_t e = pre*16; e < (pre+1)*16 ; e++){
774 dexy[l][m][e] = dexyc[l][m][e];
775 };
776 this->FindBaseRaw(l,m,pre,nst,qp);
777 }
778
779 void CaloLevel0::FindBaseRaw(Int_t l, Int_t m, Int_t pre){
780 Int_t n = 0;
781 Float_t q = 0;
782 this->FindBaseRaw(l,m,pre,n,q);
783 }
784
785 void CaloLevel0::FindBaseRaw(Int_t l, Int_t m, Int_t pre, Int_t &nst, Float_t &qp){
786 //
787 Float_t minstrip = 100000.;
788 Float_t rms = 0.;
789 Int_t process = 0;
790 Int_t onlmask[16];
791 memset(onlmask, 0, 16*sizeof(Int_t));
792 //
793 while ( process < 2 ){
794 //
795 minstrip = 100000.;
796 rms = 0.;
797 base[l][m][pre] = 0.;
798 qp = 0.;
799 //
800 Int_t spos = -1;
801 Int_t ee = 0;
802 for (Int_t e = pre*16; e < (pre+1)*16 ; e++){
803 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 ) {
804 minstrip = dexy[l][m][e]-calped[l][m][e];
805 rms = calthr[l][m][pre];
806 spos = ee;
807 };
808 ee++;
809 qp += (dexy[l][m][e]-calped[l][m][e]-sbase[l][m][e]);
810 };
811 //
812 if ( debug && l==0 ){
813 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);
814 };
815 if ( minstrip != 100000. ) {
816 Float_t strip6s = 0.;
817 for (Int_t e = pre*16; e < (pre+1)*16 ; e++){
818 if ( (dexy[l][m][e]-calped[l][m][e]) >= minstrip && (dexy[l][m][e]-calped[l][m][e]) <= (minstrip+rms) ) {
819 strip6s += 1.;
820 base[l][m][pre] += (dexy[l][m][e] - calped[l][m][e]);
821 };
822 //
823 // compression
824 //
825 // if ( abs((int)(dexy[l][m][e]-calped[l][m][e])) <= (minstrip+rms) ) {
826 // dexyc[l][m][e] = 0.;
827 // } else {
828 dexyc[l][m][e] = dexy[l][m][e];
829 // };
830 };
831 //
832 if ( strip6s == 1. && process < 1 ){
833 onlmask[spos] = 1;
834 process++;
835 continue;
836 };
837 process += 2;
838 nst = (Int_t)strip6s;
839 //
840 if ( debug && l==1 ){
841 printf(" strip6s %f \n",strip6s);
842 };
843 // if ( strip6s >= 9. ){
844 if ( (strip6s >= 2. && process == 2) || (strip6s >= 9. && process > 2) ){
845 Double_t arro = base[l][m][pre]/strip6s;
846 Float_t deci = 1000.*((float)arro - float(int(arro)));
847 if ( deci < 500. ) {
848 arro = double(int(arro));
849 } else {
850 arro = 1. + double(int(arro));
851 };
852 base[l][m][pre] = arro;
853 //
854 // 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
855 //
856 if ( debug && l==1 ) printf(" Calculated baseline: base %f sbase*1.02 %f \n",base[l][m][pre],1.02*sbase[l][m][pre]);
857 //
858 if ( strip6s < 4 && base[l][m][pre] > 1.02*sbase[l][m][pre] && sbase[l][m][pre] > 0. ){
859 if ( debug ) printf(" Suspicious calculated baseline: base %f sbase*1.02 %f strip6s %i \n",base[l][m][pre],1.02*sbase[l][m][pre],(Int_t)strip6s);
860 base[l][m][pre] = 31000.;
861 for (Int_t e = pre*16; e < (pre+1)*16 ; e++){
862 dexyc[l][m][e] = dexy[l][m][e];
863 };
864 };
865 } else {
866 base[l][m][pre] = 31000.;
867 for (Int_t e = pre*16; e < (pre+1)*16 ; e++){
868 dexyc[l][m][e] = dexy[l][m][e];
869 };
870 };
871 } else {
872 process += 2;
873 base[l][m][pre] = 31000.;
874 for (Int_t e = pre*16; e < (pre+1)*16 ; e++){
875 dexyc[l][m][e] = dexy[l][m][e];
876 };
877 };
878 };
879 }
880
881 Int_t CaloLevel0::Calibrate(Int_t ei){
882 //
883 // get entry ei
884 //
885 l0calo->GetEntry(ei);
886 //
887 // if it was not a selftrigger event, could it ever been a selftrigger event? if so trigty = 3.
888 //
889 clevel2->nsatstrip = 0.;
890 Int_t val = 0;
891 Int_t del = 1100;
892 for (Int_t sec = 0; sec < 4; sec++){
893 for (Int_t dsec = 0; dsec < 7; dsec++){
894 val = (Int_t)de->calselftrig[sec][dsec];
895 del = delay(val);
896 clevel2->selfdelay[sec][dsec] = del;
897 };
898 };
899 val = 0;
900 del = 1100;
901 if ( clevel2->trigty != 2. ){
902 Bool_t ck = false;
903 for (Int_t sec = 0; sec < 4; sec++){
904 val = (Int_t)de->calselftrig[sec][6];
905 del = delay(val);
906 if ( del < 1100 ){
907 clevel2->wartrig = 0.;
908 clevel2->trigty = 3.;
909 ck = true;
910 break;
911 };
912 };
913 if ( !ck ) clevel2->wartrig = 100.;
914 } else {
915 Bool_t ck = false;
916 for (Int_t sec = 0; sec < 4; sec++){
917 val = (Int_t)de->calselftrig[sec][6];
918 del = delay(val);
919 if ( del < 1100 ){
920 clevel2->wartrig = 0.;
921 ck = true;
922 };
923 };
924 if ( !ck ) clevel2->wartrig = 100.;
925 };
926 //
927 Int_t se = 5;
928 Int_t done = 0;
929 Int_t pre = -1;
930 Bool_t isCOMP = false;
931 Bool_t isFULL = false;
932 Bool_t isRAW = false;
933 Float_t ener;
934 Int_t doneb = 0;
935 Int_t donec = 0;
936 Int_t ck[6] = {0,0,0,0,0,0};
937 Int_t ipre = 0;
938 // Int_t ip[3] = {0};
939 Int_t ip[3] = {0,0,0};
940 Float_t base0, base1, base2;
941 base0 = 0.;
942 base1 = 0.;
943 base2 = 0.;
944 Float_t qpre[6] = {0.,0.,0.,0.,0.,0.};
945 Float_t ene[96];
946 Int_t chdone[4] = {0,0,0,0};
947 Int_t pe = 0;
948 //
949 Float_t ener0 = 0.;
950 Float_t cbase0 = 0.;
951 Bool_t pproblem = false;
952 //
953 Float_t tim = 0.;
954 Int_t plo = 0;
955 Int_t fbi = 0;
956 Int_t cle = 0;
957 //
958 // run over views and planes
959 //
960 for (Int_t l = 0; l < 2; l++){
961 for (Int_t m = 0; m < 22; m++){
962 //
963 // determine the section number
964 //
965 se = 5;
966 if (l == 0 && m%2 == 0) se = 3;
967 if (l == 0 && m%2 != 0) se = 2;
968 if (l == 1 && m%2 != 0) se = 1;
969 if (l == 1 && m%2 == 0) se = 0;
970 //
971 // determine what kind of event we are going to analyze
972 //
973 isCOMP = false;
974 isFULL = false;
975 isRAW = false;
976 if ( de->stwerr[se] & (1 << 16) ) isCOMP = true;
977 if ( de->stwerr[se] & (1 << 17) ) isFULL = true;
978 if ( de->stwerr[se] & (1 << 3) ) isRAW = true;
979 if ( !chdone[se] ){
980 //
981 // check for any error in the event
982 //
983 clevel2->crc[se] = 0;
984 if ( de->perror[se] == 132 ){
985 clevel2->crc[se] = 1;
986 pe++;
987 };
988 clevel2->perr[se] = 0;
989 if ( de->perror[se] != 0 ){
990 clevel2->perr[se] = 1;
991 pe++;
992 };
993 clevel2->swerr[se] = 0;
994 for (Int_t j = 0; j < 7 ; j++){
995 if ( (j != 3) && (de->stwerr[se] & (1 << j)) ){
996 clevel2->swerr[se] = 1;
997 pe++;
998 };
999 };
1000 chdone[se] = 1;
1001 };
1002 if ( clevel2->crc[se] == 0 && (clevel1->good2 == 1 || clevel2->trigty >= 2) ){
1003 pre = -1;
1004 //
1005 for (Int_t nn = 0; nn < 96; nn++){
1006 ene[nn] = 0.;
1007 dexy[l][m][nn] = de->dexy[l][m][nn] ;
1008 dexyc[l][m][nn] = de->dexyc[l][m][nn] ;
1009 };
1010 //
1011 // run over preamplifiers
1012 //
1013 pre = -1;
1014 cbase0 = 0.;
1015 Int_t nstt[2];
1016 Float_t rqp[2];
1017 for (Int_t i = 0; i < 3; i++){
1018 // nstt[0] = 0; // BUG
1019 // nstt[1] = 0; // BUG
1020 nstt[0] = 1000;
1021 nstt[1] = 1000;
1022 rqp[0] = 0.;
1023 rqp[1] = 0.;
1024 for (Int_t j = 0; j < 2; j++){
1025 pre = j + i*2;
1026 //
1027 // baseline check and calculation
1028 //
1029 if ( !isRAW ){
1030 //
1031 // if it is a compress event with fully transmitted pre try to calculate the baseline
1032 //
1033 if ( de->base[l][m][pre] != 0. && de->base[l][m][pre]<31000. ) {
1034 base[l][m][pre] = de->base[l][m][pre] ;
1035 } else {
1036 FindBaseCompress(l,m,pre,nstt[j],rqp[j]);
1037 };
1038 cbase0 += base[l][m][pre];
1039 } else {
1040 //
1041 // if it is a raw event calculate the baseline.
1042 //
1043 FindBaseRaw(l,m,pre,nstt[j],rqp[j]);
1044 cbase0 += base[l][m][pre];
1045 };
1046 };
1047 //
1048 // 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
1049 //
1050 if ( nstt[0] < 4 && nstt[1] >= 4 && nstt[0] != 1000 && nstt[1] != 1000 ) base[l][m][pre-1] = 31000.;
1051 if ( nstt[0] >= 4 && nstt[1] < 4 && nstt[0] != 1000 && nstt[1] != 1000 ) base[l][m][pre] = 31000.;
1052 // if ( nstt[0] < 4 && nstt[1] >= 4 ) base[l][m][pre-1] = 31000.; // BUG
1053 // if ( nstt[0] >= 4 && nstt[1] < 4 ) base[l][m][pre] = 31000.; // BUG
1054 // //
1055 // // 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
1056 // //
1057 // if ( nstt[0] < 4 && nstt[1] < 4 ){
1058 // if ( rqp[0] >= rqp[1] ) base[l][m][pre-1] = 31000.;
1059 // if ( rqp[0] < rqp[1] ) base[l][m][pre] = 31000.;
1060 // };
1061 };
1062 //
1063 // run over strips
1064 //
1065 pre = -1;
1066 ener0 = 0.;
1067 for (Int_t i = 0 ; i < 3 ; i++){
1068 ip[i] = 0;
1069 for (Int_t n = i*32 ; n < (i+1)*32 ; n++){
1070 if (n%16 == 0) {
1071 done = 0;
1072 doneb = 0;
1073 donec = 0;
1074 pre++;
1075 ck[pre] = 0;
1076 qpre[pre] = 0.;
1077 };
1078 //
1079 // baseline check and calculation
1080 //
1081 // no suitable new baseline, use old ones!
1082 //
1083 if ( !done ){
1084 if ( (base[l][m][pre] == 31000. || base[l][m][pre] == 0.) ){
1085 ck[pre] = 1;
1086 if (pre%2 == 0) {
1087 ip[i] = pre + 1;
1088 } else {
1089 ip[i] = pre - 1;
1090 };
1091 if ( (base[l][m][ip[i]] == 31000. || base[l][m][ip[i]] == 0. || !crosst ) ){
1092 //
1093 ck[pre] = 2;
1094 if ( sbase[l][m][pre] == 31000. || sbase[l][m][pre] == 0. ) {
1095 ck[pre] = 3;
1096 };
1097 };
1098 // done = 1;
1099 };
1100 done = 1;
1101 };
1102 //
1103 // CALIBRATION ALGORITHM
1104 //
1105 if ( !doneb ){
1106 if ( debug ) printf(" ck[pre] is %i \n",ck[pre]);
1107 switch (ck[pre]) {
1108 case 0:
1109 base0 = base[l][m][pre];
1110 base2 = calbase[l][m][pre];
1111 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]);
1112 break;
1113 case 1:
1114 base0 = base[l][m][ip[i]];
1115 base2 = calbase[l][m][ip[i]];
1116 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]]);
1117 break;
1118 case 2:
1119 base0 = sbase[l][m][pre];
1120 base2 = calbase[l][m][pre];
1121 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]);
1122 break;
1123 case 3:
1124 base0 = calbase[l][m][pre];
1125 base2 = calbase[l][m][pre];
1126 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]);
1127 break;
1128 };
1129 base1 = calbase[l][m][pre];
1130 doneb = 1;
1131 };
1132 ener = dexyc[l][m][n];
1133 ener0 += ener;
1134 clevel1->estrip[n][m][l] = 0.;
1135 if ( base0>0 && base0 < 30000. ){
1136 // if ( !donec && (base0 - base1 + base2) != 0. ){
1137 // sbase[l][m][pre] = base0 - base1 + base2;
1138 if ( !donec && (base0 + base1 - base2) != 0. ){
1139 sbase[l][m][pre] = base0 + base1 - base2;
1140 donec = 1;
1141 };
1142 if ( ener > 0. ){
1143 clevel1->estrip[n][m][l] = (ener - calped[l][m][n] - base0 - base1 + base2)/mip[l][m][n] ;
1144 //
1145 // 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)
1146 //
1147 qpre[pre] += clevel1->estrip[n][m][l];
1148 //
1149 //
1150 };
1151 };
1152 };
1153 if ( crosst ){
1154 if (ck[pre] == 1 || ck[pre-1] == 1){
1155 if (ck[pre] == 1){
1156 ipre = pre;
1157 ip[i] = pre - 1;
1158 } else {
1159 ipre = pre - 1;
1160 ip[i] = pre;
1161 };
1162 // if (ip[i]%2 == 0) {
1163 // ipre = ip[i] + 1;
1164 // } else {
1165 // ipre = ip[i] - 1;
1166 // };
1167 for (Int_t j = ipre*16 ; j < (ipre+1)*16 ; j++){
1168 if ( !ctground ){
1169 clevel1->estrip[j][m][l] += qpre[ipre] * ctprecor[l][m][ipre] - qpre[ip[i]] * ctprecor[l][m][ip[i]];
1170 } else {
1171 clevel1->estrip[j][m][l] += (qpre[ipre] - qpre[ip[i]]) * 0.00478;
1172 };
1173 };
1174 };
1175 if (ck[pre] == 2 && ck[pre-1] == 2){
1176 for (Int_t j = i*32 ; j < (i+1)*32 ; j++){
1177 // ipre = j/16 + 1;
1178 ipre = j/16 ;
1179 if ( !ctground ){
1180 clevel1->estrip[j][m][l] += qpre[ipre] * ctprecor[l][m][ipre];
1181 } else {
1182 clevel1->estrip[j][m][l] += qpre[ipre] * 0.00478;
1183 };
1184 };
1185 };
1186 };
1187 };
1188 //
1189 if ( ener0 == 0. && cbase0 == 0. && !pproblem && clevel2->perr[se] == 0){
1190 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);
1191 pproblem = true;
1192 pe++;
1193 };
1194 //
1195 Int_t j4 = -4;
1196 Int_t jjj = -3;
1197 Int_t jj = -2;
1198 Int_t jjpre = -1;
1199 Int_t jjjpre = -1;
1200 for (Int_t j = 0 ; j < 100 ; j++){
1201 jj++;
1202 jjj++;
1203 j4++;
1204 if ( j < 96 ) ene[j] = clevel1->estrip[j][m][l];
1205 if ( crosst ){
1206 if ( jj >= 0 && jj < 96 ){
1207 if ( !ctground ){
1208 if ( jj%16 == 0 ) jjpre++;
1209 if ( jj != 0 && jj != 32 && jj != 64 ) ene[jj-1] += -clevel1->estrip[jj][m][l] * ctneigcor[l][m][jjpre];
1210 if ( jj != 31 && jj != 63 && jj != 95 ) ene[jj+1] += -clevel1->estrip[jj][m][l] * ctneigcor[l][m][jjpre];
1211 } else {
1212 if ( jj != 0 && jj != 32 && jj != 64 ) ene[jj-1] += -clevel1->estrip[jj][m][l] * 0.01581;
1213 if ( jj != 31 && jj != 63 && jj != 95 ) ene[jj+1] += -clevel1->estrip[jj][m][l] * 0.01581;
1214 };
1215 };
1216 if ( jjj >= 0 && jjj < 96 ){
1217 if ( !ctground ){
1218 if ( jjj%16 == 0 ) jjjpre++;
1219 if ( jjj != 0 && jjj != 32 && jjj != 64 ) clevel1->estrip[jjj-1][m][l] += -ene[jjj] * ctneigcor[l][m][jjjpre];
1220 if ( jjj != 31 && jjj != 63 && jjj != 95 ) clevel1->estrip[jjj+1][m][l] += -ene[jjj] * ctneigcor[l][m][jjjpre];
1221 } else {
1222 if ( jjj != 0 && jjj != 32 && jjj != 64 ) clevel1->estrip[jjj-1][m][l] += -ene[jjj] * 0.01581;
1223 if ( jjj != 31 && jjj != 63 && jjj != 95 ) clevel1->estrip[jjj+1][m][l] += -ene[jjj] * 0.01581;
1224 };
1225 };
1226 };
1227 if ( j4 >= 0 && j4 < 96 ){
1228 //
1229 // NOTICE: THE FOLLOWING LINE EXCLUDE ALL STRIPS FOR WHICH THE RMS*4 IS GREATER THAN 26 !!! <=============== IMPORTANT! =================>
1230 //
1231 if ( obadmask[l][m][j4] == 1 || clevel1->estrip[j4][m][l] <= clevel1->emin || calrms[l][m][j4] > 26 ){
1232 clevel1->estrip[j4][m][l] = 0.;
1233 };
1234 //
1235 // code and save the energy for each strip in svstrip
1236 //
1237 if ( clevel1->estrip[j4][m][l] > clevel1->emin ){
1238 //
1239 Float_t savel1 = clevel1->estrip[j4][m][l];
1240 if ( dexyc[l][m][j4] == 32767. ){
1241 savel1 += 5000.;
1242 clevel2->nsatstrip += 1.;
1243 };
1244 //
1245 tim = 100000.;
1246 plo = m;
1247 fbi = 0;
1248 if ( savel1 > 0.99995 ){
1249 tim = 10000.;
1250 plo = m;
1251 fbi = 1;
1252 };
1253 if ( savel1 > 9.9995 ){
1254 tim = 1000.;
1255 plo = 22 + m;
1256 fbi = 1;
1257 };
1258 if ( savel1 > 99.995 ){
1259 tim = 100.;
1260 plo = 22 + m;
1261 fbi = 0;
1262 };
1263 if ( savel1 > 999.95 ){
1264 tim = 10.;
1265 plo = 44 + m;
1266 fbi = 0;
1267 };
1268 if ( savel1 > 9999.5 ){
1269 tim = 1.;
1270 plo = 66 + m;
1271 fbi = 0;
1272 };
1273 //
1274 cle = (Int_t)lroundf(tim*savel1);
1275 //
1276 if ( l == 0 ){
1277 //
1278 // +-PPSSmmmm.mmmm
1279 //
1280 svstrip[istrip] = fbi*1000000000 + plo*10000000 + j4*100000 + cle;
1281 } else {
1282 svstrip[istrip] = -(fbi*1000000000 + plo*10000000 + j4*100000 + cle);
1283 };
1284 //
1285 // if ( ei >= -770 ) printf(" j %i l %i m %i estrip %f \n",j4,l,m,clevel1->estrip[j4][m][l]);
1286 // if ( ei >= -770 ) printf(" num lim %i fbi %i tim %f plo %i cle %i \n",numeric_limits<Int_t>::max(),fbi,tim,plo,cle);
1287 // if ( ei >= -770 ) printf(" svstrip %i \n",svstrip[istrip]);
1288 //
1289 istrip++;
1290 };
1291 };
1292 };
1293 //
1294 } else {
1295 for (Int_t nn = 0; nn < 96; nn++){
1296 clevel1->estrip[nn][m][l] = 0.;
1297 };
1298 };
1299 };
1300 };
1301 if ( !pe ){
1302 clevel2->good = 1;
1303 } else {
1304 clevel2->good = 0;
1305 };
1306 return(0);
1307 }
1308
1309 void CaloLevel0::GetTrkVar(){
1310 calol2tr();
1311 }
1312
1313 void CaloLevel0::FillTrkVar(CaloLevel2 *ca, Int_t nutrk){
1314 //
1315 CaloTrkVar *t_ca = new CaloTrkVar();
1316 //
1317 t_ca->trkseqno = trkseqno;
1318 t_ca->ncore = (Int_t)clevel2->ncore;
1319 t_ca->qcore = clevel2->qcore;
1320 t_ca->noint = (Int_t)clevel2->noint;
1321 t_ca->ncyl = (Int_t)clevel2->ncyl;
1322 t_ca->qcyl = clevel2->qcyl;
1323 t_ca->qtrack = clevel2->qtrack;
1324 t_ca->qtrackx = clevel2->qtrackx;
1325 t_ca->qtracky = clevel2->qtracky;
1326 t_ca->dxtrack = clevel2->dxtrack;
1327 t_ca->dytrack = clevel2->dytrack;
1328 t_ca->qlast = clevel2->qlast;
1329 t_ca->nlast = (Int_t)clevel2->nlast;
1330 t_ca->qpre = clevel2->qpre;
1331 t_ca->npre = (Int_t)clevel2->npre;
1332 t_ca->qpresh = clevel2->qpresh;
1333 t_ca->npresh = (Int_t)clevel2->npresh;
1334 t_ca->qtr = clevel2->qtr;
1335 t_ca->ntr = (Int_t)clevel2->ntr;
1336 t_ca->planetot = (Int_t)clevel2->planetot;
1337 t_ca->qmean = clevel2->qmean;
1338 t_ca->dX0l = clevel2->dX0l;
1339 t_ca->qlow = clevel2->qlow;
1340 t_ca->nlow = (Int_t)clevel2->nlow;
1341 //
1342 if ( trkseqno == -1 ){
1343 // ca->impx = clevel2->impx;
1344 // ca->impy = clevel2->impy;
1345 ca->tanx[1] = clevel2->tanx;
1346 ca->tany[1] = clevel2->tany;
1347 ca->elen = clevel2->elen;
1348 ca->selen = clevel2->selen;
1349 // memcpy(ca->cibar,clevel2->cibar,sizeof(clevel2->cibar));
1350 // memcpy(ca->cbar,clevel2->cbar,sizeof(clevel2->cbar));
1351 memcpy(t_ca->tibar,clevel2->cibar,sizeof(clevel2->cibar));
1352 memcpy(t_ca->tbar,clevel2->cbar,sizeof(clevel2->cbar));
1353 memcpy(ca->planemax,clevel2->planemax,sizeof(clevel2->planemax));
1354 memcpy(ca->selfdelay,clevel2->selfdelay,sizeof(clevel2->selfdelay));
1355 ca->varcfit[2] = clevel2->varcfit[0];
1356 ca->varcfit[3] = clevel2->varcfit[1];
1357 ca->npcfit[2] = clevel2->npcfit[0];
1358 ca->npcfit[3] = clevel2->npcfit[1];
1359 // memcpy(ca->varcfit,clevel2->varcfit,sizeof(clevel2->varcfit));
1360 // memcpy(ca->npcfit,clevel2->npcfit,sizeof(clevel2->npcfit));
1361 } else {
1362 memcpy(t_ca->tibar,clevel2->tibar,sizeof(clevel2->tibar));
1363 memcpy(t_ca->tbar,clevel2->tbar,sizeof(clevel2->tbar));
1364 };
1365 //
1366 //
1367 if ( !(ca->CaloTrk) ) ca->CaloTrk = new TClonesArray("CaloTrkVar",1); //ELENA
1368 TClonesArray &t = *ca->CaloTrk;
1369 new(t[nutrk]) CaloTrkVar(*t_ca);
1370 //
1371 delete t_ca;
1372 //
1373 ClearTrkVar();
1374 }
1375
1376 void CaloLevel0::GetCommonVar(){
1377 calol2cm();
1378 }
1379
1380 void CaloLevel0::FillCommonVar(CaloLevel1 *c1, CaloLevel2 *ca){
1381 //
1382 ca->good = clevel2->good;
1383 if ( clevel2->trigty == 2. ){
1384 ca->selftrigger = 1;
1385 } else {
1386 ca->selftrigger = 0;
1387 };
1388 //
1389 ca->selftrigger += (Int_t)clevel2->wartrig;
1390 //
1391 memcpy(ca->perr,clevel2->perr,sizeof(clevel2->perr));
1392 memcpy(ca->swerr,clevel2->swerr,sizeof(clevel2->swerr));
1393 memcpy(ca->crc,clevel2->crc,sizeof(clevel2->crc));
1394 ca->nstrip = (Int_t)clevel2->nstrip;
1395 ca->nsatstrip = (Int_t)clevel2->nsatstrip;
1396 ca->qtot = clevel2->qtot;
1397 // ca->impx = clevel2->impx;
1398 // ca->impy = clevel2->impy;
1399 ca->tanx[0] = clevel2->tanx;
1400 ca->tany[0] = clevel2->tany;
1401 ca->nx22 = (Int_t)clevel2->nx22;
1402 ca->qx22 = clevel2->qx22;
1403 ca->qmax = clevel2->qmax;
1404 ca->elen = clevel2->elen;
1405 ca->selen = clevel2->selen;
1406 memcpy(ca->qq,clevel2->qq,sizeof(clevel2->qq));
1407 memcpy(ca->planemax,clevel2->planemax,sizeof(clevel2->planemax));
1408 memcpy(ca->selfdelay,clevel2->selfdelay,sizeof(clevel2->selfdelay));
1409 ca->varcfit[0] = clevel2->varcfit[0];
1410 ca->varcfit[1] = clevel2->varcfit[1];
1411 ca->npcfit[0] = clevel2->npcfit[0];
1412 ca->npcfit[1] = clevel2->npcfit[1];
1413 ca->fitmode[0] = clevel2->fmode[0];
1414 ca->fitmode[1] = clevel2->fmode[1];
1415 // memcpy(ca->varcfit,clevel2->varcfit,sizeof(clevel2->varcfit));
1416 // memcpy(ca->npcfit,clevel2->npcfit,sizeof(clevel2->npcfit));
1417 memcpy(ca->cibar,clevel2->cibar,sizeof(clevel2->cibar));
1418 memcpy(ca->cbar,clevel2->cbar,sizeof(clevel2->cbar));
1419 //
1420 if ( c1 ){
1421 c1->istrip = istrip;
1422 c1->estrip = TArrayI(istrip,svstrip);
1423 };
1424 //
1425 }
1426
1427 void CaloLevel0::ClearStructs(){
1428 ClearTrkVar();
1429 ClearCommonVar();
1430 }
1431
1432 void CaloLevel0::RunClose(){
1433 l0tr->Delete();
1434 ClearStructs();
1435 //
1436 memset(dexy, 0, 2*22*96*sizeof(Float_t));
1437 memset(dexyc, 0, 2*22*96*sizeof(Float_t));
1438 memset(base, 0, 2*22*6*sizeof(Float_t));
1439 memset(sbase, 0, 2*22*6*sizeof(Float_t));
1440 memset(ctprecor, 0, 2*22*6*sizeof(Float_t));
1441 memset(ctneigcor, 0, 2*22*6*sizeof(Float_t));
1442 //
1443 }
1444
1445 //
1446 // Private methods
1447 //
1448
1449 void CaloLevel0::ClearTrkVar(){
1450 clevel2->ncore = 0;
1451 clevel2->qcore = 0.;
1452 clevel2->noint = 0.;
1453 clevel2->ncyl = 0.;
1454 clevel2->qcyl = 0.;
1455 clevel2->qtrack = 0.;
1456 clevel2->qtrackx = 0.;
1457 clevel2->qtracky = 0.;
1458 clevel2->dxtrack = 0.;
1459 clevel2->dytrack = 0.;
1460 clevel2->qlast = 0.;
1461 clevel2->nlast = 0.;
1462 clevel2->qpre = 0.;
1463 clevel2->npre = 0.;
1464 clevel2->qpresh = 0.;
1465 clevel2->npresh = 0.;
1466 clevel2->qlow = 0.;
1467 clevel2->nlow = 0.;
1468 clevel2->qtr = 0.;
1469 clevel2->ntr = 0.;
1470 clevel2->planetot = 0.;
1471 clevel2->qmean = 0.;
1472 clevel2->dX0l = 0.;
1473 clevel2->elen = 0.;
1474 clevel2->selen = 0.;
1475 memset(clevel1->al_p, 0, 5*2*sizeof(Double_t));
1476 memset(clevel2->tibar, 0, 2*22*sizeof(Int_t));
1477 memset(clevel2->tbar, 0, 2*22*sizeof(Float_t));
1478 }
1479
1480 void CaloLevel0::ClearCommonVar(){
1481 istrip = 0;
1482 clevel2->trigty = -1.;
1483 clevel2->wartrig = 0.;
1484 clevel2->good = 0;
1485 clevel2->nstrip = 0.;
1486 clevel2->nsatstrip = 0.;
1487 clevel2->qtot = 0.;
1488 // clevel2->impx = 0.;
1489 // clevel2->impy = 0.;
1490 clevel2->tanx = 0.; // this is correct since it refers to the fortran structure
1491 clevel2->tany = 0.; // this is correct since it refers to the fortran structure
1492 clevel2->qmax = 0.;
1493 clevel2->nx22 = 0.;
1494 clevel2->qx22 = 0.;
1495 memset(clevel2->perr, 0, 4*sizeof(Int_t));
1496 memset(clevel2->swerr, 0, 4*sizeof(Int_t));
1497 memset(clevel2->crc, 0, 4*sizeof(Int_t));
1498 memset(clevel2->qq, 0, 4*sizeof(Int_t));
1499 memset(clevel2->varcfit, 0, 4*sizeof(Float_t));
1500 memset(clevel2->npcfit, 0, 4*sizeof(Int_t));
1501 memset(clevel2->planemax, 0, 2*sizeof(Int_t));
1502 memset(clevel2->selfdelay, 0, 4*7*sizeof(Int_t));
1503 memset(clevel2->fmode, 0, 2*sizeof(Int_t));
1504 memset(clevel2->cibar, 0, 2*22*sizeof(Int_t));
1505 memset(clevel2->cbar, 0, 2*22*sizeof(Float_t));
1506 }
1507
1508 void CaloLevel0::ClearCalibVals(Int_t s){
1509 //
1510 for ( Int_t d=0 ; d<11 ;d++ ){
1511 Int_t pre = -1;
1512 for ( Int_t j=0; j<96 ;j++){
1513 if ( j%16 == 0 ) pre++;
1514 if ( s == 2 ){
1515 calped[0][2*d+1][j] = 0.;
1516 cstwerr[3] = 0.;
1517 cperror[3] = 0.;
1518 calgood[0][2*d+1][j] = 0.;
1519 calthr[0][2*d+1][pre] = 0.;
1520 calrms[0][2*d+1][j] = 0.;
1521 calbase[0][2*d+1][pre] = 0.;
1522 calvar[0][2*d+1][pre] = 0.;
1523 };
1524 if ( s == 3 ){
1525 calped[0][2*d][j] = 0.;
1526 cstwerr[1] = 0.;
1527 cperror[1] = 0.;
1528 calgood[0][2*d][j] = 0.;
1529 calthr[0][2*d][pre] = 0.;
1530 calrms[0][2*d][j] = 0.;
1531 calbase[0][2*d][pre] = 0.;
1532 calvar[0][2*d][pre] = 0.;
1533 };
1534 if ( s == 0 ){
1535 calped[1][2*d][j] = 0.;
1536 cstwerr[0] = 0.;
1537 cperror[0] = 0.;
1538 calgood[1][2*d][j] = 0.;
1539 calthr[1][2*d][pre] = 0.;
1540 calrms[1][2*d][j] = 0.;
1541 calbase[1][2*d][pre] = 0.;
1542 calvar[1][2*d][pre] = 0.;
1543 };
1544 if ( s == 1 ){
1545 calped[1][2*d+1][j] = 0.;
1546 cstwerr[2] = 0.;
1547 cperror[2] = 0.;
1548 calgood[1][2*d+1][j] = 0.;
1549 calthr[1][2*d+1][pre] = 0.;
1550 calrms[1][2*d+1][j] = 0.;
1551 calbase[1][2*d+1][pre] = 0.;
1552 calvar[1][2*d+1][pre] = 0.;
1553 };
1554 };
1555 };
1556 return;
1557 }
1558
1559 Int_t CaloLevel0::Update(GL_TABLES *glt, UInt_t atime, Int_t s){
1560 //
1561 const TString host = glt->CGetHost();
1562 const TString user = glt->CGetUser();
1563 const TString psw = glt->CGetPsw();
1564 TSQLServer *dbc = TSQLServer::Connect(host.Data(),user.Data(),psw.Data());
1565 if ( !dbc->IsConnected() ) throw -116;
1566 stringstream myquery;
1567 myquery.str("");
1568 myquery << "SET time_zone='+0:00'";
1569 dbc->Query(myquery.str().c_str());
1570 Int_t sgnl = 0;
1571 //
1572 GL_CALO_CALIB *glcalo = new GL_CALO_CALIB();
1573 //
1574 sgnl = 0;
1575 //
1576 idcalib[s] = 0;
1577 fromtime[s] = 0;
1578 totime[s] = 0;
1579 calibno[s] = 0;
1580 ClearCalibVals(s);
1581 //
1582 UInt_t uptime = 0;
1583 //
1584 sgnl = glcalo->Query_GL_CALO_CALIB(atime,uptime,s,dbc);
1585 if ( sgnl < 0 ){
1586 if ( verbose ) printf(" CALORIMETER - ERROR: error from GLTables\n");
1587 return(sgnl);
1588 };
1589 //
1590 idcalib[s] = glcalo->ID_ROOT_L0;
1591 fromtime[s] = glcalo->FROM_TIME;
1592 if ( glcalo->TO_TIME < atime ){ // calibration is corrupted and we are using the one that preceed the good one
1593 totime[s] = uptime;
1594 } else {
1595 totime[s] = glcalo->TO_TIME;
1596 };
1597 // totime[s] = glcalo->TO_TIME;
1598 calibno[s] = glcalo->EV_ROOT;
1599 //
1600 if ( totime[s] == 0 ){
1601 if ( verbose ) printf(" CALORIMETER - WARNING: data with no associated calibration\n");
1602 ClearCalibVals(s);
1603 sgnl = 100;
1604 };
1605 //
1606 // determine path and name and entry of the calibration file
1607 //
1608 GL_ROOT *glroot = new GL_ROOT();
1609 if ( verbose ) printf("\n");
1610 if ( verbose ) printf(" ** SECTION %i **\n",s);
1611 //
1612 sgnl = glroot->Query_GL_ROOT(idcalib[s],dbc);
1613 if ( sgnl < 0 ){
1614 if ( verbose ) printf(" CALORIMETER - ERROR: error from GLTables\n");
1615 return(sgnl);
1616 };
1617 //
1618 stringstream name;
1619 name.str("");
1620 name << glroot->PATH.Data() << "/";
1621 name << glroot->NAME.Data();
1622 //
1623 fcalname[s] = (TString)name.str().c_str();
1624 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]);
1625 //
1626 sgnl = LoadCalib(s);
1627 //
1628 if ( sgnl != 0 ) return(sgnl);
1629 delete glcalo;
1630 delete glroot;
1631 //
1632 return(0);
1633 //
1634 }
1635
1636 Int_t CaloLevel0::LoadCalib(Int_t s){
1637 //
1638 ifstream myfile;
1639 myfile.open(fcalname[s].Data());
1640 if ( !myfile ){
1641 return(-107);
1642 };
1643 myfile.close();
1644 //
1645 TFile *File = new TFile(fcalname[s].Data());
1646 if ( !File ) return(-108);
1647 TTree *tr = (TTree*)File->Get("CalibCalPed");
1648 if ( !tr ) return(-109);
1649 //
1650 TBranch *calo = tr->GetBranch("CalibCalPed");
1651 //
1652 pamela::CalibCalPedEvent *ce = 0;
1653 tr->SetBranchAddress("CalibCalPed", &ce);
1654 //
1655 Long64_t ncalibs = calo->GetEntries();
1656 //
1657 if ( !ncalibs ) return(-110);
1658 //
1659 calo->GetEntry(calibno[s]);
1660 //
1661 if (ce->cstwerr[s] != 0 && ce->cperror[s] == 0 ) {
1662 for ( Int_t d=0 ; d<11 ;d++ ){
1663 Int_t pre = -1;
1664 for ( Int_t j=0; j<96 ;j++){
1665 if ( j%16 == 0 ) pre++;
1666 if ( s == 2 ){
1667 calped[0][2*d+1][j] = ce->calped[3][d][j];
1668 cstwerr[3] = ce->cstwerr[3];
1669 cperror[3] = ce->cperror[3];
1670 calgood[0][2*d+1][j] = ce->calgood[3][d][j];
1671 calthr[0][2*d+1][pre] = ce->calthr[3][d][pre];
1672 calrms[0][2*d+1][j] = ce->calrms[3][d][j];
1673 calbase[0][2*d+1][pre] = ce->calbase[3][d][pre];
1674 calvar[0][2*d+1][pre] = ce->calvar[3][d][pre];
1675 };
1676 if ( s == 3 ){
1677 calped[0][2*d][j] = ce->calped[1][d][j];
1678 cstwerr[1] = ce->cstwerr[1];
1679 cperror[1] = ce->cperror[1];
1680 calgood[0][2*d][j] = ce->calgood[1][d][j];
1681 calthr[0][2*d][pre] = ce->calthr[1][d][pre];
1682 calrms[0][2*d][j] = ce->calrms[1][d][j];
1683 calbase[0][2*d][pre] = ce->calbase[1][d][pre];
1684 calvar[0][2*d][pre] = ce->calvar[1][d][pre];
1685 };
1686 if ( s == 0 ){
1687 calped[1][2*d][j] = ce->calped[0][d][j];
1688 cstwerr[0] = ce->cstwerr[0];
1689 cperror[0] = ce->cperror[0];
1690 calgood[1][2*d][j] = ce->calgood[0][d][j];
1691 calthr[1][2*d][pre] = ce->calthr[0][d][pre];
1692 calrms[1][2*d][j] = ce->calrms[0][d][j];
1693 calbase[1][2*d][pre] = ce->calbase[0][d][pre];
1694 calvar[1][2*d][pre] = ce->calvar[0][d][pre];
1695 };
1696 if ( s == 1 ){
1697 calped[1][2*d+1][j] = ce->calped[2][d][j];
1698 cstwerr[2] = ce->cstwerr[2];
1699 cperror[2] = ce->cperror[2];
1700 calgood[1][2*d+1][j] = ce->calgood[2][d][j];
1701 calthr[1][2*d+1][pre] = ce->calthr[2][d][pre];
1702 calrms[1][2*d+1][j] = ce->calrms[2][d][j];
1703 calbase[1][2*d+1][pre] = ce->calbase[2][d][pre];
1704 calvar[1][2*d+1][pre] = ce->calvar[2][d][pre];
1705 };
1706 };
1707 };
1708 } else {
1709 if ( verbose ) printf(" CALORIMETER - ERROR: problems finding a good calibration in this file! \n\n ");
1710 return(-111);
1711 };
1712 File->Close();
1713 return(0);
1714 }

  ViewVC Help
Powered by ViewVC 1.1.23