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