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