Fixes for #86059: Install data when ALICE_ROOT does not point to source (Christian)
[u/mrichter/AliRoot.git] / MUON / AliMUONClusterFinderPeakCOG.cxx
1 /**************************************************************************
2  * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
3  *                                                                        *
4  * Author: The ALICE Off-line Project.                                    *
5  * Contributors are mentioned in the code where appropriate.              *
6  *                                                                        *
7  * Permission to use, copy, modify and distribute this software and its   *
8  * documentation strictly for non-commercial purposes is hereby granted   *
9  * without fee, provided that the above copyright notice appears in all   *
10  * copies and that both the copyright notice and this permission notice   *
11  * appear in the supporting documentation. The authors make no claims     *
12  * about the suitability of this software for any purpose. It is          *
13  * provided "as is" without express or implied warranty.                  *
14  **************************************************************************/
15
16 /* $Id$ */
17
18 //-----------------------------------------------------------------------------
19 /// \class AliMUONClusterFinderPeakCOG
20 /// 
21 /// Clusterizer class based on simple peak finder
22 ///
23 /// Pre-clustering is handled by AliMUONPreClusterFinder
24 /// From a precluster a pixel array is built, and its local maxima are used
25 /// to get pads and compute pad center of gravity.
26 ///
27 /// \author Laurent Aphecetche (for the "new" C++ structure) and 
28 /// Alexander Zinchenko, JINR Dubna, for the hardcore of it ;-)
29 //-----------------------------------------------------------------------------
30
31 #include "AliMUONClusterFinderPeakCOG.h"
32 #include "AliMUONCluster.h"
33 #include "AliMUONPad.h"
34
35 #include "AliMpPad.h"
36 #include "AliMpVSegmentation.h"
37 #include "AliMpEncodePair.h"
38
39 #include "AliLog.h"
40 #include "AliRunLoader.h"
41 //#include "AliCodeTimer.h"
42
43 #include <Riostream.h>
44 #include <TH2.h>
45 //#include <TCanvas.h>
46 #include <TMath.h>
47
48 /// \cond CLASSIMP
49 ClassImp(AliMUONClusterFinderPeakCOG)
50 /// \endcond
51  
52 const Double_t AliMUONClusterFinderPeakCOG::fgkZeroSuppression = 6; // average zero suppression value
53 //const Double_t AliMUONClusterFinderMLEM::fgkDistancePrecision = 1e-6; // (cm) used to check overlaps and so on
54 const Double_t AliMUONClusterFinderPeakCOG::fgkDistancePrecision = 1e-3; // (cm) used to check overlaps and so on
55 const TVector2 AliMUONClusterFinderPeakCOG::fgkIncreaseSize(-AliMUONClusterFinderPeakCOG::fgkDistancePrecision,-AliMUONClusterFinderPeakCOG::fgkDistancePrecision);
56 const TVector2 AliMUONClusterFinderPeakCOG::fgkDecreaseSize(AliMUONClusterFinderPeakCOG::fgkDistancePrecision,AliMUONClusterFinderPeakCOG::fgkDistancePrecision);
57
58 // Status flags for pads
59 const Int_t AliMUONClusterFinderPeakCOG::fgkZero = 0x0; ///< pad "basic" state
60 const Int_t AliMUONClusterFinderPeakCOG::fgkMustKeep = 0x1; ///< do not kill (for pixels)
61 const Int_t AliMUONClusterFinderPeakCOG::fgkUseForFit = 0x10; ///< should be used for fit
62 const Int_t AliMUONClusterFinderPeakCOG::fgkOver = 0x100; ///< processing is over
63 const Int_t AliMUONClusterFinderPeakCOG::fgkModified = 0x1000; ///< modified pad charge 
64 const Int_t AliMUONClusterFinderPeakCOG::fgkCoupled = 0x10000; ///< coupled pad  
65
66 //_____________________________________________________________________________
67 AliMUONClusterFinderPeakCOG::AliMUONClusterFinderPeakCOG(Bool_t plot, AliMUONVClusterFinder* clusterFinder)
68   : AliMUONVClusterFinder(),
69 fPreClusterFinder(clusterFinder),
70 fPreCluster(0x0),
71 fClusterList(),
72 fEventNumber(0),
73 fDetElemId(-1),
74 fClusterNumber(0),
75 fHistAnode(0x0),
76 fPixArray(new TObjArray(20)),
77 fDebug(0),
78 fPlot(plot),
79 fNClusters(0),
80 fNAddVirtualPads(0)
81 {
82   /// Constructor
83  
84   fkSegmentation[1] = fkSegmentation[0] = 0x0; 
85
86   if (fPlot) fDebug = 1;
87 }
88
89 //_____________________________________________________________________________
90 AliMUONClusterFinderPeakCOG::~AliMUONClusterFinderPeakCOG()
91 {
92 /// Destructor
93   delete fPixArray; fPixArray = 0;
94   delete fPreClusterFinder;
95   AliInfo(Form("Total clusters %d AddVirtualPad needed %d",
96                fNClusters,fNAddVirtualPads));
97 }
98
99 //_____________________________________________________________________________
100 Bool_t 
101 AliMUONClusterFinderPeakCOG::Prepare(Int_t detElemId, TObjArray* pads[2],
102                                      const AliMpArea& area, 
103                                      const AliMpVSegmentation* seg[2])
104 {
105   /// Prepare for clustering
106 //  AliCodeTimerAuto("",0)
107   
108   for ( Int_t i = 0; i < 2; ++i )
109   {
110     fkSegmentation[i] = seg[i];
111   }
112   
113   // Find out the DetElemId
114   fDetElemId = detElemId;
115   
116   // find out current event number, and reset the cluster number
117   AliRunLoader *runLoader = AliRunLoader::Instance();
118   fEventNumber = runLoader ? runLoader->GetEventNumber() : 0;
119   fClusterNumber = -1;
120   fClusterList.Delete();
121   
122   AliDebug(3,Form("EVT %d DE %d",fEventNumber,fDetElemId));
123   
124   if ( fPreClusterFinder->NeedSegmentation() )
125   {
126     return fPreClusterFinder->Prepare(detElemId,pads,area,seg);
127   }
128   else
129   {
130     return fPreClusterFinder->Prepare(detElemId,pads,area);
131   }
132 }
133
134 //_____________________________________________________________________________
135 AliMUONCluster* 
136 AliMUONClusterFinderPeakCOG::NextCluster()
137 {
138   /// Return next cluster
139 //  AliCodeTimerAuto("",0)
140   
141   // if the list of clusters is not void, pick one from there
142   TObject* o = fClusterList.At(++fClusterNumber);
143   if ( o != 0x0 ) return static_cast<AliMUONCluster*>(o);
144   
145   //FIXME : at this point, must check whether we've used all the digits
146   //from precluster : if not, let the preclustering know about those unused
147   //digits, so it can reuse them
148   
149   // if the cluster list is exhausted, we need to go to the next
150   // pre-cluster and treat it
151
152   fClusterList.Delete(); // reset the list of clusters for this pre-cluster
153   fClusterNumber = -1; 
154     
155   fPreCluster = fPreClusterFinder->NextCluster();
156   
157   if (!fPreCluster)
158   {
159     // we are done
160     return 0x0;
161   }
162     
163   WorkOnPreCluster();
164
165   // WorkOnPreCluster may have used only part of the pads, so we check that
166   // now, and let the unused pads be reused by the preclustering...
167   
168   Int_t mult = fPreCluster->Multiplicity();
169   for ( Int_t i = 0; i < mult; ++i )
170   {
171     AliMUONPad* pad = fPreCluster->Pad(i);
172     if ( !pad->IsUsed() )
173     {
174       fPreClusterFinder->UsePad(*pad);
175     }
176   }
177   
178   return NextCluster();
179 }
180
181 //_____________________________________________________________________________
182 Bool_t
183 AliMUONClusterFinderPeakCOG::WorkOnPreCluster()
184 {
185   /// Starting from a precluster, builds a pixel array, and then
186   /// extract clusters from this array
187   
188   //  AliCodeTimerAuto("",0)
189
190   if (fDebug) {
191     cout << " *** Event # " << fEventNumber 
192          << " det. elem.: " << fDetElemId << endl;
193     for (Int_t j = 0; j < fPreCluster->Multiplicity(); ++j) {
194       AliMUONPad* pad = fPreCluster->Pad(j);
195       printf(" bbb %3d %1d %8.4f %8.4f %8.4f %8.4f %6.1f %3d %3d %2d %1d %1d \n",
196              j, pad->Cathode(), pad->Coord(0), pad->Coord(1), pad->DX()*2, pad->DY()*2,
197              pad->Charge(), pad->Ix(), pad->Iy(), pad->Status(), pad->IsReal(), pad->IsSaturated());
198     }
199   }
200
201   AliMUONCluster* cluster = CheckPrecluster(*fPreCluster);
202   if (!cluster) return kFALSE;
203
204   BuildPixArray(*cluster);
205   
206   if ( fPixArray->GetLast() < 0 )
207   {
208     AliDebug(1,"No pixel for the above cluster");
209     delete cluster;
210     return kFALSE;
211   }
212   
213   Int_t nMax = 1, localMax[100], maxPos[100] = {0};
214   Double_t maxVal[100];
215   
216   nMax = FindLocalMaxima(fPixArray, localMax, maxVal); // find local maxima
217
218   if (nMax > 1) TMath::Sort(nMax, maxVal, maxPos, kTRUE); // in descending order
219   
220   for (Int_t i = 0; i < nMax; ++i) 
221   {
222     FindCluster(*cluster,localMax, maxPos[i]);
223   }
224
225   delete cluster;
226   if (fPlot == 0) {
227     delete fHistAnode;
228     fHistAnode = 0x0;
229   }
230   return kTRUE;
231 }
232
233 //_____________________________________________________________________________
234 Bool_t 
235 AliMUONClusterFinderPeakCOG::Overlap(const AliMUONPad& pad, const AliMUONPad& pixel)
236 {
237   /// Check if the pad and the pixel overlaps
238
239   // make a fake pad from the pixel
240   AliMUONPad tmp(pad.DetElemId(),pad.Cathode(),pad.Ix(),pad.Iy(),
241                  pixel.Coord(0),pixel.Coord(1),
242                  pixel.Size(0),pixel.Size(1),0);
243   
244   return AliMUONPad::AreOverlapping(pad,tmp,fgkDecreaseSize);
245 }
246
247 //_____________________________________________________________________________
248 AliMUONCluster* 
249 AliMUONClusterFinderPeakCOG::CheckPrecluster(const AliMUONCluster& origCluster)
250 {
251   /// Check precluster in order to attempt to simplify it (mostly for
252   /// two-cathode preclusters)
253     
254   //  AliCodeTimerAuto("",0)
255
256   // Disregard small clusters (leftovers from splitting or noise)
257   if ((origCluster.Multiplicity()==1 || origCluster.Multiplicity()==2) &&
258       origCluster.Charge(0)+origCluster.Charge(1) < 1.525) // JC: adc -> fc
259   { 
260     return 0x0;
261   }
262
263   AliMUONCluster* cluster = new AliMUONCluster(origCluster);
264
265   AliDebug(2,"Start of CheckPreCluster=");
266   //StdoutToAliDebug(2,cluster->Print("full"));
267
268   AliMUONCluster* rv(0x0);
269   
270   if (cluster->Multiplicity(0) && cluster->Multiplicity(1))
271   { 
272     rv = CheckPreclusterTwoCathodes(cluster);
273   }
274   else
275   {
276     rv = cluster;
277   }
278   return rv;
279 }
280
281 //_____________________________________________________________________________
282 AliMUONCluster*
283 AliMUONClusterFinderPeakCOG::CheckPreclusterTwoCathodes(AliMUONCluster* cluster)
284 {
285   /// Check two-cathode cluster
286   
287   Int_t npad = cluster->Multiplicity();
288   Int_t* flags = new Int_t[npad];
289   for (Int_t j = 0; j < npad; ++j) flags[j] = 0;
290   
291   // Check pad overlaps
292   for ( Int_t i = 0; i < npad; ++i) 
293   {
294     AliMUONPad* padi = cluster->Pad(i);
295     if ( padi->Cathode() != 0 ) continue;
296     for (Int_t j = i+1; j < npad; ++j) 
297     {
298       AliMUONPad* padj = cluster->Pad(j);
299       if ( padj->Cathode() != 1 ) continue;
300       if ( !AliMUONPad::AreOverlapping(*padi,*padj,fgkDecreaseSize) ) continue;
301       flags[i] = flags[j] = 1; // mark overlapped pads
302     } 
303   } 
304   
305   // Check if all pads overlap
306   Int_t nFlags=0;
307   for (Int_t i = 0; i < npad; ++i) 
308   {
309     if (!flags[i]) ++nFlags;
310   }
311   
312   if (nFlags > 0) 
313   {
314     // not all pads overlap.
315     if (fDebug) cout << " nFlags: " << nFlags << endl;
316     TObjArray toBeRemoved;
317     for (Int_t i = 0; i < npad; ++i) 
318     {
319       AliMUONPad* pad = cluster->Pad(i);
320       if (flags[i]) continue;
321       Int_t cath = pad->Cathode();
322       Int_t cath1 = TMath::Even(cath);
323       // Check for edge effect (missing pads on the _other_ cathode)
324       AliMpPad mpPad 
325         = fkSegmentation[cath1]->PadByPosition(pad->Position().X(),pad->Position().Y(),kFALSE);
326       if (!mpPad.IsValid()) continue;
327       //if (nFlags == 1 && pad->Charge() < fgkZeroSuppression * 3) continue;
328       if (nFlags == 1 && pad->Charge() < 3.05) continue; // JC: adc -> fc
329       AliDebug(2,Form("Releasing the following pad : de,cath,ix,iy %d,%d,%d,%d charge %e",
330                       fDetElemId,pad->Cathode(),pad->Ix(),pad->Iy(),pad->Charge()));
331       toBeRemoved.AddLast(pad);
332       fPreCluster->Pad(i)->Release();
333     }
334     Int_t nRemove = toBeRemoved.GetEntriesFast();
335     for ( Int_t i = 0; i < nRemove; ++i )
336     {
337       cluster->RemovePad(static_cast<AliMUONPad*>(toBeRemoved.UncheckedAt(i)));
338     }
339   } 
340   
341   // Check correlations of cathode charges
342   if ( !cluster->IsSaturated() && cluster->ChargeAsymmetry() > 1 )
343   {
344     // big difference
345     Int_t cathode = cluster->MaxRawChargeCathode();
346     Int_t imin(-1);
347     Int_t imax(-1);
348     Double_t cmax(0);
349     Double_t cmin(1E9);
350     
351     // get min and max pad charges on the cathode opposite to the 
352     // max pad (given by MaxRawChargeCathode())
353     //
354     Int_t mult = cluster->Multiplicity();
355     for ( Int_t i = 0; i < mult; ++i )
356     {
357       AliMUONPad* pad = cluster->Pad(i);
358       if ( pad->Cathode() != cathode || !pad->IsReal() )
359       {
360         // only consider pads in the opposite cathode, and
361         // only consider real pads (i.e. exclude the virtual ones)
362         continue;
363       }
364       if ( pad->Charge() < cmin )
365       {
366         cmin = pad->Charge();
367         imin = i;
368         if (imax < 0) {
369           imax = imin;
370           cmax = cmin;
371         }
372       }
373       else if ( pad->Charge() > cmax )
374       {
375         cmax = pad->Charge();
376         imax = i;
377       }      
378     }
379     AliDebug(2,Form("Pad imin,imax %d,%d cmin,cmax %e,%e",
380                     imin,imax,cmin,cmax));
381     //
382     // arrange pads according to their distance to the max, normalized
383     // to the pad size
384     Double_t* dist = new Double_t[mult];
385     Double_t dxMin(1E9);
386     Double_t dyMin(1E9);
387     Double_t dmin(0);
388     
389     AliMUONPad* padmax = cluster->Pad(imax);
390     
391     for ( Int_t i = 0; i < mult; ++i )
392     {
393       dist[i] = 0.0;
394       if ( i == imax) continue;
395       AliMUONPad* pad = cluster->Pad(i);
396       if ( pad->Cathode() != cathode || !pad->IsReal() ) continue;
397       Double_t dx = (pad->X()-padmax->X())/padmax->DX()/2.0;
398       Double_t dy = (pad->Y()-padmax->Y())/padmax->DY()/2.0;
399       dist[i] = TMath::Sqrt(dx*dx+dy*dy);
400       if ( i == imin )
401       {
402         dmin = dist[i] + 1E-3; // distance to the pad with minimum charge
403         dxMin = dx;
404         dyMin = dy;
405       }      
406     }
407     
408     TMath::Sort(mult,dist,flags,kFALSE); // in ascending order
409     Double_t xmax(-1), distPrev(999);
410     TObjArray toBeRemoved;
411     
412     for ( Int_t i = 0; i < mult; ++i )
413     {
414       Int_t indx = flags[i];
415       AliMUONPad* pad = cluster->Pad(indx);
416       if ( pad->Cathode() != cathode || !pad->IsReal() ) continue;
417       if ( dist[indx] > dmin )
418       {
419         // farther than the minimum pad
420         Double_t dx = (pad->X()-padmax->X())/padmax->DX()/2.0;
421         Double_t dy = (pad->Y()-padmax->Y())/padmax->DY()/2.0;
422         dx *= dxMin;
423         dy *= dyMin;
424         if (dx >= 0 && dy >= 0) continue;
425         if (TMath::Abs(dx) > TMath::Abs(dy) && dx >= 0) continue;
426         if (TMath::Abs(dy) > TMath::Abs(dx) && dy >= 0) continue;        
427       }
428       if (dist[indx] > distPrev + 1) break; // overstepping empty pads
429       if ( pad->Charge() <= cmax || TMath::Abs(dist[indx]-xmax) < 1E-3 )
430       {
431         // release pad
432         if (TMath::Abs(dist[indx]-xmax) < 1.e-3) 
433         {
434           cmax = TMath::Max(pad->Charge(),cmax);
435         }
436         else
437         {
438           cmax = pad->Charge();
439         }
440         xmax = dist[indx];
441         distPrev = dist[indx];
442         AliDebug(2,Form("Releasing the following pad : de,cath,ix,iy %d,%d,%d,%d charge %e",
443                         fDetElemId,pad->Cathode(),pad->Ix(),pad->Iy(),
444                         pad->Charge()));
445   
446         toBeRemoved.AddLast(pad);
447         fPreCluster->Pad(indx)->Release();
448       }
449     }
450     Int_t nRemove = toBeRemoved.GetEntriesFast();
451     for ( Int_t i = 0; i < nRemove; ++i )
452     {
453       cluster->RemovePad(static_cast<AliMUONPad*>(toBeRemoved.UncheckedAt(i)));
454     }    
455     delete[] dist;
456   } // if ( !cluster->IsSaturated() && 
457   
458   delete[] flags;
459   
460   AliDebug(2,"End of CheckPreClusterTwoCathodes=");
461   //StdoutToAliDebug(2,cluster->Print("full"));
462
463   return cluster;    
464 }
465
466 //_____________________________________________________________________________
467 void
468 AliMUONClusterFinderPeakCOG::CheckOverlaps()
469 {
470   /// For debug only : check if some pixels overlap...
471   
472   Int_t nPix = fPixArray->GetLast()+1;
473   Int_t dummy(0);
474   
475   for ( Int_t i = 0; i < nPix; ++i )
476   {
477     AliMUONPad* pixelI = Pixel(i);
478     AliMUONPad pi(dummy,dummy,dummy,dummy,
479                   pixelI->Coord(0),pixelI->Coord(1),
480                   pixelI->Size(0),pixelI->Size(1),0.0);
481     
482     for ( Int_t j = i+1; j < nPix; ++j )
483     {
484       AliMUONPad* pixelJ = Pixel(j);
485       AliMUONPad pj(dummy,dummy,dummy,dummy,
486                     pixelJ->Coord(0),pixelJ->Coord(1),
487                     pixelJ->Size(0),pixelJ->Size(1),0.0);  
488       AliMpArea area;
489       
490       if ( AliMUONPad::AreOverlapping(pi,pj,fgkDecreaseSize,area) )
491       {
492         AliInfo(Form("The following 2 pixels (%d and %d) overlap !",i,j));
493         /*
494         StdoutToAliInfo(pixelI->Print();
495                         cout << " Surface = " << pixelI->Size(0)*pixelI->Size(1)*4 << endl;
496                         pixelJ->Print();
497                         cout << " Surface = " << pixelJ->Size(0)*pixelJ->Size(1)*4 << endl;
498                         cout << " Area surface = " << area.Dimensions().X()*area.Dimensions().Y()*4 << endl;
499                         cout << "-------" << endl;
500                         );
501         */        
502       }
503     }    
504   }
505 }
506
507 //_____________________________________________________________________________
508 void AliMUONClusterFinderPeakCOG::BuildPixArray(AliMUONCluster& cluster)
509 {
510   /// Build pixel array 
511   
512   Int_t npad = cluster.Multiplicity();
513   if (npad<=0) 
514   {
515     AliWarning("Got no pad at all ?!");
516   }
517   
518   fPixArray->Delete();
519   BuildPixArrayOneCathode(cluster);
520   
521 //  StdoutToAliDebug(2,cout << "End of BuildPixelArray:" << endl;
522 //                   fPixArray->Print(););
523   //CheckOverlaps();//FIXME : this is for debug only. Remove it.
524 }
525
526 //_____________________________________________________________________________
527 void AliMUONClusterFinderPeakCOG::BuildPixArrayOneCathode(AliMUONCluster& cluster)
528 {
529   /// Build the pixel array
530
531 //  AliDebug(2,Form("cluster.Multiplicity=%d",cluster.Multiplicity()));
532
533   TVector2 dim = cluster.MinPadDimensions (-1, kFALSE);
534   Double_t width[2] = {dim.X(), dim.Y()}, xy0[2] = { 0.0, 0.0 };
535   Int_t found[2] = {0}, mult = cluster.Multiplicity();
536
537   for ( Int_t i = 0; i < mult; ++i) {
538     AliMUONPad* pad = cluster.Pad(i);
539     for (Int_t j = 0; j < 2; ++j) {
540       if (found[j] == 0 && TMath::Abs(pad->Size(j)-width[j]) < fgkDistancePrecision) { 
541         xy0[j] = pad->Coord(j);
542         found[j] = 1;
543       }
544     }
545     if (found[0] && found[1]) break;
546   }
547
548   Double_t min[2], max[2];
549   Int_t cath0 = 0, cath1 = 1;
550   if (cluster.Multiplicity(0) == 0) cath0 = 1;
551   else if (cluster.Multiplicity(1) == 0) cath1 = 0;
552
553   Double_t leftDownX, leftDownY;
554   cluster.Area(cath0).LeftDownCorner(leftDownX, leftDownY);
555   Double_t rightUpX, rightUpY; 
556   cluster.Area(cath0).RightUpCorner(rightUpX, rightUpY);
557   min[0] = leftDownX;
558   min[1] = leftDownY;
559   max[0] = rightUpX;
560   max[1] = rightUpY;
561   if (cath1 != cath0) {
562     cluster.Area(cath1).LeftDownCorner(leftDownX, leftDownY);
563     cluster.Area(cath1).RightUpCorner(rightUpX, rightUpY);
564     min[0] = TMath::Max (min[0], leftDownX);
565     min[1] = TMath::Max (min[1], leftDownY);
566     max[0] = TMath::Min (max[0], rightUpX);
567     max[1] = TMath::Min (max[1], rightUpY);
568   }
569
570   // Adjust limits
571   //width[0] /= 2; width[1] /= 2; // just for check
572   Int_t nbins[2];
573   for (Int_t i = 0; i < 2; ++i) {
574     Double_t dist = (min[i] - xy0[i]) / width[i] / 2;
575     if (TMath::Abs(dist) < 1.e-6) dist = -1.e-6;
576     min[i] = xy0[i] + (TMath::Nint(dist-TMath::Sign(1.e-6,dist)) 
577                        + TMath::Sign(0.5,dist)) * width[i] * 2;
578     nbins[i] = TMath::Nint ((max[i] - min[i]) / width[i] / 2);
579     if (nbins[i] == 0) ++nbins[i];
580     max[i] = min[i] + nbins[i] * width[i] * 2;
581     //cout << dist << " " << min[i] << " " << max[i] << " " << nbins[i] << endl;
582   }
583
584   // Book histogram
585   TH2D *hist1 = new TH2D ("Grid", "", nbins[0], min[0], max[0], nbins[1], min[1], max[1]);
586   TH2D *hist2 = new TH2D ("Entries", "", nbins[0], min[0], max[0], nbins[1], min[1], max[1]);
587   TAxis *xaxis = hist1->GetXaxis();
588   TAxis *yaxis = hist1->GetYaxis();
589
590   // Fill histogram
591   for ( Int_t i = 0; i < mult; ++i) {
592     AliMUONPad* pad = cluster.Pad(i);
593     Int_t ix0 = xaxis->FindBin(pad->X());
594     Int_t iy0 = yaxis->FindBin(pad->Y());
595     PadOverHist(0, ix0, iy0, pad, hist1, hist2);
596   }
597
598   // Store pixels
599   for (Int_t i = 1; i <= nbins[0]; ++i) {
600     Double_t x = xaxis->GetBinCenter(i);
601     for (Int_t j = 1; j <= nbins[1]; ++j) {
602       if (hist2->GetCellContent(i,j) < 0.01525) continue; // JC: adc -> fc
603       if (cath0 != cath1) {
604         // Two-sided cluster
605         Double_t cont = hist2->GetCellContent(i,j);
606         if (cont < 999.) continue;
607         if (cont-Int_t(cont/1000.)*1000. < 0.07625) continue; // JC: adc -> fc
608       }
609       //if (hist2->GetCellContent(i,j) < 1.1 && cluster.Multiplicity(0) && 
610       //  cluster.Multiplicity(1)) continue;
611       Double_t y = yaxis->GetBinCenter(j);
612       Double_t charge = hist1->GetCellContent(i,j);
613       AliMUONPad* pixPtr = new AliMUONPad(x, y, width[0], width[1], charge);
614       fPixArray->Add(pixPtr);
615     }  
616   }
617   /*
618   if (fPixArray->GetEntriesFast() == 1) {
619     // Split pixel into 2
620     AliMUONPad* pixPtr = static_cast<AliMUONPad*> (fPixArray->UncheckedAt(0));
621     pixPtr->SetSize(0,width[0]/2.);
622     pixPtr->Shift(0,-width[0]/4.);
623     pixPtr = new AliMUONPad(pixPtr->X()+width[0], pixPtr->Y(), width[0]/2., width[1], pixPtr->Charge());
624     fPixArray->Add(pixPtr);
625   }
626   */
627   //fPixArray->Print();
628   delete hist1;
629   delete hist2;
630 }
631
632 //_____________________________________________________________________________
633 void AliMUONClusterFinderPeakCOG::PadOverHist(Int_t idir, Int_t ix0, Int_t iy0, AliMUONPad *pad,
634                                               TH2D *hist1, TH2D *hist2)
635 {
636   /// "Span" pad over histogram in the direction idir
637
638   TAxis *axis = idir == 0 ? hist1->GetXaxis() : hist1->GetYaxis();
639   Int_t nbins = axis->GetNbins(), cath = pad->Cathode();
640   Double_t bin = axis->GetBinWidth(1), amask = TMath::Power(1000.,cath*1.);
641
642   Int_t nbinPad = (Int_t)(pad->Size(idir)/bin*2+fgkDistancePrecision) + 1; // number of bins covered by pad
643
644   for (Int_t i = 0; i < nbinPad; ++i) {
645     Int_t ixy = idir == 0 ? ix0 + i : iy0 + i;
646     if (ixy > nbins) break;
647     Double_t lowEdge = axis->GetBinLowEdge(ixy);
648     if (lowEdge + fgkDistancePrecision > pad->Coord(idir) + pad->Size(idir)) break;
649     if (idir == 0) PadOverHist(1, ixy, iy0, pad, hist1, hist2); // span in the other direction
650     else {
651       // Fill histogram
652       Double_t cont = pad->Charge();
653       if (hist2->GetCellContent(ix0, ixy) > 0.01525) // JC: adc -> fc
654         cont = TMath::Min (hist1->GetCellContent(ix0, ixy), cont) 
655              + TMath::Min (TMath::Max(hist1->GetCellContent(ix0, ixy),cont)*0.1, 1.525); // JC: adc -> fc
656       hist1->SetCellContent(ix0, ixy, cont);
657       hist2->SetCellContent(ix0, ixy, hist2->GetCellContent(ix0, ixy)+amask);
658     }
659   }
660
661   for (Int_t i = -1; i > -nbinPad; --i) {
662     Int_t ixy = idir == 0 ? ix0 + i : iy0 + i;
663     if (ixy < 1) break;
664     Double_t upEdge = axis->GetBinUpEdge(ixy);
665     if (upEdge - fgkDistancePrecision < pad->Coord(idir) - pad->Size(idir)) break;
666     if (idir == 0) PadOverHist(1, ixy, iy0, pad, hist1, hist2); // span in the other direction
667     else {
668       // Fill histogram
669       Double_t cont = pad->Charge();
670       if (hist2->GetCellContent(ix0, ixy) > 0.01525) // JC: adc -> fc
671         cont = TMath::Min (hist1->GetCellContent(ix0, ixy), cont)
672              + TMath::Min (TMath::Max(hist1->GetCellContent(ix0, ixy),cont)*0.1,1.525);  // JC: adc -> fc
673       hist1->SetCellContent(ix0, ixy, cont);
674       hist2->SetCellContent(ix0, ixy, hist2->GetCellContent(ix0, ixy)+amask);
675     }
676   }
677 }
678
679 //_____________________________________________________________________________
680 Int_t AliMUONClusterFinderPeakCOG::FindLocalMaxima(TObjArray *pixArray, Int_t *localMax, Double_t *maxVal)
681 {
682 /// Find local maxima in pixel space 
683
684   AliDebug(1,Form("nPix=%d",pixArray->GetLast()+1));
685
686   //TH2D *hist = NULL;
687   //delete ((TH2D*) gROOT->FindObject("anode"));
688   //if (pixArray == fPixArray) hist = (TH2D*) gROOT->FindObject("anode");
689   //else { hist = (TH2D*) gROOT->FindObject("anode1"); cout << hist << endl; }
690   //if (hist) hist->Delete();
691   delete fHistAnode;
692  
693   Double_t xylim[4] = {999, 999, 999, 999};
694
695   Int_t nPix = pixArray->GetEntriesFast();
696   
697   if ( nPix <= 0 ) return 0;
698   
699   AliMUONPad *pixPtr = 0;
700   for (Int_t ipix = 0; ipix < nPix; ++ipix) {
701     pixPtr = (AliMUONPad*) pixArray->UncheckedAt(ipix);
702     for (Int_t i = 0; i < 4; ++i) 
703          xylim[i] = TMath::Min (xylim[i], (i%2 ? -1 : 1)*pixPtr->Coord(i/2));
704   }
705   for (Int_t i = 0; i < 4; ++i) xylim[i] -= pixPtr->Size(i/2); 
706
707   Int_t nx = TMath::Nint ((-xylim[1]-xylim[0])/pixPtr->Size(0)/2);
708   Int_t ny = TMath::Nint ((-xylim[3]-xylim[2])/pixPtr->Size(1)/2);
709   if (pixArray == fPixArray) fHistAnode = new TH2D("anode","anode",nx,xylim[0],-xylim[1],ny,xylim[2],-xylim[3]);
710   else fHistAnode = new TH2D("anode1","anode1",nx,xylim[0],-xylim[1],ny,xylim[2],-xylim[3]);
711   for (Int_t ipix = 0; ipix < nPix; ++ipix) {
712     pixPtr = (AliMUONPad*) pixArray->UncheckedAt(ipix);
713     fHistAnode->Fill(pixPtr->Coord(0), pixPtr->Coord(1), pixPtr->Charge());
714   }
715 //  if (fDraw && pixArray == fPixArray) fDraw->DrawHist("c2", hist);
716
717   Int_t nMax = 0, indx, nxy = ny * nx;
718   Int_t *isLocalMax = new Int_t[nxy];
719   for (Int_t i = 0; i < nxy; ++i) isLocalMax[i] = 0; 
720
721   for (Int_t i = 1; i <= ny; ++i) {
722     indx = (i-1) * nx;
723     for (Int_t j = 1; j <= nx; ++j) {
724       if (fHistAnode->GetCellContent(j,i) < 0.07625) continue; // JC: adc -> fc
725       //if (isLocalMax[indx+j-1] < 0) continue;
726       if (isLocalMax[indx+j-1] != 0) continue;
727       FlagLocalMax(fHistAnode, i, j, isLocalMax);
728     }
729   }
730
731   for (Int_t i = 1; i <= ny; ++i) {
732     indx = (i-1) * nx;
733     for (Int_t j = 1; j <= nx; ++j) {
734       if (isLocalMax[indx+j-1] > 0) { 
735         localMax[nMax] = indx + j - 1; 
736         maxVal[nMax++] = fHistAnode->GetCellContent(j,i);
737         if (nMax > 99) break;
738       }
739     }
740     if (nMax > 99) {
741       AliError(" Too many local maxima !!!");
742       break;
743     }
744   }
745   if (fDebug) cout << " Local max: " << nMax << endl;
746   delete [] isLocalMax; 
747   return nMax;
748 }
749
750 //_____________________________________________________________________________
751 void AliMUONClusterFinderPeakCOG::FlagLocalMax(TH2D *hist, Int_t i, Int_t j, Int_t *isLocalMax)
752 {
753 /// Flag pixels (whether or not local maxima)
754
755   Int_t nx = hist->GetNbinsX();
756   Int_t ny = hist->GetNbinsY();
757   Int_t cont = TMath::Nint (hist->GetCellContent(j,i));
758   Int_t cont1 = 0, indx = (i-1)*nx+j-1, indx1 = 0, indx2 = 0;
759
760   Int_t ie = i + 2, je = j + 2;
761   for (Int_t i1 = i-1; i1 < ie; ++i1) {
762     if (i1 < 1 || i1 > ny) continue;
763     indx1 = (i1 - 1) * nx;
764     for (Int_t j1 = j-1; j1 < je; ++j1) {
765       if (j1 < 1 || j1 > nx) continue;
766       if (i == i1 && j == j1) continue;
767       indx2 = indx1 + j1 - 1;
768       cont1 = TMath::Nint (hist->GetCellContent(j1,i1));
769       if (cont < cont1) { isLocalMax[indx] = -1; return; }
770       else if (cont > cont1) isLocalMax[indx2] = -1;
771       else { // the same charge
772         isLocalMax[indx] = 1; 
773         if (isLocalMax[indx2] == 0) {
774           FlagLocalMax(hist, i1, j1, isLocalMax);
775           if (isLocalMax[indx2] < 0) { isLocalMax[indx] = -1; return; }
776           else isLocalMax[indx2] = -1;
777         }
778       } 
779     }
780   }
781   isLocalMax[indx] = 1; // local maximum
782 }
783
784 //_____________________________________________________________________________
785 void AliMUONClusterFinderPeakCOG::FindCluster(AliMUONCluster& cluster, 
786                                               const Int_t *localMax, Int_t iMax)
787 {
788 /// Find pixel cluster around local maximum \a iMax and pick up pads
789 /// overlapping with it
790
791   //TH2D *hist = (TH2D*) gROOT->FindObject("anode");
792   /* Just for check
793   TCanvas* c = new TCanvas("Anode","Anode",800,600);
794   c->cd();
795   hist->Draw("lego1Fb"); // debug
796   c->Update();
797   Int_t tmp;
798   cin >> tmp;
799   */
800   Int_t nx = fHistAnode->GetNbinsX();
801   //Int_t ny = hist->GetNbinsY();
802   Int_t ic = localMax[iMax] / nx + 1;
803   Int_t jc = localMax[iMax] % nx + 1;
804
805   // Get min pad dimensions for the precluster
806   Int_t nSides = 2;
807   if (cluster.Multiplicity(0) == 0 || cluster.Multiplicity(1) == 0) nSides = 1;
808   TVector2 dim0 = cluster.MinPadDimensions(0, -1, kFALSE);
809   TVector2 dim1 = cluster.MinPadDimensions(1, -1, kFALSE);
810   //Double_t width[2][2] = {{dim0.X(), dim0.Y()},{dim1.X(),dim1.Y()}};
811   Int_t nonb[2] = {1, 0}; // coordinate index vs cathode
812   if (nSides == 1 || dim0.X() < dim1.X() - fgkDistancePrecision) {
813     nonb[0] = 0;
814     nonb[1] = 1;
815   }
816   Bool_t samex = kFALSE, samey = kFALSE;
817   if (TMath::Abs(dim0.X()-dim1.X()) < fgkDistancePrecision) samex = kTRUE; // the same X pad size on both planes 
818   if (TMath::Abs(dim0.Y()-dim1.Y()) < fgkDistancePrecision) samey = kTRUE; // the same Y pad size on both planes 
819
820   // Drop all pixels from the array - pick up only the ones from the cluster
821   //fPixArray->Delete();
822
823   Double_t wx = fHistAnode->GetXaxis()->GetBinWidth(1)/2; 
824   Double_t wy = fHistAnode->GetYaxis()->GetBinWidth(1)/2;  
825   Double_t yc = fHistAnode->GetYaxis()->GetBinCenter(ic);
826   Double_t xc = fHistAnode->GetXaxis()->GetBinCenter(jc);
827   Double_t cont = fHistAnode->GetCellContent(jc,ic);
828   AliMUONPad pixel(xc, yc, wx, wy, cont);
829   if (fDebug) pixel.Print("full"); 
830
831   Int_t npad = cluster.Multiplicity();
832   
833   // Pick up pads which overlap with the maximum pixel and find pads with the max signal
834   Double_t qMax[2] = {0}; 
835   AliMUONPad *matrix[2][3] = {{0x0,0x0,0x0},{0x0,0x0,0x0}};
836   for (Int_t j = 0; j < npad; ++j) 
837   {
838     AliMUONPad* pad = cluster.Pad(j);
839     if ( Overlap(*pad,pixel) )
840     {
841       if (fDebug) { cout << j << " "; pad->Print("full"); }
842       if (pad->Charge() > qMax[pad->Cathode()]) {
843         qMax[pad->Cathode()] = pad->Charge();
844         matrix[pad->Cathode()][1] = pad;
845         if (nSides == 1) matrix[!pad->Cathode()][1] = pad;
846       }
847     }
848   }
849   //if (nSides == 2 && (matrix[0][1] == 0x0 || matrix[1][1] == 0x0)) return; // ???
850
851   // Find neighbours of maxima to have 3 pads per direction (if possible)
852   for (Int_t j = 0; j < npad; ++j) 
853   {
854     AliMUONPad* pad = cluster.Pad(j);
855     Int_t cath = pad->Cathode();
856     if (pad == matrix[cath][1]) continue;
857     Int_t nLoops = 3 - nSides;
858
859     for (Int_t k = 0; k < nLoops; ++k) {
860       Int_t cath1 = cath;
861       if (k) cath1 = !cath;
862
863       // Check the coordinate corresponding to the cathode (bending or non-bending case)
864       Double_t dist = pad->Coord(nonb[cath1]) - matrix[cath][1]->Coord(nonb[cath1]);
865       Double_t dir = TMath::Sign (1., dist);
866       dist = TMath::Abs(dist) - pad->Size(nonb[cath1]) - matrix[cath][1]->Size(nonb[cath1]);
867
868       if (TMath::Abs(dist) < fgkDistancePrecision) {
869         // Check the other coordinate
870         dist = pad->Coord(!nonb[cath1]) - matrix[cath1][1]->Coord(!nonb[cath1]);
871         if (TMath::Abs(dist) > 
872             TMath::Max(pad->Size(!nonb[cath1]), matrix[cath1][1]->Size(!nonb[cath1])) - fgkDistancePrecision) break;
873         Int_t idir = TMath::Nint (dir);
874         if (matrix[cath1][1+idir] == 0x0) matrix[cath1][1+idir] = pad;
875         else if (pad->Charge() > matrix[cath1][1+idir]->Charge()) matrix[cath1][1+idir] = pad; // diff. segmentation
876         //cout << pad->Coord(nonb[cath1]) << " " << pad->Coord(!nonb[cath1]) << " " << pad->Size(nonb[cath1]) << " " << pad->Size(!nonb[cath1]) << " " << pad->Charge() << endl ;
877         break;
878       }
879     }
880   }
881
882   Double_t coord[2] = {0.}, qAver = 0.;
883   for (Int_t i = 0; i < 2; ++i) {
884     Double_t q = 0.;
885     Double_t coordQ = 0.;
886     Int_t cath = matrix[i][1]->Cathode();
887     if (i && nSides == 1) cath = !cath;
888     for (Int_t j = 0; j < 3; ++j) {
889       if (matrix[i][j] == 0x0) continue;
890       Double_t dq = matrix[i][j]->Charge();
891       q += dq;
892       coordQ += dq * matrix[i][j]->Coord(nonb[cath]);
893       //coordQ += (matrix[i][j]->Charge() * matrix[i][j]->Coord(nonb[cath]));
894     }
895     coord[cath] = coordQ / q;
896     qAver = TMath::Max (qAver, q);
897   }
898
899   //qAver = TMath::Sqrt(qAver);
900   if ( qAver >= 2.135 ) // JC: adc -> fc
901   {
902     
903     AliMUONCluster* cluster1 = new AliMUONCluster(cluster);
904       
905     cluster1->SetCharge(qAver,qAver);
906     if (nonb[0] == 1) 
907       cluster1->SetPosition(TVector2(coord[1],coord[0]),TVector2(0.,0.));
908     else 
909       cluster1->SetPosition(TVector2(coord[0],coord[1]),TVector2(0.,0.));
910
911     cluster1->SetChi2(0.);
912       
913     // FIXME: we miss some information in this cluster, as compared to 
914     // the original AddRawCluster code.
915       
916     AliDebug(2,Form("Adding RawCluster detElemId %4d mult %2d charge %5d (xl,yl)=(%9.6g,%9.6g)",
917                     fDetElemId,cluster1->Multiplicity(),(Int_t)cluster1->Charge(),
918                     cluster1->Position().X(),cluster1->Position().Y()));
919         
920     fClusterList.Add(cluster1);
921   }
922 }
923
924 //_____________________________________________________________________________
925 AliMUONClusterFinderPeakCOG&  
926 AliMUONClusterFinderPeakCOG::operator=(const AliMUONClusterFinderPeakCOG& rhs)
927 {
928 /// Protected assignement operator
929
930   if (this == &rhs) return *this;
931
932   AliFatal("Not implemented.");
933     
934   return *this;  
935 }    
936
937 //_____________________________________________________________________________
938 void AliMUONClusterFinderPeakCOG::PadsInXandY(AliMUONCluster& cluster,
939                                            Int_t &nInX, Int_t &nInY) const
940 {
941   /// Find number of pads in X and Y-directions (excluding virtual ones and
942   /// overflows)
943
944   //Int_t statusToTest = 1;
945   Int_t statusToTest = fgkUseForFit;
946   
947   //if ( nInX < 0 ) statusToTest = 0;
948   if ( nInX < 0 ) statusToTest = fgkZero;
949        
950   Bool_t mustMatch(kTRUE);
951
952   Long_t cn = cluster.NofPads(statusToTest,mustMatch);
953   
954   nInX = AliMp::PairFirst(cn);
955   nInY = AliMp::PairSecond(cn);
956 }
957
958 //_____________________________________________________________________________
959 void AliMUONClusterFinderPeakCOG::RemovePixel(Int_t i)
960 {
961   /// Remove pixel at index i
962   AliMUONPad* pixPtr = Pixel(i);
963   fPixArray->RemoveAt(i); 
964   delete pixPtr;
965 }
966
967 //_____________________________________________________________________________
968 AliMUONPad* 
969 AliMUONClusterFinderPeakCOG::Pixel(Int_t i) const
970 {
971   /// Returns pixel at index i
972   return static_cast<AliMUONPad*>(fPixArray->UncheckedAt(i));
973 }
974
975 //_____________________________________________________________________________
976 void 
977 AliMUONClusterFinderPeakCOG::Print(Option_t* what) const
978 {
979   /// printout
980   TString swhat(what);
981   swhat.ToLower();
982   if ( swhat.Contains("precluster") )
983   {
984     if ( fPreCluster) fPreCluster->Print();
985   }
986 }
987
988