/[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.16 - (show annotations) (download)
Fri Jun 6 14:18:16 2014 UTC (10 years, 5 months ago) by mocchiut
Branch: MAIN
CVS Tags: v10REDr01, v10RED, HEAD
Changes since 1.15: +6 -2 lines
New extended algorithm track-related added to CALO, TOF and ORB

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

  ViewVC Help
Powered by ViewVC 1.1.23