Modifications needed by the HBT analysis (P.Skowronski)
[u/mrichter/AliRoot.git] / MUON / AliMUONSegmentationV01.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 //  Segmentation and Response classes version 01   //
20 /////////////////////////////////////////////////////
21
22 #include <TBox.h> 
23 #include <TTUBE.h>
24 #include <TBRIK.h>
25 #include <TNode.h>  
26 #include <TGeometry.h>  
27 #include <TF1.h> 
28 #include <TObjArray.h>
29 #include <Riostream.h>
30
31 #include "AliMUONSegmentationV01.h"
32 #include "AliMUON.h"
33 #include "AliMUONChamber.h"
34 #include "AliRun.h"
35
36
37
38 //___________________________________________
39 ClassImp(AliMUONSegmentationV01)
40
41   AliMUONSegmentationV01::AliMUONSegmentationV01(const AliMUONSegmentationV01& segmentation):AliMUONSegmentationV0(segmentation)
42 {
43 // Dummy copy constructor
44 }
45
46 AliMUONSegmentationV01::AliMUONSegmentationV01() 
47 {
48 // Default constructor
49     fRSec = 0;
50     fNDiv = 0;      
51     fDpxD = 0;
52     fCorrA = 0;
53 }
54
55 AliMUONSegmentationV01::AliMUONSegmentationV01(Int_t nsec) 
56 {
57 //  Non default constructor
58
59     fNsec = nsec;
60     fRSec = new TArrayF(fNsec);
61     fNDiv = new TArrayI(fNsec);      
62     fDpxD = new TArrayF(fNsec);      
63
64
65     (*fRSec)[0]=(*fRSec)[1]=(*fRSec)[2]=(*fRSec)[3]=0;     
66     (*fNDiv)[0]=(*fNDiv)[1]=(*fNDiv)[2]=(*fNDiv)[3]=0;     
67     (*fDpxD)[0]=(*fDpxD)[1]=(*fDpxD)[2]=(*fDpxD)[3]=0;     
68     fCorrA = new TObjArray(3);
69     fCorrA->AddAt(0,0);
70     fCorrA->AddAt(0,1);
71     fCorrA->AddAt(0,2);
72     fOffsetY=0;
73
74
75 AliMUONSegmentationV01::~AliMUONSegmentationV01() 
76 {
77 // Destructor
78     if (fRSec) delete fRSec;
79     if (fNDiv) delete fNDiv;
80     if (fDpxD) delete fDpxD;
81     if (fCorrA) {
82         fCorrA->Delete();
83         delete fCorrA;
84     }
85
86
87
88 Float_t AliMUONSegmentationV01::Dpx(Int_t isec) const
89 {
90 //
91 // Returns x-pad size for given sector isec
92    Float_t dpx = (*fDpxD)[isec];
93    return dpx;
94 }
95
96 Float_t AliMUONSegmentationV01::Dpy(Int_t /*isec*/) const
97 {
98 //
99 // Returns y-pad size for given sector isec
100    return fDpy;
101 }
102     
103 void   AliMUONSegmentationV01::SetSegRadii(Float_t  r[4])
104 {
105 //
106 // Set the radii of the segmentation zones 
107     for (Int_t i=0; i<4; i++) {
108         (*fRSec)[i]=r[i];
109     }
110 }
111
112
113 void AliMUONSegmentationV01::SetPadDivision(Int_t ndiv[4])
114 {
115 //
116 // Defines the pad size perp. to the anode wire (y) for different sectors. 
117 // Pad sizes are defined as integral fractions ndiv of a basis pad size
118 // fDpx
119 // 
120     for (Int_t i=0; i<4; i++) {
121         (*fNDiv)[i]=ndiv[i];
122     }
123     ndiv[0]=ndiv[1];
124 }
125
126
127 void AliMUONSegmentationV01::Init(Int_t chamber)
128 {
129 //
130 //  Fill the arrays fCx (x-contour) and fNpxS (ix-contour) for each sector
131 //  These arrays help in converting from real to pad co-ordinates and
132 //  vice versa.
133 //  This version approximates concentric segmentation zones
134 //
135     Int_t isec;
136     //printf("\n Initialise Segmentation V01\n");
137
138
139     fNpy=Int_t((*fRSec)[fNsec-1]/fDpy)+1;
140
141     (*fDpxD)[fNsec-1]=fDpx;
142     if (fNsec > 1) {
143         for (Int_t i=fNsec-2; i>=0; i--){
144             (*fDpxD)[i]=(*fDpxD)[fNsec-1]/(*fNDiv)[i];
145         }
146     }
147 //
148 // fill the arrays defining the pad segmentation boundaries
149     Float_t ry;
150     Int_t   dnx;
151     Int_t   add;
152 //
153 //  loop over sections
154     for(isec=0; isec<fNsec; isec++) {
155 //  
156 //  loop over pads along the aode wires
157         for (Int_t iy=1; iy<=fNpy; iy++) {
158 //
159             Float_t x=iy*fDpy-fDpy/2;
160             if (x > (*fRSec)[isec]) {
161                 fNpxS[isec][iy]=0;
162                 fCx[isec][iy]=0;
163             } else {
164                 ry=TMath::Sqrt((*fRSec)[isec]*(*fRSec)[isec]-x*x);
165                 if (isec > 1) {
166                     dnx= Int_t((ry-fCx[isec-1][iy])/(*fDpxD)[isec]);
167                     if (isec < fNsec-1) {
168                         if (TMath::Odd((Long_t)dnx)) dnx++;             
169                     }
170                     fNpxS[isec][iy]=fNpxS[isec-1][iy]+dnx;
171                     fCx[isec][iy]=fCx[isec-1][iy]+dnx*(*fDpxD)[isec];
172                 } else if (isec == 1) {
173                     dnx= Int_t((ry-fCx[isec-1][iy])/(*fDpxD)[isec]);
174                     fNpxS[isec][iy]=fNpxS[isec-1][iy]+dnx;
175                     add=4 - (fNpxS[isec][iy])%4;
176                     if (add < 4) fNpxS[isec][iy]+=add; 
177                     dnx=fNpxS[isec][iy]-fNpxS[isec-1][iy];
178                     fCx[isec][iy]=fCx[isec-1][iy]+dnx*(*fDpxD)[isec];
179                 } else {
180                     dnx=Int_t(ry/(*fDpxD)[isec]);
181                     fNpxS[isec][iy]=dnx;
182                     fCx[isec][iy]=dnx*(*fDpxD)[isec];
183                 }
184             }
185         } // y-pad loop
186     } // sector loop
187 // reference to chamber
188     AliMUON *pMUON  = (AliMUON *) gAlice->GetModule("MUON");
189     fChamber=&(pMUON->Chamber(chamber));
190     fZ = fChamber->Z();
191     fId=chamber;
192 }
193
194 Int_t AliMUONSegmentationV01::Sector(Int_t ix, Int_t iy)
195 {
196 // Returns sector number for given pad position
197 //
198     Int_t absix=TMath::Abs(ix);
199     Int_t absiy=TMath::Abs(iy);
200     Int_t isec=0;
201     for (Int_t i=0; i<fNsec; i++) {
202         if (absix<=fNpxS[i][absiy]){
203             isec=i;
204             break;
205         }
206     }
207     return isec;
208 }
209 //______________________________________________________________________
210 void AliMUONSegmentationV01::GetPadI(Float_t x, Float_t y, Int_t &ix, Int_t &iy)
211 {
212 //  Returns pad coordinates (ix,iy) for given real coordinates (x,y)
213 //
214     iy = (y-fOffsetY >0)? 
215       Int_t((y-fOffsetY)/fDpy)+1 
216       : 
217       Int_t((y-fOffsetY)/fDpy)-1;
218   
219     if (iy >  fNpy) iy= fNpy;
220     if (iy < -fNpy) iy=-fNpy;
221 //
222 //  Find sector isec
223     Int_t isec=-1;
224     Float_t absx=TMath::Abs(x);
225     Int_t absiy=TMath::Abs(iy);
226     for (Int_t i=0; i < fNsec; i++) {
227         if (absx <= fCx[i][absiy]) {
228             isec=i;
229             break;
230         }
231     }
232     if (isec>0) {
233         ix= Int_t((absx-fCx[isec-1][absiy])/(*fDpxD)[isec])
234             +fNpxS[isec-1][absiy]+1;
235     } else if (isec == 0) {
236         ix= Int_t(absx/(*fDpxD)[isec])+1;
237     } else {
238         ix=fNpxS[fNsec-1][absiy]+1;     
239     }
240     ix = (x>0) ? ix:-ix;
241 }
242 //________________________________________________________________
243 void AliMUONSegmentationV01::GetPadI(Float_t x, Float_t y , Float_t /*z*/, Int_t &ix, Int_t &iy) 
244 {
245   GetPadI(x, y, ix, iy);
246 }
247 //________________________________________________________________
248
249 void AliMUONSegmentationV01::
250 GetPadC(Int_t ix, Int_t iy, Float_t &x, Float_t &y)
251 {
252 //  Returns real coordinates (x,y) for given pad coordinates (ix,iy)
253 //
254     y = (iy>0) ?
255       Float_t(iy*fDpy)-fDpy/2.+fOffsetY
256       :
257       Float_t(iy*fDpy)+fDpy/2.+fOffsetY; 
258
259 //
260 //  Find sector isec
261     Int_t isec=AliMUONSegmentationV01::Sector(ix,iy);
262 //
263     Int_t absix=TMath::Abs(ix);
264     Int_t absiy=TMath::Abs(iy);
265     if (isec) {
266         x=fCx[isec-1][absiy]+(absix-fNpxS[isec-1][absiy])*(*fDpxD)[isec];
267         x=(ix>0) ?  x-(*fDpxD)[isec]/2 : -x+(*fDpxD)[isec]/2;
268     } else {
269         x=y=0;
270     }
271 }
272
273 void AliMUONSegmentationV01::
274 SetPad(Int_t ix, Int_t iy)
275 {
276     //
277     // Sets virtual pad coordinates, needed for evaluating pad response 
278     // outside the tracking program 
279     GetPadC(ix,iy,fX,fY);
280     fSector=Sector(ix,iy);
281 }
282
283 //______________________________________________________________________
284 void AliMUONSegmentationV01::FirstPad(Float_t xhit, Float_t yhit, Float_t dx, Float_t dy)
285 {
286 // Initialises iteration over pads for charge distribution algorithm
287 //
288     //
289     // Find the wire position (center of charge distribution)
290     Float_t x0a=GetAnod(xhit);
291     fXhit=x0a;
292     fYhit=yhit;
293     
294     //
295     // and take fNsigma*sigma around this center
296     Float_t x01=x0a  - dx;
297     Float_t x02=x0a  + dx;
298     Float_t y01=yhit - dy;
299     Float_t y02=yhit + dy;
300     //
301     // find the pads over which the charge distributes
302
303     GetPadI(x01,y01,fIxmin,fIymin);
304     GetPadI(x02,y02,fIxmax,fIymax);
305     fXmin=x01;
306     fXmax=x02;
307     fYmin=y01;
308     fYmax=y02;
309     
310     // 
311     // Set current pad to lower left corner
312     if (fIxmax < fIxmin) fIxmax=fIxmin;
313     if (fIymax < fIymin) fIymax=fIymin;    
314     fIx=fIxmin;
315     fIy=fIymin;
316     GetPadC(fIx,fIy,fX,fY);
317 }
318
319
320 void AliMUONSegmentationV01::NextPad()
321 {
322 // Stepper for the iteration over pads
323 //
324 // Step to next pad in the integration region
325   // 
326   // Step to next pad in integration region
327     Float_t xc,yc;
328     Int_t   iyc;
329     
330 //  step from left to right    
331
332     if (fX < fXmax && fX != 0) {
333         if (fIx==-1) fIx++;
334         fIx++;
335 //  step up 
336     } else if (fIy != fIymax) {
337         if (fIy==-1) fIy++;
338         fIy++;
339 //      get y-position of next row (yc), xc not used here       
340         GetPadC(fIx,fIy,xc,yc);
341 //      get x-pad coordiante for first pad in row (fIx)
342         GetPadI(fXmin,yc,fIx,iyc);
343     } else {
344         fIx=-1;
345         fIy=-1;
346     }
347     GetPadC(fIx,fIy,fX,fY);
348     fSector=Sector(fIx,fIy);
349     if (MorePads() && 
350         (fSector ==-1 || fSector==0)) 
351         NextPad();
352 }
353
354 Int_t AliMUONSegmentationV01::MorePads()
355
356 {
357 // Stopping condition for the iterator over pads
358 //
359 // Are there more pads in the integration region
360     return  (fIx != -1  || fIy != -1);
361 /*
362     if ((fX >= fXmax  && fIy >= fIymax) || fY==0) {
363         return 0;
364     } else {
365         return 1;
366     }
367 */
368 }
369 //______________________________________________________________________
370 void AliMUONSegmentationV01::FirstPad(Float_t xhit, Float_t yhit, Float_t /*zhit*/, Float_t dx, Float_t dy)
371 {
372   FirstPad(xhit, yhit, dx, dy);
373 }
374
375
376 void AliMUONSegmentationV01::
377 IntegrationLimits(Float_t& x1,Float_t& x2,Float_t& y1, Float_t& y2)
378 {
379 //  Returns integration limits for current pad
380 //
381     x1=fXhit-fX-Dpx(fSector)/2.;
382     x2=x1+Dpx(fSector);
383     y1=fYhit-fY-Dpy(fSector)/2.;
384     y2=y1+Dpy(fSector);    
385 }
386
387 void AliMUONSegmentationV01::
388 Neighbours(Int_t iX, Int_t iY, Int_t* Nlist, Int_t Xlist[10], Int_t Ylist[10])
389 {
390 // Returns list of next neighbours for given Pad (iX, iY)
391 //
392     const Float_t kEpsilon=fDpy/1000;
393     
394     Float_t x,y;
395     Int_t   ixx, iyy, isec1;
396 //
397     Int_t   isec0=AliMUONSegmentationV01::Sector(iX,iY);
398     Int_t i=0;
399 //    
400 //  step right
401     Xlist[i]=iX+1;
402     if (Xlist[i]==0) Xlist[i]++;
403     Ylist[i++]=iY;
404 //
405 //  step left    
406     Xlist[i]=iX-1;
407     if (Xlist[i]==0) Xlist[i]--;
408     Ylist[i++]=iY;
409 //
410 //  step up 
411     AliMUONSegmentationV01::GetPadC(iX,iY,x,y);
412     AliMUONSegmentationV01::GetPadI(x+kEpsilon,y+fDpy,ixx,iyy);
413     Xlist[i]=ixx;
414     Ylist[i++]=iyy;
415     isec1=AliMUONSegmentationV01::Sector(ixx,iyy);
416     if (isec1==isec0) {
417 //
418 //  no sector boundary crossing
419 //      Xlist[i]=ixx+1;
420 //      Ylist[i++]=iY+1;
421         
422 //      Xlist[i]=ixx-1;
423 //      Ylist[i++]=iY+1;
424     } else if (isec1 < isec0) {
425 // finer segmentation
426 //      Xlist[i]=ixx+1;
427 //      Ylist[i++]=iY+1;
428         
429         Xlist[i]=ixx-1;
430         Ylist[i++]=iyy;
431         
432 //      Xlist[i]=ixx-2;
433 //      Ylist[i++]=iY+1;
434     } else {
435 // coarser segmenation  
436 /*
437         if (TMath::Odd(iX-fNpxS[isec1-1][iY+1])) {
438             Xlist[i]=ixx-1;
439             Ylist[i++]=iY+1;
440         } else {
441             Xlist[i]=ixx+1;
442             Ylist[i++]=iY+1;
443         }
444 */
445     }
446
447 //
448 // step down 
449     AliMUONSegmentationV01::GetPadC(iX,iY,x,y);
450     AliMUONSegmentationV01::GetPadI(x+kEpsilon,y-fDpy,ixx,iyy);
451     Xlist[i]=ixx;
452     Ylist[i++]=iyy;
453     isec1=AliMUONSegmentationV01::Sector(ixx,iyy);
454     if (isec1==isec0) {
455 //
456 //  no sector boundary crossing
457 /*
458     Xlist[i]=ixx+1;
459     Ylist[i++]=iY-1;
460         
461     Xlist[i]=ixx-1;
462     Ylist[i++]=iY-1;
463 */
464     } else if (isec1 < isec0) {
465 // finer segmentation
466 //      Xlist[i]=ixx+1;
467 //      Ylist[i++]=iY-1;
468         
469         Xlist[i]=ixx-1;
470         Ylist[i++]=iyy;
471         
472 //      Xlist[i]=ixx-2;
473 //      Ylist[i++]=iY-1;
474     } else {
475 // coarser segmentation 
476 /*
477         if (TMath::Odd(iX-fNpxS[isec1-1][iY-1])) {
478             Xlist[i]=ixx-1;
479             Ylist[i++]=iY-1;
480         } else {
481             Xlist[i]=ixx+1;
482             Ylist[i++]=iY-1;
483         }
484 */
485     }
486     *Nlist=i;
487 }
488
489 void AliMUONSegmentationV01::GiveTestPoints(Int_t &n, Float_t *x, Float_t *y) const
490 {
491 // Returns test point on the pad plane.
492 // Used during determination of the segmoid correction of the COG-method
493
494     n=3;
495     x[0]=((*fRSec)[0]+(*fRSec)[1])/2/TMath::Sqrt(2.);
496     y[0]=x[0];
497     x[1]=((*fRSec)[1]+(*fRSec)[2])/2/TMath::Sqrt(2.);
498     y[1]=x[1];
499     x[2]=((*fRSec)[2]+(*fRSec)[3])/2/TMath::Sqrt(2.);
500     y[2]=x[2];
501 }
502
503 void AliMUONSegmentationV01::Draw(const char* opt) const
504 {
505  
506 // Draws the segmentation zones
507 //
508   if (!strcmp(opt,"eventdisplay")) { 
509     const int kColorMUON  = kBlue;
510
511     TRotMatrix* rot000 = new TRotMatrix("Rot000"," ", 90,  0, 90, 90, 0, 0);
512     TRotMatrix* rot090 = new TRotMatrix("Rot090"," ", 90, 90, 90,180, 0, 0);
513     TRotMatrix* rot180 = new TRotMatrix("Rot180"," ", 90,180, 90,270, 0, 0);
514     TRotMatrix* rot270 = new TRotMatrix("Rot270"," ", 90,270, 90,  0, 0, 0);
515
516     char nameChamber[9], nameSense[9], nameFrame[9], nameNode[9];
517     char nameSense1[9], nameSense2[9];    
518     TNode *node, *nodeF;
519  
520     sprintf(nameChamber,"C_MUON%d",fId+1);
521     sprintf(nameSense,"S_MUON%d",fId+1);
522     sprintf(nameSense1,"S1_MUON%d",fId+1);
523     sprintf(nameSense2,"S2_MUON%d",fId+1);
524     sprintf(nameFrame,"F_MUON%d",fId+1);        
525
526     TNode* top=gAlice->GetGeometry()->GetNode("alice");
527
528     Float_t rmin = (*fRSec)[0]-3;
529     Float_t rmax = (*fRSec)[3]+3;
530     new TTUBE(nameChamber,"Mother","void",rmin,rmax,0.25,1.);
531     rmin = (*fRSec)[0];
532     rmax = (*fRSec)[3];
533     new TTUBE(nameSense,"Sens. region","void",rmin,rmax,0.25, 1.);
534     Float_t dx=(rmax-rmin)/2;
535     Float_t dy=3.;
536     Float_t dz=0.25;
537     TBRIK* frMUON = new TBRIK(nameFrame,"Frame","void",dx,dy,dz);
538     top->cd();
539     sprintf(nameNode,"MUON%d",100+fId+1);
540     node = new TNode(nameNode,"ChamberNode",nameChamber,0,0,fChamber->Z(),"");
541     node->SetLineColor(kColorMUON);
542     AliMUON *pMUON  = (AliMUON *) gAlice->GetModule("MUON");
543     (pMUON->Nodes())->Add(node);
544     node->cd();
545     sprintf(nameNode,"MUON%d",200+fId+1);
546     node = new TNode(nameNode,"Sens. Region Node",nameSense,0,0,0,"");
547     node->SetLineColor(kColorMUON);
548     node->cd();
549     Float_t dr=dx+rmin;
550     sprintf(nameNode,"MUON%d",300+fId+1);
551     nodeF = new TNode(nameNode,"Frame0",frMUON,dr, 0, 0,rot000,"");
552     nodeF->SetLineColor(kColorMUON);
553     node->cd();
554     sprintf(nameNode,"MUON%d",400+fId+1);
555     nodeF = new TNode(nameNode,"Frame1",frMUON,0 ,dr,0,rot090,"");
556     nodeF->SetLineColor(kColorMUON);
557     node->cd();
558     sprintf(nameNode,"MUON%d",500+fId+1);
559     nodeF = new TNode(nameNode,"Frame2",frMUON,-dr,0,0,rot180,"");
560     nodeF->SetLineColor(kColorMUON);
561     node  ->cd();
562     sprintf(nameNode,"MUON%d",600+fId+1);   
563     nodeF = new TNode(nameNode,"Frame3",frMUON,0,-dr,0,rot270,"");
564     nodeF->SetLineColor(kColorMUON);   
565   } else {
566     TBox *box;
567     
568     Float_t dx=0.95/fCx[3][1]/2;
569     Float_t dy=0.95/(Float_t(Npy()))/2;
570     Float_t x0,y0,x1,y1;
571     Float_t xc=0.5;
572     Float_t yc=0.5;
573     
574     for (Int_t iy=1; iy<Npy(); iy++) {
575       for (Int_t isec=0; isec<4; isec++) {
576         if (isec==0) {
577           x0=0;
578           x1=fCx[isec][iy]*dx;
579         } else {
580           x0=fCx[isec-1][iy]*dx;
581           x1=fCx[isec][iy]*dx;
582         }
583         y0=Float_t(iy-1)*dy;
584         y1=y0+dy;
585         box=new TBox(x0+xc,y0+yc,x1+xc,y1+yc);
586         box->SetFillColor(isec+1);
587         box->Draw();
588         
589         box=new TBox(-x1+xc,y0+yc,-x0+xc,y1+yc);
590         box->SetFillColor(isec+1);
591         box->Draw();
592         
593         box=new TBox(x0+xc,-y1+yc,x1+xc,-y0+yc);
594         box->SetFillColor(isec+1);
595         box->Draw();
596         
597         box=new TBox(-x1+xc,-y1+yc,-x0+xc,-y0+yc);
598         box->SetFillColor(isec+1);
599         box->Draw();
600       }
601     }
602   }
603 }
604 void AliMUONSegmentationV01::SetCorrFunc(Int_t isec, TF1* func)
605 {
606 // Set the correction function
607     fCorrA->AddAt(func,isec);
608 }
609
610 TF1* AliMUONSegmentationV01::CorrFunc(Int_t isec) const
611
612 // Get correction function
613   //PH    return (TF1*) (*fCorrA)[isec];
614     return (TF1*) fCorrA->At(isec);
615 }
616
617 AliMUONSegmentationV01& AliMUONSegmentationV01::operator 
618 =(const AliMUONSegmentationV01 & /*rhs*/)
619 {
620 // Dummy assignment operator
621     return *this;
622 }
623