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

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

Parent Directory Parent Directory | Revision Log Revision Log

Revision 1.15 - (hide annotations) (download)
Tue Nov 29 13:50:11 2011 UTC (13 years ago) by pam-fi
Branch: MAIN
Changes since 1.14: +2 -0 lines 
Some fixes to prevent memory leaks or double delete.
1 mocchiut 1.1 /**
2     * \file src/CaloLevel1.cpp
3     * \author Emiliano Mocchiutti
4     *
5     **/
6     #include <CaloLevel1.h>
7     //
8     ClassImp(CaloStrip);
9     ClassImp(CaloLevel1);
10    
11 mocchiut 1.8 Float_t CaloStrip::UXal = CTX;
12     Float_t CaloStrip::UYal = CTY;
13     Float_t CaloStrip::UZal = CTZ;
14     Bool_t CaloStrip::paramload = false;
15    
16 mocchiut 1.1 /**
17     * CaloStrip default constructor
18     **/
19     CaloStrip::CaloStrip() {
20     c1 = 0;
21     this->Clear();
22 mocchiut 1.5 };
23    
24     /**
25     * CaloStrip default constructor
26     **/
27     CaloStrip::CaloStrip(Bool_t mechalig) {
28     c1 = 0;
29     this->Clear();
30     if ( mechalig ){
31     ismech = true;
32 mocchiut 1.8 paramload = true;
33 mocchiut 1.5 UXal = MECHCTX;
34     UYal = MECHCTY;
35     UZal = MECHCTZ;
36     } else {
37     ismech = false;
38     UseStandardAlig();
39     };
40     };
41    
42     /**
43     * CaloStrip default constructor
44     **/
45     CaloStrip::CaloStrip(CaloLevel1 *calo) {
46     c1 = calo->GetCaloLevel1();
47     this->Clear();
48     ismech = false;
49 mocchiut 1.4 UseStandardAlig();
50 mocchiut 1.1 };
51    
52     /**
53     * CaloStrip default constructor
54     **/
55 mocchiut 1.5 CaloStrip::CaloStrip(CaloLevel1 *calo, Bool_t mechalig) {
56 mocchiut 1.1 c1 = calo->GetCaloLevel1();
57     this->Clear();
58 mocchiut 1.5 if ( mechalig ){
59     ismech = true;
60 mocchiut 1.8 paramload = true;
61 mocchiut 1.5 UXal = MECHCTX;
62     UYal = MECHCTY;
63     UZal = MECHCTZ;
64     } else {
65     ismech = false;
66     UseStandardAlig();
67     };
68 mocchiut 1.1 };
69    
70     /**
71     * Clear variables
72     **/
73 mocchiut 1.10 void CaloStrip::Clear(Option_t *t) {
74 mocchiut 1.1 fE = 0.;
75     fX = 0.;
76     fY = 0.;
77     fZ = 0.;
78     fView = 0;
79     fPlane = 0;
80     fStrip = 0;
81     };
82    
83     /**
84 mocchiut 1.4 * Connect to the DB and retrieve alignement parameters
85     **/
86     void CaloStrip::UseStandardAlig(){
87     //
88 mocchiut 1.8 if ( !paramload ){
89 mocchiut 1.4 //
90 mocchiut 1.8 paramload = true;
91     ismech = false;
92 mocchiut 1.4 //
93 mocchiut 1.8 stringstream aligfile;
94     Int_t error = 0;
95     FILE *f = 0;
96     ifstream badfile;
97     GL_PARAM *glparam = new GL_PARAM();
98     //
99     TString host = "mysql://localhost/pamelaprod";
100     TString user = "anonymous";
101     TString psw = "";
102     const char *pamdbhost=gSystem->Getenv("PAM_DBHOST");
103     const char *pamdbuser=gSystem->Getenv("PAM_DBUSER");
104     const char *pamdbpsw=gSystem->Getenv("PAM_DBPSW");
105     if ( !pamdbhost ) pamdbhost = "";
106     if ( !pamdbuser ) pamdbuser = "";
107     if ( !pamdbpsw ) pamdbpsw = "";
108     if ( strcmp(pamdbhost,"") ) host = pamdbhost;
109     if ( strcmp(pamdbuser,"") ) user = pamdbuser;
110     if ( strcmp(pamdbpsw,"") ) psw = pamdbpsw;
111     TSQLServer *dbc = TSQLServer::Connect(host.Data(),user.Data(),psw.Data());
112     //
113     UXal = 0.;
114     UYal = 0.;
115     UZal = 0.;
116 mocchiut 1.4 //
117 mocchiut 1.8 if ( dbc ){
118     //
119     // determine where I can find calorimeter ADC to MIP conversion file
120     //
121     printf(" Querying DB for calorimeter parameters files...\n");
122 mocchiut 1.4 //
123     //
124 mocchiut 1.8 //
125     error = 0;
126     error = glparam->Query_GL_PARAM(1000,102,dbc);
127     if ( error >= 0 ){
128     //
129     aligfile.str("");
130     aligfile << glparam->PATH.Data() << "/";
131     aligfile << glparam->NAME.Data();
132 mocchiut 1.4 //
133 mocchiut 1.8 printf("\n Using parameter file: \n %s \n",aligfile.str().c_str());
134     f = fopen(aligfile.str().c_str(),"rb");
135     if ( f ){
136     //
137     fread(&UXal,sizeof(UXal),1,f);
138     fread(&UYal,sizeof(UYal),1,f);
139     fread(&UZal,sizeof(UZal),1,f);
140     //
141     fclose(f);
142     };
143 mocchiut 1.4 //
144     };
145 mocchiut 1.8 dbc->Close();
146 pam-fi 1.15 delete dbc;
147     dbc = 0;
148 mocchiut 1.8 };
149     if ( !UXal ){
150 mocchiut 1.4 //
151 mocchiut 1.8 printf(" No able to query DB for calorimeter parameters files\n Using hard-coded parameters \n");
152     UXal = CTX;
153     UYal = CTY;
154     UZal = CTZ;
155 mocchiut 1.4 };
156 mocchiut 1.8 //
157 mocchiut 1.4 };
158     //
159     };
160    
161     /**
162 mocchiut 1.1 * Given a strip returns its position in the PAMELA reference system
163     **/
164     void CaloStrip::Set(Int_t view, Int_t plane, Int_t strip) {
165     //
166     this->Clear();
167     //
168     if ( view < 0 || view > 1 ){
169     printf(" ERROR: 0 =< view =< 1 \n");
170     return;
171     };
172     if ( plane < 0 || plane > 21 ){
173     printf(" ERROR: 0 =< plane =< 21 \n");
174     return;
175     };
176     if ( strip < 0 || strip > 95 ){
177     printf(" ERROR: 0 =< strip =< 95 \n");
178     return;
179     };
180     //
181     Float_t lShift = 0.;
182     Float_t lPos = 0.;
183     extern struct shift shift_;
184     //
185     // Find MIPs for that strip
186     //
187     if ( c1 ) fE = c1->GetEstrip(view, plane, strip);
188     //
189     fView = view + 1;
190     fPlane = plane + 1;
191     fStrip = strip + 1;
192 mocchiut 1.7 //
193 mocchiut 1.1 if ( fPlane%2 ){
194 mocchiut 1.6 lShift = +0.5;
195     } else {
196 mocchiut 1.1 lShift = -0.5;
197     };
198     //
199     shift_.shift = lShift;
200     //
201     Float_t zplane[22];
202     zplane[0] = 0.;
203     Int_t ii = 0;
204     for ( Int_t i = 2; i < 23; i++){
205     ii = i-1;
206     if ( i%2 ){
207     zplane[ii] = zplane[ii-1] - 10.09;
208     } else {
209     zplane[ii] = zplane[ii-1] - 8.09;
210     };
211     };
212     //
213     millim_(&fStrip,&lPos);
214     //
215     if ( fView == 1 ){
216     //
217     // X view
218     //
219 mocchiut 1.4 fX = (lPos - UXal)/10.;
220 mocchiut 1.1 fY = 0.;
221 mocchiut 1.4 fZ = (zplane[fPlane-1] - 5.81 + UZal)/10.;
222 mocchiut 1.6 //
223     } else {
224 mocchiut 1.1 //
225     // Y view
226     //
227     fX = 0;
228 mocchiut 1.4 fY = (lPos - UYal)/10.;
229     fZ = (zplane[fPlane-1] + UZal)/10.;
230 mocchiut 1.1 };
231     //
232     };
233    
234     /**
235     * Given a point in the space (PAMELA ref system) returns the closest strip
236     **/
237     void CaloStrip::Set(Float_t X, Float_t Y, Float_t Z) {
238     //
239 mocchiut 1.2 fX = X;
240     fY = Y;
241     fZ = Z;
242     //
243     Float_t zplane[22];
244     zplane[0] = 0.;
245     Int_t ii = 0;
246     for ( Int_t i = 2; i < 23; i++){
247     ii = i-1;
248     if ( i%2 ){
249     zplane[ii] = zplane[ii-1] - 10.09;
250     } else {
251     zplane[ii] = zplane[ii-1] - 8.09;
252     };
253     };
254     //
255     Float_t dzx[22];
256     Float_t dzy[22];
257     for ( Int_t i=0; i < 22; i++){
258 mocchiut 1.4 dzx[i] = fabs(fZ*10. - (zplane[i] - 5.81 + UZal));
259     dzy[i] = fabs(fZ*10. - (zplane[i] + UZal));
260 mocchiut 1.2 };
261     //
262     Float_t minx = TMath::MinElement(22,dzx);
263     Float_t miny = TMath::MinElement(22,dzy);
264     //
265     // find view
266     //
267     if ( minx < miny ){
268     fView = 1;
269     } else {
270     fView = 2;
271     };
272     //
273     // find plane
274     //
275     Float_t pos = 0.;
276     //
277     for ( Int_t i=0; i < 22; i++){
278     if ( fView == 1 ){
279     if ( dzx[i] == minx ){
280     fPlane = i+1;
281 mocchiut 1.4 pos = fX*10. + UXal;
282 mocchiut 1.2 };
283     } else {
284     if ( dzy[i] == miny ){
285     fPlane = i+1;
286 mocchiut 1.4 pos = fY*10. + UYal;
287 mocchiut 1.2 };
288     };
289     };
290     //
291     // find strip
292     //
293     Float_t dxy[96];
294     Float_t stpos = 0.;
295     //
296     CaloStrip *ca = new CaloStrip();
297     //
298     for ( Int_t i=0; i < 96; i++){
299     ca->Set(fView-1,fPlane-1,i);
300     if ( fView == 1 ){
301 mocchiut 1.4 stpos = ca->GetX()*10. + UXal;
302 mocchiut 1.2 } else {
303 mocchiut 1.4 stpos = ca->GetY()*10. + UYal;
304 mocchiut 1.2 };
305     dxy[i] = fabs(pos - stpos);
306     };
307     //
308     delete ca;
309 mocchiut 1.1 //
310 mocchiut 1.2 Float_t mins = TMath::MinElement(96,dxy);
311     //
312     for ( Int_t i=0; i < 96; i++){
313     if ( dxy[i] == mins ) fStrip = i+1;
314     };
315 mocchiut 1.1 };
316    
317     /**
318 mocchiut 1.14 * Given a point in the space or a strip it returns the Silicon sensor number. Numbering goes like this:
319     *
320     * y ^
321     * |
322     * | 6 7 8
323     * | 3 4 5
324     * | 0 1 2
325     * | -----------> x
326     *
327     **/
328     Int_t CaloStrip::GetSiSensor() {
329     //
330     // fX fY fZ // fView fPlane
331     //
332     Float_t deadsi = 0.096;
333     Float_t dead = 0.05;
334     Float_t sidim = 8.00;
335     // Float_t stripdim = 0.244;
336     Float_t sensitarea = 7.808;
337     //
338     Float_t xoffset = 0.;
339     Float_t yoffset = 0.;
340     //
341     if ( (fView-1) == 0 && !((fPlane-1)%2) ){
342     xoffset = +0.05;
343     yoffset = 0.0;
344     };
345     if ( (fView-1) == 0 && ((fPlane-1)%2) ){
346     xoffset = -0.05;
347     yoffset = -0.20;
348     };
349     if ( (fView-1) == 1 && !((fPlane-1)%2) ){
350     xoffset = +0.10;
351     yoffset = -0.05;
352     };
353     if ( (fView-1) == 1 && ((fPlane-1)%2) ){
354     xoffset = -0.10;
355     yoffset = -0.15;
356     };
357     //
358     Int_t iind = -1;
359     Int_t jind = -1;
360     //
361     for (Int_t i = 0; i < 3; i++){
362     if ( (fX+xoffset+12.1) >= (deadsi+(sidim+dead)*i) && (fX+xoffset+12.1) <= (sensitarea+deadsi+(sidim+dead)*i) ){
363     iind = i;
364     break;
365     };
366     };
367     //
368     for (Int_t j = 0; j < 3; j++){
369     if ( (fY+yoffset+12.1) >= (deadsi+(sidim+dead)*j) && (fY+yoffset+12.1) <= (sensitarea+deadsi+(sidim+dead)*j) ){
370     jind = j;
371     break;
372     };
373     };
374     //
375     Int_t sensor = -1;
376     if ( iind != -1 && jind != -1 ){
377     sensor = iind + jind * 3;
378     };
379     //
380     // printf(" View %i Plane %i x %f y %f z %f xoffset %f yoffset %f iind %i jind %i \n",fView,fPlane,fX,fY,fZ,xoffset,yoffset,iind,jind);
381     //
382     return(sensor);
383     //
384     };
385    
386     /**
387 mocchiut 1.1 * CaloLevel1 constructor
388     **/
389     CaloLevel1::CaloLevel1() {
390     //
391     estrip = TArrayI(0,NULL);
392     //
393     this->Clear();
394     //
395     };
396    
397     /**
398     * Clear the CaloLevel1 object
399     **/
400 mocchiut 1.10 void CaloLevel1::Clear(Option_t *t) {
401 mocchiut 1.1 //
402     istrip = 0;
403     estrip.Reset();
404     //
405     };
406    
407     /**
408     * Returns the detected energy for the given strip once loaded the event
409     **/
410     Float_t CaloLevel1::GetEstrip(Int_t sview, Int_t splane, Int_t sstrip){
411     Int_t view = -1;
412     Int_t plane = -1;
413     Int_t strip = -1;
414     Float_t mip = 0.;
415     //
416     if ( istrip == 0 ) return(0.);
417     //
418     for (Int_t i = 0; i<istrip; i++ ){
419     //
420     mip = DecodeEstrip(i,view,plane,strip);
421     //
422     if ( view == sview && splane == plane && sstrip == strip ) return(mip);
423     //
424     // entry are ordered by strip, plane and view number. Go out if you pass the input strip
425     //
426     if ( view == sview && plane == splane && strip > sstrip ) return(0.);
427     if ( view == sview && plane > splane ) return(0.);
428     if ( view > sview ) return(0.);
429     //
430     };
431     return(0.);
432     };
433    
434     /**
435     * Given estrip entry returns energy plus view, plane and strip numbers
436     **/
437     Float_t CaloLevel1::DecodeEstrip(Int_t entry, Int_t &view, Int_t &plane, Int_t &strip){
438 mocchiut 1.9 Bool_t sat = false;
439     Float_t mip=this->DecodeEstrip(entry,view,plane,strip,sat);
440     return(mip);
441     };
442    
443     /**
444     * Given estrip entry returns energy plus view, plane, strip numbers and saturation info
445     **/
446     Float_t CaloLevel1::DecodeEstrip(Int_t entry, Int_t &view, Int_t &plane, Int_t &strip, Bool_t &saturated){
447 mocchiut 1.1 //
448 mocchiut 1.9 if ( entry>istrip ){
449     //
450     printf(" ERROR: problems decoding entry %i, it seems that number of entries is %i\n",entry,istrip);
451     //
452     return(0.);
453     };
454 mocchiut 1.1 //
455     // printf(" num lim %f \n",std::numeric_limits<Float_t>::max());
456     // printf(" estrip.At(%i) = %i \n",entry,estrip.At(entry));
457     //
458     Int_t eval = 0;
459     if ( estrip.At(entry) > 0. ){
460     view = 0;
461     eval = estrip.At(entry);
462     } else {
463     view = 1;
464     eval = -estrip.At(entry);
465     };
466     //
467     Int_t fbi = 0;
468     fbi = (Int_t)truncf((Float_t)(eval/1000000000));
469     //
470     Int_t plom = 0;
471     plom = (Int_t)truncf((Float_t)((eval-fbi*1000000000)/10000000));
472     //
473     Float_t tim = 100000.;
474     plane = plom;
475     if ( fbi == 1 ) tim = 10000.;
476     if ( plom > 21 ){
477     plane = plom - 22;
478     if ( fbi == 1 ){
479     tim = 1000.;
480     } else {
481     tim = 100.;
482     };
483     };
484     if ( plom > 43 ){
485     plane = plom - 44;
486     tim = 10.;
487     };
488     if ( plom > 65 ){
489     plane = plom - 66;
490     tim = 1.;
491     };
492     //
493     strip = (Int_t)truncf((Float_t)((eval - fbi*1000000000 -plom*10000000)/100000));
494     //
495 mocchiut 1.13 Double_t mip = (Double_t)(((Float_t)(eval - fbi*1000000000 -plom*10000000 -strip*100000))/tim);
496 mocchiut 1.1 //
497 mocchiut 1.9 saturated = false;
498     if ( mip > 5000. ){
499     mip -= 5000.;
500     saturated = true;
501     };
502 mocchiut 1.13 if ( mip > 0. && mip < 99999. ) return((Float_t)mip);
503 mocchiut 1.1 //
504 mocchiut 1.9 printf(" ERROR: problems decoding value %i at entry %i \n",estrip.At(entry),entry);
505 mocchiut 1.1 //
506     view = -1;
507     plane = -1;
508     strip = -1;
509     return(0.);
510     }
511    
512 mocchiut 1.3 /*
513     * Returns energy released on plane nplane (where 0<= nplane <= 43, 0 = 1Y, 1 = 1X, 2 = 2Y, 3 = 2X, etc. etc.).
514     */
515     Float_t CaloLevel1::qtotpl(Int_t nplane){
516 mocchiut 1.9 Bool_t sat = false;
517     Float_t mip = this->qtotpl(nplane,sat);
518     return(mip);
519     };
520    
521     /*
522     * Returns energy released on plane nplane (where 0<= nplane <= 43, 0 = 1Y, 1 = 1X, 2 = 2Y, 3 = 2X, etc. etc.).
523     */
524     Float_t CaloLevel1::qtotpl(Int_t nplane, Bool_t &sat){
525 mocchiut 1.3 //
526 mocchiut 1.9 sat = false;
527 mocchiut 1.3 Int_t sview = 1;
528     if ( nplane%2 ) sview = 0;
529     //
530 mocchiut 1.11 // Int_t splane = nplane-(sview+1)/2;
531 mocchiut 1.12 Int_t splane = (nplane+sview-1)/2;
532 mocchiut 1.3 //
533 mocchiut 1.9 Float_t totmip = qtotpl(sview,splane,sat);
534 mocchiut 1.3 //
535     return(totmip);
536     //
537     };
538    
539     /*
540     * Returns energy released on view "view" (0 = X, 1 = Y) and plane "plane" ( 0 <= plane <= 21 ).
541     */
542     Float_t CaloLevel1::qtotpl(Int_t sview, Int_t splane){
543 mocchiut 1.9 Bool_t sat = false;
544     Float_t mip = this->qtotpl(sview,splane,sat);
545     return(mip);
546     };
547    
548     /*
549     * Returns energy released on view "view" (0 = X, 1 = Y) and plane "plane" ( 0 <= plane <= 21 ).
550     */
551     Float_t CaloLevel1::qtotpl(Int_t sview, Int_t splane, Bool_t &sat){
552 mocchiut 1.3 //
553     Int_t view = -1;
554     Int_t plane = -1;
555     Int_t strip = -1;
556 mocchiut 1.9 Bool_t lsat = false;
557     sat = false;
558 mocchiut 1.3 //
559     Float_t mip = 0.;
560     Float_t totmip = 0.;
561     //
562     if ( istrip == 0 ) return(0.);
563     //
564     for (Int_t i = 0; i<istrip; i++ ){
565     //
566 mocchiut 1.9 mip = DecodeEstrip(i,view,plane,strip,lsat);
567 mocchiut 1.3 //
568 mocchiut 1.9 if ( view == sview && splane == plane ){
569     if ( lsat ) sat = true;
570     totmip += mip;
571 mocchiut 1.13 //printf(" totmip %f mip %f \n",totmip,mip);
572 mocchiut 1.9 };
573 mocchiut 1.3 //
574     // entry are ordered by strip, plane and view number. Go out if you pass the input strip
575     //
576     if ( view == sview && plane > splane ) return(totmip);
577     if ( view > sview ) return(totmip);
578     //
579     };
580     //
581     return(totmip);
582     //
583     };
  ViewVC Help
Powered by ViewVC 1.1.23