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