MUON/AliMUONClusterFinderVS.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 /* $Id$ */
  17
  18 #include <TMinuit.h>
  19 #include <TF1.h>
  20
  21 #include "AliMUONClusterFinderVS.h"
  22 #include "AliMUONDigit.h"
  23 #include "AliMUONRawCluster.h"
  24 #include "AliSegmentation.h"
  25 #include "AliMUONResponse.h"
  26 #include "AliMUONClusterInput.h"
  27 #include "AliMUONHitMapA1.h"
  28
  29 //_____________________________________________________________________
  30 // This function is minimized in the double-Mathieson fit
  31 void fcnS2(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
  32 void fcnS1(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
  33 void fcnCombiS1(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
  34 void fcnCombiS2(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
  35
  36 ClassImp(AliMUONClusterFinderVS)
  37
  38 AliMUONClusterFinderVS::AliMUONClusterFinderVS()
  39   : TObject()
  40 {
  41 // Default constructor
  42     fInput=AliMUONClusterInput::Instance();
  43     fHitMap[0] = 0;
  44     fHitMap[1] = 0;
  45     fTrack[0]=fTrack[1]=-1;
  46     fDebugLevel = 0; // make silent default
  47     fGhostChi2Cut = 1e6; // nothing done by default
  48     fSeg[0]    = 0;
  49     fSeg[1]    = 0;
  50     for(Int_t i=0; i<100; i++) {
  51       for (Int_t j=0; j<2; j++) {
  52         fDig[i][j] = 0;
  53       }
  54     }
  55     fRawClusters = new TClonesArray("AliMUONRawCluster",1000);
  56     fNRawClusters = 0;
  57 }
  58  //____________________________________________________________________________
  59 AliMUONClusterFinderVS::~AliMUONClusterFinderVS()
  60 {
  61   // Reset tracks information
  62    fNRawClusters = 0;
  63    if (fRawClusters) {
  64      fRawClusters->Delete();
  65      delete fRawClusters;
  66    }
  67 }
  68
  69 AliMUONClusterFinderVS::AliMUONClusterFinderVS(const AliMUONClusterFinderVS & clusterFinder):TObject(clusterFinder)
  70 {
  71 // Protected copy constructor
  72
  73   Fatal("AliMUONClusterFinderAZModule", "Not implemented.");
  74 }
  75 //____________________________________________________________________________
  76 void AliMUONClusterFinderVS::ResetRawClusters()
  77 {
  78   // Reset tracks information
  79   fNRawClusters = 0;
  80   if (fRawClusters) fRawClusters->Clear();
  81 }
  82 //____________________________________________________________________________
  83 void AliMUONClusterFinderVS::Decluster(AliMUONRawCluster *cluster)
  84 {
  85 // Decluster by local maxima
  86     SplitByLocalMaxima(cluster);
  87 }
  88 //____________________________________________________________________________
  89 void AliMUONClusterFinderVS::SplitByLocalMaxima(AliMUONRawCluster *c)
  90 {
  91 // Split complex cluster by local maxima
  92     Int_t cath, i;
  93
  94     fInput->SetCluster(c);
  95
  96     fMul[0]=c->GetMultiplicity(0);
  97     fMul[1]=c->GetMultiplicity(1);
  98
  99 //
 100 //  dump digit information into arrays
 101 //
 102
 103     Float_t qtot;
 104
 105     for (cath=0; cath<2; cath++) {
 106         qtot=0;
 107         for (i=0; i<fMul[cath]; i++)
 108         {
 109             // pointer to digit
 110             fDig[i][cath]=fInput->Digit(cath, c->GetIndex(i, cath));
 111             // pad coordinates
 112             fIx[i][cath]= fDig[i][cath]->PadX();
 113             fIy[i][cath]= fDig[i][cath]->PadY();
 114             // pad charge
 115             fQ[i][cath] = fDig[i][cath]->Signal();
 116             // pad centre coordinates
 117             fSeg[cath]->
 118                 GetPadC(fIx[i][cath], fIy[i][cath], fX[i][cath], fY[i][cath], fZ[i][cath]);
 119         } // loop over cluster digits
 120     }  // loop over cathodes
 121
 122
 123     FindLocalMaxima(c);
 124
 125 //
 126 //  Initialise and perform mathieson fits
 127     Float_t chi2, oldchi2;
 128 //  ++++++++++++++++++*************+++++++++++++++++++++
 129 //  (1) No more than one local maximum per cathode plane
 130 //  +++++++++++++++++++++++++++++++*************++++++++
 131     if ((fNLocal[0]==1 && (fNLocal[1]==0 ||  fNLocal[1]==1)) ||
 132         (fNLocal[0]==0 && fNLocal[1]==1)) {
 133 // Perform combined single Mathieson fit
 134 // Initial values for coordinates (x,y)
 135
 136         // One local maximum on cathodes 1 and 2 (X->cathode 2, Y->cathode 1)
 137         if (fNLocal[0]==1 &&  fNLocal[1]==1) {
 138             fXInit[0]=c->GetX(1);
 139             fYInit[0]=c->GetY(0);
 140             // One local maximum on cathode 1 (X,Y->cathode 1)
 141         } else if (fNLocal[0]==1) {
 142             fXInit[0]=c->GetX(0);
 143             fYInit[0]=c->GetY(0);
 144             // One local maximum on cathode 2  (X,Y->cathode 2)
 145         } else {
 146             fXInit[0]=c->GetX(1);
 147             fYInit[0]=c->GetY(1);
 148         }
 149         if (fDebugLevel)
 150             fprintf(stderr,"\n cas (1) CombiSingleMathiesonFit(c)\n");
 151         chi2=CombiSingleMathiesonFit(c);
 152 //      Int_t ndf = fgNbins[0]+fgNbins[1]-2;
 153 //      Float_t prob = TMath::Prob(Double_t(chi2),ndf);
 154 //      prob1->Fill(prob);
 155 //      chi2_1->Fill(chi2);
 156         oldchi2=chi2;
 157         if (fDebugLevel)
 158             fprintf(stderr," chi2 %f ",chi2);
 159
 160         c->SetX(0, fXFit[0]);
 161         c->SetY(0, fYFit[0]);
 162
 163         c->SetX(1,fXFit[0]);
 164         c->SetY(1,fYFit[0]);
 165         c->SetChi2(0,chi2);
 166         c->SetChi2(1,chi2);
 167         // Force on anod
 168         c->SetX(0, fSeg[0]->GetAnod(c->GetX(0)));
 169         c->SetX(1, fSeg[1]->GetAnod(c->GetX(1)));
 170
 171 // If reasonable chi^2 add result to the list of rawclusters
 172         if (chi2 < 0.3) {
 173             AddRawCluster(*c);
 174 // If not try combined double Mathieson Fit
 175         } else {
 176           if (fDebugLevel)
 177             fprintf(stderr," MAUVAIS CHI2 !!!\n");
 178             if (fNLocal[0]==1 &&  fNLocal[1]==1) {
 179                 fXInit[0]=fX[fIndLocal[0][1]][1];
 180                 fYInit[0]=fY[fIndLocal[0][0]][0];
 181                 fXInit[1]=fX[fIndLocal[0][1]][1];
 182                 fYInit[1]=fY[fIndLocal[0][0]][0];
 183             } else if (fNLocal[0]==1) {
 184                 fXInit[0]=fX[fIndLocal[0][0]][0];
 185                 fYInit[0]=fY[fIndLocal[0][0]][0];
 186                 fXInit[1]=fX[fIndLocal[0][0]][0];
 187                 fYInit[1]=fY[fIndLocal[0][0]][0];
 188             } else {
 189                 fXInit[0]=fX[fIndLocal[0][1]][1];
 190                 fYInit[0]=fY[fIndLocal[0][1]][1];
 191                 fXInit[1]=fX[fIndLocal[0][1]][1];
 192                 fYInit[1]=fY[fIndLocal[0][1]][1];
 193             }
 194
 195 //  Initial value for charge ratios
 196             fQrInit[0]=0.5;
 197             fQrInit[1]=0.5;
 198             if (fDebugLevel)
 199             fprintf(stderr,"\n cas (1) CombiDoubleMathiesonFit(c)\n");
 200             chi2=CombiDoubleMathiesonFit(c);
 201 //          Int_t ndf = fgNbins[0]+fgNbins[1]-6;
 202 //          Float_t prob = TMath::Prob(chi2,ndf);
 203 //          prob2->Fill(prob);
 204 //          chi2_2->Fill(chi2);
 205
 206 // Was this any better ??
 207             if (fDebugLevel)
 208               fprintf(stderr," Old and new chi2 %f %f ", oldchi2, chi2);
 209             if (fFitStat!=0 && chi2>0 && (2.*chi2 < oldchi2)) {
 210               if (fDebugLevel)
 211                 fprintf(stderr," Split\n");
 212                 // Split cluster into two according to fit result
 213                 Split(c);
 214             } else {
 215               if (fDebugLevel)
 216                 fprintf(stderr," Don't Split\n");
 217                 // Don't split
 218                 AddRawCluster(*c);
 219             }
 220         }
 221
 222 //  +++++++++++++++++++++++++++++++++++++++
 223 //  (2) Two local maxima per cathode plane
 224 //  +++++++++++++++++++++++++++++++++++++++
 225     } else if (fNLocal[0]==2 &&  fNLocal[1]==2) {
 226 //
 227 //  Let's look for ghosts first
 228
 229         Float_t xm[4][2], ym[4][2];
 230         Float_t dpx, dpy, dx, dy;
 231         Int_t ixm[4][2], iym[4][2];
 232         Int_t isec, im1, im2, ico;
 233 //
 234 //  Form the 2x2 combinations
 235 //  0-0, 0-1, 1-0, 1-1
 236         ico=0;
 237         for (im1=0; im1<2; im1++) {
 238             for (im2=0; im2<2; im2++) {
 239                 xm[ico][0]=fX[fIndLocal[im1][0]][0];
 240                 ym[ico][0]=fY[fIndLocal[im1][0]][0];
 241                 xm[ico][1]=fX[fIndLocal[im2][1]][1];
 242                 ym[ico][1]=fY[fIndLocal[im2][1]][1];
 243
 244                 ixm[ico][0]=fIx[fIndLocal[im1][0]][0];
 245                 iym[ico][0]=fIy[fIndLocal[im1][0]][0];
 246                 ixm[ico][1]=fIx[fIndLocal[im2][1]][1];
 247                 iym[ico][1]=fIy[fIndLocal[im2][1]][1];
 248                 ico++;
 249             }
 250         }
 251 // ico = 0 : first local maximum on cathodes 1 and 2
 252 // ico = 1 : fisrt local maximum on cathode 1 and second on cathode 2
 253 // ico = 2 : second local maximum on cathode 1 and first on cathode 1
 254 // ico = 3 : second local maximum on cathodes 1 and 2
 255
 256 // Analyse the combinations and keep those that are possible !
 257 // For each combination check consistency in x and y
 258         Int_t   iacc;
 259         Bool_t  accepted[4];
 260         Float_t dr[4] = {1.e4, 1.e4, 1.e4, 1.e4};
 261         iacc=0;
 262
 263 // In case of staggering maxima are displaced by exactly half the pad-size in y.
 264 // We have to take into account the numerical precision in the consistency check;
 265         Float_t eps = 1.e-5;
 266 //
 267         for (ico=0; ico<4; ico++) {
 268             accepted[ico]=kFALSE;
 269 // cathode one: x-coordinate
 270             isec=fSeg[0]->Sector(ixm[ico][0], iym[ico][0]);
 271             dpx=fSeg[0]->Dpx(isec)/2.;
 272             dx=TMath::Abs(xm[ico][0]-xm[ico][1]);
 273 // cathode two: y-coordinate
 274             isec=fSeg[1]->Sector(ixm[ico][1], iym[ico][1]);
 275             dpy=fSeg[1]->Dpy(isec)/2.;
 276             dy=TMath::Abs(ym[ico][0]-ym[ico][1]);
 277             if (fDebugLevel>1)
 278                 printf("\n %i %f %f %f %f %f %f \n", ico, ym[ico][0], ym[ico][1], dy, dpy, dx, dpx );
 279             if ((dx <= dpx) && (dy <= dpy+eps)) {
 280                 // consistent
 281                 accepted[ico]=kTRUE;
 282                 dr[ico] = TMath::Sqrt(dx*dx+dy*dy);
 283                 iacc++;
 284             } else {
 285                 // reject
 286                 accepted[ico]=kFALSE;
 287             }
 288         }
 289         if (fDebugLevel)
 290           printf("\n iacc= %d:\n", iacc);
 291         if (iacc == 3) {
 292             if (accepted[0] && accepted[1]) {
 293                 if (dr[0] >= dr[1]) {
 294                     accepted[0]=kFALSE;
 295                 } else {
 296                     accepted[1]=kFALSE;
 297                 }
 298             }
 299
 300             if (accepted[2] && accepted[3]) {
 301                 if (dr[2] >= dr[3]) {
 302                     accepted[2]=kFALSE;
 303                 } else {
 304                     accepted[3]=kFALSE;
 305                 }
 306             }
 307 /*
 308 // eliminate one candidate
 309             Float_t drmax = 0;
 310             Int_t icobad = -1;
 311
 312             for (ico=0; ico<4; ico++) {
 313                 if (accepted[ico] && dr[ico] > drmax) {
 314                     icobad = ico;
 315                     drmax  = dr[ico];
 316                 }
 317             }
 318
 319             accepted[icobad] = kFALSE;
 320 */
 321             iacc = 2;
 322         }
 323
 324
 325         if (fDebugLevel) {
 326           printf("\n iacc= %d:\n", iacc);
 327             if (iacc==2) {
 328                 fprintf(stderr,"\n iacc=2: No problem ! \n");
 329             } else if (iacc==4) {
 330                 fprintf(stderr,"\n iacc=4: Ok, but ghost problem !!! \n");
 331             } else if (iacc==0) {
 332                 fprintf(stderr,"\n iacc=0: I don't know what to do with this !!!!!!!!! \n");
 333             }
 334         }
 335
 336 //  Initial value for charge ratios
 337         fQrInit[0]=Float_t(fQ[fIndLocal[0][0]][0])/
 338             Float_t(fQ[fIndLocal[0][0]][0]+fQ[fIndLocal[1][0]][0]);
 339         fQrInit[1]=Float_t(fQ[fIndLocal[0][1]][1])/
 340             Float_t(fQ[fIndLocal[0][1]][1]+fQ[fIndLocal[1][1]][1]);
 341
 342 // ******* iacc = 0 *******
 343 // No combinations found between the 2 cathodes
 344 // We keep the center of gravity of the cluster
 345         if (iacc==0) {
 346             AddRawCluster(*c);
 347         }
 348
 349 // ******* iacc = 1 *******
 350 // Only one combination found between the 2 cathodes
 351         if (iacc==1) {
 352 // Initial values for the 2 maxima (x,y)
 353
 354 // 1 maximum is initialised with the maximum of the combination found (X->cathode 2, Y->cathode 1)
 355 // 1 maximum is initialised with the other maximum of the first cathode
 356             if (accepted[0]){
 357                 fprintf(stderr,"ico=0\n");
 358                 fXInit[0]=xm[0][1];
 359                 fYInit[0]=ym[0][0];
 360                 fXInit[1]=xm[3][0];
 361                 fYInit[1]=ym[3][0];
 362             } else if (accepted[1]){
 363                 fprintf(stderr,"ico=1\n");
 364                 fXInit[0]=xm[1][1];
 365                 fYInit[0]=ym[1][0];
 366                 fXInit[1]=xm[2][0];
 367                 fYInit[1]=ym[2][0];
 368             } else if (accepted[2]){
 369                 fprintf(stderr,"ico=2\n");
 370                 fXInit[0]=xm[2][1];
 371                 fYInit[0]=ym[2][0];
 372                 fXInit[1]=xm[1][0];
 373                 fYInit[1]=ym[1][0];
 374             } else if (accepted[3]){
 375                 fprintf(stderr,"ico=3\n");
 376                 fXInit[0]=xm[3][1];
 377                 fYInit[0]=ym[3][0];
 378                 fXInit[1]=xm[0][0];
 379                 fYInit[1]=ym[0][0];
 380             }
 381             if (fDebugLevel)
 382                 fprintf(stderr,"\n cas (2) CombiDoubleMathiesonFit(c)\n");
 383             chi2=CombiDoubleMathiesonFit(c);
 384 //          Int_t ndf = fgNbins[0]+fgNbins[1]-6;
 385 //          Float_t prob = TMath::Prob(chi2,ndf);
 386 //          prob2->Fill(prob);
 387 //          chi2_2->Fill(chi2);
 388             if (fDebugLevel)
 389                 fprintf(stderr," chi2 %f\n",chi2);
 390
 391 // If reasonable chi^2 add result to the list of rawclusters
 392             if (chi2<10) {
 393                 Split(c);
 394
 395             } else {
 396 // 1 maximum is initialised with the maximum of the combination found (X->cathode 2, Y->cathode 1)
 397 // 1 maximum is initialised with the other maximum of the second cathode
 398                 if (accepted[0]){
 399                     fprintf(stderr,"ico=0\n");
 400                     fXInit[0]=xm[0][1];
 401                     fYInit[0]=ym[0][0];
 402                     fXInit[1]=xm[3][1];
 403                     fYInit[1]=ym[3][1];
 404                 } else if (accepted[1]){
 405                     fprintf(stderr,"ico=1\n");
 406                     fXInit[0]=xm[1][1];
 407                     fYInit[0]=ym[1][0];
 408                     fXInit[1]=xm[2][1];
 409                     fYInit[1]=ym[2][1];
 410                 } else if (accepted[2]){
 411                     fprintf(stderr,"ico=2\n");
 412                     fXInit[0]=xm[2][1];
 413                     fYInit[0]=ym[2][0];
 414                     fXInit[1]=xm[1][1];
 415                     fYInit[1]=ym[1][1];
 416                 } else if (accepted[3]){
 417                     fprintf(stderr,"ico=3\n");
 418                     fXInit[0]=xm[3][1];
 419                     fYInit[0]=ym[3][0];
 420                     fXInit[1]=xm[0][1];
 421                     fYInit[1]=ym[0][1];
 422                 }
 423                 if (fDebugLevel)
 424                     fprintf(stderr,"\n cas (2) CombiDoubleMathiesonFit(c)\n");
 425                 chi2=CombiDoubleMathiesonFit(c);
 426 //              Int_t ndf = fgNbins[0]+fgNbins[1]-6;
 427 //              Float_t prob = TMath::Prob(chi2,ndf);
 428 //              prob2->Fill(prob);
 429 //              chi2_2->Fill(chi2);
 430                 if (fDebugLevel)
 431                     fprintf(stderr," chi2 %f\n",chi2);
 432
 433 // If reasonable chi^2 add result to the list of rawclusters
 434                 if (chi2<10) {
 435                     Split(c);
 436                 } else {
 437 //We keep only the combination found (X->cathode 2, Y->cathode 1)
 438                     for (Int_t ico=0; ico<2; ico++) {
 439                         if (accepted[ico]) {
 440                             AliMUONRawCluster cnew;
 441                             Int_t cath;
 442                             for (cath=0; cath<2; cath++) {
 443                                 cnew.SetX(cath, Float_t(xm[ico][1]));
 444                                 cnew.SetY(cath, Float_t(ym[ico][0]));
 445                                 cnew.SetZ(cath, fZPlane);
 446
 447                                 cnew.SetMultiplicity(cath,c->GetMultiplicity(cath));
 448                                 for (i=0; i<fMul[cath]; i++) {
 449                                     cnew.SetIndex(i, cath, c->GetIndex(i,cath));
 450                                     fSeg[cath]->SetPad(fIx[i][cath], fIy[i][cath]);
 451                                 }
 452                                 fprintf(stderr,"\nRawCluster %d cath %d\n",ico,cath);
 453                                 fprintf(stderr,"mult_av %d\n",c->GetMultiplicity(cath));
 454                                 FillCluster(&cnew,cath);
 455                             }
 456                             cnew.SetClusterType(cnew.PhysicsContribution());
 457                             AddRawCluster(cnew);
 458                             fNPeaks++;
 459                         }
 460                     }
 461                 }
 462             }
 463         }
 464
 465 // ******* iacc = 2 *******
 466 // Two combinations found between the 2 cathodes
 467         if (iacc==2) {
 468 // Was the same maximum taken twice
 469             if ((accepted[0]&&accepted[1]) || (accepted[2]&&accepted[3])) {
 470                 fprintf(stderr,"\n Maximum taken twice !!!\n");
 471
 472 // Have a try !! with that
 473                 if (accepted[0]&&accepted[3]) {
 474                     fXInit[0]=xm[0][1];
 475                     fYInit[0]=ym[0][0];
 476                     fXInit[1]=xm[1][1];
 477                     fYInit[1]=ym[1][0];
 478                 } else {
 479                     fXInit[0]=xm[2][1];
 480                     fYInit[0]=ym[2][0];
 481                     fXInit[1]=xm[3][1];
 482                     fYInit[1]=ym[3][0];
 483                 }
 484                 if (fDebugLevel)
 485                     fprintf(stderr,"\n cas (2) CombiDoubleMathiesonFit(c)\n");
 486                 chi2=CombiDoubleMathiesonFit(c);
 487 //                  Int_t ndf = fgNbins[0]+fgNbins[1]-6;
 488 //                  Float_t prob = TMath::Prob(chi2,ndf);
 489 //                  prob2->Fill(prob);
 490 //                  chi2_2->Fill(chi2);
 491                 Split(c);
 492
 493             } else {
 494 // No ghosts ! No Problems ! -  Perform one fit only !
 495                 if (accepted[0]&&accepted[3]) {
 496                     fXInit[0]=xm[0][1];
 497                     fYInit[0]=ym[0][0];
 498                     fXInit[1]=xm[3][1];
 499                     fYInit[1]=ym[3][0];
 500                 } else {
 501                     fXInit[0]=xm[1][1];
 502                     fYInit[0]=ym[1][0];
 503                     fXInit[1]=xm[2][1];
 504                     fYInit[1]=ym[2][0];
 505                 }
 506                 if (fDebugLevel)
 507                     fprintf(stderr,"\n cas (2) CombiDoubleMathiesonFit(c)\n");
 508                 chi2=CombiDoubleMathiesonFit(c);
 509 //                  Int_t ndf = fgNbins[0]+fgNbins[1]-6;
 510 //                  Float_t prob = TMath::Prob(chi2,ndf);
 511 //                  prob2->Fill(prob);
 512 //                  chi2_2->Fill(chi2);
 513                 if (fDebugLevel)
 514                     fprintf(stderr," chi2 %f\n",chi2);
 515                 Split(c);
 516             }
 517
 518 // ******* iacc = 4 *******
 519 // Four combinations found between the 2 cathodes
 520 // Ghost !!
 521         } else if (iacc==4) {
 522 // Perform fits for the two possibilities !!
 523 // Accept if charges are compatible on both cathodes
 524 // If none are compatible, keep everything
 525             fXInit[0]=xm[0][1];
 526             fYInit[0]=ym[0][0];
 527             fXInit[1]=xm[3][1];
 528             fYInit[1]=ym[3][0];
 529             if (fDebugLevel)
 530                 fprintf(stderr,"\n cas (2) CombiDoubleMathiesonFit(c)\n");
 531             chi2=CombiDoubleMathiesonFit(c);
 532 //              Int_t ndf = fgNbins[0]+fgNbins[1]-6;
 533 //              Float_t prob = TMath::Prob(chi2,ndf);
 534 //              prob2->Fill(prob);
 535 //              chi2_2->Fill(chi2);
 536             if (fDebugLevel)
 537                 fprintf(stderr," chi2 %f\n",chi2);
 538             // store results of fit and postpone decision
 539             Double_t sXFit[2],sYFit[2],sQrFit[2];
 540             Float_t sChi2[2];
 541             for (Int_t i=0;i<2;i++) {
 542                 sXFit[i]=fXFit[i];
 543                 sYFit[i]=fYFit[i];
 544                 sQrFit[i]=fQrFit[i];
 545                 sChi2[i]=fChi2[i];
 546             }
 547             fXInit[0]=xm[1][1];
 548             fYInit[0]=ym[1][0];
 549             fXInit[1]=xm[2][1];
 550             fYInit[1]=ym[2][0];
 551             if (fDebugLevel)
 552                 fprintf(stderr,"\n cas (2) CombiDoubleMathiesonFit(c)\n");
 553             chi2=CombiDoubleMathiesonFit(c);
 554 //              ndf = fgNbins[0]+fgNbins[1]-6;
 555 //              prob = TMath::Prob(chi2,ndf);
 556 //              prob2->Fill(prob);
 557 //              chi2_2->Fill(chi2);
 558             if (fDebugLevel)
 559                 fprintf(stderr," chi2 %f\n",chi2);
 560             // We have all informations to perform the decision
 561             // Compute the chi2 for the 2 possibilities
 562             Float_t chi2fi,chi2si,chi2f,chi2s;
 563
 564             chi2f = (TMath::Log(fInput->TotalCharge(0)*fQrFit[0]
 565                   /  (fInput->TotalCharge(1)*fQrFit[1]) )
 566                   / fInput->Response()->ChargeCorrel() );
 567             chi2f *=chi2f;
 568             chi2fi = (TMath::Log(fInput->TotalCharge(0)*(1-fQrFit[0])
 569                   /  (fInput->TotalCharge(1)*(1-fQrFit[1])) )
 570                   / fInput->Response()->ChargeCorrel() );
 571             chi2f += chi2fi*chi2fi;
 572
 573             chi2s = (TMath::Log(fInput->TotalCharge(0)*sQrFit[0]
 574                   /  (fInput->TotalCharge(1)*sQrFit[1]) )
 575                   / fInput->Response()->ChargeCorrel() );
 576             chi2s *=chi2s;
 577             chi2si = (TMath::Log(fInput->TotalCharge(0)*(1-sQrFit[0])
 578                   /  (fInput->TotalCharge(1)*(1-sQrFit[1])) )
 579                   / fInput->Response()->ChargeCorrel() );
 580             chi2s += chi2si*chi2si;
 581
 582             // usefull to store the charge matching chi2 in the cluster
 583             // fChi2[0]=sChi2[1]=chi2f;
 584             // fChi2[1]=sChi2[0]=chi2s;
 585
 586             if (chi2f<=fGhostChi2Cut && chi2s<=fGhostChi2Cut)
 587                 c->SetGhost(1);
 588             if   (chi2f>fGhostChi2Cut && chi2s>fGhostChi2Cut) {
 589                 // we keep the ghost
 590                 c->SetGhost(2);
 591                 chi2s=-1;
 592                 chi2f=-1;
 593             }
 594             if (chi2f<=fGhostChi2Cut)
 595                 Split(c);
 596             if (chi2s<=fGhostChi2Cut) {
 597                 // retreive saved values
 598                 for (Int_t i=0;i<2;i++) {
 599                     fXFit[i]=sXFit[i];
 600                     fYFit[i]=sYFit[i];
 601                     fQrFit[i]=sQrFit[i];
 602                     fChi2[i]=sChi2[i];
 603                 }
 604                 Split(c);
 605             }
 606             c->SetGhost(0);
 607         }
 608
 609     } else if (fNLocal[0]==2 &&  fNLocal[1]==1) {
 610 //  ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 611 //  (3) Two local maxima on cathode 1 and one maximum on cathode 2
 612 //  ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 613 //
 614         Float_t xm[4][2], ym[4][2];
 615         Float_t dpx, dpy, dx, dy;
 616         Int_t ixm[4][2], iym[4][2];
 617         Int_t isec, im1, ico;
 618 //
 619 //  Form the 2x2 combinations
 620 //  0-0, 0-1, 1-0, 1-1
 621         ico=0;
 622         for (im1=0; im1<2; im1++) {
 623             xm[ico][0]=fX[fIndLocal[im1][0]][0];
 624             ym[ico][0]=fY[fIndLocal[im1][0]][0];
 625             xm[ico][1]=fX[fIndLocal[0][1]][1];
 626             ym[ico][1]=fY[fIndLocal[0][1]][1];
 627
 628             ixm[ico][0]=fIx[fIndLocal[im1][0]][0];
 629             iym[ico][0]=fIy[fIndLocal[im1][0]][0];
 630             ixm[ico][1]=fIx[fIndLocal[0][1]][1];
 631             iym[ico][1]=fIy[fIndLocal[0][1]][1];
 632             ico++;
 633         }
 634 // ico = 0 : first local maximum on cathodes 1 and 2
 635 // ico = 1 : second local maximum on cathode 1 and first on cathode 2
 636
 637 // Analyse the combinations and keep those that are possible !
 638 // For each combination check consistency in x and y
 639         Int_t iacc;
 640         Bool_t accepted[4];
 641         iacc=0;
 642         // In case of staggering maxima are displaced by exactly half the pad-size in y.
 643         // We have to take into account the numerical precision in the consistency check;
 644
 645         Float_t eps = 1.e-5;
 646
 647         for (ico=0; ico<2; ico++) {
 648             accepted[ico]=kFALSE;
 649             isec=fSeg[0]->Sector(ixm[ico][0], iym[ico][0]);
 650             dpx=fSeg[0]->Dpx(isec)/2.;
 651             dx=TMath::Abs(xm[ico][0]-xm[ico][1]);
 652             isec=fSeg[1]->Sector(ixm[ico][1], iym[ico][1]);
 653             dpy=fSeg[1]->Dpy(isec)/2.;
 654             dy=TMath::Abs(ym[ico][0]-ym[ico][1]);
 655             if (fDebugLevel>1)
 656                 printf("\n %i %f %f %f %f \n", ico, ym[ico][0], ym[ico][1], dy, dpy );
 657             if ((dx <= dpx) && (dy <= dpy+eps)) {
 658                 // consistent
 659                 accepted[ico]=kTRUE;
 660                 iacc++;
 661             } else {
 662                 // reject
 663                 accepted[ico]=kFALSE;
 664             }
 665         }
 666
 667         Float_t chi21 = 100;
 668         Float_t chi22 = 100;
 669         Float_t chi23 = 100;
 670
 671         //  Initial value for charge ratios
 672         fQrInit[0]=Float_t(fQ[fIndLocal[0][0]][0])/
 673             Float_t(fQ[fIndLocal[0][0]][0]+fQ[fIndLocal[1][0]][0]);
 674         fQrInit[1]=fQrInit[0];
 675
 676         if (accepted[0] && accepted[1]) {
 677
 678             fXInit[0]=0.5*(xm[0][1]+xm[0][0]);
 679             fYInit[0]=ym[0][0];
 680             fXInit[1]=0.5*(xm[0][1]+xm[1][0]);
 681             fYInit[1]=ym[1][0];
 682             fQrInit[0]=0.5;
 683             fQrInit[1]=0.5;
 684             chi23=CombiDoubleMathiesonFit(c);
 685             if (chi23<10) {
 686                 Split(c);
 687                 Float_t yst;
 688                 yst = fYFit[0];
 689                 fYFit[0] = fYFit[1];
 690                 fYFit[1] = yst;
 691                 Split(c);
 692             }
 693         } else if (accepted[0]) {
 694             fXInit[0]=xm[0][1];
 695             fYInit[0]=ym[0][0];
 696             fXInit[1]=xm[1][0];
 697             fYInit[1]=ym[1][0];
 698             chi21=CombiDoubleMathiesonFit(c);
 699 //          Int_t ndf = fgNbins[0]+fgNbins[1]-6;
 700 //          Float_t prob = TMath::Prob(chi2,ndf);
 701 //          prob2->Fill(prob);
 702 //          chi2_2->Fill(chi21);
 703             if (fDebugLevel)
 704                 fprintf(stderr," chi2 %f\n",chi21);
 705             if (chi21<10) Split(c);
 706         } else if (accepted[1]) {
 707             fXInit[0]=xm[1][1];
 708             fYInit[0]=ym[1][0];
 709             fXInit[1]=xm[0][0];
 710             fYInit[1]=ym[0][0];
 711             chi22=CombiDoubleMathiesonFit(c);
 712 //          Int_t ndf = fgNbins[0]+fgNbins[1]-6;
 713 //          Float_t prob = TMath::Prob(chi2,ndf);
 714 //          prob2->Fill(prob);
 715 //          chi2_2->Fill(chi22);
 716             if (fDebugLevel)
 717                 fprintf(stderr," chi2 %f\n",chi22);
 718             if (chi22<10) Split(c);
 719         }
 720
 721         if (chi21 > 10 && chi22 > 10 && chi23 > 10) {
 722 // We keep only the combination found (X->cathode 2, Y->cathode 1)
 723             for (Int_t ico=0; ico<2; ico++) {
 724                 if (accepted[ico]) {
 725                     AliMUONRawCluster cnew;
 726                     Int_t cath;
 727                     for (cath=0; cath<2; cath++) {
 728                         cnew.SetX(cath, Float_t(xm[ico][1]));
 729                         cnew.SetY(cath, Float_t(ym[ico][0]));
 730                         cnew.SetZ(cath, fZPlane);
 731                         cnew.SetMultiplicity(cath, c->GetMultiplicity(cath));
 732                         for (i=0; i<fMul[cath]; i++) {
 733                             cnew.SetIndex(i, cath, c->GetIndex(i, cath));
 734                             fSeg[cath]->SetPad(fIx[i][cath], fIy[i][cath]);
 735                         }
 736                         fprintf(stderr,"\nRawCluster %d cath %d\n",ico,cath);
 737                         fprintf(stderr,"mult_av %d\n",c->GetMultiplicity(cath));
 738                         FillCluster(&cnew,cath);
 739                     }
 740                     cnew.SetClusterType(cnew.PhysicsContribution());
 741                     AddRawCluster(cnew);
 742                     fNPeaks++;
 743                 }
 744             }
 745         }
 746
 747 //  +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 748 //  (3') One local maximum on cathode 1 and two maxima on cathode 2
 749 //  +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 750     } else if (fNLocal[0]==1 && fNLocal[1]==2) {
 751         Float_t xm[4][2], ym[4][2];
 752         Float_t dpx, dpy, dx, dy;
 753         Int_t ixm[4][2], iym[4][2];
 754         Int_t isec, im1, ico;
 755 //
 756 //  Form the 2x2 combinations
 757 //  0-0, 0-1, 1-0, 1-1
 758         ico=0;
 759         for (im1=0; im1<2; im1++) {
 760             xm[ico][0]=fX[fIndLocal[0][0]][0];
 761             ym[ico][0]=fY[fIndLocal[0][0]][0];
 762             xm[ico][1]=fX[fIndLocal[im1][1]][1];
 763             ym[ico][1]=fY[fIndLocal[im1][1]][1];
 764
 765             ixm[ico][0]=fIx[fIndLocal[0][0]][0];
 766             iym[ico][0]=fIy[fIndLocal[0][0]][0];
 767             ixm[ico][1]=fIx[fIndLocal[im1][1]][1];
 768             iym[ico][1]=fIy[fIndLocal[im1][1]][1];
 769             ico++;
 770         }
 771 // ico = 0 : first local maximum on cathodes 1 and 2
 772 // ico = 1 : first local maximum on cathode 1 and second on cathode 2
 773
 774 // Analyse the combinations and keep those that are possible !
 775 // For each combination check consistency in x and y
 776         Int_t iacc;
 777         Bool_t accepted[4];
 778         iacc=0;
 779         // In case of staggering maxima are displaced by exactly half the pad-size in y.
 780         // We have to take into account the numerical precision in the consistency check;
 781         Float_t eps = 1.e-5;
 782
 783
 784         for (ico=0; ico<2; ico++) {
 785             accepted[ico]=kFALSE;
 786             isec=fSeg[0]->Sector(ixm[ico][0], iym[ico][0]);
 787             dpx=fSeg[0]->Dpx(isec)/2.;
 788             dx=TMath::Abs(xm[ico][0]-xm[ico][1]);
 789             isec=fSeg[1]->Sector(ixm[ico][1], iym[ico][1]);
 790             dpy=fSeg[1]->Dpy(isec)/2.;
 791             dy=TMath::Abs(ym[ico][0]-ym[ico][1]);
 792             if (fDebugLevel>0)
 793                 printf("\n %i %f %f %f %f \n", ico, ym[ico][0], ym[ico][1], dy, dpy );
 794             if ((dx <= dpx) && (dy <= dpy+eps)) {
 795                 // consistent
 796                 accepted[ico]=kTRUE;
 797                 fprintf(stderr,"ico %d\n",ico);
 798                 iacc++;
 799             } else {
 800                 // reject
 801                 accepted[ico]=kFALSE;
 802             }
 803         }
 804
 805         Float_t chi21 = 100;
 806         Float_t chi22 = 100;
 807         Float_t chi23 = 100;
 808
 809         fQrInit[1]=Float_t(fQ[fIndLocal[0][1]][1])/
 810             Float_t(fQ[fIndLocal[0][1]][1]+fQ[fIndLocal[1][1]][1]);
 811
 812         fQrInit[0]=fQrInit[1];
 813
 814
 815         if (accepted[0] && accepted[1]) {
 816             fXInit[0]=xm[0][1];
 817             fYInit[0]=0.5*(ym[0][0]+ym[0][1]);
 818             fXInit[1]=xm[1][1];
 819             fYInit[1]=0.5*(ym[0][0]+ym[1][1]);
 820             fQrInit[0]=0.5;
 821             fQrInit[1]=0.5;
 822             chi23=CombiDoubleMathiesonFit(c);
 823             if (chi23<10) {
 824                 Split(c);
 825                 Float_t yst;
 826                 yst = fYFit[0];
 827                 fYFit[0] = fYFit[1];
 828                 fYFit[1] = yst;
 829                 Split(c);
 830             }
 831         } else if (accepted[0]) {
 832             fXInit[0]=xm[0][0];
 833             fYInit[0]=ym[0][1];
 834             fXInit[1]=xm[1][1];
 835             fYInit[1]=ym[1][1];
 836             chi21=CombiDoubleMathiesonFit(c);
 837 //          Int_t ndf = fgNbins[0]+fgNbins[1]-6;
 838 //          Float_t prob = TMath::Prob(chi2,ndf);
 839 //          prob2->Fill(prob);
 840 //          chi2_2->Fill(chi21);
 841             if (fDebugLevel)
 842                 fprintf(stderr," chi2 %f\n",chi21);
 843             if (chi21<10) Split(c);
 844         } else if (accepted[1]) {
 845             fXInit[0]=xm[1][0];
 846             fYInit[0]=ym[1][1];
 847             fXInit[1]=xm[0][1];
 848             fYInit[1]=ym[0][1];
 849             chi22=CombiDoubleMathiesonFit(c);
 850 //          Int_t ndf = fgNbins[0]+fgNbins[1]-6;
 851 //          Float_t prob = TMath::Prob(chi2,ndf);
 852 //          prob2->Fill(prob);
 853 //          chi2_2->Fill(chi22);
 854             if (fDebugLevel)
 855                 fprintf(stderr," chi2 %f\n",chi22);
 856             if (chi22<10) Split(c);
 857         }
 858
 859         if (chi21 > 10 && chi22 > 10 && chi23 > 10) {
 860 //We keep only the combination found (X->cathode 2, Y->cathode 1)
 861             for (Int_t ico=0; ico<2; ico++) {
 862                 if (accepted[ico]) {
 863                     AliMUONRawCluster cnew;
 864                     Int_t cath;
 865                     for (cath=0; cath<2; cath++) {
 866                         cnew.SetX(cath, Float_t(xm[ico][1]));
 867                         cnew.SetY(cath, Float_t(ym[ico][0]));
 868                         cnew.SetZ(cath, fZPlane);
 869                         cnew.SetMultiplicity(cath, c->GetMultiplicity(cath));
 870                         for (i=0; i<fMul[cath]; i++) {
 871                             cnew.SetIndex(i, cath, c->GetIndex(i, cath));
 872                             fSeg[cath]->SetPad(fIx[i][cath], fIy[i][cath]);
 873                         }
 874                         fprintf(stderr,"\nRawCluster %d cath %d\n",ico,cath);
 875                         fprintf(stderr,"mult_av %d\n",c->GetMultiplicity(cath));
 876                         FillCluster(&cnew,cath);
 877                     }
 878                     cnew.SetClusterType(cnew.PhysicsContribution());
 879                     AddRawCluster(cnew);
 880                     fNPeaks++;
 881                 }
 882             }
 883         }
 884
 885 //  ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 886 //  (4) At least three local maxima on cathode 1 or on cathode 2
 887 //  ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 888     } else if (fNLocal[0]>2 || fNLocal[1]>2) {
 889         Int_t param = fNLocal[0]*fNLocal[1];
 890         Int_t ii;
 891
 892         Float_t ** xm = new Float_t * [param];
 893         for (ii=0; ii<param; ii++) xm[ii]=new Float_t [2];
 894         Float_t ** ym = new Float_t * [param];
 895         for (ii=0; ii<param; ii++) ym[ii]=new Float_t [2];
 896         Int_t ** ixm = new Int_t * [param];
 897         for (ii=0; ii<param; ii++) ixm[ii]=new Int_t [2];
 898         Int_t ** iym = new Int_t * [param];
 899         for (ii=0; ii<param; ii++) iym[ii]=new Int_t [2];
 900
 901         Int_t isec, ico;
 902         Float_t dpx, dpy, dx, dy;
 903
 904         ico=0;
 905         for (Int_t im1=0; im1<fNLocal[0]; im1++) {
 906             for (Int_t im2=0; im2<fNLocal[1]; im2++) {
 907                 xm[ico][0]=fX[fIndLocal[im1][0]][0];
 908                 ym[ico][0]=fY[fIndLocal[im1][0]][0];
 909                 xm[ico][1]=fX[fIndLocal[im2][1]][1];
 910                 ym[ico][1]=fY[fIndLocal[im2][1]][1];
 911
 912                 ixm[ico][0]=fIx[fIndLocal[im1][0]][0];
 913                 iym[ico][0]=fIy[fIndLocal[im1][0]][0];
 914                 ixm[ico][1]=fIx[fIndLocal[im2][1]][1];
 915                 iym[ico][1]=fIy[fIndLocal[im2][1]][1];
 916                 ico++;
 917             }
 918         }
 919
 920         Int_t nIco = ico;
 921         if (fDebugLevel)
 922             fprintf(stderr,"nIco %d\n",nIco);
 923         for (ico=0; ico<nIco; ico++) {
 924             if (fDebugLevel)
 925                 fprintf(stderr,"ico = %d\n",ico);
 926             isec=fSeg[0]->Sector(ixm[ico][0], iym[ico][0]);
 927             dpx=fSeg[0]->Dpx(isec)/2.;
 928             dx=TMath::Abs(xm[ico][0]-xm[ico][1]);
 929             isec=fSeg[1]->Sector(ixm[ico][1], iym[ico][1]);
 930             dpy=fSeg[1]->Dpy(isec)/2.;
 931             dy=TMath::Abs(ym[ico][0]-ym[ico][1]);
 932             if (fDebugLevel) {
 933                 fprintf(stderr,"dx %f dpx %f dy %f dpy %f\n",dx,dpx,dy,dpy);
 934                 fprintf(stderr,"  X %f Y %f\n",xm[ico][1],ym[ico][0]);
 935             }
 936             if ((dx <= dpx) && (dy <= dpy)) {
 937                 if (fDebugLevel)
 938                     fprintf(stderr,"ok\n");
 939                 Int_t cath;
 940                 AliMUONRawCluster cnew;
 941                 for (cath=0; cath<2; cath++) {
 942                     cnew.SetX(cath, Float_t(xm[ico][1]));
 943                     cnew.SetY(cath, Float_t(ym[ico][0]));
 944                     cnew.SetZ(cath, fZPlane);
 945                     cnew.SetMultiplicity(cath, c->GetMultiplicity(cath));
 946                     for (i=0; i<fMul[cath]; i++) {
 947                         cnew.SetIndex(i, cath, c->GetIndex(i, cath));
 948                         fSeg[cath]->SetPad(fIx[i][cath], fIy[i][cath]);
 949                     }
 950                     FillCluster(&cnew,cath);
 951                 }
 952                 cnew.SetClusterType(cnew.PhysicsContribution());
 953                 AddRawCluster(cnew);
 954                 fNPeaks++;
 955             }
 956         }
 957         delete [] xm;
 958         delete [] ym;
 959         delete [] ixm;
 960         delete [] iym;
 961     }
 962 }
 963
 964 void AliMUONClusterFinderVS::FindLocalMaxima(AliMUONRawCluster* /*c*/)
 965 {
 966 // Find all local maxima of a cluster
 967     if (fDebugLevel)
 968         printf("\n Find Local maxima  !");
 969
 970     AliMUONDigit* digt;
 971
 972     Int_t cath, cath1; // loops over cathodes
 973     Int_t i;           // loops over digits
 974     Int_t j;           // loops over cathodes
 975 //
 976 //  Find local maxima
 977 //
 978 //  counters for number of local maxima
 979     fNLocal[0]=fNLocal[1]=0;
 980 //  flags digits as local maximum
 981     Bool_t isLocal[100][2];
 982     for (i=0; i<100;i++) {
 983         isLocal[i][0]=isLocal[i][1]=kFALSE;
 984     }
 985 //  number of next neighbours and arrays to store them
 986     Int_t nn;
 987     Int_t x[10], y[10];
 988 // loop over cathodes
 989     for (cath=0; cath<2; cath++) {
 990 // loop over cluster digits
 991         for (i=0; i<fMul[cath]; i++) {
 992 // get neighbours for that digit and assume that it is local maximum
 993             fSeg[cath]->Neighbours(fIx[i][cath], fIy[i][cath], &nn, x, y);
 994             isLocal[i][cath]=kTRUE;
 995             Int_t isec= fSeg[cath]->Sector(fIx[i][cath], fIy[i][cath]);
 996             Float_t a0 = fSeg[cath]->Dpx(isec)*fSeg[cath]->Dpy(isec);
 997 // loop over next neighbours, if at least one neighbour has higher charger assumption
 998 // digit is not local maximum
 999             for (j=0; j<nn; j++) {
1000                 if (fHitMap[cath]->TestHit(x[j], y[j])==kEmpty) continue;
1001                 digt=(AliMUONDigit*) fHitMap[cath]->GetHit(x[j], y[j]);
1002                 isec=fSeg[cath]->Sector(x[j], y[j]);
1003                 Float_t a1 = fSeg[cath]->Dpx(isec)*fSeg[cath]->Dpy(isec);
1004                 if (digt->Signal()/a1 > fQ[i][cath]/a0) {
1005                     isLocal[i][cath]=kFALSE;
1006                     break;
1007 //
1008 // handle special case of neighbouring pads with equal signal
1009                 } else if (digt->Signal() == fQ[i][cath]) {
1010                     if (fNLocal[cath]>0) {
1011                         for (Int_t k=0; k<fNLocal[cath]; k++) {
1012                             if (x[j]==fIx[fIndLocal[k][cath]][cath]
1013                                 && y[j]==fIy[fIndLocal[k][cath]][cath])
1014                             {
1015                                 isLocal[i][cath]=kFALSE;
1016                             }
1017                         } // loop over local maxima
1018                     } // are there already local maxima
1019                 } // same charge ?
1020             } // loop over next neighbours
1021             if (isLocal[i][cath]) {
1022                 fIndLocal[fNLocal[cath]][cath]=i;
1023                 fNLocal[cath]++;
1024             }
1025         } // loop over all digits
1026     } // loop over cathodes
1027
1028     if (fDebugLevel) {
1029         printf("\n Found %d %d %d %d local Maxima\n",
1030                fNLocal[0], fNLocal[1], fMul[0], fMul[1]);
1031         fprintf(stderr,"\n Cathode 1 local Maxima %d Multiplicite %d\n",fNLocal[0], fMul[0]);
1032         fprintf(stderr," Cathode 2 local Maxima %d Multiplicite %d\n",fNLocal[1], fMul[1]);
1033     }
1034     Int_t ix, iy, isec;
1035     Float_t dpx, dpy;
1036
1037
1038     if (fNLocal[1]==2 &&  (fNLocal[0]==1 || fNLocal[0]==0)) {
1039         Int_t iback=fNLocal[0];
1040
1041 //  Two local maxima on cathode 2 and one maximum on cathode 1
1042 //  Look for local maxima considering up and down neighbours on the 1st cathode only
1043 //
1044 //  Loop over cluster digits
1045         cath=0;
1046         cath1=1;
1047
1048         for (i=0; i<fMul[cath]; i++) {
1049             isec=fSeg[cath]->Sector(fIx[i][cath],fIy[i][cath]);
1050             dpy=fSeg[cath]->Dpy(isec);
1051             dpx=fSeg[cath]->Dpx(isec);
1052             if (isLocal[i][cath]) continue;
1053 // Pad position should be consistent with position of local maxima on the opposite cathode
1054             if ((TMath::Abs(fX[i][cath]-fX[fIndLocal[0][cath1]][cath1]) > dpx/2.) &&
1055                 (TMath::Abs(fX[i][cath]-fX[fIndLocal[1][cath1]][cath1]) > dpx/2.))
1056                 continue;
1057
1058 // get neighbours for that digit and assume that it is local maximum
1059             isLocal[i][cath]=kTRUE;
1060 // compare signal to that on the two neighbours on the left and on the right
1061 // iNN counts the number of neighbours with signal, it should be 1 or 2
1062             Int_t iNN=0;
1063
1064             for (fSeg[cath]
1065                      ->FirstPad(fX[i][cath], fY[i][cath], fZPlane, 0., dpy);
1066                  fSeg[cath]
1067                      ->MorePads();
1068                  fSeg[cath]
1069                      ->NextPad())
1070             {
1071                 ix = fSeg[cath]->Ix();
1072                 iy = fSeg[cath]->Iy();
1073                 // skip the current pad
1074                 if (iy == fIy[i][cath]) continue;
1075
1076                 if (fHitMap[cath]->TestHit(ix, iy)!=kEmpty) {
1077                     iNN++;
1078                     digt=(AliMUONDigit*) fHitMap[cath]->GetHit(ix,iy);
1079                     if (digt->Signal() > fQ[i][cath]) isLocal[i][cath]=kFALSE;
1080                 }
1081             } // Loop over pad neighbours in y
1082             if (isLocal[i][cath] && iNN>0) {
1083                 fIndLocal[fNLocal[cath]][cath]=i;
1084                 fNLocal[cath]++;
1085             }
1086         } // loop over all digits
1087 // if one additional maximum has been found we are happy
1088 // if more maxima have been found restore the previous situation
1089         if (fDebugLevel) {
1090             fprintf(stderr,
1091                     "\n New search gives %d local maxima for cathode 1 \n",
1092                     fNLocal[0]);
1093             fprintf(stderr,
1094                     "                  %d local maxima for cathode 2 \n",
1095                     fNLocal[1]);
1096         }
1097         if (fNLocal[cath]>2) {
1098             fNLocal[cath]=iback;
1099         }
1100
1101     } // 1,2 local maxima
1102
1103     if (fNLocal[0]==2 &&  (fNLocal[1]==1 || fNLocal[1]==0)) {
1104         Int_t iback=fNLocal[1];
1105
1106 //  Two local maxima on cathode 1 and one maximum on cathode 2
1107 //  Look for local maxima considering left and right neighbours on the 2nd cathode only
1108         cath=1;
1109         Int_t cath1 = 0;
1110         Float_t eps = 1.e-5;
1111
1112 //
1113 //  Loop over cluster digits
1114         for (i=0; i<fMul[cath]; i++) {
1115             isec=fSeg[cath]->Sector(fIx[i][cath],fIy[i][cath]);
1116             dpx=fSeg[cath]->Dpx(isec);
1117             dpy=fSeg[cath]->Dpy(isec);
1118             if (isLocal[i][cath]) continue;
1119 // Pad position should be consistent with position of local maxima on the opposite cathode
1120             if ((TMath::Abs(fY[i][cath]-fY[fIndLocal[0][cath1]][cath1]) > dpy/2.+eps) &&
1121                 (TMath::Abs(fY[i][cath]-fY[fIndLocal[1][cath1]][cath1]) > dpy/2.+eps))
1122                 continue;
1123
1124 //
1125 // get neighbours for that digit and assume that it is local maximum
1126             isLocal[i][cath]=kTRUE;
1127 // compare signal to that on the two neighbours on the left and on the right
1128
1129 // iNN counts the number of neighbours with signal, it should be 1 or 2
1130             Int_t iNN=0;
1131             for (fSeg[cath]
1132                      ->FirstPad(fX[i][cath], fY[i][cath], fZPlane, dpx, 0.);
1133                  fSeg[cath]
1134                      ->MorePads();
1135                  fSeg[cath]
1136                      ->NextPad())
1137             {
1138
1139                 ix = fSeg[cath]->Ix();
1140                 iy = fSeg[cath]->Iy();
1141
1142                 // skip the current pad
1143                 if (ix == fIx[i][cath]) continue;
1144
1145                 if (fHitMap[cath]->TestHit(ix, iy)!=kEmpty) {
1146                     iNN++;
1147                     digt=(AliMUONDigit*) fHitMap[cath]->GetHit(ix,iy);
1148                     if (digt->Signal() > fQ[i][cath]) isLocal[i][cath]=kFALSE;
1149                 }
1150             } // Loop over pad neighbours in x
1151             if (isLocal[i][cath] && iNN>0) {
1152                 fIndLocal[fNLocal[cath]][cath]=i;
1153                 fNLocal[cath]++;
1154             }
1155         } // loop over all digits
1156 // if one additional maximum has been found we are happy
1157 // if more maxima have been found restore the previous situation
1158         if (fDebugLevel) {
1159             fprintf(stderr,"\n New search gives %d local maxima for cathode 1 \n",fNLocal[0]);
1160             fprintf(stderr,"\n                  %d local maxima for cathode 2 \n",fNLocal[1]);
1161             printf("\n New search gives %d %d \n",fNLocal[0],fNLocal[1]);
1162         }
1163         if (fNLocal[cath]>2) {
1164             fNLocal[cath]=iback;
1165         }
1166     } // 2,1 local maxima
1167 }
1168
1169
1170 void  AliMUONClusterFinderVS::FillCluster(AliMUONRawCluster* c, Int_t flag, Int_t cath)
1171 {
1172 //
1173 //  Completes cluster information starting from list of digits
1174 //
1175     AliMUONDigit* dig;
1176     Float_t x, y, z;
1177     Int_t  ix, iy;
1178
1179     if (cath==1) {
1180         c->SetPeakSignal(cath,c->GetPeakSignal(0));
1181     } else {
1182         c->SetPeakSignal(cath,0);
1183     }
1184
1185
1186     if (flag) {
1187         c->SetX(cath,0.);
1188         c->SetY(cath,0.);
1189         c->SetCharge(cath,0);
1190     }
1191
1192     if (fDebugLevel)
1193         fprintf(stderr,"\n fPeakSignal %d\n",c->GetPeakSignal(cath));
1194     for (Int_t i=0; i<c->GetMultiplicity(cath); i++)
1195     {
1196         dig= fInput->Digit(cath,c->GetIndex(i,cath));
1197         ix=dig->PadX()+c->GetOffset(i,cath);
1198         iy=dig->PadY();
1199         Int_t q=dig->Signal();
1200         if (!flag) q=Int_t(q*c->GetContrib(i,cath));
1201 //      fprintf(stderr,"q %d c->fPeakSignal[ %d ] %d\n",q,cath,c->fPeakSignal[cath]);
1202         if (dig->Physics() >= dig->Signal()) {
1203             c->SetPhysics(i,2);
1204         } else if (dig->Physics() == 0) {
1205             c->SetPhysics(i,0);
1206         } else  c->SetPhysics(i,1);
1207 //
1208 //
1209         if (fDebugLevel>1)
1210             fprintf(stderr,"q %d c->fPeakSignal[cath] %d\n",q,c->GetPeakSignal(cath));
1211 // peak signal and track list
1212         if (q>c->GetPeakSignal(cath)) {
1213             c->SetPeakSignal(cath, q);
1214             c->SetTrack(0,dig->Hit());
1215             c->SetTrack(1,dig->Track(0));
1216             c->SetTrack(2,dig->Track(1));
1217 //          fprintf(stderr," c->fTracks[0] %d c->fTracks[1] %d\n",dig->fHit,dig->fTracks[0]);
1218         }
1219 //
1220         if (flag) {
1221             fSeg[cath]->GetPadC(ix, iy, x, y, z);
1222             c->AddX(cath, q*x);
1223             c->AddY(cath, q*y);
1224             c->AddCharge(cath, q);
1225         }
1226     } // loop over digits
1227     if (fDebugLevel)
1228         fprintf(stderr," fin du cluster c\n");
1229
1230
1231     if (flag) {
1232         c->SetX(cath, c->GetX(cath)/c->GetCharge(cath));
1233 // Force on anod
1234         c->SetX(cath, fSeg[cath]->GetAnod(c->GetX(cath)));
1235         c->SetY(cath, c->GetY(cath)/c->GetCharge(cath));
1236 //
1237 //  apply correction to the coordinate along the anode wire
1238 //
1239         x=c->GetX(cath);
1240         y=c->GetY(cath);
1241         fSeg[cath]->GetPadI(x, y, fZPlane, ix, iy);
1242         fSeg[cath]->GetPadC(ix, iy, x, y, z);
1243         Int_t isec=fSeg[cath]->Sector(ix,iy);
1244         TF1* cogCorr = fSeg[cath]->CorrFunc(isec-1);
1245
1246         if (cogCorr) {
1247             Float_t yOnPad=(c->GetY(cath)-y)/fSeg[cath]->Dpy(isec);
1248             c->SetY(cath, c->GetY(cath)-cogCorr->Eval(yOnPad, 0, 0));
1249         }
1250     }
1251 }
1252
1253 void  AliMUONClusterFinderVS::FillCluster(AliMUONRawCluster* c, Int_t cath)
1254 {
1255 //
1256 //  Completes cluster information starting from list of digits
1257 //
1258     static Float_t dr0;
1259
1260     AliMUONDigit* dig;
1261
1262     if (cath==0) {
1263         dr0 = 10000;
1264     }
1265
1266     Float_t xpad, ypad, zpad;
1267     Float_t dx, dy, dr;
1268
1269     for (Int_t i=0; i<c->GetMultiplicity(cath); i++)
1270     {
1271         dig = fInput->Digit(cath,c->GetIndex(i,cath));
1272         fSeg[cath]->
1273         GetPadC(dig->PadX(),dig->PadY(),xpad,ypad, zpad);
1274         if (fDebugLevel)
1275             fprintf(stderr,"x %f y %f cx %f cy %f\n",xpad,ypad,c->GetX(0),c->GetY(0));
1276         dx = xpad - c->GetX(0);
1277         dy = ypad - c->GetY(0);
1278         dr = TMath::Sqrt(dx*dx+dy*dy);
1279
1280         if (dr < dr0) {
1281             dr0 = dr;
1282             if (fDebugLevel)
1283                 fprintf(stderr," dr %f\n",dr);
1284             Int_t q=dig->Signal();
1285             if (dig->Physics() >= dig->Signal()) {
1286                 c->SetPhysics(i,2);
1287             } else if (dig->Physics() == 0) {
1288                 c->SetPhysics(i,0);
1289             } else  c->SetPhysics(i,1);
1290             c->SetPeakSignal(cath,q);
1291             c->SetTrack(0,dig->Hit());
1292             c->SetTrack(1,dig->Track(0));
1293             c->SetTrack(2,dig->Track(1));
1294             if (fDebugLevel)
1295                 fprintf(stderr," c->fTracks[0] %d c->fTracks[1] %d\n",dig->Hit(),
1296                     dig->Track(0));
1297         }
1298 //
1299     } // loop over digits
1300
1301 //  apply correction to the coordinate along the anode wire
1302 // Force on anod
1303     c->SetX(cath,fSeg[cath]->GetAnod(c->GetX(cath)));
1304 }
1305
1306 void  AliMUONClusterFinderVS::FindCluster(Int_t i, Int_t j, Int_t cath, AliMUONRawCluster &c){
1307
1308
1309 //
1310 //  Find a super cluster on both cathodes
1311 //
1312 //
1313 //  Add i,j as element of the cluster
1314 //
1315
1316     Int_t idx = fHitMap[cath]->GetHitIndex(i,j);
1317     AliMUONDigit* dig = (AliMUONDigit*) fHitMap[cath]->GetHit(i,j);
1318     Int_t q=dig->Signal();
1319     Int_t theX=dig->PadX();
1320     Int_t theY=dig->PadY();
1321
1322     if (q > TMath::Abs(c.GetPeakSignal(0)) && q > TMath::Abs(c.GetPeakSignal(1))) {
1323         c.SetPeakSignal(cath,q);
1324         c.SetTrack(0,dig->Hit());
1325         c.SetTrack(1,dig->Track(0));
1326         c.SetTrack(2,dig->Track(1));
1327     }
1328
1329 //
1330 //  Make sure that list of digits is ordered
1331 //
1332     Int_t mu=c.GetMultiplicity(cath);
1333     c.SetIndex(mu, cath, idx);
1334
1335     if (dig->Physics() >= dig->Signal()) {
1336         c.SetPhysics(mu,2);
1337     } else if (dig->Physics() == 0) {
1338         c.SetPhysics(mu,0);
1339     } else  c.SetPhysics(mu,1);
1340
1341
1342     if (mu > 0) {
1343         for (Int_t ind = mu-1; ind >= 0; ind--) {
1344             Int_t ist=c.GetIndex(ind,cath);
1345             Int_t ql=fInput->Digit(cath, ist)->Signal();
1346             Int_t ix=fInput->Digit(cath, ist)->PadX();
1347             Int_t iy=fInput->Digit(cath, ist)->PadY();
1348
1349             if (q>ql || (q==ql && theX > ix && theY < iy)) {
1350                 c.SetIndex(ind, cath, idx);
1351                 c.SetIndex(ind+1, cath, ist);
1352             } else {
1353
1354                 break;
1355             }
1356         }
1357     }
1358
1359     c.SetMultiplicity(cath, c.GetMultiplicity(cath)+1);
1360     if (c.GetMultiplicity(cath) >= 50 ) {
1361         printf("FindCluster - multiplicity >50  %d \n",c.GetMultiplicity(0));
1362         c.SetMultiplicity(cath, 49);
1363     }
1364
1365 // Prepare center of gravity calculation
1366     Float_t x, y, z;
1367     fSeg[cath]->GetPadC(i, j, x, y, z);
1368
1369     c.AddX(cath,q*x);
1370     c.AddY(cath,q*y);
1371     c.AddCharge(cath,q);
1372 //
1373 // Flag hit as "taken"
1374     fHitMap[cath]->FlagHit(i,j);
1375 //
1376 //  Now look recursively for all neighbours and pad hit on opposite cathode
1377 //
1378 //  Loop over neighbours
1379     Int_t ix,iy;
1380     ix=iy=0;
1381     Int_t nn;
1382     Int_t xList[10], yList[10];
1383     fSeg[cath]->Neighbours(i,j,&nn,xList,yList);
1384     for (Int_t in=0; in<nn; in++) {
1385         ix=xList[in];
1386         iy=yList[in];
1387
1388         if (fHitMap[cath]->TestHit(ix,iy)==kUnused) {
1389             if (fDebugLevel>1)
1390                 printf("\n Neighbours %d %d %d", cath, ix, iy);
1391             FindCluster(ix, iy, cath, c);
1392         }
1393
1394    }
1395     Int_t nOpp=0;
1396     Int_t iXopp[50], iYopp[50];
1397
1398 //  Neighbours on opposite cathode
1399 //  Take into account that several pads can overlap with the present pad
1400     Int_t isec=fSeg[cath]->Sector(i,j);
1401     Int_t iop;
1402     Float_t dx, dy;
1403
1404     if (cath==0) {
1405         iop = 1;
1406         dx  = (fSeg[cath]->Dpx(isec))/2.;
1407         dy  = 0.;
1408     } else {
1409         iop = 0;
1410         dx  = 0.;
1411         dy  = (fSeg[cath]->Dpy(isec))/2;
1412     }
1413 // loop over pad neighbours on opposite cathode
1414     for (fSeg[iop]->FirstPad(x, y, fZPlane, dx, dy);
1415          fSeg[iop]->MorePads();
1416          fSeg[iop]->NextPad())
1417     {
1418
1419         ix = fSeg[iop]->Ix(); iy = fSeg[iop]->Iy();
1420         if (fDebugLevel > 1)
1421             printf("\n ix, iy: %f %f %f %d %d %d", x,y,z,ix, iy, fSector);
1422         if (fHitMap[iop]->TestHit(ix,iy)==kUnused){
1423             iXopp[nOpp]=ix;
1424             iYopp[nOpp++]=iy;
1425             if (fDebugLevel > 1)
1426                 printf("\n Opposite %d %d %d", iop, ix, iy);
1427         }
1428
1429     } // Loop over pad neighbours
1430 //  This had to go outside the loop since recursive calls inside the iterator are not possible
1431 //
1432     Int_t jopp;
1433     for (jopp=0; jopp<nOpp; jopp++) {
1434         if (fHitMap[iop]->TestHit(iXopp[jopp],iYopp[jopp]) == kUnused)
1435             FindCluster(iXopp[jopp], iYopp[jopp], iop, c);
1436     }
1437 }
1438
1439 //_____________________________________________________________________________
1440
1441 void AliMUONClusterFinderVS::FindRawClusters()
1442 {
1443   //
1444   // MUON cluster finder from digits -- finds neighbours on both cathodes and
1445   // fills the tree with raw clusters
1446   //
1447
1448     ResetRawClusters();
1449 //  Return if no input datad available
1450     if (!fInput->NDigits(0) && !fInput->NDigits(1)) return;
1451
1452     fSeg[0] = fInput->Segmentation(0);
1453     fSeg[1] = fInput->Segmentation(1);
1454
1455     fHitMap[0]  = new AliMUONHitMapA1(fSeg[0], fInput->Digits(0));
1456     fHitMap[1]  = new AliMUONHitMapA1(fSeg[1], fInput->Digits(1));
1457
1458
1459     AliMUONDigit *dig;
1460
1461     Int_t ndig, cath;
1462     Int_t nskip=0;
1463     Int_t ncls=0;
1464     fHitMap[0]->FillHits();
1465     fHitMap[1]->FillHits();
1466 //
1467 //  Outer Loop over Cathodes
1468     for (cath=0; cath<2; cath++) {
1469         for (ndig=0; ndig<fInput->NDigits(cath); ndig++) {
1470             dig = fInput->Digit(cath, ndig);
1471             Int_t i=dig->PadX();
1472             Int_t j=dig->PadY();
1473             if (fHitMap[cath]->TestHit(i,j)==kUsed ||fHitMap[0]->TestHit(i,j)==kEmpty) {
1474                 nskip++;
1475                 continue;
1476             }
1477             if (fDebugLevel)
1478                 fprintf(stderr,"\n CATHODE %d CLUSTER %d\n",cath,ncls);
1479             AliMUONRawCluster c;
1480             c.SetMultiplicity(0, 0);
1481             c.SetMultiplicity(1, 0);
1482             c.SetPeakSignal(cath,dig->Signal());
1483             c.SetTrack(0, dig->Hit());
1484             c.SetTrack(1, dig->Track(0));
1485             c.SetTrack(2, dig->Track(1));
1486             // tag the beginning of cluster list in a raw cluster
1487             c.SetNcluster(0,-1);
1488             Float_t xcu, ycu;
1489             fSeg[cath]->GetPadC(i,j,xcu, ycu, fZPlane);
1490             fSector= fSeg[cath]->Sector(i,j)/100;
1491             if (fDebugLevel)
1492                 printf("\n New Seed %d %d ", i,j);
1493
1494
1495             FindCluster(i,j,cath,c);
1496 //          ^^^^^^^^^^^^^^^^^^^^^^^^
1497             // center of gravity
1498             if (c.GetX(0)!=0.) c.SetX(0, c.GetX(0)/c.GetCharge(0)); // c.fX[0] /= c.fQ[0];
1499 // Force on anod
1500             c.SetX(0,fSeg[0]->GetAnod(c.GetX(0)));
1501             if (c.GetY(0)!=0.) c.SetY(0, c.GetY(0)/c.GetCharge(0)); // c.fY[0] /= c.fQ[0];
1502
1503             if(c.GetCharge(1)!=0.) c.SetX(1, c.GetX(1)/c.GetCharge(1));  // c.fX[1] /= c.fQ[1];
1504
1505            // Force on anod
1506             c.SetX(1, fSeg[0]->GetAnod(c.GetX(1)));
1507             if(c.GetCharge(1)!=0.) c.SetY(1, c.GetY(1)/c.GetCharge(1));// c.fY[1] /= c.fQ[1];
1508
1509             c.SetZ(0, fZPlane);
1510             c.SetZ(1, fZPlane);
1511
1512             if (fDebugLevel) {
1513                 fprintf(stderr,"\n Cathode 1 multiplicite %d X(CG) %f Y(CG) %f\n",
1514                         c.GetMultiplicity(0),c.GetX(0),c.GetY(0));
1515                 fprintf(stderr," Cathode 2 multiplicite %d X(CG) %f Y(CG) %f\n",
1516                         c.GetMultiplicity(1),c.GetX(1),c.GetY(1));
1517             }
1518 //      Analyse cluster and decluster if necessary
1519 //
1520         ncls++;
1521         c.SetNcluster(1,fNRawClusters);
1522         c.SetClusterType(c.PhysicsContribution());
1523
1524         fNPeaks=0;
1525 //
1526 //
1527         Decluster(&c);
1528 //
1529 //      reset Cluster object
1530         { // begin local scope
1531             for (int k=0;k<c.GetMultiplicity(0);k++) c.SetIndex(k, 0, 0);
1532         } // end local scope
1533
1534         { // begin local scope
1535             for (int k=0;k<c.GetMultiplicity(1);k++) c.SetIndex(k, 1, 0);
1536         } // end local scope
1537
1538         c.SetMultiplicity(0,0);
1539         c.SetMultiplicity(1,0);
1540
1541
1542         } // end loop ndig
1543     } // end loop cathodes
1544     delete fHitMap[0];
1545     delete fHitMap[1];
1546 }
1547
1548 Float_t AliMUONClusterFinderVS::SingleMathiesonFit(AliMUONRawCluster *c, Int_t cath)
1549 {
1550 // Performs a single Mathieson fit on one cathode
1551 //
1552     Double_t arglist[20];
1553     Int_t ierflag=0;
1554     AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
1555
1556     clusterInput.Fitter()->SetFCN(fcnS1);
1557     clusterInput.Fitter()->mninit(2,10,7);
1558     clusterInput.Fitter()->SetPrintLevel(-1+fDebugLevel);
1559     arglist[0]=-1;
1560     clusterInput.Fitter()->mnexcm("SET NOW", arglist, 0, ierflag);
1561 // Set starting values
1562     static Double_t vstart[2];
1563     vstart[0]=c->GetX(1);
1564     vstart[1]=c->GetY(0);
1565
1566
1567 // lower and upper limits
1568     static Double_t lower[2], upper[2];
1569     Int_t ix,iy;
1570     fSeg[cath]->GetPadI(c->GetX(cath), c->GetY(cath), fZPlane, ix, iy);
1571     Int_t isec=fSeg[cath]->Sector(ix, iy);
1572     lower[0]=vstart[0]-fSeg[cath]->Dpx(isec)/2;
1573     lower[1]=vstart[1]-fSeg[cath]->Dpy(isec)/2;
1574
1575     upper[0]=lower[0]+fSeg[cath]->Dpx(isec);
1576     upper[1]=lower[1]+fSeg[cath]->Dpy(isec);
1577
1578 // step sizes
1579     static Double_t step[2]={0.0005, 0.0005};
1580
1581     clusterInput.Fitter()->mnparm(0,"x1",vstart[0],step[0],lower[0],upper[0],ierflag);
1582     clusterInput.Fitter()->mnparm(1,"y1",vstart[1],step[1],lower[1],upper[1],ierflag);
1583 // ready for minimisation
1584     arglist[0]= -1;
1585     arglist[1]= 0;
1586
1587     clusterInput.Fitter()->mnexcm("SET NOGR", arglist, 0, ierflag);
1588     clusterInput.Fitter()->mnexcm("MIGRAD", arglist, 0, ierflag);
1589     //    clusterInput.Fitter()->mnexcm("EXIT" , arglist, 0, ierflag);
1590     Double_t fmin, fedm, errdef;
1591     Int_t   npari, nparx, istat;
1592
1593     clusterInput.Fitter()->mnstat(fmin, fedm, errdef, npari, nparx, istat);
1594     fFitStat=istat;
1595
1596 // Print results
1597 // Get fitted parameters
1598     Double_t xrec, yrec;
1599     TString chname;
1600     Double_t epxz, b1, b2;
1601     Int_t ierflg;
1602     clusterInput.Fitter()->mnpout(0, chname, xrec, epxz, b1, b2, ierflg);
1603     clusterInput.Fitter()->mnpout(1, chname, yrec, epxz, b1, b2, ierflg);
1604     fXFit[cath]=xrec;
1605     fYFit[cath]=yrec;
1606     return fmin;
1607 }
1608
1609 Float_t AliMUONClusterFinderVS::CombiSingleMathiesonFit(AliMUONRawCluster * /*c*/)
1610 {
1611 // Perform combined Mathieson fit on both cathode planes
1612 //
1613     Double_t arglist[20];
1614     Int_t ierflag=0;
1615     AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
1616     clusterInput.Fitter()->SetFCN(fcnCombiS1);
1617     clusterInput.Fitter()->mninit(2,10,7);
1618     clusterInput.Fitter()->SetPrintLevel(-1+fDebugLevel);
1619     arglist[0]=-1;
1620     clusterInput.Fitter()->mnexcm("SET NOW", arglist, 0, ierflag);
1621     static Double_t vstart[2];
1622     vstart[0]=fXInit[0];
1623     vstart[1]=fYInit[0];
1624
1625
1626 // lower and upper limits
1627     static Float_t lower[2], upper[2];
1628     Int_t ix,iy,isec;
1629     fSeg[0]->GetPadI(fXInit[0], fYInit[0], fZPlane, ix, iy);
1630     isec=fSeg[0]->Sector(ix, iy);
1631     Float_t dpy=fSeg[0]->Dpy(isec);
1632     fSeg[1]->GetPadI(fXInit[0], fYInit[0], fZPlane, ix, iy);
1633     isec=fSeg[1]->Sector(ix, iy);
1634     Float_t dpx=fSeg[1]->Dpx(isec);
1635
1636     Int_t icount;
1637     Float_t xdum, ydum, zdum;
1638
1639 //  Find save upper and lower limits
1640
1641     icount = 0;
1642
1643     for (fSeg[1]->FirstPad(fXInit[0], fYInit[0], fZPlane, dpx, 0.);
1644          fSeg[1]->MorePads(); fSeg[1]->NextPad())
1645     {
1646         ix=fSeg[1]->Ix(); iy=fSeg[1]->Iy();
1647         fSeg[1]->GetPadC(ix,iy, upper[0], ydum, zdum);
1648         if (icount ==0) lower[0]=upper[0];
1649         icount++;
1650     }
1651
1652     if (lower[0]>upper[0]) {xdum=lower[0]; lower[0]=upper[0]; upper[0]=xdum;}
1653
1654     icount=0;
1655     if (fDebugLevel)
1656         printf("\n single y %f %f", fXInit[0], fYInit[0]);
1657
1658     for (fSeg[0]->FirstPad(fXInit[0], fYInit[0], fZPlane, 0., dpy);
1659          fSeg[0]->MorePads(); fSeg[0]->NextPad())
1660     {
1661         ix=fSeg[0]->Ix(); iy=fSeg[0]->Iy();
1662         fSeg[0]->GetPadC(ix,iy,xdum,upper[1],zdum);
1663         if (icount ==0) lower[1]=upper[1];
1664         icount++;
1665         if (fDebugLevel)
1666             printf("\n upper lower %d %f %f", icount, upper[1], lower[1]);
1667     }
1668
1669     if (lower[1]>upper[1]) {xdum=lower[1]; lower[1]=upper[1]; upper[1]=xdum;}
1670
1671 // step sizes
1672     static Double_t step[2]={0.00001, 0.0001};
1673
1674     clusterInput.Fitter()->mnparm(0,"x1",vstart[0],step[0],lower[0],upper[0],ierflag);
1675     clusterInput.Fitter()->mnparm(1,"y1",vstart[1],step[1],lower[1],upper[1],ierflag);
1676 // ready for minimisation
1677     arglist[0]= -1;
1678     arglist[1]= 0;
1679
1680     clusterInput.Fitter()->mnexcm("SET NOGR", arglist, 0, ierflag);
1681     clusterInput.Fitter()->mnexcm("MIGRAD", arglist, 0, ierflag);
1682     //    clusterInput.Fitter()->mnexcm("EXIT" , arglist, 0, ierflag);
1683     Double_t fmin, fedm, errdef;
1684     Int_t   npari, nparx, istat;
1685
1686     clusterInput.Fitter()->mnstat(fmin, fedm, errdef, npari, nparx, istat);
1687     fFitStat=istat;
1688
1689 // Print results
1690 // Get fitted parameters
1691     Double_t xrec, yrec;
1692     TString chname;
1693     Double_t epxz, b1, b2;
1694     Int_t ierflg;
1695     clusterInput.Fitter()->mnpout(0, chname, xrec, epxz, b1, b2, ierflg);
1696     clusterInput.Fitter()->mnpout(1, chname, yrec, epxz, b1, b2, ierflg);
1697     fXFit[0]=xrec;
1698     fYFit[0]=yrec;
1699     return fmin;
1700 }
1701
1702 Bool_t AliMUONClusterFinderVS::DoubleMathiesonFit(AliMUONRawCluster * /*c*/, Int_t cath)
1703 {
1704 // Performs a double Mathieson fit on one cathode
1705 //
1706
1707 //
1708 //  Initialise global variables for fit
1709     Double_t arglist[20];
1710     Int_t ierflag=0;
1711     AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
1712     clusterInput.Fitter()->SetFCN(fcnS2);
1713     clusterInput.Fitter()->mninit(5,10,7);
1714     clusterInput.Fitter()->SetPrintLevel(-1+fDebugLevel);
1715     arglist[0]=-1;
1716     clusterInput.Fitter()->mnexcm("SET NOW", arglist, 0, ierflag);
1717 // Set starting values
1718     static Double_t vstart[5];
1719     vstart[0]=fX[fIndLocal[0][cath]][cath];
1720     vstart[1]=fY[fIndLocal[0][cath]][cath];
1721     vstart[2]=fX[fIndLocal[1][cath]][cath];
1722     vstart[3]=fY[fIndLocal[1][cath]][cath];
1723     vstart[4]=Float_t(fQ[fIndLocal[0][cath]][cath])/
1724         Float_t(fQ[fIndLocal[0][cath]][cath]+fQ[fIndLocal[1][cath]][cath]);
1725 // lower and upper limits
1726     static Float_t lower[5], upper[5];
1727     Int_t isec=fSeg[cath]->Sector(fIx[fIndLocal[0][cath]][cath], fIy[fIndLocal[0][cath]][cath]);
1728     lower[0]=vstart[0]-fSeg[cath]->Dpx(isec);
1729     lower[1]=vstart[1]-fSeg[cath]->Dpy(isec);
1730
1731     upper[0]=lower[0]+2.*fSeg[cath]->Dpx(isec);
1732     upper[1]=lower[1]+2.*fSeg[cath]->Dpy(isec);
1733
1734     isec=fSeg[cath]->Sector(fIx[fIndLocal[1][cath]][cath], fIy[fIndLocal[1][cath]][cath]);
1735     lower[2]=vstart[2]-fSeg[cath]->Dpx(isec)/2;
1736     lower[3]=vstart[3]-fSeg[cath]->Dpy(isec)/2;
1737
1738     upper[2]=lower[2]+fSeg[cath]->Dpx(isec);
1739     upper[3]=lower[3]+fSeg[cath]->Dpy(isec);
1740
1741     lower[4]=0.;
1742     upper[4]=1.;
1743 // step sizes
1744     static Double_t step[5]={0.0005, 0.0005, 0.0005, 0.0005, 0.0001};
1745
1746     clusterInput.Fitter()->mnparm(0,"x1",vstart[0],step[0],lower[0],upper[0],ierflag);
1747     clusterInput.Fitter()->mnparm(1,"y1",vstart[1],step[1],lower[1],upper[1],ierflag);
1748     clusterInput.Fitter()->mnparm(2,"x2",vstart[2],step[2],lower[2],upper[2],ierflag);
1749     clusterInput.Fitter()->mnparm(3,"y2",vstart[3],step[3],lower[3],upper[3],ierflag);
1750     clusterInput.Fitter()->mnparm(4,"a0",vstart[4],step[4],lower[4],upper[4],ierflag);
1751 // ready for minimisation
1752     arglist[0]= -1;
1753     arglist[1]= 0;
1754
1755     clusterInput.Fitter()->mnexcm("SET NOGR", arglist, 0, ierflag);
1756     clusterInput.Fitter()->mnexcm("MIGRAD", arglist, 0, ierflag);
1757     //    clusterInput.Fitter()->mnexcm("EXIT" , arglist, 0, ierflag);
1758 // Get fitted parameters
1759     Double_t xrec[2], yrec[2], qfrac;
1760     TString chname;
1761     Double_t epxz, b1, b2;
1762     Int_t ierflg;
1763     clusterInput.Fitter()->mnpout(0, chname, xrec[0], epxz, b1, b2, ierflg);
1764     clusterInput.Fitter()->mnpout(1, chname, yrec[0], epxz, b1, b2, ierflg);
1765     clusterInput.Fitter()->mnpout(2, chname, xrec[1], epxz, b1, b2, ierflg);
1766     clusterInput.Fitter()->mnpout(3, chname, yrec[1], epxz, b1, b2, ierflg);
1767     clusterInput.Fitter()->mnpout(4, chname, qfrac,   epxz, b1, b2, ierflg);
1768
1769     Double_t fmin, fedm, errdef;
1770     Int_t   npari, nparx, istat;
1771
1772     clusterInput.Fitter()->mnstat(fmin, fedm, errdef, npari, nparx, istat);
1773     fFitStat=istat;
1774     return kTRUE;
1775 }
1776
1777 Float_t AliMUONClusterFinderVS::CombiDoubleMathiesonFit(AliMUONRawCluster * /*c*/)
1778 {
1779 //
1780 // Perform combined double Mathieson fit on both cathode planes
1781 //
1782     Double_t arglist[20];
1783     Int_t ierflag=0;
1784     AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
1785     clusterInput.Fitter()->SetFCN(fcnCombiS2);
1786     clusterInput.Fitter()->mninit(6,10,7);
1787     clusterInput.Fitter()->SetPrintLevel(-1+fDebugLevel);
1788     arglist[0]=-1;
1789     clusterInput.Fitter()->mnexcm("SET NOW", arglist, 0, ierflag);
1790 // Set starting values
1791     static Double_t vstart[6];
1792     vstart[0]=fXInit[0];
1793     vstart[1]=fYInit[0];
1794     vstart[2]=fXInit[1];
1795     vstart[3]=fYInit[1];
1796     vstart[4]=fQrInit[0];
1797     vstart[5]=fQrInit[1];
1798 // lower and upper limits
1799     static Float_t lower[6], upper[6];
1800     Int_t ix,iy,isec;
1801     Float_t dpx, dpy;
1802
1803     fSeg[1]->GetPadI(fXInit[0], fYInit[0], fZPlane, ix, iy);
1804     isec=fSeg[1]->Sector(ix, iy);
1805     dpx=fSeg[1]->Dpx(isec);
1806
1807     fSeg[0]->GetPadI(fXInit[0], fYInit[0], fZPlane, ix, iy);
1808     isec=fSeg[0]->Sector(ix, iy);
1809     dpy=fSeg[0]->Dpy(isec);
1810
1811
1812     Int_t icount;
1813     Float_t xdum, ydum, zdum;
1814     if (fDebugLevel)
1815         printf("\n Cluster Finder: %f %f %f %f  ", fXInit[0], fXInit[1],fYInit[0], fYInit[1] );
1816
1817 //  Find save upper and lower limits
1818     icount = 0;
1819
1820     for (fSeg[1]->FirstPad(fXInit[0], fYInit[0], fZPlane, dpx, 0.);
1821          fSeg[1]->MorePads(); fSeg[1]->NextPad())
1822     {
1823         ix=fSeg[1]->Ix(); iy=fSeg[1]->Iy();
1824 //      if (fHitMap[1]->TestHit(ix, iy) == kEmpty) continue;
1825         fSeg[1]->GetPadC(ix,iy,upper[0],ydum,zdum);
1826         if (icount ==0) lower[0]=upper[0];
1827         icount++;
1828     }
1829     if (lower[0]>upper[0]) {xdum=lower[0]; lower[0]=upper[0]; upper[0]=xdum;}
1830 //    vstart[0] = 0.5*(lower[0]+upper[0]);
1831
1832
1833     icount=0;
1834
1835     for (fSeg[0]->FirstPad(fXInit[0], fYInit[0], fZPlane, 0., dpy);
1836          fSeg[0]->MorePads(); fSeg[0]->NextPad())
1837     {
1838         ix=fSeg[0]->Ix(); iy=fSeg[0]->Iy();
1839 //      if (fHitMap[0]->TestHit(ix, iy) == kEmpty) continue;
1840         fSeg[0]->GetPadC(ix,iy,xdum,upper[1],zdum);
1841         if (icount ==0) lower[1]=upper[1];
1842         icount++;
1843     }
1844
1845     if (lower[1]>upper[1]) {xdum=lower[1]; lower[1]=upper[1]; upper[1]=xdum;}
1846 //     vstart[1] = 0.5*(lower[1]+upper[1]);
1847
1848
1849     fSeg[1]->GetPadI(fXInit[1], fYInit[1], fZPlane, ix, iy);
1850     isec=fSeg[1]->Sector(ix, iy);
1851     dpx=fSeg[1]->Dpx(isec);
1852     fSeg[0]->GetPadI(fXInit[1], fYInit[1], fZPlane, ix, iy);
1853     isec=fSeg[0]->Sector(ix, iy);
1854     dpy=fSeg[0]->Dpy(isec);
1855
1856
1857 //  Find save upper and lower limits
1858
1859     icount=0;
1860
1861     for (fSeg[1]->FirstPad(fXInit[1], fYInit[1], fZPlane, dpx, 0);
1862          fSeg[1]->MorePads(); fSeg[1]->NextPad())
1863     {
1864         ix=fSeg[1]->Ix(); iy=fSeg[1]->Iy();
1865 //      if (fHitMap[1]->TestHit(ix, iy) == kEmpty) continue;
1866         fSeg[1]->GetPadC(ix,iy,upper[2],ydum,zdum);
1867         if (icount ==0) lower[2]=upper[2];
1868         icount++;
1869     }
1870     if (lower[2]>upper[2]) {xdum=lower[2]; lower[2]=upper[2]; upper[2]=xdum;}
1871     //    vstart[2] = 0.5*(lower[2]+upper[2]);
1872
1873     icount=0;
1874
1875     for (fSeg[0]->FirstPad(fXInit[1], fYInit[1], fZPlane, 0, dpy);
1876          fSeg[0]-> MorePads(); fSeg[0]->NextPad())
1877     {
1878         ix=fSeg[0]->Ix(); iy=fSeg[0]->Iy();
1879 //      if (fHitMap[0]->TestHit(ix, iy) != kEmpty) continue;
1880
1881         fSeg[0]->GetPadC(ix,iy,xdum,upper[3],zdum);
1882         if (icount ==0) lower[3]=upper[3];
1883         icount++;
1884
1885     }
1886     if (lower[3]>upper[3]) {xdum=lower[3]; lower[3]=upper[3]; upper[3]=xdum;}
1887
1888 //     vstart[3] = 0.5*(lower[3]+upper[3]);
1889
1890     lower[4]=0.;
1891     upper[4]=1.;
1892     lower[5]=0.;
1893     upper[5]=1.;
1894
1895 // step sizes
1896     static Double_t step[6]={0.0005, 0.0005, 0.0005, 0.0005, 0.001, 0.001};
1897     clusterInput.Fitter()->mnparm(0,"x1",vstart[0],step[0],lower[0],upper[0],ierflag);
1898     clusterInput.Fitter()->mnparm(1,"y1",vstart[1],step[1],lower[1],upper[1],ierflag);
1899     clusterInput.Fitter()->mnparm(2,"x2",vstart[2],step[2],lower[2],upper[2],ierflag);
1900     clusterInput.Fitter()->mnparm(3,"y2",vstart[3],step[3],lower[3],upper[3],ierflag);
1901     clusterInput.Fitter()->mnparm(4,"a0",vstart[4],step[4],lower[4],upper[4],ierflag);
1902     clusterInput.Fitter()->mnparm(5,"a1",vstart[5],step[5],lower[5],upper[5],ierflag);
1903 // ready for minimisation
1904     arglist[0]= -1;
1905     arglist[1]= 0;
1906
1907     clusterInput.Fitter()->mnexcm("SET NOGR", arglist, 0, ierflag);
1908     clusterInput.Fitter()->mnexcm("MIGRAD", arglist, 0, ierflag);
1909     //    clusterInput.Fitter()->mnexcm("EXIT" , arglist, 0, ierflag);
1910 // Get fitted parameters
1911     TString chname;
1912     Double_t epxz, b1, b2;
1913     Int_t ierflg;
1914     clusterInput.Fitter()->mnpout(0, chname, fXFit[0],  epxz, b1, b2, ierflg);
1915     clusterInput.Fitter()->mnpout(1, chname, fYFit[0],  epxz, b1, b2, ierflg);
1916     clusterInput.Fitter()->mnpout(2, chname, fXFit[1],  epxz, b1, b2, ierflg);
1917     clusterInput.Fitter()->mnpout(3, chname, fYFit[1],  epxz, b1, b2, ierflg);
1918     clusterInput.Fitter()->mnpout(4, chname, fQrFit[0], epxz, b1, b2, ierflg);
1919     clusterInput.Fitter()->mnpout(5, chname, fQrFit[1], epxz, b1, b2, ierflg);
1920
1921     Double_t fmin, fedm, errdef;
1922     Int_t   npari, nparx, istat;
1923
1924     clusterInput.Fitter()->mnstat(fmin, fedm, errdef, npari, nparx, istat);
1925     fFitStat=istat;
1926
1927     fChi2[0]=fmin;
1928     fChi2[1]=fmin;
1929     return fmin;
1930 }
1931
1932 void AliMUONClusterFinderVS::Split(AliMUONRawCluster* c)
1933 {
1934 //
1935 // One cluster for each maximum
1936 //
1937     Int_t i, j, cath;
1938     AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
1939     for (j=0; j<2; j++) {
1940         AliMUONRawCluster cnew;
1941         cnew.SetGhost(c->GetGhost());
1942         for (cath=0; cath<2; cath++) {
1943             cnew.SetChi2(cath,fChi2[0]);
1944             // ?? why not cnew.fChi2[cath]=fChi2[cath];
1945
1946             if (fNPeaks == 0) {
1947                 cnew.SetNcluster(0,-1);
1948                 cnew.SetNcluster(1,fNRawClusters);
1949             } else {
1950                 cnew.SetNcluster(0,fNPeaks);
1951                 cnew.SetNcluster(1,0);
1952             }
1953             cnew.SetMultiplicity(cath,0);
1954             cnew.SetX(cath, Float_t(fXFit[j]));
1955             cnew.SetY(cath, Float_t(fYFit[j]));
1956             cnew.SetZ(cath, fZPlane);
1957             if (j==0) {
1958                 cnew.SetCharge(cath, Int_t(clusterInput.TotalCharge(cath)*fQrFit[cath]));
1959             } else {
1960                 cnew.SetCharge(cath, Int_t(clusterInput.TotalCharge(cath)*(1-fQrFit[cath])));
1961             }
1962             fSeg[cath]->SetHit(fXFit[j],fYFit[j],fZPlane);
1963             for (i=0; i<fMul[cath]; i++) {
1964                 cnew.SetIndex(cnew.GetMultiplicity(cath), cath, c->GetIndex(i,cath));
1965                 fSeg[cath]->SetPad(fIx[i][cath], fIy[i][cath]);
1966                 Float_t q1=fInput->Response()->IntXY(fSeg[cath]);
1967                 cnew.SetContrib(i, cath, q1*Float_t(cnew.GetCharge(cath))/Float_t(fQ[i][cath]));
1968                 cnew.SetMultiplicity(cath, cnew.GetMultiplicity(cath)+1 );
1969             }
1970             FillCluster(&cnew,0,cath);
1971         } // cathode loop
1972
1973         cnew.SetClusterType(cnew.PhysicsContribution());
1974         if (cnew.GetCharge(0)>0 && cnew.GetCharge(1)>0) AddRawCluster(cnew);
1975         fNPeaks++;
1976     }
1977 }
1978
1979
1980 //
1981 // Minimisation functions
1982 // Single Mathieson
1983 void fcnS1(Int_t & /*npar*/, Double_t * /*gin*/, Double_t &f, Double_t *par, Int_t /*iflag*/)
1984 {
1985     AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
1986     Int_t i;
1987     Float_t delta;
1988     Float_t chisq=0;
1989     Float_t qcont=0;
1990     Float_t qtot=0;
1991
1992     for (i=0; i<clusterInput.Nmul(0); i++) {
1993         Float_t q0=clusterInput.Charge(i,0);
1994         Float_t q1=clusterInput.DiscrChargeS1(i,par);
1995         delta=(q0-q1)/q0;
1996         chisq+=delta*delta;
1997         qcont+=q1;
1998         qtot+=q0;
1999     }
2000     f=chisq;
2001 }
2002
2003 void fcnCombiS1(Int_t & /*npar*/, Double_t * /*gin*/, Double_t &f, Double_t *par, Int_t /*iflag*/)
2004 {
2005     AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
2006     Int_t i, cath;
2007     Float_t delta;
2008     Float_t chisq=0;
2009     Float_t qcont=0;
2010     Float_t qtot=0;
2011
2012     for (cath=0; cath<2; cath++) {
2013         for (i=0; i<clusterInput.Nmul(cath); i++) {
2014             Float_t q0=clusterInput.Charge(i,cath);
2015             Float_t q1=clusterInput.DiscrChargeCombiS1(i,par,cath);
2016             delta=(q0-q1)/q0;
2017             chisq+=delta*delta;
2018             qcont+=q1;
2019             qtot+=q0;
2020         }
2021     }
2022     f=chisq;
2023 }
2024
2025 // Double Mathieson
2026 void fcnS2(Int_t & /*npar*/, Double_t * /*gin*/, Double_t &f, Double_t *par, Int_t /*iflag*/)
2027 {
2028     AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
2029     Int_t i;
2030     Float_t delta;
2031     Float_t chisq=0;
2032     Float_t qcont=0;
2033     Float_t qtot=0;
2034
2035     for (i=0; i<clusterInput.Nmul(0); i++) {
2036
2037         Float_t q0=clusterInput.Charge(i,0);
2038         Float_t q1=clusterInput.DiscrChargeS2(i,par);
2039         delta=(q0-q1)/q0;
2040         chisq+=delta*delta;
2041         qcont+=q1;
2042         qtot+=q0;
2043     }
2044     f=chisq;
2045 }
2046
2047 // Double Mathieson
2048 void fcnCombiS2(Int_t & /*npar*/, Double_t * /*gin*/, Double_t &f, Double_t *par, Int_t /*iflag*/)
2049 {
2050     AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
2051     Int_t i, cath;
2052     Float_t delta;
2053     Float_t chisq=0;
2054     Float_t qcont=0;
2055     Float_t qtot=0;
2056     for (cath=0; cath<2; cath++) {
2057         for (i=0; i<clusterInput.Nmul(cath); i++) {
2058             Float_t q0=clusterInput.Charge(i,cath);
2059             Float_t q1=clusterInput.DiscrChargeCombiS2(i,par,cath);
2060             delta=(q0-q1)/q0;
2061             chisq+=delta*delta;
2062             qcont+=q1;
2063             qtot+=q0;
2064         }
2065     }
2066     f=chisq;
2067 }
2068
2069 void AliMUONClusterFinderVS::AddRawCluster(const AliMUONRawCluster& c)
2070 {
2071   //
2072   // Add a raw cluster copy to the list
2073   //
2074
2075 //     AliMUON *pMUON=(AliMUON*)gAlice->GetModule("MUON");
2076 //     pMUON->GetMUONData()->AddRawCluster(fInput->Chamber(),c);
2077 //     fNRawClusters++;
2078
2079
2080     TClonesArray &lrawcl = *fRawClusters;
2081     new(lrawcl[fNRawClusters++]) AliMUONRawCluster(c);
2082     if (fDebugLevel)
2083         fprintf(stderr,"\nfNRawClusters %d\n",fNRawClusters);
2084 }
2085
2086 Bool_t AliMUONClusterFinderVS::TestTrack(Int_t t) const {
2087 // Test if track was user selected
2088     if (fTrack[0]==-1 || fTrack[1]==-1) {
2089         return kTRUE;
2090     } else if (t==fTrack[0] || t==fTrack[1]) {
2091         return kTRUE;
2092     } else {
2093         return kFALSE;
2094     }
2095 }
2096
2097 AliMUONClusterFinderVS& AliMUONClusterFinderVS
2098 ::operator = (const AliMUONClusterFinderVS& rhs)
2099 {
2100 // Protected assignement operator
2101
2102   if (this == &rhs) return *this;
2103
2104   Fatal("operator=", "Not implemented.");
2105
2106   return *this;
2107 }
2108
2109
2110
2111
2112
2113
2114
2115
2116