/[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.7 - (hide annotations) (download)
Fri Mar 30 11:22:07 2007 UTC (17 years, 8 months ago) by mocchiut
Branch: MAIN
Changes since 1.6: +1 -18 lines
Shift inversion on X views bug 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     /**
12     * CaloStrip default constructor
13     **/
14     CaloStrip::CaloStrip() {
15     c1 = 0;
16     this->Clear();
17 mocchiut 1.5 };
18    
19     /**
20     * CaloStrip default constructor
21     **/
22     CaloStrip::CaloStrip(Bool_t mechalig) {
23     c1 = 0;
24     this->Clear();
25     if ( mechalig ){
26     ismech = true;
27     UXal = MECHCTX;
28     UYal = MECHCTY;
29     UZal = MECHCTZ;
30     } else {
31     ismech = false;
32     UseStandardAlig();
33     };
34     };
35    
36     /**
37     * CaloStrip default constructor
38     **/
39     CaloStrip::CaloStrip(CaloLevel1 *calo) {
40     c1 = calo->GetCaloLevel1();
41     this->Clear();
42     ismech = false;
43 mocchiut 1.4 UseStandardAlig();
44 mocchiut 1.1 };
45    
46     /**
47     * CaloStrip default constructor
48     **/
49 mocchiut 1.5 CaloStrip::CaloStrip(CaloLevel1 *calo, Bool_t mechalig) {
50 mocchiut 1.1 c1 = calo->GetCaloLevel1();
51     this->Clear();
52 mocchiut 1.5 if ( mechalig ){
53     ismech = true;
54     UXal = MECHCTX;
55     UYal = MECHCTY;
56     UZal = MECHCTZ;
57     } else {
58     ismech = false;
59     UseStandardAlig();
60     };
61 mocchiut 1.1 };
62    
63     /**
64     * Clear variables
65     **/
66 mocchiut 1.4 void CaloStrip::Clear() {
67 mocchiut 1.1 fE = 0.;
68     fX = 0.;
69     fY = 0.;
70     fZ = 0.;
71     fView = 0;
72     fPlane = 0;
73     fStrip = 0;
74     };
75    
76     /**
77 mocchiut 1.4 * Connect to the DB and retrieve alignement parameters
78     **/
79     void CaloStrip::UseStandardAlig(){
80     //
81     ismech = false;
82     //
83     stringstream aligfile;
84     Int_t error = 0;
85     FILE *f = 0;
86     ifstream badfile;
87     GL_PARAM *glparam = new GL_PARAM();
88     //
89     TString host = "mysql://localhost/pamelaprod";
90     TString user = "anonymous";
91     TString psw = "";
92     const char *pamdbhost=gSystem->Getenv("PAM_DBHOST");
93     const char *pamdbuser=gSystem->Getenv("PAM_DBUSER");
94     const char *pamdbpsw=gSystem->Getenv("PAM_DBPSW");
95     if ( !pamdbhost ) pamdbhost = "";
96     if ( !pamdbuser ) pamdbuser = "";
97     if ( !pamdbpsw ) pamdbpsw = "";
98     if ( strcmp(pamdbhost,"") ) host = pamdbhost;
99     if ( strcmp(pamdbuser,"") ) user = pamdbuser;
100     if ( strcmp(pamdbpsw,"") ) psw = pamdbpsw;
101     TSQLServer *dbc = TSQLServer::Connect(host.Data(),user.Data(),psw.Data());
102     //
103     UXal = 0.;
104     UYal = 0.;
105     UZal = 0.;
106     //
107     if ( dbc ){
108     //
109     // determine where I can find calorimeter ADC to MIP conversion file
110     //
111     printf(" Querying DB for calorimeter parameters files...\n");
112     //
113     //
114     //
115     error = 0;
116     error = glparam->Query_GL_PARAM(1000,102,dbc);
117     if ( error >= 0 ){
118     //
119     aligfile.str("");
120     aligfile << glparam->PATH.Data() << "/";
121     aligfile << glparam->NAME.Data();
122     //
123     printf("\n Using parameter file: \n %s \n",aligfile.str().c_str());
124     f = fopen(aligfile.str().c_str(),"rb");
125     if ( f ){
126     //
127     fread(&UXal,sizeof(UXal),1,f);
128     fread(&UYal,sizeof(UYal),1,f);
129     fread(&UZal,sizeof(UZal),1,f);
130     //
131     fclose(f);
132     };
133     //
134     };
135     dbc->Close();
136     };
137     if ( !UXal ){
138     //
139     printf(" No able to query DB for calorimeter parameters files\n Using hard-coded parameters \n");
140     UXal = CTX;
141     UYal = CTY;
142     UZal = CTZ;
143     };
144     //
145     };
146    
147     /**
148 mocchiut 1.1 * Given a strip returns its position in the PAMELA reference system
149     **/
150     void CaloStrip::Set(Int_t view, Int_t plane, Int_t strip) {
151     //
152     this->Clear();
153     //
154     if ( view < 0 || view > 1 ){
155     printf(" ERROR: 0 =< view =< 1 \n");
156     return;
157     };
158     if ( plane < 0 || plane > 21 ){
159     printf(" ERROR: 0 =< plane =< 21 \n");
160     return;
161     };
162     if ( strip < 0 || strip > 95 ){
163     printf(" ERROR: 0 =< strip =< 95 \n");
164     return;
165     };
166     //
167     Float_t lShift = 0.;
168     Float_t lPos = 0.;
169     extern struct shift shift_;
170     //
171     // Find MIPs for that strip
172     //
173     if ( c1 ) fE = c1->GetEstrip(view, plane, strip);
174     //
175     fView = view + 1;
176     fPlane = plane + 1;
177     fStrip = strip + 1;
178 mocchiut 1.7 //
179 mocchiut 1.1 if ( fPlane%2 ){
180 mocchiut 1.6 lShift = +0.5;
181     } else {
182 mocchiut 1.1 lShift = -0.5;
183     };
184     //
185     shift_.shift = lShift;
186     //
187     Float_t zplane[22];
188     zplane[0] = 0.;
189     Int_t ii = 0;
190     for ( Int_t i = 2; i < 23; i++){
191     ii = i-1;
192     if ( i%2 ){
193     zplane[ii] = zplane[ii-1] - 10.09;
194     } else {
195     zplane[ii] = zplane[ii-1] - 8.09;
196     };
197     };
198     //
199     millim_(&fStrip,&lPos);
200     //
201     if ( fView == 1 ){
202     //
203     // X view
204     //
205 mocchiut 1.4 fX = (lPos - UXal)/10.;
206 mocchiut 1.1 fY = 0.;
207 mocchiut 1.4 fZ = (zplane[fPlane-1] - 5.81 + UZal)/10.;
208 mocchiut 1.6 //
209     } else {
210 mocchiut 1.1 //
211     // Y view
212     //
213     fX = 0;
214 mocchiut 1.4 fY = (lPos - UYal)/10.;
215     fZ = (zplane[fPlane-1] + UZal)/10.;
216 mocchiut 1.1 };
217     //
218     };
219    
220     /**
221     * Given a point in the space (PAMELA ref system) returns the closest strip
222     **/
223     void CaloStrip::Set(Float_t X, Float_t Y, Float_t Z) {
224     //
225 mocchiut 1.2 fX = X;
226     fY = Y;
227     fZ = Z;
228     //
229     Float_t zplane[22];
230     zplane[0] = 0.;
231     Int_t ii = 0;
232     for ( Int_t i = 2; i < 23; i++){
233     ii = i-1;
234     if ( i%2 ){
235     zplane[ii] = zplane[ii-1] - 10.09;
236     } else {
237     zplane[ii] = zplane[ii-1] - 8.09;
238     };
239     };
240     //
241     Float_t dzx[22];
242     Float_t dzy[22];
243     for ( Int_t i=0; i < 22; i++){
244 mocchiut 1.4 dzx[i] = fabs(fZ*10. - (zplane[i] - 5.81 + UZal));
245     dzy[i] = fabs(fZ*10. - (zplane[i] + UZal));
246 mocchiut 1.2 };
247     //
248     Float_t minx = TMath::MinElement(22,dzx);
249     Float_t miny = TMath::MinElement(22,dzy);
250     //
251     // find view
252     //
253     if ( minx < miny ){
254     fView = 1;
255     } else {
256     fView = 2;
257     };
258     //
259     // find plane
260     //
261     Float_t pos = 0.;
262     //
263     for ( Int_t i=0; i < 22; i++){
264     if ( fView == 1 ){
265     if ( dzx[i] == minx ){
266     fPlane = i+1;
267 mocchiut 1.4 pos = fX*10. + UXal;
268 mocchiut 1.2 };
269     } else {
270     if ( dzy[i] == miny ){
271     fPlane = i+1;
272 mocchiut 1.4 pos = fY*10. + UYal;
273 mocchiut 1.2 };
274     };
275     };
276     //
277     // find strip
278     //
279     Float_t dxy[96];
280     Float_t stpos = 0.;
281     //
282     CaloStrip *ca = new CaloStrip();
283     //
284     for ( Int_t i=0; i < 96; i++){
285     ca->Set(fView-1,fPlane-1,i);
286     if ( fView == 1 ){
287 mocchiut 1.4 stpos = ca->GetX()*10. + UXal;
288 mocchiut 1.2 } else {
289 mocchiut 1.4 stpos = ca->GetY()*10. + UYal;
290 mocchiut 1.2 };
291     dxy[i] = fabs(pos - stpos);
292     };
293     //
294     delete ca;
295 mocchiut 1.1 //
296 mocchiut 1.2 Float_t mins = TMath::MinElement(96,dxy);
297     //
298     for ( Int_t i=0; i < 96; i++){
299     if ( dxy[i] == mins ) fStrip = i+1;
300     };
301 mocchiut 1.1 };
302    
303     /**
304     * CaloLevel1 constructor
305     **/
306     CaloLevel1::CaloLevel1() {
307     //
308     estrip = TArrayI(0,NULL);
309     //
310     this->Clear();
311     //
312     };
313    
314     /**
315     * Clear the CaloLevel1 object
316     **/
317     void CaloLevel1::Clear() {
318     //
319     istrip = 0;
320     estrip.Reset();
321     //
322     };
323    
324     /**
325     * Returns the detected energy for the given strip once loaded the event
326     **/
327     Float_t CaloLevel1::GetEstrip(Int_t sview, Int_t splane, Int_t sstrip){
328     Int_t view = -1;
329     Int_t plane = -1;
330     Int_t strip = -1;
331     Float_t mip = 0.;
332     //
333     if ( istrip == 0 ) return(0.);
334     //
335     for (Int_t i = 0; i<istrip; i++ ){
336     //
337     mip = DecodeEstrip(i,view,plane,strip);
338     //
339     if ( view == sview && splane == plane && sstrip == strip ) return(mip);
340     //
341     // entry are ordered by strip, plane and view number. Go out if you pass the input strip
342     //
343     if ( view == sview && plane == splane && strip > sstrip ) return(0.);
344     if ( view == sview && plane > splane ) return(0.);
345     if ( view > sview ) return(0.);
346     //
347     };
348     return(0.);
349     };
350    
351     /**
352     * Given estrip entry returns energy plus view, plane and strip numbers
353     **/
354     Float_t CaloLevel1::DecodeEstrip(Int_t entry, Int_t &view, Int_t &plane, Int_t &strip){
355     //
356     if ( entry>istrip ) return(0.);
357     //
358     // printf(" num lim %f \n",std::numeric_limits<Float_t>::max());
359     // printf(" estrip.At(%i) = %i \n",entry,estrip.At(entry));
360     //
361     Int_t eval = 0;
362     if ( estrip.At(entry) > 0. ){
363     view = 0;
364     eval = estrip.At(entry);
365     } else {
366     view = 1;
367     eval = -estrip.At(entry);
368     };
369     //
370     Int_t fbi = 0;
371     fbi = (Int_t)truncf((Float_t)(eval/1000000000));
372     //
373     Int_t plom = 0;
374     plom = (Int_t)truncf((Float_t)((eval-fbi*1000000000)/10000000));
375     //
376     Float_t tim = 100000.;
377     plane = plom;
378     if ( fbi == 1 ) tim = 10000.;
379     if ( plom > 21 ){
380     plane = plom - 22;
381     if ( fbi == 1 ){
382     tim = 1000.;
383     } else {
384     tim = 100.;
385     };
386     };
387     if ( plom > 43 ){
388     plane = plom - 44;
389     tim = 10.;
390     };
391     if ( plom > 65 ){
392     plane = plom - 66;
393     tim = 1.;
394     };
395     //
396     strip = (Int_t)truncf((Float_t)((eval - fbi*1000000000 -plom*10000000)/100000));
397     //
398     Float_t mip = ((Float_t)(eval - fbi*1000000000 -plom*10000000 -strip*100000))/tim;
399     //
400     if ( mip > 0. && mip < 99999. ) return(mip);
401     //
402     printf(" WARNING: problems decoding value %i at entry %i \n",estrip.At(entry),entry);
403     //
404     view = -1;
405     plane = -1;
406     strip = -1;
407     return(0.);
408     }
409    
410 mocchiut 1.3 /*
411     * Returns energy released on plane nplane (where 0<= nplane <= 43, 0 = 1Y, 1 = 1X, 2 = 2Y, 3 = 2X, etc. etc.).
412     */
413     Float_t CaloLevel1::qtotpl(Int_t nplane){
414     //
415     Int_t sview = 1;
416     if ( nplane%2 ) sview = 0;
417     //
418     Int_t splane = nplane-(sview+1)/2;
419     //
420     Float_t totmip = qtotpl(sview,splane);
421     //
422     return(totmip);
423     //
424     };
425    
426     /*
427     * Returns energy released on view "view" (0 = X, 1 = Y) and plane "plane" ( 0 <= plane <= 21 ).
428     */
429     Float_t CaloLevel1::qtotpl(Int_t sview, Int_t splane){
430     //
431     Int_t view = -1;
432     Int_t plane = -1;
433     Int_t strip = -1;
434     //
435     Float_t mip = 0.;
436     Float_t totmip = 0.;
437     //
438     if ( istrip == 0 ) return(0.);
439     //
440     for (Int_t i = 0; i<istrip; i++ ){
441     //
442     mip = DecodeEstrip(i,view,plane,strip);
443     //
444     if ( view == sview && splane == plane ) totmip += mip;
445     //
446     // entry are ordered by strip, plane and view number. Go out if you pass the input strip
447     //
448     if ( view == sview && plane > splane ) return(totmip);
449     if ( view > sview ) return(totmip);
450     //
451     };
452     //
453     return(totmip);
454     //
455     };

  ViewVC Help
Powered by ViewVC 1.1.23