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

  ViewVC Help
Powered by ViewVC 1.1.23