5a60360ea491a38d51f322bcb753ad49b471b791
[u/mrichter/AliRoot.git] / MUON / AliMUONSegResV01.cxx
1 /////////////////////////////////////////////////////
2 //  Segmentation and Response classes version 01   //
3 /////////////////////////////////////////////////////
4
5 #include <TTUBE.h>
6 #include <TNode.h> 
7 #include <TBox.h> 
8 #include <TRandom.h> 
9
10 #include "AliMUONSegResV01.h"
11 #include "AliRun.h"
12 #include "AliMC.h"
13 #include "iostream.h"
14
15 //___________________________________________
16 ClassImp(AliMUONsegmentationV01)
17
18 Float_t AliMUONsegmentationV01::Dpx(Int_t isec)
19 {
20    return fDpxD[isec];
21 }
22
23 Float_t AliMUONsegmentationV01::Dpy(Int_t isec)
24 {
25    return fDpy;
26 }
27
28 AliMUONsegmentationV01::AliMUONsegmentationV01() 
29 {
30     fNsec=4;
31     fRSec.Set(fNsec);    
32     fNDiv.Set(fNsec);      
33     fDpxD.Set(fNsec);      
34     fRSec[0]=fRSec[1]=fRSec[2]=fRSec[3]=0;     
35     fNDiv[0]=fNDiv[1]=fNDiv[2]=fNDiv[3]=0;     
36     fDpxD[0]=fDpxD[1]=fDpxD[2]=fDpxD[3]=0;     
37     fCorr = new TObjArray(3);
38     (*fCorr)[0]=0;
39     (*fCorr)[1]=0;
40     (*fCorr)[2]=0;
41
42     
43 void   AliMUONsegmentationV01::SetSegRadii(Float_t  r[4])
44 {
45     for (Int_t i=0; i<4; i++) {
46         fRSec[i]=r[i];
47         printf("\n R %d %f \n",i,fRSec[i]);
48         
49     }
50 }
51
52
53 void AliMUONsegmentationV01::SetPadDivision(Int_t ndiv[4])
54 {
55 //
56 // Defines the pad size perp. to the anode wire (y) for different sectors. 
57 //
58     for (Int_t i=0; i<4; i++) {
59         fNDiv[i]=ndiv[i];
60         printf("\n Ndiv %d %d \n",i,fNDiv[i]);
61     }
62     ndiv[0]=ndiv[1];
63 }
64
65
66 void AliMUONsegmentationV01::Init(AliMUONchamber* Chamber)
67 {
68 //
69 //  Fill the arrays fCx (x-contour) and fNpxS (ix-contour) for each sector
70 //  These arrays help in converting from real to pad co-ordinates and
71 //  vice versa
72 //
73     Int_t isec;
74     printf("\n Initialise segmentation v01 -- test !!!!!!!!!!!!!! \n");
75     fNpy=Int_t(fRSec[fNsec-1]/fDpy)+1;
76
77     fDpxD[fNsec-1]=fDpx;
78     if (fNsec > 1) {
79         for (Int_t i=fNsec-2; i>=0; i--){
80             fDpxD[i]=fDpxD[fNsec-1]/fNDiv[i];
81             printf("\n test ---dx %d %f \n",i,fDpxD[i]);
82         }
83     }
84 //
85 // fill the arrays defining the pad segmentation boundaries
86     Float_t ry;
87     Int_t   dnx;
88     Int_t   add;
89 //
90 //  loop over sections
91     for(isec=0; isec<fNsec; isec++) {
92 //  
93 //  loop over pads along the aode wires
94         for (Int_t iy=1; iy<=fNpy; iy++) {
95 //
96             Float_t x=iy*fDpy-fDpy/2;
97             if (x > fRSec[isec]) {
98                 fNpxS[isec][iy]=0;
99                  fCx[isec][iy]=0;
100             } else {
101                 ry=TMath::Sqrt(fRSec[isec]*fRSec[isec]-x*x);
102                 if (isec > 1) {
103                     dnx= Int_t((ry-fCx[isec-1][iy])/fDpxD[isec]);
104                     if (isec < fNsec-1) {
105                         if (TMath::Odd((Long_t)dnx)) dnx++;             
106                     }
107                     fNpxS[isec][iy]=fNpxS[isec-1][iy]+dnx;
108                     fCx[isec][iy]=fCx[isec-1][iy]+dnx*fDpxD[isec];
109                 } else if (isec == 1) {
110                     dnx= Int_t((ry-fCx[isec-1][iy])/fDpxD[isec]);
111                     fNpxS[isec][iy]=fNpxS[isec-1][iy]+dnx;
112                     add=4 - (fNpxS[isec][iy])%4;
113                     if (add < 4) fNpxS[isec][iy]+=add; 
114                     dnx=fNpxS[isec][iy]-fNpxS[isec-1][iy];
115                     fCx[isec][iy]=fCx[isec-1][iy]+dnx*fDpxD[isec];
116                 } else {
117                     dnx=Int_t(ry/fDpxD[isec]);
118                     fNpxS[isec][iy]=dnx;
119                     fCx[isec][iy]=dnx*fDpxD[isec];
120                 }
121             }
122         } // y-pad loop
123     } // sector loop
124     //   
125     // for debugging only 
126    
127     printf("segmentationv01 - I was here ! \n");
128 }
129
130 Int_t AliMUONsegmentationV01::Sector(Int_t ix, Int_t iy)
131 {
132     Int_t absix=TMath::Abs(ix);
133     Int_t absiy=TMath::Abs(iy);
134     Int_t isec=0;
135     for (Int_t i=0; i<fNsec; i++) {
136         if (absix<=fNpxS[i][absiy]){
137             isec=i;
138             break;
139         }
140     }
141     return isec;
142 }
143
144     void AliMUONsegmentationV01::
145     GetPadIxy(Float_t x, Float_t y, Int_t &ix, Int_t &iy)
146 {
147 //  returns pad coordinates (ix,iy) for given real coordinates (x,y)
148 //
149     iy = (y>0)? Int_t(y/fDpy)+1 : Int_t(y/fDpy)-1;
150     if (iy >  fNpy) iy= fNpy;
151     if (iy < -fNpy) iy=-fNpy;
152 //
153 //  Find sector isec
154     Int_t isec=-1;
155     Float_t absx=TMath::Abs(x);
156     Int_t absiy=TMath::Abs(iy);
157     for (Int_t i=0; i < fNsec; i++) {
158         if (absx <= fCx[i][absiy]) {
159             isec=i;
160             break;
161         }
162     }
163     if (isec>0) {
164         ix= Int_t((absx-fCx[isec-1][absiy])/fDpxD[isec])
165             +fNpxS[isec-1][absiy]+1;
166     } else if (isec == 0) {
167         ix= Int_t(absx/fDpxD[isec])+1;
168     } else {
169         ix=fNpxS[fNsec-1][absiy]+1;     
170     }
171 //    printf("\n something %d %d \n",isec,absiy);
172     
173     ix = (x>0) ? ix:-ix;
174 }
175
176 void AliMUONsegmentationV01::
177 GetPadCxy(Int_t ix, Int_t iy, Float_t &x, Float_t &y)
178 {
179 //  returns real coordinates (x,y) for given pad coordinates (ix,iy)
180 //
181     y = (iy>0) ? Float_t(iy*fDpy)-fDpy/2. : Float_t(iy*fDpy)+fDpy/2.;
182 //
183 //  Find sector isec
184     Int_t isec=AliMUONsegmentationV01::Sector(ix,iy);
185 //
186     Int_t absix=TMath::Abs(ix);
187     Int_t absiy=TMath::Abs(iy);
188     if (isec) {
189         x=fCx[isec-1][absiy]+(absix-fNpxS[isec-1][absiy])*fDpxD[isec];
190         x=(ix>0) ?  x-fDpxD[isec]/2 : -x+fDpxD[isec]/2;
191     } else {
192         x=y=0;
193     }
194 }
195
196 void AliMUONsegmentationV01::
197 SetPad(Int_t ix, Int_t iy)
198 {
199     GetPadCxy(ix,iy,fx,fy);
200     fSector=Sector(ix,iy);
201 }
202
203
204 void AliMUONsegmentationV01::
205 FirstPad(Float_t xhit, Float_t yhit, Float_t dx, Float_t dy)
206 {
207     //
208     // Find the wire position (center of charge distribution)
209     Float_t x0a=GetAnod(xhit);
210     fxhit=x0a;
211     fyhit=yhit;
212     
213     //
214     // and take fNsigma*sigma around this center
215     Float_t x01=x0a  - dx;
216     Float_t x02=x0a  + dx;
217     Float_t y01=yhit - dy;
218     Float_t y02=yhit + dy;
219     //
220     // find the pads over which the charge distributes
221     GetPadIxy(x01,y01,fixmin,fiymin);
222     GetPadIxy(x02,y02,fixmax,fiymax);
223     fxmin=x01;
224     fxmax=x02;
225     fymin=y01;
226     fymax=y02;
227     
228     // 
229     // Set current pad to lower left corner
230     if (fixmax < fixmin) fixmax=fixmin;
231     if (fiymax < fiymin) fiymax=fiymin;    
232     fix=fixmin;
233     fiy=fiymin;
234     GetPadCxy(fix,fiy,fx,fy);
235 }
236
237
238 void AliMUONsegmentationV01::NextPad()
239 {
240   // 
241   // Step to next pad in integration region
242     Float_t xc,yc;
243     Int_t   iyc;
244     
245 //  step from left to right    
246     if (fx < fxmax && fx != 0) {
247         if (fix==-1) fix++;
248         fix++;
249 //  step up 
250     } else if (fiy != fiymax) {
251         if (fiy==-1) fiy++;
252         fiy++;
253 //      get y-position of next row (yc), xc not used here       
254         GetPadCxy(fix,fiy,xc,yc);
255 //      get x-pad coordiante for 1 pad in row (fix)
256         GetPadIxy(fxmin,yc,fix,iyc);
257     } else {
258         printf("\n Error: Stepping outside integration region\n ");
259     }
260     GetPadCxy(fix,fiy,fx,fy);
261     fSector=Sector(fix,fiy);
262     if (MorePads() && 
263         (fSector ==-1 || fSector==0 || 
264          TMath::Abs(fx)<1.5 || TMath::Abs(fy)<1.5)) 
265         NextPad();
266     
267 //    printf("\n this pad %f %f %d %d \n",fx,fy,fix,fiy);
268     
269 }
270
271 Int_t AliMUONsegmentationV01::MorePads()
272 //
273 // Are there more pads in the integration region
274 {
275     if ((fx >= fxmax  && fiy >= fiymax) || fy==0) {
276         return 0;
277     } else {
278         return 1;
279     }
280 }
281
282 void AliMUONsegmentationV01::
283 IntegrationLimits(Float_t& x1,Float_t& x2,Float_t& y1, Float_t& y2)
284 {
285     x1=fxhit-fx-Dpx(fSector)/2.;
286     x2=x1+Dpx(fSector);
287     y1=fyhit-fy-Dpy(fSector)/2.;
288     y2=y1+Dpy(fSector);    
289 }
290
291 void AliMUONsegmentationV01::
292 Neighbours(Int_t iX, Int_t iY, Int_t* Nlist, Int_t Xlist[10], Int_t Ylist[10])
293 {
294     const Float_t epsilon=fDpy/1000;
295     
296     Float_t x,y;
297     Int_t   ixx, iyy, isec1;
298 //
299     Int_t   isec0=AliMUONsegmentationV01::Sector(iX,iY);
300     Int_t i=0;
301 //    
302 //  step right
303     Xlist[i]=iX+1;
304     Ylist[i++]=iY;
305 //
306 //  step left    
307     Xlist[i]=iX-1;
308     Ylist[i++]=iY;
309 //
310 //  step up 
311     AliMUONsegmentationV01::GetPadCxy(iX,iY,x,y);
312     AliMUONsegmentationV01::GetPadIxy(x+epsilon,y+fDpy,ixx,iyy);
313     Xlist[i]=ixx;
314     Ylist[i++]=iY+1;
315     isec1=AliMUONsegmentationV01::Sector(ixx,iyy);
316     if (isec1==isec0) {
317 //
318 //  no sector boundary crossing
319         Xlist[i]=ixx+1;
320         Ylist[i++]=iY+1;
321         
322         Xlist[i]=ixx-1;
323         Ylist[i++]=iY+1;
324     } else if (isec1 < isec0) {
325 // finer segmentation
326         Xlist[i]=ixx+1;
327         Ylist[i++]=iY+1;
328         
329         Xlist[i]=ixx-1;
330         Ylist[i++]=iY+1;
331         
332         Xlist[i]=ixx-2;
333         Ylist[i++]=iY+1;
334     } else {
335 // coarser segmenation  
336
337         if (TMath::Odd(iX-fNpxS[isec1-1][iY+1])) {
338             Xlist[i]=ixx-1;
339             Ylist[i++]=iY+1;
340         } else {
341             Xlist[i]=ixx+1;
342             Ylist[i++]=iY+1;
343         }
344     }
345 //
346 // step down 
347     AliMUONsegmentationV01::GetPadCxy(iX,iY,x,y);
348     AliMUONsegmentationV01::GetPadIxy(x+epsilon,y-fDpy,ixx,iyy);
349     Xlist[i]=ixx;
350     Ylist[i++]=iY-1;
351     isec1=AliMUONsegmentationV01::Sector(ixx,iyy);
352     if (isec1==isec0) {
353 //
354 //  no sector boundary crossing
355     Xlist[i]=ixx+1;
356     Ylist[i++]=iY-1;
357         
358     Xlist[i]=ixx-1;
359     Ylist[i++]=iY-1;
360     } else if (isec1 < isec0) {
361 // finer segmentation
362         Xlist[i]=ixx+1;
363         Ylist[i++]=iY-1;
364         
365         Xlist[i]=ixx-1;
366         Ylist[i++]=iY-1;
367         
368         Xlist[i]=ixx-2;
369         Ylist[i++]=iY-1;
370     } else {
371 // coarser segmentation 
372
373         if (TMath::Odd(iX-fNpxS[isec1-1][iY-1])) {
374             Xlist[i]=ixx-1;
375             Ylist[i++]=iY-1;
376         } else {
377             Xlist[i]=ixx+1;
378             Ylist[i++]=iY-1;
379         }
380
381     }
382     *Nlist=i;
383 }
384
385 //___________________________________________
386 void AliMUONsegmentationV01::
387 FitXY(AliMUONRecCluster* Cluster,TClonesArray* MUONdigits)
388     // Default : Centre of gravity method
389 {
390     Float_t x=0;
391     Float_t y=0;
392     Float_t q=0;
393     Float_t xToAdd;
394     Float_t yToAdd;
395     
396     if (gAlice->TreeD()->GetReadEvent() != Cluster->GetCathod()+1)
397         // next line warns if in the future cathod 1 is not event 2 !
398         printf("ClusterFillXY : not reading the right cathod !\n");
399     for(Int_t clusDigit=Cluster->FirstDigitIndex();
400         clusDigit!=Cluster->InvalidDigitIndex();
401         clusDigit=Cluster->NextDigitIndex()) {
402         AliMUONdigit* pDigit=(AliMUONdigit*)MUONdigits->UncheckedAt(clusDigit);
403         GetPadCxy(pDigit->fPadX,pDigit->fPadY,xToAdd,yToAdd);
404         x+= xToAdd*pDigit->fSignal;
405         y+= yToAdd*pDigit->fSignal;
406         q+= (Float_t) pDigit->fSignal;
407     }
408     Cluster->fX=x/q;
409     Cluster->fY=y/q;
410 }
411
412 void AliMUONsegmentationV01::GiveTestPoints(Int_t &n, Float_t *x, Float_t *y)
413 {
414     n=3;
415     x[0]=(fRSec[0]+fRSec[1])/2/TMath::Sqrt(2.);
416     y[0]=x[0];
417     x[1]=(fRSec[1]+fRSec[2])/2/TMath::Sqrt(2.);
418     y[1]=x[1];
419     x[2]=(fRSec[2]+fRSec[3])/2/TMath::Sqrt(2.);
420     y[2]=x[2];
421 }
422
423 void AliMUONsegmentationV01::Draw()
424 {
425     TBox *box;
426     
427     Float_t dx=0.95/fCx[3][1]/2;
428     Float_t dy=0.95/(Float_t(Npy()))/2;
429     Float_t x0,y0,x1,y1;
430     Float_t xc=0.5;
431     Float_t yc=0.5;
432     
433     for (Int_t iy=1; iy<Npy(); iy++)
434     {
435         for (Int_t isec=0; isec<4; isec++) {
436             if (isec==0) {
437                 x0=0;
438                 x1=fCx[isec][iy]*dx;
439             } else {
440                 x0=fCx[isec-1][iy]*dx;
441                 x1=fCx[isec][iy]*dx;
442             }
443             y0=Float_t(iy-1)*dy;
444             y1=y0+dy;
445             box=new TBox(x0+xc,y0+yc,x1+xc,y1+yc);
446             box->SetFillColor(isec+1);
447             box->Draw();
448
449             box=new TBox(-x1+xc,y0+yc,-x0+xc,y1+yc);
450             box->SetFillColor(isec+1);
451             box->Draw();
452
453             box=new TBox(x0+xc,-y1+yc,x1+xc,-y0+yc);
454             box->SetFillColor(isec+1);
455             box->Draw();
456
457             box=new TBox(-x1+xc,-y1+yc,-x0+xc,-y0+yc);
458             box->SetFillColor(isec+1);
459             box->Draw();
460         }
461     }
462 }
463 void AliMUONsegmentationV01::SetCorrFunc(Int_t isec, TF1* func)
464 {
465     (*fCorr)[isec]=func;
466     
467 }
468
469 TF1* AliMUONsegmentationV01::CorrFunc(Int_t isec)
470
471     return (TF1*) (*fCorr)[isec];
472 }
473