/[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.15 - (show 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 /**
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 delete dbc;
147 dbc = 0;
148 };
149 if ( !UXal ){
150 //
151 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 };
156 //
157 };
158 //
159 };
160
161 /**
162 * 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 //
193 if ( fPlane%2 ){
194 lShift = +0.5;
195 } else {
196 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 fX = (lPos - UXal)/10.;
220 fY = 0.;
221 fZ = (zplane[fPlane-1] - 5.81 + UZal)/10.;
222 //
223 } else {
224 //
225 // Y view
226 //
227 fX = 0;
228 fY = (lPos - UYal)/10.;
229 fZ = (zplane[fPlane-1] + UZal)/10.;
230 };
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 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 dzx[i] = fabs(fZ*10. - (zplane[i] - 5.81 + UZal));
259 dzy[i] = fabs(fZ*10. - (zplane[i] + UZal));
260 };
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 pos = fX*10. + UXal;
282 };
283 } else {
284 if ( dzy[i] == miny ){
285 fPlane = i+1;
286 pos = fY*10. + UYal;
287 };
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 stpos = ca->GetX()*10. + UXal;
302 } else {
303 stpos = ca->GetY()*10. + UYal;
304 };
305 dxy[i] = fabs(pos - stpos);
306 };
307 //
308 delete ca;
309 //
310 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 };
316
317 /**
318 * 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 * 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 void CaloLevel1::Clear(Option_t *t) {
401 //
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 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 //
448 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 //
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 Double_t mip = (Double_t)(((Float_t)(eval - fbi*1000000000 -plom*10000000 -strip*100000))/tim);
496 //
497 saturated = false;
498 if ( mip > 5000. ){
499 mip -= 5000.;
500 saturated = true;
501 };
502 if ( mip > 0. && mip < 99999. ) return((Float_t)mip);
503 //
504 printf(" ERROR: problems decoding value %i at entry %i \n",estrip.At(entry),entry);
505 //
506 view = -1;
507 plane = -1;
508 strip = -1;
509 return(0.);
510 }
511
512 /*
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 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 //
526 sat = false;
527 Int_t sview = 1;
528 if ( nplane%2 ) sview = 0;
529 //
530 // Int_t splane = nplane-(sview+1)/2;
531 Int_t splane = (nplane+sview-1)/2;
532 //
533 Float_t totmip = qtotpl(sview,splane,sat);
534 //
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 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 //
553 Int_t view = -1;
554 Int_t plane = -1;
555 Int_t strip = -1;
556 Bool_t lsat = false;
557 sat = false;
558 //
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 mip = DecodeEstrip(i,view,plane,strip,lsat);
567 //
568 if ( view == sview && splane == plane ){
569 if ( lsat ) sat = true;
570 totmip += mip;
571 //printf(" totmip %f mip %f \n",totmip,mip);
572 };
573 //
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