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 "AliMUONDigit.h"
25 #include "AliMUONRawCluster.h"
26 #include "AliMUONGeometrySegmentation.h"
27 #include "AliMUONMathieson.h"
28 #include "AliMUONClusterInput.h"
29 #include "AliMUONHitMapA1.h"
32 //_____________________________________________________________________
33 // This function is minimized in the double-Mathieson fit
34 void fcnS2(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
35 void fcnS1(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
36 void fcnCombiS1(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
37 void fcnCombiS2(Int_t &npar, Double_t *gin, Double_t &f, Double_t *par, Int_t iflag);
39 ClassImp(AliMUONClusterFinderVS)
41 AliMUONClusterFinderVS::AliMUONClusterFinderVS()
44 // Default constructor
45 fInput=AliMUONClusterInput::Instance();
46 // cout << " TYPE" << fSegmentationType << endl;
49 fTrack[0]=fTrack[1]=-1;
50 fDebugLevel = 0; // make silent default
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 (fHitMap[cath]->TestHit(x[j], y[j])==kEmpty) continue;
990 digt=(AliMUONDigit*) fHitMap[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 (fHitMap[cath]->TestHit(ix, iy)!=kEmpty) {
1066 digt=(AliMUONDigit*) fHitMap[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 (fHitMap[cath]->TestHit(ix, iy)!=kEmpty) {
1131 digt=(AliMUONDigit*) fHitMap[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 = fHitMap[cath]->GetHitIndex(i,j);
1303 AliMUONDigit* dig = (AliMUONDigit*) fHitMap[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 fHitMap[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 (fHitMap[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 (fHitMap[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 (fHitMap[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 if(fInput->GetSegmentationType() == 1)
1441 AliFatal("Old Segmentation no more supported.");
1443 fSeg2[0] = fInput->Segmentation2(0);
1444 fSeg2[1] = fInput->Segmentation2(1);
1446 fHitMap[0] = new AliMUONHitMapA1(fInput->DetElemId(), fSeg2[0], fInput->Digits(0));
1447 fHitMap[1] = new AliMUONHitMapA1(fInput->DetElemId(), fSeg2[1], fInput->Digits(1));
1454 fHitMap[0]->FillHits();
1455 fHitMap[1]->FillHits();
1457 // Outer Loop over Cathodes
1458 for (cath=0; cath<2; cath++) {
1460 for (ndig=0; ndig<fInput->NDigits(cath); ndig++) {
1461 dig = fInput->Digit(cath, ndig);
1462 Int_t padx = dig->PadX();
1463 Int_t pady = dig->PadY();
1464 if (fHitMap[cath]->TestHit(padx,pady)==kUsed ||fHitMap[0]->TestHit(padx,pady)==kEmpty) {
1468 AliDebug(1,Form("\n CATHODE %d CLUSTER %d\n",cath,ncls));
1469 AliMUONRawCluster clus;
1470 clus.SetMultiplicity(0, 0);
1471 clus.SetMultiplicity(1, 0);
1472 clus.SetPeakSignal(cath,dig->Signal());
1473 clus.SetTrack(0, dig->Hit());
1474 clus.SetTrack(1, dig->Track(0));
1475 clus.SetTrack(2, dig->Track(1));
1477 AliDebug(1,Form("idDE %d Padx %d Pady %d", fInput->DetElemId(), padx, pady));
1479 // tag the beginning of cluster list in a raw cluster
1480 clus.SetNcluster(0,-1);
1482 fSeg2[cath]->GetPadC(fInput->DetElemId(), padx, pady, xcu, ycu, fZPlane);
1483 fSector= fSeg2[cath]->Sector(fInput->DetElemId(), padx, pady)/100;
1488 FindCluster(padx,pady,cath,clus);
1489 //^^^^^^^^^^^^^^^^^^^^^^^^
1490 // center of gravity
1491 if (clus.GetX(0)!=0.) clus.SetX(0, clus.GetX(0)/clus.GetCharge(0)); // clus.fX[0] /= clus.fQ[0];
1494 clus.SetX(0,fSeg2[0]->GetAnod(fInput->DetElemId(), clus.GetX(0)));
1495 if (clus.GetY(0)!=0.) clus.SetY(0, clus.GetY(0)/clus.GetCharge(0)); // clus.fY[0] /= clus.fQ[0];
1497 if(clus.GetCharge(1)!=0.) clus.SetX(1, clus.GetX(1)/clus.GetCharge(1)); // clus.fX[1] /= clus.fQ[1];
1500 clus.SetX(1, fSeg2[0]->GetAnod(fInput->DetElemId(),clus.GetX(1)));
1501 if(clus.GetCharge(1)!=0.) clus.SetY(1, clus.GetY(1)/clus.GetCharge(1));// clus.fY[1] /= clus.fQ[1];
1503 clus.SetZ(0, fZPlane);
1504 clus.SetZ(1, fZPlane);
1506 AliDebug(1,Form("\n Cathode 1 multiplicite %d X(CG) %f Y(CG) %f\n",
1507 clus.GetMultiplicity(0),clus.GetX(0),clus.GetY(0)));
1508 AliDebug(1,Form(" Cathode 2 multiplicite %d X(CG) %f Y(CG) %f\n",
1509 clus.GetMultiplicity(1),clus.GetX(1),clus.GetY(1)));
1510 // Analyse cluster and decluster if necessary
1513 clus.SetNcluster(1,fNRawClusters);
1514 clus.SetClusterType(clus.PhysicsContribution());
1521 // reset Cluster object
1522 { // begin local scope
1523 for (int k=0;k<clus.GetMultiplicity(0);k++) clus.SetIndex(k, 0, 0);
1524 } // end local scope
1526 { // begin local scope
1527 for (int k=0;k<clus.GetMultiplicity(1);k++) clus.SetIndex(k, 1, 0);
1528 } // end local scope
1530 clus.SetMultiplicity(0,0);
1531 clus.SetMultiplicity(1,0);
1535 } // end loop cathodes
1540 Float_t AliMUONClusterFinderVS::SingleMathiesonFit(AliMUONRawCluster *c, Int_t cath)
1542 // Performs a single Mathieson fit on one cathode
1544 Double_t arglist[20];
1546 AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
1548 clusterInput.Fitter()->SetFCN(fcnS1);
1549 clusterInput.Fitter()->mninit(2,10,7);
1550 clusterInput.Fitter()->SetPrintLevel(-1+fDebugLevel);
1552 clusterInput.Fitter()->mnexcm("SET NOW", arglist, 0, ierflag);
1553 // Set starting values
1554 static Double_t vstart[2];
1555 vstart[0]=c->GetX(1);
1556 vstart[1]=c->GetY(0);
1559 // lower and upper limits
1560 static Double_t lower[2], upper[2];
1562 fSeg2[cath]->GetPadI(fInput->DetElemId(), c->GetX(cath), c->GetY(cath), fZPlane, ix, iy);
1563 isec=fSeg2[cath]->Sector(fInput->DetElemId(), ix, iy);
1565 lower[0]=vstart[0]-fSeg2[cath]->Dpx(fInput->DetElemId(), isec)/2;
1566 lower[1]=vstart[1]-fSeg2[cath]->Dpy(fInput->DetElemId(), isec)/2;
1568 upper[0]=lower[0]+fSeg2[cath]->Dpx(fInput->DetElemId(), isec);
1569 upper[1]=lower[1]+fSeg2[cath]->Dpy(fInput->DetElemId(), isec);
1573 static Double_t step[2]={0.0005, 0.0005};
1575 clusterInput.Fitter()->mnparm(0,"x1",vstart[0],step[0],lower[0],upper[0],ierflag);
1576 clusterInput.Fitter()->mnparm(1,"y1",vstart[1],step[1],lower[1],upper[1],ierflag);
1577 // ready for minimisation
1581 clusterInput.Fitter()->mnexcm("SET NOGR", arglist, 0, ierflag);
1582 clusterInput.Fitter()->mnexcm("MIGRAD", arglist, 0, ierflag);
1583 // clusterInput.Fitter()->mnexcm("EXIT" , arglist, 0, ierflag);
1584 Double_t fmin, fedm, errdef;
1585 Int_t npari, nparx, istat;
1587 clusterInput.Fitter()->mnstat(fmin, fedm, errdef, npari, nparx, istat);
1591 // Get fitted parameters
1592 Double_t xrec, yrec;
1594 Double_t epxz, b1, b2;
1596 clusterInput.Fitter()->mnpout(0, chname, xrec, epxz, b1, b2, ierflg);
1597 clusterInput.Fitter()->mnpout(1, chname, yrec, epxz, b1, b2, ierflg);
1603 Float_t AliMUONClusterFinderVS::CombiSingleMathiesonFit(AliMUONRawCluster * /*c*/)
1605 // Perform combined Mathieson fit on both cathode planes
1607 Double_t arglist[20];
1609 AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
1610 clusterInput.Fitter()->SetFCN(fcnCombiS1);
1611 clusterInput.Fitter()->mninit(2,10,7);
1612 clusterInput.Fitter()->SetPrintLevel(-1+fDebugLevel);
1614 clusterInput.Fitter()->mnexcm("SET NOW", arglist, 0, ierflag);
1615 static Double_t vstart[2];
1616 vstart[0]=fXInit[0];
1617 vstart[1]=fYInit[0];
1620 // lower and upper limits
1621 static Float_t lower[2], upper[2];
1625 fSeg2[0]->GetPadI(fInput->DetElemId(), fXInit[0], fYInit[0], fZPlane, ix, iy);
1626 isec=fSeg2[0]->Sector(fInput->DetElemId(), ix, iy);
1627 dpy=fSeg2[0]->Dpy(fInput->DetElemId(), isec);
1628 fSeg2[1]->GetPadI(fInput->DetElemId(), fXInit[0], fYInit[0], fZPlane, ix, iy);
1629 isec=fSeg2[1]->Sector(fInput->DetElemId(), ix, iy);
1630 dpx=fSeg2[1]->Dpx(fInput->DetElemId(), isec);
1633 Float_t xdum, ydum, zdum;
1635 // Find save upper and lower limits
1638 for (fSeg2[1]->FirstPad(fInput->DetElemId(),fXInit[0], fYInit[0], fZPlane, dpx, 0.);
1639 fSeg2[1]->MorePads(fInput->DetElemId());
1640 fSeg2[1]->NextPad(fInput->DetElemId()))
1642 ix=fSeg2[1]->Ix(); iy=fSeg2[1]->Iy();
1643 fSeg2[1]->GetPadC(fInput->DetElemId(), ix,iy, upper[0], ydum, zdum);
1644 if (icount ==0) lower[0]=upper[0];
1648 if (lower[0]>upper[0]) {xdum=lower[0]; lower[0]=upper[0]; upper[0]=xdum;}
1651 AliDebug(1,Form("\n single y %f %f", fXInit[0], fYInit[0]));
1653 for (fSeg2[0]->FirstPad(fInput->DetElemId(), fXInit[0], fYInit[0], fZPlane, 0., dpy);
1654 fSeg2[0]->MorePads(fInput->DetElemId());
1655 fSeg2[0]->NextPad(fInput->DetElemId()))
1657 ix=fSeg2[0]->Ix(); iy=fSeg2[0]->Iy();
1658 fSeg2[0]->GetPadC(fInput->DetElemId(), ix,iy,xdum,upper[1],zdum);
1659 if (icount ==0) lower[1]=upper[1];
1661 AliDebug(1,Form("\n upper lower %d %f %f", icount, upper[1], lower[1]));
1664 if (lower[1]>upper[1]) {xdum=lower[1]; lower[1]=upper[1]; upper[1]=xdum;}
1667 static Double_t step[2]={0.00001, 0.0001};
1669 clusterInput.Fitter()->mnparm(0,"x1",vstart[0],step[0],lower[0],upper[0],ierflag);
1670 clusterInput.Fitter()->mnparm(1,"y1",vstart[1],step[1],lower[1],upper[1],ierflag);
1671 // ready for minimisation
1675 clusterInput.Fitter()->mnexcm("SET NOGR", arglist, 0, ierflag);
1676 clusterInput.Fitter()->mnexcm("MIGRAD", arglist, 0, ierflag);
1677 // clusterInput.Fitter()->mnexcm("EXIT" , arglist, 0, ierflag);
1678 Double_t fmin, fedm, errdef;
1679 Int_t npari, nparx, istat;
1681 clusterInput.Fitter()->mnstat(fmin, fedm, errdef, npari, nparx, istat);
1685 // Get fitted parameters
1686 Double_t xrec, yrec;
1688 Double_t epxz, b1, b2;
1690 clusterInput.Fitter()->mnpout(0, chname, xrec, epxz, b1, b2, ierflg);
1691 clusterInput.Fitter()->mnpout(1, chname, yrec, epxz, b1, b2, ierflg);
1697 Bool_t AliMUONClusterFinderVS::DoubleMathiesonFit(AliMUONRawCluster * /*c*/, Int_t cath)
1699 // Performs a double Mathieson fit on one cathode
1703 // Initialise global variables for fit
1704 Double_t arglist[20];
1706 AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
1707 clusterInput.Fitter()->SetFCN(fcnS2);
1708 clusterInput.Fitter()->mninit(5,10,7);
1709 clusterInput.Fitter()->SetPrintLevel(-1+fDebugLevel);
1711 clusterInput.Fitter()->mnexcm("SET NOW", arglist, 0, ierflag);
1712 // Set starting values
1713 static Double_t vstart[5];
1714 vstart[0]=fX[fIndLocal[0][cath]][cath];
1715 vstart[1]=fY[fIndLocal[0][cath]][cath];
1716 vstart[2]=fX[fIndLocal[1][cath]][cath];
1717 vstart[3]=fY[fIndLocal[1][cath]][cath];
1718 vstart[4]=Float_t(fQ[fIndLocal[0][cath]][cath])/
1719 Float_t(fQ[fIndLocal[0][cath]][cath]+fQ[fIndLocal[1][cath]][cath]);
1720 // lower and upper limits
1721 static Float_t lower[5], upper[5];
1724 isec=fSeg2[cath]->Sector(fInput->DetElemId(),fIx[fIndLocal[0][cath]][cath],
1725 fIy[fIndLocal[0][cath]][cath]);
1726 lower[0]=vstart[0]-fSeg2[cath]->Dpx(fInput->DetElemId(),isec);
1727 lower[1]=vstart[1]-fSeg2[cath]->Dpy(fInput->DetElemId(),isec);
1729 upper[0]=lower[0]+2.*fSeg2[cath]->Dpx(fInput->DetElemId(),isec);
1730 upper[1]=lower[1]+2.*fSeg2[cath]->Dpy(fInput->DetElemId(),isec);
1732 isec=fSeg2[cath]->Sector(fInput->DetElemId(),fIx[fIndLocal[1][cath]][cath],
1733 fIy[fIndLocal[1][cath]][cath]);
1734 lower[2]=vstart[2]-fSeg2[cath]->Dpx(fInput->DetElemId(),isec)/2;
1735 lower[3]=vstart[3]-fSeg2[cath]->Dpy(fInput->DetElemId(),isec)/2;
1737 upper[2]=lower[2]+fSeg2[cath]->Dpx(fInput->DetElemId(),isec);
1738 upper[1]=lower[1]+2.*fSeg2[cath]->Dpy(fInput->DetElemId(),isec);
1745 static Double_t step[5]={0.0005, 0.0005, 0.0005, 0.0005, 0.0001};
1747 clusterInput.Fitter()->mnparm(0,"x1",vstart[0],step[0],lower[0],upper[0],ierflag);
1748 clusterInput.Fitter()->mnparm(1,"y1",vstart[1],step[1],lower[1],upper[1],ierflag);
1749 clusterInput.Fitter()->mnparm(2,"x2",vstart[2],step[2],lower[2],upper[2],ierflag);
1750 clusterInput.Fitter()->mnparm(3,"y2",vstart[3],step[3],lower[3],upper[3],ierflag);
1751 clusterInput.Fitter()->mnparm(4,"a0",vstart[4],step[4],lower[4],upper[4],ierflag);
1752 // ready for minimisation
1756 clusterInput.Fitter()->mnexcm("SET NOGR", arglist, 0, ierflag);
1757 clusterInput.Fitter()->mnexcm("MIGRAD", arglist, 0, ierflag);
1758 // clusterInput.Fitter()->mnexcm("EXIT" , arglist, 0, ierflag);
1759 // Get fitted parameters
1760 Double_t xrec[2], yrec[2], qfrac;
1762 Double_t epxz, b1, b2;
1764 clusterInput.Fitter()->mnpout(0, chname, xrec[0], epxz, b1, b2, ierflg);
1765 clusterInput.Fitter()->mnpout(1, chname, yrec[0], epxz, b1, b2, ierflg);
1766 clusterInput.Fitter()->mnpout(2, chname, xrec[1], epxz, b1, b2, ierflg);
1767 clusterInput.Fitter()->mnpout(3, chname, yrec[1], epxz, b1, b2, ierflg);
1768 clusterInput.Fitter()->mnpout(4, chname, qfrac, epxz, b1, b2, ierflg);
1770 Double_t fmin, fedm, errdef;
1771 Int_t npari, nparx, istat;
1773 clusterInput.Fitter()->mnstat(fmin, fedm, errdef, npari, nparx, istat);
1778 Float_t AliMUONClusterFinderVS::CombiDoubleMathiesonFit(AliMUONRawCluster * /*c*/)
1781 // Perform combined double Mathieson fit on both cathode planes
1783 Double_t arglist[20];
1785 AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
1786 clusterInput.Fitter()->SetFCN(fcnCombiS2);
1787 clusterInput.Fitter()->mninit(6,10,7);
1788 clusterInput.Fitter()->SetPrintLevel(-1+fDebugLevel);
1790 clusterInput.Fitter()->mnexcm("SET NOW", arglist, 0, ierflag);
1791 // Set starting values
1792 static Double_t vstart[6];
1793 vstart[0]=fXInit[0];
1794 vstart[1]=fYInit[0];
1795 vstart[2]=fXInit[1];
1796 vstart[3]=fYInit[1];
1797 vstart[4]=fQrInit[0];
1798 vstart[5]=fQrInit[1];
1799 // lower and upper limits
1800 static Float_t lower[6], upper[6];
1804 fSeg2[1]->GetPadI(fInput->DetElemId(),fXInit[0], fYInit[0], fZPlane, ix, iy);
1805 isec=fSeg2[1]->Sector(fInput->DetElemId(),ix, iy);
1806 dpx=fSeg2[1]->Dpx(fInput->DetElemId(), isec);
1808 fSeg2[0]->GetPadI(fInput->DetElemId(), fXInit[0], fYInit[0], fZPlane, ix, iy);
1809 isec=fSeg2[0]->Sector(fInput->DetElemId(), ix, iy);
1810 dpy=fSeg2[0]->Dpy(fInput->DetElemId(), isec);
1815 Float_t xdum, ydum, zdum;
1816 AliDebug(1,Form("\n Cluster Finder: %f %f %f %f ", fXInit[0], fXInit[1],fYInit[0], fYInit[1] ));
1818 // Find save upper and lower limits
1821 for (fSeg2[1]->FirstPad(fInput->DetElemId(),fXInit[0], fYInit[0], fZPlane, dpx, 0.);
1822 fSeg2[1]->MorePads(fInput->DetElemId());
1823 fSeg2[1]->NextPad(fInput->DetElemId()))
1825 ix=fSeg2[1]->Ix(); iy=fSeg2[1]->Iy();
1826 // if (fHitMap[1]->TestHit(ix, iy) == kEmpty) continue;
1827 fSeg2[1]->GetPadC(fInput->DetElemId(),ix,iy,upper[0],ydum,zdum);
1828 if (icount ==0) lower[0]=upper[0];
1831 if (lower[0]>upper[0]) {xdum=lower[0]; lower[0]=upper[0]; upper[0]=xdum;}
1832 // vstart[0] = 0.5*(lower[0]+upper[0]);
1837 for (fSeg2[0]->FirstPad(fInput->DetElemId(),fXInit[0], fYInit[0], fZPlane, 0., dpy);
1838 fSeg2[0]->MorePads(fInput->DetElemId());
1839 fSeg2[0]->NextPad(fInput->DetElemId()))
1841 ix=fSeg2[0]->Ix(); iy=fSeg2[0]->Iy();
1842 // if (fHitMap[0]->TestHit(ix, iy) == kEmpty) continue;
1843 fSeg2[0]->GetPadC(fInput->DetElemId(),ix,iy,xdum,upper[1],zdum);
1844 if (icount ==0) lower[1]=upper[1];
1848 if (lower[1]>upper[1]) {xdum=lower[1]; lower[1]=upper[1]; upper[1]=xdum;}
1849 // vstart[1] = 0.5*(lower[1]+upper[1]);
1852 fSeg2[1]->GetPadI(fInput->DetElemId(),fXInit[1], fYInit[1], fZPlane, ix, iy);
1853 isec=fSeg2[1]->Sector(fInput->DetElemId(),ix, iy);
1854 dpx=fSeg2[1]->Dpx(fInput->DetElemId(),isec);
1855 fSeg2[0]->GetPadI(fInput->DetElemId(),fXInit[1], fYInit[1], fZPlane, ix, iy);
1856 isec=fSeg2[0]->Sector(fInput->DetElemId(),ix, iy);
1857 dpy=fSeg2[0]->Dpy(fInput->DetElemId(),isec);
1860 // Find save upper and lower limits
1864 for (fSeg2[1]->FirstPad(fInput->DetElemId(),fXInit[1], fYInit[1], fZPlane, dpx, 0);
1865 fSeg2[1]->MorePads(fInput->DetElemId());
1866 fSeg2[1]->NextPad(fInput->DetElemId()))
1868 ix=fSeg2[1]->Ix(); iy=fSeg2[1]->Iy();
1869 // if (fHitMap[1]->TestHit(ix, iy) == kEmpty) continue;
1870 fSeg2[1]->GetPadC(fInput->DetElemId(),ix,iy,upper[2],ydum,zdum);
1871 if (icount ==0) lower[2]=upper[2];
1874 if (lower[2]>upper[2]) {xdum=lower[2]; lower[2]=upper[2]; upper[2]=xdum;}
1875 // vstart[2] = 0.5*(lower[2]+upper[2]);
1879 for (fSeg2[0]->FirstPad(fInput->DetElemId(),fXInit[1], fYInit[1], fZPlane, 0, dpy);
1880 fSeg2[0]-> MorePads(fInput->DetElemId());
1881 fSeg2[0]->NextPad(fInput->DetElemId()))
1883 ix=fSeg2[0]->Ix(); iy=fSeg2[0]->Iy();
1884 // if (fHitMap[0]->TestHit(ix, iy) != kEmpty) continue;
1886 fSeg2[0]->GetPadC(fInput->DetElemId(),ix,iy,xdum,upper[3],zdum);
1887 if (icount ==0) lower[3]=upper[3];
1891 if (lower[3]>upper[3]) {xdum=lower[3]; lower[3]=upper[3]; upper[3]=xdum;}
1899 static Double_t step[6]={0.0005, 0.0005, 0.0005, 0.0005, 0.001, 0.001};
1900 clusterInput.Fitter()->mnparm(0,"x1",vstart[0],step[0],lower[0],upper[0],ierflag);
1901 clusterInput.Fitter()->mnparm(1,"y1",vstart[1],step[1],lower[1],upper[1],ierflag);
1902 clusterInput.Fitter()->mnparm(2,"x2",vstart[2],step[2],lower[2],upper[2],ierflag);
1903 clusterInput.Fitter()->mnparm(3,"y2",vstart[3],step[3],lower[3],upper[3],ierflag);
1904 clusterInput.Fitter()->mnparm(4,"a0",vstart[4],step[4],lower[4],upper[4],ierflag);
1905 clusterInput.Fitter()->mnparm(5,"a1",vstart[5],step[5],lower[5],upper[5],ierflag);
1906 // ready for minimisation
1910 clusterInput.Fitter()->mnexcm("SET NOGR", arglist, 0, ierflag);
1911 clusterInput.Fitter()->mnexcm("MIGRAD", arglist, 0, ierflag);
1912 // clusterInput.Fitter()->mnexcm("EXIT" , arglist, 0, ierflag);
1913 // Get fitted parameters
1915 Double_t epxz, b1, b2;
1917 clusterInput.Fitter()->mnpout(0, chname, fXFit[0], epxz, b1, b2, ierflg);
1918 clusterInput.Fitter()->mnpout(1, chname, fYFit[0], epxz, b1, b2, ierflg);
1919 clusterInput.Fitter()->mnpout(2, chname, fXFit[1], epxz, b1, b2, ierflg);
1920 clusterInput.Fitter()->mnpout(3, chname, fYFit[1], epxz, b1, b2, ierflg);
1921 clusterInput.Fitter()->mnpout(4, chname, fQrFit[0], epxz, b1, b2, ierflg);
1922 clusterInput.Fitter()->mnpout(5, chname, fQrFit[1], epxz, b1, b2, ierflg);
1924 Double_t fmin, fedm, errdef;
1925 Int_t npari, nparx, istat;
1927 clusterInput.Fitter()->mnstat(fmin, fedm, errdef, npari, nparx, istat);
1935 void AliMUONClusterFinderVS::Split(AliMUONRawCluster* c)
1938 // One cluster for each maximum
1941 AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
1942 for (j=0; j<2; j++) {
1943 AliMUONRawCluster cnew;
1944 cnew.SetGhost(c->GetGhost());
1945 for (cath=0; cath<2; cath++) {
1946 cnew.SetChi2(cath,fChi2[0]);
1947 // ?? why not cnew.fChi2[cath]=fChi2[cath];
1950 cnew.SetNcluster(0,-1);
1951 cnew.SetNcluster(1,fNRawClusters);
1953 cnew.SetNcluster(0,fNPeaks);
1954 cnew.SetNcluster(1,0);
1956 cnew.SetMultiplicity(cath,0);
1957 cnew.SetX(cath, Float_t(fXFit[j]));
1958 cnew.SetY(cath, Float_t(fYFit[j]));
1959 cnew.SetZ(cath, fZPlane);
1961 cnew.SetCharge(cath, Int_t(clusterInput.TotalCharge(cath)*fQrFit[cath]));
1963 cnew.SetCharge(cath, Int_t(clusterInput.TotalCharge(cath)*(1-fQrFit[cath])));
1965 fSeg2[cath]->SetHit(fInput->DetElemId(), fXFit[j],fYFit[j],fZPlane);
1967 for (i=0; i<fMul[cath]; i++) {
1969 cnew.SetIndex(cnew.GetMultiplicity(cath), cath, c->GetIndex(i,cath));
1971 fSeg2[cath]->SetPad(fInput->DetElemId(),fIx[i][cath], fIy[i][cath]);
1972 q1 = fInput->Mathieson()->IntXY(fInput->DetElemId(),fSeg2[cath]);
1974 cnew.SetContrib(i, cath, q1*Float_t(cnew.GetCharge(cath))/Float_t(fQ[i][cath]));
1975 cnew.SetMultiplicity(cath, cnew.GetMultiplicity(cath)+1 );
1977 FillCluster(&cnew,0,cath);
1979 cnew.SetClusterType(cnew.PhysicsContribution());
1980 if (cnew.GetCharge(0)>0 && cnew.GetCharge(1)>0) AddRawCluster(cnew);
1984 void AliMUONClusterFinderVS::AddRawCluster(AliMUONRawCluster& c)
1987 // Add a raw cluster copy to the list
1989 // AliMUON *pMUON=(AliMUON*)gAlice->GetModule("MUON");
1990 // pMUON->GetMUONData()->AddRawCluster(fInput->Chamber(),c);
1993 // Setting detection element in raw cluster for alignment
1995 c.SetDetElemId(fInput->DetElemId());
1997 TClonesArray &lrawcl = *fRawClusters;
1998 new(lrawcl[fNRawClusters++]) AliMUONRawCluster(c);
1999 AliDebug(1,Form("\nfNRawClusters %d\n",fNRawClusters));
2002 AliMUONClusterFinderVS& AliMUONClusterFinderVS
2003 ::operator = (const AliMUONClusterFinderVS& rhs)
2005 // Protected assignement operator
2007 if (this == &rhs) return *this;
2009 AliFatal("Not implemented.");
2015 // Minimisation functions
2017 void fcnS1(Int_t & /*npar*/, Double_t * /*gin*/, Double_t &f, Double_t *par, Int_t /*iflag*/)
2019 AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
2026 for (i=0; i<clusterInput.Nmul(0); i++) {
2027 Float_t q0=clusterInput.Charge(i,0);
2028 Float_t q1=clusterInput.DiscrChargeS1(i,par);
2037 void fcnCombiS1(Int_t & /*npar*/, Double_t * /*gin*/, Double_t &f, Double_t *par, Int_t /*iflag*/)
2039 AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
2046 for (cath=0; cath<2; cath++) {
2047 for (i=0; i<clusterInput.Nmul(cath); i++) {
2048 Float_t q0=clusterInput.Charge(i,cath);
2049 Float_t q1=clusterInput.DiscrChargeCombiS1(i,par,cath);
2060 void fcnS2(Int_t & /*npar*/, Double_t * /*gin*/, Double_t &f, Double_t *par, Int_t /*iflag*/)
2062 AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
2069 for (i=0; i<clusterInput.Nmul(0); i++) {
2071 Float_t q0=clusterInput.Charge(i,0);
2072 Float_t q1=clusterInput.DiscrChargeS2(i,par);
2082 void fcnCombiS2(Int_t & /*npar*/, Double_t * /*gin*/, Double_t &f, Double_t *par, Int_t /*iflag*/)
2084 AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
2090 for (cath=0; cath<2; cath++) {
2091 for (i=0; i<clusterInput.Nmul(cath); i++) {
2092 Float_t q0=clusterInput.Charge(i,cath);
2093 Float_t q1=clusterInput.DiscrChargeCombiS2(i,par,cath);