3e11db9770e499d73f8442afc62c7c2c1ee20865
[u/mrichter/AliRoot.git] / MUON / AliMUONClusterFinderAZ.cxx
1 #include "AliMUONClusterFinderAZ.h"
2
3 #include <fcntl.h>
4 #include <Riostream.h>
5 #include <TROOT.h>
6 #include <TCanvas.h>
7 #include <TLine.h>
8 #include <TTree.h>
9 #include <TH2.h>
10 #include <TView.h>
11 #include <TStyle.h>
12 #include <TMinuit.h>
13 #include <TMatrixD.h>
14
15 #include "AliHeader.h"
16 #include "AliRun.h"
17 #include "AliMUON.h"
18 #include "AliMUONChamber.h"
19 #include "AliMUONDigit.h"
20 #include "AliMUONHit.h"
21 #include "AliMUONChamber.h"
22 #include "AliMUONRawCluster.h"
23 #include "AliMUONClusterInput.h"
24 #include "AliMUONPixel.h"
25
26 // This function is used for fitting
27 void fcn1(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
28
29 ClassImp(AliMUONClusterFinderAZ)
30
31 AliMUONClusterFinderAZ* AliMUONClusterFinderAZ::fgClusterFinder = NULL;
32 TMinuit* AliMUONClusterFinderAZ::fgMinuit = NULL; 
33
34 //_____________________________________________________________________________
35 AliMUONClusterFinderAZ::AliMUONClusterFinderAZ(Bool_t draw=0, Int_t iReco=0)
36 {
37 // Constructor
38   for (Int_t i=0; i<4; i++) {fHist[i] = 0;}
39   fMuonDigits = 0;
40   fSegmentation[1] = fSegmentation[0] = 0; 
41   if (!fgClusterFinder) fgClusterFinder = this;
42   if (!fgMinuit) fgMinuit = new TMinuit(8);
43   fDraw = draw;
44   fReco = iReco;
45   fPixArray = new TObjArray(20); 
46   /*
47   fPoints = 0;
48   fPhits = 0;
49   fRpoints = 0;
50   fCanvas = 0;
51   fNextCathode = kFALSE; 
52   fColPad = 0;
53   */
54 }
55
56 //_____________________________________________________________________________
57 AliMUONClusterFinderAZ::~AliMUONClusterFinderAZ()
58 {
59   // Destructor
60   delete fgMinuit; fgMinuit = 0; delete fPixArray; fPixArray = 0;
61   /*
62   // Delete space point structure
63   if (fPoints) fPoints->Delete();
64   delete fPoints;
65   fPoints     = 0;
66   //
67   if (fPhits) fPhits->Delete();
68   delete fPhits;
69   fPhits     = 0;
70   //
71   if (fRpoints) fRpoints->Delete();
72   delete fRpoints;
73   fRpoints     = 0;
74   */
75 }
76
77 //_____________________________________________________________________________
78 void AliMUONClusterFinderAZ::FindRawClusters()
79 {
80 // To provide the same interface as in AliMUONClusterFinderVS
81
82   EventLoop (gAlice->GetHeader()->GetEvent(), AliMUONClusterInput::Instance()->Chamber());
83 }
84
85 //_____________________________________________________________________________
86 void AliMUONClusterFinderAZ::EventLoop(Int_t nev=0, Int_t ch=0)
87 {
88 // Loop over events
89   
90   FILE *lun = 0;
91   TCanvas *c1 = 0;
92   TView *view = 0;
93   TH2F *hist = 0;
94   Double_t p1[3]={0}, p2[3];
95   TTree *TR = 0;
96   if (fDraw) {
97     // File
98     lun = fopen("pool.dat","w");
99     c1 = new TCanvas("c1","Clusters",0,0,600,700);
100     c1->Divide(1,2);
101     new TCanvas("c2","Mlem",700,0,600,350);
102   }
103
104 newev:
105   Int_t nparticles = 0, nent;
106   if (!fReco) nparticles = gAlice->GetEvent(nev);
107   else nparticles = gAlice->GetNtrack();
108   cout << "nev         " << nev <<endl;
109   cout << "nparticles  " << nparticles <<endl;
110   if (nparticles <= 0) return;
111   
112   TTree *TH = gAlice->TreeH();
113   Int_t ntracks = (Int_t) TH->GetEntries();
114   cout<<"ntracks "<<ntracks<<endl;
115     
116   // Get pointers to Alice detectors and Digits containers
117   AliMUON *MUON  = (AliMUON*) gAlice->GetModule("MUON");
118   if (!MUON) return;
119   //       TClonesArray *Particles = gAlice->Particles();
120   if (!fReco) {
121     TR = gAlice->TreeR();
122     if (TR) {
123       MUON->ResetRawClusters();
124       nent = (Int_t) TR->GetEntries();
125       if (nent != 1) {
126         cout << "Error in MUONdrawClust" << endl;
127         cout << " nent = " <<  nent << " not equal to 1" << endl;
128         //exit(0);
129       }
130     } // if (TR)
131   } // if (!fReco)
132
133   TTree *TD = gAlice->TreeD();
134   //MUON->ResetDigits();
135
136   TClonesArray *MUONrawclust;
137   AliMUONChamber*       iChamber = 0;
138   
139   // As default draw the first cluster of the chamber #0
140       
141 newchamber:
142   if (ch > 9) {if (fReco) return; nev++; ch = 0; goto newev;}
143   //gAlice->ResetDigits();
144   fMuonDigits  = MUON->DigitsAddress(ch);
145   if (fMuonDigits == 0) return;
146   iChamber = &(MUON->Chamber(ch));
147   fSegmentation[0] = iChamber->SegmentationModel(1);
148   fSegmentation[1] = iChamber->SegmentationModel(2);
149   fResponse = iChamber->ResponseModel();
150     
151   nent = 0;
152  
153   if (TD) {
154     nent = (Int_t) TD->GetEntries();
155     //printf(" entries %d \n", nent);
156   }
157
158   Int_t ndigits[2]={9,9}, nShown[2]={0};
159   for (Int_t i=0; i<2; i++) {
160     for (Int_t j=0; j<kDim; j++) {fUsed[i][j]=kFALSE;}
161   }
162
163 next:
164   if (ndigits[0] == nShown[0] && ndigits[1] == nShown[1]) {
165     // No more clusters
166     if (fReco) return;
167     ch++;
168     goto newchamber; // next chamber
169   }
170   Float_t xpad, ypad, zpad, zpad0;
171   TLine *line[99]={0};
172   Int_t nLine = 0;
173   Bool_t first = kTRUE;
174   cout << " *** Event # " << nev << " chamber: " << ch << endl;
175   fnPads[0] = fnPads[1] = 0;
176   for (Int_t i=0; i<kDim; i++) {fPadIJ[1][i] = 0;}
177   //for (Int_t iii = 0; iii<999; iii++) { 
178   for (Int_t iii = 0; iii<2; iii++) { 
179     Int_t cath = TMath::Odd(iii);
180     gAlice->ResetDigits();
181     TD->GetEvent(cath);
182     fMuonDigits  = MUON->DigitsAddress(ch);
183
184     ndigits[cath] = fMuonDigits->GetEntriesFast();
185     if (!ndigits[0] && !ndigits[1]) {if (fReco) return; ch++; goto newchamber;}
186     if (ndigits[cath] == 0) continue;
187     cout << " ndigits: " << ndigits[cath] << " " << cath << endl;
188
189     AliMUONDigit  *mdig;
190     Int_t digit;
191
192     Bool_t EOC = kTRUE; // end-of-cluster
193     for (digit = 0; digit < ndigits[cath]; digit++) {
194       mdig    = (AliMUONDigit*)fMuonDigits->UncheckedAt(digit);
195       if (mdig->Cathode() != cath) continue;
196       if (first) {
197         // Find first unused pad
198         if (fUsed[cath][digit]) continue;
199         fSegmentation[cath]->GetPadC(mdig->PadX(),mdig->PadY(),xpad,ypad,zpad0);
200       } else {
201         if (fUsed[cath][digit]) continue;
202         fSegmentation[cath]->GetPadC(mdig->PadX(),mdig->PadY(),xpad,ypad,zpad);
203         if (TMath::Abs(zpad-zpad0)>0.1) continue; // different slats
204         // Find a pad overlapping with the cluster
205         if (!Overlap(cath,mdig)) continue;
206       }
207       // Add pad - recursive call
208       AddPad(cath,digit);
209       EOC = kFALSE;
210       if (digit >= 0) break;
211     }
212     if (first && EOC) {
213       // No more unused pads 
214       if (cath == 0) continue; // on cathode #0 - check #1
215       else {
216         // No more clusters
217         if (fReco) return;
218         ch++;
219         goto newchamber; // next chamber
220       }
221     }
222     if (EOC) break; // cluster found
223     first = kFALSE;
224     cout << " nPads: " << fnPads[cath] << " " << nShown[cath]+fnPads[cath] << " " << cath << endl;
225   } // for (Int_t iii = 0;
226
227   
228   if (fReco) goto skip;
229   char hName[4];
230   for (Int_t cath = 0; cath<2; cath++) {
231     // Build histograms
232     if (fHist[cath*2]) {fHist[cath*2]->Delete(); fHist[cath*2] = 0;}
233     if (fHist[cath*2+1]) {fHist[cath*2+1]->Delete(); fHist[cath*2+1] = 0;}
234     if (fnPads[cath] == 0) continue; // cluster on one cathode only
235     Float_t wxMin=999, wxMax=0, wyMin=999, wyMax=0; 
236     Int_t minDx=0, maxDx=0, minDy=0, maxDy=0;
237     for (Int_t i=0; i<fnPads[0]+fnPads[1]; i++) {
238       if (fPadIJ[0][i] != cath) continue;
239       if (fXyq[3][i] < wxMin) {wxMin = fXyq[3][i]; minDx = i;}
240       if (fXyq[3][i] > wxMax) {wxMax = fXyq[3][i]; maxDx = i;}
241       if (fXyq[4][i] < wyMin) {wyMin = fXyq[4][i]; minDy = i;}
242       if (fXyq[4][i] > wyMax) {wyMax = fXyq[4][i]; maxDy = i;}
243     }
244     cout << minDx << maxDx << minDy << maxDy << endl;
245     Int_t nx, ny, padSize;
246     Float_t xmin=9999, xmax=-9999, ymin=9999, ymax=-9999;
247     if (TMath::Nint(fXyq[3][minDx]*1000) == TMath::Nint(fXyq[3][maxDx]*1000) &&
248         TMath::Nint(fXyq[4][minDy]*1000) == TMath::Nint(fXyq[4][maxDy]*1000)) {
249       // the same segmentation
250       cout << " Same" << endl;
251       cout << fXyq[3][minDx] << " " << fXyq[3][maxDx] << " " << fXyq[4][minDy] << " " << fXyq[4][maxDy] << endl;
252       for (Int_t i=0; i<fnPads[0]+fnPads[1]; i++) {
253         if (fPadIJ[0][i] != cath) continue;
254         if (fXyq[0][i] < xmin) xmin = fXyq[0][i];
255         if (fXyq[0][i] > xmax) xmax = fXyq[0][i];
256         if (fXyq[1][i] < ymin) ymin = fXyq[1][i];
257         if (fXyq[1][i] > ymax) ymax = fXyq[1][i];
258       }
259       xmin -= fXyq[3][minDx]; xmax += fXyq[3][minDx];
260       ymin -= fXyq[4][minDy]; ymax += fXyq[4][minDy];
261       nx = TMath::Nint ((xmax-xmin)/wxMin/2);
262       ny = TMath::Nint ((ymax-ymin)/wyMin/2);
263       sprintf(hName,"h%d",cath*2);
264       fHist[cath*2] = new TH2F(hName,"cluster",nx,xmin,xmax,ny,ymin,ymax);
265       cout << fHist[cath*2] << " " << fnPads[cath] << endl;
266       for (Int_t i=0; i<fnPads[0]+fnPads[1]; i++) {
267         if (fPadIJ[0][i] != cath) continue;
268         fHist[cath*2]->Fill(fXyq[0][i],fXyq[1][i],fXyq[2][i]);
269         //cout << fXyq[0][i] << fXyq[1][i] << fXyq[2][i] << endl;
270       }
271     } else {
272       // different segmentation in the cluster
273       cout << " Different" << endl;
274       cout << fXyq[3][minDx] << " " << fXyq[3][maxDx] << " " << fXyq[4][minDy] << " " << fXyq[4][maxDy] << endl;
275       Int_t nOK = 0;
276       Int_t indx, locMin, locMax;
277       if (TMath::Nint(fXyq[3][minDx]*1000) != TMath::Nint(fXyq[3][maxDx]*1000)) {
278         // different segmentation along x
279         indx = 0;
280         locMin = minDx;
281         locMax = maxDx;
282       } else {
283         // different segmentation along y
284         indx = 1;
285         locMin = minDy;
286         locMax = maxDy;
287       }
288       Int_t loc = locMin;
289       for (Int_t i=0; i<2; i++) {
290         // loop over different pad sizes
291         if (i>0) loc = locMax;
292         padSize = TMath::Nint(fXyq[indx+3][loc]*1000);
293         xmin = 9999; xmax = -9999; ymin = 9999; ymax = -9999;
294         for (Int_t j=0; j<fnPads[0]+fnPads[1]; j++) {
295           if (fPadIJ[0][j] != cath) continue;
296           if (TMath::Nint(fXyq[indx+3][j]*1000) != padSize) continue;
297           nOK++;
298           xmin = TMath::Min (xmin,fXyq[0][j]);
299           xmax = TMath::Max (xmax,fXyq[0][j]);
300           ymin = TMath::Min (ymin,fXyq[1][j]);
301           ymax = TMath::Max (ymax,fXyq[1][j]);
302         }
303         xmin -= fXyq[3][loc]; xmax += fXyq[3][loc];
304         ymin -= fXyq[4][loc]; ymax += fXyq[4][loc];
305         nx = TMath::Nint ((xmax-xmin)/fXyq[3][loc]/2);
306         ny = TMath::Nint ((ymax-ymin)/fXyq[4][loc]/2);
307         sprintf(hName,"h%d",cath*2+i);
308         fHist[cath*2+i] = new TH2F(hName,"cluster",nx,xmin,xmax,ny,ymin,ymax);
309         for (Int_t j=0; j<fnPads[0]+fnPads[1]; j++) {
310           if (fPadIJ[0][j] != cath) continue;
311           if (TMath::Nint(fXyq[indx+3][j]*1000) != padSize) continue;
312           fHist[cath*2+i]->Fill(fXyq[0][j],fXyq[1][j],fXyq[2][j]);
313         }
314       } // for (Int_t i=0;
315       if (nOK != fnPads[cath]) cout << " *** Too many segmentations: nPads, nOK " << fnPads[cath] << " " << nOK << endl;
316     } // if (TMath::Nint(fXyq[3][minDx]*1000)
317   } // for (Int_t cath = 0;
318         
319   // Draw histograms and coordinates
320   for (Int_t cath=0; cath<2; cath++) {
321     if (cath == 0) ModifyHistos();
322     if (fnPads[cath] == 0) continue; // cluster on one cathode only
323     if (fDraw) {
324       c1->cd(cath+1);
325       gPad->SetTheta(55);
326       gPad->SetPhi(30);
327       Double_t x, y, x0, y0, r1=999, r2=0;
328       if (fHist[cath*2+1]) {
329         // 
330         x0 = fHist[cath*2]->GetXaxis()->GetXmin() - 1000*TMath::Cos(30*TMath::Pi()/180);
331         y0 = fHist[cath*2]->GetYaxis()->GetXmin() - 1000*TMath::Sin(30*TMath::Pi()/180);
332         r1 = 0;
333         Int_t ihist=cath*2;
334         for (Int_t iy=1; iy<=fHist[ihist]->GetNbinsY(); iy++) {
335           y = fHist[ihist]->GetYaxis()->GetBinCenter(iy) 
336             + fHist[ihist]->GetYaxis()->GetBinWidth(iy);
337           for (Int_t ix=1; ix<=fHist[ihist]->GetNbinsX(); ix++) {
338             if (fHist[ihist]->GetCellContent(ix,iy) > 0.1) {
339               x = fHist[ihist]->GetXaxis()->GetBinCenter(ix)
340                 + fHist[ihist]->GetXaxis()->GetBinWidth(ix);
341               r1 = TMath::Max (r1,TMath::Sqrt((x-x0)*(x-x0)+(y-y0)*(y-y0)));
342             }
343           }
344         }
345         ihist = cath*2 + 1 ;
346         for (Int_t iy=1; iy<=fHist[ihist]->GetNbinsY(); iy++) {
347           y = fHist[ihist]->GetYaxis()->GetBinCenter(iy)
348             + fHist[ihist]->GetYaxis()->GetBinWidth(iy);
349           for (Int_t ix=1; ix<=fHist[ihist]->GetNbinsX(); ix++) {
350             if (fHist[ihist]->GetCellContent(ix,iy) > 0.1) {
351               x = fHist[ihist]->GetXaxis()->GetBinCenter(ix)
352                 + fHist[ihist]->GetXaxis()->GetBinWidth(ix);
353               r2 = TMath::Max (r2,TMath::Sqrt((x-x0)*(x-x0)+(y-y0)*(y-y0)));
354             }
355           }
356         }
357         cout << r1 << " " << r2 << endl;
358       } // if (fHist[cath*2+1])
359       if (r1 > r2) {
360         //fHist[cath*2]->Draw("lego1");
361         fHist[cath*2]->Draw("lego1Fb");
362         //if (fHist[cath*2+1]) fHist[cath*2+1]->Draw("lego1SameAxisBb");
363         if (fHist[cath*2+1]) fHist[cath*2+1]->Draw("lego1SameAxisBbFb");
364       } else {
365         //fHist[cath*2+1]->Draw("lego1");
366         fHist[cath*2+1]->Draw("lego1Fb");
367         //fHist[cath*2]->Draw("lego1SameAxisBb");
368         fHist[cath*2]->Draw("lego1SameAxisFbBb");
369       }
370       c1->Update();
371     } // if (fDraw)
372   } // for (Int_t cath = 0;
373
374   // Draw generated hits
375   Double_t xNDC[6];
376   hist = fHist[0] ? fHist[0] : fHist[2];
377   p2[2] = hist->GetMaximum();
378   view = 0;
379   if (c1) view = c1->Pad()->GetView();
380   cout << " *** GEANT hits *** " << endl;
381   fnMu = 0;
382   Int_t ix, iy, iok;
383   for (Int_t i=0; i<ntracks; i++) {
384     TH->GetEvent(i);
385     for (AliMUONHit* mHit=(AliMUONHit*)MUON->FirstHit(-1); 
386          mHit;
387          mHit=(AliMUONHit*)MUON->NextHit()) {
388       if (mHit->Chamber() != ch+1) continue;  // chamber number
389       if (TMath::Abs(mHit->Z()-zpad0) > 1) continue; // different slat
390       p2[0] = p1[0] = mHit->X();        // x-pos of hit
391       p2[1] = p1[1] = mHit->Y();        // y-pos
392       if (p1[0] < hist->GetXaxis()->GetXmin() || 
393           p1[0] > hist->GetXaxis()->GetXmax()) continue;
394       if (p1[1] < hist->GetYaxis()->GetXmin() || 
395           p1[1] > hist->GetYaxis()->GetXmax()) continue;
396       // Check if track comes thru pads with signal
397       iok = 0;
398       for (Int_t ihist=0; ihist<4; ihist++) {
399         if (!fHist[ihist]) continue;
400         ix = fHist[ihist]->GetXaxis()->FindBin(p1[0]);
401         iy = fHist[ihist]->GetYaxis()->FindBin(p1[1]);
402         if (fHist[ihist]->GetCellContent(ix,iy) > 0.5) {iok = 1; break;}
403       }
404       if (!iok) continue;
405       gStyle->SetLineColor(1);
406       if (TMath::Abs((Int_t)mHit->Particle()) == 13) {
407         gStyle->SetLineColor(4);
408         fnMu++;
409         if (fnMu <= 2) {
410           fxyMu[fnMu-1][0] = p1[0];
411           fxyMu[fnMu-1][1] = p1[1];
412         }
413       }     
414       printf(" X=%10.4f, Y=%10.4f, Z=%10.4f\n",p1[0],p1[1],mHit->Z());
415       if (view) {
416         view->WCtoNDC(p1, &xNDC[0]);
417         view->WCtoNDC(p2, &xNDC[3]);
418         for (Int_t ipad=1; ipad<3; ipad++) {
419           c1->cd(ipad);
420           //c1->DrawLine(xpad[0],xpad[1],xpad[3],xpad[4]);
421           line[nLine] = new TLine(xNDC[0],xNDC[1],xNDC[3],xNDC[4]);
422           line[nLine++]->Draw();
423         }
424       }
425     } // for (AliMUONHit* mHit=
426   } // for (Int_t i=0; i<ntracks;
427
428   // Draw reconstructed coordinates
429   MUONrawclust = MUON->RawClustAddress(ch);
430   TR->GetEvent(ch);
431   //cout << MUONrawclust << " " << MUONrawclust->GetEntries() << endl;
432   AliMUONRawCluster *mRaw;
433   gStyle->SetLineColor(3);
434   cout << " *** Reconstructed hits *** " << endl;
435   for (Int_t i=0; i<MUONrawclust->GetEntries(); i++) {
436     mRaw = (AliMUONRawCluster*)MUONrawclust->UncheckedAt(i);
437     if (TMath::Abs(mRaw->fZ[0]-zpad0) > 1) continue; // different slat
438     p2[0] = p1[0] = mRaw->fX[0];        // x-pos of hit
439     p2[1] = p1[1] = mRaw->fY[0];        // y-pos
440     if (p1[0] < hist->GetXaxis()->GetXmin() || 
441         p1[0] > hist->GetXaxis()->GetXmax()) continue;
442     if (p1[1] < hist->GetYaxis()->GetXmin() || 
443         p1[1] > hist->GetYaxis()->GetXmax()) continue;
444     /*
445       TD->GetEvent(cath);
446       cout << mRaw->fMultiplicity[0] << mRaw->fMultiplicity[1] << endl;
447       for (Int_t j=0; j<mRaw->fMultiplicity[cath]; j++) {
448       Int_t digit = mRaw->fIndexMap[j][cath];
449       cout << ((AliMUONDigit*)fMuonDigits->UncheckedAt(digit))->Signal() << endl;
450       }
451     */
452     // Check if track comes thru pads with signal
453     iok = 0;
454     for (Int_t ihist=0; ihist<4; ihist++) {
455       if (!fHist[ihist]) continue;
456       ix = fHist[ihist]->GetXaxis()->FindBin(p1[0]);
457       iy = fHist[ihist]->GetYaxis()->FindBin(p1[1]);
458       if (fHist[ihist]->GetCellContent(ix,iy) > 0.5) {iok = 1; break;}
459     }
460     if (!iok) continue;
461     printf(" X=%10.4f, Y=%10.4f, Z=%10.4f\n",p1[0],p1[1],mRaw->fZ[0]);
462     if (view) {
463       view->WCtoNDC(p1, &xNDC[0]);
464       view->WCtoNDC(p2, &xNDC[3]);
465       for (Int_t ipad=1; ipad<3; ipad++) {
466         c1->cd(ipad);
467         line[nLine] = new TLine(xNDC[0],xNDC[1],xNDC[3],xNDC[4]);
468         line[nLine++]->Draw();
469       }
470     }
471   } // for (Int_t i=0; i<MUONrawclust->GetEntries();
472   if (fDraw) c1->Update();
473
474 skip:
475   // Use MLEM for cluster finder
476   fZpad = zpad0;
477   Int_t nMax = 1, localMax[100], maxPos[100];
478   Double_t maxVal[100];
479   
480   if (CheckPrecluster(nShown)) {
481     BuildPixArray();
482     if (fnPads[0]+fnPads[1] > 50) nMax = FindLocalMaxima(localMax, maxVal);
483     if (nMax > 1) TMath::Sort(nMax, maxVal, maxPos, kTRUE); // in decreasing order
484     for (Int_t i=0; i<nMax; i++) {
485       if (nMax > 1) FindCluster(localMax, maxPos[i]);
486       if (!MainLoop()) cout << " MainLoop failed " << endl;
487       if (i < nMax-1) {
488         for (Int_t j=0; j<fnPads[0]+fnPads[1]; j++) {
489           if (fPadIJ[1][j] == 0) continue; // pad charge was not modified
490           fPadIJ[1][j] = 0;
491           fXyq[2][j] = fXyq[5][j]; // use backup charge value
492         }
493       }
494     }
495   }
496   if (fReco) goto next;
497
498   for (Int_t i=0; i<fnMu; i++) {
499     // Check again if muon come thru the used pads (due to extra splitting)
500     for (Int_t j=0; j<fnPads[0]+fnPads[1]; j++) {
501       if (TMath::Abs(fxyMu[i][0]-fXyq[0][j])<fXyq[3][j] && 
502           TMath::Abs(fxyMu[i][1]-fXyq[1][j])<fXyq[4][j]) {
503         printf("%12.3e %12.3e %12.3e %12.3e\n",fxyMu[i][2],fxyMu[i][3],fxyMu[i][4],fxyMu[i][5]);
504         if (lun) fprintf(lun,"%4d %2d %12.3e %12.3e %12.3e %12.3e\n",nev,ch,fxyMu[i][2],fxyMu[i][3],fxyMu[i][4],fxyMu[i][5]);
505         break;
506       }
507     }
508   } // for (Int_t i=0; i<fnMu;
509
510   // What's next?
511   char command[8];
512   cout << " What is next? " << endl;
513   command[0] = ' '; 
514   if (fDraw) gets(command);
515   if (command[0] == 'n' || command[0] == 'N') {nev++; goto newev;} // next event 
516   else if (command[0] == 'q' || command[0] == 'Q') {fclose(lun); return;} // exit display 
517   //else if (command[0] == 'r' || command[0] == 'R') goto redraw; // redraw points
518   else if (command[0] == 'c' || command[0] == 'C') {
519     // new chamber
520     sscanf(command+1,"%d",&ch);
521     goto newchamber;
522   } 
523   else if (command[0] == 'e' || command[0] == 'E') {
524     // new event
525     sscanf(command+1,"%d",&nev);
526     goto newev;
527   } 
528   else goto next; // Next cluster
529 }
530
531 //_____________________________________________________________________________
532 void AliMUONClusterFinderAZ::ModifyHistos(void)
533 {
534   // Modify histograms to bring them to the same size
535   Int_t nhist = 0;
536   Float_t hlim[4][4], hbin[4][4]; // first index - xmin, xmax, ymin, ymax
537   Float_t binMin[4] = {999,999,999,999};
538
539   for (Int_t i=0; i<4; i++) {
540     if (!fHist[i]) continue;
541     hlim[0][nhist] = fHist[i]->GetXaxis()->GetXmin(); // xmin
542     hlim[1][nhist] = fHist[i]->GetXaxis()->GetXmax(); // xmax
543     hlim[2][nhist] = fHist[i]->GetYaxis()->GetXmin(); // ymin
544     hlim[3][nhist] = fHist[i]->GetYaxis()->GetXmax(); // ymax
545     hbin[0][nhist] = hbin[1][nhist] = fHist[i]->GetXaxis()->GetBinWidth(1);
546     hbin[2][nhist] = hbin[3][nhist] = fHist[i]->GetYaxis()->GetBinWidth(1);
547     binMin[0] = TMath::Min(binMin[0],hbin[0][nhist]);
548     binMin[2] = TMath::Min(binMin[2],hbin[2][nhist]);
549     nhist++;
550   }
551   binMin[1] = binMin[0];
552   binMin[3] = binMin[2];
553   cout << " Nhist: " << nhist << endl;
554
555   Int_t imin, imax;
556   for (Int_t lim=0; lim<4; lim++) {
557     while (1) {
558       imin = TMath::LocMin(nhist,hlim[lim]);
559       imax = TMath::LocMax(nhist,hlim[lim]);
560       if (TMath::Abs(hlim[lim][imin]-hlim[lim][imax])<0.01*binMin[lim]) break;
561       if (lim == 0 || lim == 2) {
562         // find lower limit
563         hlim[lim][imax] -= hbin[lim][imax];
564       } else {
565         // find upper limit
566         hlim[lim][imin] += hbin[lim][imin];
567       }
568     } // while (1)
569   }
570     
571   // Rebuild histograms 
572   nhist = 0;
573   TH2F *hist = 0;
574   Int_t nx, ny;
575   Double_t x, y, cont, cmax=0;
576   char hName[4];
577   for (Int_t ihist=0; ihist<4; ihist++) {
578     if (!fHist[ihist]) continue;
579     nx = TMath::Nint((hlim[1][nhist]-hlim[0][nhist])/hbin[0][nhist]);
580     ny = TMath::Nint((hlim[3][nhist]-hlim[2][nhist])/hbin[2][nhist]);
581     //hist =  new TH2F("h","hist",nx,hlim[0][nhist],hlim[1][nhist],ny,hlim[2][nhist],hlim[3][nhist]);
582     sprintf(hName,"hh%d",ihist);
583     hist =  new TH2F(hName,"hist",nx,hlim[0][nhist],hlim[1][nhist],ny,hlim[2][nhist],hlim[3][nhist]);
584     for (Int_t i=1; i<=fHist[ihist]->GetNbinsX(); i++) {
585       x = fHist[ihist]->GetXaxis()->GetBinCenter(i);
586       for (Int_t j=1; j<=fHist[ihist]->GetNbinsY(); j++) {
587         y = fHist[ihist]->GetYaxis()->GetBinCenter(j);
588         cont = fHist[ihist]->GetCellContent(i,j);
589         hist->Fill(x,y,cont);
590       }
591     }
592     cmax = TMath::Max (cmax,hist->GetMaximum());
593     fHist[ihist]->Delete();
594     fHist[ihist] = new TH2F(*hist);
595     hist->Delete(); 
596     nhist++;
597   }
598   printf("%f \n",cmax);
599
600   for (Int_t ihist=0; ihist<4; ihist++) {
601     if (!fHist[ihist]) continue;
602     fHist[ihist]->SetMaximum(cmax);
603   }
604 }
605
606 //_____________________________________________________________________________
607 void AliMUONClusterFinderAZ::AddPad(Int_t cath, Int_t digit)
608 {
609   // Add pad to the cluster
610   AliMUONDigit *mdig = (AliMUONDigit*)fMuonDigits->UncheckedAt(digit);
611
612   Int_t charge = mdig->Signal();
613   // get the center of the pad
614   Float_t xpad, ypad, zpad;
615   fSegmentation[cath]->GetPadC(mdig->PadX(), mdig->PadY(), xpad, ypad, zpad);
616   
617   Int_t   isec = fSegmentation[cath]->Sector(mdig->PadX(), mdig->PadY());
618   Int_t nPads = fnPads[0] + fnPads[1];
619   fXyq[0][nPads] = xpad;
620   fXyq[1][nPads] = ypad;
621   fXyq[2][nPads] = charge;
622   fXyq[3][nPads] = fSegmentation[cath]->Dpx(isec)/2;
623   fXyq[4][nPads] = fSegmentation[cath]->Dpy(isec)/2;
624   fXyq[5][nPads] = digit;
625   fPadIJ[0][nPads] = cath;
626   fPadIJ[1][nPads] = 0;
627   fUsed[cath][digit] = kTRUE;
628   //cout << " bbb " << fXyq[cath][2][nPads] << " " << fXyq[cath][0][nPads] << " " << fXyq[cath][1][nPads] << " " << fXyq[cath][3][nPads] << " " << fXyq[cath][4][nPads] << " " << zpad << " " << nPads << endl;
629   fnPads[cath]++;
630
631   // Check neighbours
632   Int_t nn, ix, iy, xList[10], yList[10];
633   AliMUONDigit  *mdig1;
634
635   Int_t ndigits = fMuonDigits->GetEntriesFast();
636   fSegmentation[cath]->Neighbours(mdig->PadX(),mdig->PadY(),&nn,xList,yList); 
637   for (Int_t in=0; in<nn; in++) {
638     ix=xList[in];
639     iy=yList[in];
640     for (Int_t digit1 = 0; digit1 < ndigits; digit1++) {
641       if (digit1 == digit) continue;
642       mdig1 = (AliMUONDigit*)fMuonDigits->UncheckedAt(digit1);
643       if (mdig1->Cathode() != cath) continue;
644       if (!fUsed[cath][digit1] && mdig1->PadX() == ix && mdig1->PadY() == iy) {
645         fUsed[cath][digit1] = kTRUE;
646         // Add pad - recursive call
647         AddPad(cath,digit1);
648       }
649     } //for (Int_t digit1 = 0;
650   } // for (Int_t in=0;
651 }
652
653 //_____________________________________________________________________________
654 Bool_t AliMUONClusterFinderAZ::Overlap(Int_t cath, TObject *dig)
655 {
656   // Check if the pad from one cathode overlaps with a pad 
657   // in the precluster on the other cathode
658
659   AliMUONDigit *mdig = (AliMUONDigit*) dig;
660
661   Float_t xpad, ypad, zpad;
662   fSegmentation[cath]->GetPadC(mdig->PadX(), mdig->PadY(), xpad, ypad, zpad);
663   Int_t   isec = fSegmentation[cath]->Sector(mdig->PadX(), mdig->PadY());
664
665   Float_t xy1[4], xy12[4];
666   xy1[0] = xpad - fSegmentation[cath]->Dpx(isec)/2;
667   xy1[1] = xy1[0] + fSegmentation[cath]->Dpx(isec);
668   xy1[2] = ypad - fSegmentation[cath]->Dpy(isec)/2;
669   xy1[3] = xy1[2] + fSegmentation[cath]->Dpy(isec);
670   //cout << " ok " << fnPads[0]+fnPads[1] << xy1[0] << xy1[1] << xy1[2] << xy1[3] << endl;
671
672   Int_t cath1 = TMath::Even(cath);
673   for (Int_t i=0; i<fnPads[0]+fnPads[1]; i++) {
674     if (fPadIJ[0][i] != cath1) continue;
675     if (Overlap(xy1, i, xy12, 0)) return kTRUE;
676   }
677   return kFALSE;
678 }
679
680 //_____________________________________________________________________________
681 Bool_t AliMUONClusterFinderAZ::Overlap(Float_t *xy1, Int_t iPad, Float_t *xy12, Int_t iSkip)
682 {
683   // Check if the pads xy1 and iPad overlap and return overlap area
684
685   Float_t xy2[4];
686   xy2[0] = fXyq[0][iPad] - fXyq[3][iPad];
687   xy2[1] = fXyq[0][iPad] + fXyq[3][iPad];
688   if (xy1[0] > xy2[1]-1.e-4 || xy1[1] < xy2[0]+1.e-4) return kFALSE;
689   xy2[2] = fXyq[1][iPad] - fXyq[4][iPad];
690   xy2[3] = fXyq[1][iPad] + fXyq[4][iPad];
691   if (xy1[2] > xy2[3]-1.e-4 || xy1[3] < xy2[2]+1.e-4) return kFALSE;
692   if (!iSkip) return kTRUE; // just check overlap (w/out computing the area)
693   xy12[0] = TMath::Max (xy1[0],xy2[0]);
694   xy12[1] = TMath::Min (xy1[1],xy2[1]);
695   xy12[2] = TMath::Max (xy1[2],xy2[2]);
696   xy12[3] = TMath::Min (xy1[3],xy2[3]);
697   return kTRUE;
698 }
699
700 //_____________________________________________________________________________
701 /*
702 Bool_t AliMUONClusterFinderAZ::Overlap(Int_t i, Int_t j, Float_t *xy12, Int_t iSkip)
703 {
704   // Check if the pads i and j overlap and return overlap area
705
706   Float_t xy1[4], xy2[4];
707   return Overlap(xy1, xy2, xy12, iSkip);
708 }
709 */
710 //_____________________________________________________________________________
711 Bool_t AliMUONClusterFinderAZ::CheckPrecluster(Int_t *nShown)
712 {
713   // Check precluster in order to attempt to simplify it (mostly for
714   // two-cathode preclusters)
715
716   Int_t i1, i2;
717   Float_t xy1[4], xy12[4];
718   
719   Int_t npad = fnPads[0] + fnPads[1];
720
721   // If pads have the same size take average of pads on both cathodes 
722   Int_t sameSize = (fnPads[0] && fnPads[1]) ? 1 : 0;
723   if (sameSize) {
724     Double_t xSize = -1, ySize = 0;
725     for (Int_t i=0; i<npad; i++) {
726       if (fXyq[2][i] < 0) continue;
727       if (xSize < 0) { xSize = fXyq[3][i]; ySize = fXyq[4][i]; }
728       if (TMath::Abs(xSize-fXyq[3][i]) > 1.e-4 ||  TMath::Abs(ySize-fXyq[4][i]) > 1.e-4) { sameSize = 0; break; }
729     }
730   } // if (sameSize)
731   if (sameSize && (fnPads[0] > 2 || fnPads[1] > 2)) {
732     nShown[0] += fnPads[0];
733     nShown[1] += fnPads[1];
734     fnPads[0] = fnPads[1] = 0;
735     Int_t div;
736     for (Int_t i=0; i<npad; i++) {
737       if (fXyq[2][i] < 0) continue; // used pad
738       fXyq[2][fnPads[0]] = fXyq[2][i];
739       div = 1;
740       for (Int_t j=i+1; j<npad; j++) {
741         if (fPadIJ[0][j] == fPadIJ[0][i]) continue; // same cathode
742         if (TMath::Abs(fXyq[0][j]-fXyq[0][i]) > 1.e-4) continue;
743         if (TMath::Abs(fXyq[1][j]-fXyq[1][i]) > 1.e-4) continue;
744         fXyq[2][fnPads[0]] += fXyq[2][j];
745         div = 2;
746         fXyq[2][j] = -2;
747         break;
748       }
749       fXyq[2][fnPads[0]] /= div;
750       fXyq[0][fnPads[0]] = fXyq[0][i];
751       fXyq[1][fnPads[0]] = fXyq[1][i];
752       fPadIJ[0][fnPads[0]++] = 0;
753     }
754   } // if (sameSize)
755
756   // Check if one-cathode precluster
757   i1 = fnPads[0]!=0 ? 0 : 1;
758   i2 = fnPads[1]!=0 ? 1 : 0;
759
760   if (i1 != i2) { // two-cathode 
761
762     Int_t *flags = new Int_t[npad];
763     for (Int_t i=0; i<npad; i++) { flags[i] = 0; }
764
765     // Check pad overlaps
766     for (Int_t i=0; i<npad; i++) {
767       if (fPadIJ[0][i] != i1) continue;
768       xy1[0] = fXyq[0][i] - fXyq[3][i];
769       xy1[1] = fXyq[0][i] + fXyq[3][i];
770       xy1[2] = fXyq[1][i] - fXyq[4][i];
771       xy1[3] = fXyq[1][i] + fXyq[4][i];
772       for (Int_t j=0; j<npad; j++) {
773         if (fPadIJ[0][j] != i2) continue;
774         if (!Overlap(xy1, j, xy12, 0)) continue;
775         flags[i] = flags[j] = 1; // mark overlapped pads
776       } // for (Int_t j=0;
777     } // for (Int_t i=0;
778
779     // Check if all pads overlap
780     Int_t digit=0, cath, nFlags=0;
781     for (Int_t i=0; i<npad; i++) {nFlags += !flags[i];}
782     if (nFlags) cout << " nFlags = " << nFlags << endl;
783     //if (nFlags > 2 || (Float_t)nFlags / npad > 0.2) { // why 2 ??? - empirical choice
784     if (nFlags > 0) {
785       for (Int_t i=0; i<npad; i++) {
786         if (flags[i]) continue;
787         digit = TMath::Nint (fXyq[5][i]);
788         cath = fPadIJ[0][i];
789         fUsed[cath][digit] = kFALSE; // release pad
790         fXyq[2][i] = -2;
791         fnPads[cath]--;
792       }
793     } // if (nFlags > 2)
794
795     // Check correlations of cathode charges
796     if (fnPads[0] && fnPads[1]) { // two-cathode
797       Double_t sum[2]={0};
798       Int_t over[2] = {1, 1};
799       for (Int_t i=0; i<npad; i++) {
800         cath = fPadIJ[0][i];
801         if (fXyq[2][i] > 0) sum[cath] += fXyq[2][i];
802         if (fXyq[2][i] > fResponse->MaxAdc()-1) over[cath] = 0;
803       }
804       cout << " Total charge: " << sum[0] << " " << sum[1] << endl;
805       if ((over[0] || over[1]) && TMath::Abs(sum[0]-sum[1])/(sum[0]+sum[1])*2 > 1) { // 3 times difference
806         cout << " Release " << endl;
807         // Big difference
808         cath = sum[0]>sum[1] ? 0 : 1;
809         Int_t imax = 0;
810         Double_t cmax=-1;
811         Double_t *dist = new Double_t[npad];
812         for (Int_t i=0; i<npad; i++) {
813           if (fPadIJ[0][i] != cath) continue;
814           if (fXyq[2][i] < cmax) continue;
815           cmax = fXyq[2][i];
816           imax = i;
817         }
818         // Arrange pads according to their distance to the max, 
819         // normalized to the pad size
820         for (Int_t i=0; i<npad; i++) {
821           dist[i] = 0;
822           if (fPadIJ[0][i] != cath) continue;
823           if (i == imax) continue; 
824           if (fXyq[2][i] < 0) continue;
825           dist[i] = (fXyq[0][i]-fXyq[0][imax])*(fXyq[0][i]-fXyq[0][imax])/
826                      fXyq[3][imax]/fXyq[3][imax]/4;
827           dist[i] += (fXyq[1][i]-fXyq[1][imax])*(fXyq[1][i]-fXyq[1][imax])/
828                       fXyq[4][imax]/fXyq[4][imax]/4;
829           dist[i] = TMath::Sqrt (dist[i]);
830         }
831         TMath::Sort(npad, dist, flags, kFALSE); // in increasing order
832         Int_t indx;
833         Double_t xmax = -1;
834         for (Int_t i=0; i<npad; i++) {
835           indx = flags[i];
836           if (fPadIJ[0][indx] != cath) continue;
837           if (fXyq[2][indx] < 0) continue;
838           if (fXyq[2][indx] <= cmax || TMath::Abs(dist[indx]-xmax)<1.e-3) {
839             // Release pads
840             if (TMath::Abs(dist[indx]-xmax)<1.e-3) 
841                 cmax = TMath::Max((Double_t)(fXyq[2][indx]),cmax);
842             else cmax = fXyq[2][indx];
843             xmax = dist[indx];
844             digit = TMath::Nint (fXyq[5][indx]);
845             fUsed[cath][digit] = kFALSE; 
846             fXyq[2][indx] = -2;
847             fnPads[cath]--;
848             // xmax = dist[i]; // Bug?
849           }
850           else break;
851         } 
852         delete [] dist; dist = 0;
853       } // TMath::Abs(sum[0]-sum[1])...
854     } // if (fnPads[0] && fnPads[1])
855     delete [] flags; flags = 0;
856   } // if (i1 != i2) 
857
858   if (!sameSize) { nShown[0] += fnPads[0]; nShown[1] += fnPads[1]; }
859
860   // Move released pads to the right
861   Int_t beg = 0, end = npad-1, padij;
862   Double_t xyq;
863   while (beg < end) {
864     if (fXyq[2][beg] > 0) { beg++; continue; }
865     for (Int_t j=end; j>beg; j--) {
866       if (fXyq[2][j] < 0) continue;
867       end = j - 1;
868       for (Int_t j1=0; j1<2; j1++) {
869         padij = fPadIJ[j1][beg]; 
870         fPadIJ[j1][beg] = fPadIJ[j1][j];
871         fPadIJ[j1][j] = padij;
872       }
873       for (Int_t j1=0; j1<6; j1++) {
874         xyq = fXyq[j1][beg]; 
875         fXyq[j1][beg] = fXyq[j1][j];
876         fXyq[j1][j] = xyq;
877       }
878       break;
879     } // for (Int_t j=end;
880     beg++;
881   } // while
882   npad = fnPads[0] + fnPads[1];
883   if (npad > 500) { cout << " ***** Too large cluster. Give up. " << npad << endl; return kFALSE; }
884   // Back up charge value
885   for (Int_t j=0; j<npad; j++) fXyq[5][j] = fXyq[2][j];
886
887   return kTRUE;
888 }
889
890 //_____________________________________________________________________________
891 void AliMUONClusterFinderAZ::BuildPixArray()
892 {
893   // Build pixel array for MLEM method
894   
895   Int_t nPix=0, i1, i2;
896   Float_t xy1[4], xy12[4];
897   AliMUONPixel *pixPtr=0;
898
899   Int_t npad = fnPads[0] + fnPads[1];
900
901   // One cathode is empty
902   i1 = fnPads[0]!=0 ? 0 : 1;
903   i2 = fnPads[1]!=0 ? 1 : 0;
904
905   // Build array of pixels on anode plane
906   if (i1 == i2) { // one-cathode precluster
907     for (Int_t j=0; j<npad; j++) {
908       pixPtr = new AliMUONPixel();
909       for (Int_t i=0; i<2; i++) {
910         pixPtr->SetCoord(i, fXyq[i][j]); // pixel coordinates
911         pixPtr->SetSize(i, fXyq[i+3][j]); // pixel size
912       }
913       pixPtr->SetCharge(fXyq[2][j]); // charge
914       fPixArray->Add((TObject*)pixPtr);
915       nPix++;
916     }
917   } else { // two-cathode precluster    
918     for (Int_t i=0; i<npad; i++) {
919       if (fPadIJ[0][i] != i1) continue;
920       xy1[0] = fXyq[0][i] - fXyq[3][i];
921       xy1[1] = fXyq[0][i] + fXyq[3][i];
922       xy1[2] = fXyq[1][i] - fXyq[4][i];
923       xy1[3] = fXyq[1][i] + fXyq[4][i];
924       for (Int_t j=0; j<npad; j++) {
925         if (fPadIJ[0][j] != i2) continue;
926         if (!Overlap(xy1, j, xy12, 1)) continue;
927         pixPtr = new AliMUONPixel();
928         for (Int_t k=0; k<2; k++) {
929           pixPtr->SetCoord(k, (xy12[2*k]+xy12[2*k+1])/2); // pixel coordinates
930           pixPtr->SetSize(k, xy12[2*k+1]-pixPtr->Coord(k)); // size
931         }
932         pixPtr->SetCharge(TMath::Min (fXyq[2][i],fXyq[2][j])); //charge
933         fPixArray->Add((TObject*)pixPtr);
934         nPix++;
935       } // for (Int_t j=0;
936     } // for (Int_t i=0;
937   } // else
938
939   Float_t wxmin=999, wymin=999;
940   for (Int_t i=0; i<npad; i++) {
941     if (fPadIJ[0][i] == i1) wymin = TMath::Min (wymin,fXyq[4][i]);
942     if (fPadIJ[0][i] == i2) wxmin = TMath::Min (wxmin,fXyq[3][i]);
943   }
944   cout << wxmin << " " << wymin << endl;
945
946   // Check if small pixel X-size
947   AjustPixel(wxmin, 0);
948   // Check if small pixel Y-size
949   AjustPixel(wymin, 1);
950   // Check if large pixel size
951   AjustPixel(wxmin, wymin);
952
953   // Remove discarded pixels
954   for (Int_t i=0; i<nPix; i++) {
955     pixPtr = (AliMUONPixel*) fPixArray->UncheckedAt(i);
956     //pixPtr->Print();
957     if (pixPtr->Charge() < 1) { fPixArray->RemoveAt(i); delete pixPtr; }// discarded pixel
958   }
959   fPixArray->Compress();
960   nPix = fPixArray->GetEntriesFast();
961
962   if (nPix > npad) {
963     cout << nPix << endl;
964     // Too many pixels - sort and remove pixels with the lowest signal
965     fPixArray->Sort();
966     for (Int_t i=npad; i<nPix; i++) {
967       pixPtr = (AliMUONPixel*) fPixArray->UncheckedAt(i);
968       //pixPtr->Print();
969       fPixArray->RemoveAt(i);
970       delete pixPtr;
971     }
972     nPix = npad;
973   } // if (nPix > npad)
974
975   // Set pixel charges to the same value (for MLEM)
976   for (Int_t i=0; i<nPix; i++) {
977     pixPtr = (AliMUONPixel*) fPixArray->UncheckedAt(i);
978     //pixPtr->SetCharge(10);
979     cout << i+1 << " " << pixPtr->Coord(0) << " " << pixPtr->Coord(1) << " " << pixPtr->Size(0) << " " << pixPtr->Size(1) << endl;
980   }
981 }
982
983 //_____________________________________________________________________________
984 void AliMUONClusterFinderAZ::AjustPixel(Float_t width, Int_t ixy)
985 {
986   // Check if some pixels have small size (ajust if necessary)
987
988   AliMUONPixel *pixPtr, *pixPtr1 = 0;
989   Int_t ixy1 = TMath::Even(ixy);
990   Int_t nPix = fPixArray->GetEntriesFast();
991
992   for (Int_t i=0; i<nPix; i++) {
993     pixPtr = (AliMUONPixel*) fPixArray->UncheckedAt(i);
994     if (pixPtr->Charge() < 1) continue; // discarded pixel
995     if (pixPtr->Size(ixy)-width < -1.e-4) {
996       // try to merge 
997       cout << " Small X or Y: " << ixy << " " << pixPtr->Size(ixy) << " " << width << " " << pixPtr->Coord(0) << " " << pixPtr->Coord(1) << endl;
998       for (Int_t j=i+1; j<nPix; j++) {
999         pixPtr1 = (AliMUONPixel*) fPixArray->UncheckedAt(j);
1000         if (pixPtr1->Charge() < 1) continue; // discarded pixel
1001         if (TMath::Abs(pixPtr1->Size(ixy)-width) < 1.e-4) continue; // right size 
1002         if (TMath::Abs(pixPtr1->Coord(ixy1)-pixPtr->Coord(ixy1)) > 1.e-4) continue; // different rows/columns
1003         if (TMath::Abs(pixPtr1->Coord(ixy)-pixPtr->Coord(ixy)) < 2*width) {
1004           // merge
1005           pixPtr->SetSize(ixy, width);
1006           pixPtr->SetCoord(ixy, (pixPtr->Coord(ixy)+pixPtr1->Coord(ixy))/2);
1007           pixPtr->SetCharge(TMath::Min (pixPtr->Charge(),pixPtr1->Charge()));
1008           pixPtr1->SetCharge(0);
1009           pixPtr1 = 0;
1010           break;
1011         }
1012       } // for (Int_t j=i+1;
1013       //if (!pixPtr1) { cout << " I am here!" << endl; pixPtr->SetSize(ixy, width); } // ???
1014       //else if (pixPtr1->Charge() > 0.5 || i == nPix-1) {
1015       if (pixPtr1 || i == nPix-1) {
1016         // edge pixel - just increase its size
1017         cout << " Edge ..." << endl; 
1018         for (Int_t j=0; j<fnPads[0]+fnPads[1]; j++) {
1019           // ???if (fPadIJ[0][j] != i1) continue;
1020           if (TMath::Abs(pixPtr->Coord(ixy1)-fXyq[ixy1][j]) > 1.e-4) continue;
1021           if (pixPtr->Coord(ixy) < fXyq[ixy][j]) 
1022             pixPtr->Shift(ixy, -pixPtr->Size(ixy));
1023           else pixPtr->Shift(ixy, pixPtr->Size(ixy));
1024           pixPtr->SetSize(ixy, width);
1025           break;
1026         }
1027       }
1028     } // if (pixPtr->Size(ixy)-width < -1.e-4)
1029   } // for (Int_t i=0; i<nPix;
1030   return;
1031 }
1032   
1033 //_____________________________________________________________________________
1034 void AliMUONClusterFinderAZ::AjustPixel(Float_t wxmin, Float_t wymin)
1035 {
1036   // Check if some pixels have large size (ajust if necessary)
1037
1038   Int_t nx, ny;
1039   Int_t nPix = fPixArray->GetEntriesFast();
1040   AliMUONPixel *pixPtr, *pixPtr1, pix;
1041
1042   // Check if large pixel size
1043   for (Int_t i=0; i<nPix; i++) {
1044     pixPtr = (AliMUONPixel*) fPixArray->UncheckedAt(i);
1045     if (pixPtr->Charge() < 1) continue; // discarded pixel
1046     if (pixPtr->Size(0)-wxmin > 1.e-4 || pixPtr->Size(1)-wymin > 1.e-4) {
1047       cout << " Different " << pixPtr->Size(0) << " " << wxmin << " " << pixPtr->Size(1) << " " << wymin << endl;
1048       pix = *pixPtr;
1049       nx = TMath::Nint (pix.Size(0)/wxmin);
1050       ny = TMath::Nint (pix.Size(1)/wymin);
1051       pix.Shift(0, -pix.Size(0)-wxmin);
1052       pix.Shift(1, -pix.Size(1)-wymin);
1053       pix.SetSize(0, wxmin);
1054       pix.SetSize(1, wymin);
1055       for (Int_t ii=0; ii<nx; ii++) {
1056         pix.Shift(0, wxmin*2);
1057         for (Int_t jj=0; jj<ny; jj++) {
1058           pix.Shift(1, wymin*2);
1059           pixPtr1 = new AliMUONPixel(pix);
1060           fPixArray->Add((TObject*)pixPtr1);
1061         }
1062       }
1063       pixPtr->SetCharge(0);
1064     }
1065   } // for (Int_t i=0; i<nPix;
1066   return;
1067 }
1068
1069 //_____________________________________________________________________________
1070 Bool_t AliMUONClusterFinderAZ::MainLoop()
1071 {
1072   // Repeat MLEM algorithm until pixel size becomes sufficiently small
1073   
1074   TH2D *mlem;
1075
1076   Int_t ix, iy;
1077   //Int_t nn, xList[10], yList[10];
1078   Int_t nPix = fPixArray->GetEntriesFast();
1079   Int_t npadTot = fnPads[0] + fnPads[1], npadOK = 0;
1080   AliMUONPixel *pixPtr = 0;
1081   Double_t *coef = 0, *probi = 0; 
1082   for (Int_t i=0; i<npadTot; i++) if (fPadIJ[1][i] == 0) npadOK++;
1083
1084   while (1) {
1085
1086     mlem = (TH2D*) gROOT->FindObject("mlem");
1087     if (mlem) mlem->Delete();
1088     // Calculate coefficients
1089     cout << " nPix, npadTot, npadOK " << nPix << " " << npadTot << " " << npadOK << endl;
1090
1091     // Calculate coefficients and pixel visibilities
1092     coef = new Double_t [npadTot*nPix];
1093     probi = new Double_t [nPix];
1094     Int_t indx = 0, cath;
1095     for (Int_t ipix=0; ipix<nPix; ipix++) {
1096       pixPtr = (AliMUONPixel*) fPixArray->UncheckedAt(ipix);
1097       probi[ipix] = 0;
1098       for (Int_t j=0; j<npadTot; j++) {
1099         if (fPadIJ[1][j] < 0) { coef[j*nPix+ipix] = 0; continue; }
1100         cath = fPadIJ[0][j];
1101         fSegmentation[cath]->GetPadI(fXyq[0][j],fXyq[1][j],fZpad,ix,iy);
1102         fSegmentation[cath]->SetPad(ix,iy);
1103         /*
1104           fSegmentation[cath]->Neighbours(ix,iy,&nn,xList,yList); 
1105           if (nn != 4) {
1106           cout << nn << ": ";
1107           for (Int_t i=0; i<nn; i++) {cout << xList[i] << " " << yList[i] << ", ";}
1108           cout << endl;
1109           }
1110         */
1111         Double_t sum = 0;
1112         fSegmentation[cath]->SetHit(pixPtr->Coord(0),pixPtr->Coord(1),fZpad);
1113         sum += fResponse->IntXY(fSegmentation[cath]);
1114         indx = j*nPix + ipix;
1115         coef[indx] = sum; 
1116         probi[ipix] += coef[indx];
1117         //cout << j << " " << ipix << " " << coef[indx] << endl;
1118       } // for (Int_t j=0;
1119       //cout << " prob: " << probi[ipix] << endl;
1120       if (probi[ipix] < 0.01) pixPtr->SetCharge(0); // "invisible" pixel
1121     } // for (Int_t ipix=0;
1122
1123     // MLEM algorithm
1124     Mlem(coef, probi);
1125
1126     Double_t xylim[4] = {999, 999, 999, 999};
1127     for (Int_t ipix=0; ipix<nPix; ipix++) {
1128       pixPtr = (AliMUONPixel*) fPixArray->UncheckedAt(ipix);
1129       for (Int_t i=0; i<4; i++) 
1130         xylim[i] = TMath::Min (xylim[i], (i%2 ? -1 : 1)*pixPtr->Coord(i/2));
1131       //cout << ipix+1; pixPtr->Print();
1132     }
1133     for (Int_t i=0; i<4; i++) {
1134       xylim[i] -= pixPtr->Size(i/2); cout << (i%2 ? -1 : 1)*xylim[i] << " "; }
1135     cout << endl;
1136
1137     // Ajust histogram to approximately the same limits as for the pads
1138     // (for good presentation)
1139     //*
1140     Float_t xypads[4];
1141     if (fHist[0]) {
1142       xypads[0] = fHist[0]->GetXaxis()->GetXmin();
1143       xypads[1] = -fHist[0]->GetXaxis()->GetXmax();
1144       xypads[2] = fHist[0]->GetYaxis()->GetXmin();
1145       xypads[3] = -fHist[0]->GetYaxis()->GetXmax();
1146       for (Int_t i=0; i<4; i++) {
1147         while(1) {
1148           if (xylim[i] < xypads[i]) break;
1149           xylim[i] -= 2*pixPtr->Size(i/2);
1150         }
1151       }
1152     } // if (fHist[0])
1153     //*/
1154
1155     Int_t nx = TMath::Nint ((-xylim[1]-xylim[0])/pixPtr->Size(0)/2);
1156     Int_t ny = TMath::Nint ((-xylim[3]-xylim[2])/pixPtr->Size(1)/2);
1157     mlem = new TH2D("mlem","mlem",nx,xylim[0],-xylim[1],ny,xylim[2],-xylim[3]);
1158     for (Int_t ipix=0; ipix<nPix; ipix++) {
1159       pixPtr = (AliMUONPixel*) fPixArray->UncheckedAt(ipix);
1160       mlem->Fill(pixPtr->Coord(0),pixPtr->Coord(1),pixPtr->Charge());
1161     }
1162     //gPad->GetCanvas()->cd(3);
1163     if (fDraw) {
1164       ((TCanvas*)gROOT->FindObject("c2"))->cd();
1165       gPad->SetTheta(55);
1166       gPad->SetPhi(30);
1167       mlem->Draw("lego1Fb");
1168       gPad->Update();
1169       gets((char*)&ix);
1170     }
1171
1172     // Check if the total charge of pixels is too low
1173     Double_t qTot = 0;
1174     for (Int_t i=0; i<nPix; i++) {
1175       pixPtr = (AliMUONPixel*) fPixArray->UncheckedAt(i);
1176       qTot += pixPtr->Charge();
1177     }
1178     if (qTot < 1.e-4 || npadOK < 3 && qTot < 50) {
1179       delete [] coef; delete [] probi; coef = 0; probi = 0;
1180       fPixArray->Delete(); 
1181       return kFALSE; 
1182     }
1183
1184     // Plot data - expectation
1185     /*
1186     Double_t x, y, cont;
1187     for (Int_t j=0; j<npadTot; j++) {
1188       Double_t sum1 = 0;
1189       for (Int_t i=0; i<nPix; i++) {
1190         // Caculate expectation
1191         pixPtr = (AliMUONPixel*) fPixArray->UncheckedAt(i);
1192         sum1 += pixPtr->Charge()*coef[j*nPix+i];
1193       }
1194       sum1 = TMath::Min (sum1,(Double_t)fResponse->MaxAdc());
1195       x = fXyq[0][j];
1196       y = fXyq[1][j];
1197       cath = fPadIJ[0][j];
1198       Int_t ihist = cath*2;
1199       ix = fHist[ihist]->GetXaxis()->FindBin(x);
1200       iy = fHist[ihist]->GetYaxis()->FindBin(y);
1201       cont = fHist[ihist]->GetCellContent(ix,iy);
1202       if (cont == 0 && fHist[ihist+1]) {
1203         ihist += 1;
1204         ix = fHist[ihist]->GetXaxis()->FindBin(x);
1205         iy = fHist[ihist]->GetYaxis()->FindBin(y);
1206       }
1207       fHist[ihist]->SetBinContent(ix,iy,fXyq[2][j]-sum1);
1208     }
1209     ((TCanvas*)gROOT->FindObject("c1"))->cd(1);
1210     //gPad->SetTheta(55);
1211     //gPad->SetPhi(30);
1212     //mlem->Draw("lego1");
1213     gPad->Modified();
1214     ((TCanvas*)gROOT->FindObject("c1"))->cd(2);
1215     gPad->Modified();
1216     */
1217
1218     // Calculate position of the center-of-gravity around the maximum pixel
1219     Double_t xyCOG[2];
1220     FindCOG(mlem, xyCOG);
1221
1222     if (TMath::Min(pixPtr->Size(0),pixPtr->Size(1)) < 0.07 && pixPtr->Size(0) > pixPtr->Size(1)) break;
1223     //if (TMath::Min(pixPtr->Size(0),pixPtr->Size(1)) >= 0.07 || pixPtr->Size(0) < pixPtr->Size(1)) {
1224     // Sort pixels according to the charge
1225     fPixArray->Sort();
1226     /*
1227     for (Int_t i=0; i<nPix; i++) {
1228       pixPtr = (AliMUONPixel*) fPixArray->UncheckedAt(i);
1229       cout << i+1; pixPtr->Print();
1230     }
1231     */
1232     Double_t pixMin = 0.01*((AliMUONPixel*)fPixArray->UncheckedAt(0))->Charge();
1233     pixMin = TMath::Min (pixMin,50.);
1234
1235     // Decrease pixel size and shift pixels to make them centered at 
1236     // the maximum one
1237     indx = (pixPtr->Size(0)>pixPtr->Size(1)) ? 0 : 1;
1238     Double_t width = 0, shift[2]={0};
1239     ix = 1;
1240     for (Int_t i=0; i<4; i++) xylim[i] = 999;
1241     Int_t nPix1 = nPix; nPix = 0;
1242     for (Int_t ipix=0; ipix<nPix1; ipix++) {
1243       pixPtr = (AliMUONPixel*) fPixArray->UncheckedAt(ipix);
1244       if (nPix >= npadOK) { // too many pixels already
1245         fPixArray->RemoveAt(ipix); 
1246         delete pixPtr; 
1247         continue;
1248       }
1249       if (pixPtr->Charge() < pixMin) { // low charge
1250         fPixArray->RemoveAt(ipix); 
1251         delete pixPtr; 
1252         continue;
1253       }
1254       for (Int_t i=0; i<2; i++) {
1255         if (!i) {
1256           pixPtr->SetCharge(10);
1257           pixPtr->SetSize(indx, pixPtr->Size(indx)/2);
1258           width = -pixPtr->Size(indx);
1259           pixPtr->Shift(indx, width);
1260           // Shift pixel position
1261           if (ix) {
1262             ix = 0;
1263             for (Int_t j=0; j<2; j++) {
1264               shift[j] = pixPtr->Coord(j) - xyCOG[j];
1265               shift[j] -= ((Int_t)(shift[j]/pixPtr->Size(j)/2))*pixPtr->Size(j)*2;
1266             }
1267             //cout << ipix << " " << i << " " << shift[0] << " " << shift[1] << endl;
1268           } // if (ix)
1269           pixPtr->Shift(0, -shift[0]);
1270           pixPtr->Shift(1, -shift[1]);
1271         } else {
1272           pixPtr = new AliMUONPixel(*pixPtr);
1273           pixPtr->Shift(indx, -2*width);
1274           fPixArray->Add((TObject*)pixPtr);
1275         } // else
1276         //pixPtr->Print();
1277         for (Int_t i=0; i<4; i++) 
1278           xylim[i] = TMath::Min (xylim[i], (i%2 ? -1 : 1)*pixPtr->Coord(i/2));
1279       } // for (Int_t i=0; i<2;
1280       nPix += 2;
1281     } // for (Int_t ipix=0;
1282
1283     fPixArray->Compress();
1284     nPix = fPixArray->GetEntriesFast();
1285
1286     // Remove excessive pixels
1287     if (nPix > npadOK) {
1288       for (Int_t ipix=npadOK; ipix<nPix; ipix++) { 
1289         pixPtr = (AliMUONPixel*) fPixArray->UncheckedAt(ipix);
1290         fPixArray->RemoveAt(ipix); 
1291         delete pixPtr;
1292       }
1293     } else {
1294       pixPtr = (AliMUONPixel*) fPixArray->UncheckedAt(0);
1295       // add pixels if the maximum is at the limit of pixel area
1296       // start from Y-direction
1297       Int_t j = 0;
1298       for (Int_t i=3; i>-1; i--) {
1299         if (nPix < npadOK && 
1300             TMath::Abs((i%2 ? -1 : 1)*xylim[i]-xyCOG[i/2]) < pixPtr->Size(i/2)) {
1301           pixPtr = new AliMUONPixel(*pixPtr);
1302           pixPtr->SetCoord(i/2, xyCOG[i/2]+(i%2 ? 2:-2)*pixPtr->Size(i/2));
1303           j = TMath::Even (i/2);
1304           pixPtr->SetCoord(j, xyCOG[j]);
1305           fPixArray->Add((TObject*)pixPtr);
1306           nPix++;
1307         }
1308       }
1309     } // else    
1310
1311     fPixArray->Compress();
1312     nPix = fPixArray->GetEntriesFast();
1313     delete [] coef; delete [] probi; coef = 0; probi = 0;
1314   } // while (1)
1315
1316   // remove pixels with low signal or low visibility
1317   // Cuts are empirical !!!
1318   Double_t thresh = TMath::Max (mlem->GetMaximum()/100.,1.);
1319   thresh = TMath::Min (thresh,50.);
1320   Double_t cmax = -1, charge = 0;
1321   for (Int_t i=0; i<nPix; i++) cmax = TMath::Max (cmax,probi[i]); 
1322   // Mark pixels which should be removed
1323   for (Int_t i=0; i<nPix; i++) {
1324     pixPtr = (AliMUONPixel*) fPixArray->UncheckedAt(i);
1325     charge = pixPtr->Charge();
1326     if (charge < thresh) pixPtr->SetCharge(-charge);
1327     else if (cmax > 1.91) {
1328       if (probi[i] < 1.9) pixPtr->SetCharge(-charge);
1329     }
1330     else if (probi[i] < cmax*0.9) pixPtr->SetCharge(-charge);
1331   }
1332   // Move charge of removed pixels to their nearest neighbour (to keep total charge the same)
1333   Int_t near = 0;
1334   for (Int_t i=0; i<nPix; i++) {
1335     pixPtr = (AliMUONPixel*) fPixArray->UncheckedAt(i);
1336     charge = pixPtr->Charge();
1337     if (charge > 0) continue;
1338     near = FindNearest(pixPtr);
1339     pixPtr->SetCharge(0);
1340     pixPtr = (AliMUONPixel*) fPixArray->UncheckedAt(near);
1341     pixPtr->SetCharge(pixPtr->Charge() - charge);
1342   }
1343   // Update histogram
1344   for (Int_t i=0; i<nPix; i++) {
1345     pixPtr = (AliMUONPixel*) fPixArray->UncheckedAt(i);
1346     ix = mlem->GetXaxis()->FindBin(pixPtr->Coord(0));
1347     iy = mlem->GetYaxis()->FindBin(pixPtr->Coord(1));
1348     mlem->SetBinContent(ix, iy, pixPtr->Charge());
1349   }
1350   if (fDraw) {
1351     ((TCanvas*)gROOT->FindObject("c2"))->cd();
1352     gPad->SetTheta(55);
1353     gPad->SetPhi(30);
1354     mlem->Draw("lego1Fb");
1355     gPad->Update();
1356   }
1357
1358   fxyMu[0][6] = fxyMu[1][6] = 9999;
1359   // Try to split into clusters
1360   Bool_t ok = kTRUE;
1361   if (mlem->GetSum() < 1) ok = kFALSE;
1362   else Split(mlem, coef);
1363   delete [] coef; delete [] probi; coef = 0; probi = 0;
1364   fPixArray->Delete(); 
1365   return ok;
1366 }
1367
1368 //_____________________________________________________________________________
1369 void AliMUONClusterFinderAZ::Mlem(Double_t *coef, Double_t *probi)
1370 {
1371   // Use MLEM to find pixel charges
1372   
1373   Int_t nPix = fPixArray->GetEntriesFast();
1374   Int_t npad = fnPads[0] + fnPads[1];
1375   Double_t *probi1 = new Double_t [nPix];
1376   Int_t indx, indx1;
1377   AliMUONPixel *pixPtr;
1378
1379   for (Int_t iter=0; iter<15; iter++) {
1380     // Do iterations
1381     for (Int_t ipix=0; ipix<nPix; ipix++) {
1382       // Correct each pixel
1383       if (probi[ipix] < 0.01) continue; // skip "invisible" pixel
1384       Double_t sum = 0;
1385       probi1[ipix] = probi[ipix];
1386       for (Int_t j=0; j<npad; j++) {
1387         if (fPadIJ[1][j] < 0) continue; 
1388         Double_t sum1 = 0;
1389         indx1 = j*nPix;
1390         indx = indx1 + ipix;
1391         for (Int_t i=0; i<nPix; i++) {
1392           // Caculate expectation
1393           pixPtr = (AliMUONPixel*) fPixArray->UncheckedAt(i);
1394           sum1 += pixPtr->Charge()*coef[indx1+i];
1395         } // for (Int_t i=0;
1396         if (fXyq[2][j] > fResponse->MaxAdc()-1 && sum1 > fResponse->MaxAdc()) { probi1[ipix] -= coef[indx]; continue; } // correct for pad charge overflows
1397         //cout << sum1 << " " << fXyq[2][j] << " " << coef[j*nPix+ipix] << endl;
1398         if (coef[indx] > 1.e-6) sum += fXyq[2][j]*coef[indx]/sum1;
1399       } // for (Int_t j=0;
1400       pixPtr = (AliMUONPixel*) fPixArray->UncheckedAt(ipix);
1401       if (probi1[ipix] > 1.e-6) pixPtr->SetCharge(pixPtr->Charge()*sum/probi1[ipix]);
1402     } // for (Int_t ipix=0;
1403   } // for (Int_t iter=0;
1404   delete [] probi1;
1405   return;
1406 }
1407
1408 //_____________________________________________________________________________
1409 void AliMUONClusterFinderAZ::FindCOG(TH2D *mlem, Double_t *xyc)
1410 {
1411   // Calculate position of the center-of-gravity around the maximum pixel
1412
1413   Int_t ixmax, iymax, ix, nsumx=0, nsumy=0, nsum=0;
1414   Int_t i1 = -9, j1 = -9;
1415   mlem->GetMaximumBin(ixmax,iymax,ix);
1416   Int_t nx = mlem->GetNbinsX();
1417   Int_t ny = mlem->GetNbinsY();
1418   Double_t thresh = mlem->GetMaximum()/10;
1419   Double_t x, y, cont, xq=0, yq=0, qq=0;
1420
1421   for (Int_t i=TMath::Max(1,iymax-1); i<=TMath::Min(ny,iymax+1); i++) {
1422     y = mlem->GetYaxis()->GetBinCenter(i);
1423     for (Int_t j=TMath::Max(1,ixmax-1); j<=TMath::Min(nx,ixmax+1); j++) {
1424       cont = mlem->GetCellContent(j,i);
1425       if (cont < thresh) continue;
1426       if (i != i1) {i1 = i; nsumy++;}
1427       if (j != j1) {j1 = j; nsumx++;}
1428       x = mlem->GetXaxis()->GetBinCenter(j);
1429       xq += x*cont;
1430       yq += y*cont;
1431       qq += cont;
1432       nsum++;
1433     }
1434   }
1435
1436   Double_t cmax = 0;
1437   Int_t i2 = 0, j2 = 0;
1438   x = y = 0;
1439   if (nsumy == 1) {
1440     // one bin in Y - add one more (with the largest signal)
1441     for (Int_t i=TMath::Max(1,iymax-1); i<=TMath::Min(ny,iymax+1); i++) {
1442       if (i == iymax) continue;
1443       for (Int_t j=TMath::Max(1,ixmax-1); j<=TMath::Min(nx,ixmax+1); j++) {
1444         cont = mlem->GetCellContent(j,i);
1445         if (cont > cmax) {
1446           cmax = cont;
1447           x = mlem->GetXaxis()->GetBinCenter(j);
1448           y = mlem->GetYaxis()->GetBinCenter(i);
1449           i2 = i;
1450           j2 = j;
1451         }
1452       }
1453     }
1454     xq += x*cmax;
1455     yq += y*cmax;
1456     qq += cmax;
1457     if (i2 != i1) nsumy++;
1458     if (j2 != j1) nsumx++;
1459     nsum++;
1460   } // if (nsumy == 1)
1461
1462   if (nsumx == 1) {
1463     // one bin in X - add one more (with the largest signal)
1464     cmax = x = y = 0;
1465     for (Int_t j=TMath::Max(1,ixmax-1); j<=TMath::Min(nx,ixmax+1); j++) {
1466       if (j == ixmax) continue;
1467       for (Int_t i=TMath::Max(1,iymax-1); i<=TMath::Min(ny,iymax+1); i++) {
1468         cont = mlem->GetCellContent(j,i);
1469         if (cont > cmax) {
1470           cmax = cont;
1471           x = mlem->GetXaxis()->GetBinCenter(j);
1472           y = mlem->GetYaxis()->GetBinCenter(i);
1473           i2 = i;
1474           j2 = j;
1475         }
1476       }
1477     }
1478     xq += x*cmax;
1479     yq += y*cmax;
1480     qq += cmax;
1481     if (i2 != i1) nsumy++;
1482     if (j2 != j1) nsumx++;
1483     nsum++;
1484   } // if (nsumx == 1)
1485
1486   xyc[0] = xq/qq; xyc[1] = yq/qq;
1487   cout << xyc[0] << " " << xyc[1] << " " << qq << " " << nsum << " " << nsumx << " " << nsumy << endl;
1488   return;
1489 }
1490
1491 //_____________________________________________________________________________
1492 Int_t AliMUONClusterFinderAZ::FindNearest(AliMUONPixel *pixPtr0)
1493 {
1494   // Find the pixel nearest to the given one
1495   // (algorithm may be not very efficient)
1496
1497   Int_t nPix = fPixArray->GetEntriesFast(), imin = 0;
1498   Double_t rmin = 99999, dx = 0, dy = 0, r = 0;
1499   Double_t xc = pixPtr0->Coord(0), yc = pixPtr0->Coord(1);
1500   AliMUONPixel *pixPtr;
1501
1502   for (Int_t i=0; i<nPix; i++) {
1503     pixPtr = (AliMUONPixel*) fPixArray->UncheckedAt(i);
1504     if (pixPtr->Charge() < 0.5) continue;
1505     dx = (xc - pixPtr->Coord(0)) / pixPtr->Size(0);
1506     dy = (yc - pixPtr->Coord(1)) / pixPtr->Size(1);
1507     r = dx *dx + dy * dy;
1508     if (r < rmin) { rmin = r; imin = i; }
1509   }
1510   return imin;
1511 }
1512
1513 //_____________________________________________________________________________
1514 void AliMUONClusterFinderAZ::Split(TH2D *mlem, Double_t *coef)
1515 {
1516   // The main steering function to work with clusters of pixels in anode
1517   // plane (find clusters, decouple them from each other, merge them (if
1518   // necessary), pick up coupled pads, call the fitting function)
1519   
1520   Int_t nx = mlem->GetNbinsX();
1521   Int_t ny = mlem->GetNbinsY();
1522   Int_t nPix = fPixArray->GetEntriesFast();
1523
1524   Bool_t *used = new Bool_t[ny*nx];
1525   Double_t cont;
1526   Int_t nclust = 0, indx, indx1;
1527
1528   for (Int_t i=0; i<ny*nx; i++) used[i] = kFALSE; 
1529
1530   TObjArray *clusters[200]={0};
1531   TObjArray *pix;
1532
1533   // Find clusters of histogram bins (easier to work in 2-D space)
1534   for (Int_t i=1; i<=ny; i++) {
1535     for (Int_t j=1; j<=nx; j++) {
1536       indx = (i-1)*nx + j - 1;
1537       if (used[indx]) continue;
1538       cont = mlem->GetCellContent(j,i);
1539       if (cont < 0.5) continue;
1540       pix = new TObjArray(20);
1541       used[indx] = 1;
1542       pix->Add(BinToPix(mlem,j,i));
1543       AddBin(mlem, i, j, 0, used, pix); // recursive call
1544       clusters[nclust++] = pix;
1545       if (nclust > 200) { cout << " Too many clusters " << endl; ::exit(0); }
1546     } // for (Int_t j=1; j<=nx; j++) {
1547   } // for (Int_t i=1; i<=ny;
1548   cout << nclust << endl;
1549   delete [] used; used = 0;
1550   
1551   // Compute couplings between clusters and clusters to pads
1552   Int_t npad = fnPads[0] + fnPads[1];
1553
1554   // Exclude pads with overflows
1555   for (Int_t j=0; j<npad; j++) {
1556     if (fXyq[2][j] > fResponse->MaxAdc()-1) fPadIJ[1][j] = -9;
1557     else fPadIJ[1][j] = 0;
1558   }
1559
1560   // Compute couplings of clusters to pads
1561   TMatrixD *aij_clu_pad = new TMatrixD(nclust,npad);
1562   *aij_clu_pad = 0;
1563   Int_t npxclu;
1564   for (Int_t iclust=0; iclust<nclust; iclust++) {
1565     pix = clusters[iclust];
1566     npxclu = pix->GetEntriesFast();
1567     for (Int_t i=0; i<npxclu; i++) {
1568       indx = fPixArray->IndexOf(pix->UncheckedAt(i));
1569       for (Int_t j=0; j<npad; j++) {
1570         // Exclude overflows
1571         if (fPadIJ[1][j] < 0) continue;
1572         if (coef[j*nPix+indx] < kCouplMin) continue;
1573         (*aij_clu_pad)(iclust,j) += coef[j*nPix+indx];
1574       }
1575     }
1576   }
1577   // Compute couplings between clusters
1578   TMatrixD *aij_clu_clu = new TMatrixD(nclust,nclust);
1579   *aij_clu_clu = 0;
1580   for (Int_t iclust=0; iclust<nclust; iclust++) {
1581     for (Int_t j=0; j<npad; j++) {
1582       // Exclude overflows
1583       if (fPadIJ[1][j] < 0) continue;
1584       if ((*aij_clu_pad)(iclust,j) < kCouplMin) continue;
1585       for (Int_t iclust1=iclust+1; iclust1<nclust; iclust1++) {
1586         if ((*aij_clu_pad)(iclust1,j) < kCouplMin) continue;
1587         (*aij_clu_clu)(iclust,iclust1) += 
1588           TMath::Sqrt ((*aij_clu_pad)(iclust,j)*(*aij_clu_pad)(iclust1,j));
1589       }
1590     }
1591   }
1592   for (Int_t iclust=0; iclust<nclust; iclust++) {
1593     for (Int_t iclust1=iclust+1; iclust1<nclust; iclust1++) {
1594       (*aij_clu_clu)(iclust1,iclust) = (*aij_clu_clu)(iclust,iclust1);
1595     }
1596   }
1597
1598   if (nclust > 1) aij_clu_clu->Print();
1599
1600   // Find groups of coupled clusters
1601   used = new Bool_t[nclust];
1602   for (Int_t i=0; i<nclust; i++) used[i] = kFALSE;
1603   Int_t *clustNumb = new Int_t[nclust];
1604   Int_t nCoupled, nForFit, minGroup[3], clustFit[3], nfit = 0;
1605   Double_t parOk[8];
1606
1607   for (Int_t igroup=0; igroup<nclust; igroup++) {
1608     if (used[igroup]) continue;
1609     used[igroup] = kTRUE;
1610     clustNumb[0] = igroup;
1611     nCoupled = 1;
1612     // Find group of coupled clusters
1613     AddCluster(igroup, nclust, aij_clu_clu, used, clustNumb, nCoupled); // recursive
1614     cout << " nCoupled: " << nCoupled << endl;
1615     for (Int_t i=0; i<nCoupled; i++) cout << clustNumb[i] << " "; cout << endl; 
1616
1617     while (nCoupled > 0) {
1618
1619       if (nCoupled < 4) {
1620         nForFit = nCoupled;
1621         for (Int_t i=0; i<nCoupled; i++) clustFit[i] = clustNumb[i];
1622       } else {
1623         // Too many coupled clusters to fit - try to decouple them
1624         // Find the lowest coupling of 1, 2, min(3,nLinks/2) pixels with 
1625         // all the others in the group 
1626         for (Int_t j=0; j<3; j++) minGroup[j] = -1;
1627         Double_t coupl = MinGroupCoupl(nCoupled, clustNumb, aij_clu_clu, minGroup);
1628
1629         // Flag clusters for fit
1630         nForFit = 0;
1631         while (minGroup[nForFit] >= 0 && nForFit < 3) {
1632           cout << clustNumb[minGroup[nForFit]] << " ";
1633           clustFit[nForFit] = clustNumb[minGroup[nForFit]];
1634           clustNumb[minGroup[nForFit]] -= 999;
1635           nForFit++;
1636         }
1637         cout << nForFit << " " << coupl << endl;
1638       } // else
1639
1640       // Select pads for fit. 
1641       if (SelectPad(nCoupled, nForFit, clustNumb, clustFit, aij_clu_pad) < 3 && nCoupled > 1) {
1642         // Deselect pads
1643         for (Int_t j=0; j<npad; j++) if (TMath::Abs(fPadIJ[1][j]) == 1) fPadIJ[1][j] = 0;
1644         // Merge the failed cluster candidates (with too few pads to fit) with 
1645         // the one with the strongest coupling
1646         Merge(nForFit, nCoupled, clustNumb, clustFit, clusters, aij_clu_clu, aij_clu_pad);
1647       } else {
1648         // Do the fit
1649         nfit = Fit(nForFit, clustFit, clusters, parOk);
1650       }
1651
1652       // Subtract the fitted charges from pads with strong coupling and/or
1653       // return pads for further use
1654       UpdatePads(nfit, parOk);
1655
1656       // Mark used pads
1657       for (Int_t j=0; j<npad; j++) {if (fPadIJ[1][j] == 1) fPadIJ[1][j] = -1;}
1658
1659       // Sort the clusters (move to the right the used ones)
1660       Int_t beg = 0, end = nCoupled - 1;
1661       while (beg < end) {
1662         if (clustNumb[beg] >= 0) { beg++; continue; }
1663         for (Int_t j=end; j>beg; j--) {
1664           if (clustNumb[j] < 0) continue;
1665           end = j - 1;
1666           indx = clustNumb[beg];
1667           clustNumb[beg] = clustNumb[j];
1668           clustNumb[j] = indx;
1669           break;
1670         }
1671         beg++;
1672       }
1673
1674       nCoupled -= nForFit;
1675       if (nCoupled > 3) {
1676         // Remove couplings of used clusters
1677         for (Int_t iclust=nCoupled; iclust<nCoupled+nForFit; iclust++) {
1678           indx = clustNumb[iclust] + 999;
1679           for (Int_t iclust1=0; iclust1<nCoupled; iclust1++) {
1680             indx1 = clustNumb[iclust1];
1681             (*aij_clu_clu)(indx,indx1) = (*aij_clu_clu)(indx1,indx) = 0;
1682           }
1683         }
1684
1685         // Update the remaining clusters couplings (exclude couplings from 
1686         // the used pads)
1687         for (Int_t j=0; j<npad; j++) {
1688           if (fPadIJ[1][j] != -1) continue;
1689           for (Int_t iclust=0; iclust<nCoupled; iclust++) {
1690             indx = clustNumb[iclust];
1691             if ((*aij_clu_pad)(indx,j) < kCouplMin) continue;
1692             for (Int_t iclust1=iclust+1; iclust1<nCoupled; iclust1++) {
1693               indx1 = clustNumb[iclust1];
1694               if ((*aij_clu_pad)(indx1,j) < kCouplMin) continue;
1695               // Check this
1696               (*aij_clu_clu)(indx,indx1) -= 
1697                 TMath::Sqrt ((*aij_clu_pad)(indx,j)*(*aij_clu_pad)(indx1,j));
1698               (*aij_clu_clu)(indx1,indx) = (*aij_clu_clu)(indx,indx1);
1699             }
1700           }
1701           fPadIJ[1][j] = -9;
1702         } // for (Int_t j=0; j<npad;
1703       } // if (nCoupled > 3)
1704     } // while (nCoupled > 0)
1705   } // for (Int_t igroup=0; igroup<nclust;
1706
1707   //delete aij_clu; aij_clu = 0; delete aij_clu_pad; aij_clu_pad = 0;
1708   aij_clu_clu->Delete(); aij_clu_pad->Delete();
1709   for (Int_t iclust=0; iclust<nclust; iclust++) {
1710     pix = clusters[iclust]; 
1711     pix->Clear();
1712     delete pix; pix = 0;
1713   }
1714   delete [] clustNumb; clustNumb = 0; delete [] used; used = 0;
1715 }
1716
1717 //_____________________________________________________________________________
1718 void AliMUONClusterFinderAZ::AddBin(TH2D *mlem, Int_t ic, Int_t jc, Int_t mode, Bool_t *used, TObjArray *pix)
1719 {
1720   // Add a bin to the cluster
1721
1722   Int_t nx = mlem->GetNbinsX();
1723   Int_t ny = mlem->GetNbinsY();
1724   Double_t cont1, cont = mlem->GetCellContent(jc,ic);
1725   AliMUONPixel *pixPtr = 0;
1726
1727   for (Int_t i=TMath::Max(ic-1,1); i<=TMath::Min(ic+1,ny); i++) {
1728     for (Int_t j=TMath::Max(jc-1,1); j<=TMath::Min(jc+1,nx); j++) {
1729       if (i != ic && j != jc) continue;
1730       if (used[(i-1)*nx+j-1]) continue;
1731       cont1 = mlem->GetCellContent(j,i);
1732       if (mode && cont1 > cont) continue;
1733       used[(i-1)*nx+j-1] = kTRUE;
1734       if (cont1 < 0.5) continue;
1735       if (pix) pix->Add(BinToPix(mlem,j,i)); 
1736       else {
1737         pixPtr = new AliMUONPixel (mlem->GetXaxis()->GetBinCenter(j), 
1738                                    mlem->GetYaxis()->GetBinCenter(i), 0, 0, cont1);
1739         fPixArray->Add((TObject*)pixPtr);
1740       }
1741       AddBin(mlem, i, j, mode, used, pix); // recursive call
1742     }
1743   }
1744 }
1745
1746 //_____________________________________________________________________________
1747 TObject* AliMUONClusterFinderAZ::BinToPix(TH2D *mlem, Int_t jc, Int_t ic)
1748 {
1749   // Translate histogram bin to pixel 
1750   
1751   Double_t yc = mlem->GetYaxis()->GetBinCenter(ic);
1752   Double_t xc = mlem->GetXaxis()->GetBinCenter(jc);
1753   
1754   Int_t nPix = fPixArray->GetEntriesFast();
1755   AliMUONPixel *pixPtr;
1756
1757   // Compare pixel and bin positions
1758   for (Int_t i=0; i<nPix; i++) {
1759     pixPtr = (AliMUONPixel*) fPixArray->UncheckedAt(i);
1760     if (pixPtr->Charge() < 0.5) continue;
1761     if (TMath::Abs(pixPtr->Coord(0)-xc)<1.e-4 && TMath::Abs(pixPtr->Coord(1)-yc)<1.e-4) return (TObject*) pixPtr;
1762   }
1763   cout << " Something wrong ??? " << endl;
1764   return NULL;
1765 }
1766
1767 //_____________________________________________________________________________
1768 void AliMUONClusterFinderAZ::AddCluster(Int_t ic, Int_t nclust, TMatrixD *aij_clu_clu, Bool_t *used, Int_t *clustNumb, Int_t &nCoupled)
1769 {
1770   // Add a cluster to the group of coupled clusters
1771
1772   for (Int_t i=0; i<nclust; i++) {
1773     if (used[i]) continue;
1774     if ((*aij_clu_clu)(i,ic) < kCouplMin) continue;
1775     used[i] = kTRUE;
1776     clustNumb[nCoupled++] = i;
1777     AddCluster(i, nclust, aij_clu_clu, used, clustNumb, nCoupled);
1778   }
1779 }
1780
1781 //_____________________________________________________________________________
1782 Double_t AliMUONClusterFinderAZ::MinGroupCoupl(Int_t nCoupled, Int_t *clustNumb, TMatrixD *aij_clu_clu, Int_t *minGroup)
1783 {
1784   // Find group of clusters with minimum coupling to all the others
1785
1786   Int_t i123max = TMath::Min(3,nCoupled/2); 
1787   Int_t indx, indx1, indx2, indx3, nTot = 0;
1788   Double_t *coupl1 = 0, *coupl2 = 0, *coupl3 = 0;
1789
1790   for (Int_t i123=1; i123<=i123max; i123++) {
1791
1792     if (i123 == 1) {
1793       coupl1 = new Double_t [nCoupled];
1794       for (Int_t i=0; i<nCoupled; i++) coupl1[i] = 0;
1795     }
1796     else if (i123 == 2) {
1797       nTot = nCoupled*nCoupled;
1798       coupl2 = new Double_t [nTot];
1799       for (Int_t i=0; i<nTot; i++) coupl2[i] = 9999;
1800     } else {
1801       nTot = nTot*nCoupled;
1802       coupl3 = new Double_t [nTot];
1803       for (Int_t i=0; i<nTot; i++) coupl3[i] = 9999;
1804     } // else
1805
1806     for (Int_t i=0; i<nCoupled; i++) {
1807       indx1 = clustNumb[i];
1808       for (Int_t j=i+1; j<nCoupled; j++) {
1809         indx2 = clustNumb[j];
1810         if (i123 == 1) {
1811           coupl1[i] += (*aij_clu_clu)(indx1,indx2);
1812           coupl1[j] += (*aij_clu_clu)(indx1,indx2);
1813         } 
1814         else if (i123 == 2) {
1815           indx = i*nCoupled + j;
1816           coupl2[indx] = coupl1[i] + coupl1[j];
1817           coupl2[indx] -= 2 * ((*aij_clu_clu)(indx1,indx2));
1818         } else {
1819           for (Int_t k=j+1; k<nCoupled; k++) {
1820             indx3 = clustNumb[k];
1821             indx = i*nCoupled*nCoupled + j*nCoupled + k;
1822             coupl3[indx] = coupl2[i*nCoupled+j] + coupl1[k];
1823             coupl3[indx] -= 2 * ((*aij_clu_clu)(indx1,indx3)+(*aij_clu_clu)(indx2,indx3));
1824           }
1825         } // else
1826       } // for (Int_t j=i+1;
1827     } // for (Int_t i=0;
1828   } // for (Int_t i123=1;
1829
1830   // Find minimum coupling
1831   Double_t couplMin = 9999;
1832   Int_t locMin = 0;
1833
1834   for (Int_t i123=1; i123<=i123max; i123++) {
1835     if (i123 == 1) {
1836       locMin = TMath::LocMin(nCoupled, coupl1);
1837       couplMin = coupl1[locMin];
1838       minGroup[0] = locMin;
1839       delete [] coupl1; coupl1 = 0;
1840     } 
1841     else if (i123 == 2) {
1842       locMin = TMath::LocMin(nCoupled*nCoupled, coupl2);
1843       if (coupl2[locMin] < couplMin) {
1844         couplMin = coupl2[locMin];
1845         minGroup[0] = locMin/nCoupled;
1846         minGroup[1] = locMin%nCoupled;
1847       }
1848       delete [] coupl2; coupl2 = 0;
1849     } else {
1850       locMin = TMath::LocMin(nTot, coupl3);
1851       if (coupl3[locMin] < couplMin) {
1852         couplMin = coupl3[locMin];
1853         minGroup[0] = locMin/nCoupled/nCoupled;
1854         minGroup[1] = locMin%(nCoupled*nCoupled)/nCoupled;
1855         minGroup[2] = locMin%nCoupled;
1856       }
1857       delete [] coupl3; coupl3 = 0;
1858     } // else
1859   } // for (Int_t i123=1;
1860   return couplMin;
1861 }
1862
1863 //_____________________________________________________________________________
1864 Int_t AliMUONClusterFinderAZ::SelectPad(Int_t nCoupled, Int_t nForFit, Int_t *clustNumb, Int_t *clustFit, TMatrixD *aij_clu_pad)
1865 {
1866   // Select pads for fit. If too many coupled clusters, find pads giving 
1867   // the strongest coupling with the rest of clusters and exclude them from the fit.
1868
1869   Int_t npad = fnPads[0] + fnPads[1];
1870   Double_t *pad_pix = 0;
1871
1872   if (nCoupled > 3) {
1873     pad_pix = new Double_t[npad];
1874     for (Int_t i=0; i<npad; i++) pad_pix[i] = 0; 
1875   }
1876
1877   Int_t nOK = 0, indx, indx1;
1878   for (Int_t iclust=0; iclust<nForFit; iclust++) {
1879     indx = clustFit[iclust];
1880     for (Int_t j=0; j<npad; j++) {
1881       if (fPadIJ[1][j] < 0) continue; // exclude overflows and used pads
1882       if ((*aij_clu_pad)(indx,j) < kCouplMin) continue;
1883       fPadIJ[1][j] = 1; // pad to be used in fit
1884       nOK++;
1885       if (nCoupled > 3) {
1886         // Check other clusters
1887         for (Int_t iclust1=0; iclust1<nCoupled; iclust1++) {
1888           indx1 = clustNumb[iclust1];
1889           if (indx1 < 0) continue;
1890           if ((*aij_clu_pad)(indx1,j) < kCouplMin) continue;
1891           pad_pix[j] += (*aij_clu_pad)(indx1,j);
1892         }
1893       } // if (nCoupled > 3)
1894     } // for (Int_t j=0; j<npad;
1895   } // for (Int_t iclust=0; iclust<nForFit
1896   if (nCoupled < 4) return nOK;
1897
1898   Double_t aaa = 0;
1899   for (Int_t j=0; j<npad; j++) {
1900     if (pad_pix[j] < kCouplMin) continue;
1901     cout << j << " " << pad_pix[j] << " "; 
1902     cout << fXyq[0][j] << " " << fXyq[1][j] << endl;
1903     aaa += pad_pix[j];
1904     fPadIJ[1][j] = -1; // exclude pads with strong coupling to the other clusters
1905     nOK--;
1906   }
1907   delete [] pad_pix; pad_pix = 0;
1908   return nOK;
1909 }
1910   
1911 //_____________________________________________________________________________
1912 void AliMUONClusterFinderAZ::Merge(Int_t nForFit, Int_t nCoupled, Int_t *clustNumb, Int_t *clustFit, TObjArray **clusters, TMatrixD *aij_clu_clu, TMatrixD *aij_clu_pad)
1913 {
1914   // Merge the group of clusters with the one having the strongest coupling with them
1915
1916   Int_t indx, indx1, npxclu, npxclu1, imax=0;
1917   TObjArray *pix, *pix1;
1918   Double_t couplMax;
1919
1920   for (Int_t icl=0; icl<nForFit; icl++) {
1921     indx = clustFit[icl];
1922     pix = clusters[indx];
1923     npxclu = pix->GetEntriesFast();
1924     couplMax = -1;
1925     for (Int_t icl1=0; icl1<nCoupled; icl1++) {
1926       indx1 = clustNumb[icl1];
1927       if (indx1 < 0) continue;
1928       if ((*aij_clu_clu)(indx,indx1) > couplMax) {
1929         couplMax = (*aij_clu_clu)(indx,indx1);
1930         imax = indx1;
1931       }
1932     } // for (Int_t icl1=0;
1933     /*if (couplMax < kCouplMin) {
1934       cout << " Oops " << couplMax << endl;
1935       aij_clu_clu->Print();
1936       cout << icl << " " << indx << " " << npxclu << " " << nLinks << endl;
1937       ::exit(0);
1938       }*/
1939     // Add to it
1940     pix1 = clusters[imax];
1941     npxclu1 = pix1->GetEntriesFast();
1942     // Add pixels 
1943     for (Int_t i=0; i<npxclu; i++) { pix1->Add(pix->UncheckedAt(i)); pix->RemoveAt(i); }
1944     cout << " New number of pixels: " << npxclu1 << " " << pix1->GetEntriesFast() << endl;
1945     //Add cluster-to-cluster couplings
1946     //aij_clu_clu->Print();
1947     for (Int_t icl1=0; icl1<nCoupled; icl1++) {
1948       indx1 = clustNumb[icl1];
1949       if (indx1 < 0 || indx1 == imax) continue;
1950       (*aij_clu_clu)(indx1,imax) += (*aij_clu_clu)(indx,indx1);
1951       (*aij_clu_clu)(imax,indx1) = (*aij_clu_clu)(indx1,imax);
1952     }
1953     (*aij_clu_clu)(indx,imax) = (*aij_clu_clu)(imax,indx) = 0;
1954     //aij_clu_clu->Print();
1955     //Add cluster-to-pad couplings
1956     for (Int_t j=0; j<fnPads[0]+fnPads[1]; j++) {
1957       if (fPadIJ[1][j] < 0) continue; // exclude overflows and used pads
1958       (*aij_clu_pad)(imax,j) += (*aij_clu_pad)(indx,j);
1959       (*aij_clu_pad)(indx,j) = 0;
1960     }
1961   } // for (Int_t icl=0; icl<nForFit;
1962 }
1963
1964 //_____________________________________________________________________________
1965 Int_t AliMUONClusterFinderAZ::Fit(Int_t nfit, Int_t *clustFit, TObjArray **clusters, Double_t *parOk)
1966 {
1967   // Find selected clusters to selected pad charges
1968   
1969   TH2D *mlem = (TH2D*) gROOT->FindObject("mlem");
1970   //Int_t nx = mlem->GetNbinsX();
1971   //Int_t ny = mlem->GetNbinsY();
1972   Double_t xmin = mlem->GetXaxis()->GetXmin() - mlem->GetXaxis()->GetBinWidth(1);
1973   Double_t xmax = mlem->GetXaxis()->GetXmax() + mlem->GetXaxis()->GetBinWidth(1);
1974   Double_t ymin = mlem->GetYaxis()->GetXmin() - mlem->GetYaxis()->GetBinWidth(1);
1975   Double_t ymax = mlem->GetYaxis()->GetXmax() + mlem->GetYaxis()->GetBinWidth(1);
1976   //Double_t qmin = 0, qmax = 1;
1977   Double_t step[3]={0.01,0.002,0.02};
1978
1979   Double_t cont, cmax = 0, xseed = 0, yseed = 0, errOk[8];
1980   TObjArray *pix;
1981   Int_t npxclu;
1982
1983   // Number of pads to use
1984   Int_t npads = 0;
1985   for (Int_t i=0; i<fnPads[0]+fnPads[1]; i++) {if (fPadIJ[1][i] == 1) npads++;}
1986   for (Int_t i=0; i<nfit; i++) {cout << i+1 << " " << clustFit[i] << " ";}
1987   cout << nfit << endl;
1988   cout << " Number of pads to fit: " << npads << endl;
1989   fNpar = 0;
1990   fQtot = 0;
1991   if (npads < 2) return 0; 
1992
1993   // Take cluster maxima as fitting seeds
1994   AliMUONPixel *pixPtr;
1995   Double_t xyseed[3][2], qseed[3];
1996   for (Int_t ifit=1; ifit<=nfit; ifit++) {
1997     cmax = 0;
1998     pix = clusters[clustFit[ifit-1]];
1999     npxclu = pix->GetEntriesFast();
2000     for (Int_t clu=0; clu<npxclu; clu++) {
2001       pixPtr = (AliMUONPixel*) pix->UncheckedAt(clu);
2002       cont = pixPtr->Charge();
2003       fQtot += cont;
2004       if (cont > cmax) { 
2005         cmax = cont; 
2006         xseed = pixPtr->Coord(0);
2007         yseed = pixPtr->Coord(1);
2008       }
2009     }
2010     xyseed[ifit-1][0] = xseed;
2011     xyseed[ifit-1][1] = yseed;
2012     qseed[ifit-1] = cmax;
2013   } // for (Int_t ifit=1;
2014
2015   Int_t nDof, maxSeed[3];
2016   Double_t fmin, chi2o = 9999, chi2n;
2017
2018   // Try to fit with one-track hypothesis, then 2-track. If chi2/dof is 
2019   // lower, try 3-track (if number of pads is sufficient).
2020   
2021   TMath::Sort(nfit, qseed, maxSeed, kTRUE); // in decreasing order
2022   nfit = TMath::Min (nfit, (npads + 1) / 3);
2023
2024   Double_t *gin = 0, func0, func1, param[8], param0[2][8], deriv[2][8], step0[8];
2025   Double_t shift[8], stepMax, derMax, parmin[8], parmax[8], func2[2], shift0;
2026   Double_t delta[8], scMax, dder[8], estim, shiftSave = 0;
2027   Int_t min, max, nCall = 0, memory[8] = {0}, nLoop, idMax = 0, iestMax = 0, nFail;
2028
2029   for (Int_t iseed=0; iseed<nfit; iseed++) {
2030
2031     for (Int_t j=0; j<3; j++) step0[fNpar+j] = shift[fNpar+j] = step[j];
2032     param[fNpar] = xyseed[maxSeed[iseed]][0];
2033     parmin[fNpar] = xmin; 
2034     parmax[fNpar++] = xmax; 
2035     param[fNpar] = xyseed[maxSeed[iseed]][1];
2036     parmin[fNpar] = ymin; 
2037     parmax[fNpar++] = ymax; 
2038     if (fNpar > 2) {
2039       param[fNpar] = fNpar == 4 ? 0.5 : 0.3;
2040       parmin[fNpar] = 0; 
2041       parmax[fNpar++] = 1; 
2042     }
2043
2044     // Try new algorithm
2045     min = nLoop = 1; stepMax = func2[1] = derMax = 999999; nFail = 0;
2046
2047     while (1) {
2048       max = !min;
2049       fcn1(fNpar, gin, func0, param, 1); nCall++;
2050       //cout << " Func: " << func0 << endl;
2051
2052       func2[max] = func0;
2053       for (Int_t j=0; j<fNpar; j++) {
2054         param0[max][j] = param[j];
2055         delta[j] = step0[j];
2056         param[j] += delta[j] / 10;
2057         if (j > 0) param[j-1] -= delta[j-1] / 10;
2058         fcn1(fNpar, gin, func1, param, 1); nCall++;
2059         deriv[max][j] = (func1 - func0) / delta[j] * 10; // first derivative
2060         //cout << j << " " << deriv[max][j] << endl;
2061         dder[j] = param0[0][j] != param0[1][j] ? (deriv[0][j] - deriv[1][j]) / 
2062                                                  (param0[0][j] - param0[1][j]) : 0; // second derivative
2063       }
2064       param[fNpar-1] -= delta[fNpar-1] / 10;
2065       if (nCall > 2000) ::exit(0);
2066
2067       min = func2[0] < func2[1] ? 0 : 1;
2068       nFail = min == max ? 0 : nFail + 1;
2069
2070       stepMax = derMax = estim = 0;
2071       for (Int_t j=0; j<fNpar; j++) { 
2072         // Estimated distance to minimum
2073         shift0 = shift[j];
2074         if (nLoop == 1) shift[j] = TMath::Sign (step0[j], -deriv[max][j]); // first step
2075         else if (TMath::Abs(deriv[0][j]) < 1.e-3 && TMath::Abs(deriv[1][j]) < 1.e-3) shift[j] = 0;
2076         else if (deriv[min][j]*deriv[!min][j] > 0 && TMath::Abs(deriv[min][j]) > TMath::Abs(deriv[!min][j])
2077               || TMath::Abs(deriv[0][j]-deriv[1][j]) < 1.e-3) {
2078           shift[j] = -TMath::Sign (shift[j], (func2[0]-func2[1]) * (param0[0][j]-param0[1][j]));
2079           if (min == max) { 
2080             if (memory[j] > 1) { shift[j] *= 2; } //cout << " Memory " << memory[j] << " " << shift[j] << endl; }
2081             memory[j]++;
2082           }
2083         } else {
2084           shift[j] = -deriv[min][j] / dder[j];
2085           memory[j] = 0;
2086         }
2087         if (TMath::Abs(shift[j])/step0[j] > estim) { 
2088           estim = TMath::Abs(shift[j])/step0[j];
2089           iestMax = j;
2090         }
2091
2092         // Too big step
2093         if (TMath::Abs(shift[j])/step0[j] > 10) shift[j] = TMath::Sign(10.,shift[j]) * step0[j]; // 
2094
2095         // Failed to improve minimum
2096         if (min != max) {
2097           memory[j] = 0;
2098           param[j] = param0[min][j];
2099           if (TMath::Abs(shift[j]+shift0) > 0.1*step0[j]) shift[j] = (shift[j] + shift0) / 2;
2100           else shift[j] /= -2;
2101         } 
2102
2103         // Too big step
2104         if (TMath::Abs(shift[j]*deriv[min][j]) > func2[min]) 
2105           shift[j] = TMath::Sign (func2[min]/deriv[min][j], shift[j]);
2106
2107         // Introduce step relaxation factor
2108         if (memory[j] < 3) {
2109           scMax = 1 + 4 / TMath::Max(nLoop/2.,1.);
2110           if (TMath::Abs(shift0) > 0 && TMath::Abs(shift[j]/shift0) > scMax) 
2111             shift[j] = TMath::Sign (shift0*scMax, shift[j]);
2112         }
2113         param[j] += shift[j]; 
2114           
2115         //cout << " xxx " << j << " " << shift[j] << " " << param[j] << endl;
2116         stepMax = TMath::Max (stepMax, TMath::Abs(shift[j]/step0[j]));
2117         if (TMath::Abs(deriv[min][j]) > derMax) {
2118           idMax = j;
2119           derMax = TMath::Abs (deriv[min][j]);
2120         }
2121       } // for (Int_t j=0; j<fNpar;
2122       //cout << max << " " << func2[min] << " " << derMax << " " << stepMax << " " << estim << " " << iestMax << " " << nCall << endl;
2123       if (estim < 1 && derMax < 2 || nLoop > 100) break; // minimum was found
2124
2125       nLoop++;
2126       // Check for small step
2127       if (shift[idMax] == 0) { shift[idMax] = step0[idMax]/10; param[idMax] += shift[idMax]; continue; }
2128       if (!memory[idMax] && derMax > 0.5 && nLoop > 10) {
2129         //cout << " ok " << deriv[min][idMax] << " " << deriv[!min][idMax] << " " << dder[idMax]*shift[idMax] << " " << shift[idMax] << endl;
2130         if (dder[idMax] != 0 && TMath::Abs(deriv[min][idMax]/dder[idMax]/shift[idMax]) > 10) {
2131           if (min == max) dder[idMax] = -dder[idMax];
2132           shift[idMax] = -deriv[min][idMax] / dder[idMax] / 10; 
2133           param[idMax] += shift[idMax];
2134           stepMax = TMath::Max (stepMax, TMath::Abs(shift[idMax])/step0[idMax]);
2135           //cout << shift[idMax] << " " << param[idMax] << endl;
2136           if (min == max) shiftSave = shift[idMax];
2137         }
2138         if (nFail > 10) {
2139           param[idMax] -= shift[idMax];
2140           shift[idMax] = 4 * shiftSave * (gRandom->Rndm(0) - 0.5);
2141           param[idMax] += shift[idMax];
2142           //cout << shift[idMax] << endl;
2143         }
2144       }      
2145     } // while (1)
2146     fmin = func2[min];
2147
2148     nDof = npads - fNpar;
2149     chi2n = nDof ? fmin/nDof : 0;
2150
2151     if (chi2n*1.2+1.e-6 > chi2o ) { fNpar -= 3; break; }
2152     // Save parameters and errors
2153     for (Int_t i=0; i<fNpar; i++) {
2154       parOk[i] = param0[min][i];
2155       errOk[i] = fmin;
2156     }
2157
2158     cout << chi2o << " " << chi2n << endl;
2159     chi2o = chi2n;
2160     if (fmin < 0.1) break; // !!!???
2161   } // for (Int_t iseed=0; 
2162
2163   for (Int_t i=0; i<fNpar; i++) {
2164     if (i == 4 || i == 7) continue;
2165     cout << parOk[i] << " " << errOk[i] << endl;
2166   }
2167   nfit = (fNpar + 1) / 3;
2168   Double_t rad;
2169   Int_t indx, imax;
2170   if (fReco) {
2171     for (Int_t j=0; j<nfit; j++) {
2172       indx = j<2 ? j*2 : j*2+1;  
2173       AddRawCluster (parOk[indx], parOk[indx+1], errOk[indx]);
2174     }
2175     return nfit;
2176   } 
2177   for (Int_t i=0; i<fnMu; i++) {
2178     cmax = fxyMu[i][6];
2179     for (Int_t j=0; j<nfit; j++) {
2180       indx = j<2 ? j*2 : j*2+1;  
2181       rad = (fxyMu[i][0]-parOk[indx])*(fxyMu[i][0]-parOk[indx]) +
2182             (fxyMu[i][1]-parOk[indx+1])*(fxyMu[i][1]-parOk[indx+1]);
2183       if (rad < cmax) {
2184         cmax = rad; 
2185         imax = indx;
2186         fxyMu[i][6] = cmax;
2187         fxyMu[i][2] = parOk[imax] - fxyMu[i][0];
2188         fxyMu[i][4] = parOk[imax+1] - fxyMu[i][1];
2189         fxyMu[i][3] = errOk[imax];
2190         fxyMu[i][5] = errOk[imax+1];
2191       }
2192     }      
2193   }
2194   return nfit;
2195 }  
2196
2197 //_____________________________________________________________________________
2198 void fcn1(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag)
2199 {
2200   // Fit for one track
2201   AliMUONClusterFinderAZ& c = *(AliMUONClusterFinderAZ::fgClusterFinder);    
2202   
2203   Int_t cath, ix, iy, indx, npads=0;
2204   Double_t charge, delta, coef=0, chi2=0;
2205   for (Int_t j=0; j<c.fnPads[0]+c.fnPads[1]; j++) {
2206     if (c.fPadIJ[1][j] != 1) continue;
2207     cath = c.fPadIJ[0][j];
2208     npads++;
2209     c.fSegmentation[cath]->GetPadI(c.fXyq[0][j],c.fXyq[1][j],c.fZpad,ix,iy);
2210     c.fSegmentation[cath]->SetPad(ix,iy);
2211     charge = 0;
2212     for (Int_t i=c.fNpar/3; i>=0; i--) { // sum over tracks
2213       indx = i<2 ? 2*i : 2*i+1;
2214       c.fSegmentation[cath]->SetHit(par[indx],par[indx+1],c.fZpad);
2215       //charge += c.fResponse->IntXY(c.fSegmentation[cath])*par[icl*3+2];
2216       if (c.fNpar == 2) coef = 1;
2217       else coef = i==c.fNpar/3 ? par[indx+2] : 1-coef;
2218       //coef = TMath::Max (coef, 0.);
2219       if (c.fNpar == 8 && i < 2) coef = i==1 ? coef*par[indx+2] : coef - par[7];
2220       //coef = TMath::Max (coef, 0.);
2221       charge += c.fResponse->IntXY(c.fSegmentation[cath])*coef;
2222     }
2223     charge *= c.fQtot;
2224     //if (c.fXyq[2][j] > c.fResponse->MaxAdc()-1 && charge > 
2225     //  c.fResponse->MaxAdc()) charge = c.fResponse->MaxAdc(); 
2226     delta = charge - c.fXyq[2][j];
2227     delta /= TMath::Sqrt ((Double_t)c.fXyq[2][j]);
2228     //chi2 += TMath::Abs(delta);
2229     chi2 += delta*delta;
2230   } // for (Int_t j=0;
2231   f = chi2; 
2232   Double_t qAver = c.fQtot/npads; //(c.fnPads[0]+c.fnPads[1]);
2233   f = chi2/qAver;
2234 }
2235
2236 //_____________________________________________________________________________
2237 void AliMUONClusterFinderAZ::UpdatePads(Int_t nfit, Double_t *par)
2238 {
2239   // Subtract the fitted charges from pads with strong coupling
2240
2241   Int_t cath, ix, iy, indx;
2242   Double_t charge, coef=0;
2243   for (Int_t j=0; j<fnPads[0]+fnPads[1]; j++) {
2244     if (fPadIJ[1][j] != -1) continue;
2245     if (fNpar != 0) {
2246       cath = fPadIJ[0][j];
2247       fSegmentation[cath]->GetPadI(fXyq[0][j],fXyq[1][j],fZpad,ix,iy);
2248       fSegmentation[cath]->SetPad(ix,iy);
2249       charge = 0;
2250       for (Int_t i=fNpar/3; i>=0; i--) { // sum over tracks
2251         indx = i<2 ? 2*i : 2*i+1;
2252         fSegmentation[cath]->SetHit(par[indx],par[indx+1],fZpad);
2253         if (fNpar == 2) coef = 1;
2254         else coef = i==fNpar/3 ? par[indx+2] : 1-coef;
2255         if (fNpar == 8 && i < 2) coef = i==1 ? coef*par[indx+2] : coef - par[7];
2256         charge += fResponse->IntXY(fSegmentation[cath])*coef;
2257       }
2258       charge *= fQtot;
2259       fXyq[2][j] -= charge;
2260     } // if (fNpar != 0)
2261     if (fXyq[2][j] > fResponse->ZeroSuppression()) fPadIJ[1][j] = 0; // return pad for further using
2262   } // for (Int_t j=0;
2263 }  
2264
2265 //_____________________________________________________________________________
2266 Bool_t AliMUONClusterFinderAZ::TestTrack(Int_t t) {
2267 // Test if track was user selected
2268   return kTRUE;
2269   /*
2270     if (fTrack[0]==-1 || fTrack[1]==-1) {
2271         return kTRUE;
2272     } else if (t==fTrack[0] || t==fTrack[1]) {
2273         return kTRUE;
2274     } else {
2275         return kFALSE;
2276     }
2277   */
2278 }
2279
2280 //_____________________________________________________________________________
2281 void AliMUONClusterFinderAZ::AddRawCluster(Double_t x, Double_t y, Double_t fmin)
2282 {
2283   //
2284   // Add a raw cluster copy to the list
2285   //
2286   AliMUONRawCluster cnew;
2287   AliMUON *pMUON=(AliMUON*)gAlice->GetModule("MUON");
2288   //pMUON->AddRawCluster(fInput->Chamber(),c); 
2289
2290   Int_t cath;    
2291   for (cath=0; cath<2; cath++) {
2292     cnew.fX[cath] = x;
2293     cnew.fY[cath] = y;
2294     cnew.fZ[cath] = fZpad;
2295     cnew.fQ[cath] = 100;
2296     cnew.fPeakSignal[cath] = 20;
2297     cnew.fMultiplicity[cath] = 5;
2298     cnew.fNcluster[cath] = 1;
2299     cnew.fChi2[cath] = fmin; //0.1;
2300     /*
2301     cnew.fMultiplicity[cath]=c->fMultiplicity[cath];
2302     for (i=0; i<fMul[cath]; i++) {
2303       cnew.fIndexMap[i][cath]=c->fIndexMap[i][cath];
2304       fSeg[cath]->SetPad(fIx[i][cath], fIy[i][cath]);
2305     }
2306     fprintf(stderr,"\nRawCluster %d cath %d\n",ico,cath);
2307     fprintf(stderr,"mult_av %d\n",c->fMultiplicity[cath]);
2308     FillCluster(&cnew,cath);
2309     */
2310   } 
2311   //cnew.fClusterType=cnew.PhysicsContribution();
2312   pMUON->AddRawCluster(AliMUONClusterInput::Instance()->Chamber(),cnew); 
2313   //fNPeaks++;
2314 }
2315
2316 //_____________________________________________________________________________
2317 Int_t AliMUONClusterFinderAZ::FindLocalMaxima(Int_t *localMax, Double_t *maxVal)
2318 {
2319   // Find local maxima in pixel space for large preclusters in order to
2320   // try to split them into smaller pieces (to speed up the MLEM procedure)
2321
2322   TH2D *hist = (TH2D*) gROOT->FindObject("anode");
2323   if (hist) hist->Delete();
2324
2325   Double_t xylim[4] = {999, 999, 999, 999};
2326   Int_t nPix = fPixArray->GetEntriesFast();
2327   AliMUONPixel *pixPtr = 0;
2328   for (Int_t ipix=0; ipix<nPix; ipix++) {
2329     pixPtr = (AliMUONPixel*) fPixArray->UncheckedAt(ipix);
2330     for (Int_t i=0; i<4; i++) 
2331          xylim[i] = TMath::Min (xylim[i], (i%2 ? -1 : 1)*pixPtr->Coord(i/2));
2332   }
2333   for (Int_t i=0; i<4; i++) xylim[i] -= pixPtr->Size(i/2); 
2334
2335   Int_t nx = TMath::Nint ((-xylim[1]-xylim[0])/pixPtr->Size(0)/2);
2336   Int_t ny = TMath::Nint ((-xylim[3]-xylim[2])/pixPtr->Size(1)/2);
2337   hist = new TH2D("anode","anode",nx,xylim[0],-xylim[1],ny,xylim[2],-xylim[3]);
2338   for (Int_t ipix=0; ipix<nPix; ipix++) {
2339     pixPtr = (AliMUONPixel*) fPixArray->UncheckedAt(ipix);
2340     hist->Fill(pixPtr->Coord(0), pixPtr->Coord(1), pixPtr->Charge());
2341   }
2342   if (fDraw) {
2343     ((TCanvas*)gROOT->FindObject("c2"))->cd();
2344     gPad->SetTheta(55);
2345     gPad->SetPhi(30);
2346     hist->Draw("lego1Fb");
2347     gPad->Update();
2348     int ia;
2349     cin >> ia;
2350   }
2351
2352   Int_t nMax = 0, indx;
2353   Int_t *isLocalMax = new Int_t[ny*nx];
2354   for (Int_t i=0; i<ny*nx; i++) isLocalMax[i] = 0;
2355
2356   for (Int_t i=1; i<=ny; i++) {
2357     indx = (i-1) * nx;
2358     for (Int_t j=1; j<=nx; j++) {
2359       if (hist->GetCellContent(j,i) < 0.5) continue;
2360       //if (isLocalMax[indx+j-1] < 0) continue;
2361       if (isLocalMax[indx+j-1] != 0) continue;
2362       FlagLocalMax(hist, i, j, isLocalMax);
2363     }
2364   }
2365
2366   for (Int_t i=1; i<=ny; i++) {
2367     indx = (i-1) * nx;
2368     for (Int_t j=1; j<=nx; j++) {
2369       if (isLocalMax[indx+j-1] > 0) { 
2370         localMax[nMax] = indx + j - 1; 
2371         maxVal[nMax++] = hist->GetCellContent(j,i);
2372       }
2373       if (nMax > 99) { cout << " Too many local maxima !!!" << endl; ::exit(0); }
2374     }
2375   }
2376   cout << " Local max: " << nMax << endl;
2377   delete [] isLocalMax; isLocalMax = 0;
2378   return nMax;
2379 }
2380
2381 //_____________________________________________________________________________
2382 void AliMUONClusterFinderAZ::FlagLocalMax(TH2D *hist, Int_t i, Int_t j, Int_t *isLocalMax)
2383 {
2384   // Flag pixels (whether or not local maxima)
2385
2386   Int_t nx = hist->GetNbinsX();
2387   Int_t ny = hist->GetNbinsY();
2388   Int_t cont = TMath::Nint (hist->GetCellContent(j,i));
2389   Int_t cont1 = 0;
2390
2391   for (Int_t i1=i-1; i1<i+2; i1++) {
2392     if (i1 < 1 || i1 > ny) continue;
2393     for (Int_t j1=j-1; j1<j+2; j1++) {
2394       if (j1 < 1 || j1 > nx) continue;
2395       if (i == i1 && j == j1) continue;
2396       cont1 = TMath::Nint (hist->GetCellContent(j1,i1));
2397       if (cont < cont1) { isLocalMax[(i-1)*nx+j-1] = -1; return; }
2398       else if (cont > cont1) isLocalMax[(i1-1)*nx+j1-1] = -1;
2399       else { // the same charge
2400         isLocalMax[(i-1)*nx+j-1] = 1; 
2401         if (isLocalMax[(i1-1)*nx+j1-1] == 0) {
2402           FlagLocalMax(hist, i1, j1, isLocalMax);
2403           if (isLocalMax[(i1-1)*nx+j1-1] < 0) { isLocalMax[(i-1)*nx+j-1] = -1; return; }
2404           else isLocalMax[(i1-1)*nx+j1-1] = -1;
2405         }
2406       } 
2407     }
2408   }
2409   isLocalMax[(i-1)*nx+j-1] = 1; // local maximum
2410 }
2411
2412 //_____________________________________________________________________________
2413 void AliMUONClusterFinderAZ::FindCluster(Int_t *localMax, Int_t iMax)
2414 {
2415   // Find pixel cluster around local maximum #iMax and pick up pads
2416   // overlapping with it
2417
2418   TH2D *hist = (TH2D*) gROOT->FindObject("anode");
2419   Int_t nx = hist->GetNbinsX();
2420   Int_t ny = hist->GetNbinsY();
2421   Int_t ic = localMax[iMax] / nx + 1;
2422   Int_t jc = localMax[iMax] % nx + 1;
2423   Bool_t *used = new Bool_t[ny*nx];
2424   for (Int_t i=0; i<ny*nx; i++) used[i] = kFALSE;
2425
2426   // Drop all pixels from the array - pick up only the ones from the cluster
2427   fPixArray->Delete();
2428
2429   Double_t wx = hist->GetXaxis()->GetBinWidth(1)/2; 
2430   Double_t wy = hist->GetYaxis()->GetBinWidth(1)/2;  
2431   Double_t yc = hist->GetYaxis()->GetBinCenter(ic);
2432   Double_t xc = hist->GetXaxis()->GetBinCenter(jc);
2433   Double_t cont = hist->GetCellContent(jc,ic);
2434   AliMUONPixel *pixPtr = new AliMUONPixel (xc, yc, wx, wy, cont);
2435   fPixArray->Add((TObject*)pixPtr);
2436   used[(ic-1)*nx+jc-1] = kTRUE;
2437   AddBin(hist, ic, jc, 1, used, (TObjArray*)0); // recursive call
2438
2439   Int_t nPix = fPixArray->GetEntriesFast(), npad = fnPads[0] + fnPads[1];
2440   for (Int_t i=0; i<nPix; i++) {
2441     ((AliMUONPixel*)fPixArray->UncheckedAt(i))->SetSize(0,wx); 
2442     ((AliMUONPixel*)fPixArray->UncheckedAt(i))->SetSize(1,wy); 
2443   }
2444   cout << iMax << " " << nPix << endl;
2445
2446   Float_t xy[4], xy12[4];
2447   // Pick up pads which overlap with found pixels
2448   for (Int_t i=0; i<npad; i++) fPadIJ[1][i] = -1;
2449   for (Int_t i=0; i<nPix; i++) {
2450     pixPtr = (AliMUONPixel*) fPixArray->UncheckedAt(i);
2451     for (Int_t j=0; j<4; j++) 
2452       xy[j] = pixPtr->Coord(j/2) + (j%2 ? 1 : -1)*pixPtr->Size(j/2);
2453     for (Int_t j=0; j<npad; j++) 
2454       if (Overlap(xy, j, xy12, 0)) fPadIJ[1][j] = 0; // flag for use
2455   }
2456
2457   delete [] used; used = 0;
2458 }