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