MUON/AliMUONClusterFinderAZ.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 AliMUONClusterFinderAZ
  20 // -------------------------------
  21 // Clusterizer class based on the Expectation-Maximization algorithm
  22 // Author: Alexander Zinchenko, JINR Dubna
  23
  24 #include <stdlib.h>
  25 #include <Riostream.h>
  26 #include <TH2.h>
  27 #include <TMinuit.h>
  28 #include <TMatrixD.h>
  29 #include <TRandom.h>
  30 #include <TROOT.h>
  31 #include <TMath.h>
  32
  33 #include "AliMUONClusterFinderAZ.h"
  34 #include "AliMpVSegmentation.h"
  35 #include "AliMUONGeometryModuleTransformer.h"
  36 #include "AliMUONVDigit.h"
  37 #include "AliMUONCluster.h"
  38 #include "AliMUONPixel.h"
  39 #include "AliMUONMathieson.h"
  40 #include "AliLog.h"
  41 #include <TClonesArray.h>
  42 #include "AliMpDEManager.h"
  43 #include "AliMUONVDigitStore.h"
  44 #include "AliMUONConstants.h"
  45
  46 /// \cond CLASSIMP
  47 ClassImp(AliMUONClusterFinderAZ)
  48 /// \endcond
  49
  50  const Double_t AliMUONClusterFinderAZ::fgkCouplMin = 1.e-3; // threshold on coupling
  51  const Double_t AliMUONClusterFinderAZ::fgkZeroSuppression = 6; // average zero suppression value
  52  const Double_t AliMUONClusterFinderAZ::fgkSaturation = 3000; // average saturation level
  53  AliMUONClusterFinderAZ* AliMUONClusterFinderAZ::fgClusterFinder = 0x0;
  54  TMinuit* AliMUONClusterFinderAZ::fgMinuit = 0x0;
  55 //FILE *lun1 = fopen("nxny.dat","w");
  56
  57 //_____________________________________________________________________________
  58 AliMUONClusterFinderAZ::AliMUONClusterFinderAZ(Bool_t draw)
  59   : AliMUONVClusterFinder(),
  60     fNpar(0),
  61     fQtot(0),
  62     fReco(1),
  63     fCathBeg(0),
  64 //    fDraw(0x0),
  65     fPixArray(0x0),
  66     fnCoupled(0),
  67     fDebug(0),
  68 fRawClusters(new TClonesArray("AliMUONCluster",100)),
  69 fDigitStore(0x0),
  70 fDetElemId(-1),
  71 fChamberId(-1),
  72 fMathieson(0x0),
  73 fCurrentCluster(-1)
  74 {
  75 /// Constructor
  76   fnPads[0]=fnPads[1]=0;
  77
  78   for (Int_t i=0; i<7; i++)
  79     for (Int_t j=0; j<fgkDim; j++)
  80       fXyq[i][j]= 9999.;
  81
  82   for (Int_t i=0; i<4; i++)
  83     for (Int_t j=0; j<fgkDim; j++)
  84       fPadIJ[i][j]=-1;
  85
  86   for (Int_t i=0; i<2; i++)
  87     for (Int_t j=0; j<fgkDim; j++)
  88       fUsed[i][j] = 0;
  89
  90   fSegmentation[1] = fSegmentation[0] = 0x0;
  91
  92   fPadBeg[0] = fPadBeg[1] = 0;
  93
  94   if (!fgMinuit) fgMinuit = new TMinuit(8);
  95   if (!fgClusterFinder) fgClusterFinder = this;
  96   fPixArray = new TObjArray(20);
  97
  98   if (draw) {
  99     fDebug = 1;
 100     fReco = 0;
 101     //    fDraw = new AliMUONClusterDrawAZ(this);
 102   }
 103   AliInfo(" *** Running AZ cluster finder *** ");
 104 }
 105
 106 //_____________________________________________________________________________
 107 AliMUONClusterFinderAZ::~AliMUONClusterFinderAZ()
 108 {
 109 /// Destructor
 110   delete fgMinuit; fgMinuit = 0; delete fPixArray; fPixArray = 0;
 111   delete fMathieson;
 112 //  delete fDraw;
 113 }
 114
 115 //_____________________________________________________________________________
 116 Bool_t
 117 AliMUONClusterFinderAZ::Prepare(const AliMpVSegmentation* segmentations[2],
 118                                 const AliMUONVDigitStore& digitStore)
 119 {
 120   /// Prepare for the clusterization of one detection element, which digits
 121   /// are in digitStore
 122
 123   fSegmentation[0] = segmentations[0];
 124   fSegmentation[1] = segmentations[1];
 125   fDigitStore = &digitStore;
 126   fDetElemId = -1;
 127   TIter next(digitStore.CreateIterator());
 128   AliMUONVDigit* digit = static_cast<AliMUONVDigit*>(next());
 129   if (digit) fDetElemId = digit->DetElemId();
 130   fCurrentCluster = -1;
 131   if ( fDetElemId > 0 )
 132   {
 133     fChamberId = AliMpDEManager::GetChamberId(fDetElemId);
 134     AliMp::StationType stationType = AliMpDEManager::GetStationType(fDetElemId);
 135
 136     Float_t kx3 = AliMUONConstants::SqrtKx3();
 137     Float_t ky3 = AliMUONConstants::SqrtKy3();
 138     Float_t pitch = AliMUONConstants::Pitch();
 139
 140     if ( stationType == AliMp::kStation1 )
 141     {
 142       kx3 = AliMUONConstants::SqrtKx3St1();
 143       ky3 = AliMUONConstants::SqrtKy3St1();
 144       pitch = AliMUONConstants::PitchSt1();
 145     }
 146
 147     delete fMathieson;
 148     fMathieson = new AliMUONMathieson;
 149
 150     fMathieson->SetPitch(pitch);
 151     fMathieson->SetSqrtKx3AndDeriveKx2Kx4(kx3);
 152     fMathieson->SetSqrtKy3AndDeriveKy2Ky4(ky3);
 153
 154     return kTRUE;
 155   }
 156   return kFALSE;
 157 }
 158
 159 //_____________________________________________________________________________
 160 AliMUONCluster*
 161 AliMUONClusterFinderAZ::NextCluster()
 162 {
 163   /// Return the next cluster in the iteration
 164   if ( fCurrentCluster == -1 )
 165   {
 166     FindRawClusters(0);
 167   }
 168
 169   ++fCurrentCluster;
 170   if ( fCurrentCluster <= fRawClusters->GetLast() )
 171   {
 172     return static_cast<AliMUONCluster*>(fRawClusters->At(fCurrentCluster));
 173   }
 174   return 0x0;
 175 }
 176
 177 //_____________________________________________________________________________
 178 void AliMUONClusterFinderAZ::FindRawClusters(Int_t ch)
 179 {
 180   /// To comply with old old old interface...
 181   ResetRawClusters();
 182   EventLoop (ch);
 183 }
 184
 185 //_____________________________________________________________________________
 186 void AliMUONClusterFinderAZ::EventLoop(Int_t)
 187 {
 188 /// Loop over digits
 189
 190  // if (fDraw && !fDraw->FindEvCh(nev, ch)) return;
 191
 192 //  AliInfo("");
 193 //  fDigitStore->Print();
 194
 195   Int_t ndigits[2] = {9,9}, nShown[2] = {0};
 196   if (fReco != 2) { // skip initialization for the combined cluster / track
 197     fCathBeg = fPadBeg[0] = fPadBeg[1] = 0;
 198     for (Int_t i = 0; i < 2; i++) {
 199       for (Int_t j = 0; j < fgkDim; j++) { fUsed[i][j] = kFALSE; }
 200     }
 201   }
 202
 203 next:
 204   if (fReco == 2 && (nShown[0] || nShown[1])) return; // only one precluster for the combined finder
 205   if (ndigits[0] == nShown[0] && ndigits[1] == nShown[1]) return;
 206
 207   Bool_t first = kTRUE;
 208   fnPads[0] = fnPads[1] = 0;
 209   for (Int_t i = 0; i < fgkDim; i++)
 210   {
 211     fPadIJ[1][i] = 0;
 212     fDigitId[i] = 0;
 213   }
 214
 215   for (Int_t iii = fCathBeg; iii < 2; iii++)
 216   {
 217     Int_t cath = TMath::Odd(iii);
 218     TIter next(fDigitStore->CreateIterator(fDetElemId,fDetElemId,cath));
 219
 220     AliMUONVDigit  *mdig;
 221     Bool_t eEOC = kTRUE; // end-of-cluster
 222
 223     while ( ( mdig = static_cast<AliMUONVDigit*>(next()) ) )
 224     {
 225       if (first)
 226       {
 227         // Find first unused pad
 228         if (mdig->IsUsed()) continue;
 229       }
 230       else
 231       {
 232         if (mdig->IsUsed()) continue;
 233         // Find a pad overlapping with the cluster
 234         if (!Overlap(cath,*mdig)) continue;
 235       }
 236       // Add pad - recursive call
 237       AddPad(cath,*mdig);
 238       eEOC = kFALSE;
 239       break;
 240     }
 241     if (first && eEOC)
 242     {
 243       // No more unused pads
 244       if (cath == 0) continue; // on cathode #0 - check #1
 245       else return; // No more clusters
 246     }
 247     if (eEOC) break; // cluster found
 248     first = kFALSE;
 249     if (fDebug) cout << " nPads: " << fnPads[cath] << " " << nShown[cath]+fnPads[cath] << " " << cath << endl;
 250   } // for (Int_t iii = 0;
 251
 252 //  if (fDraw) fDraw->DrawCluster();
 253
 254   // Use MLEM for cluster finder
 255   Int_t nMax = 1, localMax[100], maxPos[100];
 256   Double_t maxVal[100];
 257
 258   if (CheckPrecluster(nShown)) {
 259 //    AliInfo("After CheckPrecluster");
 260 //    Print();
 261     BuildPixArray();
 262 //    AliInfo("PixArray");
 263 //    fPixArray->Print();
 264     //*
 265     if (fnPads[0]+fnPads[1] > 50) nMax = FindLocalMaxima(fPixArray, localMax, maxVal);
 266     if (nMax > 1) TMath::Sort(nMax, maxVal, maxPos, kTRUE); // in decreasing order
 267     Int_t iSimple = 0, nInX = -1, nInY;
 268     PadsInXandY(nInX, nInY);
 269     if (fDebug) cout << "Pads in X and Y: " << nInX << " " << nInY << endl;
 270     if (nMax == 1 && nInX < 4 && nInY < 4) iSimple = 1; //1; // simple cluster
 271     //*/
 272     /* For test
 273     Int_t iSimple = 0, nInX = -1, nInY;
 274     PadsInXandY(nInX, nInY);
 275     if (fDebug) cout << "Pads in X and Y: " << nInX << " " << nInY << endl;
 276     if (nMax == 1 && nInX < 4 && nInY < 4) iSimple = 1; //1; // simple cluster
 277     if (!iSimple) nMax = FindLocalMaxima(fPixArray, localMax, maxVal);
 278     nMax = 1;
 279     if (nMax > 1) TMath::Sort(nMax, maxVal, maxPos, kTRUE); // in decreasing order
 280     */
 281     for (Int_t i=0; i<nMax; i++) {
 282       if (nMax > 1) FindCluster(localMax, maxPos[i]);
 283       MainLoop(iSimple);
 284       if (i < nMax-1) {
 285         for (Int_t j=0; j<fnPads[0]+fnPads[1]; j++) {
 286           if (fPadIJ[1][j] == 0) continue; // pad charge was not modified
 287           fPadIJ[1][j] = 0;
 288           fXyq[2][j] = fXyq[6][j]; // use backup charge value
 289         }
 290       }
 291     } // for (Int_t i=0; i<nMax;
 292     if (nMax > 1) ((TH2D*) gROOT->FindObject("anode"))->Delete();
 293     TH2D *mlem = (TH2D*) gROOT->FindObject("mlem");
 294     if (mlem) mlem->Delete();
 295   }
 296 //  if (!fDraw || fDraw->Next()) goto next;
 297      goto next;
 298 }
 299
 300 //_____________________________________________________________________________
 301 void AliMUONClusterFinderAZ::AddPad(Int_t cath, AliMUONVDigit& mdig)
 302 {
 303 /// Add pad to the cluster
 304
 305 //  AliInfo("");
 306 //  StdoutToAliWarning(mdig.Print(););
 307
 308   Float_t charge = mdig.Charge();
 309
 310   AliMpPad pad = fSegmentation[cath]->PadByIndices(AliMpIntPair(mdig.PadX(),mdig.PadY()));
 311
 312   // get the center of the pad
 313   Float_t xpad = pad.Position().X();
 314   Float_t ypad = pad.Position().Y();
 315
 316 //  Int_t isec = fSegmentation[cath]->Sector(mdig.PadX(), mdig.PadY());
 317   Int_t nPads = fnPads[0] + fnPads[1];
 318   fXyq[0][nPads] = xpad;
 319   fXyq[1][nPads] = ypad;
 320   fXyq[2][nPads] = charge;
 321   fXyq[3][nPads] = pad.Dimensions().X();
 322   fXyq[4][nPads] = pad.Dimensions().Y();
 323   fXyq[5][nPads] = -1;
 324   fDigitId[nPads] = mdig.GetUniqueID();
 325   fXyq[6][nPads] = 0;
 326   fPadIJ[0][nPads] = cath;
 327   fPadIJ[1][nPads] = 0;
 328   fPadIJ[2][nPads] = mdig.PadX();
 329   fPadIJ[3][nPads] = mdig.PadY();
 330   mdig.Used(kTRUE);
 331   if (fDebug) printf(" bbb %d %d %f %f %f %f %f %3d %3d \n", nPads, cath,
 332                      xpad, ypad, fXyq[3][nPads]*2, fXyq[4][nPads]*2,
 333                      charge, mdig.PadX(), mdig.PadY());
 334   fnPads[cath]++;
 335
 336   // Check neighbours
 337   TObjArray neighbours;
 338   Int_t nn = fSegmentation[cath]->GetNeighbours(pad,neighbours);
 339   for (Int_t in = 0; in < nn; ++in)
 340   {
 341     AliMpPad* p = static_cast<AliMpPad*>(neighbours.At(in));
 342     AliMUONVDigit* mdig1 = static_cast<AliMUONVDigit*>
 343       (fDigitStore->FindObject(fDetElemId,p->GetLocation().GetFirst(),p->GetLocation().GetSecond()));
 344     if ( mdig1 && !mdig1->IsUsed() )
 345     {
 346       AddPad(cath,*mdig1);
 347     }
 348   } // for (Int_t in = 0;
 349 }
 350
 351 //_____________________________________________________________________________
 352 Bool_t AliMUONClusterFinderAZ::Overlap(Int_t cath, const AliMUONVDigit& mdig)
 353 {
 354 /// Check if the pad from one cathode overlaps with a pad
 355 /// in the precluster on the other cathode
 356
 357   AliMpPad pad = fSegmentation[cath]->PadByIndices(AliMpIntPair(mdig.PadX(), mdig.PadY()));
 358
 359   Float_t xpad = pad.Position().X();
 360   Float_t ypad = pad.Position().Y();
 361
 362   Float_t dx = pad.Dimensions().X();
 363   Float_t dy = pad.Dimensions().Y();
 364
 365   Float_t xy1[4], xy12[4];
 366   xy1[0] = xpad - dx;
 367   xy1[1] = xy1[0] + dx*2;
 368   xy1[2] = ypad - dy;
 369   xy1[3] = xy1[2] + dy*2;
 370
 371   Int_t cath1 = TMath::Even(cath);
 372   for (Int_t i=0; i<fnPads[0]+fnPads[1]; i++) {
 373     if (fPadIJ[0][i] != cath1) continue;
 374     if (Overlap(xy1, i, xy12, 0)) return kTRUE;
 375   }
 376   return kFALSE;
 377 }
 378
 379 //_____________________________________________________________________________
 380 Bool_t AliMUONClusterFinderAZ::Overlap(Float_t *xy1, Int_t iPad, Float_t *xy12, Int_t iSkip)
 381 {
 382 /// Check if the pads xy1 and iPad overlap and return overlap area
 383
 384   Float_t xy2[4];
 385   xy2[0] = fXyq[0][iPad] - fXyq[3][iPad];
 386   xy2[1] = fXyq[0][iPad] + fXyq[3][iPad];
 387   if (xy1[0] > xy2[1]-1.e-4 || xy1[1] < xy2[0]+1.e-4) return kFALSE;
 388   xy2[2] = fXyq[1][iPad] - fXyq[4][iPad];
 389   xy2[3] = fXyq[1][iPad] + fXyq[4][iPad];
 390   if (xy1[2] > xy2[3]-1.e-4 || xy1[3] < xy2[2]+1.e-4) return kFALSE;
 391   if (!iSkip) return kTRUE; // just check overlap (w/out computing the area)
 392   xy12[0] = TMath::Max (xy1[0],xy2[0]);
 393   xy12[1] = TMath::Min (xy1[1],xy2[1]);
 394   xy12[2] = TMath::Max (xy1[2],xy2[2]);
 395   xy12[3] = TMath::Min (xy1[3],xy2[3]);
 396   return kTRUE;
 397 }
 398
 399 //_____________________________________________________________________________
 400 void
 401 AliMUONClusterFinderAZ::Used(Int_t indx, Bool_t value)
 402 {
 403   /// Change the Used status of the pad at index indx
 404   AliMUONVDigit* digit = static_cast<AliMUONVDigit*>
 405   (fDigitStore->FindObject(fDigitId[indx]));
 406   if (!digit)
 407   {
 408     AliError(Form("Did not find digit %d",fDigitId[indx]));
 409   }
 410   else
 411   {
 412     digit->Used(value);
 413   }
 414 }
 415
 416 //_____________________________________________________________________________
 417 void
 418 AliMUONClusterFinderAZ::PrintPixel(Int_t i) const
 419 {
 420   /// Printout one pixel
 421   AliMUONPixel* pixel = static_cast<AliMUONPixel*>(fPixArray->UncheckedAt(i));
 422   if (pixel) pixel->Print();
 423 }
 424
 425 //_____________________________________________________________________________
 426 void
 427 AliMUONClusterFinderAZ::PrintPad(Int_t i) const
 428 {
 429   /// Printout one pad
 430   Int_t cathode = fPadIJ[0][i];
 431   UInt_t index = fDigitId[i];
 432   Int_t ix = fPadIJ[2][i];
 433   Int_t iy = fPadIJ[3][i];
 434
 435   cout << Form("i=%4d status %1d cathode %1d index %u ix %3d iy %3d (x,y)=(%7.2f,%7.2f) (dx,dy)=(%7.2f,%7.2f) Q=%7.2f",
 436                i,fPadIJ[1][i],cathode,index,ix,iy,fXyq[0][i],fXyq[1][i],
 437                fXyq[3][i],fXyq[4][i],
 438                fXyq[2][i]) << endl;
 439 }
 440
 441 //_____________________________________________________________________________
 442 void
 443 AliMUONClusterFinderAZ::Print(Option_t*) const
 444 {
 445   /// Print current state
 446   Int_t nPads = fnPads[0] + fnPads[1];
 447   cout << "PreCluster npads=" << nPads << "(" << fnPads[0] << ","
 448     << fnPads[1] << ")" << endl;
 449   for ( Int_t i = 0; i < nPads; ++i )
 450   {
 451     PrintPad(i);
 452   }
 453 }
 454
 455 //_____________________________________________________________________________
 456 Bool_t AliMUONClusterFinderAZ::CheckPrecluster(Int_t *nShown)
 457 {
 458 /// Check precluster in order to attempt to simplify it (mostly for
 459 /// two-cathode preclusters)
 460
 461 //  AliInfo("CheckPrecluster");
 462 //  Print();
 463
 464   Int_t i1, i2, cath=0;
 465   Float_t xy1[4], xy12[4];
 466
 467   Int_t npad = fnPads[0] + fnPads[1];
 468   if (npad == 1) {
 469     // Disregard one-pad clusters (leftovers from splitting)
 470     nShown[0] += fnPads[0];
 471     nShown[1] += fnPads[1];
 472     return kFALSE;
 473   }
 474
 475   // If pads have the same size take average of pads on both cathodes
 476   //Int_t sameSize = (fnPads[0] && fnPads[1]) ? 1 : 0;
 477   Int_t sameSize = 0; //AZ - 17-01-06
 478
 479   if (sameSize) {
 480     Double_t xSize = -1, ySize = 0;
 481     for (Int_t i=0; i<npad; i++) {
 482       if (fXyq[2][i] < 0) continue;
 483       if (xSize < 0) { xSize = fXyq[3][i]; ySize = fXyq[4][i]; }
 484       if (TMath::Abs(xSize-fXyq[3][i]) > 1.e-4 ||  TMath::Abs(ySize-fXyq[4][i]) > 1.e-4) { sameSize = 0; break; }
 485     }
 486   } // if (sameSize)
 487   if (sameSize && fnPads[0] == 1 && fnPads[1] == 1) sameSize = 0; //AZ
 488   // Handle shift by half a pad in Station 1
 489   if (sameSize) {
 490     Int_t cath0 = fPadIJ[0][0];
 491     for (Int_t i = 1; i < npad; i++) {
 492       if (fPadIJ[0][i] == cath0) continue;
 493       Double_t dx = TMath::Abs ((fXyq[0][i] - fXyq[0][0]) / fXyq[3][i] / 2);
 494       Int_t idx = (Int_t) TMath::Abs ((fXyq[0][i] - fXyq[0][0]) / fXyq[3][i] / 2);
 495       if (TMath::Abs (dx - idx) > 0.001) sameSize = 0;
 496       break;
 497     }
 498   } // if (sameSize)
 499
 500   if (sameSize && (fnPads[0] >= 2 || fnPads[1] >= 2)) {
 501     nShown[0] += fnPads[0];
 502     nShown[1] += fnPads[1];
 503     fnPads[0] = fnPads[1] = 0;
 504     Int_t div;
 505     for (Int_t i=0; i<npad; i++) {
 506       if (fXyq[2][i] < 0) continue; // used pad
 507       fXyq[2][fnPads[0]] = fXyq[2][i];
 508       div = 1;
 509       cath = fPadIJ[0][i];
 510       for (Int_t j=i+1; j<npad; j++) {
 511         if (fPadIJ[0][j] == fPadIJ[0][i]) continue; // same cathode
 512         if (TMath::Abs(fXyq[0][j]-fXyq[0][i]) > 1.e-4) continue;
 513         if (TMath::Abs(fXyq[1][j]-fXyq[1][i]) > 1.e-4) continue;
 514         fXyq[2][fnPads[0]] += fXyq[2][j];
 515         div = 2;
 516         fXyq[2][j] = -2;
 517         if (cath) fDigitId[fnPads[0]] = fDigitId[j]; // save digit number for cath 0
 518         break;
 519       }
 520       // Flag that the digit from the other cathode
 521       // LA commented  if (cath && div == 1) fXyq[5][fnPads[0]] = -fXyq[5][i] - 1;
 522       if (cath && div == 1) fDigitId[fnPads[0]] = fDigitId[i];
 523       // If low pad charge take the other equal to 0
 524       //if (div == 1 && fXyq[2][fnPads[0]] < fgkZeroSuppression + 1.5*3) div = 2;
 525       fXyq[2][fnPads[0]] /= div;
 526       fXyq[0][fnPads[0]] = fXyq[0][i];
 527       fXyq[1][fnPads[0]] = fXyq[1][i];
 528       fPadIJ[2][fnPads[0]] = fPadIJ[2][i];
 529       fPadIJ[3][fnPads[0]] = fPadIJ[3][i];
 530       fPadIJ[0][fnPads[0]++] = 0;
 531     }
 532   } // if (sameSize)
 533
 534   // Check if one-cathode precluster
 535   i1 = fnPads[0]!=0 ? 0 : 1;
 536   i2 = fnPads[1]!=0 ? 1 : 0;
 537
 538   if (i1 != i2) { // two-cathode
 539
 540     Int_t *flags = new Int_t[npad];
 541     for (Int_t i=0; i<npad; i++) { flags[i] = 0; }
 542
 543     // Check pad overlaps
 544     for (Int_t i=0; i<npad; i++) {
 545       if (fPadIJ[0][i] != i1) continue;
 546       xy1[0] = fXyq[0][i] - fXyq[3][i];
 547       xy1[1] = fXyq[0][i] + fXyq[3][i];
 548       xy1[2] = fXyq[1][i] - fXyq[4][i];
 549       xy1[3] = fXyq[1][i] + fXyq[4][i];
 550       for (Int_t j=0; j<npad; j++) {
 551         if (fPadIJ[0][j] != i2) continue;
 552         if (!Overlap(xy1, j, xy12, 0)) continue;
 553         flags[i] = flags[j] = 1; // mark overlapped pads
 554       } // for (Int_t j=0;
 555     } // for (Int_t i=0;
 556
 557     // Check if all pads overlap
 558     Int_t nFlags=0;
 559     for (Int_t i=0; i<npad; i++) {
 560       if (flags[i]) continue;
 561       nFlags ++;
 562       if (fDebug) cout << i << " " << fPadIJ[0][i] << " " << fXyq[0][i] << " " << fXyq[1][i] << endl;
 563     }
 564     if (fDebug && nFlags) cout << " nFlags = " << nFlags << endl;
 565     //if (nFlags > 2 || (Float_t)nFlags / npad > 0.2) { // why 2 ??? - empirical choice
 566     if (nFlags > 0) {
 567       for (Int_t i=0; i<npad; i++) {
 568         if (flags[i]) continue;
 569         cath = fPadIJ[0][i];
 570         // Check for edge effect (missing pads on the other cathode)
 571         Int_t cath1 = TMath::Even(cath), ix, iy;
 572         ix = iy = 0;
 573         AliMpPad pad = fSegmentation[cath1]->PadByPosition(TVector2(fXyq[0][i], fXyq[1][i]));
 574         if (!pad.IsValid()) continue;
 575         if (nFlags == 1 && fXyq[2][i] < fgkZeroSuppression * 3) continue;
 576         Used(i,kFALSE);
 577         fXyq[2][i] = -2;
 578         fnPads[cath]--;
 579       }
 580       //      if (fDraw) fDraw->UpdateCluster(npad);
 581     } // if (nFlags > 2)
 582
 583     // Check correlations of cathode charges
 584     if (fnPads[0] && fnPads[1]) { // two-cathode
 585       Double_t sum[2]={0};
 586       Int_t over[2] = {1, 1};
 587       for (Int_t i=0; i<npad; i++) {
 588         cath = fPadIJ[0][i];
 589         if (fXyq[2][i] > 0) sum[cath] += fXyq[2][i];
 590         if (fXyq[2][i] > fgkSaturation-1) over[cath] = 0;
 591       }
 592       if (fDebug) cout << " Total charge: " << sum[0] << " " << sum[1] << endl;
 593       if ((over[0] || over[1]) && TMath::Abs(sum[0]-sum[1])/(sum[0]+sum[1])*2 > 1) { // 3 times difference
 594         if (fDebug) cout << " Release " << endl;
 595         // Big difference
 596         cath = sum[0] > sum[1] ? 0 : 1;
 597         Int_t imax = 0, imin = 0;
 598         Double_t cmax = -1, cmin = 9999, dxMin = 0, dyMin = 0;
 599         Double_t *dist = new Double_t[npad];
 600         for (Int_t i = 0; i < npad; i++) {
 601           if (fPadIJ[0][i] != cath || fXyq[2][i] < 0) continue;
 602           if (fXyq[2][i] < cmin) {
 603             cmin = fXyq[2][i];
 604             imin = i;
 605           }
 606           if (fXyq[2][i] < cmax) continue;
 607           cmax = fXyq[2][i];
 608           imax = i;
 609         }
 610         // Arrange pads according to their distance to the max,
 611         // normalized to the pad size
 612         for (Int_t i = 0; i < npad; i++) {
 613           dist[i] = 0;
 614           if (fPadIJ[0][i] != cath || fXyq[2][i] < 0) continue;
 615           if (i == imax) continue;
 616           Double_t dx = (fXyq[0][i] - fXyq[0][imax]) / fXyq[3][imax] / 2;
 617           Double_t dy = (fXyq[1][i] - fXyq[1][imax]) / fXyq[4][imax] / 2;
 618           dist[i] = TMath::Sqrt (dx * dx + dy * dy);
 619           if (i == imin) {
 620             cmin = dist[i] + 0.001; // distance to the pad with minimum charge
 621             dxMin = dx;
 622             dyMin = dy;
 623           }
 624         }
 625         TMath::Sort(npad, dist, flags, kFALSE); // in increasing order
 626         Int_t indx;
 627         Double_t xmax = -1;
 628         for (Int_t i = 0; i < npad; i++) {
 629           indx = flags[i];
 630           if (fPadIJ[0][indx] != cath || fXyq[2][indx] < 0) continue;
 631           if (dist[indx] > cmin) {
 632             // Farther than the minimum pad
 633             Double_t dx = (fXyq[0][indx] - fXyq[0][imax]) / fXyq[3][imax] / 2;
 634             Double_t dy = (fXyq[1][indx] - fXyq[1][imax]) / fXyq[4][imax] / 2;
 635             dx *= dxMin;
 636             dy *= dyMin;
 637             if (dx >= 0 && dy >= 0) continue;
 638             if (TMath::Abs(dx) > TMath::Abs(dy) && dx >= 0) continue;
 639             if (TMath::Abs(dy) > TMath::Abs(dx) && dy >= 0) continue;
 640           }
 641           if (fXyq[2][indx] <= cmax || TMath::Abs(dist[indx]-xmax) < 1.e-3) {
 642             // Release pads
 643             if (TMath::Abs(dist[indx]-xmax) < 1.e-3)
 644                 cmax = TMath::Max((Double_t)(fXyq[2][indx]),cmax);
 645             else cmax = fXyq[2][indx];
 646             xmax = dist[indx];
 647       Used(indx,kFALSE);
 648             fXyq[2][indx] = -2;
 649             fnPads[cath]--;
 650           }
 651         } // for (Int_t i = 0; i < npad;
 652
 653         // Check pad overlaps once more
 654         for (Int_t j = 0; j < npad; j++) flags[j] = 0;
 655         for (Int_t k = 0; k < npad; k++) {
 656           if (fXyq[2][k] < 0 || fPadIJ[0][k] != i1) continue;
 657           xy1[0] = fXyq[0][k] - fXyq[3][k];
 658           xy1[1] = fXyq[0][k] + fXyq[3][k];
 659           xy1[2] = fXyq[1][k] - fXyq[4][k];
 660           xy1[3] = fXyq[1][k] + fXyq[4][k];
 661           for (Int_t j = 0; j < npad; j++) {
 662             if (fXyq[2][j] < 0) continue;
 663             if (fPadIJ[0][j] != i2) continue;
 664             if (!Overlap(xy1, j, xy12, 0)) continue;
 665             flags[k] = flags[j] = 1; // mark overlapped pads
 666           } // for (Int_t j = 0;
 667         } // for (Int_t k = 0;
 668         nFlags = 0;
 669         for (Int_t j = 0; j < npad; j++) {
 670           if (fXyq[2][j] < 0 || flags[j]) continue;
 671           nFlags ++;
 672         }
 673         if (nFlags == fnPads[0] + fnPads[1]) {
 674           // No overlap
 675           for (Int_t j = 0; j < npad; j++) {
 676             if (fXyq[2][j] < 0 || fPadIJ[0][j] != cath) continue;
 677             fXyq[2][j] = -2;
 678             fnPads[cath]--;
 679           }
 680         }
 681         delete [] dist; dist = 0;
 682 //      if (fDraw) fDraw->UpdateCluster(npad);
 683       } // TMath::Abs(sum[0]-sum[1])...
 684     } // if (fnPads[0] && fnPads[1])
 685     delete [] flags; flags = 0;
 686   } // if (i1 != i2)
 687
 688   if (!sameSize) { nShown[0] += fnPads[0]; nShown[1] += fnPads[1]; }
 689
 690   // Move released pads to the right
 691   Int_t beg = 0, end = npad-1, padij;
 692   Double_t xyq;
 693   while (beg < end) {
 694     if (fXyq[2][beg] > 0) { beg++; continue; }
 695     for (Int_t j=end; j>beg; j--) {
 696       if (fXyq[2][j] < 0) continue;
 697       end = j - 1;
 698       for (Int_t j1=0; j1<4; j1++) {
 699         padij = fPadIJ[j1][beg];
 700         fPadIJ[j1][beg] = fPadIJ[j1][j];
 701         fPadIJ[j1][j] = padij;
 702       }
 703       for (Int_t j1=0; j1<6; j1++) {
 704         xyq = fXyq[j1][beg];
 705         fXyq[j1][beg] = fXyq[j1][j];
 706         fXyq[j1][j] = xyq;
 707       }
 708       break;
 709     } // for (Int_t j=end;
 710     beg++;
 711   } // while
 712   npad = fnPads[0] + fnPads[1];
 713   if (npad > 500) {
 714     AliWarning(Form(" *** Too large cluster. Give up. %d ", npad));
 715     return kFALSE;
 716   }
 717   // Back up charge value
 718   for (Int_t j = 0; j < npad; j++) fXyq[6][j] = fXyq[2][j];
 719
 720   return kTRUE;
 721 }
 722
 723 //_____________________________________________________________________________
 724 void AliMUONClusterFinderAZ::BuildPixArray()
 725 {
 726 /// Build pixel array for MLEM method
 727
 728   Int_t nPix=0, i1, i2;
 729   Float_t xy1[4], xy12[4];
 730   AliMUONPixel *pixPtr=0;
 731
 732   Int_t npad = fnPads[0] + fnPads[1];
 733
 734   // One cathode is empty
 735   i1 = fnPads[0]!=0 ? 0 : 1;
 736   i2 = fnPads[1]!=0 ? 1 : 0;
 737
 738   // Build array of pixels on anode plane
 739   if (i1 == i2) { // one-cathode precluster
 740     for (Int_t j=0; j<npad; j++) {
 741       pixPtr = new AliMUONPixel();
 742       for (Int_t i=0; i<2; i++) {
 743         pixPtr->SetCoord(i, fXyq[i][j]); // pixel coordinates
 744         pixPtr->SetSize(i, fXyq[i+3][j]); // pixel size
 745       }
 746       pixPtr->SetCharge(fXyq[2][j]); // charge
 747       fPixArray->Add((TObject*)pixPtr);
 748       nPix++;
 749     }
 750   } else { // two-cathode precluster
 751     i1 = fPadIJ[0][0];
 752     i2 = TMath::Even (i1);
 753     for (Int_t i = 0; i < npad; i++) {
 754       if (fPadIJ[0][i] != i1) continue;
 755       xy1[0] = fXyq[0][i] - fXyq[3][i];
 756       xy1[1] = fXyq[0][i] + fXyq[3][i];
 757       xy1[2] = fXyq[1][i] - fXyq[4][i];
 758       xy1[3] = fXyq[1][i] + fXyq[4][i];
 759       for (Int_t j = 1; j < npad; j++) {
 760         if (fPadIJ[0][j] != i2) continue;
 761         if (!Overlap(xy1, j, xy12, 1)) continue;
 762         pixPtr = new AliMUONPixel();
 763         for (Int_t k=0; k<2; k++) {
 764           pixPtr->SetCoord(k, (xy12[2*k]+xy12[2*k+1])/2); // pixel coordinates
 765           pixPtr->SetSize(k, xy12[2*k+1]-pixPtr->Coord(k)); // size
 766         }
 767         pixPtr->SetCharge(TMath::Min (fXyq[2][i],fXyq[2][j])); //charge
 768         fPixArray->Add((TObject*)pixPtr);
 769         //cout << nPix << " " << pixPtr->Coord(0) << " " << pixPtr->Size(0) << " " << pixPtr->Coord(1) << " " << pixPtr->Size(1) << " " << pixPtr->Charge() << endl;
 770         nPix++;
 771       } // for (Int_t j=0;
 772     } // for (Int_t i=0;
 773   } // else
 774
 775   Float_t xPadMin = 999, yPadMin = 999;
 776   for (Int_t i = 0; i < npad; i++) {
 777     xPadMin = TMath::Min (xPadMin, fXyq[3][i]);
 778     yPadMin = TMath::Min (yPadMin, fXyq[4][i]);
 779   }
 780   if (fDebug) cout << xPadMin << " " << yPadMin << endl;
 781
 782   Float_t wxmin = 999, wymin = 999;
 783   for (Int_t i = 0; i < nPix; i++) {
 784     pixPtr = (AliMUONPixel*) fPixArray->UncheckedAt(i);
 785     wxmin = TMath::Min ((Double_t)wxmin, pixPtr->Size(0));
 786     wymin = TMath::Min ((Double_t)wymin, pixPtr->Size(1));
 787   }
 788   if (fDebug) cout << wxmin << " " << wymin << endl;
 789   wxmin = TMath::Abs (wxmin - xPadMin/2) > 0.001 ? xPadMin : xPadMin / 2;
 790   wymin = TMath::Abs (wymin - yPadMin/2) > 0.001 ? yPadMin : yPadMin / 2;
 791   //wxmin = xPadMin; wymin = yPadMin;
 792
 793   // Check if small pixel X-size
 794   AdjustPixel(wxmin, 0);
 795   // Check if small pixel Y-size
 796   AdjustPixel(wymin, 1);
 797   // Check if large pixel size
 798   AdjustPixel(wxmin, wymin);
 799
 800   // Remove discarded pixels
 801   for (Int_t i=0; i<nPix; i++) {
 802     pixPtr = (AliMUONPixel*) fPixArray->UncheckedAt(i);
 803     //pixPtr->Print();
 804     if (pixPtr->Charge() < 1) { fPixArray->RemoveAt(i); delete pixPtr; }// discarded pixel
 805   }
 806   fPixArray->Compress();
 807   nPix = fPixArray->GetEntriesFast();
 808
 809   if (nPix > npad) {
 810     if (fDebug) cout << nPix << endl;
 811     // Too many pixels - sort and remove pixels with the lowest signal
 812     fPixArray->Sort();
 813     for (Int_t i=npad; i<nPix; i++) {
 814       pixPtr = (AliMUONPixel*) fPixArray->UncheckedAt(i);
 815       //pixPtr->Print();
 816       fPixArray->RemoveAt(i);
 817       delete pixPtr;
 818     }
 819     nPix = npad;
 820   } // if (nPix > npad)
 821
 822   // Set pixel charges to the same value (for MLEM)
 823   for (Int_t i=0; i<nPix; i++) {
 824     pixPtr = (AliMUONPixel*) fPixArray->UncheckedAt(i);
 825     //pixPtr->SetCharge(10);
 826     if (fDebug) cout << i+1 << " " << pixPtr->Coord(0) << " " << pixPtr->Coord(1) << " " << pixPtr->Size(0) << " " << pixPtr->Size(1) << endl;
 827   }
 828 }
 829
 830 //_____________________________________________________________________________
 831 void AliMUONClusterFinderAZ::AdjustPixel(Float_t width, Int_t ixy)
 832 {
 833 /// Check if some pixels have small size (adjust if necessary)
 834
 835   AliMUONPixel *pixPtr, *pixPtr1 = 0;
 836   Int_t ixy1 = TMath::Even(ixy);
 837   Int_t nPix = fPixArray->GetEntriesFast();
 838
 839   for (Int_t i=0; i<nPix; i++) {
 840     pixPtr = (AliMUONPixel*) fPixArray->UncheckedAt(i);
 841     if (pixPtr->Charge() < 1) continue; // discarded pixel
 842     if (pixPtr->Size(ixy)-width < -1.e-4) {
 843       // try to merge
 844       if (fDebug) cout << i << " Small X or Y: " << ixy << " " << pixPtr->Size(ixy) << " " << width << " " << pixPtr->Coord(0) << " " << pixPtr->Coord(1) << endl;
 845       for (Int_t j=i+1; j<nPix; j++) {
 846         pixPtr1 = (AliMUONPixel*) fPixArray->UncheckedAt(j);
 847         if (pixPtr1->Charge() < 1) continue; // discarded pixel
 848         if (TMath::Abs(pixPtr1->Size(ixy)-width) < 1.e-4) continue; // right size
 849         if (TMath::Abs(pixPtr1->Coord(ixy1)-pixPtr->Coord(ixy1)) > 1.e-4) continue; // different rows/columns
 850         if (TMath::Abs(pixPtr1->Coord(ixy)-pixPtr->Coord(ixy)) < 2*width) {
 851           // merge
 852           Double_t tmp = pixPtr->Coord(ixy) + pixPtr1->Size(ixy) *
 853             TMath::Sign (1., pixPtr1->Coord(ixy) - pixPtr->Coord(ixy));
 854           pixPtr->SetCoord(ixy, tmp);
 855           pixPtr->SetSize(ixy, width);
 856           pixPtr->SetCharge(TMath::Min (pixPtr->Charge(),pixPtr1->Charge()));
 857           pixPtr1->SetCharge(0);
 858           pixPtr1 = 0;
 859           break;
 860         }
 861       } // for (Int_t j=i+1;
 862       //if (!pixPtr1) { cout << " I am here!" << endl; pixPtr->SetSize(ixy, width); } // ???
 863       //else if (pixPtr1->Charge() > 0.5 || i == nPix-1) {
 864       if (pixPtr1 || i == nPix-1) {
 865         // edge pixel - just increase its size
 866         if (fDebug) cout << " Edge ..." << endl;
 867         for (Int_t j=0; j<fnPads[0]+fnPads[1]; j++) {
 868           //if (fPadIJ[0][j] != ixy1) continue;
 869           //???-check if (TMath::Abs(pixPtr->Coord(ixy1)-fXyq[ixy1][j]) > 1.e-4) continue;
 870           if (pixPtr->Coord(ixy) < fXyq[ixy][j])
 871             //pixPtr->Shift(ixy, -pixPtr->Size(ixy));
 872             pixPtr->Shift(ixy, pixPtr->Size(ixy)-width);
 873           //else pixPtr->Shift(ixy, pixPtr->Size(ixy));
 874           else pixPtr->Shift(ixy, -pixPtr->Size(ixy)+width);
 875           pixPtr->SetSize(ixy, width);
 876           break;
 877         }
 878       }
 879     } // if (pixPtr->Size(ixy)-width < -1.e-4)
 880   } // for (Int_t i=0; i<nPix;
 881   return;
 882 }
 883
 884 //_____________________________________________________________________________
 885 void AliMUONClusterFinderAZ::AdjustPixel(Float_t wxmin, Float_t wymin)
 886 {
 887 /// Check if some pixels have large size (adjust if necessary)
 888
 889   Int_t n1[2], n2[2], iOK = 1, nPix = fPixArray->GetEntriesFast();
 890   AliMUONPixel *pixPtr, pix;
 891   Double_t xy0[2] = {9999, 9999}, wxy[2], dist[2] = {0};
 892
 893   // Check if large pixel size
 894   for (Int_t i = 0; i < nPix; i++) {
 895     pixPtr = (AliMUONPixel*) fPixArray->UncheckedAt(i);
 896     if (pixPtr->Charge() < 1) continue; // discarded pixel
 897     if (pixPtr->Size(0) - wxmin < 1.e-4) {
 898       if (xy0[0] > 9998) xy0[0] = pixPtr->Coord(0); // position of a "normal" pixel
 899       if (pixPtr->Size(1) - wymin < 1.e-4) {
 900         if (xy0[1] > 9998) xy0[1] = pixPtr->Coord(1); // position of a "normal" pixel
 901         continue;
 902       } else iOK = 0; // large pixel
 903     } else {
 904       iOK = 0; // large pixel
 905       if (xy0[1] > 9998 && pixPtr->Size(1) - wymin < 1.e-4) xy0[1] = pixPtr->Coord(1); // "normal" pixel
 906     }
 907     if (xy0[0] < 9998 && xy0[1] < 9998) break;
 908   }
 909   if (iOK) return;
 910
 911   wxy[0] = wxmin;
 912   wxy[1] = wymin;
 913   //cout << xy0[0] << " " << xy0[1] << endl;
 914   for (Int_t i = 0; i < nPix; i++) {
 915     pixPtr = (AliMUONPixel*) fPixArray->UncheckedAt(i);
 916     if (pixPtr->Charge() < 1) continue; // discarded pixel
 917     n1[0] = n1[1] = 999;
 918     n2[0] = n2[1] = 1;
 919     for (Int_t j = 0; j < 2; j++) {
 920       if (pixPtr->Size(j) - wxy[j] < 1.e-4) continue;
 921       dist[j] = (pixPtr->Coord(j) - xy0[j]) / wxy[j] / 2; // normalized distance to "normal" pixel
 922       n2[j] = TMath::Nint (pixPtr->Size(j) / wxy[j]);
 923       n1[j] = n2[j] == 1 ? TMath::Nint(dist[j]) : (Int_t)dist[j];
 924     }
 925     if (n1[0] > 998 && n1[1] > 998) continue;
 926     if (fDebug) cout << " Different " << pixPtr->Size(0) << " " << wxy[0] << " "
 927                      << pixPtr->Size(1) << " " << wxy[1] <<endl;
 928
 929     if (n2[0] > 2 || n2[1] > 2) {
 930       //cout << n2[0] << " " << n2[1] << endl;
 931       if (n2[0] > 2 && n1[0] < 999) n1[0]--;
 932       if (n2[1] > 2 && n1[1] < 999) n1[1]--;
 933     }
 934     //cout << n1[0] << " " << n2[0] << " " << n1[1] << " " << n2[1] << endl;
 935     pix = *pixPtr;
 936     pix.SetSize(0, wxy[0]); pix.SetSize(1, wxy[1]);
 937     //pixPtr->Print();
 938     for (Int_t ii = 0; ii < n2[0]; ii++) {
 939       if (n1[0] < 999) pix.SetCoord(0, xy0[0] + (n1[0] + TMath::Sign(1.,dist[0]) * ii) * 2 * wxy[0]);
 940       for (Int_t jj = 0; jj < n2[1]; jj++) {
 941         if (n1[1] < 999) pix.SetCoord(1, xy0[1] + (n1[1] + TMath::Sign(1.,dist[1]) * jj) * 2 * wxy[1]);
 942         fPixArray->Add(new AliMUONPixel(pix));
 943         //pix.Print();
 944       }
 945     }
 946     pixPtr->SetCharge(0);
 947   } // for (Int_t i = 0; i < nPix;
 948 }
 949
 950 //_____________________________________________________________________________
 951 Float_t
 952 AliMUONClusterFinderAZ::ChargeIntegration(Double_t x, Double_t y,
 953                                           Double_t padX, Double_t padY,
 954                                           Double_t padDX, Double_t padDY)
 955 {
 956   /// Compute the Mathieson integral on pad area, assuming the center
 957   /// of the Mathieson is at (x,y)
 958
 959   Double_t llx = x - padX - padDX;
 960   Double_t lly = y - padY - padDY;
 961   Double_t urx = llx + 2.0*padDX;
 962   Double_t ury = lly + 2.0*padDY;
 963
 964   return fMathieson->IntXY(llx,lly,urx,ury);
 965 }
 966
 967 //_____________________________________________________________________________
 968 Bool_t AliMUONClusterFinderAZ::MainLoop(Int_t iSimple)
 969 {
 970 /// Repeat MLEM algorithm until pixel size becomes sufficiently small
 971
 972   TH2D *mlem;
 973
 974   Int_t ix, iy;
 975   //Int_t nn, xList[10], yList[10];
 976   Int_t nPix = fPixArray->GetEntriesFast();
 977   AliMUONPixel *pixPtr = 0;
 978   Double_t *coef = 0, *probi = 0;
 979   AddVirtualPad(); // add virtual pads if necessary
 980   Int_t npadTot = fnPads[0] + fnPads[1], npadOK = 0;
 981   for (Int_t i = 0; i < npadTot; i++) if (fPadIJ[1][i] == 0) npadOK++;
 982 //  if (fDraw) fDraw->ResetMuon();
 983
 984   while (1) {
 985
 986     mlem = (TH2D*) gROOT->FindObject("mlem");
 987     if (mlem) mlem->Delete();
 988     // Calculate coefficients
 989     if (fDebug) cout << " nPix, npadTot, npadOK " << nPix << " " << npadTot << " " << npadOK << endl;
 990
 991     // Calculate coefficients and pixel visibilities
 992     coef = new Double_t [npadTot*nPix];
 993     probi = new Double_t [nPix];
 994     for (Int_t ipix=0; ipix<nPix; ipix++) probi[ipix] = 0;
 995     Int_t indx = 0, indx1 = 0;
 996
 997     for (Int_t j=0; j<npadTot; j++)
 998     {
 999       indx = j*nPix;
1000
1001 //      if (fPadIJ[1][j] == 0)
1002 //      {
1003 //        cath = fPadIJ[0][j];
1004 //        ix = fPadIJ[2][j];
1005 //        iy = fPadIJ[3][j];
1006 //        fSegmentation[cath]->SetPad(ix, iy);
1007 //      }
1008
1009       for (Int_t ipix=0; ipix<nPix; ipix++)
1010       {
1011         indx1 = indx + ipix;
1012         if (fPadIJ[1][j] < 0) { coef[indx1] = 0; continue; }
1013         pixPtr = (AliMUONPixel*) fPixArray->UncheckedAt(ipix);
1014
1015         Float_t q = ChargeIntegration(pixPtr->Coord(0),pixPtr->Coord(1),
1016                                       fXyq[0][j],fXyq[1][j],
1017                                       fXyq[3][j],fXyq[4][j]);
1018
1019 //        AliInfo(Form("pad %d pixel %d",j,ipix));
1020 //        PrintPad(j);
1021 //        PrintPixel(ipix);
1022
1023         coef[indx1] = q;
1024         probi[ipix] += coef[indx1];
1025
1026 //        AliInfo(Form("indx1=%d q=%e",indx1,q));
1027
1028       } // for (Int_t ipix=0;
1029     } // for (Int_t j=0;
1030
1031     for (Int_t ipix=0; ipix<nPix; ipix++) if (probi[ipix] < 0.01) pixPtr->SetCharge(0); // "invisible" pixel
1032
1033     // MLEM algorithm
1034     Mlem(coef, probi, 15);
1035
1036     Double_t xylim[4] = {999, 999, 999, 999};
1037     for (Int_t ipix=0; ipix<nPix; ipix++) {
1038       pixPtr = (AliMUONPixel*) fPixArray->UncheckedAt(ipix);
1039       //cout << ipix+1; pixPtr->Print();
1040       for (Int_t i=0; i<4; i++)
1041         xylim[i] = TMath::Min (xylim[i], (i%2 ? -1 : 1)*pixPtr->Coord(i/2));
1042     }
1043     for (Int_t i=0; i<4; i++) {
1044       xylim[i] -= pixPtr->Size(i/2); if (fDebug) cout << (i%2 ? -1 : 1)*xylim[i] << " "; }
1045     if (fDebug) cout << endl;
1046
1047     // Adjust histogram to approximately the same limits as for the pads
1048     // (for good presentation)
1049 //    if (fDraw) fDraw->AdjustHist(xylim, pixPtr);
1050
1051     Int_t nx = TMath::Nint ((-xylim[1]-xylim[0])/pixPtr->Size(0)/2);
1052     Int_t ny = TMath::Nint ((-xylim[3]-xylim[2])/pixPtr->Size(1)/2);
1053
1054     mlem = new TH2D("mlem","mlem",nx,xylim[0],-xylim[1],ny,xylim[2],-xylim[3]);
1055     for (Int_t ipix=0; ipix<nPix; ipix++) {
1056       pixPtr = (AliMUONPixel*) fPixArray->UncheckedAt(ipix);
1057       mlem->Fill(pixPtr->Coord(0),pixPtr->Coord(1),pixPtr->Charge());
1058     }
1059 //    if (fDraw) fDraw->DrawHist("c2", mlem);
1060
1061     // Check if the total charge of pixels is too low
1062     Double_t qTot = 0;
1063     for (Int_t i=0; i<nPix; i++) {
1064       pixPtr = (AliMUONPixel*) fPixArray->UncheckedAt(i);
1065       qTot += pixPtr->Charge();
1066     }
1067     if (qTot < 1.e-4 || npadOK < 3 && qTot < 7) {
1068       delete [] coef; delete [] probi; coef = 0; probi = 0;
1069       fPixArray->Delete();
1070       for (Int_t i=0; i<npadTot; i++) if (fPadIJ[1][i] == 0) fPadIJ[1][i] = -1;
1071       return kFALSE;
1072     }
1073
1074     // Plot data - expectation
1075     /*
1076     Double_t x, y, cont;
1077     for (Int_t j=0; j<npadTot; j++) {
1078       Double_t sum1 = 0;
1079       for (Int_t i=0; i<nPix; i++) {
1080         // Caculate expectation
1081         pixPtr = (AliMUONPixel*) fPixArray->UncheckedAt(i);
1082         sum1 += pixPtr->Charge()*coef[j*nPix+i];
1083       }
1084       sum1 = TMath::Min (sum1,fgkSaturation);
1085       x = fXyq[0][j];
1086       y = fXyq[1][j];
1087       cath = fPadIJ[0][j];
1088       Int_t ihist = cath*2;
1089       ix = fHist[ihist]->GetXaxis()->FindBin(x);
1090       iy = fHist[ihist]->GetYaxis()->FindBin(y);
1091       cont = fHist[ihist]->GetCellContent(ix,iy);
1092       if (cont == 0 && fHist[ihist+1]) {
1093         ihist += 1;
1094         ix = fHist[ihist]->GetXaxis()->FindBin(x);
1095         iy = fHist[ihist]->GetYaxis()->FindBin(y);
1096       }
1097       fHist[ihist]->SetBinContent(ix,iy,fXyq[2][j]-sum1);
1098     }
1099     ((TCanvas*)gROOT->FindObject("c1"))->cd(1);
1100     //gPad->SetTheta(55);
1101     //gPad->SetPhi(30);
1102     //mlem->Draw("lego1");
1103     gPad->Modified();
1104     ((TCanvas*)gROOT->FindObject("c1"))->cd(2);
1105     gPad->Modified();
1106     */
1107
1108     if (iSimple) {
1109       // Simple cluster - skip further passes thru EM-procedure
1110       Simple();
1111       delete [] coef; delete [] probi; coef = 0; probi = 0;
1112       fPixArray->Delete();
1113       return kTRUE;
1114     }
1115
1116     // Calculate position of the center-of-gravity around the maximum pixel
1117     Double_t xyCOG[2];
1118     FindCOG(mlem, xyCOG);
1119
1120     if (TMath::Min(pixPtr->Size(0),pixPtr->Size(1)) < 0.07 && pixPtr->Size(0) > pixPtr->Size(1)) break;
1121     //if (TMath::Min(pixPtr->Size(0),pixPtr->Size(1)) < 0.007 && pixPtr->Size(0) > pixPtr->Size(1)) break;
1122     //if (TMath::Min(pixPtr->Size(0),pixPtr->Size(1)) >= 0.07 || pixPtr->Size(0) < pixPtr->Size(1)) {
1123     // Sort pixels according to the charge
1124     fPixArray->Sort();
1125     /*
1126     for (Int_t i=0; i<nPix; i++) {
1127       pixPtr = (AliMUONPixel*) fPixArray->UncheckedAt(i);
1128       cout << i+1; pixPtr->Print();
1129     }
1130     */
1131     Double_t pixMin = 0.01*((AliMUONPixel*)fPixArray->UncheckedAt(0))->Charge();
1132     pixMin = TMath::Min (pixMin,50.);
1133
1134     // Decrease pixel size and shift pixels to make them centered at
1135     // the maximum one
1136     indx = (pixPtr->Size(0)>pixPtr->Size(1)) ? 0 : 1;
1137     Double_t width = 0, shift[2]={0};
1138     ix = 1;
1139     for (Int_t i=0; i<4; i++) xylim[i] = 999;
1140     Int_t nPix1 = nPix; nPix = 0;
1141     for (Int_t ipix=0; ipix<nPix1; ipix++) {
1142       pixPtr = (AliMUONPixel*) fPixArray->UncheckedAt(ipix);
1143       if (nPix >= npadOK) { // too many pixels already
1144         fPixArray->RemoveAt(ipix);
1145         delete pixPtr;
1146         continue;
1147       }
1148       if (pixPtr->Charge() < pixMin) { // low charge
1149         fPixArray->RemoveAt(ipix);
1150         delete pixPtr;
1151         continue;
1152       }
1153       for (Int_t i=0; i<2; i++) {
1154         if (!i) {
1155           pixPtr->SetCharge(10);
1156           pixPtr->SetSize(indx, pixPtr->Size(indx)/2);
1157           width = -pixPtr->Size(indx);
1158           pixPtr->Shift(indx, width);
1159           // Shift pixel position
1160           if (ix) {
1161             ix = 0;
1162             for (Int_t j=0; j<2; j++) {
1163               shift[j] = pixPtr->Coord(j) - xyCOG[j];
1164               shift[j] -= ((Int_t)(shift[j]/pixPtr->Size(j)/2))*pixPtr->Size(j)*2;
1165             }
1166             //cout << ipix << " " << i << " " << shift[0] << " " << shift[1] << endl;
1167           } // if (ix)
1168           pixPtr->Shift(0, -shift[0]);
1169           pixPtr->Shift(1, -shift[1]);
1170         } else {
1171           pixPtr = new AliMUONPixel(*pixPtr);
1172           pixPtr->Shift(indx, -2*width);
1173           fPixArray->Add((TObject*)pixPtr);
1174         } // else
1175         //pixPtr->Print();
1176         for (Int_t i=0; i<4; i++)
1177           xylim[i] = TMath::Min (xylim[i], (i%2 ? -1 : 1)*pixPtr->Coord(i/2));
1178       } // for (Int_t i=0; i<2;
1179       nPix += 2;
1180     } // for (Int_t ipix=0;
1181
1182     fPixArray->Compress();
1183     nPix = fPixArray->GetEntriesFast();
1184
1185     // Remove excessive pixels
1186     if (nPix > npadOK) {
1187       for (Int_t ipix=npadOK; ipix<nPix; ipix++) {
1188         pixPtr = (AliMUONPixel*) fPixArray->UncheckedAt(ipix);
1189         fPixArray->RemoveAt(ipix);
1190         delete pixPtr;
1191       }
1192     } else {
1193       pixPtr = (AliMUONPixel*) fPixArray->UncheckedAt(0);
1194       // add pixels if the maximum is at the limit of pixel area
1195       // start from Y-direction
1196       Int_t j = 0;
1197       for (Int_t i=3; i>-1; i--) {
1198         if (nPix < npadOK &&
1199             TMath::Abs((i%2 ? -1 : 1)*xylim[i]-xyCOG[i/2]) < pixPtr->Size(i/2)) {
1200           pixPtr = new AliMUONPixel(*pixPtr);
1201           pixPtr->SetCoord(i/2, xyCOG[i/2]+(i%2 ? 2:-2)*pixPtr->Size(i/2));
1202           j = TMath::Even (i/2);
1203           pixPtr->SetCoord(j, xyCOG[j]);
1204           fPixArray->Add((TObject*)pixPtr);
1205           nPix++;
1206         }
1207       }
1208     } // else
1209
1210     fPixArray->Compress();
1211     nPix = fPixArray->GetEntriesFast();
1212     delete [] coef; delete [] probi; coef = 0; probi = 0;
1213   } // while (1)
1214
1215   // remove pixels with low signal or low visibility
1216   // Cuts are empirical !!!
1217   Double_t thresh = TMath::Max (mlem->GetMaximum()/100.,1.);
1218   thresh = TMath::Min (thresh,50.);
1219   Double_t cmax = -1, charge = 0;
1220   for (Int_t i=0; i<nPix; i++) cmax = TMath::Max (cmax,probi[i]);
1221   //cout << thresh << " " << cmax << " " << cmax*0.9 << endl;
1222   // Mark pixels which should be removed
1223   for (Int_t i=0; i<nPix; i++) {
1224     pixPtr = (AliMUONPixel*) fPixArray->UncheckedAt(i);
1225     charge = pixPtr->Charge();
1226     if (charge < thresh) pixPtr->SetCharge(-charge);
1227     //else if (cmax > 1.91) {
1228     //  if (probi[i] < 1.9) pixPtr->SetCharge(-charge);
1229     //}
1230     //AZ else if (probi[i] < cmax*0.9) pixPtr->SetCharge(-charge);
1231     //18-01-06 else if (probi[i] < cmax*0.8) pixPtr->SetCharge(-charge);
1232     //cout << i << " " << pixPtr->Coord(0) << " " << pixPtr->Coord(1) << " " << charge << " " << probi[i] << endl;
1233   }
1234   // Move charge of removed pixels to their nearest neighbour (to keep total charge the same)
1235   Int_t near = 0;
1236   for (Int_t i=0; i<nPix; i++) {
1237     pixPtr = (AliMUONPixel*) fPixArray->UncheckedAt(i);
1238     charge = pixPtr->Charge();
1239     if (charge > 0) continue;
1240     near = FindNearest(pixPtr);
1241     pixPtr->SetCharge(0);
1242     probi[i] = 0; // make it "invisible"
1243     pixPtr = (AliMUONPixel*) fPixArray->UncheckedAt(near);
1244     pixPtr->SetCharge(pixPtr->Charge() + (-charge));
1245   }
1246   Mlem(coef,probi,2);
1247   // Update histogram
1248   for (Int_t i=0; i<nPix; i++) {
1249     pixPtr = (AliMUONPixel*) fPixArray->UncheckedAt(i);
1250     ix = mlem->GetXaxis()->FindBin(pixPtr->Coord(0));
1251     iy = mlem->GetYaxis()->FindBin(pixPtr->Coord(1));
1252     mlem->SetBinContent(ix, iy, pixPtr->Charge());
1253   }
1254 //  if (fDraw) fDraw->DrawHist("c2", mlem);
1255
1256   // Try to split into clusters
1257   Bool_t ok = kTRUE;
1258   if (mlem->GetSum() < 1) ok = kFALSE;
1259   else Split(mlem, coef);
1260   delete [] coef; delete [] probi; coef = 0; probi = 0;
1261   fPixArray->Delete();
1262   return ok;
1263 }
1264
1265 //_____________________________________________________________________________
1266 void AliMUONClusterFinderAZ::Mlem(Double_t *coef, Double_t *probi, Int_t nIter)
1267 {
1268 /// Use MLEM to find pixel charges
1269
1270   Int_t nPix = fPixArray->GetEntriesFast();
1271   Int_t npad = fnPads[0] + fnPads[1];
1272   Double_t *probi1 = new Double_t [nPix];
1273   Double_t probMax = 0;
1274   Int_t indx, indx1;
1275   AliMUONPixel *pixPtr;
1276
1277   for (Int_t ipix=0; ipix<nPix; ipix++) if (probi[ipix] > probMax) probMax = probi[ipix];
1278   for (Int_t iter=0; iter<nIter; iter++) {
1279     // Do iterations
1280     for (Int_t ipix=0; ipix<nPix; ipix++) {
1281       // Correct each pixel
1282       if (probi[ipix] < 0.01) continue; // skip "invisible" pixel
1283       Double_t sum = 0;
1284       //probi1[ipix] = probi[ipix];
1285       probi1[ipix] = probMax;
1286       for (Int_t j=0; j<npad; j++) {
1287         if (fPadIJ[1][j] < 0) continue;
1288         Double_t sum1 = 0;
1289         indx1 = j*nPix;
1290         indx = indx1 + ipix;
1291         for (Int_t i=0; i<nPix; i++) {
1292           // Caculate expectation
1293           pixPtr = (AliMUONPixel*) fPixArray->UncheckedAt(i);
1294           sum1 += pixPtr->Charge()*coef[indx1+i];
1295         } // for (Int_t i=0;
1296         if (fXyq[2][j] > fgkSaturation-1 && sum1 > fXyq[2][j]) { probi1[ipix] -= coef[indx]; continue; } // correct for pad charge overflows
1297         //cout << sum1 << " " << fXyq[2][j] << " " << coef[j*nPix+ipix] << endl;
1298         if (coef[indx] > 1.e-6) sum += fXyq[2][j]*coef[indx]/sum1;
1299       } // for (Int_t j=0;
1300       pixPtr = (AliMUONPixel*) fPixArray->UncheckedAt(ipix);
1301       if (probi1[ipix] > 1.e-6) pixPtr->SetCharge(pixPtr->Charge()*sum/probi1[ipix]);
1302     } // for (Int_t ipix=0;
1303   } // for (Int_t iter=0;
1304   delete [] probi1;
1305   return;
1306 }
1307
1308 //_____________________________________________________________________________
1309 void AliMUONClusterFinderAZ::FindCOG(TH2D *mlem, Double_t *xyc)
1310 {
1311 /// Calculate position of the center-of-gravity around the maximum pixel
1312
1313   Int_t ixmax, iymax, ix, nsumx=0, nsumy=0, nsum=0;
1314   Int_t i1 = -9, j1 = -9;
1315   mlem->GetMaximumBin(ixmax,iymax,ix);
1316   Int_t nx = mlem->GetNbinsX();
1317   Int_t ny = mlem->GetNbinsY();
1318   Double_t thresh = mlem->GetMaximum()/10;
1319   Double_t x, y, cont, xq=0, yq=0, qq=0;
1320
1321   for (Int_t i=TMath::Max(1,iymax-1); i<=TMath::Min(ny,iymax+1); i++) {
1322     y = mlem->GetYaxis()->GetBinCenter(i);
1323     for (Int_t j=TMath::Max(1,ixmax-1); j<=TMath::Min(nx,ixmax+1); j++) {
1324       cont = mlem->GetCellContent(j,i);
1325       if (cont < thresh) continue;
1326       if (i != i1) {i1 = i; nsumy++;}
1327       if (j != j1) {j1 = j; nsumx++;}
1328       x = mlem->GetXaxis()->GetBinCenter(j);
1329       xq += x*cont;
1330       yq += y*cont;
1331       qq += cont;
1332       nsum++;
1333     }
1334   }
1335
1336   Double_t cmax = 0;
1337   Int_t i2 = 0, j2 = 0;
1338   x = y = 0;
1339   if (nsumy == 1) {
1340     // one bin in Y - add one more (with the largest signal)
1341     for (Int_t i=TMath::Max(1,iymax-1); i<=TMath::Min(ny,iymax+1); i++) {
1342       if (i == iymax) continue;
1343       for (Int_t j=TMath::Max(1,ixmax-1); j<=TMath::Min(nx,ixmax+1); j++) {
1344         cont = mlem->GetCellContent(j,i);
1345         if (cont > cmax) {
1346           cmax = cont;
1347           x = mlem->GetXaxis()->GetBinCenter(j);
1348           y = mlem->GetYaxis()->GetBinCenter(i);
1349           i2 = i;
1350           j2 = j;
1351         }
1352       }
1353     }
1354     xq += x*cmax;
1355     yq += y*cmax;
1356     qq += cmax;
1357     if (i2 != i1) nsumy++;
1358     if (j2 != j1) nsumx++;
1359     nsum++;
1360   } // if (nsumy == 1)
1361
1362   if (nsumx == 1) {
1363     // one bin in X - add one more (with the largest signal)
1364     cmax = x = y = 0;
1365     for (Int_t j=TMath::Max(1,ixmax-1); j<=TMath::Min(nx,ixmax+1); j++) {
1366       if (j == ixmax) continue;
1367       for (Int_t i=TMath::Max(1,iymax-1); i<=TMath::Min(ny,iymax+1); i++) {
1368         cont = mlem->GetCellContent(j,i);
1369         if (cont > cmax) {
1370           cmax = cont;
1371           x = mlem->GetXaxis()->GetBinCenter(j);
1372           y = mlem->GetYaxis()->GetBinCenter(i);
1373           i2 = i;
1374           j2 = j;
1375         }
1376       }
1377     }
1378     xq += x*cmax;
1379     yq += y*cmax;
1380     qq += cmax;
1381     if (i2 != i1) nsumy++;
1382     if (j2 != j1) nsumx++;
1383     nsum++;
1384   } // if (nsumx == 1)
1385
1386   xyc[0] = xq/qq; xyc[1] = yq/qq;
1387   if (fDebug) cout << xyc[0] << " " << xyc[1] << " " << qq << " " << nsum << " " << nsumx << " " << nsumy << endl;
1388   return;
1389 }
1390
1391 //_____________________________________________________________________________
1392 Int_t AliMUONClusterFinderAZ::FindNearest(AliMUONPixel *pixPtr0)
1393 {
1394 /// Find the pixel nearest to the given one
1395 /// (algorithm may be not very efficient)
1396
1397   Int_t nPix = fPixArray->GetEntriesFast(), imin = 0;
1398   Double_t rmin = 99999, dx = 0, dy = 0, r = 0;
1399   Double_t xc = pixPtr0->Coord(0), yc = pixPtr0->Coord(1);
1400   AliMUONPixel *pixPtr;
1401
1402   for (Int_t i=0; i<nPix; i++) {
1403     pixPtr = (AliMUONPixel*) fPixArray->UncheckedAt(i);
1404     if (pixPtr->Charge() < 0.5) continue;
1405     dx = (xc - pixPtr->Coord(0)) / pixPtr->Size(0);
1406     dy = (yc - pixPtr->Coord(1)) / pixPtr->Size(1);
1407     r = dx *dx + dy * dy;
1408     if (r < rmin) { rmin = r; imin = i; }
1409   }
1410   return imin;
1411 }
1412
1413 //_____________________________________________________________________________
1414 void AliMUONClusterFinderAZ::Split(TH2D *mlem, Double_t *coef)
1415 {
1416 /// The main steering function to work with clusters of pixels in anode
1417 /// plane (find clusters, decouple them from each other, merge them (if
1418 /// necessary), pick up coupled pads, call the fitting function)
1419
1420   Int_t nx = mlem->GetNbinsX();
1421   Int_t ny = mlem->GetNbinsY();
1422   Int_t nPix = fPixArray->GetEntriesFast();
1423
1424   Bool_t *used = new Bool_t[ny*nx];
1425   Double_t cont;
1426   Int_t nclust = 0, indx, indx1;
1427
1428   for (Int_t i=0; i<ny*nx; i++) used[i] = kFALSE;
1429
1430   TObjArray *clusters[200]={0};
1431   TObjArray *pix;
1432
1433   // Find clusters of histogram bins (easier to work in 2-D space)
1434   for (Int_t i=1; i<=ny; i++) {
1435     for (Int_t j=1; j<=nx; j++) {
1436       indx = (i-1)*nx + j - 1;
1437       if (used[indx]) continue;
1438       cont = mlem->GetCellContent(j,i);
1439       if (cont < 0.5) continue;
1440       pix = new TObjArray(20);
1441       used[indx] = 1;
1442       pix->Add(BinToPix(mlem,j,i));
1443       AddBin(mlem, i, j, 0, used, pix); // recursive call
1444       if (nclust >= 200) AliFatal(" Too many clusters !!!");
1445       clusters[nclust++] = pix;
1446     } // for (Int_t j=1; j<=nx; j++) {
1447   } // for (Int_t i=1; i<=ny;
1448   if (fDebug) cout << nclust << endl;
1449   delete [] used; used = 0;
1450
1451   // Compute couplings between clusters and clusters to pads
1452   Int_t npad = fnPads[0] + fnPads[1];
1453
1454   // Write out some information for algorithm development
1455   Int_t cath=0, npadx[2]={0}, npady[2]={0};
1456   Double_t xlow[2]={9999,9999}, xhig[2]={-9999,-9999};
1457   Double_t ylow[2]={9999,9999}, yhig[2]={-9999,-9999};
1458   for (Int_t j=0; j<npad; j++) {
1459     if (fXyq[3][j] < 0) continue; // exclude virtual pads
1460     cath = fPadIJ[0][j];
1461     if (fXyq[0][j] < xlow[cath]-0.001) {
1462       if (fXyq[0][j]+fXyq[3][j] <= xlow[cath] && npadx[cath]) npadx[cath]++;
1463       xlow[cath] = fXyq[0][j];
1464     }
1465     if (fXyq[0][j] > xhig[cath]+0.001) {
1466       if (fXyq[0][j]-fXyq[3][j] >= xhig[cath]) npadx[cath]++;
1467       xhig[cath] = fXyq[0][j];
1468     }
1469     if (fXyq[1][j] < ylow[cath]-0.001) {
1470       if (fXyq[1][j]+fXyq[4][j] <= ylow[cath] && npady[cath]) npady[cath]++;
1471       ylow[cath] = fXyq[1][j];
1472     }
1473     if (fXyq[1][j] > yhig[cath]+0.001) {
1474       if (fXyq[1][j]-fXyq[4][j] >= yhig[cath]) npady[cath]++;
1475       yhig[cath] = fXyq[1][j];
1476     }
1477   }
1478
1479   // Exclude pads with overflows
1480   for (Int_t j=0; j<npad; j++) {
1481     if (fXyq[2][j] > fgkSaturation-1) fPadIJ[1][j] = -5;
1482     else fPadIJ[1][j] = 0;
1483   }
1484
1485   // Compute couplings of clusters to pads
1486   TMatrixD *aijclupad = new TMatrixD(nclust,npad);
1487   *aijclupad = 0;
1488   Int_t npxclu;
1489   for (Int_t iclust=0; iclust<nclust; iclust++) {
1490     pix = clusters[iclust];
1491     npxclu = pix->GetEntriesFast();
1492     for (Int_t i=0; i<npxclu; i++) {
1493       indx = fPixArray->IndexOf(pix->UncheckedAt(i));
1494       for (Int_t j=0; j<npad; j++) {
1495         if (fPadIJ[1][j] < 0 && fPadIJ[1][j] != -5) continue;
1496         if (coef[j*nPix+indx] < fgkCouplMin) continue;
1497         (*aijclupad)(iclust,j) += coef[j*nPix+indx];
1498       }
1499     }
1500   }
1501   // Compute couplings between clusters
1502   TMatrixD *aijcluclu = new TMatrixD(nclust,nclust);
1503   *aijcluclu = 0;
1504   for (Int_t iclust=0; iclust<nclust; iclust++) {
1505     for (Int_t j=0; j<npad; j++) {
1506       // Exclude overflows
1507       if (fPadIJ[1][j] < 0) continue;
1508       if ((*aijclupad)(iclust,j) < fgkCouplMin) continue;
1509       for (Int_t iclust1=iclust+1; iclust1<nclust; iclust1++) {
1510         if ((*aijclupad)(iclust1,j) < fgkCouplMin) continue;
1511         (*aijcluclu)(iclust,iclust1) +=
1512           TMath::Sqrt ((*aijclupad)(iclust,j)*(*aijclupad)(iclust1,j));
1513       }
1514     }
1515   }
1516   for (Int_t iclust=0; iclust<nclust; iclust++) {
1517     for (Int_t iclust1=iclust+1; iclust1<nclust; iclust1++) {
1518       (*aijcluclu)(iclust1,iclust) = (*aijcluclu)(iclust,iclust1);
1519     }
1520   }
1521
1522   if (fDebug && nclust > 1) aijcluclu->Print();
1523
1524   // Find groups of coupled clusters
1525   used = new Bool_t[nclust];
1526   for (Int_t i=0; i<nclust; i++) used[i] = kFALSE;
1527   Int_t *clustNumb = new Int_t[nclust];
1528   Int_t nCoupled, nForFit, minGroup[3], clustFit[3], nfit = 0;
1529   Double_t parOk[8];
1530
1531   for (Int_t igroup=0; igroup<nclust; igroup++) {
1532     if (used[igroup]) continue;
1533     used[igroup] = kTRUE;
1534     clustNumb[0] = igroup;
1535     nCoupled = 1;
1536     // Find group of coupled clusters
1537     AddCluster(igroup, nclust, aijcluclu, used, clustNumb, nCoupled); // recursive
1538     if (fDebug) {
1539       cout << " nCoupled: " << nCoupled << endl;
1540       for (Int_t i=0; i<nCoupled; i++) cout << clustNumb[i] << " "; cout << endl;
1541     }
1542     fnCoupled = nCoupled;
1543
1544     while (nCoupled > 0) {
1545
1546       if (nCoupled < 4) {
1547         nForFit = nCoupled;
1548         for (Int_t i=0; i<nCoupled; i++) clustFit[i] = clustNumb[i];
1549       } else {
1550         // Too many coupled clusters to fit - try to decouple them
1551         // Find the lowest coupling of 1, 2, min(3,nLinks/2) pixels with
1552         // all the others in the group
1553         for (Int_t j=0; j<3; j++) minGroup[j] = -1;
1554         Double_t coupl = MinGroupCoupl(nCoupled, clustNumb, aijcluclu, minGroup);
1555
1556         // Flag clusters for fit
1557         nForFit = 0;
1558         while (minGroup[nForFit] >= 0 && nForFit < 3) {
1559           if (fDebug) cout << clustNumb[minGroup[nForFit]] << " ";
1560           clustFit[nForFit] = clustNumb[minGroup[nForFit]];
1561           clustNumb[minGroup[nForFit]] -= 999;
1562           nForFit++;
1563         }
1564         if (fDebug) cout << nForFit << " " << coupl << endl;
1565       } // else
1566
1567       // Select pads for fit.
1568       if (SelectPad(nCoupled, nForFit, clustNumb, clustFit, aijclupad) < 3 && nCoupled > 1) {
1569         // Deselect pads
1570         for (Int_t j=0; j<npad; j++) {
1571           if (TMath::Abs(fPadIJ[1][j]) == 1) fPadIJ[1][j] = 0;
1572           if (TMath::Abs(fPadIJ[1][j]) == -9) fPadIJ[1][j] = -5;
1573         }
1574         // Merge the failed cluster candidates (with too few pads to fit) with
1575         // the one with the strongest coupling
1576         Merge(nForFit, nCoupled, clustNumb, clustFit, clusters, aijcluclu, aijclupad);
1577       } else {
1578         // Do the fit
1579         nfit = Fit(0, nForFit, clustFit, clusters, parOk);
1580       }
1581
1582       // Subtract the fitted charges from pads with strong coupling and/or
1583       // return pads for further use
1584       UpdatePads(nfit, parOk);
1585
1586       // Mark used pads
1587       for (Int_t j=0; j<npad; j++) {
1588         if (fPadIJ[1][j] == 1) fPadIJ[1][j] = -1;
1589         if (fPadIJ[1][j] == -9) fPadIJ[1][j] = -5;
1590       }
1591
1592       // Sort the clusters (move to the right the used ones)
1593       Int_t beg = 0, end = nCoupled - 1;
1594       while (beg < end) {
1595         if (clustNumb[beg] >= 0) { beg++; continue; }
1596         for (Int_t j=end; j>beg; j--) {
1597           if (clustNumb[j] < 0) continue;
1598           end = j - 1;
1599           indx = clustNumb[beg];
1600           clustNumb[beg] = clustNumb[j];
1601           clustNumb[j] = indx;
1602           break;
1603         }
1604         beg++;
1605       }
1606
1607       nCoupled -= nForFit;
1608       if (nCoupled > 3) {
1609         // Remove couplings of used clusters
1610         for (Int_t iclust=nCoupled; iclust<nCoupled+nForFit; iclust++) {
1611           indx = clustNumb[iclust] + 999;
1612           for (Int_t iclust1=0; iclust1<nCoupled; iclust1++) {
1613             indx1 = clustNumb[iclust1];
1614             (*aijcluclu)(indx,indx1) = (*aijcluclu)(indx1,indx) = 0;
1615           }
1616         }
1617
1618         // Update the remaining clusters couplings (exclude couplings from
1619         // the used pads)
1620         for (Int_t j=0; j<npad; j++) {
1621           if (fPadIJ[1][j] != -1) continue;
1622           for (Int_t iclust=0; iclust<nCoupled; iclust++) {
1623             indx = clustNumb[iclust];
1624             if ((*aijclupad)(indx,j) < fgkCouplMin) continue;
1625             for (Int_t iclust1=iclust+1; iclust1<nCoupled; iclust1++) {
1626               indx1 = clustNumb[iclust1];
1627               if ((*aijclupad)(indx1,j) < fgkCouplMin) continue;
1628               // Check this
1629               (*aijcluclu)(indx,indx1) -=
1630                 TMath::Sqrt ((*aijclupad)(indx,j)*(*aijclupad)(indx1,j));
1631               (*aijcluclu)(indx1,indx) = (*aijcluclu)(indx,indx1);
1632             }
1633           }
1634           fPadIJ[1][j] = -8;
1635         } // for (Int_t j=0; j<npad;
1636       } // if (nCoupled > 3)
1637     } // while (nCoupled > 0)
1638   } // for (Int_t igroup=0; igroup<nclust;
1639
1640   aijcluclu->Delete(); aijclupad->Delete();
1641   for (Int_t iclust=0; iclust<nclust; iclust++) {
1642     pix = clusters[iclust];
1643     pix->Clear();
1644     delete pix; pix = 0;
1645   }
1646   delete [] clustNumb; clustNumb = 0; delete [] used; used = 0;
1647 }
1648
1649 //_____________________________________________________________________________
1650 void AliMUONClusterFinderAZ::AddBin(TH2D *mlem, Int_t ic, Int_t jc, Int_t mode, Bool_t *used, TObjArray *pix)
1651 {
1652 /// Add a bin to the cluster
1653
1654   Int_t nx = mlem->GetNbinsX();
1655   Int_t ny = mlem->GetNbinsY();
1656   Double_t cont1, cont = mlem->GetCellContent(jc,ic);
1657   AliMUONPixel *pixPtr = 0;
1658
1659   for (Int_t i=TMath::Max(ic-1,1); i<=TMath::Min(ic+1,ny); i++) {
1660     for (Int_t j=TMath::Max(jc-1,1); j<=TMath::Min(jc+1,nx); j++) {
1661       if (i != ic && j != jc) continue;
1662       if (used[(i-1)*nx+j-1]) continue;
1663       cont1 = mlem->GetCellContent(j,i);
1664       if (mode && cont1 > cont) continue;
1665       used[(i-1)*nx+j-1] = kTRUE;
1666       if (cont1 < 0.5) continue;
1667       if (pix) pix->Add(BinToPix(mlem,j,i));
1668       else {
1669         pixPtr = new AliMUONPixel (mlem->GetXaxis()->GetBinCenter(j),
1670                                    mlem->GetYaxis()->GetBinCenter(i), 0, 0, cont1);
1671         fPixArray->Add((TObject*)pixPtr);
1672       }
1673       AddBin(mlem, i, j, mode, used, pix); // recursive call
1674     }
1675   }
1676 }
1677
1678 //_____________________________________________________________________________
1679 TObject* AliMUONClusterFinderAZ::BinToPix(TH2D *mlem, Int_t jc, Int_t ic)
1680 {
1681 /// Translate histogram bin to pixel
1682
1683   Double_t yc = mlem->GetYaxis()->GetBinCenter(ic);
1684   Double_t xc = mlem->GetXaxis()->GetBinCenter(jc);
1685
1686   Int_t nPix = fPixArray->GetEntriesFast();
1687   AliMUONPixel *pixPtr = NULL;
1688
1689   // Compare pixel and bin positions
1690   for (Int_t i=0; i<nPix; i++) {
1691     pixPtr = (AliMUONPixel*) fPixArray->UncheckedAt(i);
1692     if (pixPtr->Charge() < 0.5) continue;
1693     if (TMath::Abs(pixPtr->Coord(0)-xc)<1.e-4 && TMath::Abs(pixPtr->Coord(1)-yc)<1.e-4) return (TObject*) pixPtr;
1694   }
1695   AliError(Form(" Something wrong ??? %f %f ", xc, yc));
1696   return NULL;
1697 }
1698
1699 //_____________________________________________________________________________
1700 void AliMUONClusterFinderAZ::AddCluster(Int_t ic, Int_t nclust, TMatrixD *aijcluclu, Bool_t *used, Int_t *clustNumb, Int_t &nCoupled)
1701 {
1702 /// Add a cluster to the group of coupled clusters
1703
1704   for (Int_t i=0; i<nclust; i++) {
1705     if (used[i]) continue;
1706     if ((*aijcluclu)(i,ic) < fgkCouplMin) continue;
1707     used[i] = kTRUE;
1708     clustNumb[nCoupled++] = i;
1709     AddCluster(i, nclust, aijcluclu, used, clustNumb, nCoupled);
1710   }
1711 }
1712
1713 //_____________________________________________________________________________
1714 Double_t AliMUONClusterFinderAZ::MinGroupCoupl(Int_t nCoupled, Int_t *clustNumb, TMatrixD *aijcluclu, Int_t *minGroup)
1715 {
1716 /// Find group of clusters with minimum coupling to all the others
1717
1718   Int_t i123max = TMath::Min(3,nCoupled/2);
1719   Int_t indx, indx1, indx2, indx3, nTot = 0;
1720   Double_t *coupl1 = 0, *coupl2 = 0, *coupl3 = 0;
1721
1722   for (Int_t i123=1; i123<=i123max; i123++) {
1723
1724     if (i123 == 1) {
1725       coupl1 = new Double_t [nCoupled];
1726       for (Int_t i=0; i<nCoupled; i++) coupl1[i] = 0;
1727     }
1728     else if (i123 == 2) {
1729       nTot = nCoupled*nCoupled;
1730       coupl2 = new Double_t [nTot];
1731       for (Int_t i=0; i<nTot; i++) coupl2[i] = 9999;
1732     } else {
1733       nTot = nTot*nCoupled;
1734       coupl3 = new Double_t [nTot];
1735       for (Int_t i=0; i<nTot; i++) coupl3[i] = 9999;
1736     } // else
1737
1738     for (Int_t i=0; i<nCoupled; i++) {
1739       indx1 = clustNumb[i];
1740       for (Int_t j=i+1; j<nCoupled; j++) {
1741         indx2 = clustNumb[j];
1742         if (i123 == 1) {
1743           coupl1[i] += (*aijcluclu)(indx1,indx2);
1744           coupl1[j] += (*aijcluclu)(indx1,indx2);
1745         }
1746         else if (i123 == 2) {
1747           indx = i*nCoupled + j;
1748           coupl2[indx] = coupl1[i] + coupl1[j];
1749           coupl2[indx] -= 2 * ((*aijcluclu)(indx1,indx2));
1750         } else {
1751           for (Int_t k=j+1; k<nCoupled; k++) {
1752             indx3 = clustNumb[k];
1753             indx = i*nCoupled*nCoupled + j*nCoupled + k;
1754             coupl3[indx] = coupl2[i*nCoupled+j] + coupl1[k];
1755             coupl3[indx] -= 2 * ((*aijcluclu)(indx1,indx3)+(*aijcluclu)(indx2,indx3));
1756           }
1757         } // else
1758       } // for (Int_t j=i+1;
1759     } // for (Int_t i=0;
1760   } // for (Int_t i123=1;
1761
1762   // Find minimum coupling
1763   Double_t couplMin = 9999;
1764   Int_t locMin = 0;
1765
1766   for (Int_t i123=1; i123<=i123max; i123++) {
1767     if (i123 == 1) {
1768       locMin = TMath::LocMin(nCoupled, coupl1);
1769       couplMin = coupl1[locMin];
1770       minGroup[0] = locMin;
1771       delete [] coupl1; coupl1 = 0;
1772     }
1773     else if (i123 == 2) {
1774       locMin = TMath::LocMin(nCoupled*nCoupled, coupl2);
1775       if (coupl2[locMin] < couplMin) {
1776         couplMin = coupl2[locMin];
1777         minGroup[0] = locMin/nCoupled;
1778         minGroup[1] = locMin%nCoupled;
1779       }
1780       delete [] coupl2; coupl2 = 0;
1781     } else {
1782       locMin = TMath::LocMin(nTot, coupl3);
1783       if (coupl3[locMin] < couplMin) {
1784         couplMin = coupl3[locMin];
1785         minGroup[0] = locMin/nCoupled/nCoupled;
1786         minGroup[1] = locMin%(nCoupled*nCoupled)/nCoupled;
1787         minGroup[2] = locMin%nCoupled;
1788       }
1789       delete [] coupl3; coupl3 = 0;
1790     } // else
1791   } // for (Int_t i123=1;
1792   return couplMin;
1793 }
1794
1795 //_____________________________________________________________________________
1796 Int_t AliMUONClusterFinderAZ::SelectPad(Int_t nCoupled, Int_t nForFit, Int_t *clustNumb, Int_t *clustFit, TMatrixD *aijclupad)
1797 {
1798 /// Select pads for fit. If too many coupled clusters, find pads giving
1799 /// the strongest coupling with the rest of clusters and exclude them from the fit.
1800
1801   Int_t npad = fnPads[0] + fnPads[1];
1802   Double_t *padpix = 0;
1803
1804   if (nCoupled > 3) {
1805     padpix = new Double_t[npad];
1806     for (Int_t i=0; i<npad; i++) padpix[i] = 0;
1807   }
1808
1809   Int_t nOK = 0, indx, indx1;
1810   for (Int_t iclust=0; iclust<nForFit; iclust++) {
1811     indx = clustFit[iclust];
1812     for (Int_t j=0; j<npad; j++) {
1813       if ((*aijclupad)(indx,j) < fgkCouplMin) continue;
1814       if (fPadIJ[1][j] == -5) fPadIJ[1][j] = -9; // flag overflow
1815       if (fPadIJ[1][j] < 0) continue; // exclude overflows and used pads
1816       if (!fPadIJ[1][j]) { fPadIJ[1][j] = 1; nOK++; } // pad to be used in fit
1817       if (nCoupled > 3) {
1818         // Check other clusters
1819         for (Int_t iclust1=0; iclust1<nCoupled; iclust1++) {
1820           indx1 = clustNumb[iclust1];
1821           if (indx1 < 0) continue;
1822           if ((*aijclupad)(indx1,j) < fgkCouplMin) continue;
1823           padpix[j] += (*aijclupad)(indx1,j);
1824         }
1825       } // if (nCoupled > 3)
1826     } // for (Int_t j=0; j<npad;
1827   } // for (Int_t iclust=0; iclust<nForFit
1828   if (nCoupled < 4) return nOK;
1829
1830   Double_t aaa = 0;
1831   for (Int_t j=0; j<npad; j++) {
1832     if (padpix[j] < fgkCouplMin) continue;
1833     if (fDebug) cout << j << " " << padpix[j] << " " << fXyq[0][j] << " " << fXyq[1][j] << endl;
1834     aaa += padpix[j];
1835     fPadIJ[1][j] = -1; // exclude pads with strong coupling to the other clusters
1836     nOK--;
1837   }
1838   delete [] padpix; padpix = 0;
1839   return nOK;
1840 }
1841
1842 //_____________________________________________________________________________
1843 void AliMUONClusterFinderAZ::Merge(Int_t nForFit, Int_t nCoupled, Int_t *clustNumb, Int_t *clustFit, TObjArray **clusters, TMatrixD *aijcluclu, TMatrixD *aijclupad)
1844 {
1845 /// Merge the group of clusters with the one having the strongest coupling with them
1846
1847   Int_t indx, indx1, npxclu, npxclu1, imax=0;
1848   TObjArray *pix, *pix1;
1849   Double_t couplMax;
1850
1851   for (Int_t icl=0; icl<nForFit; icl++) {
1852     indx = clustFit[icl];
1853     pix = clusters[indx];
1854     npxclu = pix->GetEntriesFast();
1855     couplMax = -1;
1856     for (Int_t icl1=0; icl1<nCoupled; icl1++) {
1857       indx1 = clustNumb[icl1];
1858       if (indx1 < 0) continue;
1859       if ((*aijcluclu)(indx,indx1) > couplMax) {
1860         couplMax = (*aijcluclu)(indx,indx1);
1861         imax = indx1;
1862       }
1863     } // for (Int_t icl1=0;
1864     /*if (couplMax < fgkCouplMin) {
1865       cout << " Oops " << couplMax << endl;
1866       aijcluclu->Print();
1867       cout << icl << " " << indx << " " << npxclu << " " << nLinks << endl;
1868       ::exit(0);
1869       }*/
1870     // Add to it
1871     pix1 = clusters[imax];
1872     npxclu1 = pix1->GetEntriesFast();
1873     // Add pixels
1874     for (Int_t i=0; i<npxclu; i++) { pix1->Add(pix->UncheckedAt(i)); pix->RemoveAt(i); }
1875     if (fDebug) cout << " New number of pixels: " << npxclu1 << " " << pix1->GetEntriesFast() << endl;
1876     //Add cluster-to-cluster couplings
1877     //aijcluclu->Print();
1878     for (Int_t icl1=0; icl1<nCoupled; icl1++) {
1879       indx1 = clustNumb[icl1];
1880       if (indx1 < 0 || indx1 == imax) continue;
1881       (*aijcluclu)(indx1,imax) += (*aijcluclu)(indx,indx1);
1882       (*aijcluclu)(imax,indx1) = (*aijcluclu)(indx1,imax);
1883     }
1884     (*aijcluclu)(indx,imax) = (*aijcluclu)(imax,indx) = 0;
1885     //aijcluclu->Print();
1886     //Add cluster-to-pad couplings
1887     for (Int_t j=0; j<fnPads[0]+fnPads[1]; j++) {
1888       if (fPadIJ[1][j] < 0 && fPadIJ[1][j] != -5) continue; // exclude used pads
1889       (*aijclupad)(imax,j) += (*aijclupad)(indx,j);
1890       (*aijclupad)(indx,j) = 0;
1891     }
1892   } // for (Int_t icl=0; icl<nForFit;
1893 }
1894
1895 //_____________________________________________________________________________
1896 Int_t AliMUONClusterFinderAZ::Fit(Int_t iSimple, Int_t nfit, Int_t *clustFit, TObjArray **clusters, Double_t *parOk)
1897 {
1898 /// Find selected clusters to selected pad charges
1899
1900   TH2D *mlem = (TH2D*) gROOT->FindObject("mlem");
1901   Double_t xmin = mlem->GetXaxis()->GetXmin() - mlem->GetXaxis()->GetBinWidth(1);
1902   Double_t xmax = mlem->GetXaxis()->GetXmax() + mlem->GetXaxis()->GetBinWidth(1);
1903   Double_t ymin = mlem->GetYaxis()->GetXmin() - mlem->GetYaxis()->GetBinWidth(1);
1904   Double_t ymax = mlem->GetYaxis()->GetXmax() + mlem->GetYaxis()->GetBinWidth(1);
1905   Double_t step[3]={0.01,0.002,0.02}, xPad = 0, yPad = 99999;
1906
1907   // Number of pads to use and number of virtual pads
1908   Int_t npads = 0, nVirtual = 0, nfit0 = nfit;
1909   for (Int_t i=0; i<fnPads[0]+fnPads[1]; i++) {
1910     if (fXyq[3][i] < 0) nVirtual++;
1911     if (fPadIJ[1][i] != 1) continue;
1912     if (fXyq[3][i] > 0) {
1913       npads++;
1914       if (yPad > 9999) {
1915         xPad = fXyq[0][i];
1916         yPad = fXyq[1][i];
1917       } else {
1918         if (fXyq[4][i] < fXyq[3][i]) yPad = fXyq[1][i];
1919         else xPad = fXyq[0][i];
1920       }
1921     }
1922   }
1923   if (fDebug) {
1924     for (Int_t i=0; i<nfit; i++) {cout << i+1 << " " << clustFit[i] << " ";}
1925     cout << nfit << endl;
1926     cout << " Number of pads to fit: " << npads << endl;
1927   }
1928   fNpar = 0;
1929   fQtot = 0;
1930   if (npads < 2) return 0;
1931
1932   AliMUONVDigit *mdig = 0;
1933   Int_t tracks[3] = {-1, -1, -1};
1934   for (Int_t cath=0; cath<2; cath++) {
1935     for (Int_t i=0; i<fnPads[0]+fnPads[1]; i++) {
1936       if (fPadIJ[0][i] != cath) continue;
1937       if (fPadIJ[1][i] != 1) continue;
1938       if (fXyq[3][i] < 0) continue; // exclude virtual pads
1939       UInt_t digit = fDigitId[i];
1940       mdig = static_cast<AliMUONVDigit*>(fDigitStore->FindObject(digit));
1941       if (!mdig) continue; // protection for cluster display
1942       if (mdig->Hit() >= 0) {
1943         if (tracks[0] < 0) {
1944           tracks[0] = mdig->Hit();
1945           tracks[1] = mdig->Track(0);
1946         } else if (mdig->Track(0) < tracks[1]) {
1947           tracks[0] = mdig->Hit();
1948           tracks[1] = mdig->Track(0);
1949         }
1950       }
1951       if (mdig->Track(1) >= 0 && mdig->Track(1) != tracks[1]) {
1952         if (tracks[2] < 0) tracks[2] = mdig->Track(1);
1953         else tracks[2] = TMath::Min (tracks[2], mdig->Track(1));
1954       }
1955       //if (!mdig) break;
1956       //cout << mdig->Hit() << " " << mdig->Track(0) << " " << mdig->Track(1) <<endl;
1957     } // for (Int_t i=0;
1958   } // for (Int_t cath=0;
1959   //cout << tracks[0] << " " << tracks[1] << " " << tracks[2] <<endl;
1960
1961   // Get number of pads in X and Y
1962   Int_t nInX = 0, nInY;
1963   PadsInXandY(nInX, nInY);
1964   //cout << " nInX and Y: " << nInX << " " << nInY << endl;
1965
1966   Int_t nfitMax = 3;
1967   nfitMax = TMath::Min (nfitMax, (npads + 1) / 3);
1968   if (nfitMax > 1) {
1969     if (nInX < 3 && nInY < 3 || nInX == 3 && nInY < 3 || nInX < 3 && nInY == 3) nfitMax = 1; // not enough pads in each direction
1970   }
1971   if (nfit > nfitMax) nfit = nfitMax;
1972
1973   // Take cluster maxima as fitting seeds
1974   TObjArray *pix;
1975   AliMUONPixel *pixPtr;
1976   Int_t npxclu;
1977   Double_t cont, cmax = 0, xseed = 0, yseed = 0, errOk[8], qq = 0;
1978   Double_t xyseed[3][2], qseed[3], xyCand[3][2] = {{0},{0}}, sigCand[3][2] = {{0},{0}};
1979
1980   for (Int_t ifit=1; ifit<=nfit0; ifit++) {
1981     cmax = 0;
1982     pix = clusters[clustFit[ifit-1]];
1983     npxclu = pix->GetEntriesFast();
1984     //qq = 0;
1985     for (Int_t clu=0; clu<npxclu; clu++) {
1986       pixPtr = (AliMUONPixel*) pix->UncheckedAt(clu);
1987       cont = pixPtr->Charge();
1988       fQtot += cont;
1989       if (cont > cmax) {
1990         cmax = cont;
1991         xseed = pixPtr->Coord(0);
1992         yseed = pixPtr->Coord(1);
1993       }
1994       qq += cont;
1995       /*
1996       xyCand[ifit-1][0] += pixPtr->Coord(0) * cont;
1997       xyCand[ifit-1][1] += pixPtr->Coord(1) * cont;
1998       sigCand[ifit-1][0] += pixPtr->Coord(0) * pixPtr->Coord(0) * cont;
1999       sigCand[ifit-1][1] += pixPtr->Coord(1) * pixPtr->Coord(1) * cont;
2000       */
2001       xyCand[0][0] += pixPtr->Coord(0) * cont;
2002       xyCand[0][1] += pixPtr->Coord(1) * cont;
2003       sigCand[0][0] += pixPtr->Coord(0) * pixPtr->Coord(0) * cont;
2004       sigCand[0][1] += pixPtr->Coord(1) * pixPtr->Coord(1) * cont;
2005     }
2006     xyseed[ifit-1][0] = xseed;
2007     xyseed[ifit-1][1] = yseed;
2008     qseed[ifit-1] = cmax;
2009     /*
2010     xyCand[ifit-1][0] /= qq; // <x>
2011     xyCand[ifit-1][1] /= qq; // <y>
2012     sigCand[ifit-1][0] = sigCand[ifit-1][0]/qq - xyCand[ifit-1][0]*xyCand[ifit-1][0]; // <x^2> - <x>^2
2013     sigCand[ifit-1][0] = sigCand[ifit-1][0] > 0 ? TMath::Sqrt (sigCand[ifit-1][0]) : 0;
2014     sigCand[ifit-1][1] = sigCand[ifit-1][1]/qq - xyCand[ifit-1][1]*xyCand[ifit-1][1]; // <y^2> - <y>^2
2015     sigCand[ifit-1][1] = sigCand[ifit-1][1] > 0 ? TMath::Sqrt (sigCand[ifit-1][1]) : 0;
2016     cout << xyCand[ifit-1][0] << " " << xyCand[ifit-1][1] << " " << sigCand[ifit-1][0] << " " << sigCand[ifit-1][1] << endl;
2017     */
2018   } // for (Int_t ifit=1;
2019
2020   xyCand[0][0] /= qq; // <x>
2021   xyCand[0][1] /= qq; // <y>
2022   sigCand[0][0] = sigCand[0][0]/qq - xyCand[0][0]*xyCand[0][0]; // <x^2> - <x>^2
2023   sigCand[0][0] = sigCand[0][0] > 0 ? TMath::Sqrt (sigCand[0][0]) : 0;
2024   sigCand[0][1] = sigCand[0][1]/qq - xyCand[0][1]*xyCand[0][1]; // <y^2> - <y>^2
2025   sigCand[0][1] = sigCand[0][1] > 0 ? TMath::Sqrt (sigCand[0][1]) : 0;
2026   if (fDebug) cout << xyCand[0][0] << " " << xyCand[0][1] << " " << sigCand[0][0] << " " << sigCand[0][1] << endl;
2027
2028   Int_t nDof, maxSeed[3], nMax = 0;
2029   Double_t fmin, chi2o = 9999, chi2n;
2030
2031   TMath::Sort(nfit0, qseed, maxSeed, kTRUE); // in decreasing order
2032   /*
2033   Int_t itmp[100], localMax[100];
2034   Double_t maxVal[100];
2035   if (!iSimple && nfit < nfitMax) {
2036     // Try to split pixel cluster according to local maxima
2037     Int_t nfit1 = nfit;
2038     for (Int_t iclus = 0; iclus < nfit1; iclus++) {
2039       nMax = FindLocalMaxima (clusters[clustFit[maxSeed[iclus]]], localMax, maxVal);
2040       TH2D *hist = (TH2D*) gROOT->FindObject("anode1");
2041       if (nMax == 1) { hist->Delete(); continue; }
2042       // Add extra fitting seeds from local maxima
2043       Int_t ixseed = hist->GetXaxis()->FindBin(xyseed[maxSeed[iclus]][0]);
2044       Int_t iyseed = hist->GetYaxis()->FindBin(xyseed[maxSeed[iclus]][1]);
2045       Int_t nx = hist->GetNbinsX();
2046       TMath::Sort(nMax, maxVal, itmp, kTRUE); // in decreasing order
2047       for (Int_t j = 0; j < nMax; j++) {
2048         Int_t iyc = localMax[itmp[j]] / nx + 1;
2049         Int_t ixc = localMax[itmp[j]] % nx + 1;
2050         if (ixc == ixseed && iyc == iyseed) continue; // local max already taken for seeding
2051         xyseed[nfit][0] = hist->GetXaxis()->GetBinCenter(ixc);
2052         xyseed[nfit][1] = hist->GetYaxis()->GetBinCenter(iyc);
2053         qseed[nfit] = maxVal[itmp[j]];
2054         maxSeed[nfit] = nfit++;
2055         if (nfit >= nfitMax) break;
2056       }
2057       hist->Delete();
2058       if (nfit >= nfitMax) break;
2059     } // for (Int_t iclus = 0;
2060     //nfit0 = nfit;
2061     //TMath::Sort(nfit0, qseed, maxSeed, kTRUE); // in decreasing order
2062   } //if (!iSimple && nfit < nfitMax)
2063   */
2064
2065   Double_t *gin = 0, func0, func1, param[8], step0[8];
2066   Double_t param0[2][8]={{0},{0}}, deriv[2][8]={{0},{0}};
2067   Double_t shift[8], stepMax, derMax, parmin[8], parmax[8], func2[2], shift0;
2068   Double_t delta[8], scMax, dder[8], estim, shiftSave = 0;
2069   Int_t min, max, nCall = 0, memory[8] = {0}, nLoop, idMax = 0, iestMax = 0, nFail;
2070   Double_t rad, dist[3] = {0};
2071
2072   // Try to fit with one-track hypothesis, then 2-track. If chi2/dof is
2073   // lower, try 3-track (if number of pads is sufficient).
2074   for (Int_t iseed=0; iseed<nfit; iseed++) {
2075
2076     if (iseed) { for (Int_t j=0; j<fNpar; j++) param[j] = parOk[j]; } // for bounded params
2077     for (Int_t j=0; j<3; j++) step0[fNpar+j] = shift[fNpar+j] = step[j];
2078     if (nfit == 1) param[fNpar] = xyCand[0][0]; // take COG
2079     else param[fNpar] = xyseed[maxSeed[iseed]][0];
2080     parmin[fNpar] = xmin;
2081     parmax[fNpar++] = xmax;
2082     if (nfit == 1) param[fNpar] = xyCand[0][1]; // take COG
2083     else param[fNpar] = xyseed[maxSeed[iseed]][1];
2084     parmin[fNpar] = ymin;
2085     parmax[fNpar++] = ymax;
2086     if (fNpar > 2) {
2087       param[fNpar] = fNpar == 4 ? 0.5 : 0.3;
2088       parmin[fNpar] = 0;
2089       parmax[fNpar++] = 1;
2090     }
2091     if (iseed) { for (Int_t j=0; j<fNpar; j++) param0[1][j] = 0; }
2092
2093     // Try new algorithm
2094     min = nLoop = 1; stepMax = func2[1] = derMax = 999999; nFail = 0;
2095
2096     while (1) {
2097       max = !min;
2098       Fcn1(fNpar, gin, func0, param, 1); nCall++;
2099       //cout << " Func: " << func0 << endl;
2100
2101       func2[max] = func0;
2102       for (Int_t j=0; j<fNpar; j++) {
2103         param0[max][j] = param[j];
2104         delta[j] = step0[j];
2105         param[j] += delta[j] / 10;
2106         if (j > 0) param[j-1] -= delta[j-1] / 10;
2107         Fcn1(fNpar, gin, func1, param, 1); nCall++;
2108         deriv[max][j] = (func1 - func0) / delta[j] * 10; // first derivative
2109         //cout << j << " " << deriv[max][j] << endl;
2110         dder[j] = param0[0][j] != param0[1][j] ? (deriv[0][j] - deriv[1][j]) /
2111           (param0[0][j] - param0[1][j]) : 0; // second derivative
2112       }
2113       param[fNpar-1] -= delta[fNpar-1] / 10;
2114       if (nCall > 2000) break;
2115
2116       min = func2[0] < func2[1] ? 0 : 1;
2117       nFail = min == max ? 0 : nFail + 1;
2118
2119       stepMax = derMax = estim = 0;
2120       for (Int_t j=0; j<fNpar; j++) {
2121         // Estimated distance to minimum
2122         shift0 = shift[j];
2123         if (nLoop == 1) shift[j] = TMath::Sign (step0[j], -deriv[max][j]); // first step
2124         else if (TMath::Abs(deriv[0][j]) < 1.e-3 && TMath::Abs(deriv[1][j]) < 1.e-3) shift[j] = 0;
2125         else if (deriv[min][j]*deriv[!min][j] > 0 && TMath::Abs(deriv[min][j]) > TMath::Abs(deriv[!min][j])
2126                  //|| TMath::Abs(deriv[0][j]-deriv[1][j]) < 1.e-3) {
2127           || TMath::Abs(deriv[0][j]-deriv[1][j]) < 1.e-3 || TMath::Abs(dder[j]) < 1.e-6) {
2128           shift[j] = -TMath::Sign (shift[j], (func2[0]-func2[1]) * (param0[0][j]-param0[1][j]));
2129           if (min == max) {
2130             if (memory[j] > 1) { shift[j] *= 2; } //cout << " Memory " << memory[j] << " " << shift[j] << endl; }
2131             memory[j]++;
2132           }
2133         } else {
2134           shift[j] = dder[j] != 0 ? -deriv[min][j] / dder[j] : 0;
2135           memory[j] = 0;
2136         }
2137         if (TMath::Abs(shift[j])/step0[j] > estim) {
2138           estim = TMath::Abs(shift[j])/step0[j];
2139           iestMax = j;
2140         }
2141
2142         // Too big step
2143         if (TMath::Abs(shift[j])/step0[j] > 10) shift[j] = TMath::Sign(10.,shift[j]) * step0[j]; //
2144
2145         // Failed to improve minimum
2146         if (min != max) {
2147           memory[j] = 0;
2148           param[j] = param0[min][j];
2149           if (TMath::Abs(shift[j]+shift0) > 0.1*step0[j]) shift[j] = (shift[j] + shift0) / 2;
2150           else shift[j] /= -2;
2151         }
2152
2153         // Too big step
2154         if (TMath::Abs(shift[j]*deriv[min][j]) > func2[min])
2155           shift[j] = TMath::Sign (func2[min]/deriv[min][j], shift[j]);
2156
2157         // Introduce step relaxation factor
2158         if (memory[j] < 3) {
2159           scMax = 1 + 4 / TMath::Max(nLoop/2.,1.);
2160           if (TMath::Abs(shift0) > 0 && TMath::Abs(shift[j]/shift0) > scMax)
2161             shift[j] = TMath::Sign (shift0*scMax, shift[j]);
2162         }
2163         param[j] += shift[j];
2164         //AZ Check parameter limits 27-12-2004
2165         if (param[j] < parmin[j]) {
2166           shift[j] = parmin[j] - param[j];
2167           param[j] = parmin[j];
2168         } else if (param[j] > parmax[j]) {
2169           shift[j] = parmax[j] - param[j];
2170           param[j] = parmax[j];
2171         }
2172         //cout << " xxx " << j << " " << shift[j] << " " << param[j] << endl;
2173         stepMax = TMath::Max (stepMax, TMath::Abs(shift[j]/step0[j]));
2174         if (TMath::Abs(deriv[min][j]) > derMax) {
2175           idMax = j;
2176           derMax = TMath::Abs (deriv[min][j]);
2177         }
2178       } // for (Int_t j=0; j<fNpar;
2179       //cout << max << " " << func2[min] << " " << derMax << " " << stepMax << " " << estim << " " << iestMax << " " << nCall << endl;
2180       if (estim < 1 && derMax < 2 || nLoop > 150) break; // minimum was found
2181
2182       nLoop++;
2183       // Check for small step
2184       if (shift[idMax] == 0) { shift[idMax] = step0[idMax]/10; param[idMax] += shift[idMax]; continue; }
2185       if (!memory[idMax] && derMax > 0.5 && nLoop > 10) {
2186         //cout << " ok " << deriv[min][idMax] << " " << deriv[!min][idMax] << " " << dder[idMax]*shift[idMax] << " " << shift[idMax] << endl;
2187         if (dder[idMax] != 0 && TMath::Abs(deriv[min][idMax]/dder[idMax]/shift[idMax]) > 10) {
2188           if (min == max) dder[idMax] = -dder[idMax];
2189           shift[idMax] = -deriv[min][idMax] / dder[idMax] / 10;
2190           param[idMax] += shift[idMax];
2191           stepMax = TMath::Max (stepMax, TMath::Abs(shift[idMax])/step0[idMax]);
2192           //cout << shift[idMax] << " " << param[idMax] << endl;
2193           if (min == max) shiftSave = shift[idMax];
2194         }
2195         if (nFail > 10) {
2196           param[idMax] -= shift[idMax];
2197           shift[idMax] = 4 * shiftSave * (gRandom->Rndm(0) - 0.5);
2198           param[idMax] += shift[idMax];
2199           //cout << shift[idMax] << endl;
2200         }
2201       }
2202     } // while (1)
2203     fmin = func2[min];
2204
2205     nDof = npads - fNpar + nVirtual;
2206     if (!nDof) nDof++;
2207     chi2n = fmin / nDof;
2208     if (fDebug) cout << " Chi2 " << chi2n << " " << fNpar << endl;
2209
2210     if (chi2n*1.2+1.e-6 > chi2o ) { fNpar -= 3; break; }
2211
2212     // Save parameters and errors
2213
2214     if (nInX == 1) {
2215       // One pad per direction
2216       for (Int_t i=0; i<fNpar; i++) if (i == 0 || i == 2 || i == 5) param0[min][i] = xPad;
2217     }
2218     if (nInY == 1) {
2219       // One pad per direction
2220       for (Int_t i=0; i<fNpar; i++) if (i == 1 || i == 3 || i == 6) param0[min][i] = yPad;
2221     }
2222
2223     /*
2224     if (iseed > 0) {
2225       // Find distance to the nearest neighbour
2226       dist[0] = dist[1] = TMath::Sqrt ((param0[min][0]-param0[min][2])*
2227                                        (param0[min][0]-param0[min][2])
2228                                       +(param0[min][1]-param0[min][3])*
2229                                        (param0[min][1]-param0[min][3]));
2230       if (iseed > 1) {
2231         dist[2] = TMath::Sqrt ((param0[min][0]-param0[min][5])*
2232                                (param0[min][0]-param0[min][5])
2233                               +(param0[min][1]-param0[min][6])*
2234                                (param0[min][1]-param0[min][6]));
2235         rad = TMath::Sqrt ((param0[min][2]-param0[min][5])*
2236                            (param0[min][2]-param0[min][5])
2237                           +(param0[min][3]-param0[min][6])*
2238                            (param0[min][3]-param0[min][6]));
2239         if (dist[2] < dist[0]) dist[0] = dist[2];
2240         if (rad < dist[1]) dist[1] = rad;
2241         if (rad < dist[2]) dist[2] = rad;
2242       }
2243       cout << dist[0] << " " << dist[1] << " " << dist[2] << endl;
2244       if (dist[TMath::LocMin(iseed+1,dist)] < 1.) { fNpar -= 3; break; }
2245     }
2246     */
2247
2248     for (Int_t i=0; i<fNpar; i++) {
2249       parOk[i] = param0[min][i];
2250       //errOk[i] = fmin;
2251       errOk[i] = chi2n;
2252       // Bounded params
2253       parOk[i] = TMath::Max (parOk[i], parmin[i]);
2254       parOk[i] = TMath::Min (parOk[i], parmax[i]);
2255     }
2256
2257     chi2o = chi2n;
2258     if (fmin < 0.1) break; // !!!???
2259   } // for (Int_t iseed=0;
2260
2261   if (fDebug) {
2262     for (Int_t i=0; i<fNpar; i++) {
2263       if (i == 4 || i == 7) {
2264         if (i == 7 || i == 4 && fNpar < 7) cout << parOk[i] << endl;
2265         else cout << parOk[i] * (1-parOk[7]) << endl;
2266         continue;
2267       }
2268       cout << parOk[i] << " " << errOk[i] << endl;
2269     }
2270   }
2271   nfit = (fNpar + 1) / 3;
2272   dist[0] = dist[1] = dist[2] = 0;
2273
2274   if (nfit > 1) {
2275     // Find distance to the nearest neighbour
2276     dist[0] = dist[1] = TMath::Sqrt ((parOk[0]-parOk[2])*
2277                                      (parOk[0]-parOk[2])
2278                                     +(parOk[1]-parOk[3])*
2279                                      (parOk[1]-parOk[3]));
2280     if (nfit > 2) {
2281       dist[2] = TMath::Sqrt ((parOk[0]-parOk[5])*
2282                              (parOk[0]-parOk[5])
2283                             +(parOk[1]-parOk[6])*
2284                              (parOk[1]-parOk[6]));
2285       rad = TMath::Sqrt ((parOk[2]-parOk[5])*
2286                          (parOk[2]-parOk[5])
2287                         +(parOk[3]-parOk[6])*
2288                          (parOk[3]-parOk[6]));
2289       if (dist[2] < dist[0]) dist[0] = dist[2];
2290       if (rad < dist[1]) dist[1] = rad;
2291       if (rad < dist[2]) dist[2] = rad;
2292     }
2293   }
2294
2295   Int_t indx;
2296   fnPads[1] -= nVirtual;
2297 //  if (!fDraw) {
2298     Double_t coef = 0;
2299     if (iSimple) fnCoupled = 0;
2300     //for (Int_t j=0; j<nfit; j++) {
2301     for (Int_t j=nfit-1; j>=0; j--) {
2302       indx = j<2 ? j*2 : j*2+1;
2303       if (nfit == 1) coef = 1;
2304       else coef = j==nfit-1 ? parOk[indx+2] : 1-coef;
2305       coef = TMath::Max (coef, 0.);
2306       if (nfit == 3 && j < 2) coef = j==1 ? coef*parOk[indx+2] : coef - parOk[7];
2307       coef = TMath::Max (coef, 0.);
2308       AddRawCluster (parOk[indx], parOk[indx+1], coef*fQtot, errOk[indx], nfit0+10*nfit+100*nMax+10000*fnCoupled, tracks,
2309                      //sigCand[maxSeed[j]][0], sigCand[maxSeed[j]][1]);
2310                      //sigCand[0][0], sigCand[0][1], dist[j]);
2311                      sigCand[0][0], sigCand[0][1], dist[TMath::LocMin(nfit,dist)]);
2312     }
2313 //  } else fDraw->FillMuon(nfit, parOk, errOk);
2314   return nfit;
2315 }
2316
2317 //_____________________________________________________________________________
2318 void AliMUONClusterFinderAZ::Fcn1(Int_t & /*npar*/, Double_t * /*gin*/, Double_t &f, Double_t *par, Int_t /*iflag*/)
2319 {
2320 /// Fit for one track
2321 /// AZ for Muinuit AliMUONClusterFinderAZ& c = *(AliMUONClusterFinderAZ::fgClusterFinder);
2322
2323   AliMUONClusterFinderAZ& c = *this; //AZ
2324
2325   Int_t cath, ix, iy, indx, npads=0;
2326   Double_t charge, delta, coef=0, chi2=0, qTot = 0;
2327   for (Int_t j=0; j<c.fnPads[0]+c.fnPads[1]; j++) {
2328     if (c.fPadIJ[1][j] != 1) continue;
2329     cath = c.fPadIJ[0][j];
2330     if (c.fXyq[3][j] > 0) npads++; // exclude virtual pads
2331     qTot += c.fXyq[2][j];
2332     ix = c.fPadIJ[2][j];
2333     iy = c.fPadIJ[3][j];
2334 //    c.fSegmentation[cath]->SetPad(ix, iy);
2335     charge = 0;
2336     for (Int_t i=c.fNpar/3; i>=0; i--) { // sum over tracks
2337       indx = i<2 ? 2*i : 2*i+1;
2338 //      c.fSegmentation[cath]->SetHit(par[indx], par[indx+1], c.fZpad);
2339       if (c.fNpar == 2) coef = 1;
2340       else coef = i==c.fNpar/3 ? par[indx+2] : 1-coef;
2341       coef = TMath::Max (coef, 0.);
2342       if (c.fNpar == 8 && i < 2) coef = i==1 ? coef*par[indx+2] : coef - par[7];
2343       coef = TMath::Max (coef, 0.);
2344 //      charge += fMathieson->IntXY(fDetElemId, fSegmentation[cath])*coef;
2345       charge += ChargeIntegration(par[indx],par[indx+1],
2346                                   c.fXyq[0][j],c.fXyq[1][j],
2347                                   c.fXyq[3][j],c.fXyq[4][j]);
2348     }
2349     charge *= c.fQtot;
2350     delta = charge - c.fXyq[2][j];
2351     delta *= delta;
2352     delta /= c.fXyq[2][j];
2353     //if (cath) delta /= 5; // just for test
2354     chi2 += delta;
2355   } // for (Int_t j=0;
2356   f = chi2;
2357   Double_t qAver = qTot/npads; //(c.fnPads[0]+c.fnPads[1]);
2358   f = chi2/qAver;
2359 }
2360
2361 //_____________________________________________________________________________
2362 void AliMUONClusterFinderAZ::UpdatePads(Int_t /*nfit*/, Double_t *par)
2363 {
2364 /// Subtract the fitted charges from pads with strong coupling
2365
2366   Int_t cath, ix, iy, indx;
2367   Double_t charge, coef=0;
2368   for (Int_t j=0; j<fnPads[0]+fnPads[1]; j++) {
2369     if (fPadIJ[1][j] != -1) continue;
2370     if (fNpar != 0) {
2371       cath = fPadIJ[0][j];
2372       ix = fPadIJ[2][j];
2373       iy = fPadIJ[3][j];
2374       //      fSegmentation[cath]->SetPad(ix, iy);
2375       charge = 0;
2376       for (Int_t i=fNpar/3; i>=0; i--) { // sum over tracks
2377         indx = i<2 ? 2*i : 2*i+1;
2378 //        fSegmentation[cath]->SetHit(par[indx], par[indx+1], fZpad);
2379         if (fNpar == 2) coef = 1;
2380         else coef = i==fNpar/3 ? par[indx+2] : 1-coef;
2381         coef = TMath::Max (coef, 0.);
2382         if (fNpar == 8 && i < 2) coef = i==1 ? coef*par[indx+2] : coef - par[7];
2383         coef = TMath::Max (coef, 0.);
2384 //        charge += fMathieson->IntXY(fDetElemId,fSegmentation[cath])*coef;
2385         charge += ChargeIntegration(par[indx],par[indx+1],
2386                                     fXyq[0][j],fXyq[1][j],
2387                                     fXyq[3][j],fXyq[4][j]);
2388       }
2389       charge *= fQtot;
2390       fXyq[2][j] -= charge;
2391     } // if (fNpar != 0)
2392     if (fXyq[2][j] > fgkZeroSuppression) fPadIJ[1][j] = 0; // return pad for further using
2393   } // for (Int_t j=0;
2394 }
2395
2396 //_____________________________________________________________________________
2397 Bool_t AliMUONClusterFinderAZ::TestTrack(Int_t /*t*/) const
2398 {
2399 /// Test if track was user selected
2400
2401   return kTRUE;
2402   /*
2403     if (fTrack[0]==-1 || fTrack[1]==-1) {
2404         return kTRUE;
2405     } else if (t==fTrack[0] || t==fTrack[1]) {
2406         return kTRUE;
2407     } else {
2408         return kFALSE;
2409     }
2410   */
2411 }
2412
2413 //_____________________________________________________________________________
2414 void AliMUONClusterFinderAZ::AddRawCluster(Double_t x, Double_t y,
2415                                            Double_t qTot,
2416                                            Double_t /*fmin*/, Int_t /*nfit*/,
2417                                            Int_t* /*tracks*/,
2418                                            Double_t /*sigx*/, Double_t /*sigy*/,
2419                                            Double_t /*dist*/)
2420 {
2421 /// Add a raw cluster copy to the list
2422
2423   if (qTot <= 0.501) return;
2424
2425 //  Int_t cath, npads[2] = {0}, nover[2] = {0};
2426 //  for (Int_t j=0; j<fnPads[0]+fnPads[1]; j++)
2427 //  {
2428 //    cath = fPadIJ[0][j];
2429 //    // There was an overflow
2430 //    if (fPadIJ[1][j] == -9) nover[cath]++;
2431 //    if (fPadIJ[1][j] != 1 && fPadIJ[1][j] != -9) continue;
2432 //    cnew.SetMultiplicity(cath,cnew.GetMultiplicity(cath)+1);
2433 //    if (fXyq[2][j] > cnew.GetPeakSignal(cath)) cnew.SetPeakSignal(cath,fXyq[2][j]);
2434 //    //cnew.SetCharge(cath,cnew.GetCharge(cath) + TMath::Nint (fXyq[2][j]));
2435 //    cnew.SetContrib(npads[cath],cath,fXyq[2][j]);
2436 //    cnew.SetIndex(npads[cath],cath,TMath::Nint (fXyq[5][j]));
2437 //    cnew.SetDetElemId(fDetElemId);
2438 //    npads[cath]++;
2439 //  }
2440
2441 //  cnew.SetClusterType(nover[0] + nover[1] * 100);
2442 //  for (Int_t j=0; j<3; j++) cnew.SetTrack(j,tracks[j]);
2443
2444 //  Double_t xg, yg, zg;
2445 //  for (cath=0; cath<2; cath++)
2446 //  {
2447 //    // Perform local-to-global transformation
2448 //    cnew.SetX(cath, xg);
2449 //    cnew.SetY(cath, yg);
2450 //    cnew.SetZ(cath, zg);
2451 //    cnew.SetCharge(cath, TMath::Nint(qTot));
2452 //    //cnew.SetPeakSignal(cath,20);
2453 //    //cnew.SetMultiplicity(cath, 5);
2454 //    cnew.SetNcluster(cath, nfit);
2455 //    cnew.SetChi2(cath, fmin); //0.;1
2456 //  }
2457   // Evaluate measurement errors
2458   //AZ Errors(&cnew);
2459
2460   AliMUONCluster cnew;
2461
2462   cnew.SetCharge(qTot,qTot);
2463   cnew.SetPosition(TVector2(x,y),TVector2(0.0,0.0));
2464
2465 //  cnew.SetGhost(nfit); //cnew.SetX(1,sigx); cnew.SetY(1,sigy); cnew.SetZ(1,dist);
2466   //cnew.fClusterType=cnew.PhysicsContribution();
2467   new((*fRawClusters)[fRawClusters->GetLast()+1]) AliMUONCluster(cnew);
2468 //  if (fDebug) cout << fNRawClusters << " " << fChamberId << endl;
2469   //fNPeaks++;
2470 }
2471
2472 //_____________________________________________________________________________
2473 Int_t AliMUONClusterFinderAZ::FindLocalMaxima(TObjArray *pixArray, Int_t *localMax, Double_t *maxVal)
2474 {
2475 /// Find local maxima in pixel space for large preclusters in order to
2476 /// try to split them into smaller pieces (to speed up the MLEM procedure)
2477 /// or to find additional fitting seeds if clusters were not completely resolved
2478
2479   TH2D *hist = NULL;
2480   //if (pixArray == fPixArray) hist = (TH2D*) gROOT->FindObject("anode");
2481   //else { hist = (TH2D*) gROOT->FindObject("anode1"); cout << hist << endl; }
2482   //if (hist) hist->Delete();
2483
2484   Double_t xylim[4] = {999, 999, 999, 999};
2485   Int_t nPix = pixArray->GetEntriesFast();
2486   AliMUONPixel *pixPtr = 0;
2487   for (Int_t ipix=0; ipix<nPix; ipix++) {
2488     pixPtr = (AliMUONPixel*) pixArray->UncheckedAt(ipix);
2489     for (Int_t i=0; i<4; i++)
2490          xylim[i] = TMath::Min (xylim[i], (i%2 ? -1 : 1)*pixPtr->Coord(i/2));
2491   }
2492   for (Int_t i=0; i<4; i++) xylim[i] -= pixPtr->Size(i/2);
2493
2494   Int_t nx = TMath::Nint ((-xylim[1]-xylim[0])/pixPtr->Size(0)/2);
2495   Int_t ny = TMath::Nint ((-xylim[3]-xylim[2])/pixPtr->Size(1)/2);
2496   if (pixArray == fPixArray) hist = new TH2D("anode","anode",nx,xylim[0],-xylim[1],ny,xylim[2],-xylim[3]);
2497   else hist = new TH2D("anode1","anode1",nx,xylim[0],-xylim[1],ny,xylim[2],-xylim[3]);
2498   for (Int_t ipix=0; ipix<nPix; ipix++) {
2499     pixPtr = (AliMUONPixel*) pixArray->UncheckedAt(ipix);
2500     hist->Fill(pixPtr->Coord(0), pixPtr->Coord(1), pixPtr->Charge());
2501   }
2502 //  if (fDraw && pixArray == fPixArray) fDraw->DrawHist("c2", hist);
2503
2504   Int_t nMax = 0, indx;
2505   Int_t *isLocalMax = new Int_t[ny*nx];
2506   for (Int_t i=0; i<ny*nx; i++) isLocalMax[i] = 0;
2507
2508   for (Int_t i=1; i<=ny; i++) {
2509     indx = (i-1) * nx;
2510     for (Int_t j=1; j<=nx; j++) {
2511       if (hist->GetCellContent(j,i) < 0.5) continue;
2512       //if (isLocalMax[indx+j-1] < 0) continue;
2513       if (isLocalMax[indx+j-1] != 0) continue;
2514       FlagLocalMax(hist, i, j, isLocalMax);
2515     }
2516   }
2517
2518   for (Int_t i=1; i<=ny; i++) {
2519     indx = (i-1) * nx;
2520     for (Int_t j=1; j<=nx; j++) {
2521       if (isLocalMax[indx+j-1] > 0) {
2522         localMax[nMax] = indx + j - 1;
2523         maxVal[nMax++] = hist->GetCellContent(j,i);
2524         if (nMax > 99) AliFatal(" Too many local maxima !!!");
2525       }
2526     }
2527   }
2528   if (fDebug) cout << " Local max: " << nMax << endl;
2529   delete [] isLocalMax; isLocalMax = 0;
2530   return nMax;
2531 }
2532
2533 //_____________________________________________________________________________
2534 void AliMUONClusterFinderAZ::FlagLocalMax(TH2D *hist, Int_t i, Int_t j, Int_t *isLocalMax)
2535 {
2536 /// Flag pixels (whether or not local maxima)
2537
2538   Int_t nx = hist->GetNbinsX();
2539   Int_t ny = hist->GetNbinsY();
2540   Int_t cont = TMath::Nint (hist->GetCellContent(j,i));
2541   Int_t cont1 = 0, indx = (i-1)*nx+j-1, indx1 = 0, indx2 = 0;
2542
2543   for (Int_t i1=i-1; i1<i+2; i1++) {
2544     if (i1 < 1 || i1 > ny) continue;
2545     indx1 = (i1 - 1) * nx;
2546     for (Int_t j1=j-1; j1<j+2; j1++) {
2547       if (j1 < 1 || j1 > nx) continue;
2548       if (i == i1 && j == j1) continue;
2549       indx2 = indx1 + j1 - 1;
2550       cont1 = TMath::Nint (hist->GetCellContent(j1,i1));
2551       if (cont < cont1) { isLocalMax[indx] = -1; return; }
2552       else if (cont > cont1) isLocalMax[indx2] = -1;
2553       else { // the same charge
2554         isLocalMax[indx] = 1;
2555         if (isLocalMax[indx2] == 0) {
2556           FlagLocalMax(hist, i1, j1, isLocalMax);
2557           if (isLocalMax[indx2] < 0) { isLocalMax[indx] = -1; return; }
2558           else isLocalMax[indx2] = -1;
2559         }
2560       }
2561     }
2562   }
2563   isLocalMax[indx] = 1; // local maximum
2564 }
2565
2566 //_____________________________________________________________________________
2567 void AliMUONClusterFinderAZ::FindCluster(Int_t *localMax, Int_t iMax)
2568 {
2569 /// Find pixel cluster around local maximum \a iMax and pick up pads
2570 /// overlapping with it
2571
2572   TH2D *hist = (TH2D*) gROOT->FindObject("anode");
2573   Int_t nx = hist->GetNbinsX();
2574   Int_t ny = hist->GetNbinsY();
2575   Int_t ic = localMax[iMax] / nx + 1;
2576   Int_t jc = localMax[iMax] % nx + 1;
2577   Bool_t *used = new Bool_t[ny*nx];
2578   for (Int_t i=0; i<ny*nx; i++) used[i] = kFALSE;
2579
2580   // Drop all pixels from the array - pick up only the ones from the cluster
2581   fPixArray->Delete();
2582
2583   Double_t wx = hist->GetXaxis()->GetBinWidth(1)/2;
2584   Double_t wy = hist->GetYaxis()->GetBinWidth(1)/2;
2585   Double_t yc = hist->GetYaxis()->GetBinCenter(ic);
2586   Double_t xc = hist->GetXaxis()->GetBinCenter(jc);
2587   Double_t cont = hist->GetCellContent(jc,ic);
2588   AliMUONPixel *pixPtr = new AliMUONPixel (xc, yc, wx, wy, cont);
2589   fPixArray->Add((TObject*)pixPtr);
2590   used[(ic-1)*nx+jc-1] = kTRUE;
2591   AddBin(hist, ic, jc, 1, used, (TObjArray*)0); // recursive call
2592
2593   Int_t nPix = fPixArray->GetEntriesFast(), npad = fnPads[0] + fnPads[1];
2594   for (Int_t i=0; i<nPix; i++) {
2595     ((AliMUONPixel*)fPixArray->UncheckedAt(i))->SetSize(0,wx);
2596     ((AliMUONPixel*)fPixArray->UncheckedAt(i))->SetSize(1,wy);
2597   }
2598   if (fDebug) cout << iMax << " " << nPix << endl;
2599
2600   Float_t xy[4], xy12[4];
2601   // Pick up pads which overlap with found pixels
2602   for (Int_t i=0; i<npad; i++) fPadIJ[1][i] = -1;
2603   for (Int_t i=0; i<nPix; i++) {
2604     pixPtr = (AliMUONPixel*) fPixArray->UncheckedAt(i);
2605     for (Int_t j=0; j<4; j++)
2606       xy[j] = pixPtr->Coord(j/2) + (j%2 ? 1 : -1)*pixPtr->Size(j/2);
2607     for (Int_t j=0; j<npad; j++)
2608       if (Overlap(xy, j, xy12, 0)) fPadIJ[1][j] = 0; // flag for use
2609   }
2610
2611   delete [] used; used = 0;
2612 }
2613
2614 //_____________________________________________________________________________
2615 void AliMUONClusterFinderAZ::AddVirtualPad()
2616 {
2617 /// Add virtual pad (with small charge) to improve fit for some
2618 /// clusters (when pad with max charge is at the extreme of the cluster)
2619
2620   // Get number of pads in X and Y-directions
2621   Int_t nInX = -1, nInY;
2622   PadsInXandY(nInX, nInY);
2623   //return;
2624
2625   // Add virtual pad only if number of pads per direction == 2
2626   if (nInX != 2 && nInY != 2) return;
2627
2628   // Find pads with max charge
2629   Int_t maxpad[2][2] = {{-1, -1}, {-1, -1}}, cath;
2630   Double_t sigmax[2] = {0}, aamax[2] = {0};
2631   for (Int_t j=0; j<fnPads[0]+fnPads[1]; j++) {
2632     if (fPadIJ[1][j] != 0) continue;
2633     cath = fPadIJ[0][j];
2634     if (fXyq[2][j] > sigmax[cath]) {
2635       maxpad[cath][1] = maxpad[cath][0];
2636       aamax[cath] = sigmax[cath];
2637       sigmax[cath] = fXyq[2][j];
2638       maxpad[cath][0] = j;
2639     }
2640   }
2641   if (maxpad[0][0] >= 0 && maxpad[0][1] < 0 || maxpad[1][0] >= 0 && maxpad[1][1] < 0) {
2642     for (Int_t j=0; j<fnPads[0]+fnPads[1]; j++) {
2643       if (fPadIJ[1][j] != 0) continue;
2644       cath = fPadIJ[0][j];
2645       if (j == maxpad[cath][0] || j == maxpad[cath][1]) continue;
2646       if (fXyq[2][j] > aamax[cath]) {
2647         aamax[cath] = fXyq[2][j];
2648         maxpad[cath][1] = j;
2649       }
2650     }
2651   }
2652   // Check for mirrors (side X on cathode 0)
2653   Bool_t mirror = kFALSE;
2654   if (maxpad[0][0] >= 0 && maxpad[1][0] >= 0) {
2655     mirror = fXyq[3][maxpad[0][0]] < fXyq[4][maxpad[0][0]];
2656     if (!mirror && TMath::Abs(fXyq[3][maxpad[0][0]]-fXyq[3][maxpad[1][0]]) < 0.001) {
2657       // Special case when pads on both cathodes have the same size
2658       Int_t yud[2] = {0};
2659       for (Int_t j = 0; j < fnPads[0]+fnPads[1]; j++) {
2660         cath = fPadIJ[0][j];
2661         if (j == maxpad[cath][0]) continue;
2662         if (fPadIJ[2][j] != fPadIJ[2][maxpad[cath][0]]) continue;
2663         if (fPadIJ[3][j] + 1 == fPadIJ[3][maxpad[cath][0]] ||
2664             fPadIJ[3][j] - 1 == fPadIJ[3][maxpad[cath][0]]) yud[cath]++;
2665       }
2666       if (!yud[0]) mirror = kTRUE; // take the other cathode
2667     } // if (!mirror &&...
2668   } // if (maxpad[0][0] >= 0 && maxpad[1][0] >= 0)
2669
2670   // Find neughbours of pads with max charges
2671   Int_t xList[10], yList[10], ix0, iy0, ix, iy, neighb;
2672   for (cath=0; cath<2; cath++) {
2673     if (!cath && maxpad[0][0] < 0) continue; // one-sided cluster - cathode 1
2674     if (cath && maxpad[1][0] < 0) break; // one-sided cluster - cathode 0
2675     if (maxpad[1][0] >= 0) {
2676       if (!mirror) {
2677         if (!cath && nInY != 2) continue;
2678         if (cath && nInX != 2 && (maxpad[0][0] >= 0 || nInY != 2)) continue;
2679       } else {
2680         if (!cath && nInX != 2) continue;
2681         if (cath && nInY != 2 && (maxpad[0][0] >= 0 || nInX != 2)) continue;
2682       }
2683     }
2684
2685     Int_t iAddX = 0, iAddY = 0, ix1 = 0, iy1 = 0, iPad = 0;
2686     if (maxpad[0][0] < 0) iPad = 1;
2687
2688     for (iPad=0; iPad<2; iPad++) {
2689       if (maxpad[cath][iPad] < 0) continue;
2690       if (iPad && !iAddX && !iAddY) break;
2691       if (iPad && fXyq[2][maxpad[cath][1]] / sigmax[cath] < 0.5) break;
2692
2693       Int_t neighbx = 0, neighby = 0;
2694       ix0 = fPadIJ[2][maxpad[cath][iPad]];
2695       iy0 = fPadIJ[3][maxpad[cath][iPad]];
2696       TObjArray neighbours;
2697       AliMpPad pad = fSegmentation[cath]->PadByIndices(AliMpIntPair(ix0, iy0));
2698       Int_t nn = fSegmentation[cath]->GetNeighbours(pad,neighbours);
2699       for (Int_t j=0; j<nn; j++) {
2700         AliMpPad* pad = static_cast<AliMpPad*>(neighbours.At(j));
2701         Int_t xx = pad->GetIndices().GetFirst();
2702         Int_t yy = pad->GetIndices().GetSecond();
2703         if (TMath::Abs(xx-ix0) == 1 || xx*ix0 == -1) neighbx++;
2704         if (TMath::Abs(yy-iy0) == 1 || yy*iy0 == -1) neighby++;
2705       }
2706       if (!mirror) {
2707         if (cath) neighb = neighbx;
2708         else neighb = neighby;
2709         if (maxpad[0][0] < 0) neighb += neighby;
2710         else if (maxpad[1][0] < 0) neighb += neighbx;
2711       } else {
2712         if (!cath) neighb = neighbx;
2713         else neighb = neighby;
2714         if (maxpad[0][0] < 0) neighb += neighbx;
2715         else if (maxpad[1][0] < 0) neighb += neighby;
2716       }
2717
2718       for (Int_t j=0; j<fnPads[0]+fnPads[1]; j++) {
2719         if (fPadIJ[0][j] != cath) continue;
2720         ix = fPadIJ[2][j];
2721         iy = fPadIJ[3][j];
2722         if (iy == iy0 && ix == ix0) continue;
2723         for (Int_t k=0; k<nn; k++) {
2724           if (xList[k] != ix || yList[k] != iy) continue;
2725           if (!mirror) {
2726             if ((!cath || maxpad[0][0] < 0) &&
2727                 (TMath::Abs(iy-iy0) == 1 || iy*iy0 == -1)) {
2728               if (!iPad && TMath::Abs(ix-ix0) == 1 || ix*ix0 == -1) ix1 = xList[k]; //19-12-05
2729               xList[k] = yList[k] = 0;
2730               neighb--;
2731               break;
2732             }
2733             if ((cath || maxpad[1][0] < 0) &&
2734                 (TMath::Abs(ix-ix0) == 1 || ix*ix0 == -1)) {
2735               if (!iPad) ix1 = xList[k]; //19-12-05
2736               xList[k] = yList[k] = 0;
2737               neighb--;
2738             }
2739           } else {
2740             if ((!cath || maxpad[0][0] < 0) &&
2741                 (TMath::Abs(ix-ix0) == 1 || ix*ix0 == -1)) {
2742               if (!iPad) ix1 = xList[k]; //19-12-05
2743               xList[k] = yList[k] = 0;
2744               neighb--;
2745               break;
2746             }
2747             if ((cath || maxpad[1][0] < 0) &&
2748                 (TMath::Abs(iy-iy0) == 1 || iy*iy0 == -1)) {
2749               xList[k] = yList[k] = 0;
2750               neighb--;
2751             }
2752           }
2753           break;
2754         } // for (Int_t k=0; k<nn;
2755         if (!neighb) break;
2756       } // for (Int_t j=0; j<fnPads[0]+fnPads[1];
2757       if (!neighb) continue;
2758
2759       // Add virtual pad
2760       Int_t npads, isec;
2761       isec = 0;
2762       for (Int_t j=0; j<nn; j++) {
2763         if (xList[j] == 0 && yList[j] == 0) continue;
2764         npads = fnPads[0] + fnPads[1];
2765         fPadIJ[0][npads] = cath;
2766         fPadIJ[1][npads] = 0;
2767         ix = xList[j];
2768         iy = yList[j];
2769         if (TMath::Abs(ix-ix0) == 1 || ix*ix0 == -1) {
2770           if (iy != iy0) continue; // new segmentation - check
2771           if (nInX != 2) continue; // new
2772           if (!mirror) {
2773             if (!cath && maxpad[1][0] >= 0) continue;
2774           } else {
2775             if (cath && maxpad[0][0] >= 0) continue;
2776           }
2777           if (iPad && !iAddX) continue;
2778     AliMpPad pad = fSegmentation[cath]->PadByIndices(AliMpIntPair(ix,iy));
2779     fXyq[0][npads] = pad.Position().X();
2780     fXyq[1][npads] = pad.Position().Y();
2781           if (fXyq[0][npads] > 1.e+5) continue; // temporary fix
2782           if (ix == ix1) continue; //19-12-05
2783           if (ix1 == ix0) continue;
2784           if (maxpad[1][0] < 0 || mirror && maxpad[0][0] >= 0) {
2785             if (!iPad) fXyq[2][npads] = TMath::Min (sigmax[0]/100, 5.);
2786             else fXyq[2][npads] = TMath::Min (aamax[0]/100, 5.);
2787           }
2788           else {
2789             if (!iPad) fXyq[2][npads] = TMath::Min (sigmax[1]/100, 5.);
2790             else fXyq[2][npads] = TMath::Min (aamax[1]/100, 5.);
2791           }
2792           fXyq[2][npads] = TMath::Max (fXyq[2][npads], (float)1);
2793           fXyq[3][npads] = -2; // flag
2794           fPadIJ[2][npads] = ix;
2795           fPadIJ[3][npads] = iy;
2796           fnPads[1]++;
2797           iAddX = npads;
2798           if (fDebug) printf(" ***** Add virtual pad in X ***** %f %f %f %3d %3d \n", fXyq[2][npads],
2799                              fXyq[0][npads], fXyq[1][npads], ix, iy);
2800           ix1 = ix0;
2801           continue;
2802         }
2803         if (nInY != 2) continue;
2804         if (!mirror && cath && maxpad[0][0] >= 0) continue;
2805         if (mirror && !cath && maxpad[1][0] >= 0) continue;
2806         if (TMath::Abs(iy-iy0) == 1 || TMath::Abs(iy*iy0) == 1) {
2807           if (ix != ix0) continue; // new segmentation - check
2808           if (iPad && !iAddY) continue;
2809     AliMpPad pad = fSegmentation[cath]->PadByIndices(AliMpIntPair(ix,iy));
2810     fXyq[0][npads] = pad.Position().X();
2811     fXyq[1][npads] = pad.Position().Y();
2812           if (iy1 == iy0) continue;
2813           //if (iPad && iy1 == iy0) continue;
2814           if (maxpad[0][0] < 0 || mirror && maxpad[1][0] >= 0) {
2815             if (!iPad) fXyq[2][npads] = TMath::Min (sigmax[1]/15, fgkZeroSuppression);
2816             else fXyq[2][npads] = TMath::Min (aamax[1]/15, fgkZeroSuppression);
2817           }
2818           else {
2819             if (!iPad) fXyq[2][npads] = TMath::Min (sigmax[0]/15, fgkZeroSuppression);
2820             else fXyq[2][npads] = TMath::Min (aamax[0]/15, fgkZeroSuppression);
2821           }
2822           fXyq[2][npads] = TMath::Max (fXyq[2][npads], (float)1);
2823           fXyq[3][npads] = -2; // flag
2824           fPadIJ[2][npads] = ix;
2825           fPadIJ[3][npads] = iy;
2826           fnPads[1]++;
2827           iAddY = npads;
2828           if (fDebug) printf(" ***** Add virtual pad in Y ***** %f %f %f %3d %3d \n", fXyq[2][npads],
2829                              fXyq[0][npads], fXyq[1][npads], ix, iy);
2830           iy1 = iy0;
2831         }
2832       } // for (Int_t j=0; j<nn;
2833     } // for (Int_t iPad=0;
2834   } // for (cath=0; cath<2;
2835   return;
2836 }
2837
2838 //_____________________________________________________________________________
2839 void AliMUONClusterFinderAZ::PadsInXandY(Int_t &nInX, Int_t &nInY)
2840 {
2841 /// Find number of pads in X and Y-directions (excluding virtual ones and
2842 /// overflows)
2843
2844   static Int_t nXsaved = 0, nYsaved = 0;
2845   nXsaved = nYsaved = 0;
2846   //if (nInX >= 0) {nInX = nXsaved; nInY = nYsaved; return; }
2847   Float_t *xPad0 = NULL, *yPad0 = NULL, *xPad1 = NULL, *yPad1 = NULL;
2848   Float_t wMinX[2] = {99, 99}, wMinY[2] = {99, 99};
2849   Int_t *nPad0 = NULL, *nPad1 = NULL;
2850   Int_t nPads = fnPads[0] + fnPads[1];
2851   if (fnPads[0]) {
2852     xPad0 = new Float_t[nPads];
2853     yPad0 = new Float_t[nPads];
2854     nPad0 = new Int_t[nPads];
2855   }
2856   if (fnPads[1]) {
2857     xPad1 = new Float_t[nPads];
2858     yPad1 = new Float_t[nPads];
2859     nPad1 = new Int_t[nPads];
2860   }
2861   Int_t n0 = 0, n1 = 0, cath, npadx[2] = {1, 1}, npady[2] = {1, 1};
2862   for (Int_t j = 0; j < nPads; j++) {
2863     if (nInX < 0 && fPadIJ[1][j] != 0) continue; // before fit
2864     else if (nInX == 0 && fPadIJ[1][j] != 1) continue; // fit - exclude overflows
2865     else if (nInX > 0 && fPadIJ[1][j] != 1 && fPadIJ[1][j] != -9) continue; // exclude non-marked
2866     if (nInX <= 0 && fXyq[2][j] > fgkSaturation-1) continue; // skip overflows
2867     cath = fPadIJ[0][j];
2868     if (fXyq[3][j] > 0) { // exclude virtual pads
2869       wMinX[cath] = TMath::Min (wMinX[cath], fXyq[3][j]);
2870       wMinY[cath] = TMath::Min (wMinY[cath], fXyq[4][j]);
2871       //20-12-05 }
2872       if (cath) { xPad1[n1] = fXyq[0][j]; yPad1[n1++] = fXyq[1][j]; }
2873       else { xPad0[n0] = fXyq[0][j]; yPad0[n0++] = fXyq[1][j]; }
2874     }
2875   }
2876
2877   // Sort
2878   if (n0) {
2879     TMath::Sort (n0, xPad0, nPad0); // in X
2880     for (Int_t i = 1; i < n0; i++)
2881       if (xPad0[nPad0[i]] - xPad0[nPad0[i-1]] < -0.01) npadx[0]++;
2882     TMath::Sort (n0, yPad0, nPad0); // in Y
2883     for (Int_t i = 1; i < n0; i++)
2884       if (yPad0[nPad0[i]] - yPad0[nPad0[i-1]] < -0.01) npady[0]++;
2885   }
2886
2887   if (n1) {
2888     TMath::Sort (n1, xPad1, nPad1); // in X
2889     for (Int_t i = 1; i < n1; i++)
2890       if (xPad1[nPad1[i]] - xPad1[nPad1[i-1]] < -0.01) npadx[1]++;
2891     TMath::Sort (n1, yPad1, nPad1); // in Y
2892     for (Int_t i = 1; i < n1; i++)
2893       if (yPad1[nPad1[i]] - yPad1[nPad1[i-1]] < -0.01) npady[1]++;
2894   }
2895   if (fnPads[0]) { delete [] xPad0; delete [] yPad0; delete [] nPad0; }
2896   if (fnPads[1]) { delete [] xPad1; delete [] yPad1; delete [] nPad1; }
2897   if (TMath::Abs (wMinY[0] - wMinY[1]) < 1.e-3) nInY = TMath::Max (npady[0], npady[1]);
2898   else nInY = wMinY[0] < wMinY[1] ? npady[0] : npady[1];
2899   if (TMath::Abs (wMinX[0] - wMinX[1]) < 1.e-3) nInX = TMath::Max (npadx[0], npadx[1]);
2900   else nInX = wMinX[0] < wMinX[1] ? npadx[0] : npadx[1];
2901 }
2902
2903 //_____________________________________________________________________________
2904 void AliMUONClusterFinderAZ::Simple()
2905 {
2906 /// Process simple cluster (small number of pads) without EM-procedure
2907
2908   Int_t nForFit = 1, clustFit[1] = {0}, nfit;
2909   Double_t parOk[3] = {0.};
2910   TObjArray *clusters[1];
2911   clusters[0] = fPixArray;
2912   for (Int_t i = 0; i < fnPads[0]+fnPads[1]; i++) {
2913     if (fXyq[2][i] > fgkSaturation-1) fPadIJ[1][i] = -9;
2914     else fPadIJ[1][i] = 1;
2915   }
2916   nfit = Fit(1, nForFit, clustFit, clusters, parOk);
2917 }
2918
2919 //_____________________________________________________________________________
2920 void AliMUONClusterFinderAZ::Errors(AliMUONRawCluster* /*clus*/)
2921 {
2922 /// Correct reconstructed coordinates for some clusters and evaluate errors
2923
2924   AliWarning("Reimplement me!");
2925
2926 //  Double_t qTot = clus->GetCharge(0), fmin = clus->GetChi2(0);
2927 //  Double_t xreco = clus->GetX(0), yreco = clus->GetY(0), zreco = clus->GetZ(0);
2928 //  Double_t sigmax[2] = {0};
2929 //
2930 //  Int_t nInX = 1, nInY, maxdig[2] ={-1, -1}, digit, cath1, isec;
2931 //  PadsInXandY(nInX, nInY);
2932 //
2933 //  // Find pad with maximum signal
2934 //  for (Int_t cath = 0; cath < 2; cath++) {
2935 //    for (Int_t j = 0; j < clus->GetMultiplicity(cath); j++) {
2936 //      cath1 = cath;
2937 //      digit = clus->GetIndex(j, cath);
2938 //      if (digit < 0) { cath1 = TMath::Even(cath); digit = -digit - 1; } // from the other cathode
2939 //
2940 //      if (clus->GetContrib(j,cath) > sigmax[cath1]) {
2941 //      sigmax[cath1] = clus->GetContrib(j,cath);
2942 //      maxdig[cath1] = digit;
2943 //      }
2944 //    }
2945 //  }
2946 //
2947 //  // Size of pad with maximum signal and reco coordinate distance from the pad center
2948 //  AliMUONVDigit *mdig = 0;
2949 //  Double_t wx[2], wy[2], dxc[2], dyc[2];
2950 //  Float_t xpad, ypad, zpad;
2951 //  Int_t ix, iy;
2952 //  for (Int_t cath = 0; cath < 2; cath++) {
2953 //    if (maxdig[cath] < 0) continue;
2954 //    mdig = fDigitStore->Find(maxdig[cath]);
2955 //    isec = fSegmentation[cath]->Sector(mdig->PadX(), mdig->PadY());
2956 //    wx[cath] = fSegmentation[cath]->Dpx(isec);
2957 //    wy[cath] = fSegmentation[cath]->Dpy(isec);
2958 //    fSegmentation[cath]->GetPadI(xreco, yreco, zreco, ix, iy);
2959 //    isec = fSegmentation[cath]->Sector(ix, iy);
2960 //    if (isec > 0) {
2961 //      fSegmentation[cath]->GetPadC(ix, iy, xpad, ypad, zpad);
2962 //      dxc[cath] = xreco - xpad;
2963 //      dyc[cath] = yreco - ypad;
2964 //    }
2965 //  }
2966 //
2967 //  // Check if pad with max charge at the edge (number of neughbours)
2968 //  Int_t nn, xList[10], yList[10], neighbx[2][2] = {{0,0}, {0,0}}, neighby[2][2]= {{0,0}, {0,0}};
2969 //  for (Int_t cath = 0; cath < 2; cath++) {
2970 //    if (maxdig[cath] < 0) continue;
2971 //    mdig = fDigitStore->FindObject(maxdig[cath]);
2972 //    fSegmentation[cath]->Neighbours(mdig->PadX(), mdig->PadY(), &nn, xList, yList);
2973 //    isec = fSegmentation[cath]->Sector(mdig->PadX(), mdig->PadY());
2974 //    for (Int_t j=0; j<nn; j++) {
2975 //      fSegmentation[cath]->GetPadC(xList[j], yList[j], xpad, ypad, zpad);
2976 //      if (TMath::Abs(xpad) < 1 && TMath::Abs(ypad) < 1) continue;
2977 //      if (xList[j] == mdig->PadX()-1 || mdig->PadX() == 1 &&
2978 //        xList[j] == -1) neighbx[cath][0] = 1;
2979 //      else if (xList[j] == mdig->PadX()+1 || mdig->PadX() == -1 &&
2980 //             xList[j] == 1) neighbx[cath][1] = 1;
2981 //      if (yList[j] == mdig->PadY()-1 || mdig->PadY() == 1 &&
2982 //        yList[j] == -1) neighby[cath][0] = 1;
2983 //      else if (yList[j] == mdig->PadY()+1 || mdig->PadY() == -1 &&
2984 //             yList[j] == 1) neighby[cath][1] = 1;
2985 //    } // for (Int_t j=0; j<nn;
2986 //    if (neighbx[cath][0] && neighbx[cath][1]) neighbx[cath][0] = 0;
2987 //    else if (neighbx[cath][1]) neighbx[cath][0] = -1;
2988 //    else neighbx[cath][0] = 1;
2989 //    if (neighby[cath][0] && neighby[cath][1]) neighby[cath][0] = 0;
2990 //    else if (neighby[cath][1]) neighby[cath][0] = -1;
2991 //    else neighby[cath][0] = 1;
2992 //  }
2993 //
2994 //  Int_t iOver = clus->GetClusterType();
2995 //  // One-sided cluster
2996 //  if (!clus->GetMultiplicity(0)) {
2997 //    neighby[0][0] = neighby[1][0];
2998 //    wy[0] = wy[1];
2999 //    if (iOver < 99) iOver += 100 * iOver;
3000 //    dyc[0] = dyc[1];
3001 //  } else if (!clus->GetMultiplicity(1)) {
3002 //    neighbx[1][0] = neighbx[0][0];
3003 //    wx[1] = wx[0];
3004 //    if (iOver < 99) iOver += 100 * iOver;
3005 //    dxc[1] = dxc[0];
3006 //  }
3007 //
3008 //  // Apply corrections and evaluate errors
3009 //  Double_t errY, errX;
3010 //  Errors(nInY, nInX, neighby[0][0],neighbx[1][0], fmin, wy[0]*10, wx[1]*10, iOver,
3011 //       dyc[0], dxc[1], qTot, yreco, xreco, errY, errX);
3012 //  errY = TMath::Max (errY, 0.01);
3013 //  //errY = 0.01;
3014 //  //errX = TMath::Max (errX, 0.144);
3015 //  clus->SetX(0, xreco); clus->SetY(0, yreco);
3016 //  clus->SetErrX(errX); clus->SetErrY(errY);
3017 }
3018
3019 //_____________________________________________________________________________
3020 void AliMUONClusterFinderAZ::Errors(Int_t ny, Int_t nx, Int_t iby, Int_t ibx, Double_t fmin,
3021                                     Double_t wy, Double_t wx, Int_t iover,
3022                                     Double_t dyc, Double_t /*dxc*/, Double_t qtot,
3023                                     Double_t &yrec, Double_t &xrec, Double_t &erry, Double_t &errx)
3024 {
3025 /// Correct reconstructed coordinates for some clusters and evaluate errors
3026
3027     erry = 0.01;
3028     errx = 0.144;
3029     Int_t iovery = iover % 100;
3030     Double_t corr = 0;
3031
3032 /* ---> Ny = 1 */
3033     if (ny == 1) {
3034       if (iby != 0) {
3035         // edge effect
3036         yrec += iby * (0.1823+0.2008)/2;
3037         erry = 0.04587;
3038       } else {
3039         // Find "effective pad width"
3040         Double_t width = 0.218 / (1.31e-4 * TMath::Exp (2.688 * TMath::Log(qtot)) + 1) * 2;
3041         width = TMath::Min (width, 0.4);
3042         erry = width / TMath::Sqrt(12.);
3043         erry = TMath::Max (erry, 0.01293);
3044       }
3045       goto x; //return;
3046     }
3047
3048 /* ---> "Bad" fit */
3049     if (fmin > 0.4) {
3050       erry = 0.1556;
3051       if (ny == 5) erry = 0.06481;
3052       goto x; //return;
3053     }
3054
3055 /* ---> By != 0 */
3056     if (iby != 0) {
3057       if (ny > 2) {
3058         erry = 0.00417; //0.01010
3059       } else {
3060         // ny = 2
3061         if (dyc * iby > -0.05) {
3062           Double_t dyc2 = dyc * dyc;
3063           if (iby < 0) {
3064             corr = 0.019 - 0.602 * dyc + 8.739 * dyc2 - 44.209 * dyc2 * dyc;
3065             corr = TMath::Min (corr, TMath::Abs(-0.25-dyc));
3066             yrec -= corr;
3067             //dyc -= corr;
3068             erry = 0.00814;
3069           } else {
3070             corr = 0.006 + 0.300 * dyc + 6.147 * dyc2 + 42.039 * dyc2 * dyc;
3071             corr = TMath::Min (corr, 0.25-dyc);
3072             yrec += corr;
3073             //dyc += corr;
3074             erry = 0.01582;
3075           }
3076         } else {
3077           erry = (0.00303 + 0.00296) / 2;
3078         }
3079       }
3080       goto x; //return;
3081     }
3082
3083 /* ---> Overflows */
3084     if (iovery != 0) {
3085       if (qtot < 3000) {
3086         erry = 0.0671;
3087       } else {
3088         if (iovery > 1) {
3089           erry = 0.09214;
3090         } else if (TMath::Abs(wy - 5) < 0.1) {
3091           erry = 0.061; //0.06622
3092         } else {
3093           erry = 0.00812; // 0.01073
3094         }
3095       }
3096       goto x; //return;
3097     }
3098
3099 /* ---> "Good" but very high signal */
3100     if (qtot > 4000) {
3101       if (TMath::Abs(wy - 4) < 0.1) {
3102         erry = 0.00117;
3103       } else if (fmin < 0.03 && qtot < 6000) {
3104         erry = 0.01003;
3105       } else {
3106         erry = 0.1931;
3107       }
3108       goto x; //return;
3109     }
3110
3111 /* ---> "Good" clusters */
3112     if (ny > 3) {
3113       if (TMath::Abs(wy - 5) < 0.1) {
3114         erry = 0.0011; //0.00304
3115       } else if (qtot < 400.) {
3116         erry = 0.0165;
3117       } else {
3118         erry = 0.00135; // 0.00358
3119       }
3120     } else if (ny == 3) {
3121       if (TMath::Abs(wy - 4) < 0.1) {
3122         erry = 35.407 / (1 + TMath::Exp(5.511*TMath::Log(qtot/265.51))) + 11.564;
3123         //erry = 83.512 / (1 + TMath::Exp(3.344*TMath::Log(qtot/211.58))) + 12.260;
3124       } else {
3125         erry = 147.03 / (1 + TMath::Exp(1.713*TMath::Log(qtot/73.151))) + 9.575;
3126         //erry = 91.743 / (1 + TMath::Exp(2.332*TMath::Log(qtot/151.67))) + 11.453;
3127       }
3128       erry *= 1.e-4;
3129     } else {
3130       // ny = 2
3131       if (TMath::Abs(wy - 4) < 0.1) {
3132         erry = 60.800 / (1 + TMath::Exp(3.305*TMath::Log(qtot/104.53))) + 11.702;
3133         //erry = 73.128 / (1 + TMath::Exp(5.676*TMath::Log(qtot/120.93))) + 17.839;
3134       } else {
3135         erry = 117.98 / (1 + TMath::Exp(2.005*TMath::Log(qtot/37.649))) + 21.431;
3136         //erry = 99.066 / (1 + TMath::Exp(4.900*TMath::Log(qtot/107.57))) + 25.315;
3137       }
3138       erry *= 1.e-4;
3139     }
3140     //return;
3141
3142  x:
3143 /* ---> X-coordinate */
3144 /* ---> Y-side */
3145     if (wx > 11) {
3146       errx = 0.0036;
3147       xrec -= 0.1385;
3148       return;
3149     }
3150 /* ---> Nx = 1 */
3151     if (nx == 1) {
3152       if (TMath::Abs(wx - 6) < 0.1) {
3153         if (qtot < 40) errx = 0.1693;
3154         else errx = 0.06241;
3155       } else if (TMath::Abs(wx - 7.5) < 0.1) {
3156         if (qtot < 40) errx = 0.2173;
3157         else errx = 0.07703;
3158       } else if (TMath::Abs(wx - 10) < 0.1) {
3159         if (ibx == 0) {
3160           if (qtot < 40) errx = 0.2316;
3161           else errx = 0.1426;
3162         } else {
3163           xrec += (0.2115 + 0.1942) / 2 * ibx;
3164           errx = 0.1921;
3165         }
3166       }
3167       return;
3168     }
3169 /* ---> "Bad" fit */
3170     if (fmin > 0.5) {
3171       errx = 0.1591;
3172       return;
3173     }
3174 /* ---> Bx != 0 */
3175     if (ibx != 0) {
3176       if (ibx > 0) { errx = 0.06761; xrec -= 0.03832; }
3177       else { errx = 0.06653; xrec += 0.02581; }
3178       return;
3179     }
3180 /* ---> Overflows */
3181     if (iover != 0) {
3182       if (TMath::Abs(wx - 6) < 0.1) errx = 0.06979;
3183       else if (TMath::Abs(wx - 7.5) < 0.1) errx = 0.1089;
3184       else if (TMath::Abs(wx - 10) < 0.1) errx = 0.09847;
3185       return;
3186     }
3187 /* ---> Good */
3188     if (TMath::Abs(wx - 6) < 0.1) errx = 0.06022;
3189     else if (TMath::Abs(wx - 7.5) < 0.1) errx = 0.07247;
3190     else if (TMath::Abs(wx - 10) < 0.1) errx = 0.07359;
3191 }
3192
3193 //___________________________________________________________________________
3194 void AliMUONClusterFinderAZ::ResetRawClusters()
3195 {
3196   /// Reset tracks information
3197   if (fRawClusters) fRawClusters->Clear("C");
3198 }