/[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.3 - (hide annotations) (download)
Thu Jan 11 16:08:02 2007 UTC (17 years, 11 months ago) by mocchiut
Branch: MAIN
CVS Tags: v3r00
Changes since 1.2: +46 -0 lines
Calorimeter: added methods to get released energy on any plane

1 mocchiut 1.1 /**
2     * \file src/CaloLevel1.cpp
3     * \author Emiliano Mocchiutti
4     *
5     **/
6     #include <TObject.h>
7 mocchiut 1.2 #include <TMath.h>
8 mocchiut 1.1 #include <CaloLevel1.h>
9     //
10     ClassImp(CaloStrip);
11     ClassImp(CaloLevel1);
12    
13     /**
14     * CaloStrip default constructor
15     **/
16     CaloStrip::CaloStrip() {
17     c1 = 0;
18     this->Clear();
19     };
20    
21     /**
22     * CaloStrip default constructor
23     **/
24     CaloStrip::CaloStrip(CaloLevel1 *calo) {
25     c1 = calo->GetCaloLevel1();
26     this->Clear();
27     };
28    
29     /**
30     * Clear variables
31     **/
32     void CaloStrip::Clear() {
33     fE = 0.;
34     fX = 0.;
35     fY = 0.;
36     fZ = 0.;
37     fView = 0;
38     fPlane = 0;
39     fStrip = 0;
40     };
41    
42     /**
43     * Given a strip returns its position in the PAMELA reference system
44     **/
45     void CaloStrip::Set(Int_t view, Int_t plane, Int_t strip) {
46     //
47     this->Clear();
48     //
49     if ( view < 0 || view > 1 ){
50     printf(" ERROR: 0 =< view =< 1 \n");
51     return;
52     };
53     if ( plane < 0 || plane > 21 ){
54     printf(" ERROR: 0 =< plane =< 21 \n");
55     return;
56     };
57     if ( strip < 0 || strip > 95 ){
58     printf(" ERROR: 0 =< strip =< 95 \n");
59     return;
60     };
61     //
62     Float_t lShift = 0.;
63     Float_t lPos = 0.;
64     extern struct shift shift_;
65     //
66     // Find MIPs for that strip
67     //
68     if ( c1 ) fE = c1->GetEstrip(view, plane, strip);
69     //
70     fView = view + 1;
71     fPlane = plane + 1;
72     fStrip = strip + 1;
73     if ( fPlane%2 ){
74     lShift = -0.5;
75     } else {
76     lShift = 0.5;
77     };
78     //
79     shift_.shift = lShift;
80     //
81     Float_t zplane[22];
82     zplane[0] = 0.;
83     Int_t ii = 0;
84     for ( Int_t i = 2; i < 23; i++){
85     ii = i-1;
86     if ( i%2 ){
87     zplane[ii] = zplane[ii-1] - 10.09;
88     } else {
89     zplane[ii] = zplane[ii-1] - 8.09;
90     };
91     };
92     //
93     millim_(&fStrip,&lPos);
94     //
95     if ( fView == 1 ){
96     //
97     // X view
98     //
99     fX = (lPos - CTX)/10.;
100     fY = 0.;
101     fZ = (zplane[fPlane-1] - 5.81 + CTZ)/10.;
102     } else {
103     //
104     // Y view
105     //
106     fX = 0;
107     fY = (lPos - CTY)/10.;
108     fZ = (zplane[fPlane-1] + CTZ)/10.;
109     };
110     //
111     };
112    
113     /**
114     * Given a point in the space (PAMELA ref system) returns the closest strip
115     **/
116     void CaloStrip::Set(Float_t X, Float_t Y, Float_t Z) {
117     //
118 mocchiut 1.2 fX = X;
119     fY = Y;
120     fZ = Z;
121     //
122     Float_t zplane[22];
123     zplane[0] = 0.;
124     Int_t ii = 0;
125     for ( Int_t i = 2; i < 23; i++){
126     ii = i-1;
127     if ( i%2 ){
128     zplane[ii] = zplane[ii-1] - 10.09;
129     } else {
130     zplane[ii] = zplane[ii-1] - 8.09;
131     };
132     };
133     //
134     Float_t dzx[22];
135     Float_t dzy[22];
136     for ( Int_t i=0; i < 22; i++){
137     dzx[i] = fabs(fZ*10. - (zplane[i] - 5.81 + CTZ));
138     dzy[i] = fabs(fZ*10. - (zplane[i] + CTZ));
139     };
140     //
141     Float_t minx = TMath::MinElement(22,dzx);
142     Float_t miny = TMath::MinElement(22,dzy);
143     //
144     // find view
145     //
146     if ( minx < miny ){
147     fView = 1;
148     } else {
149     fView = 2;
150     };
151     //
152     // find plane
153     //
154     Float_t pos = 0.;
155     //
156     for ( Int_t i=0; i < 22; i++){
157     if ( fView == 1 ){
158     if ( dzx[i] == minx ){
159     fPlane = i+1;
160     pos = fX*10. + CTX;
161     };
162     } else {
163     if ( dzy[i] == miny ){
164     fPlane = i+1;
165     pos = fY*10. + CTY;
166     };
167     };
168     };
169     //
170     // find strip
171     //
172     Float_t dxy[96];
173     Float_t stpos = 0.;
174     //
175     CaloStrip *ca = new CaloStrip();
176     //
177     for ( Int_t i=0; i < 96; i++){
178     ca->Set(fView-1,fPlane-1,i);
179     if ( fView == 1 ){
180     stpos = ca->GetX()*10. + CTX;
181     } else {
182     stpos = ca->GetY()*10. + CTY;
183     };
184     dxy[i] = fabs(pos - stpos);
185     };
186     //
187     delete ca;
188 mocchiut 1.1 //
189 mocchiut 1.2 Float_t mins = TMath::MinElement(96,dxy);
190     //
191     for ( Int_t i=0; i < 96; i++){
192     if ( dxy[i] == mins ) fStrip = i+1;
193     };
194 mocchiut 1.1 };
195    
196     /**
197     * CaloLevel1 constructor
198     **/
199     CaloLevel1::CaloLevel1() {
200     //
201     estrip = TArrayI(0,NULL);
202     //
203     this->Clear();
204     //
205     };
206    
207     /**
208     * Clear the CaloLevel1 object
209     **/
210     void CaloLevel1::Clear() {
211     //
212     istrip = 0;
213     estrip.Reset();
214     //
215     };
216    
217     /**
218     * Returns the detected energy for the given strip once loaded the event
219     **/
220     Float_t CaloLevel1::GetEstrip(Int_t sview, Int_t splane, Int_t sstrip){
221     Int_t view = -1;
222     Int_t plane = -1;
223     Int_t strip = -1;
224     Float_t mip = 0.;
225     //
226     if ( istrip == 0 ) return(0.);
227     //
228     for (Int_t i = 0; i<istrip; i++ ){
229     //
230     mip = DecodeEstrip(i,view,plane,strip);
231     //
232     if ( view == sview && splane == plane && sstrip == strip ) return(mip);
233     //
234     // entry are ordered by strip, plane and view number. Go out if you pass the input strip
235     //
236     if ( view == sview && plane == splane && strip > sstrip ) return(0.);
237     if ( view == sview && plane > splane ) return(0.);
238     if ( view > sview ) return(0.);
239     //
240     };
241     return(0.);
242     };
243    
244     /**
245     * Given estrip entry returns energy plus view, plane and strip numbers
246     **/
247     Float_t CaloLevel1::DecodeEstrip(Int_t entry, Int_t &view, Int_t &plane, Int_t &strip){
248     //
249     if ( entry>istrip ) return(0.);
250     //
251     // printf(" num lim %f \n",std::numeric_limits<Float_t>::max());
252     // printf(" estrip.At(%i) = %i \n",entry,estrip.At(entry));
253     //
254     Int_t eval = 0;
255     if ( estrip.At(entry) > 0. ){
256     view = 0;
257     eval = estrip.At(entry);
258     } else {
259     view = 1;
260     eval = -estrip.At(entry);
261     };
262     //
263     Int_t fbi = 0;
264     fbi = (Int_t)truncf((Float_t)(eval/1000000000));
265     //
266     Int_t plom = 0;
267     plom = (Int_t)truncf((Float_t)((eval-fbi*1000000000)/10000000));
268     //
269     Float_t tim = 100000.;
270     plane = plom;
271     if ( fbi == 1 ) tim = 10000.;
272     if ( plom > 21 ){
273     plane = plom - 22;
274     if ( fbi == 1 ){
275     tim = 1000.;
276     } else {
277     tim = 100.;
278     };
279     };
280     if ( plom > 43 ){
281     plane = plom - 44;
282     tim = 10.;
283     };
284     if ( plom > 65 ){
285     plane = plom - 66;
286     tim = 1.;
287     };
288     //
289     strip = (Int_t)truncf((Float_t)((eval - fbi*1000000000 -plom*10000000)/100000));
290     //
291     Float_t mip = ((Float_t)(eval - fbi*1000000000 -plom*10000000 -strip*100000))/tim;
292     //
293     if ( mip > 0. && mip < 99999. ) return(mip);
294     //
295     printf(" WARNING: problems decoding value %i at entry %i \n",estrip.At(entry),entry);
296     //
297     view = -1;
298     plane = -1;
299     strip = -1;
300     return(0.);
301     }
302    
303 mocchiut 1.3 /*
304     * Returns energy released on plane nplane (where 0<= nplane <= 43, 0 = 1Y, 1 = 1X, 2 = 2Y, 3 = 2X, etc. etc.).
305     */
306     Float_t CaloLevel1::qtotpl(Int_t nplane){
307     //
308     Int_t sview = 1;
309     if ( nplane%2 ) sview = 0;
310     //
311     Int_t splane = nplane-(sview+1)/2;
312     //
313     Float_t totmip = qtotpl(sview,splane);
314     //
315     return(totmip);
316     //
317     };
318    
319     /*
320     * Returns energy released on view "view" (0 = X, 1 = Y) and plane "plane" ( 0 <= plane <= 21 ).
321     */
322     Float_t CaloLevel1::qtotpl(Int_t sview, Int_t splane){
323     //
324     Int_t view = -1;
325     Int_t plane = -1;
326     Int_t strip = -1;
327     //
328     Float_t mip = 0.;
329     Float_t totmip = 0.;
330     //
331     if ( istrip == 0 ) return(0.);
332     //
333     for (Int_t i = 0; i<istrip; i++ ){
334     //
335     mip = DecodeEstrip(i,view,plane,strip);
336     //
337     if ( view == sview && splane == plane ) totmip += mip;
338     //
339     // entry are ordered by strip, plane and view number. Go out if you pass the input strip
340     //
341     if ( view == sview && plane > splane ) return(totmip);
342     if ( view > sview ) return(totmip);
343     //
344     };
345     //
346     return(totmip);
347     //
348     };

  ViewVC Help
Powered by ViewVC 1.1.23