/[PAMELA software]/DarthVader/CalorimeterLevel2/src/CaloLevel0.cpp
ViewVC logotype

Annotation of /DarthVader/CalorimeterLevel2/src/CaloLevel0.cpp

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.25 - (hide annotations) (download)
Fri Dec 12 16:18:17 2008 UTC (16 years, 2 months ago) by mocchiut
Branch: MAIN
CVS Tags: v6r01, v6r00
Changes since 1.24: +25 -0 lines
Calo optimization + bug fixed, ToF errors -315/6/7 fixed

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

  ViewVC Help
Powered by ViewVC 1.1.23