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