Code from MUON-dev joined
[u/mrichter/AliRoot.git] / MUON / AliMUONSegmentationV01.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /*
17 $Log$
18 Revision 1.1.2.1  2000/06/09 21:37:30  morsch
19 AliMUONSegmentationV01 code  from  AliMUONSegResV01.cxx
20
21 */
22
23
24 /////////////////////////////////////////////////////
25 //  Segmentation and Response classes version 01   //
26 /////////////////////////////////////////////////////
27
28 #include <TBox.h> 
29 #include <TF1.h> 
30 #include <TObjArray.h>
31 #include <iostream.h>
32
33 #include "AliMUONSegmentationV01.h"
34 #include "AliMUON.h"
35
36
37
38 //___________________________________________
39 ClassImp(AliMUONSegmentationV01)
40
41 AliMUONSegmentationV01::AliMUONSegmentationV01(const AliMUONSegmentationV01& segmentation)
42 {
43 // Dummy copy constructor
44 }
45 AliMUONSegmentationV01::AliMUONSegmentationV01() 
46 {
47 // Default constructor
48     fNsec=4;
49     fRSec.Set(fNsec);    
50     fNDiv.Set(fNsec);      
51     fDpxD.Set(fNsec);      
52     fRSec[0]=fRSec[1]=fRSec[2]=fRSec[3]=0;     
53     fNDiv[0]=fNDiv[1]=fNDiv[2]=fNDiv[3]=0;     
54     fDpxD[0]=fDpxD[1]=fDpxD[2]=fDpxD[3]=0;     
55     fCorr = new TObjArray(3);
56     (*fCorr)[0]=0;
57     (*fCorr)[1]=0;
58     (*fCorr)[2]=0;
59
60
61 Float_t AliMUONSegmentationV01::Dpx(Int_t isec)
62 {
63 //
64 // Returns x-pad size for given sector isec
65    return fDpxD[isec];
66 }
67
68 Float_t AliMUONSegmentationV01::Dpy(Int_t isec)
69 {
70 //
71 // Returns y-pad size for given sector isec
72    return fDpy;
73 }
74     
75 void   AliMUONSegmentationV01::SetSegRadii(Float_t  r[4])
76 {
77 //
78 // Set the radii of the segmentation zones 
79     for (Int_t i=0; i<4; i++) {
80         fRSec[i]=r[i];
81         printf("\n R %d %f \n",i,fRSec[i]);
82         
83     }
84 }
85
86
87 void AliMUONSegmentationV01::SetPadDivision(Int_t ndiv[4])
88 {
89 //
90 // Defines the pad size perp. to the anode wire (y) for different sectors. 
91 // Pad sizes are defined as integral fractions ndiv of a basis pad size
92 // fDpx
93 // 
94     for (Int_t i=0; i<4; i++) {
95         fNDiv[i]=ndiv[i];
96         printf("\n Ndiv %d %d \n",i,fNDiv[i]);
97     }
98     ndiv[0]=ndiv[1];
99 }
100
101
102 void AliMUONSegmentationV01::Init(AliMUONChamber* Chamber)
103 {
104 //
105 //  Fill the arrays fCx (x-contour) and fNpxS (ix-contour) for each sector
106 //  These arrays help in converting from real to pad co-ordinates and
107 //  vice versa.
108 //  This version approximates concentric segmentation zones
109 //
110     Int_t isec;
111     printf("\n Initialise segmentation v01 -- test !!!!!!!!!!!!!! \n");
112     fNpy=Int_t(fRSec[fNsec-1]/fDpy)+1;
113
114     fDpxD[fNsec-1]=fDpx;
115     if (fNsec > 1) {
116         for (Int_t i=fNsec-2; i>=0; i--){
117             fDpxD[i]=fDpxD[fNsec-1]/fNDiv[i];
118             printf("\n test ---dx %d %f \n",i,fDpxD[i]);
119         }
120     }
121 //
122 // fill the arrays defining the pad segmentation boundaries
123     Float_t ry;
124     Int_t   dnx;
125     Int_t   add;
126 //
127 //  loop over sections
128     for(isec=0; isec<fNsec; isec++) {
129 //  
130 //  loop over pads along the aode wires
131         for (Int_t iy=1; iy<=fNpy; iy++) {
132 //
133             Float_t x=iy*fDpy-fDpy/2;
134             if (x > fRSec[isec]) {
135                 fNpxS[isec][iy]=0;
136                 fCx[isec][iy]=0;
137             } else {
138                 ry=TMath::Sqrt(fRSec[isec]*fRSec[isec]-x*x);
139                 if (isec > 1) {
140                     dnx= Int_t((ry-fCx[isec-1][iy])/fDpxD[isec]);
141                     if (isec < fNsec-1) {
142                         if (TMath::Odd((Long_t)dnx)) dnx++;             
143                     }
144                     fNpxS[isec][iy]=fNpxS[isec-1][iy]+dnx;
145                     fCx[isec][iy]=fCx[isec-1][iy]+dnx*fDpxD[isec];
146                 } else if (isec == 1) {
147                     dnx= Int_t((ry-fCx[isec-1][iy])/fDpxD[isec]);
148                     fNpxS[isec][iy]=fNpxS[isec-1][iy]+dnx;
149                     add=4 - (fNpxS[isec][iy])%4;
150                     if (add < 4) fNpxS[isec][iy]+=add; 
151                     dnx=fNpxS[isec][iy]-fNpxS[isec-1][iy];
152                     fCx[isec][iy]=fCx[isec-1][iy]+dnx*fDpxD[isec];
153                 } else {
154                     dnx=Int_t(ry/fDpxD[isec]);
155                     fNpxS[isec][iy]=dnx;
156                     fCx[isec][iy]=dnx*fDpxD[isec];
157                 }
158             }
159         } // y-pad loop
160     } // sector loop
161 }
162
163 Int_t AliMUONSegmentationV01::Sector(Int_t ix, Int_t iy)
164 {
165 // Returns sector number for given pad position
166 //
167     Int_t absix=TMath::Abs(ix);
168     Int_t absiy=TMath::Abs(iy);
169     Int_t isec=0;
170     for (Int_t i=0; i<fNsec; i++) {
171         if (absix<=fNpxS[i][absiy]){
172             isec=i;
173             break;
174         }
175     }
176     return isec;
177 }
178
179 void AliMUONSegmentationV01::
180 GetPadIxy(Float_t x, Float_t y, Int_t &ix, Int_t &iy)
181 {
182 //  Returns pad coordinates (ix,iy) for given real coordinates (x,y)
183 //
184     iy = (y>0)? Int_t(y/fDpy)+1 : Int_t(y/fDpy)-1;
185     if (iy >  fNpy) iy= fNpy;
186     if (iy < -fNpy) iy=-fNpy;
187 //
188 //  Find sector isec
189     Int_t isec=-1;
190     Float_t absx=TMath::Abs(x);
191     Int_t absiy=TMath::Abs(iy);
192     for (Int_t i=0; i < fNsec; i++) {
193         if (absx <= fCx[i][absiy]) {
194             isec=i;
195             break;
196         }
197     }
198     if (isec>0) {
199         ix= Int_t((absx-fCx[isec-1][absiy])/fDpxD[isec])
200             +fNpxS[isec-1][absiy]+1;
201     } else if (isec == 0) {
202         ix= Int_t(absx/fDpxD[isec])+1;
203     } else {
204         ix=fNpxS[fNsec-1][absiy]+1;     
205     }
206     ix = (x>0) ? ix:-ix;
207 }
208
209 void AliMUONSegmentationV01::
210 GetPadCxy(Int_t ix, Int_t iy, Float_t &x, Float_t &y)
211 {
212 //  Returns real coordinates (x,y) for given pad coordinates (ix,iy)
213 //
214     y = (iy>0) ? Float_t(iy*fDpy)-fDpy/2. : Float_t(iy*fDpy)+fDpy/2.;
215 //
216 //  Find sector isec
217     Int_t isec=AliMUONSegmentationV01::Sector(ix,iy);
218 //
219     Int_t absix=TMath::Abs(ix);
220     Int_t absiy=TMath::Abs(iy);
221     if (isec) {
222         x=fCx[isec-1][absiy]+(absix-fNpxS[isec-1][absiy])*fDpxD[isec];
223         x=(ix>0) ?  x-fDpxD[isec]/2 : -x+fDpxD[isec]/2;
224     } else {
225         x=y=0;
226     }
227 }
228
229 void AliMUONSegmentationV01::
230 SetPad(Int_t ix, Int_t iy)
231 {
232     //
233     // Sets virtual pad coordinates, needed for evaluating pad response 
234     // outside the tracking program 
235     GetPadCxy(ix,iy,fx,fy);
236     fSector=Sector(ix,iy);
237 }
238
239
240 void AliMUONSegmentationV01::
241 FirstPad(Float_t xhit, Float_t yhit, Float_t dx, Float_t dy)
242 {
243 // Initialises iteration over pads for charge distribution algorithm
244 //
245     //
246     // Find the wire position (center of charge distribution)
247     Float_t x0a=GetAnod(xhit);
248     fxhit=x0a;
249     fyhit=yhit;
250     
251     //
252     // and take fNsigma*sigma around this center
253     Float_t x01=x0a  - dx;
254     Float_t x02=x0a  + dx;
255     Float_t y01=yhit - dy;
256     Float_t y02=yhit + dy;
257     //
258     // find the pads over which the charge distributes
259     GetPadIxy(x01,y01,fixmin,fiymin);
260     GetPadIxy(x02,y02,fixmax,fiymax);
261     fxmin=x01;
262     fxmax=x02;
263     fymin=y01;
264     fymax=y02;
265     
266     // 
267     // Set current pad to lower left corner
268     if (fixmax < fixmin) fixmax=fixmin;
269     if (fiymax < fiymin) fiymax=fiymin;    
270     fix=fixmin;
271     fiy=fiymin;
272     GetPadCxy(fix,fiy,fx,fy);
273 }
274
275
276 void AliMUONSegmentationV01::NextPad()
277 {
278 // Stepper for the iteration over pads
279 //
280 // Step to next pad in the integration region
281   // 
282   // Step to next pad in integration region
283     Float_t xc,yc;
284     Int_t   iyc;
285     
286 //  step from left to right    
287     if (fx < fxmax && fx != 0) {
288         if (fix==-1) fix++;
289         fix++;
290 //  step up 
291     } else if (fiy != fiymax) {
292         if (fiy==-1) fiy++;
293         fiy++;
294 //      get y-position of next row (yc), xc not used here       
295         GetPadCxy(fix,fiy,xc,yc);
296 //      get x-pad coordiante for first pad in row (fix)
297         GetPadIxy(fxmin,yc,fix,iyc);
298     } else {
299         printf("\n Error: Stepping outside integration region\n ");
300     }
301     GetPadCxy(fix,fiy,fx,fy);
302     fSector=Sector(fix,fiy);
303     if (MorePads() && 
304         (fSector ==-1 || fSector==0)) 
305         NextPad();
306 }
307
308 Int_t AliMUONSegmentationV01::MorePads()
309 // Stopping condition for the iterator over pads
310 //
311 // Are there more pads in the integration region
312 {
313     if ((fx >= fxmax  && fiy >= fiymax) || fy==0) {
314         return 0;
315     } else {
316         return 1;
317     }
318 }
319
320 void AliMUONSegmentationV01::
321 IntegrationLimits(Float_t& x1,Float_t& x2,Float_t& y1, Float_t& y2)
322 {
323 //  Returns integration limits for current pad
324 //
325     x1=fxhit-fx-Dpx(fSector)/2.;
326     x2=x1+Dpx(fSector);
327     y1=fyhit-fy-Dpy(fSector)/2.;
328     y2=y1+Dpy(fSector);    
329 }
330
331 void AliMUONSegmentationV01::
332 Neighbours(Int_t iX, Int_t iY, Int_t* Nlist, Int_t Xlist[10], Int_t Ylist[10])
333 {
334 // Returns list of next neighbours for given Pad (iX, iY)
335 //
336     const Float_t kEpsilon=fDpy/1000;
337     
338     Float_t x,y;
339     Int_t   ixx, iyy, isec1;
340 //
341     Int_t   isec0=AliMUONSegmentationV01::Sector(iX,iY);
342     Int_t i=0;
343 //    
344 //  step right
345     Xlist[i]=iX+1;
346     if (Xlist[i]==0) Xlist[i]++;
347     Ylist[i++]=iY;
348 //
349 //  step left    
350     Xlist[i]=iX-1;
351     if (Xlist[i]==0) Xlist[i]--;
352     Ylist[i++]=iY;
353 //
354 //  step up 
355     AliMUONSegmentationV01::GetPadCxy(iX,iY,x,y);
356     AliMUONSegmentationV01::GetPadIxy(x+kEpsilon,y+fDpy,ixx,iyy);
357     Xlist[i]=ixx;
358     Ylist[i++]=iyy;
359     isec1=AliMUONSegmentationV01::Sector(ixx,iyy);
360     if (isec1==isec0) {
361 //
362 //  no sector boundary crossing
363 //      Xlist[i]=ixx+1;
364 //      Ylist[i++]=iY+1;
365         
366 //      Xlist[i]=ixx-1;
367 //      Ylist[i++]=iY+1;
368     } else if (isec1 < isec0) {
369 // finer segmentation
370 //      Xlist[i]=ixx+1;
371 //      Ylist[i++]=iY+1;
372         
373         Xlist[i]=ixx-1;
374         Ylist[i++]=iyy;
375         
376 //      Xlist[i]=ixx-2;
377 //      Ylist[i++]=iY+1;
378     } else {
379 // coarser segmenation  
380 /*
381         if (TMath::Odd(iX-fNpxS[isec1-1][iY+1])) {
382             Xlist[i]=ixx-1;
383             Ylist[i++]=iY+1;
384         } else {
385             Xlist[i]=ixx+1;
386             Ylist[i++]=iY+1;
387         }
388 */
389     }
390
391 //
392 // step down 
393     AliMUONSegmentationV01::GetPadCxy(iX,iY,x,y);
394     AliMUONSegmentationV01::GetPadIxy(x+kEpsilon,y-fDpy,ixx,iyy);
395     Xlist[i]=ixx;
396     Ylist[i++]=iyy;
397     isec1=AliMUONSegmentationV01::Sector(ixx,iyy);
398     if (isec1==isec0) {
399 //
400 //  no sector boundary crossing
401 /*
402     Xlist[i]=ixx+1;
403     Ylist[i++]=iY-1;
404         
405     Xlist[i]=ixx-1;
406     Ylist[i++]=iY-1;
407 */
408     } else if (isec1 < isec0) {
409 // finer segmentation
410 //      Xlist[i]=ixx+1;
411 //      Ylist[i++]=iY-1;
412         
413         Xlist[i]=ixx-1;
414         Ylist[i++]=iyy;
415         
416 //      Xlist[i]=ixx-2;
417 //      Ylist[i++]=iY-1;
418     } else {
419 // coarser segmentation 
420 /*
421         if (TMath::Odd(iX-fNpxS[isec1-1][iY-1])) {
422             Xlist[i]=ixx-1;
423             Ylist[i++]=iY-1;
424         } else {
425             Xlist[i]=ixx+1;
426             Ylist[i++]=iY-1;
427         }
428 */
429     }
430     *Nlist=i;
431 }
432
433 void AliMUONSegmentationV01::GiveTestPoints(Int_t &n, Float_t *x, Float_t *y)
434 {
435 // Returns test point on the pad plane.
436 // Used during determination of the segmoid correction of the COG-method
437
438     n=3;
439     x[0]=(fRSec[0]+fRSec[1])/2/TMath::Sqrt(2.);
440     y[0]=x[0];
441     x[1]=(fRSec[1]+fRSec[2])/2/TMath::Sqrt(2.);
442     y[1]=x[1];
443     x[2]=(fRSec[2]+fRSec[3])/2/TMath::Sqrt(2.);
444     y[2]=x[2];
445 }
446
447 void AliMUONSegmentationV01::Draw()
448 {
449 // Draws the segmentation zones
450 //
451     TBox *box;
452     
453     Float_t dx=0.95/fCx[3][1]/2;
454     Float_t dy=0.95/(Float_t(Npy()))/2;
455     Float_t x0,y0,x1,y1;
456     Float_t xc=0.5;
457     Float_t yc=0.5;
458     
459     for (Int_t iy=1; iy<Npy(); iy++)
460     {
461         for (Int_t isec=0; isec<4; isec++) {
462             if (isec==0) {
463                 x0=0;
464                 x1=fCx[isec][iy]*dx;
465             } else {
466                 x0=fCx[isec-1][iy]*dx;
467                 x1=fCx[isec][iy]*dx;
468             }
469             y0=Float_t(iy-1)*dy;
470             y1=y0+dy;
471             box=new TBox(x0+xc,y0+yc,x1+xc,y1+yc);
472             box->SetFillColor(isec+1);
473             box->Draw();
474
475             box=new TBox(-x1+xc,y0+yc,-x0+xc,y1+yc);
476             box->SetFillColor(isec+1);
477             box->Draw();
478
479             box=new TBox(x0+xc,-y1+yc,x1+xc,-y0+yc);
480             box->SetFillColor(isec+1);
481             box->Draw();
482
483             box=new TBox(-x1+xc,-y1+yc,-x0+xc,-y0+yc);
484             box->SetFillColor(isec+1);
485             box->Draw();
486         }
487     }
488 }
489 void AliMUONSegmentationV01::SetCorrFunc(Int_t isec, TF1* func)
490 {
491     (*fCorr)[isec]=func;
492 }
493
494 TF1* AliMUONSegmentationV01::CorrFunc(Int_t isec)
495
496     return (TF1*) (*fCorr)[isec];
497 }
498
499 AliMUONSegmentationV01& AliMUONSegmentationV01::operator 
500 =(const AliMUONSegmentationV01 & rhs)
501 {
502 // Dummy assignment operator
503     return *this;
504 }