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 #include "AliLog.h"
  29
  30 //_____________________________________________________________________
  31 // This function is minimized in the double-Mathieson fit
  32 void fcnS2(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
  33 void fcnS1(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
  34 void fcnCombiS1(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
  35 void fcnCombiS2(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
  36
  37 ClassImp(AliMUONClusterFinderVS)
  38
  39 AliMUONClusterFinderVS::AliMUONClusterFinderVS()
  40   : TObject()
  41 {
  42 // Default constructor
  43     fInput=AliMUONClusterInput::Instance();
  44     fHitMap[0] = 0;
  45     fHitMap[1] = 0;
  46     fTrack[0]=fTrack[1]=-1;
  47     fDebugLevel = 0; // make silent default
  48     fGhostChi2Cut = 1e6; // nothing done by default
  49     fSeg[0]    = 0;
  50     fSeg[1]    = 0;
  51     for(Int_t i=0; i<100; i++) {
  52       for (Int_t j=0; j<2; j++) {
  53         fDig[i][j] = 0;
  54       }
  55     }
  56     fRawClusters = new TClonesArray("AliMUONRawCluster",1000);
  57     fNRawClusters = 0;
  58 }
  59  //____________________________________________________________________________
  60 AliMUONClusterFinderVS::~AliMUONClusterFinderVS()
  61 {
  62   // Reset tracks information
  63    fNRawClusters = 0;
  64    if (fRawClusters) {
  65      fRawClusters->Delete();
  66      delete fRawClusters;
  67    }
  68 }
  69
  70 AliMUONClusterFinderVS::AliMUONClusterFinderVS(const AliMUONClusterFinderVS & clusterFinder):TObject(clusterFinder)
  71 {
  72 // Protected copy constructor
  73
  74   AliFatal("Not implemented.");
  75 }
  76 //____________________________________________________________________________
  77 void AliMUONClusterFinderVS::ResetRawClusters()
  78 {
  79   // Reset tracks information
  80   fNRawClusters = 0;
  81   if (fRawClusters) fRawClusters->Clear();
  82 }
  83 //____________________________________________________________________________
  84 void AliMUONClusterFinderVS::Decluster(AliMUONRawCluster *cluster)
  85 {
  86 // Decluster by local maxima
  87     SplitByLocalMaxima(cluster);
  88 }
  89 //____________________________________________________________________________
  90 void AliMUONClusterFinderVS::SplitByLocalMaxima(AliMUONRawCluster *c)
  91 {
  92 // Split complex cluster by local maxima
  93     Int_t cath, i;
  94
  95     fInput->SetCluster(c);
  96
  97     fMul[0]=c->GetMultiplicity(0);
  98     fMul[1]=c->GetMultiplicity(1);
  99
 100 //
 101 //  dump digit information into arrays
 102 //
 103
 104     Float_t qtot;
 105
 106     for (cath=0; cath<2; cath++) {
 107         qtot=0;
 108         for (i=0; i<fMul[cath]; i++)
 109         {
 110             // pointer to digit
 111             fDig[i][cath]=fInput->Digit(cath, c->GetIndex(i, cath));
 112             // pad coordinates
 113             fIx[i][cath]= fDig[i][cath]->PadX();
 114             fIy[i][cath]= fDig[i][cath]->PadY();
 115             // pad charge
 116             fQ[i][cath] = fDig[i][cath]->Signal();
 117             // pad centre coordinates
 118             fSeg[cath]->
 119                 GetPadC(fIx[i][cath], fIy[i][cath], fX[i][cath], fY[i][cath], fZ[i][cath]);
 120         } // loop over cluster digits
 121     }  // loop over cathodes
 122
 123
 124     FindLocalMaxima(c);
 125
 126 //
 127 //  Initialise and perform mathieson fits
 128     Float_t chi2, oldchi2;
 129 //  ++++++++++++++++++*************+++++++++++++++++++++
 130 //  (1) No more than one local maximum per cathode plane
 131 //  +++++++++++++++++++++++++++++++*************++++++++
 132     if ((fNLocal[0]==1 && (fNLocal[1]==0 ||  fNLocal[1]==1)) ||
 133         (fNLocal[0]==0 && fNLocal[1]==1)) {
 134 // Perform combined single Mathieson fit
 135 // Initial values for coordinates (x,y)
 136
 137         // One local maximum on cathodes 1 and 2 (X->cathode 2, Y->cathode 1)
 138         if (fNLocal[0]==1 &&  fNLocal[1]==1) {
 139             fXInit[0]=c->GetX(1);
 140             fYInit[0]=c->GetY(0);
 141             // One local maximum on cathode 1 (X,Y->cathode 1)
 142         } else if (fNLocal[0]==1) {
 143             fXInit[0]=c->GetX(0);
 144             fYInit[0]=c->GetY(0);
 145             // One local maximum on cathode 2  (X,Y->cathode 2)
 146         } else {
 147             fXInit[0]=c->GetX(1);
 148             fYInit[0]=c->GetY(1);
 149         }
 150         AliDebug(1,"cas (1) CombiSingleMathiesonFit(c)");
 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         AliDebug(1,Form(" chi2 %f ",chi2));
 158
 159         c->SetX(0, fXFit[0]);
 160         c->SetY(0, fYFit[0]);
 161
 162         c->SetX(1,fXFit[0]);
 163         c->SetY(1,fYFit[0]);
 164         c->SetChi2(0,chi2);
 165         c->SetChi2(1,chi2);
 166         // Force on anod
 167         c->SetX(0, fSeg[0]->GetAnod(c->GetX(0)));
 168         c->SetX(1, fSeg[1]->GetAnod(c->GetX(1)));
 169
 170 // If reasonable chi^2 add result to the list of rawclusters
 171         if (chi2 < 0.3) {
 172             AddRawCluster(*c);
 173 // If not try combined double Mathieson Fit
 174         } else {
 175           if (fDebugLevel)
 176             fprintf(stderr," MAUVAIS CHI2 !!!\n");
 177             if (fNLocal[0]==1 &&  fNLocal[1]==1) {
 178                 fXInit[0]=fX[fIndLocal[0][1]][1];
 179                 fYInit[0]=fY[fIndLocal[0][0]][0];
 180                 fXInit[1]=fX[fIndLocal[0][1]][1];
 181                 fYInit[1]=fY[fIndLocal[0][0]][0];
 182             } else if (fNLocal[0]==1) {
 183                 fXInit[0]=fX[fIndLocal[0][0]][0];
 184                 fYInit[0]=fY[fIndLocal[0][0]][0];
 185                 fXInit[1]=fX[fIndLocal[0][0]][0];
 186                 fYInit[1]=fY[fIndLocal[0][0]][0];
 187             } else {
 188                 fXInit[0]=fX[fIndLocal[0][1]][1];
 189                 fYInit[0]=fY[fIndLocal[0][1]][1];
 190                 fXInit[1]=fX[fIndLocal[0][1]][1];
 191                 fYInit[1]=fY[fIndLocal[0][1]][1];
 192             }
 193
 194 //  Initial value for charge ratios
 195             fQrInit[0]=0.5;
 196             fQrInit[1]=0.5;
 197             if (fDebugLevel)
 198             fprintf(stderr,"\n cas (1) CombiDoubleMathiesonFit(c)\n");
 199             chi2=CombiDoubleMathiesonFit(c);
 200 //          Int_t ndf = fgNbins[0]+fgNbins[1]-6;
 201 //          Float_t prob = TMath::Prob(chi2,ndf);
 202 //          prob2->Fill(prob);
 203 //          chi2_2->Fill(chi2);
 204
 205 // Was this any better ??
 206             AliDebug(1,Form(" Old and new chi2 %f %f ", oldchi2, chi2));
 207             if (fFitStat!=0 && chi2>0 && (2.*chi2 < oldchi2)) {
 208               AliDebug(1,"Split");
 209                 // Split cluster into two according to fit result
 210                 Split(c);
 211             } else {
 212               AliDebug(1,"Do not Split");
 213                 // Don't split
 214                 AddRawCluster(*c);
 215             }
 216         }
 217
 218 //  +++++++++++++++++++++++++++++++++++++++
 219 //  (2) Two local maxima per cathode plane
 220 //  +++++++++++++++++++++++++++++++++++++++
 221     } else if (fNLocal[0]==2 &&  fNLocal[1]==2) {
 222 //
 223 //  Let's look for ghosts first
 224
 225         Float_t xm[4][2], ym[4][2];
 226         Float_t dpx, dpy, dx, dy;
 227         Int_t ixm[4][2], iym[4][2];
 228         Int_t isec, im1, im2, ico;
 229 //
 230 //  Form the 2x2 combinations
 231 //  0-0, 0-1, 1-0, 1-1
 232         ico=0;
 233         for (im1=0; im1<2; im1++) {
 234             for (im2=0; im2<2; im2++) {
 235                 xm[ico][0]=fX[fIndLocal[im1][0]][0];
 236                 ym[ico][0]=fY[fIndLocal[im1][0]][0];
 237                 xm[ico][1]=fX[fIndLocal[im2][1]][1];
 238                 ym[ico][1]=fY[fIndLocal[im2][1]][1];
 239
 240                 ixm[ico][0]=fIx[fIndLocal[im1][0]][0];
 241                 iym[ico][0]=fIy[fIndLocal[im1][0]][0];
 242                 ixm[ico][1]=fIx[fIndLocal[im2][1]][1];
 243                 iym[ico][1]=fIy[fIndLocal[im2][1]][1];
 244                 ico++;
 245             }
 246         }
 247 // ico = 0 : first local maximum on cathodes 1 and 2
 248 // ico = 1 : fisrt local maximum on cathode 1 and second on cathode 2
 249 // ico = 2 : second local maximum on cathode 1 and first on cathode 1
 250 // ico = 3 : second local maximum on cathodes 1 and 2
 251
 252 // Analyse the combinations and keep those that are possible !
 253 // For each combination check consistency in x and y
 254         Int_t   iacc;
 255         Bool_t  accepted[4];
 256         Float_t dr[4] = {1.e4, 1.e4, 1.e4, 1.e4};
 257         iacc=0;
 258
 259 // In case of staggering maxima are displaced by exactly half the pad-size in y.
 260 // We have to take into account the numerical precision in the consistency check;
 261         Float_t eps = 1.e-5;
 262 //
 263         for (ico=0; ico<4; ico++) {
 264             accepted[ico]=kFALSE;
 265 // cathode one: x-coordinate
 266             isec=fSeg[0]->Sector(ixm[ico][0], iym[ico][0]);
 267             dpx=fSeg[0]->Dpx(isec)/2.;
 268             dx=TMath::Abs(xm[ico][0]-xm[ico][1]);
 269 // cathode two: y-coordinate
 270             isec=fSeg[1]->Sector(ixm[ico][1], iym[ico][1]);
 271             dpy=fSeg[1]->Dpy(isec)/2.;
 272             dy=TMath::Abs(ym[ico][0]-ym[ico][1]);
 273             if (fDebugLevel>1)
 274                 printf("\n %i %f %f %f %f %f %f \n", ico, ym[ico][0], ym[ico][1], dy, dpy, dx, dpx );
 275             if ((dx <= dpx) && (dy <= dpy+eps)) {
 276                 // consistent
 277                 accepted[ico]=kTRUE;
 278                 dr[ico] = TMath::Sqrt(dx*dx+dy*dy);
 279                 iacc++;
 280             } else {
 281                 // reject
 282                 accepted[ico]=kFALSE;
 283             }
 284         }
 285         if (fDebugLevel)
 286           printf("\n iacc= %d:\n", iacc);
 287         if (iacc == 3) {
 288             if (accepted[0] && accepted[1]) {
 289                 if (dr[0] >= dr[1]) {
 290                     accepted[0]=kFALSE;
 291                 } else {
 292                     accepted[1]=kFALSE;
 293                 }
 294             }
 295
 296             if (accepted[2] && accepted[3]) {
 297                 if (dr[2] >= dr[3]) {
 298                     accepted[2]=kFALSE;
 299                 } else {
 300                     accepted[3]=kFALSE;
 301                 }
 302             }
 303 /*
 304 // eliminate one candidate
 305             Float_t drmax = 0;
 306             Int_t icobad = -1;
 307
 308             for (ico=0; ico<4; ico++) {
 309                 if (accepted[ico] && dr[ico] > drmax) {
 310                     icobad = ico;
 311                     drmax  = dr[ico];
 312                 }
 313             }
 314
 315             accepted[icobad] = kFALSE;
 316 */
 317             iacc = 2;
 318         }
 319
 320
 321         if (fDebugLevel) {
 322           printf("\n iacc= %d:\n", iacc);
 323             if (iacc==2) {
 324                 fprintf(stderr,"\n iacc=2: No problem ! \n");
 325             } else if (iacc==4) {
 326                 fprintf(stderr,"\n iacc=4: Ok, but ghost problem !!! \n");
 327             } else if (iacc==0) {
 328                 fprintf(stderr,"\n iacc=0: I don't know what to do with this !!!!!!!!! \n");
 329             }
 330         }
 331
 332 //  Initial value for charge ratios
 333         fQrInit[0]=Float_t(fQ[fIndLocal[0][0]][0])/
 334             Float_t(fQ[fIndLocal[0][0]][0]+fQ[fIndLocal[1][0]][0]);
 335         fQrInit[1]=Float_t(fQ[fIndLocal[0][1]][1])/
 336             Float_t(fQ[fIndLocal[0][1]][1]+fQ[fIndLocal[1][1]][1]);
 337
 338 // ******* iacc = 0 *******
 339 // No combinations found between the 2 cathodes
 340 // We keep the center of gravity of the cluster
 341         if (iacc==0) {
 342             AddRawCluster(*c);
 343         }
 344
 345 // ******* iacc = 1 *******
 346 // Only one combination found between the 2 cathodes
 347         if (iacc==1) {
 348 // Initial values for the 2 maxima (x,y)
 349
 350 // 1 maximum is initialised with the maximum of the combination found (X->cathode 2, Y->cathode 1)
 351 // 1 maximum is initialised with the other maximum of the first cathode
 352             if (accepted[0]){
 353                 AliDebug(1,"ico=0");
 354                 fXInit[0]=xm[0][1];
 355                 fYInit[0]=ym[0][0];
 356                 fXInit[1]=xm[3][0];
 357                 fYInit[1]=ym[3][0];
 358             } else if (accepted[1]){
 359                 AliDebug(1,"ico=1");
 360                 fXInit[0]=xm[1][1];
 361                 fYInit[0]=ym[1][0];
 362                 fXInit[1]=xm[2][0];
 363                 fYInit[1]=ym[2][0];
 364             } else if (accepted[2]){
 365                 AliDebug(1,"ico=2");
 366                 fXInit[0]=xm[2][1];
 367                 fYInit[0]=ym[2][0];
 368                 fXInit[1]=xm[1][0];
 369                 fYInit[1]=ym[1][0];
 370             } else if (accepted[3]){
 371                 AliDebug(1,"ico=3");
 372                 fXInit[0]=xm[3][1];
 373                 fYInit[0]=ym[3][0];
 374                 fXInit[1]=xm[0][0];
 375                 fYInit[1]=ym[0][0];
 376             }
 377                 AliDebug(1,"cas (2) CombiDoubleMathiesonFit(c)");
 378             chi2=CombiDoubleMathiesonFit(c);
 379 //          Int_t ndf = fgNbins[0]+fgNbins[1]-6;
 380 //          Float_t prob = TMath::Prob(chi2,ndf);
 381 //          prob2->Fill(prob);
 382 //          chi2_2->Fill(chi2);
 383             AliDebug(1,Form(" chi2 %f\n",chi2));
 384
 385 // If reasonable chi^2 add result to the list of rawclusters
 386             if (chi2<10) {
 387                 Split(c);
 388
 389             } else {
 390 // 1 maximum is initialised with the maximum of the combination found (X->cathode 2, Y->cathode 1)
 391 // 1 maximum is initialised with the other maximum of the second cathode
 392                 if (accepted[0]){
 393                         AliDebug(1,"ico=0");
 394                     fXInit[0]=xm[0][1];
 395                     fYInit[0]=ym[0][0];
 396                     fXInit[1]=xm[3][1];
 397                     fYInit[1]=ym[3][1];
 398                 } else if (accepted[1]){
 399                         AliDebug(1,"ico=1");
 400                     fXInit[0]=xm[1][1];
 401                     fYInit[0]=ym[1][0];
 402                     fXInit[1]=xm[2][1];
 403                     fYInit[1]=ym[2][1];
 404                 } else if (accepted[2]){
 405                         AliDebug(1,"ico=2");
 406                     fXInit[0]=xm[2][1];
 407                     fYInit[0]=ym[2][0];
 408                     fXInit[1]=xm[1][1];
 409                     fYInit[1]=ym[1][1];
 410                 } else if (accepted[3]){
 411                         AliDebug(1,"ico=3");
 412                     fXInit[0]=xm[3][1];
 413                     fYInit[0]=ym[3][0];
 414                     fXInit[1]=xm[0][1];
 415                     fYInit[1]=ym[0][1];
 416                 }
 417                 AliDebug(1,"\n cas (2) CombiDoubleMathiesonFit(c)\n");
 418                 chi2=CombiDoubleMathiesonFit(c);
 419 //              Int_t ndf = fgNbins[0]+fgNbins[1]-6;
 420 //              Float_t prob = TMath::Prob(chi2,ndf);
 421 //              prob2->Fill(prob);
 422 //              chi2_2->Fill(chi2);
 423                 if (fDebugLevel)
 424                     fprintf(stderr," chi2 %f\n",chi2);
 425
 426 // If reasonable chi^2 add result to the list of rawclusters
 427                 if (chi2<10) {
 428                     Split(c);
 429                 } else {
 430 //We keep only the combination found (X->cathode 2, Y->cathode 1)
 431                     for (Int_t ico=0; ico<2; ico++) {
 432                         if (accepted[ico]) {
 433                             AliMUONRawCluster cnew;
 434                             Int_t cath;
 435                             for (cath=0; cath<2; cath++) {
 436                                 cnew.SetX(cath, Float_t(xm[ico][1]));
 437                                 cnew.SetY(cath, Float_t(ym[ico][0]));
 438                                 cnew.SetZ(cath, fZPlane);
 439
 440                                 cnew.SetMultiplicity(cath,c->GetMultiplicity(cath));
 441                                 for (i=0; i<fMul[cath]; i++) {
 442                                     cnew.SetIndex(i, cath, c->GetIndex(i,cath));
 443                                     fSeg[cath]->SetPad(fIx[i][cath], fIy[i][cath]);
 444                                 }
 445                                 if (fDebugLevel) {
 446                                   fprintf(stderr,"\nRawCluster %d cath %d\n",ico,cath);
 447                                   fprintf(stderr,"mult_av %d\n",c->GetMultiplicity(cath));
 448                                 }
 449                                 FillCluster(&cnew,cath);
 450                             }
 451                             cnew.SetClusterType(cnew.PhysicsContribution());
 452                             AddRawCluster(cnew);
 453                             fNPeaks++;
 454                         }
 455                     }
 456                 }
 457             }
 458         }
 459
 460 // ******* iacc = 2 *******
 461 // Two combinations found between the 2 cathodes
 462         if (iacc==2) {
 463 // Was the same maximum taken twice
 464             if ((accepted[0]&&accepted[1]) || (accepted[2]&&accepted[3])) {
 465                 if (fDebugLevel) fprintf(stderr,"\n Maximum taken twice !!!\n");
 466
 467 // Have a try !! with that
 468                 if (accepted[0]&&accepted[3]) {
 469                     fXInit[0]=xm[0][1];
 470                     fYInit[0]=ym[0][0];
 471                     fXInit[1]=xm[1][1];
 472                     fYInit[1]=ym[1][0];
 473                 } else {
 474                     fXInit[0]=xm[2][1];
 475                     fYInit[0]=ym[2][0];
 476                     fXInit[1]=xm[3][1];
 477                     fYInit[1]=ym[3][0];
 478                 }
 479                 if (fDebugLevel)
 480                     fprintf(stderr,"\n cas (2) CombiDoubleMathiesonFit(c)\n");
 481                 chi2=CombiDoubleMathiesonFit(c);
 482 //                  Int_t ndf = fgNbins[0]+fgNbins[1]-6;
 483 //                  Float_t prob = TMath::Prob(chi2,ndf);
 484 //                  prob2->Fill(prob);
 485 //                  chi2_2->Fill(chi2);
 486                 Split(c);
 487
 488             } else {
 489 // No ghosts ! No Problems ! -  Perform one fit only !
 490                 if (accepted[0]&&accepted[3]) {
 491                     fXInit[0]=xm[0][1];
 492                     fYInit[0]=ym[0][0];
 493                     fXInit[1]=xm[3][1];
 494                     fYInit[1]=ym[3][0];
 495                 } else {
 496                     fXInit[0]=xm[1][1];
 497                     fYInit[0]=ym[1][0];
 498                     fXInit[1]=xm[2][1];
 499                     fYInit[1]=ym[2][0];
 500                 }
 501                 if (fDebugLevel)
 502                     fprintf(stderr,"\n cas (2) CombiDoubleMathiesonFit(c)\n");
 503                 chi2=CombiDoubleMathiesonFit(c);
 504 //                  Int_t ndf = fgNbins[0]+fgNbins[1]-6;
 505 //                  Float_t prob = TMath::Prob(chi2,ndf);
 506 //                  prob2->Fill(prob);
 507 //                  chi2_2->Fill(chi2);
 508                 if (fDebugLevel)
 509                     fprintf(stderr," chi2 %f\n",chi2);
 510                 Split(c);
 511             }
 512
 513 // ******* iacc = 4 *******
 514 // Four combinations found between the 2 cathodes
 515 // Ghost !!
 516         } else if (iacc==4) {
 517 // Perform fits for the two possibilities !!
 518 // Accept if charges are compatible on both cathodes
 519 // If none are compatible, keep everything
 520             fXInit[0]=xm[0][1];
 521             fYInit[0]=ym[0][0];
 522             fXInit[1]=xm[3][1];
 523             fYInit[1]=ym[3][0];
 524             if (fDebugLevel)
 525                 fprintf(stderr,"\n cas (2) CombiDoubleMathiesonFit(c)\n");
 526             chi2=CombiDoubleMathiesonFit(c);
 527 //              Int_t ndf = fgNbins[0]+fgNbins[1]-6;
 528 //              Float_t prob = TMath::Prob(chi2,ndf);
 529 //              prob2->Fill(prob);
 530 //              chi2_2->Fill(chi2);
 531             if (fDebugLevel)
 532                 fprintf(stderr," chi2 %f\n",chi2);
 533             // store results of fit and postpone decision
 534             Double_t sXFit[2],sYFit[2],sQrFit[2];
 535             Float_t sChi2[2];
 536             for (Int_t i=0;i<2;i++) {
 537                 sXFit[i]=fXFit[i];
 538                 sYFit[i]=fYFit[i];
 539                 sQrFit[i]=fQrFit[i];
 540                 sChi2[i]=fChi2[i];
 541             }
 542             fXInit[0]=xm[1][1];
 543             fYInit[0]=ym[1][0];
 544             fXInit[1]=xm[2][1];
 545             fYInit[1]=ym[2][0];
 546             if (fDebugLevel)
 547                 fprintf(stderr,"\n cas (2) CombiDoubleMathiesonFit(c)\n");
 548             chi2=CombiDoubleMathiesonFit(c);
 549 //              ndf = fgNbins[0]+fgNbins[1]-6;
 550 //              prob = TMath::Prob(chi2,ndf);
 551 //              prob2->Fill(prob);
 552 //              chi2_2->Fill(chi2);
 553             if (fDebugLevel)
 554                 fprintf(stderr," chi2 %f\n",chi2);
 555             // We have all informations to perform the decision
 556             // Compute the chi2 for the 2 possibilities
 557             Float_t chi2fi,chi2si,chi2f,chi2s;
 558
 559             chi2f = (TMath::Log(fInput->TotalCharge(0)*fQrFit[0]
 560                   /  (fInput->TotalCharge(1)*fQrFit[1]) )
 561                   / fInput->Response()->ChargeCorrel() );
 562             chi2f *=chi2f;
 563             chi2fi = (TMath::Log(fInput->TotalCharge(0)*(1-fQrFit[0])
 564                   /  (fInput->TotalCharge(1)*(1-fQrFit[1])) )
 565                   / fInput->Response()->ChargeCorrel() );
 566             chi2f += chi2fi*chi2fi;
 567
 568             chi2s = (TMath::Log(fInput->TotalCharge(0)*sQrFit[0]
 569                   /  (fInput->TotalCharge(1)*sQrFit[1]) )
 570                   / fInput->Response()->ChargeCorrel() );
 571             chi2s *=chi2s;
 572             chi2si = (TMath::Log(fInput->TotalCharge(0)*(1-sQrFit[0])
 573                   /  (fInput->TotalCharge(1)*(1-sQrFit[1])) )
 574                   / fInput->Response()->ChargeCorrel() );
 575             chi2s += chi2si*chi2si;
 576
 577             // usefull to store the charge matching chi2 in the cluster
 578             // fChi2[0]=sChi2[1]=chi2f;
 579             // fChi2[1]=sChi2[0]=chi2s;
 580
 581             if (chi2f<=fGhostChi2Cut && chi2s<=fGhostChi2Cut)
 582                 c->SetGhost(1);
 583             if   (chi2f>fGhostChi2Cut && chi2s>fGhostChi2Cut) {
 584                 // we keep the ghost
 585                 c->SetGhost(2);
 586                 chi2s=-1;
 587                 chi2f=-1;
 588             }
 589             if (chi2f<=fGhostChi2Cut)
 590                 Split(c);
 591             if (chi2s<=fGhostChi2Cut) {
 592                 // retreive saved values
 593                 for (Int_t i=0;i<2;i++) {
 594                     fXFit[i]=sXFit[i];
 595                     fYFit[i]=sYFit[i];
 596                     fQrFit[i]=sQrFit[i];
 597                     fChi2[i]=sChi2[i];
 598                 }
 599                 Split(c);
 600             }
 601             c->SetGhost(0);
 602         }
 603
 604     } else if (fNLocal[0]==2 &&  fNLocal[1]==1) {
 605 //  ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 606 //  (3) Two local maxima on cathode 1 and one maximum on cathode 2
 607 //  ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 608 //
 609         Float_t xm[4][2], ym[4][2];
 610         Float_t dpx, dpy, dx, dy;
 611         Int_t ixm[4][2], iym[4][2];
 612         Int_t isec, im1, ico;
 613 //
 614 //  Form the 2x2 combinations
 615 //  0-0, 0-1, 1-0, 1-1
 616         ico=0;
 617         for (im1=0; im1<2; im1++) {
 618             xm[ico][0]=fX[fIndLocal[im1][0]][0];
 619             ym[ico][0]=fY[fIndLocal[im1][0]][0];
 620             xm[ico][1]=fX[fIndLocal[0][1]][1];
 621             ym[ico][1]=fY[fIndLocal[0][1]][1];
 622
 623             ixm[ico][0]=fIx[fIndLocal[im1][0]][0];
 624             iym[ico][0]=fIy[fIndLocal[im1][0]][0];
 625             ixm[ico][1]=fIx[fIndLocal[0][1]][1];
 626             iym[ico][1]=fIy[fIndLocal[0][1]][1];
 627             ico++;
 628         }
 629 // ico = 0 : first local maximum on cathodes 1 and 2
 630 // ico = 1 : second local maximum on cathode 1 and first on cathode 2
 631
 632 // Analyse the combinations and keep those that are possible !
 633 // For each combination check consistency in x and y
 634         Int_t iacc;
 635         Bool_t accepted[4];
 636         iacc=0;
 637         // In case of staggering maxima are displaced by exactly half the pad-size in y.
 638         // We have to take into account the numerical precision in the consistency check;
 639
 640         Float_t eps = 1.e-5;
 641
 642         for (ico=0; ico<2; ico++) {
 643             accepted[ico]=kFALSE;
 644             isec=fSeg[0]->Sector(ixm[ico][0], iym[ico][0]);
 645             dpx=fSeg[0]->Dpx(isec)/2.;
 646             dx=TMath::Abs(xm[ico][0]-xm[ico][1]);
 647             isec=fSeg[1]->Sector(ixm[ico][1], iym[ico][1]);
 648             dpy=fSeg[1]->Dpy(isec)/2.;
 649             dy=TMath::Abs(ym[ico][0]-ym[ico][1]);
 650             if (fDebugLevel>1)
 651                 printf("\n %i %f %f %f %f \n", ico, ym[ico][0], ym[ico][1], dy, dpy );
 652             if ((dx <= dpx) && (dy <= dpy+eps)) {
 653                 // consistent
 654                 accepted[ico]=kTRUE;
 655                 iacc++;
 656             } else {
 657                 // reject
 658                 accepted[ico]=kFALSE;
 659             }
 660         }
 661
 662         Float_t chi21 = 100;
 663         Float_t chi22 = 100;
 664         Float_t chi23 = 100;
 665
 666         //  Initial value for charge ratios
 667         fQrInit[0]=Float_t(fQ[fIndLocal[0][0]][0])/
 668             Float_t(fQ[fIndLocal[0][0]][0]+fQ[fIndLocal[1][0]][0]);
 669         fQrInit[1]=fQrInit[0];
 670
 671         if (accepted[0] && accepted[1]) {
 672
 673             fXInit[0]=0.5*(xm[0][1]+xm[0][0]);
 674             fYInit[0]=ym[0][0];
 675             fXInit[1]=0.5*(xm[0][1]+xm[1][0]);
 676             fYInit[1]=ym[1][0];
 677             fQrInit[0]=0.5;
 678             fQrInit[1]=0.5;
 679             chi23=CombiDoubleMathiesonFit(c);
 680             if (chi23<10) {
 681                 Split(c);
 682                 Float_t yst;
 683                 yst = fYFit[0];
 684                 fYFit[0] = fYFit[1];
 685                 fYFit[1] = yst;
 686                 Split(c);
 687             }
 688         } else if (accepted[0]) {
 689             fXInit[0]=xm[0][1];
 690             fYInit[0]=ym[0][0];
 691             fXInit[1]=xm[1][0];
 692             fYInit[1]=ym[1][0];
 693             chi21=CombiDoubleMathiesonFit(c);
 694 //          Int_t ndf = fgNbins[0]+fgNbins[1]-6;
 695 //          Float_t prob = TMath::Prob(chi2,ndf);
 696 //          prob2->Fill(prob);
 697 //          chi2_2->Fill(chi21);
 698             if (fDebugLevel)
 699                 fprintf(stderr," chi2 %f\n",chi21);
 700             if (chi21<10) Split(c);
 701         } else if (accepted[1]) {
 702             fXInit[0]=xm[1][1];
 703             fYInit[0]=ym[1][0];
 704             fXInit[1]=xm[0][0];
 705             fYInit[1]=ym[0][0];
 706             chi22=CombiDoubleMathiesonFit(c);
 707 //          Int_t ndf = fgNbins[0]+fgNbins[1]-6;
 708 //          Float_t prob = TMath::Prob(chi2,ndf);
 709 //          prob2->Fill(prob);
 710 //          chi2_2->Fill(chi22);
 711             if (fDebugLevel)
 712                 fprintf(stderr," chi2 %f\n",chi22);
 713             if (chi22<10) Split(c);
 714         }
 715
 716         if (chi21 > 10 && chi22 > 10 && chi23 > 10) {
 717 // We keep only the combination found (X->cathode 2, Y->cathode 1)
 718             for (Int_t ico=0; ico<2; ico++) {
 719                 if (accepted[ico]) {
 720                     AliMUONRawCluster cnew;
 721                     Int_t cath;
 722                     for (cath=0; cath<2; cath++) {
 723                         cnew.SetX(cath, Float_t(xm[ico][1]));
 724                         cnew.SetY(cath, Float_t(ym[ico][0]));
 725                         cnew.SetZ(cath, fZPlane);
 726                         cnew.SetMultiplicity(cath, c->GetMultiplicity(cath));
 727                         for (i=0; i<fMul[cath]; i++) {
 728                             cnew.SetIndex(i, cath, c->GetIndex(i, cath));
 729                             fSeg[cath]->SetPad(fIx[i][cath], fIy[i][cath]);
 730                         }
 731                         if (fDebugLevel){
 732                           fprintf(stderr,"\nRawCluster %d cath %d\n",ico,cath);
 733                           fprintf(stderr,"mult_av %d\n",c->GetMultiplicity(cath));
 734                         }
 735                         FillCluster(&cnew,cath);
 736                     }
 737                     cnew.SetClusterType(cnew.PhysicsContribution());
 738                     AddRawCluster(cnew);
 739                     fNPeaks++;
 740                 }
 741             }
 742         }
 743
 744 //  +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 745 //  (3') One local maximum on cathode 1 and two maxima on cathode 2
 746 //  +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 747     } else if (fNLocal[0]==1 && fNLocal[1]==2) {
 748         Float_t xm[4][2], ym[4][2];
 749         Float_t dpx, dpy, dx, dy;
 750         Int_t ixm[4][2], iym[4][2];
 751         Int_t isec, im1, ico;
 752 //
 753 //  Form the 2x2 combinations
 754 //  0-0, 0-1, 1-0, 1-1
 755         ico=0;
 756         for (im1=0; im1<2; im1++) {
 757             xm[ico][0]=fX[fIndLocal[0][0]][0];
 758             ym[ico][0]=fY[fIndLocal[0][0]][0];
 759             xm[ico][1]=fX[fIndLocal[im1][1]][1];
 760             ym[ico][1]=fY[fIndLocal[im1][1]][1];
 761
 762             ixm[ico][0]=fIx[fIndLocal[0][0]][0];
 763             iym[ico][0]=fIy[fIndLocal[0][0]][0];
 764             ixm[ico][1]=fIx[fIndLocal[im1][1]][1];
 765             iym[ico][1]=fIy[fIndLocal[im1][1]][1];
 766             ico++;
 767         }
 768 // ico = 0 : first local maximum on cathodes 1 and 2
 769 // ico = 1 : first local maximum on cathode 1 and second on cathode 2
 770
 771 // Analyse the combinations and keep those that are possible !
 772 // For each combination check consistency in x and y
 773         Int_t iacc;
 774         Bool_t accepted[4];
 775         iacc=0;
 776         // In case of staggering maxima are displaced by exactly half the pad-size in y.
 777         // We have to take into account the numerical precision in the consistency check;
 778         Float_t eps = 1.e-5;
 779
 780
 781         for (ico=0; ico<2; ico++) {
 782             accepted[ico]=kFALSE;
 783             isec=fSeg[0]->Sector(ixm[ico][0], iym[ico][0]);
 784             dpx=fSeg[0]->Dpx(isec)/2.;
 785             dx=TMath::Abs(xm[ico][0]-xm[ico][1]);
 786             isec=fSeg[1]->Sector(ixm[ico][1], iym[ico][1]);
 787             dpy=fSeg[1]->Dpy(isec)/2.;
 788             dy=TMath::Abs(ym[ico][0]-ym[ico][1]);
 789             if (fDebugLevel>0)
 790                 printf("\n %i %f %f %f %f \n", ico, ym[ico][0], ym[ico][1], dy, dpy );
 791             if ((dx <= dpx) && (dy <= dpy+eps)) {
 792                 // consistent
 793                 accepted[ico]=kTRUE;
 794                 if (fDebugLevel) fprintf(stderr,"ico %d\n",ico);
 795                 iacc++;
 796             } else {
 797                 // reject
 798                 accepted[ico]=kFALSE;
 799             }
 800         }
 801
 802         Float_t chi21 = 100;
 803         Float_t chi22 = 100;
 804         Float_t chi23 = 100;
 805
 806         fQrInit[1]=Float_t(fQ[fIndLocal[0][1]][1])/
 807             Float_t(fQ[fIndLocal[0][1]][1]+fQ[fIndLocal[1][1]][1]);
 808
 809         fQrInit[0]=fQrInit[1];
 810
 811
 812         if (accepted[0] && accepted[1]) {
 813             fXInit[0]=xm[0][1];
 814             fYInit[0]=0.5*(ym[0][0]+ym[0][1]);
 815             fXInit[1]=xm[1][1];
 816             fYInit[1]=0.5*(ym[0][0]+ym[1][1]);
 817             fQrInit[0]=0.5;
 818             fQrInit[1]=0.5;
 819             chi23=CombiDoubleMathiesonFit(c);
 820             if (chi23<10) {
 821                 Split(c);
 822                 Float_t yst;
 823                 yst = fYFit[0];
 824                 fYFit[0] = fYFit[1];
 825                 fYFit[1] = yst;
 826                 Split(c);
 827             }
 828         } else if (accepted[0]) {
 829             fXInit[0]=xm[0][0];
 830             fYInit[0]=ym[0][1];
 831             fXInit[1]=xm[1][1];
 832             fYInit[1]=ym[1][1];
 833             chi21=CombiDoubleMathiesonFit(c);
 834 //          Int_t ndf = fgNbins[0]+fgNbins[1]-6;
 835 //          Float_t prob = TMath::Prob(chi2,ndf);
 836 //          prob2->Fill(prob);
 837 //          chi2_2->Fill(chi21);
 838             if (fDebugLevel)
 839                 fprintf(stderr," chi2 %f\n",chi21);
 840             if (chi21<10) Split(c);
 841         } else if (accepted[1]) {
 842             fXInit[0]=xm[1][0];
 843             fYInit[0]=ym[1][1];
 844             fXInit[1]=xm[0][1];
 845             fYInit[1]=ym[0][1];
 846             chi22=CombiDoubleMathiesonFit(c);
 847 //          Int_t ndf = fgNbins[0]+fgNbins[1]-6;
 848 //          Float_t prob = TMath::Prob(chi2,ndf);
 849 //          prob2->Fill(prob);
 850 //          chi2_2->Fill(chi22);
 851             if (fDebugLevel)
 852                 fprintf(stderr," chi2 %f\n",chi22);
 853             if (chi22<10) Split(c);
 854         }
 855
 856         if (chi21 > 10 && chi22 > 10 && chi23 > 10) {
 857 //We keep only the combination found (X->cathode 2, Y->cathode 1)
 858             for (Int_t ico=0; ico<2; ico++) {
 859                 if (accepted[ico]) {
 860                     AliMUONRawCluster cnew;
 861                     Int_t cath;
 862                     for (cath=0; cath<2; cath++) {
 863                         cnew.SetX(cath, Float_t(xm[ico][1]));
 864                         cnew.SetY(cath, Float_t(ym[ico][0]));
 865                         cnew.SetZ(cath, fZPlane);
 866                         cnew.SetMultiplicity(cath, c->GetMultiplicity(cath));
 867                         for (i=0; i<fMul[cath]; i++) {
 868                             cnew.SetIndex(i, cath, c->GetIndex(i, cath));
 869                             fSeg[cath]->SetPad(fIx[i][cath], fIy[i][cath]);
 870                         }
 871                         if (fDebugLevel) {
 872                           fprintf(stderr,"\nRawCluster %d cath %d\n",ico,cath);
 873                           fprintf(stderr,"mult_av %d\n",c->GetMultiplicity(cath));
 874                         }
 875                         FillCluster(&cnew,cath);
 876                     }
 877                     cnew.SetClusterType(cnew.PhysicsContribution());
 878                     AddRawCluster(cnew);
 879                     fNPeaks++;
 880                 }
 881             }
 882         }
 883
 884 //  ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 885 //  (4) At least three local maxima on cathode 1 or on cathode 2
 886 //  ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
 887     } else if (fNLocal[0]>2 || fNLocal[1]>2) {
 888         Int_t param = fNLocal[0]*fNLocal[1];
 889         Int_t ii;
 890
 891         Float_t ** xm = new Float_t * [param];
 892         for (ii=0; ii<param; ii++) xm[ii]=new Float_t [2];
 893         Float_t ** ym = new Float_t * [param];
 894         for (ii=0; ii<param; ii++) ym[ii]=new Float_t [2];
 895         Int_t ** ixm = new Int_t * [param];
 896         for (ii=0; ii<param; ii++) ixm[ii]=new Int_t [2];
 897         Int_t ** iym = new Int_t * [param];
 898         for (ii=0; ii<param; ii++) iym[ii]=new Int_t [2];
 899
 900         Int_t isec, ico;
 901         Float_t dpx, dpy, dx, dy;
 902
 903         ico=0;
 904         for (Int_t im1=0; im1<fNLocal[0]; im1++) {
 905             for (Int_t im2=0; im2<fNLocal[1]; im2++) {
 906                 xm[ico][0]=fX[fIndLocal[im1][0]][0];
 907                 ym[ico][0]=fY[fIndLocal[im1][0]][0];
 908                 xm[ico][1]=fX[fIndLocal[im2][1]][1];
 909                 ym[ico][1]=fY[fIndLocal[im2][1]][1];
 910
 911                 ixm[ico][0]=fIx[fIndLocal[im1][0]][0];
 912                 iym[ico][0]=fIy[fIndLocal[im1][0]][0];
 913                 ixm[ico][1]=fIx[fIndLocal[im2][1]][1];
 914                 iym[ico][1]=fIy[fIndLocal[im2][1]][1];
 915                 ico++;
 916             }
 917         }
 918
 919         Int_t nIco = ico;
 920         if (fDebugLevel)
 921             fprintf(stderr,"nIco %d\n",nIco);
 922         for (ico=0; ico<nIco; ico++) {
 923             if (fDebugLevel)
 924                 fprintf(stderr,"ico = %d\n",ico);
 925             isec=fSeg[0]->Sector(ixm[ico][0], iym[ico][0]);
 926             dpx=fSeg[0]->Dpx(isec)/2.;
 927             dx=TMath::Abs(xm[ico][0]-xm[ico][1]);
 928             isec=fSeg[1]->Sector(ixm[ico][1], iym[ico][1]);
 929             dpy=fSeg[1]->Dpy(isec)/2.;
 930             dy=TMath::Abs(ym[ico][0]-ym[ico][1]);
 931             if (fDebugLevel) {
 932                 fprintf(stderr,"dx %f dpx %f dy %f dpy %f\n",dx,dpx,dy,dpy);
 933                 fprintf(stderr,"  X %f Y %f\n",xm[ico][1],ym[ico][0]);
 934             }
 935             if ((dx <= dpx) && (dy <= dpy)) {
 936                 if (fDebugLevel)
 937                     fprintf(stderr,"ok\n");
 938                 Int_t cath;
 939                 AliMUONRawCluster cnew;
 940                 for (cath=0; cath<2; cath++) {
 941                     cnew.SetX(cath, Float_t(xm[ico][1]));
 942                     cnew.SetY(cath, Float_t(ym[ico][0]));
 943                     cnew.SetZ(cath, fZPlane);
 944                     cnew.SetMultiplicity(cath, c->GetMultiplicity(cath));
 945                     for (i=0; i<fMul[cath]; i++) {
 946                         cnew.SetIndex(i, cath, c->GetIndex(i, cath));
 947                         fSeg[cath]->SetPad(fIx[i][cath], fIy[i][cath]);
 948                     }
 949                     FillCluster(&cnew,cath);
 950                 }
 951                 cnew.SetClusterType(cnew.PhysicsContribution());
 952                 AddRawCluster(cnew);
 953                 fNPeaks++;
 954             }
 955         }
 956         delete [] xm;
 957         delete [] ym;
 958         delete [] ixm;
 959         delete [] iym;
 960     }
 961 }
 962
 963 void AliMUONClusterFinderVS::FindLocalMaxima(AliMUONRawCluster* /*c*/)
 964 {
 965 // Find all local maxima of a cluster
 966     if (fDebugLevel)
 967         printf("\n Find Local maxima  !");
 968
 969     AliMUONDigit* digt;
 970
 971     Int_t cath, cath1; // loops over cathodes
 972     Int_t i;           // loops over digits
 973     Int_t j;           // loops over cathodes
 974 //
 975 //  Find local maxima
 976 //
 977 //  counters for number of local maxima
 978     fNLocal[0]=fNLocal[1]=0;
 979 //  flags digits as local maximum
 980     Bool_t isLocal[100][2];
 981     for (i=0; i<100;i++) {
 982         isLocal[i][0]=isLocal[i][1]=kFALSE;
 983     }
 984 //  number of next neighbours and arrays to store them
 985     Int_t nn;
 986     Int_t x[10], y[10];
 987 // loop over cathodes
 988     for (cath=0; cath<2; cath++) {
 989 // loop over cluster digits
 990         for (i=0; i<fMul[cath]; i++) {
 991 // get neighbours for that digit and assume that it is local maximum
 992             fSeg[cath]->Neighbours(fIx[i][cath], fIy[i][cath], &nn, x, y);
 993             isLocal[i][cath]=kTRUE;
 994             Int_t isec= fSeg[cath]->Sector(fIx[i][cath], fIy[i][cath]);
 995             Float_t a0 = fSeg[cath]->Dpx(isec)*fSeg[cath]->Dpy(isec);
 996 // loop over next neighbours, if at least one neighbour has higher charger assumption
 997 // digit is not local maximum
 998             for (j=0; j<nn; j++) {
 999                 if (fHitMap[cath]->TestHit(x[j], y[j])==kEmpty) continue;
1000                 digt=(AliMUONDigit*) fHitMap[cath]->GetHit(x[j], y[j]);
1001                 isec=fSeg[cath]->Sector(x[j], y[j]);
1002                 Float_t a1 = fSeg[cath]->Dpx(isec)*fSeg[cath]->Dpy(isec);
1003                 if (digt->Signal()/a1 > fQ[i][cath]/a0) {
1004                     isLocal[i][cath]=kFALSE;
1005                     break;
1006 //
1007 // handle special case of neighbouring pads with equal signal
1008                 } else if (digt->Signal() == fQ[i][cath]) {
1009                     if (fNLocal[cath]>0) {
1010                         for (Int_t k=0; k<fNLocal[cath]; k++) {
1011                             if (x[j]==fIx[fIndLocal[k][cath]][cath]
1012                                 && y[j]==fIy[fIndLocal[k][cath]][cath])
1013                             {
1014                                 isLocal[i][cath]=kFALSE;
1015                             }
1016                         } // loop over local maxima
1017                     } // are there already local maxima
1018                 } // same charge ?
1019             } // loop over next neighbours
1020             if (isLocal[i][cath]) {
1021                 fIndLocal[fNLocal[cath]][cath]=i;
1022                 fNLocal[cath]++;
1023             }
1024         } // loop over all digits
1025     } // loop over cathodes
1026
1027     if (fDebugLevel) {
1028         printf("\n Found %d %d %d %d local Maxima\n",
1029                fNLocal[0], fNLocal[1], fMul[0], fMul[1]);
1030         fprintf(stderr,"\n Cathode 1 local Maxima %d Multiplicite %d\n",fNLocal[0], fMul[0]);
1031         fprintf(stderr," Cathode 2 local Maxima %d Multiplicite %d\n",fNLocal[1], fMul[1]);
1032     }
1033     Int_t ix, iy, isec;
1034     Float_t dpx, dpy;
1035
1036
1037     if (fNLocal[1]==2 &&  (fNLocal[0]==1 || fNLocal[0]==0)) {
1038         Int_t iback=fNLocal[0];
1039
1040 //  Two local maxima on cathode 2 and one maximum on cathode 1
1041 //  Look for local maxima considering up and down neighbours on the 1st cathode only
1042 //
1043 //  Loop over cluster digits
1044         cath=0;
1045         cath1=1;
1046
1047         for (i=0; i<fMul[cath]; i++) {
1048             isec=fSeg[cath]->Sector(fIx[i][cath],fIy[i][cath]);
1049             dpy=fSeg[cath]->Dpy(isec);
1050             dpx=fSeg[cath]->Dpx(isec);
1051             if (isLocal[i][cath]) continue;
1052 // Pad position should be consistent with position of local maxima on the opposite cathode
1053             if ((TMath::Abs(fX[i][cath]-fX[fIndLocal[0][cath1]][cath1]) > dpx/2.) &&
1054                 (TMath::Abs(fX[i][cath]-fX[fIndLocal[1][cath1]][cath1]) > dpx/2.))
1055                 continue;
1056
1057 // get neighbours for that digit and assume that it is local maximum
1058             isLocal[i][cath]=kTRUE;
1059 // compare signal to that on the two neighbours on the left and on the right
1060 // iNN counts the number of neighbours with signal, it should be 1 or 2
1061             Int_t iNN=0;
1062
1063             for (fSeg[cath]
1064                      ->FirstPad(fX[i][cath], fY[i][cath], fZPlane, 0., dpy);
1065                  fSeg[cath]
1066                      ->MorePads();
1067                  fSeg[cath]
1068                      ->NextPad())
1069             {
1070                 ix = fSeg[cath]->Ix();
1071                 iy = fSeg[cath]->Iy();
1072                 // skip the current pad
1073                 if (iy == fIy[i][cath]) continue;
1074
1075                 if (fHitMap[cath]->TestHit(ix, iy)!=kEmpty) {
1076                     iNN++;
1077                     digt=(AliMUONDigit*) fHitMap[cath]->GetHit(ix,iy);
1078                     if (digt->Signal() > fQ[i][cath]) isLocal[i][cath]=kFALSE;
1079                 }
1080             } // Loop over pad neighbours in y
1081             if (isLocal[i][cath] && iNN>0) {
1082                 fIndLocal[fNLocal[cath]][cath]=i;
1083                 fNLocal[cath]++;
1084             }
1085         } // loop over all digits
1086 // if one additional maximum has been found we are happy
1087 // if more maxima have been found restore the previous situation
1088         if (fDebugLevel) {
1089             fprintf(stderr,
1090                     "\n New search gives %d local maxima for cathode 1 \n",
1091                     fNLocal[0]);
1092             fprintf(stderr,
1093                     "                  %d local maxima for cathode 2 \n",
1094                     fNLocal[1]);
1095         }
1096         if (fNLocal[cath]>2) {
1097             fNLocal[cath]=iback;
1098         }
1099
1100     } // 1,2 local maxima
1101
1102     if (fNLocal[0]==2 &&  (fNLocal[1]==1 || fNLocal[1]==0)) {
1103         Int_t iback=fNLocal[1];
1104
1105 //  Two local maxima on cathode 1 and one maximum on cathode 2
1106 //  Look for local maxima considering left and right neighbours on the 2nd cathode only
1107         cath=1;
1108         Int_t cath1 = 0;
1109         Float_t eps = 1.e-5;
1110
1111 //
1112 //  Loop over cluster digits
1113         for (i=0; i<fMul[cath]; i++) {
1114             isec=fSeg[cath]->Sector(fIx[i][cath],fIy[i][cath]);
1115             dpx=fSeg[cath]->Dpx(isec);
1116             dpy=fSeg[cath]->Dpy(isec);
1117             if (isLocal[i][cath]) continue;
1118 // Pad position should be consistent with position of local maxima on the opposite cathode
1119             if ((TMath::Abs(fY[i][cath]-fY[fIndLocal[0][cath1]][cath1]) > dpy/2.+eps) &&
1120                 (TMath::Abs(fY[i][cath]-fY[fIndLocal[1][cath1]][cath1]) > dpy/2.+eps))
1121                 continue;
1122
1123 //
1124 // get neighbours for that digit and assume that it is local maximum
1125             isLocal[i][cath]=kTRUE;
1126 // compare signal to that on the two neighbours on the left and on the right
1127
1128 // iNN counts the number of neighbours with signal, it should be 1 or 2
1129             Int_t iNN=0;
1130             for (fSeg[cath]
1131                      ->FirstPad(fX[i][cath], fY[i][cath], fZPlane, dpx, 0.);
1132                  fSeg[cath]
1133                      ->MorePads();
1134                  fSeg[cath]
1135                      ->NextPad())
1136             {
1137
1138                 ix = fSeg[cath]->Ix();
1139                 iy = fSeg[cath]->Iy();
1140
1141                 // skip the current pad
1142                 if (ix == fIx[i][cath]) continue;
1143
1144                 if (fHitMap[cath]->TestHit(ix, iy)!=kEmpty) {
1145                     iNN++;
1146                     digt=(AliMUONDigit*) fHitMap[cath]->GetHit(ix,iy);
1147                     if (digt->Signal() > fQ[i][cath]) isLocal[i][cath]=kFALSE;
1148                 }
1149             } // Loop over pad neighbours in x
1150             if (isLocal[i][cath] && iNN>0) {
1151                 fIndLocal[fNLocal[cath]][cath]=i;
1152                 fNLocal[cath]++;
1153             }
1154         } // loop over all digits
1155 // if one additional maximum has been found we are happy
1156 // if more maxima have been found restore the previous situation
1157         if (fDebugLevel) {
1158             fprintf(stderr,"\n New search gives %d local maxima for cathode 1 \n",fNLocal[0]);
1159             fprintf(stderr,"\n                  %d local maxima for cathode 2 \n",fNLocal[1]);
1160             printf("\n New search gives %d %d \n",fNLocal[0],fNLocal[1]);
1161         }
1162         if (fNLocal[cath]>2) {
1163             fNLocal[cath]=iback;
1164         }
1165     } // 2,1 local maxima
1166 }
1167
1168
1169 void  AliMUONClusterFinderVS::FillCluster(AliMUONRawCluster* c, Int_t flag, Int_t cath)
1170 {
1171 //
1172 //  Completes cluster information starting from list of digits
1173 //
1174     AliMUONDigit* dig;
1175     Float_t x, y, z;
1176     Int_t  ix, iy;
1177
1178     if (cath==1) {
1179         c->SetPeakSignal(cath,c->GetPeakSignal(0));
1180     } else {
1181         c->SetPeakSignal(cath,0);
1182     }
1183
1184
1185     if (flag) {
1186         c->SetX(cath,0.);
1187         c->SetY(cath,0.);
1188         c->SetCharge(cath,0);
1189     }
1190
1191     if (fDebugLevel)
1192         fprintf(stderr,"\n fPeakSignal %d\n",c->GetPeakSignal(cath));
1193     for (Int_t i=0; i<c->GetMultiplicity(cath); i++)
1194     {
1195         dig= fInput->Digit(cath,c->GetIndex(i,cath));
1196         ix=dig->PadX()+c->GetOffset(i,cath);
1197         iy=dig->PadY();
1198         Int_t q=dig->Signal();
1199         if (!flag) q=Int_t(q*c->GetContrib(i,cath));
1200 //      fprintf(stderr,"q %d c->fPeakSignal[ %d ] %d\n",q,cath,c->fPeakSignal[cath]);
1201         if (dig->Physics() >= dig->Signal()) {
1202             c->SetPhysics(i,2);
1203         } else if (dig->Physics() == 0) {
1204             c->SetPhysics(i,0);
1205         } else  c->SetPhysics(i,1);
1206 //
1207 //
1208         if (fDebugLevel>1)
1209             fprintf(stderr,"q %d c->fPeakSignal[cath] %d\n",q,c->GetPeakSignal(cath));
1210 // peak signal and track list
1211         if (q>c->GetPeakSignal(cath)) {
1212             c->SetPeakSignal(cath, q);
1213             c->SetTrack(0,dig->Hit());
1214             c->SetTrack(1,dig->Track(0));
1215             c->SetTrack(2,dig->Track(1));
1216 //          fprintf(stderr," c->fTracks[0] %d c->fTracks[1] %d\n",dig->fHit,dig->fTracks[0]);
1217         }
1218 //
1219         if (flag) {
1220             fSeg[cath]->GetPadC(ix, iy, x, y, z);
1221             c->AddX(cath, q*x);
1222             c->AddY(cath, q*y);
1223             c->AddCharge(cath, q);
1224         }
1225     } // loop over digits
1226     if (fDebugLevel)
1227         fprintf(stderr," fin du cluster c\n");
1228
1229
1230     if (flag) {
1231         c->SetX(cath, c->GetX(cath)/c->GetCharge(cath));
1232 // Force on anod
1233         c->SetX(cath, fSeg[cath]->GetAnod(c->GetX(cath)));
1234         c->SetY(cath, c->GetY(cath)/c->GetCharge(cath));
1235 //
1236 //  apply correction to the coordinate along the anode wire
1237 //
1238         x=c->GetX(cath);
1239         y=c->GetY(cath);
1240         fSeg[cath]->GetPadI(x, y, fZPlane, ix, iy);
1241         fSeg[cath]->GetPadC(ix, iy, x, y, z);
1242         Int_t isec=fSeg[cath]->Sector(ix,iy);
1243         TF1* cogCorr = fSeg[cath]->CorrFunc(isec-1);
1244
1245         if (cogCorr) {
1246             Float_t yOnPad=(c->GetY(cath)-y)/fSeg[cath]->Dpy(isec);
1247             c->SetY(cath, c->GetY(cath)-cogCorr->Eval(yOnPad, 0, 0));
1248         }
1249     }
1250 }
1251
1252 void  AliMUONClusterFinderVS::FillCluster(AliMUONRawCluster* c, Int_t cath)
1253 {
1254 //
1255 //  Completes cluster information starting from list of digits
1256 //
1257     static Float_t dr0;
1258
1259     AliMUONDigit* dig;
1260
1261     if (cath==0) {
1262         dr0 = 10000;
1263     }
1264
1265     Float_t xpad, ypad, zpad;
1266     Float_t dx, dy, dr;
1267
1268     for (Int_t i=0; i<c->GetMultiplicity(cath); i++)
1269     {
1270         dig = fInput->Digit(cath,c->GetIndex(i,cath));
1271         fSeg[cath]->
1272         GetPadC(dig->PadX(),dig->PadY(),xpad,ypad, zpad);
1273         if (fDebugLevel)
1274             fprintf(stderr,"x %f y %f cx %f cy %f\n",xpad,ypad,c->GetX(0),c->GetY(0));
1275         dx = xpad - c->GetX(0);
1276         dy = ypad - c->GetY(0);
1277         dr = TMath::Sqrt(dx*dx+dy*dy);
1278
1279         if (dr < dr0) {
1280             dr0 = dr;
1281             if (fDebugLevel)
1282                 fprintf(stderr," dr %f\n",dr);
1283             Int_t q=dig->Signal();
1284             if (dig->Physics() >= dig->Signal()) {
1285                 c->SetPhysics(i,2);
1286             } else if (dig->Physics() == 0) {
1287                 c->SetPhysics(i,0);
1288             } else  c->SetPhysics(i,1);
1289             c->SetPeakSignal(cath,q);
1290             c->SetTrack(0,dig->Hit());
1291             c->SetTrack(1,dig->Track(0));
1292             c->SetTrack(2,dig->Track(1));
1293             if (fDebugLevel)
1294                 fprintf(stderr," c->fTracks[0] %d c->fTracks[1] %d\n",dig->Hit(),
1295                     dig->Track(0));
1296         }
1297 //
1298     } // loop over digits
1299
1300 //  apply correction to the coordinate along the anode wire
1301 // Force on anod
1302     c->SetX(cath,fSeg[cath]->GetAnod(c->GetX(cath)));
1303 }
1304
1305 void  AliMUONClusterFinderVS::FindCluster(Int_t i, Int_t j, Int_t cath, AliMUONRawCluster &c){
1306
1307
1308 //
1309 //  Find a super cluster on both cathodes
1310 //
1311 //
1312 //  Add i,j as element of the cluster
1313 //
1314
1315     Int_t idx = fHitMap[cath]->GetHitIndex(i,j);
1316     AliMUONDigit* dig = (AliMUONDigit*) fHitMap[cath]->GetHit(i,j);
1317     Int_t q=dig->Signal();
1318     Int_t theX=dig->PadX();
1319     Int_t theY=dig->PadY();
1320
1321     if (q > TMath::Abs(c.GetPeakSignal(0)) && q > TMath::Abs(c.GetPeakSignal(1))) {
1322         c.SetPeakSignal(cath,q);
1323         c.SetTrack(0,dig->Hit());
1324         c.SetTrack(1,dig->Track(0));
1325         c.SetTrack(2,dig->Track(1));
1326     }
1327
1328 //
1329 //  Make sure that list of digits is ordered
1330 //
1331     Int_t mu=c.GetMultiplicity(cath);
1332     c.SetIndex(mu, cath, idx);
1333
1334     if (dig->Physics() >= dig->Signal()) {
1335         c.SetPhysics(mu,2);
1336     } else if (dig->Physics() == 0) {
1337         c.SetPhysics(mu,0);
1338     } else  c.SetPhysics(mu,1);
1339
1340
1341     if (mu > 0) {
1342         for (Int_t ind = mu-1; ind >= 0; ind--) {
1343             Int_t ist=c.GetIndex(ind,cath);
1344             Int_t ql=fInput->Digit(cath, ist)->Signal();
1345             Int_t ix=fInput->Digit(cath, ist)->PadX();
1346             Int_t iy=fInput->Digit(cath, ist)->PadY();
1347
1348             if (q>ql || (q==ql && theX > ix && theY < iy)) {
1349                 c.SetIndex(ind, cath, idx);
1350                 c.SetIndex(ind+1, cath, ist);
1351             } else {
1352
1353                 break;
1354             }
1355         }
1356     }
1357
1358     c.SetMultiplicity(cath, c.GetMultiplicity(cath)+1);
1359     if (c.GetMultiplicity(cath) >= 50 ) {
1360       if (fDebugLevel)
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   AliFatal("Not implemented.");
2105
2106   return *this;
2107 }
2108
2109
2110
2111
2112
2113
2114
2115
2116