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