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