/[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.33 - (hide annotations) (download)
Tue Nov 5 14:02:07 2013 UTC (11 years ago) by mocchiut
Branch: MAIN
Changes since 1.32: +7 -0 lines
New code for retracking

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

  ViewVC Help
Powered by ViewVC 1.1.23