AliRoot Core  v5-06-30 (35d6c57)
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
AliMUONClusterFinderPeakFit.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 //-----------------------------------------------------------------------------
30 //-----------------------------------------------------------------------------
31 
33 #include "AliMUONCluster.h"
34 #include "AliMUONConstants.h"
35 #include "AliMUONPad.h"
36 #include "AliMUONMathieson.h"
37 
38 #include "AliMpDEManager.h"
39 #include "AliMpPad.h"
40 #include "AliMpVSegmentation.h"
41 #include "AliMpEncodePair.h"
42 
43 #include "AliLog.h"
44 #include "AliRunLoader.h"
45 //#include "AliCodeTimer.h"
46 
47 #include <Riostream.h>
48 #include <TH2.h>
49 #include <TVirtualFitter.h>
50 #include <TMath.h>
51 //#include <TCanvas.h>
52 
53 using std::endl;
54 using std::cout;
58 
59 const Double_t AliMUONClusterFinderPeakFit::fgkZeroSuppression = 6; // average zero suppression value
60 //const Double_t AliMUONClusterFinderMLEM::fgkDistancePrecision = 1e-6; // (cm) used to check overlaps and so on
61 const Double_t AliMUONClusterFinderPeakFit::fgkDistancePrecision = 1e-3; // (cm) used to check overlaps and so on
62 const TVector2 AliMUONClusterFinderPeakFit::fgkIncreaseSize(-AliMUONClusterFinderPeakFit::fgkDistancePrecision,-AliMUONClusterFinderPeakFit::fgkDistancePrecision);
63 const TVector2 AliMUONClusterFinderPeakFit::fgkDecreaseSize(AliMUONClusterFinderPeakFit::fgkDistancePrecision,AliMUONClusterFinderPeakFit::fgkDistancePrecision);
64 
65 // Status flags for pads
66 const Int_t AliMUONClusterFinderPeakFit::fgkZero = 0x0;
67 const Int_t AliMUONClusterFinderPeakFit::fgkMustKeep = 0x1;
68 const Int_t AliMUONClusterFinderPeakFit::fgkUseForFit = 0x10;
69 const Int_t AliMUONClusterFinderPeakFit::fgkOver = 0x100;
70 const Int_t AliMUONClusterFinderPeakFit::fgkModified = 0x1000;
71 const Int_t AliMUONClusterFinderPeakFit::fgkCoupled = 0x10000;
72 
73 namespace
74 {
75  //_____________________________________________________________________________
76  Double_t Param2Coef(Int_t icand, Double_t coef, Double_t *par, Int_t nHits)
77  {
79 
80  //Int_t nHits = TMath::Nint(par[8]);
81  if (nHits == 1) return 1.;
82  if (nHits == 2) return icand==0 ? par[2] : TMath::Max(1.-par[2],0.);
83  if (icand == 0) return par[2];
84  if (icand == 1) return TMath::Max((1.-par[2])*par[5], 0.);
85  return TMath::Max(1.-par[2]-coef,0.);
86  }
87 
88  //___________________________________________________________________________
89  void
90  FitFunction(Int_t& /*notused*/, Double_t* /*notused*/,
91  Double_t& f, Double_t* par,
92  Int_t /*notused*/)
93  {
95 
96  TObjArray* userObjects = static_cast<TObjArray*>(TVirtualFitter::GetFitter()->GetObjectFit());
97 
98  AliMUONCluster* cluster = static_cast<AliMUONCluster*>(userObjects->At(0));
99  AliMUONMathieson* mathieson = static_cast<AliMUONMathieson*>(userObjects->At(1));
100  AliMUONClusterFinderPeakFit* finder =
101  static_cast<AliMUONClusterFinderPeakFit*>(userObjects->At(2));
102 
103  f = 0.0;
104  Int_t nHits = finder->GetNMax(), npads = cluster->Multiplicity();
105  Double_t qTot = cluster->Charge(), coef = 0;
106  //if (cluster->Multiplicity(0) == 0 || cluster->Multiplicity(1) == 0) qTot *= 2.;
107 
108  for ( Int_t i = 0 ; i < npads; ++i )
109  {
110  AliMUONPad* pad = cluster->Pad(i);
111  // skip pads w/ saturation or other problem(s)
112  //if ( pad->Status() ) continue;
113  if ( pad->IsSaturated() ) continue;
114  Double_t charge = 0.;
115  for (Int_t j = 0; j < nHits; ++j) {
116  // Sum over hits
117  Int_t indx = 3 * j;
118  TVector2 lowerLeft = TVector2(par[indx],par[indx+1]) - pad->Position() - pad->Dimensions();
119  TVector2 upperRight(lowerLeft + pad->Dimensions()*2.0);
120  Double_t estimatedCharge = mathieson->IntXY(lowerLeft.X(),lowerLeft.Y(),
121  upperRight.X(),upperRight.Y());
122  coef = Param2Coef(j, coef, par, nHits);
123  charge += estimatedCharge * coef;
124  }
125  charge *= qTot;
126 
127  Double_t delta = charge - pad->Charge();
128  delta *= delta;
129  delta /= pad->Charge();
130  f += delta;
131  }
132  f /= (qTot/npads);
133  //cout << qTot << " " << par[0] << " " << par[1] << " " << f << endl;
134  }
135 }
136 
137 //_____________________________________________________________________________
140 fPreClusterFinder(clusterFinder),
141 fPreCluster(0x0),
142 fClusterList(),
143 fMathieson(0x0),
144 fEventNumber(0),
145 fDetElemId(-1),
146 fClusterNumber(0),
147 fNMax(0),
148 fHistAnode(0x0),
149 fPixArray(new TObjArray(20)),
150 fDebug(0),
151 fPlot(plot),
152 fNClusters(0),
153 fNAddVirtualPads(0)
154 {
156 
157  fkSegmentation[1] = fkSegmentation[0] = 0x0;
158 
159  if (fPlot) fDebug = 1;
160 }
161 
162 //_____________________________________________________________________________
164 {
166  delete fPixArray; fPixArray = 0;
167  delete fPreClusterFinder;
168  AliInfo(Form("Total clusters %d AddVirtualPad needed %d",
170  delete fMathieson;
171 }
172 
173 //_____________________________________________________________________________
174 Bool_t
176  const AliMpArea& area,
177  const AliMpVSegmentation* seg[2])
178 {
180 // AliCodeTimerAuto("",0)
181 
182  for ( Int_t i = 0; i < 2; ++i )
183  {
184  fkSegmentation[i] = seg[i];
185  }
186 
187  // Find out the DetElemId
188  fDetElemId = detElemId;
189 
190  // find out current event number, and reset the cluster number
191  AliRunLoader *runLoader = AliRunLoader::Instance();
192  fEventNumber = runLoader ? runLoader->GetEventNumber() : 0;
193  fClusterNumber = -1;
194  fClusterList.Delete();
195 
196  AliDebug(3,Form("EVT %d DE %d",fEventNumber,fDetElemId));
197 
199 
200  Float_t kx3 = AliMUONConstants::SqrtKx3();
201  Float_t ky3 = AliMUONConstants::SqrtKy3();
202  Float_t pitch = AliMUONConstants::Pitch();
203 
204  if ( stationType == AliMq::kStation1 )
205  {
208  pitch = AliMUONConstants::PitchSt1();
209  }
210 
211  delete fMathieson;
213 
214  fMathieson->SetPitch(pitch);
217 
219  {
220  return fPreClusterFinder->Prepare(detElemId,pads,area,seg);
221  }
222  else
223  {
224  return fPreClusterFinder->Prepare(detElemId,pads,area);
225  }
226 }
227 
228 //_____________________________________________________________________________
231 {
233 // AliCodeTimerAuto("",0)
234 
235  // if the list of clusters is not void, pick one from there
236  TObject* o = fClusterList.At(++fClusterNumber);
237  if ( o != 0x0 ) return static_cast<AliMUONCluster*>(o);
238 
239  //FIXME : at this point, must check whether we've used all the digits
240  //from precluster : if not, let the preclustering know about those unused
241  //digits, so it can reuse them
242 
243  // if the cluster list is exhausted, we need to go to the next
244  // pre-cluster and treat it
245 
246  fClusterList.Delete(); // reset the list of clusters for this pre-cluster
247  fClusterNumber = -1;
248 
250 
251  if (!fPreCluster)
252  {
253  // we are done
254  return 0x0;
255  }
256 
258 
259  // WorkOnPreCluster may have used only part of the pads, so we check that
260  // now, and let the unused pads be reused by the preclustering...
261 
262  Int_t mult = fPreCluster->Multiplicity();
263  for ( Int_t i = 0; i < mult; ++i )
264  {
265  AliMUONPad* pad = fPreCluster->Pad(i);
266  if ( !pad->IsUsed() )
267  {
268  fPreClusterFinder->UsePad(*pad);
269  }
270  }
271 
272  return NextCluster();
273 }
274 
275 //_____________________________________________________________________________
276 Bool_t
278 {
281 
282  // AliCodeTimerAuto("",0)
283 
284  if (fDebug) {
285  cout << " *** Event # " << fEventNumber
286  << " det. elem.: " << fDetElemId << endl;
287  for (Int_t j = 0; j < fPreCluster->Multiplicity(); ++j) {
288  AliMUONPad* pad = fPreCluster->Pad(j);
289  printf(" bbb %3d %1d %8.4f %8.4f %8.4f %8.4f %6.1f %3d %3d %2d %1d %1d \n",
290  j, pad->Cathode(), pad->Coord(0), pad->Coord(1), pad->DX()*2, pad->DY()*2,
291  pad->Charge(), pad->Ix(), pad->Iy(), pad->Status(), pad->IsReal(), pad->IsSaturated());
292  }
293  }
294 
296  if (!cluster) return kFALSE;
297 
298  BuildPixArray(*cluster);
299 
300  if ( fPixArray->GetLast() < 0 )
301  {
302  AliDebug(1,"No pixel for the above cluster");
303  delete cluster;
304  return kFALSE;
305  }
306 
307  Int_t nMax = 1, localMax[100], maxPos[100] = {0};
308  Double_t maxVal[100];
309 
310  nMax = FindLocalMaxima(fPixArray, localMax, maxVal); // find local maxima
311 
312  if (nMax > 1) TMath::Sort(nMax, maxVal, maxPos, kTRUE); // in descending order
313 
314  if (nMax <= 3) {
315  FindClusterFit(*cluster, localMax, maxPos, nMax);
316  } else {
317  for (Int_t i = 0; i < nMax; ++i)
318  {
319  FindClusterCOG(*cluster, localMax, maxPos[i]);
320  }
321  }
322 
323  delete cluster;
324  if (fPlot == 0) {
325  delete fHistAnode;
326  fHistAnode = 0x0;
327  }
328  return kTRUE;
329 }
330 
331 //_____________________________________________________________________________
332 Bool_t
334 {
336 
337  // make a fake pad from the pixel
338  AliMUONPad tmp(pad.DetElemId(),pad.Cathode(),pad.Ix(),pad.Iy(),
339  pixel.Coord(0),pixel.Coord(1),
340  pixel.Size(0),pixel.Size(1),0);
341 
343 }
344 
345 //_____________________________________________________________________________
348 {
351 
352  // AliCodeTimerAuto("",0)
353 
354  // Disregard small clusters (leftovers from splitting or noise)
355  if ((origCluster.Multiplicity()==1 || origCluster.Multiplicity()==2) &&
356  origCluster.Charge(0)+origCluster.Charge(1) < 1.525) // JC: adc -> fc
357  {
358  return 0x0;
359  }
360 
361  AliMUONCluster* cluster = new AliMUONCluster(origCluster);
362 
363  AliDebug(2,"Start of CheckPreCluster=");
364  //StdoutToAliDebug(2,cluster->Print("full"));
365 
366  AliMUONCluster* rv(0x0);
367 
368  if (cluster->Multiplicity(0) && cluster->Multiplicity(1))
369  {
370  rv = CheckPreclusterTwoCathodes(cluster);
371  }
372  else
373  {
374  rv = cluster;
375  }
376  return rv;
377 }
378 
379 //_____________________________________________________________________________
382 {
384 
385  Int_t npad = cluster->Multiplicity();
386  Int_t* flags = new Int_t[npad];
387  for (Int_t j = 0; j < npad; ++j) flags[j] = 0;
388 
389  // Check pad overlaps
390  for ( Int_t i = 0; i < npad; ++i)
391  {
392  AliMUONPad* padi = cluster->Pad(i);
393  if ( padi->Cathode() != 0 ) continue;
394  for (Int_t j = i+1; j < npad; ++j)
395  {
396  AliMUONPad* padj = cluster->Pad(j);
397  if ( padj->Cathode() != 1 ) continue;
398  if ( !AliMUONPad::AreOverlapping(*padi,*padj,fgkDecreaseSize) ) continue;
399  flags[i] = flags[j] = 1; // mark overlapped pads
400  }
401  }
402 
403  // Check if all pads overlap
404  Int_t nFlags=0;
405  for (Int_t i = 0; i < npad; ++i)
406  {
407  if (!flags[i]) ++nFlags;
408  }
409 
410  if (nFlags > 0)
411  {
412  // not all pads overlap.
413  if (fDebug) cout << " nFlags: " << nFlags << endl;
414  TObjArray toBeRemoved;
415  for (Int_t i = 0; i < npad; ++i)
416  {
417  AliMUONPad* pad = cluster->Pad(i);
418  if (flags[i]) continue;
419  Int_t cath = pad->Cathode();
420  Int_t cath1 = TMath::Even(cath);
421  // Check for edge effect (missing pads on the _other_ cathode)
422  AliMpPad mpPad
423  = fkSegmentation[cath1]->PadByPosition(pad->Position().X(),pad->Position().Y(),kFALSE);
424  if (!mpPad.IsValid()) continue;
425  //if (nFlags == 1 && pad->Charge() < fgkZeroSuppression * 3) continue;
426  if (nFlags == 1 && pad->Charge() < 3.05) continue; // JC: adc -> fc
427  AliDebug(2,Form("Releasing the following pad : de,cath,ix,iy %d,%d,%d,%d charge %e",
428  fDetElemId,pad->Cathode(),pad->Ix(),pad->Iy(),pad->Charge()));
429  toBeRemoved.AddLast(pad);
430  fPreCluster->Pad(i)->Release();
431  }
432  Int_t nRemove = toBeRemoved.GetEntriesFast();
433  for ( Int_t i = 0; i < nRemove; ++i )
434  {
435  cluster->RemovePad(static_cast<AliMUONPad*>(toBeRemoved.UncheckedAt(i)));
436  }
437  }
438 
439  // Check correlations of cathode charges
440  if ( !cluster->IsSaturated() && cluster->ChargeAsymmetry() > 1 )
441  {
442  // big difference
443  Int_t cathode = cluster->MaxRawChargeCathode();
444  Int_t imin(-1);
445  Int_t imax(-1);
446  Double_t cmax(0);
447  Double_t cmin(1E9);
448 
449  // get min and max pad charges on the cathode opposite to the
450  // max pad (given by MaxRawChargeCathode())
451  //
452  Int_t mult = cluster->Multiplicity();
453  for ( Int_t i = 0; i < mult; ++i )
454  {
455  AliMUONPad* pad = cluster->Pad(i);
456  if ( pad->Cathode() != cathode || !pad->IsReal() )
457  {
458  // only consider pads in the opposite cathode, and
459  // only consider real pads (i.e. exclude the virtual ones)
460  continue;
461  }
462  if ( pad->Charge() < cmin )
463  {
464  cmin = pad->Charge();
465  imin = i;
466  if (imax < 0) {
467  imax = imin;
468  cmax = cmin;
469  }
470  }
471  else if ( pad->Charge() > cmax )
472  {
473  cmax = pad->Charge();
474  imax = i;
475  }
476  }
477  AliDebug(2,Form("Pad imin,imax %d,%d cmin,cmax %e,%e",
478  imin,imax,cmin,cmax));
479  //
480  // arrange pads according to their distance to the max, normalized
481  // to the pad size
482  Double_t* dist = new Double_t[mult];
483  Double_t dxMin(1E9);
484  Double_t dyMin(1E9);
485  Double_t dmin(0);
486 
487  AliMUONPad* padmax = cluster->Pad(imax);
488 
489  for ( Int_t i = 0; i < mult; ++i )
490  {
491  dist[i] = 0.0;
492  if ( i == imax) continue;
493  AliMUONPad* pad = cluster->Pad(i);
494  if ( pad->Cathode() != cathode || !pad->IsReal() ) continue;
495  Double_t dx = (pad->X()-padmax->X())/padmax->DX()/2.0;
496  Double_t dy = (pad->Y()-padmax->Y())/padmax->DY()/2.0;
497  dist[i] = TMath::Sqrt(dx*dx+dy*dy);
498  if ( i == imin )
499  {
500  dmin = dist[i] + 1E-3; // distance to the pad with minimum charge
501  dxMin = dx;
502  dyMin = dy;
503  }
504  }
505 
506  TMath::Sort(mult,dist,flags,kFALSE); // in ascending order
507  Double_t xmax(-1), distPrev(999);
508  TObjArray toBeRemoved;
509 
510  for ( Int_t i = 0; i < mult; ++i )
511  {
512  Int_t indx = flags[i];
513  AliMUONPad* pad = cluster->Pad(indx);
514  if ( pad->Cathode() != cathode || !pad->IsReal() ) continue;
515  if ( dist[indx] > dmin )
516  {
517  // farther than the minimum pad
518  Double_t dx = (pad->X()-padmax->X())/padmax->DX()/2.0;
519  Double_t dy = (pad->Y()-padmax->Y())/padmax->DY()/2.0;
520  dx *= dxMin;
521  dy *= dyMin;
522  if (dx >= 0 && dy >= 0) continue;
523  if (TMath::Abs(dx) > TMath::Abs(dy) && dx >= 0) continue;
524  if (TMath::Abs(dy) > TMath::Abs(dx) && dy >= 0) continue;
525  }
526  if (dist[indx] > distPrev + 1) break; // overstepping empty pads
527  if ( pad->Charge() <= cmax || TMath::Abs(dist[indx]-xmax) < 1E-3 )
528  {
529  // release pad
530  if (TMath::Abs(dist[indx]-xmax) < 1.e-3)
531  {
532  cmax = TMath::Max(pad->Charge(),cmax);
533  }
534  else
535  {
536  cmax = pad->Charge();
537  }
538  xmax = dist[indx];
539  distPrev = dist[indx];
540  AliDebug(2,Form("Releasing the following pad : de,cath,ix,iy %d,%d,%d,%d charge %e",
541  fDetElemId,pad->Cathode(),pad->Ix(),pad->Iy(),
542  pad->Charge()));
543 
544  toBeRemoved.AddLast(pad);
545  fPreCluster->Pad(indx)->Release();
546  }
547  }
548  Int_t nRemove = toBeRemoved.GetEntriesFast();
549  for ( Int_t i = 0; i < nRemove; ++i )
550  {
551  cluster->RemovePad(static_cast<AliMUONPad*>(toBeRemoved.UncheckedAt(i)));
552  }
553  delete[] dist;
554  } // if ( !cluster->IsSaturated() &&
555 
556  delete[] flags;
557 
558  AliDebug(2,"End of CheckPreClusterTwoCathodes=");
559  //StdoutToAliDebug(2,cluster->Print("full"));
560 
561  return cluster;
562 }
563 
564 //_____________________________________________________________________________
565 void
567 {
569 
570  Int_t nPix = fPixArray->GetLast()+1;
571  Int_t dummy(0);
572 
573  for ( Int_t i = 0; i < nPix; ++i )
574  {
575  AliMUONPad* pixelI = Pixel(i);
576  AliMUONPad pi(dummy,dummy,dummy,dummy,
577  pixelI->Coord(0),pixelI->Coord(1),
578  pixelI->Size(0),pixelI->Size(1),0.0);
579 
580  for ( Int_t j = i+1; j < nPix; ++j )
581  {
582  AliMUONPad* pixelJ = Pixel(j);
583  AliMUONPad pj(dummy,dummy,dummy,dummy,
584  pixelJ->Coord(0),pixelJ->Coord(1),
585  pixelJ->Size(0),pixelJ->Size(1),0.0);
586  AliMpArea area;
587 
589  {
590  AliInfo(Form("The following 2 pixels (%d and %d) overlap !",i,j));
591  /*
592  StdoutToAliInfo(pixelI->Print();
593  cout << " Surface = " << pixelI->Size(0)*pixelI->Size(1)*4 << endl;
594  pixelJ->Print();
595  cout << " Surface = " << pixelJ->Size(0)*pixelJ->Size(1)*4 << endl;
596  cout << " Area surface = " << area.Dimensions().X()*area.Dimensions().Y()*4 << endl;
597  cout << "-------" << endl;
598  );
599  */
600  }
601  }
602  }
603 }
604 
605 //_____________________________________________________________________________
607 {
609 
610  Int_t npad = cluster.Multiplicity();
611  if (npad<=0)
612  {
613  AliWarning("Got no pad at all ?!");
614  }
615 
616  fPixArray->Delete();
617  BuildPixArrayOneCathode(cluster);
618 
619 // StdoutToAliDebug(2,cout << "End of BuildPixelArray:" << endl;
620 // fPixArray->Print(););
621  //CheckOverlaps();//FIXME : this is for debug only. Remove it.
622 }
623 
624 //_____________________________________________________________________________
626 {
628 
629 // AliDebug(2,Form("cluster.Multiplicity=%d",cluster.Multiplicity()));
630 
631  TVector2 dim = cluster.MinPadDimensions (-1, kFALSE);
632  Double_t width[2] = {dim.X(), dim.Y()}, xy0[2] = { 0.0, 0.0 };
633  Int_t found[2] = {0}, mult = cluster.Multiplicity();
634 
635  for ( Int_t i = 0; i < mult; ++i) {
636  AliMUONPad* pad = cluster.Pad(i);
637  for (Int_t j = 0; j < 2; ++j) {
638  if (found[j] == 0 && TMath::Abs(pad->Size(j)-width[j]) < fgkDistancePrecision) {
639  xy0[j] = pad->Coord(j);
640  found[j] = 1;
641  }
642  }
643  if (found[0] && found[1]) break;
644  }
645 
646  Double_t min[2], max[2];
647  Int_t cath0 = 0, cath1 = 1;
648  if (cluster.Multiplicity(0) == 0) cath0 = 1;
649  else if (cluster.Multiplicity(1) == 0) cath1 = 0;
650 
651  Double_t leftDownX, leftDownY;
652  cluster.Area(cath0).LeftDownCorner(leftDownX, leftDownY);
653  Double_t rightUpX, rightUpY;
654  cluster.Area(cath0).RightUpCorner(rightUpX, rightUpY);
655  min[0] = leftDownX;
656  min[1] = leftDownY;
657  max[0] = rightUpX;
658  max[1] = rightUpY;
659  if (cath1 != cath0) {
660  cluster.Area(cath1).LeftDownCorner(leftDownX, leftDownY);
661  cluster.Area(cath1).RightUpCorner(rightUpX, rightUpY);
662  min[0] = TMath::Max (min[0], leftDownX);
663  min[1] = TMath::Max (min[1], leftDownY);
664  max[0] = TMath::Min (max[0], rightUpX);
665  max[1] = TMath::Min (max[1], rightUpY);
666  }
667 
668  // Adjust limits
669  if (cath0 != cath1) {
670  TVector2 dim0 = cluster.MinPadDimensions (0, -1, kFALSE);
671  TVector2 dim1 = cluster.MinPadDimensions (1, -1, kFALSE);
672  if (TMath::Abs(dim0.Y()-dim1.Y()) < fgkDistancePrecision) {
673  // The same size of pads on both cathodes - check position
674  AliMUONPad* pad0 = cluster.Pad(0);
675  for ( Int_t i = 1; i < mult; ++i) {
676  AliMUONPad* pad = cluster.Pad(i);
677  if (pad->Cathode() == pad0->Cathode()) continue;
678  Double_t dist = TMath::Abs (pad0->Coord(1) - pad->Coord(1));
679  Double_t dd = dist - Int_t(dist/width[1]/2.) * width[1] * 2.;
680  if (TMath::Abs(dd/width[1]/2.-0.5) < fgkDistancePrecision) {
681  // Half pad shift between cathodes
682  width[0] /= 2.;
683  width[1] /= 2.;
684  }
685  break;
686  }
687  }
688  }
689 
690  Int_t nbins[2];
691  for (Int_t i = 0; i < 2; ++i) {
692  Double_t dist = (min[i] - xy0[i]) / width[i] / 2;
693  if (TMath::Abs(dist) < 1.e-6) dist = -1.e-6;
694  min[i] = xy0[i] + (TMath::Nint(dist-TMath::Sign(1.e-6,dist))
695  + TMath::Sign(0.5,dist)) * width[i] * 2;
696  nbins[i] = TMath::Nint ((max[i] - min[i]) / width[i] / 2);
697  if (nbins[i] == 0) ++nbins[i];
698  max[i] = min[i] + nbins[i] * width[i] * 2;
699  //cout << dist << " " << min[i] << " " << max[i] << " " << nbins[i] << endl;
700  }
701 
702  // Book histogram
703  TH2D *hist1 = new TH2D ("Grid", "", nbins[0], min[0], max[0], nbins[1], min[1], max[1]);
704  TH2D *hist2 = new TH2D ("Entries", "", nbins[0], min[0], max[0], nbins[1], min[1], max[1]);
705  TAxis *xaxis = hist1->GetXaxis();
706  TAxis *yaxis = hist1->GetYaxis();
707 
708  // Fill histogram
709  for ( Int_t i = 0; i < mult; ++i) {
710  AliMUONPad* pad = cluster.Pad(i);
711  Int_t ix0 = xaxis->FindBin(pad->X());
712  Int_t iy0 = yaxis->FindBin(pad->Y());
713  PadOverHist(0, ix0, iy0, pad, hist1, hist2);
714  }
715 
716  // Store pixels
717  for (Int_t i = 1; i <= nbins[0]; ++i) {
718  Double_t x = xaxis->GetBinCenter(i);
719  for (Int_t j = 1; j <= nbins[1]; ++j) {
720  if (hist2->GetBinContent(hist2->GetBin(i,j)) < 0.01525) continue; // JC: adc -> fc
721  if (cath0 != cath1) {
722  // Two-sided cluster
723  Double_t cont = hist2->GetBinContent(hist2->GetBin(i,j));
724  if (cont < 999.) continue;
725  if (cont-Int_t(cont/1000.)*1000. < 0.07625) continue; // JC: adc -> fc
726  }
727  //if (hist2->GetBinContent(hist2->GetBin(i,j)) < 1.1 && cluster.Multiplicity(0) &&
728  // cluster.Multiplicity(1)) continue;
729  Double_t y = yaxis->GetBinCenter(j);
730  Double_t charge = hist1->GetBinContent(hist1->GetBin(i,j));
731  AliMUONPad* pixPtr = new AliMUONPad(x, y, width[0], width[1], charge);
732  fPixArray->Add(pixPtr);
733  }
734  }
735  /*
736  if (fPixArray->GetEntriesFast() == 1) {
737  // Split pixel into 2
738  AliMUONPad* pixPtr = static_cast<AliMUONPad*> (fPixArray->UncheckedAt(0));
739  pixPtr->SetSize(0,width[0]/2.);
740  pixPtr->Shift(0,-width[0]/4.);
741  pixPtr = new AliMUONPad(pixPtr->X()+width[0], pixPtr->Y(), width[0]/2., width[1], pixPtr->Charge());
742  fPixArray->Add(pixPtr);
743  }
744  */
745  //fPixArray->Print();
746  delete hist1;
747  delete hist2;
748 }
749 
750 //_____________________________________________________________________________
751 void AliMUONClusterFinderPeakFit::PadOverHist(Int_t idir, Int_t ix0, Int_t iy0, AliMUONPad *pad,
752  TH2D *hist1, TH2D *hist2)
753 {
755 
756  TAxis *axis = idir == 0 ? hist1->GetXaxis() : hist1->GetYaxis();
757  Int_t nbins = axis->GetNbins(), cath = pad->Cathode();
758  Double_t bin = axis->GetBinWidth(1), amask = TMath::Power(1000.,cath*1.);
759 
760  Int_t nbinPad = (Int_t)(pad->Size(idir)/bin*2+fgkDistancePrecision) + 1; // number of bins covered by pad
761 
762  for (Int_t i = 0; i < nbinPad; ++i) {
763  Int_t ixy = idir == 0 ? ix0 + i : iy0 + i;
764  if (ixy > nbins) break;
765  Double_t lowEdge = axis->GetBinLowEdge(ixy);
766  if (lowEdge + fgkDistancePrecision > pad->Coord(idir) + pad->Size(idir)) break;
767  if (idir == 0) PadOverHist(1, ixy, iy0, pad, hist1, hist2); // span in the other direction
768  else {
769  // Fill histogram
770  Double_t cont = pad->Charge();
771  if (hist2->GetBinContent(hist2->GetBin(ix0, ixy)) > 0.01525) // JC: adc -> fc
772  cont = TMath::Min (hist1->GetBinContent(hist1->GetBin(ix0, ixy)), cont)
773  + TMath::Min (TMath::Max(hist1->GetBinContent(hist1->GetBin(ix0, ixy)),cont)*0.1, 1.525); // JC: adc -> fc
774  hist1->SetBinContent(hist1->GetBin(ix0, ixy), cont);
775  hist2->SetBinContent(hist2->GetBin(ix0, ixy), hist2->GetBinContent( hist2->GetBin(ix0, ixy))+amask);
776  }
777  }
778 
779  for (Int_t i = -1; i > -nbinPad; --i) {
780  Int_t ixy = idir == 0 ? ix0 + i : iy0 + i;
781  if (ixy < 1) break;
782  Double_t upEdge = axis->GetBinUpEdge(ixy);
783  if (upEdge - fgkDistancePrecision < pad->Coord(idir) - pad->Size(idir)) break;
784  if (idir == 0) PadOverHist(1, ixy, iy0, pad, hist1, hist2); // span in the other direction
785  else {
786  // Fill histogram
787  Double_t cont = pad->Charge();
788  if (hist2->GetBinContent(hist2->GetBin(ix0, ixy)) > 0.01525) // JC: adc -> fc
789  cont = TMath::Min (hist1->GetBinContent(hist1->GetBin(ix0, ixy)), cont)
790  + TMath::Min (TMath::Max(hist1->GetBinContent(hist1->GetBin(ix0, ixy)),cont)*0.1,1.525); // JC: adc -> fc
791  hist1->SetBinContent(hist1->GetBin(ix0, ixy), cont);
792  hist2->SetBinContent(hist2->GetBin(ix0, ixy), hist2->GetBinContent( hist2->GetBin(ix0, ixy))+amask);
793  }
794  }
795 }
796 
797 //_____________________________________________________________________________
798 Int_t AliMUONClusterFinderPeakFit::FindLocalMaxima(TObjArray *pixArray, Int_t *localMax, Double_t *maxVal)
799 {
801 
802  AliDebug(1,Form("nPix=%d",pixArray->GetLast()+1));
803 
804  //TH2D *hist = NULL;
805  //delete ((TH2D*) gROOT->FindObject("anode"));
806  //if (pixArray == fPixArray) hist = (TH2D*) gROOT->FindObject("anode");
807  //else { hist = (TH2D*) gROOT->FindObject("anode1"); cout << hist << endl; }
808  //if (hist) hist->Delete();
809  delete fHistAnode;
810 
811  Double_t xylim[4] = {999, 999, 999, 999};
812 
813  Int_t nPix = pixArray->GetEntriesFast();
814 
815  if ( nPix <= 0 ) return 0;
816 
817  AliMUONPad *pixPtr = 0;
818  for (Int_t ipix = 0; ipix < nPix; ++ipix) {
819  pixPtr = (AliMUONPad*) pixArray->UncheckedAt(ipix);
820  for (Int_t i = 0; i < 4; ++i)
821  xylim[i] = TMath::Min (xylim[i], (i%2 ? -1 : 1)*pixPtr->Coord(i/2));
822  }
823  for (Int_t i = 0; i < 4; ++i) xylim[i] -= pixPtr->Size(i/2);
824 
825  Int_t nx = TMath::Nint ((-xylim[1]-xylim[0])/pixPtr->Size(0)/2);
826  Int_t ny = TMath::Nint ((-xylim[3]-xylim[2])/pixPtr->Size(1)/2);
827  if (pixArray == fPixArray) fHistAnode = new TH2D("anode","anode",nx,xylim[0],-xylim[1],ny,xylim[2],-xylim[3]);
828  else fHistAnode = new TH2D("anode1","anode1",nx,xylim[0],-xylim[1],ny,xylim[2],-xylim[3]);
829  for (Int_t ipix = 0; ipix < nPix; ++ipix) {
830  pixPtr = (AliMUONPad*) pixArray->UncheckedAt(ipix);
831  fHistAnode->Fill(pixPtr->Coord(0), pixPtr->Coord(1), pixPtr->Charge());
832  }
833 // if (fDraw && pixArray == fPixArray) fDraw->DrawHist("c2", hist);
834 
835  Int_t nMax = 0, indx, nxy = ny * nx;
836  Int_t *isLocalMax = new Int_t[nxy];
837  for (Int_t i = 0; i < nxy; ++i) isLocalMax[i] = 0;
838 
839  for (Int_t i = 1; i <= ny; ++i) {
840  indx = (i-1) * nx;
841  for (Int_t j = 1; j <= nx; ++j) {
842  if (fHistAnode->GetBinContent(fHistAnode->GetBin(j,i)) < 0.07625) continue; // JC: adc -> fc
843  //if (isLocalMax[indx+j-1] < 0) continue;
844  if (isLocalMax[indx+j-1] != 0) continue;
845  FlagLocalMax(fHistAnode, i, j, isLocalMax);
846  }
847  }
848 
849  for (Int_t i = 1; i <= ny; ++i) {
850  indx = (i-1) * nx;
851  for (Int_t j = 1; j <= nx; ++j) {
852  if (isLocalMax[indx+j-1] > 0) {
853  localMax[nMax] = indx + j - 1;
854  maxVal[nMax++] = fHistAnode->GetBinContent(fHistAnode->GetBin(j,i));
855  if (nMax > 99) break;
856  }
857  }
858  if (nMax > 99) {
859  AliError(" Too many local maxima !!!");
860  break;
861  }
862  }
863  if (fDebug) cout << " Local max: " << nMax << endl;
864  delete [] isLocalMax;
865  return nMax;
866 }
867 
868 //_____________________________________________________________________________
869 void AliMUONClusterFinderPeakFit::FlagLocalMax(TH2D *hist, Int_t i, Int_t j, Int_t *isLocalMax)
870 {
872 
873  Int_t nx = hist->GetNbinsX();
874  Int_t ny = hist->GetNbinsY();
875  Int_t cont = TMath::Nint (hist->GetBinContent(hist->GetBin(j,i)));
876  Int_t cont1 = 0, indx = (i-1)*nx+j-1, indx1 = 0, indx2 = 0;
877 
878  Int_t ie = i + 2, je = j + 2;
879  for (Int_t i1 = i-1; i1 < ie; ++i1) {
880  if (i1 < 1 || i1 > ny) continue;
881  indx1 = (i1 - 1) * nx;
882  for (Int_t j1 = j-1; j1 < je; ++j1) {
883  if (j1 < 1 || j1 > nx) continue;
884  if (i == i1 && j == j1) continue;
885  indx2 = indx1 + j1 - 1;
886  cont1 = TMath::Nint (hist->GetBinContent(hist->GetBin(j1,i1)));
887  if (cont < cont1) { isLocalMax[indx] = -1; return; }
888  else if (cont > cont1) isLocalMax[indx2] = -1;
889  else { // the same charge
890  isLocalMax[indx] = 1;
891  if (isLocalMax[indx2] == 0) {
892  FlagLocalMax(hist, i1, j1, isLocalMax);
893  if (isLocalMax[indx2] < 0) { isLocalMax[indx] = -1; return; }
894  else isLocalMax[indx2] = -1;
895  }
896  }
897  }
898  }
899  isLocalMax[indx] = 1; // local maximum
900 }
901 
902 //_____________________________________________________________________________
903 void AliMUONClusterFinderPeakFit::FindClusterFit(AliMUONCluster& cluster, const Int_t *localMax,
904  const Int_t *maxPos, Int_t nMax)
905 {
907 
908  //if (cluster.Multiplicity(0) == 0 || cluster.Multiplicity(1) == 0) cout << cluster.Multiplicity(0) << " " << cluster.Multiplicity(1) << " " << cluster.Charge() << " " << cluster.Charge(0) << " " << cluster.Charge(1) << " " << endl;
909 
910  //TH2D *hist = (TH2D*) gROOT->FindObject("anode");
911  Int_t nx = fHistAnode->GetNbinsX();
912  //Int_t ny = hist->GetNbinsY();
913  Double_t xmin = fHistAnode->GetXaxis()->GetXmin(); //- hist->GetXaxis()->GetBinWidth(1);
914  Double_t xmax = fHistAnode->GetXaxis()->GetXmax(); //+ hist->GetXaxis()->GetBinWidth(1);
915  Double_t ymin = fHistAnode->GetYaxis()->GetXmin(); //- hist->GetYaxis()->GetBinWidth(1);
916  Double_t ymax = fHistAnode->GetYaxis()->GetXmax(); //+ hist->GetYaxis()->GetBinWidth(1);
917 
918  TVirtualFitter* fitter = TVirtualFitter::Fitter(0,nMax*3);
919  fitter->Clear("");
920  fitter->SetFCN(FitFunction);
921 
922  Float_t stepX = 0.01; // cm
923  Float_t stepY = 0.01; // cm
924  Float_t stepQ = 0.01; //
925 
926  Double_t args[10] = {-1.}; // disable printout
927 
928  fitter->ExecuteCommand("SET PRINT",args,1);
929  fitter->ExecuteCommand("SET NOW",args,0); // no warnings
930 
931  Int_t indx = 0;
932  fNMax = nMax;
933  for (Int_t i = 0; i < nMax; ++i) {
934  Int_t ic = localMax[maxPos[i]] / nx + 1;
935  Int_t jc = localMax[maxPos[i]] % nx + 1;
936  Double_t yc = fHistAnode->GetYaxis()->GetBinCenter(ic);
937  Double_t xc = fHistAnode->GetXaxis()->GetBinCenter(jc);
938  indx = 3 * i;
939  fitter->SetParameter(indx,"Hit X position",xc,stepX,xmin,xmax);
940  fitter->SetParameter(indx+1,"Hit Y position",yc,stepY,ymin,ymax);
941  fitter->SetParameter(indx+2,"Hit contribution",0.6,stepQ,0.,1.);
942  }
943  fitter->SetParameter(indx+2,"Hit contribution",0.,0.,0,0);
944  //fitter->SetParameter(8,"Number of hits",nMax,0.,0,0);
945 
946  TObjArray userObjects;
947 
948  userObjects.Add(&cluster);
949  userObjects.Add(fMathieson);
950  userObjects.Add(this);
951 
952  fitter->SetObjectFit(&userObjects);
953 
954  args[0] = 500.;
955  args[1] = 1.;
956  /*Int_t stat =*/ fitter->ExecuteCommand("MIGRAD",args,2);
957  //if (stat) { cout << " stat = " << stat << " " << fDetElemId << endl; /*exit(0);*/ }
958  //Int_t nvpar, nparx;
959  //Double_t amin, edm, errdef;
960  //fitter->GetStats(amin, edm, errdef, nvpar, nparx);
961  //cout << amin << endl;
962 
963  Double_t qTot = cluster.Charge(), par[9] = {0.}, err[9] = {0.}, coef = 0.;
964  //par[8] = nMax;
965  for (Int_t j = 0; j < nMax; ++j) {
966  indx = 3 * j;
967  par[indx+2] = fitter->GetParameter(indx+2);
968  coef = Param2Coef(j, coef, par, nMax);
969  par[indx] = fitter->GetParameter(indx);
970  par[indx+1] = fitter->GetParameter(indx+1);
971  err[indx] = fitter->GetParError(indx);
972  err[indx+1] = fitter->GetParError(indx+1);
973 
974  if ( coef*qTot >= 2.135 ) // JC: adc -> fc
975  {
976  AliMUONCluster* cluster1 = new AliMUONCluster(cluster);
977 
978  cluster1->SetCharge(coef*qTot,coef*qTot);
979 
980  cluster1->SetPosition(TVector2(par[indx],par[indx+1]),TVector2(err[indx],err[indx+1]));
981  cluster1->SetChi2(0.);
982 
983  // FIXME: we miss some information in this cluster, as compared to
984  // the original AddRawCluster code.
985 
986  AliDebug(2,Form("Adding RawCluster detElemId %4d mult %2d charge %5d (xl,yl)=(%9.6g,%9.6g)",
987  fDetElemId,cluster1->Multiplicity(),(Int_t)cluster1->Charge(),
988  cluster1->Position().X(),cluster1->Position().Y()));
989 
990  fClusterList.Add(cluster1);
991  }
992  }
993 }
994 
995 //_____________________________________________________________________________
997  const Int_t *localMax, Int_t iMax)
998 {
1000 
1001  //TH2D *hist = (TH2D*) gROOT->FindObject("anode");
1002  /* Just for check
1003  TCanvas* c = new TCanvas("Anode","Anode",800,600);
1004  c->cd();
1005  hist->Draw("lego1Fb"); // debug
1006  c->Update();
1007  Int_t tmp;
1008  cin >> tmp;
1009  */
1010  Int_t nx = fHistAnode->GetNbinsX();
1011  //Int_t ny = hist->GetNbinsY();
1012  Int_t ic = localMax[iMax] / nx + 1;
1013  Int_t jc = localMax[iMax] % nx + 1;
1014 
1015  // Get min pad dimensions for the precluster
1016  Int_t nSides = 2;
1017  if (cluster.Multiplicity(0) == 0 || cluster.Multiplicity(1) == 0) nSides = 1;
1018  TVector2 dim0 = cluster.MinPadDimensions(0, -1, kFALSE);
1019  TVector2 dim1 = cluster.MinPadDimensions(1, -1, kFALSE);
1020  //Double_t width[2][2] = {{dim0.X(), dim0.Y()},{dim1.X(),dim1.Y()}};
1021  Int_t nonb[2] = {1, 0}; // coordinate index vs cathode
1022  if (nSides == 1 || dim0.X() < dim1.X() - fgkDistancePrecision) {
1023  nonb[0] = 0;
1024  nonb[1] = 1;
1025  }
1026 
1027  // Drop all pixels from the array - pick up only the ones from the cluster
1028  //fPixArray->Delete();
1029 
1030  Double_t wx = fHistAnode->GetXaxis()->GetBinWidth(1)/2;
1031  Double_t wy = fHistAnode->GetYaxis()->GetBinWidth(1)/2;
1032  Double_t yc = fHistAnode->GetYaxis()->GetBinCenter(ic);
1033  Double_t xc = fHistAnode->GetXaxis()->GetBinCenter(jc);
1034  Double_t cont = fHistAnode->GetBinContent(fHistAnode->GetBin(jc,ic));
1035  AliMUONPad pixel(xc, yc, wx, wy, cont);
1036  if (fDebug) pixel.Print("full");
1037 
1038  Int_t npad = cluster.Multiplicity();
1039 
1040  // Pick up pads which overlap with the maximum pixel and find pads with the max signal
1041  Double_t qMax[2] = {0};
1042  AliMUONPad *matrix[2][3] = {{0x0,0x0,0x0},{0x0,0x0,0x0}};
1043  for (Int_t j = 0; j < npad; ++j)
1044  {
1045  AliMUONPad* pad = cluster.Pad(j);
1046  if ( Overlap(*pad,pixel) )
1047  {
1048  if (fDebug) { cout << j << " "; pad->Print("full"); }
1049  if (pad->Charge() > qMax[pad->Cathode()]) {
1050  qMax[pad->Cathode()] = pad->Charge();
1051  matrix[pad->Cathode()][1] = pad;
1052  if (nSides == 1) matrix[!pad->Cathode()][1] = pad;
1053  }
1054  }
1055  }
1056  //if (nSides == 2 && (matrix[0][1] == 0x0 || matrix[1][1] == 0x0)) return; // ???
1057 
1058  // Find neighbours of maxima to have 3 pads per direction (if possible)
1059  for (Int_t j = 0; j < npad; ++j)
1060  {
1061  AliMUONPad* pad = cluster.Pad(j);
1062  Int_t cath = pad->Cathode();
1063  if (pad == matrix[cath][1]) continue;
1064  Int_t nLoops = 3 - nSides;
1065 
1066  for (Int_t k = 0; k < nLoops; ++k) {
1067  Int_t cath1 = cath;
1068  if (k) cath1 = !cath;
1069 
1070  // Check the coordinate corresponding to the cathode (bending or non-bending case)
1071  Double_t dist = pad->Coord(nonb[cath1]) - matrix[cath][1]->Coord(nonb[cath1]);
1072  Double_t dir = TMath::Sign (1., dist);
1073  dist = TMath::Abs(dist) - pad->Size(nonb[cath1]) - matrix[cath][1]->Size(nonb[cath1]);
1074 
1075  if (TMath::Abs(dist) < fgkDistancePrecision) {
1076  // Check the other coordinate
1077  dist = pad->Coord(!nonb[cath1]) - matrix[cath1][1]->Coord(!nonb[cath1]);
1078  if (TMath::Abs(dist) >
1079  TMath::Max(pad->Size(!nonb[cath1]), matrix[cath1][1]->Size(!nonb[cath1])) - fgkDistancePrecision) break;
1080  Int_t idir = TMath::Nint (dir);
1081  if (matrix[cath1][1+idir] == 0x0) matrix[cath1][1+idir] = pad;
1082  else if (pad->Charge() > matrix[cath1][1+idir]->Charge()) matrix[cath1][1+idir] = pad; // diff. segmentation
1083  //cout << pad->Coord(nonb[cath1]) << " " << pad->Coord(!nonb[cath1]) << " " << pad->Size(nonb[cath1]) << " " << pad->Size(!nonb[cath1]) << " " << pad->Charge() << endl ;
1084  break;
1085  }
1086  }
1087  }
1088 
1089  Double_t coord[2] = {0.}, qAver = 0.;
1090  for (Int_t i = 0; i < 2; ++i) {
1091  Double_t q = 0.;
1092  Double_t coordQ = 0.;
1093  Int_t cath = matrix[i][1]->Cathode();
1094  if (i && nSides == 1) cath = !cath;
1095  for (Int_t j = 0; j < 3; ++j) {
1096  if (matrix[i][j] == 0x0) continue;
1097  Double_t dq = matrix[i][j]->Charge();
1098  q += dq;
1099  coordQ += dq * matrix[i][j]->Coord(nonb[cath]);
1100  //coordQ += (matrix[i][j]->Charge() * matrix[i][j]->Coord(nonb[cath]));
1101  }
1102  coord[cath] = coordQ / q;
1103  qAver = TMath::Max (qAver, q);
1104  }
1105 
1106  //qAver = TMath::Sqrt(qAver);
1107  if ( qAver >= 2.135 ) // JC: adc -> fc
1108  {
1109 
1110  AliMUONCluster* cluster1 = new AliMUONCluster(cluster);
1111 
1112  cluster1->SetCharge(qAver,qAver);
1113  if (nonb[0] == 1)
1114  cluster1->SetPosition(TVector2(coord[1],coord[0]),TVector2(0.,0.));
1115  else
1116  cluster1->SetPosition(TVector2(coord[0],coord[1]),TVector2(0.,0.));
1117 
1118  cluster1->SetChi2(0.);
1119 
1120  // FIXME: we miss some information in this cluster, as compared to
1121  // the original AddRawCluster code.
1122 
1123  AliDebug(2,Form("Adding RawCluster detElemId %4d mult %2d charge %5d (xl,yl)=(%9.6g,%9.6g)",
1124  fDetElemId,cluster1->Multiplicity(),(Int_t)cluster1->Charge(),
1125  cluster1->Position().X(),cluster1->Position().Y()));
1126 
1127  fClusterList.Add(cluster1);
1128  }
1129 }
1130 
1131 //_____________________________________________________________________________
1134 {
1136 
1137  if (this == &rhs) return *this;
1138 
1139  AliFatal("Not implemented.");
1140 
1141  return *this;
1142 }
1143 
1144 //_____________________________________________________________________________
1146  Int_t &nInX, Int_t &nInY) const
1147 {
1150 
1151  //Int_t statusToTest = 1;
1152  Int_t statusToTest = fgkUseForFit;
1153 
1154  //if ( nInX < 0 ) statusToTest = 0;
1155  if ( nInX < 0 ) statusToTest = fgkZero;
1156 
1157  Bool_t mustMatch(kTRUE);
1158 
1159  Long_t cn = cluster.NofPads(statusToTest,mustMatch);
1160 
1161  nInX = AliMp::PairFirst(cn);
1162  nInY = AliMp::PairSecond(cn);
1163 }
1164 
1165 //_____________________________________________________________________________
1167 {
1169  AliMUONPad* pixPtr = Pixel(i);
1170  fPixArray->RemoveAt(i);
1171  delete pixPtr;
1172 }
1173 
1174 //_____________________________________________________________________________
1175 AliMUONPad*
1177 {
1179  return static_cast<AliMUONPad*>(fPixArray->UncheckedAt(i));
1180 }
1181 
1182 //_____________________________________________________________________________
1183 void
1185 {
1187  TString swhat(what);
1188  swhat.ToLower();
1189  if ( swhat.Contains("precluster") )
1190  {
1191  if ( fPreCluster) fPreCluster->Print();
1192  }
1193 }
1194 
1195 
static Bool_t AreOverlapping(const AliMUONPad &d1, const AliMUONPad &d2, const TVector2 &precision, AliMpArea &overlapArea)
Definition: AliMUONPad.cxx:221
static AliMq::Station12Type GetStation12Type(Int_t detElemId)
AliMUONCluster * CheckPrecluster(const AliMUONCluster &cluster)
Check precluster to simplify it (if possible), and return the simplified cluster. ...
void PadsInXandY(AliMUONCluster &cluster, Int_t &nInX, Int_t &nInY) const
printf("Chi2/npoints = %f\n", TMath::Sqrt(chi2/npoints))
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.
static const Int_t fgkUseForFit
should be used for fit
void SetSqrtKx3AndDeriveKx2Kx4(Float_t SqrtKx3)
Mathieson sqrt{Kx3} and derived Kx2 and Kx4.
static Float_t Pitch()
Return wire pitch.
AliMUONClusterFinderPeakFit & operator=(const AliMUONClusterFinderPeakFit &rhs)
Not implemented.
Bool_t Overlap(const AliMUONPad &pad, const AliMUONPad &pixel)
Checks whether a pad and a pixel have an overlapping area.
Bool_t IsSaturated() const
Return info whether this pad is saturated or not.
Definition: AliMUONPad.h:67
void BuildPixArrayOneCathode(AliMUONCluster &cluster)
Float_t IntXY(Float_t xi1, Float_t yi1, Float_t xi2, Float_t yi2) const
Charge integration on region (x1,y1,x2,y2).
#define TObjArray
A group of adjacent pads.
static Float_t SqrtKy3St1()
Return SqrtKy3 for Station 1 & 2.
Int_t fClusterNumber
! current cluster number
AliMUONCluster * CheckPreclusterTwoCathodes(AliMUONCluster *cluster)
static Float_t SqrtKx3()
Return SqrtKx3 for Slat.
TFile f("CalibObjects.root")
Implementation of Mathieson response.
TVector2 Dimensions() const
Return half dimensions in x and y (cm)
Definition: AliMUONPad.h:56
virtual void Print(Option_t *opt="") const
Int_t fNAddVirtualPads
! number of clusters for which we added virtual pads
A rectangle area positioned in plane..
Definition: AliMpArea.h:20
TH2D * fHistAnode
! histo for local maxima search
virtual void Print(Option_t *opt="") const
static const Double_t fgkDistancePrecision
used to check overlaps and so on
Int_t fEventNumber
! current event being processed
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
AliMUONCluster * fPreCluster
! current pre-cluster
void SetPosition(const TVector2 &pos, const TVector2 &errorOnPos)
Set (x,y) of that cluster and errors.
static Float_t PitchSt1()
Return wire pitch for Station 1 & 2.
static Float_t SqrtKy3()
Return SqrtKy3 for Slat.
virtual Bool_t UsePad(const AliMUONPad &pad)
ClassImp(TPCGenInfo)
Definition: AliTPCCmpNG.C:254
void PadOverHist(Int_t idir, Int_t ix0, Int_t iy0, AliMUONPad *pad, TH2D *h1, TH2D *h2)
Double_t Y() const
Return position in y (cm)
Definition: AliMUONPad.h:128
static Float_t SqrtKx3St1()
Return SqrtKx3 for Station 1 & 2.
return clusterFinder
virtual AliMUONCluster * NextCluster()
Int_t Iy() const
Return y-index.
Definition: AliMUONPad.h:80
static const TVector2 fgkDecreaseSize
idem
void Release()
Detach this pad from a cluster.
Definition: AliMUONPad.h:90
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
void SetPitch(Float_t p1)
Int_t Cathode() const
Return cathode number.
Definition: AliMUONPad.h:43
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.
static const Int_t fgkZero
pad "basic" state
Double_t Size(Int_t ixy) const
Definition: AliMUONPad.cxx:508
void RemovePad(AliMUONPad *pad)
void FindClusterFit(AliMUONCluster &cluster, const Int_t *localMax, const Int_t *maxPos, Int_t nMax)
TVector2 Position() const
Return (x,y) of that cluster.
void SetChi2(Float_t chi2)
Set chi2 of the RawCharge fit.
Bool_t IsValid() const
Return validity.
Definition: AliMpPad.h:89
Int_t fDetElemId
! current DE being processed
TVector2 Position() const
Return positions in x and y (cm)
Definition: AliMUONPad.h:85
virtual Bool_t Prepare(Int_t detElemId, TObjArray *pads[2], const AliMpArea &area, const AliMpVSegmentation *seg[2])
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 fNClusters
! total number of clusters
TObjArray * fPixArray
! collection of pixels
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.
void FlagLocalMax(TH2D *hist, Int_t i, Int_t j, Int_t *isLocalMax)
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.
AliMUONClusterFinderPeakFit(Bool_t plot, AliMUONVClusterFinder *clusterFinder)
const AliMpVSegmentation * fkSegmentation[2]
! new segmentation
void Print(Option_t *opt="") const
Definition: AliMUONPad.cxx:420
Double_t Charge() const
Return pad charge.
Definition: AliMUONPad.h:48
TObjArray fClusterList
! clusters corresponding to the current pre-cluster
Double_t Coord(Int_t ixy) const
Definition: AliMUONPad.cxx:354
Class which encapsuate all information about a pad.
Definition: AliMpPad.h:22
AliMUONMathieson * fMathieson
! Mathieson to compute the charge repartition
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.
Int_t FindLocalMaxima(TObjArray *pixArray, Int_t *localMax, Double_t *maxVal)
Cluster finder in MUON arm of ALICE.
Int_t PairSecond(MpPair_t pair)
Decode the second integer from encoded pair.
Bool_t fPlot
! whether we should plot thing (for debug only, quite slow!)
void FindClusterCOG(AliMUONCluster &cluster, const Int_t *localMax, Int_t iMax)
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
Int_t MaxRawChargeCathode() const
Return the max raw charge on the chathod.
AliMUONPad * Pad(Int_t index) const
void BuildPixArray(AliMUONCluster &cluster)
build array of pixels
void SetSqrtKy3AndDeriveKy2Ky4(Float_t SqrtKy3)
Mathieson sqrt{Ky3} and derived Ky2 and Ky4.
AliMUONVClusterFinder * fPreClusterFinder
! the pre-clustering worker
Combination of digit and mppad informations.
Definition: AliMUONPad.h:25
Int_t GetNMax() const
Return the number of local maxima.
Double_t DX() const
Return half dimensions in x (cm)
Definition: AliMUONPad.h:59
Int_t fNMax
! number of local maxima