/[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.14 - (hide annotations) (download)
Mon May 12 14:36:09 2008 UTC (17 years, 5 months ago) by mocchiut
Branch: MAIN
CVS Tags: v6r01, v6r00, v9r00, v9r01
Changes since 1.13: +69 -0 lines 
New method in CaloStrip + cross-talk bugs fixed
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     };
147     if ( !UXal ){
148 mocchiut 1.4 //
149 mocchiut 1.8 printf(" No able to query DB for calorimeter parameters files\n Using hard-coded parameters \n");
150     UXal = CTX;
151     UYal = CTY;
152     UZal = CTZ;
153 mocchiut 1.4 };
154 mocchiut 1.8 //
155 mocchiut 1.4 };
156     //
157     };
158    
159     /**
160 mocchiut 1.1 * Given a strip returns its position in the PAMELA reference system
161     **/
162     void CaloStrip::Set(Int_t view, Int_t plane, Int_t strip) {
163     //
164     this->Clear();
165     //
166     if ( view < 0 || view > 1 ){
167     printf(" ERROR: 0 =< view =< 1 \n");
168     return;
169     };
170     if ( plane < 0 || plane > 21 ){
171     printf(" ERROR: 0 =< plane =< 21 \n");
172     return;
173     };
174     if ( strip < 0 || strip > 95 ){
175     printf(" ERROR: 0 =< strip =< 95 \n");
176     return;
177     };
178     //
179     Float_t lShift = 0.;
180     Float_t lPos = 0.;
181     extern struct shift shift_;
182     //
183     // Find MIPs for that strip
184     //
185     if ( c1 ) fE = c1->GetEstrip(view, plane, strip);
186     //
187     fView = view + 1;
188     fPlane = plane + 1;
189     fStrip = strip + 1;
190 mocchiut 1.7 //
191 mocchiut 1.1 if ( fPlane%2 ){
192 mocchiut 1.6 lShift = +0.5;
193     } else {
194 mocchiut 1.1 lShift = -0.5;
195     };
196     //
197     shift_.shift = lShift;
198     //
199     Float_t zplane[22];
200     zplane[0] = 0.;
201     Int_t ii = 0;
202     for ( Int_t i = 2; i < 23; i++){
203     ii = i-1;
204     if ( i%2 ){
205     zplane[ii] = zplane[ii-1] - 10.09;
206     } else {
207     zplane[ii] = zplane[ii-1] - 8.09;
208     };
209     };
210     //
211     millim_(&fStrip,&lPos);
212     //
213     if ( fView == 1 ){
214     //
215     // X view
216     //
217 mocchiut 1.4 fX = (lPos - UXal)/10.;
218 mocchiut 1.1 fY = 0.;
219 mocchiut 1.4 fZ = (zplane[fPlane-1] - 5.81 + UZal)/10.;
220 mocchiut 1.6 //
221     } else {
222 mocchiut 1.1 //
223     // Y view
224     //
225     fX = 0;
226 mocchiut 1.4 fY = (lPos - UYal)/10.;
227     fZ = (zplane[fPlane-1] + UZal)/10.;
228 mocchiut 1.1 };
229     //
230     };
231    
232     /**
233     * Given a point in the space (PAMELA ref system) returns the closest strip
234     **/
235     void CaloStrip::Set(Float_t X, Float_t Y, Float_t Z) {
236     //
237 mocchiut 1.2 fX = X;
238     fY = Y;
239     fZ = Z;
240     //
241     Float_t zplane[22];
242     zplane[0] = 0.;
243     Int_t ii = 0;
244     for ( Int_t i = 2; i < 23; i++){
245     ii = i-1;
246     if ( i%2 ){
247     zplane[ii] = zplane[ii-1] - 10.09;
248     } else {
249     zplane[ii] = zplane[ii-1] - 8.09;
250     };
251     };
252     //
253     Float_t dzx[22];
254     Float_t dzy[22];
255     for ( Int_t i=0; i < 22; i++){
256 mocchiut 1.4 dzx[i] = fabs(fZ*10. - (zplane[i] - 5.81 + UZal));
257     dzy[i] = fabs(fZ*10. - (zplane[i] + UZal));
258 mocchiut 1.2 };
259     //
260     Float_t minx = TMath::MinElement(22,dzx);
261     Float_t miny = TMath::MinElement(22,dzy);
262     //
263     // find view
264     //
265     if ( minx < miny ){
266     fView = 1;
267     } else {
268     fView = 2;
269     };
270     //
271     // find plane
272     //
273     Float_t pos = 0.;
274     //
275     for ( Int_t i=0; i < 22; i++){
276     if ( fView == 1 ){
277     if ( dzx[i] == minx ){
278     fPlane = i+1;
279 mocchiut 1.4 pos = fX*10. + UXal;
280 mocchiut 1.2 };
281     } else {
282     if ( dzy[i] == miny ){
283     fPlane = i+1;
284 mocchiut 1.4 pos = fY*10. + UYal;
285 mocchiut 1.2 };
286     };
287     };
288     //
289     // find strip
290     //
291     Float_t dxy[96];
292     Float_t stpos = 0.;
293     //
294     CaloStrip *ca = new CaloStrip();
295     //
296     for ( Int_t i=0; i < 96; i++){
297     ca->Set(fView-1,fPlane-1,i);
298     if ( fView == 1 ){
299 mocchiut 1.4 stpos = ca->GetX()*10. + UXal;
300 mocchiut 1.2 } else {
301 mocchiut 1.4 stpos = ca->GetY()*10. + UYal;
302 mocchiut 1.2 };
303     dxy[i] = fabs(pos - stpos);
304     };
305     //
306     delete ca;
307 mocchiut 1.1 //
308 mocchiut 1.2 Float_t mins = TMath::MinElement(96,dxy);
309     //
310     for ( Int_t i=0; i < 96; i++){
311     if ( dxy[i] == mins ) fStrip = i+1;
312     };
313 mocchiut 1.1 };
314    
315     /**
316 mocchiut 1.14 * Given a point in the space or a strip it returns the Silicon sensor number. Numbering goes like this:
317     *
318     * y ^
319     * |
320     * | 6 7 8
321     * | 3 4 5
322     * | 0 1 2
323     * | -----------> x
324     *
325     **/
326     Int_t CaloStrip::GetSiSensor() {
327     //
328     // fX fY fZ // fView fPlane
329     //
330     Float_t deadsi = 0.096;
331     Float_t dead = 0.05;
332     Float_t sidim = 8.00;
333     // Float_t stripdim = 0.244;
334     Float_t sensitarea = 7.808;
335     //
336     Float_t xoffset = 0.;
337     Float_t yoffset = 0.;
338     //
339     if ( (fView-1) == 0 && !((fPlane-1)%2) ){
340     xoffset = +0.05;
341     yoffset = 0.0;
342     };
343     if ( (fView-1) == 0 && ((fPlane-1)%2) ){
344     xoffset = -0.05;
345     yoffset = -0.20;
346     };
347     if ( (fView-1) == 1 && !((fPlane-1)%2) ){
348     xoffset = +0.10;
349     yoffset = -0.05;
350     };
351     if ( (fView-1) == 1 && ((fPlane-1)%2) ){
352     xoffset = -0.10;
353     yoffset = -0.15;
354     };
355     //
356     Int_t iind = -1;
357     Int_t jind = -1;
358     //
359     for (Int_t i = 0; i < 3; i++){
360     if ( (fX+xoffset+12.1) >= (deadsi+(sidim+dead)*i) && (fX+xoffset+12.1) <= (sensitarea+deadsi+(sidim+dead)*i) ){
361     iind = i;
362     break;
363     };
364     };
365     //
366     for (Int_t j = 0; j < 3; j++){
367     if ( (fY+yoffset+12.1) >= (deadsi+(sidim+dead)*j) && (fY+yoffset+12.1) <= (sensitarea+deadsi+(sidim+dead)*j) ){
368     jind = j;
369     break;
370     };
371     };
372     //
373     Int_t sensor = -1;
374     if ( iind != -1 && jind != -1 ){
375     sensor = iind + jind * 3;
376     };
377     //
378     // 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);
379     //
380     return(sensor);
381     //
382     };
383    
384     /**
385 mocchiut 1.1 * CaloLevel1 constructor
386     **/
387     CaloLevel1::CaloLevel1() {
388     //
389     estrip = TArrayI(0,NULL);
390     //
391     this->Clear();
392     //
393     };
394    
395     /**
396     * Clear the CaloLevel1 object
397     **/
398 mocchiut 1.10 void CaloLevel1::Clear(Option_t *t) {
399 mocchiut 1.1 //
400     istrip = 0;
401     estrip.Reset();
402     //
403     };
404    
405     /**
406     * Returns the detected energy for the given strip once loaded the event
407     **/
408     Float_t CaloLevel1::GetEstrip(Int_t sview, Int_t splane, Int_t sstrip){
409     Int_t view = -1;
410     Int_t plane = -1;
411     Int_t strip = -1;
412     Float_t mip = 0.;
413     //
414     if ( istrip == 0 ) return(0.);
415     //
416     for (Int_t i = 0; i<istrip; i++ ){
417     //
418     mip = DecodeEstrip(i,view,plane,strip);
419     //
420     if ( view == sview && splane == plane && sstrip == strip ) return(mip);
421     //
422     // entry are ordered by strip, plane and view number. Go out if you pass the input strip
423     //
424     if ( view == sview && plane == splane && strip > sstrip ) return(0.);
425     if ( view == sview && plane > splane ) return(0.);
426     if ( view > sview ) return(0.);
427     //
428     };
429     return(0.);
430     };
431    
432     /**
433     * Given estrip entry returns energy plus view, plane and strip numbers
434     **/
435     Float_t CaloLevel1::DecodeEstrip(Int_t entry, Int_t &view, Int_t &plane, Int_t &strip){
436 mocchiut 1.9 Bool_t sat = false;
437     Float_t mip=this->DecodeEstrip(entry,view,plane,strip,sat);
438     return(mip);
439     };
440    
441     /**
442     * Given estrip entry returns energy plus view, plane, strip numbers and saturation info
443     **/
444     Float_t CaloLevel1::DecodeEstrip(Int_t entry, Int_t &view, Int_t &plane, Int_t &strip, Bool_t &saturated){
445 mocchiut 1.1 //
446 mocchiut 1.9 if ( entry>istrip ){
447     //
448     printf(" ERROR: problems decoding entry %i, it seems that number of entries is %i\n",entry,istrip);
449     //
450     return(0.);
451     };
452 mocchiut 1.1 //
453     // printf(" num lim %f \n",std::numeric_limits<Float_t>::max());
454     // printf(" estrip.At(%i) = %i \n",entry,estrip.At(entry));
455     //
456     Int_t eval = 0;
457     if ( estrip.At(entry) > 0. ){
458     view = 0;
459     eval = estrip.At(entry);
460     } else {
461     view = 1;
462     eval = -estrip.At(entry);
463     };
464     //
465     Int_t fbi = 0;
466     fbi = (Int_t)truncf((Float_t)(eval/1000000000));
467     //
468     Int_t plom = 0;
469     plom = (Int_t)truncf((Float_t)((eval-fbi*1000000000)/10000000));
470     //
471     Float_t tim = 100000.;
472     plane = plom;
473     if ( fbi == 1 ) tim = 10000.;
474     if ( plom > 21 ){
475     plane = plom - 22;
476     if ( fbi == 1 ){
477     tim = 1000.;
478     } else {
479     tim = 100.;
480     };
481     };
482     if ( plom > 43 ){
483     plane = plom - 44;
484     tim = 10.;
485     };
486     if ( plom > 65 ){
487     plane = plom - 66;
488     tim = 1.;
489     };
490     //
491     strip = (Int_t)truncf((Float_t)((eval - fbi*1000000000 -plom*10000000)/100000));
492     //
493 mocchiut 1.13 Double_t mip = (Double_t)(((Float_t)(eval - fbi*1000000000 -plom*10000000 -strip*100000))/tim);
494 mocchiut 1.1 //
495 mocchiut 1.9 saturated = false;
496     if ( mip > 5000. ){
497     mip -= 5000.;
498     saturated = true;
499     };
500 mocchiut 1.13 if ( mip > 0. && mip < 99999. ) return((Float_t)mip);
501 mocchiut 1.1 //
502 mocchiut 1.9 printf(" ERROR: problems decoding value %i at entry %i \n",estrip.At(entry),entry);
503 mocchiut 1.1 //
504     view = -1;
505     plane = -1;
506     strip = -1;
507     return(0.);
508     }
509    
510 mocchiut 1.3 /*
511     * Returns energy released on plane nplane (where 0<= nplane <= 43, 0 = 1Y, 1 = 1X, 2 = 2Y, 3 = 2X, etc. etc.).
512     */
513     Float_t CaloLevel1::qtotpl(Int_t nplane){
514 mocchiut 1.9 Bool_t sat = false;
515     Float_t mip = this->qtotpl(nplane,sat);
516     return(mip);
517     };
518    
519     /*
520     * Returns energy released on plane nplane (where 0<= nplane <= 43, 0 = 1Y, 1 = 1X, 2 = 2Y, 3 = 2X, etc. etc.).
521     */
522     Float_t CaloLevel1::qtotpl(Int_t nplane, Bool_t &sat){
523 mocchiut 1.3 //
524 mocchiut 1.9 sat = false;
525 mocchiut 1.3 Int_t sview = 1;
526     if ( nplane%2 ) sview = 0;
527     //
528 mocchiut 1.11 // Int_t splane = nplane-(sview+1)/2;
529 mocchiut 1.12 Int_t splane = (nplane+sview-1)/2;
530 mocchiut 1.3 //
531 mocchiut 1.9 Float_t totmip = qtotpl(sview,splane,sat);
532 mocchiut 1.3 //
533     return(totmip);
534     //
535     };
536    
537     /*
538     * Returns energy released on view "view" (0 = X, 1 = Y) and plane "plane" ( 0 <= plane <= 21 ).
539     */
540     Float_t CaloLevel1::qtotpl(Int_t sview, Int_t splane){
541 mocchiut 1.9 Bool_t sat = false;
542     Float_t mip = this->qtotpl(sview,splane,sat);
543     return(mip);
544     };
545    
546     /*
547     * Returns energy released on view "view" (0 = X, 1 = Y) and plane "plane" ( 0 <= plane <= 21 ).
548     */
549     Float_t CaloLevel1::qtotpl(Int_t sview, Int_t splane, Bool_t &sat){
550 mocchiut 1.3 //
551     Int_t view = -1;
552     Int_t plane = -1;
553     Int_t strip = -1;
554 mocchiut 1.9 Bool_t lsat = false;
555     sat = false;
556 mocchiut 1.3 //
557     Float_t mip = 0.;
558     Float_t totmip = 0.;
559     //
560     if ( istrip == 0 ) return(0.);
561     //
562     for (Int_t i = 0; i<istrip; i++ ){
563     //
564 mocchiut 1.9 mip = DecodeEstrip(i,view,plane,strip,lsat);
565 mocchiut 1.3 //
566 mocchiut 1.9 if ( view == sview && splane == plane ){
567     if ( lsat ) sat = true;
568     totmip += mip;
569 mocchiut 1.13 //printf(" totmip %f mip %f \n",totmip,mip);
570 mocchiut 1.9 };
571 mocchiut 1.3 //
572     // entry are ordered by strip, plane and view number. Go out if you pass the input strip
573     //
574     if ( view == sview && plane > splane ) return(totmip);
575     if ( view > sview ) return(totmip);
576     //
577     };
578     //
579     return(totmip);
580     //
581     };
  ViewVC Help
Powered by ViewVC 1.1.23