MUON/AliMUONClusterFinder.cxx

   1 #include "AliMUONClusterFinder.h"
   2 #include "TTree.h"
   3 #include "AliRun.h"
   4 #include <TCanvas.h>
   5 #include <TH1.h>
   6 #include <TPad.h>
   7 #include <TGraph.h>
   8 #include <TPostScript.h>
   9 #include <TMinuit.h>
  10
  11 //----------------------------------------------------------
  12 static AliMUONsegmentation* gSegmentation;
  13 static AliMUONresponse*     gResponse;
  14 static Int_t                gix[500];
  15 static Int_t                giy[500];
  16 static Float_t              gCharge[500];
  17 static Int_t                gNbins;
  18 static Int_t                gFirst=kTRUE;
  19 static TMinuit *gMyMinuit ;
  20 void fcn(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
  21 static Int_t                gChargeTot;
  22
  23 //----------------------------------------------------------
  24
  25 ClassImp(AliMUONClusterFinder)
  26
  27     AliMUONClusterFinder::AliMUONClusterFinder
  28 (AliMUONsegmentation *segmentation, AliMUONresponse *response,
  29  TClonesArray *digits, Int_t chamber)
  30 {
  31     fSegmentation=segmentation;
  32     fResponse=response;
  33
  34     fDigits=digits;
  35     fNdigits = fDigits->GetEntriesFast();
  36     fChamber=chamber;
  37     fRawClusters=new TClonesArray("AliMUONRawCluster",10000);
  38     fNRawClusters=0;
  39     fCogCorr = 0;
  40     SetNperMax();
  41     SetClusterSize();
  42     SetDeclusterFlag();
  43     fNPeaks=-1;
  44 }
  45
  46     AliMUONClusterFinder::AliMUONClusterFinder()
  47 {
  48     fSegmentation=0;
  49     fResponse=0;
  50
  51     fDigits=0;
  52     fNdigits = 0;
  53     fChamber=-1;
  54     fRawClusters=new TClonesArray("AliMUONRawCluster",10000);
  55     fNRawClusters=0;
  56     fHitMap = 0;
  57     fCogCorr = 0;
  58     SetNperMax();
  59     SetClusterSize();
  60     SetDeclusterFlag();
  61     fNPeaks=-1;
  62 }
  63
  64 void AliMUONClusterFinder::AddRawCluster(const AliMUONRawCluster c)
  65 {
  66   //
  67   // Add a raw cluster copy to the list
  68   //
  69     AliMUON *MUON=(AliMUON*)gAlice->GetModule("MUON");
  70     MUON->AddRawCluster(fChamber,c);
  71     fNRawClusters++;
  72 }
  73
  74
  75
  76 void AliMUONClusterFinder::Decluster(AliMUONRawCluster *cluster)
  77 {
  78 //    AliMUONdigit *dig;
  79 //    Int_t q;
  80
  81
  82     Int_t mul = cluster->fMultiplicity;
  83 //    printf("Decluster - multiplicity   %d \n",mul);
  84
  85     if (mul == 1 || mul ==2) {
  86 //      printf("\n Nothing special for 1- and 2-clusters \n");
  87 //
  88 // Nothing special for 1- and 2-clusters
  89         if (fNPeaks != 0) {
  90             cluster->fNcluster[0]=fNPeaks;
  91             cluster->fNcluster[1]=0;
  92         }
  93         AddRawCluster(*cluster);
  94         fNPeaks++;
  95     } else if (mul ==3) {
  96 //
  97 // 3-cluster, check topology
  98 //      printf("\n 3-cluster, check topology \n");
  99         if (fDeclusterFlag) {
 100           if (Centered(cluster)) {
 101             // ok, cluster is centered
 102             //      printf("\n ok, cluster is centered \n");
 103           } else {
 104             // cluster is not centered, split into 2+1
 105             //      printf("\n cluster is not centered, split into 2+1 \n");
 106           }
 107         } else {
 108           if (fNPeaks != 0) {
 109             cluster->fNcluster[0]=fNPeaks;
 110             cluster->fNcluster[1]=0;
 111           }
 112           AddRawCluster(*cluster);
 113           fNPeaks++;
 114         }
 115     } else {
 116       //if (mul < 12) {
 117           //      printf("Decluster - multiplicity > 45   %d \n",mul);
 118           //printf("Decluster - multiplicity < 25   %d \n",mul);
 119 //
 120 // 4-and more-pad clusters
 121 //
 122       if (mul <= fClusterSize) {
 123         if (fDeclusterFlag) {
 124             SplitByLocalMaxima(cluster);
 125         } else {
 126             if (fNPeaks != 0) {
 127                 cluster->fNcluster[0]=fNPeaks;
 128                 cluster->fNcluster[1]=0;
 129             }
 130             AddRawCluster(*cluster);
 131             fNPeaks++;
 132         }
 133       }
 134         //}
 135     } // multiplicity
 136 }
 137
 138
 139 Bool_t AliMUONClusterFinder::Centered(AliMUONRawCluster *cluster)
 140 {
 141     AliMUONdigit* dig;
 142     dig= (AliMUONdigit*)fDigits->UncheckedAt(cluster->fIndexMap[0]);
 143     Int_t ix=dig->fPadX;
 144     Int_t iy=dig->fPadY;
 145     Int_t nn;
 146     Int_t X[kMaxNeighbours], Y[kMaxNeighbours], XN[kMaxNeighbours], YN[kMaxNeighbours];
 147
 148     fSegmentation->Neighbours(ix,iy,&nn,X,Y);
 149     Int_t nd=0;
 150     for (Int_t i=0; i<nn; i++) {
 151         if (fHitMap->TestHit(X[i],Y[i]) == used) {
 152             XN[nd]=X[i];
 153             YN[nd]=Y[i];
 154             nd++;
 155         }
 156     }
 157     if (nd==2) {
 158 //
 159 // cluster is centered !
 160         if (fNPeaks != 0) {
 161             cluster->fNcluster[0]=fNPeaks;
 162             cluster->fNcluster[1]=0;
 163         }
 164         AddRawCluster(*cluster);
 165         fNPeaks++;
 166         return kTRUE;
 167     } else if (nd ==1) {
 168 //
 169 // Highest signal on an edge, split cluster into 2+1
 170 //
 171 // who is the neighbour ?
 172         Int_t nind=fHitMap->GetHitIndex(XN[0], YN[0]);
 173         Int_t i1= (nind==cluster->fIndexMap[1]) ? 1:2;
 174         Int_t i2= (nind==cluster->fIndexMap[1]) ? 2:1;
 175 //
 176 // 2-cluster
 177         AliMUONRawCluster cnew;
 178         if (fNPeaks == 0) {
 179             cnew.fNcluster[0]=-1;
 180             cnew.fNcluster[1]=fNRawClusters;
 181         } else {
 182             cnew.fNcluster[0]=fNPeaks;
 183             cnew.fNcluster[1]=0;
 184         }
 185         cnew.fMultiplicity=2;
 186         cnew.fIndexMap[0]=cluster->fIndexMap[0];
 187         cnew.fIndexMap[1]=cluster->fIndexMap[i1];
 188         FillCluster(&cnew);
 189         cnew.fClusterType=cnew.PhysicsContribution();
 190         AddRawCluster(cnew);
 191         fNPeaks++;
 192 //
 193 // 1-cluster
 194         cluster->fMultiplicity=1;
 195         cluster->fIndexMap[0]=cluster->fIndexMap[i2];
 196         cluster->fIndexMap[1]=0;
 197         cluster->fIndexMap[2]=0;
 198         FillCluster(cluster);
 199         if (fNPeaks != 0) {
 200             cluster->fNcluster[0]=fNPeaks;
 201             cluster->fNcluster[1]=0;
 202         }
 203         cluster->fClusterType=cluster->PhysicsContribution();
 204         AddRawCluster(*cluster);
 205         fNPeaks++;
 206         return kFALSE;
 207     } else {
 208         printf("\n Completely screwed up %d !! \n",nd);
 209
 210     }
 211
 212         return kFALSE;
 213 }
 214 void AliMUONClusterFinder::SplitByLocalMaxima(AliMUONRawCluster *c)
 215 {
 216     AliMUONdigit* dig[100], *digt;
 217     Int_t ix[100], iy[100], q[100];
 218     Float_t x[100], y[100];
 219     Int_t i; // loops over digits
 220     Int_t j; // loops over local maxima
 221     //    Float_t xPeak[2];
 222     //    Float_t yPeak[2];
 223     //    Int_t threshold=500;
 224     Int_t mul=c->fMultiplicity;
 225 //
 226 //  dump digit information into arrays
 227 //
 228     for (i=0; i<mul; i++)
 229     {
 230         dig[i]= (AliMUONdigit*)fDigits->UncheckedAt(c->fIndexMap[i]);
 231         ix[i]= dig[i]->fPadX;
 232         iy[i]= dig[i]->fPadY;
 233         q[i] = dig[i]->fSignal;
 234         fSegmentation->GetPadCxy(ix[i], iy[i], x[i], y[i]);
 235     }
 236 //
 237 //  Find local maxima
 238 //
 239     Bool_t IsLocal[100];
 240     Int_t NLocal=0;
 241     Int_t AssocPeak[100];
 242     Int_t IndLocal[100];
 243     Int_t nn;
 244     Int_t X[kMaxNeighbours], Y[kMaxNeighbours];
 245     for (i=0; i<mul; i++) {
 246         fSegmentation->Neighbours(ix[i], iy[i], &nn, X, Y);
 247         IsLocal[i]=kTRUE;
 248         for (j=0; j<nn; j++) {
 249             if (fHitMap->TestHit(X[j],Y[j])==empty) continue;
 250             digt=(AliMUONdigit*) fHitMap->GetHit(X[j], Y[j]);
 251             if (digt->fSignal > q[i]) {
 252                 IsLocal[i]=kFALSE;
 253                 break;
 254 //
 255 // handle special case of neighbouring pads with equal signal
 256             } else if (digt->fSignal == q[i]) {
 257                 if (NLocal >0) {
 258                     for (Int_t k=0; k<NLocal; k++) {
 259                         if (X[j]==ix[IndLocal[k]] && Y[j]==iy[IndLocal[k]]){
 260                             IsLocal[i]=kFALSE;
 261                         }
 262                     }
 263                 }
 264             }
 265         } // loop over next neighbours
 266         // Maxima should not be on the edge
 267         if (IsLocal[i]) {
 268             IndLocal[NLocal]=i;
 269             NLocal++;
 270         }
 271     } // loop over all digits
 272 //    printf("Found %d local Maxima",NLocal);
 273 //
 274 // If only one local maximum found but multiplicity is high
 275 // take global maximum from the list of digits.
 276     if (NLocal==1 && mul>12) {
 277         Int_t nnew=0;
 278         for (i=0; i<mul; i++) {
 279             if (!IsLocal[i]) {
 280                 IndLocal[NLocal]=i;
 281                 IsLocal[i]=kTRUE;
 282                 NLocal++;
 283                 nnew++;
 284             }
 285             if (nnew==1) break;
 286         }
 287     }
 288
 289 // If number of local maxima is 2 try to fit a double gaussian
 290     if (NLocal==2) {
 291 //
 292 //  Initialise global variables for fit
 293         gFirst=1;
 294         gSegmentation=fSegmentation;
 295         gResponse    =fResponse;
 296         gNbins=mul;
 297
 298         for (i=0; i<mul; i++) {
 299             gix[i]=ix[i];
 300             giy[i]=iy[i];
 301             gCharge[i]=Float_t(q[i]);
 302         }
 303 //
 304         if (gFirst) {
 305           gFirst=kFALSE;
 306           gMyMinuit = new TMinuit(5);
 307         }
 308         gMyMinuit->SetFCN(fcn);
 309         gMyMinuit->mninit(5,10,7);
 310         Double_t arglist[20];
 311         Int_t ierflag=0;
 312         arglist[0]=1;
 313 //      gMyMinuit->mnexcm("SET ERR",arglist,1,ierflag);
 314 // Set starting values
 315         static Double_t vstart[5];
 316         vstart[0]=x[IndLocal[0]];
 317         vstart[1]=y[IndLocal[0]];
 318         vstart[2]=x[IndLocal[1]];
 319         vstart[3]=y[IndLocal[1]];
 320         vstart[4]=Float_t(q[IndLocal[0]])/
 321             Float_t(q[IndLocal[0]]+q[IndLocal[1]]);
 322 // lower and upper limits
 323         static Double_t lower[5], upper[5];
 324         Int_t isec=fSegmentation->Sector(ix[IndLocal[0]], iy[IndLocal[0]]);
 325         lower[0]=vstart[0]-fSegmentation->Dpx(isec)/2;
 326         lower[1]=vstart[1]-fSegmentation->Dpy(isec)/2;
 327 //      lower[1]=vstart[1];
 328
 329         upper[0]=lower[0]+fSegmentation->Dpx(isec);
 330         upper[1]=lower[1]+fSegmentation->Dpy(isec);
 331 //      upper[1]=vstart[1];
 332
 333         isec=fSegmentation->Sector(ix[IndLocal[1]], iy[IndLocal[1]]);
 334         lower[2]=vstart[2]-fSegmentation->Dpx(isec)/2;
 335         lower[3]=vstart[3]-fSegmentation->Dpy(isec)/2;
 336 //      lower[3]=vstart[3];
 337
 338         upper[2]=lower[2]+fSegmentation->Dpx(isec);
 339         upper[3]=lower[3]+fSegmentation->Dpy(isec);
 340 //      upper[3]=vstart[3];
 341
 342         lower[4]=0.;
 343         upper[4]=1.;
 344 // step sizes
 345         static Double_t step[5]={0.005, 0.03, 0.005, 0.03, 0.01};
 346
 347         gMyMinuit->mnparm(0,"x1",vstart[0],step[0],lower[0],upper[0],ierflag);
 348         gMyMinuit->mnparm(1,"y1",vstart[1],step[1],lower[1],upper[1],ierflag);
 349         gMyMinuit->mnparm(2,"x2",vstart[2],step[2],lower[2],upper[2],ierflag);
 350         gMyMinuit->mnparm(3,"y2",vstart[3],step[3],lower[3],upper[3],ierflag);
 351         gMyMinuit->mnparm(4,"a0",vstart[4],step[4],lower[4],upper[4],ierflag);
 352 // ready for minimisation
 353         gMyMinuit->SetPrintLevel(-1);
 354         gMyMinuit->mnexcm("SET OUT", arglist, 0, ierflag);
 355         arglist[0]= -1;
 356         arglist[1]= 0;
 357
 358         gMyMinuit->mnexcm("SET NOGR", arglist, 0, ierflag);
 359         gMyMinuit->mnexcm("SCAN", arglist, 0, ierflag);
 360         gMyMinuit->mnexcm("EXIT" , arglist, 0, ierflag);
 361 // Print results
 362 //      Double_t amin,edm,errdef;
 363 //      Int_t nvpar,nparx,icstat;
 364 //      gMyMinuit->mnstat(amin,edm,errdef,nvpar,nparx,icstat);
 365 //      gMyMinuit->mnprin(3,amin);
 366 // Get fitted parameters
 367
 368         Double_t xrec[2], yrec[2], qfrac;
 369         TString chname;
 370         Double_t epxz, b1, b2;
 371         Int_t ierflg;
 372         gMyMinuit->mnpout(0, chname, xrec[0], epxz, b1, b2, ierflg);
 373         gMyMinuit->mnpout(1, chname, yrec[0], epxz, b1, b2, ierflg);
 374         gMyMinuit->mnpout(2, chname, xrec[1], epxz, b1, b2, ierflg);
 375         gMyMinuit->mnpout(3, chname, yrec[1], epxz, b1, b2, ierflg);
 376         gMyMinuit->mnpout(4, chname, qfrac,   epxz, b1, b2, ierflg);
 377         printf("\n %f %f %f %f %f\n", xrec[0], yrec[0], xrec[1], yrec[1],qfrac);
 378 //      delete gMyMinuit;
 379
 380
 381  //
 382  // One cluster for each maximum
 383  //
 384          for (j=0; j<2; j++) {
 385              AliMUONRawCluster cnew;
 386              if (fNPeaks == 0) {
 387                  cnew.fNcluster[0]=-1;
 388                  cnew.fNcluster[1]=fNRawClusters;
 389              } else {
 390                  cnew.fNcluster[0]=fNPeaks;
 391                  cnew.fNcluster[1]=0;
 392              }
 393              cnew.fMultiplicity=0;
 394              cnew.fX=Float_t(xrec[j]);
 395              cnew.fY=Float_t(yrec[j]);
 396              if (j==0) {
 397                  cnew.fQ=Int_t(gChargeTot*qfrac);
 398              } else {
 399                  cnew.fQ=Int_t(gChargeTot*(1-qfrac));
 400              }
 401              gSegmentation->SetHit(xrec[j],yrec[j]);
 402              for (i=0; i<mul; i++) {
 403                  cnew.fIndexMap[cnew.fMultiplicity]=c->fIndexMap[i];
 404                  gSegmentation->SetPad(gix[i], giy[i]);
 405                  Float_t q1=gResponse->IntXY(gSegmentation);
 406                  cnew.fContMap[cnew.fMultiplicity]=Float_t(q[i])/(q1*cnew.fQ);
 407                  cnew.fMultiplicity++;
 408              }
 409              FillCluster(&cnew,0);
 410              //printf("\n x,y %f %f ", cnew.fX, cnew.fY);
 411              cnew.fClusterType=cnew.PhysicsContribution();
 412              AddRawCluster(cnew);
 413              fNPeaks++;
 414          }
 415      }
 416
 417      Bool_t fitted=kTRUE;
 418
 419      if (NLocal !=2 || !fitted) {
 420  // Check if enough local clusters have been found,
 421  // if not add global maxima to the list
 422  //
 423          Int_t nPerMax;
 424          if (NLocal!=0) {
 425              nPerMax=mul/NLocal;
 426          } else {
 427              printf("\n Warning, no local maximum found \n");
 428              nPerMax=fNperMax+1;
 429          }
 430
 431          if (nPerMax > fNperMax) {
 432              Int_t nGlob=mul/fNperMax-NLocal+1;
 433              if (nGlob > 0) {
 434                  Int_t nnew=0;
 435                  for (i=0; i<mul; i++) {
 436                      if (!IsLocal[i]) {
 437                          IndLocal[NLocal]=i;
 438                          IsLocal[i]=kTRUE;
 439                          NLocal++;
 440                          nnew++;
 441                      }
 442                      if (nnew==nGlob) break;
 443                  }
 444              }
 445          }
 446  //
 447  // Associate hits to peaks
 448  //
 449          for (i=0; i<mul; i++) {
 450              Float_t dmin=1.E10;
 451              Float_t qmax=0;
 452              if (IsLocal[i]) continue;
 453              for (j=0; j<NLocal; j++) {
 454                  Int_t il=IndLocal[j];
 455                  Float_t d=TMath::Sqrt((x[i]-x[il])*(x[i]-x[il])
 456                                    +(y[i]-y[il])*(y[i]-y[il]));
 457                  Float_t ql=q[il];
 458  //
 459  // Select nearest peak
 460  //
 461                  if (d<dmin) {
 462                      dmin=d;
 463                      qmax=ql;
 464                      AssocPeak[i]=j;
 465                  } else if (d==dmin) {
 466  //
 467  // If more than one take highest peak
 468  //
 469                      if (ql>qmax) {
 470                          dmin=d;
 471                          qmax=ql;
 472                          AssocPeak[i]=j;
 473                      }
 474                  }
 475              }
 476          }
 477
 478
 479  //
 480  // One cluster for each maximum
 481  //
 482          for (j=0; j<NLocal; j++) {
 483              AliMUONRawCluster cnew;
 484              if (fNPeaks == 0) {
 485                  cnew.fNcluster[0]=-1;
 486                  cnew.fNcluster[1]=fNRawClusters;
 487              } else {
 488                  cnew.fNcluster[0]=fNPeaks;
 489                  cnew.fNcluster[1]=0;
 490              }
 491              cnew.fIndexMap[0]=c->fIndexMap[IndLocal[j]];
 492              cnew.fMultiplicity=1;
 493              for (i=0; i<mul; i++) {
 494                  if (IsLocal[i]) continue;
 495                  if (AssocPeak[i]==j) {
 496                      cnew.fIndexMap[cnew.fMultiplicity]=c->fIndexMap[i];
 497                      cnew.fMultiplicity++;
 498                  }
 499              }
 500              FillCluster(&cnew);
 501              cnew.fClusterType=cnew.PhysicsContribution();
 502              AddRawCluster(cnew);
 503              fNPeaks++;
 504          }
 505      }
 506
 507 }
 508
 509
 510 void  AliMUONClusterFinder::FillCluster(AliMUONRawCluster* c, Int_t flag)
 511 {
 512 //
 513 //  Completes cluster information starting from list of digits
 514 //
 515     AliMUONdigit* dig;
 516     Float_t x, y;
 517     Int_t  ix, iy;
 518     Float_t frac=0;
 519
 520     c->fPeakSignal=0;
 521     if (flag) {
 522         c->fX=0;
 523         c->fY=0;
 524         c->fQ=0;
 525     }
 526     //c->fQ=0;
 527
 528
 529     for (Int_t i=0; i<c->fMultiplicity; i++)
 530     {
 531         dig= (AliMUONdigit*)fDigits->UncheckedAt(c->fIndexMap[i]);
 532         ix=dig->fPadX+c->fOffsetMap[i];
 533         iy=dig->fPadY;
 534         Int_t q=dig->fSignal;
 535         if (dig->fPhysics >= dig->fSignal) {
 536           c->fPhysicsMap[i]=2;
 537         } else if (dig->fPhysics == 0) {
 538           c->fPhysicsMap[i]=0;
 539         } else  c->fPhysicsMap[i]=1;
 540 //
 541 //
 542 // peak signal and track list
 543         if (flag) {
 544            if (q>c->fPeakSignal) {
 545               c->fPeakSignal=q;
 546 /*
 547             c->fTracks[0]=dig->fTracks[0];
 548             c->fTracks[1]=dig->fTracks[1];
 549             c->fTracks[2]=dig->fTracks[2];
 550 */
 551               //c->fTracks[0]=dig->fTrack;
 552             c->fTracks[0]=dig->fHit;
 553             c->fTracks[1]=dig->fTracks[0];
 554             c->fTracks[2]=dig->fTracks[1];
 555            }
 556         } else {
 557            if (c->fContMap[i] > frac) {
 558               frac=c->fContMap[i];
 559               c->fPeakSignal=q;
 560 /*
 561             c->fTracks[0]=dig->fTracks[0];
 562             c->fTracks[1]=dig->fTracks[1];
 563             c->fTracks[2]=dig->fTracks[2];
 564 */
 565               //c->fTracks[0]=dig->fTrack;
 566             c->fTracks[0]=dig->fHit;
 567             c->fTracks[1]=dig->fTracks[0];
 568             c->fTracks[2]=dig->fTracks[1];
 569            }
 570         }
 571 //
 572         if (flag) {
 573             fSegmentation->GetPadCxy(ix, iy, x, y);
 574             c->fX += q*x;
 575             c->fY += q*y;
 576             c->fQ += q;
 577         }
 578
 579     } // loop over digits
 580
 581  if (flag) {
 582
 583      c->fX/=c->fQ;
 584      c->fX=fSegmentation->GetAnod(c->fX);
 585      c->fY/=c->fQ;
 586 //
 587 //  apply correction to the coordinate along the anode wire
 588 //
 589      x=c->fX;
 590      y=c->fY;
 591      fSegmentation->GetPadIxy(x, y, ix, iy);
 592      fSegmentation->GetPadCxy(ix, iy, x, y);
 593      Int_t isec=fSegmentation->Sector(ix,iy);
 594      TF1* CogCorr = fSegmentation->CorrFunc(isec-1);
 595
 596      if (CogCorr) {
 597          Float_t YonPad=(c->fY-y)/fSegmentation->Dpy(isec);
 598          c->fY=c->fY-CogCorr->Eval(YonPad, 0, 0);
 599      }
 600  }
 601 }
 602
 603
 604 void  AliMUONClusterFinder::FindCluster(Int_t i, Int_t j, AliMUONRawCluster &c){
 605 //
 606 //  Find clusters
 607 //
 608 //
 609 //  Add i,j as element of the cluster
 610 //
 611
 612     Int_t idx = fHitMap->GetHitIndex(i,j);
 613     AliMUONdigit* dig = (AliMUONdigit*) fHitMap->GetHit(i,j);
 614     Int_t q=dig->fSignal;
 615     if (q > TMath::Abs(c.fPeakSignal)) {
 616         c.fPeakSignal=q;
 617 /*
 618         c.fTracks[0]=dig->fTracks[0];
 619         c.fTracks[1]=dig->fTracks[1];
 620         c.fTracks[2]=dig->fTracks[2];
 621 */
 622         //c.fTracks[0]=dig->fTrack;
 623         c.fTracks[0]=dig->fHit;
 624         c.fTracks[1]=dig->fTracks[0];
 625         c.fTracks[2]=dig->fTracks[1];
 626     }
 627 //
 628 //  Make sure that list of digits is ordered
 629 //
 630     Int_t mu=c.fMultiplicity;
 631     c.fIndexMap[mu]=idx;
 632
 633     if (dig->fPhysics >= dig->fSignal) {
 634         c.fPhysicsMap[mu]=2;
 635     } else if (dig->fPhysics == 0) {
 636         c.fPhysicsMap[mu]=0;
 637     } else  c.fPhysicsMap[mu]=1;
 638
 639     if (mu > 0) {
 640         for (Int_t ind=mu-1; ind>=0; ind--) {
 641             Int_t ist=(c.fIndexMap)[ind];
 642             Int_t ql=((AliMUONdigit*)fDigits
 643                       ->UncheckedAt(ist))->fSignal;
 644             if (q>ql) {
 645                 c.fIndexMap[ind]=idx;
 646                 c.fIndexMap[ind+1]=ist;
 647             } else {
 648                 break;
 649             }
 650         }
 651     }
 652
 653     c.fMultiplicity++;
 654
 655     if (c.fMultiplicity >= 50 ) {
 656         printf("FindCluster - multiplicity >50  %d \n",c.fMultiplicity);
 657         c.fMultiplicity=49;
 658     }
 659
 660 // Prepare center of gravity calculation
 661     Float_t x, y;
 662     fSegmentation->GetPadCxy(i, j, x, y);
 663     c.fX += q*x;
 664     c.fY += q*y;
 665     c.fQ += q;
 666 // Flag hit as taken
 667     fHitMap->FlagHit(i,j);
 668 //
 669 //  Now look recursively for all neighbours
 670 //
 671     Int_t nn;
 672     Int_t Xlist[kMaxNeighbours], Ylist[kMaxNeighbours];
 673     fSegmentation->Neighbours(i,j,&nn,Xlist,Ylist);
 674     for (Int_t in=0; in<nn; in++) {
 675         Int_t ix=Xlist[in];
 676         Int_t iy=Ylist[in];
 677         if (fHitMap->TestHit(ix,iy)==unused) FindCluster(ix, iy, c);
 678     }
 679 }
 680
 681 //_____________________________________________________________________________
 682
 683 void AliMUONClusterFinder::FindRawClusters()
 684 {
 685   //
 686   // simple MUON cluster finder from digits -- finds neighbours and
 687   // fill the tree with raw clusters
 688   //
 689
 690     if (!fNdigits) return;
 691
 692     fHitMap = new AliMUONHitMapA1(fSegmentation, fDigits);
 693
 694
 695     AliMUONdigit *dig;
 696
 697     Int_t ndig;
 698     Int_t nskip=0;
 699     Int_t ncls=0;
 700     fHitMap->FillHits();
 701     for (ndig=0; ndig<fNdigits; ndig++) {
 702         dig = (AliMUONdigit*)fDigits->UncheckedAt(ndig);
 703         Int_t i=dig->fPadX;
 704         Int_t j=dig->fPadY;
 705         if (fHitMap->TestHit(i,j)==used ||fHitMap->TestHit(i,j)==empty) {
 706             nskip++;
 707             continue;
 708         }
 709         AliMUONRawCluster c;
 710         c.fMultiplicity=0;
 711         c.fPeakSignal=dig->fSignal;
 712 /*
 713         c.fTracks[0]=dig->fTracks[0];
 714         c.fTracks[1]=dig->fTracks[1];
 715         c.fTracks[2]=dig->fTracks[2];
 716 */
 717         //c.fTracks[0]=dig->fTrack;
 718         c.fTracks[0]=dig->fHit;
 719         c.fTracks[1]=dig->fTracks[0];
 720         c.fTracks[2]=dig->fTracks[1];
 721         // tag the beginning of cluster list in a raw cluster
 722         c.fNcluster[0]=-1;
 723         FindCluster(i,j, c);
 724         // center of gravity
 725         c.fX /= c.fQ;
 726         c.fX=fSegmentation->GetAnod(c.fX);
 727         c.fY /= c.fQ;
 728 //
 729 //  apply correction to the coordinate along the anode wire
 730 //
 731         Int_t ix,iy;
 732         Float_t x=c.fX;
 733         Float_t y=c.fY;
 734         fSegmentation->GetPadIxy(x, y, ix, iy);
 735         fSegmentation->GetPadCxy(ix, iy, x, y);
 736         Int_t isec=fSegmentation->Sector(ix,iy);
 737         TF1* CogCorr=fSegmentation->CorrFunc(isec-1);
 738         if (CogCorr) {
 739             Float_t YonPad=(c.fY-y)/fSegmentation->Dpy(isec);
 740             c.fY=c.fY-CogCorr->Eval(YonPad,0,0);
 741         }
 742
 743 //
 744 //      Analyse cluster and decluster if necessary
 745 //
 746     ncls++;
 747     c.fNcluster[1]=fNRawClusters;
 748     c.fClusterType=c.PhysicsContribution();
 749     Decluster(&c);
 750     fNPeaks=0;
 751 //
 752 //
 753 //
 754 //      reset Cluster object
 755         for (int k=0;k<c.fMultiplicity;k++) {
 756             c.fIndexMap[k]=0;
 757         }
 758         c.fMultiplicity=0;
 759     } // end loop ndig
 760     delete fHitMap;
 761 }
 762
 763 void AliMUONClusterFinder::
 764 CalibrateCOG()
 765 {
 766     Float_t x[5];
 767     Float_t y[5];
 768     Int_t n, i;
 769     TF1 func;
 770     if (fSegmentation) {
 771         fSegmentation->GiveTestPoints(n, x, y);
 772         for (i=0; i<n; i++) {
 773             Float_t xtest=x[i];
 774             Float_t ytest=y[i];
 775             SinoidalFit(xtest, ytest, func);
 776             fSegmentation->SetCorrFunc(i, new TF1(func));
 777         }
 778     }
 779 }
 780
 781
 782 void AliMUONClusterFinder::
 783 SinoidalFit(Float_t x, Float_t y, TF1 &func)
 784 {
 785 //
 786     static Int_t count=0;
 787     char canvasname[3];
 788     count++;
 789     sprintf(canvasname,"c%d",count);
 790
 791     const Int_t ns=101;
 792     Float_t xg[ns], yg[ns], xrg[ns], yrg[ns];
 793     Float_t xsig[ns], ysig[ns];
 794
 795     AliMUONsegmentation *segmentation=fSegmentation;
 796
 797     Int_t ix,iy;
 798     segmentation->GetPadIxy(x,y,ix,iy);
 799     segmentation->GetPadCxy(ix,iy,x,y);
 800     Int_t isec=segmentation->Sector(ix,iy);
 801 // Pad Limits
 802     Float_t xmin = x-segmentation->Dpx(isec)/2;
 803     Float_t ymin = y-segmentation->Dpy(isec)/2;
 804 //
 805 //      Integration Limits
 806     Float_t dxI=fResponse->SigmaIntegration()*fResponse->ChargeSpreadX();
 807     Float_t dyI=fResponse->SigmaIntegration()*fResponse->ChargeSpreadY();
 808
 809 //
 810 //  Scanning
 811 //
 812     Int_t i;
 813     Float_t qp;
 814 //
 815 //  y-position
 816     Float_t yscan=ymin;
 817     Float_t dy=segmentation->Dpy(isec)/(ns-1);
 818
 819     for (i=0; i<ns; i++) {
 820 //
 821 //      Pad Loop
 822 //
 823         Float_t sum=0;
 824         Float_t qcheck=0;
 825         segmentation->SigGenInit(x, yscan, 0);
 826
 827         for (segmentation->FirstPad(x, yscan, dxI, dyI);
 828              segmentation->MorePads();
 829              segmentation->NextPad())
 830         {
 831             qp=fResponse->IntXY(segmentation);
 832             qp=TMath::Abs(qp);
 833 //
 834 //
 835             if (qp > 1.e-4) {
 836                 qcheck+=qp;
 837                 Int_t ixs=segmentation->Ix();
 838                 Int_t iys=segmentation->Iy();
 839                 Float_t xs,ys;
 840                 segmentation->GetPadCxy(ixs,iys,xs,ys);
 841                 sum+=qp*ys;
 842             }
 843         } // Pad loop
 844         Float_t ycog=sum/qcheck;
 845         yg[i]=(yscan-y)/segmentation->Dpy(isec);
 846         yrg[i]=(ycog-y)/segmentation->Dpy(isec);
 847         ysig[i]=ycog-yscan;
 848         yscan+=dy;
 849     } // scan loop
 850 //
 851 //  x-position
 852     Float_t xscan=xmin;
 853     Float_t dx=segmentation->Dpx(isec)/(ns-1);
 854
 855     for (i=0; i<ns; i++) {
 856 //
 857 //      Pad Loop
 858 //
 859         Float_t sum=0;
 860         Float_t qcheck=0;
 861         segmentation->SigGenInit(xscan, y, 0);
 862
 863         for (segmentation->FirstPad(xscan, y, dxI, dyI);
 864              segmentation->MorePads();
 865              segmentation->NextPad())
 866         {
 867             qp=fResponse->IntXY(segmentation);
 868             qp=TMath::Abs(qp);
 869 //
 870 //
 871             if (qp > 1.e-2) {
 872                 qcheck+=qp;
 873                 Int_t ixs=segmentation->Ix();
 874                 Int_t iys=segmentation->Iy();
 875                 Float_t xs,ys;
 876                 segmentation->GetPadCxy(ixs,iys,xs,ys);
 877                 sum+=qp*xs;
 878             }
 879         } // Pad loop
 880         Float_t xcog=sum/qcheck;
 881         xcog=segmentation->GetAnod(xcog);
 882
 883         xg[i]=(xscan-x)/segmentation->Dpx(isec);
 884         xrg[i]=(xcog-x)/segmentation->Dpx(isec);
 885         xsig[i]=xcog-xscan;
 886         xscan+=dx;
 887     }
 888 //
 889 // Creates a Root function based on function sinoid above
 890 // and perform the fit
 891 //
 892     //    TGraph *graphx = new TGraph(ns,xg ,xsig);
 893     //    TGraph *graphxr= new TGraph(ns,xrg,xsig);
 894     //    TGraph *graphy = new TGraph(ns,yg ,ysig);
 895     TGraph *graphyr= new TGraph(ns,yrg,ysig);
 896
 897     Double_t sinoid(Double_t *x, Double_t *par);
 898     new TF1("sinoidf",sinoid,0.5,0.5,5);
 899     graphyr->Fit("sinoidf","Q");
 900     func = *((TF1*)((graphyr->GetListOfFunctions())->At(0)));
 901 /*
 902
 903     TCanvas *c1=new TCanvas(canvasname,canvasname,400,10,600,700);
 904     TPad* pad11 = new TPad("pad11"," ",0.01,0.51,0.49,0.99);
 905     TPad* pad12 = new TPad("pad12"," ",0.51,0.51,0.99,0.99);
 906     TPad* pad13 = new TPad("pad13"," ",0.01,0.01,0.49,0.49);
 907     TPad* pad14 = new TPad("pad14"," ",0.51,0.01,0.99,0.49);
 908     pad11->SetFillColor(11);
 909     pad12->SetFillColor(11);
 910     pad13->SetFillColor(11);
 911     pad14->SetFillColor(11);
 912     pad11->Draw();
 913     pad12->Draw();
 914     pad13->Draw();
 915     pad14->Draw();
 916
 917 //
 918     pad11->cd();
 919     graphx->SetFillColor(42);
 920     graphx->SetMarkerColor(4);
 921     graphx->SetMarkerStyle(21);
 922     graphx->Draw("AC");
 923     graphx->GetHistogram()->SetXTitle("x on pad");
 924     graphx->GetHistogram()->SetYTitle("xcog-x");
 925
 926
 927     pad12->cd();
 928     graphxr->SetFillColor(42);
 929     graphxr->SetMarkerColor(4);
 930     graphxr->SetMarkerStyle(21);
 931     graphxr->Draw("AP");
 932     graphxr->GetHistogram()->SetXTitle("xcog on pad");
 933     graphxr->GetHistogram()->SetYTitle("xcog-x");
 934
 935
 936     pad13->cd();
 937     graphy->SetFillColor(42);
 938     graphy->SetMarkerColor(4);
 939     graphy->SetMarkerStyle(21);
 940     graphy->Draw("AF");
 941     graphy->GetHistogram()->SetXTitle("y on pad");
 942     graphy->GetHistogram()->SetYTitle("ycog-y");
 943
 944
 945
 946     pad14->cd();
 947     graphyr->SetFillColor(42);
 948     graphyr->SetMarkerColor(4);
 949     graphyr->SetMarkerStyle(21);
 950     graphyr->Draw("AF");
 951     graphyr->GetHistogram()->SetXTitle("ycog on pad");
 952     graphyr->GetHistogram()->SetYTitle("ycog-y");
 953
 954     c1->Update();
 955 */
 956 }
 957
 958 Double_t sinoid(Double_t *x, Double_t *par)
 959 {
 960     Double_t arg = -2*TMath::Pi()*x[0];
 961     Double_t fitval= par[0]*TMath::Sin(arg)+
 962         par[1]*TMath::Sin(2*arg)+
 963         par[2]*TMath::Sin(3*arg)+
 964         par[3]*TMath::Sin(4*arg)+
 965         par[4]*TMath::Sin(5*arg);
 966     return fitval;
 967  }
 968
 969
 970 Double_t DoubleGauss(Double_t *x, Double_t *par)
 971 {
 972     Double_t arg1 = (x[0]-par[1])/0.18;
 973     Double_t arg2 = (x[0]-par[3])/0.18;
 974     Double_t fitval= par[0]*TMath::Exp(-arg1*arg1/2)
 975         +par[2]*TMath::Exp(-arg2*arg2/2);
 976     return fitval;
 977  }
 978
 979 Float_t DiscrCharge(Int_t i,Double_t *par)
 980 {
 981 // par[0]    x-position of first  cluster
 982 // par[1]    y-position of first  cluster
 983 // par[2]    x-position of second cluster
 984 // par[3]    y-position of second cluster
 985 // par[4]    charge fraction of first  cluster
 986 // 1-par[4]  charge fraction of second cluster
 987
 988     static Float_t qtot;
 989     if (gFirst) {
 990         qtot=0;
 991         for (Int_t jbin=0; jbin<gNbins; jbin++) {
 992             qtot+=gCharge[jbin];
 993         }
 994         gFirst=0;
 995         //printf("\n sum of charge from DiscrCharge %f\n", qtot);
 996         gChargeTot=Int_t(qtot);
 997
 998     }
 999     gSegmentation->SetPad(gix[i], giy[i]);
1000 //  First Cluster
1001     gSegmentation->SetHit(par[0],par[1]);
1002     Float_t q1=gResponse->IntXY(gSegmentation);
1003
1004 //  Second Cluster
1005     gSegmentation->SetHit(par[2],par[3]);
1006     Float_t q2=gResponse->IntXY(gSegmentation);
1007
1008     Float_t value = qtot*(par[4]*q1+(1.-par[4])*q2);
1009     return value;
1010 }
1011
1012 //
1013 // Minimisation function
1014 void fcn(Int_t &, Double_t *, Double_t &f, Double_t *par, Int_t )
1015 {
1016     Int_t i;
1017     Float_t delta;
1018     Float_t chisq=0;
1019     Float_t qcont=0;
1020     Float_t qtot=0;
1021
1022     for (i=0; i<gNbins; i++) {
1023         Float_t q0=gCharge[i];
1024         Float_t q1=DiscrCharge(i,par);
1025         delta=(q0-q1)/TMath::Sqrt(q0);
1026         chisq+=delta*delta;
1027         qcont+=q1;
1028         qtot+=q0;
1029     }
1030     chisq=chisq+=(qtot-qcont)*(qtot-qcont)*0.5;
1031     f=chisq;
1032 }
1033
1034
1035
1036
1037
1038
1039
1040
1041
1042
1043