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