/[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.4 - (hide annotations) (download)
Mon Mar 26 14:02:06 2007 UTC (17 years, 8 months ago) by mocchiut
Branch: MAIN
Changes since 1.3: +91 -13 lines
Bug in CaloStrip/Level0 fixed + new alignement parameters added

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

  ViewVC Help
Powered by ViewVC 1.1.23