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