New version of MUON from M.Bondila & A.Morsch
[u/mrichter/AliRoot.git] / MUON / AliMUONSegResV01.cxx
CommitLineData
a897a37a 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//___________________________________________
16ClassImp(AliMUONsegmentationV01)
17
18Float_t AliMUONsegmentationV01::Dpx(Int_t isec)
19{
20 return fDpxD[isec];
21}
22
23Float_t AliMUONsegmentationV01::Dpy(Int_t isec)
24{
25 return fDpy;
26}
27
28AliMUONsegmentationV01::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
43void 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
53void 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
66void 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
130Int_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
176void AliMUONsegmentationV01::
177GetPadCxy(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
196void AliMUONsegmentationV01::
197SetPad(Int_t ix, Int_t iy)
198{
199 GetPadCxy(ix,iy,fx,fy);
200 fSector=Sector(ix,iy);
201}
202
203
204void AliMUONsegmentationV01::
205FirstPad(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
238void 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
271Int_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
282void AliMUONsegmentationV01::
283IntegrationLimits(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
291void AliMUONsegmentationV01::
292Neighbours(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//___________________________________________
386void AliMUONsegmentationV01::
387FitXY(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
412void 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
423void 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}
463void AliMUONsegmentationV01::SetCorrFunc(Int_t isec, TF1* func)
464{
465 (*fCorr)[isec]=func;
466
467}
468
469TF1* AliMUONsegmentationV01::CorrFunc(Int_t isec)
470{
471 return (TF1*) (*fCorr)[isec];
472}
473