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