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