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