/[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.20 - (hide annotations) (download)
Fri Apr 18 18:55:47 2008 UTC (16 years, 7 months ago) by mocchiut
Branch: MAIN
Changes since 1.19: +3 -2 lines
ToF upgrade and bugs fixed, Calorimeter crosstalk bugs fixed, new ToF calibrations

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

  ViewVC Help
Powered by ViewVC 1.1.23