MUON/AliMUONv01.cxx

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