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