AliRoot Core  v5-06-30 (35d6c57)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliMUONClusterFinderMLEM.cxx
Go to the documentation of this file.
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 //-----------------------------------------------------------------------------
29 //-----------------------------------------------------------------------------
30 
32 #include "AliLog.h"
33 #include "AliMUONCluster.h"
35 #include "AliMUONVDigit.h"
36 #include "AliMUONPad.h"
38 #include "AliMpPad.h"
39 #include "AliMpVPadIterator.h"
40 #include "AliMpVSegmentation.h"
41 #include "AliRunLoader.h"
42 #include "AliMUONVDigitStore.h"
43 #include <Riostream.h>
44 #include <TH2.h>
45 #include <TMinuit.h>
46 #include <TCanvas.h>
47 #include <TMath.h>
48 #include "AliCodeTimer.h"
49 
50 using std::endl;
51 using std::cout;
55 
56 const Double_t AliMUONClusterFinderMLEM::fgkDistancePrecision = 1e-3; // (cm) used to check overlaps and so on
57 const TVector2 AliMUONClusterFinderMLEM::fgkIncreaseSize(-AliMUONClusterFinderMLEM::fgkDistancePrecision,-AliMUONClusterFinderMLEM::fgkDistancePrecision);
58 const TVector2 AliMUONClusterFinderMLEM::fgkDecreaseSize(AliMUONClusterFinderMLEM::fgkDistancePrecision,AliMUONClusterFinderMLEM::fgkDistancePrecision);
59 
60 // Status flags for pads
61 const Int_t AliMUONClusterFinderMLEM::fgkZero = 0x0;
62 const Int_t AliMUONClusterFinderMLEM::fgkMustKeep = 0x1;
63 const Int_t AliMUONClusterFinderMLEM::fgkUseForFit = 0x10;
64 const Int_t AliMUONClusterFinderMLEM::fgkOver = 0x100;
65 const Int_t AliMUONClusterFinderMLEM::fgkModified = 0x1000;
66 const Int_t AliMUONClusterFinderMLEM::fgkCoupled = 0x10000;
67 
68 //_____________________________________________________________________________
71 fPreClusterFinder(clusterFinder),
72 fPreCluster(0x0),
73 fClusterList(),
74 fEventNumber(0),
75 fDetElemId(-1),
76 fClusterNumber(0),
77 fHistMlem(0x0),
78 fHistAnode(0x0),
79 fPixArray(new TObjArray(20)),
80 fDebug(0),
81 fPlot(plot),
82 fSplitter(0x0),
83 fNClusters(0),
84 fNAddVirtualPads(0),
85 fLowestPixelCharge(0),
86 fLowestPadCharge(0),
87 fLowestClusterCharge(0)
88 {
90 
91  fkSegmentation[1] = fkSegmentation[0] = 0x0;
92 
93  if (fPlot) fDebug = 1;
94 }
95 
96 //_____________________________________________________________________________
98 {
100  delete fPixArray; fPixArray = 0;
101 // delete fDraw;
102  delete fPreClusterFinder;
103  delete fSplitter;
104  AliInfo(Form("Total clusters %d AddVirtualPad needed %d",
106 }
107 
108 //_____________________________________________________________________________
109 Bool_t
111  TObjArray* pads[2],
112  const AliMpArea& area,
113  const AliMpVSegmentation* seg[2])
114 {
116 // AliCodeTimerAuto("",0)
117 
118  for ( Int_t i = 0; i < 2; ++i )
119  {
120  fkSegmentation[i] = seg[i];
121  }
122 
123  // Find out the DetElemId
124  fDetElemId = detElemId;
125 
126  delete fSplitter;
128  fPixArray,
133 
134  // find out current event number, and reset the cluster number
135  AliRunLoader *runLoader = AliRunLoader::Instance();
136  fEventNumber = runLoader ? runLoader->GetEventNumber() : 0;
137  fClusterNumber = -1;
138  fClusterList.Delete();
139  fPixArray->Delete();
140 
141  AliDebug(3,Form("EVT %d DE %d",fEventNumber,fDetElemId));
142 
144  {
145  return fPreClusterFinder->Prepare(detElemId,pads,area,seg);
146  }
147  else
148  {
149  return fPreClusterFinder->Prepare(detElemId,pads,area);
150  }
151 }
152 
153 //_____________________________________________________________________________
156 {
158 // AliCodeTimerAuto("",0)
159 
160  // if the list of clusters is not void, pick one from there
161  TObject* o(0x0);
162 
163  // do we have clusters in our list ?
164  if ( fClusterNumber < fClusterList.GetLast() )
165  {
166  o = fClusterList.At(++fClusterNumber);
167  }
168 
169  if ( o != 0x0 ) return static_cast<AliMUONCluster*>(o);
170 
171  //FIXME : at this point, must check whether we've used all the digits
172  //from precluster : if not, let the preclustering know about those unused
173  //digits, so it can reuse them
174 
175  // if the cluster list is exhausted, we need to go to the next
176  // pre-cluster and treat it
177 
179 
180  fPixArray->Delete();
181  fClusterList.Delete(); // reset the list of clusters for this pre-cluster
182  fClusterNumber = -1; //AZ
183 
184  if (!fPreCluster)
185  {
186  // we are done
187  return 0x0;
188  }
189 
191 
192  // WorkOnPreCluster may have used only part of the pads, so we check that
193  // now, and let the unused pads be reused by the preclustering...
194 
195  Int_t mult = fPreCluster->Multiplicity();
196  for ( Int_t i = 0; i < mult; ++i )
197  {
198  AliMUONPad* pad = fPreCluster->Pad(i);
199  if ( !pad->IsUsed() )
200  {
201  fPreClusterFinder->UsePad(*pad);
202  }
203  }
204 
205  return NextCluster();
206 }
207 
208 //_____________________________________________________________________________
209 Bool_t
211 {
214 
215  // AliCodeTimerAuto("",0)
216 
217  if (fDebug) {
218  cout << " *** Event # " << fEventNumber
219  << " det. elem.: " << fDetElemId << endl;
220  for (Int_t j = 0; j < fPreCluster->Multiplicity(); ++j) {
221  AliMUONPad* pad = fPreCluster->Pad(j);
222  printf(" bbb %3d %1d %8.4f %8.4f %8.4f %8.4f %6.1f %3d %3d %2d %1d %1d \n",
223  j, pad->Cathode(), pad->Coord(0), pad->Coord(1), pad->DX()*2, pad->DY()*2,
224  pad->Charge(), pad->Ix(), pad->Iy(), pad->Status(), pad->IsReal(), pad->IsSaturated());
225  }
226  }
227 
229  if (!cluster) return kFALSE;
230 
231  BuildPixArray(*cluster);
232 
233  if ( fPixArray->GetLast() < 0 )
234  {
235  AliDebug(1,"No pixel for the above cluster");
236  delete cluster;
237  return kFALSE;
238  }
239 
240  // Use MLEM for cluster finder
241  Int_t nMax = 1, localMax[100], maxPos[100];
242  Double_t maxVal[100];
243 
244  Int_t iSimple = 0, nInX = -1, nInY;
245 
246  PadsInXandY(*cluster,nInX, nInY);
247 
248  if (nInX < 4 && nInY < 4)
249  {
250  iSimple = 1; // simple cluster
251  }
252  else
253  {
254  nMax = FindLocalMaxima(fPixArray, localMax, maxVal); // for small clusters just to tag pixels
255  if (nMax > 1) {
256  if (cluster->Multiplicity() <= 50) nMax = 1; // for small clusters
257  if (nMax > 1) TMath::Sort(nMax, maxVal, maxPos, kTRUE); // in descending order
258  }
259  }
260 
261  for (Int_t i = 0; i < nMax; ++i)
262  {
263  if (nMax > 1)
264  {
265  FindCluster(*cluster,localMax, maxPos[i]);
266  }
267 
268  MainLoop(*cluster,iSimple);
269 
270  if (i < nMax-1)
271  {
272  Int_t mult = cluster->Multiplicity();
273  for (Int_t j = 0; j < mult; ++j)
274  {
275  AliMUONPad* pad = cluster->Pad(j);
276  //if ( pad->Status() == 0 ) continue; // pad charge was not modified
277  if ( pad->Status() != fgkOver ) continue; // pad was not used
278  //pad->SetStatus(0);
279  pad->SetStatus(fgkZero);
280  pad->RevertCharge(); // use backup charge value
281  }
282  }
283  } // for (Int_t i=0; i<nMax;
284  delete fHistMlem;
285  delete fHistAnode;
286  fHistMlem = fHistAnode = 0x0;
287  delete cluster;
288  return kTRUE;
289 }
290 
291 //_____________________________________________________________________________
292 Bool_t
294 {
296 
297  // make a fake pad from the pixel
298  AliMUONPad tmp(pad.DetElemId(),pad.Cathode(),pad.Ix(),pad.Iy(),
299  pixel.Coord(0),pixel.Coord(1),
300  pixel.Size(0),pixel.Size(1),0);
301 
303 }
304 
305 //_____________________________________________________________________________
308 {
311 
312  AliCodeTimerAuto("",0)
313 
314  // Disregard small clusters (leftovers from splitting or noise)
315  if ((origCluster.Multiplicity()==1 || origCluster.Multiplicity()==2) &&
316  origCluster.Charge(0)+origCluster.Charge(1) < fLowestClusterCharge )
317  {
318  return 0x0;
319  }
320 
321  AliMUONCluster* cluster = new AliMUONCluster(origCluster);
322 
323  AliDebug(2,"Start of CheckPreCluster=");
324  //StdoutToAliDebug(2,cluster->Print("full"));
325 
326  AliMUONCluster* rv(0x0);
327 
328  if (cluster->Multiplicity(0) && cluster->Multiplicity(1))
329  {
330  rv = CheckPreclusterTwoCathodes(cluster);
331  }
332  else
333  {
334  rv = cluster;
335  }
336  return rv;
337 }
338 
339 //_____________________________________________________________________________
342 {
344 
345  Int_t npad = cluster->Multiplicity();
346  Int_t* flags = new Int_t[npad];
347  for (Int_t j = 0; j < npad; ++j) flags[j] = 0;
348 
349  // Check pad overlaps
350  for ( Int_t i = 0; i < npad; ++i)
351  {
352  AliMUONPad* padi = cluster->Pad(i);
353  if ( padi->Cathode() != 0 ) continue;
354  for (Int_t j = i+1; j < npad; ++j)
355  {
356  AliMUONPad* padj = cluster->Pad(j);
357  if ( padj->Cathode() != 1 ) continue;
358  if ( !AliMUONPad::AreOverlapping(*padi,*padj,fgkDecreaseSize) ) continue;
359  flags[i] = flags[j] = 1; // mark overlapped pads
360  }
361  }
362 
363  // Check if all pads overlap
364  Int_t nFlags=0;
365  for (Int_t i = 0; i < npad; ++i)
366  {
367  if (!flags[i]) ++nFlags;
368  }
369 
370  if (nFlags > 0)
371  {
372  // not all pads overlap.
373  if (fDebug) cout << " nFlags: " << nFlags << endl;
374  TObjArray toBeRemoved;
375  for (Int_t i = 0; i < npad; ++i)
376  {
377  AliMUONPad* pad = cluster->Pad(i);
378  if (flags[i]) continue;
379  Int_t cath = pad->Cathode();
380  Int_t cath1 = TMath::Even(cath);
381  // Check for edge effect (missing pads on the _other_ cathode)
382  AliMpPad mpPad =
383  fkSegmentation[cath1]->PadByPosition(pad->Position().X(),
384  pad->Position().Y(),kFALSE);
385  if (!mpPad.IsValid()) continue;
386  if (nFlags == 1 && pad->Charge() < fLowestPadCharge) continue;
387  AliDebug(2,Form("Releasing the following pad : de,cath,ix,iy %d,%d,%d,%d charge %e",
388  fDetElemId,pad->Cathode(),pad->Ix(),pad->Iy(),pad->Charge()));
389  toBeRemoved.AddLast(pad);
390  fPreCluster->Pad(i)->Release();
391  }
392  Int_t nRemove = toBeRemoved.GetEntriesFast();
393  for ( Int_t i = 0; i < nRemove; ++i )
394  {
395  cluster->RemovePad(static_cast<AliMUONPad*>(toBeRemoved.UncheckedAt(i)));
396  }
397  }
398 
399  // Check correlations of cathode charges
400  if ( !cluster->IsSaturated() && cluster->ChargeAsymmetry() > 1 )
401  {
402  // big difference
403  Int_t cathode = cluster->MaxRawChargeCathode();
404  Int_t imin(-1);
405  Int_t imax(-1);
406  Double_t cmax(0);
407  Double_t cmin(1E9);
408 
409  // get min and max pad charges on the cathode opposite to the
410  // max pad (given by MaxRawChargeCathode())
411  //
412  Int_t mult = cluster->Multiplicity();
413  for ( Int_t i = 0; i < mult; ++i )
414  {
415  AliMUONPad* pad = cluster->Pad(i);
416  if ( pad->Cathode() != cathode || !pad->IsReal() )
417  {
418  // only consider pads in the opposite cathode, and
419  // only consider real pads (i.e. exclude the virtual ones)
420  continue;
421  }
422  if ( pad->Charge() < cmin )
423  {
424  cmin = pad->Charge();
425  imin = i;
426  if (imax < 0) {
427  imax = imin;
428  cmax = cmin;
429  }
430  }
431  else if ( pad->Charge() > cmax )
432  {
433  cmax = pad->Charge();
434  imax = i;
435  }
436  }
437  AliDebug(2,Form("Pad imin,imax %d,%d cmin,cmax %e,%e",
438  imin,imax,cmin,cmax));
439  //
440  // arrange pads according to their distance to the max, normalized
441  // to the pad size
442  Double_t* dist = new Double_t[mult];
443  Double_t dxMin(1E9);
444  Double_t dyMin(1E9);
445  Double_t dmin(0);
446 
447  AliMUONPad* padmax = cluster->Pad(imax);
448 
449  for ( Int_t i = 0; i < mult; ++i )
450  {
451  dist[i] = 0.0;
452  if ( i == imax) continue;
453  AliMUONPad* pad = cluster->Pad(i);
454  if ( pad->Cathode() != cathode || !pad->IsReal() ) continue;
455  Double_t dx = (pad->X()-padmax->X())/padmax->DX()/2.0;
456  Double_t dy = (pad->Y()-padmax->Y())/padmax->DY()/2.0;
457  dist[i] = TMath::Sqrt(dx*dx+dy*dy);
458  if ( i == imin )
459  {
460  dmin = dist[i] + 1E-3; // distance to the pad with minimum charge
461  dxMin = dx;
462  dyMin = dy;
463  }
464  }
465 
466  TMath::Sort(mult,dist,flags,kFALSE); // in ascending order
467  Double_t xmax(-1), distPrev(999);
468  TObjArray toBeRemoved;
469 
470  for ( Int_t i = 0; i < mult; ++i )
471  {
472  Int_t indx = flags[i];
473  AliMUONPad* pad = cluster->Pad(indx);
474  if ( pad->Cathode() != cathode || !pad->IsReal() ) continue;
475  if ( dist[indx] > dmin )
476  {
477  // farther than the minimum pad
478  Double_t dx = (pad->X()-padmax->X())/padmax->DX()/2.0;
479  Double_t dy = (pad->Y()-padmax->Y())/padmax->DY()/2.0;
480  dx *= dxMin;
481  dy *= dyMin;
482  if (dx >= 0 && dy >= 0) continue;
483  if (TMath::Abs(dx) > TMath::Abs(dy) && dx >= 0) continue;
484  if (TMath::Abs(dy) > TMath::Abs(dx) && dy >= 0) continue;
485  }
486  if (dist[indx] > distPrev + 1) break; // overstepping empty pads
487  if ( pad->Charge() <= cmax || TMath::Abs(dist[indx]-xmax) < 1E-3 )
488  {
489  // release pad
490  if (TMath::Abs(dist[indx]-xmax) < 1.e-3)
491  {
492  cmax = TMath::Max(pad->Charge(),cmax);
493  }
494  else
495  {
496  cmax = pad->Charge();
497  }
498  xmax = dist[indx];
499  distPrev = dist[indx];
500  AliDebug(2,Form("Releasing the following pad : de,cath,ix,iy %d,%d,%d,%d charge %e",
501  fDetElemId,pad->Cathode(),pad->Ix(),pad->Iy(),
502  pad->Charge()));
503 
504  toBeRemoved.AddLast(pad);
505  fPreCluster->Pad(indx)->Release();
506  }
507  }
508  Int_t nRemove = toBeRemoved.GetEntriesFast();
509  for ( Int_t i = 0; i < nRemove; ++i )
510  {
511  cluster->RemovePad(static_cast<AliMUONPad*>(toBeRemoved.UncheckedAt(i)));
512  }
513  delete[] dist;
514  }
515 
516  delete[] flags;
517 
518  AliDebug(2,"End of CheckPreClusterTwoCathodes=");
519  //StdoutToAliDebug(2,cluster->Print("full"));
520 
521  return cluster;
522 }
523 
524 //_____________________________________________________________________________
525 void
527 {
529 
530  Int_t nPix = fPixArray->GetLast()+1;
531  Int_t dummy(0);
532 
533  for ( Int_t i = 0; i < nPix; ++i )
534  {
535  AliMUONPad* pixelI = Pixel(i);
536  AliMUONPad pi(dummy,dummy,dummy,dummy,
537  pixelI->Coord(0),pixelI->Coord(1),
538  pixelI->Size(0),pixelI->Size(1),0.0);
539 
540  for ( Int_t j = i+1; j < nPix; ++j )
541  {
542  AliMUONPad* pixelJ = Pixel(j);
543  AliMUONPad pj(dummy,dummy,dummy,dummy,
544  pixelJ->Coord(0),pixelJ->Coord(1),
545  pixelJ->Size(0),pixelJ->Size(1),0.0);
546  AliMpArea area;
547 
549  {
550  AliInfo(Form("The following 2 pixels (%d and %d) overlap !",i,j));
551  /*
552  StdoutToAliInfo(pixelI->Print();
553  cout << " Surface = " << pixelI->Size(0)*pixelI->Size(1)*4 << endl;
554  pixelJ->Print();
555  cout << " Surface = " << pixelJ->Size(0)*pixelJ->Size(1)*4 << endl;
556  cout << " Area surface = " << area.GetDimensionX()*area.GetDimensionY()*4 << endl;
557  cout << "-------" << endl;
558  );
559  */
560  }
561  }
562  }
563 }
564 
565 //_____________________________________________________________________________
567 {
569 
570  Int_t npad = cluster.Multiplicity();
571  if (npad<=0)
572  {
573  AliWarning("Got no pad at all ?!");
574  }
575 
576  fPixArray->Delete();
577  BuildPixArrayOneCathode(cluster);
578 
579  Int_t nPix = fPixArray->GetLast()+1;
580 
581 // AliDebug(2,Form("nPix after BuildPixArray=%d",nPix));
582 
583  if ( nPix > npad )
584  {
585 // AliDebug(2,Form("Will trim number of pixels to number of pads"));
586 
587  // Too many pixels - sort and remove pixels with the lowest signal
588  fPixArray->Sort();
589  for ( Int_t i = npad; i < nPix; ++i )
590  {
591  RemovePixel(i);
592  }
593  fPixArray->Compress();
594  } // if (nPix > npad)
595 
596 // StdoutToAliDebug(2,cout << "End of BuildPixelArray:" << endl;
597 // fPixArray->Print(););
598  //CheckOverlaps();//FIXME : this is for debug only. Remove it.
599 }
600 
601 //_____________________________________________________________________________
603 {
605 
606 // AliDebug(2,Form("cluster.Multiplicity=%d",cluster.Multiplicity()));
607 
608  TVector2 dim = cluster.MinPadDimensions (-1, kFALSE);
609  Double_t width[2] = {dim.X(), dim.Y()}, xy0[2]={99999,99999};
610  Int_t found[2] = {0,0}, mult = cluster.Multiplicity();
611 
612  for ( Int_t i = 0; i < mult; ++i) {
613  AliMUONPad* pad = cluster.Pad(i);
614  for (Int_t j = 0; j < 2; ++j) {
615  if (found[j] == 0 && TMath::Abs(pad->Size(j)-width[j]) < fgkDistancePrecision) {
616  xy0[j] = pad->Coord(j);
617  found[j] = 1;
618  }
619  }
620  if (found[0] && found[1]) break;
621  }
622 
623  Double_t min[2], max[2];
624  Int_t cath0 = 0, cath1 = 1;
625  if (cluster.Multiplicity(0) == 0) cath0 = 1;
626  else if (cluster.Multiplicity(1) == 0) cath1 = 0;
627 
628 
629  Double_t leftDownX, leftDownY;
630  cluster.Area(cath0).LeftDownCorner(leftDownX, leftDownY);
631  Double_t rightUpX, rightUpY;
632  cluster.Area(cath0).RightUpCorner(rightUpX, rightUpY);
633  min[0] = leftDownX;
634  min[1] = leftDownY;
635  max[0] = rightUpX;
636  max[1] = rightUpY;;
637  if (cath1 != cath0) {
638  cluster.Area(cath1).LeftDownCorner(leftDownX, leftDownY);
639  cluster.Area(cath1).RightUpCorner(rightUpX, rightUpY);
640  min[0] = TMath::Max (min[0], leftDownX);
641  min[1] = TMath::Max (min[1], leftDownY);
642  max[0] = TMath::Min (max[0], rightUpX);
643  max[1] = TMath::Min (max[1], rightUpY);
644  }
645 
646  // Adjust limits
647  //width[0] /= 2; width[1] /= 2; // just for check
648  Int_t nbins[2]={0,0};
649  for (Int_t i = 0; i < 2; ++i) {
650  Double_t dist = (min[i] - xy0[i]) / width[i] / 2;
651  if (TMath::Abs(dist) < 1.e-6) dist = -1.e-6;
652  min[i] = xy0[i] + (TMath::Nint(dist-TMath::Sign(1.e-6,dist))
653  + TMath::Sign(0.5,dist)) * width[i] * 2;
654  nbins[i] = TMath::Nint ((max[i] - min[i]) / width[i] / 2);
655  if (nbins[i] == 0) ++nbins[i];
656  max[i] = min[i] + nbins[i] * width[i] * 2;
657  //cout << dist << " " << min[i] << " " << max[i] << " " << nbins[i] << endl;
658  }
659 
660  // Book histogram
661  TH2D *hist1 = new TH2D ("Grid", "", nbins[0], min[0], max[0], nbins[1], min[1], max[1]);
662  TH2D *hist2 = new TH2D ("Entries", "", nbins[0], min[0], max[0], nbins[1], min[1], max[1]);
663  TAxis *xaxis = hist1->GetXaxis();
664  TAxis *yaxis = hist1->GetYaxis();
665 
666  // Fill histogram
667  for ( Int_t i = 0; i < mult; ++i) {
668  AliMUONPad* pad = cluster.Pad(i);
669  Int_t ix0 = xaxis->FindBin(pad->X());
670  Int_t iy0 = yaxis->FindBin(pad->Y());
671  PadOverHist(0, ix0, iy0, pad, hist1, hist2);
672  }
673 
674  // Store pixels
675  for (Int_t i = 1; i <= nbins[0]; ++i) {
676  Double_t x = xaxis->GetBinCenter(i);
677  for (Int_t j = 1; j <= nbins[1]; ++j) {
678  if (hist2->GetBinContent(hist2->GetBin(i,j)) < 0.1) continue;
679  //if (hist2->GetBinContent(hist2->GetBin(i,j)) < 1.1 && cluster.Multiplicity(0) &&
680  // cluster.Multiplicity(1)) continue;
681  if (cath0 != cath1) {
682  // Two-sided cluster
683  Double_t cont = hist2->GetBinContent(hist2->GetBin(i,j));
684  if (cont < 999.) continue;
685  if (cont-Int_t(cont/1000.)*1000. < 0.5) continue;
686  }
687  Double_t y = yaxis->GetBinCenter(j);
688  Double_t charge = hist1->GetBinContent(hist1->GetBin(i,j));
689  AliMUONPad* pixPtr = new AliMUONPad(x, y, width[0], width[1], charge);
690  fPixArray->Add(pixPtr);
691  }
692  }
693  //*
694  if (fPixArray->GetEntriesFast() == 1) {
695  // Split pixel into 2
696  AliMUONPad* pixPtr = static_cast<AliMUONPad*> (fPixArray->UncheckedAt(0));
697  pixPtr->SetSize(0,width[0]/2.);
698  pixPtr->Shift(0,-width[0]/4.);
699  pixPtr = new AliMUONPad(pixPtr->X()+width[0], pixPtr->Y(), width[0]/2., width[1], pixPtr->Charge());
700  fPixArray->Add(pixPtr);
701  }
702  //*/
703  //fPixArray->Print();
704  delete hist1;
705  delete hist2;
706 }
707 
708 //_____________________________________________________________________________
709 void AliMUONClusterFinderMLEM::PadOverHist(Int_t idir, Int_t ix0, Int_t iy0, AliMUONPad *pad,
710  TH2D *hist1, TH2D *hist2)
711 {
713 
714  TAxis *axis = idir == 0 ? hist1->GetXaxis() : hist1->GetYaxis();
715  Int_t nbins = axis->GetNbins(), cath = pad->Cathode();
716  Double_t bin = axis->GetBinWidth(1), amask = TMath::Power(1000.,cath*1.);
717 
718  Int_t nbinPad = (Int_t)(pad->Size(idir)/bin*2+fgkDistancePrecision) + 1; // number of bins covered by pad
719 
720  for (Int_t i = 0; i < nbinPad; ++i) {
721  Int_t ixy = idir == 0 ? ix0 + i : iy0 + i;
722  if (ixy > nbins) break;
723  Double_t lowEdge = axis->GetBinLowEdge(ixy);
724  if (lowEdge + fgkDistancePrecision > pad->Coord(idir) + pad->Size(idir)) break;
725  if (idir == 0) PadOverHist(1, ixy, iy0, pad, hist1, hist2); // span in the other direction
726  else {
727  // Fill histogram
728  Double_t cont = pad->Charge();
729  if (hist2->GetBinContent(hist2->GetBin(ix0, ixy)) > 0.1)
730  cont = TMath::Min (hist1->GetBinContent(hist1->GetBin(ix0, ixy)), cont);
731  hist1->SetBinContent(hist1->GetBin(ix0, ixy), cont);
732  //hist2->SetBinContent(hist2->GetBin(ix0, ixy), hist2->GetBinContent(hist2->GetBin(ix0, ixy))+1);
733  hist2->SetBinContent(hist2->GetBin(ix0, ixy), hist2->GetBinContent(hist2->GetBin(ix0, ixy))+amask);
734  }
735  }
736 
737  for (Int_t i = -1; i > -nbinPad; --i) {
738  Int_t ixy = idir == 0 ? ix0 + i : iy0 + i;
739  if (ixy < 1) break;
740  Double_t upEdge = axis->GetBinUpEdge(ixy);
741  if (upEdge - fgkDistancePrecision < pad->Coord(idir) - pad->Size(idir)) break;
742  if (idir == 0) PadOverHist(1, ixy, iy0, pad, hist1, hist2); // span in the other direction
743  else {
744  // Fill histogram
745  Double_t cont = pad->Charge();
746  if (hist2->GetBinContent(hist2->GetBin(ix0, ixy)) > 0.1)
747  cont = TMath::Min (hist1->GetBinContent(hist1->GetBin(ix0, ixy)), cont);
748  hist1->SetBinContent(hist1->GetBin(ix0, ixy), cont);
749  //hist2->SetBinContent(hist1->GetBin(ix0, ixy), hist2->GetBinContent(hist2->GetBin(ix0, ixy))+1);
750  hist2->SetBinContent(hist2->GetBin(ix0, ixy), hist2->GetBinContent(hist2->GetBin(ix0, ixy))+amask);
751  }
752  }
753 }
754 
755 //_____________________________________________________________________________
756 void
757 AliMUONClusterFinderMLEM::Plot(const char* /*basename*/)
758 {
760 
761  return; //AZ
762 // if (!fPlot) return;
763 //
764 // TCanvas* c = new TCanvas("MLEM","MLEM",800,600);
765 // c->Draw();
766 // Draw();
767 // c->Modified();
768 // c->Update();
769 // c->Print(Form("%s.EVT%d.DE%d.CLU%d.png",basename,fEventNumber,
770 // fDetElemId,fClusterNumber));
771 }
772 
773 //_____________________________________________________________________________
774 void
776  Double_t* coef,
777  Double_t* probi)
778 {
780 
781  Int_t npadTot = cluster.Multiplicity();
782  Int_t nPix = fPixArray->GetLast()+1;
783 
784  //memset(probi,0,nPix*sizeof(Double_t));
785  for (Int_t j = 0; j < npadTot*nPix; ++j) coef[j] = 0.;
786  for (Int_t j = 0; j < nPix; ++j) probi[j] = 0.;
787 
788  Int_t mult = cluster.Multiplicity();
789  for ( Int_t j = 0; j < mult; ++j )
790  {
791  AliMUONPad* pad = cluster.Pad(j);
792  Int_t indx = j*nPix;
793 
794  for ( Int_t ipix = 0; ipix < nPix; ++ipix )
795  {
796  Int_t indx1 = indx + ipix;
797  //if (pad->Status() < 0)
798  if (pad->Status() != fgkZero)
799  {
800  coef[indx1] = 0;
801  continue;
802  }
803  AliMUONPad* pixPtr = Pixel(ipix);
804  // coef is the charge (given by Mathieson integral) on pad, assuming
805  // the Mathieson is center at pixel.
806  coef[indx1] = fSplitter->ChargeIntegration(pixPtr->Coord(0), pixPtr->Coord(1), *pad);
807 // AliDebug(2,Form("pad=(%d,%d,%e,%e,%e,%e) pix=(%e,%e,%e,%e) coef %e",
808 // pad->Ix(),pad->Iy(),
809 // pad->X(),pad->Y(),
810 // pad->DX(),pad->DY(),
811 // pixPtr->Coord(0),pixPtr->Coord(1),
812 // pixPtr->Size(0),pixPtr->Size(1),
813 // coef[indx1]));
814 
815  probi[ipix] += coef[indx1];
816  }
817  }
818 }
819 
820 //_____________________________________________________________________________
821 Bool_t AliMUONClusterFinderMLEM::MainLoop(AliMUONCluster& cluster, Int_t iSimple)
822 {
824 
825  // AliCodeTimerAuto("",0)
826 
827  Int_t nPix = fPixArray->GetLast()+1;
828 
829  AliDebug(2,Form("nPix=%d iSimple=%d, precluster=",nPix,iSimple));
830  //StdoutToAliDebug(2,cluster.Print("full"););
831 
832  if ( nPix < 0 )
833  {
834  AliDebug(1,"No pixels, why am I here then ?");
835  }
836 
837  AddVirtualPad(cluster); // add virtual pads if necessary
838 
839  Int_t npadTot = cluster.Multiplicity();
840  Int_t npadOK = 0;
841  for (Int_t i = 0; i < npadTot; ++i)
842  {
843  //if (cluster.Pad(i)->Status() == 0) ++npadOK;
844  if (cluster.Pad(i)->Status() == fgkZero) ++npadOK;
845  }
846 
847  Double_t* coef(0x0);
848  Double_t* probi(0x0);
849  Int_t lc(0); // loop counter
850 
851  //Plot("mlem.start");
852  AliMUONPad* pixPtr = Pixel(0);
853  Double_t xylim[4] = {pixPtr->X(), -pixPtr->X(), pixPtr->Y(), -pixPtr->Y()};
854 
855  while (1)
856  {
857  ++lc;
858  delete fHistMlem;
859  fHistMlem = 0x0;
860 
861  AliDebug(2,Form("lc %d nPix %d(%d) npadTot %d npadOK %d",lc,nPix,fPixArray->GetLast()+1,npadTot,npadOK));
862  AliDebug(2,Form("EVT%d PixArray=",fEventNumber));
863  //StdoutToAliDebug(2,fPixArray->Print("full"));
864 
865  coef = new Double_t [npadTot*nPix];
866  probi = new Double_t [nPix];
867 
868  // Calculate coefficients and pixel visibilities
869  ComputeCoefficients(cluster,coef,probi);
870 
871  for (Int_t ipix = 0; ipix < nPix; ++ipix)
872  {
873  if (probi[ipix] < 0.01)
874  {
875  AliMUONPad* pixel = Pixel(ipix);
876  AliDebug(2,Form("Setting the following pixel to invisible as its probi<0.01:"));
877  //StdoutToAliDebug(2,cout << Form(" -- ipix %3d --- "); pixel->Print(););
878  pixel->SetCharge(0); // "invisible" pixel
879  }
880  }
881 
882  // MLEM algorithm
883  Mlem(cluster,coef, probi, 15);
884 
885  // Find histogram limits for the 1'st pass only - for others computed below
886  if (lc == 1) {
887  for ( Int_t ipix = 1; ipix < nPix; ++ipix )
888  {
889  pixPtr = Pixel(ipix);
890  for ( Int_t i = 0; i < 2; ++i )
891  {
892  Int_t indx = i * 2;
893  if (pixPtr->Coord(i) < xylim[indx]) xylim[indx] = pixPtr->Coord(i);
894  else if (-pixPtr->Coord(i) < xylim[indx+1]) xylim[indx+1] = -pixPtr->Coord(i);
895  }
896  }
897  } else pixPtr = Pixel(0);
898 
899  for (Int_t i = 0; i < 4; i++)
900  {
901  xylim[i] -= pixPtr->Size(i/2);
902  }
903 
904  Int_t nx = TMath::Nint ((-xylim[1]-xylim[0])/pixPtr->Size(0)/2);
905  Int_t ny = TMath::Nint ((-xylim[3]-xylim[2])/pixPtr->Size(1)/2);
906 
907  //StdoutToAliDebug(2,cout << "pixel used for nx,ny computation : "; pixPtr->Print(););
908  AliDebug(2,Form("lc %d pixPtr size = %e,%e nx,ny=%d,%d xylim=%e,%e,%e,%e",
909  lc,pixPtr->Size(0),pixPtr->Size(1),nx,ny,
910  xylim[0],-xylim[1],xylim[2],-xylim[3]
911  ));
912 
913  AliDebug(2,Form("LowestPadCharge=%e",fLowestPadCharge));
914 
915  delete fHistMlem;
916 
917  fHistMlem = new TH2D("mlem","mlem",nx,xylim[0],-xylim[1],ny,xylim[2],-xylim[3]);
918 
919  for (Int_t ipix = 0; ipix < nPix; ++ipix)
920  {
921  AliMUONPad* pixPtr2 = Pixel(ipix);
922  fHistMlem->Fill(pixPtr2->Coord(0),pixPtr2->Coord(1),pixPtr2->Charge());
923  }
924 
925  // Check if the total charge of pixels is too low
926  Double_t qTot = 0;
927  for ( Int_t i = 0; i < nPix; ++i)
928  {
929  qTot += Pixel(i)->Charge();
930  }
931 
932  if ( qTot < 1.e-4 || ( npadOK < 3 && qTot < fLowestClusterCharge ) )
933  {
934  AliDebug(1,Form("Deleting the above cluster (charge %e too low, npadOK=%d)",qTot,npadOK));
935  delete [] coef;
936  delete [] probi;
937  fPixArray->Delete();
938  for ( Int_t i = 0; i < npadTot; ++i)
939  {
940  AliMUONPad* pad = cluster.Pad(i);
941  //if ( pad->Status() == 0) pad->SetStatus(-1);
942  if ( pad->Status() == fgkZero) pad->SetStatus(fgkOver);
943  }
944  return kFALSE;
945  }
946 
947  if (iSimple)
948  {
949  // Simple cluster - skip further passes thru EM-procedure
950  Simple(cluster);
951  delete [] coef;
952  delete [] probi;
953  fPixArray->Delete();
954  return kTRUE;
955  }
956 
957  // Calculate position of the center-of-gravity around the maximum pixel
958  Double_t xyCOG[2];
959  FindCOG(xyCOG);
960 
961  if (TMath::Min(pixPtr->Size(0),pixPtr->Size(1)) < 0.07 &&
962  pixPtr->Size(0) > pixPtr->Size(1)) break;
963 
964  // Sort pixels according to the charge
965  MaskPeaks(1); // mask local maxima
966  fPixArray->Sort();
967  MaskPeaks(0); // unmask local maxima
968  Double_t pixMin = 0.01*Pixel(0)->Charge();
969  pixMin = TMath::Min(pixMin,100*fLowestPixelCharge);
970 
971  // Decrease pixel size and shift pixels to make them centered at
972  // the maximum one
973  Int_t indx = (pixPtr->Size(0)>pixPtr->Size(1)) ? 0 : 1;
974  Int_t ix(1);
975  Double_t width = 0;
976  Double_t shift[2] = { 0.0, 0.0 };
977  for (Int_t i = 0; i < 4; ++i) xylim[i] = 999;
978  Int_t nPix1 = nPix;
979  nPix = 0;
980  for (Int_t ipix = 0; ipix < nPix1; ++ipix)
981  {
982  AliMUONPad* pixPtr2 = Pixel(ipix);
983  if ( nPix >= npadOK // too many pixels already
984  ||
985  ((pixPtr2->Charge() < pixMin) && (pixPtr2->Status() != fgkMustKeep)) // too low charge
986  )
987  {
988  RemovePixel(ipix);
989  continue;
990  }
991  for (Int_t i = 0; i < 2; ++i)
992  {
993  if (!i)
994  {
995  pixPtr2->SetCharge(fLowestPadCharge);
996  pixPtr2->SetSize(indx, pixPtr2->Size(indx)/2);
997  width = -pixPtr2->Size(indx);
998  pixPtr2->Shift(indx, width);
999  // Shift pixel position
1000  if (ix)
1001  {
1002  ix = 0;
1003  for (Int_t j = 0; j < 2; ++j)
1004  {
1005  shift[j] = pixPtr2->Coord(j) - xyCOG[j];
1006  shift[j] -= ((Int_t)(shift[j]/pixPtr2->Size(j)/2))*pixPtr2->Size(j)*2;
1007  }
1008  } // if (ix)
1009  pixPtr2->Shift(0, -shift[0]);
1010  pixPtr2->Shift(1, -shift[1]);
1011  ++nPix;
1012  }
1013  else if (nPix < npadOK)
1014  {
1015  pixPtr2 = new AliMUONPad(*pixPtr2);
1016  pixPtr2->Shift(indx, -2*width);
1017  pixPtr2->SetStatus(fgkZero);
1018  fPixArray->Add(pixPtr2);
1019  ++nPix;
1020  }
1021  else continue; // skip adjustment of histo limits
1022  for (Int_t j = 0; j < 4; ++j)
1023  {
1024  xylim[j] = TMath::Min (xylim[j], (j%2 ? -1 : 1)*pixPtr2->Coord(j/2));
1025  }
1026  } // for (Int_t i=0; i<2;
1027  } // for (Int_t ipix=0;
1028 
1029  fPixArray->Compress();
1030 
1031  AliDebug(2,Form("After shift:"));
1032  //StdoutToAliDebug(2,fPixArray->Print("","full"););
1033  //Plot(Form("mlem.lc%d",lc+1));
1034 
1035  AliDebug(2,Form(" xyCOG=%9.6f %9.6f xylim=%9.6f,%9.6f,%9.6f,%9.6f",
1036  xyCOG[0],xyCOG[1],
1037  xylim[0],xylim[1],
1038  xylim[2],xylim[3]));
1039 
1040  if (nPix < npadOK)
1041  {
1042  AliMUONPad* pixPtr2 = Pixel(0);
1043  // add pixels if the maximum is at the limit of pixel area:
1044  // start from Y-direction
1045  Int_t j = 0;
1046  for (Int_t i = 3; i > -1; --i)
1047  {
1048  if (nPix < npadOK &&
1049  TMath::Abs((i%2 ? -1 : 1)*xylim[i]-xyCOG[i/2]) < pixPtr2->Size(i/2))
1050  {
1051  //AliMUONPad* p = static_cast<AliMUONPad*>(pixPtr->Clone());
1052  AliMUONPad* p = new AliMUONPad(*pixPtr2);
1053  p->SetCoord(i/2, xyCOG[i/2]+(i%2 ? 2:-2)*pixPtr2->Size(i/2));
1054  xylim[i] = p->Coord(i/2) * (i%2 ? -1 : 1); // update histo limits
1055  j = TMath::Even (i/2);
1056  p->SetCoord(j, xyCOG[j]);
1057  AliDebug(2,Form("Adding pixel on the edge (i=%d) ",i));
1058  //StdoutToAliDebug(2,cout << " ---- ";
1059  // p->Print("corners"););
1060  fPixArray->Add(p);
1061  ++nPix;
1062  }
1063  }
1064  }
1065  nPix = fPixArray->GetEntriesFast();
1066  delete [] coef;
1067  delete [] probi;
1068  } // while (1)
1069 
1070  AliDebug(2,Form("At the end of while loop nPix=%d : ",fPixArray->GetLast()+1));
1071  //StdoutToAliDebug(2,fPixArray->Print("","full"););
1072 
1073  // remove pixels with low signal or low visibility
1074  // Cuts are empirical !!!
1075  Double_t thresh = TMath::Max (fHistMlem->GetMaximum()/100.,2.0*fLowestPixelCharge);
1076  thresh = TMath::Min (thresh,100.0*fLowestPixelCharge);
1077  Double_t charge = 0;
1078 
1079  // Mark pixels which should be removed
1080  for (Int_t i = 0; i < nPix; ++i)
1081  {
1082  AliMUONPad* pixPtr2 = Pixel(i);
1083  charge = pixPtr2->Charge();
1084  if (charge < thresh)
1085  {
1086  pixPtr2->SetCharge(-charge);
1087  }
1088  }
1089 
1090  // Move charge of removed pixels to their nearest neighbour (to keep total charge the same)
1091  Int_t near = 0;
1092  for (Int_t i = 0; i < nPix; ++i)
1093  {
1094  AliMUONPad* pixPtr2 = Pixel(i);
1095  charge = pixPtr2->Charge();
1096  if (charge > 0) continue;
1097  near = FindNearest(pixPtr2);
1098  pixPtr2->SetCharge(0);
1099  probi[i] = 0; // make it "invisible"
1100  AliMUONPad* pnear = Pixel(near);
1101  pnear->SetCharge(pnear->Charge() + (-charge));
1102  }
1103  Mlem(cluster,coef,probi,2);
1104 
1105  AliDebug(2,Form("Before splitting nPix=%d EVT %d DE %d",fPixArray->GetLast()+1,fEventNumber,fDetElemId));
1106  //StdoutToAliDebug(2,fPixArray->Print("","full"););
1107  //Plot("mlem.beforesplit");
1108 
1109  // Update histogram
1110  for (Int_t i = 0; i < nPix; ++i)
1111  {
1112  AliMUONPad* pixPtr2 = Pixel(i);
1113  Int_t ix = fHistMlem->GetXaxis()->FindBin(pixPtr2->Coord(0));
1114  Int_t iy = fHistMlem->GetYaxis()->FindBin(pixPtr2->Coord(1));
1115  fHistMlem->SetBinContent(ix, iy, pixPtr2->Charge());
1116  }
1117 
1118  // Try to split into clusters
1119  Bool_t ok = kTRUE;
1120  if (fHistMlem->GetSum() < 2.0*fLowestPixelCharge)
1121  {
1122  ok = kFALSE;
1123  }
1124  else
1125  {
1126  fSplitter->Split(cluster,fHistMlem,coef,fClusterList);
1127  }
1128 
1129  delete [] coef;
1130  delete [] probi;
1131  fPixArray->Delete();
1132 
1133  return ok;
1134 }
1135 
1136 //_____________________________________________________________________________
1138 {
1141 
1142  for (Int_t i = 0; i < fPixArray->GetEntriesFast(); ++i) {
1143  AliMUONPad* pix = Pixel(i);
1144  if (pix->Status() == fgkMustKeep) {
1145  if (mask == 1) pix->SetCharge(pix->Charge()+10000.);
1146  else pix->SetCharge(pix->Charge()-10000.);
1147  }
1148  }
1149 }
1150 
1151 //_____________________________________________________________________________
1153  const Double_t* coef, Double_t* probi,
1154  Int_t nIter)
1155 {
1157 
1158  Int_t nPix = fPixArray->GetEntriesFast();
1159 
1160  Int_t npad = cluster.Multiplicity();
1161 
1162  Double_t* probi1 = new Double_t[nPix];
1163  Double_t probMax = TMath::MaxElement(nPix,probi);
1164 
1165  for (Int_t iter = 0; iter < nIter; ++iter)
1166  {
1167  // Do iterations
1168  for (Int_t ipix = 0; ipix < nPix; ++ipix)
1169  {
1170  Pixel(ipix)->SetChargeBackup(0);
1171  // Correct each pixel
1172  probi1[ipix] = 0;
1173  if (probi[ipix] < 0.01) continue; // skip "invisible" pixel
1174  Double_t sum = 0;
1175  probi1[ipix] = probMax;
1176  for (Int_t j = 0; j < npad; ++j)
1177  {
1178  AliMUONPad* pad = cluster.Pad(j);
1179  //if (pad->Status() < 0) continue;
1180  if (pad->Status() != fgkZero) continue;
1181  Double_t sum1 = 0;
1182  Int_t indx1 = j*nPix;
1183  Int_t indx = indx1 + ipix;
1184  // Calculate expectation
1185  for (Int_t i = 0; i < nPix; ++i)
1186  {
1187  sum1 += Pixel(i)->Charge()*coef[indx1+i];
1188  //cout << i << " " << Pixel(i)->Charge() << " " << coef[indx1+i] << endl;
1189  }
1190  if ( pad->IsSaturated() && sum1 > pad->Charge() )
1191  {
1192  // correct for pad charge overflows
1193  probi1[ipix] -= coef[indx];
1194  continue;
1195  //sum1 = pad->Charge();
1196  }
1197 
1198  if (sum1 > 1.e-6) sum += pad->Charge()*coef[indx]/sum1;
1199  } // for (Int_t j=0;
1200  AliMUONPad* pixPtr = Pixel(ipix);
1201  if (probi1[ipix] > 1.e-6)
1202  {
1203  //AZ pixPtr->SetCharge(pixPtr->Charge()*sum/probi1[ipix]);
1204  pixPtr->SetChargeBackup(pixPtr->Charge()*sum/probi1[ipix]);
1205  }
1206  //cout << " xxx " << ipix << " " << pixPtr->Charge() << " " << pixPtr->ChargeBackup() << " " << sum << " " << probi1[ipix] << endl;
1207  } // for (Int_t ipix=0;
1208  Double_t qTot = 0;
1209  for (Int_t i = 0; i < nPix; ++i) {
1210  AliMUONPad* pixPtr = Pixel(i);
1211  pixPtr->RevertCharge();
1212  qTot += pixPtr->Charge();
1213  }
1214  if (qTot < 1.e-6) {
1215  // Can happen in clusters with large number of overflows - speeding up
1216  delete [] probi1;
1217  return;
1218  }
1219  } // for (Int_t iter=0;
1220  delete [] probi1;
1221 }
1222 
1223 //_____________________________________________________________________________
1225 {
1227 
1228  Int_t ixmax, iymax, ix, nsumx=0, nsumy=0, nsum=0;
1229  Int_t i1 = -9, j1 = -9;
1230  fHistMlem->GetMaximumBin(ixmax,iymax,ix);
1231  Int_t nx = fHistMlem->GetNbinsX();
1232  Int_t ny = fHistMlem->GetNbinsY();
1233  Double_t thresh = fHistMlem->GetMaximum()/10;
1234  Double_t x, y, cont, xq=0, yq=0, qq=0;
1235 
1236  Int_t ie = TMath::Min(ny,iymax+1), je = TMath::Min(nx,ixmax+1);
1237  for (Int_t i = TMath::Max(1,iymax-1); i <= ie; ++i) {
1238  y = fHistMlem->GetYaxis()->GetBinCenter(i);
1239  for (Int_t j = TMath::Max(1,ixmax-1); j <= je; ++j) {
1240  cont = fHistMlem->GetBinContent(fHistMlem->GetBin(j,i));
1241  if (cont < thresh) continue;
1242  if (i != i1) {i1 = i; nsumy++;}
1243  if (j != j1) {j1 = j; nsumx++;}
1244  x = fHistMlem->GetXaxis()->GetBinCenter(j);
1245  xq += x*cont;
1246  yq += y*cont;
1247  qq += cont;
1248  nsum++;
1249  }
1250  }
1251 
1252  Double_t cmax = 0;
1253  Int_t i2 = 0, j2 = 0;
1254  x = y = 0;
1255  if (nsumy == 1) {
1256  // one bin in Y - add one more (with the largest signal)
1257  for (Int_t i = TMath::Max(1,iymax-1); i <= ie; ++i) {
1258  if (i == iymax) continue;
1259  for (Int_t j = TMath::Max(1,ixmax-1); j <= je; ++j) {
1260  cont = fHistMlem->GetBinContent(fHistMlem->GetBin(j,i));
1261  if (cont > cmax) {
1262  cmax = cont;
1263  x = fHistMlem->GetXaxis()->GetBinCenter(j);
1264  y = fHistMlem->GetYaxis()->GetBinCenter(i);
1265  i2 = i;
1266  j2 = j;
1267  }
1268  }
1269  }
1270  xq += x*cmax;
1271  yq += y*cmax;
1272  qq += cmax;
1273  if (i2 != i1) nsumy++;
1274  if (j2 != j1) nsumx++;
1275  nsum++;
1276  } // if (nsumy == 1)
1277 
1278  if (nsumx == 1) {
1279  // one bin in X - add one more (with the largest signal)
1280  cmax = x = y = 0;
1281  for (Int_t j = TMath::Max(1,ixmax-1); j <= je; ++j) {
1282  if (j == ixmax) continue;
1283  for (Int_t i = TMath::Max(1,iymax-1); i <= ie; ++i) {
1284  cont = fHistMlem->GetBinContent(fHistMlem->GetBin(j,i));
1285  if (cont > cmax) {
1286  cmax = cont;
1287  x = fHistMlem->GetXaxis()->GetBinCenter(j);
1288  y = fHistMlem->GetYaxis()->GetBinCenter(i);
1289  i2 = i;
1290  j2 = j;
1291  }
1292  }
1293  }
1294  xq += x*cmax;
1295  yq += y*cmax;
1296  qq += cmax;
1297  if (i2 != i1) nsumy++;
1298  if (j2 != j1) nsumx++;
1299  nsum++;
1300  } // if (nsumx == 1)
1301 
1302  xyc[0] = xq/qq; xyc[1] = yq/qq;
1303  AliDebug(2,Form("x,y COG = %e,%e",xyc[0],xyc[1]));
1304 }
1305 
1306 //_____________________________________________________________________________
1308 {
1311 
1312  Int_t nPix = fPixArray->GetEntriesFast(), imin = 0;
1313  Double_t rmin = 99999, dx = 0, dy = 0, r = 0;
1314  Double_t xc = pixPtr0->Coord(0), yc = pixPtr0->Coord(1);
1315  AliMUONPad *pixPtr;
1316 
1317  for (Int_t i = 0; i < nPix; ++i) {
1318  pixPtr = (AliMUONPad*) fPixArray->UncheckedAt(i);
1319  if (pixPtr == pixPtr0 || pixPtr->Charge() < fLowestPixelCharge) continue;
1320  dx = (xc - pixPtr->Coord(0)) / pixPtr->Size(0);
1321  dy = (yc - pixPtr->Coord(1)) / pixPtr->Size(1);
1322  r = dx *dx + dy * dy;
1323  if (r < rmin) { rmin = r; imin = i; }
1324  }
1325  return imin;
1326 }
1327 
1328 //_____________________________________________________________________________
1329 void
1331 {
1333 
1334  AliMpArea area(fPreCluster->Area());
1335 
1336  gPad->Range(area.LeftBorder(),area.DownBorder(),area.RightBorder(),area.UpBorder());
1337 
1338  gVirtualX->SetFillStyle(1001);
1339  gVirtualX->SetFillColor(3);
1340  gVirtualX->SetLineColor(3);
1341 
1342  Double_t s(1.0);
1343 
1344  for ( Int_t i = 0; i <= fPixArray->GetLast(); ++i)
1345  {
1346  AliMUONPad* pixel = Pixel(i);
1347 
1348  gPad->PaintBox(pixel->Coord(0)-pixel->Size(0)*s,
1349  pixel->Coord(1)-pixel->Size(1)*s,
1350  pixel->Coord(0)+pixel->Size(0)*s,
1351  pixel->Coord(1)+pixel->Size(1)*s);
1352 
1353 // for ( Int_t sign = -1; sign < 2; sign +=2 )
1354 // {
1355 // gPad->PaintLine(pixel->Coord(0) - pixel->Size(0),
1356 // pixel->Coord(1) + sign*pixel->Size(1),
1357 // pixel->Coord(0) + pixel->Size(0),
1358 // pixel->Coord(1) - sign*pixel->Size(1)
1359 // );
1360 // }
1361  }
1362 
1363 
1364  gVirtualX->SetFillStyle(0);
1365 
1366  fPreCluster->Paint();
1367 
1368  gVirtualX->SetLineColor(1);
1369  gVirtualX->SetLineWidth(2);
1370  gVirtualX->SetFillStyle(0);
1371  gVirtualX->SetTextColor(1);
1372  gVirtualX->SetTextAlign(22);
1373 
1374  for ( Int_t i = 0; i <= fPixArray->GetLast(); ++i)
1375  {
1376  AliMUONPad* pixel = Pixel(i);
1377  gPad->PaintBox(pixel->Coord(0)-pixel->Size(0),
1378  pixel->Coord(1)-pixel->Size(1),
1379  pixel->Coord(0)+pixel->Size(0),
1380  pixel->Coord(1)+pixel->Size(1));
1381  gVirtualX->SetTextSize(pixel->Size(0)*60);
1382 
1383  gPad->PaintText(pixel->Coord(0),pixel->Coord(1),Form("%d",(Int_t)(pixel->Charge())));
1384  }
1385 }
1386 
1387 //_____________________________________________________________________________
1388 Int_t AliMUONClusterFinderMLEM::FindLocalMaxima(TObjArray *pixArray, Int_t *localMax, Double_t *maxVal)
1389 {
1393 
1394  AliDebug(1,Form("nPix=%d",pixArray->GetLast()+1));
1395 
1396  Double_t xylim[4] = {999, 999, 999, 999};
1397 
1398  Int_t nPix = pixArray->GetEntriesFast();
1399 
1400  if ( nPix <= 0 ) return 0;
1401 
1402  AliMUONPad *pixPtr = 0;
1403  for (Int_t ipix = 0; ipix < nPix; ++ipix) {
1404  pixPtr = (AliMUONPad*) pixArray->UncheckedAt(ipix);
1405  for (Int_t i = 0; i < 4; ++i)
1406  xylim[i] = TMath::Min (xylim[i], (i%2 ? -1 : 1)*pixPtr->Coord(i/2));
1407  }
1408  for (Int_t i = 0; i < 4; ++i) xylim[i] -= pixPtr->Size(i/2);
1409 
1410  Int_t nx = TMath::Nint ((-xylim[1]-xylim[0])/pixPtr->Size(0)/2);
1411  Int_t ny = TMath::Nint ((-xylim[3]-xylim[2])/pixPtr->Size(1)/2);
1412  if (pixArray == fPixArray) fHistAnode = new TH2D("anode","anode",nx,xylim[0],-xylim[1],ny,xylim[2],-xylim[3]);
1413  else fHistAnode = new TH2D("anode1","anode1",nx,xylim[0],-xylim[1],ny,xylim[2],-xylim[3]);
1414  for (Int_t ipix = 0; ipix < nPix; ++ipix) {
1415  pixPtr = (AliMUONPad*) pixArray->UncheckedAt(ipix);
1416  fHistAnode->Fill(pixPtr->Coord(0), pixPtr->Coord(1), pixPtr->Charge());
1417  }
1418 // if (fDraw && pixArray == fPixArray) fDraw->DrawHist("c2", hist);
1419 
1420  Int_t nMax = 0, indx, nxy = ny * nx;
1421  Int_t *isLocalMax = new Int_t[nxy];
1422  for (Int_t i = 0; i < nxy; ++i) isLocalMax[i] = 0;
1423 
1424  for (Int_t i = 1; i <= ny; ++i) {
1425  indx = (i-1) * nx;
1426  for (Int_t j = 1; j <= nx; ++j) {
1427  if (fHistAnode->GetBinContent(fHistAnode->GetBin(j,i)) < fLowestPixelCharge) continue;
1428  //if (isLocalMax[indx+j-1] < 0) continue;
1429  if (isLocalMax[indx+j-1] != 0) continue;
1430  FlagLocalMax(fHistAnode, i, j, isLocalMax);
1431  }
1432  }
1433 
1434  for (Int_t i = 1; i <= ny; ++i) {
1435  indx = (i-1) * nx;
1436  for (Int_t j = 1; j <= nx; ++j) {
1437  if (isLocalMax[indx+j-1] > 0) {
1438  localMax[nMax] = indx + j - 1;
1439  maxVal[nMax++] = fHistAnode->GetBinContent(fHistAnode->GetBin(j,i));
1440  ((AliMUONPad*)fSplitter->BinToPix(fHistAnode, j, i))->SetStatus(fgkMustKeep);
1441  if (nMax > 99) break;
1442  }
1443  }
1444  if (nMax > 99) {
1445  AliError(" Too many local maxima !!!");
1446  break;
1447  }
1448  }
1449  if (fDebug) cout << " Local max: " << nMax << endl;
1450  delete [] isLocalMax;
1451  return nMax;
1452 }
1453 
1454 //_____________________________________________________________________________
1455 void AliMUONClusterFinderMLEM::FlagLocalMax(TH2D *hist, Int_t i, Int_t j, Int_t *isLocalMax)
1456 {
1458 
1459  Int_t nx = hist->GetNbinsX();
1460  Int_t ny = hist->GetNbinsY();
1461  Int_t cont = TMath::Nint (hist->GetBinContent(hist->GetBin(j,i)));
1462  Int_t cont1 = 0, indx = (i-1)*nx+j-1, indx1 = 0, indx2 = 0;
1463 
1464  Int_t ie = i + 2, je = j + 2;
1465  for (Int_t i1 = i-1; i1 < ie; ++i1) {
1466  if (i1 < 1 || i1 > ny) continue;
1467  indx1 = (i1 - 1) * nx;
1468  for (Int_t j1 = j-1; j1 < je; ++j1) {
1469  if (j1 < 1 || j1 > nx) continue;
1470  if (i == i1 && j == j1) continue;
1471  indx2 = indx1 + j1 - 1;
1472  cont1 = TMath::Nint (hist->GetBinContent(hist->GetBin(j1,i1)));
1473  if (cont < cont1) { isLocalMax[indx] = -1; return; }
1474  else if (cont > cont1) isLocalMax[indx2] = -1;
1475  else { // the same charge
1476  isLocalMax[indx] = 1;
1477  if (isLocalMax[indx2] == 0) {
1478  FlagLocalMax(hist, i1, j1, isLocalMax);
1479  if (isLocalMax[indx2] < 0) { isLocalMax[indx] = -1; return; }
1480  else isLocalMax[indx2] = -1;
1481  }
1482  }
1483  }
1484  }
1485  isLocalMax[indx] = 1; // local maximum
1486 }
1487 
1488 //_____________________________________________________________________________
1490  const Int_t *localMax, Int_t iMax)
1491 {
1494 
1495  /* Just for check
1496  TCanvas* c = new TCanvas("Anode","Anode",800,600);
1497  c->cd();
1498  hist->Draw("lego1Fb"); // debug
1499  c->Update();
1500  Int_t tmp;
1501  cin >> tmp;
1502  */
1503  Int_t nx = fHistAnode->GetNbinsX();
1504  Int_t ny = fHistAnode->GetNbinsY();
1505  Int_t ic = localMax[iMax] / nx + 1;
1506  Int_t jc = localMax[iMax] % nx + 1;
1507  Int_t nxy = ny * nx;
1508  Bool_t *used = new Bool_t[nxy];
1509  for (Int_t i = 0; i < nxy; ++i) used[i] = kFALSE;
1510 
1511  // Drop all pixels from the array - pick up only the ones from the cluster
1512  fPixArray->Delete();
1513 
1514  Double_t wx = fHistAnode->GetXaxis()->GetBinWidth(1)/2;
1515  Double_t wy = fHistAnode->GetYaxis()->GetBinWidth(1)/2;
1516  Double_t yc = fHistAnode->GetYaxis()->GetBinCenter(ic);
1517  Double_t xc = fHistAnode->GetXaxis()->GetBinCenter(jc);
1518  Double_t cont = fHistAnode->GetBinContent( fHistAnode->GetBin(jc,ic));
1519  fPixArray->Add(new AliMUONPad (xc, yc, wx, wy, cont));
1520  used[(ic-1)*nx+jc-1] = kTRUE;
1521  AddBinSimple(fHistAnode, ic, jc);
1522  //fSplitter->AddBin(hist, ic, jc, 1, used, (TObjArray*)0); // recursive call
1523 
1524  Int_t nPix = fPixArray->GetEntriesFast();
1525  Int_t npad = cluster.Multiplicity();
1526 
1527  for (Int_t i = 0; i < nPix; ++i)
1528  {
1529  AliMUONPad* pixPtr = Pixel(i);
1530  pixPtr->SetSize(0,wx);
1531  pixPtr->SetSize(1,wy);
1532  }
1533 
1534  // Pick up pads which overlap with found pixels
1535  for (Int_t i = 0; i < npad; ++i)
1536  {
1537  //cluster.Pad(i)->SetStatus(-1);
1538  cluster.Pad(i)->SetStatus(fgkOver); // just the dirty trick
1539  }
1540 
1541  for (Int_t i = 0; i < nPix; ++i)
1542  {
1543  AliMUONPad* pixPtr = Pixel(i);
1544  for (Int_t j = 0; j < npad; ++j)
1545  {
1546  AliMUONPad* pad = cluster.Pad(j);
1547  //if (pad->Status() == 0) continue;
1548  if (pad->Status() == fgkZero) continue;
1549  if ( Overlap(*pad,*pixPtr) )
1550  {
1551  //pad->SetStatus(0);
1552  pad->SetStatus(fgkZero);
1553  if (fDebug) { cout << j << " "; pad->Print("full"); }
1554  }
1555  }
1556  }
1557 
1558  delete [] used;
1559 }
1560 
1561 //_____________________________________________________________________________
1562 void
1563 AliMUONClusterFinderMLEM::AddBinSimple(TH2D *hist, Int_t ic, Int_t jc)
1564 {
1566 
1567  Int_t nx = hist->GetNbinsX();
1568  Int_t ny = hist->GetNbinsY();
1569  Double_t cont1, cont = hist->GetBinContent(hist->GetBin(jc,ic));
1570  AliMUONPad *pixPtr = 0;
1571 
1572  Int_t ie = TMath::Min(ic+1,ny), je = TMath::Min(jc+1,nx);
1573  for (Int_t i = TMath::Max(ic-1,1); i <= ie; ++i) {
1574  for (Int_t j = TMath::Max(jc-1,1); j <= je; ++j) {
1575  cont1 = hist->GetBinContent(hist->GetBin(j,i));
1576  if (cont1 > cont) continue;
1577  if (cont1 < fLowestPixelCharge) continue;
1578  pixPtr = new AliMUONPad (hist->GetXaxis()->GetBinCenter(j),
1579  hist->GetYaxis()->GetBinCenter(i), 0, 0, cont1);
1580  fPixArray->Add(pixPtr);
1581  }
1582  }
1583 }
1584 
1585 //_____________________________________________________________________________
1588 {
1590 
1591  if (this == &rhs) return *this;
1592 
1593  AliFatal("Not implemented.");
1594 
1595  return *this;
1596 }
1597 
1598 //_____________________________________________________________________________
1600 {
1603 
1604  // Find out non-bending and bending planes
1605  Int_t nonb[2] = {1, 0}; // non-bending and bending cathodes
1606 
1607  TVector2 dim0 = cluster.MinPadDimensions(0, 0, kTRUE);
1608  TVector2 dim1 = cluster.MinPadDimensions(1, 0, kTRUE);
1609  if (dim0.X() < dim1.X() - fgkDistancePrecision) {
1610  nonb[0] = 0;
1611  nonb[1] = 1;
1612  }
1613 
1614  Bool_t same = kFALSE;
1615  if (TMath::Abs(dim0.Y()-dim1.Y()) < fgkDistancePrecision) same = kTRUE; // the same pad size on both planes
1616 
1617  Long_t cn;
1618  Bool_t check[2] = {kFALSE, kFALSE};
1619  Int_t nxy[2];
1620  nxy[0] = nxy[1] = 0;
1621  for (Int_t inb = 0; inb < 2; ++inb) {
1622  cn = cluster.NofPads(nonb[inb], 0, kTRUE);
1623  if (inb == 0 && AliMp::PairFirst(cn) == 2) check[inb] = kTRUE; // check non-bending plane
1624  else if (inb == 1 && AliMp::PairSecond(cn) == 2) check[inb] = kTRUE; // check bending plane
1625  if (same) {
1626  nxy[0] = TMath::Max (nxy[0], AliMp::PairFirst(cn));
1627  nxy[1] = TMath::Max (nxy[1], AliMp::PairSecond(cn));
1628  if (inb == 0 && nxy[0] < 2) nonb[inb] = !nonb[inb];
1629  else if (inb == 1 && AliMp::PairSecond(cn) < 2) nonb[inb] = !nonb[inb];
1630  }
1631  }
1632  if (same) {
1633  if (nxy[0] > 2) check[0] = kFALSE;
1634  if (nxy[1] > 2) check[1] = kFALSE;
1635  }
1636  if (!check[0] && !check[1]) return;
1637 
1638  for (Int_t inb = 0; inb < 2; ++inb) {
1639  if (!check[inb]) continue;
1640 
1641  // Find pads with maximum and next to maximum charges
1642  Int_t maxPads[2] = {-1, -1};
1643  Double_t amax[2] = {0};
1644  Int_t mult = cluster.Multiplicity();
1645  for (Int_t j = 0; j < mult; ++j) {
1646  AliMUONPad *pad = cluster.Pad(j);
1647  if (pad->Cathode() != nonb[inb]) continue;
1648  for (Int_t j2 = 0; j2 < 2; ++j2) {
1649  if (pad->Charge() > amax[j2]) {
1650  if (j2 == 0) { amax[1] = amax[0]; maxPads[1] = maxPads[0]; }
1651  amax[j2] = pad->Charge();
1652  maxPads[j2] = j;
1653  break;
1654  }
1655  }
1656  }
1657 
1658  // Find min and max dimensions of the cluster
1659  Double_t limits[2] = {9999, -9999};
1660  for (Int_t j = 0; j < mult; ++j) {
1661  AliMUONPad *pad = cluster.Pad(j);
1662  if (pad->Cathode() != nonb[inb]) continue;
1663  if (pad->Coord(inb) < limits[0]) limits[0] = pad->Coord(inb);
1664  if (pad->Coord(inb) > limits[1]) limits[1] = pad->Coord(inb);
1665  }
1666 
1667  // Loop over max and next to max pads
1668  Bool_t add = kFALSE;
1669  Int_t idirAdd = 0;
1670  for (Int_t j = 0; j < 2; ++j) {
1671  if (j == 1) {
1672  if (maxPads[j] < 0) continue;
1673  if (!add) break;
1674  if (amax[1] / amax[0] < 0.5) break;
1675  }
1676  // Check if pad at the cluster limit
1677  AliMUONPad *pad = cluster.Pad(maxPads[j]);
1678  Int_t idir = 0;
1679  if (TMath::Abs(pad->Coord(inb)-limits[0]) < fgkDistancePrecision) idir = -1;
1680  else if (TMath::Abs(pad->Coord(inb)-limits[1]) < fgkDistancePrecision) idir = 1;
1681  else {
1682  //cout << " *** Pad not at the cluster limit: " << j << endl;
1683  break;
1684  }
1685  if (j == 1 && idir == idirAdd) break; // this pad has been already added
1686 
1687  // Add pad (if it exists)
1688  TVector2 pos;
1689  if (inb == 0) pos.Set (pad->X() + idir * (pad->DX()+fgkDistancePrecision), pad->Y());
1690  else pos.Set (pad->X(), pad->Y() + idir * (pad->DY()+fgkDistancePrecision));
1691  //AliMpPad mppad = fkSegmentation[nonb[inb]]->PadByPosition(pos,kTRUE);
1692  AliMpPad mppad = fkSegmentation[nonb[inb]]->PadByPosition(pos.X(), pos.Y(),kFALSE);
1693  if (!mppad.IsValid()) continue; // non-existing pad
1694  AliMUONPad muonPad(fDetElemId, nonb[inb], mppad.GetIx(), mppad.GetIy(),
1695  mppad.GetPositionX(), mppad.GetPositionY(),
1696  mppad.GetDimensionX(), mppad.GetDimensionY(), 0);
1697  if (inb == 0) muonPad.SetCharge(TMath::Min (amax[j]/100, fLowestPadCharge));
1698  //else muonPad.SetCharge(TMath::Min (amax[j]/15, fgkZeroSuppression));
1699  else muonPad.SetCharge(TMath::Min (amax[j]/15, fLowestPadCharge));
1700  if (muonPad.Charge() < 2.0*fLowestPixelCharge) muonPad.SetCharge(2.0*fLowestPixelCharge);
1701  muonPad.SetReal(kFALSE);
1702  if (fDebug) printf(" ***** Add virtual pad in %d direction ***** %f %f %f %3d %3d %f %f \n",
1703  inb, muonPad.Charge(), muonPad.X(), muonPad.Y(), muonPad.Ix(),
1704  muonPad.Iy(), muonPad.DX(), muonPad.DY());
1705  cluster.AddPad(muonPad); // add pad to the cluster
1706  add = kTRUE;
1707  idirAdd = idir;
1708  }
1709  }
1710 }
1711 
1712 //_____________________________________________________________________________
1714  Int_t &nInX, Int_t &nInY) const
1715 {
1718 
1719  //Int_t statusToTest = 1;
1720  Int_t statusToTest = fgkUseForFit;
1721 
1722  //if ( nInX < 0 ) statusToTest = 0;
1723  if ( nInX < 0 ) statusToTest = fgkZero;
1724 
1725  Bool_t mustMatch(kTRUE);
1726 
1727  Long_t cn = cluster.NofPads(statusToTest,mustMatch);
1728 
1729  nInX = AliMp::PairFirst(cn);
1730  nInY = AliMp::PairSecond(cn);
1731 }
1732 
1733 //_____________________________________________________________________________
1735 {
1737  AliMUONPad* pixPtr = Pixel(i);
1738  fPixArray->RemoveAt(i);
1739  delete pixPtr;
1740 }
1741 
1742 //_____________________________________________________________________________
1744 {
1746 
1747  Int_t nForFit = 1, clustFit[1] = {0};
1748  Double_t parOk[3] = {0.};
1749  TObjArray *clusters[1];
1750  clusters[0] = fPixArray;
1751 
1752  AliDebug(1,Form("nPix=%d",fPixArray->GetLast()+1));
1753 
1754  Int_t mult = cluster.Multiplicity();
1755  for (Int_t i = 0; i < mult; ++i)
1756  {
1757  AliMUONPad* pad = cluster.Pad(i);
1758  /*
1759  if ( pad->IsSaturated())
1760  {
1761  pad->SetStatus(-9);
1762  }
1763  else
1764  {
1765  pad->SetStatus(1);
1766  }
1767  */
1768  if (!pad->IsSaturated()) pad->SetStatus(fgkUseForFit);
1769  }
1770  fSplitter->Fit(cluster,1, nForFit, clustFit, clusters, parOk, fClusterList, fHistMlem);
1771 }
1772 
1773 //_____________________________________________________________________________
1774 AliMUONPad*
1776 {
1778  return static_cast<AliMUONPad*>(fPixArray->UncheckedAt(i));
1779 }
1780 
1781 //_____________________________________________________________________________
1782 void
1783 AliMUONClusterFinderMLEM::Print(Option_t* what) const
1784 {
1786  TString swhat(what);
1787  swhat.ToLower();
1788  if ( swhat.Contains("precluster") )
1789  {
1790  if ( fPreCluster) fPreCluster->Print();
1791  }
1792 }
1793 
1794 //_____________________________________________________________________________
1795 void
1796 AliMUONClusterFinderMLEM::SetChargeHints(Double_t lowestPadCharge, Double_t lowestClusterCharge)
1797 {
1799  fLowestPadCharge=lowestPadCharge;
1800  fLowestClusterCharge=lowestClusterCharge;
1802 }
1803 
1804 
static Bool_t AreOverlapping(const AliMUONPad &d1, const AliMUONPad &d2, const TVector2 &precision, AliMpArea &overlapArea)
Definition: AliMUONPad.cxx:221
AliMUONVClusterFinder * fPreClusterFinder
! the pre-clustering worker
Int_t FindNearest(const AliMUONPad *pixPtr0)
Int_t fNClusters
! total number of clusters
printf("Chi2/npoints = %f\n", TMath::Sqrt(chi2/npoints))
void BuildPixArrayOneCathode(AliMUONCluster &cluster)
Bool_t IsReal() const
Return info whether this is a real pad or a virtual one.
Definition: AliMUONPad.h:64
Interface of a cluster finder.
void SetStatus(Int_t status)
Set status word.
Definition: AliMUONPad.h:104
Double_t fLowestPixelCharge
! see AliMUONRecoParam
Bool_t IsSaturated() const
Return info whether this pad is saturated or not.
Definition: AliMUONPad.h:67
AliMUONPad * Pixel(Int_t i) const
#define TObjArray
Int_t GetIy() const
Definition: AliMpPad.cxx:280
TObject * BinToPix(TH2 *mlem, Int_t jc, Int_t ic)
Int_t FindLocalMaxima(TObjArray *pixArray, Int_t *localMax, Double_t *maxVal)
A group of adjacent pads.
void RevertCharge()
Set charge to value in backup charge.
Definition: AliMUONPad.h:93
virtual void Print(Option_t *opt="") const
void Split(const AliMUONCluster &cluster, TH2 *mlem, Double_t *coef, TObjArray &clusterList)
Splitter class for the MLEM algorithm.
A rectangle area positioned in plane..
Definition: AliMpArea.h:20
void PadsInXandY(AliMUONCluster &cluster, Int_t &nInX, Int_t &nInY) const
virtual Bool_t Prepare(Int_t detElemId, TObjArray *pads[2], const AliMpArea &area, const AliMpVSegmentation *segmentations[2])
Double_t DY() const
Return half dimensions in y (cm)
Definition: AliMUONPad.h:61
Bool_t MainLoop(AliMUONCluster &cluster, Int_t iSimple)
Int_t Multiplicity() const
virtual Bool_t Prepare(Int_t detElemId, TObjArray *pads[2], const AliMpArea &area)
void Shift(Int_t ixy, Double_t shift)
Definition: AliMUONPad.cxx:500
virtual Bool_t NeedSegmentation() const
TH2D * fHistMlem
! histogram for MLEM procedure
void AddBinSimple(TH2D *mlem, Int_t ic, Int_t jc)
void BuildPixArray(AliMUONCluster &cluster)
build array of pixels
void SetSize(Int_t ixy, Double_t Size)
Definition: AliMUONPad.cxx:481
AliMUONCluster * fPreCluster
! current pre-cluster
virtual Bool_t UsePad(const AliMUONPad &pad)
ClassImp(TPCGenInfo)
Definition: AliTPCCmpNG.C:254
AliMUONCluster * CheckPrecluster(const AliMUONCluster &cluster)
Check precluster to simplify it (if possible), and return the simplified cluster. ...
static const Int_t fgkMustKeep
do not kill (for pixels)
void ComputeCoefficients(AliMUONCluster &cluster, Double_t *coef, Double_t *probi)
void sum()
Double_t Y() const
Return position in y (cm)
Definition: AliMUONPad.h:128
void SetChargeBackup(Double_t charge)
Set charge backup.
Definition: AliMUONPad.h:99
Int_t fClusterNumber
! current cluster number
Bool_t Overlap(const AliMUONPad &pad, const AliMUONPad &pixel)
Checks whether a pad and a pixel have an overlapping area.
void FlagLocalMax(TH2D *hist, Int_t i, Int_t j, Int_t *isLocalMax)
void SetCharge(Double_t charge)
Set charge.
Definition: AliMUONPad.h:96
return clusterFinder
static const TVector2 fgkDecreaseSize
idem
Int_t Iy() const
Return y-index.
Definition: AliMUONPad.h:80
void Release()
Detach this pad from a cluster.
Definition: AliMUONPad.h:90
static const Int_t fgkZero
pad "basic" state
TObjArray fClusterList
! clusters corresponding to the current pre-cluster
void Simple(AliMUONCluster &cluster)
Process simple cluster.
void LeftDownCorner(Double_t &x, Double_t &y) const
Definition: AliMpArea.cxx:156
Double_t X() const
Return position in x (cm)
Definition: AliMUONPad.h:126
Double_t GetPositionY() const
Return the pad x position (in cm)
Definition: AliMpPad.h:81
void SetDebug(Int_t debug)
Set debug level.
virtual void Print(Option_t *opt="") const
void AddVirtualPad(AliMUONCluster &cluster)
Int_t Cathode() const
Return cathode number.
Definition: AliMUONPad.h:43
AliMUONClusterFinderMLEM & operator=(const AliMUONClusterFinderMLEM &rhs)
Not implemented.
Bool_t IsUsed() const
Return true if is used.
Definition: AliMUONPad.h:75
AliMUONCluster * CheckPreclusterTwoCathodes(AliMUONCluster *cluster)
Double_t Size(Int_t ixy) const
Definition: AliMUONPad.cxx:508
void RemovePad(AliMUONPad *pad)
Int_t fDetElemId
! current DE being processed
AliMUONPad * AddPad(const AliMUONPad &pad)
Int_t fEventNumber
! current event being processed
TObjArray * fPixArray
! collection of pixels
Bool_t IsValid() const
Return validity.
Definition: AliMpPad.h:89
TVector2 Position() const
Return positions in x and y (cm)
Definition: AliMUONPad.h:85
Float_t Charge() const
virtual AliMpPad PadByPosition(Double_t x, Double_t y, Bool_t warning=true) const =0
Find pad by position.
Int_t GetIx() const
Definition: AliMpPad.cxx:272
static const Int_t fgkOver
processing is over
static const Int_t fgkUseForFit
should be used for fit
void RightUpCorner(Double_t &x, Double_t &y) const
Definition: AliMpArea.cxx:184
virtual AliMUONCluster * NextCluster()=0
void Plot(const char *outputfile)
void SetCoord(Int_t ixy, Double_t Coord)
Definition: AliMUONPad.cxx:462
void Mlem(AliMUONCluster &cluster, const Double_t *coef, Double_t *probi, Int_t nIter)
Long_t NofPads(Int_t cathode, Int_t statusMask, Bool_t matchMask) const
Compute number of pads in X and Y direction for a given cathode.
The abstract base class for the segmentation.
Int_t Ix() const
Return x-index.
Definition: AliMUONPad.h:78
AliMpArea Area() const
Area that contains all the pads (whatever the cathode)
Double_t GetDimensionY() const
Return the y pad dimension - half length (in cm)
Definition: AliMpPad.h:86
virtual void Paint(Option_t *opt="")
Int_t PairFirst(MpPair_t pair)
Decode the first integer from encoded pair.
void Print(Option_t *opt="") const
Definition: AliMUONPad.cxx:420
Double_t Charge() const
Return pad charge.
Definition: AliMUONPad.h:48
Int_t Fit(const AliMUONCluster &cluster, Int_t iSimple, Int_t nfit, const Int_t *clustFit, TObjArray **clusters, Double_t *parOk, TObjArray &clusterList, TH2 *mlem)
virtual AliMUONCluster * NextCluster()
Double_t Coord(Int_t ixy) const
Definition: AliMUONPad.cxx:354
Class which encapsuate all information about a pad.
Definition: AliMpPad.h:22
Int_t fNAddVirtualPads
! number of clusters for which we added virtual pads
Bool_t IsSaturated(Int_t cathode) const
Whether we have at least one saturated pad in a given cathode.
Double_t fLowestClusterCharge
! see AliMUONRecoParam
void FindCluster(AliMUONCluster &cluster, const Int_t *localMax, Int_t iMax)
Double_t GetPositionX() const
Return the pad x position (in cm)
Definition: AliMpPad.h:79
Int_t DetElemId() const
Return detection element id.
Definition: AliMUONPad.h:53
AliMUONClusterSplitterMLEM * fSplitter
! helper class to go from pixel arrays to clusters
Float_t ChargeAsymmetry() const
Return the cathode's charges asymmetry.
Float_t ChargeIntegration(Double_t x, Double_t y, const AliMUONPad &pad)
virtual void Paint(Option_t *opt="")
Int_t PairSecond(MpPair_t pair)
Decode the second integer from encoded pair.
virtual void SetChargeHints(Double_t lowestPadCharge, Double_t lowestClusterCharge)
TVector2 MinPadDimensions(Int_t cathode, Int_t statusMask, Bool_t matchMask) const
Return the smallest pad dimensions for a given cathode.
Int_t Status() const
Return status word.
Definition: AliMUONPad.h:123
TH2D * fHistAnode
! histogram for local maxima search
Double_t fLowestPadCharge
! see AliMUONRecoParam
static const Double_t fgkDistancePrecision
used to check overlaps and so on
Int_t MaxRawChargeCathode() const
Return the max raw charge on the chathod.
AliMUONPad * Pad(Int_t index) const
Double_t GetDimensionX() const
Return the x pad dimension - half length (in cm)
Definition: AliMpPad.h:84
Cluster finder in MUON arm of ALICE.
Combination of digit and mppad informations.
Definition: AliMUONPad.h:25
void PadOverHist(Int_t idir, Int_t ix0, Int_t iy0, AliMUONPad *pad, TH2D *hist1, TH2D *hist2)
const AliMpVSegmentation * fkSegmentation[2]
! new segmentation
Double_t DX() const
Return half dimensions in x (cm)
Definition: AliMUONPad.h:59