MUON/AliMUONClusterSplitterMLEM.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 AliMUONClusterSplitterMLEM
  20 ///
  21 /// Splitter class for the MLEM algorithm. Performs fitting procedure
  22 /// with up to 3 hit candidates and tries to split clusters if the number
  23 /// of candidates exceeds 3.
  24 ///
  25 /// \author Laurent Aphecetche (for the "new" C++ structure) and
  26 /// Alexander Zinchenko, JINR Dubna, for the hardcore of it ;-)
  27 //-----------------------------------------------------------------------------
  28
  29 #include "AliMUONClusterSplitterMLEM.h"
  30 #include "AliMUONClusterFinderMLEM.h" // for status flag constants
  31
  32 #include "AliMUONCluster.h"
  33 #include "AliMUONPad.h"
  34 #include "AliMUONPad.h"
  35 #include "AliMUONConstants.h"
  36 #include "AliMpDEManager.h"
  37 #include "AliMUONMathieson.h"
  38
  39 #include "AliLog.h"
  40
  41 #include <TClonesArray.h>
  42 #include <TH2.h>
  43 #include <TMath.h>
  44 #include <TMatrixD.h>
  45 #include <TObjArray.h>
  46 #include <TRandom.h>
  47 #include <Riostream.h>
  48
  49 /// \cond CLASSIMP
  50 ClassImp(AliMUONClusterSplitterMLEM)
  51 /// \endcond
  52
  53 //const Double_t AliMUONClusterSplitterMLEM::fgkCouplMin = 1.e-3; // threshold on coupling
  54 const Double_t AliMUONClusterSplitterMLEM::fgkCouplMin = 1.e-2; // threshold on coupling
  55
  56 //_____________________________________________________________________________
  57 AliMUONClusterSplitterMLEM::AliMUONClusterSplitterMLEM(Int_t detElemId,
  58                                                        TObjArray* pixArray)
  59 : TObject(),
  60 fPixArray(pixArray),
  61 fMathieson(0x0),
  62 fDetElemId(detElemId),
  63 fNpar(0),
  64 fQtot(0),
  65 fnCoupled(0),
  66 fDebug(0)
  67 {
  68   /// Constructor
  69
  70   AliMq::Station12Type stationType = AliMpDEManager::GetStation12Type(fDetElemId);
  71
  72   Float_t kx3 = AliMUONConstants::SqrtKx3();
  73   Float_t ky3 = AliMUONConstants::SqrtKy3();
  74   Float_t pitch = AliMUONConstants::Pitch();
  75
  76   if ( stationType == AliMq::kStation1 )
  77   {
  78     kx3 = AliMUONConstants::SqrtKx3St1();
  79     ky3 = AliMUONConstants::SqrtKy3St1();
  80     pitch = AliMUONConstants::PitchSt1();
  81   }
  82
  83   fMathieson = new AliMUONMathieson;
  84
  85   fMathieson->SetPitch(pitch);
  86   fMathieson->SetSqrtKx3AndDeriveKx2Kx4(kx3);
  87   fMathieson->SetSqrtKy3AndDeriveKy2Ky4(ky3);
  88
  89 }
  90
  91 //_____________________________________________________________________________
  92 AliMUONClusterSplitterMLEM::~AliMUONClusterSplitterMLEM()
  93 {
  94   /// Destructor
  95
  96   delete fMathieson;
  97 }
  98
  99 //_____________________________________________________________________________
 100 void
 101 AliMUONClusterSplitterMLEM::AddBin(TH2 *mlem,
 102                                    Int_t ic, Int_t jc, Int_t mode,
 103                                    Bool_t *used, TObjArray *pix)
 104 {
 105   /// Add a bin to the cluster
 106
 107   Int_t nx = mlem->GetNbinsX();
 108   Int_t ny = mlem->GetNbinsY();
 109   Double_t cont1, cont = mlem->GetCellContent(jc,ic);
 110   AliMUONPad *pixPtr = 0;
 111
 112   Int_t ie = TMath::Min(ic+1,ny), je = TMath::Min(jc+1,nx);
 113   for (Int_t i = TMath::Max(ic-1,1); i <= ie; ++i) {
 114     for (Int_t j = TMath::Max(jc-1,1); j <= je; ++j) {
 115       if (i != ic && j != jc) continue;
 116       if (used[(i-1)*nx+j-1]) continue;
 117       cont1 = mlem->GetCellContent(j,i);
 118       if (mode && cont1 > cont) continue;
 119       used[(i-1)*nx+j-1] = kTRUE;
 120       if (cont1 < 0.5) continue;
 121       if (pix) pix->Add(BinToPix(mlem,j,i));
 122       else {
 123         pixPtr = new AliMUONPad (mlem->GetXaxis()->GetBinCenter(j),
 124                                  mlem->GetYaxis()->GetBinCenter(i), 0, 0, cont1);
 125         fPixArray->Add(pixPtr);
 126       }
 127       AddBin(mlem, i, j, mode, used, pix); // recursive call
 128     }
 129   }
 130 }
 131
 132 //_____________________________________________________________________________
 133 void
 134 AliMUONClusterSplitterMLEM::AddCluster(Int_t ic, Int_t nclust,
 135                                        TMatrixD& aijcluclu,
 136                                        Bool_t *used, Int_t *clustNumb, Int_t &nCoupled)
 137 {
 138   /// Add a cluster to the group of coupled clusters
 139
 140   for (Int_t i = 0; i < nclust; ++i) {
 141     if (used[i]) continue;
 142     if (aijcluclu(i,ic) < fgkCouplMin) continue;
 143     used[i] = kTRUE;
 144     clustNumb[nCoupled++] = i;
 145     AddCluster(i, nclust, aijcluclu, used, clustNumb, nCoupled);
 146   }
 147 }
 148
 149 //_____________________________________________________________________________
 150 TObject*
 151 AliMUONClusterSplitterMLEM::BinToPix(TH2 *mlem,
 152                                      Int_t jc, Int_t ic)
 153 {
 154   /// Translate histogram bin to pixel
 155
 156   Double_t yc = mlem->GetYaxis()->GetBinCenter(ic);
 157   Double_t xc = mlem->GetXaxis()->GetBinCenter(jc);
 158
 159   Int_t nPix = fPixArray->GetEntriesFast();
 160   AliMUONPad *pixPtr = NULL;
 161
 162   // Compare pixel and bin positions
 163   for (Int_t i = 0; i < nPix; ++i) {
 164     pixPtr = (AliMUONPad*) fPixArray->UncheckedAt(i);
 165     if (pixPtr->Charge() < 0.5) continue;
 166     if (TMath::Abs(pixPtr->Coord(0)-xc)<1.e-4 && TMath::Abs(pixPtr->Coord(1)-yc)<1.e-4)
 167     {
 168       //return (TObject*) pixPtr;
 169       return pixPtr;
 170     }
 171   }
 172   AliError(Form(" Something wrong ??? %f %f ", xc, yc));
 173   return NULL;
 174 }
 175
 176 //_____________________________________________________________________________
 177 Float_t
 178 AliMUONClusterSplitterMLEM::ChargeIntegration(Double_t x, Double_t y,
 179                                               const AliMUONPad& pad)
 180 {
 181   /// Compute the Mathieson integral on pad area, assuming the center
 182   /// of the Mathieson is at (x,y)
 183
 184   TVector2 lowerLeft(TVector2(x,y)-pad.Position()-pad.Dimensions());
 185   TVector2 upperRight(lowerLeft + pad.Dimensions()*2.0);
 186
 187         return fMathieson->IntXY(lowerLeft.X(),lowerLeft.Y(),
 188                            upperRight.X(),upperRight.Y());
 189 }
 190
 191 //_____________________________________________________________________________
 192 void
 193 AliMUONClusterSplitterMLEM::Fcn1(const AliMUONCluster& cluster,
 194                                     Int_t & /*fNpar*/, Double_t * /*gin*/,
 195                                     Double_t &f, Double_t *par, Int_t iflag)
 196 {
 197   /// Computes the functional to be minimized
 198
 199   Int_t indx, npads=0;
 200   Double_t charge, delta, coef=0, chi2=0, qTot = 0;
 201   static Double_t qAver = 0;
 202
 203   Int_t mult = cluster.Multiplicity(), iend = fNpar / 3;
 204   for (Int_t j = 0; j < mult; ++j)
 205   {
 206     AliMUONPad* pad = cluster.Pad(j);
 207     //if ( pad->Status() !=1 || pad->IsSaturated() ) continue;
 208     if ( pad->Status() != AliMUONClusterFinderMLEM::GetUseForFitFlag()) continue;
 209     if (iflag == 0) {
 210       if ( pad->IsReal() ) npads++; // exclude virtual pads
 211       qTot += pad->Charge();
 212     }
 213     charge = 0;
 214     for (Int_t i = 0; i <= iend; ++i)
 215     {
 216       // sum over hits
 217       indx = 3 * i;
 218       coef = Param2Coef(i, coef, par);
 219       charge += ChargeIntegration(par[indx],par[indx+1],*pad) * coef;
 220     }
 221     charge *= fQtot;
 222     delta = charge - pad->Charge();
 223     delta *= delta;
 224     delta /= pad->Charge();
 225     chi2 += delta;
 226   } // for (Int_t j=0;
 227   if (iflag == 0) qAver = qTot / npads;
 228   f = chi2 / qAver;
 229 }
 230
 231 //_____________________________________________________________________________
 232 Double_t AliMUONClusterSplitterMLEM::Param2Coef(Int_t icand, Double_t coef, Double_t *par)
 233 {
 234   /// Extract hit contribution scale factor from fit parameters
 235
 236   if (fNpar == 2) return 1.;
 237   if (fNpar == 5) return icand==0 ? par[2] : TMath::Max(1.-par[2],0.);
 238   if (icand == 0) return par[2];
 239   if (icand == 1) return TMath::Max((1.-par[2])*par[5], 0.);
 240   return TMath::Max(1.-par[2]-coef,0.);
 241 }
 242
 243 //_____________________________________________________________________________
 244 Int_t
 245 AliMUONClusterSplitterMLEM::Fit(const AliMUONCluster& cluster,
 246                                 Int_t iSimple, Int_t nfit,
 247                                 Int_t *clustFit, TObjArray **clusters,
 248                                 Double_t *parOk,
 249                                 TObjArray& clusterList, TH2 *mlem)
 250 {
 251   /// Steering function and fitting procedure for the fit of pad charge distribution
 252
 253   //  AliDebug(2,Form("iSimple=%d nfit=%d",iSimple,nfit));
 254
 255   Double_t xmin = mlem->GetXaxis()->GetXmin() - mlem->GetXaxis()->GetBinWidth(1);
 256   Double_t xmax = mlem->GetXaxis()->GetXmax() + mlem->GetXaxis()->GetBinWidth(1);
 257   Double_t ymin = mlem->GetYaxis()->GetXmin() - mlem->GetYaxis()->GetBinWidth(1);
 258   Double_t ymax = mlem->GetYaxis()->GetXmax() + mlem->GetYaxis()->GetBinWidth(1);
 259   Double_t xPad = 0, yPad = 99999;
 260
 261   // Number of pads to use and number of virtual pads
 262   Int_t npads = 0, nVirtual = 0, nfit0 = nfit;
 263   //cluster.Print("full");
 264   Int_t mult = cluster.Multiplicity();
 265   for (Int_t i = 0; i < mult; ++i )
 266   {
 267     AliMUONPad* pad = cluster.Pad(i);
 268     if ( !pad->IsReal() ) ++nVirtual;
 269     //if ( pad->Status() !=1 || pad->IsSaturated() ) continue;
 270     if ( pad->Status() != AliMUONClusterFinderMLEM::GetUseForFitFlag() ) continue;
 271     if ( pad->IsReal() )
 272     {
 273       ++npads;
 274       if (yPad > 9999)
 275       {
 276         xPad = pad->X();
 277         yPad = pad->Y();
 278       }
 279       else
 280       {
 281         if (pad->DY() < pad->DX() )
 282         {
 283           yPad = pad->Y();
 284         }
 285         else
 286         {
 287           xPad = pad->X();
 288         }
 289       }
 290     }
 291   }
 292
 293   fNpar = 0;
 294   fQtot = 0;
 295
 296   if (npads < 2) return 0;
 297
 298   // FIXME : AliWarning("Reconnect the following code for hit/track passing ?");
 299
 300   //  Int_t tracks[3] = {-1, -1, -1};
 301
 302   /*
 303    Int_t digit = 0;
 304    AliMUONDigit *mdig = 0;
 305    for (Int_t cath=0; cath<2; cath++) {
 306      for (Int_t i=0; i<fnPads[0]+fnPads[1]; i++) {
 307        if (fPadIJ[0][i] != cath) continue;
 308        if (fPadIJ[1][i] != 1) continue;
 309        if (fXyq[3][i] < 0) continue; // exclude virtual pads
 310        digit = TMath::Nint (fXyq[5][i]);
 311        if (digit >= 0) mdig = fInput->Digit(cath,digit);
 312        else mdig = fInput->Digit(TMath::Even(cath),-digit-1);
 313        //if (!mdig) mdig = fInput->Digit(TMath::Even(cath),digit);
 314        if (!mdig) continue; // protection for cluster display
 315        if (mdig->Hit() >= 0) {
 316          if (tracks[0] < 0) {
 317            tracks[0] = mdig->Hit();
 318            tracks[1] = mdig->Track(0);
 319          } else if (mdig->Track(0) < tracks[1]) {
 320            tracks[0] = mdig->Hit();
 321            tracks[1] = mdig->Track(0);
 322          }
 323        }
 324        if (mdig->Track(1) >= 0 && mdig->Track(1) != tracks[1]) {
 325          if (tracks[2] < 0) tracks[2] = mdig->Track(1);
 326          else tracks[2] = TMath::Min (tracks[2], mdig->Track(1));
 327        }
 328      } // for (Int_t i=0;
 329   } // for (Int_t cath=0;
 330    */
 331
 332   // Get number of pads in X and Y
 333   //const Int_t kStatusToTest(1);
 334   const Int_t kStatusToTest(AliMUONClusterFinderMLEM::GetUseForFitFlag());
 335
 336   AliMpIntPair nofPads = cluster.NofPads(kStatusToTest);
 337   Int_t nInX = nofPads.GetFirst();
 338   Int_t nInY = nofPads.GetSecond();
 339
 340   if (fDebug) {
 341     Int_t npadOK = 0;
 342     for (Int_t j = 0; j < cluster.Multiplicity(); ++j) {
 343       AliMUONPad *pad = cluster.Pad(j);
 344       //if (pad->Status() == 1 && !pad->IsSaturated()) npadOK++;
 345       if (pad->Status() == AliMUONClusterFinderMLEM::GetUseForFitFlag() && !pad->IsSaturated()) npadOK++;
 346     }
 347     cout << " Number of pads to fit: " << npadOK << endl;
 348     cout << " nInX and Y: " << nInX << " " << nInY << endl;
 349   }
 350
 351   Int_t nfitMax = 3;
 352   nfitMax = TMath::Min (nfitMax, (npads + 1) / 3);
 353   if (nfitMax > 1) {
 354     if (((nInX < 3) && (nInY < 3)) || ((nInX == 3) && (nInY < 3)) || ((nInX < 3) && (nInY == 3))) nfitMax = 1; // not enough pads in each direction
 355   }
 356   if (nfit > nfitMax) nfit = nfitMax;
 357
 358   // Take cluster maxima as fitting seeds
 359   TObjArray *pix;
 360   AliMUONPad *pixPtr;
 361   Int_t npxclu;
 362   Double_t cont, cmax = 0, xseed = 0, yseed = 0, errOk[8], qq = 0;
 363   Double_t xyseed[3][2], qseed[3], xyCand[3][2] = {{0},{0}}, sigCand[3][2] = {{0},{0}};
 364
 365   for (Int_t ifit = 1; ifit <= nfit0; ++ifit)
 366   {
 367     cmax = 0;
 368     pix = clusters[clustFit[ifit-1]];
 369     npxclu = pix->GetEntriesFast();
 370     //qq = 0;
 371     for (Int_t clu = 0; clu < npxclu; ++clu)
 372     {
 373       pixPtr = (AliMUONPad*) pix->UncheckedAt(clu);
 374       cont = pixPtr->Charge();
 375       fQtot += cont;
 376       if (cont > cmax)
 377       {
 378         cmax = cont;
 379         xseed = pixPtr->Coord(0);
 380         yseed = pixPtr->Coord(1);
 381       }
 382       qq += cont;
 383       xyCand[0][0] += pixPtr->Coord(0) * cont;
 384       xyCand[0][1] += pixPtr->Coord(1) * cont;
 385       sigCand[0][0] += pixPtr->Coord(0) * pixPtr->Coord(0) * cont;
 386       sigCand[0][1] += pixPtr->Coord(1) * pixPtr->Coord(1) * cont;
 387     }
 388     xyseed[ifit-1][0] = xseed;
 389     xyseed[ifit-1][1] = yseed;
 390     qseed[ifit-1] = cmax;
 391   } // for (Int_t ifit=1;
 392
 393   xyCand[0][0] /= qq; // <x>
 394   xyCand[0][1] /= qq; // <y>
 395   sigCand[0][0] = sigCand[0][0]/qq - xyCand[0][0]*xyCand[0][0]; // <x^2> - <x>^2
 396   sigCand[0][0] = sigCand[0][0] > 0 ? TMath::Sqrt (sigCand[0][0]) : 0;
 397   sigCand[0][1] = sigCand[0][1]/qq - xyCand[0][1]*xyCand[0][1]; // <y^2> - <y>^2
 398   sigCand[0][1] = sigCand[0][1] > 0 ? TMath::Sqrt (sigCand[0][1]) : 0;
 399   if (fDebug) cout << xyCand[0][0] << " " << xyCand[0][1] << " " << sigCand[0][0] << " " << sigCand[0][1] << endl;
 400
 401   Int_t nDof, maxSeed[3];//, nMax = 0;
 402
 403   TMath::Sort(nfit0, qseed, maxSeed, kTRUE); // in decreasing order
 404
 405   Double_t step[3]={0.01,0.002,0.02}, fmin, chi2o = 9999, chi2n;
 406   Double_t *gin = 0, func0, func1, param[8], step0[8];
 407   Double_t param0[2][8]={{0},{0}}, deriv[2][8]={{0},{0}};
 408   Double_t shift[8], stepMax, derMax, parmin[8], parmax[8], func2[2], shift0;
 409   Double_t delta[8], scMax, dder[8], estim, shiftSave = 0;
 410   Int_t min, max, nCall = 0, nLoop, idMax = 0, iestMax = 0, nFail;
 411   Double_t rad, dist[3] = {0};
 412
 413   // Try to fit with one-track hypothesis, then 2-track. If chi2/dof is
 414   // lower, try 3-track (if number of pads is sufficient).
 415   Int_t iflag = 0; // for the first call of fcn1
 416   for (Int_t iseed = 0; iseed < nfit; ++iseed)
 417   {
 418
 419     Int_t memory[8] = {0};
 420     if (iseed)
 421     {
 422       for (Int_t j = 0; j < fNpar; ++j)
 423       {
 424         param[j] = parOk[j];
 425       }
 426       param[fNpar] = 0.6;
 427       parmin[fNpar] = 0;
 428       parmax[fNpar++] = 1;
 429     }
 430
 431     if (nfit == 1)
 432     {
 433       param[fNpar] = xyCand[0][0]; // take COG
 434     }
 435     else
 436     {
 437       param[fNpar] = xyseed[maxSeed[iseed]][0];
 438       //param[fNpar] = fNpar==0 ? -16.1651 : -15.2761;
 439     }
 440     parmin[fNpar] = xmin;
 441     parmax[fNpar++] = xmax;
 442     if (nfit == 1)
 443     {
 444       param[fNpar] = xyCand[0][1]; // take COG
 445     }
 446     else
 447     {
 448       param[fNpar] = xyseed[maxSeed[iseed]][1];
 449       //param[fNpar] = fNpar==1 ? -15.1737 : -15.8487;
 450     }
 451     parmin[fNpar] = ymin;
 452     parmax[fNpar++] = ymax;
 453
 454     for (Int_t j = 0; j < fNpar; ++j)
 455     {
 456       step0[j] = shift[j] = step[j%3];
 457     }
 458
 459     if (iseed)
 460     {
 461       for (Int_t j = 0; j < fNpar; ++j)
 462       {
 463         param0[1][j] = 0;
 464       }
 465     }
 466     if (fDebug) {
 467       for (Int_t j = 0; j < fNpar; ++j) cout << param[j] << " ";
 468       cout << endl;
 469     }
 470
 471     // Try new algorithm
 472     min = nLoop = 1; stepMax = func2[1] = derMax = 999999; nFail = 0;
 473
 474     while (1)
 475     {
 476       max = !min;
 477       Fcn1(cluster,fNpar, gin, func0, param, iflag); nCall++;
 478       iflag = 1;
 479       //cout << " Func: " << func0 << endl;
 480
 481       func2[max] = func0;
 482       for (Int_t j = 0; j < fNpar; ++j)
 483       {
 484         param0[max][j] = param[j];
 485         delta[j] = step0[j];
 486         param[j] += delta[j] / 10;
 487         if (j > 0) param[j-1] -= delta[j-1] / 10;
 488         Fcn1(cluster,fNpar, gin, func1, param, iflag); nCall++;
 489         deriv[max][j] = (func1 - func0) / delta[j] * 10; // first derivative
 490         //cout << j << " " << deriv[max][j] << endl;
 491         dder[j] = param0[0][j] != param0[1][j] ? (deriv[0][j] - deriv[1][j]) /
 492           (param0[0][j] - param0[1][j]) : 0; // second derivative
 493       }
 494       param[fNpar-1] -= delta[fNpar-1] / 10;
 495       if (nCall > 2000) break;
 496
 497       min = func2[0] < func2[1] ? 0 : 1;
 498       nFail = min == max ? 0 : nFail + 1;
 499
 500       stepMax = derMax = estim = 0;
 501       for (Int_t j = 0; j < fNpar; ++j)
 502       {
 503         // Estimated distance to minimum
 504         shift0 = shift[j];
 505         if (nLoop == 1)
 506         {
 507           shift[j] = TMath::Sign (step0[j], -deriv[max][j]); // first step
 508         }
 509         else if (TMath::Abs(deriv[0][j]) < 1.e-3 && TMath::Abs(deriv[1][j]) < 1.e-3)
 510         {
 511           shift[j] = 0;
 512         }
 513         else if (((deriv[min][j]*deriv[!min][j] > 0) && (TMath::Abs(deriv[min][j]) > TMath::Abs(deriv[!min][j])))
 514                  || (TMath::Abs(deriv[0][j]-deriv[1][j]) < 1.e-3) || (TMath::Abs(dder[j]) < 1.e-6))
 515         {
 516           shift[j] = -TMath::Sign (shift[j], (func2[0]-func2[1]) * (param0[0][j]-param0[1][j]));
 517           if (min == max)
 518           {
 519             if (memory[j] > 1)
 520             {
 521               shift[j] *= 2;
 522             }
 523             memory[j]++;
 524           }
 525         }
 526         else
 527         {
 528           shift[j] = dder[j] != 0 ? -deriv[min][j] / dder[j] : 0;
 529           memory[j] = 0;
 530         }
 531
 532         Double_t es = TMath::Abs(shift[j]) / step0[j];
 533         if (es > estim)
 534         {
 535           estim = es;
 536           iestMax = j;
 537         }
 538
 539         // Too big step
 540         if (TMath::Abs(shift[j])/step0[j] > 10) shift[j] = TMath::Sign(10.,shift[j]) * step0[j]; //
 541
 542         // Failed to improve minimum
 543         if (min != max)
 544         {
 545           memory[j] = 0;
 546           param[j] = param0[min][j];
 547           if (TMath::Abs(shift[j]+shift0) > 0.1*step0[j])
 548           {
 549             shift[j] = (shift[j] + shift0) / 2;
 550           }
 551           else
 552           {
 553             shift[j] /= -2;
 554           }
 555         }
 556
 557         // Too big step
 558         if (TMath::Abs(shift[j]*deriv[min][j]) > func2[min])
 559         {
 560           shift[j] = TMath::Sign (func2[min]/deriv[min][j], shift[j]);
 561         }
 562
 563         // Introduce step relaxation factor
 564         if (memory[j] < 3)
 565         {
 566           scMax = 1 + 4 / TMath::Max(nLoop/2.,1.);
 567           if (TMath::Abs(shift0) > 0 && TMath::Abs(shift[j]/shift0) > scMax)
 568           {
 569             shift[j] = TMath::Sign (shift0*scMax, shift[j]);
 570           }
 571         }
 572         param[j] += shift[j];
 573         // Check parameter limits
 574         if (param[j] < parmin[j])
 575         {
 576           shift[j] = parmin[j] - param[j];
 577           param[j] = parmin[j];
 578         }
 579         else if (param[j] > parmax[j])
 580         {
 581           shift[j] = parmax[j] - param[j];
 582           param[j] = parmax[j];
 583         }
 584         //cout << " xxx " << j << " " << shift[j] << " " << param[j] << endl;
 585         stepMax = TMath::Max (stepMax, TMath::Abs(shift[j]/step0[j]));
 586         if (TMath::Abs(deriv[min][j]) > derMax)
 587         {
 588           idMax = j;
 589           derMax = TMath::Abs (deriv[min][j]);
 590         }
 591       } // for (Int_t j=0; j<fNpar;
 592
 593       if (((estim < 1) && (derMax < 2)) || nLoop > 150) break; // minimum was found
 594
 595       nLoop++;
 596
 597       // Check for small step
 598       if (shift[idMax] == 0)
 599       {
 600         shift[idMax] = step0[idMax]/10;
 601         param[idMax] += shift[idMax];
 602         continue;
 603       }
 604
 605       if (!memory[idMax] && derMax > 0.5 && nLoop > 10)
 606       {
 607         if (dder[idMax] != 0 && TMath::Abs(deriv[min][idMax]/dder[idMax]/shift[idMax]) > 10)
 608         {
 609           if (min == max) dder[idMax] = -dder[idMax];
 610           shift[idMax] = -deriv[min][idMax] / dder[idMax] / 10;
 611           param[idMax] += shift[idMax];
 612           stepMax = TMath::Max (stepMax, TMath::Abs(shift[idMax])/step0[idMax]);
 613           if (min == max) shiftSave = shift[idMax];
 614         }
 615         if (nFail > 10)
 616         {
 617           param[idMax] -= shift[idMax];
 618           shift[idMax] = 4 * shiftSave * (gRandom->Rndm(0) - 0.5);
 619           param[idMax] += shift[idMax];
 620         }
 621       }
 622     } // while (1)
 623
 624     fmin = func2[min];
 625
 626     nDof = npads - fNpar + nVirtual;
 627     if (!nDof) nDof++;
 628     chi2n = fmin / nDof;
 629     if (fDebug) cout << " Chi2 " << chi2n << " " << fNpar << endl;
 630
 631     //if (fNpar > 2) cout << param0[min][fNpar-3] << " " << chi2n * (1+TMath::Min(1-param0[min][fNpar-3],0.25)) << endl;
 632     //if (chi2n*1.2+1.e-6 > chi2o )
 633     if (fNpar > 2 && (chi2n > chi2o || ((iseed == nfit-1)
 634                                         && (chi2n * (1+TMath::Min(1-param0[min][fNpar-3],0.25)) > chi2o))))
 635       { fNpar -= 3; break; }
 636
 637     // Save parameters and errors
 638
 639     if (nInX == 1) {
 640       // One pad per direction
 641       //for (Int_t i=0; i<fNpar; ++i) if (i == 0 || i == 2 || i == 5) param0[min][i] = xPad;
 642       for (Int_t i=0; i<fNpar; ++i) if (i == 0 || i == 2 || i == 5)
 643         param0[min][i] = xyCand[0][0];
 644     }
 645     if (nInY == 1) {
 646       // One pad per direction
 647       //for (Int_t i=0; i<fNpar; ++i) if (i == 1 || i == 3 || i == 6) param0[min][i] = yPad;
 648       for (Int_t i=0; i<fNpar; ++i) if (i == 1 || i == 3 || i == 6)
 649         param0[min][i] = xyCand[0][1];
 650     }
 651
 652     /*
 653       if (iseed > 0) {
 654       // Find distance to the nearest neighbour
 655       dist[0] = dist[1] = TMath::Sqrt ((param0[min][0]-param0[min][2])*
 656       (param0[min][0]-param0[min][2])
 657       +(param0[min][1]-param0[min][3])*
 658       (param0[min][1]-param0[min][3]));
 659       if (iseed > 1) {
 660       dist[2] = TMath::Sqrt ((param0[min][0]-param0[min][5])*
 661       (param0[min][0]-param0[min][5])
 662       +(param0[min][1]-param0[min][6])*
 663       (param0[min][1]-param0[min][6]));
 664       rad = TMath::Sqrt ((param0[min][2]-param0[min][5])*
 665       (param0[min][2]-param0[min][5])
 666       +(param0[min][3]-param0[min][6])*
 667       (param0[min][3]-param0[min][6]));
 668       if (dist[2] < dist[0]) dist[0] = dist[2];
 669       if (rad < dist[1]) dist[1] = rad;
 670       if (rad < dist[2]) dist[2] = rad;
 671       }
 672       cout << dist[0] << " " << dist[1] << " " << dist[2] << endl;
 673       if (dist[TMath::LocMin(iseed+1,dist)] < 1.) { fNpar -= 3; break; }
 674       }
 675     */
 676
 677     for (Int_t i = 0; i < fNpar; ++i) {
 678       parOk[i] = param0[min][i];
 679       //errOk[i] = fmin;
 680       errOk[i] = chi2n;
 681       // Bounded params
 682       parOk[i] = TMath::Max (parOk[i], parmin[i]);
 683       parOk[i] = TMath::Min (parOk[i], parmax[i]);
 684     }
 685
 686     chi2o = chi2n;
 687     if (fmin < 0.1) break; // !!!???
 688   } // for (Int_t iseed=0;
 689
 690   if (fDebug) {
 691     for (Int_t i=0; i<fNpar; ++i) {
 692       if (i == 4 || i == 7) {
 693         if ((i == 7) || ((i == 4) && (fNpar < 7))) cout << parOk[i] << endl;
 694         else cout << parOk[i] * (1-parOk[7]) << endl;
 695         continue;
 696       }
 697       cout << parOk[i] << " " << errOk[i] << endl;
 698     }
 699   }
 700   nfit = (fNpar + 1) / 3;
 701   dist[0] = dist[1] = dist[2] = 0;
 702
 703   if (nfit > 1) {
 704     // Find distance to the nearest neighbour
 705     dist[0] = dist[1] = TMath::Sqrt ((parOk[0]-parOk[2])*
 706                                      (parOk[0]-parOk[2])
 707                                      +(parOk[1]-parOk[3])*
 708                                      (parOk[1]-parOk[3]));
 709     if (nfit > 2) {
 710       dist[2] = TMath::Sqrt ((parOk[0]-parOk[5])*
 711                              (parOk[0]-parOk[5])
 712                              +(parOk[1]-parOk[6])*
 713                              (parOk[1]-parOk[6]));
 714       rad = TMath::Sqrt ((parOk[2]-parOk[5])*
 715                          (parOk[2]-parOk[5])
 716                          +(parOk[3]-parOk[6])*
 717                          (parOk[3]-parOk[6]));
 718       if (dist[2] < dist[0]) dist[0] = dist[2];
 719       if (rad < dist[1]) dist[1] = rad;
 720       if (rad < dist[2]) dist[2] = rad;
 721     }
 722   }
 723
 724   Int_t indx;
 725
 726   Double_t coef = 0;
 727   if (iSimple) fnCoupled = 0;
 728   for (Int_t j = 0; j < nfit; ++j) {
 729     indx = 3 * j;
 730     coef = Param2Coef(j, coef, parOk);
 731
 732     //void AliMUONClusterFinderMLEM::AddRawCluster(Double_t x, Double_t y,
 733     //                                             Double_t qTot, Double_t fmin,
 734     //                                             Int_t nfit, Int_t *tracks,
 735     //                                             Double_t /*sigx*/,
 736     //                                             Double_t /*sigy*/,
 737     //                                             Double_t /*dist*/)
 738
 739     if ( coef*fQtot >= 14 )
 740     {
 741       //AZ AliMUONCluster* cluster1 = new AliMUONCluster();
 742       AliMUONCluster* cluster1 = new AliMUONCluster(cluster);
 743
 744       cluster1->SetCharge(coef*fQtot,coef*fQtot);
 745       cluster1->SetPosition(TVector2(parOk[indx],parOk[indx+1]),TVector2(sigCand[0][0],sigCand[0][1]));
 746       cluster1->SetChi2(dist[TMath::LocMin(nfit,dist)]);
 747
 748       // FIXME: we miss some information in this cluster, as compared to
 749       // the original AddRawCluster code.
 750
 751       AliDebug(2,Form("Adding RawCluster detElemId %4d mult %2d charge %5d (xl,yl)=(%9.6g,%9.6g)",
 752                       fDetElemId,cluster1->Multiplicity(),(Int_t)cluster1->Charge(),
 753                       cluster1->Position().X(),cluster1->Position().Y()));
 754
 755       clusterList.Add(cluster1);
 756     }
 757     //      AddRawCluster (parOk[indx], // double x
 758     //                     parOk[indx+1], // double y
 759     //                     coef*qTot, // double charge
 760     //                     errOk[indx], // double fmin
 761     //                     nfit0+10*nfit+100*nMax+10000*fnCoupled, // int nfit
 762     //                     tracks, // int* tracks
 763     //                     sigCand[0][0], // double sigx
 764     //                     sigCand[0][1], // double sigy
 765     //                     dist[TMath::LocMin(nfit,dist)] // double dist
 766     //                     );
 767   }
 768   return nfit;
 769 }
 770
 771
 772 //_____________________________________________________________________________
 773 void
 774 AliMUONClusterSplitterMLEM::Split(const AliMUONCluster& cluster,
 775                                   TH2 *mlem, Double_t *coef,
 776                                   TObjArray& clusterList)
 777 {
 778   /// The main steering function to work with clusters of pixels in anode
 779   /// plane (find clusters, decouple them from each other, merge them (if
 780   /// necessary), pick up coupled pads, call the fitting function)
 781
 782   Int_t nx = mlem->GetNbinsX();
 783   Int_t ny = mlem->GetNbinsY();
 784   Int_t nPix = fPixArray->GetEntriesFast();
 785
 786   Double_t cont;
 787   Int_t nclust = 0, indx, indx1, nxy = ny * nx;
 788   Bool_t *used = new Bool_t[nxy];
 789
 790   for (Int_t j = 0; j < nxy; ++j) used[j] = kFALSE;
 791
 792   TObjArray *clusters[200]={0};
 793   TObjArray *pix;
 794
 795   // Find clusters of histogram bins (easier to work in 2-D space)
 796   for (Int_t i = 1; i <= ny; ++i)
 797   {
 798     for (Int_t j = 1; j <= nx; ++j)
 799     {
 800       indx = (i-1)*nx + j - 1;
 801       if (used[indx]) continue;
 802       cont = mlem->GetCellContent(j,i);
 803       if (cont < 0.5) continue;
 804       pix = new TObjArray(20);
 805       used[indx] = 1;
 806       pix->Add(BinToPix(mlem,j,i));
 807       AddBin(mlem, i, j, 0, used, pix); // recursive call
 808       if (nclust >= 200) AliFatal(" Too many clusters !!!");
 809       clusters[nclust++] = pix;
 810     } // for (Int_t j=1; j<=nx; j++) {
 811   } // for (Int_t i=1; i<=ny;
 812   if (fDebug) cout << nclust << endl;
 813   delete [] used;
 814
 815   // Compute couplings between clusters and clusters to pads
 816   Int_t npad = cluster.Multiplicity();
 817
 818   // Exclude pads with overflows
 819   /*
 820   for (Int_t j = 0; j < npad; ++j)
 821   {
 822     AliMUONPad* pad = cluster.Pad(j);
 823     if ( pad->IsSaturated() )
 824     {
 825       pad->SetStatus(-5);
 826     }
 827     else
 828     {
 829       pad->SetStatus(0);
 830     }
 831   }
 832   */
 833
 834   // Compute couplings of clusters to pads (including overflows)
 835   TMatrixD aijclupad(nclust,npad);
 836   aijclupad = 0;
 837   Int_t npxclu;
 838   for (Int_t iclust = 0; iclust < nclust; ++iclust)
 839   {
 840     pix = clusters[iclust];
 841     npxclu = pix->GetEntriesFast();
 842     for (Int_t i = 0; i < npxclu; ++i)
 843     {
 844       indx = fPixArray->IndexOf(pix->UncheckedAt(i));
 845       for (Int_t j = 0; j < npad; ++j)
 846       {
 847         //AliMUONPad* pad = cluster.Pad(j);
 848         //if ( pad->Status() < 0 && pad->Status() != -5) continue;
 849         if (coef[j*nPix+indx] < fgkCouplMin) continue;
 850         aijclupad(iclust,j) += coef[j*nPix+indx];
 851       }
 852     }
 853   }
 854
 855   // Compute couplings between clusters (exclude overflows)
 856   TMatrixD aijcluclu(nclust,nclust);
 857   aijcluclu = 0;
 858   for (Int_t iclust = 0; iclust < nclust; ++iclust)
 859   {
 860     for (Int_t j = 0; j < npad; ++j)
 861     {
 862       // Exclude overflows
 863       //if ( cluster.Pad(j)->Status() < 0) continue;
 864       if ( cluster.Pad(j)->IsSaturated()) continue;
 865       if (aijclupad(iclust,j) < fgkCouplMin) continue;
 866       for (Int_t iclust1=iclust+1; iclust1<nclust; iclust1++)
 867       {
 868         if (aijclupad(iclust1,j) < fgkCouplMin) continue;
 869         aijcluclu(iclust,iclust1) +=
 870           TMath::Sqrt (aijclupad(iclust,j)*aijclupad(iclust1,j));
 871       }
 872     }
 873   }
 874   for (Int_t iclust = 0; iclust < nclust; ++iclust)
 875   {
 876     for (Int_t iclust1 = iclust+1; iclust1 < nclust; ++iclust1)
 877     {
 878       aijcluclu(iclust1,iclust) = aijcluclu(iclust,iclust1);
 879     }
 880   }
 881
 882   if (fDebug && nclust > 1) aijcluclu.Print();
 883
 884   // Find groups of coupled clusters
 885   used = new Bool_t[nclust];
 886   for (Int_t j = 0; j < nclust; ++j) used[j] = kFALSE;
 887
 888   Int_t *clustNumb = new Int_t[nclust];
 889   Int_t nCoupled, nForFit, minGroup[3], clustFit[3], nfit = 0;
 890   //Double_t parOk[8];
 891   Double_t parOk[8] = {0}; //AZ
 892
 893   for (Int_t igroup = 0; igroup < nclust; ++igroup)
 894   {
 895     if (used[igroup]) continue;
 896     used[igroup] = kTRUE;
 897     clustNumb[0] = igroup;
 898     nCoupled = 1;
 899     // Find group of coupled clusters
 900     AddCluster(igroup, nclust, aijcluclu, used, clustNumb, nCoupled); // recursive
 901
 902     if (fDebug) {
 903       cout << " nCoupled: " << nCoupled << endl;
 904       for (Int_t i=0; i<nCoupled; ++i) cout << clustNumb[i] << " "; cout << endl;
 905     }
 906
 907     fnCoupled = nCoupled;
 908
 909     while (nCoupled > 0)
 910     {
 911       if (nCoupled < 4)
 912       {
 913         nForFit = nCoupled;
 914         for (Int_t i = 0; i < nCoupled; ++i) clustFit[i] = clustNumb[i];
 915       }
 916       else
 917       {
 918         // Too many coupled clusters to fit - try to decouple them
 919         // Find the lowest coupling of 1, 2, min(3,nLinks/2) pixels with
 920         // all the others in the group
 921         for (Int_t j = 0; j < 3; ++j) minGroup[j] = -1;
 922         Double_t coupl = MinGroupCoupl(nCoupled, clustNumb, aijcluclu, minGroup);
 923
 924         // Flag clusters for fit
 925         nForFit = 0;
 926         while (minGroup[nForFit] >= 0 && nForFit < 3)
 927         {
 928           if (fDebug) cout << clustNumb[minGroup[nForFit]] << " ";
 929           clustFit[nForFit] = clustNumb[minGroup[nForFit]];
 930           clustNumb[minGroup[nForFit]] -= 999;
 931           nForFit++;
 932         }
 933         if (fDebug) cout << " nForFit " << nForFit << " " << coupl << endl;
 934       } // else
 935
 936       // Select pads for fit.
 937       if (SelectPad(cluster,nCoupled, nForFit, clustNumb, clustFit, aijclupad) < 3 && nCoupled > 1)
 938       {
 939         // Deselect pads
 940         for (Int_t j = 0; j < npad; ++j)
 941         {
 942           AliMUONPad* pad = cluster.Pad(j);
 943           //if ( pad->Status()==1 ) pad->SetStatus(0);
 944           //if ( pad->Status()==-9) pad->SetStatus(-5);
 945           if ( pad->Status() == AliMUONClusterFinderMLEM::GetUseForFitFlag() ||
 946                pad->Status() == AliMUONClusterFinderMLEM::GetCoupledFlag())
 947             pad->SetStatus(AliMUONClusterFinderMLEM::GetZeroFlag());
 948         }
 949         // Merge the failed cluster candidates (with too few pads to fit) with
 950         // the one with the strongest coupling
 951         Merge(cluster,nForFit, nCoupled, clustNumb, clustFit, clusters, aijcluclu, aijclupad);
 952       }
 953       else
 954       {
 955         // Do the fit
 956         nfit = Fit(cluster,0, nForFit, clustFit, clusters, parOk, clusterList, mlem);
 957         if (nfit == 0) {
 958           //cout << " (nfit == 0) " << fNpar << " " << cluster.Multiplicity() << endl;
 959           fNpar = 0; // should be 0 by itself but just in case ...
 960         }
 961       }
 962
 963       // Subtract the fitted charges from pads with strong coupling and/or
 964       // return pads for further use
 965       UpdatePads(cluster,nfit, parOk);
 966
 967       // Mark used pads
 968       for (Int_t j = 0; j < npad; ++j)
 969       {
 970         AliMUONPad* pad = cluster.Pad(j);
 971         //if ( pad->Status()==1 ) pad->SetStatus(-2);
 972         //if ( pad->Status()==-9) pad->SetStatus(-5);
 973         if ( pad->Status() == AliMUONClusterFinderMLEM::GetUseForFitFlag() )
 974           pad->SetStatus(AliMUONClusterFinderMLEM::GetModifiedFlag());
 975       }
 976
 977       // Sort the clusters (move to the right the used ones)
 978       Int_t beg = 0, end = nCoupled - 1;
 979       while (beg < end)
 980       {
 981         if (clustNumb[beg] >= 0) { ++beg; continue; }
 982         for (Int_t j = end; j > beg; --j)
 983         {
 984           if (clustNumb[j] < 0) continue;
 985           end = j - 1;
 986           indx = clustNumb[beg];
 987           clustNumb[beg] = clustNumb[j];
 988           clustNumb[j] = indx;
 989           break;
 990         }
 991         ++beg;
 992       }
 993
 994       nCoupled -= nForFit;
 995       if (nCoupled > 3)
 996       {
 997         // Remove couplings of used clusters
 998         for (Int_t iclust = nCoupled; iclust < nCoupled+nForFit; ++iclust)
 999         {
1000           indx = clustNumb[iclust] + 999;
1001           for (Int_t iclust1 = 0; iclust1 < nCoupled; ++iclust1)
1002           {
1003             indx1 = clustNumb[iclust1];
1004             aijcluclu(indx,indx1) = aijcluclu(indx1,indx) = 0;
1005           }
1006         }
1007
1008         // Update the remaining clusters couplings (subtract couplings from
1009         // the used pads) - overflows excluded
1010         for (Int_t j = 0; j < npad; ++j)
1011         {
1012           AliMUONPad* pad = cluster.Pad(j);
1013           //if ( pad->Status() != -2) continue;
1014           if ( pad->Status() != AliMUONClusterFinderMLEM::GetModifiedFlag()) continue;
1015           for (Int_t iclust=0; iclust<nCoupled; ++iclust)
1016           {
1017             indx = clustNumb[iclust];
1018             if (aijclupad(indx,j) < fgkCouplMin) continue;
1019             for (Int_t iclust1 = iclust+1; iclust1 < nCoupled; ++iclust1)
1020             {
1021               indx1 = clustNumb[iclust1];
1022               if (aijclupad(indx1,j) < fgkCouplMin) continue;
1023               // Check this
1024               aijcluclu(indx,indx1) -=
1025                 TMath::Sqrt (aijclupad(indx,j)*aijclupad(indx1,j));
1026               aijcluclu(indx1,indx) = aijcluclu(indx,indx1);
1027             }
1028           }
1029           //pad->SetStatus(-8);
1030           pad->SetStatus(AliMUONClusterFinderMLEM::GetOverFlag());
1031         } // for (Int_t j=0; j<npad;
1032       } // if (nCoupled > 3)
1033     } // while (nCoupled > 0)
1034   } // for (Int_t igroup=0; igroup<nclust;
1035
1036   for (Int_t iclust = 0; iclust < nclust; ++iclust)
1037   {
1038     pix = clusters[iclust];
1039     pix->Clear();
1040     delete pix;
1041   }
1042   delete [] clustNumb;
1043   delete [] used;
1044
1045 }
1046
1047 //_____________________________________________________________________________
1048 void
1049 AliMUONClusterSplitterMLEM::Merge(const AliMUONCluster& cluster,
1050                                      Int_t nForFit, Int_t nCoupled,
1051                                      Int_t *clustNumb, Int_t *clustFit,
1052                                      TObjArray **clusters,
1053                                      TMatrixD& aijcluclu, TMatrixD& aijclupad)
1054 {
1055   /// Merge the group of clusters with the one having the strongest coupling with them
1056
1057   Int_t indx, indx1, npxclu, npxclu1, imax=0;
1058   TObjArray *pix, *pix1;
1059   Double_t couplMax;
1060
1061   for (Int_t icl = 0; icl < nForFit; ++icl)
1062   {
1063     indx = clustFit[icl];
1064     pix = clusters[indx];
1065     npxclu = pix->GetEntriesFast();
1066     couplMax = -1;
1067     for (Int_t icl1 = 0; icl1 < nCoupled; ++icl1)
1068     {
1069       indx1 = clustNumb[icl1];
1070       if (indx1 < 0) continue;
1071       if ( aijcluclu(indx,indx1) > couplMax)
1072       {
1073         couplMax = aijcluclu(indx,indx1);
1074         imax = indx1;
1075       }
1076     } // for (Int_t icl1=0;
1077       // Add to it
1078     pix1 = clusters[imax];
1079     npxclu1 = pix1->GetEntriesFast();
1080     // Add pixels
1081     for (Int_t i = 0; i < npxclu; ++i)
1082     {
1083       pix1->Add(pix->UncheckedAt(i));
1084       pix->RemoveAt(i);
1085     }
1086
1087     //Add cluster-to-cluster couplings
1088     for (Int_t icl1 = 0; icl1 < nCoupled; ++icl1)
1089     {
1090       indx1 = clustNumb[icl1];
1091       if (indx1 < 0 || indx1 == imax) continue;
1092       aijcluclu(indx1,imax) += aijcluclu(indx,indx1);
1093       aijcluclu(imax,indx1) = aijcluclu(indx1,imax);
1094     }
1095     aijcluclu(indx,imax) = aijcluclu(imax,indx) = 0;
1096
1097     //Add cluster-to-pad couplings
1098     Int_t mult = cluster.Multiplicity();
1099     for (Int_t j = 0; j < mult; ++j)
1100     {
1101       AliMUONPad* pad = cluster.Pad(j);
1102       //if ( pad->Status() < 0 && pad->Status() != -5 ) continue;// exclude used pads
1103       if ( pad->Status() != AliMUONClusterFinderMLEM::GetZeroFlag()) continue;// exclude used pads
1104         aijclupad(imax,j) += aijclupad(indx,j);
1105         aijclupad(indx,j) = 0;
1106     }
1107   } // for (Int_t icl=0; icl<nForFit;
1108 }
1109
1110
1111 //_____________________________________________________________________________
1112 Double_t
1113 AliMUONClusterSplitterMLEM::MinGroupCoupl(Int_t nCoupled, Int_t *clustNumb,
1114                                           TMatrixD& aijcluclu, Int_t *minGroup)
1115 {
1116   /// Find group of clusters with minimum coupling to all the others
1117
1118   Int_t i123max = TMath::Min(3,nCoupled/2);
1119   Int_t indx, indx1, indx2, indx3, nTot = 0;
1120   Double_t *coupl1 = 0, *coupl2 = 0, *coupl3 = 0;
1121
1122   for (Int_t i123 = 1; i123 <= i123max; ++i123) {
1123
1124     if (i123 == 1) {
1125       coupl1 = new Double_t [nCoupled];
1126       for (Int_t i = 0; i < nCoupled; ++i) coupl1[i] = 0;
1127     }
1128     else if (i123 == 2) {
1129       nTot = nCoupled*nCoupled;
1130       coupl2 = new Double_t [nTot];
1131       for (Int_t i = 0; i < nTot; ++i) coupl2[i] = 9999;
1132     } else {
1133       nTot = nTot*nCoupled;
1134       coupl3 = new Double_t [nTot];
1135       for (Int_t i = 0; i < nTot; ++i) coupl3[i] = 9999;
1136     } // else
1137
1138     for (Int_t i = 0; i < nCoupled; ++i) {
1139       indx1 = clustNumb[i];
1140       for (Int_t j = i+1; j < nCoupled; ++j) {
1141         indx2 = clustNumb[j];
1142         if (i123 == 1) {
1143           coupl1[i] += aijcluclu(indx1,indx2);
1144           coupl1[j] += aijcluclu(indx1,indx2);
1145         }
1146         else if (i123 == 2) {
1147           indx = i*nCoupled + j;
1148           coupl2[indx] = coupl1[i] + coupl1[j];
1149           coupl2[indx] -= 2 * (aijcluclu(indx1,indx2));
1150         } else {
1151           for (Int_t k = j+1; k < nCoupled; ++k) {
1152             indx3 = clustNumb[k];
1153             indx = i*nCoupled*nCoupled + j*nCoupled + k;
1154             coupl3[indx] = coupl2[i*nCoupled+j] + coupl1[k];
1155             coupl3[indx] -= 2 * (aijcluclu(indx1,indx3)+aijcluclu(indx2,indx3));
1156           }
1157         } // else
1158       } // for (Int_t j=i+1;
1159     } // for (Int_t i=0;
1160   } // for (Int_t i123=1;
1161
1162   // Find minimum coupling
1163   Double_t couplMin = 9999;
1164   Int_t locMin = 0;
1165
1166   for (Int_t i123 = 1; i123 <= i123max; ++i123) {
1167     if (i123 == 1) {
1168       locMin = TMath::LocMin(nCoupled, coupl1);
1169       couplMin = coupl1[locMin];
1170       minGroup[0] = locMin;
1171       delete [] coupl1;
1172     }
1173     else if (i123 == 2) {
1174       locMin = TMath::LocMin(nCoupled*nCoupled, coupl2);
1175       if (coupl2[locMin] < couplMin) {
1176         couplMin = coupl2[locMin];
1177         minGroup[0] = locMin/nCoupled;
1178         minGroup[1] = locMin%nCoupled;
1179       }
1180       delete [] coupl2;
1181     } else {
1182       locMin = TMath::LocMin(nTot, coupl3);
1183       if (coupl3[locMin] < couplMin) {
1184         couplMin = coupl3[locMin];
1185         minGroup[0] = locMin/nCoupled/nCoupled;
1186         minGroup[1] = locMin%(nCoupled*nCoupled)/nCoupled;
1187         minGroup[2] = locMin%nCoupled;
1188       }
1189       delete [] coupl3;
1190     } // else
1191   } // for (Int_t i123=1;
1192   return couplMin;
1193 }
1194
1195 //_____________________________________________________________________________
1196 Int_t
1197 AliMUONClusterSplitterMLEM::SelectPad(const AliMUONCluster& cluster,
1198                                           Int_t nCoupled, Int_t nForFit,
1199                                           Int_t *clustNumb, Int_t *clustFit,
1200                                           TMatrixD& aijclupad)
1201 {
1202   /// Select pads for fit. If too many coupled clusters, find pads giving
1203   /// the strongest coupling with the rest of clusters and exclude them from the fit.
1204
1205   Int_t npad = cluster.Multiplicity();
1206   Double_t *padpix = 0;
1207
1208   if (nCoupled > 3)
1209   {
1210     padpix = new Double_t[npad];
1211     for (Int_t i = 0; i < npad; ++i) padpix[i] = 0.;
1212   }
1213
1214   Int_t nOK = 0, indx, indx1;
1215   for (Int_t iclust = 0; iclust < nForFit; ++iclust)
1216   {
1217     indx = clustFit[iclust];
1218     for (Int_t j = 0; j < npad; ++j)
1219     {
1220       if ( aijclupad(indx,j) < fgkCouplMin) continue;
1221       AliMUONPad* pad = cluster.Pad(j);
1222       /*
1223       if ( pad->Status() == -5 ) pad->SetStatus(-9); // flag overflow
1224       if ( pad->Status() < 0 ) continue; // exclude overflows and used pads
1225       if ( !pad->Status() )
1226       {
1227         pad->SetStatus(1);
1228         ++nOK; // pad to be used in fit
1229       }
1230       */
1231       if ( pad->Status() != AliMUONClusterFinderMLEM::GetZeroFlag()
1232            || pad->IsSaturated() ) continue; // used pads and overflows
1233       pad->SetStatus(AliMUONClusterFinderMLEM::GetUseForFitFlag());
1234       ++nOK; // pad to be used in fit
1235
1236       if (nCoupled > 3)
1237       {
1238         // Check other clusters
1239         for (Int_t iclust1 = 0; iclust1 < nCoupled; ++iclust1)
1240         {
1241           indx1 = clustNumb[iclust1];
1242           if (indx1 < 0) continue;
1243           if ( aijclupad(indx1,j) < fgkCouplMin ) continue;
1244           padpix[j] += aijclupad(indx1,j);
1245         }
1246       } // if (nCoupled > 3)
1247     } // for (Int_t j=0; j<npad;
1248   } // for (Int_t iclust=0; iclust<nForFit
1249   if (nCoupled < 4) return nOK;
1250
1251   Double_t aaa = 0;
1252   for (Int_t j = 0; j < npad; ++j)
1253   {
1254     if (padpix[j] < fgkCouplMin) continue;
1255     aaa += padpix[j];
1256     //cluster.Pad(j)->SetStatus(-1); // exclude pads with strong coupling to the other clusters
1257     cluster.Pad(j)->SetStatus(AliMUONClusterFinderMLEM::GetCoupledFlag()); // exclude pads with strong coupling to the other clusters
1258     nOK--;
1259   }
1260   delete [] padpix;
1261   return nOK;
1262 }
1263
1264 //_____________________________________________________________________________
1265 void
1266 AliMUONClusterSplitterMLEM::UpdatePads(const AliMUONCluster& cluster,
1267                                           Int_t /*nfit*/, Double_t *par)
1268 {
1269   /// Subtract the fitted charges from pads with strong coupling
1270
1271   Int_t indx, mult = cluster.Multiplicity(), iend = fNpar/3;
1272   Double_t charge, coef=0;
1273
1274   for (Int_t j = 0; j < mult; ++j)
1275   {
1276     AliMUONPad* pad = cluster.Pad(j);
1277     //if ( pad->Status() != -1 ) continue;
1278     if ( pad->Status() != AliMUONClusterFinderMLEM::GetCoupledFlag() ) continue;
1279     if (fNpar != 0)
1280     {
1281       charge = 0;
1282       for (Int_t i = 0; i <= iend; ++i)
1283       {
1284         // sum over hits
1285         indx = 3 * i;
1286         coef = Param2Coef(i, coef, par);
1287         charge += ChargeIntegration(par[indx],par[indx+1],*pad) * coef;
1288       }
1289       charge *= fQtot;
1290       pad->SetCharge(pad->Charge()-charge);
1291     } // if (fNpar != 0)
1292
1293     //if (pad->Charge() > 6 /*fgkZeroSuppression*/) pad->SetStatus(0);
1294     if (pad->Charge() > 6 /*fgkZeroSuppression*/) pad->SetStatus(AliMUONClusterFinderMLEM::GetZeroFlag());
1295     // return pad for further using // FIXME: remove usage of zerosuppression here
1296     else pad->SetStatus(AliMUONClusterFinderMLEM::GetOverFlag()); // do not use anymore
1297
1298   } // for (Int_t j=0;
1299 }
1300
1301