/[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.6 - (hide annotations) (download)
Mon Sep 3 08:42:13 2007 UTC (17 years, 3 months ago) by mocchiut
Branch: MAIN
Changes since 1.5: +12 -0 lines
Calorimeter: added selfdelay variable

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     calopar1 = true;
68     calopar2 = true;
69     calopar3 = true;
70     crosst = true;
71     ftcalopar1 = 0;
72     ttcalopar1 = 0;
73     ftcalopar2 = 0;
74     ttcalopar2 = 0;
75     ftcalopar3 = 0;
76     ttcalopar3 = 0;
77     }
78    
79     void CaloLevel0::SetCrossTalk(Bool_t ct){
80     crosst = ct;
81     };
82    
83     void CaloLevel0::SetVerbose(Bool_t ct){
84     verbose = ct;
85     };
86    
87     /**
88     * Initialize CaloLevel0 object
89     **/
90     void CaloLevel0::ProcessingInit(TSQLServer *dbc, UInt_t hs, Int_t &sgnl, TTree *l0tree, Bool_t isdeb, Bool_t isverb){
91     //
92     debug = isdeb;
93     verbose = isverb;
94     //
95     l0tr=(TTree*)l0tree;
96     de = new pamela::calorimeter::CalorimeterEvent();
97     l0calo = (TBranch*)l0tr->GetBranch("Calorimeter");
98     l0tr->SetBranchAddress("Calorimeter", &de);
99     //
100     trkseqno = 0;
101     ClearStructs();
102     //
103     GL_CALO_CALIB *glcalo = new GL_CALO_CALIB();
104     //
105     sgnl = 0;
106     UInt_t uptime = 0;
107     //
108     for (Int_t s = 0; s < 4; s++){
109     idcalib[s] = 0;
110     fromtime[s] = 0;
111     totime[s] = 0;
112     calibno[s] = 0;
113     ClearCalibVals(s);
114     //
115     sgnl = glcalo->Query_GL_CALO_CALIB(hs,uptime,s,dbc);
116     if ( sgnl < 0 ){
117     if ( verbose ) printf(" CALORIMETER - ERROR: error from GLTables\n");
118     return;
119     };
120     //
121     idcalib[s] = glcalo->ID_ROOT_L0;
122     fromtime[s] = glcalo->FROM_TIME;
123     if ( glcalo->TO_TIME < hs ){ // calibration is corrupted and we are using the one that preceed the good one
124     totime[s] = uptime;
125     } else {
126     totime[s] = glcalo->TO_TIME;
127     };
128     calibno[s] = glcalo->EV_ROOT;
129     //
130     if ( totime[s] == 0 ){
131     if ( verbose ) printf(" CALORIMETER - WARNING: data with no associated calibration\n");
132     ClearCalibVals(s);
133     sgnl = 100;
134     };
135     };
136     //
137     // determine path and name and entry of the calibration file
138     //
139     GL_ROOT *glroot = new GL_ROOT();
140     if ( verbose ) printf("\n");
141     for (Int_t s = 0; s < 4; s++){
142     if ( verbose ) printf(" ** SECTION %i **\n",s);
143     if ( totime[s] > 0 ){
144     //
145     sgnl = glroot->Query_GL_ROOT(idcalib[s],dbc);
146     if ( sgnl < 0 ){
147     if ( verbose ) printf(" CALORIMETER - ERROR: error from GLTables\n");
148     return;
149     };
150     //
151     stringstream name;
152     name.str("");
153     name << glroot->PATH.Data() << "/";
154     name << glroot->NAME.Data();
155     //
156     fcalname[s] = (TString)name.str().c_str();
157     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]);
158     } else {
159     if ( verbose ) printf(" - runheader at time %u. NO CALIBRATION INCLUDE THE RUNHEADER! ",hs);
160     };
161     sgnl = LoadCalib(s);
162     if ( sgnl ) break;
163     };
164     //
165     delete glcalo;
166     delete glroot;
167     //
168     return;
169     //
170     }
171    
172     Int_t CaloLevel0::ChkCalib(TSQLServer *dbc, UInt_t atime){
173     Int_t sgnl = 0;
174     for ( Int_t s = 0; s < 4; s++){
175     if ( atime > totime[s] ){
176     if ( !dbc->IsConnected() ) throw -116;
177     sgnl = Update(dbc,atime,s);
178     if ( sgnl < 0 ) return(sgnl);
179     };
180     };
181     return(sgnl);
182     }
183    
184 mocchiut 1.2 Int_t CaloLevel0::ChkParam(TSQLServer *dbc, UInt_t runheader, Bool_t mechal){
185 mocchiut 1.1 stringstream calfile;
186     stringstream bmfile;
187     stringstream aligfile;
188     Int_t error = 0;
189     FILE *f = 0;
190     ifstream badfile;
191     GL_PARAM *glparam = new GL_PARAM();
192     //
193     if ( calopar1 || ( ttcalopar1 != 0 && ttcalopar1 < runheader ) ){
194 mocchiut 1.2 //
195     //
196     //
197     if ( debug ) printf(" calopar1 %i ftcalopar1 %u ttcalopar1 %u runheader %u \n",calopar1,ftcalopar1,ttcalopar1,runheader);
198     //
199 mocchiut 1.1 calopar1 = false;
200     //
201     // determine where I can find calorimeter ADC to MIP conversion file
202     //
203     if ( verbose ) printf(" Querying DB for calorimeter parameters files...\n");
204     //
205     error = 0;
206     error = glparam->Query_GL_PARAM(runheader,101,dbc);
207     if ( error < 0 ) return(error);
208     //
209     calfile.str("");
210     calfile << glparam->PATH.Data() << "/";
211     calfile << glparam->NAME.Data();
212     ftcalopar1 = glparam->FROM_TIME;
213     ttcalopar1 = glparam->TO_TIME;
214     //
215     if ( verbose ) printf("\n Using ADC to MIP conversion file: \n %s \n",calfile.str().c_str());
216     f = fopen(calfile.str().c_str(),"rb");
217     if ( !f ){
218     if ( verbose ) printf(" CALORIMETER - ERROR: no ADC to MIP file!\n");
219     return(-105);
220     };
221     //
222     for (Int_t m = 0; m < 2 ; m++ ){
223     for (Int_t k = 0; k < 22; k++ ){
224     for (Int_t l = 0; l < 96; l++ ){
225     fread(&mip[m][k][l],sizeof(mip[m][k][l]),1,f);
226 mocchiut 1.5 if ( debug ) printf(" %f \n",mip[m][k][l]);
227 mocchiut 1.1 };
228     };
229     };
230     fclose(f);
231     };
232     //
233     if ( calopar2 || ( ttcalopar2 != 0 && ttcalopar2 < runheader ) ){
234 mocchiut 1.2 if ( debug ) printf(" calopar2 %i ftcalopar2 %u ttcalopar2 %u runheader %u \n",calopar2,ftcalopar2,ttcalopar2,runheader);
235 mocchiut 1.1 calopar2 = false;
236     //
237     // determine where I can find calorimeter alignment file
238     //
239     //
240     error = 0;
241     error = glparam->Query_GL_PARAM(runheader,102,dbc);
242     if ( error < 0 ) return(error);
243     //
244     aligfile.str("");
245     aligfile << glparam->PATH.Data() << "/";
246     aligfile << glparam->NAME.Data();
247     ftcalopar2 = glparam->FROM_TIME;
248     ttcalopar2 = glparam->TO_TIME;
249     //
250 mocchiut 1.2 if ( verbose ) printf("\n Using parameter file: \n %s \n",aligfile.str().c_str());
251 mocchiut 1.1 f = fopen(aligfile.str().c_str(),"rb");
252     if ( !f ){
253 mocchiut 1.2 if ( verbose ) printf(" CALORIMETER - ERROR: no parameter file!\n");
254 mocchiut 1.1 return(-106);
255     };
256     //
257 mocchiut 1.2 if ( !mechal ){
258     //
259     fread(&clevel1->xalig,sizeof(clevel1->xalig),1,f);
260     if ( debug ) printf(" xalig = %f \n",clevel1->xalig);
261     fread(&clevel1->yalig,sizeof(clevel1->yalig),1,f);
262     if ( debug ) printf(" yalig = %f \n",clevel1->yalig);
263     fread(&clevel1->zalig,sizeof(clevel1->zalig),1,f);
264     if ( debug ) printf(" zalig = %f \n",clevel1->zalig);
265     } else {
266     if ( verbose ) printf("\n Using MECHANICAL alignement parameters \n");
267     //
268     CaloStrip cs = CaloStrip();
269     cs.UseMechanicalAlig();
270     clevel1->xalig = cs.GetXalig();
271     if ( debug ) printf(" xalig = %f \n",clevel1->xalig);
272     clevel1->yalig = cs.GetYalig();
273     if ( debug ) printf(" yalig = %f \n",clevel1->yalig);
274     clevel1->zalig = cs.GetZalig();
275     if ( debug ) printf(" zalig = %f \n",clevel1->zalig);
276     //
277     Float_t tmp = 0;
278     fread(&tmp,sizeof(clevel1->xalig),1,f);
279     fread(&tmp,sizeof(clevel1->yalig),1,f);
280     fread(&tmp,sizeof(clevel1->zalig),1,f);
281     //
282     };
283 mocchiut 1.1 fread(&clevel1->emin,sizeof(clevel1->emin),1,f);
284     if ( debug ) printf(" signal threshold = %f \n",clevel1->emin);
285     //
286     fclose(f);
287     };
288     //
289     // Load offline bad strip mask
290     //
291     if ( calopar3 || ( ttcalopar3 != 0 && ttcalopar3 < runheader ) ){
292 mocchiut 1.2 if ( debug ) printf(" calopar3 %i ftcalopar3 %u ttcalopar3 %u runheader %u \n",calopar3,ftcalopar3,ttcalopar3,runheader);
293 mocchiut 1.1 calopar3 = false;
294     //
295     // determine where I can find calorimeter alignment file
296     //
297     //
298     error = 0;
299     error = glparam->Query_GL_PARAM(runheader,103,dbc);
300     if ( error < 0 ) return(error);
301     //
302     bmfile.str("");
303     bmfile << glparam->PATH.Data() << "/";
304     bmfile << glparam->NAME.Data();
305     ftcalopar3 = glparam->FROM_TIME;
306     ttcalopar3 = glparam->TO_TIME;
307     //
308     if ( verbose ) printf("\n Using bad strip offline mask file: \n %s \n\n",bmfile.str().c_str());
309     badfile.open(bmfile.str().c_str());
310     if ( !badfile ){
311     if ( verbose ) printf(" CALORIMETER - ERROR: no bad strip offline mask file!\n");
312     return(-115);
313     };
314     //
315     Bool_t isdone = false;
316     Int_t bad = 0;
317     Int_t view = 1;
318     Int_t strip = 0;
319     Int_t plane = 21;
320     while ( !isdone ) {
321     badfile >> bad;
322     obadmask[view][plane][strip] = bad;
323     if ( debug && bad ) printf(" SETTING view %i plane %i strip %i BAD = %i \n",view,plane,strip,bad);
324     strip++;
325     if ( strip > 95 ){
326     strip = 0;
327     plane--;
328     if ( plane < 0 ){
329     plane = 21;
330     view--;
331     };
332     if ( view < 0 ) isdone = true;
333     };
334     };
335     //
336     badfile.close();
337     };
338     //
339     delete glparam;
340     //
341     return(0);
342     }
343    
344    
345    
346     void CaloLevel0::FindBaseRaw(Int_t l, Int_t m, Int_t pre){
347     Float_t minstrip = 100000.;
348     Float_t rms = 0.;
349     base[l][m][pre] = 0.;
350     for (Int_t e = pre*16; e < (pre+1)*16 ; e++){
351     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.) {
352     minstrip = dexy[l][m][e]-calped[l][m][e];
353     rms = calthr[l][m][pre];
354     };
355     };
356     if ( minstrip != 100000. ) {
357     Float_t strip6s = 0.;
358     for (Int_t e = pre*16; e < (pre+1)*16 ; e++){
359     if ( (dexy[l][m][e]-calped[l][m][e]) >= minstrip && (dexy[l][m][e]-calped[l][m][e]) <= (minstrip+rms) ) {
360     strip6s += 1.;
361     base[l][m][pre] += (dexy[l][m][e] - calped[l][m][e]);
362     };
363     //
364     // compression
365     //
366     if ( abs((int)(dexy[l][m][e]-calped[l][m][e])) <= (minstrip+rms) ) {
367     dexyc[l][m][e] = 0.;
368     } else {
369     dexyc[l][m][e] = dexy[l][m][e];
370     };
371     };
372     if ( strip6s >= 9. ){
373     Double_t arro = base[l][m][pre]/strip6s;
374     Float_t deci = 1000.*((float)arro - float(int(arro)));
375     if ( deci < 500. ) {
376     arro = double(int(arro));
377     } else {
378     arro = 1. + double(int(arro));
379     };
380     base[l][m][pre] = arro;
381     } else {
382     base[l][m][pre] = 31000.;
383     for (Int_t e = pre*16; e < (pre+1)*16 ; e++){
384     dexyc[l][m][e] = dexy[l][m][e];
385     };
386     };
387     } else {
388     base[l][m][pre] = 31000.;
389     };
390     }
391    
392     Int_t CaloLevel0::Calibrate(Int_t ei){
393     //
394     // get entry ei
395     //
396     l0calo->GetEntry(ei);
397     //
398     // if it was not a selftrigger event, could it ever been a selftrigger event? if so trigty = 3.
399     //
400     Int_t val = 0;
401     Int_t del = 1100;
402 mocchiut 1.6 for (Int_t sec = 0; sec < 4; sec++){
403     for (Int_t dsec = 0; dsec < 7; dsec++){
404     val = (Int_t)de->calselftrig[sec][dsec];
405     del = delay(val);
406     clevel2->selfdelay[sec][dsec] = del;
407     };
408     };
409     val = 0;
410     del = 1100;
411 mocchiut 1.1 if ( clevel2->trigty != 2. ){
412     Bool_t ck = false;
413     for (Int_t sec = 0; sec < 4; sec++){
414     val = (Int_t)de->calselftrig[sec][6];
415     del = delay(val);
416     if ( del < 1100 ){
417     clevel2->wartrig = 0.;
418     clevel2->trigty = 3.;
419     ck = true;
420     break;
421     };
422     };
423     if ( !ck ) clevel2->wartrig = 100.;
424     } else {
425     Bool_t ck = false;
426     for (Int_t sec = 0; sec < 4; sec++){
427     val = (Int_t)de->calselftrig[sec][6];
428     del = delay(val);
429     if ( del < 1100 ){
430     clevel2->wartrig = 0.;
431     ck = true;
432     };
433     };
434     if ( !ck ) clevel2->wartrig = 100.;
435     };
436     //
437     Int_t se = 5;
438     Int_t done = 0;
439     Int_t pre = -1;
440     Bool_t isCOMP = false;
441     Bool_t isFULL = false;
442     Bool_t isRAW = false;
443     Float_t ener;
444     Int_t doneb = 0;
445     Int_t donec = 0;
446     Int_t ck = 0;
447     Int_t ipre = 0;
448     Int_t ip[3] = {0};
449     Float_t base0, base1, base2;
450     base0 = 0.;
451     base1 = 0.;
452     base2 = 0.;
453     Float_t qpre[6] = {0.,0.,0.,0.,0.,0.};
454     Float_t ene[96];
455     Int_t chdone[4] = {0,0,0,0};
456     Int_t pe = 0;
457     //
458     Float_t ener0 = 0.;
459     Float_t cbase0 = 0.;
460     Bool_t pproblem = false;
461     //
462     Float_t tim = 0.;
463     Int_t plo = 0;
464     Int_t fbi = 0;
465     Int_t cle = 0;
466     //
467     // run over views and planes
468     //
469     for (Int_t l = 0; l < 2; l++){
470     for (Int_t m = 0; m < 22; m++){
471     //
472 mocchiut 1.2 // determine the section number
473 mocchiut 1.1 //
474     se = 5;
475 mocchiut 1.2 if (l == 0 && m%2 == 0) se = 3;
476 mocchiut 1.1 if (l == 0 && m%2 != 0) se = 2;
477 mocchiut 1.2 if (l == 1 && m%2 != 0) se = 1;
478     if (l == 1 && m%2 == 0) se = 0;
479 mocchiut 1.1 //
480     // determine what kind of event we are going to analyze
481     //
482     isCOMP = false;
483     isFULL = false;
484     isRAW = false;
485     if ( de->stwerr[se] & (1 << 16) ) isCOMP = true;
486     if ( de->stwerr[se] & (1 << 17) ) isFULL = true;
487     if ( de->stwerr[se] & (1 << 3) ) isRAW = true;
488     if ( !chdone[se] ){
489     //
490     // check for any error in the event
491     //
492     clevel2->crc[se] = 0;
493     if ( de->perror[se] == 132 ){
494     clevel2->crc[se] = 1;
495     pe++;
496     };
497     clevel2->perr[se] = 0;
498     if ( de->perror[se] != 0 ){
499     clevel2->perr[se] = 1;
500     pe++;
501     };
502     clevel2->swerr[se] = 0;
503     for (Int_t j = 0; j < 7 ; j++){
504     if ( (j != 3) && (de->stwerr[se] & (1 << j)) ){
505     clevel2->swerr[se] = 1;
506     pe++;
507     };
508     };
509     chdone[se] = 1;
510     };
511     if ( clevel2->crc[se] == 0 && (clevel1->good2 == 1 || clevel2->trigty >= 2) ){
512     pre = -1;
513     //
514     for (Int_t nn = 0; nn < 96; nn++){
515     ene[nn] = 0.;
516     dexy[l][m][nn] = de->dexy[l][m][nn] ;
517     dexyc[l][m][nn] = de->dexyc[l][m][nn] ;
518     };
519     //
520     // run over preamplifiers
521     //
522     pre = -1;
523     cbase0 = 0.;
524     for (Int_t i = 0; i < 3; i++){
525     for (Int_t j = 0; j < 2; j++){
526     pre = j + i*2;
527     //
528     // baseline check and calculation
529     //
530     if ( !isRAW ) {
531     base[l][m][pre] = de->base[l][m][pre] ;
532     cbase0 += base[l][m][pre];
533     } else {
534     //
535     // if it is a raw event and we haven't checked
536     // yet, calculate the baseline.
537     //
538     FindBaseRaw(l,m,pre);
539     cbase0 += base[l][m][pre];
540     };
541     };
542     };
543     //
544     // run over strips
545     //
546     pre = -1;
547     ener0 = 0.;
548     for (Int_t i = 0 ; i < 3 ; i++){
549     ip[i] = 0;
550     for (Int_t n = i*32 ; n < (i+1)*32 ; n++){
551     if (n%16 == 0) {
552     ck = 0;
553     done = 0;
554     doneb = 0;
555     donec = 0;
556     pre++;
557     qpre[pre] = 0.;
558     };
559     //
560     // baseline check and calculation
561     //
562     // no suitable new baseline, use old ones!
563     //
564     if ( !done ){
565     if ( (base[l][m][pre] == 31000. || base[l][m][pre] == 0.) ){
566     ck = 1;
567     if (pre%2 == 0) {
568     ip[i] = pre + 1;
569     } else {
570     ip[i] = pre - 1;
571     };
572 mocchiut 1.3 if ( (base[l][m][ip[i]] == 31000. || base[l][m][ip[i]] == 0. || !crosst ) ){
573 mocchiut 1.1 //
574     ck = 2;
575     if ( sbase[l][m][pre] == 31000. || sbase[l][m][pre] == 0. ) {
576     ck = 3;
577     };
578     };
579     done = 1;
580     };
581     };
582     //
583     // CALIBRATION ALGORITHM
584     //
585     if ( !doneb ){
586 mocchiut 1.3 if ( debug ) printf(" ck is %i \n",ck);
587 mocchiut 1.1 switch (ck) {
588     case 0:
589     base0 = base[l][m][pre];
590     base2 = calbase[l][m][pre];
591 mocchiut 1.3 if ( debug ) printf(" base0 = base l m pre = %f base2 = calbase l m pre = %f \n",base[l][m][pre],calbase[l][m][pre]);
592 mocchiut 1.1 break;
593     case 1:
594     base0 = base[l][m][ip[i]];
595     base2 = calbase[l][m][ip[i]];
596 mocchiut 1.3 if ( debug ) printf(" base0 = base l m ip(i) = %f base2 = calbase l m ip(i) = %f \n",base[l][m][ip[i]],calbase[l][m][ip[i]]);
597 mocchiut 1.1 break;
598     case 2:
599     base0 = sbase[l][m][pre];
600 mocchiut 1.3 base2 = calbase[l][m][pre];
601     if ( debug ) printf(" base0 = sbase l m pre = %f base2 = calbase l m pre = %f \n",sbase[l][m][pre],calbase[l][m][pre]);
602 mocchiut 1.1 break;
603     case 3:
604     base0 = calbase[l][m][pre];
605     base2 = calbase[l][m][pre];
606 mocchiut 1.3 if ( debug ) printf(" base0 = calbase l m pre = %f base2 = calbase l m pre = %f \n",calbase[l][m][pre],calbase[l][m][pre]);
607 mocchiut 1.1 break;
608     };
609     base1 = calbase[l][m][pre];
610     doneb = 1;
611     };
612     ener = dexyc[l][m][n];
613     ener0 += ener;
614     clevel1->estrip[n][m][l] = 0.;
615     if ( base0>0 && base0 < 30000. ){
616 mocchiut 1.3 // if ( !donec && (base0 - base1 + base2) != 0. ){
617     // sbase[l][m][pre] = base0 - base1 + base2;
618     if ( !donec && (base0 + base1 - base2) != 0. ){
619     sbase[l][m][pre] = base0 + base1 - base2;
620 mocchiut 1.1 donec = 1;
621     };
622     if ( ener > 0. ){
623     clevel1->estrip[n][m][l] = (ener - calped[l][m][n] - base0 - base1 + base2)/mip[l][m][n] ;
624     //
625     // 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)
626     //
627     qpre[pre] += clevel1->estrip[n][m][l];
628 mocchiut 1.3 //
629     //
630 mocchiut 1.1 };
631     };
632     };
633     if ( crosst ){
634     if (ck == 1){
635     if (ip[i]%2 == 0) {
636     ipre = ip[i] + 1;
637     } else {
638     ipre = ip[i] - 1;
639     };
640     for (Int_t j = ipre*16 ; j < (ipre+1)*16 ; j++){
641     clevel1->estrip[j][m][l] += (qpre[ipre] - qpre[ip[i]]) * 0.00478;
642     };
643     };
644     if (ck == 2){
645     for (Int_t j = i*32 ; j < (i+1)*32 ; j++){
646     ipre = j/16 + 1;
647     clevel1->estrip[j][m][l] += qpre[ipre] * 0.00478;
648     };
649     };
650     };
651     };
652     //
653     if ( ener0 == 0. && cbase0 == 0. && !pproblem && clevel2->perr[se] == 0){
654     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);
655     pproblem = true;
656     pe++;
657     };
658     //
659     Int_t j4 = -4;
660     Int_t jjj = -3;
661     Int_t jj = -2;
662     for (Int_t j = 0 ; j < 100 ; j++){
663     jj++;
664     jjj++;
665     j4++;
666     if ( j < 96 ) ene[j] = clevel1->estrip[j][m][l];
667     if ( crosst ){
668     if ( jj >= 0 && jj < 96 ){
669     if ( jj != 0 && jj != 32 && jj != 64 ) ene[jj-1] += -clevel1->estrip[jj][m][l] * 0.01581;
670     if ( jj != 31 && jj != 63 && jj != 95 ) ene[jj+1] += -clevel1->estrip[jj][m][l] * 0.01581;
671     };
672     if ( jjj >= 0 && jjj < 96 ){
673     if ( jjj != 0 && jjj != 32 && jjj != 64 ) clevel1->estrip[jjj-1][m][l] += -ene[jjj] * 0.01581;
674     if ( jjj != 31 && jjj != 63 && jjj != 95 ) clevel1->estrip[jjj+1][m][l] += -ene[jjj] * 0.01581;
675     };
676     };
677     if ( j4 >= 0 && j4 < 96 ){
678     //
679     // NOTICE: THE FOLLOWING LINE EXCLUDE ALL STRIPS FOR WHICH THE RMS*4 IS GREATER THAN 26 !!! <=============== IMPORTANT! =================>
680     //
681     if ( obadmask[l][m][j4] == 1 || clevel1->estrip[j4][m][l] <= clevel1->emin || calrms[l][m][j4] > 26 ){
682     clevel1->estrip[j4][m][l] = 0.;
683     };
684     //
685     // code and save the energy for each strip in svstrip
686     //
687     if ( clevel1->estrip[j4][m][l] > clevel1->emin ){
688     //
689     tim = 100000.;
690     plo = m;
691     fbi = 0;
692     if ( clevel1->estrip[j4][m][l] > 0.99995 ){
693     tim = 10000.;
694     plo = m;
695     fbi = 1;
696     };
697     if ( clevel1->estrip[j4][m][l] > 9.9995 ){
698     tim = 1000.;
699     plo = 22 + m;
700     fbi = 1;
701     };
702     if ( clevel1->estrip[j4][m][l] > 99.995 ){
703     tim = 100.;
704     plo = 22 + m;
705     fbi = 0;
706     };
707     if ( clevel1->estrip[j4][m][l] > 999.95 ){
708     tim = 10.;
709     plo = 44 + m;
710     fbi = 0;
711     };
712     if ( clevel1->estrip[j4][m][l] > 9999.5 ){
713     tim = 1.;
714     plo = 66 + m;
715     fbi = 0;
716     };
717     //
718     cle = (Int_t)lroundf(tim*clevel1->estrip[j4][m][l]);
719     //
720     if ( l == 0 ){
721     //
722     // +-PPSSmmmm.mmmm
723     //
724     svstrip[istrip] = fbi*1000000000 + plo*10000000 + j4*100000 + cle;
725     } else {
726     svstrip[istrip] = -(fbi*1000000000 + plo*10000000 + j4*100000 + cle);
727     };
728     //
729     // if ( ei >= -770 ) printf(" j %i l %i m %i estrip %f \n",j4,l,m,clevel1->estrip[j4][m][l]);
730     // 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);
731     // if ( ei >= -770 ) printf(" svstrip %i \n",svstrip[istrip]);
732     //
733     istrip++;
734     };
735     };
736     };
737     //
738     } else {
739     for (Int_t nn = 0; nn < 96; nn++){
740     clevel1->estrip[nn][m][l] = 0.;
741     };
742     };
743     };
744     };
745     if ( !pe ){
746     clevel2->good = 1;
747     } else {
748     clevel2->good = 0;
749     };
750     return(0);
751     }
752    
753     void CaloLevel0::GetTrkVar(){
754     calol2tr();
755     }
756    
757     void CaloLevel0::FillTrkVar(CaloLevel2 *ca, Int_t nutrk){
758     //
759     CaloTrkVar *t_ca = new CaloTrkVar();
760     //
761     t_ca->trkseqno = trkseqno;
762     t_ca->ncore = (Int_t)clevel2->ncore;
763     t_ca->qcore = clevel2->qcore;
764     t_ca->noint = (Int_t)clevel2->noint;
765     t_ca->ncyl = (Int_t)clevel2->ncyl;
766     t_ca->qcyl = clevel2->qcyl;
767     t_ca->qtrack = clevel2->qtrack;
768     t_ca->qtrackx = clevel2->qtrackx;
769     t_ca->qtracky = clevel2->qtracky;
770     t_ca->dxtrack = clevel2->dxtrack;
771     t_ca->dytrack = clevel2->dytrack;
772     t_ca->qlast = clevel2->qlast;
773     t_ca->nlast = (Int_t)clevel2->nlast;
774     t_ca->qpre = clevel2->qpre;
775     t_ca->npre = (Int_t)clevel2->npre;
776     t_ca->qpresh = clevel2->qpresh;
777     t_ca->npresh = (Int_t)clevel2->npresh;
778     t_ca->qtr = clevel2->qtr;
779     t_ca->ntr = (Int_t)clevel2->ntr;
780     t_ca->planetot = (Int_t)clevel2->planetot;
781     t_ca->qmean = clevel2->qmean;
782     t_ca->dX0l = clevel2->dX0l;
783     t_ca->qlow = clevel2->qlow;
784     t_ca->nlow = (Int_t)clevel2->nlow;
785     //
786     if ( trkseqno == -1 ){
787     // ca->impx = clevel2->impx;
788     // ca->impy = clevel2->impy;
789     ca->tanx[1] = clevel2->tanx;
790     ca->tany[1] = clevel2->tany;
791     ca->elen = clevel2->elen;
792     ca->selen = clevel2->selen;
793     // memcpy(ca->cibar,clevel2->cibar,sizeof(clevel2->cibar));
794     // memcpy(ca->cbar,clevel2->cbar,sizeof(clevel2->cbar));
795     memcpy(t_ca->tibar,clevel2->cibar,sizeof(clevel2->cibar));
796     memcpy(t_ca->tbar,clevel2->cbar,sizeof(clevel2->cbar));
797     memcpy(ca->planemax,clevel2->planemax,sizeof(clevel2->planemax));
798 mocchiut 1.6 memcpy(ca->selfdelay,clevel2->selfdelay,sizeof(clevel2->selfdelay));
799 mocchiut 1.1 ca->varcfit[2] = clevel2->varcfit[0];
800     ca->varcfit[3] = clevel2->varcfit[1];
801     ca->npcfit[2] = clevel2->npcfit[0];
802     ca->npcfit[3] = clevel2->npcfit[1];
803     // memcpy(ca->varcfit,clevel2->varcfit,sizeof(clevel2->varcfit));
804     // memcpy(ca->npcfit,clevel2->npcfit,sizeof(clevel2->npcfit));
805     } else {
806     memcpy(t_ca->tibar,clevel2->tibar,sizeof(clevel2->tibar));
807     memcpy(t_ca->tbar,clevel2->tbar,sizeof(clevel2->tbar));
808     };
809     //
810     //
811     if ( !(ca->CaloTrk) ) ca->CaloTrk = new TClonesArray("CaloTrkVar",1); //ELENA
812     TClonesArray &t = *ca->CaloTrk;
813     new(t[nutrk]) CaloTrkVar(*t_ca);
814     //
815     delete t_ca;
816     //
817     ClearTrkVar();
818     }
819    
820     void CaloLevel0::GetCommonVar(){
821     calol2cm();
822     }
823    
824     void CaloLevel0::FillCommonVar(CaloLevel1 *c1, CaloLevel2 *ca){
825     //
826     ca->good = clevel2->good;
827     if ( clevel2->trigty == 2. ){
828     ca->selftrigger = 1;
829     } else {
830     ca->selftrigger = 0;
831     };
832     //
833     ca->selftrigger += (Int_t)clevel2->wartrig;
834     //
835     memcpy(ca->perr,clevel2->perr,sizeof(clevel2->perr));
836     memcpy(ca->swerr,clevel2->swerr,sizeof(clevel2->swerr));
837     memcpy(ca->crc,clevel2->crc,sizeof(clevel2->crc));
838     ca->nstrip = (Int_t)clevel2->nstrip;
839     ca->qtot = clevel2->qtot;
840     // ca->impx = clevel2->impx;
841     // ca->impy = clevel2->impy;
842     ca->tanx[0] = clevel2->tanx;
843     ca->tany[0] = clevel2->tany;
844     ca->nx22 = (Int_t)clevel2->nx22;
845     ca->qx22 = clevel2->qx22;
846     ca->qmax = clevel2->qmax;
847     ca->elen = clevel2->elen;
848     ca->selen = clevel2->selen;
849     memcpy(ca->qq,clevel2->qq,sizeof(clevel2->qq));
850     memcpy(ca->planemax,clevel2->planemax,sizeof(clevel2->planemax));
851 mocchiut 1.6 memcpy(ca->selfdelay,clevel2->selfdelay,sizeof(clevel2->selfdelay));
852 mocchiut 1.1 ca->varcfit[0] = clevel2->varcfit[0];
853     ca->varcfit[1] = clevel2->varcfit[1];
854     ca->npcfit[0] = clevel2->npcfit[0];
855     ca->npcfit[1] = clevel2->npcfit[1];
856     ca->fitmode[0] = clevel2->fmode[0];
857     ca->fitmode[1] = clevel2->fmode[1];
858     // memcpy(ca->varcfit,clevel2->varcfit,sizeof(clevel2->varcfit));
859     // memcpy(ca->npcfit,clevel2->npcfit,sizeof(clevel2->npcfit));
860     memcpy(ca->cibar,clevel2->cibar,sizeof(clevel2->cibar));
861     memcpy(ca->cbar,clevel2->cbar,sizeof(clevel2->cbar));
862     //
863     if ( c1 ){
864     c1->istrip = istrip;
865     c1->estrip = TArrayI(istrip,svstrip);
866     };
867     //
868     }
869    
870     void CaloLevel0::ClearStructs(){
871     ClearTrkVar();
872     ClearCommonVar();
873     }
874    
875     void CaloLevel0::RunClose(){
876     l0tr->Delete();
877     ClearStructs();
878     //
879     memset(dexy, 0, 2*22*96*sizeof(Float_t));
880     memset(dexyc, 0, 2*22*96*sizeof(Float_t));
881     memset(base, 0, 2*22*6*sizeof(Float_t));
882     memset(sbase, 0, 2*22*6*sizeof(Float_t));
883     //
884     }
885    
886     //
887     // Private methods
888     //
889    
890     void CaloLevel0::ClearTrkVar(){
891     clevel2->ncore = 0;
892     clevel2->qcore = 0.;
893     clevel2->noint = 0.;
894     clevel2->ncyl = 0.;
895     clevel2->qcyl = 0.;
896     clevel2->qtrack = 0.;
897     clevel2->qtrackx = 0.;
898     clevel2->qtracky = 0.;
899     clevel2->dxtrack = 0.;
900     clevel2->dytrack = 0.;
901     clevel2->qlast = 0.;
902     clevel2->nlast = 0.;
903     clevel2->qpre = 0.;
904     clevel2->npre = 0.;
905     clevel2->qpresh = 0.;
906     clevel2->npresh = 0.;
907     clevel2->qlow = 0.;
908     clevel2->nlow = 0.;
909     clevel2->qtr = 0.;
910     clevel2->ntr = 0.;
911     clevel2->planetot = 0.;
912     clevel2->qmean = 0.;
913     clevel2->dX0l = 0.;
914     clevel2->elen = 0.;
915     clevel2->selen = 0.;
916     memset(clevel1->al_p, 0, 5*2*sizeof(Double_t));
917     memset(clevel2->tibar, 0, 2*22*sizeof(Int_t));
918     memset(clevel2->tbar, 0, 2*22*sizeof(Float_t));
919     }
920    
921     void CaloLevel0::ClearCommonVar(){
922     istrip = 0;
923     clevel2->trigty = -1.;
924     clevel2->wartrig = 0.;
925     clevel2->good = 0;
926     clevel2->nstrip = 0.;
927     clevel2->qtot = 0.;
928     // clevel2->impx = 0.;
929     // clevel2->impy = 0.;
930     clevel2->tanx = 0.; // this is correct since it refers to the fortran structure
931     clevel2->tany = 0.; // this is correct since it refers to the fortran structure
932     clevel2->qmax = 0.;
933     clevel2->nx22 = 0.;
934     clevel2->qx22 = 0.;
935     memset(clevel2->perr, 0, 4*sizeof(Int_t));
936     memset(clevel2->swerr, 0, 4*sizeof(Int_t));
937     memset(clevel2->crc, 0, 4*sizeof(Int_t));
938     memset(clevel2->qq, 0, 4*sizeof(Int_t));
939     memset(clevel2->varcfit, 0, 4*sizeof(Float_t));
940     memset(clevel2->npcfit, 0, 4*sizeof(Int_t));
941     memset(clevel2->planemax, 0, 2*sizeof(Int_t));
942 mocchiut 1.6 memset(clevel2->selfdelay, 0, 4*7*sizeof(Int_t));
943 mocchiut 1.1 memset(clevel2->fmode, 0, 2*sizeof(Int_t));
944     memset(clevel2->cibar, 0, 2*22*sizeof(Int_t));
945     memset(clevel2->cbar, 0, 2*22*sizeof(Float_t));
946     }
947    
948     void CaloLevel0::ClearCalibVals(Int_t s){
949     //
950     for ( Int_t d=0 ; d<11 ;d++ ){
951     Int_t pre = -1;
952     for ( Int_t j=0; j<96 ;j++){
953     if ( j%16 == 0 ) pre++;
954     if ( s == 2 ){
955     calped[0][2*d+1][j] = 0.;
956     cstwerr[3] = 0.;
957     cperror[3] = 0.;
958     calgood[0][2*d+1][j] = 0.;
959     calthr[0][2*d+1][pre] = 0.;
960     calrms[0][2*d+1][j] = 0.;
961     calbase[0][2*d+1][pre] = 0.;
962     calvar[0][2*d+1][pre] = 0.;
963     };
964     if ( s == 3 ){
965     calped[0][2*d][j] = 0.;
966     cstwerr[1] = 0.;
967     cperror[1] = 0.;
968     calgood[0][2*d][j] = 0.;
969     calthr[0][2*d][pre] = 0.;
970     calrms[0][2*d][j] = 0.;
971     calbase[0][2*d][pre] = 0.;
972     calvar[0][2*d][pre] = 0.;
973     };
974     if ( s == 0 ){
975     calped[1][2*d][j] = 0.;
976     cstwerr[0] = 0.;
977     cperror[0] = 0.;
978     calgood[1][2*d][j] = 0.;
979     calthr[1][2*d][pre] = 0.;
980     calrms[1][2*d][j] = 0.;
981     calbase[1][2*d][pre] = 0.;
982     calvar[1][2*d][pre] = 0.;
983     };
984     if ( s == 1 ){
985     calped[1][2*d+1][j] = 0.;
986     cstwerr[2] = 0.;
987     cperror[2] = 0.;
988     calgood[1][2*d+1][j] = 0.;
989     calthr[1][2*d+1][pre] = 0.;
990     calrms[1][2*d+1][j] = 0.;
991     calbase[1][2*d+1][pre] = 0.;
992     calvar[1][2*d+1][pre] = 0.;
993     };
994     };
995     };
996     return;
997     }
998    
999     Int_t CaloLevel0::Update(TSQLServer *dbc, UInt_t atime, Int_t s){
1000     //
1001     Int_t sgnl = 0;
1002     //
1003     GL_CALO_CALIB *glcalo = new GL_CALO_CALIB();
1004     //
1005     sgnl = 0;
1006     //
1007     idcalib[s] = 0;
1008     fromtime[s] = 0;
1009     totime[s] = 0;
1010     calibno[s] = 0;
1011     ClearCalibVals(s);
1012     //
1013     UInt_t uptime = 0;
1014     //
1015     sgnl = glcalo->Query_GL_CALO_CALIB(atime,uptime,s,dbc);
1016     if ( sgnl < 0 ){
1017     if ( verbose ) printf(" CALORIMETER - ERROR: error from GLTables\n");
1018     return(sgnl);
1019     };
1020     //
1021     idcalib[s] = glcalo->ID_ROOT_L0;
1022     fromtime[s] = glcalo->FROM_TIME;
1023     if ( glcalo->TO_TIME < atime ){ // calibration is corrupted and we are using the one that preceed the good one
1024     totime[s] = uptime;
1025     } else {
1026     totime[s] = glcalo->TO_TIME;
1027     };
1028     // totime[s] = glcalo->TO_TIME;
1029     calibno[s] = glcalo->EV_ROOT;
1030     //
1031     if ( totime[s] == 0 ){
1032     if ( verbose ) printf(" CALORIMETER - WARNING: data with no associated calibration\n");
1033     ClearCalibVals(s);
1034     sgnl = 100;
1035     };
1036     //
1037     // determine path and name and entry of the calibration file
1038     //
1039     GL_ROOT *glroot = new GL_ROOT();
1040     if ( verbose ) printf("\n");
1041     if ( verbose ) printf(" ** SECTION %i **\n",s);
1042     //
1043     sgnl = glroot->Query_GL_ROOT(idcalib[s],dbc);
1044     if ( sgnl < 0 ){
1045     if ( verbose ) printf(" CALORIMETER - ERROR: error from GLTables\n");
1046     return(sgnl);
1047     };
1048     //
1049     stringstream name;
1050     name.str("");
1051     name << glroot->PATH.Data() << "/";
1052     name << glroot->NAME.Data();
1053     //
1054     fcalname[s] = (TString)name.str().c_str();
1055     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]);
1056     //
1057     sgnl = LoadCalib(s);
1058     //
1059     if ( sgnl != 0 ) return(sgnl);
1060     delete glcalo;
1061     delete glroot;
1062     //
1063     return(0);
1064     //
1065     }
1066    
1067     Int_t CaloLevel0::LoadCalib(Int_t s){
1068     //
1069     ifstream myfile;
1070     myfile.open(fcalname[s].Data());
1071     if ( !myfile ){
1072     return(-107);
1073     };
1074     myfile.close();
1075     //
1076     TFile *File = new TFile(fcalname[s].Data());
1077     if ( !File ) return(-108);
1078     TTree *tr = (TTree*)File->Get("CalibCalPed");
1079     if ( !tr ) return(-109);
1080     //
1081     TBranch *calo = tr->GetBranch("CalibCalPed");
1082     //
1083     pamela::CalibCalPedEvent *ce = 0;
1084     tr->SetBranchAddress("CalibCalPed", &ce);
1085     //
1086     Long64_t ncalibs = calo->GetEntries();
1087     //
1088     if ( !ncalibs ) return(-110);
1089     //
1090     calo->GetEntry(calibno[s]);
1091     //
1092     if (ce->cstwerr[s] != 0 && ce->cperror[s] == 0 ) {
1093     for ( Int_t d=0 ; d<11 ;d++ ){
1094     Int_t pre = -1;
1095     for ( Int_t j=0; j<96 ;j++){
1096     if ( j%16 == 0 ) pre++;
1097     if ( s == 2 ){
1098     calped[0][2*d+1][j] = ce->calped[3][d][j];
1099     cstwerr[3] = ce->cstwerr[3];
1100     cperror[3] = ce->cperror[3];
1101     calgood[0][2*d+1][j] = ce->calgood[3][d][j];
1102     calthr[0][2*d+1][pre] = ce->calthr[3][d][pre];
1103     calrms[0][2*d+1][j] = ce->calrms[3][d][j];
1104     calbase[0][2*d+1][pre] = ce->calbase[3][d][pre];
1105     calvar[0][2*d+1][pre] = ce->calvar[3][d][pre];
1106     };
1107     if ( s == 3 ){
1108     calped[0][2*d][j] = ce->calped[1][d][j];
1109     cstwerr[1] = ce->cstwerr[1];
1110     cperror[1] = ce->cperror[1];
1111     calgood[0][2*d][j] = ce->calgood[1][d][j];
1112     calthr[0][2*d][pre] = ce->calthr[1][d][pre];
1113     calrms[0][2*d][j] = ce->calrms[1][d][j];
1114     calbase[0][2*d][pre] = ce->calbase[1][d][pre];
1115     calvar[0][2*d][pre] = ce->calvar[1][d][pre];
1116     };
1117     if ( s == 0 ){
1118     calped[1][2*d][j] = ce->calped[0][d][j];
1119     cstwerr[0] = ce->cstwerr[0];
1120     cperror[0] = ce->cperror[0];
1121     calgood[1][2*d][j] = ce->calgood[0][d][j];
1122     calthr[1][2*d][pre] = ce->calthr[0][d][pre];
1123     calrms[1][2*d][j] = ce->calrms[0][d][j];
1124     calbase[1][2*d][pre] = ce->calbase[0][d][pre];
1125     calvar[1][2*d][pre] = ce->calvar[0][d][pre];
1126     };
1127     if ( s == 1 ){
1128     calped[1][2*d+1][j] = ce->calped[2][d][j];
1129     cstwerr[2] = ce->cstwerr[2];
1130     cperror[2] = ce->cperror[2];
1131     calgood[1][2*d+1][j] = ce->calgood[2][d][j];
1132     calthr[1][2*d+1][pre] = ce->calthr[2][d][pre];
1133     calrms[1][2*d+1][j] = ce->calrms[2][d][j];
1134     calbase[1][2*d+1][pre] = ce->calbase[2][d][pre];
1135     calvar[1][2*d+1][pre] = ce->calvar[2][d][pre];
1136     };
1137     };
1138     };
1139     } else {
1140     if ( verbose ) printf(" CALORIMETER - ERROR: problems finding a good calibration in this file! \n\n ");
1141     return(-111);
1142     };
1143     File->Close();
1144     return(0);
1145     }

  ViewVC Help
Powered by ViewVC 1.1.23