Modifications needed by the HBT analysis (P.Skowronski)
[u/mrichter/AliRoot.git] / MUON / AliMUONSegmentationV1.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 /* $Id$ */
17
18 /////////////////////////////////////////////////////////
19 //  Manager and hits classes for set:MUON version LYON //
20 /////////////////////////////////////////////////////////
21
22 #include <TMath.h>
23 #include "AliMUONChamber.h"
24 #include "AliMUONSegmentationV1.h"
25 #include "AliRun.h"
26 #include "AliMUON.h"
27 #include "AliMUONChamber.h"
28
29 //___________________________________________
30 ClassImp(AliMUONSegmentationV1)
31
32   AliMUONSegmentationV1::AliMUONSegmentationV1(const AliMUONSegmentationV1& segmentation):AliSegmentation(segmentation)
33 {
34 // Dummy copy constructor
35 }
36
37
38 AliMUONSegmentationV1::AliMUONSegmentationV1()
39
40 {
41   // initizalize the class with default settings
42     fNzone=1;
43     fDAnod=0.0;
44     fDpx=0.0;
45     fDpx=0.0; // forces crash if not initialized by user
46     fNZoneCut[0]=0;
47     fSensOffset=0;
48     fCorr = 0;
49 }
50
51
52 void AliMUONSegmentationV1::Init(Int_t chamber)
53 {
54     // valid only for T5/6
55     // beware : frMin is SENSITIVE radius by definition.
56     AliMUON *pMUON  = (AliMUON *) gAlice->GetModule("MUON");
57     AliMUONChamber* iChamber=&(pMUON->Chamber(chamber));
58
59     frSensMin2 = (iChamber->RInner())*(iChamber->RInner());
60     frSensMax2 = (iChamber->ROuter())*(iChamber->ROuter());
61     fNpx=(Int_t) (iChamber->ROuter()/fDpx) + 1;
62     fNpy=(Int_t) (iChamber->ROuter()/fDpy) + 1;
63     //    fNwire=3;
64     DefaultCut();
65     fCorr=0;
66
67     fZ = iChamber->Z();
68
69 }
70
71 void AliMUONSegmentationV1::DefaultCut(void)
72 {
73 // Set the default cuts
74     SetNzone(3);
75     AddCut(0,5*6,18*8);
76     AddCut(0,9*6,15*8);
77     AddCut(0,11*6,12*8);
78     AddCut(0,12*6,9*8);
79     AddCut(0,13*6,6*8);
80     AddCut(1,6*6,20*12);
81     AddCut(1,12*6,18*12);
82     AddCut(1,15*6,15*12);
83     AddCut(1,18*6,12*12);
84     AddCut(1,21*6,9*12);
85     SetSensOffset(3.0);
86     SetDAnod(0.325);
87 }
88
89 Int_t AliMUONSegmentationV1::GetiAnod(Float_t xhit)
90 {
91 // Get anode number
92     Int_t kwire=Int_t((TMath::Abs(xhit)-fSensOffset)/fDAnod)+1;
93     return (xhit>0) ? kwire : -kwire ;
94 }
95
96 Float_t AliMUONSegmentationV1::GetAnod(Float_t xhit) const
97 {
98 // Get anode position
99     Int_t kwire=Int_t((TMath::Abs(xhit)-fSensOffset)/fDAnod)+1; // to be compatible ...
100     return (xhit>0) ? fDAnod*(kwire-0.5)+fSensOffset : -fDAnod*(kwire-0.5)-fSensOffset ;
101 }
102
103
104 void AliMUONSegmentationV1::SetPadSize(Float_t p1, Float_t p2)
105 {
106 // For chamber T5/6 p1 and p2 should be same for each zone
107     fDpx=p1;
108     fDpy=p2;
109 }
110
111 void AliMUONSegmentationV1::
112 GetPadI(Float_t x, Float_t y, Int_t &ix, Int_t &iy)
113 {
114 //  returns pad coordinates (ix,iy) for given real coordinates (x,y)
115 //
116     ix = (x>0)? Int_t((x-fSensOffset)/fDpx)+1 : Int_t((x+fSensOffset)/fDpx)-1;
117     iy = (y>0)? Int_t((y-fSensOffset)/fDpy)+1 : Int_t((y+fSensOffset)/fDpy)-1;
118 }
119
120 void AliMUONSegmentationV1::
121 GetPadC(Int_t ix, Int_t iy, Float_t &x, Float_t &y)
122 {
123 //  returns real coordinates (x,y) for given pad coordinates (ix,iy)
124 //
125     x = (ix>0) ? (Float_t(ix)-0.5)*fDpx+fSensOffset : (Float_t(ix)+0.5)*fDpx-fSensOffset;
126     y = (iy>0) ? (Float_t(iy)-0.5)*fDpy+fSensOffset : (Float_t(iy)+0.5)*fDpy-fSensOffset;
127 }
128
129 void AliMUONSegmentationV1::AddCut(Int_t Zone, Int_t nX, Int_t nY)
130 {
131 // the pad nX,nY is last INSIDE zone Zone. First pad is labelled 1 and not 0
132     if (Zone+1>=fNzone) 
133 // no cut for last Zone : it is the natural boundary of the chamber
134         printf("AliMUONSegmentationV1::AddCut ==> Zone %d not allowed !\n",Zone);
135     fZoneX[Zone][fNZoneCut[Zone]] = nX;
136     fZoneY[Zone][fNZoneCut[Zone]] = nY;
137     fNZoneCut[Zone]++;
138 }
139
140 Int_t AliMUONSegmentationV1::GetZone(Float_t X, Float_t Y)
141 {
142 // Get segmentation zone
143     Int_t iX, iY;
144     GetPadI(X,Y,iX,iY);
145     return GetZone( iX , iY );
146 }
147
148 Int_t AliMUONSegmentationV1::GetZone(Int_t nX, Int_t nY)
149 {
150 // Beware : first pad begins at 1 !!
151     Int_t aX =  TMath::Abs(nX);
152     Int_t aY =  TMath::Abs(nY);
153     Int_t zone=fNzone-1;
154     for (Int_t iZone=fNzone-2;iZone>=0;iZone--) 
155     {
156         for (Int_t iCut=0;iCut<fNZoneCut[iZone];iCut++)
157             if ( aY<=fZoneY[iZone][iCut] && aX<=fZoneX[iZone][iCut] )
158             {
159                 zone=iZone;
160                 break;
161             } 
162     }
163     return zone;
164 }
165
166 void AliMUONSegmentationV1::
167 SetHit(Float_t xhit, Float_t yhit)
168 {
169 // Find the wire position (center of charge distribution)
170     fXhit=xhit;
171     fYhit=yhit;
172 }
173
174 void AliMUONSegmentationV1::
175 SetPad(Int_t ix, Int_t iy)
176 {
177 // Set current pad position
178     GetPadC(ix,iy,fX,fY);
179 }
180
181
182 void AliMUONSegmentationV1::SetPadCoord(Int_t iX, Int_t iY)
183 {    
184 // Set current pad coordinates
185 GetPadC(iX,iY,fX,fY);
186  Float_t radius2;
187  if ( ( (radius2=fX*fX+fY*fY) > frSensMax2 || radius2 < frSensMin2 ) 
188       && MorePads() )
189      NextPad();
190 }
191
192 void AliMUONSegmentationV1::FirstPad(Float_t xhit, Float_t yhit, Float_t dx, Float_t dy)
193 {
194     //
195     // Find the wire position (center of charge distribution)
196     Float_t x0a=GetAnod(xhit);
197     fXhit=x0a;
198     fYhit=yhit;
199     //
200     // and take fNsigma*sigma around this center
201     Float_t x01=x0a  - dx;
202     Float_t x02=x0a  + dx;
203     Float_t y01=yhit - dy;
204     Float_t y02=yhit + dy;
205
206     // Do not cross over frames...
207     if (x01 * x0a < 0) 
208       x01 = TMath::Sign(fSensOffset, x0a);
209     if (x02 * x0a < 0) 
210       x02 = TMath::Sign(fSensOffset, x0a);
211     if (y01 * yhit < 0) 
212       y01 = TMath::Sign(fSensOffset, yhit);
213     if (y02 * yhit < 0) 
214       y02 = TMath::Sign(fSensOffset, yhit);
215     //
216     // find the pads over which the charge distributes
217     GetPadI(x01,y01,fIxmin,fIymin);
218     GetPadI(x02,y02,fIxmax,fIymax);
219     // 
220     // Set current pad to lower left corner
221     fIx=fIxmin;
222     fIy=fIymin;
223     SetPadCoord(fIx,fIy);
224 }
225
226 void AliMUONSegmentationV1::NextPad()
227 {
228   // 
229   // Step to next pad in integration region
230     if (fIx != fIxmax) {
231         fIx++;
232     } else if (fIy != fIymax) {
233         fIx=fIxmin;
234         fIy++;
235     } else 
236         printf("\n Error: Stepping outside integration region\n ");
237     SetPadCoord(fIx,fIy);
238 }
239
240 Int_t AliMUONSegmentationV1::MorePads()
241 {
242 //
243 // Are there more pads in the integration region
244
245     if (fIx == fIxmax && fIy == fIymax) {
246         return 0;
247     } else {
248         return 1;       
249     }
250 }
251
252 Int_t AliMUONSegmentationV1::IsParallel2(Int_t iX, Int_t iY)
253 {
254 // test if the pad is read in parallel for zone 2
255 // iX and iY are assumed to be positive and starting at 0 numbering (cF. iX)
256 // returns 1 or 2 if read in parallel, 
257 // according to the actual number in the chain, 0 else
258 //
259 // chainage is        result is 
260 // 1  2  3  1  2  3   1 1 1 2 2 2     y
261 // 7  8  9 10 11 12   0 0 0 0 0 0     ^
262 // 4  5  6  4  5  6   1 1 1 2 2 2     +->x
263 //
264
265     if (iY%3==1) return 0;
266     return (iX%6)/3+1;
267 }
268
269 Int_t AliMUONSegmentationV1::IsParallel3(Int_t iX, Int_t iY)
270 {
271 // test if the pad is read in parallel for zone 3
272 // iX and iY are assumed to be positive and starting at 0 numbering (cF. iX)
273 // returns 1,2 or 3 if read in parallel, 
274 // according to the actual number in the chain, 0 else
275 //
276 // chainage is                 result is
277 //16  2  3  1  2  3  1  2  3   0 1 1 1 2 2 2 3 3
278 // 7  8  9 10 11 12 13 14 15   0 0 0 0 0 0 0 0 0
279 // 4  5  6  4  5  6  4  5  6   1 1 1 2 2 2 3 3 3
280 //
281
282     if (iY%3==1) return 0;
283     return (iX%9)/3+1 - (iY%3==2 && iX%3==0);
284 }
285
286 Int_t AliMUONSegmentationV1::NParallel2(Int_t /*iX*/, Int_t iY)
287 {
288 // returns the number of pads connected in parallel for zone 2
289 // iX and iY are assumed to be positive and starting at 0 numbering (cF. iX)
290 //
291 // result is
292 // 2 2 2 2 2 2
293 // 1 1 1 1 1 1
294 // 2 2 2 2 2 2
295 //
296
297     if (iY%3==1) return 1;
298     return 2;
299 }
300
301 Int_t AliMUONSegmentationV1::NParallel3(Int_t iX, Int_t iY)
302 {
303 // test if the pad is read in parallel for zone 3
304 // iX and iY are assumed to be positive and starting at 0 numbering (cF. iX)
305 // returns 1,2 or 3 if read in parallel, 
306 // according to the actual number in the chain, 0 else
307 //
308 // result is 
309 // 1 3 3 2 3 3 2 3 3 
310 // 1 1 1 1 1 1 1 1 1
311 // 3 3 3 3 3 3 3 3 3
312 //
313
314     if (iY%3==1) return 1;
315     if (iY%3==2 && iX%9==0) return 1;
316     return 3 - (iY%3==2 && iX%3==0);
317 }
318
319
320 Int_t AliMUONSegmentationV1::Ix(Int_t trueX, Int_t trueY)
321 {
322 // returns the X number of pad which corresponds to the logical 
323 // channel, expressed in x and y.
324
325     Int_t wix = TMath::Abs(trueX)-1;
326     Int_t wiy = TMath::Abs(trueY)-1;
327     Int_t zone = GetZone(trueX,trueY);
328     Int_t par3;
329     switch (zone) {
330     case 0: return trueX;
331     case 1:
332         if (IsParallel2(wix,wiy) == 2)
333             return (trueX>0)? trueX-3 : trueX+3 ;
334         return trueX;
335     case 2:
336         if ( (par3= IsParallel3(wix,wiy)) )
337             return (trueX>0) ? trueX-3*(par3-1) : trueX+3*(par3-1) ;
338         return trueX ;
339     default :
340         printf("Couille dans AliMUONSegmentationV1::ix\n");
341     }
342     return -1;
343 }
344
345 Int_t AliMUONSegmentationV1::Ix() 
346 {
347 // returns the X number of pad which has to increment charge
348 // due to parallel read-out
349     return Ix(fIx,fIy);
350 }
351
352 Int_t AliMUONSegmentationV1::ISector() 
353 {
354 // This function is of no use for this kind of segmentation.
355     return GetZone(fIx,fIy);
356 }
357
358 void AliMUONSegmentationV1::SigGenInit(Float_t x,Float_t y,Float_t /*z*/)
359 {
360 //
361 //  Initialises pad and wire position during stepping
362     fXt =x;
363     fYt =y;
364     GetPadI(x,y,fIxt,fIyt);
365     fIwt= GetiAnod(x);
366
367 }
368
369 Int_t AliMUONSegmentationV1::SigGenCond(Float_t x,Float_t y,Float_t /*z*/)
370 {
371 //
372 //  Signal will be generated if particle crosses pad boundary or
373 //  boundary between two wires. 
374     Int_t ixt;
375     Int_t iyt;
376     GetPadI(x,y,ixt,iyt);
377     Int_t iwt= GetiAnod(x);
378     
379     if ((ixt != fIxt) || (iyt !=fIyt) || (iwt != fIwt)) {
380         return 1;
381     } else {
382         return 0;
383     }
384 }
385
386 void AliMUONSegmentationV1::
387 IntegrationLimits(Float_t& x1,Float_t& x2,Float_t& y1, Float_t& y2)
388 {
389 // Get integration limits
390     x1=fXhit-fX-fDpx/2.;
391     x2=x1+fDpx;
392     y1=fYhit-fY-fDpy/2.;
393     y2=y1+fDpy;    
394 }
395
396 void AliMUONSegmentationV1::GetNParallelAndOffset(Int_t iX, Int_t iY,Int_t
397 *Nparallel, Int_t* Offset)
398 {
399 // Get parallel pad
400     Int_t wix = TMath::Abs(iX)-1;
401     Int_t wiy = TMath::Abs(iY)-1;
402     Int_t zone = GetZone(iX,iY);
403     switch (zone) {
404     case 0: 
405         *Nparallel=1;
406         *Offset=0;
407         break;
408     case 1:
409         *Nparallel = NParallel2(wix,wiy);
410         (iX>0) ? *Offset =3 : *Offset = -3;
411         if (IsParallel2(wix,wiy)>1)
412             printf("GetNParallelAndOffset called for existing channel -> answer is crazy\n");
413         break;
414     case 2:
415         *Nparallel = NParallel3(wix,wiy);
416         (iX>0) ? *Offset =3 : *Offset = -3;
417         if (IsParallel3(wix,wiy)>1)
418             printf("GetNParallelAndOffset called for existing channel -> answer is crazy\n");
419         break;
420     }
421 }
422
423
424 Float_t AliMUONSegmentationV1::Distance2AndOffset(Int_t iX, Int_t iY, Float_t X, Float_t Y, Int_t *Offset)
425 {
426 //
427 // Computes the offset for which the physical pad has the minimum distance squared
428 // (returned value) to the given coordinates
429
430     Int_t nPara,offset;
431     GetNParallelAndOffset(iX,iY,&nPara,&offset);
432     Float_t d2min=1E10;
433     for (Int_t i=0;i<nPara; i++)
434     {
435         Float_t x,y;
436         GetPadC(iX+i*offset,iY,x,y);
437         Float_t d2=(x-X)*(x-X) + (y-Y)*(y-Y);
438         if ( d2min > d2)
439         {
440             d2min = d2;
441             *Offset = i*offset;
442         }
443     }
444     return d2min; 
445 }
446
447 void AliMUONSegmentationV1::CleanNeighbours(Int_t* Nlist, Int_t *Xlist, 
448         Int_t *Ylist)
449 {
450 // In the raw neighbours list, some pads do not exist 
451 // and some others are read in parallel ...
452 // So we prune non-existing neighbours from the list (event if this should be
453 // at last not be a problem due to the clustering algorithm...)
454
455     Int_t nTot=0;
456     for (Int_t nList=0;nList<*Nlist;nList++)
457     {
458         // prune if it does not exist
459         if ( Xlist[nList]==0 || Ylist[nList]==0 )
460             continue;
461         // compute true position
462         Xlist[nTot] = Ix(Xlist[nList],Ylist[nList]) ;
463         Ylist[nTot] = Ylist[nList] ;
464         // and prune if it does already exist
465         Int_t nTest;
466         for (nTest=0;nTest<nTot; nTest++)
467         {
468             if ( Xlist[nTest]==Xlist[nTot] && Ylist[nTest]==Ylist[nTot])
469                 // we found it
470                 break ;
471         }
472         if (nTest==nTot)
473             nTot++;
474     }
475     *Nlist = nTot;
476 }
477
478 void AliMUONSegmentationV1::
479 NeighboursNonDiag(Int_t iX, Int_t iY, Int_t* Nlist, Int_t Xlist[12], Int_t Ylist[12])
480 {
481 // returns the X number of pad which has to increment charge
482 // due to parallel read-out
483
484     Int_t nParallel, offset;
485     GetNParallelAndOffset(iX,iY,&nParallel,&offset);
486 //
487 // now fill raw list of neighbours
488     *Nlist=4*nParallel;
489     Xlist[0]=Xlist[1]=iX;Xlist[2]=iX-1;Xlist[3]=iX+1;
490     Ylist[0]=iY-1;Ylist[1]=iY+1;Ylist[2]=Ylist[3]=iY;
491     if (nParallel>1) {
492         Xlist[4]=Xlist[5]=iX+offset;Xlist[6]=iX+offset-1;Xlist[7]=iX+offset+1;
493         Ylist[4]=iY-1;Ylist[5]=iY+1;Ylist[6]=Ylist[7]=iY;
494         if (nParallel>2) {
495             Xlist[8]=Xlist[9]=iX+2*offset;Xlist[10]=iX+2*offset-1;Xlist[11]=iX+2*offset+1;
496             Ylist[8]=iY-1;Ylist[9]=iY+1;Ylist[10]=Ylist[11]=iY;
497         }
498     }
499     CleanNeighbours(Nlist,Xlist,Ylist);
500 }
501
502 void AliMUONSegmentationV1::
503 NeighboursDiag(Int_t iX, Int_t iY, Int_t* Nlist, Int_t Xlist[24], Int_t Ylist[24])
504 {
505 // returns the X number of pad which has to increment charge
506 // due to parallel read-out
507
508     Int_t nParallel, offset;
509     GetNParallelAndOffset(iX,iY,&nParallel,&offset);
510 //
511 // now fill raw list of neighbours
512     *Nlist=0;
513     for (Int_t i=0;i<nParallel;i++)
514         for (Int_t dx=-1;dx<2;dx++)
515             for (Int_t dy=-1;dy<2;dy++)
516             {
517                 if (dx==dy && dy==0)
518                     continue; 
519                 Xlist[*Nlist] = iX + dx + i*offset;
520                 Ylist[*Nlist] = iY + dy;
521                 (*Nlist)++;
522             }
523     CleanNeighbours(Nlist,Xlist,Ylist);
524 }
525
526 void AliMUONSegmentationV1::Neighbours(Int_t iX, Int_t iY, Int_t* Nlist, 
527                                        Int_t Xlist[24], Int_t Ylist[24])
528 {
529 // Get neighbours
530 NeighboursDiag(iX,iY,Nlist,Xlist,Ylist);
531 }
532
533
534 void AliMUONSegmentationV1::GiveTestPoints(Int_t &n, Float_t *x, Float_t *y) const
535 {
536 // Return a test point
537     n=1;
538     x[0]=(TMath::Sqrt(frSensMax2)-TMath::Sqrt(frSensMin2))/2/TMath::Sqrt(2.);
539     y[0]=x[0];
540 }
541
542 AliMUONSegmentationV1& AliMUONSegmentationV1::operator =(const AliMUONSegmentationV1 & /*rhs*/)
543 {
544 // Dummy assignment operator
545     return *this;
546 }