Code from MUON-dev joined
[u/mrichter/AliRoot.git] / MUON / AliMUONSegResV01.cxx
CommitLineData
4c039060 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
a897a37a 16/////////////////////////////////////////////////////
17// Segmentation and Response classes version 01 //
18/////////////////////////////////////////////////////
19
20#include <TTUBE.h>
21#include <TNode.h>
22#include <TBox.h>
23#include <TRandom.h>
24
25#include "AliMUONSegResV01.h"
26#include "AliRun.h"
27#include "AliMC.h"
28#include "iostream.h"
29
30//___________________________________________
a9e2aefa 31ClassImp(AliMUONSegmentationV01)
a897a37a 32
a897a37a 33
a9e2aefa 34AliMUONSegmentationV01::AliMUONSegmentationV01()
a897a37a 35{
36 fNsec=4;
37 fRSec.Set(fNsec);
38 fNDiv.Set(fNsec);
39 fDpxD.Set(fNsec);
40 fRSec[0]=fRSec[1]=fRSec[2]=fRSec[3]=0;
41 fNDiv[0]=fNDiv[1]=fNDiv[2]=fNDiv[3]=0;
42 fDpxD[0]=fDpxD[1]=fDpxD[2]=fDpxD[3]=0;
43 fCorr = new TObjArray(3);
44 (*fCorr)[0]=0;
45 (*fCorr)[1]=0;
46 (*fCorr)[2]=0;
47}
a9e2aefa 48
49Float_t AliMUONSegmentationV01::Dpx(Int_t isec)
50{
51//
52// Returns x-pad size for given sector isec
53 return fDpxD[isec];
54}
55
56Float_t AliMUONSegmentationV01::Dpy(Int_t isec)
57{
58//
59// Returns y-pad size for given sector isec
60 return fDpy;
61}
a897a37a 62
a9e2aefa 63void AliMUONSegmentationV01::SetSegRadii(Float_t r[4])
a897a37a 64{
a9e2aefa 65//
66// Set the radii of the segmentation zones
a897a37a 67 for (Int_t i=0; i<4; i++) {
68 fRSec[i]=r[i];
a9e2aefa 69 printf("\n R %d %f \n",i,fRSec[i]);
a897a37a 70
71 }
72}
73
74
a9e2aefa 75void AliMUONSegmentationV01::SetPadDivision(Int_t ndiv[4])
a897a37a 76{
77//
78// Defines the pad size perp. to the anode wire (y) for different sectors.
a9e2aefa 79// Pad sizes are defined as integral fractions ndiv of a basis pad size
80// fDpx
81//
a897a37a 82 for (Int_t i=0; i<4; i++) {
83 fNDiv[i]=ndiv[i];
a9e2aefa 84 printf("\n Ndiv %d %d \n",i,fNDiv[i]);
a897a37a 85 }
86 ndiv[0]=ndiv[1];
87}
88
89
a9e2aefa 90void AliMUONSegmentationV01::Init(AliMUONChamber* Chamber)
a897a37a 91{
92//
93// Fill the arrays fCx (x-contour) and fNpxS (ix-contour) for each sector
94// These arrays help in converting from real to pad co-ordinates and
a9e2aefa 95// vice versa.
96// This version approximates concentric segmentation zones
a897a37a 97//
98 Int_t isec;
a9e2aefa 99 printf("\n Initialise segmentation v01 -- test !!!!!!!!!!!!!! \n");
a897a37a 100 fNpy=Int_t(fRSec[fNsec-1]/fDpy)+1;
101
102 fDpxD[fNsec-1]=fDpx;
103 if (fNsec > 1) {
104 for (Int_t i=fNsec-2; i>=0; i--){
105 fDpxD[i]=fDpxD[fNsec-1]/fNDiv[i];
a9e2aefa 106 printf("\n test ---dx %d %f \n",i,fDpxD[i]);
a897a37a 107 }
108 }
109//
110// fill the arrays defining the pad segmentation boundaries
111 Float_t ry;
112 Int_t dnx;
113 Int_t add;
114//
115// loop over sections
116 for(isec=0; isec<fNsec; isec++) {
117//
118// loop over pads along the aode wires
119 for (Int_t iy=1; iy<=fNpy; iy++) {
120//
121 Float_t x=iy*fDpy-fDpy/2;
122 if (x > fRSec[isec]) {
123 fNpxS[isec][iy]=0;
a9e2aefa 124 fCx[isec][iy]=0;
a897a37a 125 } else {
126 ry=TMath::Sqrt(fRSec[isec]*fRSec[isec]-x*x);
127 if (isec > 1) {
128 dnx= Int_t((ry-fCx[isec-1][iy])/fDpxD[isec]);
129 if (isec < fNsec-1) {
130 if (TMath::Odd((Long_t)dnx)) dnx++;
131 }
132 fNpxS[isec][iy]=fNpxS[isec-1][iy]+dnx;
133 fCx[isec][iy]=fCx[isec-1][iy]+dnx*fDpxD[isec];
134 } else if (isec == 1) {
135 dnx= Int_t((ry-fCx[isec-1][iy])/fDpxD[isec]);
136 fNpxS[isec][iy]=fNpxS[isec-1][iy]+dnx;
137 add=4 - (fNpxS[isec][iy])%4;
138 if (add < 4) fNpxS[isec][iy]+=add;
139 dnx=fNpxS[isec][iy]-fNpxS[isec-1][iy];
140 fCx[isec][iy]=fCx[isec-1][iy]+dnx*fDpxD[isec];
141 } else {
142 dnx=Int_t(ry/fDpxD[isec]);
143 fNpxS[isec][iy]=dnx;
144 fCx[isec][iy]=dnx*fDpxD[isec];
145 }
146 }
147 } // y-pad loop
148 } // sector loop
a897a37a 149}
150
a9e2aefa 151Int_t AliMUONSegmentationV01::Sector(Int_t ix, Int_t iy)
a897a37a 152{
a9e2aefa 153// Returns sector number for given pad position
154//
a897a37a 155 Int_t absix=TMath::Abs(ix);
156 Int_t absiy=TMath::Abs(iy);
157 Int_t isec=0;
158 for (Int_t i=0; i<fNsec; i++) {
159 if (absix<=fNpxS[i][absiy]){
160 isec=i;
161 break;
162 }
163 }
164 return isec;
165}
166
a9e2aefa 167void AliMUONSegmentationV01::
168GetPadIxy(Float_t x, Float_t y, Int_t &ix, Int_t &iy)
a897a37a 169{
a9e2aefa 170// Returns pad coordinates (ix,iy) for given real coordinates (x,y)
a897a37a 171//
172 iy = (y>0)? Int_t(y/fDpy)+1 : Int_t(y/fDpy)-1;
173 if (iy > fNpy) iy= fNpy;
174 if (iy < -fNpy) iy=-fNpy;
175//
176// Find sector isec
177 Int_t isec=-1;
178 Float_t absx=TMath::Abs(x);
179 Int_t absiy=TMath::Abs(iy);
180 for (Int_t i=0; i < fNsec; i++) {
181 if (absx <= fCx[i][absiy]) {
182 isec=i;
183 break;
184 }
185 }
186 if (isec>0) {
187 ix= Int_t((absx-fCx[isec-1][absiy])/fDpxD[isec])
188 +fNpxS[isec-1][absiy]+1;
189 } else if (isec == 0) {
190 ix= Int_t(absx/fDpxD[isec])+1;
191 } else {
192 ix=fNpxS[fNsec-1][absiy]+1;
193 }
a897a37a 194 ix = (x>0) ? ix:-ix;
195}
196
a9e2aefa 197void AliMUONSegmentationV01::
a897a37a 198GetPadCxy(Int_t ix, Int_t iy, Float_t &x, Float_t &y)
199{
a9e2aefa 200// Returns real coordinates (x,y) for given pad coordinates (ix,iy)
a897a37a 201//
202 y = (iy>0) ? Float_t(iy*fDpy)-fDpy/2. : Float_t(iy*fDpy)+fDpy/2.;
203//
204// Find sector isec
a9e2aefa 205 Int_t isec=AliMUONSegmentationV01::Sector(ix,iy);
a897a37a 206//
207 Int_t absix=TMath::Abs(ix);
208 Int_t absiy=TMath::Abs(iy);
209 if (isec) {
210 x=fCx[isec-1][absiy]+(absix-fNpxS[isec-1][absiy])*fDpxD[isec];
211 x=(ix>0) ? x-fDpxD[isec]/2 : -x+fDpxD[isec]/2;
212 } else {
213 x=y=0;
214 }
215}
216
a9e2aefa 217void AliMUONSegmentationV01::
a897a37a 218SetPad(Int_t ix, Int_t iy)
219{
a9e2aefa 220 //
221 // Sets virtual pad coordinates, needed for evaluating pad response
222 // outside the tracking program
a897a37a 223 GetPadCxy(ix,iy,fx,fy);
224 fSector=Sector(ix,iy);
225}
226
227
a9e2aefa 228void AliMUONSegmentationV01::
a897a37a 229FirstPad(Float_t xhit, Float_t yhit, Float_t dx, Float_t dy)
230{
a9e2aefa 231// Initialises iteration over pads for charge distribution algorithm
232//
a897a37a 233 //
234 // Find the wire position (center of charge distribution)
235 Float_t x0a=GetAnod(xhit);
236 fxhit=x0a;
237 fyhit=yhit;
238
239 //
240 // and take fNsigma*sigma around this center
241 Float_t x01=x0a - dx;
242 Float_t x02=x0a + dx;
243 Float_t y01=yhit - dy;
244 Float_t y02=yhit + dy;
245 //
246 // find the pads over which the charge distributes
247 GetPadIxy(x01,y01,fixmin,fiymin);
248 GetPadIxy(x02,y02,fixmax,fiymax);
249 fxmin=x01;
250 fxmax=x02;
251 fymin=y01;
252 fymax=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
a9e2aefa 264void AliMUONSegmentationV01::NextPad()
a897a37a 265{
a9e2aefa 266// Stepper for the iteration over pads
267//
268// Step to next pad in the integration region
a897a37a 269 //
270 // Step to next pad in integration region
271 Float_t xc,yc;
272 Int_t iyc;
273
274// step from left to right
275 if (fx < fxmax && fx != 0) {
276 if (fix==-1) fix++;
277 fix++;
278// step up
279 } else if (fiy != fiymax) {
280 if (fiy==-1) fiy++;
281 fiy++;
282// get y-position of next row (yc), xc not used here
283 GetPadCxy(fix,fiy,xc,yc);
a9e2aefa 284// get x-pad coordiante for first pad in row (fix)
a897a37a 285 GetPadIxy(fxmin,yc,fix,iyc);
286 } else {
287 printf("\n Error: Stepping outside integration region\n ");
288 }
289 GetPadCxy(fix,fiy,fx,fy);
290 fSector=Sector(fix,fiy);
291 if (MorePads() &&
a9e2aefa 292 (fSector ==-1 || fSector==0))
a897a37a 293 NextPad();
a897a37a 294}
295
a9e2aefa 296Int_t AliMUONSegmentationV01::MorePads()
297// Stopping condition for the iterator over pads
a897a37a 298//
299// Are there more pads in the integration region
300{
301 if ((fx >= fxmax && fiy >= fiymax) || fy==0) {
302 return 0;
303 } else {
304 return 1;
305 }
306}
307
a9e2aefa 308void AliMUONSegmentationV01::
a897a37a 309IntegrationLimits(Float_t& x1,Float_t& x2,Float_t& y1, Float_t& y2)
310{
a9e2aefa 311// Returns integration limits for current pad
312//
a897a37a 313 x1=fxhit-fx-Dpx(fSector)/2.;
314 x2=x1+Dpx(fSector);
315 y1=fyhit-fy-Dpy(fSector)/2.;
316 y2=y1+Dpy(fSector);
317}
318
a9e2aefa 319void AliMUONSegmentationV01::
a897a37a 320Neighbours(Int_t iX, Int_t iY, Int_t* Nlist, Int_t Xlist[10], Int_t Ylist[10])
321{
a9e2aefa 322// Returns list of next neighbours for given Pad (iX, iY)
323//
a897a37a 324 const Float_t epsilon=fDpy/1000;
325
326 Float_t x,y;
327 Int_t ixx, iyy, isec1;
328//
a9e2aefa 329 Int_t isec0=AliMUONSegmentationV01::Sector(iX,iY);
a897a37a 330 Int_t i=0;
331//
332// step right
333 Xlist[i]=iX+1;
a9e2aefa 334 if (Xlist[i]==0) Xlist[i]++;
a897a37a 335 Ylist[i++]=iY;
336//
337// step left
338 Xlist[i]=iX-1;
a9e2aefa 339 if (Xlist[i]==0) Xlist[i]--;
a897a37a 340 Ylist[i++]=iY;
341//
342// step up
a9e2aefa 343 AliMUONSegmentationV01::GetPadCxy(iX,iY,x,y);
344 AliMUONSegmentationV01::GetPadIxy(x+epsilon,y+fDpy,ixx,iyy);
a897a37a 345 Xlist[i]=ixx;
a9e2aefa 346 Ylist[i++]=iyy;
347 isec1=AliMUONSegmentationV01::Sector(ixx,iyy);
a897a37a 348 if (isec1==isec0) {
349//
350// no sector boundary crossing
a9e2aefa 351// Xlist[i]=ixx+1;
352// Ylist[i++]=iY+1;
a897a37a 353
a9e2aefa 354// Xlist[i]=ixx-1;
355// Ylist[i++]=iY+1;
a897a37a 356 } else if (isec1 < isec0) {
357// finer segmentation
a9e2aefa 358// Xlist[i]=ixx+1;
359// Ylist[i++]=iY+1;
a897a37a 360
361 Xlist[i]=ixx-1;
a9e2aefa 362 Ylist[i++]=iyy;
a897a37a 363
a9e2aefa 364// Xlist[i]=ixx-2;
365// Ylist[i++]=iY+1;
a897a37a 366 } else {
367// coarser segmenation
a9e2aefa 368/*
a897a37a 369 if (TMath::Odd(iX-fNpxS[isec1-1][iY+1])) {
370 Xlist[i]=ixx-1;
371 Ylist[i++]=iY+1;
372 } else {
373 Xlist[i]=ixx+1;
374 Ylist[i++]=iY+1;
375 }
a9e2aefa 376*/
a897a37a 377 }
a9e2aefa 378
a897a37a 379//
380// step down
a9e2aefa 381 AliMUONSegmentationV01::GetPadCxy(iX,iY,x,y);
382 AliMUONSegmentationV01::GetPadIxy(x+epsilon,y-fDpy,ixx,iyy);
a897a37a 383 Xlist[i]=ixx;
a9e2aefa 384 Ylist[i++]=iyy;
385 isec1=AliMUONSegmentationV01::Sector(ixx,iyy);
a897a37a 386 if (isec1==isec0) {
387//
388// no sector boundary crossing
a9e2aefa 389/*
a897a37a 390 Xlist[i]=ixx+1;
391 Ylist[i++]=iY-1;
392
393 Xlist[i]=ixx-1;
394 Ylist[i++]=iY-1;
a9e2aefa 395*/
a897a37a 396 } else if (isec1 < isec0) {
397// finer segmentation
a9e2aefa 398// Xlist[i]=ixx+1;
399// Ylist[i++]=iY-1;
a897a37a 400
401 Xlist[i]=ixx-1;
a9e2aefa 402 Ylist[i++]=iyy;
a897a37a 403
a9e2aefa 404// Xlist[i]=ixx-2;
405// Ylist[i++]=iY-1;
a897a37a 406 } else {
407// coarser segmentation
a9e2aefa 408/*
a897a37a 409 if (TMath::Odd(iX-fNpxS[isec1-1][iY-1])) {
410 Xlist[i]=ixx-1;
411 Ylist[i++]=iY-1;
412 } else {
413 Xlist[i]=ixx+1;
414 Ylist[i++]=iY-1;
415 }
a9e2aefa 416*/
a897a37a 417 }
418 *Nlist=i;
419}
420
a9e2aefa 421void AliMUONSegmentationV01::GiveTestPoints(Int_t &n, Float_t *x, Float_t *y)
a897a37a 422{
a9e2aefa 423// Returns test point on the pad plane.
424// Used during determination of the segmoid correction of the COG-method
425
a897a37a 426 n=3;
427 x[0]=(fRSec[0]+fRSec[1])/2/TMath::Sqrt(2.);
428 y[0]=x[0];
429 x[1]=(fRSec[1]+fRSec[2])/2/TMath::Sqrt(2.);
430 y[1]=x[1];
431 x[2]=(fRSec[2]+fRSec[3])/2/TMath::Sqrt(2.);
432 y[2]=x[2];
433}
434
a9e2aefa 435void AliMUONSegmentationV01::Draw()
a897a37a 436{
a9e2aefa 437// Draws the segmentation zones
438//
a897a37a 439 TBox *box;
440
441 Float_t dx=0.95/fCx[3][1]/2;
442 Float_t dy=0.95/(Float_t(Npy()))/2;
443 Float_t x0,y0,x1,y1;
444 Float_t xc=0.5;
445 Float_t yc=0.5;
446
447 for (Int_t iy=1; iy<Npy(); iy++)
448 {
449 for (Int_t isec=0; isec<4; isec++) {
450 if (isec==0) {
451 x0=0;
452 x1=fCx[isec][iy]*dx;
453 } else {
454 x0=fCx[isec-1][iy]*dx;
455 x1=fCx[isec][iy]*dx;
456 }
457 y0=Float_t(iy-1)*dy;
458 y1=y0+dy;
459 box=new TBox(x0+xc,y0+yc,x1+xc,y1+yc);
460 box->SetFillColor(isec+1);
461 box->Draw();
462
463 box=new TBox(-x1+xc,y0+yc,-x0+xc,y1+yc);
464 box->SetFillColor(isec+1);
465 box->Draw();
466
467 box=new TBox(x0+xc,-y1+yc,x1+xc,-y0+yc);
468 box->SetFillColor(isec+1);
469 box->Draw();
470
471 box=new TBox(-x1+xc,-y1+yc,-x0+xc,-y0+yc);
472 box->SetFillColor(isec+1);
473 box->Draw();
474 }
475 }
476}
a9e2aefa 477void AliMUONSegmentationV01::SetCorrFunc(Int_t isec, TF1* func)
a897a37a 478{
479 (*fCorr)[isec]=func;
a897a37a 480}
481
a9e2aefa 482TF1* AliMUONSegmentationV01::CorrFunc(Int_t isec)
a897a37a 483{
484 return (TF1*) (*fCorr)[isec];
485}
486