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