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

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

Parent Directory Parent Directory | Revision Log Revision Log


Revision 1.14 - (show annotations) (download)
Mon May 12 14:36:09 2008 UTC (16 years, 7 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 /**
2 * \file src/CaloLevel1.cpp
3 * \author Emiliano Mocchiutti
4 *
5 **/
6 #include <CaloLevel1.h>
7 //
8 ClassImp(CaloStrip);
9 ClassImp(CaloLevel1);
10
11 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 /**
17 * CaloStrip default constructor
18 **/
19 CaloStrip::CaloStrip() {
20 c1 = 0;
21 this->Clear();
22 };
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 paramload = true;
33 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 UseStandardAlig();
50 };
51
52 /**
53 * CaloStrip default constructor
54 **/
55 CaloStrip::CaloStrip(CaloLevel1 *calo, Bool_t mechalig) {
56 c1 = calo->GetCaloLevel1();
57 this->Clear();
58 if ( mechalig ){
59 ismech = true;
60 paramload = true;
61 UXal = MECHCTX;
62 UYal = MECHCTY;
63 UZal = MECHCTZ;
64 } else {
65 ismech = false;
66 UseStandardAlig();
67 };
68 };
69
70 /**
71 * Clear variables
72 **/
73 void CaloStrip::Clear(Option_t *t) {
74 fE = 0.;
75 fX = 0.;
76 fY = 0.;
77 fZ = 0.;
78 fView = 0;
79 fPlane = 0;
80 fStrip = 0;
81 };
82
83 /**
84 * Connect to the DB and retrieve alignement parameters
85 **/
86 void CaloStrip::UseStandardAlig(){
87 //
88 if ( !paramload ){
89 //
90 paramload = true;
91 ismech = false;
92 //
93 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 //
117 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 //
123 //
124 //
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 //
133 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 //
144 };
145 dbc->Close();
146 };
147 if ( !UXal ){
148 //
149 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 };
154 //
155 };
156 //
157 };
158
159 /**
160 * 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 //
191 if ( fPlane%2 ){
192 lShift = +0.5;
193 } else {
194 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 fX = (lPos - UXal)/10.;
218 fY = 0.;
219 fZ = (zplane[fPlane-1] - 5.81 + UZal)/10.;
220 //
221 } else {
222 //
223 // Y view
224 //
225 fX = 0;
226 fY = (lPos - UYal)/10.;
227 fZ = (zplane[fPlane-1] + UZal)/10.;
228 };
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 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 dzx[i] = fabs(fZ*10. - (zplane[i] - 5.81 + UZal));
257 dzy[i] = fabs(fZ*10. - (zplane[i] + UZal));
258 };
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 pos = fX*10. + UXal;
280 };
281 } else {
282 if ( dzy[i] == miny ){
283 fPlane = i+1;
284 pos = fY*10. + UYal;
285 };
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 stpos = ca->GetX()*10. + UXal;
300 } else {
301 stpos = ca->GetY()*10. + UYal;
302 };
303 dxy[i] = fabs(pos - stpos);
304 };
305 //
306 delete ca;
307 //
308 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 };
314
315 /**
316 * 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 * 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 void CaloLevel1::Clear(Option_t *t) {
399 //
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 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 //
446 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 //
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 Double_t mip = (Double_t)(((Float_t)(eval - fbi*1000000000 -plom*10000000 -strip*100000))/tim);
494 //
495 saturated = false;
496 if ( mip > 5000. ){
497 mip -= 5000.;
498 saturated = true;
499 };
500 if ( mip > 0. && mip < 99999. ) return((Float_t)mip);
501 //
502 printf(" ERROR: problems decoding value %i at entry %i \n",estrip.At(entry),entry);
503 //
504 view = -1;
505 plane = -1;
506 strip = -1;
507 return(0.);
508 }
509
510 /*
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 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 //
524 sat = false;
525 Int_t sview = 1;
526 if ( nplane%2 ) sview = 0;
527 //
528 // Int_t splane = nplane-(sview+1)/2;
529 Int_t splane = (nplane+sview-1)/2;
530 //
531 Float_t totmip = qtotpl(sview,splane,sat);
532 //
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 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 //
551 Int_t view = -1;
552 Int_t plane = -1;
553 Int_t strip = -1;
554 Bool_t lsat = false;
555 sat = false;
556 //
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 mip = DecodeEstrip(i,view,plane,strip,lsat);
565 //
566 if ( view == sview && splane == plane ){
567 if ( lsat ) sat = true;
568 totmip += mip;
569 //printf(" totmip %f mip %f \n",totmip,mip);
570 };
571 //
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