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 **************************************************************************/
18 #include "AliMUONClusterFinderVS.h"
19 #include "AliMUONDigit.h"
20 #include "AliMUONRawCluster.h"
21 #include "AliMUONGeometrySegmentation.h"
22 #include "AliMUONMathieson.h"
23 #include "AliMUONClusterInput.h"
24 #include "AliMUONDigitMapA1.h"
31 #include <Riostream.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 fSeg2[0] = fInput->Segmentation2(0);
1441 fSeg2[1] = fInput->Segmentation2(1);
1443 Int_t detElemId = fInput->DetElemId();
1445 Int_t npx0 = fSeg2[0]->Npx(detElemId)+1;
1446 Int_t npy0 = fSeg2[0]->Npy(detElemId)+1;
1447 fDigitMap[0] = new AliMUONDigitMapA1(detElemId, npx0, npy0);
1449 Int_t npx1 = fSeg2[0]->Npx(detElemId)+1;
1450 Int_t npy1 = fSeg2[0]->Npy(detElemId)+1;
1451 fDigitMap[1] = new AliMUONDigitMapA1(detElemId, npx1, npy1);
1459 fDigitMap[0]->FillHits(fInput->Digits(0));
1460 fDigitMap[1]->FillHits(fInput->Digits(1));
1462 // Outer Loop over Cathodes
1463 for (cath = 0; cath < 2; cath++) {
1465 for (ndig=0; ndig<fInput->NDigits(cath); ndig++) {
1466 dig = fInput->Digit(cath, ndig);
1467 Int_t padx = dig->PadX();
1468 Int_t pady = dig->PadY();
1469 if (fDigitMap[cath]->TestHit(padx,pady)==kUsed ||fDigitMap[0]->TestHit(padx,pady)==kEmpty) {
1473 AliDebug(1,Form("\n CATHODE %d CLUSTER %d\n",cath,ncls));
1474 AliMUONRawCluster clus;
1475 clus.SetMultiplicity(0, 0);
1476 clus.SetMultiplicity(1, 0);
1477 clus.SetPeakSignal(cath,dig->Signal());
1478 clus.SetTrack(0, dig->Hit());
1479 clus.SetTrack(1, dig->Track(0));
1480 clus.SetTrack(2, dig->Track(1));
1482 AliDebug(1,Form("idDE %d Padx %d Pady %d", fInput->DetElemId(), padx, pady));
1484 // tag the beginning of cluster list in a raw cluster
1485 clus.SetNcluster(0,-1);
1487 fSeg2[cath]->GetPadC(fInput->DetElemId(), padx, pady, xcu, ycu, fZPlane);
1488 fSector= fSeg2[cath]->Sector(fInput->DetElemId(), padx, pady)/100;
1493 FindCluster(padx,pady,cath,clus);
1494 //^^^^^^^^^^^^^^^^^^^^^^^^
1495 // center of gravity
1496 if (clus.GetX(0)!=0.) clus.SetX(0, clus.GetX(0)/clus.GetCharge(0)); // clus.fX[0] /= clus.fQ[0];
1499 clus.SetX(0,fSeg2[0]->GetAnod(fInput->DetElemId(), clus.GetX(0)));
1500 if (clus.GetY(0)!=0.) clus.SetY(0, clus.GetY(0)/clus.GetCharge(0)); // clus.fY[0] /= clus.fQ[0];
1502 if(clus.GetCharge(1)!=0.) clus.SetX(1, clus.GetX(1)/clus.GetCharge(1)); // clus.fX[1] /= clus.fQ[1];
1505 clus.SetX(1, fSeg2[0]->GetAnod(fInput->DetElemId(),clus.GetX(1)));
1506 if(clus.GetCharge(1)!=0.) clus.SetY(1, clus.GetY(1)/clus.GetCharge(1));// clus.fY[1] /= clus.fQ[1];
1508 clus.SetZ(0, fZPlane);
1509 clus.SetZ(1, fZPlane);
1511 AliDebug(1,Form("\n Cathode 1 multiplicite %d X(CG) %f Y(CG) %f\n",
1512 clus.GetMultiplicity(0),clus.GetX(0),clus.GetY(0)));
1513 AliDebug(1,Form(" Cathode 2 multiplicite %d X(CG) %f Y(CG) %f\n",
1514 clus.GetMultiplicity(1),clus.GetX(1),clus.GetY(1)));
1515 // Analyse cluster and decluster if necessary
1518 clus.SetNcluster(1,fNRawClusters);
1519 clus.SetClusterType(clus.PhysicsContribution());
1526 // reset Cluster object
1527 { // begin local scope
1528 for (int k=0;k<clus.GetMultiplicity(0);k++) clus.SetIndex(k, 0, 0);
1529 } // end local scope
1531 { // begin local scope
1532 for (int k=0;k<clus.GetMultiplicity(1);k++) clus.SetIndex(k, 1, 0);
1533 } // end local scope
1535 clus.SetMultiplicity(0,0);
1536 clus.SetMultiplicity(1,0);
1540 } // end loop cathodes
1541 delete fDigitMap[0];
1542 delete fDigitMap[1];
1545 Float_t AliMUONClusterFinderVS::SingleMathiesonFit(AliMUONRawCluster *c, Int_t cath)
1547 // Performs a single Mathieson fit on one cathode
1549 Double_t arglist[20];
1551 AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
1553 clusterInput.Fitter()->SetFCN(fcnS1);
1554 clusterInput.Fitter()->mninit(2,10,7);
1555 clusterInput.Fitter()->SetPrintLevel(-1 + AliLog::GetGlobalDebugLevel());
1557 clusterInput.Fitter()->mnexcm("SET NOW", arglist, 0, ierflag);
1558 // Set starting values
1559 static Double_t vstart[2];
1560 vstart[0]=c->GetX(1);
1561 vstart[1]=c->GetY(0);
1564 // lower and upper limits
1565 static Double_t lower[2], upper[2];
1567 fSeg2[cath]->GetPadI(fInput->DetElemId(), c->GetX(cath), c->GetY(cath), fZPlane, ix, iy);
1568 isec=fSeg2[cath]->Sector(fInput->DetElemId(), ix, iy);
1570 lower[0]=vstart[0]-fSeg2[cath]->Dpx(fInput->DetElemId(), isec)/2;
1571 lower[1]=vstart[1]-fSeg2[cath]->Dpy(fInput->DetElemId(), isec)/2;
1573 upper[0]=lower[0]+fSeg2[cath]->Dpx(fInput->DetElemId(), isec);
1574 upper[1]=lower[1]+fSeg2[cath]->Dpy(fInput->DetElemId(), isec);
1578 static Double_t step[2]={0.0005, 0.0005};
1580 clusterInput.Fitter()->mnparm(0,"x1",vstart[0],step[0],lower[0],upper[0],ierflag);
1581 clusterInput.Fitter()->mnparm(1,"y1",vstart[1],step[1],lower[1],upper[1],ierflag);
1582 // ready for minimisation
1586 clusterInput.Fitter()->mnexcm("SET NOGR", arglist, 0, ierflag);
1587 clusterInput.Fitter()->mnexcm("MIGRAD", arglist, 0, ierflag);
1588 // clusterInput.Fitter()->mnexcm("EXIT" , arglist, 0, ierflag);
1589 Double_t fmin, fedm, errdef;
1590 Int_t npari, nparx, istat;
1592 clusterInput.Fitter()->mnstat(fmin, fedm, errdef, npari, nparx, istat);
1596 // Get fitted parameters
1597 Double_t xrec, yrec;
1599 Double_t epxz, b1, b2;
1601 clusterInput.Fitter()->mnpout(0, chname, xrec, epxz, b1, b2, ierflg);
1602 clusterInput.Fitter()->mnpout(1, chname, yrec, epxz, b1, b2, ierflg);
1608 Float_t AliMUONClusterFinderVS::CombiSingleMathiesonFit(AliMUONRawCluster * /*c*/)
1610 // Perform combined Mathieson fit on both cathode planes
1612 Double_t arglist[20];
1614 AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
1615 clusterInput.Fitter()->SetFCN(fcnCombiS1);
1616 clusterInput.Fitter()->mninit(2,10,7);
1617 clusterInput.Fitter()->SetPrintLevel(-1 + AliLog::GetGlobalDebugLevel());
1619 clusterInput.Fitter()->mnexcm("SET NOW", arglist, 0, ierflag);
1620 static Double_t vstart[2];
1621 vstart[0]=fXInit[0];
1622 vstart[1]=fYInit[0];
1625 // lower and upper limits
1626 static Float_t lower[2], upper[2];
1630 fSeg2[0]->GetPadI(fInput->DetElemId(), fXInit[0], fYInit[0], fZPlane, ix, iy);
1631 isec=fSeg2[0]->Sector(fInput->DetElemId(), ix, iy);
1632 dpy=fSeg2[0]->Dpy(fInput->DetElemId(), isec);
1633 fSeg2[1]->GetPadI(fInput->DetElemId(), fXInit[0], fYInit[0], fZPlane, ix, iy);
1634 isec=fSeg2[1]->Sector(fInput->DetElemId(), ix, iy);
1635 dpx=fSeg2[1]->Dpx(fInput->DetElemId(), isec);
1638 Float_t xdum, ydum, zdum;
1640 // Find save upper and lower limits
1643 for (fSeg2[1]->FirstPad(fInput->DetElemId(),fXInit[0], fYInit[0], fZPlane, dpx, 0.);
1644 fSeg2[1]->MorePads(fInput->DetElemId());
1645 fSeg2[1]->NextPad(fInput->DetElemId()))
1647 ix=fSeg2[1]->Ix(); iy=fSeg2[1]->Iy();
1648 fSeg2[1]->GetPadC(fInput->DetElemId(), ix,iy, upper[0], ydum, zdum);
1649 if (icount ==0) lower[0]=upper[0];
1653 if (lower[0]>upper[0]) {xdum=lower[0]; lower[0]=upper[0]; upper[0]=xdum;}
1656 AliDebug(1,Form("\n single y %f %f", fXInit[0], fYInit[0]));
1658 for (fSeg2[0]->FirstPad(fInput->DetElemId(), fXInit[0], fYInit[0], fZPlane, 0., dpy);
1659 fSeg2[0]->MorePads(fInput->DetElemId());
1660 fSeg2[0]->NextPad(fInput->DetElemId()))
1662 ix=fSeg2[0]->Ix(); iy=fSeg2[0]->Iy();
1663 fSeg2[0]->GetPadC(fInput->DetElemId(), ix,iy,xdum,upper[1],zdum);
1664 if (icount ==0) lower[1]=upper[1];
1666 AliDebug(1,Form("\n upper lower %d %f %f", icount, upper[1], lower[1]));
1669 if (lower[1]>upper[1]) {xdum=lower[1]; lower[1]=upper[1]; upper[1]=xdum;}
1672 static Double_t step[2]={0.00001, 0.0001};
1674 clusterInput.Fitter()->mnparm(0,"x1",vstart[0],step[0],lower[0],upper[0],ierflag);
1675 clusterInput.Fitter()->mnparm(1,"y1",vstart[1],step[1],lower[1],upper[1],ierflag);
1676 // ready for minimisation
1680 clusterInput.Fitter()->mnexcm("SET NOGR", arglist, 0, ierflag);
1681 clusterInput.Fitter()->mnexcm("MIGRAD", arglist, 0, ierflag);
1682 // clusterInput.Fitter()->mnexcm("EXIT" , arglist, 0, ierflag);
1683 Double_t fmin, fedm, errdef;
1684 Int_t npari, nparx, istat;
1686 clusterInput.Fitter()->mnstat(fmin, fedm, errdef, npari, nparx, istat);
1690 // Get fitted parameters
1691 Double_t xrec, yrec;
1693 Double_t epxz, b1, b2;
1695 clusterInput.Fitter()->mnpout(0, chname, xrec, epxz, b1, b2, ierflg);
1696 clusterInput.Fitter()->mnpout(1, chname, yrec, epxz, b1, b2, ierflg);
1702 Bool_t AliMUONClusterFinderVS::DoubleMathiesonFit(AliMUONRawCluster * /*c*/, Int_t cath)
1704 // Performs a double Mathieson fit on one cathode
1708 // Initialise global variables for fit
1709 Double_t arglist[20];
1711 AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
1712 clusterInput.Fitter()->SetFCN(fcnS2);
1713 clusterInput.Fitter()->mninit(5,10,7);
1714 clusterInput.Fitter()->SetPrintLevel(-1 + AliLog::GetGlobalDebugLevel());
1716 clusterInput.Fitter()->mnexcm("SET NOW", arglist, 0, ierflag);
1717 // Set starting values
1718 static Double_t vstart[5];
1719 vstart[0]=fX[fIndLocal[0][cath]][cath];
1720 vstart[1]=fY[fIndLocal[0][cath]][cath];
1721 vstart[2]=fX[fIndLocal[1][cath]][cath];
1722 vstart[3]=fY[fIndLocal[1][cath]][cath];
1723 vstart[4]=Float_t(fQ[fIndLocal[0][cath]][cath])/
1724 Float_t(fQ[fIndLocal[0][cath]][cath]+fQ[fIndLocal[1][cath]][cath]);
1725 // lower and upper limits
1726 static Float_t lower[5], upper[5];
1729 isec=fSeg2[cath]->Sector(fInput->DetElemId(),fIx[fIndLocal[0][cath]][cath],
1730 fIy[fIndLocal[0][cath]][cath]);
1731 lower[0]=vstart[0]-fSeg2[cath]->Dpx(fInput->DetElemId(),isec);
1732 lower[1]=vstart[1]-fSeg2[cath]->Dpy(fInput->DetElemId(),isec);
1734 upper[0]=lower[0]+2.*fSeg2[cath]->Dpx(fInput->DetElemId(),isec);
1735 upper[1]=lower[1]+2.*fSeg2[cath]->Dpy(fInput->DetElemId(),isec);
1737 isec=fSeg2[cath]->Sector(fInput->DetElemId(),fIx[fIndLocal[1][cath]][cath],
1738 fIy[fIndLocal[1][cath]][cath]);
1739 lower[2]=vstart[2]-fSeg2[cath]->Dpx(fInput->DetElemId(),isec)/2;
1740 lower[3]=vstart[3]-fSeg2[cath]->Dpy(fInput->DetElemId(),isec)/2;
1742 upper[2]=lower[2]+fSeg2[cath]->Dpx(fInput->DetElemId(),isec);
1743 upper[1]=lower[1]+2.*fSeg2[cath]->Dpy(fInput->DetElemId(),isec);
1750 static Double_t step[5]={0.0005, 0.0005, 0.0005, 0.0005, 0.0001};
1752 clusterInput.Fitter()->mnparm(0,"x1",vstart[0],step[0],lower[0],upper[0],ierflag);
1753 clusterInput.Fitter()->mnparm(1,"y1",vstart[1],step[1],lower[1],upper[1],ierflag);
1754 clusterInput.Fitter()->mnparm(2,"x2",vstart[2],step[2],lower[2],upper[2],ierflag);
1755 clusterInput.Fitter()->mnparm(3,"y2",vstart[3],step[3],lower[3],upper[3],ierflag);
1756 clusterInput.Fitter()->mnparm(4,"a0",vstart[4],step[4],lower[4],upper[4],ierflag);
1757 // ready for minimisation
1761 clusterInput.Fitter()->mnexcm("SET NOGR", arglist, 0, ierflag);
1762 clusterInput.Fitter()->mnexcm("MIGRAD", arglist, 0, ierflag);
1763 // clusterInput.Fitter()->mnexcm("EXIT" , arglist, 0, ierflag);
1764 // Get fitted parameters
1765 Double_t xrec[2], yrec[2], qfrac;
1767 Double_t epxz, b1, b2;
1769 clusterInput.Fitter()->mnpout(0, chname, xrec[0], epxz, b1, b2, ierflg);
1770 clusterInput.Fitter()->mnpout(1, chname, yrec[0], epxz, b1, b2, ierflg);
1771 clusterInput.Fitter()->mnpout(2, chname, xrec[1], epxz, b1, b2, ierflg);
1772 clusterInput.Fitter()->mnpout(3, chname, yrec[1], epxz, b1, b2, ierflg);
1773 clusterInput.Fitter()->mnpout(4, chname, qfrac, epxz, b1, b2, ierflg);
1775 Double_t fmin, fedm, errdef;
1776 Int_t npari, nparx, istat;
1778 clusterInput.Fitter()->mnstat(fmin, fedm, errdef, npari, nparx, istat);
1783 Float_t AliMUONClusterFinderVS::CombiDoubleMathiesonFit(AliMUONRawCluster * /*c*/)
1786 // Perform combined double Mathieson fit on both cathode planes
1788 Double_t arglist[20];
1790 AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
1791 clusterInput.Fitter()->SetFCN(fcnCombiS2);
1792 clusterInput.Fitter()->mninit(6,10,7);
1793 clusterInput.Fitter()->SetPrintLevel(-1 + AliLog::GetGlobalDebugLevel());
1795 clusterInput.Fitter()->mnexcm("SET NOW", arglist, 0, ierflag);
1796 // Set starting values
1797 static Double_t vstart[6];
1798 vstart[0]=fXInit[0];
1799 vstart[1]=fYInit[0];
1800 vstart[2]=fXInit[1];
1801 vstart[3]=fYInit[1];
1802 vstart[4]=fQrInit[0];
1803 vstart[5]=fQrInit[1];
1804 // lower and upper limits
1805 static Float_t lower[6], upper[6];
1809 fSeg2[1]->GetPadI(fInput->DetElemId(),fXInit[0], fYInit[0], fZPlane, ix, iy);
1810 isec=fSeg2[1]->Sector(fInput->DetElemId(),ix, iy);
1811 dpx=fSeg2[1]->Dpx(fInput->DetElemId(), isec);
1813 fSeg2[0]->GetPadI(fInput->DetElemId(), fXInit[0], fYInit[0], fZPlane, ix, iy);
1814 isec=fSeg2[0]->Sector(fInput->DetElemId(), ix, iy);
1815 dpy=fSeg2[0]->Dpy(fInput->DetElemId(), isec);
1820 Float_t xdum, ydum, zdum;
1821 AliDebug(1,Form("\n Cluster Finder: %f %f %f %f ", fXInit[0], fXInit[1],fYInit[0], fYInit[1] ));
1823 // Find save upper and lower limits
1826 for (fSeg2[1]->FirstPad(fInput->DetElemId(),fXInit[0], fYInit[0], fZPlane, dpx, 0.);
1827 fSeg2[1]->MorePads(fInput->DetElemId());
1828 fSeg2[1]->NextPad(fInput->DetElemId()))
1830 ix=fSeg2[1]->Ix(); iy=fSeg2[1]->Iy();
1831 // if (fDigitMap[1]->TestHit(ix, iy) == kEmpty) continue;
1832 fSeg2[1]->GetPadC(fInput->DetElemId(),ix,iy,upper[0],ydum,zdum);
1833 if (icount ==0) lower[0]=upper[0];
1836 if (lower[0]>upper[0]) {xdum=lower[0]; lower[0]=upper[0]; upper[0]=xdum;}
1837 // vstart[0] = 0.5*(lower[0]+upper[0]);
1842 for (fSeg2[0]->FirstPad(fInput->DetElemId(),fXInit[0], fYInit[0], fZPlane, 0., dpy);
1843 fSeg2[0]->MorePads(fInput->DetElemId());
1844 fSeg2[0]->NextPad(fInput->DetElemId()))
1846 ix=fSeg2[0]->Ix(); iy=fSeg2[0]->Iy();
1847 // if (fDigitMap[0]->TestHit(ix, iy) == kEmpty) continue;
1848 fSeg2[0]->GetPadC(fInput->DetElemId(),ix,iy,xdum,upper[1],zdum);
1849 if (icount ==0) lower[1]=upper[1];
1853 if (lower[1]>upper[1]) {xdum=lower[1]; lower[1]=upper[1]; upper[1]=xdum;}
1854 // vstart[1] = 0.5*(lower[1]+upper[1]);
1857 fSeg2[1]->GetPadI(fInput->DetElemId(),fXInit[1], fYInit[1], fZPlane, ix, iy);
1858 isec=fSeg2[1]->Sector(fInput->DetElemId(),ix, iy);
1859 dpx=fSeg2[1]->Dpx(fInput->DetElemId(),isec);
1860 fSeg2[0]->GetPadI(fInput->DetElemId(),fXInit[1], fYInit[1], fZPlane, ix, iy);
1861 isec=fSeg2[0]->Sector(fInput->DetElemId(),ix, iy);
1862 dpy=fSeg2[0]->Dpy(fInput->DetElemId(),isec);
1865 // Find save upper and lower limits
1869 for (fSeg2[1]->FirstPad(fInput->DetElemId(),fXInit[1], fYInit[1], fZPlane, dpx, 0);
1870 fSeg2[1]->MorePads(fInput->DetElemId());
1871 fSeg2[1]->NextPad(fInput->DetElemId()))
1873 ix=fSeg2[1]->Ix(); iy=fSeg2[1]->Iy();
1874 // if (fDigitMap[1]->TestHit(ix, iy) == kEmpty) continue;
1875 fSeg2[1]->GetPadC(fInput->DetElemId(),ix,iy,upper[2],ydum,zdum);
1876 if (icount ==0) lower[2]=upper[2];
1879 if (lower[2]>upper[2]) {xdum=lower[2]; lower[2]=upper[2]; upper[2]=xdum;}
1880 // vstart[2] = 0.5*(lower[2]+upper[2]);
1884 for (fSeg2[0]->FirstPad(fInput->DetElemId(),fXInit[1], fYInit[1], fZPlane, 0, dpy);
1885 fSeg2[0]-> MorePads(fInput->DetElemId());
1886 fSeg2[0]->NextPad(fInput->DetElemId()))
1888 ix=fSeg2[0]->Ix(); iy=fSeg2[0]->Iy();
1889 // if (fDigitMap[0]->TestHit(ix, iy) != kEmpty) continue;
1891 fSeg2[0]->GetPadC(fInput->DetElemId(),ix,iy,xdum,upper[3],zdum);
1892 if (icount ==0) lower[3]=upper[3];
1896 if (lower[3]>upper[3]) {xdum=lower[3]; lower[3]=upper[3]; upper[3]=xdum;}
1904 static Double_t step[6]={0.0005, 0.0005, 0.0005, 0.0005, 0.001, 0.001};
1905 clusterInput.Fitter()->mnparm(0,"x1",vstart[0],step[0],lower[0],upper[0],ierflag);
1906 clusterInput.Fitter()->mnparm(1,"y1",vstart[1],step[1],lower[1],upper[1],ierflag);
1907 clusterInput.Fitter()->mnparm(2,"x2",vstart[2],step[2],lower[2],upper[2],ierflag);
1908 clusterInput.Fitter()->mnparm(3,"y2",vstart[3],step[3],lower[3],upper[3],ierflag);
1909 clusterInput.Fitter()->mnparm(4,"a0",vstart[4],step[4],lower[4],upper[4],ierflag);
1910 clusterInput.Fitter()->mnparm(5,"a1",vstart[5],step[5],lower[5],upper[5],ierflag);
1911 // ready for minimisation
1915 clusterInput.Fitter()->mnexcm("SET NOGR", arglist, 0, ierflag);
1916 clusterInput.Fitter()->mnexcm("MIGRAD", arglist, 0, ierflag);
1917 // clusterInput.Fitter()->mnexcm("EXIT" , arglist, 0, ierflag);
1918 // Get fitted parameters
1920 Double_t epxz, b1, b2;
1922 clusterInput.Fitter()->mnpout(0, chname, fXFit[0], epxz, b1, b2, ierflg);
1923 clusterInput.Fitter()->mnpout(1, chname, fYFit[0], epxz, b1, b2, ierflg);
1924 clusterInput.Fitter()->mnpout(2, chname, fXFit[1], epxz, b1, b2, ierflg);
1925 clusterInput.Fitter()->mnpout(3, chname, fYFit[1], epxz, b1, b2, ierflg);
1926 clusterInput.Fitter()->mnpout(4, chname, fQrFit[0], epxz, b1, b2, ierflg);
1927 clusterInput.Fitter()->mnpout(5, chname, fQrFit[1], epxz, b1, b2, ierflg);
1929 Double_t fmin, fedm, errdef;
1930 Int_t npari, nparx, istat;
1932 clusterInput.Fitter()->mnstat(fmin, fedm, errdef, npari, nparx, istat);
1940 void AliMUONClusterFinderVS::Split(AliMUONRawCluster* c)
1943 // One cluster for each maximum
1946 AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
1947 for (j=0; j<2; j++) {
1948 AliMUONRawCluster cnew;
1949 cnew.SetGhost(c->GetGhost());
1950 for (cath=0; cath<2; cath++) {
1951 cnew.SetChi2(cath,fChi2[0]);
1952 // ?? why not cnew.fChi2[cath]=fChi2[cath];
1955 cnew.SetNcluster(0,-1);
1956 cnew.SetNcluster(1,fNRawClusters);
1958 cnew.SetNcluster(0,fNPeaks);
1959 cnew.SetNcluster(1,0);
1961 cnew.SetMultiplicity(cath,0);
1962 cnew.SetX(cath, Float_t(fXFit[j]));
1963 cnew.SetY(cath, Float_t(fYFit[j]));
1964 cnew.SetZ(cath, fZPlane);
1966 cnew.SetCharge(cath, Int_t(clusterInput.TotalCharge(cath)*fQrFit[cath]));
1968 cnew.SetCharge(cath, Int_t(clusterInput.TotalCharge(cath)*(1-fQrFit[cath])));
1970 fSeg2[cath]->SetHit(fInput->DetElemId(), fXFit[j],fYFit[j],fZPlane);
1972 for (i=0; i<fMul[cath]; i++) {
1974 cnew.SetIndex(cnew.GetMultiplicity(cath), cath, c->GetIndex(i,cath));
1976 fSeg2[cath]->SetPad(fInput->DetElemId(),fIx[i][cath], fIy[i][cath]);
1977 q1 = fInput->Mathieson()->IntXY(fInput->DetElemId(),fSeg2[cath]);
1979 cnew.SetContrib(i, cath, q1*Float_t(cnew.GetCharge(cath))/Float_t(fQ[i][cath]));
1980 cnew.SetMultiplicity(cath, cnew.GetMultiplicity(cath)+1 );
1982 FillCluster(&cnew,0,cath);
1984 cnew.SetClusterType(cnew.PhysicsContribution());
1985 if (cnew.GetCharge(0)>0 && cnew.GetCharge(1)>0) AddRawCluster(cnew);
1989 void AliMUONClusterFinderVS::AddRawCluster(AliMUONRawCluster& c)
1992 // Add a raw cluster copy to the list
1994 // AliMUON *pMUON=(AliMUON*)gAlice->GetModule("MUON");
1995 // pMUON->GetMUONData()->AddRawCluster(fInput->Chamber(),c);
1998 // Setting detection element in raw cluster for alignment
2000 c.SetDetElemId(fInput->DetElemId());
2002 TClonesArray &lrawcl = *fRawClusters;
2003 new(lrawcl[fNRawClusters++]) AliMUONRawCluster(c);
2004 AliDebug(1,Form("\nfNRawClusters %d\n",fNRawClusters));
2007 AliMUONClusterFinderVS& AliMUONClusterFinderVS
2008 ::operator = (const AliMUONClusterFinderVS& rhs)
2010 // Protected assignement operator
2012 if (this == &rhs) return *this;
2014 AliFatal("Not implemented.");
2020 // Minimisation functions
2022 void fcnS1(Int_t & /*npar*/, Double_t * /*gin*/, Double_t &f, Double_t *par, Int_t /*iflag*/)
2024 AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
2031 for (i=0; i<clusterInput.Nmul(0); i++) {
2032 Float_t q0=clusterInput.Charge(i,0);
2033 Float_t q1=clusterInput.DiscrChargeS1(i,par);
2042 void fcnCombiS1(Int_t & /*npar*/, Double_t * /*gin*/, Double_t &f, Double_t *par, Int_t /*iflag*/)
2044 AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
2051 for (cath=0; cath<2; cath++) {
2052 for (i=0; i<clusterInput.Nmul(cath); i++) {
2053 Float_t q0=clusterInput.Charge(i,cath);
2054 Float_t q1=clusterInput.DiscrChargeCombiS1(i,par,cath);
2065 void fcnS2(Int_t & /*npar*/, Double_t * /*gin*/, Double_t &f, Double_t *par, Int_t /*iflag*/)
2067 AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
2074 for (i=0; i<clusterInput.Nmul(0); i++) {
2076 Float_t q0=clusterInput.Charge(i,0);
2077 Float_t q1=clusterInput.DiscrChargeS2(i,par);
2087 void fcnCombiS2(Int_t & /*npar*/, Double_t * /*gin*/, Double_t &f, Double_t *par, Int_t /*iflag*/)
2089 AliMUONClusterInput& clusterInput = *(AliMUONClusterInput::Instance());
2095 for (cath=0; cath<2; cath++) {
2096 for (i=0; i<clusterInput.Nmul(cath); i++) {
2097 Float_t q0=clusterInput.Charge(i,cath);
2098 Float_t q1=clusterInput.DiscrChargeCombiS2(i,par,cath);