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 AliMUONClusterSplitterMLEM
20 : ///
21 : /// Splitter class for the MLEM algorithm. Performs fitting procedure
22 : /// with up to 3 hit candidates and tries to split clusters if the number
23 : /// of candidates exceeds 3.
24 : ///
25 : /// \author Laurent Aphecetche (for the "new" C++ structure) and
26 : /// Alexander Zinchenko, JINR Dubna, for the hardcore of it ;-)
27 : //-----------------------------------------------------------------------------
28 :
29 : #include "AliMUONClusterSplitterMLEM.h"
30 : #include "AliMUONClusterFinderMLEM.h" // for status flag constants
31 :
32 : #include "AliMUONCluster.h"
33 : #include "AliMUONPad.h"
34 : #include "AliMUONPad.h"
35 : #include "AliMUONConstants.h"
36 : #include "AliMpDEManager.h"
37 : #include "AliMUONMathieson.h"
38 :
39 : #include "AliMpEncodePair.h"
40 :
41 : #include "AliLog.h"
42 :
43 : #include <TClonesArray.h>
44 : #include <TH2.h>
45 : #include <TMath.h>
46 : #include <TMatrixD.h>
47 : #include <TObjArray.h>
48 : #include <TRandom.h>
49 : #include <Riostream.h>
50 :
51 : using std::endl;
52 : using std::cout;
53 : /// \cond CLASSIMP
54 18 : ClassImp(AliMUONClusterSplitterMLEM)
55 : /// \endcond
56 :
57 : //const Double_t AliMUONClusterSplitterMLEM::fgkCouplMin = 1.e-3; // threshold on coupling
58 : const Double_t AliMUONClusterSplitterMLEM::fgkCouplMin = 1.e-2; // threshold on coupling
59 :
60 : //_____________________________________________________________________________
61 : AliMUONClusterSplitterMLEM::AliMUONClusterSplitterMLEM(Int_t detElemId,
62 : TObjArray* pixArray,
63 : Double_t lowestPixelCharge,
64 : Double_t lowestPadCharge,
65 : Double_t lowestClusterCharge)
66 144 : : TObject(),
67 144 : fPixArray(pixArray),
68 144 : fMathieson(0x0),
69 144 : fDetElemId(detElemId),
70 144 : fNpar(0),
71 144 : fQtot(0),
72 144 : fnCoupled(0),
73 144 : fDebug(0),
74 144 : fLowestPixelCharge(lowestPixelCharge),
75 144 : fLowestPadCharge(lowestPadCharge),
76 144 : fLowestClusterCharge(lowestClusterCharge)
77 720 : {
78 : /// Constructor
79 :
80 144 : AliMq::Station12Type stationType = AliMpDEManager::GetStation12Type(fDetElemId);
81 :
82 144 : Float_t kx3 = AliMUONConstants::SqrtKx3();
83 144 : Float_t ky3 = AliMUONConstants::SqrtKy3();
84 144 : Float_t pitch = AliMUONConstants::Pitch();
85 :
86 144 : if ( stationType == AliMq::kStation1 )
87 : {
88 28 : kx3 = AliMUONConstants::SqrtKx3St1();
89 28 : ky3 = AliMUONConstants::SqrtKy3St1();
90 28 : pitch = AliMUONConstants::PitchSt1();
91 28 : }
92 :
93 432 : fMathieson = new AliMUONMathieson;
94 :
95 144 : fMathieson->SetPitch(pitch);
96 144 : fMathieson->SetSqrtKx3AndDeriveKx2Kx4(kx3);
97 144 : fMathieson->SetSqrtKy3AndDeriveKy2Ky4(ky3);
98 :
99 288 : }
100 :
101 : //_____________________________________________________________________________
102 : AliMUONClusterSplitterMLEM::~AliMUONClusterSplitterMLEM()
103 864 : {
104 : /// Destructor
105 :
106 288 : delete fMathieson;
107 432 : }
108 :
109 : //_____________________________________________________________________________
110 : void
111 : AliMUONClusterSplitterMLEM::AddBin(TH2 *mlem,
112 : Int_t ic, Int_t jc, Int_t mode,
113 : Bool_t *used, TObjArray *pix)
114 : {
115 : /// Add a bin to the cluster
116 :
117 814 : Int_t nx = mlem->GetNbinsX();
118 407 : Int_t ny = mlem->GetNbinsY();
119 407 : Double_t cont1, cont = mlem->GetBinContent(mlem->GetBin(jc,ic));
120 : AliMUONPad *pixPtr = 0;
121 :
122 407 : Int_t ie = TMath::Min(ic+1,ny), je = TMath::Min(jc+1,nx);
123 3018 : for (Int_t i = TMath::Max(ic-1,1); i <= ie; ++i) {
124 7980 : for (Int_t j = TMath::Max(jc-1,1); j <= je; ++j) {
125 4709 : if (i != ic && j != jc) continue;
126 1762 : if (used[(i-1)*nx+j-1]) continue;
127 706 : cont1 = mlem->GetBinContent(mlem->GetBin(j,i));
128 706 : if (mode && cont1 > cont) continue;
129 706 : used[(i-1)*nx+j-1] = kTRUE;
130 706 : if (cont1 < fLowestPixelCharge) continue;
131 688 : if (pix) pix->Add(BinToPix(mlem,j,i));
132 : else {
133 0 : pixPtr = new AliMUONPad (mlem->GetXaxis()->GetBinCenter(j),
134 0 : mlem->GetYaxis()->GetBinCenter(i), 0, 0, cont1);
135 0 : fPixArray->Add(pixPtr);
136 : }
137 344 : AddBin(mlem, i, j, mode, used, pix); // recursive call
138 344 : }
139 : }
140 407 : }
141 :
142 : //_____________________________________________________________________________
143 : void
144 : AliMUONClusterSplitterMLEM::AddCluster(Int_t ic, Int_t nclust,
145 : TMatrixD& aijcluclu,
146 : Bool_t *used, Int_t *clustNumb, Int_t &nCoupled)
147 : {
148 : /// Add a cluster to the group of coupled clusters
149 :
150 431 : for (Int_t i = 0; i < nclust; ++i) {
151 121 : if (used[i]) continue;
152 21 : if (aijcluclu(i,ic) < fgkCouplMin) continue;
153 21 : used[i] = kTRUE;
154 21 : clustNumb[nCoupled++] = i;
155 21 : AddCluster(i, nclust, aijcluclu, used, clustNumb, nCoupled);
156 21 : }
157 63 : }
158 :
159 : //_____________________________________________________________________________
160 : TObject*
161 : AliMUONClusterSplitterMLEM::BinToPix(TH2 *mlem,
162 : Int_t jc, Int_t ic)
163 : {
164 : /// Translate histogram bin to pixel
165 :
166 900 : Double_t yc = mlem->GetYaxis()->GetBinCenter(ic);
167 450 : Double_t xc = mlem->GetXaxis()->GetBinCenter(jc);
168 :
169 450 : Int_t nPix = fPixArray->GetEntriesFast();
170 : AliMUONPad *pixPtr = NULL;
171 :
172 : // Compare pixel and bin positions
173 5414 : for (Int_t i = 0; i < nPix; ++i) {
174 2707 : pixPtr = (AliMUONPad*) fPixArray->UncheckedAt(i);
175 2707 : if (pixPtr->Charge() < fLowestPixelCharge) continue;
176 3401 : if (TMath::Abs(pixPtr->Coord(0)-xc)<1.e-4 && TMath::Abs(pixPtr->Coord(1)-yc)<1.e-4)
177 : {
178 : //return (TObject*) pixPtr;
179 450 : return pixPtr;
180 : }
181 : }
182 0 : AliError(Form(" Something wrong ??? %f %f ", xc, yc));
183 0 : return NULL;
184 450 : }
185 :
186 : //_____________________________________________________________________________
187 : Float_t
188 : AliMUONClusterSplitterMLEM::ChargeIntegration(Double_t x, Double_t y,
189 : const AliMUONPad& pad)
190 : {
191 : /// Compute the Mathieson integral on pad area, assuming the center
192 : /// of the Mathieson is at (x,y)
193 :
194 280992 : TVector2 lowerLeft(TVector2(x,y)-pad.Position()-pad.Dimensions());
195 187328 : TVector2 upperRight(lowerLeft + pad.Dimensions()*2.0);
196 :
197 93664 : return fMathieson->IntXY(lowerLeft.X(),lowerLeft.Y(),
198 46832 : upperRight.X(),upperRight.Y());
199 46832 : }
200 :
201 : //_____________________________________________________________________________
202 : void
203 : AliMUONClusterSplitterMLEM::Fcn1(const AliMUONCluster& cluster,
204 : Int_t & /*fNpar*/, Double_t * /*gin*/,
205 : Double_t &f, Double_t *par, Int_t iflag)
206 : {
207 : /// Computes the functional to be minimized
208 :
209 : Int_t indx, npads=0;
210 : Double_t charge, delta, coef=0, chi2=0, qTot = 0;
211 : static Double_t qAver = 0;
212 :
213 4710 : Int_t mult = cluster.Multiplicity(), iend = fNpar / 3;
214 44316 : for (Int_t j = 0; j < mult; ++j)
215 : {
216 19803 : AliMUONPad* pad = cluster.Pad(j);
217 : //if ( pad->Status() !=1 || pad->IsSaturated() ) continue;
218 39450 : if ( pad->Status() != AliMUONClusterFinderMLEM::GetUseForFitFlag() ||
219 19803 : pad->Charge() == 0 ) continue;
220 19647 : if (iflag == 0) {
221 2738 : if ( pad->IsReal() ) npads++; // exclude virtual pads
222 1422 : qTot += pad->Charge();
223 1422 : }
224 : charge = 0;
225 78588 : for (Int_t i = 0; i <= iend; ++i)
226 : {
227 : // sum over hits
228 19647 : indx = 3 * i;
229 19647 : coef = Param2Coef(i, coef, par);
230 19647 : charge += ChargeIntegration(par[indx],par[indx+1],*pad) * coef;
231 : }
232 19647 : charge *= fQtot;
233 19647 : delta = charge - pad->Charge();
234 19647 : delta *= delta;
235 19647 : delta /= pad->Charge();
236 19647 : chi2 += delta;
237 19647 : } // for (Int_t j=0;
238 2519 : if (iflag == 0 && npads) qAver = qTot / npads;
239 2355 : if (!npads && iflag==0)
240 : {
241 0 : AliError(Form("Got npads=0. Please check"));
242 0 : }
243 2355 : f = chi2 / qAver;
244 2355 : }
245 :
246 : //_____________________________________________________________________________
247 : Double_t AliMUONClusterSplitterMLEM::Param2Coef(Int_t icand, Double_t coef, Double_t *par) const
248 : {
249 : /// Extract hit contribution scale factor from fit parameters
250 :
251 59433 : if (fNpar == 2) return 1.;
252 0 : if (fNpar == 5) return icand==0 ? par[2] : TMath::Max(1.-par[2],0.);
253 0 : if (icand == 0) return par[2];
254 0 : if (icand == 1) return TMath::Max((1.-par[2])*par[5], 0.);
255 0 : return TMath::Max(1.-par[2]-coef,0.);
256 19811 : }
257 :
258 : //_____________________________________________________________________________
259 : Int_t
260 : AliMUONClusterSplitterMLEM::Fit(const AliMUONCluster& cluster,
261 : Int_t iSimple, Int_t nfit,
262 : const Int_t *clustFit, TObjArray **clusters,
263 : Double_t *parOk,
264 : TObjArray& clusterList, TH2 *mlem)
265 : {
266 : /// Steering function and fitting procedure for the fit of pad charge distribution
267 :
268 : // AliDebug(2,Form("iSimple=%d nfit=%d",iSimple,nfit));
269 :
270 328 : Double_t xmin = mlem->GetXaxis()->GetXmin() - mlem->GetXaxis()->GetBinWidth(1);
271 164 : Double_t xmax = mlem->GetXaxis()->GetXmax() + mlem->GetXaxis()->GetBinWidth(1);
272 164 : Double_t ymin = mlem->GetYaxis()->GetXmin() - mlem->GetYaxis()->GetBinWidth(1);
273 164 : Double_t ymax = mlem->GetYaxis()->GetXmax() + mlem->GetYaxis()->GetBinWidth(1);
274 :
275 : // Number of pads to use and number of virtual pads
276 : Int_t npads = 0, nVirtual = 0, nfit0 = nfit;
277 : //cluster.Print("full");
278 164 : Int_t mult = cluster.Multiplicity();
279 3188 : for (Int_t i = 0; i < mult; ++i )
280 : {
281 1430 : AliMUONPad* pad = cluster.Pad(i);
282 1536 : if ( !pad->IsReal() ) ++nVirtual;
283 : //if ( pad->Status() !=1 || pad->IsSaturated() ) continue;
284 1438 : if ( pad->Status() != AliMUONClusterFinderMLEM::GetUseForFitFlag() ) continue;
285 1422 : if ( pad->IsReal() )
286 : {
287 1316 : ++npads;
288 1316 : }
289 1422 : }
290 :
291 164 : fNpar = 0;
292 164 : fQtot = 0;
293 :
294 164 : if (npads < 2) return 0;
295 :
296 : // FIXME : AliWarning("Reconnect the following code for hit/track passing ?");
297 :
298 : // Int_t tracks[3] = {-1, -1, -1};
299 :
300 : /*
301 : Int_t digit = 0;
302 : AliMUONDigit *mdig = 0;
303 : for (Int_t cath=0; cath<2; cath++) {
304 : for (Int_t i=0; i<fnPads[0]+fnPads[1]; i++) {
305 : if (fPadIJ[0][i] != cath) continue;
306 : if (fPadIJ[1][i] != 1) continue;
307 : if (fXyq[3][i] < 0) continue; // exclude virtual pads
308 : digit = TMath::Nint (fXyq[5][i]);
309 : if (digit >= 0) mdig = fInput->Digit(cath,digit);
310 : else mdig = fInput->Digit(TMath::Even(cath),-digit-1);
311 : //if (!mdig) mdig = fInput->Digit(TMath::Even(cath),digit);
312 : if (!mdig) continue; // protection for cluster display
313 : if (mdig->Hit() >= 0) {
314 : if (tracks[0] < 0) {
315 : tracks[0] = mdig->Hit();
316 : tracks[1] = mdig->Track(0);
317 : } else if (mdig->Track(0) < tracks[1]) {
318 : tracks[0] = mdig->Hit();
319 : tracks[1] = mdig->Track(0);
320 : }
321 : }
322 : if (mdig->Track(1) >= 0 && mdig->Track(1) != tracks[1]) {
323 : if (tracks[2] < 0) tracks[2] = mdig->Track(1);
324 : else tracks[2] = TMath::Min (tracks[2], mdig->Track(1));
325 : }
326 : } // for (Int_t i=0;
327 : } // for (Int_t cath=0;
328 : */
329 :
330 : // Get number of pads in X and Y
331 : //const Int_t kStatusToTest(1);
332 164 : const Int_t kStatusToTest(AliMUONClusterFinderMLEM::GetUseForFitFlag());
333 :
334 164 : Long_t nofPads = cluster.NofPads(kStatusToTest);
335 164 : Int_t nInX = AliMp::PairFirst(nofPads);
336 164 : Int_t nInY = AliMp::PairSecond(nofPads);
337 :
338 164 : if (fDebug) {
339 : Int_t npadOK = 0;
340 0 : for (Int_t j = 0; j < cluster.Multiplicity(); ++j) {
341 0 : AliMUONPad *pad = cluster.Pad(j);
342 : //if (pad->Status() == 1 && !pad->IsSaturated()) npadOK++;
343 0 : if (pad->Status() == AliMUONClusterFinderMLEM::GetUseForFitFlag() && !pad->IsSaturated()) npadOK++;
344 : }
345 0 : cout << " Number of pads to fit: " << npadOK << endl;
346 0 : cout << " nInX and Y: " << nInX << " " << nInY << endl;
347 0 : }
348 :
349 : Int_t nfitMax = 3;
350 164 : nfitMax = TMath::Min (nfitMax, (npads + 1) / 3);
351 164 : if (nfitMax > 1) {
352 552 : if (((nInX < 3) && (nInY < 3)) || ((nInX == 3) && (nInY < 3)) || ((nInX < 3) && (nInY == 3))) nfitMax = 1; // not enough pads in each direction
353 : }
354 164 : if (nfit > nfitMax) nfit = nfitMax;
355 :
356 : // Take cluster maxima as fitting seeds
357 : TObjArray *pix;
358 : AliMUONPad *pixPtr;
359 : Int_t npxclu;
360 164 : Double_t cont, cmax = 0, xseed = 0, yseed = 0, errOk[8], qq = 0;
361 :
362 2952 : for ( int i = 0; i < 8; ++i ) errOk[i]=0.0;
363 :
364 164 : Double_t xyseed[3][2], qseed[3], xyCand[3][2] = {{0},{0}}, sigCand[3][2] = {{0},{0}};
365 :
366 694 : for (Int_t ifit = 1; ifit <= nfit0; ++ifit)
367 : {
368 : cmax = 0;
369 183 : pix = clusters[clustFit[ifit-1]];
370 183 : npxclu = pix->GetEntriesFast();
371 : //qq = 0;
372 2616 : for (Int_t clu = 0; clu < npxclu; ++clu)
373 : {
374 1125 : pixPtr = (AliMUONPad*) pix->UncheckedAt(clu);
375 1125 : cont = pixPtr->Charge();
376 1125 : fQtot += cont;
377 1125 : if (cont > cmax)
378 : {
379 : cmax = cont;
380 439 : xseed = pixPtr->Coord(0);
381 439 : yseed = pixPtr->Coord(1);
382 439 : }
383 1125 : qq += cont;
384 1125 : xyCand[0][0] += pixPtr->Coord(0) * cont;
385 1125 : xyCand[0][1] += pixPtr->Coord(1) * cont;
386 1125 : sigCand[0][0] += pixPtr->Coord(0) * pixPtr->Coord(0) * cont;
387 1125 : sigCand[0][1] += pixPtr->Coord(1) * pixPtr->Coord(1) * cont;
388 : }
389 183 : xyseed[ifit-1][0] = xseed;
390 183 : xyseed[ifit-1][1] = yseed;
391 183 : qseed[ifit-1] = cmax;
392 : } // for (Int_t ifit=1;
393 :
394 164 : xyCand[0][0] /= qq; // <x>
395 164 : xyCand[0][1] /= qq; // <y>
396 164 : sigCand[0][0] = sigCand[0][0]/qq - xyCand[0][0]*xyCand[0][0]; // <x^2> - <x>^2
397 492 : sigCand[0][0] = sigCand[0][0] > 0 ? TMath::Sqrt (sigCand[0][0]) : 0;
398 164 : sigCand[0][1] = sigCand[0][1]/qq - xyCand[0][1]*xyCand[0][1]; // <y^2> - <y>^2
399 492 : sigCand[0][1] = sigCand[0][1] > 0 ? TMath::Sqrt (sigCand[0][1]) : 0;
400 164 : if (fDebug) cout << xyCand[0][0] << " " << xyCand[0][1] << " " << sigCand[0][0] << " " << sigCand[0][1] << endl;
401 :
402 164 : Int_t nDof, maxSeed[3];//, nMax = 0;
403 :
404 164 : if ( nfit0 < 0 || nfit0 > 3 ) {
405 0 : AliErrorStream() << "Wrong nfit0 value: " << nfit0 << endl;
406 0 : return nfit;
407 : }
408 164 : TMath::Sort(nfit0, qseed, maxSeed, kTRUE); // in decreasing order
409 :
410 164 : Double_t step[3]={0.01,0.002,0.02}, fmin, chi2o = 9999, chi2n;
411 164 : Double_t *gin = 0, func0, func1, param[8]={0}, step0[8]={0};
412 164 : Double_t param0[2][8]={{0},{0}}, deriv[2][8]={{0},{0}};
413 164 : Double_t shift[8]={0}, stepMax, derMax, parmin[8]={0}, parmax[8]={0}, func2[2]={0}, shift0;
414 164 : Double_t delta[8]={0}, scMax, dder[8], estim, shiftSave = 0;
415 : Int_t min, max, nCall = 0, nLoop, idMax = 0, nFail;
416 164 : Double_t rad, dist[3] = {0};
417 :
418 : // Try to fit with one-track hypothesis, then 2-track. If chi2/dof is
419 : // lower, try 3-track (if number of pads is sufficient).
420 : Int_t iflag = 0; // for the first call of fcn1
421 358 : for (Int_t iseed = 0; iseed < nfit; ++iseed)
422 : {
423 :
424 164 : Int_t memory[8] = {0};
425 164 : if (iseed)
426 : {
427 0 : for (Int_t j = 0; j < fNpar; ++j)
428 : {
429 0 : param[j] = parOk[j];
430 : }
431 0 : param[fNpar] = 0.6;
432 0 : parmin[fNpar] = 1E-9;
433 0 : parmax[fNpar++] = 1;
434 0 : }
435 :
436 164 : if (nfit == 1)
437 : {
438 149 : param[fNpar] = xyCand[0][0]; // take COG
439 149 : }
440 : else
441 : {
442 15 : param[fNpar] = xyseed[maxSeed[iseed]][0];
443 : //param[fNpar] = fNpar==0 ? -16.1651 : -15.2761;
444 : }
445 164 : parmin[fNpar] = xmin;
446 164 : parmax[fNpar++] = xmax;
447 164 : if (nfit == 1)
448 : {
449 149 : param[fNpar] = xyCand[0][1]; // take COG
450 149 : }
451 : else
452 : {
453 15 : param[fNpar] = xyseed[maxSeed[iseed]][1];
454 : //param[fNpar] = fNpar==1 ? -15.1737 : -15.8487;
455 : }
456 164 : parmin[fNpar] = ymin;
457 164 : parmax[fNpar++] = ymax;
458 :
459 984 : for (Int_t j = 0; j < fNpar; ++j)
460 : {
461 328 : step0[j] = shift[j] = step[j%3];
462 : }
463 :
464 164 : if (iseed)
465 : {
466 0 : for (Int_t j = 0; j < fNpar; ++j)
467 : {
468 0 : param0[1][j] = 0;
469 : }
470 0 : }
471 164 : if (fDebug) {
472 0 : for (Int_t j = 0; j < fNpar; ++j) cout << param[j] << " ";
473 0 : cout << endl;
474 0 : }
475 :
476 : // Try new algorithm
477 164 : min = nLoop = 1; stepMax = func2[1] = derMax = 999999; nFail = 0;
478 :
479 164 : while (1)
480 : {
481 785 : max = !min;
482 785 : Fcn1(cluster,fNpar, gin, func0, param, iflag); nCall++;
483 : iflag = 1;
484 : //cout << " Func: " << func0 << endl;
485 :
486 785 : func2[max] = func0;
487 4710 : for (Int_t j = 0; j < fNpar; ++j)
488 : {
489 1570 : param0[max][j] = param[j];
490 1570 : delta[j] = step0[j];
491 1570 : param[j] += delta[j] / 10;
492 2355 : if (j > 0) param[j-1] -= delta[j-1] / 10;
493 1570 : Fcn1(cluster,fNpar, gin, func1, param, iflag); nCall++;
494 1570 : deriv[max][j] = (func1 - func0) / delta[j] * 10; // first derivative
495 : //cout << j << " " << deriv[max][j] << endl;
496 6248 : dder[j] = param0[0][j] != param0[1][j] ? (deriv[0][j] - deriv[1][j]) /
497 1554 : (param0[0][j] - param0[1][j]) : 0; // second derivative
498 : }
499 785 : param[fNpar-1] -= delta[fNpar-1] / 10;
500 785 : if (nCall > 2000) break;
501 :
502 785 : min = func2[0] < func2[1] ? 0 : 1;
503 785 : nFail = min == max ? 0 : nFail + 1;
504 :
505 : stepMax = derMax = estim = 0;
506 4710 : for (Int_t j = 0; j < fNpar; ++j)
507 : {
508 : // Estimated distance to minimum
509 1570 : shift0 = shift[j];
510 1570 : if (nLoop == 1)
511 : {
512 328 : shift[j] = TMath::Sign (step0[j], -deriv[max][j]); // first step
513 328 : }
514 1274 : else if (TMath::Abs(deriv[0][j]) < 1.e-3 && TMath::Abs(deriv[1][j]) < 1.e-3)
515 : {
516 22 : shift[j] = 0;
517 22 : }
518 3234 : else if (((deriv[min][j]*deriv[!min][j] > 0) && (TMath::Abs(deriv[min][j]) > TMath::Abs(deriv[!min][j])))
519 3158 : || (TMath::Abs(deriv[0][j]-deriv[1][j]) < 1.e-3) || (TMath::Abs(dder[j]) < 1.e-6))
520 : {
521 80 : shift[j] = -TMath::Sign (shift[j], (func2[0]-func2[1]) * (param0[0][j]-param0[1][j]));
522 80 : if (min == max)
523 : {
524 66 : if (memory[j] > 1)
525 : {
526 14 : shift[j] *= 2;
527 14 : }
528 66 : memory[j]++;
529 66 : }
530 : }
531 : else
532 : {
533 3420 : shift[j] = dder[j] != 0 ? -deriv[min][j] / dder[j] : 0;
534 1140 : memory[j] = 0;
535 : }
536 :
537 1570 : Double_t es = TMath::Abs(shift[j]) / step0[j];
538 1570 : if (es > estim)
539 : {
540 : estim = es;
541 1208 : }
542 :
543 : // Too big step
544 1842 : if (TMath::Abs(shift[j])/step0[j] > 10) shift[j] = TMath::Sign(10.,shift[j]) * step0[j]; //
545 :
546 : // Failed to improve minimum
547 1570 : if (min != max)
548 : {
549 42 : memory[j] = 0;
550 42 : param[j] = param0[min][j];
551 84 : if (TMath::Abs(shift[j]+shift0) > 0.1*step0[j])
552 : {
553 70 : shift[j] = (shift[j] + shift0) / 2;
554 28 : }
555 : else
556 : {
557 14 : shift[j] /= -2;
558 : }
559 : }
560 :
561 : // Too big step
562 1570 : if (TMath::Abs(shift[j]*deriv[min][j]) > func2[min])
563 : {
564 211 : shift[j] = TMath::Sign (func2[min]/deriv[min][j], shift[j]);
565 211 : }
566 :
567 : // Introduce step relaxation factor
568 1570 : if (memory[j] < 3)
569 : {
570 1556 : scMax = 1 + 4 / TMath::Max(nLoop/2.,1.);
571 3100 : if (TMath::Abs(shift0) > 0 && TMath::Abs(shift[j]/shift0) > scMax)
572 : {
573 199 : shift[j] = TMath::Sign (shift0*scMax, shift[j]);
574 199 : }
575 : }
576 1570 : param[j] += shift[j];
577 : // Check parameter limits
578 1570 : if (param[j] < parmin[j])
579 : {
580 0 : shift[j] = parmin[j] - param[j];
581 0 : param[j] = parmin[j];
582 0 : }
583 1570 : else if (param[j] > parmax[j])
584 : {
585 0 : shift[j] = parmax[j] - param[j];
586 0 : param[j] = parmax[j];
587 0 : }
588 : //cout << " xxx " << j << " " << shift[j] << " " << param[j] << endl;
589 1570 : stepMax = TMath::Max (stepMax, TMath::Abs(shift[j]/step0[j]));
590 1570 : if (TMath::Abs(deriv[min][j]) > derMax)
591 : {
592 : idMax = j;
593 1198 : derMax = TMath::Abs (deriv[min][j]);
594 1198 : }
595 : } // for (Int_t j=0; j<fNpar;
596 :
597 785 : if (((estim < 1) && (derMax < 2)) || nLoop > 150) break; // minimum was found
598 :
599 621 : nLoop++;
600 :
601 : // Check for small step
602 621 : if (shift[idMax] == 0)
603 : {
604 0 : shift[idMax] = step0[idMax]/10;
605 0 : param[idMax] += shift[idMax];
606 0 : continue;
607 : }
608 :
609 621 : if (!memory[idMax] && derMax > 0.5 && nLoop > 10)
610 : {
611 0 : if (dder[idMax] != 0 && TMath::Abs(deriv[min][idMax]/dder[idMax]/shift[idMax]) > 10)
612 : {
613 0 : if (min == max) dder[idMax] = -dder[idMax];
614 0 : shift[idMax] = -deriv[min][idMax] / dder[idMax] / 10;
615 0 : param[idMax] += shift[idMax];
616 0 : stepMax = TMath::Max (stepMax, TMath::Abs(shift[idMax])/step0[idMax]);
617 0 : if (min == max) shiftSave = shift[idMax];
618 : }
619 0 : if (nFail > 10)
620 : {
621 0 : param[idMax] -= shift[idMax];
622 0 : shift[idMax] = 4 * shiftSave * (gRandom->Rndm(0) - 0.5);
623 0 : param[idMax] += shift[idMax];
624 0 : }
625 : }
626 : } // while (1)
627 :
628 164 : fmin = func2[min];
629 :
630 164 : nDof = npads - fNpar + nVirtual;
631 164 : if (!nDof) nDof++;
632 164 : chi2n = fmin / nDof;
633 164 : if (fDebug) cout << " Chi2 " << chi2n << " " << fNpar << endl;
634 :
635 : //if (fNpar > 2) cout << param0[min][fNpar-3] << " " << chi2n * (1+TMath::Min(1-param0[min][fNpar-3],0.25)) << endl;
636 : //if (chi2n*1.2+1.e-6 > chi2o )
637 164 : if (fNpar > 2 && (chi2n > chi2o || ((iseed == nfit-1)
638 0 : && (chi2n * (1+TMath::Min(1-param0[min][fNpar-3],0.25)) > chi2o))))
639 0 : { fNpar -= 3; break; }
640 :
641 : // Save parameters and errors
642 :
643 164 : if (nInX == 1) {
644 : // One pad per direction
645 : //for (Int_t i=0; i<fNpar; ++i) if (i == 0 || i == 2 || i == 5) param0[min][i] = xPad;
646 16 : for (Int_t i=0; i<fNpar; ++i) if (i == 0 || i == 2 || i == 5)
647 2 : param0[min][i] = xyCand[0][0];
648 2 : }
649 164 : if (nInY == 1) {
650 : // One pad per direction
651 : //for (Int_t i=0; i<fNpar; ++i) if (i == 1 || i == 3 || i == 6) param0[min][i] = yPad;
652 0 : for (Int_t i=0; i<fNpar; ++i) if (i == 1 || i == 3 || i == 6)
653 0 : param0[min][i] = xyCand[0][1];
654 0 : }
655 :
656 : /*
657 : if (iseed > 0) {
658 : // Find distance to the nearest neighbour
659 : dist[0] = dist[1] = TMath::Sqrt ((param0[min][0]-param0[min][2])*
660 : (param0[min][0]-param0[min][2])
661 : +(param0[min][1]-param0[min][3])*
662 : (param0[min][1]-param0[min][3]));
663 : if (iseed > 1) {
664 : dist[2] = TMath::Sqrt ((param0[min][0]-param0[min][5])*
665 : (param0[min][0]-param0[min][5])
666 : +(param0[min][1]-param0[min][6])*
667 : (param0[min][1]-param0[min][6]));
668 : rad = TMath::Sqrt ((param0[min][2]-param0[min][5])*
669 : (param0[min][2]-param0[min][5])
670 : +(param0[min][3]-param0[min][6])*
671 : (param0[min][3]-param0[min][6]));
672 : if (dist[2] < dist[0]) dist[0] = dist[2];
673 : if (rad < dist[1]) dist[1] = rad;
674 : if (rad < dist[2]) dist[2] = rad;
675 : }
676 : cout << dist[0] << " " << dist[1] << " " << dist[2] << endl;
677 : if (dist[TMath::LocMin(iseed+1,dist)] < 1.) { fNpar -= 3; break; }
678 : }
679 : */
680 :
681 984 : for (Int_t i = 0; i < fNpar; ++i) {
682 328 : parOk[i] = param0[min][i];
683 : //errOk[i] = fmin;
684 328 : errOk[i] = chi2n;
685 : // Bounded params
686 328 : parOk[i] = TMath::Max (parOk[i], parmin[i]);
687 328 : parOk[i] = TMath::Min (parOk[i], parmax[i]);
688 : }
689 :
690 : chi2o = chi2n;
691 318 : if (fmin < 0.1) break; // !!!???
692 174 : } // for (Int_t iseed=0;
693 :
694 164 : if (fDebug) {
695 0 : for (Int_t i=0; i<fNpar; ++i) {
696 0 : if (i == 4 || i == 7) {
697 0 : if ((i == 7) || ((i == 4) && (fNpar < 7))) cout << parOk[i] << endl;
698 0 : else cout << parOk[i] * (1-parOk[7]) << endl;
699 : continue;
700 : }
701 0 : cout << parOk[i] << " " << errOk[i] << endl;
702 0 : }
703 0 : }
704 164 : nfit = (fNpar + 1) / 3;
705 164 : dist[0] = dist[1] = dist[2] = 0;
706 :
707 164 : if (nfit > 1) {
708 : // Find distance to the nearest neighbour
709 0 : dist[0] = dist[1] = TMath::Sqrt ((parOk[0]-parOk[2])*
710 : (parOk[0]-parOk[2])
711 0 : +(parOk[1]-parOk[3])*
712 : (parOk[1]-parOk[3]));
713 0 : if (nfit > 2) {
714 0 : dist[2] = TMath::Sqrt ((parOk[0]-parOk[5])*
715 : (parOk[0]-parOk[5])
716 0 : +(parOk[1]-parOk[6])*
717 : (parOk[1]-parOk[6]));
718 0 : rad = TMath::Sqrt ((parOk[2]-parOk[5])*
719 : (parOk[2]-parOk[5])
720 0 : +(parOk[3]-parOk[6])*
721 : (parOk[3]-parOk[6]));
722 0 : if (dist[2] < dist[0]) dist[0] = dist[2];
723 0 : if (rad < dist[1]) dist[1] = rad;
724 0 : if (rad < dist[2]) dist[2] = rad;
725 : }
726 : }
727 :
728 : Int_t indx;
729 :
730 : Double_t coef = 0;
731 286 : if (iSimple) fnCoupled = 0;
732 820 : for (Int_t j = 0; j < nfit; ++j) {
733 164 : indx = 3 * j;
734 164 : coef = Param2Coef(j, coef, parOk);
735 :
736 : //void AliMUONClusterFinderMLEM::AddRawCluster(Double_t x, Double_t y,
737 : // Double_t qTot, Double_t fmin,
738 : // Int_t nfit, Int_t *tracks,
739 : // Double_t /*sigx*/,
740 : // Double_t /*sigy*/,
741 : // Double_t /*dist*/)
742 :
743 164 : if ( coef*fQtot >= fLowestClusterCharge )
744 : {
745 : //AZ AliMUONCluster* cluster1 = new AliMUONCluster();
746 164 : AliMUONCluster* cluster1 = new AliMUONCluster(cluster);
747 :
748 164 : cluster1->SetCharge(coef*fQtot,coef*fQtot);
749 492 : cluster1->SetPosition(TVector2(parOk[indx],parOk[indx+1]),TVector2(sigCand[0][0],sigCand[0][1]));
750 : //cluster1->SetChi2(dist[TMath::LocMin(nfit,dist)]);
751 164 : Int_t idx = TMath::LocMin(nfit,dist);
752 164 : if ( idx < 0 || idx > 2 ) {
753 0 : AliErrorStream() << "Wrong index value: " << idx << endl;
754 0 : return nfit;
755 : }
756 164 : cluster1->SetChi2(dist[idx]);
757 :
758 : // FIXME: we miss some information in this cluster, as compared to
759 : // the original AddRawCluster code.
760 :
761 492 : AliDebug(2,Form("Adding RawCluster detElemId %4d mult %2d charge %5d (xl,yl)=(%9.6g,%9.6g)",
762 : fDetElemId,cluster1->Multiplicity(),(Int_t)cluster1->Charge(),
763 : cluster1->Position().X(),cluster1->Position().Y()));
764 :
765 164 : clusterList.Add(cluster1);
766 164 : }
767 : // AddRawCluster (parOk[indx], // double x
768 : // parOk[indx+1], // double y
769 : // coef*qTot, // double charge
770 : // errOk[indx], // double fmin
771 : // nfit0+10*nfit+100*nMax+10000*fnCoupled, // int nfit
772 : // tracks, // int* tracks
773 : // sigCand[0][0], // double sigx
774 : // sigCand[0][1], // double sigy
775 : // dist[TMath::LocMin(nfit,dist)] // double dist
776 : // );
777 : }
778 164 : return nfit;
779 492 : }
780 :
781 :
782 : //_____________________________________________________________________________
783 : void
784 : AliMUONClusterSplitterMLEM::Split(const AliMUONCluster& cluster,
785 : TH2 *mlem, Double_t *coef,
786 : TObjArray& clusterList)
787 : {
788 : /// The main steering function to work with clusters of pixels in anode
789 : /// plane (find clusters, decouple them from each other, merge them (if
790 : /// necessary), pick up coupled pads, call the fitting function)
791 :
792 84 : Int_t nx = mlem->GetNbinsX();
793 42 : Int_t ny = mlem->GetNbinsY();
794 42 : Int_t nPix = fPixArray->GetEntriesFast();
795 :
796 : Double_t cont;
797 42 : Int_t nclust = 0, indx, indx1, nxy = ny * nx;
798 42 : Bool_t *used = new Bool_t[nxy];
799 :
800 2266 : for (Int_t j = 0; j < nxy; ++j) used[j] = kFALSE;
801 :
802 42 : TObjArray *clusters[200]={0};
803 : TObjArray *pix;
804 :
805 : // Find clusters of histogram bins (easier to work in 2-D space)
806 598 : for (Int_t i = 1; i <= ny; ++i)
807 : {
808 2696 : for (Int_t j = 1; j <= nx; ++j)
809 : {
810 1091 : indx = (i-1)*nx + j - 1;
811 1091 : if (used[indx]) continue;
812 477 : cont = mlem->GetBinContent(mlem->GetBin(j,i));
813 477 : if (cont < fLowestPixelCharge) continue;
814 63 : pix = new TObjArray(20);
815 63 : used[indx] = 1;
816 63 : pix->Add(BinToPix(mlem,j,i));
817 63 : AddBin(mlem, i, j, 0, used, pix); // recursive call
818 63 : if (nclust >= 200) AliFatal(" Too many clusters !!!");
819 63 : clusters[nclust++] = pix;
820 63 : } // for (Int_t j=1; j<=nx; j++) {
821 : } // for (Int_t i=1; i<=ny;
822 42 : if (fDebug) cout << nclust << endl;
823 84 : delete [] used;
824 :
825 : // Compute couplings between clusters and clusters to pads
826 42 : Int_t npad = cluster.Multiplicity();
827 :
828 : // Exclude pads with overflows
829 : /*
830 : for (Int_t j = 0; j < npad; ++j)
831 : {
832 : AliMUONPad* pad = cluster.Pad(j);
833 : if ( pad->IsSaturated() )
834 : {
835 : pad->SetStatus(-5);
836 : }
837 : else
838 : {
839 : pad->SetStatus(0);
840 : }
841 : }
842 : */
843 :
844 : // Compute couplings of clusters to pads (including overflows)
845 42 : TMatrixD aijclupad(nclust,npad);
846 42 : aijclupad = 0;
847 : Int_t npxclu;
848 210 : for (Int_t iclust = 0; iclust < nclust; ++iclust)
849 : {
850 63 : pix = clusters[iclust];
851 63 : npxclu = pix->GetEntriesFast();
852 940 : for (Int_t i = 0; i < npxclu; ++i)
853 : {
854 407 : indx = fPixArray->IndexOf(pix->UncheckedAt(i));
855 10168 : for (Int_t j = 0; j < npad; ++j)
856 : {
857 : //AliMUONPad* pad = cluster.Pad(j);
858 : //if ( pad->Status() < 0 && pad->Status() != -5) continue;
859 4677 : if (coef[j*nPix+indx] < fgkCouplMin) continue;
860 6854 : aijclupad(iclust,j) += coef[j*nPix+indx];
861 3427 : }
862 : }
863 : }
864 :
865 : // Compute couplings between clusters (exclude overflows)
866 42 : TMatrixD aijcluclu(nclust,nclust);
867 42 : aijcluclu = 0;
868 210 : for (Int_t iclust = 0; iclust < nclust; ++iclust)
869 : {
870 1444 : for (Int_t j = 0; j < npad; ++j)
871 : {
872 : // Exclude overflows
873 : //if ( cluster.Pad(j)->Status() < 0) continue;
874 1318 : if ( cluster.Pad(j)->IsSaturated()) continue;
875 1318 : if (aijclupad(iclust,j) < fgkCouplMin) continue;
876 1600 : for (Int_t iclust1=iclust+1; iclust1<nclust; iclust1++)
877 : {
878 420 : if (aijclupad(iclust1,j) < fgkCouplMin) continue;
879 356 : aijcluclu(iclust,iclust1) +=
880 534 : TMath::Sqrt (aijclupad(iclust,j)*aijclupad(iclust1,j));
881 178 : }
882 590 : }
883 : }
884 210 : for (Int_t iclust = 0; iclust < nclust; ++iclust)
885 : {
886 184 : for (Int_t iclust1 = iclust+1; iclust1 < nclust; ++iclust1)
887 : {
888 87 : aijcluclu(iclust1,iclust) = aijcluclu(iclust,iclust1);
889 : }
890 : }
891 :
892 42 : if (fDebug && nclust > 1) aijcluclu.Print();
893 :
894 : // Find groups of coupled clusters
895 42 : used = new Bool_t[nclust];
896 210 : for (Int_t j = 0; j < nclust; ++j) used[j] = kFALSE;
897 :
898 84 : Int_t *clustNumb = new Int_t[nclust];
899 42 : Int_t nCoupled, nForFit, minGroup[3], clustFit[3], nfit = 0;
900 : //Double_t parOk[8];
901 42 : Double_t parOk[8] = {0}; //AZ
902 :
903 210 : for (Int_t igroup = 0; igroup < nclust; ++igroup)
904 : {
905 63 : if (used[igroup]) continue;
906 42 : used[igroup] = kTRUE;
907 42 : clustNumb[0] = igroup;
908 42 : nCoupled = 1;
909 : // Find group of coupled clusters
910 42 : AddCluster(igroup, nclust, aijcluclu, used, clustNumb, nCoupled); // recursive
911 :
912 42 : if (fDebug) {
913 0 : cout << " nCoupled: " << nCoupled << endl;
914 0 : for (Int_t i=0; i<nCoupled; ++i) cout << clustNumb[i] << " "; cout << endl;
915 : }
916 :
917 42 : fnCoupled = nCoupled;
918 :
919 172 : while (nCoupled > 0)
920 : {
921 44 : if (nCoupled < 4)
922 : {
923 : nForFit = nCoupled;
924 206 : for (Int_t i = 0; i < nCoupled; ++i) clustFit[i] = clustNumb[i];
925 42 : }
926 : else
927 : {
928 : // Too many coupled clusters to fit - try to decouple them
929 : // Find the lowest coupling of 1, 2, min(3,nLinks/2) pixels with
930 : // all the others in the group
931 16 : for (Int_t j = 0; j < 3; ++j) minGroup[j] = -1;
932 2 : Double_t coupl = MinGroupCoupl(nCoupled, clustNumb, aijcluclu, minGroup);
933 :
934 : // Flag clusters for fit
935 : nForFit = 0;
936 16 : while (nForFit < 3 && minGroup[nForFit] >= 0)
937 : {
938 6 : if (fDebug) cout << clustNumb[minGroup[nForFit]] << " ";
939 2 : clustFit[nForFit] = clustNumb[minGroup[nForFit]];
940 2 : clustNumb[minGroup[nForFit]] -= 999;
941 2 : nForFit++;
942 : }
943 2 : if (fDebug) cout << " nForFit " << nForFit << " " << coupl << endl;
944 : } // else
945 :
946 : // Select pads for fit.
947 88 : if (SelectPad(cluster,nCoupled, nForFit, clustNumb, clustFit, aijclupad) < 3 && nCoupled > 1)
948 : {
949 : // Deselect pads
950 40 : for (Int_t j = 0; j < npad; ++j)
951 : {
952 18 : AliMUONPad* pad = cluster.Pad(j);
953 : //if ( pad->Status()==1 ) pad->SetStatus(0);
954 : //if ( pad->Status()==-9) pad->SetStatus(-5);
955 36 : if ( pad->Status() == AliMUONClusterFinderMLEM::GetUseForFitFlag() ||
956 18 : pad->Status() == AliMUONClusterFinderMLEM::GetCoupledFlag())
957 10 : pad->SetStatus(AliMUONClusterFinderMLEM::GetZeroFlag());
958 : }
959 : // Merge the failed cluster candidates (with too few pads to fit) with
960 : // the one with the strongest coupling
961 2 : Merge(cluster,nForFit, nCoupled, clustNumb, clustFit, clusters, aijcluclu, aijclupad);
962 : }
963 : else
964 : {
965 : // Do the fit
966 42 : nfit = Fit(cluster,0, nForFit, clustFit, clusters, parOk, clusterList, mlem);
967 42 : if (nfit == 0) {
968 : //cout << " (nfit == 0) " << fNpar << " " << cluster.Multiplicity() << endl;
969 0 : fNpar = 0; // should be 0 by itself but just in case ...
970 0 : }
971 : }
972 :
973 : // Subtract the fitted charges from pads with strong coupling and/or
974 : // return pads for further use
975 44 : UpdatePads(cluster,nfit, parOk);
976 :
977 : // Mark used pads
978 1028 : for (Int_t j = 0; j < npad; ++j)
979 : {
980 470 : AliMUONPad* pad = cluster.Pad(j);
981 : //if ( pad->Status()==1 ) pad->SetStatus(-2);
982 : //if ( pad->Status()==-9) pad->SetStatus(-5);
983 470 : if ( pad->Status() == AliMUONClusterFinderMLEM::GetUseForFitFlag() )
984 444 : pad->SetStatus(AliMUONClusterFinderMLEM::GetModifiedFlag());
985 : }
986 :
987 : // Sort the clusters (move to the right the used ones)
988 44 : Int_t beg = 0, end = nCoupled - 1;
989 113 : while (beg < end)
990 : {
991 50 : if (clustNumb[beg] >= 0) { ++beg; continue; }
992 0 : for (Int_t j = end; j > beg; --j)
993 : {
994 0 : if (clustNumb[j] < 0) continue;
995 0 : end = j - 1;
996 0 : indx = clustNumb[beg];
997 0 : clustNumb[beg] = clustNumb[j];
998 0 : clustNumb[j] = indx;
999 0 : break;
1000 : }
1001 0 : ++beg;
1002 : }
1003 :
1004 44 : nCoupled -= nForFit;
1005 44 : if (nCoupled > 3)
1006 : {
1007 : // Remove couplings of used clusters
1008 0 : for (Int_t iclust = nCoupled; iclust < nCoupled+nForFit; ++iclust)
1009 : {
1010 0 : indx = clustNumb[iclust] + 999;
1011 0 : for (Int_t iclust1 = 0; iclust1 < nCoupled; ++iclust1)
1012 : {
1013 0 : indx1 = clustNumb[iclust1];
1014 0 : aijcluclu(indx,indx1) = aijcluclu(indx1,indx) = 0;
1015 : }
1016 : }
1017 :
1018 : // Update the remaining clusters couplings (subtract couplings from
1019 : // the used pads) - overflows excluded
1020 0 : for (Int_t j = 0; j < npad; ++j)
1021 : {
1022 0 : AliMUONPad* pad = cluster.Pad(j);
1023 : //if ( pad->Status() != -2) continue;
1024 0 : if ( pad->Status() != AliMUONClusterFinderMLEM::GetModifiedFlag()) continue;
1025 0 : for (Int_t iclust=0; iclust<nCoupled; ++iclust)
1026 : {
1027 0 : indx = clustNumb[iclust];
1028 0 : if (aijclupad(indx,j) < fgkCouplMin) continue;
1029 0 : for (Int_t iclust1 = iclust+1; iclust1 < nCoupled; ++iclust1)
1030 : {
1031 0 : indx1 = clustNumb[iclust1];
1032 0 : if (aijclupad(indx1,j) < fgkCouplMin) continue;
1033 : // Check this
1034 0 : aijcluclu(indx,indx1) -=
1035 0 : TMath::Sqrt (aijclupad(indx,j)*aijclupad(indx1,j));
1036 0 : aijcluclu(indx1,indx) = aijcluclu(indx,indx1);
1037 0 : }
1038 0 : }
1039 : //pad->SetStatus(-8);
1040 0 : pad->SetStatus(AliMUONClusterFinderMLEM::GetOverFlag());
1041 0 : } // for (Int_t j=0; j<npad;
1042 0 : } // if (nCoupled > 3)
1043 : } // while (nCoupled > 0)
1044 : } // for (Int_t igroup=0; igroup<nclust;
1045 :
1046 210 : for (Int_t iclust = 0; iclust < nclust; ++iclust)
1047 : {
1048 63 : pix = clusters[iclust];
1049 63 : pix->Clear();
1050 126 : delete pix;
1051 : }
1052 84 : delete [] clustNumb;
1053 84 : delete [] used;
1054 :
1055 42 : }
1056 :
1057 : //_____________________________________________________________________________
1058 : void
1059 : AliMUONClusterSplitterMLEM::Merge(const AliMUONCluster& cluster,
1060 : Int_t nForFit, Int_t nCoupled,
1061 : const Int_t *clustNumb, const Int_t *clustFit,
1062 : TObjArray **clusters,
1063 : TMatrixD& aijcluclu, TMatrixD& aijclupad)
1064 : {
1065 : /// Merge the group of clusters with the one having the strongest coupling with them
1066 :
1067 : Int_t indx, indx1, npxclu, imax=0;
1068 : TObjArray *pix, *pix1;
1069 : Double_t couplMax;
1070 :
1071 10 : for (Int_t icl = 0; icl < nForFit; ++icl)
1072 : {
1073 2 : indx = clustFit[icl];
1074 2 : pix = clusters[indx];
1075 2 : npxclu = pix->GetEntriesFast();
1076 : couplMax = -1;
1077 20 : for (Int_t icl1 = 0; icl1 < nCoupled; ++icl1)
1078 : {
1079 8 : indx1 = clustNumb[icl1];
1080 8 : if (indx1 < 0) continue;
1081 6 : if ( aijcluclu(indx,indx1) > couplMax)
1082 : {
1083 6 : couplMax = aijcluclu(indx,indx1);
1084 : imax = indx1;
1085 6 : }
1086 : } // for (Int_t icl1=0;
1087 : // Add to it
1088 2 : pix1 = clusters[imax];
1089 : // Add pixels
1090 8 : for (Int_t i = 0; i < npxclu; ++i)
1091 : {
1092 2 : pix1->Add(pix->UncheckedAt(i));
1093 2 : pix->RemoveAt(i);
1094 : }
1095 :
1096 : //Add cluster-to-cluster couplings
1097 20 : for (Int_t icl1 = 0; icl1 < nCoupled; ++icl1)
1098 : {
1099 8 : indx1 = clustNumb[icl1];
1100 14 : if (indx1 < 0 || indx1 == imax) continue;
1101 4 : aijcluclu(indx1,imax) += aijcluclu(indx,indx1);
1102 4 : aijcluclu(imax,indx1) = aijcluclu(indx1,imax);
1103 4 : }
1104 2 : aijcluclu(indx,imax) = aijcluclu(imax,indx) = 0;
1105 :
1106 : //Add cluster-to-pad couplings
1107 2 : Int_t mult = cluster.Multiplicity();
1108 40 : for (Int_t j = 0; j < mult; ++j)
1109 : {
1110 18 : AliMUONPad* pad = cluster.Pad(j);
1111 : //if ( pad->Status() < 0 && pad->Status() != -5 ) continue;// exclude used pads
1112 18 : if ( pad->Status() != AliMUONClusterFinderMLEM::GetZeroFlag()) continue;// exclude used pads
1113 18 : aijclupad(imax,j) += aijclupad(indx,j);
1114 18 : aijclupad(indx,j) = 0;
1115 18 : }
1116 : } // for (Int_t icl=0; icl<nForFit;
1117 2 : }
1118 :
1119 :
1120 : //_____________________________________________________________________________
1121 : Double_t
1122 : AliMUONClusterSplitterMLEM::MinGroupCoupl(Int_t nCoupled, const Int_t *clustNumb,
1123 : const TMatrixD& aijcluclu, Int_t *minGroup)
1124 : {
1125 : /// Find group of clusters with minimum coupling to all the others
1126 :
1127 4 : Int_t i123max = TMath::Min(3,nCoupled/2);
1128 : Int_t indx, indx1, indx2, indx3, nTot = 0;
1129 : Double_t *coupl1 = 0, *coupl2 = 0, *coupl3 = 0;
1130 :
1131 12 : for (Int_t i123 = 1; i123 <= i123max; ++i123) {
1132 :
1133 4 : if (i123 == 1) {
1134 2 : coupl1 = new Double_t [nCoupled];
1135 20 : for (Int_t i = 0; i < nCoupled; ++i) coupl1[i] = 0;
1136 2 : }
1137 2 : else if (i123 == 2) {
1138 2 : nTot = nCoupled*nCoupled;
1139 2 : coupl2 = new Double_t [nTot];
1140 68 : for (Int_t i = 0; i < nTot; ++i) coupl2[i] = 9999;
1141 2 : } else {
1142 0 : nTot = nTot*nCoupled;
1143 0 : coupl3 = new Double_t [nTot];
1144 0 : for (Int_t i = 0; i < nTot; ++i) coupl3[i] = 9999;
1145 : } // else
1146 :
1147 40 : for (Int_t i = 0; i < nCoupled; ++i) {
1148 16 : indx1 = clustNumb[i];
1149 80 : for (Int_t j = i+1; j < nCoupled; ++j) {
1150 24 : indx2 = clustNumb[j];
1151 24 : if (i123 == 1) {
1152 12 : coupl1[i] += aijcluclu(indx1,indx2);
1153 12 : coupl1[j] += aijcluclu(indx1,indx2);
1154 12 : }
1155 12 : else if (i123 == 2) {
1156 12 : indx = i*nCoupled + j;
1157 12 : coupl2[indx] = coupl1[i] + coupl1[j];
1158 12 : coupl2[indx] -= 2 * (aijcluclu(indx1,indx2));
1159 12 : } else {
1160 0 : for (Int_t k = j+1; k < nCoupled; ++k) {
1161 0 : indx3 = clustNumb[k];
1162 0 : indx = i*nCoupled*nCoupled + j*nCoupled + k;
1163 0 : coupl3[indx] = coupl2[i*nCoupled+j] + coupl1[k];
1164 0 : coupl3[indx] -= 2 * (aijcluclu(indx1,indx3)+aijcluclu(indx2,indx3));
1165 : }
1166 : } // else
1167 : } // for (Int_t j=i+1;
1168 : } // for (Int_t i=0;
1169 : } // for (Int_t i123=1;
1170 :
1171 : // Find minimum coupling
1172 : Double_t couplMin = 9999;
1173 : Int_t locMin = 0;
1174 :
1175 12 : for (Int_t i123 = 1; i123 <= i123max; ++i123) {
1176 4 : if (i123 == 1) {
1177 2 : locMin = TMath::LocMin(nCoupled, coupl1);
1178 2 : couplMin = coupl1[locMin];
1179 2 : minGroup[0] = locMin;
1180 4 : delete [] coupl1;
1181 : }
1182 2 : else if (i123 == 2) {
1183 2 : locMin = TMath::LocMin(nCoupled*nCoupled, coupl2);
1184 2 : if (coupl2[locMin] < couplMin) {
1185 : couplMin = coupl2[locMin];
1186 0 : minGroup[0] = locMin/nCoupled;
1187 0 : minGroup[1] = locMin%nCoupled;
1188 0 : }
1189 4 : delete [] coupl2;
1190 : } else {
1191 0 : locMin = TMath::LocMin(nTot, coupl3);
1192 0 : if (coupl3[locMin] < couplMin) {
1193 : couplMin = coupl3[locMin];
1194 0 : minGroup[0] = locMin/nCoupled/nCoupled;
1195 0 : minGroup[1] = locMin%(nCoupled*nCoupled)/nCoupled;
1196 0 : minGroup[2] = locMin%nCoupled;
1197 0 : }
1198 0 : delete [] coupl3;
1199 : } // else
1200 : } // for (Int_t i123=1;
1201 2 : return couplMin;
1202 : }
1203 :
1204 : //_____________________________________________________________________________
1205 : Int_t
1206 : AliMUONClusterSplitterMLEM::SelectPad(const AliMUONCluster& cluster,
1207 : Int_t nCoupled, Int_t nForFit,
1208 : const Int_t *clustNumb, const Int_t *clustFit,
1209 : const TMatrixD& aijclupad)
1210 : {
1211 : /// Select pads for fit. If too many coupled clusters, find pads giving
1212 : /// the strongest coupling with the rest of clusters and exclude them from the fit.
1213 :
1214 88 : Int_t npad = cluster.Multiplicity();
1215 : Double_t *padpix = 0;
1216 :
1217 44 : if (nCoupled > 3)
1218 : {
1219 2 : padpix = new Double_t[npad];
1220 40 : for (Int_t i = 0; i < npad; ++i) padpix[i] = 0.;
1221 2 : }
1222 :
1223 : Int_t nOK = 0, indx, indx1;
1224 214 : for (Int_t iclust = 0; iclust < nForFit; ++iclust)
1225 : {
1226 63 : indx = clustFit[iclust];
1227 1444 : for (Int_t j = 0; j < npad; ++j)
1228 : {
1229 659 : if ( aijclupad(indx,j) < fgkCouplMin) continue;
1230 590 : AliMUONPad* pad = cluster.Pad(j);
1231 : /*
1232 : if ( pad->Status() == -5 ) pad->SetStatus(-9); // flag overflow
1233 : if ( pad->Status() < 0 ) continue; // exclude overflows and used pads
1234 : if ( !pad->Status() )
1235 : {
1236 : pad->SetStatus(1);
1237 : ++nOK; // pad to be used in fit
1238 : }
1239 : */
1240 1044 : if ( pad->Status() != AliMUONClusterFinderMLEM::GetZeroFlag()
1241 1180 : || pad->IsSaturated() ) continue; // used pads and overflows
1242 454 : pad->SetStatus(AliMUONClusterFinderMLEM::GetUseForFitFlag());
1243 454 : ++nOK; // pad to be used in fit
1244 :
1245 454 : if (nCoupled > 3)
1246 : {
1247 : // Check other clusters
1248 100 : for (Int_t iclust1 = 0; iclust1 < nCoupled; ++iclust1)
1249 : {
1250 40 : indx1 = clustNumb[iclust1];
1251 40 : if (indx1 < 0) continue;
1252 30 : if ( aijclupad(indx1,j) < fgkCouplMin ) continue;
1253 22 : padpix[j] += aijclupad(indx1,j);
1254 22 : }
1255 10 : } // if (nCoupled > 3)
1256 454 : } // for (Int_t j=0; j<npad;
1257 : } // for (Int_t iclust=0; iclust<nForFit
1258 86 : if (nCoupled < 4) return nOK;
1259 :
1260 : Double_t aaa = 0;
1261 40 : for (Int_t j = 0; j < npad; ++j)
1262 : {
1263 18 : if (padpix[j] < fgkCouplMin) continue;
1264 10 : aaa += padpix[j];
1265 : //cluster.Pad(j)->SetStatus(-1); // exclude pads with strong coupling to the other clusters
1266 10 : cluster.Pad(j)->SetStatus(AliMUONClusterFinderMLEM::GetCoupledFlag()); // exclude pads with strong coupling to the other clusters
1267 10 : nOK--;
1268 10 : }
1269 4 : delete [] padpix;
1270 : return nOK;
1271 44 : }
1272 :
1273 : //_____________________________________________________________________________
1274 : void
1275 : AliMUONClusterSplitterMLEM::UpdatePads(const AliMUONCluster& cluster,
1276 : Int_t /*nfit*/, Double_t *par)
1277 : {
1278 : /// Subtract the fitted charges from pads with strong coupling
1279 :
1280 88 : Int_t indx, mult = cluster.Multiplicity(), iend = fNpar/3;
1281 : Double_t charge, coef=0;
1282 :
1283 1028 : for (Int_t j = 0; j < mult; ++j)
1284 : {
1285 470 : AliMUONPad* pad = cluster.Pad(j);
1286 : //if ( pad->Status() != -1 ) continue;
1287 940 : if ( pad->Status() != AliMUONClusterFinderMLEM::GetCoupledFlag() ) continue;
1288 0 : if (fNpar != 0)
1289 : {
1290 : charge = 0;
1291 0 : for (Int_t i = 0; i <= iend; ++i)
1292 : {
1293 : // sum over hits
1294 0 : indx = 3 * i;
1295 0 : coef = Param2Coef(i, coef, par);
1296 0 : charge += ChargeIntegration(par[indx],par[indx+1],*pad) * coef;
1297 : }
1298 0 : charge *= fQtot;
1299 0 : pad->SetCharge(pad->Charge()-charge);
1300 0 : } // if (fNpar != 0)
1301 :
1302 : //if (pad->Charge() > 6 /*fgkZeroSuppression*/) pad->SetStatus(0);
1303 0 : if (pad->Charge() > fLowestPadCharge) pad->SetStatus(AliMUONClusterFinderMLEM::GetZeroFlag());
1304 : // return pad for further using // FIXME: remove usage of zerosuppression here
1305 0 : else pad->SetStatus(AliMUONClusterFinderMLEM::GetOverFlag()); // do not use anymore
1306 :
1307 0 : } // for (Int_t j=0;
1308 44 : }
1309 :
1310 :
|
