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