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