/[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.16 - (hide 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 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 mocchiut 1.16 debug = false;
22 mocchiut 1.1 this->Clear();
23 mocchiut 1.5 };
24    
25     /**
26     * CaloStrip default constructor
27     **/
28     CaloStrip::CaloStrip(Bool_t mechalig) {
29     c1 = 0;
30 mocchiut 1.16 debug = false;
31 mocchiut 1.5 this->Clear();
32     if ( mechalig ){
33     ismech = true;
34 mocchiut 1.8 paramload = true;
35 mocchiut 1.5 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 mocchiut 1.16 debug = false;
50 mocchiut 1.5 this->Clear();
51     ismech = false;
52 mocchiut 1.4 UseStandardAlig();
53 mocchiut 1.1 };
54    
55     /**
56     * CaloStrip default constructor
57     **/
58 mocchiut 1.5 CaloStrip::CaloStrip(CaloLevel1 *calo, Bool_t mechalig) {
59 mocchiut 1.1 c1 = calo->GetCaloLevel1();
60 mocchiut 1.16 debug = false;
61 mocchiut 1.1 this->Clear();
62 mocchiut 1.5 if ( mechalig ){
63     ismech = true;
64 mocchiut 1.8 paramload = true;
65 mocchiut 1.5 UXal = MECHCTX;
66     UYal = MECHCTY;
67     UZal = MECHCTZ;
68     } else {
69     ismech = false;
70     UseStandardAlig();
71     };
72 mocchiut 1.1 };
73    
74     /**
75     * Clear variables
76     **/
77 mocchiut 1.10 void CaloStrip::Clear(Option_t *t) {
78 mocchiut 1.1 fE = 0.;
79     fX = 0.;
80     fY = 0.;
81     fZ = 0.;
82     fView = 0;
83     fPlane = 0;
84     fStrip = 0;
85     };
86    
87     /**
88 mocchiut 1.4 * Connect to the DB and retrieve alignement parameters
89     **/
90     void CaloStrip::UseStandardAlig(){
91     //
92 mocchiut 1.8 if ( !paramload ){
93 mocchiut 1.4 //
94 mocchiut 1.8 paramload = true;
95     ismech = false;
96 mocchiut 1.4 //
97 mocchiut 1.8 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 mocchiut 1.4 //
121 mocchiut 1.8 if ( dbc ){
122     //
123     // determine where I can find calorimeter ADC to MIP conversion file
124     //
125 mocchiut 1.16 if ( debug ) printf(" Querying DB for calorimeter parameters files...\n");
126 mocchiut 1.4 //
127     //
128 mocchiut 1.8 //
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 mocchiut 1.4 //
137 mocchiut 1.16 if (debug) printf("\n Using parameter file: \n %s \n",aligfile.str().c_str());
138 mocchiut 1.8 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 mocchiut 1.4 //
148     };
149 mocchiut 1.8 dbc->Close();
150 pam-fi 1.15 delete dbc;
151     dbc = 0;
152 mocchiut 1.8 };
153     if ( !UXal ){
154 mocchiut 1.4 //
155 mocchiut 1.8 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 mocchiut 1.4 };
160 mocchiut 1.8 //
161 mocchiut 1.4 };
162     //
163     };
164    
165     /**
166 mocchiut 1.1 * 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 mocchiut 1.7 //
197 mocchiut 1.1 if ( fPlane%2 ){
198 mocchiut 1.6 lShift = +0.5;
199     } else {
200 mocchiut 1.1 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 mocchiut 1.4 fX = (lPos - UXal)/10.;
224 mocchiut 1.1 fY = 0.;
225 mocchiut 1.4 fZ = (zplane[fPlane-1] - 5.81 + UZal)/10.;
226 mocchiut 1.6 //
227     } else {
228 mocchiut 1.1 //
229     // Y view
230     //
231     fX = 0;
232 mocchiut 1.4 fY = (lPos - UYal)/10.;
233     fZ = (zplane[fPlane-1] + UZal)/10.;
234 mocchiut 1.1 };
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 mocchiut 1.2 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 mocchiut 1.4 dzx[i] = fabs(fZ*10. - (zplane[i] - 5.81 + UZal));
263     dzy[i] = fabs(fZ*10. - (zplane[i] + UZal));
264 mocchiut 1.2 };
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 mocchiut 1.4 pos = fX*10. + UXal;
286 mocchiut 1.2 };
287     } else {
288     if ( dzy[i] == miny ){
289     fPlane = i+1;
290 mocchiut 1.4 pos = fY*10. + UYal;
291 mocchiut 1.2 };
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 mocchiut 1.4 stpos = ca->GetX()*10. + UXal;
306 mocchiut 1.2 } else {
307 mocchiut 1.4 stpos = ca->GetY()*10. + UYal;
308 mocchiut 1.2 };
309     dxy[i] = fabs(pos - stpos);
310     };
311     //
312     delete ca;
313 mocchiut 1.1 //
314 mocchiut 1.2 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 mocchiut 1.1 };
320    
321     /**
322 mocchiut 1.14 * 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 mocchiut 1.1 * 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 mocchiut 1.10 void CaloLevel1::Clear(Option_t *t) {
405 mocchiut 1.1 //
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 mocchiut 1.9 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 mocchiut 1.1 //
452 mocchiut 1.9 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 mocchiut 1.1 //
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 mocchiut 1.13 Double_t mip = (Double_t)(((Float_t)(eval - fbi*1000000000 -plom*10000000 -strip*100000))/tim);
500 mocchiut 1.1 //
501 mocchiut 1.9 saturated = false;
502     if ( mip > 5000. ){
503     mip -= 5000.;
504     saturated = true;
505     };
506 mocchiut 1.13 if ( mip > 0. && mip < 99999. ) return((Float_t)mip);
507 mocchiut 1.1 //
508 mocchiut 1.9 printf(" ERROR: problems decoding value %i at entry %i \n",estrip.At(entry),entry);
509 mocchiut 1.1 //
510     view = -1;
511     plane = -1;
512     strip = -1;
513     return(0.);
514     }
515    
516 mocchiut 1.3 /*
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 mocchiut 1.9 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 mocchiut 1.3 //
530 mocchiut 1.9 sat = false;
531 mocchiut 1.3 Int_t sview = 1;
532     if ( nplane%2 ) sview = 0;
533     //
534 mocchiut 1.11 // Int_t splane = nplane-(sview+1)/2;
535 mocchiut 1.12 Int_t splane = (nplane+sview-1)/2;
536 mocchiut 1.3 //
537 mocchiut 1.9 Float_t totmip = qtotpl(sview,splane,sat);
538 mocchiut 1.3 //
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 mocchiut 1.9 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 mocchiut 1.3 //
557     Int_t view = -1;
558     Int_t plane = -1;
559     Int_t strip = -1;
560 mocchiut 1.9 Bool_t lsat = false;
561     sat = false;
562 mocchiut 1.3 //
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 mocchiut 1.9 mip = DecodeEstrip(i,view,plane,strip,lsat);
571 mocchiut 1.3 //
572 mocchiut 1.9 if ( view == sview && splane == plane ){
573     if ( lsat ) sat = true;
574     totmip += mip;
575 mocchiut 1.13 //printf(" totmip %f mip %f \n",totmip,mip);
576 mocchiut 1.9 };
577 mocchiut 1.3 //
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