1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
18 //-----------------------------------------------------------------------------
19 /// \class AliMUONClusterFinderMLEM
21 /// Clusterizer class based on the Expectation-Maximization algorithm
23 /// Pre-clustering is handled by AliMUONPreClusterFinder
24 /// From a precluster a pixel array is built, and from this pixel array
25 /// a list of clusters is output (using AliMUONClusterSplitterMLEM).
27 /// \author Laurent Aphecetche (for the "new" C++ structure) and
28 /// Alexander Zinchenko, JINR Dubna, for the hardcore of it ;-)
29 //-----------------------------------------------------------------------------
31 #include "AliMUONClusterFinderMLEM.h"
33 #include "AliMUONCluster.h"
34 #include "AliMUONClusterSplitterMLEM.h"
35 #include "AliMUONVDigit.h"
36 #include "AliMUONPad.h"
37 #include "AliMUONPreClusterFinder.h"
39 #include "AliMpVPadIterator.h"
40 #include "AliMpVSegmentation.h"
41 #include "AliRunLoader.h"
42 #include "AliMUONVDigitStore.h"
43 #include <Riostream.h>
49 #include "AliCodeTimer.h"
52 ClassImp(AliMUONClusterFinderMLEM)
55 const Double_t AliMUONClusterFinderMLEM::fgkZeroSuppression = 6; // average zero suppression value
56 //const Double_t AliMUONClusterFinderMLEM::fgkDistancePrecision = 1e-6; // (cm) used to check overlaps and so on
57 const Double_t AliMUONClusterFinderMLEM::fgkDistancePrecision = 1e-3; // (cm) used to check overlaps and so on
58 const TVector2 AliMUONClusterFinderMLEM::fgkIncreaseSize(-AliMUONClusterFinderMLEM::fgkDistancePrecision,-AliMUONClusterFinderMLEM::fgkDistancePrecision);
59 const TVector2 AliMUONClusterFinderMLEM::fgkDecreaseSize(AliMUONClusterFinderMLEM::fgkDistancePrecision,AliMUONClusterFinderMLEM::fgkDistancePrecision);
61 TMinuit* AliMUONClusterFinderMLEM::fgMinuit = 0x0;
62 // Status flags for pads
63 const Int_t AliMUONClusterFinderMLEM::fgkZero = 0x0; ///< pad "basic" state
64 const Int_t AliMUONClusterFinderMLEM::fgkMustKeep = 0x1; ///< do not kill (for pixels)
65 const Int_t AliMUONClusterFinderMLEM::fgkUseForFit = 0x10; ///< should be used for fit
66 const Int_t AliMUONClusterFinderMLEM::fgkOver = 0x100; ///< processing is over
67 const Int_t AliMUONClusterFinderMLEM::fgkModified = 0x1000; ///< modified pad charge
68 const Int_t AliMUONClusterFinderMLEM::fgkCoupled = 0x10000; ///< coupled pad
70 //_____________________________________________________________________________
71 AliMUONClusterFinderMLEM::AliMUONClusterFinderMLEM(Bool_t plot, AliMUONVClusterFinder* clusterFinder)
72 : AliMUONVClusterFinder(),
73 fPreClusterFinder(clusterFinder),
80 fPixArray(new TObjArray(20)),
89 fSegmentation[1] = fSegmentation[0] = 0x0;
91 fPadBeg[0] = fPadBeg[1] = fCathBeg;
93 if (!fgMinuit) fgMinuit = new TMinuit(8);
95 if (fPlot) fDebug = 1;
98 //_____________________________________________________________________________
99 AliMUONClusterFinderMLEM::~AliMUONClusterFinderMLEM()
102 delete fgMinuit; fgMinuit = 0; delete fPixArray; fPixArray = 0;
104 delete fPreClusterFinder;
106 AliInfo(Form("Total clusters %d AddVirtualPad needed %d",
107 fNClusters,fNAddVirtualPads));
110 //_____________________________________________________________________________
112 AliMUONClusterFinderMLEM::Prepare(Int_t detElemId,
113 TClonesArray* pads[2],
114 const AliMpArea& area,
115 const AliMpVSegmentation* seg[2])
117 /// Prepare for clustering
118 // AliCodeTimerAuto("")
120 for ( Int_t i = 0; i < 2; ++i )
122 fSegmentation[i] = seg[i];
125 // Find out the DetElemId
126 fDetElemId = detElemId;
129 fSplitter = new AliMUONClusterSplitterMLEM(fDetElemId,fPixArray);
130 fSplitter->SetDebug(fDebug);
132 // find out current event number, and reset the cluster number
133 fEventNumber = AliRunLoader::GetRunLoader()->GetEventNumber();
135 fClusterList.Delete();
137 AliDebug(3,Form("EVT %d DE %d",fEventNumber,fDetElemId));
139 if ( fPreClusterFinder->NeedSegmentation() )
141 return fPreClusterFinder->Prepare(detElemId,pads,area,seg);
145 return fPreClusterFinder->Prepare(detElemId,pads,area);
149 //_____________________________________________________________________________
151 AliMUONClusterFinderMLEM::NextCluster()
153 /// Return next cluster
154 // AliCodeTimerAuto("")
156 // if the list of clusters is not void, pick one from there
157 TObject* o = fClusterList.At(++fClusterNumber);
158 if ( o != 0x0 ) return static_cast<AliMUONCluster*>(o);
160 //FIXME : at this point, must check whether we've used all the digits
161 //from precluster : if not, let the preclustering know about those unused
162 //digits, so it can reuse them
164 // if the cluster list is exhausted, we need to go to the next
165 // pre-cluster and treat it
167 fPreCluster = fPreClusterFinder->NextCluster();
175 fClusterList.Delete(); // reset the list of clusters for this pre-cluster
176 fClusterNumber = -1; //AZ
180 // WorkOnPreCluster may have used only part of the pads, so we check that
181 // now, and let the unused pads be reused by the preclustering...
183 Int_t mult = fPreCluster->Multiplicity();
184 for ( Int_t i = 0; i < mult; ++i )
186 AliMUONPad* pad = fPreCluster->Pad(i);
187 if ( !pad->IsUsed() )
189 fPreClusterFinder->UsePad(*pad);
193 return NextCluster();
196 //_____________________________________________________________________________
198 AliMUONClusterFinderMLEM::WorkOnPreCluster()
200 /// Starting from a precluster, builds a pixel array, and then
201 /// extract clusters from this array
203 // AliCodeTimerAuto("")
206 cout << " *** Event # " << AliRunLoader::GetRunLoader()->GetEventNumber()
207 << " det. elem.: " << fDetElemId << endl;
208 for (Int_t j = 0; j < fPreCluster->Multiplicity(); ++j) {
209 AliMUONPad* pad = fPreCluster->Pad(j);
210 printf(" bbb %3d %1d %8.4f %8.4f %8.4f %8.4f %6.1f %3d %3d %2d %1d %1d \n",
211 j, pad->Cathode(), pad->Coord(0), pad->Coord(1), pad->DX()*2, pad->DY()*2,
212 pad->Charge(), pad->Ix(), pad->Iy(), pad->Status(), pad->IsReal(), pad->IsSaturated());
216 AliMUONCluster* cluster = CheckPrecluster(*fPreCluster);
217 if (!cluster) return kFALSE;
219 BuildPixArray(*cluster);
221 if ( fPixArray->GetLast() < 0 )
223 AliDebug(1,"No pixel for the above cluster");
228 // Use MLEM for cluster finder
229 Int_t nMax = 1, localMax[100], maxPos[100];
230 Double_t maxVal[100];
232 Int_t iSimple = 0, nInX = -1, nInY;
234 PadsInXandY(*cluster,nInX, nInY);
236 if (nInX < 4 && nInY < 4)
238 iSimple = 1; // simple cluster
242 nMax = FindLocalMaxima(fPixArray, localMax, maxVal); // for small clusters just to tag pixels
244 if (cluster->Multiplicity() <= 50) nMax = 1; // for small clusters
245 if (nMax > 1) TMath::Sort(nMax, maxVal, maxPos, kTRUE); // in descending order
249 for (Int_t i = 0; i < nMax; ++i)
253 FindCluster(*cluster,localMax, maxPos[i]);
256 MainLoop(*cluster,iSimple);
260 Int_t mult = cluster->Multiplicity();
261 for (Int_t j = 0; j < mult; ++j)
263 AliMUONPad* pad = cluster->Pad(j);
264 //if ( pad->Status() == 0 ) continue; // pad charge was not modified
265 if ( pad->Status() != fgkOver ) continue; // pad was not used
267 pad->SetStatus(fgkZero);
268 pad->RevertCharge(); // use backup charge value
271 } // for (Int_t i=0; i<nMax;
272 //if (nMax > 1) ((TH2D*) gROOT->FindObject("anode"))->Delete();
273 delete ((TH2D*) gROOT->FindObject("anode"));
274 //TH2D *mlem = (TH2D*) gROOT->FindObject("mlem");
275 //if (mlem) mlem->Delete();
280 //_____________________________________________________________________________
282 AliMUONClusterFinderMLEM::Overlap(const AliMUONPad& pad, const AliMUONPad& pixel)
284 /// Check if the pad and the pixel overlaps
286 // make a fake pad from the pixel
287 AliMUONPad tmp(pad.DetElemId(),pad.Cathode(),pad.Ix(),pad.Iy(),
288 pixel.Coord(0),pixel.Coord(1),
289 pixel.Size(0),pixel.Size(1),0);
291 return AliMUONPad::AreOverlapping(pad,tmp,fgkDecreaseSize);
294 //_____________________________________________________________________________
296 AliMUONClusterFinderMLEM::CheckPrecluster(const AliMUONCluster& origCluster)
298 /// Check precluster in order to attempt to simplify it (mostly for
299 /// two-cathode preclusters)
303 // Disregard small clusters (leftovers from splitting or noise)
304 if ((origCluster.Multiplicity()==1 || origCluster.Multiplicity()==2) &&
305 origCluster.Charge(0)+origCluster.Charge(1) < 10)
310 AliMUONCluster* cluster = new AliMUONCluster(origCluster);
312 AliDebug(2,"Start of CheckPreCluster=");
313 //StdoutToAliDebug(2,cluster->Print("full"));
315 AliMUONCluster* rv(0x0);
317 if (cluster->Multiplicity(0) && cluster->Multiplicity(1))
319 rv = CheckPreclusterTwoCathodes(cluster);
328 //_____________________________________________________________________________
330 AliMUONClusterFinderMLEM::CheckPreclusterTwoCathodes(AliMUONCluster* cluster)
332 /// Check two-cathode cluster
334 Int_t npad = cluster->Multiplicity();
335 Int_t* flags = new Int_t[npad];
336 for (Int_t j = 0; j < npad; ++j) flags[j] = 0;
338 // Check pad overlaps
339 for ( Int_t i = 0; i < npad; ++i)
341 AliMUONPad* padi = cluster->Pad(i);
342 if ( padi->Cathode() != 0 ) continue;
343 for (Int_t j = i+1; j < npad; ++j)
345 AliMUONPad* padj = cluster->Pad(j);
346 if ( padj->Cathode() != 1 ) continue;
347 if ( !AliMUONPad::AreOverlapping(*padi,*padj,fgkDecreaseSize) ) continue;
348 flags[i] = flags[j] = 1; // mark overlapped pads
352 // Check if all pads overlap
354 for (Int_t i = 0; i < npad; ++i)
356 if (!flags[i]) ++nFlags;
361 // not all pads overlap.
362 if (fDebug) cout << " nFlags: " << nFlags << endl;
363 TObjArray toBeRemoved;
364 for (Int_t i = 0; i < npad; ++i)
366 AliMUONPad* pad = cluster->Pad(i);
367 if (flags[i]) continue;
368 Int_t cath = pad->Cathode();
369 Int_t cath1 = TMath::Even(cath);
370 // Check for edge effect (missing pads on the _other_ cathode)
371 AliMpPad mpPad = fSegmentation[cath1]->PadByPosition(pad->Position(),kFALSE);
372 if (!mpPad.IsValid()) continue;
373 //if (nFlags == 1 && pad->Charge() < fgkZeroSuppression * 3) continue;
374 if (nFlags == 1 && pad->Charge() < 20) continue;
375 AliDebug(2,Form("Releasing the following pad : de,cath,ix,iy %d,%d,%d,%d charge %e",
376 fDetElemId,pad->Cathode(),pad->Ix(),pad->Iy(),pad->Charge()));
377 toBeRemoved.AddLast(pad);
378 fPreCluster->Pad(i)->Release();
380 Int_t nRemove = toBeRemoved.GetEntriesFast();
381 for ( Int_t i = 0; i < nRemove; ++i )
383 cluster->RemovePad(static_cast<AliMUONPad*>(toBeRemoved.UncheckedAt(i)));
387 // Check correlations of cathode charges
388 if ( !cluster->IsSaturated() && cluster->ChargeAsymmetry() > 1 )
391 Int_t cathode = cluster->MaxRawChargeCathode();
397 // get min and max pad charges on the cathode opposite to the
398 // max pad (given by MaxRawChargeCathode())
400 Int_t mult = cluster->Multiplicity();
401 for ( Int_t i = 0; i < mult; ++i )
403 AliMUONPad* pad = cluster->Pad(i);
404 if ( pad->Cathode() != cathode || !pad->IsReal() )
406 // only consider pads in the opposite cathode, and
407 // only consider real pads (i.e. exclude the virtual ones)
410 if ( pad->Charge() < cmin )
412 cmin = pad->Charge();
419 else if ( pad->Charge() > cmax )
421 cmax = pad->Charge();
425 AliDebug(2,Form("Pad imin,imax %d,%d cmin,cmax %e,%e",
426 imin,imax,cmin,cmax));
428 // arrange pads according to their distance to the max, normalized
430 Double_t* dist = new Double_t[mult];
435 AliMUONPad* padmax = cluster->Pad(imax);
437 for ( Int_t i = 0; i < mult; ++i )
440 if ( i == imax) continue;
441 AliMUONPad* pad = cluster->Pad(i);
442 if ( pad->Cathode() != cathode || !pad->IsReal() ) continue;
443 Double_t dx = (pad->X()-padmax->X())/padmax->DX()/2.0;
444 Double_t dy = (pad->Y()-padmax->Y())/padmax->DY()/2.0;
445 dist[i] = TMath::Sqrt(dx*dx+dy*dy);
448 dmin = dist[i] + 1E-3; // distance to the pad with minimum charge
454 TMath::Sort(mult,dist,flags,kFALSE); // in ascending order
455 Double_t xmax(-1), distPrev(999);
456 TObjArray toBeRemoved;
458 for ( Int_t i = 0; i < mult; ++i )
460 Int_t indx = flags[i];
461 AliMUONPad* pad = cluster->Pad(indx);
462 if ( pad->Cathode() != cathode || !pad->IsReal() ) continue;
463 if ( dist[indx] > dmin )
465 // farther than the minimum pad
466 Double_t dx = (pad->X()-padmax->X())/padmax->DX()/2.0;
467 Double_t dy = (pad->Y()-padmax->Y())/padmax->DY()/2.0;
470 if (dx >= 0 && dy >= 0) continue;
471 if (TMath::Abs(dx) > TMath::Abs(dy) && dx >= 0) continue;
472 if (TMath::Abs(dy) > TMath::Abs(dx) && dy >= 0) continue;
474 if (dist[indx] > distPrev + 1) break; // overstepping empty pads
475 if ( pad->Charge() <= cmax || TMath::Abs(dist[indx]-xmax) < 1E-3 )
478 if (TMath::Abs(dist[indx]-xmax) < 1.e-3)
480 cmax = TMath::Max(pad->Charge(),cmax);
484 cmax = pad->Charge();
487 distPrev = dist[indx];
488 AliDebug(2,Form("Releasing the following pad : de,cath,ix,iy %d,%d,%d,%d charge %e",
489 fDetElemId,pad->Cathode(),pad->Ix(),pad->Iy(),
492 toBeRemoved.AddLast(pad);
493 fPreCluster->Pad(indx)->Release();
496 Int_t nRemove = toBeRemoved.GetEntriesFast();
497 for ( Int_t i = 0; i < nRemove; ++i )
499 cluster->RemovePad(static_cast<AliMUONPad*>(toBeRemoved.UncheckedAt(i)));
506 AliDebug(2,"End of CheckPreClusterTwoCathodes=");
507 //StdoutToAliDebug(2,cluster->Print("full"));
512 //_____________________________________________________________________________
514 AliMUONClusterFinderMLEM::CheckOverlaps()
516 /// For debug only : check if some pixels overlap...
518 Int_t nPix = fPixArray->GetLast()+1;
521 for ( Int_t i = 0; i < nPix; ++i )
523 AliMUONPad* pixelI = Pixel(i);
524 AliMUONPad pi(dummy,dummy,dummy,dummy,
525 pixelI->Coord(0),pixelI->Coord(1),
526 pixelI->Size(0),pixelI->Size(1),0.0);
528 for ( Int_t j = i+1; j < nPix; ++j )
530 AliMUONPad* pixelJ = Pixel(j);
531 AliMUONPad pj(dummy,dummy,dummy,dummy,
532 pixelJ->Coord(0),pixelJ->Coord(1),
533 pixelJ->Size(0),pixelJ->Size(1),0.0);
536 if ( AliMUONPad::AreOverlapping(pi,pj,fgkDecreaseSize,area) )
538 AliInfo(Form("The following 2 pixels (%d and %d) overlap !",i,j));
540 StdoutToAliInfo(pixelI->Print();
541 cout << " Surface = " << pixelI->Size(0)*pixelI->Size(1)*4 << endl;
543 cout << " Surface = " << pixelJ->Size(0)*pixelJ->Size(1)*4 << endl;
544 cout << " Area surface = " << area.Dimensions().X()*area.Dimensions().Y()*4 << endl;
545 cout << "-------" << endl;
553 //_____________________________________________________________________________
554 void AliMUONClusterFinderMLEM::BuildPixArray(AliMUONCluster& cluster)
556 /// Build pixel array for MLEM method
558 Int_t npad = cluster.Multiplicity();
561 AliWarning("Got no pad at all ?!");
565 BuildPixArrayOneCathode(cluster);
567 Int_t nPix = fPixArray->GetLast()+1;
569 // AliDebug(2,Form("nPix after BuildPixArray=%d",nPix));
573 // AliDebug(2,Form("Will trim number of pixels to number of pads"));
575 // Too many pixels - sort and remove pixels with the lowest signal
577 for ( Int_t i = npad; i < nPix; ++i )
581 fPixArray->Compress();
582 } // if (nPix > npad)
584 // StdoutToAliDebug(2,cout << "End of BuildPixelArray:" << endl;
585 // fPixArray->Print(););
586 //CheckOverlaps();//FIXME : this is for debug only. Remove it.
589 //_____________________________________________________________________________
590 void AliMUONClusterFinderMLEM::BuildPixArrayOneCathode(AliMUONCluster& cluster)
592 /// Build the pixel array
594 // AliDebug(2,Form("cluster.Multiplicity=%d",cluster.Multiplicity()));
596 TVector2 dim = cluster.MinPadDimensions (-1, kFALSE);
597 Double_t width[2] = {dim.X(), dim.Y()}, xy0[2];
598 Int_t found[2] = {0}, mult = cluster.Multiplicity();
600 for ( Int_t i = 0; i < mult; ++i) {
601 AliMUONPad* pad = cluster.Pad(i);
602 for (Int_t j = 0; j < 2; ++j) {
603 if (found[j] == 0 && TMath::Abs(pad->Size(j)-width[j]) < fgkDistancePrecision) {
604 xy0[j] = pad->Coord(j);
608 if (found[0] && found[1]) break;
611 Double_t min[2], max[2];
612 Int_t cath0 = 0, cath1 = 1;
613 if (cluster.Multiplicity(0) == 0) cath0 = 1;
614 else if (cluster.Multiplicity(1) == 0) cath1 = 0;
616 TVector2 leftDown = cluster.Area(cath0).LeftDownCorner();
617 TVector2 rightUp = cluster.Area(cath0).RightUpCorner();
618 min[0] = leftDown.X();
619 min[1] = leftDown.Y();
620 max[0] = rightUp.X();
621 max[1] = rightUp.Y();
622 if (cath1 != cath0) {
623 leftDown = cluster.Area(cath1).LeftDownCorner();
624 rightUp = cluster.Area(cath1).RightUpCorner();
625 min[0] = TMath::Max (min[0], leftDown.X());
626 min[1] = TMath::Max (min[1], leftDown.Y());
627 max[0] = TMath::Min (max[0], rightUp.X());
628 max[1] = TMath::Min (max[1], rightUp.Y());
632 //width[0] /= 2; width[1] /= 2; // just for check
634 for (Int_t i = 0; i < 2; ++i) {
635 Double_t dist = (min[i] - xy0[i]) / width[i] / 2;
636 if (TMath::Abs(dist) < 1.e-6) dist = -1.e-6;
637 min[i] = xy0[i] + (TMath::Nint(dist-TMath::Sign(1.e-6,dist))
638 + TMath::Sign(0.5,dist)) * width[i] * 2;
639 nbins[i] = TMath::Nint ((max[i] - min[i]) / width[i] / 2);
640 if (nbins[i] == 0) ++nbins[i];
641 max[i] = min[i] + nbins[i] * width[i] * 2;
642 //cout << dist << " " << min[i] << " " << max[i] << " " << nbins[i] << endl;
646 TH2D *hist1 = new TH2D ("Grid", "", nbins[0], min[0], max[0], nbins[1], min[1], max[1]);
647 TH2D *hist2 = new TH2D ("Entries", "", nbins[0], min[0], max[0], nbins[1], min[1], max[1]);
648 TAxis *xaxis = hist1->GetXaxis();
649 TAxis *yaxis = hist1->GetYaxis();
652 for ( Int_t i = 0; i < mult; ++i) {
653 AliMUONPad* pad = cluster.Pad(i);
654 Int_t ix0 = xaxis->FindBin(pad->X());
655 Int_t iy0 = yaxis->FindBin(pad->Y());
656 PadOverHist(0, ix0, iy0, pad);
660 for (Int_t i = 1; i <= nbins[0]; ++i) {
661 Double_t x = xaxis->GetBinCenter(i);
662 for (Int_t j = 1; j <= nbins[1]; ++j) {
663 if (hist2->GetCellContent(i,j) < 0.1) continue;
664 //if (hist2->GetCellContent(i,j) < 1.1 && cluster.Multiplicity(0) &&
665 // cluster.Multiplicity(1)) continue;
666 if (cath0 != cath1) {
668 Double_t cont = hist2->GetCellContent(i,j);
669 if (cont < 999.) continue;
670 if (cont-Int_t(cont/1000.)*1000. < 0.5) continue;
672 Double_t y = yaxis->GetBinCenter(j);
673 Double_t charge = hist1->GetCellContent(i,j);
674 AliMUONPad* pixPtr = new AliMUONPad(x, y, width[0], width[1], charge);
675 fPixArray->Add(pixPtr);
679 if (fPixArray->GetEntriesFast() == 1) {
680 // Split pixel into 2
681 AliMUONPad* pixPtr = static_cast<AliMUONPad*> (fPixArray->UncheckedAt(0));
682 pixPtr->SetSize(0,width[0]/2.);
683 pixPtr->Shift(0,-width[0]/4.);
684 pixPtr = new AliMUONPad(pixPtr->X()+width[0], pixPtr->Y(), width[0]/2., width[1], pixPtr->Charge());
685 fPixArray->Add(pixPtr);
688 //fPixArray->Print();
693 //_____________________________________________________________________________
694 void AliMUONClusterFinderMLEM::PadOverHist(Int_t idir, Int_t ix0, Int_t iy0, AliMUONPad *pad)
696 /// "Span" pad over histogram in the direction idir
698 TH2D *hist1 = static_cast<TH2D*> (gROOT->FindObject("Grid"));
699 TH2D *hist2 = static_cast<TH2D*> (gROOT->FindObject("Entries"));
700 TAxis *axis = idir == 0 ? hist1->GetXaxis() : hist1->GetYaxis();
701 Int_t nbins = axis->GetNbins(), cath = pad->Cathode();
702 Double_t bin = axis->GetBinWidth(1), amask = TMath::Power(1000.,cath*1.);
704 Int_t nbinPad = (Int_t)(pad->Size(idir)/bin*2+fgkDistancePrecision) + 1; // number of bins covered by pad
706 for (Int_t i = 0; i < nbinPad; ++i) {
707 Int_t ixy = idir == 0 ? ix0 + i : iy0 + i;
708 if (ixy > nbins) break;
709 Double_t lowEdge = axis->GetBinLowEdge(ixy);
710 if (lowEdge + fgkDistancePrecision > pad->Coord(idir) + pad->Size(idir)) break;
711 if (idir == 0) PadOverHist(1, ixy, iy0, pad); // span in the other direction
714 Double_t cont = pad->Charge();
715 if (hist2->GetCellContent(ix0, ixy) > 0.1)
716 cont = TMath::Min (hist1->GetCellContent(ix0, ixy), cont);
717 hist1->SetCellContent(ix0, ixy, cont);
718 //hist2->SetCellContent(ix0, ixy, hist2->GetCellContent(ix0, ixy)+1);
719 hist2->SetCellContent(ix0, ixy, hist2->GetCellContent(ix0, ixy)+amask);
723 for (Int_t i = -1; i > -nbinPad; --i) {
724 Int_t ixy = idir == 0 ? ix0 + i : iy0 + i;
726 Double_t upEdge = axis->GetBinUpEdge(ixy);
727 if (upEdge - fgkDistancePrecision < pad->Coord(idir) - pad->Size(idir)) break;
728 if (idir == 0) PadOverHist(1, ixy, iy0, pad); // span in the other direction
731 Double_t cont = pad->Charge();
732 if (hist2->GetCellContent(ix0, ixy) > 0.1)
733 cont = TMath::Min (hist1->GetCellContent(ix0, ixy), cont);
734 hist1->SetCellContent(ix0, ixy, cont);
735 //hist2->SetCellContent(ix0, ixy, hist2->GetCellContent(ix0, ixy)+1);
736 hist2->SetCellContent(ix0, ixy, hist2->GetCellContent(ix0, ixy)+amask);
741 //_____________________________________________________________________________
743 AliMUONClusterFinderMLEM::Plot(const char* basename)
745 /// Make a plot and save it as png
750 TCanvas* c = new TCanvas("MLEM","MLEM",800,600);
755 c->Print(Form("%s.EVT%d.DE%d.CLU%d.png",basename,fEventNumber,
756 fDetElemId,fClusterNumber));
759 //_____________________________________________________________________________
761 AliMUONClusterFinderMLEM::ComputeCoefficients(AliMUONCluster& cluster,
765 /// Compute coefficients needed for MLEM algorithm
767 Int_t nPix = fPixArray->GetLast()+1;
769 //memset(probi,0,nPix*sizeof(Double_t));
770 for (Int_t j = 0; j < nPix; ++j) probi[j] = 0.;
772 Int_t mult = cluster.Multiplicity();
773 for ( Int_t j = 0; j < mult; ++j )
775 AliMUONPad* pad = cluster.Pad(j);
778 for ( Int_t ipix = 0; ipix < nPix; ++ipix )
780 Int_t indx1 = indx + ipix;
781 //if (pad->Status() < 0)
782 if (pad->Status() != fgkZero)
787 AliMUONPad* pixPtr = Pixel(ipix);
788 // coef is the charge (given by Mathieson integral) on pad, assuming
789 // the Mathieson is center at pixel.
790 coef[indx1] = fSplitter->ChargeIntegration(pixPtr->Coord(0), pixPtr->Coord(1), *pad);
791 // AliDebug(2,Form("pad=(%d,%d,%e,%e,%e,%e) pix=(%e,%e,%e,%e) coef %e",
792 // pad->Ix(),pad->Iy(),
793 // pad->X(),pad->Y(),
794 // pad->DX(),pad->DY(),
795 // pixPtr->Coord(0),pixPtr->Coord(1),
796 // pixPtr->Size(0),pixPtr->Size(1),
799 probi[ipix] += coef[indx1];
804 //_____________________________________________________________________________
805 Bool_t AliMUONClusterFinderMLEM::MainLoop(AliMUONCluster& cluster, Int_t iSimple)
807 /// Repeat MLEM algorithm until pixel size becomes sufficiently small
809 // AliCodeTimerAuto("")
811 Int_t nPix = fPixArray->GetLast()+1;
813 AliDebug(2,Form("nPix=%d iSimple=%d, precluster=",nPix,iSimple));
814 //StdoutToAliDebug(2,cluster.Print("full"););
818 AliDebug(1,"No pixels, why am I here then ?");
821 AddVirtualPad(cluster); // add virtual pads if necessary
823 Int_t npadTot = cluster.Multiplicity();
825 for (Int_t i = 0; i < npadTot; ++i)
827 //if (cluster.Pad(i)->Status() == 0) ++npadOK;
828 if (cluster.Pad(i)->Status() == fgkZero) ++npadOK;
833 Double_t* probi(0x0);
834 Int_t lc(0); // loop counter
836 //Plot("mlem.start");
837 AliMUONPad* pixPtr = Pixel(0);
838 Double_t xylim[4] = {pixPtr->X(), -pixPtr->X(), pixPtr->Y(), -pixPtr->Y()};
843 mlem = (TH2D*) gROOT->FindObject("mlem");
846 AliDebug(2,Form("lc %d nPix %d(%d) npadTot %d npadOK %d",lc,nPix,fPixArray->GetLast()+1,npadTot,npadOK));
847 AliDebug(2,Form("EVT%d PixArray=",fEventNumber));
848 //StdoutToAliDebug(2,fPixArray->Print("","full"));
850 coef = new Double_t [npadTot*nPix];
851 probi = new Double_t [nPix];
853 // Calculate coefficients and pixel visibilities
854 ComputeCoefficients(cluster,coef,probi);
856 for (Int_t ipix = 0; ipix < nPix; ++ipix)
858 if (probi[ipix] < 0.01)
860 AliMUONPad* pixel = Pixel(ipix);
861 AliDebug(2,Form("Setting the following pixel to invisible as its probi<0.01:"));
862 //StdoutToAliDebug(2,cout << Form(" -- ipix %3d --- "); pixel->Print(););
863 pixel->SetCharge(0); // "invisible" pixel
868 Mlem(cluster,coef, probi, 15);
870 // Find histogram limits for the 1'st pass only - for others computed below
872 for ( Int_t ipix = 1; ipix < nPix; ++ipix )
874 pixPtr = Pixel(ipix);
875 for ( Int_t i = 0; i < 2; ++i )
878 if (pixPtr->Coord(i) < xylim[indx]) xylim[indx] = pixPtr->Coord(i);
879 else if (-pixPtr->Coord(i) < xylim[indx+1]) xylim[indx+1] = -pixPtr->Coord(i);
882 } else pixPtr = Pixel(0);
884 for (Int_t i = 0; i < 4; i++)
886 xylim[i] -= pixPtr->Size(i/2);
889 Int_t nx = TMath::Nint ((-xylim[1]-xylim[0])/pixPtr->Size(0)/2);
890 Int_t ny = TMath::Nint ((-xylim[3]-xylim[2])/pixPtr->Size(1)/2);
892 //StdoutToAliDebug(2,cout << "pixel used for nx,ny computation : "; pixPtr->Print(););
893 AliDebug(2,Form("lc %d pixPtr size = %e,%e nx,ny=%d,%d xylim=%e,%e,%e,%e",
894 lc,pixPtr->Size(0),pixPtr->Size(1),nx,ny,
895 xylim[0],-xylim[1],xylim[2],-xylim[3]
898 mlem = new TH2D("mlem","mlem",nx,xylim[0],-xylim[1],ny,xylim[2],-xylim[3]);
900 for (Int_t ipix = 0; ipix < nPix; ++ipix)
902 AliMUONPad* pixPtr = Pixel(ipix);
903 mlem->Fill(pixPtr->Coord(0),pixPtr->Coord(1),pixPtr->Charge());
906 // Check if the total charge of pixels is too low
908 for ( Int_t i = 0; i < nPix; ++i)
910 qTot += Pixel(i)->Charge();
913 if ( qTot < 1.e-4 || ( npadOK < 3 && qTot < 7 ) )
915 AliDebug(1,Form("Deleting the above cluster (charge %e too low, npadOK=%d)",qTot,npadOK));
919 for ( Int_t i = 0; i < npadTot; ++i)
921 AliMUONPad* pad = cluster.Pad(i);
922 //if ( pad->Status() == 0) pad->SetStatus(-1);
923 if ( pad->Status() == fgkZero) pad->SetStatus(fgkOver);
930 // Simple cluster - skip further passes thru EM-procedure
938 // Calculate position of the center-of-gravity around the maximum pixel
940 FindCOG(mlem, xyCOG);
942 if (TMath::Min(pixPtr->Size(0),pixPtr->Size(1)) < 0.07 &&
943 pixPtr->Size(0) > pixPtr->Size(1)) break;
945 // Sort pixels according to the charge
946 MaskPeaks(1); // mask local maxima
948 MaskPeaks(0); // unmask local maxima
949 Double_t pixMin = 0.01*Pixel(0)->Charge();
950 pixMin = TMath::Min(pixMin,50.);
952 // Decrease pixel size and shift pixels to make them centered at
954 Int_t indx = (pixPtr->Size(0)>pixPtr->Size(1)) ? 0 : 1;
957 Double_t shift[2] = { 0.0, 0.0 };
958 for (Int_t i = 0; i < 4; ++i) xylim[i] = 999;
961 for (Int_t ipix = 0; ipix < nPix1; ++ipix)
963 AliMUONPad* pixPtr = Pixel(ipix);
964 if ( nPix >= npadOK // too many pixels already
966 pixPtr->Charge() < pixMin && pixPtr->Status() != fgkMustKeep // too low charge
972 for (Int_t i = 0; i < 2; ++i)
976 pixPtr->SetCharge(10);
977 pixPtr->SetSize(indx, pixPtr->Size(indx)/2);
978 width = -pixPtr->Size(indx);
979 pixPtr->Shift(indx, width);
980 // Shift pixel position
984 for (Int_t j = 0; j < 2; ++j)
986 shift[j] = pixPtr->Coord(j) - xyCOG[j];
987 shift[j] -= ((Int_t)(shift[j]/pixPtr->Size(j)/2))*pixPtr->Size(j)*2;
990 pixPtr->Shift(0, -shift[0]);
991 pixPtr->Shift(1, -shift[1]);
994 else if (nPix < npadOK)
996 pixPtr = new AliMUONPad(*pixPtr);
997 pixPtr->Shift(indx, -2*width);
998 pixPtr->SetStatus(fgkZero);
999 fPixArray->Add(pixPtr);
1002 else continue; // skip adjustment of histo limits
1003 for (Int_t j = 0; j < 4; ++j)
1005 xylim[j] = TMath::Min (xylim[j], (j%2 ? -1 : 1)*pixPtr->Coord(j/2));
1007 } // for (Int_t i=0; i<2;
1008 } // for (Int_t ipix=0;
1010 fPixArray->Compress();
1012 AliDebug(2,Form("After shift:"));
1013 //StdoutToAliDebug(2,fPixArray->Print("","full"););
1014 //Plot(Form("mlem.lc%d",lc+1));
1016 AliDebug(2,Form(" xyCOG=%9.6f %9.6f xylim=%9.6f,%9.6f,%9.6f,%9.6f",
1019 xylim[2],xylim[3]));
1023 AliMUONPad* pixPtr = Pixel(0);
1024 // add pixels if the maximum is at the limit of pixel area:
1025 // start from Y-direction
1027 for (Int_t i = 3; i > -1; --i)
1029 if (nPix < npadOK &&
1030 TMath::Abs((i%2 ? -1 : 1)*xylim[i]-xyCOG[i/2]) < pixPtr->Size(i/2))
1032 //AliMUONPad* p = static_cast<AliMUONPad*>(pixPtr->Clone());
1033 AliMUONPad* p = new AliMUONPad(*pixPtr);
1034 p->SetCoord(i/2, xyCOG[i/2]+(i%2 ? 2:-2)*pixPtr->Size(i/2));
1035 xylim[i] = p->Coord(i/2) * (i%2 ? -1 : 1); // update histo limits
1036 j = TMath::Even (i/2);
1037 p->SetCoord(j, xyCOG[j]);
1038 AliDebug(2,Form("Adding pixel on the edge (i=%d) ",i));
1039 //StdoutToAliDebug(2,cout << " ---- ";
1040 // p->Print("corners"););
1046 nPix = fPixArray->GetEntriesFast();
1051 AliDebug(2,Form("At the end of while loop nPix=%d : ",fPixArray->GetLast()+1));
1052 //StdoutToAliDebug(2,fPixArray->Print("","full"););
1054 // remove pixels with low signal or low visibility
1055 // Cuts are empirical !!!
1056 Double_t thresh = TMath::Max (mlem->GetMaximum()/100.,1.);
1057 thresh = TMath::Min (thresh,50.);
1058 Double_t charge = 0;
1060 // Mark pixels which should be removed
1061 for (Int_t i = 0; i < nPix; ++i)
1063 AliMUONPad* pixPtr = Pixel(i);
1064 charge = pixPtr->Charge();
1065 if (charge < thresh)
1067 pixPtr->SetCharge(-charge);
1071 // Move charge of removed pixels to their nearest neighbour (to keep total charge the same)
1073 for (Int_t i = 0; i < nPix; ++i)
1075 AliMUONPad* pixPtr = Pixel(i);
1076 charge = pixPtr->Charge();
1077 if (charge > 0) continue;
1078 near = FindNearest(pixPtr);
1079 pixPtr->SetCharge(0);
1080 probi[i] = 0; // make it "invisible"
1081 AliMUONPad* pnear = Pixel(near);
1082 pnear->SetCharge(pnear->Charge() + (-charge));
1084 Mlem(cluster,coef,probi,2);
1086 AliDebug(2,Form("Before splitting nPix=%d EVT %d DE %d",fPixArray->GetLast()+1,fEventNumber,fDetElemId));
1087 //StdoutToAliDebug(2,fPixArray->Print("","full"););
1088 //Plot("mlem.beforesplit");
1091 for (Int_t i = 0; i < nPix; ++i)
1093 AliMUONPad* pixPtr = Pixel(i);
1094 Int_t ix = mlem->GetXaxis()->FindBin(pixPtr->Coord(0));
1095 Int_t iy = mlem->GetYaxis()->FindBin(pixPtr->Coord(1));
1096 mlem->SetBinContent(ix, iy, pixPtr->Charge());
1099 // Try to split into clusters
1101 if (mlem->GetSum() < 1)
1107 fSplitter->Split(cluster,mlem,coef,fClusterList);
1112 fPixArray->Delete();
1117 //_____________________________________________________________________________
1118 void AliMUONClusterFinderMLEM::MaskPeaks(Int_t mask)
1120 /// Mask/unmask pixels corresponding to local maxima (add/subtract 10000 to their charge
1121 /// - to avoid loosing low charge pixels after sorting)
1123 for (Int_t i = 0; i < fPixArray->GetEntriesFast(); ++i) {
1124 AliMUONPad* pix = Pixel(i);
1125 if (pix->Status() == fgkMustKeep) {
1126 if (mask == 1) pix->SetCharge(pix->Charge()+10000.);
1127 else pix->SetCharge(pix->Charge()-10000.);
1132 //_____________________________________________________________________________
1133 void AliMUONClusterFinderMLEM::Mlem(AliMUONCluster& cluster,
1134 Double_t* coef, Double_t* probi,
1137 /// Use MLEM to find pixel charges
1139 Int_t nPix = fPixArray->GetEntriesFast();
1141 Int_t npad = cluster.Multiplicity();
1143 Double_t* probi1 = new Double_t[nPix];
1144 Double_t probMax = TMath::MaxElement(nPix,probi);
1146 for (Int_t iter = 0; iter < nIter; ++iter)
1149 for (Int_t ipix = 0; ipix < nPix; ++ipix)
1151 Pixel(ipix)->SetChargeBackup(0);
1152 // Correct each pixel
1153 if (probi[ipix] < 0.01) continue; // skip "invisible" pixel
1155 probi1[ipix] = probMax;
1156 for (Int_t j = 0; j < npad; ++j)
1158 AliMUONPad* pad = cluster.Pad(j);
1159 //if (pad->Status() < 0) continue;
1160 if (pad->Status() != fgkZero) continue;
1162 Int_t indx1 = j*nPix;
1163 Int_t indx = indx1 + ipix;
1164 // Calculate expectation
1165 for (Int_t i = 0; i < nPix; ++i)
1167 sum1 += Pixel(i)->Charge()*coef[indx1+i];
1168 //cout << i << " " << Pixel(i)->Charge() << " " << coef[indx1+i] << endl;
1170 if ( pad->IsSaturated() && sum1 > pad->Charge() )
1172 // correct for pad charge overflows
1173 probi1[ipix] -= coef[indx];
1175 //sum1 = pad->Charge();
1178 if (sum1 > 1.e-6) sum += pad->Charge()*coef[indx]/sum1;
1179 //if (coef[indx] > 1.e-6) sum += pad->Charge()*coef[indx]/sum1;
1180 } // for (Int_t j=0;
1181 AliMUONPad* pixPtr = Pixel(ipix);
1182 if (probi1[ipix] > 1.e-6)
1184 //AZ pixPtr->SetCharge(pixPtr->Charge()*sum/probi1[ipix]);
1185 pixPtr->SetChargeBackup(pixPtr->Charge()*sum/probi1[ipix]);
1187 //cout << " xxx " << ipix << " " << pixPtr->Charge() << " " << pixPtr->ChargeBackup() << " " << sum << " " << probi1[ipix] << endl;
1188 } // for (Int_t ipix=0;
1190 for (Int_t i = 0; i < nPix; ++i) {
1191 AliMUONPad* pixPtr = Pixel(i);
1192 pixPtr->RevertCharge();
1193 qTot += pixPtr->Charge();
1196 // Can happen in clusters with large number of overflows - speeding up
1200 } // for (Int_t iter=0;
1204 //_____________________________________________________________________________
1205 void AliMUONClusterFinderMLEM::FindCOG(TH2D *mlem, Double_t *xyc)
1207 /// Calculate position of the center-of-gravity around the maximum pixel
1209 Int_t ixmax, iymax, ix, nsumx=0, nsumy=0, nsum=0;
1210 Int_t i1 = -9, j1 = -9;
1211 mlem->GetMaximumBin(ixmax,iymax,ix);
1212 Int_t nx = mlem->GetNbinsX();
1213 Int_t ny = mlem->GetNbinsY();
1214 Double_t thresh = mlem->GetMaximum()/10;
1215 Double_t x, y, cont, xq=0, yq=0, qq=0;
1217 Int_t ie = TMath::Min(ny,iymax+1), je = TMath::Min(nx,ixmax+1);
1218 for (Int_t i = TMath::Max(1,iymax-1); i <= ie; ++i) {
1219 y = mlem->GetYaxis()->GetBinCenter(i);
1220 for (Int_t j = TMath::Max(1,ixmax-1); j <= je; ++j) {
1221 cont = mlem->GetCellContent(j,i);
1222 if (cont < thresh) continue;
1223 if (i != i1) {i1 = i; nsumy++;}
1224 if (j != j1) {j1 = j; nsumx++;}
1225 x = mlem->GetXaxis()->GetBinCenter(j);
1234 Int_t i2 = 0, j2 = 0;
1237 // one bin in Y - add one more (with the largest signal)
1238 for (Int_t i = TMath::Max(1,iymax-1); i <= ie; ++i) {
1239 if (i == iymax) continue;
1240 for (Int_t j = TMath::Max(1,ixmax-1); j <= je; ++j) {
1241 cont = mlem->GetCellContent(j,i);
1244 x = mlem->GetXaxis()->GetBinCenter(j);
1245 y = mlem->GetYaxis()->GetBinCenter(i);
1254 if (i2 != i1) nsumy++;
1255 if (j2 != j1) nsumx++;
1257 } // if (nsumy == 1)
1260 // one bin in X - add one more (with the largest signal)
1262 for (Int_t j = TMath::Max(1,ixmax-1); j <= je; ++j) {
1263 if (j == ixmax) continue;
1264 for (Int_t i = TMath::Max(1,iymax-1); i <= ie; ++i) {
1265 cont = mlem->GetCellContent(j,i);
1268 x = mlem->GetXaxis()->GetBinCenter(j);
1269 y = mlem->GetYaxis()->GetBinCenter(i);
1278 if (i2 != i1) nsumy++;
1279 if (j2 != j1) nsumx++;
1281 } // if (nsumx == 1)
1283 xyc[0] = xq/qq; xyc[1] = yq/qq;
1284 AliDebug(2,Form("x,y COG = %e,%e",xyc[0],xyc[1]));
1287 //_____________________________________________________________________________
1288 Int_t AliMUONClusterFinderMLEM::FindNearest(AliMUONPad *pixPtr0)
1290 /// Find the pixel nearest to the given one
1291 /// (algorithm may be not very efficient)
1293 Int_t nPix = fPixArray->GetEntriesFast(), imin = 0;
1294 Double_t rmin = 99999, dx = 0, dy = 0, r = 0;
1295 Double_t xc = pixPtr0->Coord(0), yc = pixPtr0->Coord(1);
1298 for (Int_t i = 0; i < nPix; ++i) {
1299 pixPtr = (AliMUONPad*) fPixArray->UncheckedAt(i);
1300 if (pixPtr == pixPtr0 || pixPtr->Charge() < 0.5) continue;
1301 dx = (xc - pixPtr->Coord(0)) / pixPtr->Size(0);
1302 dy = (yc - pixPtr->Coord(1)) / pixPtr->Size(1);
1303 r = dx *dx + dy * dy;
1304 if (r < rmin) { rmin = r; imin = i; }
1309 //_____________________________________________________________________________
1311 AliMUONClusterFinderMLEM::Paint(Option_t*)
1313 /// Paint cluster and pixels
1315 AliMpArea area(fPreCluster->Area());
1317 gPad->Range(area.LeftBorder(),area.DownBorder(),area.RightBorder(),area.UpBorder());
1319 gVirtualX->SetFillStyle(1001);
1320 gVirtualX->SetFillColor(3);
1321 gVirtualX->SetLineColor(3);
1325 for ( Int_t i = 0; i <= fPixArray->GetLast(); ++i)
1327 AliMUONPad* pixel = Pixel(i);
1329 gPad->PaintBox(pixel->Coord(0)-pixel->Size(0)*s,
1330 pixel->Coord(1)-pixel->Size(1)*s,
1331 pixel->Coord(0)+pixel->Size(0)*s,
1332 pixel->Coord(1)+pixel->Size(1)*s);
1334 // for ( Int_t sign = -1; sign < 2; sign +=2 )
1336 // gPad->PaintLine(pixel->Coord(0) - pixel->Size(0),
1337 // pixel->Coord(1) + sign*pixel->Size(1),
1338 // pixel->Coord(0) + pixel->Size(0),
1339 // pixel->Coord(1) - sign*pixel->Size(1)
1345 gVirtualX->SetFillStyle(0);
1347 fPreCluster->Paint();
1349 gVirtualX->SetLineColor(1);
1350 gVirtualX->SetLineWidth(2);
1351 gVirtualX->SetFillStyle(0);
1352 gVirtualX->SetTextColor(1);
1353 gVirtualX->SetTextAlign(22);
1355 for ( Int_t i = 0; i <= fPixArray->GetLast(); ++i)
1357 AliMUONPad* pixel = Pixel(i);
1358 gPad->PaintBox(pixel->Coord(0)-pixel->Size(0),
1359 pixel->Coord(1)-pixel->Size(1),
1360 pixel->Coord(0)+pixel->Size(0),
1361 pixel->Coord(1)+pixel->Size(1));
1362 gVirtualX->SetTextSize(pixel->Size(0)*60);
1364 gPad->PaintText(pixel->Coord(0),pixel->Coord(1),Form("%d",(Int_t)(pixel->Charge())));
1368 //_____________________________________________________________________________
1369 Int_t AliMUONClusterFinderMLEM::FindLocalMaxima(TObjArray *pixArray, Int_t *localMax, Double_t *maxVal)
1371 /// Find local maxima in pixel space for large preclusters in order to
1372 /// try to split them into smaller pieces (to speed up the MLEM procedure)
1373 /// or to find additional fitting seeds if clusters were not completely resolved
1375 AliDebug(1,Form("nPix=%d",pixArray->GetLast()+1));
1378 //if (pixArray == fPixArray) hist = (TH2D*) gROOT->FindObject("anode");
1379 //else { hist = (TH2D*) gROOT->FindObject("anode1"); cout << hist << endl; }
1380 //if (hist) hist->Delete();
1382 Double_t xylim[4] = {999, 999, 999, 999};
1384 Int_t nPix = pixArray->GetEntriesFast();
1385 AliMUONPad *pixPtr = 0;
1386 for (Int_t ipix = 0; ipix < nPix; ++ipix) {
1387 pixPtr = (AliMUONPad*) pixArray->UncheckedAt(ipix);
1388 for (Int_t i = 0; i < 4; ++i)
1389 xylim[i] = TMath::Min (xylim[i], (i%2 ? -1 : 1)*pixPtr->Coord(i/2));
1391 for (Int_t i = 0; i < 4; ++i) xylim[i] -= pixPtr->Size(i/2);
1393 Int_t nx = TMath::Nint ((-xylim[1]-xylim[0])/pixPtr->Size(0)/2);
1394 Int_t ny = TMath::Nint ((-xylim[3]-xylim[2])/pixPtr->Size(1)/2);
1395 if (pixArray == fPixArray) hist = new TH2D("anode","anode",nx,xylim[0],-xylim[1],ny,xylim[2],-xylim[3]);
1396 else hist = new TH2D("anode1","anode1",nx,xylim[0],-xylim[1],ny,xylim[2],-xylim[3]);
1397 for (Int_t ipix = 0; ipix < nPix; ++ipix) {
1398 pixPtr = (AliMUONPad*) pixArray->UncheckedAt(ipix);
1399 hist->Fill(pixPtr->Coord(0), pixPtr->Coord(1), pixPtr->Charge());
1401 // if (fDraw && pixArray == fPixArray) fDraw->DrawHist("c2", hist);
1403 Int_t nMax = 0, indx, nxy = ny * nx;
1404 Int_t *isLocalMax = new Int_t[nxy];
1405 for (Int_t i = 0; i < nxy; ++i) isLocalMax[i] = 0;
1407 for (Int_t i = 1; i <= ny; ++i) {
1409 for (Int_t j = 1; j <= nx; ++j) {
1410 if (hist->GetCellContent(j,i) < 0.5) continue;
1411 //if (isLocalMax[indx+j-1] < 0) continue;
1412 if (isLocalMax[indx+j-1] != 0) continue;
1413 FlagLocalMax(hist, i, j, isLocalMax);
1417 for (Int_t i = 1; i <= ny; ++i) {
1419 for (Int_t j = 1; j <= nx; ++j) {
1420 if (isLocalMax[indx+j-1] > 0) {
1421 localMax[nMax] = indx + j - 1;
1422 maxVal[nMax++] = hist->GetCellContent(j,i);
1423 ((AliMUONPad*)fSplitter->BinToPix(hist, j, i))->SetStatus(fgkMustKeep);
1424 if (nMax > 99) AliFatal(" Too many local maxima !!!");
1428 if (fDebug) cout << " Local max: " << nMax << endl;
1429 delete [] isLocalMax;
1433 //_____________________________________________________________________________
1434 void AliMUONClusterFinderMLEM::FlagLocalMax(TH2D *hist, Int_t i, Int_t j, Int_t *isLocalMax)
1436 /// Flag pixels (whether or not local maxima)
1438 Int_t nx = hist->GetNbinsX();
1439 Int_t ny = hist->GetNbinsY();
1440 Int_t cont = TMath::Nint (hist->GetCellContent(j,i));
1441 Int_t cont1 = 0, indx = (i-1)*nx+j-1, indx1 = 0, indx2 = 0;
1443 Int_t ie = i + 2, je = j + 2;
1444 for (Int_t i1 = i-1; i1 < ie; ++i1) {
1445 if (i1 < 1 || i1 > ny) continue;
1446 indx1 = (i1 - 1) * nx;
1447 for (Int_t j1 = j-1; j1 < je; ++j1) {
1448 if (j1 < 1 || j1 > nx) continue;
1449 if (i == i1 && j == j1) continue;
1450 indx2 = indx1 + j1 - 1;
1451 cont1 = TMath::Nint (hist->GetCellContent(j1,i1));
1452 if (cont < cont1) { isLocalMax[indx] = -1; return; }
1453 else if (cont > cont1) isLocalMax[indx2] = -1;
1454 else { // the same charge
1455 isLocalMax[indx] = 1;
1456 if (isLocalMax[indx2] == 0) {
1457 FlagLocalMax(hist, i1, j1, isLocalMax);
1458 if (isLocalMax[indx2] < 0) { isLocalMax[indx] = -1; return; }
1459 else isLocalMax[indx2] = -1;
1464 isLocalMax[indx] = 1; // local maximum
1467 //_____________________________________________________________________________
1468 void AliMUONClusterFinderMLEM::FindCluster(AliMUONCluster& cluster,
1469 Int_t *localMax, Int_t iMax)
1471 /// Find pixel cluster around local maximum \a iMax and pick up pads
1472 /// overlapping with it
1474 TH2D *hist = (TH2D*) gROOT->FindObject("anode");
1476 TCanvas* c = new TCanvas("Anode","Anode",800,600);
1478 hist->Draw("lego1Fb"); // debug
1483 Int_t nx = hist->GetNbinsX();
1484 Int_t ny = hist->GetNbinsY();
1485 Int_t ic = localMax[iMax] / nx + 1;
1486 Int_t jc = localMax[iMax] % nx + 1;
1487 Int_t nxy = ny * nx;
1488 Bool_t *used = new Bool_t[nxy];
1489 for (Int_t i = 0; i < nxy; ++i) used[i] = kFALSE;
1491 // Drop all pixels from the array - pick up only the ones from the cluster
1492 fPixArray->Delete();
1494 Double_t wx = hist->GetXaxis()->GetBinWidth(1)/2;
1495 Double_t wy = hist->GetYaxis()->GetBinWidth(1)/2;
1496 Double_t yc = hist->GetYaxis()->GetBinCenter(ic);
1497 Double_t xc = hist->GetXaxis()->GetBinCenter(jc);
1498 Double_t cont = hist->GetCellContent(jc,ic);
1499 fPixArray->Add(new AliMUONPad (xc, yc, wx, wy, cont));
1500 used[(ic-1)*nx+jc-1] = kTRUE;
1501 AddBinSimple(hist, ic, jc);
1502 //fSplitter->AddBin(hist, ic, jc, 1, used, (TObjArray*)0); // recursive call
1504 Int_t nPix = fPixArray->GetEntriesFast();
1505 Int_t npad = cluster.Multiplicity();
1507 for (Int_t i = 0; i < nPix; ++i)
1509 AliMUONPad* pixPtr = Pixel(i);
1510 pixPtr->SetSize(0,wx);
1511 pixPtr->SetSize(1,wy);
1514 // Pick up pads which overlap with found pixels
1515 for (Int_t i = 0; i < npad; ++i)
1517 //cluster.Pad(i)->SetStatus(-1);
1518 cluster.Pad(i)->SetStatus(fgkOver); // just the dirty trick
1521 for (Int_t i = 0; i < nPix; ++i)
1523 AliMUONPad* pixPtr = Pixel(i);
1524 for (Int_t j = 0; j < npad; ++j)
1526 AliMUONPad* pad = cluster.Pad(j);
1527 //if (pad->Status() == 0) continue;
1528 if (pad->Status() == fgkZero) continue;
1529 if ( Overlap(*pad,*pixPtr) )
1531 //pad->SetStatus(0);
1532 pad->SetStatus(fgkZero);
1533 if (fDebug) { cout << j << " "; pad->Print("full"); }
1541 //_____________________________________________________________________________
1543 AliMUONClusterFinderMLEM::AddBinSimple(TH2D *mlem, Int_t ic, Int_t jc)
1545 /// Add adjacent bins (+-1 in X and Y) to the cluster
1547 Int_t nx = mlem->GetNbinsX();
1548 Int_t ny = mlem->GetNbinsY();
1549 Double_t cont1, cont = mlem->GetCellContent(jc,ic);
1550 AliMUONPad *pixPtr = 0;
1552 Int_t ie = TMath::Min(ic+1,ny), je = TMath::Min(jc+1,nx);
1553 for (Int_t i = TMath::Max(ic-1,1); i <= ie; ++i) {
1554 for (Int_t j = TMath::Max(jc-1,1); j <= je; ++j) {
1555 cont1 = mlem->GetCellContent(j,i);
1556 if (cont1 > cont) continue;
1557 if (cont1 < 0.5) continue;
1558 pixPtr = new AliMUONPad (mlem->GetXaxis()->GetBinCenter(j),
1559 mlem->GetYaxis()->GetBinCenter(i), 0, 0, cont1);
1560 fPixArray->Add(pixPtr);
1565 //_____________________________________________________________________________
1566 AliMUONClusterFinderMLEM&
1567 AliMUONClusterFinderMLEM::operator=(const AliMUONClusterFinderMLEM& rhs)
1569 /// Protected assignement operator
1571 if (this == &rhs) return *this;
1573 AliFatal("Not implemented.");
1578 //_____________________________________________________________________________
1579 void AliMUONClusterFinderMLEM::AddVirtualPad(AliMUONCluster& cluster)
1581 /// Add virtual pad (with small charge) to improve fit for clusters
1582 /// with number of pads == 2 per direction
1584 // Find out non-bending and bending planes
1585 Int_t nonb[2] = {1, 0}; // non-bending and bending cathodes
1587 TVector2 dim0 = cluster.MinPadDimensions(0, 0, kTRUE);
1588 TVector2 dim1 = cluster.MinPadDimensions(1, 0, kTRUE);
1589 if (dim0.X() < dim1.X() - fgkDistancePrecision) {
1594 Bool_t same = kFALSE;
1595 if (TMath::Abs(dim0.Y()-dim1.Y()) < fgkDistancePrecision) same = kTRUE; // the same pad size on both planes
1598 Bool_t check[2] = {kFALSE, kFALSE};
1600 nxy[0] = nxy[1] = 0;
1601 for (Int_t inb = 0; inb < 2; ++inb) {
1602 cn = cluster.NofPads(nonb[inb], 0, kTRUE);
1603 if (inb == 0 && cn.GetFirst() == 2) check[inb] = kTRUE; // check non-bending plane
1604 else if (inb == 1 && cn.GetSecond() == 2) check[inb] = kTRUE; // check bending plane
1606 nxy[0] = TMath::Max (nxy[0], cn.GetFirst());
1607 nxy[1] = TMath::Max (nxy[1], cn.GetSecond());
1608 if (inb == 0 && nxy[0] < 2) nonb[inb] = !nonb[inb];
1609 else if (inb == 1 && cn.GetSecond() < 2) nonb[inb] = !nonb[inb];
1613 if (nxy[0] > 2) check[0] = kFALSE;
1614 if (nxy[1] > 2) check[1] = kFALSE;
1616 if (!check[0] && !check[1]) return;
1618 for (Int_t inb = 0; inb < 2; ++inb) {
1619 if (!check[inb]) continue;
1621 // Find pads with maximum and next to maximum charges
1622 Int_t maxPads[2] = {-1, -1};
1623 Double_t amax[2] = {0};
1624 Int_t mult = cluster.Multiplicity();
1625 for (Int_t j = 0; j < mult; ++j) {
1626 AliMUONPad *pad = cluster.Pad(j);
1627 if (pad->Cathode() != nonb[inb]) continue;
1628 for (Int_t j2 = 0; j2 < 2; ++j2) {
1629 if (pad->Charge() > amax[j2]) {
1630 if (j2 == 0) { amax[1] = amax[0]; maxPads[1] = maxPads[0]; }
1631 amax[j2] = pad->Charge();
1638 // Find min and max dimensions of the cluster
1639 Double_t limits[2] = {9999, -9999};
1640 for (Int_t j = 0; j < mult; ++j) {
1641 AliMUONPad *pad = cluster.Pad(j);
1642 if (pad->Cathode() != nonb[inb]) continue;
1643 if (pad->Coord(inb) < limits[0]) limits[0] = pad->Coord(inb);
1644 if (pad->Coord(inb) > limits[1]) limits[1] = pad->Coord(inb);
1647 // Loop over max and next to max pads
1648 Bool_t add = kFALSE;
1650 for (Int_t j = 0; j < 2; ++j) {
1652 if (maxPads[j] < 0) continue;
1654 if (amax[1] / amax[0] < 0.5) break;
1656 // Check if pad at the cluster limit
1657 AliMUONPad *pad = cluster.Pad(maxPads[j]);
1659 if (TMath::Abs(pad->Coord(inb)-limits[0]) < fgkDistancePrecision) idir = -1;
1660 else if (TMath::Abs(pad->Coord(inb)-limits[1]) < fgkDistancePrecision) idir = 1;
1662 //cout << " *** Pad not at the cluster limit: " << j << endl;
1665 if (j == 1 && idir == idirAdd) break; // this pad has been already added
1667 // Add pad (if it exists)
1669 if (inb == 0) pos.Set (pad->X() + idir * (pad->DX()+fgkDistancePrecision), pad->Y());
1670 else pos.Set (pad->X(), pad->Y() + idir * (pad->DY()+fgkDistancePrecision));
1671 //AliMpPad mppad = fSegmentation[nonb[inb]]->PadByPosition(pos,kTRUE);
1672 AliMpPad mppad = fSegmentation[nonb[inb]]->PadByPosition(pos,kFALSE);
1673 if (!mppad.IsValid()) continue; // non-existing pad
1674 cn = mppad.GetIndices();
1675 AliMUONPad muonPad(fDetElemId, nonb[inb], cn.GetFirst(), cn.GetSecond(),
1676 mppad.Position().X(), mppad.Position().Y(),
1677 mppad.Dimensions().X(), mppad.Dimensions().Y(), 0);
1678 if (inb == 0) muonPad.SetCharge(TMath::Min (amax[j]/100, 5.));
1679 //else muonPad.SetCharge(TMath::Min (amax[j]/15, fgkZeroSuppression));
1680 else muonPad.SetCharge(TMath::Min (amax[j]/15, 6.));
1681 if (muonPad.Charge() < 1.) muonPad.SetCharge(1.);
1682 muonPad.SetReal(kFALSE);
1683 if (fDebug) printf(" ***** Add virtual pad in %d direction ***** %f %f %f %3d %3d %f %f \n",
1684 inb, muonPad.Charge(), muonPad.X(), muonPad.Y(), muonPad.Ix(),
1685 muonPad.Iy(), muonPad.DX(), muonPad.DY());
1686 cluster.AddPad(muonPad); // add pad to the cluster
1693 //_____________________________________________________________________________
1694 void AliMUONClusterFinderMLEM::PadsInXandY(AliMUONCluster& cluster,
1695 Int_t &nInX, Int_t &nInY) const
1697 /// Find number of pads in X and Y-directions (excluding virtual ones and
1700 //Int_t statusToTest = 1;
1701 Int_t statusToTest = fgkUseForFit;
1703 //if ( nInX < 0 ) statusToTest = 0;
1704 if ( nInX < 0 ) statusToTest = fgkZero;
1706 Bool_t mustMatch(kTRUE);
1708 AliMpIntPair cn = cluster.NofPads(statusToTest,mustMatch);
1710 nInX = cn.GetFirst();
1711 nInY = cn.GetSecond();
1714 //_____________________________________________________________________________
1715 void AliMUONClusterFinderMLEM::RemovePixel(Int_t i)
1717 /// Remove pixel at index i
1718 AliMUONPad* pixPtr = Pixel(i);
1719 fPixArray->RemoveAt(i);
1723 //_____________________________________________________________________________
1724 void AliMUONClusterFinderMLEM::Simple(AliMUONCluster& cluster)
1726 /// Process simple cluster (small number of pads) without EM-procedure
1728 Int_t nForFit = 1, clustFit[1] = {0}, nfit;
1729 Double_t parOk[3] = {0.};
1730 TObjArray *clusters[1];
1731 clusters[0] = fPixArray;
1733 AliDebug(1,Form("nPix=%d",fPixArray->GetLast()+1));
1735 Int_t mult = cluster.Multiplicity();
1736 for (Int_t i = 0; i < mult; ++i)
1738 AliMUONPad* pad = cluster.Pad(i);
1740 if ( pad->IsSaturated())
1749 if (!pad->IsSaturated()) pad->SetStatus(fgkUseForFit);
1751 nfit = fSplitter->Fit(cluster,1, nForFit, clustFit, clusters, parOk, fClusterList);
1754 //_____________________________________________________________________________
1756 AliMUONClusterFinderMLEM::Pixel(Int_t i) const
1758 /// Returns pixel at index i
1759 return static_cast<AliMUONPad*>(fPixArray->UncheckedAt(i));
1762 //_____________________________________________________________________________
1764 AliMUONClusterFinderMLEM::Print(Option_t* what) const
1767 TString swhat(what);
1769 if ( swhat.Contains("precluster") )
1771 if ( fPreCluster) fPreCluster->Print();