MUON/AliMUONSegResV1.cxx

   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 //___________________________________________
  17 ClassImp(AliMUONsegmentationV1)
  18
  19 AliMUONsegmentationV1::AliMUONsegmentationV1()
  20   // initizalize the class with default settings
  21 {
  22 fNzone=1;
  23 fDAnod=0.0;
  24 fDpx=0.0;
  25 fDpx=0.0; // forces crash if not initialized by user
  26 fNZoneCut[0]=0;
  27 fSensOffset=0;
  28 }
  29
  30
  31 void 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
  44 void AliMUONsegmentationV1::DefaultCut(void)
  45 {
  46 SetNzone(3);
  47 AddCut(0,5*6,18*8);
  48 AddCut(0,9*6,15*8);
  49 AddCut(0,11*6,12*8);
  50 AddCut(0,12*6,9*8);
  51 AddCut(0,13*6,6*8);
  52 AddCut(1,6*6,20*12);
  53 AddCut(1,12*6,18*12);
  54 AddCut(1,15*6,15*12);
  55 AddCut(1,18*6,12*12);
  56 AddCut(1,21*6,9*12);
  57 SetSensOffset(3.0);
  58 SetDAnod(0.325);
  59 }
  60
  61 Int_t AliMUONsegmentationV1::GetiAnod(Float_t xhit)
  62 {
  63 Int_t kwire=Int_t((TMath::Abs(xhit)-fSensOffset)/fDAnod)+1;
  64 return (xhit>0) ? kwire : -kwire ;
  65 }
  66
  67 Float_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
  74 void AliMUONsegmentationV1::SetPADSIZ(Float_t p1, Float_t p2)
  75 {
  76   fDpx=p1;
  77   fDpy=p2;
  78 }
  79
  80 void 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
  89 void AliMUONsegmentationV1::
  90 GetPadCxy(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
  98 void 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
 100 if (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);
 102 fZoneX[Zone][fNZoneCut[Zone]] = nX;
 103 fZoneY[Zone][fNZoneCut[Zone]] = nY;
 104 fNZoneCut[Zone]++;
 105 }
 106
 107 Int_t AliMUONsegmentationV1::GetZone(Float_t X, Float_t Y)
 108 {
 109 Int_t iX, iY;
 110 GetPadIxy(X,Y,iX,iY);
 111 return GetZone( iX , iY );
 112 }
 113
 114 Int_t AliMUONsegmentationV1::GetZone(Int_t nX, Int_t nY)
 115 {// Beware : first pad begins at 1 !!
 116 Int_t aX =  TMath::Abs(nX);
 117 Int_t aY =  TMath::Abs(nY);
 118 Int_t zone=fNzone-1;
 119 for (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   }
 128 return zone;
 129 }
 130
 131 void AliMUONsegmentationV1::
 132 SetHit(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
 141 void AliMUONsegmentationV1::
 142 SetPad(Int_t ix, Int_t iy)
 143 {
 144     GetPadCxy(ix,iy,fx,fy);
 145 }
 146
 147
 148 void AliMUONsegmentationV1::SetPadCoord(Int_t iX, Int_t iY)
 149 {
 150 GetPadCxy(iX,iY,fx,fy);
 151 Float_t radius2;
 152 if ( ( (radius2=fx*fx+fy*fy) > frSensMax2 || radius2 < frSensMin2 )
 153      && MorePads() )
 154   NextPad();
 155 }
 156
 157 void 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
 192 void 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
 206 Int_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
 217 Int_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 {
 229 if (iY%3==1) return 0;
 230 return (iX%6)/3+1;
 231 }
 232
 233 Int_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 {
 245 if (iY%3==1) return 0;
 246 return (iX%9)/3+1 - (iY%3==2 && iX%3==0);
 247 }
 248
 249 Int_t AliMUONsegmentationV1::NParallel2(Int_t , 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 {
 259 if (iY%3==1) return 1;
 260 return 2;
 261 }
 262
 263 Int_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 {
 275 if (iY%3==1) return 1;
 276 if (iY%3==2 && iX%9==0) return 1;
 277 return 3 - (iY%3==2 && iX%3==0);
 278 }
 279
 280
 281 Int_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 {
 285 Int_t wix = TMath::Abs(trueX)-1;
 286 Int_t wiy = TMath::Abs(trueY)-1;
 287 Int_t zone = GetZone(trueX,trueY);
 288 Int_t par3;
 289 switch (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  }
 302 return -1;
 303 }
 304
 305 Int_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
 310 Int_t AliMUONsegmentationV1::ISector()
 311 // This function is of no use for this kind of segmentation.
 312 {
 313 return GetZone(fix,fiy);
 314 }
 315
 316 void AliMUONsegmentationV1::SigGenInit(Float_t x,Float_t y,Float_t )
 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
 327 Int_t AliMUONsegmentationV1::SigGenCond(Float_t x,Float_t y,Float_t )
 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
 344 void AliMUONsegmentationV1::
 345 IntegrationLimits(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
 353 void AliMUONsegmentationV1::GetNParallelAndOffset(Int_t iX, Int_t iY,Int_t
 354 *Nparallel, Int_t* Offset)
 355 {
 356 Int_t wix = TMath::Abs(iX)-1;
 357 Int_t wiy = TMath::Abs(iY)-1;
 358 Int_t zone = GetZone(iX,iY);
 359 switch (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
 380 Float_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 {
 385 Int_t nPara,offset;
 386 GetNParallelAndOffset(iX,iY,&nPara,&offset);
 387 Float_t d2min=1E10;
 388 for (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   }
 399 return d2min;
 400 }
 401
 402 void 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 {
 409 Int_t nTot=0;
 410 for (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
 432 void AliMUONsegmentationV1::
 433 NeighboursNonDiag(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 {
 437 Int_t nParallel, offset;
 438 GetNParallelAndOffset(iX,iY,&nParallel,&offset);
 439 //
 440 // now fill raw list of neighbours
 441 *Nlist=4*nParallel;
 442 Xlist[0]=Xlist[1]=iX;Xlist[2]=iX-1;Xlist[3]=iX+1;
 443 Ylist[0]=iY-1;Ylist[1]=iY+1;Ylist[2]=Ylist[3]=iY;
 444 if (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   }
 452 CleanNeighbours(Nlist,Xlist,Ylist);
 453 }
 454
 455 void AliMUONsegmentationV1::
 456 NeighboursDiag(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 {
 460 Int_t nParallel, offset;
 461 GetNParallelAndOffset(iX,iY,&nParallel,&offset);
 462 //
 463 // now fill raw list of neighbours
 464 *Nlist=0;
 465 for (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       }
 475 CleanNeighbours(Nlist,Xlist,Ylist);
 476 }
 477
 478 void 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
 483
 484 void AliMUONsegmentationV1::GiveTestPoints(Int_t &n, Float_t *x, Float_t *y)
 485 {
 486     n=1;
 487     x[0]=(TMath::Sqrt(frSensMax2)-TMath::Sqrt(frSensMin2))/2/TMath::Sqrt(2.);
 488     y[0]=x[0];
 489 }
 490