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