1 /**************************************************************************
2 * Copyright(c) 1998-1999, ALICE Experiment at CERN, All rights reserved. *
4 * Author: The ALICE Off-line Project. *
5 * Contributors are mentioned in the code where appropriate. *
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 **************************************************************************/
21 #include <Riostream.h>
23 #include "AliMUONClusterFinderVS.h"
24 #include "AliMpPlaneType.h"
25 #include "AliMpVSegmentation.h"
26 #include "AliMUONDigit.h"
27 #include "AliMUONRawCluster.h"
28 #include "AliMUONGeometrySegmentation.h"
29 #include "AliMUONMathieson.h"
30 #include "AliMUONClusterInput.h"
31 #include "AliMUONDigitMapA1.h"
34 //_____________________________________________________________________
35 // This function is minimized in the double-Mathieson fit
36 void fcnS2(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
37 void fcnS1(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
38 void fcnCombiS1(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
39 void fcnCombiS2(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
41 ClassImp(AliMUONClusterFinderVS)
43 AliMUONClusterFinderVS::AliMUONClusterFinderVS()
46 // Default constructor
47 fInput=AliMUONClusterInput::Instance();
50 fTrack[0]=fTrack[1]=-1;
51 fGhostChi2Cut = 1e6; // nothing done by default
55 for(Int_t i=0; i<100; i++) {
56 for (Int_t j=0; j<2; j++) {
60 fRawClusters = new TClonesArray("AliMUONRawCluster",1000);
63 //____________________________________________________________________________
64 AliMUONClusterFinderVS::~AliMUONClusterFinderVS()
66 // Reset tracks information
69 fRawClusters->Delete();
74 AliMUONClusterFinderVS::AliMUONClusterFinderVS(const AliMUONClusterFinderVS & clusterFinder):TObject(clusterFinder)
76 // Protected copy constructor
78 AliFatal("Not implemented.");
80 //____________________________________________________________________________
81 void AliMUONClusterFinderVS::ResetRawClusters()
83 // Reset tracks information
85 if (fRawClusters) fRawClusters->Clear();
87 //____________________________________________________________________________
88 void AliMUONClusterFinderVS::Decluster(AliMUONRawCluster *cluster)
90 // Decluster by local maxima
91 SplitByLocalMaxima(cluster);
93 //____________________________________________________________________________
94 void AliMUONClusterFinderVS::SplitByLocalMaxima(AliMUONRawCluster *c)
96 // Split complex cluster by local maxima
99 fInput->SetCluster(c);
101 fMul[0]=c->GetMultiplicity(0);
102 fMul[1]=c->GetMultiplicity(1);
105 // dump digit information into arrays
110 for (cath=0; cath<2; cath++) {
113 for (i=0; i<fMul[cath]; i++) {
115 fDig[i][cath]=fInput->Digit(cath, c->GetIndex(i, cath));
117 fIx[i][cath]= fDig[i][cath]->PadX();
118 fIy[i][cath]= fDig[i][cath]->PadY();
120 fQ[i][cath] = fDig[i][cath]->Signal();
121 // pad centre coordinates
123 GetPadC(fInput->DetElemId(), fIx[i][cath], fIy[i][cath], fX[i][cath], fY[i][cath], fZ[i][cath]);
124 } // loop over cluster digits
126 } // loop over cathodes
132 // Initialise and perform mathieson fits
133 Float_t chi2, oldchi2;
134 // ++++++++++++++++++*************+++++++++++++++++++++
135 // (1) No more than one local maximum per cathode plane
136 // +++++++++++++++++++++++++++++++*************++++++++
137 if ((fNLocal[0]==1 && (fNLocal[1]==0 || fNLocal[1]==1)) ||
138 (fNLocal[0]==0 && fNLocal[1]==1)) {
139 // Perform combined single Mathieson fit
140 // Initial values for coordinates (x,y)
142 // One local maximum on cathodes 1 and 2 (X->cathode 2, Y->cathode 1)
143 if (fNLocal[0]==1 && fNLocal[1]==1) {
144 fXInit[0]=c->GetX(1);
145 fYInit[0]=c->GetY(0);
146 // One local maximum on cathode 1 (X,Y->cathode 1)
147 } else if (fNLocal[0]==1) {
148 fXInit[0]=c->GetX(0);
149 fYInit[0]=c->GetY(0);
150 // One local maximum on cathode 2 (X,Y->cathode 2)
152 fXInit[0]=c->GetX(1);
153 fYInit[0]=c->GetY(1);
155 AliDebug(1,"cas (1) CombiSingleMathiesonFit(c)");
156 chi2=CombiSingleMathiesonFit(c);
157 // Int_t ndf = fgNbins[0]+fgNbins[1]-2;
158 // Float_t prob = TMath::Prob(Double_t(chi2),ndf);
159 // prob1->Fill(prob);
160 // chi2_1->Fill(chi2);
162 AliDebug(1,Form(" chi2 %f ",chi2));
164 c->SetX(0, fXFit[0]);
165 c->SetY(0, fYFit[0]);
173 c->SetX(0, fSeg2[0]->GetAnod(fInput->DetElemId(), c->GetX(0)));
174 c->SetX(1, fSeg2[1]->GetAnod(fInput->DetElemId(), c->GetX(1)));
176 // c->SetDetElemId(fInput->DetElemId());
177 // If reasonable chi^2 add result to the list of rawclusters
180 // If not try combined double Mathieson Fit
182 AliDebug(1," MAUVAIS CHI2 !!!\n");
183 if (fNLocal[0]==1 && fNLocal[1]==1) {
184 fXInit[0]=fX[fIndLocal[0][1]][1];
185 fYInit[0]=fY[fIndLocal[0][0]][0];
186 fXInit[1]=fX[fIndLocal[0][1]][1];
187 fYInit[1]=fY[fIndLocal[0][0]][0];
188 } else if (fNLocal[0]==1) {
189 fXInit[0]=fX[fIndLocal[0][0]][0];
190 fYInit[0]=fY[fIndLocal[0][0]][0];
191 fXInit[1]=fX[fIndLocal[0][0]][0];
192 fYInit[1]=fY[fIndLocal[0][0]][0];
194 fXInit[0]=fX[fIndLocal[0][1]][1];
195 fYInit[0]=fY[fIndLocal[0][1]][1];
196 fXInit[1]=fX[fIndLocal[0][1]][1];
197 fYInit[1]=fY[fIndLocal[0][1]][1];
200 // Initial value for charge ratios
203 AliDebug(1,"\n cas (1) CombiDoubleMathiesonFit(c)\n");
204 chi2=CombiDoubleMathiesonFit(c);
205 // Int_t ndf = fgNbins[0]+fgNbins[1]-6;
206 // Float_t prob = TMath::Prob(chi2,ndf);
207 // prob2->Fill(prob);
208 // chi2_2->Fill(chi2);
210 // Was this any better ??
211 AliDebug(1,Form(" Old and new chi2 %f %f ", oldchi2, chi2));
212 if (fFitStat!=0 && chi2>0 && (2.*chi2 < oldchi2)) {
214 // Split cluster into two according to fit result
217 AliDebug(1,"Do not Split");
223 // +++++++++++++++++++++++++++++++++++++++
224 // (2) Two local maxima per cathode plane
225 // +++++++++++++++++++++++++++++++++++++++
226 } else if (fNLocal[0]==2 && fNLocal[1]==2) {
228 // Let's look for ghosts first
230 Float_t xm[4][2], ym[4][2];
231 Float_t dpx, dpy, dx, dy;
232 Int_t ixm[4][2], iym[4][2];
233 Int_t isec, im1, im2, ico;
235 // Form the 2x2 combinations
236 // 0-0, 0-1, 1-0, 1-1
238 for (im1=0; im1<2; im1++) {
239 for (im2=0; im2<2; im2++) {
240 xm[ico][0]=fX[fIndLocal[im1][0]][0];
241 ym[ico][0]=fY[fIndLocal[im1][0]][0];
242 xm[ico][1]=fX[fIndLocal[im2][1]][1];
243 ym[ico][1]=fY[fIndLocal[im2][1]][1];
245 ixm[ico][0]=fIx[fIndLocal[im1][0]][0];
246 iym[ico][0]=fIy[fIndLocal[im1][0]][0];
247 ixm[ico][1]=fIx[fIndLocal[im2][1]][1];
248 iym[ico][1]=fIy[fIndLocal[im2][1]][1];
252 // ico = 0 : first local maximum on cathodes 1 and 2
253 // ico = 1 : fisrt local maximum on cathode 1 and second on cathode 2
254 // ico = 2 : second local maximum on cathode 1 and first on cathode 1
255 // ico = 3 : second local maximum on cathodes 1 and 2
257 // Analyse the combinations and keep those that are possible !
258 // For each combination check consistency in x and y
261 Float_t dr[4] = {1.e4, 1.e4, 1.e4, 1.e4};
264 // In case of staggering maxima are displaced by exactly half the pad-size in y.
265 // We have to take into account the numerical precision in the consistency check;
268 for (ico=0; ico<4; ico++) {
269 accepted[ico]=kFALSE;
270 // cathode one: x-coordinate
271 isec=fSeg2[0]->Sector(fInput->DetElemId(), ixm[ico][0], iym[ico][0]);
272 dpx=fSeg2[0]->Dpx(fInput->DetElemId(), isec)/2.;
274 dx=TMath::Abs(xm[ico][0]-xm[ico][1]);
275 // cathode two: y-coordinate
277 isec=fSeg2[1]->Sector(fInput->DetElemId(), ixm[ico][1], iym[ico][1]);
278 dpy=fSeg2[1]->Dpy(fInput->DetElemId(), isec)/2.;
280 dy=TMath::Abs(ym[ico][0]-ym[ico][1]);
281 AliDebug(2,Form("\n %i %f %f %f %f %f %f \n", ico, ym[ico][0], ym[ico][1], dy, dpy, dx, dpx ));
282 if ((dx <= dpx) && (dy <= dpy+eps)) {
285 dr[ico] = TMath::Sqrt(dx*dx+dy*dy);
289 accepted[ico]=kFALSE;
292 AliDebug(1,Form("\n iacc= %d:\n", iacc));
294 if (accepted[0] && accepted[1]) {
295 if (dr[0] >= dr[1]) {
302 if (accepted[2] && accepted[3]) {
303 if (dr[2] >= dr[3]) {
310 // eliminate one candidate
314 for (ico=0; ico<4; ico++) {
315 if (accepted[ico] && dr[ico] > drmax) {
321 accepted[icobad] = kFALSE;
327 AliDebug(1,Form("\n iacc= %d:\n", iacc));
329 AliDebug(1,"\n iacc=2: No problem ! \n");
330 } else if (iacc==4) {
331 AliDebug(1,"\n iacc=4: Ok, but ghost problem !!! \n");
332 } else if (iacc==0) {
333 AliDebug(1,"\n iacc=0: I don't know what to do with this !!!!!!!!! \n");
336 // Initial value for charge ratios
337 fQrInit[0]=Float_t(fQ[fIndLocal[0][0]][0])/
338 Float_t(fQ[fIndLocal[0][0]][0]+fQ[fIndLocal[1][0]][0]);
339 fQrInit[1]=Float_t(fQ[fIndLocal[0][1]][1])/
340 Float_t(fQ[fIndLocal[0][1]][1]+fQ[fIndLocal[1][1]][1]);
342 // ******* iacc = 0 *******
343 // No combinations found between the 2 cathodes
344 // We keep the center of gravity of the cluster
349 // ******* iacc = 1 *******
350 // Only one combination found between the 2 cathodes
352 // Initial values for the 2 maxima (x,y)
354 // 1 maximum is initialised with the maximum of the combination found (X->cathode 2, Y->cathode 1)
355 // 1 maximum is initialised with the other maximum of the first cathode
362 } else if (accepted[1]){
368 } else if (accepted[2]){
374 } else if (accepted[3]){
381 AliDebug(1,"cas (2) CombiDoubleMathiesonFit(c)");
382 chi2=CombiDoubleMathiesonFit(c);
383 // Int_t ndf = fgNbins[0]+fgNbins[1]-6;
384 // Float_t prob = TMath::Prob(chi2,ndf);
385 // prob2->Fill(prob);
386 // chi2_2->Fill(chi2);
387 AliDebug(1,Form(" chi2 %f\n",chi2));
389 // If reasonable chi^2 add result to the list of rawclusters
394 // 1 maximum is initialised with the maximum of the combination found (X->cathode 2, Y->cathode 1)
395 // 1 maximum is initialised with the other maximum of the second cathode
402 } else if (accepted[1]){
408 } else if (accepted[2]){
414 } else if (accepted[3]){
421 AliDebug(1,"\n cas (2) CombiDoubleMathiesonFit(c)\n");
422 chi2=CombiDoubleMathiesonFit(c);
423 // Int_t ndf = fgNbins[0]+fgNbins[1]-6;
424 // Float_t prob = TMath::Prob(chi2,ndf);
425 // prob2->Fill(prob);
426 // chi2_2->Fill(chi2);
427 AliDebug(1,Form(" chi2 %f\n",chi2));
429 // If reasonable chi^2 add result to the list of rawclusters
433 //We keep only the combination found (X->cathode 2, Y->cathode 1)
434 for (Int_t ico=0; ico<2; ico++) {
436 AliMUONRawCluster cnew;
438 for (cath=0; cath<2; cath++) {
439 cnew.SetX(cath, Float_t(xm[ico][1]));
440 cnew.SetY(cath, Float_t(ym[ico][0]));
441 cnew.SetZ(cath, fZPlane);
442 cnew.SetMultiplicity(cath,c->GetMultiplicity(cath));
443 for (i=0; i<fMul[cath]; i++) {
444 cnew.SetIndex(i, cath, c->GetIndex(i,cath));
445 fSeg2[cath]->SetPad(fInput->DetElemId(), fIx[i][cath], fIy[i][cath]);
447 AliDebug(1,Form("\nRawCluster %d cath %d\n",ico,cath));
448 AliDebug(1,Form("mult_av %d\n",c->GetMultiplicity(cath)));
449 FillCluster(&cnew,cath);
451 cnew.SetClusterType(cnew.PhysicsContribution());
460 // ******* iacc = 2 *******
461 // Two combinations found between the 2 cathodes
463 // Was the same maximum taken twice
464 if ((accepted[0]&&accepted[1]) || (accepted[2]&&accepted[3])) {
465 AliDebug(1,"\n Maximum taken twice !!!\n");
467 // Have a try !! with that
468 if (accepted[0]&&accepted[3]) {
479 AliDebug(1,"\n cas (2) CombiDoubleMathiesonFit(c)\n");
480 chi2=CombiDoubleMathiesonFit(c);
481 // Int_t ndf = fgNbins[0]+fgNbins[1]-6;
482 // Float_t prob = TMath::Prob(chi2,ndf);
483 // prob2->Fill(prob);
484 // chi2_2->Fill(chi2);
488 // No ghosts ! No Problems ! - Perform one fit only !
489 if (accepted[0]&&accepted[3]) {
500 AliDebug(1,"\n cas (2) CombiDoubleMathiesonFit(c)\n");
501 chi2=CombiDoubleMathiesonFit(c);
502 // Int_t ndf = fgNbins[0]+fgNbins[1]-6;
503 // Float_t prob = TMath::Prob(chi2,ndf);
504 // prob2->Fill(prob);
505 // chi2_2->Fill(chi2);
506 AliDebug(1,Form(" chi2 %f\n",chi2));
510 // ******* iacc = 4 *******
511 // Four combinations found between the 2 cathodes
513 } else if (iacc==4) {
514 // Perform fits for the two possibilities !!
515 // Accept if charges are compatible on both cathodes
516 // If none are compatible, keep everything
521 AliDebug(1,"\n cas (2) CombiDoubleMathiesonFit(c)\n");
522 chi2=CombiDoubleMathiesonFit(c);
523 // Int_t ndf = fgNbins[0]+fgNbins[1]-6;
524 // Float_t prob = TMath::Prob(chi2,ndf);
525 // prob2->Fill(prob);
526 // chi2_2->Fill(chi2);
527 AliDebug(1,Form(" chi2 %f\n",chi2));
528 // store results of fit and postpone decision
529 Double_t sXFit[2],sYFit[2],sQrFit[2];
531 for (Int_t i=0;i<2;i++) {
541 AliDebug(1,"\n cas (2) CombiDoubleMathiesonFit(c)\n");
542 chi2=CombiDoubleMathiesonFit(c);
543 // ndf = fgNbins[0]+fgNbins[1]-6;
544 // prob = TMath::Prob(chi2,ndf);
545 // prob2->Fill(prob);
546 // chi2_2->Fill(chi2);
547 AliDebug(1,Form(" chi2 %f\n",chi2));
548 // We have all informations to perform the decision
549 // Compute the chi2 for the 2 possibilities
550 Float_t chi2fi,chi2si,chi2f,chi2s;
552 chi2f = (TMath::Log(fInput->TotalCharge(0)*fQrFit[0]
553 / (fInput->TotalCharge(1)*fQrFit[1]) )
554 / fInput->ChargeCorrel() );
556 chi2fi = (TMath::Log(fInput->TotalCharge(0)*(1-fQrFit[0])
557 / (fInput->TotalCharge(1)*(1-fQrFit[1])) )
558 / fInput->ChargeCorrel() );
559 chi2f += chi2fi*chi2fi;
561 chi2s = (TMath::Log(fInput->TotalCharge(0)*sQrFit[0]
562 / (fInput->TotalCharge(1)*sQrFit[1]) )
563 / fInput->ChargeCorrel() );
565 chi2si = (TMath::Log(fInput->TotalCharge(0)*(1-sQrFit[0])
566 / (fInput->TotalCharge(1)*(1-sQrFit[1])) )
567 / fInput->ChargeCorrel() );
568 chi2s += chi2si*chi2si;
570 // usefull to store the charge matching chi2 in the cluster
571 // fChi2[0]=sChi2[1]=chi2f;
572 // fChi2[1]=sChi2[0]=chi2s;
574 if (chi2f<=fGhostChi2Cut && chi2s<=fGhostChi2Cut)
576 if (chi2f>fGhostChi2Cut && chi2s>fGhostChi2Cut) {
582 if (chi2f<=fGhostChi2Cut)
584 if (chi2s<=fGhostChi2Cut) {
585 // retreive saved values
586 for (Int_t i=0;i<2;i++) {
597 } else if (fNLocal[0]==2 && fNLocal[1]==1) {
598 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
599 // (3) Two local maxima on cathode 1 and one maximum on cathode 2
600 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
602 Float_t xm[4][2], ym[4][2];
603 Float_t dpx, dpy, dx, dy;
604 Int_t ixm[4][2], iym[4][2];
605 Int_t isec, im1, ico;
607 // Form the 2x2 combinations
608 // 0-0, 0-1, 1-0, 1-1
610 for (im1=0; im1<2; im1++) {
611 xm[ico][0]=fX[fIndLocal[im1][0]][0];
612 ym[ico][0]=fY[fIndLocal[im1][0]][0];
613 xm[ico][1]=fX[fIndLocal[0][1]][1];
614 ym[ico][1]=fY[fIndLocal[0][1]][1];
616 ixm[ico][0]=fIx[fIndLocal[im1][0]][0];
617 iym[ico][0]=fIy[fIndLocal[im1][0]][0];
618 ixm[ico][1]=fIx[fIndLocal[0][1]][1];
619 iym[ico][1]=fIy[fIndLocal[0][1]][1];
622 // ico = 0 : first local maximum on cathodes 1 and 2
623 // ico = 1 : second local maximum on cathode 1 and first on cathode 2
625 // Analyse the combinations and keep those that are possible !
626 // For each combination check consistency in x and y
630 // In case of staggering maxima are displaced by exactly half the pad-size in y.
631 // We have to take into account the numerical precision in the consistency check;
635 for (ico=0; ico<2; ico++) {
636 isec=fSeg2[0]->Sector(fInput->DetElemId(), ixm[ico][0], iym[ico][0]);
637 dpx=fSeg2[0]->Dpx(fInput->DetElemId(), isec)/2.;
639 dx=TMath::Abs(xm[ico][0]-xm[ico][1]);
640 isec=fSeg2[1]->Sector(fInput->DetElemId(), ixm[ico][1], iym[ico][1]);
641 dpy=fSeg2[1]->Dpy(fInput->DetElemId(), isec)/2.;
643 dy=TMath::Abs(ym[ico][0]-ym[ico][1]);
644 AliDebug(2,Form("\n %i %f %f %f %f \n", ico, ym[ico][0], ym[ico][1], dy, dpy ));
645 if ((dx <= dpx) && (dy <= dpy+eps)) {
651 accepted[ico]=kFALSE;
659 // Initial value for charge ratios
660 fQrInit[0]=Float_t(fQ[fIndLocal[0][0]][0])/
661 Float_t(fQ[fIndLocal[0][0]][0]+fQ[fIndLocal[1][0]][0]);
662 fQrInit[1]=fQrInit[0];
664 if (accepted[0] && accepted[1]) {
666 fXInit[0]=0.5*(xm[0][1]+xm[0][0]);
668 fXInit[1]=0.5*(xm[0][1]+xm[1][0]);
672 chi23=CombiDoubleMathiesonFit(c);
681 } else if (accepted[0]) {
686 chi21=CombiDoubleMathiesonFit(c);
687 // Int_t ndf = fgNbins[0]+fgNbins[1]-6;
688 // Float_t prob = TMath::Prob(chi2,ndf);
689 // prob2->Fill(prob);
690 // chi2_2->Fill(chi21);
691 AliDebug(1,Form(" chi2 %f\n",chi21));
692 if (chi21<10) Split(c);
693 } else if (accepted[1]) {
698 chi22=CombiDoubleMathiesonFit(c);
699 // Int_t ndf = fgNbins[0]+fgNbins[1]-6;
700 // Float_t prob = TMath::Prob(chi2,ndf);
701 // prob2->Fill(prob);
702 // chi2_2->Fill(chi22);
703 AliDebug(1,Form(" chi2 %f\n",chi22));
704 if (chi22<10) Split(c);
707 if (chi21 > 10 && chi22 > 10 && chi23 > 10) {
708 // We keep only the combination found (X->cathode 2, Y->cathode 1)
709 for (Int_t ico=0; ico<2; ico++) {
711 AliMUONRawCluster cnew;
713 for (cath=0; cath<2; cath++) {
714 cnew.SetX(cath, Float_t(xm[ico][1]));
715 cnew.SetY(cath, Float_t(ym[ico][0]));
716 cnew.SetZ(cath, fZPlane);
717 cnew.SetMultiplicity(cath, c->GetMultiplicity(cath));
718 for (i=0; i<fMul[cath]; i++) {
719 cnew.SetIndex(i, cath, c->GetIndex(i, cath));
720 fSeg2[cath]->SetPad(fInput->DetElemId(), fIx[i][cath], fIy[i][cath]);
723 AliDebug(1,Form("\nRawCluster %d cath %d\n",ico,cath));
724 AliDebug(1,Form("mult_av %d\n",c->GetMultiplicity(cath)));
726 FillCluster(&cnew,cath);
728 cnew.SetClusterType(cnew.PhysicsContribution());
735 // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
736 // (3') One local maximum on cathode 1 and two maxima on cathode 2
737 // +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
738 } else if (fNLocal[0]==1 && fNLocal[1]==2) {
739 Float_t xm[4][2], ym[4][2];
740 Float_t dpx, dpy, dx, dy;
741 Int_t ixm[4][2], iym[4][2];
742 Int_t isec, im1, ico;
744 // Form the 2x2 combinations
745 // 0-0, 0-1, 1-0, 1-1
747 for (im1=0; im1<2; im1++) {
748 xm[ico][0]=fX[fIndLocal[0][0]][0];
749 ym[ico][0]=fY[fIndLocal[0][0]][0];
750 xm[ico][1]=fX[fIndLocal[im1][1]][1];
751 ym[ico][1]=fY[fIndLocal[im1][1]][1];
753 ixm[ico][0]=fIx[fIndLocal[0][0]][0];
754 iym[ico][0]=fIy[fIndLocal[0][0]][0];
755 ixm[ico][1]=fIx[fIndLocal[im1][1]][1];
756 iym[ico][1]=fIy[fIndLocal[im1][1]][1];
759 // ico = 0 : first local maximum on cathodes 1 and 2
760 // ico = 1 : first local maximum on cathode 1 and second on cathode 2
762 // Analyse the combinations and keep those that are possible !
763 // For each combination check consistency in x and y
767 // In case of staggering maxima are displaced by exactly half the pad-size in y.
768 // We have to take into account the numerical precision in the consistency check;
772 for (ico=0; ico<2; ico++) {
773 accepted[ico]=kFALSE;
774 isec=fSeg2[0]->Sector(fInput->DetElemId(), ixm[ico][0], iym[ico][0]);
775 dpx=fSeg2[0]->Dpx(fInput->DetElemId(), isec)/2.;
777 dx=TMath::Abs(xm[ico][0]-xm[ico][1]);
778 isec=fSeg2[1]->Sector(fInput->DetElemId(), ixm[ico][1], iym[ico][1]);
779 dpy=fSeg2[1]->Dpy(fInput->DetElemId(), isec)/2.;
781 dy=TMath::Abs(ym[ico][0]-ym[ico][1]);
782 AliDebug(1,Form("\n %i %f %f %f %f \n", ico, ym[ico][0], ym[ico][1], dy, dpy ));
783 if ((dx <= dpx) && (dy <= dpy+eps)) {
786 AliDebug(1,Form("ico %d\n",ico));
790 accepted[ico]=kFALSE;
798 fQrInit[1]=Float_t(fQ[fIndLocal[0][1]][1])/
799 Float_t(fQ[fIndLocal[0][1]][1]+fQ[fIndLocal[1][1]][1]);
801 fQrInit[0]=fQrInit[1];
804 if (accepted[0] && accepted[1]) {
806 fYInit[0]=0.5*(ym[0][0]+ym[0][1]);
808 fYInit[1]=0.5*(ym[0][0]+ym[1][1]);
811 chi23=CombiDoubleMathiesonFit(c);
820 } else if (accepted[0]) {
825 chi21=CombiDoubleMathiesonFit(c);
826 // Int_t ndf = fgNbins[0]+fgNbins[1]-6;
827 // Float_t prob = TMath::Prob(chi2,ndf);
828 // prob2->Fill(prob);
829 // chi2_2->Fill(chi21);
830 AliDebug(1,Form(" chi2 %f\n",chi21));
831 if (chi21<10) Split(c);
832 } else if (accepted[1]) {
837 chi22=CombiDoubleMathiesonFit(c);
838 // Int_t ndf = fgNbins[0]+fgNbins[1]-6;
839 // Float_t prob = TMath::Prob(chi2,ndf);
840 // prob2->Fill(prob);
841 // chi2_2->Fill(chi22);
842 AliDebug(1,Form(" chi2 %f\n",chi22));
843 if (chi22<10) Split(c);
846 if (chi21 > 10 && chi22 > 10 && chi23 > 10) {
847 //We keep only the combination found (X->cathode 2, Y->cathode 1)
848 for (Int_t ico=0; ico<2; ico++) {
850 AliMUONRawCluster cnew;
852 for (cath=0; cath<2; cath++) {
853 cnew.SetX(cath, Float_t(xm[ico][1]));
854 cnew.SetY(cath, Float_t(ym[ico][0]));
855 cnew.SetZ(cath, fZPlane);
856 cnew.SetMultiplicity(cath, c->GetMultiplicity(cath));
857 for (i=0; i<fMul[cath]; i++) {
858 cnew.SetIndex(i, cath, c->GetIndex(i, cath));
859 fSeg2[cath]->SetPad(fInput->DetElemId(), fIx[i][cath], fIy[i][cath]);
861 AliDebug(1,Form("\nRawCluster %d cath %d\n",ico,cath));
862 AliDebug(1,Form("mult_av %d\n",c->GetMultiplicity(cath)));
863 FillCluster(&cnew,cath);
865 cnew.SetClusterType(cnew.PhysicsContribution());
872 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
873 // (4) At least three local maxima on cathode 1 or on cathode 2
874 // ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
875 } else if (fNLocal[0]>2 || fNLocal[1]>2) {
876 Int_t param = fNLocal[0]*fNLocal[1];
879 Float_t ** xm = new Float_t * [param];
880 for (ii=0; ii<param; ii++) xm[ii]=new Float_t [2];
881 Float_t ** ym = new Float_t * [param];
882 for (ii=0; ii<param; ii++) ym[ii]=new Float_t [2];
883 Int_t ** ixm = new Int_t * [param];
884 for (ii=0; ii<param; ii++) ixm[ii]=new Int_t [2];
885 Int_t ** iym = new Int_t * [param];
886 for (ii=0; ii<param; ii++) iym[ii]=new Int_t [2];
889 Float_t dpx, dpy, dx, dy;
892 for (Int_t im1=0; im1<fNLocal[0]; im1++) {
893 for (Int_t im2=0; im2<fNLocal[1]; im2++) {
894 xm[ico][0]=fX[fIndLocal[im1][0]][0];
895 ym[ico][0]=fY[fIndLocal[im1][0]][0];
896 xm[ico][1]=fX[fIndLocal[im2][1]][1];
897 ym[ico][1]=fY[fIndLocal[im2][1]][1];
899 ixm[ico][0]=fIx[fIndLocal[im1][0]][0];
900 iym[ico][0]=fIy[fIndLocal[im1][0]][0];
901 ixm[ico][1]=fIx[fIndLocal[im2][1]][1];
902 iym[ico][1]=fIy[fIndLocal[im2][1]][1];
908 AliDebug(1,Form("nIco %d\n",nIco));
909 for (ico=0; ico<nIco; ico++) {
910 AliDebug(1,Form("ico = %d\n",ico));
911 isec=fSeg2[0]->Sector(fInput->DetElemId(), ixm[ico][0], iym[ico][0]);
912 dpx=fSeg2[0]->Dpx(fInput->DetElemId(), isec)/2.;
914 dx=TMath::Abs(xm[ico][0]-xm[ico][1]);
915 isec=fSeg2[1]->Sector(fInput->DetElemId(), ixm[ico][1], iym[ico][1]);
916 dpy=fSeg2[1]->Dpy(fInput->DetElemId(), isec)/2.;
918 dy=TMath::Abs(ym[ico][0]-ym[ico][1]);
919 AliDebug(1,Form("dx %f dpx %f dy %f dpy %f\n",dx,dpx,dy,dpy));
920 AliDebug(1,Form(" X %f Y %f\n",xm[ico][1],ym[ico][0]));
921 if ((dx <= dpx) && (dy <= dpy)) {
924 AliMUONRawCluster cnew;
925 for (cath=0; cath<2; cath++) {
926 cnew.SetX(cath, Float_t(xm[ico][1]));
927 cnew.SetY(cath, Float_t(ym[ico][0]));
928 cnew.SetZ(cath, fZPlane);
929 cnew.SetMultiplicity(cath, c->GetMultiplicity(cath));
930 for (i=0; i<fMul[cath]; i++) {
931 cnew.SetIndex(i, cath, c->GetIndex(i, cath));
932 fSeg2[cath]->SetPad(fInput->DetElemId(), fIx[i][cath], fIy[i][cath]);
934 FillCluster(&cnew,cath);
936 cnew.SetClusterType(cnew.PhysicsContribution());
937 // cnew.SetDetElemId(fInput->DetElemId());
949 void AliMUONClusterFinderVS::FindLocalMaxima(AliMUONRawCluster* /*c*/)
951 // Find all local maxima of a cluster
952 AliDebug(1,"\n Find Local maxima !");
956 Int_t cath, cath1; // loops over cathodes
957 Int_t i; // loops over digits
958 Int_t j; // loops over cathodes
962 // counters for number of local maxima
963 fNLocal[0]=fNLocal[1]=0;
964 // flags digits as local maximum
965 Bool_t isLocal[100][2];
966 for (i=0; i<100;i++) {
967 isLocal[i][0]=isLocal[i][1]=kFALSE;
969 // number of next neighbours and arrays to store them
972 // loop over cathodes
973 for (cath=0; cath<2; cath++) {
974 // loop over cluster digits
975 for (i=0; i<fMul[cath]; i++) {
976 // get neighbours for that digit and assume that it is local maximum
980 fSeg2[cath]->Neighbours(fInput->DetElemId(), fIx[i][cath], fIy[i][cath], &nn, x, y);
982 isLocal[i][cath]=kTRUE;
983 isec = fSeg2[cath]->Sector(fInput->DetElemId(), fIx[i][cath], fIy[i][cath]);
984 a0 = fSeg2[cath]->Dpx(fInput->DetElemId(), isec)*fSeg2[cath]->Dpy(fInput->DetElemId(), isec);
986 // loop over next neighbours, if at least one neighbour has higher charger assumption
987 // digit is not local maximum
988 for (j=0; j<nn; j++) {
989 if (fDigitMap[cath]->TestHit(x[j], y[j])==kEmpty) continue;
990 digt=(AliMUONDigit*) fDigitMap[cath]->GetHit(x[j], y[j]);
992 isec=fSeg2[cath]->Sector(fInput->DetElemId(), x[j], y[j]);
993 a1 = fSeg2[cath]->Dpx(fInput->DetElemId(),isec)*fSeg2[cath]->Dpy(fInput->DetElemId(), isec);
995 if (digt->Signal()/a1 > fQ[i][cath]/a0) {
996 isLocal[i][cath]=kFALSE;
999 // handle special case of neighbouring pads with equal signal
1000 } else if (digt->Signal() == fQ[i][cath]) {
1001 if (fNLocal[cath]>0) {
1002 for (Int_t k=0; k<fNLocal[cath]; k++) {
1003 if (x[j]==fIx[fIndLocal[k][cath]][cath]
1004 && y[j]==fIy[fIndLocal[k][cath]][cath])
1006 isLocal[i][cath]=kFALSE;
1008 } // loop over local maxima
1009 } // are there already local maxima
1011 } // loop over next neighbours
1012 if (isLocal[i][cath]) {
1013 fIndLocal[fNLocal[cath]][cath]=i;
1016 } // loop over all digits
1017 } // loop over cathodes
1019 AliDebug(1,Form("\n Found %d %d %d %d local Maxima\n",
1020 fNLocal[0], fNLocal[1], fMul[0], fMul[1]));
1021 AliDebug(1,Form("\n Cathode 1 local Maxima %d Multiplicite %d\n",fNLocal[0], fMul[0]));
1022 AliDebug(1,Form(" Cathode 2 local Maxima %d Multiplicite %d\n",fNLocal[1], fMul[1]));
1027 if (fNLocal[1]==2 && (fNLocal[0]==1 || fNLocal[0]==0)) {
1028 Int_t iback=fNLocal[0];
1030 // Two local maxima on cathode 2 and one maximum on cathode 1
1031 // Look for local maxima considering up and down neighbours on the 1st cathode only
1033 // Loop over cluster digits
1037 for (i=0; i<fMul[cath]; i++) {
1038 isec=fSeg2[cath]->Sector(fInput->DetElemId(), fIx[i][cath],fIy[i][cath]);
1039 dpy=fSeg2[cath]->Dpy(fInput->DetElemId(), isec);
1040 dpx=fSeg2[cath]->Dpx(fInput->DetElemId(), isec);
1042 if (isLocal[i][cath]) continue;
1043 // Pad position should be consistent with position of local maxima on the opposite cathode
1044 if ((TMath::Abs(fX[i][cath]-fX[fIndLocal[0][cath1]][cath1]) > dpx/2.) &&
1045 (TMath::Abs(fX[i][cath]-fX[fIndLocal[1][cath1]][cath1]) > dpx/2.))
1048 // get neighbours for that digit and assume that it is local maximum
1049 isLocal[i][cath]=kTRUE;
1050 // compare signal to that on the two neighbours on the left and on the right
1051 // iNN counts the number of neighbours with signal, it should be 1 or 2
1055 for (fSeg2[cath]->FirstPad(fInput->DetElemId(), fX[i][cath], fY[i][cath], fZPlane, 0., dpy);
1056 fSeg2[cath]->MorePads(fInput->DetElemId());
1057 fSeg2[cath]->NextPad(fInput->DetElemId()))
1059 ix = fSeg2[cath]->Ix();
1060 iy = fSeg2[cath]->Iy();
1061 // skip the current pad
1062 if (iy == fIy[i][cath]) continue;
1064 if (fDigitMap[cath]->TestHit(ix, iy)!=kEmpty) {
1066 digt=(AliMUONDigit*) fDigitMap[cath]->GetHit(ix,iy);
1067 if (digt->Signal() > fQ[i][cath]) isLocal[i][cath]=kFALSE;
1069 } // Loop over pad neighbours in y
1071 if (isLocal[i][cath] && iNN>0) {
1072 fIndLocal[fNLocal[cath]][cath]=i;
1075 } // loop over all digits
1076 // if one additional maximum has been found we are happy
1077 // if more maxima have been found restore the previous situation
1078 AliDebug(1,Form("\n New search gives %d local maxima for cathode 1 \n",
1080 AliDebug(1,Form(" %d local maxima for cathode 2 \n",
1082 if (fNLocal[cath]>2) {
1083 fNLocal[cath]=iback;
1086 } // 1,2 local maxima
1088 if (fNLocal[0]==2 && (fNLocal[1]==1 || fNLocal[1]==0)) {
1089 Int_t iback=fNLocal[1];
1091 // Two local maxima on cathode 1 and one maximum on cathode 2
1092 // Look for local maxima considering left and right neighbours on the 2nd cathode only
1095 Float_t eps = 1.e-5;
1098 // Loop over cluster digits
1099 for (i=0; i<fMul[cath]; i++) {
1100 isec=fSeg2[cath]->Sector(fInput->DetElemId(), fIx[i][cath],fIy[i][cath]);
1101 dpx=fSeg2[cath]->Dpx(fInput->DetElemId(), isec);
1102 dpy=fSeg2[cath]->Dpy(fInput->DetElemId(), isec);
1105 if (isLocal[i][cath]) continue;
1106 // Pad position should be consistent with position of local maxima on the opposite cathode
1107 if ((TMath::Abs(fY[i][cath]-fY[fIndLocal[0][cath1]][cath1]) > dpy/2.+eps) &&
1108 (TMath::Abs(fY[i][cath]-fY[fIndLocal[1][cath1]][cath1]) > dpy/2.+eps))
1112 // get neighbours for that digit and assume that it is local maximum
1113 isLocal[i][cath]=kTRUE;
1114 // compare signal to that on the two neighbours on the left and on the right
1116 // iNN counts the number of neighbours with signal, it should be 1 or 2
1118 for (fSeg2[cath]->FirstPad(fInput->DetElemId(), fX[i][cath], fY[i][cath], fZPlane, dpx, 0.);
1119 fSeg2[cath]->MorePads(fInput->DetElemId());
1120 fSeg2[cath]->NextPad(fInput->DetElemId()))
1123 ix = fSeg2[cath]->Ix();
1124 iy = fSeg2[cath]->Iy();
1126 // skip the current pad
1127 if (ix == fIx[i][cath]) continue;
1129 if (fDigitMap[cath]->TestHit(ix, iy)!=kEmpty) {
1131 digt=(AliMUONDigit*) fDigitMap[cath]->GetHit(ix,iy);
1132 if (digt->Signal() > fQ[i][cath]) isLocal[i][cath]=kFALSE;
1134 } // Loop over pad neighbours in x
1136 if (isLocal[i][cath] && iNN>0) {
1137 fIndLocal[fNLocal[cath]][cath]=i;
1140 } // loop over all digits
1141 // if one additional maximum has been found we are happy
1142 // if more maxima have been found restore the previous situation
1143 AliDebug(1,Form("\n New search gives %d local maxima for cathode 1 \n",fNLocal[0]));
1144 AliDebug(1,Form("\n %d local maxima for cathode 2 \n",fNLocal[1]));
1145 AliDebug(1,Form("\n New search gives %d %d \n",fNLocal[0],fNLocal[1]));
1146 if (fNLocal[cath]>2) {
1147 fNLocal[cath]=iback;
1149 } // 2,1 local maxima
1153 void AliMUONClusterFinderVS::FillCluster(AliMUONRawCluster* c, Int_t flag, Int_t cath)
1156 // Completes cluster information starting from list of digits
1163 c->SetPeakSignal(cath,c->GetPeakSignal(0));
1165 c->SetPeakSignal(cath,0);
1172 c->SetCharge(cath,0);
1175 AliDebug(1,Form("\n fPeakSignal %d\n",c->GetPeakSignal(cath)));
1176 for (Int_t i=0; i<c->GetMultiplicity(cath); i++)
1178 dig= fInput->Digit(cath,c->GetIndex(i,cath));
1179 ix=dig->PadX()+c->GetOffset(i,cath);
1181 Int_t q=dig->Signal();
1182 if (!flag) q=Int_t(q*c->GetContrib(i,cath));
1183 // fprintf(stderr,"q %d c->fPeakSignal[ %d ] %d\n",q,cath,c->fPeakSignal[cath]);
1184 if (dig->Physics() >= dig->Signal()) {
1186 } else if (dig->Physics() == 0) {
1188 } else c->SetPhysics(i,1);
1191 AliDebug(2,Form("q %d c->fPeakSignal[cath] %d\n",q,c->GetPeakSignal(cath)));
1192 // peak signal and track list
1193 if (q>c->GetPeakSignal(cath)) {
1194 c->SetPeakSignal(cath, q);
1195 c->SetTrack(0,dig->Hit());
1196 c->SetTrack(1,dig->Track(0));
1197 c->SetTrack(2,dig->Track(1));
1198 // fprintf(stderr," c->fTracks[0] %d c->fTracks[1] %d\n",dig->fHit,dig->fTracks[0]);
1202 fSeg2[cath]->GetPadC(fInput->DetElemId(), ix, iy, x, y, z);
1206 c->AddCharge(cath, q);
1208 } // loop over digits
1209 AliDebug(1," fin du cluster c\n");
1213 c->SetX(cath, c->GetX(cath)/c->GetCharge(cath));
1215 c->SetX(cath, fSeg2[cath]->GetAnod(fInput->DetElemId(), c->GetX(cath)));
1216 c->SetY(cath, c->GetY(cath)/c->GetCharge(cath));
1218 // apply correction to the coordinate along the anode wire
1224 fSeg2[cath]->GetPadI(fInput->DetElemId(), x, y, fZPlane, ix, iy);
1225 fSeg2[cath]->GetPadC(fInput->DetElemId(), ix, iy, x, y, z);
1226 isec=fSeg2[cath]->Sector(fInput->DetElemId(), ix,iy);
1227 cogCorr = fSeg2[cath]->CorrFunc(fInput->DetElemId(), isec-1);
1232 yOnPad=(c->GetY(cath)-y)/fSeg2[cath]->Dpy(fInput->DetElemId(), isec);
1234 c->SetY(cath, c->GetY(cath)-cogCorr->Eval(yOnPad, 0, 0));
1235 // slat ID from digit
1241 void AliMUONClusterFinderVS::FillCluster(AliMUONRawCluster* c, Int_t cath)
1244 // Completes cluster information starting from list of digits
1254 Float_t xpad, ypad, zpad;
1257 for (Int_t i=0; i<c->GetMultiplicity(cath); i++)
1259 dig = fInput->Digit(cath,c->GetIndex(i,cath));
1261 GetPadC(fInput->DetElemId(),dig->PadX(),dig->PadY(),xpad,ypad, zpad);
1262 AliDebug(1,Form("x %f y %f cx %f cy %f\n",xpad,ypad,c->GetX(0),c->GetY(0)));
1263 dx = xpad - c->GetX(0);
1264 dy = ypad - c->GetY(0);
1265 dr = TMath::Sqrt(dx*dx+dy*dy);
1269 AliDebug(1,Form(" dr %f\n",dr));
1270 Int_t q=dig->Signal();
1271 if (dig->Physics() >= dig->Signal()) {
1273 } else if (dig->Physics() == 0) {
1275 } else c->SetPhysics(i,1);
1276 c->SetPeakSignal(cath,q);
1277 c->SetTrack(0,dig->Hit());
1278 c->SetTrack(1,dig->Track(0));
1279 c->SetTrack(2,dig->Track(1));
1281 AliDebug(1,Form(" c->fTracks[0] %d c->fTracks[1] %d\n",dig->Hit(),
1285 } // loop over digits
1287 // apply correction to the coordinate along the anode wire
1289 c->SetX(cath,fSeg2[cath]->GetAnod(fInput->DetElemId(), c->GetX(cath)));
1292 void AliMUONClusterFinderVS::FindCluster(Int_t i, Int_t j, Int_t cath, AliMUONRawCluster &c){
1296 // Find a super cluster on both cathodes
1299 // Add i,j as element of the cluster
1302 Int_t idx = fDigitMap[cath]->GetHitIndex(i,j);
1303 AliMUONDigit* dig = (AliMUONDigit*) fDigitMap[cath]->GetHit(i,j);
1304 Int_t q=dig->Signal();
1305 Int_t theX=dig->PadX();
1306 Int_t theY=dig->PadY();
1308 if (q > TMath::Abs(c.GetPeakSignal(0)) && q > TMath::Abs(c.GetPeakSignal(1))) {
1309 c.SetPeakSignal(cath,q);
1310 c.SetTrack(0,dig->Hit());
1311 c.SetTrack(1,dig->Track(0));
1312 c.SetTrack(2,dig->Track(1));
1316 // Make sure that list of digits is ordered
1318 Int_t mu=c.GetMultiplicity(cath);
1319 c.SetIndex(mu, cath, idx);
1321 if (dig->Physics() >= dig->Signal()) {
1323 } else if (dig->Physics() == 0) {
1325 } else c.SetPhysics(mu,1);
1329 for (Int_t ind = mu-1; ind >= 0; ind--) {
1330 Int_t ist=c.GetIndex(ind,cath);
1331 Int_t ql=fInput->Digit(cath, ist)->Signal();
1332 Int_t ix=fInput->Digit(cath, ist)->PadX();
1333 Int_t iy=fInput->Digit(cath, ist)->PadY();
1335 if (q>ql || (q==ql && theX > ix && theY < iy)) {
1336 c.SetIndex(ind, cath, idx);
1337 c.SetIndex(ind+1, cath, ist);
1345 c.SetMultiplicity(cath, c.GetMultiplicity(cath)+1);
1346 if (c.GetMultiplicity(cath) >= 50 ) {
1347 AliDebug(1,Form("FindCluster - multiplicity >50 %d \n",c.GetMultiplicity(0)));
1348 c.SetMultiplicity(cath, 49);
1351 // Prepare center of gravity calculation
1353 fSeg2[cath]->GetPadC(fInput->DetElemId(), i, j, x, y, z);
1356 c.AddCharge(cath,q);
1358 // Flag hit as "taken"
1359 fDigitMap[cath]->FlagHit(i,j);
1361 // Now look recursively for all neighbours and pad hit on opposite cathode
1363 // Loop over neighbours
1367 Int_t xList[10], yList[10];
1368 fSeg2[cath]->Neighbours(fInput->DetElemId(), i,j,&nn,xList,yList);
1369 for (Int_t in=0; in<nn; in++) {
1373 if (fDigitMap[cath]->TestHit(ix,iy)==kUnused) {
1374 AliDebug(2,Form("\n Neighbours %d %d %d", cath, ix, iy));
1375 FindCluster(ix, iy, cath, c);
1380 Int_t iXopp[50], iYopp[50];
1382 // Neighbours on opposite cathode
1383 // Take into account that several pads can overlap with the present pad
1385 isec=fSeg2[cath]->Sector(fInput->DetElemId(), i,j);
1392 dx = (fSeg2[cath]->Dpx(fInput->DetElemId(), isec))/2.;
1397 dy = (fSeg2[cath]->Dpy(fInput->DetElemId(), isec))/2;
1402 // loop over pad neighbours on opposite cathode
1403 for (fSeg2[iop]->FirstPad(fInput->DetElemId(), x, y, fZPlane, dx, dy);
1404 fSeg2[iop]->MorePads(fInput->DetElemId());
1405 fSeg2[iop]->NextPad(fInput->DetElemId()))
1408 ix = fSeg2[iop]->Ix(); iy = fSeg2[iop]->Iy();
1409 AliDebug(2,Form("\n ix, iy: %f %f %f %d %d %d", x,y,z,ix, iy, fSector));
1410 if (fDigitMap[iop]->TestHit(ix,iy)==kUnused){
1413 AliDebug(2,Form("\n Opposite %d %d %d", iop, ix, iy));
1416 } // Loop over pad neighbours
1417 // This had to go outside the loop since recursive calls inside the iterator are not possible
1420 for (jopp=0; jopp<nOpp; jopp++) {
1421 if (fDigitMap[iop]->TestHit(iXopp[jopp],iYopp[jopp]) == kUnused)
1422 FindCluster(iXopp[jopp], iYopp[jopp], iop, c);
1427 //_____________________________________________________________________________
1429 void AliMUONClusterFinderVS::FindRawClusters()
1432 // MUON cluster finder from digits -- finds neighbours on both cathodes and
1433 // fills the tree with raw clusters
1437 // Return if no input datad available
1438 if (!fInput->NDigits(0) && !fInput->NDigits(1)) return;
1440 AliMpPlaneType plane1, plane2;
1442 fSeg2[0] = fInput->Segmentation2(0);
1443 fSeg2[1] = fInput->Segmentation2(1);
1445 if ( fSeg2[0]->GetDirection(fInput->DetElemId()) == kDirY ||
1446 fSeg2[0]->GetDirection(fInput->DetElemId()) == kDirUndefined ) {
1447 plane1 = kBendingPlane;
1448 plane2 = kNonBendingPlane;
1450 plane2 = kBendingPlane;
1451 plane1 = kNonBendingPlane;
1454 fDigitMap[0] = new AliMUONDigitMapA1(fInput->DetElemId(), plane1);
1455 fDigitMap[1] = new AliMUONDigitMapA1(fInput->DetElemId(), plane2);
1463 fDigitMap[0]->FillHits(fInput->Digits(0));
1464 fDigitMap[1]->FillHits(fInput->Digits(1));
1466 // Outer Loop over Cathodes
1467 for (cath = 0; cath < 2; cath++) {
1469 for (ndig=0; ndig<fInput->NDigits(cath); ndig++) {
1470 dig = fInput->Digit(cath, ndig);
1471 Int_t padx = dig->PadX();
1472 Int_t pady = dig->PadY();
1473 if (fDigitMap[cath]->TestHit(padx,pady)==kUsed ||fDigitMap[0]->TestHit(padx,pady)==kEmpty) {
1477 AliDebug(1,Form("\n CATHODE %d CLUSTER %d\n",cath,ncls));
1478 AliMUONRawCluster clus;
1479 clus.SetMultiplicity(0, 0);
1480 clus.SetMultiplicity(1, 0);
1481 clus.SetPeakSignal(cath,dig->Signal());
1482 clus.SetTrack(0, dig->Hit());
1483 clus.SetTrack(1, dig->Track(0));
1484 clus.SetTrack(2, dig->Track(1));
1486 AliDebug(1,Form("idDE %d Padx %d Pady %d", fInput->DetElemId(), padx, pady));
1488 // tag the beginning of cluster list in a raw cluster
1489 clus.SetNcluster(0,-1);
1491 fSeg2[cath]->GetPadC(fInput->DetElemId(), padx, pady, xcu, ycu, fZPlane);
1492 fSector= fSeg2[cath]->Sector(fInput->DetElemId(), padx, pady)/100;
1497 FindCluster(padx,pady,cath,clus);
1498 //^^^^^^^^^^^^^^^^^^^^^^^^
1499 // center of gravity
1500 if (clus.GetX(0)!=0.) clus.SetX(0, clus.GetX(0)/clus.GetCharge(0)); // clus.fX[0] /= clus.fQ[0];
1503 clus.SetX(0,fSeg2[0]->GetAnod(fInput->DetElemId(), clus.GetX(0)));
1504 if (clus.GetY(0)!=0.) clus.SetY(0, clus.GetY(0)/clus.GetCharge(0)); // clus.fY[0] /= clus.fQ[0];
1506 if(clus.GetCharge(1)!=0.) clus.SetX(1, clus.GetX(1)/clus.GetCharge(1)); // clus.fX[1] /= clus.fQ[1];
1509 clus.SetX(1, fSeg2[0]->GetAnod(fInput->DetElemId(),clus.GetX(1)));
1510 if(clus.GetCharge(1)!=0.) clus.SetY(1, clus.GetY(1)/clus.GetCharge(1));// clus.fY[1] /= clus.fQ[1];
1512 clus.SetZ(0, fZPlane);
1513 clus.SetZ(1, fZPlane);
1515 AliDebug(1,Form("\n Cathode 1 multiplicite %d X(CG) %f Y(CG) %f\n",
1516 clus.GetMultiplicity(0),clus.GetX(0),clus.GetY(0)));
1517 AliDebug(1,Form(" Cathode 2 multiplicite %d X(CG) %f Y(CG) %f\n",
1518 clus.GetMultiplicity(1),clus.GetX(1),clus.GetY(1)));
1519 // Analyse cluster and decluster if necessary
1522 clus.SetNcluster(1,fNRawClusters);
1523 clus.SetClusterType(clus.PhysicsContribution());
1530 // reset Cluster object
1531 { // begin local scope
1532 for (int k=0;k<clus.GetMultiplicity(0);k++) clus.SetIndex(k, 0, 0);
1533 } // end local scope
1535 { // begin local scope
1536 for (int k=0;k<clus.GetMultiplicity(1);k++) clus.SetIndex(k, 1, 0);
1537 } // end local scope
1539 clus.SetMultiplicity(0,0);
1540 clus.SetMultiplicity(1,0);
1544 } // end loop cathodes
1545 delete fDigitMap[0];
1546 delete fDigitMap[1];
1549 Float_t AliMUONClusterFinderVS::SingleMathiesonFit(AliMUONRawCluster *c, Int_t cath)
1551 // Performs a single Mathieson fit on one cathode
1553 Double_t arglist[20];
1555 AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
1557 clusterInput.Fitter()->SetFCN(fcnS1);
1558 clusterInput.Fitter()->mninit(2,10,7);
1559 clusterInput.Fitter()->SetPrintLevel(-1 + AliLog::GetGlobalDebugLevel());
1561 clusterInput.Fitter()->mnexcm("SET NOW", arglist, 0, ierflag);
1562 // Set starting values
1563 static Double_t vstart[2];
1564 vstart[0]=c->GetX(1);
1565 vstart[1]=c->GetY(0);
1568 // lower and upper limits
1569 static Double_t lower[2], upper[2];
1571 fSeg2[cath]->GetPadI(fInput->DetElemId(), c->GetX(cath), c->GetY(cath), fZPlane, ix, iy);
1572 isec=fSeg2[cath]->Sector(fInput->DetElemId(), ix, iy);
1574 lower[0]=vstart[0]-fSeg2[cath]->Dpx(fInput->DetElemId(), isec)/2;
1575 lower[1]=vstart[1]-fSeg2[cath]->Dpy(fInput->DetElemId(), isec)/2;
1577 upper[0]=lower[0]+fSeg2[cath]->Dpx(fInput->DetElemId(), isec);
1578 upper[1]=lower[1]+fSeg2[cath]->Dpy(fInput->DetElemId(), isec);
1582 static Double_t step[2]={0.0005, 0.0005};
1584 clusterInput.Fitter()->mnparm(0,"x1",vstart[0],step[0],lower[0],upper[0],ierflag);
1585 clusterInput.Fitter()->mnparm(1,"y1",vstart[1],step[1],lower[1],upper[1],ierflag);
1586 // ready for minimisation
1590 clusterInput.Fitter()->mnexcm("SET NOGR", arglist, 0, ierflag);
1591 clusterInput.Fitter()->mnexcm("MIGRAD", arglist, 0, ierflag);
1592 // clusterInput.Fitter()->mnexcm("EXIT" , arglist, 0, ierflag);
1593 Double_t fmin, fedm, errdef;
1594 Int_t npari, nparx, istat;
1596 clusterInput.Fitter()->mnstat(fmin, fedm, errdef, npari, nparx, istat);
1600 // Get fitted parameters
1601 Double_t xrec, yrec;
1603 Double_t epxz, b1, b2;
1605 clusterInput.Fitter()->mnpout(0, chname, xrec, epxz, b1, b2, ierflg);
1606 clusterInput.Fitter()->mnpout(1, chname, yrec, epxz, b1, b2, ierflg);
1612 Float_t AliMUONClusterFinderVS::CombiSingleMathiesonFit(AliMUONRawCluster * /*c*/)
1614 // Perform combined Mathieson fit on both cathode planes
1616 Double_t arglist[20];
1618 AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
1619 clusterInput.Fitter()->SetFCN(fcnCombiS1);
1620 clusterInput.Fitter()->mninit(2,10,7);
1621 clusterInput.Fitter()->SetPrintLevel(-1 + AliLog::GetGlobalDebugLevel());
1623 clusterInput.Fitter()->mnexcm("SET NOW", arglist, 0, ierflag);
1624 static Double_t vstart[2];
1625 vstart[0]=fXInit[0];
1626 vstart[1]=fYInit[0];
1629 // lower and upper limits
1630 static Float_t lower[2], upper[2];
1634 fSeg2[0]->GetPadI(fInput->DetElemId(), fXInit[0], fYInit[0], fZPlane, ix, iy);
1635 isec=fSeg2[0]->Sector(fInput->DetElemId(), ix, iy);
1636 dpy=fSeg2[0]->Dpy(fInput->DetElemId(), isec);
1637 fSeg2[1]->GetPadI(fInput->DetElemId(), fXInit[0], fYInit[0], fZPlane, ix, iy);
1638 isec=fSeg2[1]->Sector(fInput->DetElemId(), ix, iy);
1639 dpx=fSeg2[1]->Dpx(fInput->DetElemId(), isec);
1642 Float_t xdum, ydum, zdum;
1644 // Find save upper and lower limits
1647 for (fSeg2[1]->FirstPad(fInput->DetElemId(),fXInit[0], fYInit[0], fZPlane, dpx, 0.);
1648 fSeg2[1]->MorePads(fInput->DetElemId());
1649 fSeg2[1]->NextPad(fInput->DetElemId()))
1651 ix=fSeg2[1]->Ix(); iy=fSeg2[1]->Iy();
1652 fSeg2[1]->GetPadC(fInput->DetElemId(), ix,iy, upper[0], ydum, zdum);
1653 if (icount ==0) lower[0]=upper[0];
1657 if (lower[0]>upper[0]) {xdum=lower[0]; lower[0]=upper[0]; upper[0]=xdum;}
1660 AliDebug(1,Form("\n single y %f %f", fXInit[0], fYInit[0]));
1662 for (fSeg2[0]->FirstPad(fInput->DetElemId(), fXInit[0], fYInit[0], fZPlane, 0., dpy);
1663 fSeg2[0]->MorePads(fInput->DetElemId());
1664 fSeg2[0]->NextPad(fInput->DetElemId()))
1666 ix=fSeg2[0]->Ix(); iy=fSeg2[0]->Iy();
1667 fSeg2[0]->GetPadC(fInput->DetElemId(), ix,iy,xdum,upper[1],zdum);
1668 if (icount ==0) lower[1]=upper[1];
1670 AliDebug(1,Form("\n upper lower %d %f %f", icount, upper[1], lower[1]));
1673 if (lower[1]>upper[1]) {xdum=lower[1]; lower[1]=upper[1]; upper[1]=xdum;}
1676 static Double_t step[2]={0.00001, 0.0001};
1678 clusterInput.Fitter()->mnparm(0,"x1",vstart[0],step[0],lower[0],upper[0],ierflag);
1679 clusterInput.Fitter()->mnparm(1,"y1",vstart[1],step[1],lower[1],upper[1],ierflag);
1680 // ready for minimisation
1684 clusterInput.Fitter()->mnexcm("SET NOGR", arglist, 0, ierflag);
1685 clusterInput.Fitter()->mnexcm("MIGRAD", arglist, 0, ierflag);
1686 // clusterInput.Fitter()->mnexcm("EXIT" , arglist, 0, ierflag);
1687 Double_t fmin, fedm, errdef;
1688 Int_t npari, nparx, istat;
1690 clusterInput.Fitter()->mnstat(fmin, fedm, errdef, npari, nparx, istat);
1694 // Get fitted parameters
1695 Double_t xrec, yrec;
1697 Double_t epxz, b1, b2;
1699 clusterInput.Fitter()->mnpout(0, chname, xrec, epxz, b1, b2, ierflg);
1700 clusterInput.Fitter()->mnpout(1, chname, yrec, epxz, b1, b2, ierflg);
1706 Bool_t AliMUONClusterFinderVS::DoubleMathiesonFit(AliMUONRawCluster * /*c*/, Int_t cath)
1708 // Performs a double Mathieson fit on one cathode
1712 // Initialise global variables for fit
1713 Double_t arglist[20];
1715 AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
1716 clusterInput.Fitter()->SetFCN(fcnS2);
1717 clusterInput.Fitter()->mninit(5,10,7);
1718 clusterInput.Fitter()->SetPrintLevel(-1 + AliLog::GetGlobalDebugLevel());
1720 clusterInput.Fitter()->mnexcm("SET NOW", arglist, 0, ierflag);
1721 // Set starting values
1722 static Double_t vstart[5];
1723 vstart[0]=fX[fIndLocal[0][cath]][cath];
1724 vstart[1]=fY[fIndLocal[0][cath]][cath];
1725 vstart[2]=fX[fIndLocal[1][cath]][cath];
1726 vstart[3]=fY[fIndLocal[1][cath]][cath];
1727 vstart[4]=Float_t(fQ[fIndLocal[0][cath]][cath])/
1728 Float_t(fQ[fIndLocal[0][cath]][cath]+fQ[fIndLocal[1][cath]][cath]);
1729 // lower and upper limits
1730 static Float_t lower[5], upper[5];
1733 isec=fSeg2[cath]->Sector(fInput->DetElemId(),fIx[fIndLocal[0][cath]][cath],
1734 fIy[fIndLocal[0][cath]][cath]);
1735 lower[0]=vstart[0]-fSeg2[cath]->Dpx(fInput->DetElemId(),isec);
1736 lower[1]=vstart[1]-fSeg2[cath]->Dpy(fInput->DetElemId(),isec);
1738 upper[0]=lower[0]+2.*fSeg2[cath]->Dpx(fInput->DetElemId(),isec);
1739 upper[1]=lower[1]+2.*fSeg2[cath]->Dpy(fInput->DetElemId(),isec);
1741 isec=fSeg2[cath]->Sector(fInput->DetElemId(),fIx[fIndLocal[1][cath]][cath],
1742 fIy[fIndLocal[1][cath]][cath]);
1743 lower[2]=vstart[2]-fSeg2[cath]->Dpx(fInput->DetElemId(),isec)/2;
1744 lower[3]=vstart[3]-fSeg2[cath]->Dpy(fInput->DetElemId(),isec)/2;
1746 upper[2]=lower[2]+fSeg2[cath]->Dpx(fInput->DetElemId(),isec);
1747 upper[1]=lower[1]+2.*fSeg2[cath]->Dpy(fInput->DetElemId(),isec);
1754 static Double_t step[5]={0.0005, 0.0005, 0.0005, 0.0005, 0.0001};
1756 clusterInput.Fitter()->mnparm(0,"x1",vstart[0],step[0],lower[0],upper[0],ierflag);
1757 clusterInput.Fitter()->mnparm(1,"y1",vstart[1],step[1],lower[1],upper[1],ierflag);
1758 clusterInput.Fitter()->mnparm(2,"x2",vstart[2],step[2],lower[2],upper[2],ierflag);
1759 clusterInput.Fitter()->mnparm(3,"y2",vstart[3],step[3],lower[3],upper[3],ierflag);
1760 clusterInput.Fitter()->mnparm(4,"a0",vstart[4],step[4],lower[4],upper[4],ierflag);
1761 // ready for minimisation
1765 clusterInput.Fitter()->mnexcm("SET NOGR", arglist, 0, ierflag);
1766 clusterInput.Fitter()->mnexcm("MIGRAD", arglist, 0, ierflag);
1767 // clusterInput.Fitter()->mnexcm("EXIT" , arglist, 0, ierflag);
1768 // Get fitted parameters
1769 Double_t xrec[2], yrec[2], qfrac;
1771 Double_t epxz, b1, b2;
1773 clusterInput.Fitter()->mnpout(0, chname, xrec[0], epxz, b1, b2, ierflg);
1774 clusterInput.Fitter()->mnpout(1, chname, yrec[0], epxz, b1, b2, ierflg);
1775 clusterInput.Fitter()->mnpout(2, chname, xrec[1], epxz, b1, b2, ierflg);
1776 clusterInput.Fitter()->mnpout(3, chname, yrec[1], epxz, b1, b2, ierflg);
1777 clusterInput.Fitter()->mnpout(4, chname, qfrac, epxz, b1, b2, ierflg);
1779 Double_t fmin, fedm, errdef;
1780 Int_t npari, nparx, istat;
1782 clusterInput.Fitter()->mnstat(fmin, fedm, errdef, npari, nparx, istat);
1787 Float_t AliMUONClusterFinderVS::CombiDoubleMathiesonFit(AliMUONRawCluster * /*c*/)
1790 // Perform combined double Mathieson fit on both cathode planes
1792 Double_t arglist[20];
1794 AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
1795 clusterInput.Fitter()->SetFCN(fcnCombiS2);
1796 clusterInput.Fitter()->mninit(6,10,7);
1797 clusterInput.Fitter()->SetPrintLevel(-1 + AliLog::GetGlobalDebugLevel());
1799 clusterInput.Fitter()->mnexcm("SET NOW", arglist, 0, ierflag);
1800 // Set starting values
1801 static Double_t vstart[6];
1802 vstart[0]=fXInit[0];
1803 vstart[1]=fYInit[0];
1804 vstart[2]=fXInit[1];
1805 vstart[3]=fYInit[1];
1806 vstart[4]=fQrInit[0];
1807 vstart[5]=fQrInit[1];
1808 // lower and upper limits
1809 static Float_t lower[6], upper[6];
1813 fSeg2[1]->GetPadI(fInput->DetElemId(),fXInit[0], fYInit[0], fZPlane, ix, iy);
1814 isec=fSeg2[1]->Sector(fInput->DetElemId(),ix, iy);
1815 dpx=fSeg2[1]->Dpx(fInput->DetElemId(), isec);
1817 fSeg2[0]->GetPadI(fInput->DetElemId(), fXInit[0], fYInit[0], fZPlane, ix, iy);
1818 isec=fSeg2[0]->Sector(fInput->DetElemId(), ix, iy);
1819 dpy=fSeg2[0]->Dpy(fInput->DetElemId(), isec);
1824 Float_t xdum, ydum, zdum;
1825 AliDebug(1,Form("\n Cluster Finder: %f %f %f %f ", fXInit[0], fXInit[1],fYInit[0], fYInit[1] ));
1827 // Find save upper and lower limits
1830 for (fSeg2[1]->FirstPad(fInput->DetElemId(),fXInit[0], fYInit[0], fZPlane, dpx, 0.);
1831 fSeg2[1]->MorePads(fInput->DetElemId());
1832 fSeg2[1]->NextPad(fInput->DetElemId()))
1834 ix=fSeg2[1]->Ix(); iy=fSeg2[1]->Iy();
1835 // if (fDigitMap[1]->TestHit(ix, iy) == kEmpty) continue;
1836 fSeg2[1]->GetPadC(fInput->DetElemId(),ix,iy,upper[0],ydum,zdum);
1837 if (icount ==0) lower[0]=upper[0];
1840 if (lower[0]>upper[0]) {xdum=lower[0]; lower[0]=upper[0]; upper[0]=xdum;}
1841 // vstart[0] = 0.5*(lower[0]+upper[0]);
1846 for (fSeg2[0]->FirstPad(fInput->DetElemId(),fXInit[0], fYInit[0], fZPlane, 0., dpy);
1847 fSeg2[0]->MorePads(fInput->DetElemId());
1848 fSeg2[0]->NextPad(fInput->DetElemId()))
1850 ix=fSeg2[0]->Ix(); iy=fSeg2[0]->Iy();
1851 // if (fDigitMap[0]->TestHit(ix, iy) == kEmpty) continue;
1852 fSeg2[0]->GetPadC(fInput->DetElemId(),ix,iy,xdum,upper[1],zdum);
1853 if (icount ==0) lower[1]=upper[1];
1857 if (lower[1]>upper[1]) {xdum=lower[1]; lower[1]=upper[1]; upper[1]=xdum;}
1858 // vstart[1] = 0.5*(lower[1]+upper[1]);
1861 fSeg2[1]->GetPadI(fInput->DetElemId(),fXInit[1], fYInit[1], fZPlane, ix, iy);
1862 isec=fSeg2[1]->Sector(fInput->DetElemId(),ix, iy);
1863 dpx=fSeg2[1]->Dpx(fInput->DetElemId(),isec);
1864 fSeg2[0]->GetPadI(fInput->DetElemId(),fXInit[1], fYInit[1], fZPlane, ix, iy);
1865 isec=fSeg2[0]->Sector(fInput->DetElemId(),ix, iy);
1866 dpy=fSeg2[0]->Dpy(fInput->DetElemId(),isec);
1869 // Find save upper and lower limits
1873 for (fSeg2[1]->FirstPad(fInput->DetElemId(),fXInit[1], fYInit[1], fZPlane, dpx, 0);
1874 fSeg2[1]->MorePads(fInput->DetElemId());
1875 fSeg2[1]->NextPad(fInput->DetElemId()))
1877 ix=fSeg2[1]->Ix(); iy=fSeg2[1]->Iy();
1878 // if (fDigitMap[1]->TestHit(ix, iy) == kEmpty) continue;
1879 fSeg2[1]->GetPadC(fInput->DetElemId(),ix,iy,upper[2],ydum,zdum);
1880 if (icount ==0) lower[2]=upper[2];
1883 if (lower[2]>upper[2]) {xdum=lower[2]; lower[2]=upper[2]; upper[2]=xdum;}
1884 // vstart[2] = 0.5*(lower[2]+upper[2]);
1888 for (fSeg2[0]->FirstPad(fInput->DetElemId(),fXInit[1], fYInit[1], fZPlane, 0, dpy);
1889 fSeg2[0]-> MorePads(fInput->DetElemId());
1890 fSeg2[0]->NextPad(fInput->DetElemId()))
1892 ix=fSeg2[0]->Ix(); iy=fSeg2[0]->Iy();
1893 // if (fDigitMap[0]->TestHit(ix, iy) != kEmpty) continue;
1895 fSeg2[0]->GetPadC(fInput->DetElemId(),ix,iy,xdum,upper[3],zdum);
1896 if (icount ==0) lower[3]=upper[3];
1900 if (lower[3]>upper[3]) {xdum=lower[3]; lower[3]=upper[3]; upper[3]=xdum;}
1908 static Double_t step[6]={0.0005, 0.0005, 0.0005, 0.0005, 0.001, 0.001};
1909 clusterInput.Fitter()->mnparm(0,"x1",vstart[0],step[0],lower[0],upper[0],ierflag);
1910 clusterInput.Fitter()->mnparm(1,"y1",vstart[1],step[1],lower[1],upper[1],ierflag);
1911 clusterInput.Fitter()->mnparm(2,"x2",vstart[2],step[2],lower[2],upper[2],ierflag);
1912 clusterInput.Fitter()->mnparm(3,"y2",vstart[3],step[3],lower[3],upper[3],ierflag);
1913 clusterInput.Fitter()->mnparm(4,"a0",vstart[4],step[4],lower[4],upper[4],ierflag);
1914 clusterInput.Fitter()->mnparm(5,"a1",vstart[5],step[5],lower[5],upper[5],ierflag);
1915 // ready for minimisation
1919 clusterInput.Fitter()->mnexcm("SET NOGR", arglist, 0, ierflag);
1920 clusterInput.Fitter()->mnexcm("MIGRAD", arglist, 0, ierflag);
1921 // clusterInput.Fitter()->mnexcm("EXIT" , arglist, 0, ierflag);
1922 // Get fitted parameters
1924 Double_t epxz, b1, b2;
1926 clusterInput.Fitter()->mnpout(0, chname, fXFit[0], epxz, b1, b2, ierflg);
1927 clusterInput.Fitter()->mnpout(1, chname, fYFit[0], epxz, b1, b2, ierflg);
1928 clusterInput.Fitter()->mnpout(2, chname, fXFit[1], epxz, b1, b2, ierflg);
1929 clusterInput.Fitter()->mnpout(3, chname, fYFit[1], epxz, b1, b2, ierflg);
1930 clusterInput.Fitter()->mnpout(4, chname, fQrFit[0], epxz, b1, b2, ierflg);
1931 clusterInput.Fitter()->mnpout(5, chname, fQrFit[1], epxz, b1, b2, ierflg);
1933 Double_t fmin, fedm, errdef;
1934 Int_t npari, nparx, istat;
1936 clusterInput.Fitter()->mnstat(fmin, fedm, errdef, npari, nparx, istat);
1944 void AliMUONClusterFinderVS::Split(AliMUONRawCluster* c)
1947 // One cluster for each maximum
1950 AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
1951 for (j=0; j<2; j++) {
1952 AliMUONRawCluster cnew;
1953 cnew.SetGhost(c->GetGhost());
1954 for (cath=0; cath<2; cath++) {
1955 cnew.SetChi2(cath,fChi2[0]);
1956 // ?? why not cnew.fChi2[cath]=fChi2[cath];
1959 cnew.SetNcluster(0,-1);
1960 cnew.SetNcluster(1,fNRawClusters);
1962 cnew.SetNcluster(0,fNPeaks);
1963 cnew.SetNcluster(1,0);
1965 cnew.SetMultiplicity(cath,0);
1966 cnew.SetX(cath, Float_t(fXFit[j]));
1967 cnew.SetY(cath, Float_t(fYFit[j]));
1968 cnew.SetZ(cath, fZPlane);
1970 cnew.SetCharge(cath, Int_t(clusterInput.TotalCharge(cath)*fQrFit[cath]));
1972 cnew.SetCharge(cath, Int_t(clusterInput.TotalCharge(cath)*(1-fQrFit[cath])));
1974 fSeg2[cath]->SetHit(fInput->DetElemId(), fXFit[j],fYFit[j],fZPlane);
1976 for (i=0; i<fMul[cath]; i++) {
1978 cnew.SetIndex(cnew.GetMultiplicity(cath), cath, c->GetIndex(i,cath));
1980 fSeg2[cath]->SetPad(fInput->DetElemId(),fIx[i][cath], fIy[i][cath]);
1981 q1 = fInput->Mathieson()->IntXY(fInput->DetElemId(),fSeg2[cath]);
1983 cnew.SetContrib(i, cath, q1*Float_t(cnew.GetCharge(cath))/Float_t(fQ[i][cath]));
1984 cnew.SetMultiplicity(cath, cnew.GetMultiplicity(cath)+1 );
1986 FillCluster(&cnew,0,cath);
1988 cnew.SetClusterType(cnew.PhysicsContribution());
1989 if (cnew.GetCharge(0)>0 && cnew.GetCharge(1)>0) AddRawCluster(cnew);
1993 void AliMUONClusterFinderVS::AddRawCluster(AliMUONRawCluster& c)
1996 // Add a raw cluster copy to the list
1998 // AliMUON *pMUON=(AliMUON*)gAlice->GetModule("MUON");
1999 // pMUON->GetMUONData()->AddRawCluster(fInput->Chamber(),c);
2002 // Setting detection element in raw cluster for alignment
2004 c.SetDetElemId(fInput->DetElemId());
2006 TClonesArray &lrawcl = *fRawClusters;
2007 new(lrawcl[fNRawClusters++]) AliMUONRawCluster(c);
2008 AliDebug(1,Form("\nfNRawClusters %d\n",fNRawClusters));
2011 AliMUONClusterFinderVS& AliMUONClusterFinderVS
2012 ::operator = (const AliMUONClusterFinderVS& rhs)
2014 // Protected assignement operator
2016 if (this == &rhs) return *this;
2018 AliFatal("Not implemented.");
2024 // Minimisation functions
2026 void fcnS1(Int_t & /*npar*/, Double_t * /*gin*/, Double_t &f, Double_t *par, Int_t /*iflag*/)
2028 AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
2035 for (i=0; i<clusterInput.Nmul(0); i++) {
2036 Float_t q0=clusterInput.Charge(i,0);
2037 Float_t q1=clusterInput.DiscrChargeS1(i,par);
2046 void fcnCombiS1(Int_t & /*npar*/, Double_t * /*gin*/, Double_t &f, Double_t *par, Int_t /*iflag*/)
2048 AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
2055 for (cath=0; cath<2; cath++) {
2056 for (i=0; i<clusterInput.Nmul(cath); i++) {
2057 Float_t q0=clusterInput.Charge(i,cath);
2058 Float_t q1=clusterInput.DiscrChargeCombiS1(i,par,cath);
2069 void fcnS2(Int_t & /*npar*/, Double_t * /*gin*/, Double_t &f, Double_t *par, Int_t /*iflag*/)
2071 AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
2078 for (i=0; i<clusterInput.Nmul(0); i++) {
2080 Float_t q0=clusterInput.Charge(i,0);
2081 Float_t q1=clusterInput.DiscrChargeS2(i,par);
2091 void fcnCombiS2(Int_t & /*npar*/, Double_t * /*gin*/, Double_t &f, Double_t *par, Int_t /*iflag*/)
2093 AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
2099 for (cath=0; cath<2; cath++) {
2100 for (i=0; i<clusterInput.Nmul(cath); i++) {
2101 Float_t q0=clusterInput.Charge(i,cath);
2102 Float_t q1=clusterInput.DiscrChargeCombiS2(i,par,cath);